Coverage for pygeodyn/inout/config.py: 77%
245 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-22 13:43 +0000
1import numpy as np
2import logging
3from pygeodyn.inout import reads
4import pygeodyn
5import os
6from os.path import isfile, join
8class ComputationConfig(object):
9 """
10 Defines the configuration of the computation, namely:
11 - All core state degrees and number of coefs
12 - Physical constants
13 - Forecast and Analysis times
14 - Paths of the prior and observation data
15 - Output file features
16 """
18 # Initialization methods
19 def __init__(self, do_shear, config_file=None):
20 """
21 Loads the config according to parameters. First tries to load config file, then config dict if the config_file is None.
23 :param do_shear: control if we add shear config parameters to config dict
24 :type do_shear: int
25 :param config_file: filename containing the config (see format hereafter)
26 :type config_file: str (default None)
28 config file format:
29 name type value
31 Example :
32 Lb int 10
34 supported format: int/float/str/bool
35 lines starting by # are not considered.
36 """
37 self.Lb = None
38 self.Lsv = None
39 self.Lu = None
40 self.Ly = None
42 self.do_shear = do_shear
44 # data is a string, assuming the filename of a config file
45 if config_file is not None:
46 self.read_from_file(config_file)
48 if self.Lb is None:
49 raise ValueError("No maximal degree was read for the magnetic field !")
50 if self.Lsv is None:
51 raise ValueError("No maximal degree was read for the secular variation !")
52 if self.Lu is None:
53 raise ValueError("No maximal degree was read for the flow !")
55 # Once the data was loaded, build all the relevant variables
56 self.init_bools()
57 self.init_out()
58 self.init_physical_constants()
59 self.init_paths()
60 self.init_times()
61 self.init_corestate_inits()
62 self.init_observations()
63 if self.do_shear:
64 if self.Ly is None:
65 raise ValueError("No maximal degree was read for the constraint on the shear !")
67 self.init_shear()
69 @property
70 def Nb(self):
71 return self.Lb * (self.Lb + 2)
73 @property
74 def Nsv(self):
75 return self.Lsv * (self.Lsv + 2)
77 @property
78 def Nu(self):
79 return self.Lu * (self.Lu + 2)
81 @property
82 def Ny(self):
83 return self.Ly * (self.Ly + 2)
85 @property
86 def Nu2(self):
87 return 2*self.Nu
89 @property
90 def Ny2(self):
91 return 2*self.Ny
93 @property
94 def Nz(self):
95 if self.pca:
96 return self.N_pca_u + self.Nsv
97 else:
98 return self.Nu2 + self.Nsv
100 @property
101 def Nuz(self):
102 if self.pca:
103 return self.N_pca_u
104 else:
105 return self.Nu2
107 def init_bools(self):
108 """
109 Init computation booleans.
110 """
111 if not hasattr(self, 'compute_e'):
112 logging.warning("Compute ER was not set in the config file! Assuming True.")
113 self.compute_e = 1
115 if not hasattr(self, 'kalman_norm'):
116 logging.warning("Kalman_norm was not set in the config file! Assuming l2 norm.")
117 self.kalman_norm = 'l2'
118 else:
119 self.kalman_norm = self.kalman_norm.lower() # tolerance with caps
121 if not hasattr(self, 'N_pca_u'):
122 logging.warning("Number of coefficients for PCA of U was not set in the config file! Assuming no PCA.")
123 self.N_pca_u = 0
124 self.pca = 0
125 else:
126 if self.N_pca_u == 0:
127 self.pca = 0
128 else:
129 self.pca = 1
131 if not hasattr(self, 'remove_spurious'):
132 logging.warning("Value to discard spurious correlations was not set! Assuming diagonal.")
133 self.remove_spurious = np.inf
135 if not hasattr(self, 'AR_type'):
136 logging.warning("Type of AR process was not set in the config file! Assuming diagonal.")
137 self.AR_type = 'diag'
139 if self.N_pca_u == 0 and self.AR_type == "AR3":
140 raise ValueError("AR3 process requires N_pca_u to be non-zero")
142 if not hasattr(self, 'combined_U_ER_forecast'):
143 logging.warning("U ER forecast dependancy was not set in the config file! Assuming independant.")
144 self.combined_U_ER_forecast = 0
146 if self.combined_U_ER_forecast == 1 and self.AR_type == 'diag':
147 raise ValueError("U ER must be independant for diag forecast")
150 def init_out(self):
151 if not hasattr(self, 'out_computed'):
152 logging.warning("Out_computed was not set in the config file ! Assuming computed states not saved.")
153 self.out_computed = 0
155 if not hasattr(self, 'out_analysis'):
156 logging.warning("Out_analysis was not set in the config file ! Assuming analysis states saved.")
157 self.out_analysis = 1
159 if not hasattr(self, 'out_forecast'):
160 logging.warning("Out_forecast was not set in the config file ! Assuming forecast states not saved.")
161 self.out_forecast = 0
163 if not hasattr(self, 'out_misfits'):
164 logging.warning("Out_misfits was not set in the config file ! Assuming misfits not saved.")
165 self.out_misfits = 0
167 if not hasattr(self, 'out_format'):
168 logging.warning("Out_format was not set in the config file ! Assuming single precision.")
169 self.out_format = 'float32'
173 def init_corestate_inits(self):
174 if not hasattr(self, 'core_state_init'):
175 logging.warning("Type of CoreState initialisation was not set in the config file! Assuming normal draw around average priors.")
176 self.core_state_init = 'normal'
178 if self.core_state_init == 'from_file':
179 if not (hasattr(self, 'init_file') and hasattr(self, 'init_date')):
180 logging.warning(
181 "CoreState initialisation was set to 'from_file' but no file or date were given. Falling back to normal draw initialisation.")
182 self.core_state_init = 'normal'
184 def init_paths(self):
185 """
186 Builds the full path according to the path given in conf or set default.
187 """
188 if not hasattr(self, 'prior_dir'):
189 logging.warning("Prior directory was not set, assuming 71path.")
190 self.prior_dir = 'data/priors/71path'
191 self.prior_type = '71path'
193 if self.prior_type == '0path':
194 if not self.AR_type == "diag":
195 raise ValueError('AR_type must be diag for 0path prior')
198 if not hasattr(self, 'prior_type'):
199 raise ValueError('Prior directory was set ({}) without prior type ! Please set the prior type.'.format(self.prior_dir))
201 # Default behaviour is to read the COVOBS-x2 model as observations
202 if not hasattr(self, 'obs_dir'):
203 self.obs_dir = 'data/observations/COVOBS-x2_maglat'
204 self.obs_type = 'COVOBS_hdf5'
206 paths = ['prior_dir', 'obs_dir']
207 if hasattr(self, 'init_file'):
208 paths.append('init_file')
210 for path_key in paths:
211 # If not absolute, build the absolute path from the given relative path
212 given_path = self.__dict__[path_key]
213 if not os.path.isabs(given_path):
214 given_path = os.path.normpath(os.path.join(pygeodyn._package_directory, given_path))
216 if not os.path.exists(given_path):
217 # Try to correct the path
218 corrected_path = os.path.join(pygeodyn._package_directory, given_path.split('pygeodyn/')[-1])
219 # If the corrected path does not exist either, throw an error
220 if not os.path.exists(corrected_path):
221 raise IOError('{} given in config does not exist ! Did you supply a relative path instead of an absolute one ?'.format(given_path))
222 given_path = corrected_path
223 self.__dict__[path_key] = given_path
226 def init_observations(self):
227 if self.obs_type == 'GO_VO':
228 if not hasattr(self, 'obs_select'):
229 logging.warning('No values of govo_obs_list found, default take all obs files in observation directory')
230 self.obs_select = all
231 if not hasattr(self, 'discard_high_lat'):
232 logging.warning("Suppression of high latitude data was not set, assuming no suppression")
233 self.discard_high_lat = False
234 if not hasattr(self, 'SW_err'): # FOR 4MONTHS SWARM DATA ONLY!!!
235 logging.warning('SW_err was not set, assuming diag error matrix for SWARM')
236 self.SW_err = "diag"
237 else:
238 self.obs_select = None
242 def init_physical_constants(self):
243 """
244 Builds the physical constants of the computation.
245 """
247 # Check for theta steps for Legendre polynomial
248 if "Nth_legendre" not in self.__dict__:
249 logging.warning("No number of theta steps was read ! Using default value 64.")
250 self.Nth_legendre = 64
252 # Check the presence of times
253 if "TauU" not in self.__dict__:
254 logging.warning("No characteristic time was read for the core flow ! Using default value 30 yrs.")
255 self.TauU = np.timedelta64(30, 'Y')
256 if "TauE" not in self.__dict__:
257 logging.warning("No characteristic time was read for the subgrid errors ! Using default value 10 yrs.")
258 self.TauE = np.timedelta64(10, 'Y')
260 # Convert TauU and TauE in decimal values (in years)
261 year_timedelta64 = np.timedelta64(1, 'Y')
262 self.TauU = self.TauU / year_timedelta64
263 self.TauE = self.TauE / year_timedelta64
265 def init_times(self):
266 """ Builds the variables linked to the times of computation (forecasts and analysis). """
268 if "t_start_analysis" not in self.__dict__:
269 raise ValueError("No analysis starting time was read !")
271 if "t_end_analysis" not in self.__dict__:
272 raise ValueError("No analysis end time was read !")
274 if "dt_f" not in self.__dict__:
275 raise ValueError("No forecast timestep was read !")
277 if self.dt_f == 0:
278 raise ValueError("Forecast timestep cannot be zero !")
280 if "dt_a_f_ratio" not in self.__dict__:
281 raise ValueError("No dt_analysis/dt_forecast ratio was read !")
283 if int(self.dt_a_f_ratio) != self.dt_a_f_ratio:
284 raise ValueError("dt_a_f_ratio must be an integer!")
285 else:
286 self.dt_a = self.dt_a_f_ratio * self.dt_f
288 if "N_dta_start_forecast" not in self.__dict__:
289 self.N_dta_start_forecast = 1
290 else:
291 if self.N_dta_start_forecast < 1:
292 raise ValueError("N_dta_start_forecast must be greater or equal to 1")
294 if self.AR_type == "AR3":
295 if "N_dta_end_forecast" not in self.__dict__:
296 self.N_dta_end_forecast = 1
297 else:
298 if self.N_dta_end_forecast < 1:
299 raise ValueError("N_dta_end_forecast must be greater or equal to 1")
300 else:
301 if "N_dta_end_forecast" not in self.__dict__:
302 self.N_dta_end_forecast = 0
304 if "dt_smoothing" not in self.__dict__:
305 logging.warning("No prior dt smoothing was read ! Using default value 3.2 years.")
306 self.dt_smoothing = 3.2
308 if "dt_sampling" not in self.__dict__:
309 logging.warning("No prior dt_sampling was read ! Using default value 5 years.")
310 self.dt_sampling = 5
312 t_start_forecast = self.t_start_analysis - self.dt_a * self.N_dta_start_forecast
313 t_end_forecast = self.t_end_analysis + self.dt_a * self.N_dta_end_forecast + self.dt_f
315 # Build time arrays
316 # First forecast at t_start + dt_f
317 self.t_forecasts = np.arange(t_start_forecast, t_end_forecast, self.dt_f)
318 # First analysis at t_start
319 self.t_analyses = np.arange(self.t_start_analysis, self.t_end_analysis + self.dt_f, self.dt_a)
320 if self.AR_type == "AR3":
321 self.t_analyses_full = np.arange(self.t_start_analysis - self.dt_a, self.t_end_analysis + self.dt_a + self.dt_f, self.dt_a)
322 else:
323 self.t_analyses_full = np.arange(self.t_start_analysis - self.dt_a, self.t_end_analysis + self.dt_f, self.dt_a)
325 def init_shear(self):
327 if not hasattr(self, 'TauG'):
328 logging.warning("TauG was not set in the config file! Using default value 20 yrs.")
329 self.TauG = np.timedelta64(20, 'Y')
331 year_timedelta64 = np.timedelta64(1, 'Y')
332 self.TauG = self.TauG / year_timedelta64
334 if not hasattr(self, 'remove_spurious_shear_u'):
335 logging.warning("Value to discard shear U spurious correlations was not set! Assuming diagonal.")
336 self.remove_spurious_shear_u = np.inf
338 if not hasattr(self, 'remove_spurious_shear_err'):
339 logging.warning("Value to discard shear error spurious correlations was not set! Assuming diagonal.")
340 self.remove_spurious_shear_err = np.inf
342 if not hasattr(self, 'prior_dir_shear'):
343 logging.warning("Shear prior directory was not set, assuming same prior as prior_dir.")
344 self.prior_dir_shear = self.prior_dir
345 self.prior_type_shear = self.prior_type
347 if not hasattr(self, 'prior_type_shear'):
348 raise ValueError('Shear prior directory was set ({}) without shear prior type ! Please set the shear prior type.'.format(self.prior_dir_shear))
350 paths = ['prior_dir_shear']
352 for path_key in paths:
353 # If not absolute, build the absolute path from the given relative path
354 given_path = self.__dict__[path_key]
355 if not os.path.isabs(given_path):
356 given_path = os.path.normpath(os.path.join(pygeodyn._package_directory, given_path))
358 if not os.path.exists(given_path):
359 # Try to correct the path
360 corrected_path = os.path.join(pygeodyn._package_directory, given_path.split('pygeodyn/')[-1])
361 # If the corrected path does not exist either, throw an error
362 if not os.path.exists(corrected_path):
363 raise IOError('{} given in config does not exist ! Did you supply a relative path instead of an absolute one ?'.format(given_path))
364 given_path = corrected_path
365 self.__dict__[path_key] = given_path
367 @property
368 def nb_forecasts(self):
369 return len(self.t_forecasts)
371 @property
372 def nb_analyses(self):
373 return len(self.t_analyses)
375 # Config methods
376 def read_from_file(self, filename):
377 """ Reads the config file from a filename and sets the read dictonary as object dict.
379 :param filename: the file from which we load the data (see format at end)
380 :type filename: str
381 """
382 file_conf = reads.parse_config_file(filename)
383 self.update(file_conf)
385 def update(self, config_dict):
386 """ Updates the internal dictionary of the object with the config_dict
388 :param config_dict: Configuration dictionary
389 :type config_dict: dict
390 """
391 self.__dict__.update(config_dict)
393 def save_hdf5(self, hdf5file):
394 """ save necessary attributes for webgeodyn reading
395 """
396 hdf5file.attrs["Lb"] = self.Lb
397 hdf5file.attrs["Lu"] = self.Lu
398 hdf5file.attrs["Lsv"] = self.Lsv