Tutorials

The code of these examples can be found in the examples package. The first three examples are meant to illustrate the basics of the EMA workbench. How to implement a model, specify its uncertainties and outcomes, and run it. The fourth example is a more extensive illustration based on Pruyt & Hamarat (2010). It shows some more advanced possibilities of the EMA workbench, including one way of handling policies.

A simple model in Python

The simplest case is where we have a model available through a python function. For example, imagine we have the simple model.

def some_model(x1=None, x2=None, x3=None):
    return {'y':x1*x2+x3}

In order to control this model from the workbench, we can make use of the Model. We can instantiate a model object, by passing it a name, and the function.

model = Model('simpleModel', function=some_model) #instantiate the model

Next, we need to specify the uncertainties and the outcomes of the model. In this case, the uncertainties are x1, x2, and x3, while the outcome is y. Both uncertainties and outcomes are attributes of the model object, so we can say

1#specify uncertainties
2model.uncertainties = [RealParameter("x1", 0.1, 10),
3                       RealParameter("x2", -0.01,0.01),
4                       RealParameter("x3", -0.01,0.01)]
5#specify outcomes
6model.outcomes = [ScalarOutcome('y')]

Here, we specify that x1 is some value between 0.1, and 10, while both x2 and x3 are somewhere between -0.01 and 0.01. Having implemented this model, we can now investigate the model behavior over the set of uncertainties by simply calling

results = perform_experiments(model, 100)

The function perform_experiments() takes the model we just specified and will execute 100 experiments. By default, these experiments are generated using a Latin Hypercube sampling, but Monte Carlo sampling and Full factorial sampling are also readily available. Read the documentation for perform_experiments() for more details.

The complete code:

 1"""
 2Created on 20 dec. 2010
 3
 4This file illustrated the use the EMA classes for a contrived example
 5It's main purpose has been to test the parallel processing functionality
 6
 7.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
 8"""
 9
10from ema_workbench import Model, RealParameter, ScalarOutcome, ema_logging, perform_experiments
11
12
13def some_model(x1=None, x2=None, x3=None):
14    return {"y": x1 * x2 + x3}
15
16
17if __name__ == "__main__":
18    ema_logging.LOG_FORMAT = "[%(name)s/%(levelname)s/%(processName)s] %(message)s"
19    ema_logging.log_to_stderr(ema_logging.INFO)
20
21    model = Model("simpleModel", function=some_model)  # instantiate the model
22
23    # specify uncertainties
24    model.uncertainties = [
25        RealParameter("x1", 0.1, 10),
26        RealParameter("x2", -0.01, 0.01),
27        RealParameter("x3", -0.01, 0.01),
28    ]
29    # specify outcomes
30    model.outcomes = [ScalarOutcome("y")]
31
32    results = perform_experiments(model, 100)

A simple model in Vensim

Imagine we have a very simple Vensim model:

_images/simpleVensimModel.png

For this example, we assume that ‘x11’ and ‘x12’ are uncertain. The state variable ‘a’ is the outcome of interest. Similar to the previous example, we have to first instantiate a vensim model object, in this case VensimModel. To this end, we need to specify the directory in which the vensim file resides, the name of the vensim file and the name of the model.

wd = r'./models/vensim example'
model = VensimModel("simpleModel", wd=wd, model_file=r'\model.vpm')

Next, we can specify the uncertainties and the outcomes.

1model.uncertainties = [RealParameter("x11", 0, 2.5),
2                       RealParameter("x12", -2.5, 2.5)]
3
4
5model.outcomes = [TimeSeriesOutcome('a')]

Note that we are using a TimeSeriesOutcome, because vensim results are time series. We can now simply run this model by calling perform_experiments().

with MultiprocessingEvaluator(model) as evaluator:
results = evaluator.perform_experiments(1000)

We now use a evaluator, which ensures that the code is executed in parallel.

Is it generally good practice to first run a model a small number of times sequentially prior to running in parallel. In this way, bugs etc. can be spotted more easily. To further help with keeping track of what is going on, it is also good practice to make use of the logging functionality provided by the workbench

