Skip to content

Commit f739486

Browse files
JohnVolkJohnVolk
committed
* basic functionality to OptimizerResult objects
 - table of top_n best results based on 4 obj functions - utilizes output json metadata files * improved srad, pet optimization methods to accept arbitrary number of parameters - new naming system based on first resampled param's mean value - improved resample and mod_params functions in optimizer.py * add new obj functions in util.py
1 parent ccb57a4 commit f739486

2 files changed

Lines changed: 199 additions & 50 deletions

File tree

‎prms_python/optimizer.py‎

Lines changed: 173 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
from .data import Data
1616
from .parameters import Parameters
1717
from .simulation import Simulation, SimulationSeries
18-
from .util import load_statvar
18+
from .util import load_statvar, nash_sutcliffe, percent_bias, rmse
1919

2020
OPJ = os.path.join
2121

@@ -43,8 +43,8 @@ class Optimizer:
4343
#dic for min/max of parameter allowable ranges, add more when needed
4444
param_ranges = {'dday_intcp': (-60.0, 10.0), 'dday_slope': (0.2, 0.9),\
4545
'jh_coef': (0.005, 0.06), 'pt_alpha': (1.0, 2.0), \
46-
'potet_coef_hru_mo': (1.0, 2.0)\
47-
}
46+
'potet_coef_hru_mo': (1.0, 1.6)\
47+
} #changed potet max
4848

