4-sa_practice_wolf_sheep.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Practice SA on the Wolf-Sheep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the second practical exercise for the participants of the workshop on *sensitivity analysis* (SA) at the Social Simulation Festival 2021.\n",
    "\n",
    "The task remains the same: is to perform a more comperhansive SA than we did previously, but the model is different - [Wolf-Sheep](https://ccl.northwestern.edu/netlogo/models/WolfSheepPredation). Again, explore more parameters, their ranges and add more scenarios. We kept the structure of the notebook as previous, so you can always read something up. Let's go!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 0. Installations and imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fatal: destination path 'SSF2021' already exists and is not an empty directory.\r\n"
     ]
    }
   ],
   "source": [
    "# Clone the repo to make its file available for Google Colab\n",
    "!git clone https://github.com/BROSE-Uninc/SSF2021.git"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "^C\r\n"
     ]
    }
   ],
   "source": [
    "# Install necessary packages\n",
    "!pip install ema_workbench mesa ipyparallel SALib &> /dev/null"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import necessary packages\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import datetime\n",
    "import random\n",
    "\n",
    "# Import EMA Workbench modules\n",
    "from ema_workbench import ReplicatorModel, RealParameter, BooleanParameter, IntegerParameter, Constant, TimeSeriesOutcome, perform_experiments, save_results, ema_logging\n",
    "\n",
    "# Initialize logger to keep track of experiments run\n",
    "ema_logging.log_to_stderr(ema_logging.INFO)\n",
    "\n",
    "# Import Mesa Wolf-Sheep model\n",
    "from SSF2021.wolf_sheep import model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Load the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The very first step of SA with EMA Workbench is to define or \"load\" the model as a function. That is, EMA Workbench treats all models as functions (read *black box*). They are supposed to have **inputs** (parameters, constants, uncertainties and policy levers) and **outputs** (outcomes, KPIs). The model structure is not interesting for EMA Workbench. It may be something simple as `def func(x)` which just returns x + 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting up the wolf-sheep model simulation as a function\n",
    "def model_wolf_sheep(height=20,\n",
    "                     width=20,\n",
    "                     initial_sheep=100,\n",
    "                     initial_wolves=50,\n",
    "                     sheep_reproduce=0.04,\n",
    "                     wolf_reproduce=0.05,\n",
    "                     wolf_gain_from_food=20,\n",
    "                     grass=False,\n",
    "                     grass_regrowth_time=30,\n",
    "                     sheep_gain_from_food=4,\n",
    "                     steps=200):\n",
    "    \n",
    "    from SSF2021.wolf_sheep import model\n",
    "    \n",
    "    # Initialising the model\n",
    "    wolf_sheep = model.WolfSheep(height=20,\n",
    "                                   width=20,\n",
    "                                   initial_sheep=100,\n",
    "                                   initial_wolves=50,\n",
    "                                   sheep_reproduce=0.04,\n",
    "                                   wolf_reproduce=0.05,\n",
    "                                   wolf_gain_from_food=20,\n",
    "                                   grass=False,\n",
    "                                   grass_regrowth_time=30,\n",
    "                                   sheep_gain_from_food=4)\n",
    "                \n",
    "    # Run the model steps times\n",
    "    wolf_sheep.run_model(steps)\n",
    "    \n",
    "    # Get model outcomes\n",
    "    outcomes = wolf_sheep.datacollector.get_model_vars_dataframe()\n",
    "    \n",
    "    # Return model outcomes\n",
    "    # below to be changed!\n",
    "    return {'TIME' : list(range(steps + 1)),\n",
    "            \"Wolves\" : outcomes[\"Wolves\"].tolist(),\n",
    "            \"Sheep\" : outcomes[\"Sheep\"].tolist()}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's parameterize and test out our Mesa model. What is supposed to happen? First, we shouldn't get any error 😅. Second, after we run `model_wolf_sheep` function it has to give us a set of model outcomes. Let's try."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'TIME': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10],\n",
       " 'Wolves': [50, 52, 53, 55, 58, 54, 53, 57, 60, 63, 62],\n",
       " 'Sheep': [100, 96, 93, 89, 82, 79, 76, 68, 63, 60, 54]}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Parametrize the model\n",
    "height=20\n",
    "width=20\n",
    "initial_sheep=100\n",
    "initial_wolves=50\n",
    "sheep_reproduce=0.04\n",
    "wolf_reproduce=0.05\n",
    "wolf_gain_from_food=20\n",
    "grass=False\n",
    "grass_regrowth_time=30\n",
    "sheep_gain_from_food=4\n",
    "steps=10\n",
    "\n",
    "model_wolf_sheep(height=height, \n",
    "                 width=width, \n",
    "                 initial_sheep=initial_sheep, \n",
    "                 initial_wolves=initial_wolves, \n",
    "                 sheep_reproduce=sheep_reproduce, \n",
    "                 wolf_reproduce=wolf_reproduce, \n",
    "                 wolf_gain_from_food=wolf_gain_from_food,\n",
    "                 grass=grass,\n",
    "                 grass_regrowth_time=grass_regrowth_time,\n",
    "                 sheep_gain_from_food=sheep_gain_from_food,\n",
    "                 steps=steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Design experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's time to design experiments. What does it mean? Well, we have to specify:\n",
    "\n",
    "* **which model parameters** aka *inputs* are we going to sample, what are their **ranges**, and random **distributions**,\n",
    "* what we will keep as **constants** and do not change over the model run,\n",
    "* and finally which **outcomes** we want to observe.\n",
    "\n",
    "It's an important step in SA workflow and we have to be careful. Because if parameter ranges are too narrow or they're sampled from e.g. a Normal distribution, there is a chance that you'll overlook import model behavior. This is why model parameters are named **uncertainties** in the EMA Workbench. We often do not know parameter vales and how to explore many plausible options."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's talk about \"tech\" part. First we have to initialize an instance of EMA Workbench called `ReplicatorModel`. This is how we \"connect\" EMA Workbench to our Python model. We have to pass a name of our model to `ReplicatorModel`, and also pass the function that we defined previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate and pass the model \n",
    "model = ReplicatorModel(..., function=...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define model parameters and their ranges to be sampled\n",
    "model.uncertainties = [RealParameter(..., ..., ...)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define model parameters that will remain constant\n",
    "model.constants = [Constant(\"parameter_A\", parameter_A)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define model outcomes\n",
    "model.outcomes = [TimeSeriesOutcome('TIME'),\n",
    "                  TimeSeriesOutcome('Wolves'),\n",
    "                  TimeSeriesOutcome('Sheep')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the number of replications\n",
    "model.replications = 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Run the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run experiments with the aforementioned parameters and outputs\n",
    "results = perform_experiments(model, scenarios=100, uncertainty_sampling='sobol')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the results\n",
    "experiments, outcomes = results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiments.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outcomes.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ema_workbench.util.utilities import save_results, load_results\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creaet a directory to store the results\n",
    "directory = 'results/wolf_sheep'\n",
    "if not os.path.exists(directory):\n",
    "    os.makedirs(directory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the results\n",
    "save_results(results, 'results/wolf_sheep/results.tar.gz')a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the results\n",
    "results = load_results('results/wolf_sheep/results.tar.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiments, outcomes = results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A bit of preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(random.choice(list(outcomes)))\n",
    "outcomes[random.choice(list(outcomes))].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_outcomes = {key:np.mean(outcomes[key],axis=1) for key in outcomes.keys()}\n",
    "mean_results = (experiments.copy(), mean_outcomes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now the shape of this array doesn't have 10 in it  \n",
    "mean_outcomes[random.choice(list(outcomes))].shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visuals!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ema_workbench.analysis.plotting import lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting all of the results\n",
    "plt.rcParams['figure.figsize'] = [10, 12]\n",
    "\n",
    "figure = lines(experiments, outcomes_2D) #show lines, and end state density\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Sensitivity analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from SALib.analyze import sobol\n",
    "from ema_workbench.em_framework.salib_samplers import get_SALib_problem\n",
    "from SSF2021.src.plot import plot_sobol_indices\n",
    "sns.set_style('white')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the problem\n",
    "problem = get_SALib_problem(model.uncertainties)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select and normalize an outcome\n",
    "normalized_outcomes = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform Sobol SA\n",
    "Si = sobol.analyze(problem=problem, Y=normalized_resistant,\n",
    "                   calc_second_order=True, print_to_console=False)\n",
    "\n",
    "# Get scores by type \n",
    "Si_filter = {k:Si[k] for k in ['ST', 'ST_conf', 'S1', 'S1_conf']}\n",
    "\n",
    "# Create a DataFrame out of them\n",
    "Si_df = pd.DataFrame(Si_filter, index=problem['names'])\n",
    "\n",
    "# Get indices and error bars\n",
    "indices = Si_df[['S1','ST']]\n",
    "err = Si_df[['S1_conf','ST_conf']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the results\n",
    "fig, ax = plt.subplots(1)\n",
    "indices.plot.bar(yerr=err.values.T,ax=ax)\n",
    "fig.set_size_inches(8,6)\n",
    "fig.subplots_adjust(bottom=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set_style('whitegrid')\n",
    "fig = plot_sobol_indices(Si, problem, criterion='ST', threshold=0.005)\n",
    "fig.set_size_inches(7,7)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}