<jupyter_start><jupyter_text>Open exploration with the Exploratory Modelling WorkbenchIn this blog, I will continue to showcase the functionality of the exploratory modelling. In the [previous blog](https://wordpress.com/post/waterprogramming.wordpress.com/14714), 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.  first a short backgroundIn 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.  open explorationTo 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).<jupyter_code>%%px --local

def process_p(values):
    values = np.asarray(values)
    values = np.mean(values, axis=0)
    return np.max(values)
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)]
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)<jupyter_output>[MainProcess/INFO] performing experiments using ipyparallel
[MainProcess/INFO] performing 100 scenarios * 10 policies * 1 model(s) = 1000 experiments
[MainProcess/INFO] 100 cases completed
[MainProcess/INFO] 200 cases completed
[MainProcess/INFO] 300 cases completed
[MainProcess/INFO] 400 cases completed
[MainProcess/INFO] 500 cases completed
[MainProcess/INFO] 600 cases completed
[MainProcess/INFO] 700 cases completed
[MainProcess/INFO] 800 cases completed
[MainProcess/INFO] 900 cases completed
[MainProcess/INFO] 1000 cases completed
[MainProcess/INFO] experiments finished<jupyter_text>Visual analysis<jupyter_code>experiments, outcomes = results
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()<jupyter_output><empty_output><jupyter_text>advanced analysis* scenario discovery* feature scoring<jupyter_code>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()
import mpld3
box1.show_tradeoff()
plt.show()<jupyter_output><empty_output><jupyter_text>linear regression<jupyter_code>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()<jupyter_output><empty_output><jupyter_text>feature scoring<jupyter_code>from ema_workbench.analysis import feature_scoring

x = experiments
y = outcomes

fs = feature_scoring.get_feature_scores_all(x, y)

sns.heatmap(fs, cmap='viridis', annot=True)
plt.show()<jupyter_output><empty_output><jupyter_text>More advanced sampling techniques<jupyter_code>with MultiprocessingEvaluator(model) as evaluator:
    sa_results = evaluator.perform_experiments(scenarios=1000, 
                                               uncertainty_sampling='sobol')
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)
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<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08579164885278609, 'fraction_numerical': 0.019951546244833973, 'mean_word_length': 3.6507620941020544, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 24, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '296026', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2124, 'n_tokens_neox': 1925, 'n_words': 623, 'path': 'docs/source/indepth_tutorial/blog2_ipyparallel.ipynb', 'repo_name': 'brodderickrodriguez/EMAworkbench', 'size': '7017'}