<jupyter_start><jupyter_text>Adaptive - IjsselmeerAuthor: <NAME>, Delft University of Technology, Delft, The Netherlands. Section of Policy Analysis  email: <EMAIL><jupyter_code>cd ..
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
import seaborn as sns
import ast
from pathways.pathways import get_policy_def, find_pathway
from ema_workbench import load_results
import collections

plt.style.use('default')
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"<jupyter_output><empty_output><jupyter_text>Load data<jupyter_code># load data from EMA
results_df = pd.read_csv('data/Ijsselmeer/2018_11_05_results.csv',index_col=0)
len(results_df)
results_df.head()
results_df[results_df['low_flow_shift']==0][results_df['peak_inflow_multiplier']==1][results_df['sea_level_shift']==0]<jupyter_output>/Users/lraso/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:1: UserWarning: Boolean Series key will be reindexed to match DataFrame index.
  if __name__ == '__main__':<jupyter_text>Specify policy dimensions<jupyter_code>actions_name = ['pump_capacity','sluices_widening','summer_target']
critical_uncertainties = ['low_flow_shift', 'peak_inflow_multiplier', 'sea_level_shift']
objectives = ['dike failure frequency','demand deficit relative']
color_list = ['clMaroon','clNavy','clOlive','clGreen','clPurple','clOrange','clYellow']
actions = []
for i, action_name in enumerate(actions_name):
    actions.append({'name':action_name, 
                    'caption': action_name, 
                    'color': color_list[i],
                    'predecessor':'current_situation'}) 

decreasing_actions = []#['winter_target'] # list actions whose effectiveness increase with decreasing of action value, default increasing
for action in actions:
    direction = -1 if action in decreasing_actions else 1
    action['values'] = np.sort(results_df[action['name']].unique())[::direction].tolist()<jupyter_output><empty_output><jupyter_text>2-action Pathways  We consider only pathways made of two actions (because the Deltares Pathway Generator does not include pathways made of more than two actions)<jupyter_code>pathways = []
for i, action in enumerate(actions):
    left_actions = actions[i+1:]
    for combined_action in left_actions:
        pathway = {'name': action['name'],
                   'caption': action['name'] + ' + ' + combined_action['name'],
                   'combined_to': combined_action['name'],
                   'predecessor': action['name'],
                   'color': action['color'], 
                   'values': [action['values'][-1]]}
        pathways.append(pathway)
# var new_var correspondence
name_crit_unc ={'low_flow_shift':'low flow reduction $[m^3/s]$',
      'peak_inflow_multiplier': 'max discharge increase $[\%]$ ',
      'sea_level_shift': 'sea level rise $[m]$',
      'wind_multiplier' : 'change wind speed, $[\%]$ '}
     # 'winter_target': 'winter target level $[mNAP]$',}
    
    
# rename
old_var_name = results_df.columns.copy(deep=True)
new_var_name = [ name_crit_unc[var] if var in name_crit_unc else var for var in old_var_name  ]
# results_df.columns = new_var_name<jupyter_output><empty_output><jupyter_text>define scenarios data<jupyter_code>scenarios_df = results_df.iloc[results_df['scenario_id'].unique(),0:len(critical_uncertainties)]
scenarios_df.shape<jupyter_output><empty_output><jupyter_text>Define policy success<jupyter_code>policy_success_threshold = {}
policy_success_threshold['dike failure frequency'] = 2500
policy_success_threshold['demand deficit relative'] = 0.015 # in *100 %

sign = {} # negative is the objective to be minimized, positive if to be maximized 
sign['dike failure frequency'] = -1
sign['demand deficit relative'] = +1<jupyter_output><empty_output><jupyter_text>Identify policy success scenarios<jupyter_code>policy_success = pd.DataFrame(results_df[['policy','scenario_id']])
for objective in objectives:
    policy_success[objective] = sign[objective] * results_df[objective] <= sign [objective] * policy_success_threshold[objective]<jupyter_output><empty_output><jupyter_text>overall success of policy success is if all criteria are satisfied<jupyter_code>policy_success['overall']= policy_success[objectives].all(axis=1)<jupyter_output><empty_output><jupyter_text>Pivot policy-success results in scenarios (rows) vs policy (columns) table<jupyter_code>success_per_policy = policy_success.pivot(columns='policy', values='overall', index='scenario_id')   