ema_logging.log_to_stderr(ema_logging.INFO)

Typically, this line appears at the start of the script. When executing the code, messages on progress or on errors will be shown.

The complete code

 1"""
 2Created on 3 Jan. 2011
 3
 4This file illustrated the use the EMA classes for a contrived vensim
 5example
 6
 7
 8.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
 9                chamarat <c.hamarat (at) tudelft (dot) nl>
10"""
11
12from ema_workbench import TimeSeriesOutcome, perform_experiments, RealParameter, ema_logging
13
14from ema_workbench.connectors.vensim import VensimModel
15
16if __name__ == "__main__":
17    # turn on logging
18    ema_logging.log_to_stderr(ema_logging.INFO)
19
20    # instantiate a model
21    wd = "./models/vensim example"
22    vensimModel = VensimModel("simpleModel", wd=wd, model_file="model.vpm")
23    vensimModel.uncertainties = [RealParameter("x11", 0, 2.5), RealParameter("x12", -2.5, 2.5)]
24
25    vensimModel.outcomes = [TimeSeriesOutcome("a")]
26
27    results = perform_experiments(vensimModel, 1000)

A simple model in Excel

In order to perform EMA on an Excel model, one can use the ExcelModel. This base class makes uses of naming cells in Excel to refer to them directly. That is, we can assume that the names of the uncertainties correspond to named cells in Excel, and similarly, that the names of the outcomes correspond to named cells or ranges of cells in Excel. When using this class, make sure that the decimal separator and thousands separator are set correctly in Excel. This can be checked via file > options > advanced. These separators should follow the anglo saxon convention.

 1"""
 2Created on 27 Jul. 2011
 3
 4This file illustrated the use the EMA classes for a model in Excel.
 5
 6It used the excel file provided by
 7`A. Sharov <https://home.comcast.net/~sharov/PopEcol/lec10/fullmod.html>`_
 8
 9This excel file implements a simple predator prey model.
10
11.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
12"""
13
14from ema_workbench import RealParameter, TimeSeriesOutcome, ema_logging, perform_experiments
15
16from ema_workbench.connectors.excel import ExcelModel
17from ema_workbench.em_framework.evaluators import MultiprocessingEvaluator
18
19if __name__ == "__main__":
20    ema_logging.log_to_stderr(level=ema_logging.INFO)
21
22    model = ExcelModel("predatorPrey", wd="./models/excelModel", model_file="excel example.xlsx")
23    model.uncertainties = [
24        RealParameter("K2", 0.01, 0.2),
25        # we can refer to a cell in the normal way
26        # we can also use named cells
27        RealParameter("KKK", 450, 550),
28        RealParameter("rP", 0.05, 0.15),
29        RealParameter("aaa", 0.00001, 0.25),
30        RealParameter("tH", 0.45, 0.55),
31        RealParameter("kk", 0.1, 0.3),
32    ]
33
34    # specification of the outcomes
35    model.outcomes = [
36        TimeSeriesOutcome("B4:B1076"),
37        # we can refer to a range in the normal way
38        TimeSeriesOutcome("P_t"),
39    ]  # we can also use named range
40
41    # name of the sheet
42    model.default_sheet = "Sheet1"
43
44    with MultiprocessingEvaluator(model) as evaluator:
45        results = perform_experiments(model, 100, reporting_interval=1, evaluator=evaluator)

The example is relatively straight forward. We instantiate an excel model, we specify the uncertainties and the outcomes. We also need to specify the sheet in excel on which the model resides. Next we can call perform_experiments().

Warning

when using named cells. Make sure that the names are defined at the sheet level and not at the workbook level

A more elaborate example: Mexican Flu

This example is derived from Pruyt & Hamarat (2010). This paper presents a small exploratory System Dynamics model related to the dynamics of the 2009 flu pandemic, also known as the Mexican flu, swine flu, or A(H1N1)v. The model was developed in May 2009 in order to quickly foster understanding about the possible dynamics of this new flu variant and to perform rough-cut policy explorations. Later, the model was also used to further develop and illustrate Exploratory Modelling and Analysis.

