Skip to content

Commit 3cf2ab9

Browse files
JohnVolkJohnVolk
committed
add example workflow for PAWN SA of PRMS parameters
1 parent 6f76b99 commit 3cf2ab9

1 file changed

Lines changed: 212 additions & 0 deletions

File tree

‎notebooks/PAWN_SA_workflow.ipynb‎

Lines changed: 212 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Example workflow of the PAWN global sensitivity analysis method using PRMS-Python on three arbitrary PRMS parameters \n",
8+
"\n",
9+
"The PAWN global, moment-independent, sensitivity analysis (SA) method is relatively straightforward to implement for any PRMS parameter using PRMS-Python objects. This example template can be modified to create the data necessary (emprical CDFs) for PAWN SA on any number of physical PRMS paraemters for any PRMS model. More information on the PAWN method can be found in the manuscript [here](http://dx.doi.org/10.1016/j.envsoft.2015.01.004). A case study using a slightly modified version of this script was used to conduct PAWN SA on 8 parameters that comprise the degree day solar radiation method in PRMS, details can be found in the PRMS-Python manuscript [here](). The control parameters and variable names for the experimental setup of PAWN as stated the literature are used within this code template, similarly this example code is heavily commented for clarity."
10+
]
11+
},
12+
{
13+
"cell_type": "code",
14+
"execution_count": 2,
15+
"metadata": {},
16+
"outputs": [],
17+
"source": [
18+
"import numpy as np\n",
19+
"import pandas as pd\n",
20+
"import os, json\n",
21+
"from prms_python import Data, Optimizer, Parameters, util\n",
22+
"from prms_python.optimizer import resample_param as resample\n",
23+
"from prms_python.optimizer import OptimizationResult\n",
24+
"import matplotlib.pyplot as plt\n",
25+
"import matplotlib\n",
26+
"matplotlib.style.use('ggplot')\n",
27+
"%matplotlib inline"
28+
]
29+
},
30+
{
31+
"cell_type": "markdown",
32+
"metadata": {},
33+
"source": [
34+
"## Define paths to initial model inputs, measured data, and initialize input objects"
35+
]
36+
},
37+
{
38+
"cell_type": "code",
39+
"execution_count": 3,
40+
"metadata": {},
41+
"outputs": [],
42+
"source": [
43+
"param_path 'path/to/model/parameters'\n",
44+
"data = Data('path/to/model/data') # this example does not modify climate forcing so a single instance is fine\n",
45+
"control_path = 'path/to/model/control'\n",
46+
"work_directory = '/path/to/output/directory'\n",
47+
"measrd_path = 'path/to/measured/output.csv'\n",
48+
"PRMS_output_variable = 'name_of_PRMS_output_variable_for_SA' # e.g. \"basin runoff 1\" or Y in PAWN"
49+
]
50+
},
51+
{
52+
"cell_type": "markdown",
53+
"metadata": {},
54+
"source": [
55+
"Note, to avoid platform specific path errors use os.path.join() as opposed to strings with platform path separators. For example if the path to the control file is $HOME/prms/control then use:\n",
56+
"\n",
57+
"```python\n",
58+
"home = os.path.expanduser(\"~\")\n",
59+
"control_path = os.path.join(home, 'prms', 'control')\n",
60+
"```\n"
61+
]
62+
},
63+
{
64+
"cell_type": "markdown",
65+
"metadata": {},
66+
"source": [
67+
"## Define sampling method, optimization title, and experimental setup control parameters"
68+
]
69+
},
70+
{
71+
"cell_type": "code",
72+
"execution_count": 11,
73+
"metadata": {},
74+
"outputs": [],
75+
"source": [
76+
"resample_meth = 'uniform' # sampling parameter values from a uniform distribution is one way for global SA\n",
77+
"title = 'ddsolrad_PAWN' # title optional\n",
78+
"nprocs = 8 # number of physical or logical processing cores to use\n",
79+
"\n",
80+
"# PAWN related experimental setup variables\n",
81+
"M = 3 # number of input factors for SA (parameters in this case)\n",
82+
"Nuc = 4000 # number of simulations to build the unconditional CDF\n",
83+
"Nc = 20 # number of times to resample each conditioning parameter\n",
84+
"n = 100 # number of simulations for buidling each conditional CDF\n",
85+
"\n",
86+
"# stages are used by the Optimizer object, \n",
87+
"# in PAWN these are essentially the names of each input factor 1,2,...,M\n",
88+
"stage_names = ['unconditional', 'cond_p1', 'cond_p2', 'cond_p3'] \n",
89+
"# each conditional stage is the exclusion of 1 of the 3 input parameters which is held 'constant' Nc times\n",
90+
"param_names_for_each_stage = [\n",
91+
" ['p1_name', 'p2_name', 'p3_name'], \n",
92+
" ['p2_name', 'p3_name'], \n",
93+
" ['p1_name', 'p3_name'], \n",
94+
" ['p1_name', 'p2_name']\n",
95+
" ]"
96+
]
97+
},
98+
{
99+
"cell_type": "markdown",
100+
"metadata": {},
101+
"source": [
102+
"## Conduct Nuc simulations on all M parameters to build the unconditional CDF"
103+
]
104+
},
105+
{
106+
"cell_type": "code",
107+
"execution_count": null,
108+
"metadata": {},
109+
"outputs": [],
110+
"source": [
111+
"stage_name = stage_names[0] # \"unconditional\"\n",
112+
"archive_dir = os.path.join(work_directory,\"{}_archived\".format(stage_name))\n",
113+
"# make archive directory if it doesn't exist\n",
114+
"if not os.path.isdir(archive_dir): \n",
115+
" os.mkdir(archive_dir)\n",
116+
"# create an Optimizer instance and call monte_carlo method\n",
117+
"optr = Optimizer(Parameters(param_path), data, control_path, work_directory, title=title)\n",
118+
"optr.monte_carlo(\n",
119+
" measured_path, \n",
120+
" param_names_for_each_stage[0], \n",
121+
" PRMS_output_variable, \n",
122+
" method=resample_meth,\n",
123+
" n_sims=Nuc, \n",
124+
" nproc=nprocs, \n",
125+
" stage=stage_name\n",
126+
" )\n",
127+
"# optionally archive output to reduce disk space\n",
128+
"result = OptimizationResult(work_directory, stage=stage_name)\n",
129+
"result.archive()"
130+
]
131+
},
132+
{
133+
"cell_type": "markdown",
134+
"metadata": {},
135+
"source": [
136+
"## Create conditional CDFs Nc times for each parameter, storing conditioning parameter values"
137+
]
138+
},
139+
{
140+
"cell_type": "code",
141+
"execution_count": null,
142+
"metadata": {},
143+
"outputs": [],
144+
"source": [
145+
"for nc in range(Nc): # number of conditional param bootstrap resamples- distributable to multiple machines/nodes\n",
146+
" for i, stage in enumerate(stage_names): # xi param of M total params\n",
147+
" if stage=='unconditional': # already built unconditional CDF- skip\n",
148+
" continue\n",
149+
" # use unconditional param list set minus conditioning param list to find conditioning param name\n",
150+
" conditional_param = list(set(stage_names['unconditional']) - set(stage_names[stage]))[0]\n",
151+
" # make a Parameter instance to resample the conditioning param and others \n",
152+
" params = Parameters(param_path)\n",
153+
" xi_values = resample(params, conditional_param)\n",
154+
" params[conditional_param] = xi_values # assign resampled conditioning parameter (xi) to Parameter object\n",
155+
" stage_name = stage+str(nc) # add nc for nc'th bootstrap resampling round\n",
156+
" # create output archive directory to hold info on xi conditioning values and correspinding model output\n",
157+
" archive_dir = os.path.join(work_directory,\"{}_archived\".format(stage_name))\n",
158+
" if not os.path.isdir(archive_dir):\n",
159+
" os.mkdir(archive_dir)\n",
160+
" if xi_values.shape == (): # some parameters may be single valued, this worked for me\n",
161+
" with open(os.path.join(archive_dir,'{}_Nc_{}.txt'.format(conditional_param,i)), 'w') as outf:\n",
162+
" outf.write(str(xi_values))\n",
163+
" else: # ndarrays can be dumped to a text file using numpy\n",
164+
" np.savetxt(os.path.join(archive_dir,'{}_Nc_{}.txt'.format(conditional_param,i)), xi_values, fmt=\"%f\")\n",
165+
" # create an Optimizer instance and call monte_carlo method with the correct parameter modifications\n",
166+
" optr = Optimizer(params, data, control_path, work_directory, title=title)\n",
167+
" optr.monte_carlo(\n",
168+
" measured_path, \n",
169+
" param_names_for_each_stage[i], \n",
170+
" PRMS_output_variable, \n",
171+
" method=resample_meth,\n",
172+
" n_sims=n, \n",
173+
" nproc=nprocs, \n",
174+
" stage=stage_name\n",
175+
" )\n",
176+
" # optionally archive output to reduce disk space\n",
177+
" result = OptimizationResult(work_directory, stage=stage_name)\n",
178+
" result.archive()"
179+
]
180+
},
181+
{
182+
"cell_type": "markdown",
183+
"metadata": {},
184+
"source": [
185+
"## That's it, now analyze output to calculate sensitivity indices for each input parameter\n",
186+
"\n",
187+
"Although this example does not include analysis of results, it is straightforward to build output CDFs from the archived JSON files. For example of accessing output from these files please refer to the Jupyter Notebooks that packs with PRMS-Python for the `OptimizationResult` object. For convenience we added Python functions that calculate emprical CDFs and the Kolmogorov-Smirnov distance between two CDFs in the `prms_python.util` module. These functions can be used to easily calculate the PAWN sensitivity analysis from the results produced using the template above. "
188+
]
189+
}
190+
],
191+
"metadata": {
192+
"kernelspec": {
193+
"display_name": "Python 3",
194+
"language": "python",
195+
"name": "python3"
196+
},
197+
"language_info": {
198+
"codemirror_mode": {
199+
"name": "ipython",
200+
"version": 3
201+
},
202+
"file_extension": ".py",
203+
"mimetype": "text/x-python",
204+
"name": "python",
205+
"nbconvert_exporter": "python",
206+
"pygments_lexer": "ipython3",
207+
"version": "3.5.2"
208+
}
209+
},
210+
"nbformat": 4,
211+
"nbformat_minor": 2
212+
}

0 commit comments

Comments
 (0)