policies = success_per_policy.columns
success_per_policy.sum()<jupyter_output><empty_output><jupyter_text>Include observational uncertainty  Observational uncertainty Signposts<jupyter_code>forcings = pd.read_csv('data/Ijsselmeer/forcings.csv',index_col=0,parse_dates=[0])
window_size = 10 # window size of ten years<jupyter_output><empty_output><jupyter_text>Signpost $S_{sea}$<jupyter_code>yearly_wl_mean = forcings['sea level'].groupby(forcings.index.year).mean()
yearly_wl_mean.plot(marker='.')
sigma_Ssea = yearly_wl_mean.std() / np.sqrt( window_size )
sigma_Ssea<jupyter_output><empty_output><jupyter_text>Signpost $S_{qmax}$<jupyter_code>yearly_q_max = forcings['inflow ijssel'].groupby(forcings.index.year).max()
yearly_q_max.plot(marker='.',ylim=[0,2500])
mu_Sqmax = yearly_q_max.mean()
sigma_Sqmax = yearly_q_max.std() / np.sqrt( window_size )
mu_Sqmax, sigma_Sqmax<jupyter_output><empty_output><jupyter_text>Signpost $S_{qmin}$<jupyter_code>yearly_q_min = forcings['inflow ijssel'].groupby(forcings.index.year).min()
yearly_q_min.plot(marker='.')
mu_Sqmin = yearly_q_min.mean()
sigma_Sqmin = yearly_q_min.std() / np.sqrt( window_size )
mu_Sqmin, sigma_Sqmin<jupyter_output><empty_output><jupyter_text>Include Observational Uncertainty<jupyter_code>N_sampling = 30
obs_unc = {'S_sea':{'uncertainty':'sea_level_shift','sigma':sigma_Ssea},
           'S_qmax':{'uncertainty':'peak_inflow_multiplier','mult_factor':mu_Sqmax,'sigma':sigma_Sqmax},
           'S_qmin':{'uncertainty':'low_flow_shift','mult_factor':-1,'mu':mu_Sqmin,'sigma':sigma_Sqmin}}
scenarios_indicators = pd.concat([scenarios_df]*(N_sampling*len(obs_unc)))
success_per_policy_indicators = pd.concat([success_per_policy]*(N_sampling*len(obs_unc)))
size_df = len(scenarios_indicators)
for signpost in obs_unc:
    mult_factor = obs_unc[signpost]['mult_factor'] if 'mult_factor' in obs_unc[signpost] else 1
    mu = obs_unc[signpost]['mu'] if 'mu' in obs_unc[signpost] else 0
    scenarios_indicators[signpost] = mult_factor * scenarios_indicators[obs_unc[signpost]['uncertainty']] + np.random.normal(mu,obs_unc[signpost]['sigma'],size_df)
scenarios_indicators.plot(kind='scatter',x='peak_inflow_multiplier',y='S_qmax')
scenarios_indicators.plot(kind='scatter',x='sea_level_shift',y='S_sea')
signposts = [signpost for signpost in obs_unc]
signposts<jupyter_output><empty_output><jupyter_text>Identify Adaptation Tipping Points by scenario discovery  Find ATP on the space of critical uncertainties  We consider two criteria* Parsimony* Fixed level of riskThe parsimonious approach identify a monitoring system made of one signpost for each action, selecting the signpost that better identify policy success, and select the trigger value that maximize density (min risk). The fixed level of risk approach selects one or more the signposts that offer a level of density not smaller than 90%, equivalent to a risk of missing the adaptation signal with 10% of chance.       Use of Range Search (PRIM)<jupyter_code>from ema_workbench.analysis import prim
signpost_based_monitor = False
observables = signposts if signpost_based_monitor else critical_uncertainties
boxes = {}

