source: issm/trunk/src/m/consistency/checkfield.py@ 26744

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 3 years ago

merged trunk-jpl and trunk for revision 26742
File size: 14.0 KB
Line 
1import numpy as np
2import os
3from re import findall, split
4from pairoptions import pairoptions
5from operator import attrgetter
6import MatlabFuncs as m
7
8
9def checkfield(md, *args):
10 """CHECKFIELD - check field consistency
11
12 Used to check model consistency
13 Requires:
14 'field' or 'fieldname' option. If 'fieldname' is provided, it will
15 retrieve it from the model md. (md.(fieldname))
16 If 'field' is provided, it will assume the argument following 'field'
17 is a numeric array.
18
19 Available options:
20 - NaN: 1 if check that there is no NaN
21 - size: [lines cols], NaN for non checked dimensions, or 'universal'
22 for any input type (nodal, element, time series, etc)
23 - > : greater than provided value
24 - >= : greater or equal to provided value
25 - < : smallerthan provided value
26 - <=: smaller or equal to provided value
27 - < vec: smallerthan provided values on each vertex
28 - timeseries: 1 if check time series consistency (size and time)
29 - values: list of acceptable values
30 - numel: list of acceptable number of elements
31 - cell: 1 if check that is cell
32 - empty: 1 if check that non empty
33 - message: overloaded error message
34
35 Usage:
36 md = checkfield(md, fieldname, options)
37 """
38
39 #get options
40 options = pairoptions(*args)
41
42 #get field from model
43 if options.exist('field'):
44 field = options.getfieldvalue('field')
45 fieldname = options.getfieldvalue('fieldname', 'no fieldname')
46 else:
47 fieldname = options.getfieldvalue('fieldname')
48 fieldprefix = split(r'\[(.*?)\]', fieldname)[0]
49 fieldindexes = findall(r'\[(.*?)\]', fieldname)
50 field = attrgetter(fieldprefix)(md)
51 for index in fieldindexes:
52 try:
53 field = field[index.strip("\'")]
54 except TypeError:
55 field = field[int(index)] #looking for an index and not a key
56
57 # that works for py2
58 # exec("field = md.{}".format(fieldname))
59 # exec("field = md.{}".format(fieldname), namespace)
60
61 if isinstance(field, (bool, int, float)):
62 field = np.array([field])
63
64 #check empty
65 if options.exist('empty'):
66 if not field:
67 md = md.checkmessage(options.getfieldvalue('message', "field '{}' is empty".format(fieldname)))
68
69 #Check size
70 if options.exist('size'):
71 fieldsize = options.getfieldvalue('size')
72 if type(fieldsize) == str:
73 if fieldsize == 'universal':
74 #Check that vector size will not be confusing for ModelProcessorx
75 if (md.mesh.numberofvertices == md.mesh.numberofelements):
76 raise Exception('number of vertices is the same as number of elements')
77 elif (md.mesh.numberofvertices + 1 == md.mesh.numberofelements):
78 raise Exception('number of vertices + 1 is the same as number of elements')
79 elif (md.mesh.numberofvertices == md.mesh.numberofelements + 1):
80 raise Exception('number of vertices is the same as number of elements + 1')
81
82 #Uniform field
83 if (field.shape[0] == 1):
84 if (np.ndim(field) > 1 and np.shape(field)[1] != 1):
85 md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))
86
87 #vertex oriented input, only one column allowed
88 elif (np.shape(field)[0] == md.mesh.numberofvertices):
89 if (np.ndim(field) > 1 and np.shape(field)[1] != 1):
90 md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))
91
92 #element oriented input, one or more column (patch) is ok
93 elif (np.shape(field)[0] == md.mesh.numberofelements):
94 #nothing to do here (either constant per element, or defined on nodes)
95 pass
96
97 #vertex time series
98 elif (np.shape(field)[0] == md.mesh.numberofvertices + 1):
99 if (np.ndim(field) > 1 and np.shape(field)[1] <= 1):
100 md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))
101
102 #element time series
103 elif (np.shape(field)[0] == md.mesh.numberofelements + 1):
104 if (np.ndim(field) > 1 and np.shape(field)[1] <= 1):
105 md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))
106
107 #else not supported
108 else:
109 md = md.checkmessage(options.getfieldvalue('message', "field '{}' is not supported".format(fieldname)))
110
111 else:
112 raise Exception("fieldsize '{}' not supported yet".format(fieldsize))
113
114 else:
115 if len(np.shape(field)) < len(fieldsize):
116 md = md.checkmessage(options.getfieldvalue('message', "field {} has size {} but should be size {}".format(fieldname, np.shape(field), fieldsize)))
117 else:
118 for i in range(np.size(fieldsize)):
119 if (not np.isnan(fieldsize[i])) and (np.shape(field)[i] != fieldsize[i]):
120 md = md.checkmessage(options.getfieldvalue('message', "field {} dimension # {} should be of size {}".format(fieldname, i, fieldsize[i])))
121
122 #Check numel
123 if options.exist('numel'):
124 fieldnumel = options.getfieldvalue('numel')
125 if (type(fieldnumel) == int and np.size(field) != fieldnumel) or (type(fieldnumel) == list and np.size(field) not in fieldnumel):
126 if len(fieldnumel) == 1:
127 md = md.checkmessage(options.getfieldvalue('message', "field '{}' size should be {}".format(fieldname, fieldnumel)))
128 elif len(fieldnumel) == 2:
129 md = md.checkmessage(options.getfieldvalue('message', "field '{}' size should be {} or {}".format(fieldname, fieldnumel[0], fieldnumel[1])))
130 else:
131 md = md.checkmessage(options.getfieldvalue('message', "field '{}' size should be {}".format(fieldname, fieldnumel)))
132
133 #check NaN
134 if options.getfieldvalue('NaN', 0):
135 if np.any(np.isnan(field)):
136 md = md.checkmessage(options.getfieldvalue('message', "NaN values found in field '{}'".format(fieldname)))
137
138 #check Inf
139 if options.getfieldvalue('Inf', 0):
140 if np.any(np.isinf(field)):
141 md = md.checkmessage(options.getfieldvalue('message', "Inf values found in field '{}'".format(fieldname)))
142
143 #check cell
144 if options.getfieldvalue('cell', 0):
145 if not isinstance(field, (tuple, list, dict)):
146 md = md.checkmessage(options.getfieldvalue('message', "field '{}' should be a tuple, list, or dict".format(fieldname)))
147
148 #check values
149 if options.exist('values'):
150 fieldvalues = options.getfieldvalue('values')
151 if False in m.ismember(field, fieldvalues):
152 if len(fieldvalues) == 1:
153 md = md.checkmessage(options.getfieldvalue('message', "field '{}' value should be '{}'".format(fieldname, fieldvalues[0])))
154 elif len(fieldvalues) == 2:
155 md = md.checkmessage(options.getfieldvalue('message', "field '{}' values should be '{}' '{}' or '{}'".format(fieldname, fieldvalues[0], fieldvalues[1])))
156 else:
157 md = md.checkmessage(options.getfieldvalue('message', "field '{}' should have values in {}".format(fieldname, fieldvalues)))
158
159 #check greater
160 if options.exist('>='):
161 lowerbound = options.getfieldvalue('>=')
162 if type(lowerbound) is str:
163 lowerbound = attrgetter(lowerbound)(md)
164 if np.size(lowerbound) > 1: #checking elementwise
165 if any(field < lowerbound):
166 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
167 else:
168 minval = np.nanmin(field)
169 if options.getfieldvalue('timeseries', 0):
170 minval = np.nanmin(field[:-1])
171 elif options.getfieldvalue('singletimeseries', 0):
172 if np.size(field) == 1: #some singletimeseries are just one value
173 minval = field
174 else:
175 minval = np.nanmin(field[0])
176
177 if minval < lowerbound:
178 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
179
180 if options.exist('>'):
181 lowerbound = options.getfieldvalue('>')
182 if type(lowerbound) is str:
183 lowerbound = attrgetter(lowerbound)(md)
184 if np.size(lowerbound) > 1: #checking elementwise
185 if any(field <= lowerbound):
186 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
187 else:
188 minval = np.nanmin(field)
189 if options.getfieldvalue('timeseries', 0):
190 minval = np.nanmin(field[:-1])
191 elif options.getfieldvalue('singletimeseries', 0):
192 if np.size(field) == 1: #some singletimeseries are just one value
193 minval = field
194 else:
195 minval = np.nanmin(field[0])
196
197 if minval <= lowerbound:
198 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values above {}".format(fieldname, lowerbound)))
199
200 #check smaller
201 if options.exist('<='):
202 upperbound = options.getfieldvalue('<=')
203 if type(upperbound) is str:
204 upperbound = attrgetter(upperbound)(md)
205 if np.size(upperbound) > 1: #checking elementwise
206 if any(field > upperbound):
207 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
208 else:
209 maxval = np.nanmax(field)
210 if options.getfieldvalue('timeseries', 0):
211 maxval = np.nanmax(field[:-1])
212 elif options.getfieldvalue('singletimeseries', 0):
213 if np.size(field) == 1: #some singletimeseries are just one value
214 maxval = field
215 else:
216 maxval = np.nanmax(field[0])
217 elif hasattr(field, 'fov_forward_indices'):
218 maxval = field.fov_forward_indices[0]
219 if maxval > upperbound:
220 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
221
222 if options.exist('<'):
223 upperbound = options.getfieldvalue('<')
224 if type(upperbound) is str:
225 upperbound = attrgetter(upperbound)(md)
226 if np.size(upperbound) > 1: #checking elementwise
227 if any(field >= upperbound):
228 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
229
230 else:
231 maxval = np.nanmax(field)
232 if options.getfieldvalue('timeseries', 0):
233 maxval = np.nanmax(field[:-1])
234 elif options.getfieldvalue('singletimeseries', 0):
235 if np.size(field) == 1: #some singletimeseries are just one value
236 maxval = field.copy()
237 else:
238 maxval = np.nanmax(field[0])
239
240 if maxval >= upperbound:
241 md = md.checkmessage(options.getfieldvalue('message', "field {} should have values below {}".format(fieldname, upperbound)))
242
243 #check file
244 if options.getfieldvalue('file', 0):
245 if not os.path.exists(field):
246 md = md.checkmessage("file provided in {}: {} does not exist".format(fieldname, field))
247
248 #Check row of strings
249 if options.exist('stringrow'):
250 if not isinstance(field, list):
251 md = md.checkmessage(options.getfieldvalue('message', "field {} should be a list".format(fieldname)))
252
253 #Check forcings (size and times)
254 if options.getfieldvalue('timeseries', 0):
255 if field.shape[0] == md.mesh.numberofvertices or field.shape[0] == md.mesh.numberofelements:
256 if np.ndim(field) > 1 and not np.size(field, 1) == 1:
257 md = md.checkmessage(options.getfieldvalue('message', "field {} should have only one column as there are md.mesh.numberofvertices lines".format(fieldname)))
258 elif field.shape[0] == md.mesh.numberofvertices + 1 or field.shape[0] == md.mesh.numberofelements + 1:
259 if np.ndim(field) > 1 and not all(field[-1, :] == np.sort(field[-1, :])):
260 md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
261 if np.ndim(field) > 1 and any(field[-1, 0:-1] == field[-1, 1:]):
262 md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
263 else:
264 md = md.checkmessage(options.getfieldvalue('message', "field {} should have md.mesh.numberofvertices or md.mesh.numberofvertices + 1 lines".format(fieldname)))
265
266 #Check single value forcings (size and times)
267 if options.getfieldvalue('singletimeseries', 0):
268 if field.shape[0] == 2:
269 if not all(field[-1, :] == np.sort(field[-1, :])):
270 md = md.checkmessage(options.getfieldvalue('message', "field {} columns should be sorted chronologically".format(fieldname)))
271 if any(field[-1, 0:-1] == field[-1, 1:]):
272 md = md.checkmessage(options.getfieldvalue('message', "field {} columns must not contain duplicate timesteps".format(fieldname)))
273 elif field.shape[0] == 1:
274 if np.ndim(field) > 1 and not np.size(field, 1) == 1:
275 md = md.checkmessage(options.getfieldvalue('message', "field {} should be either a scalar or have 2 lines".format(fieldname)))
276 else:
277 md = md.checkmessage(options.getfieldvalue('message', "field {} should have 2 lines or be a scalar".format(fieldname)))
278
279 return md
Note: See TracBrowser for help on using the repository browser.