Mexican Flu: the basic model

In the first days, weeks and months after the first reports about the outbreak of a new flu variant in Mexico and the USA, much remained unknown about the possible dynamics and consequences of the at the time plausible/imminent epidemic/pandemic of the new flu variant, first known as Swine or Mexican flu and known today as Influenza A(H1N1)v.

The exploratory model presented here is small, simple, high-level, data-poor (no complex/special structures nor detailed data beyond crude guestimates), and history-poor.

The modelled world is divided in three regions: the Western World, the densely populated Developing World, and the scarcely populated Developing World. Only the two first regions are included in the model because it is assumed that the scarcely populated regions are causally less important for dynamics of flu pandemics. Below, the figure shows the basic stock-and-flow structure. For a more elaborate description of the model, see Pruyt & Hamarat (2010).

_images/flu-model.png

Given the various uncertainties about the exact characteristics of the flu, including its fatality rate, the contact rate, the susceptibility of the population, etc. the flu case is an ideal candidate for EMA. One can use EMA to explore the kinds of dynamics that can occur, identify undesirable dynamic, and develop policies targeted at the undesirable dynamics.

In the original paper, Pruyt & Hamarat (2010). recoded the model in Python and performed the analysis in that way. Here we show how the EMA workbench can be connected to Vensim directly.

The flu model was build in Vensim. We can thus use VensimModelS as a base class.

We are interested in two outcomes:

  • deceased population region 1: the total number of deaths over the duration of the simulation.

  • peak infected fraction: the fraction of the population that is infected.

These are added to self.outcomes, using the TimeSeriesOutcome class.

The table below is adapted from Pruyt & Hamarat (2010). It shows the uncertainties, and their bounds. These are added to self.uncertainties as ParameterUncertainty instances.

Parameter

Lower Limit

Upper Limit

additional seasonal immune population fraction region 1

0.0

0.5

additional seasonal immune population fraction region 2

0.0

0.5

fatality ratio region 1

0.0001

0.1

fatality ratio region 2

0.0001

0.1

initial immune fraction of the population of region 1

0.0

0.5

initial immune fraction of the population of region 2

0.0

0.5

normal interregional contact rate

0.0

0.9

permanent immune population fraction region 1

0.0

0.5

permanent immune population fraction region 2

0.0

0.5

recovery time region 1

0.2

0.8

recovery time region 2

0.2

0.8

root contact rate region 1

1.0

10.0

root contact rate region 2

1.0

10.0

infection ratio region 1

0.0

0.1

infection ratio region 2

0.0

0.1

normal contact rate region 1

10

200

normal contact rate region 2

10

200

