Skip to content

Commit 66bb693

Browse files
JohnVolkJohnVolk
committed
Renaming system for optimization metadata files, plotting options for optimization
1 parent 85916a5 commit 66bb693

1 file changed

Lines changed: 155 additions & 62 deletions

File tree

‎prms_python/optimizer.py‎

Lines changed: 155 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,12 @@
22
optimizer.py -- Optimization routines for PRMS parameters and data.
33
'''
44
from __future__ import print_function
5+
from itertools import tee, chain
56
import pandas as pd
67
import numpy as np
8+
import matplotlib.pyplot as plt
79
import datetime as dt
8-
import os, sys, json
10+
import os, sys, json, re
911

1012
from copy import deepcopy
1113
from numpy import log10
@@ -68,6 +70,9 @@ def __init__(self, parameters, data, control_file, working_dir,
6870
self.working_dir = working_dir
6971
self.title = title
7072
self.description = description
73+
self.srad_outputs = []
74+
self.measured_srad = None
75+
self.srad_hru = None
7176

7277
def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
7378
nproc=None):
@@ -90,6 +95,12 @@ def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
9095
Returns:
9196
(SradOptimizationResult)
9297
'''
98+
# assign the optimization object a copy of measured srad for plots
99+
self.measured_srad = pd.Series.from_csv(
100+
reference_srad_path, parse_dates=True
101+
)
102+
103+
self.srad_hru = station_nhru
93104

94105
srad_start_time = dt.datetime.now()
95106
srad_start_time = srad_start_time.replace(second=0, microsecond=0)
@@ -118,7 +129,8 @@ def srad(self, reference_srad_path, station_nhru, n_sims=10, method='',\
118129
)
119130

120131
# run all scenarios
121-
outputs = series.run(nproc=nproc).outputs_iter()
132+
outputs = list(series.run(nproc=nproc).outputs_iter())
133+
self.srad_outputs.extend(outputs)
122134

123135
srad_end_time = dt.datetime.now()
124136
srad_end_time = srad_end_time.replace(second=0, microsecond=0)
@@ -128,72 +140,149 @@ def _error(x, y):
128140
ret = sum(ret)
129141
return ret
130142

131-
measured_srad = pd.Series.from_csv(
132-
reference_srad_path, parse_dates=True
133-
)
134-
135143
srad_meta = {'stage' : 'swrad',
136-
'hru_id' : station_nhru,
144+
'swrad_hru_id' : self.srad_hru,
137145
'optimization_title' : self.title,
138146
'optimization_description' : self.description,
139147
'start_time' : str(srad_start_time),
140148
'end_time' : str(srad_end_time),
141149
'measured_swrad' : reference_srad_path,
142150
'sim_dirs' : [],
143-
'original_params' : self.parameters.base_file
151+
'original_params' : self.parameters.base_file,
152+
'n_sims' : n_sims
144153
}
145154

146155
for output in outputs:
147156
srad_meta['sim_dirs'].append(output['simulation_dir'])
148-
149-
json_outfile = OPJ(self.working_dir, '{0}_swrad_opt.json'.format(self.title))
150-
157+
158+
json_outfile = OPJ(self.working_dir, _create_metafile_name(\
159+
self.working_dir, self.title, 'swrad'))
160+
151161
with open(json_outfile, 'w') as outf:
152-
json.dump(srad_meta, outf, sort_keys = True, indent = 4, ensure_ascii = False)
162+
json.dump(srad_meta, outf, sort_keys = True, indent = 4,\
163+
ensure_ascii = False)
153164

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

156-
# # calculate the top performing
157-
# errors = (
158-
# (
159-
# output['simulation_dir'],
160-
# _error(measured_srad,
161-
# output['statvar']['swrad_' + str(station_nhru)])
162-
# )
163-
#
164-
# for output in outputs
165-
# )
166-
#
167-
# print("directory, error (swrad)")
168-
# for directory, error in errors:
169-
# print(directory, error) #
170-
#
171-
# monthly_errors = {str(mo): [] for mo in range(12)}
172-
#
173-
# for directory, error in errors:
174-
#
175-
# month, intcp, slope = (el.split(':')[1] for el in
176-
# directory.split(os.sep)[-1].split('_'))
177-
#
178-
# monthly_errors[month].append((intcp, slope, error))
179-
#
180-
# rankings = {
181-
# str(mo): list(sorted(monthly_errors[str(mo)], key=lambda x: x[-1]))
182-
# for mo in range(12)
183-
# }
184-
#
185-
# tops = [(mo, rankings[str(mo)][0]) for mo in range(12)]
186-
#
187-
# # update internal parameters
188-
# for top in tops:
189-
# mo = top[0]
190-
# self.parameters['dday_intcp'][mo] = tops[mo][1][0]
191-
# self.parameters['dday_slope'][mo] = tops[mo][1][1]
192-
#
193-
# return {
194-
# 'best': tops,
195-
# 'all': rankings
196-
# }
167+
def plot_srad_optimization(self, freq='daily', method='time_series'):
168+
"""
169+
Basic plotting of current srad optimization results with
170+
limited options for quick viewing, measured, original,
171+
and simulated swrad at the correspinding HRU is plotted
172+
either as time series (all three) or scatter (measured
173+
versus simulated). Not recommended for plotting results
174+
when n_sims is very high, use plotting options from
175+
an OptimizationResult object
176+
177+
Kwargs:
178+
freq (str): frequency of time series plots, value can be 'daily'
179+
or 'monthly' for solar radiation !!!need to finish monthly!!!
180+
method (str): 'time_series' for time series sub plot of each
181+
simulation alongside measured radiation. Other choice is
182+
'correlation' which plots each measured daily solar radiation
183+
value versus the corresponding simulated variable as subplots
184+
one for each simulation in the optimization. With coefficients
185+
of determiniationi i.e. square of pearson correlation coef.
186+
"""
187+
if not self.srad_outputs:
188+
raise ValueError('You have not run any srad optimizations')
189+
190+
# indices that measured and simulated swrad share (the intersection)
191+
X = self.measured_srad
192+
idx = X.index.intersection(self.srad_outputs[0]['statvar']\
193+
['swrad_{}'.format(self.srad_hru)].index)
194+
X = X[idx]
195+
n = len(self.srad_outputs) # number of simulations to plot
196+
197+
if freq == 'daily' and method == 'time_series':
198+
fig, ax = plt.subplots(n, sharex=True, sharey=True,\
199+
figsize=(12,n*3.5))
200+
axs = ax.ravel()
201+
for i,out in enumerate(self.srad_outputs):
202+
axs[i].plot(out['statvar']['swrad_{}'.format(self.srad_hru)]\
203+
[idx], 'r.', markersize=3, label='Simulated')
204+
axs[i].plot(self.measured_srad[idx], 'k.', markersize=3,\
205+
label='Measured')
206+
axs[i].set_ylabel('sim: {}'.format(out['simulation_dir'].\
207+
split(os.sep)[-1].replace('_', ' ')))
208+
if i == 0: axs[i].legend(markerscale=5, loc='best')
209+
fig.subplots_adjust(hspace=0)
210+
fig.autofmt_xdate()
211+
plt.show()
212+
213+
elif method == 'correlation':
214+
## number of subplots and rows (two plots per row)
215+
nrow = n//2 # round down if odd n
216+
ncol = 2
217+
odd_n = False
218+
if n/2. - nrow == 0.5:
219+
nrow+=1 # odd number need extra row
220+
odd_n = True
221+
## figure
222+
fig, ax = plt.subplots(nrows=nrow, ncols=ncol, figsize=(12,n*3))
223+
axs = ax.ravel()
224+
## subplot dimensions
225+
meas_min = min(X)
226+
meas_max = max(X)
227+
228+
for i, out in enumerate(self.srad_outputs):
229+
Y = out['statvar']['swrad_{}'.format(self.srad_hru)][idx]
230+
sim_max = max(Y)
231+
sim_min = min(Y)
232+
m = max(meas_max,sim_max)
233+
axs[i].plot([0, m], [0, m], 'k--', lw=2) ## one to one line
234+
axs[i].set_xlim(meas_min,meas_max)
235+
axs[i].set_ylim(sim_min, sim_max)
236+
axs[i].scatter(X, Y, facecolors='none', edgecolor='r', s=3)
237+
axs[i].set_ylabel('sim: {}'.format(out['simulation_dir']\
238+
.split(os.sep)[-1].replace('_', ' ')))
239+
axs[i].set_xlabel('Measured shortwave radiation')
240+
axs[i].text(0.05, 0.95,r'$R^2 = {0:.2f}$'.format(\
241+
X.corr(Y)**2), fontsize=16,\
242+
ha='left', va='center', transform=axs[i].transAxes)
243+
if odd_n: # empty subplot if odd number of simulations
244+
fig.delaxes(axs[n])
245+
246+
def _create_metafile_name(out_dir, opt_title, stage):
247+
"""
248+
Search through output directory where simulations are conducted
249+
look for all metadata simulation json files and find out if the
250+
current simulation is a replicate. Then use that information to
251+
build the correct file name for the output json file. The series
252+
are typically run in parallelel that is why this step has to be
253+
done after running multiple simulations from an optimization stage.
254+
255+
Args:
256+
out_dir (str): path to directory with model results, i.e.
257+
location where simulation series outputs and optimization
258+
json files are located, aka Optimizer.working_dir
259+
opt_title (str): optimization instance title for file search
260+
stage (str): stage of optimization, e.g. 'swrad', 'pet'
261+
262+
Returns:
263+
name (str): file name for the current optimization simulation series
264+
metadata json file. E.g 'dry_creek_swrad_opt.json', or if
265+
this is the second time you have run an optimization titled
266+
'dry_creek' the next json file will be returned as
267+
'dry_creek_swrad_opt1.json' and so on with integer increments
268+
"""
269+
swrad_meta_re = re.compile(r'^{}_{}_opt(\d*)\.json'.format(opt_title, stage))
270+
reps = []
271+
for f in os.listdir(out_dir):
272+
if swrad_meta_re.match(f):
273+
nrep = swrad_meta_re.match(f).group(1)
274+
if nrep == '':
275+
reps.append(0)
276+
else:
277+
reps.append(nrep)
278+
279+
if not reps:
280+
name = '{}_{}_opt.json'.format(opt_title, stage)
281+
else:
282+
# this is the nth optimization done under the same title
283+
n = max(map(int, reps)) + 1
284+
name = '{}_{}_opt{}.json'.format(opt_title, stage, n)
285+
return name
197286