cols = ['UpperBound','LowerBound','RelativeATP','Density','Coverage']
ATP_parsimony = pd.DataFrame(index=policies, columns=cols)
ATP_risk_averse = pd.DataFrame(index=policies, columns=cols)

lower_lim = 0
upper_lim = 1 # select the second value (i.e. the max value) in the box, because the signposts are all activated by larger values

if signpost_based_monitor:
    lower_bound = np.array([58.0,817.8327073939638, -0.049]) #mu_Sqmin - 100 - 3 * sigma_Sqmin, mu_Sqmax - 3 * sigma_Sqmax ,0.0 - 3 * sigma_Ssea])
    upper_bound = np.array([254.70, 3208., 1.05])  #[mu_Sqmin - 3 * sigma_Sqmin, 2 * mu_Sqmax + 3 * sigma_Sqmax,1.0 + 3* sigma_Ssea])
else:
    lower_bound = np.array([0.,1.,0.])
    upper_bound = np.array([100.,2.,1.])

 

with open('ATP.txt','w') as file:
    for policy in policies:
        x = scenarios_indicators[observables].to_records(index=False)
        y = success_per_policy_indicators[policy].values
        prim_alg = prim.Prim(x, y, threshold=0.05,peel_alpha=0.025,paste_alpha=0.025,mass_min=0.025)
        box = prim_alg.find_box()
        box.peeling_trajectory['risk missed adaptation'] = (1 - box.peeling_trajectory['density']) * (box.peeling_trajectory['mass'] / box.peeling_trajectory['mean'])
        box.peeling_trajectory['risk regretful action'] = 1 - box.peeling_trajectory['coverage']
        boxes[policy] = box
      
        file.write(policy +'\n')
        file.write(box.peeling_trajectory.to_string()+'\n')    
# criteria 
citerium = {}
citerium['parsimony'] = 
criterium['risk'] = 
         
for policy in policies:        # criteria : single sigpost with the largest density
        i = max(box.peeling_trajectory.index[box.peeling_trajectory['res dim'] == 1])
        upper_bound_box_rec = box.box_lims[i][:][upper_lim]
        upper_bound_box = np.array([val.astype(float) for val in upper_bound_box_rec])
        ATP_parsimony['UpperBound'].loc[policy] = upper_bound_box
        
        lower_bound_box_rec = box.box_lims[i][:][lower_lim]
        lower_bound_box = np.array([val.astype(float) for val in lower_bound_box_rec])
        ATP_parsimony['LowerBound'].loc[policy] = lower_bound_box
        
        print( (upper_bound_box - lower_bound_box) / (upper_bound - lower_bound)) 
        
        ATP_parsimony['RelativeATP'].loc[policy] = np.mean( (upper_bound_box - lower_bound_box ) / (upper_bound - lower_bound) ) #/  (upper_bound - lower_bound) )
        ATP_parsimony['Density'].loc[policy] = box.peeling_trajectory.density[i]
        ATP_parsimony['Coverage'].loc[policy] = box.peeling_trajectory.coverage[i]

        
        # criteria : density not smaller than threshold %
        if signpost_based_monitor:
            th = density['uncertainties'][policy]
        else:
            th = 0.75
        
        i = min(box.peeling_trajectory.index[box.peeling_trajectory['density'] > th])
        upper_bound_box_rec = box.box_lims[i][:][upper_lim]
        upper_bound_box = np.array([val.astype(float) for val in upper_bound_box_rec])
        ATP_risk_averse['UpperBound'].loc[policy] = upper_bound_box
        
        lower_bound_box_rec = box.box_lims[i][:][lower_lim]
        lower_bound_box = np.array([val.astype(float) for val in lower_bound_box_rec])
        ATP_risk_averse['LowerBound'].loc[policy] = lower_bound_box

        ATP_risk_averse['RelativeATP'].loc[policy] = np.mean( (upper_bound_box - lower_bound_box) /  (upper_bound - lower_bound) )
        ATP_risk_averse['Density'].loc[policy] = box.peeling_trajectory.density[i]
        ATP_risk_averse['Coverage'].loc[policy] = box.peeling_trajectory.coverage[i]