Together, this results in the following code:

  1"""
  2Created on 20 May, 2011
  3
  4This module shows how you can use vensim models directly
  5instead of coding the model in Python. The underlying case
  6is the same as used in fluExample
  7
  8.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
  9                epruyt <e.pruyt (at) tudelft (dot) nl>
 10"""
 11
 12from ema_workbench import (
 13    RealParameter,
 14    TimeSeriesOutcome,
 15    ema_logging,
 16    perform_experiments,
 17    MultiprocessingEvaluator,
 18    save_results,
 19)
 20
 21from ema_workbench.connectors.vensim import VensimModel
 22
 23if __name__ == "__main__":
 24    ema_logging.log_to_stderr(ema_logging.INFO)
 25
 26    model = VensimModel("fluCase", wd=r"./models/flu", model_file=r"FLUvensimV1basecase.vpm")
 27
 28    # outcomes
 29    model.outcomes = [
 30        TimeSeriesOutcome(
 31            "deceased_population_region_1", variable_name="deceased population region 1"
 32        ),
 33        TimeSeriesOutcome("infected_fraction_R1", variable_name="infected fraction R1"),
 34    ]
 35
 36    # Plain Parametric Uncertainties
 37    model.uncertainties = [
 38        RealParameter(
 39            "additional_seasonal_immune_population_fraction_R1",
 40            0,
 41            0.5,
 42            variable_name="additional seasonal immune population fraction R1",
 43        ),
 44        RealParameter(
 45            "additional_seasonal_immune_population_fraction_R2",
 46            0,
 47            0.5,
 48            variable_name="additional seasonal immune population fraction R2",
 49        ),
 50        RealParameter(
 51            "fatality_ratio_region_1", 0.0001, 0.1, variable_name="fatality ratio region 1"
 52        ),
 53        RealParameter(
 54            "fatality_rate_region_2", 0.0001, 0.1, variable_name="fatality rate region 2"
 55        ),
 56        RealParameter(
 57            "initial_immune_fraction_of_the_population_of_region_1",
 58            0,
 59            0.5,
 60            variable_name="initial immune fraction of the population of region 1",
 61        ),
 62        RealParameter(
 63            "initial_immune_fraction_of_the_population_of_region_2",
 64            0,
 65            0.5,
 66            variable_name="initial immune fraction of the population of region 2",
 67        ),
 68        RealParameter(
 69            "normal_interregional_contact_rate",
 70            0,
 71            0.9,
 72            variable_name="normal interregional contact rate",
 73        ),
 74        RealParameter(
 75            "permanent_immune_population_fraction_R1",
 76            0,
 77            0.5,
 78            variable_name="permanent immune population fraction R1",
 79        ),
 80        RealParameter(
 81            "permanent_immune_population_fraction_R2",
 82            0,
 83            0.5,
 84            variable_name="permanent immune population fraction R2",
 85        ),
 86        RealParameter("recovery_time_region_1", 0.1, 0.75, variable_name="recovery time region 1"),
 87        RealParameter("recovery_time_region_2", 0.1, 0.75, variable_name="recovery time region 2"),
 88        RealParameter(
 89            "susceptible_to_immune_population_delay_time_region_1",
 90            0.5,
 91            2,
 92            variable_name="susceptible to immune population delay time region 1",
 93        ),
 94        RealParameter(
 95            "susceptible_to_immune_population_delay_time_region_2",
 96            0.5,
 97            2,
 98            variable_name="susceptible to immune population delay time region 2",
 99        ),
100        RealParameter(
101            "root_contact_rate_region_1", 0.01, 5, variable_name="root contact rate region 1"
102        ),
103        RealParameter(
104            "root_contact_ratio_region_2", 0.01, 5, variable_name="root contact ratio region 2"
105        ),
106        RealParameter(
107            "infection_ratio_region_1", 0, 0.15, variable_name="infection ratio region 1"
108        ),
109        RealParameter("infection_rate_region_2", 0, 0.15, variable_name="infection rate region 2"),
110        RealParameter(
111            "normal_contact_rate_region_1", 10, 100, variable_name="normal contact rate region 1"
112        ),
113        RealParameter(
114            "normal_contact_rate_region_2", 10, 200, variable_name="normal contact rate region 2"
115        ),
116    ]
117
118    nr_experiments = 1000
119    with MultiprocessingEvaluator(model) as evaluator:
120        results = perform_experiments(model, nr_experiments, evaluator=evaluator)
121
122    save_results(results, "./data/1000 flu cases no policy.tar.gz")

We have now instantiated the model, specified the uncertain factors and outcomes and run the model. We now have generated a dataset of results and can proceed to analyse the results using various analysis scripts. As a first step, one can look at the individual runs using a line plot using lines(). See plotting for some more visualizations using results from performing EMA on FluModel.

1import matplotlib.pyplot as plt
2from ema_workbench.analysis.plotting import lines
3
4figure = lines(results, density=True) #show lines, and end state density
5plt.show() #show figure

generates the following figure:

_images/tutorial-lines.png

From this figure, one can deduce that across the ensemble of possible futures, there is a subset of runs with a substantial amount of deaths. We can zoom in on those cases, identify their conditions for occurring, and use this insight for policy design.

For further analysis, it is generally convenient, to generate the results for a series of experiments and save these results. One can then use these saved results in various analysis scripts.

