Skip to content

Commit 699fa9d

Browse files
mt-digitalmt-digital
committed
simulation series test pass
1 parent b3c0f49 commit 699fa9d

5 files changed

Lines changed: 270 additions & 77 deletions

File tree

‎prms_python/__init__.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
from .data import Data
22
from .parameters import Parameters, modify_params
3-
from .simulation import Simulation
3+
from .simulation import Simulation, SimulationSeries
44
from .scenario import Scenario, ScenarioSeries
55
from .util import load_statvar, load_data_file, nash_sutcliffe

‎prms_python/optimizer.py‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
'''
2+
optimizer.py -- Optimization routines for PRMS parameters and data.
3+
'''
4+
import pandas as pd
5+
import numpy as np
6+
import os
7+
8+
from .data import Data
9+
from .parameters import Parameters
10+
from .scenario import ScenarioSeries
11+
12+
13+
class Optimizer:
14+
'''
15+
Container for a PRMS parameter optimization routine consisting of the
16+
four stages as described in Hay, et al, 2006
17+
(ftp://brrftp.cr.usgs.gov/pub/mows/software/luca_s/jawraHay.pdf).
18+
19+
Example:
20+
21+
>>> from prms_python import Data, Optimizer, Parameters
22+
>>> params = Parameters('path/to/parameters')
23+
>>> data = Data('path/to/data')
24+
>>> optr = Optimizer(params, data, title='the title', description='desc')
25+
>>> optr.srad('path/to/reference_data/measured_srad.csv')
26+
27+
'''
28+
29+
def __init__(self, parameters, data, working_dir,
30+
title=None, description=None):
31+
32+
if isinstance(parameters, Parameters):
33+
self.parameters = parameters
34+
else:
35+
raise TypeError('parameters must be instance of Parameters')
36+
37+
if isinstance(data, Data):
38+
self.data = data
39+
else:
40+
raise TypeError('data must be instance of Data')
41+
42+
self.working_dir = working_dir
43+
self.title = title
44+
self.description = description
45+
46+
def srad(self, reference_srad_path, station_nhru, method='',
47+
dday_intcp_range=None, dday_slope_range=None,
48+
intcp_delta=None, slope_delta=None):
49+
'''
50+
Optimize the monthly dday_intcp and dday_slope parameters by one of
51+
two methods: 'uniform' or 'random' for uniform sampling
52+
53+
Args:
54+
reference_srad_path (str): path to measured solar radiation data
55+
Kwargs:
56+
method (str): 'uniform' or 'random'; if 'random',
57+
intcp_delta and slope_delta are ignored, if provided
58+
dday_intcp_range ((float, float)): two-tuple of minimum and
59+
maximum value to consider for the dday_intcp parameter
60+
dday_slope_range ((float, float)): two-tuple of minimum and
61+
maximum value to consider for the dday_slope parameter
62+
intcp_delta (float): resolution of grid to test in intcp dimension
63+
slope_delta (float): resolution of grid to test in slope dimension
64+
65+
Returns:
66+
(SradOptimizationResult)
67+
'''
68+
if dday_intcp_range is None:
69+
dday_intcp_range = (-60.0, 10.0)
70+
intcp_delta = 10.0
71+
elif intcp_delta is None:
72+
intcp_delta = (dday_intcp_range[1] - dday_intcp_range[0]) / 4.0
73+
74+
if dday_slope_range is None:
75+
dday_slope_range = (0.2, 0.9)
76+
slope_delta = .05
77+
elif slope_delta is None:
78+
slope_delta = (dday_slope_range[1] - dday_slope_range[0]) / 4.0
79+
80+
# create parameters
81+
ir = dday_intcp_range
82+
sr = dday_slope_range
83+
84+
intcps = np.arange(ir[0], ir[1], intcp_delta)
85+
slopes = np.arange(sr[0], sr[1], slope_delta)
86+
87+
param_grid = np.meshgrid(intcps, slopes)
88+
89+
def _mod_params(parameters, month, intcp, slope):
90+
91+
parameters['dday_intcp'][month] = intcp
92+
parameters['dday_slope'][month] = slope
93+
94+
parameters_iter = (
95+
{
96+
'parameters':
97+
_mod_params(self.parameters, month, intcp, slope),
98+
99+
'title': '"month":{0},"dday_intcp":{1:.3f},'
100+
'"dday_slope":{2:.3f}'.format(month, intcp, slope),
101+
}
102+
for month in range(12)
103+
for intcp, slope in param_grid
104+
)
105+
106+
# create ScenarioSeries from parameters
107+
108+
# XXX TODO XXX TODO
109+
series = ScenarioSeries.from_params_iter(
110+
self.working_dir,
111+
parameters_iter=parameters_iter,
112+
title=self.title,
113+
description=self.description
114+
)
115+
116+
# run all scenarios
117+
series.run()
118+
119+
def _error(x, y):
120+
return float(abs(x - y))/float(len(x))
121+
122+
# calculate the top performing
123+
modeled_srads = (
124+
(output['title'], output['statvar']['swrad_' + str(station_nhru)])
125+
126+
# XXX TODO XXX TODO
127+
for output in series.outputs
128+
)
129+
130+
measured_srad = pd.read_csv(reference_srad_path, parse_dates=True)
131+
132+
errors = (
133+
(modeled_srads[0], _error(measured_srad, modeled_srads[1]))
134+
for modeled_srad in modeled_srads
135+
)
136+
rankings = list(sorted(errors, key=lambda x: x[1]))
137+
138+
# update internal parameters
139+
self.parameters = series.outputs[rankings[0]]['parameters']
140+
141+
return rankings
142+
143+
144+
class OptimizationResult:
145+
146+
pass
147+
148+
149+
class SradOptimizationResult(OptimizationResult):
150+
151+
pass