4949
def __init__(self, parameters, data, control_file, working_dir,
5050
title, description=None):
@@ -84,7 +84,7 @@ def __init__(self, parameters, data, control_file, working_dir,
8484
self.pet_hru = None
8585

8686
def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
87-
srad_mod='ddsolrad', nproc=None):
87+
noise_factor=0.1, srad_mod='ddsolrad', nproc=None):
8888
'''
8989
Optimize the monthly dday_intcp and dday_slope parameters
9090
(two key parameters in the ddsolrad module in PRMS) by one of
@@ -118,22 +118,27 @@ def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
118118
## resampling degree day slope and intercept simoultaneously
119119
if srad_mod == 'ddsolrad':
120120
param_names = ['dday_intcp', 'dday_slope']
121-
intcps = [resample_param(self.parameters, 'dday_intcp') for i in range(n_sims)]
122-
slopes = [resample_param(self.parameters, 'dday_slope') for i in range(n_sims)]
123-
params = [intcps, slopes]
121+
122+
params = []
123+
for name in param_names: # list of lists of resampled params
124+
tmp = []
125+
for idx in range(n_sims):
126+
tmp.append(resample_param(self.parameters, name, how=method,\
127+
noise_factor=noise_factor))
128+
params.append(list(tmp))
124129

125130
self.data.write(OPJ(self.working_dir, 'data'))
126131

127-
# TODO: better file naming for output folders
132+
# sim dirs named after first resampled param name and mean value
128133
series = SimulationSeries(
129134
Simulation.from_data(
130-
self.data, _mod_params(self.parameters, [intcps[i], slopes[i]],\
131-
'swrad', srad_mod),
135+
self.data, _mod_params(self.parameters,\
136+
[params[n][i] for n in range(len(params))],\
137+
param_names),
132138
self.control_file,
133139
OPJ(
134140
self.working_dir,
135-
'intcp:{0:.3f}_slope:{1:.3f}'.format(np.mean(params[0][i]),\
136-
np.mean(params[1][i]))
141+
'{0}:{1:.6f}'.format(param_names[0], np.mean(params[0][i]))
137142
)
138143
)
139144
for i in range(n_sims)
@@ -177,7 +182,7 @@ def _error(x, y):
177182
print('{0}\nOutput information sent to {1}\n'.format('-' * 80, json_outfile))
178183

179184
def pet(self, reference_pet_path, station_nhru, n_sims=10,\
180-
pet_mod='potet_pt', method='', nproc=None):
185+
pet_mod='potet_pt', method='uniform', noise_factor=0.1, nproc=None):
181186
'''
182187
Optimize the monthly coefficients (depending on module) parameters
183188
by one of multiple methods (in development): Monte Carlo default method
@@ -203,32 +208,37 @@ def pet(self, reference_pet_path, station_nhru, n_sims=10,\
203208
else: # get variable at a specific hru
204209
self.pet_hru = 'potet_{}'.format(station_hru)
205210

211+
pet_start_time = dt.datetime.now()
212+
pet_start_time = pet_start_time.replace(second=0, microsecond=0)
213+
206214
if pet_mod == 'potet_pt':
207215
param_names = ['potet_coef_hru_mo']
208-
params = [resample_param(self.parameters, 'potet_coef_hru_mo')\
209-
for i in range(n_sims)]
210216
elif pet_mod == 'potet_jh':
211217
param_names = ['jh_coef']
212-
params = [resample_param(self.parameters, 'jh_coef') for i in range(n_sims)]
213-
214-
215-
pet_start_time = dt.datetime.now()
216-
pet_start_time = pet_start_time.replace(second=0, microsecond=0)
217-
218-
## generate parameter set for each simulation
219-
#self.data.write(OPJ(self.working_dir, 'data'))
220-
218+
# list of lists of resampled params
219+
params = []
220+
for name in param_names:
221+
tmp = []
222+
for idx in range(n_sims):
223+
tmp.append(resample_param(self.parameters, name, how=method,\
224+
noise_factor=noise_factor))
225+
params.append(list(tmp))
226+
227+
# sim dirs named after first resampled param name and mean value
221228
series = SimulationSeries(
222229
Simulation.from_data(
223-
self.data, _mod_params(self.parameters, [params[i]], 'pet', pet_mod),
230+
self.data, _mod_params(self.parameters,\
231+
[params[n][i] for n in range(len(params))],\
232+
param_names),
224233
self.control_file,
225234
OPJ(
226235
self.working_dir,
227-
'{0}:{1:.5f}'.format('_'.join(param_names), np.mean(params[i]))
236+
'{0}:{1:.6f}'.format(param_names[0], np.mean(params[0][i]))
228237
)
229238
)
230239
for i in range(n_sims)
231240
)
241+
232242

233243
# run all scenarios
234244
outputs = list(series.run(nproc=nproc).outputs_iter())
@@ -263,7 +273,8 @@ def pet(self, reference_pet_path, station_nhru, n_sims=10,\
263273
print('{0}\nOutput information sent to {1}\n'.format('-' * 80, json_outfile))
264274

265275
def plot_optimization(self, stage, freq='daily', method='time_series',\
266-
plot_vars='both', plot_1to1=True, return_fig=False):
276+
plot_vars='both', plot_1to1=True, return_fig=False,\
277+
n_plots=4):
267278
"""
268279
Basic plotting of current optimization results with limited options.
269280
Plots measured, original simluated, and optimization simulated variabes
@@ -332,7 +343,8 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\
332343
n = len(self.pet_outputs) # number of simulations to plot
333344
else:
334345
raise ValueError('{} is not a valid optimization stage.'.format(stage))
335-
346+
# user defined number of subplots from first n_plots results
347+
if (n > n_plots): n = n_plots
336348
# styles for each plot
337349
ms = 4 # markersize for all points
338350
orig_sty = dict(linestyle='none',markersize=ms,\
@@ -358,7 +370,7 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\
358370
fig, ax = plt.subplots(n, sharex=True, sharey=True,\
359371
figsize=(12,n*3.5))
360372
axs = ax.ravel()
361-
for i,sim in enumerate(sims):
373+
for i,sim in enumerate(sims[:n]):
362374
if plot_vars in ('meas', 'both'):
363375
axs[i].plot(meas, label='Measured', **meas_sty)
364376
if plot_vars in ('orig', 'both'):
@@ -380,7 +392,7 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\
380392

381393
fig, ax = plt.subplots(nrows=nrow, ncols=ncol, figsize=(12,n*3.5))
382394
axs = ax.ravel()
383-
for i,sim in enumerate(sims):
395+
for i,sim in enumerate(sims[:n]):
384396
if plot_vars in ('meas', 'both'):
385397
axs[i].plot(meas, label='Measured', **meas_sty)
386398
if plot_vars in ('orig', 'both'):
@@ -403,7 +415,7 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\
403415
meas_min = min(X)
404416
meas_max = max(X)
405417

406-
for i, sim in enumerate(sims):
418+
for i, sim in enumerate(sims[:n]):
407419
Y = sim
408420
sim_max = max(Y)
409421
sim_min = min(Y)
@@ -465,7 +477,7 @@ def _create_metafile_name(out_dir, opt_title, stage):
465477
name = '{}_{}_opt{}.json'.format(opt_title, stage, n)
466478
return name
467479

468-
def resample_param(params, param_name, noise_factor=0.1):
480+
def resample_param(params, param_name, how='uniform', noise_factor=0.1):
469481
"""
470482
Resample PRMS parameter by shifting all values by a constant that is
471483
taken from a uniform distribution, where the range of the uniform
@@ -479,6 +491,9 @@ def resample_param(params, param_name, noise_factor=0.1):
479491
params (parameters.Parameters): parameter object
480492
param_name (str): name of PRMS parameter to resample
481493
Kwargs:
494+
how (str): distribution to resample parameters from in the case
495+
that each parameter element can be resampled (len <=366)
496+
Currently works for uniform and normal distributions.
482497
noise_factor (float): factor to multiply parameter range by,
483498
use the result as the standard deviation for the normal rand.
484499
variable used to add element wise noise. i.e. higher
@@ -539,20 +554,35 @@ def resample_param(params, param_name, noise_factor=0.1):
539554

540555
low_bnd = p_min - np.min(param) # lowest param value minus allowable min
541556
up_bnd = p_max - np.max(param)
542-
s = (p_max - p_min) * noise_factor # stddev noise, default: range*(1/10)
557+
s = (p_max - p_min) * noise_factor # variance noise, default: range*(1/10)
543558
#do resampling differently based on param dimensions
544559
if dim_case == 'resample_all_values_once':
545560
#uniform RV for shifting all values once
546561
shifted_param = np.random.uniform(low=low_bnd, high=up_bnd) + param
547562
ret = shifted_param
548563
elif dim_case == 'resample_each_value':
549-
shifted_param = np.random.uniform(low=low_bnd, high=up_bnd) + param
550-
while True:
551-
## add noise to each value from ~N(0,s)
552-
tmp = shifted_param + np.random.normal(0,s,size=(np.shape(param)))
553-
if np.max(tmp) <= p_max and np.min(tmp) >= p_min:
554-
ret = tmp
555-
break
564+
ret = copy(param)
565+
if how == 'uniform':
566+
for i, el in enumerate(param):
567+
while True:
568+
low_bnd = p_min - np.min(el)
569+
up_bnd = p_max - np.max(el)
570+
571+
tmp = el + np.random.uniform(low=low_bnd, high=up_bnd)
572+
if np.max(tmp) <= p_max and np.min(tmp) >= p_min:
573+
ret[i] = tmp
574+
break
575+
elif how == 'normal':
576+
for i, el in enumerate(param):
577+
while True:
578+
low_bnd = p_min - np.min(el)
579+
up_bnd = p_max - np.max(el)
580+
581+
tmp = el + np.random.normal(0, s)
582+
if np.max(tmp) <= p_max and np.min(tmp) >= p_min:
583+
ret[i] = tmp
584+
break
585+
556586
elif dim_case == 'nhru_nmonths':
557587
ret = copy(param)
558588
rvs = [np.random.uniform(low=low_bnd, high=up_bnd) for i in range(12)]
@@ -561,18 +591,12 @@ def resample_param(params, param_name, noise_factor=0.1):
561591

562592
return ret
563593

564-
def _mod_params(parameters, params, stage, module):
594+
def _mod_params(parameters, params, param_names):
565595
# deepcopy was crashing, raising:
566596
# TypeError: cannot serialize '_io.TextIOWrapper' object
567597
ret = copy(parameters)
568-
#print (intcp, slope)
569-
if ((stage == 'swrad') and (module=='ddsolrad')):
570-
ret['dday_intcp'] = params[0]
571-
ret['dday_slope'] = params[1]
572-
elif ((stage == 'pet') and (module=='potet_jh')):
573-
ret['jh_coef'] = params[0]
574-
elif ((stage == 'pet') and (module=='potet_pt')):
575-
ret['potet_coef_hru_mo'] = params[0]
598+
for idx, param in enumerate(params):
599+
ret[param_names[idx]] = param
576600
return ret
577601

578602

@@ -582,6 +606,10 @@ def __init__(self, working_dir, stage='all'):
582606
self.working_dir = working_dir
583607
self.stage = stage
584608
self.metadata_json_paths = self.get_optr_jsons(working_dir, stage)
609+
self.meas_swrad = self.get_measured('swrad')
610+
self.meas_pet = self.get_measured('pet')
611+
self.swrad_statvar_name = self.get_statvar_name('swrad')
612+
self.pet_statvar_name = self.get_statvar_name('pet')
585613

586614
def get_optr_jsons(self, work_dir, stage):
587615
"""
@@ -617,15 +645,110 @@ def get_optr_jsons(self, work_dir, stage):
617645
optr_metafile_re.match(f) ]
618646
return ret
619647

648+
def get_sim_dirs(self, stage):
649+
jsons = self.metadata_json_paths[stage]
650+
json_files = []
651+
sim_dirs = []
652+
for inf in jsons:
653+
with open(inf) as json_file:
654+
json_files.append(json.load(json_file))
655+
for json_file in json_files:
656+
sim_dirs.extend(json_file['sim_dirs'])
657+
# list of all simulation directory paths for stage
658+
return sim_dirs
659+
660+
def get_measured(self, stage):
661+
# only need to open one json file to get this information
662+
if not self.metadata_json_paths.get(stage):
663+
return # no optimization json files exist for given stage
664+
first_json = self.metadata_json_paths[stage][0]
665+
with open(first_json) as json_file:
666+
json_data = json.load(json_file)
667+
measured_series = pd.Series.from_csv(json_data.get('measured_{}'.\
668+
format(stage)), parse_dates=True)
669+
return measured_series
670+
671+
def get_statvar_name(self, stage):
672+
# only need to open one json file to get this information
673+
if not self.metadata_json_paths.get(stage):
674+
return # no optimization json files exist for given stage
675+
first_json = self.metadata_json_paths[stage][0]
676+
with open(first_json) as json_file:
677+
json_data = json.load(json_file)
678+
var_name = json_data.get('{}_hru_id'.format(stage))
679+
680+
return var_name
681+
682+
def result_table(self, stage, freq='daily', top_n=5, latex=False):
683+
##TODO: add stats for freq options monthly, annual (means or sum)
684+
685+
sim_dirs = self.get_sim_dirs(stage)
686+
if top_n >= len(sim_dirs): top_n = len(sim_dirs) - 1
687+
sim_names = [path.split(os.sep)[-1] for path in sim_dirs]
688+
meas_var = self.get_measured(stage)
689+
statvar_name = self.get_statvar_name(stage)
690+
result_df = pd.DataFrame(columns=\
691+
['NSE','RMSE','PBIAS','COEF_DET','ABS(PBIAS)'])
692+
for i, sim in enumerate(sim_dirs):
693+
sim_out = load_statvar(OPJ(sim, 'outputs', 'statvar.dat'))\
694+
['{}'.format(statvar_name)]
695+
idx = meas_var.index.intersection(sim_out.index)
696+
meas_var = copy(meas_var[idx])
697+
sim_out = sim_out[idx]
698+
if freq == 'daily':
699+
result_df.loc[sim_names[i]] = [nash_sutcliffe(meas_var, sim_out),\
700+
rmse(meas_var, sim_out),\
701+
percent_bias(meas_var,sim_out),\
702+
meas_var.corr(sim_out)**2,\
703+
np.abs(percent_bias(meas_var, sim_out)) ]
704+
elif freq == 'monthly':
705+
meas_mo = meas_var.groupby(meas_var.index.month).mean()
706+
sim_out = sim_out.groupby(sim_out.index.month).mean()
707+
result_df.loc[sim_names[i]] = [nash_sutcliffe(meas_mo, sim_out),\
708+
rmse(meas_mo, sim_out),\
709+
percent_bias(meas_mo, sim_out),\
710+
meas_mo.corr(sim_out)**2,\
711+
np.abs(percent_bias(meas_mo, sim_out)) ]
712+
713+
sorted_result = result_df.sort_values(by=['NSE','RMSE','ABS(PBIAS)',\
714+
'COEF_DET'], ascending=[False,True,True,False])
715+
sorted_result.columns.name = '{} parameters'.format(stage)
716+
sorted_result = sorted_result[['NSE','RMSE','PBIAS','COEF_DET']]
717+
718+
if latex: return sorted_result[:top_n].to_latex(escape=False)
719+
else: return sorted_result[:top_n]
720+
721+
def get_top_ranked_sims(self, sorted_df):
722+
## use result table to make dic with best param and statvar paths
723+
# index of table is the simulation directory names
724+
ret = {
725+
'dir_name' : [],
726+
'param_path' : [],
727+
'statvar_path' : []
728+
}
729+
730+
for i,el in enumerate(sorted_df.index):
731+
ret['dir_name'].append(el)
732+
ret['param_path'].append(OPJ(self.working_dir,el,'inputs','parameters'))
733+
ret['statvar_path'].append(OPJ(self.working_dir,el,'outputs','statvar.dat'))
734+
735+
return ret
736+
737+
620738
class SradOptimizationResult(OptimizationResult):
621739

622740
def __init__(self, working_dir, stage='swrad' ):
623741
OptimizationResult.__init__(self, working_dir, stage)
624-
742+
self.swrad_sim_dirs = self.get_sim_dirs(stage)
743+
self.meas_swrad = self.get_measured(stage)
744+
self.swrad_statvar_name = self.get_statvar_name(stage)
625745

626746

627747
class PetOptimizationResult(OptimizationResult):
628748

629749
def __init__(self, working_dir, stage='pet'):
630750
OptimizationResult.__init__(self, working_dir, stage)
751+
self.pet_sim_dirs = self.get_sim_dirs(stage)
752+
self.meas_pet = self.get_measured(stage)
753+
self.pet_statvar_name = self.get_statvar_name(stage)
631754

0 commit comments

Comments
 (0)