+# -*- coding: utf-8 -*-
+'''
+optimizer.py -- holds ``Optimizer`` and ``OptimizationResult`` classes for
+optimization routines and management conducted on PRMS parameters.
+'''
+from __future__ import print_function
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import os, sys, json, re, shutil
+import multiprocessing as mp
+from copy import copy
+from copy import deepcopy
+from numpy import log10
+from datetime import datetime
+from .data import Data
+from .parameters import Parameters
+from .simulation import Simulation, SimulationSeries
+from .util import load_statvar, nash_sutcliffe, percent_bias, rmse
+
+OPJ = os.path.join
+
+
+[docs]class Optimizer:
+
'''
+
Container for PRMS parameter optimization and related routines.
+
+
Currently the ``monte_carlo`` method provides random parameter
+
resampling routines using uniform and normal random variables.
+
+
Example:
+
>>> from prms_python import Data, Optimizer, Parameters
+
>>> params = Parameters('path/to/parameters')
+
>>> data = Data('path/to/data')
+
>>> control = 'path/to/control'
+
>>> work_directory = 'path/to/create/simulations'
+
>>> optr = Optimizer(
+
params,
+
data,
+
control,
+
work_directory,
+
title='the title',
+
description='desc')
+
>>> measured = 'path/to/measured/csv'
+
>>> statvar_name = 'basin_cfs' # or any other valid statvar
+
>>> params_to_resample = ['dday_intcp', 'dday_slope'] # list of params
+
>>> optr.monte_carlo(measured, params_to_resample, statvar_name)
+
+
'''
+
+
#dic for min/max of parameter allowable ranges, add more when needed
+
param_ranges = {'dday_intcp': (-60.0, 10.0),
+
'dday_slope': (0.2, 0.9),
+
'jh_coef': (0.005, 0.06),
+
'pt_alpha': (1.0, 2.0),
+
'potet_coef_hru_mo': (1.0, 2.0),
+
'tmax_index': (-10.0, 110.0),
+
'tmin_lapse': (-10.0, 10.0),
+
'soil_moist_max': (0.001, 10.0),
+
'rain_adj': (0.5, 2.0),
+
'ppt_rad_adj': (0.0, 0.5),
+
'radadj_intcp': (0.0, 1.0),
+
'radadj_slope': (0.0, 1.0),
+
'radj_sppt': (0.0, 1.0),
+
'radj_wppt': (0.0, 1.0),
+
'radmax': (0.1, 1.0)
+
}
+
+
def __init__(self, parameters, data, control_file, working_dir,
+
title, description=None):
+
+
if isinstance(parameters, Parameters):
+
self.parameters = parameters
+
else:
+
raise TypeError('parameters must be instance of Parameters')
+
+
if isinstance(data, Data):
+
self.data = data
+
+
else:
+
raise TypeError('data must be instance of Data')
+
+
input_dir = '{}'.format(os.sep).join(control_file.split(os.sep)[:-1])
+
if not os.path.isfile(OPJ(input_dir, 'statvar.dat')):
+
print('You have no statvar.dat file in your model directory')
+
print('Running PRMS on original data in {} for later comparison'\
+
.format(input_dir))
+
sim = Simulation(input_dir)
+
sim.run()
+
+
if not os.path.isdir(working_dir):
+
os.mkdir(working_dir)
+
+
self.input_dir = input_dir
+
self.control_file = control_file
+
self.working_dir = working_dir
+
self.title = title
+
self.description = description
+
self.measured_arb = None # for arbitrary output methods
+
self.statvar_name = None
+
self.arb_outputs = []
+
+
[docs] def monte_carlo(self, reference_path, param_names, statvar_name, \
+
stage, n_sims=10, method='uniform', mu_factor=1,\
+
noise_factor=0.1, nproc=None):
+
'''
+
The ``monte_carlo`` method of ``Optimizer`` performs parameter
+
random resampling techniques to a set of PRMS parameters and
+
executes and manages the corresponding simulations.
+
+
Arguments:
+
reference_path (str): path to measured data for optimization
+
param_names (list): list of parameter names to resample
+
statvar_name (str): name of statisical variable output name for
+
optimization
+
stage (str): custom name of optimization stage e.g. 'ddsolrad'
+
+
Keyword Arguments:
+
n_sims (int): number of simulations to conduct
+
parameter optimization/uncertaitnty analysis.
+
method (str): resampling method for parameters (normal or uniform)
+
mu_factor (float): coefficient to scale mean of the parameter(s)
+
to resample from when using the normal distribution to resample
+
i.e. a value of 1.5 will sample from a normal rv with mean
+
50% higher than the original parameter mean
+
noise_factor (float): scales the variance of noise to add to
+
parameter values when using normal rv (method='normal')
+
nproc (int): number of processors available to run PRMS simulations
+
+
Returns:
+
None
+
'''
+
if '_' in stage:
+
raise ValueError('stage name cannot contain an underscore')
+
# assign the optimization object a copy of measured data for plots
+
self.measured_arb = pd.Series.from_csv(reference_path,\
+
parse_dates=True)
+
# statistical variable output name
+
self.statvar_name = statvar_name
+
+
start_time = datetime.now().isoformat()
+
+
# resample params for all simulations- potential place to serialize
+
params = []
+
for name in param_names: # create list of lists of resampled params
+
tmp = []
+
for idx in range(n_sims):
+
tmp.append(resample_param(self.parameters, name, how=method,\
+
mu_factor=mu_factor, noise_factor=noise_factor))
+
params.append(list(tmp))
+
+
# SimulationSeries comprised of each resampled param set
+
series = SimulationSeries(
+
Simulation.from_data(
+
self.data, _mod_params(self.parameters,\
+
[params[n][i] for n in range(len(params))],\
+
param_names),
+
self.control_file,
+
OPJ(
+
self.working_dir, # name of sim: first param and mean value
+
'{0}_{1:.10f}'.format(param_names[0], np.mean(params[0][i]))
+
)
+
)
+
for i in range(n_sims)
+
)
+
+
if not nproc:
+
nproc = mp.cpu_count() // 2
+
+
# run
+
outputs = list(series.run(nproc=nproc).outputs_iter())
+
self.arb_outputs.extend(outputs) # for current instance- add outputs
+
+
end_time = datetime.now().isoformat()
+
+
# json metadata for Monte Carlo run
+
meta = { 'params_adjusted' : param_names,
+
'statvar_name' : self.statvar_name,
+
'optimization_title' : self.title,
+
'optimization_description' : self.description,
+
'start_datetime' : start_time,
+
'end_datetime' : end_time,
+
'measured' : reference_path,
+
'method' : 'Monte Carlo',
+
'mu_factor' : mu_factor,
+
'noise_factor' : noise_factor,
+
'resample': method,
+
'sim_dirs' : [],
+
'stage': stage,
+
'original_params' : self.parameters.base_file,
+
'nproc': nproc,
+
'n_sims' : n_sims
+
}
+
+
for output in outputs:
+
meta['sim_dirs'].append(output['simulation_dir'])
+
+
json_outfile = OPJ(self.working_dir, _create_metafile_name(\
+
self.working_dir, self.title, stage))
+
+
with open(json_outfile, 'w') as outf:
+
json.dump(meta, outf, sort_keys=True, indent=4,\
+
ensure_ascii=False)
+
+
print('{0}\nOutput information sent to {1}\n'.\
+
format('-' * 80, json_outfile))
+
+
[docs] def plot_optimization(self, freq='daily', method='time_series',\
+
plot_vars='both', plot_1to1=True, return_fig=False,\
+
n_plots=4):
+
"""
+
Basic plotting of current optimization results with limited options.
+
Plots measured, original simluated, and optimization simulated variabes
+
Not recommended for plotting results when n_sims is very large, instead
+
use options from an OptimizationResult object, or employ a user-defined
+
method using the result data.
+
+
Keyword Arguments:
+
freq (str): frequency of time series plots, value can be 'daily'
+
or 'monthly' for solar radiation
+
method (str): 'time_series' for time series sub plot of each
+
simulation alongside measured radiation. Another choice is
+
'correlation' which plots each measured daily solar radiation
+
value versus the corresponding simulated variable as subplots
+
one for each simulation in the optimization. With coefficients
+
of determiniationi i.e. square of pearson correlation coef.
+
plot_vars (str): what to plot alongside simulated srad:
+
'meas': plot simulated along with measured swrad
+
'orig': plot simulated along with the original simulated swrad
+
'both': plot simulated, with original simulation and measured
+
plot_1to1 (bool): if True plot one to one line on correlation
+
scatter plot, otherwise exclude.
+
return_fig (bool): flag whether to return matplotlib figure
+
+
Returns:
+
f (:obj:`matplotlib.figure.Figure`): If kwarg return_fig=True, then return
+
copy of the figure that is generated to the user.
+
"""
+
if not self.arb_outputs:
+
raise ValueError('You have not run any optimizations')
+
var_name = self.statvar_name
+
X = self.measured_arb
+
idx = X.index.intersection(load_statvar(self.arb_outputs[0]\
+
['statvar'])['{}'.format(var_name)].index)
+
X = X[idx]
+
orig = load_statvar(OPJ(self.input_dir, 'statvar.dat'))['{}'\
+
.format(var_name)][idx]
+
meas = self.measured_arb[idx]
+
sims = [load_statvar(out['statvar'])['{}'.format(var_name)][idx] for \
+
out in self.arb_outputs]
+
simdirs = [out['simulation_dir'].split(os.sep)[-1].\
+
replace('_', ' ') for out in self.arb_outputs]
+
var_name = '{}'.format(self.statvar_name)
+
n = len(self.arb_outputs) # number of simulations to plot
+
# user defined number of subplots from first n_plots results
+
if (n > n_plots): n = n_plots
+
# styles for each plot
+
ms = 4 # markersize for all points
+
orig_sty = dict(linestyle='none',markersize=ms,\
+
markerfacecolor='none', marker='s',\
+
markeredgecolor='royalblue', color='royalblue')
+
meas_sty = dict(linestyle='none',markersize=ms+1,\
+
markerfacecolor='none', marker='1',\
+
markeredgecolor='k', color='k')
+
sim_sty = dict(linestyle='none',markersize=ms,\
+
markerfacecolor='none', marker='o',\
+
markeredgecolor='r', color='r')
+
## number of subplots and rows (two plots per row)
+
nrow = n//2 # round down if odd n
+
ncol = 2
+
odd_n = False
+
if n/2. - nrow == 0.5:
+
nrow+=1 # odd number need extra row
+
odd_n = True
+
########
+
## Start plots depnding on key word arguments
+
########
+
if freq == 'daily' and method == 'time_series':
+
fig, ax = plt.subplots(n, sharex=True, sharey=True,\
+
figsize=(12,n*3.5))
+
axs = ax.ravel()
+
for i,sim in enumerate(sims[:n]):
+
if plot_vars in ('meas', 'both'):
+
axs[i].plot(meas, label='Measured', **meas_sty)
+
if plot_vars in ('orig', 'both'):
+
axs[i].plot(orig, label='Original sim.', **orig_sty)
+
axs[i].plot(sim, **sim_sty)
+
axs[i].set_ylabel('sim: {}'.format(simdirs[i]), fontsize=10)
+
if i == 0: axs[i].legend(markerscale=5, loc='best')
+
fig.subplots_adjust(hspace=0)
+
fig.autofmt_xdate()
+
#monthly means
+
elif freq == 'monthly' and method == 'time_series':
+
# compute monthly means
+
meas = meas.groupby(meas.index.month).mean()
+
orig = orig.groupby(orig.index.month).mean()
+
# change line styles for monthly plots to lines not points
+
for d in (orig_sty, meas_sty, sim_sty):
+
d['linestyle'] = '-'
+
d['marker'] = None
+
+
fig, ax = plt.subplots(nrows=nrow, ncols=ncol, figsize=(12,n*3.5))
+
axs = ax.ravel()
+
for i,sim in enumerate(sims[:n]):
+
if plot_vars in ('meas', 'both'):
+
axs[i].plot(meas, label='Measured', **meas_sty)
+
if plot_vars in ('orig', 'both'):
+
axs[i].plot(orig, label='Original sim.', **orig_sty)
+
sim = sim.groupby(sim.index.month).mean()
+
axs[i].plot(sim, **sim_sty)
+
axs[i].set_ylabel('sim: {}\nmean'.format(simdirs[i]),\
+
fontsize=10)
+
axs[i].set_xlim(0.5,12.5)
+
if i == 0: axs[i].legend(markerscale=5, loc='best')
+
if odd_n: # empty subplot if odd number of simulations
+
fig.delaxes(axs[n])
+
fig.text(0.5, 0.1, 'month')
+
#x-y scatter
+
elif method == 'correlation':
+
## figure
+
fig, ax = plt.subplots(nrows=nrow, ncols=ncol, figsize=(12,n*3))
+
axs = ax.ravel()
+
## subplot dimensions
+
meas_min = min(X)
+
meas_max = max(X)
+
+
for i, sim in enumerate(sims[:n]):
+
Y = sim
+
sim_max = max(Y)
+
sim_min = min(Y)
+
m = max(meas_max,sim_max)
+
if plot_1to1:
+
axs[i].plot([0, m], [0, m], 'k--', lw=2) ## one to one line
+
axs[i].set_xlim(meas_min,meas_max)
+
axs[i].set_ylim(sim_min, sim_max)
+
axs[i].plot(X, Y, **sim_sty)
+
axs[i].set_ylabel('sim: {}'.format(simdirs[i]))
+
axs[i].set_xlabel('Measured {}'.format(var_name))
+
axs[i].text(0.05, 0.95,r'$R^2 = {0:.2f}$'.format(\
+
X.corr(Y)**2), fontsize=16,\
+
ha='left', va='center', transform=axs[i].transAxes)
+
if odd_n: # empty subplot if odd number of simulations
+
fig.delaxes(axs[n])
+
+
if return_fig:
+
return fig
+
+def _create_metafile_name(out_dir, opt_title, stage):
+ """
+ Search through output directory where simulations are conducted
+ look for all metadata simulation json files and find out if the
+ current simulation is a replicate. Then use that information to
+ build the correct file name for the output json file. The series
+ are typically run in parallel that is why this step has to be
+ done after running multiple simulations from an optimization stage.
+
+ Arguments:
+ out_dir (str): path to directory with model results, i.e.
+ location where simulation series outputs and optimization
+ json files are located, aka Optimizer.working_dir
+ opt_title (str): optimization instance title for file search
+ stage (str): stage of optimization, e.g. 'swrad', 'pet'
+
+ Returns:
+ name (str): file name for the current optimization simulation series
+ metadata json file. E.g 'dry_creek_swrad_opt.json', or if
+ this is the second time you have run an optimization titled
+ 'dry_creek' the next json file will be returned as
+ 'dry_creek_swrad_opt1.json' and so on with integer increments
+ """
+ meta_re = re.compile(r'^{}_{}_opt(\d*)\.json'.format(opt_title, stage))
+ reps = []
+ for f in os.listdir(out_dir):
+ if meta_re.match(f):
+ nrep = meta_re.match(f).group(1)
+ if nrep == '':
+ reps.append(0)
+ else:
+ reps.append(nrep)
+
+ if not reps:
+ name = '{}_{}_opt.json'.format(opt_title, stage)
+ else:
+ # this is the nth optimization done under the same title
+ n = max(map(int, reps)) + 1
+ name = '{}_{}_opt{}.json'.format(opt_title, stage, n)
+ return name
+
+[docs]def resample_param(params, param_name, how='uniform', mu_factor=1,\
+
noise_factor=0.1):
+
"""
+
Resample PRMS parameter by shifting all values by a constant that is
+
taken from a uniform distribution, where the range of the uniform
+
values is equal to the difference between the min and max of the allowable
+
range. The parameter min and max are set in ``Optimizer.param_ranges``.
+
If the resampling method (``how`` argument) is set to 'normal', randomly
+
sample a normal distribution with mean = mean(parameter) X ``mu_factor`` and
+
sigma = param allowable range multiplied by ``noise_factor``. If parameters have
+
array length <= 366 then individual parameter values are resampled otherwise
+
resample all param values at once, e.g. by taking a single random value
+
from the uniform distribution. If they are taking all at once using the
+
normal method then the original values are scaled by mu_factor and a normal
+
random variable with mean=0 and std dev = parameter range X ``noise_factor``.
+
+
Arguments:
+
params (:class:`prms_python.Parameters`): ``Parameters`` object
+
param_name (str): name of PRMS parameter to resample
+
+
Keyword Arguments:
+
how (str): distribution to resample parameters from in the case
+
that each parameter element can be resampled (len <=366)
+
Currently works for uniform and normal distributions.
+
noise_factor (float): factor to multiply parameter range by,
+
use the result as the standard deviation for the normal rand.
+
variable used to add element wise noise. i.e. higher
+
noise_factor will result in higher variance. Must be > 0.
+
+
Returns:
+
ret (:obj:`numpy.ndarry`): ndarray of param after resampling
+
+
Raises:
+
KeyError: if ``param_name`` not a valid parameter name
+
ValueError: if the parameter range has not been set in
+
``Optimizer.param_ranges``
+
"""
+
p_min, p_max = Optimizer.param_ranges.get(param_name,(-1,-1))
+
+
# create dictionary of parameter basic info (not values)
+
param_dic = {param['name']: param for param in params.base_params}
+
if not param_dic.get(param_name):
+
raise KeyError('{} is not a valid parameter'.format(param_name))
+
+
if p_min == p_max == -1:
+
raise ValueError("""{} has not been added to the dictionary of
+
parameters to resample, add it's allowable min and max value
+
to the Optimizer.param_ranges attribute in
+
Optimizer.py""".format(param_name))
+
+
dim_case = None
+
nhru = params.dimensions['nhru']
+
ndims = param_dic.get(param_name)['ndims']
+
dimnames = param_dic.get(param_name)['dimnames']
+
length = param_dic.get(param_name)['length']
+
param = deepcopy(params[param_name])
+
+
# could expand list and check parameter name also e.g. cascade_flg
+
# is a parameter that should not be changed
+
dims_to_not_change = set(['ncascade','ncascdgw','nreach',\
+
'nsegment'])
+
if (len(set.intersection(dims_to_not_change, set(dimnames))) > 0):
+
raise ValueError("""{} should not be resampled as
+
it relates to the location of cascade flow
+
parameters.""".format(param_name))
+
+
# use param info to get dimension info- e.g. if multidimensional
+
if (ndims == 1 and length <= 366):
+
dim_case = 'resample_each_value' # for smaller param dimensions
+
elif (ndims == 1 and length > 366):
+
dim_case = 'resample_all_values_once' # covers nssr, ngw, etc.
+
elif (ndims == 2 and dimnames[1] == 'nmonths' and \
+
nhru == params.dimensions[dimnames[0]]):
+
dim_case = 'nhru_nmonths'
+
elif not dim_case:
+
raise ValueError('The {} parameter is not set for resampling'.\
+
format(param_name))
+
# #testing purposes
+
# print('name: ', param_name)
+
# print('max_val: ', p_max)
+
# print('min_val: ', p_min)
+
# print('ndims: ', ndims)
+
# print('dimnames: ', dimnames)
+
# print('length: ', length)
+
# print('resample_method: ', dim_case)
+
+
s = (p_max - p_min) * noise_factor # std_dev (s) default: param_range/10
+
#do resampling based on param dimensions and sampling distribution
+
if dim_case == 'resample_all_values_once':
+
if how == 'uniform':
+
ret = np.random.uniform(low=p_min, high=p_max, size=param.shape)
+
elif how == 'normal': # scale parameter mean if mu_factor given
+
param *= mu_factor
+
tmp = np.random.normal(0, s, size=param.shape)
+
ret = tmp + param
+
+
elif dim_case == 'resample_each_value':
+
ret = param
+
if how == 'uniform':
+
ret = np.random.uniform(low=p_min, high=p_max, size=param.shape)
+
elif how == 'normal': # the original value is considered the mean
+
if len(ret.shape) != 0:
+
for i, el in enumerate(param):
+
mu = el * mu_factor
+
ret[i] = np.random.normal(mu, s)
+
else: # single value parameter
+
mu = param * mu_factor
+
ret = np.random.normal(mu, s)
+
+
# nhru by nmonth dimensional params
+
elif dim_case == 'nhru_nmonths':
+
ret = param
+
if how == 'uniform':
+
for month in range(12):
+
ret[month] = np.random.uniform(low=p_min, high=p_max,\
+
size=param[0].shape)
+
elif how == 'normal':
+
for i, el in enumerate(param):
+
el *= mu_factor
+
tmp = np.random.normal(0, s, size=el.shape)
+
ret[i] = tmp + el
+
+
return ret
+
+def _mod_params(parameters, params, param_names):
+ # loop through list of params and assign their values to Parameter instance
+ ret = copy(parameters)
+ for idx, param in enumerate(params):
+ ret[param_names[idx]] = np.array(param)
+ return ret
+
+
+[docs]class OptimizationResult:
+
"""
+
The ``OptimizationResult`` object serves to collect and manage output
+
from an ``Optimizer`` method. Upon initialization and a given optimization
+
stage that was used when running the Optimizer method, e.g. ``monte_carlo``,
+
the class gathers all JSON metadata that was produced for the given stage.
+
The ``OptimizationResult`` has three main user methods: first ``result_table``
+
which returns the top n simulations according to four model performance
+
metrics (Nash-Sutcliffe efficiency (NSE), root-mean squared-error (RMSE),
+
percent bias (PBIAS), and the coefficient of determination (COEF_DET) as
+
calculated against measured data. For example the table may look like:
+
+
>>> ddsolrad_res = OptimizationResult(work_directory, stage=stage)
+
>>> top10 = ddsolrad_res.result_table(freq='monthly',top_n=10)
+
>>> top10
+
======================== ======== ======= ========= ========
+
ddsolrad parameters NSE RMSE PBIAS COEF_DET
+
======================== ======== ======= ========= ========
+
orig_params 0.956267 39.4725 -0.885715 0.963116
+
tmax_index_54.2224631748 0.921626 47.6092 -0.849256 0.94402
+
tmax_index_44.8823940703 0.879965 58.9194 5.79603 0.922021
+
tmax_index_47.6835387480 0.764133 82.5918 -4.78896 0.837582
+
======================== ======== ======= ========= ========
+
+
Second, the ``get_top_ranked_sims`` which returns a dictionary that map
+
key information about the top n ranked simulations, an example returned
+
dictionary may look like:
+
+
>>> {
+
'dir_name' : ['pathToSim1', 'pathToSim2'],
+
'param_path' : ['pathToSim1/input/parameters', 'pathToSim2/input/parameters'],
+
'statvar_path' : ['pathToSim1/output/statvar.dat', 'pathToSim2/output/statvar.dat'],
+
'params_adjusted' : [[param_names_sim1], [param_names_sim2]]
+
}
+
+
The third method of ``OptimizationResult`` is ``archive`` which essentially
+
opens all parameter and statvar files from each simulation of the given
+
stage and archives the parameters that were modified and their modified values
+
and the statistical variable (PRMS time series output) that is associated with
+
the optimization stage. Other ``Optimizer`` simulation metadata is also gathered
+
and new JSON metadata containing only this information is created and written
+
within a newly created "archived" subdirectory within the same directory that
+
the ``Optimizer`` routine managed simulations. The ``OptimizationResult.archive``
+
method then recursively deletes the simulation data for each of the given stage.
+
"""
+
+
def __init__(self, working_dir, stage):
+
"""
+
Create an ``OptimizationResult`` instance to manage output and analyse parameter-
+
output relationships as produced by the use of an ``Optimizer`` method of a user
+
defined optimization stage.
+
"""
+
self.working_dir = working_dir
+
self.stage = stage
+
self.metadata_json_paths = self._get_optr_jsons(working_dir, stage)
+
self.total_sims = self._count_total_sims()
+
self.statvar_name = self._get_statvar_name(stage)
+
self.measured = self._get_measured(stage)
+
self.input_dir = self._get_input_dir(stage)
+
self.input_params = self._get_input_params(stage)
+
+
# if there are more than one input param for given stage
+
if len(self.input_params) > 1:
+
print(
+
"""Warning: there were more than one initial parameter sets used for the
+
optimization for stage: {}. Make sure to compare the the correct input
+
params with their corresponding output sims.""".format(stage))
+
print('\nThis optimization stage used the\
+
following input parameter files:\n{}'.format('\n'.join(self.input_params)))
+
+
def _count_total_sims(self):
+
# total number of simulations of given stage in working directory
+
tracked_dirs = []
+
for f in self.metadata_json_paths[self.stage]:
+
with open(f) as fh:
+
json_data = json.load(fh)
+
tracked_dirs.extend(json_data.get('sim_dirs'))
+
return len(tracked_dirs)
+
+
def _get_optr_jsons(self, work_dir, stage):
+
"""
+
Retrieve locations of optimization output jsons which contain
+
metadata needed to understand optimization results.
+
Create dictionary of each optimization with stage as key and lists
+
of corresponding json file paths as values.
+
+
Arguments:
+
work_dir (str): path to directory with model results, i.e.
+
location where simulation series outputs and optimization
+
json files are located, aka Optimizer.working_dir
+
stage (str): the stage ('ddsolrad', 'jhpet', 'flow', etc.) of
+
the optimization in which to gather the jsons
+
+
Returns:
+
ret (dict): dictionary of stage (keys) and lists of
+
json file paths for that stage (values).
+
"""
+
+
ret = {}
+
optr_metafile_re = re.compile(r'^.*_{}_opt(\d*)\.json'.\
+
format(stage))
+
ret[stage] = [OPJ(work_dir, f) for f in\
+
os.listdir(work_dir) if\
+
optr_metafile_re.match(f)]
+
return ret
+
+
def _get_input_dir(self, stage):
+
# retrieves the input directory from the first json file of given stage
+
json_file = self.metadata_json_paths[stage][0]
+
with open(json_file) as jf:
+
meta_dic = json.load(jf)
+
return '{}'.format(os.sep).join(meta_dic['original_params'].\
+
split(os.sep)[:-1])
+
+
def _get_input_params(self, stage):
+
json_files = self.metadata_json_paths[stage]
+
param_paths = []
+
for json_file in json_files:
+
with open(json_file) as jf:
+
meta_dic = json.load(jf)
+
param_paths.append(meta_dic['original_params'])
+
return list(set(param_paths))
+
+
def _get_sim_dirs(self, stage):
+
jsons = self.metadata_json_paths[stage]
+
json_files = []
+
sim_dirs = []
+
for inf in jsons:
+
with open(inf) as json_file:
+
json_files.append(json.load(json_file))
+
for json_file in json_files:
+
sim_dirs.extend(json_file['sim_dirs'])
+
# list of all simulation directory paths for stage
+
return sim_dirs
+
+
def _get_measured(self, stage):
+
# only need to open one json file to get this information
+
if not self.metadata_json_paths.get(stage):
+
return # no optimization json files exist for given stage
+
first_json = self.metadata_json_paths[stage][0]
+
with open(first_json) as json_file:
+
json_data = json.load(json_file)
+
measured_series = pd.Series.from_csv(json_data.get('measured'),\
+
parse_dates=True)
+
return measured_series
+
+
def _get_statvar_name(self, stage):
+
# only need to open one json file to get this information
+
try:
+
first_json = self.metadata_json_paths[stage][0]
+
except:
+
raise ValueError("""No optimization has been run for
+
stage: {}""".format(stage))
+
with open(first_json) as json_file:
+
json_data = json.load(json_file)
+
var_name = json_data.get('statvar_name')
+
+
return var_name
+
+
def result_table(self, freq='daily', top_n=5, latex=False):
+
##TODO: add stats for freq options annual (means or sum)
+
+
sim_dirs = self._get_sim_dirs(self.stage)
+
if top_n >= len(sim_dirs):
+
top_n = len(sim_dirs) + 1 # for returning inclusive last sim
+
sim_names = [path.split(os.sep)[-1] for path in sim_dirs]
+
meas_var = self._get_measured(self.stage)
+
statvar_name = self._get_statvar_name(self.stage)
+
orig_statvar = load_statvar(OPJ(self.input_dir,'statvar.dat'))\
+
[statvar_name]
+
+
result_df = pd.DataFrame(columns=\
+
['NSE','RMSE','PBIAS','COEF_DET','ABS(PBIAS)'])
+
orig_results = pd.DataFrame(index=['orig_params'],\
+
columns=['NSE','RMSE','PBIAS','COEF_DET'])
+
# get datetime indices that overlap from measured and simulated
+
sim_out = load_statvar(OPJ(sim_dirs[0], 'outputs', 'statvar.dat'))\
+
[statvar_name]
+
idx = meas_var.index.intersection(sim_out.index)
+
meas_var = copy(meas_var[idx])
+
#sim_out = sim_out[idx]
+
orig_statvar = orig_statvar[idx]
+
+
if freq == 'monthly':
+
meas_mo = meas_var.groupby(meas_var.index.month).mean()
+
orig_mo = orig_statvar.groupby(orig_statvar.index.month).mean()
+
+
for i, sim in enumerate(sim_dirs):
+
try:
+
sim_out = load_statvar(OPJ(sim, 'outputs', 'statvar.dat'))\
+
['{}'.format(statvar_name)]
+
except: # simulation might have been removed or missing
+
pass
+
+
sim_out = sim_out[idx]
+
if freq == 'daily':
+
result_df.loc[sim_names[i]] = [\
+
nash_sutcliffe(meas_var, sim_out),\
+
rmse(meas_var, sim_out),\
+
percent_bias(meas_var,sim_out),\
+
meas_var.corr(sim_out)**2,\
+
np.abs(percent_bias(meas_var, sim_out)) ]
+
orig_results.loc['orig_params'] = [\
+
nash_sutcliffe(orig_statvar,meas_var),\
+
rmse(orig_statvar,meas_var),\
+
percent_bias(orig_statvar,meas_var),\
+
orig_statvar.corr(meas_var)**2]
+
+
elif freq == 'monthly':
+
sim_out = sim_out.groupby(sim_out.index.month).mean()
+
result_df.loc[sim_names[i]] = [\
+
nash_sutcliffe(meas_mo, sim_out),\
+
rmse(meas_mo, sim_out),\
+
percent_bias(meas_mo, sim_out),\
+
meas_mo.corr(sim_out)**2,\
+
np.abs(percent_bias(meas_mo, sim_out)) ]
+
orig_results.loc['orig_params'] = [\
+
nash_sutcliffe(orig_mo,meas_mo),\
+
rmse(orig_mo,meas_mo),\
+
percent_bias(orig_mo,meas_mo),\
+
orig_mo.corr(meas_mo)**2]
+
+
sorted_result = result_df.sort_values(by=['NSE','RMSE','ABS(PBIAS)',\
+
'COEF_DET'], ascending=[False,True,True,False])
+
sorted_result.columns.name = '{} parameters'.format(self.stage)
+
sorted_result = sorted_result[['NSE','RMSE','PBIAS','COEF_DET']]
+
sorted_result = pd.concat([orig_results,sorted_result])
+
+
if latex: return sorted_result[:top_n].to_latex(escape=False)
+
else: return sorted_result[:top_n]
+
+
def get_top_ranked_sims(self, sorted_df):
+
# use result table to make dic with best param and statvar paths
+
# index of table are simulation directory names
+
ret = {
+
'dir_name' : [],
+
'param_path' : [],
+
'statvar_path' : [],
+
'params_adjusted' : []
+
}
+
+
json_paths = self.metadata_json_paths[self.stage]
+
+
for el in sorted_df.drop('orig_params').index:
+
ret['dir_name'].append(el)
+
ret['param_path'].append(OPJ(self.working_dir,el,'inputs',\
+
'parameters'))
+
ret['statvar_path'].append(OPJ(self.working_dir,el,'outputs',\
+
'statvar.dat'))
+
for f in json_paths:
+
with open(f) as fh:
+
json_data = json.load(fh)
+
if OPJ(self.working_dir, el) in json_data.get('sim_dirs'):
+
ret['params_adjusted'].append(\
+
json_data.get('params_adjusted'))
+
+
return ret
+
+
[docs] def archive(self, remove_sims=True, remove_meta=False, metric_freq='daily'):
+
"""
+
Create archive directory to hold json files that contain
+
information of adjusted parameters, model output, and performance
+
metrics for each Optimizer simulation of the
+
OptimizationResult.stage in the OptimizationResult.working_dir.
+
+
Keyword Arguments:
+
remove_sims (bool): If True recursively delete all folders
+
and files associated with original simulations of the
+
OptimizationResult.stage in the
+
OptimizationResult.working_dir, if False do not delete
+
simulations.
+
remove_meta (bool): Whether to delete original Optimizer
+
JSON metadata files, default is False.
+
metric_freq (Str): Frequency of output metric computation
+
for recording of model performance. Can be 'daily'
+
(default) or 'monthly'. Note, other results can be computed
+
later with archived results.
+
+
Returns:
+
None
+
"""
+
# create output archive directory
+
archive_dir = OPJ(self.working_dir,"{}_archived".format(self.stage))
+
if not os.path.isdir(archive_dir):
+
os.mkdir(archive_dir)
+
+
# create table and use to make mapping dic
+
table = self.result_table(freq=metric_freq,\
+
top_n=self.total_sims, latex=False)
+
map_dic = self.get_top_ranked_sims(table)
+
+
metadata_json_paths = self.metadata_json_paths[self.stage]
+
# get measured optimization variable path
+
first_json = metadata_json_paths[0]
+
+
with open(first_json) as json_file:
+
json_data = json.load(json_file)
+
measured_path = json_data.get('measured')
+
+
# record info for each simulation and archive to JSONs
+
# pandas series and numpy arrays are converted to Python lists
+
# for JSON serialization
+
for i, sim in enumerate(map_dic.get('dir_name')):
+
json_path = OPJ(archive_dir, '{sim}.json'.format(sim=sim))
+
try:
+
output_series = load_statvar(OPJ(self.working_dir, sim,\
+
'outputs', 'statvar.dat'))\
+
[self.statvar_name]
+
except: # simulation directory was already removed
+
continue
+
+
# look for resampling method info for the particular simulation
+
for f in metadata_json_paths:
+
with open(f) as tmp_file:
+
tmp = json.load(tmp_file)
+
if OPJ(self.working_dir, sim) in tmp.get('sim_dirs'):
+
resample = tmp.get('resample')
+
noise_factor = tmp.get('noise_factor')
+
mu_factor = tmp.get('mu_factor')
+
+
json_data = {
+
'param_names' : [],
+
'param_values' : [],
+
'original_param_path' : self.input_params,
+
'measured_path' : measured_path,
+
'output_name' : self.statvar_name,
+
'output_date_index' : output_series.\
+
index.astype(str).tolist(),
+
'output_values' : output_series.values.tolist(),
+
'metric_freq' : metric_freq,
+
'resample' : resample,
+
'stage' : self.stage,
+
'mu_factor' : mu_factor,
+
'noise_factor' : noise_factor,
+
'NSE' : table.loc[sim, 'NSE'],
+
'RMSE' : table.loc[sim, 'RMSE'],
+
'PBIAS' : table.loc[sim, 'PBIAS'],
+
'COEF_DET' : table.loc[sim, 'COEF_DET']
+
}
+
+
for param in map_dic.get('params_adjusted')[i]:
+
json_data['param_names'].append(param)
+
json_data['param_values'].append(Parameters(\
+
map_dic.get('param_path')[i])[param]\
+
.tolist())
+
# save JSON file into archive directory
+
with open(json_path, 'w') as outf:
+
json.dump(json_data, outf, sort_keys = True, indent = 4,\
+
ensure_ascii = False)
+
+
# recursively delete all simulation directories after archiving
+
if remove_sims:
+
path = OPJ(self.working_dir, sim)
+
for dirpath, dirnames, filenames in os.walk(path,\
+
topdown=False):
+
shutil.rmtree(dirpath, ignore_errors=True)
+
else:
+
continue
+
+
# optional delete the original JSON metadata
+
if remove_meta:
+
for meta_file in self.metadata_json_paths[self.stage]:
+
try:
+
os.remove(meta_file)
+
except:
+
continue
+
+
+
+
+
+
+
+
+
+