‎prms_python/scenario.py‎

Lines changed: 0 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -193,11 +193,6 @@ def run(self, prms_exec='prms', nproc=None):
193193
pool = mp.Pool(processes=nproc)
194194
pool.map(_scenario_runner, self.scenarios)
195195

196-
# self.outputs = [
197-
# ScenarioOutput(uu, os.path.join(os.curdir(), d))
198-
# for uu, d in self.metadata['uuid_title_map'].items()
199-
# ]
200-
201196

202197
# multiprocessing req the function be def'd at root scope so it's picklable
203198
def _scenario_runner(scenario, prms_exec='prms'):
@@ -296,16 +291,3 @@ def __setitem__(self, key, value):
296291
def write(self, output_path):
297292
with open(output_path, 'w') as f:
298293
f.write(json.dumps(self.metadata_dict))
299-
300-
301-
class ScenarioOutput:
302-
303-
def __init__(self, scenario_uu, scenario_directory, title=None):
304-
opj = os.path.join
305-
self.uuid = scenario_uu
306-
self.scenario_directory = scenario_directory
307-
self.title = title
308-
self.data = Data(opj(scenario_directory, 'data'))
309-
self.parameters = Parameters(opj(scenario_directory, 'parameters'))
310-
self.statvar = load_statvar(opj(scenario_directory, 'statvar.dat'))
311-
self.control = open(opj(scenario_directory, 'control')).read()

‎prms_python/simulation.py‎

Lines changed: 50 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,10 @@
66

77
from .data import Data
88
from .parameters import Parameters
9+
from .util import load_statvar
10+
11+
12+
OPJ = os.path.join
913

1014

1115
class SimulationSeries(object):
@@ -14,7 +18,9 @@ class SimulationSeries(object):
1418
'''
1519

1620
def __init__(self, simulations):
17-
self.series = simulations
21+
# XXX TODO would love to not have to use list here, but otherwise can't
22+
# access the simulations after they have run through map
23+
self.series = list(simulations)
1824

1925
def run(self, prms_exec='prms', nproc=None):
2026

@@ -24,6 +30,45 @@ def run(self, prms_exec='prms', nproc=None):
2430
pool = mp.Pool(processes=nproc)
2531
pool.map(_simulation_runner, self.series)
2632

33+
return self
34+
35+
def outputs_iter(self):
36+
'''
37+
Return an iterator of directories with the path to the simulation_dir
38+
as well as a pandas.DataFrame of the statvar output, and the Data and
39+
Parameters representations used in the simulation.
40+
41+
Example:
42+
>>> ser = SimulationSeries(simulations)
43+
>>> ser.run()
44+
>>> g = ser.outputs_iter()
45+
>>> print(g.next())
46+
47+
Would return something like
48+
49+
{'simulation_dir': 'path/to/sim/', 'statvar': <pandas.DataFrame>,
50+
'data': <data.Data>, 'parameters': <parameters.Parameters>}
51+
52+
Returns:
53+
(generator(dict)):
54+
'''
55+
dirs = (s.simulation_dir for s in self.series)
56+
57+
print(self.series)
58+
59+
return (
60+
{
61+
'simulation_dir': d,
62+
'statvar': load_statvar(OPJ(d, 'outputs', 'statvar.dat')),
63+
'data': Data(OPJ(d, 'inputs', 'data')),
64+
'parameters': Parameters(OPJ(d, 'inputs', 'parameters'))
65+
}
66+
for d in dirs
67+
)
68+
69+
def __len__(self):
70+
return len(list(self.outputs_iter()))
71+
2772

2873
def _simulation_runner(sim):
2974
sim.run(prms_exec='prms')
@@ -118,7 +163,7 @@ def from_data(cls, data, parameters, control_path, simulation_dir):
118163
raise TypeError('data must be instance of Data')
119164

120165
if not isinstance(parameters, Parameters):
121-
raise TypeError('parameters must be instance of Parameters')
166+
raise TypeError('parameters must be instance of Parameters, not ' + str(type(parameters)))
122167

123168
if os.path.exists(simulation_dir):
124169
shutil.rmtree(simulation_dir)
@@ -129,17 +174,15 @@ def from_data(cls, data, parameters, control_path, simulation_dir):
129174
sim.simulation_dir = simulation_dir
130175

131176
sd = simulation_dir
132-
opj = os.path.join
133177

134-
data_path = opj(sd, 'data')
178+
data_path = OPJ(sd, 'data')
135179
data.write(data_path)
136-
params_path = opj(sd, 'parameters')
180+
params_path = OPJ(sd, 'parameters')
137181
parameters.write(params_path)
138-
shutil.copy(control_path, opj(sd, 'control'))
182+
shutil.copy(control_path, OPJ(sd, 'control'))
139183

140184
return sim
141185

142-
143186
def run(self, prms_exec='prms'):
144187

145188
cwd = os.getcwd()

0 commit comments

Comments
 (0)