lower_bound,upper_bound
ATP_parsimony
ATP_risk_averse<jupyter_output><empty_output><jupyter_text>ZOOM IN: tradeoff density-coverage<jupyter_code>coverage = {}
density = {}
if signpost_based_monitor:
    coverage['signposts'] = ATP_risk_averse['Coverage']
    density['signposts'] = ATP_risk_averse['Density']
else: 
    coverage['uncertainties'] = ATP_risk_averse['Coverage']
    density['uncertainties'] = ATP_risk_averse['Density']
density['uncertainties'][policy]
coverage['diff'] = coverage['uncertainties'] - coverage['signposts']
density['diff'] = density['uncertainties'] - coverage['signposts']
coverage['diff'], density['diff']
b = boxes["{'sluices_widening': 1, 'summer_target': -0.2, 'pump_capacity': 0}"]
b1 = boxes["{'sluices_widening': 2, 'summer_target': 0.2, 'pump_capacity': 500}"]
boxes
b.peeling_trajectory.head(15)
b.box_lims[0:10]
box.box_lims[i][:][lower_lim]
p1 = b.show_tradeoff()
b.show_pairs_scatter()
b1.show_tradeoff()
plt.show()
#b1.peeling_trajectory
#f = (b.peeling_trajectory.coverage + b.peeling_trajectory.density)/2
b1.show_pairs_scatter()
b.inspect(int(f.idxmax()))
b
b.inspect(int(f.idxmax()),style='graph')
b.inspect(11,style='graph')
box_id=None
box_res = b.inspect(4)
box_res2 = b.inspect(box_id,style='graph')<jupyter_output><empty_output><jupyter_text>Rendering adaptation pathways Creating xml to be loaded within the [Pathways Generator](https://publicwiki.deltares.nl/display/AP/Pathways+Generator) by [Deltares](www.deltares/nl)<jupyter_code>from pathways.pathways import create_patwhays_map
# find correspondence policy as in table to dictionary of actions
correspondence_table = get_policy_def(results_df,actions_name)
# auxiliary functions 
def find_policy_id(policy_dict,correspondence_table):
    for policy_id in correspondence_table:
        if correspondence_table[policy_id]==policy_dict:
            return(policy_id)

def find_policy(elements,precedent_policy,correspondence_table):
    
    for action in elements:
        print (action)
        
        new_policy = precedent_policy.copy()
        
        if 'values' in action: # not current situation
            new_policy[action['name']]=action['values'][-1]

        if 'combined_to' in action: # is a pathway
            value_combined = [combined_action['values'][-1] for combined_action in elements if combined_action['name']==action['combined_to']][0]
            new_policy[action['combined_to']] = value_combined
            
        new_policy_id = find_policy_id(new_policy,correspondence_table)
        action['policy_id'] = new_policy_id
        
        
# add action for current situation

current_situation = {}
for action_name in actions_name:
    action_value = [action['values'][0] for action in actions if action['name']==action_name][0]
    current_situation[action_name] = action_value 

action_current_situation = {'caption':'current_situation', 'name': 'current_situation','color':'clGray','predecessor_id':'-1'}

#action_current_situation['policy_id'] = find_policy_id(current_situation,correspondence_table)
# add current action to list of actins
elements = [action_current_situation] + actions + pathways
find_policy(elements,current_situation,correspondence_table)<jupyter_output><empty_output><jupyter_text>select critical uncertainty to render<jupyter_code># select critical uncertainty to render
rendered_uncertainty = 'peak_inflow_multiplier'#'sea_level_shift' # 
pos_rend_unc = critical_uncertainties.index(rendered_uncertainty)
elements
# assign tippingpoint value for the selected critical uncertainty
for action in elements:
    action['tippingpointvalue'] = ATP_risk_averse['RelativeATP'].loc[action['policy_id']]
ATP_parsimony.loc[action['policy_id']]
# check results
for el in elements:
    print(el['caption'])
    print(el['tippingpointvalue'])<jupyter_output><empty_output><jupyter_text>Generate pathway file<jupyter_code>graph_settings = {'xaxistitle':'relative change','maxtippingpoint':1, 'fontsize': 20}
xml_str = create_patwhays_map(elements,graph_settings = graph_settings)

pathway_file = open('figures/Ijsselmeer/risk_averse_density90_ATP_08_Nov.pathways', 'w')
pathway_file.write(xml_str)
pathway_file.close()
pathway<jupyter_output><empty_output><jupyter_text>Identify Adaptation Tipping Points, including observational uncertainty  Find ATP on the space of observable variables<jupyter_code>obs_unc_low_flow = 10
obs_unc_peak_discharge = 100 
obs_unc_sea_level = 0.1



obs_unc_comatrix = np.matrix([obs_unc_low_flow,0,0],[0,obs_unc_peak_discharge,0],[0,0,obs_unc_sea_level])
actions
import graphviz 
from graphviz_tree import *
test_new = False
if test_new == True:
    export_graphviz_func = export_graphviz_temp
else:
    export_graphviz_func = tree.export_graphviz
policies_names = []
for policy in policy_success_ordinal.columns:
    policies_names.append(policy.replace(',','\n'))
class_names = list(policy_success_ordered.columns)
class_names.append('None')
class_names = [policy_name.replace(',','\n') for policy_name in class_names ]
pathways = export_graphviz_func(clf, out_file=None, #max_depth=4,
                                feature_names=scenarios_df.columns, 
                                class_names=class_names, 
                                filled=True,
                                label = 'all',
                                leaves_parallel =False,
                                #node_ids=True,
                                #rounded=True,
                                #impurity = True, 
                                #proportion = True,
                                rotate = True,
                                precision = 3
                                #,minimal_info=True
                               )  
pathways_rendered = graphviz.Source(pathways)
pathways_rendered

cmap_name = 'my_list'
for n_bin, ax in zip(n_bins, axs.ravel()):
    # Create the colormap
    cm = LinearSegmentedColormap.from_list(
        cmap_name, colors, N=n_bin)
    # Fewer bins will result in "coarser" colomap interpolation
    im = ax.imshow(Z, interpolation='nearest', origin='lower', cmap=cm)
    ax.set_title("N bins: %s" % n_bin)
    fig.colorbar(im)
import pylab as pl
from matplotlib.colors import LinearSegmentedColormap


max_cost = 340 * 1.3
a = np.array([[0,max_cost]])
pl.figure(figsize=(3, 10))

colors = [(0, 0.7, 0),(0.7, 0, 0)]  # G -> R
n_bins = 30
color_map = LinearSegmentedColormap.from_list(cmap_name, colors, N=n_bins)
img = pl.imshow(a, cmap=color_map)
pl.gca().set_visible(False)
cax = pl.axes([0, 1, 0.1, 1])
pl.colorbar(cax=cax)
plt.ylabel('Cost $[M- euro]$')
pl.savefig("colorbar.pdf")

cost_policy = 100
color_for_policy = color_map(cost_policy/max_cost)[0:3]
color_for_policy
from sklearn.datasets import load_iris
iris = load_iris()
clf_iris = tree.DecisionTreeClassifier()
clf_iris.fit(iris.data, iris.target)

dot_data = tree.export_graphviz(clf, out_file=None, 
                         feature_names=iris.feature_names,  
                         class_names=iris.target_names,  
                         filled=True, rounded=True,  
                         special_characters=True)  
graph = graphviz.Source(dot_data)  
graph 
iris.target_names
# Parameters
n_classes = 3
plot_colors = "ryb"
plot_step = 0.02

# Load data
from sklearn.datasets import load_iris


iris = load_iris()

plt.figure(figsize=(10,7))

for pairidx, pair in enumerate([[0, 1], [0, 2], [0, 3],
                                [1, 2], [1, 3], [2, 3]]):
    # We only take the two corresponding features
    X = iris.data[:, pair]
    y = iris.target

    # Train
    clf = tree.DecisionTreeClassifier().fit(X, y)

    # Plot the decision boundary
    plt.subplot(2, 3, pairidx + 1)

    x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx, yy = np.meshgrid(np.arange(x_min, x_max, plot_step),
                         np.arange(y_min, y_max, plot_step))
    plt.tight_layout(h_pad=0.5, w_pad=0.5, pad=2.5)

    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])
    Z = Z.reshape(xx.shape)
    cs = plt.contourf(xx, yy, Z, cmap=plt.cm.RdYlBu)

    plt.xlabel(iris.feature_names[pair[0]])
    plt.ylabel(iris.feature_names[pair[1]])

    # Plot the training points
    for i, color in zip(range(n_classes), plot_colors):
        idx = np.where(y == i)
        plt.scatter(X[idx, 0], X[idx, 1], c=color, label=iris.target_names[i],
                    cmap=plt.cm.RdYlBu, edgecolor='black', s=15)
        
