Skip to content

Commit 7fcb6ea

Browse files
JohnVolkJohnVolk
committed
fix bugs in optimizer, improve resampling routine, add archive method, add delete util functions
1 parent 41aeb2b commit 7fcb6ea

2 files changed

Lines changed: 101 additions & 64 deletions

File tree

‎prms_python/optimizer.py‎

Lines changed: 69 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -46,8 +46,9 @@ class Optimizer:
4646
#dic for min/max of parameter allowable ranges, add more when needed
4747
param_ranges = {'dday_intcp': (-60.0, 10.0), 'dday_slope': (0.2, 0.9),\
4848
'jh_coef': (0.005, 0.06), 'pt_alpha': (1.0, 2.0), \
49-
'potet_coef_hru_mo': (1.0, 1.6), 'tmax_index': \
50-
(-10.0, 110.0)\
49+
'potet_coef_hru_mo': (1.0, 2.0), 'tmax_index': \
50+
(-10.0, 110.0), 'tmin_lapse': (-10.0, 10.0), \
51+
'soil_moist_max': (0.001, 10.0), 'rain_adj': (0.5, 2.0)\
5152
} #changed potet max
5253

5354
def __init__(self, parameters, data, control_file, working_dir,
@@ -85,7 +86,7 @@ def __init__(self, parameters, data, control_file, working_dir,
8586
self.arb_outputs = []
8687

8788
def monte_carlo(self, reference_path, param_names, statvar_name, \
88-
stage, n_sims=10, method='uniform',\
89+
stage, n_sims=10, method='uniform', mu_factor=1,\
8990
noise_factor=0.1, nproc=None):
9091
'''
9192
Optimize the monthly dday_intcp and dday_slope parameters
@@ -102,8 +103,12 @@ def monte_carlo(self, reference_path, param_names, statvar_name, \
102103
n_sims (int): number of simulations to conduct
103104
parameter optimization/uncertaitnty analysis.
104105
method (str): resampling method for parameters (normal or uniform)
106+
mu_factor (float): coefficient to scale mean of the parameter(s)
107+
to resample from when using the normal distribution to resample
108+
i.e. a value of 1.5 will sample from a normal rv with mean
109+
50% higher than the original parameter mean
105110
noise_factor (float): scales the variance of noise to add to
106-
parameter values when adding normal rv (method='normal')
111+
parameter values when using normal rv (method='normal')
107112
nproc (int): number of processors available to run PRMS simulations
108113
'''
109114
if '_' in stage:
@@ -122,7 +127,7 @@ def monte_carlo(self, reference_path, param_names, statvar_name, \
122127
tmp = []
123128
for idx in range(n_sims):
124129
tmp.append(resample_param(self.parameters, name, how=method,\
125-
noise_factor=noise_factor))
130+
mu_factor=mu_factor, noise_factor=noise_factor))
126131
params.append(list(tmp))
127132

128133
# SimulationSeries comprised of each resampled param set
@@ -157,6 +162,7 @@ def monte_carlo(self, reference_path, param_names, statvar_name, \
157162
'end_time' : str(end_time),
158163
'measured' : reference_path,
159164
'method' : 'Monte Carlo',
165+
'mu_factor' : mu_factor,
160166
'noise_factor' : noise_factor,
161167
'resample': method,
162168
'sim_dirs' : [],
@@ -359,16 +365,21 @@ def _create_metafile_name(out_dir, opt_title, stage):
359365
name = '{}_{}_opt{}.json'.format(opt_title, stage, n)
360366
return name
361367

362-
def resample_param(params, param_name, how='uniform', noise_factor=0.1):
368+
def resample_param(params, param_name, how='uniform', mu_factor=1,\
369+
noise_factor=0.1):
363370
"""
364371
Resample PRMS parameter by shifting all values by a constant that is
365372
taken from a uniform distribution, where the range of the uniform
366373
values is equal to the difference between the min(max) of the parameter
367-
set and the min(max) of the allowable range from PRMS. For parameters
368-
that have array length <= 366 add noise to each parameter element by
369-
adding a RV from a normal distribution with mean 0, sigma = param
370-
allowable range / 10.
371-
374+
set and the min(max) of the allowable range from PRMS. If the resampling
375+
method ("how" argument) is set to 'normal', randomly sample a normal
376+
distribution with mean 0 and sigma = param allowable range multiplied by
377+
noise_factor, then add this random value to original param value. If
378+
parameters have array length <= 366 then individual parameter values are
379+
resampled otherwise resample all param values at once, i.e. by taking
380+
a single random value from the uniform distribution or by adding a single
381+
random value samplied from the normal distribution to all param values.
382+
372383
Args:
373384
params (parameters.Parameters): parameter object
374385
param_name (str): name of PRMS parameter to resample
@@ -379,8 +390,7 @@ def resample_param(params, param_name, how='uniform', noise_factor=0.1):
379390
noise_factor (float): factor to multiply parameter range by,
380391
use the result as the standard deviation for the normal rand.
381392
variable used to add element wise noise. i.e. higher
382-
noise facter will result in higher noise added to each param
383-
element. Must be > 0.
393+
noise facter will result in higher variance. Must be > 0.
384394
Returns:
385395
ret (numpy.ndarry): ndarray of param after uniform random mean
386396
shift or element-wise noise addition (normal r.v.)
@@ -395,15 +405,15 @@ def resample_param(params, param_name, how='uniform', noise_factor=0.1):
395405
if p_min == p_max == -1:
396406
raise ValueError("""{} has not been added to the dictionary of
397407
parameters to resample, add it's allowable min and max value
398-
to the param_ranges dictionary in the resample function in
408+
to the Optimizer.param_ranges attribute in
399409
Optimizer.py""".format(param_name))
400410

401411
dim_case = None
402412
nhru = params.dimensions['nhru']
403413
ndims = param_dic.get(param_name)['ndims']
404414
dimnames = param_dic.get(param_name)['dimnames']
405415
length = param_dic.get(param_name)['length']
406-
param = params[param_name]
416+
param = deepcopy(params[param_name])
407417

408418
# could expand list and check parameter name also e.g. cascade_flg
409419
# is a parameter that should not be changed
@@ -436,52 +446,48 @@ def resample_param(params, param_name, how='uniform', noise_factor=0.1):
436446

437447
low_bnd = p_min - np.min(param) # lowest param value minus allowable min
438448
up_bnd = p_max - np.max(param)
439-
s = (p_max - p_min) * noise_factor # variance noise, default: range*(1/10)
449+
s = (p_max - p_min) * noise_factor # std_dev (s) default: param_range/10
440450
#do resampling differently based on param dimensions
441451
if dim_case == 'resample_all_values_once':
442452
if how == 'uniform':
443-
#uniform RV for shifting all values once
444453
shifted_param = np.random.uniform(low=low_bnd, high=up_bnd) + param
445454
ret = shifted_param
446-
elif how == 'normal':
447-
while True:
448-
tmp = np.random.normal(0, s, size=param.shape) + param
449-
if np.max(tmp) <= p_max and np.min(tmp) >= p_min:
450-
ret = tmp
451-
break
455+
elif how == 'normal': # scale parameter mean if mu_factor given
456+
mu = np.mean(param) * mu_factor
457+
if mu_factor != 1:
458+
tmp = np.random.normal(mu, s, size=param.shape)
459+
ret = tmp + param
460+
else: # if default mu_factor, add noise from N(0,s)
461+
ret = np.random.normal(0, s, size=param.shape) + param
462+
452463
elif dim_case == 'resample_each_value':
453-
ret = copy(param)
464+
ret = param
454465
if how == 'uniform':
455466
for i, el in enumerate(param):
456467
low_bnd = p_min - np.min(el) # a,b for uniform RV to add
457468
up_bnd = p_max - np.max(el)
458469
ret[i] = el + np.random.uniform(low=low_bnd, high=up_bnd)
459-
elif how == 'normal':
470+
elif how == 'normal': # the original value is the mean
460471
for i, el in enumerate(param):
461-
while True:
462-
low_bnd = p_min - np.min(el)
463-
up_bnd = p_max - np.max(el)
464-
tmp = el + np.random.normal(0, s)
465-
if np.max(tmp) <= p_max and np.min(tmp) >= p_min:
466-
ret[i] = tmp
467-
break
472+
mu = el * mu_factor
473+
ret[i] = np.random.normal(mu, s)
474+
468475
# nhru by nmonth dimensional params
469476
elif dim_case == 'nhru_nmonths':
470-
ret = copy(param)
477+
ret = param
471478
if how == 'uniform':
472479
rvs = [np.random.uniform(low=low_bnd, high=up_bnd)\
473480
for i in range(12)]
474481
for month in range(12):
475482
ret[month] += rvs[month]
476483
elif how == 'normal':
477484
for i, el in enumerate(param):
478-
while True:
479-
low_bnd = p_min - np.min(el)
480-
up_bnd = p_max - np.max(el)
481-
tmp = el + np.random.normal(0, s)
482-
if np.max(tmp) <= p_max and np.min(tmp) >= p_min:
483-
ret[i] = tmp
484-
break
485+
mu = np.mean(el) * mu_factor
486+
if mu_factor != 1:
487+
tmp = np.random.normal(mu, s, size=el.shape)
488+
ret[i] = tmp + el
489+
else:
490+
ret[i] = np.random.normal(0, s, size=el.shape) + el
485491

486492
return ret
487493

@@ -529,7 +535,7 @@ def _get_optr_jsons(self, work_dir, stage):
529535
Create dictionary of each optimization with stage as key and lists
530536
of corresponding json file paths as values.
531537
532-
Arguments:
538+
Args:
533539
work_dir (str): path to directory with model results, i.e.
534540
location where simulation series outputs and optimization
535541
json files are located, aka Optimizer.working_dir
@@ -605,7 +611,8 @@ def result_table(self, freq='daily', top_n=5, latex=False):
605611
##TODO: add stats for freq options annual (means or sum)
606612

607613
sim_dirs = self._get_sim_dirs(self.stage)
608-
if top_n >= len(sim_dirs): top_n = len(sim_dirs)
614+
if top_n >= len(sim_dirs):
615+
top_n = len(sim_dirs) + 1 # for returning inclusive last sim
609616
sim_names = [path.split(os.sep)[-1] for path in sim_dirs]
610617
meas_var = self._get_measured(self.stage)
611618
statvar_name = self._get_statvar_name(self.stage)
@@ -633,7 +640,7 @@ def result_table(self, freq='daily', top_n=5, latex=False):
633640
sim_out = load_statvar(OPJ(sim, 'outputs', 'statvar.dat'))\
634641
['{}'.format(statvar_name)]
635642
except: # simulation might have been removed or missing
636-
continue
643+
pass
637644

638645
sim_out = sim_out[idx]
639646
if freq == 'daily':
@@ -674,7 +681,7 @@ def result_table(self, freq='daily', top_n=5, latex=False):
674681

675682
def get_top_ranked_sims(self, sorted_df):
676683
# use result table to make dic with best param and statvar paths
677-
# index of table is the simulation directory names
684+
# index of table are simulation directory names
678685
ret = {
679686
'dir_name' : [],
680687
'param_path' : [],
@@ -706,15 +713,15 @@ def archive(self, remove_sims=True, metric_freq='daily'):
706713
metrics for each Optimizer simulation of the
707714
OptimizationResult.stage in the OptimizationResult.working_dir.
708715
709-
Arguments:
716+
Kwargs:
710717
remove_sims (bool) : If True recursively delete all folders
711718
and files associated with original simulations of the
712719
OptimizationResult.stage in the
713720
OptimizationResult.working_dir, if False do not delete
714721
simulations.
715722
metric_freq (Str) : Frequency of output metric computation
716723
for recording of model performance. Can be 'daily'
717-
(default) or 'monthly'. Note, new results can be computed
724+
(default) or 'monthly'. Note, other results can be computed
718725
later with archived results.
719726
Returns:
720727
None
@@ -724,13 +731,14 @@ def archive(self, remove_sims=True, metric_freq='daily'):
724731
if not os.path.isdir(archive_dir):
725732
os.mkdir(archive_dir)
726733

727-
# create table and used to make mapping dic
734+
# create table and use to make mapping dic
728735
table = self.result_table(freq=metric_freq,\
729736
top_n=self.total_sims, latex=False)
730737
map_dic = self.get_top_ranked_sims(table)
731738

739+
metadata_json_paths = self.metadata_json_paths[self.stage]
732740
# get measured optimization variable path
733-
first_json = self.metadata_json_paths[self.stage][0]
741+
first_json = metadata_json_paths[0]
734742

735743
with open(first_json) as json_file:
736744
json_data = json.load(json_file)
@@ -747,6 +755,16 @@ def archive(self, remove_sims=True, metric_freq='daily'):
747755
[self.statvar_name]
748756
except: # simulation directory was already removed
749757
continue
758+
759+
# look for resampling method info for the particular simulation
760+
for f in metadata_json_paths:
761+
with open(f) as tmp_file:
762+
tmp = json.load(tmp_file)
763+
if OPJ(self.working_dir, sim) in tmp.get('sim_dirs'):
764+
resample = tmp.get('resample')
765+
noise_factor = tmp.get('noise_factor')
766+
mu_factor = tmp.get('mu_factor')
767+
750768
json_data = {
751769
'param_names' : [],
752770
'param_values' : [],
@@ -757,6 +775,9 @@ def archive(self, remove_sims=True, metric_freq='daily'):
757775
index.astype(str).tolist(),
758776
'output_values' : output_series.values.tolist(),
759777
'metric_freq' : metric_freq,
778+
'resample' : resample,
779+
'mu_factor' : mu_factor,
780+
'noise_factor' : noise_factor,
760781
'NSE' : table.loc[sim, 'NSE'],
761782
'RMSE' : table.loc[sim, 'RMSE'],
762783
'PBIAS' : table.loc[sim, 'PBIAS'],

‎prms_python/util.py‎

Lines changed: 32 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -24,25 +24,41 @@ def remove_all_optimization_sims_of_other_stage(work_directory, stage):
2424
None
2525
"""
2626
from .optimizer import OptimizationResult # avoid circular import
27-
result = OptimizationResult(work_directory,stage=stage)
28-
tracked_dirs = []
29-
for f in result.metadata_json_paths[stage]:
30-
with open(f) as fh:
31-
json_data = json.load(fh)
32-
tracked_dirs.extend(json_data.get('sim_dirs'))
33-
count = 0
34-
for d in os.listdir(result.working_dir):
35-
path = os.path.join(result.working_dir, d)
36-
if path in tracked_dirs:
37-
continue
38-
elif os.path.isdir(path):
39-
count+=1
40-
for dirpath, dirnames, filenames in os.walk(path, topdown=False):
41-
shutil.rmtree(dirpath, ignore_errors=True)
27+
28+
try:
29+
result = OptimizationResult(work_directory,stage=stage)
30+
tracked_dirs = []
31+
for f in result.metadata_json_paths[stage]:
32+
with open(f) as fh:
33+
json_data = json.load(fh)
34+
tracked_dirs.extend(json_data.get('sim_dirs'))
35+
count = 0
36+
for d in os.listdir(result.working_dir):
37+
path = os.path.join(result.working_dir, d)
38+
if path in tracked_dirs:
39+
continue
40+
elif os.path.isdir(path) and '_archived' not in path:
41+
count+=1
42+
for dirpath, dirnames, filenames in os.walk(path,\
43+
topdown=False):
44+
shutil.rmtree(dirpath, ignore_errors=True)
45+
46+
# if no json file in working dir for given stage, delete any other sim dirs
47+
except:
48+
count = 0
49+
for d in os.listdir(work_directory):
50+
path = os.path.join(work_directory, d)
51+
if os.path.isdir(path) and '_archived' not in path:
52+
count+=1
53+
for dirpath, dirnames, filenames in os.walk(path,\
54+
topdown=False):
55+
shutil.rmtree(dirpath, ignore_errors=True)
56+
4257
print('deleted {} simulations that were either not tracked by a JSON file'\
4358
.format(count) + ' or were not part of {} optimization stage'\
4459
.format(stage))
45-
60+
61+
4662
def delete_files(work_directory, file_name=''):
4763
"""
4864
Recursively delete all files of a certain name from multiple PRMS

0 commit comments

Comments
 (0)