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

Last change on this file since 25836 was 25836, checked in by Mathieu Morlighem, 4 years ago

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