198287
def _resample_param(param, p_min, p_max, noise_factor=0.1 ):
199288
"""
@@ -204,10 +293,11 @@ def _resample_param(param, p_min, p_max, noise_factor=0.1 ):
204293
to each parameter element by adding a RV from a normal distribution
205294
with mean 0, sigma = param allowable range / 10.
206295
207-
Arguments:
296+
Args:
208297
param (numpy.ndarray): ndarray of parameter to be resampled
209298
p_min (float): lower bound of PRMS allowable range for param
210299
p_max (float): upper bound of PRMS allowable range for param
300+
Kwargs:
211301
noise_factor (float): factor to multiply parameter range by,
212302
use the result as the standard deviation for the normal rand.
213303
variable used to add element wise noise. i.e. higher
@@ -220,7 +310,7 @@ def _resample_param(param, p_min, p_max, noise_factor=0.1 ):
220310

221311
low_bnd = p_min - np.min(param) # lowest param value minus allowable min
222312
up_bnd = p_max - np.max(param)
223-
s = (p_max - p_min) * noise_factor # default is one tenth rangee (0.1)
313+
s = (p_max - p_min) * noise_factor # stddev noise, default: range*(1/10)
224314

225315
shifted_param = np.random.uniform(low=low_bnd, high=up_bnd) + param
226316
## add noise to each point keeping result within allowable range
@@ -255,27 +345,30 @@ def get_optr_jsons(self, work_dir, stage):
255345
of corresponding json file paths for each stage as values.
256346
257347
Arguments:
258-
work_dir (str): path to simulation directory where
259-
optimization simulations were conducted and where
260-
corresponding json files should exist.
348+
work_dir (str): path to directory with model results, i.e.
349+
location where simulation series outputs and optimization
350+
json files are located, aka Optimizer.working_dir
261351
stage (str): the stage ('swrad', 'pet', 'flow', etc.) of
262352
the optimization in which to gather the jsons, if
263353
stage is 'all' then each stage will be gathered.
264354
Returns:
265355
ret (dict): dictionary of stage (keys) and lists of
266356
json file paths for that stage (values).
267357
"""
268-
358+
269359
ret = {}
270360
if stage != 'all':
361+
optr_metafile_re = re.compile(r'^.*_{}_opt(\d*)\.json'.format(stage))
271362
ret[stage] = [OPJ(work_dir, f) for f in\
272363
os.listdir(work_dir) if\
273-
f.endswith('_{0}_opt.json'.format(stage)) ]
364+
optr_metafile_re.match(f) ]
274365
else:
275366
stages = ['swrad', 'pet', 'flow']
276367
for s in stages:
277-
ret[s] = self.get_optr_jsons(work_dir, s)
278-
368+
optr_metafile_re = re.compile(r'^.*_{}_opt(\d*)\.json'.format(s))
369+
ret[s] = [OPJ(work_dir, f) for f in\
370+
os.listdir(work_dir) if\
371+
optr_metafile_re.match(f) ]
279372
return ret
280373

281374
class SradOptimizationResult(OptimizationResult):

0 commit comments

Comments
 (0)