<jupyter_start><jupyter_text>.. py:currentmodule:: emat<jupyter_code># Road Test<jupyter_output><empty_output><jupyter_text>import emat, os, numpy, pandas, functools, asyncioemat.versions() logger = emat.util.loggers.log_to_stderr(30, True)<jupyter_code>## Defining the Exploratory Scope<jupyter_output><empty_output><jupyter_text>The model scope is defined in a YAML file.  For this Road Test example, the scope file is named :ref:`road_test.yaml ` and is included in the model/tests directory.<jupyter_code>road_test_scope_file = emat.package_file('model','tests','road_test.yaml')<jupyter_output><empty_output><jupyter_text>The filename for the YAML file is the first argument when creating a :class:`Scope`object, which will load and process the content of the file.<jupyter_code>road_scope = emat.Scope(road_test_scope_file)
road_scope<jupyter_output><empty_output><jupyter_text>A short summary of the scope can be reviewed using the `info` method.<jupyter_code>road_scope.info()<jupyter_output><empty_output><jupyter_text>Alternatively, more detailed information about each part of the scope can beaccessed in four list attributes:<jupyter_code>road_scope.get_constants()
road_scope.get_uncertainties()
road_scope.get_levers()
road_scope.get_measures()<jupyter_output><empty_output><jupyter_text>Using a DatabaseThe exploratory modeling process will typically generate many different sets of outputs,for different explored modeling scopes, or for different applications.  It is convenientto organize these outputs in a database structure, so they are stored consistently and readily available for subsequent analysis.The `SQLiteDB` object will create a database to store results.  When instantiated withno arguments, the database is initialized in-memory, which will not store anything todisk (which is convenient for this example, but in practice you will generally want tostore data to disk so that it can persist after this Python session ends).<jupyter_code>emat_db = emat.SQLiteDB('tempfile')<jupyter_output><empty_output><jupyter_text>An EMAT Scope can be stored in the database, to provide needed information about what the various inputs and outputs represent.<jupyter_code>road_scope.store_scope(emat_db)<jupyter_output><empty_output><jupyter_text>Trying to store another scope with the same name (or the same scope) raises a KeyError.<jupyter_code>try:
    road_scope.store_scope(emat_db)
except KeyError as err:
    print(err)<jupyter_output><empty_output><jupyter_text>We can review the names of scopes already stored in the database using the `read_scope_names` method.<jupyter_code>emat_db.read_scope_names()<jupyter_output><empty_output><jupyter_text>Experimental DesignActually running the model can be done by the user on an *ad hoc* basis (i.e., manually defining every combination of inputs that will be evaluated) but the real power of EMAT comes from runnning the modelusing algorithm-created experimental designs.An important experimental design used in exploratory modeling is the Latin Hypercube.  This design selectsa random set of experiments across multiple input dimensions, to ensure "good" coverage of the multi-dimensional modeling space.The `design_latin_hypercube` function creates such a design based on a `Scope`, and optionallystores the design of experiments in a database.<jupyter_code>from emat.experiment.experimental_design import design_experiments
