Skip to content

Commit 4837ac4

Browse files
JohnVolkJohnVolk
committed
* add Optimizer.pet method, tested and working
* adjust function _mod_params and srad method (function) to reuse functionality in Optimizer.pet
1 parent cf7b1af commit 4837ac4

1 file changed

Lines changed: 102 additions & 19 deletions

File tree

‎prms_python/optimizer.py‎

Lines changed: 102 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -40,9 +40,10 @@ class Optimizer:
4040
4141
'''
4242

43-
## constant attributes for allowable range of solrad parameters
43+
## constant attributes for allowable range of select PRMS parameters
4444
ir = (-60.0, 10.0) # PRMS dday_intcp range
4545
sr = (0.2, 0.9) # dday_slope range
46+
jh = (0.005, 0.06) # jh_coef range
4647

4748
def __init__(self, parameters, data, control_file, working_dir,
4849
title, description=None):
@@ -76,13 +77,17 @@ def __init__(self, parameters, data, control_file, working_dir,
7677
self.description = description
7778
self.srad_outputs = []
7879
self.measured_srad = None
79-
self.srad_hru = None
80+
self.srad_hru = None
81+
self.pet_outputs = []
82+
self.measured_pet = None
83+
self.pet_hru = None
8084

8185
def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
8286
nproc=None):
8387
'''
84-
Optimize the monthly dday_intcp and dday_slope parameters by one of
85-
two methods: 'uniform' or 'random' for uniform sampling
88+
Optimize the monthly dday_intcp and dday_slope parameters
89+
(two key parameters in the ddsolrad module in PRMS) by one of
90+
multiple methods (in development): Monte Carlo default method
8691
8792
Args:
8893
reference_srad_path (str): path to measured solar radiation data
@@ -91,13 +96,9 @@ def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
9196
have swrad 'nhru' listed as a statvar output in your
9297
control file.
9398
Kwargs:
94-
method (str): XXX not yet implemented- all uniform now
95-
'uniform' or 'random'; if 'random',
96-
intcp_delta and slope_delta are ignored, if provided
99+
method (str): XXX not yet implemented-
97100
n_sims (int): number of simulations to conduct
98101
parameter optimization/uncertaitnty analysis.
99-
Returns:
100-
(SradOptimizationResult)
101102
'''
102103
# assign the optimization object a copy of measured srad for plots
103104
self.measured_srad = pd.Series.from_csv(
@@ -121,7 +122,8 @@ def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
121122

122123
series = SimulationSeries(
123124
Simulation.from_data(
124-
self.data, _mod_params(self.parameters, intcps[i], slopes[i]),
125+
self.data, _mod_params(self.parameters, [intcps[i], slopes[i]],\
126+
'swrad'),
125127
self.control_file,
126128
OPJ(
127129
self.working_dir,
@@ -168,6 +170,86 @@ def _error(x, y):
168170

169171
print('{0}\nOutput information sent to {1}\n'.format('-' * 80, json_outfile))
170172

173+
def pet(self, reference_pet_path, station_nhru, n_sims=10, method='',\
174+
nproc=None):
175+
'''
176+
Optimize the monthly Jenson Haise coefficients (jh_coef) parameters
177+
(a key parameter in the potet_jh et module in PRMS) by one of
178+
multiple methods (in development): Monte Carlo default method
179+
Future development- implement for other PRMS et modules (e.g. Penmen Monteith)
180+
181+
Args:
182+
reference_pet_path (str): path to measured pet data
183+
station_nhru (int or str): hru index in PRMS that is geographically
184+
near the measured/estimated pet location. TODO: if value
185+
if 'basin' then you wish to compare to area-weighted
186+
simulated pet (`basin_potet_1` in PRMS).
187+
Kwargs:
188+
method (str): XXX not yet implemented- (default-Monte Carlo)
189+
n_sims (int): number of simulations to conduct
190+
parameter optimization/uncertaitnty analysis.
191+
'''
192+
# assign the optimization object a copy of measured pet
193+
self.measured_pet = pd.Series.from_csv(
194+
reference_pet_path, parse_dates=True
195+
)
196+
197+
self.pet_hru = station_nhru
198+
199+
pet_start_time = dt.datetime.now()
200+
pet_start_time = pet_start_time.replace(second=0, microsecond=0)
201+
202+
## shifting all parameter values by small random amount from normal distribution
203+
## generate parameter set for each simulation
204+
jh_coefs = [_resample_param(self.parameters['jh_coef'], Optimizer.jh[0],\
205+
Optimizer.jh[1]) for i in range(n_sims)]
206+
207+
self.data.write(OPJ(self.working_dir, 'data'))
208+
209+
series = SimulationSeries(
210+
Simulation.from_data(
211+
self.data, _mod_params(self.parameters, [jh_coefs[i]], 'pet'),
212+
self.control_file,
213+
OPJ(
214+
self.working_dir,
215+
'jh_coef:{0:.3f}'.format(np.mean(jh_coefs[i]))
216+
)
217+
)
218+
for i in range(n_sims)
219+
)
220+
221+
# run all scenarios
222+
outputs = list(series.run(nproc=nproc).outputs_iter())
223+
self.pet_outputs.extend(outputs)
224+
225+
pet_end_time = dt.datetime.now()
226+
pet_end_time = pet_end_time.replace(second=0, microsecond=0)
227+
228+
pet_meta = {'stage' : 'pet',
229+
'pet_hru_id' : self.pet_hru,
230+
'optimization_title' : self.title,
231+
'optimization_description' : self.description,
232+
'start_time' : str(pet_start_time),
233+
'end_time' : str(pet_end_time),
234+
'measured_pet' : reference_pet_path,
235+
'sim_dirs' : [],
236+
'original_params' : self.parameters.base_file,
237+
'n_sims' : n_sims
238+
}
239+
240+
for output in outputs:
241+
pet_meta['sim_dirs'].append(output['simulation_dir'])
242+
243+
json_outfile = OPJ(self.working_dir, _create_metafile_name(\
244+
self.working_dir, self.title, 'pet'))
245+
246+
with open(json_outfile, 'w') as outf:
247+
json.dump(pet_meta, outf, sort_keys = True, indent = 4,\
248+
ensure_ascii = False)
249+
250+
print('{0}\nOutput information sent to {1}\n'.format('-' * 80, json_outfile))
251+
252+
#TODO: make single plot function that takes stage as input (swrad, pet, flow...)
171253
def plot_srad_optimization(self, freq='daily', method='time_series',\
172254
plot_vars='both', return_fig=False):
173255
"""
@@ -308,7 +390,7 @@ def _create_metafile_name(out_dir, opt_title, stage):
308390
look for all metadata simulation json files and find out if the
309391
current simulation is a replicate. Then use that information to
310392
build the correct file name for the output json file. The series
311-
are typically run in parallelel that is why this step has to be
393+
are typically run in parallel that is why this step has to be
312394
done after running multiple simulations from an optimization stage.
313395
314396
Args:
@@ -325,11 +407,11 @@ def _create_metafile_name(out_dir, opt_title, stage):
325407
'dry_creek' the next json file will be returned as
326408
'dry_creek_swrad_opt1.json' and so on with integer increments
327409
"""
328-
swrad_meta_re = re.compile(r'^{}_{}_opt(\d*)\.json'.format(opt_title, stage))
410+
meta_re = re.compile(r'^{}_{}_opt(\d*)\.json'.format(opt_title, stage))
329411
reps = []
330412
for f in os.listdir(out_dir):
331-
if swrad_meta_re.match(f):
332-
nrep = swrad_meta_re.match(f).group(1)
413+
if meta_re.match(f):
414+
nrep = meta_re.match(f).group(1)
333415
if nrep == '':
334416
reps.append(0)
335417
else:
@@ -378,15 +460,16 @@ def _resample_param(param, p_min, p_max, noise_factor=0.1):
378460
if np.max(tmp) <= p_max and np.min(tmp) >= p_min:
379461
return tmp
380462

381-
def _mod_params(parameters, intcp, slope):
463+
def _mod_params(parameters, params, stage):
382464
# deepcopy was crashing, raising:
383465
# TypeError: cannot serialize '_io.TextIOWrapper' object
384466
ret = copy(parameters)
385467
#print (intcp, slope)
386-
387-
ret['dday_intcp'] = intcp
388-
ret['dday_slope'] = slope
389-
468+
if stage == 'swrad':
469+
ret['dday_intcp'] = params[0]
470+
ret['dday_slope'] = params[1]
471+
elif stage == 'pet':
472+
ret['jh_coef'] = params[0]
390473
return ret
391474

392475

0 commit comments

Comments
 (0)