In [1]:
import ipyparallel
client = ipyparallel.Client()
In [2]:
client.ids
Out[2]:
In [3]:
%%px --local
from __future__ import division # python2
import math
import numpy as np
from scipy.optimize import brentq
def get_antropogenic_release(xt, c1, c2, r1, r2, w1):
'''
Parameters
----------
xt : float
polution in lake at time t
c1 : float
center rbf 1
c2 : float
center rbf 2
r1 : float
ratius rbf 1
r2 : float
ratius rbf 2
w1 : float
weight of rbf 1
note:: w2 = 1 - w1
'''
rule = w1*(abs(xt-c1/r1))**3+(1-w1)*(abs(xt-c2/r2))**3
at = min(max(rule, 0.01), 0.1)
return at
def lake_model(b=0.42, q=2.0, mean=0.02, stdev=0.001, alpha=0.4,
delta=0.98, c1=0.25, c2=0.25, r1=0.5, r2=0.5,
w1=0.5, nsamples=100, steps=100, seed=None):
'''runs the lake model for 1 stochastic realisation using specified
random seed.
Parameters
----------
b : float
decay rate for P in lake (0.42 = irreversible)
q : float
recycling exponent
mean : float
mean of natural inflows
stdev : float
standard deviation of natural inflows
alpha : float
utility from pollution
delta : float
future utility discount rate
c1 : float
c2 : float
r1 : float
r2 : float
w1 : float
steps : int
the number of time steps (e.g., days)
seed : int, optional
seed for the random number generator
'''
np.random.seed(seed)
Pcrit = brentq(lambda x: x**q/(1+x**q) - b*x, 0.01, 1.5)
X = np.zeros((steps,))
decisions = np.zeros((steps,))
X[0] = 0.0
natural_inflows = np.random.lognormal(
math.log(mean**2 / math.sqrt(stdev**2 + mean**2)),
math.sqrt(math.log(1.0 + stdev**2 / mean**2)),
size=steps)
for t in range(steps-1):
decisions[t] = get_antropogenic_release(X[t], c1, c2, r1, r2, w1)
X[t+1] = (1-b)*X[t] + X[t]**q/(1+X[t]**q) + decisions[t] + natural_inflows[t]
reliability = np.sum(X < Pcrit)/steps
utility = np.sum(alpha*decisions*np.power(delta,np.arange(steps)))
# note that I have slightly changed this formulation to retain
# consistency with the equations in the papers
inertia = np.sum(np.abs(np.diff(decisions)) < 0.01)/(steps-1)
return X, utility, inertia, reliability
In this blog, I will continue to showcase the functionality of the exploratory modelling. In the previous blog, I have given a general introduction to the workbench, and showed how the DPS example that comes with Rhodium can be adapted for use with the workbench. In this blogpost, I will showcase how the workbench can be used for open exploration.
In exploratory modeling, we are interested in understanding how regions in the uncertainty space and/or the decision lever space map to the outcome space, or partitions thereof. There are two general approaches for investigating this mapping. The first one is through systematic sampling of the uncertainty or lever space. This is sometimes also known as open exploration. The second one is to search through the space in a directed manner using some type of optimization approach. This is sometimes also known as directed search.
The workbench support both open exploration and directed search. Both can be applied to investigate the mapping of the uncertainty space or the decision lever space to the outcome space. In most applications, search is used for finding promising mappings from the lever space to the outcome space, while exploration is used to stress test these mappings under a whole range of possible resolutions to the various uncertainties. This need not be the case however. Optimization can be used to discover the worst possible scenario, while sampling can be used to get insight into the sensitivity of outcomes to the various decision levers.
To showcase the functionality, let's start with a simple example. We are going to simultaneously sample over uncertainties and levers. We are going to generate 100 scenarios, 10 policies, and see how they jointly affect the outcomes. By default, the workbench will use Latin Hypercube sampling for generating both the scenarios and the policies. Each policy will be always evaluated over all scenarios (i.e. a full factorial over scenarios and policies).
In [4]:
%%px --local
def process_p(values):
values = np.asarray(values)
values = np.mean(values, axis=0)
return np.max(values)
In [5]:
from ema_workbench import (RealParameter, ScalarOutcome, Constant,
ReplicatorModel)
model = ReplicatorModel('lakeproblem', function=lake_model)
model.replications = 150
#specify uncertainties
model.uncertainties = [RealParameter('b', 0.1, 0.45),
RealParameter('q', 2.0, 4.5),
RealParameter('mean', 0.01, 0.05),
RealParameter('stdev', 0.001, 0.005),
RealParameter('delta', 0.93, 0.99)]
# set levers
model.levers = [RealParameter("c1", -2, 2),
RealParameter("c2", -2, 2),
RealParameter("r1", 0, 2),
RealParameter("r2", 0, 2),
RealParameter("w1", 0, 1)]
#specify outcomes
model.outcomes = [ScalarOutcome('max_P', kind=ScalarOutcome.MINIMIZE,
function=process_p),
ScalarOutcome('utility', kind=ScalarOutcome.MAXIMIZE,
function=np.mean),
ScalarOutcome('inertia', kind=ScalarOutcome.MINIMIZE,
function=np.mean),
ScalarOutcome('reliability', kind=ScalarOutcome.MAXIMIZE,
function=np.mean)]
# override some of the defaults of the model
model.constants = [Constant('alpha', 0.41),
Constant('steps', 100)]
In [6]:
from ema_workbench import (IpyparallelEvaluator, ema_logging,
perform_experiments)
ema_logging.log_to_stderr(ema_logging.INFO)
with IpyparallelEvaluator(model, client) as evaluator:
results = evaluator.perform_experiments(scenarios=100, policies=10)
In [ ]:
experiments, outcomes = results
In [ ]:
from ema_workbench.analysis import pairs_plotting
fig, axes = pairs_plotting.pairs_scatter(results, group_by='policy',
legend=False)
fig.set_size_inches(8,8)
plt.show()
In [ ]:
In [ ]:
from ema_workbench.analysis import prim
x = experiments
y = outcomes['max_P'] <0.8
prim_alg = prim.Prim(x, y, threshold=0.8)
box1 = prim_alg.find_box()
In [ ]:
import mpld3
box1.show_tradeoff()
plt.show()
In [ ]:
experiments, outcomes = results
for key, value in outcomes.items():
params = model.uncertainties #+ model.levers[:]
fig, axes = plt.subplots(ncols=len(params), sharey=True)
y = value
for i, param in enumerate(params):
ax = axes[i]
ax.set_xlabel(param.name)
pearson = sp.stats.pearsonr(x, y)
ax.annotate("r: {:6.3f}".format(pearson[0]), xy=(0.15, 0.85),
xycoords='axes fraction',fontsize=13)
x = experiments[param.name]
sns.regplot(x, y, ax=ax, ci=None, color='k',
scatter_kws={'alpha':0.2, 's':8, 'color':'gray'})
ax.set_xlim(param.lower_bound, param.upper_bound)
axes[0].set_ylabel(key)
plt.show()
In [ ]:
from ema_workbench.analysis import feature_scoring
x = experiments
y = outcomes
fs = feature_scoring.get_feature_scores_all(x, y)
In [ ]:
sns.heatmap(fs, cmap='viridis', annot=True)
plt.show()
In [ ]:
with MultiprocessingEvaluator(model) as evaluator:
sa_results = evaluator.perform_experiments(scenarios=1000,
uncertainty_sampling='sobol')
In [ ]:
from SALib.analyze import sobol
from ema_workbench.em_framework.salib_samplers import get_SALib_problem
experiments, outcomes = sa_results
problem = get_SALib_problem(model.uncertainties)
Si = sobol.analyze(problem, outcomes['max_P'],
calc_second_order=True, print_to_console=False)
In [ ]:
Si_filter = {k:Si[k] for k in ['ST','ST_conf','S1','S1_conf']}
Si_df = pd.DataFrame(Si_filter, index=problem['names'])
Si_df
In [ ]: