Plotting a Suite of Simulations
===============================

When we run a suite of simulations, such as when performing sensitivity
tests or uncertainty propagation, it can be handy to plot that suite of
runs together, and give the viewer a sense for how the uncertainty is
manifest over the course of the simulation.

Libraries
~~~~~~~~~

In addition to the standard plotting, simulation and data handling
libraries, we’ll use the ``seaborn`` plotting library to handle density
plots.

.. code:: ipython3

    %pylab inline
    import pysd
    import numpy as np
    import pandas as pd
    import seaborn


.. parsed-literal::

    Populating the interactive namespace from numpy and matplotlib


Loading Model
~~~~~~~~~~~~~

The model we’ll use in this example is a basic, 1-stock carbon bathtub
model, in which ``Emissions`` contribute to a stock of
``Excess Atmospheric Carbon``, which is slowly depleted through a
process of ``Natural Removal``.

|Atmospheric Bathtub|.

.. |Atmospheric Bathtub| image:: ../../../source/models/Climate/Atmospheric_Bathtub.png

.. code:: ipython3

    model = pysd.read_vensim('../../models/Climate/Atmospheric_Bathtub.mdl')

Creating a suite of run parameters
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We’ll generate a suite of simulations by drawing 1000 constant values
for the ``Emissions`` parameter from an exponential distribution.

.. code:: ipython3

    n_runs = 1000
    runs = pd.DataFrame({'Emissions': np.random.exponential(scale=10000, size=n_runs)})
    runs.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>Emissions</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>10619.175953</td>
        </tr>
        <tr>
          <th>1</th>
          <td>10854.493356</td>
        </tr>
        <tr>
          <th>2</th>
          <td>17405.296238</td>
        </tr>
        <tr>
          <th>3</th>
          <td>6649.127703</td>
        </tr>
        <tr>
          <th>4</th>
          <td>4454.821863</td>
        </tr>
      </tbody>
    </table>
    </div>



Run the model with the various parameters
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Next we’ll run the model with our various values for emissions, and
collect the resulting timeseries values of the stock of
``Excess Atmospheric Carbon``.

The resulting dataframe ``result`` contains a column for each value
simulation run, and these will form the traces for our plot.

.. code:: ipython3

    result = runs.apply(lambda p: model.run(params=dict(p))['Excess Atmospheric Carbon'],
                        axis=1).T
    result.head()




.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }
    
        .dataframe tbody tr th {
            vertical-align: top;
        }
    
        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>0</th>
          <th>1</th>
          <th>2</th>
          <th>3</th>
          <th>4</th>
          <th>5</th>
          <th>6</th>
          <th>7</th>
          <th>8</th>
          <th>9</th>
          <th>...</th>
          <th>990</th>
          <th>991</th>
          <th>992</th>
          <th>993</th>
          <th>994</th>
          <th>995</th>
          <th>996</th>
          <th>997</th>
          <th>998</th>
          <th>999</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>...</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
          <td>0.000000</td>
        </tr>
        <tr>
          <th>1</th>
          <td>10619.175953</td>
          <td>10854.493356</td>
          <td>17405.296238</td>
          <td>6649.127703</td>
          <td>4454.821863</td>
          <td>3560.003419</td>
          <td>6640.764213</td>
          <td>23972.691959</td>
          <td>5816.878223</td>
          <td>4160.733879</td>
          <td>...</td>
          <td>5354.787840</td>
          <td>24506.090733</td>
          <td>5391.310572</td>
          <td>5042.198805</td>
          <td>7717.170753</td>
          <td>13855.131882</td>
          <td>1254.482621</td>
          <td>3288.548916</td>
          <td>3995.330497</td>
          <td>689.274481</td>
        </tr>
        <tr>
          <th>2</th>
          <td>21132.160146</td>
          <td>21600.441778</td>
          <td>34636.539514</td>
          <td>13231.764129</td>
          <td>8865.095508</td>
          <td>7084.406803</td>
          <td>13215.120785</td>
          <td>47705.656999</td>
          <td>11575.587665</td>
          <td>8279.860420</td>
          <td>...</td>
          <td>10656.027801</td>
          <td>48767.120558</td>
          <td>10728.708037</td>
          <td>10033.975622</td>
          <td>15357.169798</td>
          <td>27571.712446</td>
          <td>2496.420415</td>
          <td>6544.212344</td>
          <td>7950.707690</td>
          <td>1371.656217</td>
        </tr>
        <tr>
          <th>3</th>
          <td>31540.014497</td>
          <td>32238.930716</td>
          <td>51695.470357</td>
          <td>19748.574190</td>
          <td>13231.266416</td>
          <td>10573.566154</td>
          <td>19723.733790</td>
          <td>71201.292388</td>
          <td>17276.710011</td>
          <td>12357.795695</td>
          <td>...</td>
          <td>15904.255362</td>
          <td>72785.540086</td>
          <td>16012.731529</td>
          <td>14975.834671</td>
          <td>22920.768853</td>
          <td>41151.127204</td>
          <td>3725.938832</td>
          <td>9767.319137</td>
          <td>11866.531110</td>
          <td>2047.214135</td>
        </tr>
        <tr>
          <th>4</th>
          <td>41843.790305</td>
          <td>42771.034765</td>
          <td>68583.811892</td>
          <td>26200.216151</td>
          <td>17553.775615</td>
          <td>14027.833911</td>
          <td>26167.260665</td>
          <td>94461.971424</td>
          <td>22920.821135</td>
          <td>16394.951618</td>
          <td>...</td>
          <td>21100.000648</td>
          <td>96563.775418</td>
          <td>21243.914785</td>
          <td>19868.275130</td>
          <td>30408.731918</td>
          <td>54594.747814</td>
          <td>4943.162064</td>
          <td>12958.194862</td>
          <td>15743.196296</td>
          <td>2716.016475</td>
        </tr>
      </tbody>
    </table>
    <p>5 rows × 1000 columns</p>
    </div>



