|
15 | 15 | from .data import Data |
16 | 16 | from .parameters import Parameters |
17 | 17 | from .simulation import Simulation, SimulationSeries |
18 | | -from .util import load_statvar |
| 18 | +from .util import load_statvar, nash_sutcliffe, percent_bias, rmse |
19 | 19 |
|
20 | 20 | OPJ = os.path.join |
21 | 21 |
|
@@ -43,8 +43,8 @@ class Optimizer: |
43 | 43 | #dic for min/max of parameter allowable ranges, add more when needed |
44 | 44 | param_ranges = {'dday_intcp': (-60.0, 10.0), 'dday_slope': (0.2, 0.9),\ |
45 | 45 | '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 |
48 | 48 |
|
49 | 49 | def __init__(self, parameters, data, control_file, working_dir, |
50 | 50 | title, description=None): |
@@ -84,7 +84,7 @@ def __init__(self, parameters, data, control_file, working_dir, |
84 | 84 | self.pet_hru = None |
85 | 85 |
|
86 | 86 | 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): |
88 | 88 | ''' |
89 | 89 | Optimize the monthly dday_intcp and dday_slope parameters |
90 | 90 | (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='',\ |
118 | 118 | ## resampling degree day slope and intercept simoultaneously |
119 | 119 | if srad_mod == 'ddsolrad': |
120 | 120 | 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)) |
124 | 129 |
|
125 | 130 | self.data.write(OPJ(self.working_dir, 'data')) |
126 | 131 |
|
127 | | - # TODO: better file naming for output folders |
| 132 | + # sim dirs named after first resampled param name and mean value |
128 | 133 | series = SimulationSeries( |
129 | 134 | 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), |
132 | 138 | self.control_file, |
133 | 139 | OPJ( |
134 | 140 | 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])) |
137 | 142 | ) |
138 | 143 | ) |
139 | 144 | for i in range(n_sims) |
@@ -177,7 +182,7 @@ def _error(x, y): |
177 | 182 | print('{0}\nOutput information sent to {1}\n'.format('-' * 80, json_outfile)) |
178 | 183 |
|
179 | 184 | 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): |
181 | 186 | ''' |
182 | 187 | Optimize the monthly coefficients (depending on module) parameters |
183 | 188 | 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,\ |
203 | 208 | else: # get variable at a specific hru |
204 | 209 | self.pet_hru = 'potet_{}'.format(station_hru) |
205 | 210 |
|
| 211 | + pet_start_time = dt.datetime.now() |
| 212 | + pet_start_time = pet_start_time.replace(second=0, microsecond=0) |
| 213 | + |
206 | 214 | if pet_mod == 'potet_pt': |
207 | 215 | param_names = ['potet_coef_hru_mo'] |
208 | | - params = [resample_param(self.parameters, 'potet_coef_hru_mo')\ |
209 | | - for i in range(n_sims)] |
210 | 216 | elif pet_mod == 'potet_jh': |
211 | 217 | 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 |
221 | 228 | series = SimulationSeries( |
222 | 229 | 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), |
224 | 233 | self.control_file, |
225 | 234 | OPJ( |
226 | 235 | 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])) |
228 | 237 | ) |
229 | 238 | ) |
230 | 239 | for i in range(n_sims) |
231 | 240 | ) |
| 241 | + |
232 | 242 |
|
233 | 243 | # run all scenarios |
234 | 244 | outputs = list(series.run(nproc=nproc).outputs_iter()) |
@@ -263,7 +273,8 @@ def pet(self, reference_pet_path, station_nhru, n_sims=10,\ |
263 | 273 | print('{0}\nOutput information sent to {1}\n'.format('-' * 80, json_outfile)) |
264 | 274 |
|
265 | 275 | 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): |
267 | 278 | """ |
268 | 279 | Basic plotting of current optimization results with limited options. |
269 | 280 | Plots measured, original simluated, and optimization simulated variabes |
@@ -332,7 +343,8 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\ |
332 | 343 | n = len(self.pet_outputs) # number of simulations to plot |
333 | 344 | else: |
334 | 345 | 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 |
336 | 348 | # styles for each plot |
337 | 349 | ms = 4 # markersize for all points |
338 | 350 | orig_sty = dict(linestyle='none',markersize=ms,\ |
@@ -358,7 +370,7 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\ |
358 | 370 | fig, ax = plt.subplots(n, sharex=True, sharey=True,\ |
359 | 371 | figsize=(12,n*3.5)) |
360 | 372 | axs = ax.ravel() |
361 | | - for i,sim in enumerate(sims): |
| 373 | + for i,sim in enumerate(sims[:n]): |
362 | 374 | if plot_vars in ('meas', 'both'): |
363 | 375 | axs[i].plot(meas, label='Measured', **meas_sty) |
364 | 376 | if plot_vars in ('orig', 'both'): |
@@ -380,7 +392,7 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\ |
380 | 392 |
|
381 | 393 | fig, ax = plt.subplots(nrows=nrow, ncols=ncol, figsize=(12,n*3.5)) |
382 | 394 | axs = ax.ravel() |
383 | | - for i,sim in enumerate(sims): |
| 395 | + for i,sim in enumerate(sims[:n]): |
384 | 396 | if plot_vars in ('meas', 'both'): |
385 | 397 | axs[i].plot(meas, label='Measured', **meas_sty) |
386 | 398 | if plot_vars in ('orig', 'both'): |
@@ -403,7 +415,7 @@ def plot_optimization(self, stage, freq='daily', method='time_series',\ |
403 | 415 | meas_min = min(X) |
404 | 416 | meas_max = max(X) |
405 | 417 |
|
406 | | - for i, sim in enumerate(sims): |
| 418 | + for i, sim in enumerate(sims[:n]): |
407 | 419 | Y = sim |
408 | 420 | sim_max = max(Y) |
409 | 421 | sim_min = min(Y) |
@@ -465,7 +477,7 @@ def _create_metafile_name(out_dir, opt_title, stage): |
465 | 477 | name = '{}_{}_opt{}.json'.format(opt_title, stage, n) |
466 | 478 | return name |
467 | 479 |
|
468 | | -def resample_param(params, param_name, noise_factor=0.1): |
| 480 | +def resample_param(params, param_name, how='uniform', noise_factor=0.1): |
469 | 481 | """ |
470 | 482 | Resample PRMS parameter by shifting all values by a constant that is |
471 | 483 | 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): |
479 | 491 | params (parameters.Parameters): parameter object |
480 | 492 | param_name (str): name of PRMS parameter to resample |
481 | 493 | 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. |
482 | 497 | noise_factor (float): factor to multiply parameter range by, |
483 | 498 | use the result as the standard deviation for the normal rand. |
484 | 499 | variable used to add element wise noise. i.e. higher |
@@ -539,20 +554,35 @@ def resample_param(params, param_name, noise_factor=0.1): |
539 | 554 |
|
540 | 555 | low_bnd = p_min - np.min(param) # lowest param value minus allowable min |
541 | 556 | 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) |
543 | 558 | #do resampling differently based on param dimensions |
544 | 559 | if dim_case == 'resample_all_values_once': |
545 | 560 | #uniform RV for shifting all values once |
546 | 561 | shifted_param = np.random.uniform(low=low_bnd, high=up_bnd) + param |
547 | 562 | ret = shifted_param |
548 | 563 | 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 | + |
556 | 586 | elif dim_case == 'nhru_nmonths': |
557 | 587 | ret = copy(param) |
558 | 588 | 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): |
561 | 591 |
|
562 | 592 | return ret |
563 | 593 |
|
564 | | -def _mod_params(parameters, params, stage, module): |
| 594 | +def _mod_params(parameters, params, param_names): |
565 | 595 | # deepcopy was crashing, raising: |
566 | 596 | # TypeError: cannot serialize '_io.TextIOWrapper' object |
567 | 597 | 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 |
576 | 600 | return ret |
577 | 601 |
|
578 | 602 |
|
@@ -582,6 +606,10 @@ def __init__(self, working_dir, stage='all'): |
582 | 606 | self.working_dir = working_dir |
583 | 607 | self.stage = stage |
584 | 608 | 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') |
585 | 613 |
|
586 | 614 | def get_optr_jsons(self, work_dir, stage): |
587 | 615 | """ |
@@ -617,15 +645,110 @@ def get_optr_jsons(self, work_dir, stage): |
617 | 645 | optr_metafile_re.match(f) ] |
618 | 646 | return ret |
619 | 647 |
|
| 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 | + |
620 | 738 | class SradOptimizationResult(OptimizationResult): |
621 | 739 |
|
622 | 740 | def __init__(self, working_dir, stage='swrad' ): |
623 | 741 | 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) |
625 | 745 |
|
626 | 746 |
|
627 | 747 | class PetOptimizationResult(OptimizationResult): |
628 | 748 |
|
629 | 749 | def __init__(self, working_dir, stage='pet'): |
630 | 750 | 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) |
631 | 754 |
|