plt.suptitle("Decision surface of a decision tree using paired features")
plt.legend(loc=0)
#plt.axis("tight")
actions

all_actions<jupyter_output><empty_output><jupyter_text>Exploratory data analysis<jupyter_code># plot_area_success
pos = -1
print(policies[pos])
policy_under_analysis = pd.DataFrame(success_per_policy.iloc[:,pos].rename(columns={policies[pos]:'policy'}))
                                                              
scenarios_df_success = scenarios_df.join(policy_under_analysis)
    
sns.pairplot(scenarios_df_success,hue=0)
plt.title(policies[pos])
policies
var1 = 'sea_level_shift'
var2 = 'peak_inflow_multiplier'

policy = "{'winter_target': -0.4, 'sluices_widening': 1}"

#for policy in policies:
for policy in [policy]:
    
    scenarios_vs_success = scenarios_df.join(success_per_policy[policy].apply(int))
    pivoted_success = scenarios_vs_success.pivot_table(index=var1,columns=var2,values=policy,aggfunc=np.sum) 
    pivoted_total = scenarios_vs_success.pivot_table(index=var1,columns=var2,values=policy,aggfunc='count')
    frequency_success = pivoted_success / pivoted_total
    # plotting
    x,y = frequency_success.index, frequency_success.columns
    X,Y = np.meshgrid(x,y)
    levelsf=[0.2,1]
    levels=[0.2,0.5,0.8]
    linestyles = ['dotted','solid','dotted']
    #plt1 = plt.contourf(X, Y, frequency_success.values,levelsf,label=policy,colors=policy_color[policy],linestyles = linestyles, linewidth = 5 ,alpha=0.4)#, cmap=plt.cm.jet);
    plt2 = plt.contour(X, Y, frequency_success.values,levels,label=policy,colors=policy_color[policy],linestyles = linestyles, linewidth = 5 )#, cmap=plt.cm.jet);
    
plt.xlabel(var1)
plt.ylabel(var2)

scenarios_df.head()
scenarios_df,success_per_policy.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12368701124427375, 'fraction_numerical': 0.013280273934570357, 'mean_word_length': 4.654631083202512, 'pattern_counts': {'":': 0, '<': 87, '<?xml version=': 0, '>': 88, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 3}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '76859', 'license': 'MIT', 'n_tokens_mistral': 7703, 'n_tokens_neox': 7310, 'n_words': 1323, 'path': '3_Ijsselmeer/adaptive_ijsselmeer.ipynb', 'repo_name': 'luciofaso/Monitoring_DAP', 'size': '21627'}