Draw a static plot showing the results, and a marginal density plot
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The code below is what we might use for making a static graphic for
publication in a print environment. The result is an image, and we have
programmatic control over how we want the image displayed and saved.

In the lefthand side of the plot, we draw all traces from the suite of
simulation runs. Plotting each line in the same color, and setting a low
value for ``alpha``, the opacity of each line, we can see the regions of
the plot in which a large number of simulations agree on the values the
system will take, despite the parametric uncertainty.

In the righthand plot, we use a gaussian `Kernel Density
Estimator <https://seaborn.pydata.org/generated/seaborn.kdeplot.html>`__
provided by the `seaborn <https://seaborn.pydata.org/index.html>`__
library. The KDE gives an indication of the regions in which the density
of simulation results is highest, at a specific time point in the
simulation, which we refer to here as ``density_time``.

To indicate the simulation time for which we are displaying a density
estimate, we’ll add a vertical line at the point on the lefthand plot at
which the density curve is calculated.

.. code:: ipython3

    # define when to show the density
    density_time = 85
    
    # left side: plot all traces, slightly transparent
    plt.subplot2grid((1,4), loc=(0,0), colspan=3)
    [plt.plot(result.index, result[i], 'b', alpha=.02) for i in result.columns]
    ymax = result.max().max()
    plt.ylim(0, ymax)
    
    # left side: add marker of density location
    plt.vlines(density_time, 0, ymax, 'k')
    plt.text(density_time, ymax, 'Density Calculation', ha='right', va='top', rotation=90)
    
    # right side: gaussian KDE on selected timestamp
    plt.subplot2grid((1,4), loc=(0,3))
    seaborn.kdeplot(y=result.loc[density_time])
    plt.ylim(0, ymax)
    plt.yticks([])
    plt.xticks([])
    plt.xlabel('Density')
    
    plt.suptitle('Emissions scenarios under uncertainty', fontsize=16);



.. image:: plotting_suite_of_simulations_files/plotting_suite_of_simulations_10_0.png


Static density plot with selector
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

For purposes of lightweight exploration, we can add a slider to the
chart. In this case, whenever the slider is moved, the figure is
regenerated, making this a suitable method for exploring results before
feeding in to print graphics.

.. code:: ipython3

    import matplotlib as mpl
    from ipywidgets import interact, IntSlider
    sim_time = 200
    slider_time = IntSlider(description = 'Select Time for plotting Density',
                            min=0, max=result.index[-1], value=1)

.. code:: ipython3

    @interact(density_time=slider_time)
    def update(density_time): 
        ax1 = plt.subplot2grid((1,4), loc=(0,0), colspan=3)
        [ax1.plot(result.index, result[i], 'b', alpha=.02) for i in result.columns]
        ymax = result.max().max()
        ax1.set_ylim(0, ymax)
    
        # left side: add marker of density location
        ax1.vlines(density_time, 0, ymax, 'k')
        ax1.text(density_time, ymax, 'Density Calculation', ha='right', va='top', rotation=90)
    
        # right side: gaussian KDE on selected timestamp
        ax2 = plt.subplot2grid((1,4), loc=(0,3))
        seaborn.kdeplot(y=result.loc[density_time], ax=ax2)
        ax2.set_ylim(0, ymax)
        ax2.set_yticks([])
        ax2.set_xticks([])
        ax2.set_xlabel('Density')
    
        plt.suptitle('Emissions scenarios under uncertainty', fontsize=16);
        



.. parsed-literal::

    interactive(children=(IntSlider(value=1, description='Select Time for plotting Density'), Output()), _dom_clas…