from ema_workbench import save_results
save_results(results, './1000 runs.tar.gz')

The above code snippet shows how we can use save_results() for saving the results of our experiments. save_results() stores the as csv files in a tarball.

Mexican Flu: policies

For this paper, policies were developed by using the system understanding of the analysts.

static policy

adaptive policy

running the policies

In order to be able to run the models with the policies and to compare their results with the no policy case, we need to specify the policies

1#add policies
2policies = [Policy('no policy',
3                   model_file=r'/FLUvensimV1basecase.vpm'),
4            Policy('static policy',
5                   model_file=r'/FLUvensimV1static.vpm'),
6            Policy('adaptive policy',
7                   model_file=r'/FLUvensimV1dynamic.vpm')
8            ]

In this case, we have chosen to have the policies implemented in separate vensim files. Policies require a name, and can take any other keyword arguments you like. If the keyword matches an attribute on the model object, it will be updated, so model_file is an attribute on the vensim model. When executing the policies, we update this attribute for each policy. We can pass these policies to perform_experiment() as an additional keyword argument

results = perform_experiments(model, 1000, policies=policies)

We can now proceed in the same way as before, and perform a series of experiments. Together, this results in the following code:

  1"""
  2Created on 20 May, 2011
  3
  4This module shows how you can use vensim models directly
  5instead of coding the model in Python. The underlying case
  6is the same as used in fluExample
  7
  8.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
  9                epruyt <e.pruyt (at) tudelft (dot) nl>
 10"""
 11
 12import numpy as np
 13
 14from ema_workbench import (
 15    RealParameter,
 16    TimeSeriesOutcome,
 17    ema_logging,
 18    ScalarOutcome,
 19    perform_experiments,
 20    Policy,
 21    save_results,
 22)
 23from ema_workbench.connectors.vensim import VensimModel
 24
 25if __name__ == "__main__":
 26    ema_logging.log_to_stderr(ema_logging.INFO)
 27
 28    model = VensimModel("fluCase", wd=r"./models/flu", model_file=r"FLUvensimV1basecase.vpm")
 29
 30    # outcomes
 31    model.outcomes = [
 32        TimeSeriesOutcome(
 33            "deceased_population_region_1", variable_name="deceased population region 1"
 34        ),
 35        TimeSeriesOutcome("infected_fraction_R1", variable_name="infected fraction R1"),
 36        ScalarOutcome(
 37            "max_infection_fraction", variable_name="infected fraction R1", function=np.max
 38        ),
 39    ]
 40
 41    # Plain Parametric Uncertainties
 42    model.uncertainties = [
 43        RealParameter(
 44            "additional_seasonal_immune_population_fraction_R1",
 45            0,
 46            0.5,
 47            variable_name="additional seasonal immune population fraction R1",
 48        ),
 49        RealParameter(
 50            "additional_seasonal_immune_population_fraction_R2",
 51            0,
 52            0.5,
 53            variable_name="additional seasonal immune population fraction R2",
 54        ),
 55        RealParameter(
 56            "fatality_ratio_region_1", 0.0001, 0.1, variable_name="fatality ratio region 1"
 57        ),
 58        RealParameter(
 59            "fatality_rate_region_2", 0.0001, 0.1, variable_name="fatality rate region 2"
 60        ),
 61        RealParameter(
 62            "initial_immune_fraction_of_the_population_of_region_1",
 63            0,
 64            0.5,
 65            variable_name="initial immune fraction of the population of region 1",
 66        ),
 67        RealParameter(
 68            "initial_immune_fraction_of_the_population_of_region_2",
 69            0,
 70            0.5,
 71            variable_name="initial immune fraction of the population of region 2",
 72        ),
 73        RealParameter(
 74            "normal_interregional_contact_rate",
 75            0,
 76            0.9,
 77            variable_name="normal interregional contact rate",
 78        ),
 79        RealParameter(
 80            "permanent_immune_population_fraction_R1",
 81            0,
 82            0.5,
 83            variable_name="permanent immune population fraction R1",
 84        ),
 85        RealParameter(
 86            "permanent_immune_population_fraction_R2",
 87            0,
 88            0.5,
 89            variable_name="permanent immune population fraction R2",
 90        ),
 91        RealParameter("recovery_time_region_1", 0.1, 0.75, variable_name="recovery time region 1"),
 92        RealParameter("recovery_time_region_2", 0.1, 0.75, variable_name="recovery time region 2"),
 93        RealParameter(
 94            "susceptible_to_immune_population_delay_time_region_1",
 95            0.5,
 96            2,
 97            variable_name="susceptible to immune population delay time region 1",
 98        ),
 99        RealParameter(
100            "susceptible_to_immune_population_delay_time_region_2",
101            0.5,
102            2,
103            variable_name="susceptible to immune population delay time region 2",
104        ),
105        RealParameter(
106            "root_contact_rate_region_1", 0.01, 5, variable_name="root contact rate region 1"
107        ),
108        RealParameter(
109            "root_contact_ratio_region_2", 0.01, 5, variable_name="root contact ratio region 2"
110        ),
111        RealParameter(
112            "infection_ratio_region_1", 0, 0.15, variable_name="infection ratio region 1"
113        ),
114        RealParameter("infection_rate_region_2", 0, 0.15, variable_name="infection rate region 2"),
115        RealParameter(
116            "normal_contact_rate_region_1", 10, 100, variable_name="normal contact rate region 1"
117        ),
118        RealParameter(
119            "normal_contact_rate_region_2", 10, 200, variable_name="normal contact rate region 2"
120        ),
121    ]
122
123    # add policies
124    policies = [
125        Policy("no policy", model_file=r"FLUvensimV1basecase.vpm"),
126        Policy("static policy", model_file=r"FLUvensimV1static.vpm"),
127        Policy("adaptive policy", model_file=r"FLUvensimV1dynamic.vpm"),
128    ]
129
130    results = perform_experiments(model, 1000, policies=policies)
131    save_results(results, "./data/1000 flu cases with policies.tar.gz")