design = design_experiments(road_scope, db=emat_db, n_samples_per_factor=10, sampler='lhs')
design.head()
large_design = design_experiments(road_scope, db=emat_db, n_samples=5000, sampler='lhs', design_name='lhs_large')
large_design.head()<jupyter_output><empty_output><jupyter_text>We can review what experimental designs have already been stored in the database using the `read_design_names` method of the `Database` object.<jupyter_code>emat_db.read_design_names('EMAT Road Test')<jupyter_output><empty_output><jupyter_text>Core Model in Python<jupyter_code>Up until this point, we have been considering a model in the abstract, defining in the :class:`Scope` what the inputs 
and outputs will be, and designing some experiments we would like to run with the model.  Now we will actually 
interface with the model itself.<jupyter_output><empty_output><jupyter_text>Model DefinitionIn the simplest approach for EMAT, a model can be defined as a basic Python function, which accepts allinputs (exogenous uncertainties, policy levers, and externally defined constants) as named keywordarguments, and returns a dictionary where the dictionary keys are names of performace measures, and the mapped values are the computed values for those performance measures.  The `Road_Capacity_Investment`function provided in EMAT is an example of such a function.  This made-up example considers the investment in capacity expansion for a single roadway link.  The inputs to this function are describedabove in the Scope, including uncertain parameters in the volume-delay function,traffic volumes, value of travel time savings, unit construction costs, and interest rates, and policy levers including the amount of capacity expansion and amortization period.<jupyter_code>from emat.model.core_python import PythonCoreModel, Road_Capacity_Investment<jupyter_output><empty_output><jupyter_text>The :class:`PythonCoreModel ` object provides an interface that links the basic Python function that represents the model, the :class:`Scope `, and optionally the :class:`Database ` used to manage data storage.<jupyter_code>m = PythonCoreModel(Road_Capacity_Investment, scope=road_scope, db=emat_db)<jupyter_output><empty_output><jupyter_text>From the :class:`PythonCoreModel`, which links the model, scope, design, and database, we can run the design of experiments.  This will systematically run the core model with each set of input parameters in the design, store the results inthe database, and return a pandas.DataFrame containing the results.<jupyter_code>### Model Execution<jupyter_output><empty_output><jupyter_text>lhs_results = m.run_experiments(design_name='lhs')lhs_results.head()<jupyter_code>If running a large number of experiments, it may be valuable to parallelize the 
processing using a DistributedEvaluator instead of the default SequentialEvaluator.
The DistributedEvaluator uses dask.distributed to distribute the workload to
a cluster of processes, which can all be on the same machine or distributed
over multiple networked computers. (The details of using dask.distributed in 
more complex environments are beyond this scope of this example, but interested
users can refer to that package's [documentation](https://distributed.dask.org/).)<jupyter_output><empty_output><jupyter_text>lhs_large_async = m.async_experiments(large_design, max_n_workers=8, batch_size=157) lhs_large_results = await lhs_large_async.final_results()<jupyter_code>Once a particular design has been run once, the results can be recovered from the database without re-running the model itself.<jupyter_output><empty_output><jupyter_text>reload_results = m.read_experiments(design_name='lhs')reload_results.head()<jupyter_code>It is also possible to load only the parameters, or only the performance meausures.<jupyter_output><empty_output><jupyter_text>lhs_params = m.read_experiment_parameters(design_name='lhs')lhs_params.head() lhs_outcomes = m.read_experiment_measures(design_name='lhs')lhs_outcomes.head()<jupyter_code>## Feature Scoring<jupyter_output><empty_output><jupyter_text>m.get_feature_scores('lhs')<jupyter_code>## Visualization<jupyter_output><empty_output><jupyter_text>from emat.analysis import display_experimentsdisplay_experiments(road_scope, lhs_results, rows=['time_savings', 'net_benefits', 'input_flow'])<jupyter_code>## Scenario Discovery
Scenario discovery in exploratory modeling is focused on finding scenarios that are interesting to the user.  
The process generally begins through the identification of particular outcomes that are "of interest",
and the discovery process that can seek out what factor or combination of factors can result in
those outcomes.

There are a variety of methods to use for scenario discovery.  We illustrate a few here.

### PRIM

Patient rule induction method (PRIM) is an algorithm that operates on a set of data with inputs and outputs.  
It is used for locating areas of an outcome space that are of particular interest, which it does by reducing 
the data size incrementally by small amounts in an iterative process as follows:
    
- Candidate boxes are generated.  These boxes represent incrementally smaller sets of the data.  
  Each box removes a portion of the data based on the levels of a single input variable.
  * For categorical input variables, there is one box generated for each category with each box 
    removing one category from the data set.
  * For integer and continuous variables, two boxes are generated – one box that removes a 
    portion of data representing the smallest set of values for that input variable and another 
    box that removes a portion of data representing the largest set of values for that input.  
    The step size for these variables is controlled by the analyst.
- For each candidate box, the relative improvement in the number of outcomes of interest inside 
  the box is calculated and the candidate box with the highest improvement is selected.
- The data in the selected candidate box replaces the starting data and the process is repeated.

The process ends based on a stopping criteria.  For more details on the algorithm, 
see [<NAME> (1999)](http://statweb.stanford.edu/~jhf/ftp/prim.pdf) or 
[Kwakkel and Jaxa-Rozen (2016)](https://www.sciencedirect.com/science/article/pii/S1364815215301092).

The PRIM algorithm is particularly useful for scenario discovery, which broadly is the process of 
identifying particular scenarios of interest in a large and deeply uncertain dataset.   
In the context of exploratory modeling, scenario discovery is often used to obtain a better understanding 
of areas of the uncertainty space where a policy or collection of policies performs poorly because it is 
often used in tandem with robust search methods for identifying policies that perform well 
([Kwakkel and Jaxa-Rozen (2016)](https://www.sciencedirect.com/science/article/pii/S1364815215301092)).<jupyter_output><empty_output><jupyter_text>from emat.analysis.prim import Prim of_interest = lhs_large_results['net_benefits']>0 discovery = Prim(    m.read_experiment_parameters(design_name='lhs_large'),    of_interest,    scope=road_scope,) box1 = discovery.find_box() box1.tradeoff_selector() box1.inspect(45) box1.select(45) box1.splom()<jupyter_code>### CART

Classification and Regression Trees (CART) can also be used for scenario discovery. 
They partition the explored space (i.e., the scope) into a number of sections, with each partition
being added in such a way as to maximize the difference between observations on each 
side of the newly added partition divider, subject to some constraints.<jupyter_output><empty_output><jupyter_text>from emat.workbench.analysis import cartcart_alg = cart.CART(    m.read_experiment_parameters(design_name='lhs_large'),    of_interest,)cart_alg.build_tree() from emat.util.xmle import ShowShow(cart_alg.show_tree(format='svg')) cart_alg.boxes_to_dataframe(include_stats=True)<jupyter_code>## Robust Search<jupyter_output><empty_output><jupyter_text>from emat import MeasureMAXIMIZE = Measure.MAXIMIZEMINIMIZE = Measure.MINIMIZErobustness_functions = [    Measure(        'Expected Net Benefit',         kind=Measure.INFO,         variable_name='net_benefits',         function=numpy.mean,    ),        Measure(        'Probability of Net Loss',         kind=MINIMIZE,         variable_name='net_benefits',         function=lambda x: numpy.mean(x<0),    ),    Measure(        '95%ile Travel Time',         kind=MINIMIZE,         variable_name='build_travel_time',        function=functools.partial(numpy.percentile, q=95),    ),    Measure(        '99%ile Present Cost',         kind=Measure.INFO,         variable_name='present_cost_expansion',         function=functools.partial(numpy.percentile, q=99),    ),    Measure(        'Expected Present Cost',         kind=Measure.INFO,         variable_name='present_cost_expansion',         function=numpy.mean,    ),]<jupyter_code>### Constraints

The robust optimization process solutions can be constrained to only include solutions that 
satisfy certain constraints.  These constraints can be based on the policy lever parameters
that are contained in the core model, the aggregate performance measures identified in 
the list of robustness functions, or some combination of levers and aggregate measures.<jupyter_output><empty_output><jupyter_text>from emat import Constraint<jupyter_code>The common use case for constraints in robust optimation is imposing requirements
on solution outcomes. For example, we may want to limit our robust search only
to solutions where the expected present cost of the capacity expansion is less
than some particular value (in our example here, 4000).<jupyter_output><empty_output><jupyter_text>constraint_1 = Constraint(    "Maximum Log Expected Present Cost",     outcome_names="Expected Present Cost",    function=Constraint.must_be_less_than(4000),)<jupyter_code>Our second constraint is based exclusively on an input: the capacity expansion
must be at least 10.  We could also achieve this kind of constraint by changing
the exploratory scope, but we don't necessarily want to change the scope to 
conduct a single robust optimization analysis with a constraint on a policy lever.<jupyter_output><empty_output><jupyter_text>constraint_2 = Constraint(    "Minimum Capacity Expansion",     parameter_names="expand_capacity",    function=Constraint.must_be_greater_than(10),)<jupyter_code>It is also possible to impose constraints based on a combination of inputs and outputs.
For example, suppose that the total funds available for pay-as-you-go financing are
only 1500.  We may thus want to restrict the robust search to only solutions that
are almost certainly within the available funds at 99% confidence (a model output) but only 
if the Paygo financing option is used (a model input).  This kind of constraint can
be created by giving both `parameter_names` and `outcomes_names`, and writing a constraint
function that takes two arguments.<jupyter_output><empty_output><jupyter_text>constraint_3 = Constraint(    "Maximum Paygo",     parameter_names='debt_type',    outcome_names='99%ile Present Cost',    function=lambda i,j: max(0, j-1500) if i=='Paygo' else 0,) robust_results = m.robust_optimize(    robustness_functions,      scenarios=200,     nfe=2500,     constraints=[        constraint_1,        constraint_2,        constraint_3,    ],    evaluator=get_client(),    cache_file="./cache_road_test_robust_opt")<jupyter_code>The robust_results include all the non-dominated solutions which satisfy the contraints that are found.
Note that robust optimization is generally a hard problem, and the algorithm may need to run for
a very large number of iterations in order to arrive at a good set of robust solutions.<jupyter_output><empty_output><jupyter_text>robust_results.par_coords()<jupyter_code>## Creating a Meta-Model<jupyter_output><empty_output><jupyter_text>Creating a meta-model requires an existing model, plus a set of experiments (including inputs and outputs) run with that modelto support the estimation of the meta-model.  After having completedsufficient initial runs of the core model, instantiating a meta-modelis as simple as calling a `create_metamodel_*` method on the coremodel, either giving a design of experiments stored in the databasewith :meth:`create_metamodel_from_design `or passing the experimental results directly with:meth:`create_metamodel_from_data `.<jupyter_code>mm = m.create_metamodel_from_design('lhs')
mm<jupyter_output><empty_output><jupyter_text>To demonstrate the performance of the meta-model, we can create analternate design of experiments.  Note that to get different random values,we set the `random_seed` argument to something other than the default value.<jupyter_code>design2 = design_experiments(
    scope=road_scope, 
    db=emat_db, 
    n_samples_per_factor=10, 
    sampler='lhs', 
    random_seed=2,
)
design2_results = mm.run_experiments(design2)
design2_results.head()
mm.cross_val_scores()<jupyter_output><empty_output><jupyter_text>Compare Core vs Meta Model ResultsWe can generate a variety of plots to compare the distribution of meta-model outcomeson the new design against the original model's results.<jupyter_code>from emat.analysis import contrast_experiments
contrast_experiments(road_scope, lhs_results, design2_results)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0654144305307096, 'fraction_numerical': 0.00727489564698867, 'mean_word_length': 6.037767519932857, 'pattern_counts': {'":': 0, '<': 135, '<?xml version=': 0, '>': 135, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '33', 'dir': 'jupyter-structured-clean-dedup', 'id': '661080', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4694, 'n_tokens_neox': 4489, 'n_words': 1736, 'path': 'docs/source/emat.examples/RoadTest/road_test_direct.ipynb', 'repo_name': 'jinsanity07git/tmip-emat', 'size': '16783'}