comparison of results

Using the following script, we reproduce figures similar to the 3D figures in Pruyt & Hamarat (2010). But using pairs_scatter(). It shows for the three different policies their behavior on the total number of deaths, the height of the highest peak of the pandemic, and the point in time at which this peak was reached.

 1"""
 2Created on 20 sep. 2011
 3
 4.. codeauthor:: jhkwakkel <j.h.kwakkel (at) tudelft (dot) nl>
 5"""
 6
 7import matplotlib.pyplot as plt
 8import numpy as np
 9
10from ema_workbench import load_results, ema_logging
11from ema_workbench.analysis.pairs_plotting import pairs_lines, pairs_scatter, pairs_density
12
13ema_logging.log_to_stderr(level=ema_logging.DEFAULT_LEVEL)
14
15# load the data
16fh = "./data/1000 flu cases no policy.tar.gz"
17experiments, outcomes = load_results(fh)
18
19# transform the results to the required format
20# that is, we want to know the max peak and the casualties at the end of the
21# run
22tr = {}
23
24# get time and remove it from the dict
25time = outcomes.pop("TIME")
26
27for key, value in outcomes.items():
28    if key == "deceased_population_region_1":
29        tr[key] = value[:, -1]  # we want the end value
30    else:
31        # we want the maximum value of the peak
32        max_peak = np.max(value, axis=1)
33        tr["max peak"] = max_peak
34
35        # we want the time at which the maximum occurred
36        # the code here is a bit obscure, I don't know why the transpose
37        # of value is needed. This however does produce the appropriate results
38        logical = value.T == np.max(value, axis=1)
39        tr["time of max"] = time[logical.T]
40
41pairs_scatter(experiments, tr, filter_scalar=False)
42pairs_lines(experiments, outcomes)
43pairs_density(experiments, tr, filter_scalar=False)
44plt.show()

no policy

_images/multiplot-flu-no-policy.png

static policy

_images/multiplot-flu-static-policy.png

adaptive policy

_images/multiplot-flu-adaptive-policy.png