{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# cadCAD Tutorials: The Robot and the Marbles, part 4\n",
    "In parts [1](../robot-marbles-part-1/robot-marbles-part-1.ipynb) and [2](../robot-marbles-part-2/robot-marbles-part-2.ipynb) we introduced the 'language' in which a system must be described in order for it to be interpretable by cadCAD and some of the basic concepts of the library:\n",
    "* State Variables\n",
    "* Timestep\n",
    "* State Update Functions\n",
    "* Partial State Update Blocks\n",
    "* Simulation Configuration Parameters\n",
    "* Policies\n",
    "\n",
    "In [part 3](../robot-marbles-part-3/robot-marbles-part-3.ipynb) we covered how to describe the presence of asynchronous subsystems within the system being modeled in cadCAD.\n",
    "\n",
    "So far, all the examples referred to deterministic systems: no matter how many times you ran one of those simulations, the results would be the same. However, systems are more commonly non-deterministic, and modelling them as deterministic might be an oversimplification sometimes. \n",
    "\n",
    "In this notebook, we'll cover cadCAD's support for modelling non-deterministic systems and Monte Carlo simulations. But first let's copy the base configuration with which we ended Part 3. Here's the description of that system:\n",
    "\n",
    "__The robot and the marbles__  \n",
    "* Picture a box (`box_A`) with ten marbles in it; an empty box (`box_B`) next to the first one; and __two__ robot arms capable of taking a marble from any one of the boxes and dropping it into the other one. \n",
    "* The robots are programmed to take one marble at a time from the box containing the largest number of marbles and drop it in the other box. They repeat that process until the boxes contain an equal number of marbles.\n",
    "* The robots act __asynchronously__; robot 1 acts once every two timesteps, and robot 2 acts once every three timesteps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "\n",
    "\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# List of all the state variables in the system and their initial values\n",
    "initial_conditions = {\n",
    "    'box_A': 10, # as per the description of the example, box_A starts out with 10 marbles in it\n",
    "    'box_B': 0 # as per the description of the example, box_B starts out empty\n",
    "}\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "\n",
    "\n",
    "\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# Settings of general simulation parameters, unrelated to the system itself\n",
    "# `T` is a range with the number of discrete units of time the simulation will run for;\n",
    "# `N` is the number of times the simulation will be run (Monte Carlo runs)\n",
    "# In this example, we'll run the simulation once (N=1) and its duration will be of 10 timesteps\n",
    "# We'll cover the `M` key in a future article. For now, let's leave it empty\n",
    "simulation_parameters = {\n",
    "    'T': range(10),\n",
    "    'N': 1,\n",
    "    'M': {}\n",
    "}\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# We specify the robot arm's logic in a Policy Function\n",
    "def robot_arm(params, step, sL, s):\n",
    "    add_to_A = 0\n",
    "    if (s['box_A'] > s['box_B']):\n",
    "        add_to_A = -1\n",
    "    elif (s['box_A'] < s['box_B']):\n",
    "        add_to_A = 1\n",
    "    return({'add_to_A': add_to_A, 'add_to_B': -add_to_A})\n",
    "    \n",
    "robots_periods = [2,3] # Robot 1 acts once every 2 timesteps; Robot 2 acts once every 3 timesteps\n",
    "\n",
    "def robot_arm_1(params, step, sL, s):\n",
    "    _robotId = 1\n",
    "    if s['timestep']%robots_periods[_robotId-1]==0: # on timesteps that are multiple of 2, Robot 1 acts\n",
    "        return robot_arm(params, step, sL, s)\n",
    "    else:\n",
    "        return({'add_to_A': 0, 'add_to_B': 0}) # for all other timesteps, Robot 1 doesn't interfere with the system\n",
    "\n",
    "def robot_arm_2(params, step, sL, s):\n",
    "    _robotId = 2\n",
    "    if s['timestep']%robots_periods[_robotId-1]==0: # on timesteps that are multiple of 3, Robot 2 acts\n",
    "        return robot_arm(params, step, sL, s)\n",
    "    else:\n",
    "        return({'add_to_A': 0, 'add_to_B': 0}) # for all other timesteps, Robot 2 doesn't interfere with the system\n",
    "\n",
    "    \n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# We make the state update functions less \"intelligent\",\n",
    "# ie. they simply add the number of marbles specified in _input \n",
    "# (which, per the policy function definition, may be negative)\n",
    "def increment_A(params, step, sL, s, _input):\n",
    "    y = 'box_A'\n",
    "    x = s['box_A'] + _input['add_to_A']\n",
    "    return (y, x)\n",
    "\n",
    "def increment_B(params, step, sL, s, _input):\n",
    "    y = 'box_B'\n",
    "    x = s['box_B'] + _input['add_to_B']\n",
    "    return (y, x)\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "\n",
    "\n",
    "\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# In the Partial State Update Blocks, \n",
    "# the user specifies if state update functions will be run in series or in parallel\n",
    "# and the policy functions that will be evaluated in that block\n",
    "partial_state_update_blocks = [\n",
    "    { \n",
    "        'policies': { # The following policy functions will be evaluated and their returns will be passed to the state update functions\n",
    "            'robot_arm_1': robot_arm_1,\n",
    "            'robot_arm_2': robot_arm_2\n",
    "        },\n",
    "        'variables': { # The following state variables will be updated simultaneously\n",
    "            'box_A': increment_A,\n",
    "            'box_B': increment_B\n",
    "        }\n",
    "    }\n",
    "]\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "\n",
    "\n",
    "\n",
    "from cadCAD.configuration import Configuration\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# The configurations above are then packaged into a `Configuration` object\n",
    "config = Configuration(initial_state=initial_conditions, #dict containing variable names and initial values\n",
    "                       partial_state_update_blocks=partial_state_update_blocks, #dict containing state update functions\n",
    "                       sim_config=simulation_parameters #dict containing simulation parameters\n",
    "                      )\n",
    "\n",
    "from cadCAD.engine import ExecutionMode, ExecutionContext, Executor\n",
    "exec_mode = ExecutionMode()\n",
    "exec_context = ExecutionContext(exec_mode.single_proc)\n",
    "executor = Executor(exec_context, [config]) # Pass the configuration object inside an array\n",
    "raw_result, tensor = executor.execute() # The `execute()` method returns a tuple; its first elements contains the raw results\n",
    "\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "df = pd.DataFrame(raw_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14676 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot('timestep', ['box_A', 'box_B'], grid=True, \n",
    "        xticks=list(df['timestep'].drop_duplicates()), \n",
    "        colormap = 'RdYlGn',\n",
    "        yticks=list(range(1+(df['box_A']+df['box_B']).max())));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Non-determinism\n",
    "Non-deterministic systems exhibit different behaviors on different runs for the same input. The order of heads and tails in a series of 3 coin tosses, for example, is non deterministic. \n",
    "\n",
    "Our robots and marbles system is currently modelled as a deterministic system. Meaning that every time we run the simulation: none of the robots act on timestep 1; robot 1 acts on timestep 2; robot 2 acts on timestep 3; an so on. \n",
    "\n",
    "If however we were to define that at every timestep each robot would act with a probability P, then we would have a non-deterministic (probabilistic) system. Let's make the following changes to our system.\n",
    "* Robot 1: instead of acting once every two timesteps, there's a 50% chance it will act in any given timestep\n",
    "* Robot 2: instead of acting once every three timesteps, there's a 33.33% chance it will act in any given timestep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "from numpy.random import rand\n",
    "\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# We specify each of the robots logic in a Policy Function\n",
    "robots_probabilities = [0.5,1/3] # Robot 1 acts with a 50% probability; Robot 2, 33.33%\n",
    "\n",
    "def robot_arm_1(params, step, sL, s):\n",
    "    _robotId = 1\n",
    "    if rand()<robots_probabilities[_robotId-1]: # draw a random number between 0 and 1; if it's smaller than the robot's parameter, it acts\n",
    "        return robot_arm(params, step, sL, s)\n",
    "    else:\n",
    "        return({'add_to_A': 0, 'add_to_B': 0}) # otherwise, the robot doesn't interfere with the system\n",
    "\n",
    "def robot_arm_2(params, step, sL, s):\n",
    "    _robotId = 2\n",
    "    if rand()<robots_probabilities[_robotId-1]: # draw a random number between 0 and 1; if it's smaller than the robot's parameter, it acts\n",
    "        return robot_arm(params, step, sL, s)\n",
    "    else:\n",
    "        return({'add_to_A': 0, 'add_to_B': 0}) # otherwise, the robot doesn't interfere with the system\n",
    "\n",
    "\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# In the Partial State Update Blocks, \n",
    "# the user specifies if state update functions will be run in series or in parallel\n",
    "# and the policy functions that will be evaluated in that block\n",
    "partial_state_update_blocks = [\n",
    "    { \n",
    "        'policies': { # The following policy functions will be evaluated and their returns will be passed to the state update functions\n",
    "            'robot_arm_1': robot_arm_1,\n",
    "            'robot_arm_2': robot_arm_2\n",
    "        },\n",
    "        'variables': { # The following state variables will be updated simultaneously\n",
    "            'box_A': increment_A,\n",
    "            'box_B': increment_B\n",
    "        }\n",
    "    }\n",
    "]\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run the simulation with those configurations, the system is unlikely to exhibit the same behavior as it did in its deterministic version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# The configurations above are then packaged into a `Configuration` object\n",
    "config = Configuration(initial_state=initial_conditions, #dict containing variable names and initial values\n",
    "                       partial_state_update_blocks=partial_state_update_blocks, #dict containing state update functions\n",
    "                       sim_config=simulation_parameters #dict containing simulation parameters\n",
    "                      )\n",
    "\n",
    "exec_mode = ExecutionMode()\n",
    "exec_context = ExecutionContext(exec_mode.single_proc)\n",
    "executor = Executor(exec_context, [config]) # Pass the configuration object inside an array\n",
    "raw_result, tensor = executor.execute() # The `execute()` method returns a tuple; its first elements contains the raw results\n",
    "\n",
    "df = pd.DataFrame(raw_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15456 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot('timestep', ['box_A', 'box_B'], grid=True, \n",
    "        xticks=list(df['timestep'].drop_duplicates()), \n",
    "        colormap = 'RdYlGn',\n",
    "        yticks=list(range(1+(df['box_A']+df['box_B']).max())));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And if we run it again, it returns yet another result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# The configurations above are then packaged into a `Configuration` object\n",
    "config = Configuration(initial_state=initial_conditions, #dict containing variable names and initial values\n",
    "                       partial_state_update_blocks=partial_state_update_blocks, #dict containing state update functions\n",
    "                       sim_config=simulation_parameters #dict containing simulation parameters\n",
    "                      )\n",
    "\n",
    "exec_mode = ExecutionMode()\n",
    "exec_context = ExecutionContext(exec_mode.single_proc)\n",
    "executor = Executor(exec_context, [config]) # Pass the configuration object inside an array\n",
    "raw_result, tensor = executor.execute() # The `execute()` method returns a tuple; its first elements contains the raw results\n",
    "\n",
    "df = pd.DataFrame(raw_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 15616 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot('timestep', ['box_A', 'box_B'], grid=True, \n",
    "        xticks=list(df['timestep'].drop_duplicates()), \n",
    "        colormap = 'RdYlGn',\n",
    "        yticks=list(range(1+(df['box_A']+df['box_B']).max())));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to take advantage of cadCAD's Monte Carlo simulation features, we should modify the configuration file so as to define the number of times we want the same simulation to be run. This is done in the `N` key of the `simulation_parameters` dict."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# Settings of general simulation parameters, unrelated to the system itself\n",
    "# `T` is a range with the number of discrete units of time the simulation will run for;\n",
    "# `N` is the number of times the simulation will be run (Monte Carlo runs)\n",
    "# In this example, we'll run the simulation once (N=1) and its duration will be of 10 timesteps\n",
    "# We'll cover the `M` key in a future article. For now, let's leave it empty\n",
    "simulation_parameters = {\n",
    "    'T': range(10),\n",
    "    'N': 50, # We'll run the same simulation 50 times; the random events in each simulation are independent\n",
    "    'M': {}\n",
    "}\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "\n",
    "# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # \n",
    "# The configurations above are then packaged into a `Configuration` object\n",
    "config = Configuration(initial_state=initial_conditions, #dict containing variable names and initial values\n",
    "                       partial_state_update_blocks=partial_state_update_blocks, #dict containing state update functions\n",
    "                       sim_config=simulation_parameters #dict containing simulation parameters\n",
    "                      )\n",
    "\n",
    "exec_mode = ExecutionMode()\n",
    "exec_context = ExecutionContext(exec_mode.single_proc)\n",
    "executor = Executor(exec_context, [config]) # Pass the configuration object inside an array\n",
    "raw_result, tensor = executor.execute() # The `execute()` method returns a tuple; its first elements contains the raw results\n",
    "\n",
    "df = pd.DataFrame(raw_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>box_A</th>\n",
       "      <th>box_B</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>run</th>\n",
       "      <th>timestep</th>\n",
       "      <th>substep</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">1</th>\n",
       "      <th>0</th>\n",
       "      <th>0</th>\n",
       "      <td>10</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <th>1</th>\n",
       "      <td>9</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <th>1</th>\n",
       "      <td>8</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <th>1</th>\n",
       "      <td>8</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <th>1</th>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <th>...</th>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"5\" valign=\"top\">50</th>\n",
       "      <th>6</th>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <th>1</th>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>550 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                      box_A  box_B\n",
       "run timestep substep              \n",
       "1   0        0           10      0\n",
       "    1        1            9      1\n",
       "    2        1            8      2\n",
       "    3        1            8      2\n",
       "    4        1            6      4\n",
       "...                     ...    ...\n",
       "50  6        1            5      5\n",
       "    7        1            5      5\n",
       "    8        1            5      5\n",
       "    9        1            5      5\n",
       "    10       1            5      5\n",
       "\n",
       "[550 rows x 2 columns]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display\n",
    "tmp_rows = pd.options.display.max_rows\n",
    "pd.options.display.max_rows = 10\n",
    "display(df.set_index(['run', 'timestep', 'substep']))\n",
    "pd.options.display.max_rows = tmp_rows"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting two of those runs allows us to see the different behaviors over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 14916 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 14760 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df[df['run']==1].plot('timestep', ['box_A', 'box_B'], grid=True,\n",
    "                          xticks=list(df['timestep'].drop_duplicates()), \n",
    "                          yticks=list(range(11)),\n",
    "                          colormap = 'RdYlGn');\n",
    "df[df['run']==9].plot('timestep', ['box_A', 'box_B'], grid=True,\n",
    "                            xticks=list(df['timestep'].drop_duplicates()), \n",
    "                            yticks=list(range(11)),\n",
    "                            colormap = 'RdYlGn');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we plot all those runs onto a single chart, we can see every possible trajectory for the number of marbles in each box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 69484 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = None\n",
    "for i in range(simulation_parameters['N']):\n",
    "    ax = df[df['run']==i+1].plot('timestep', ['box_A', 'box_B'],\n",
    "                                     grid=True,\n",
    "                                     xticks=list(df['timestep'].drop_duplicates()), \n",
    "                                     yticks=list(range(1+max(df['box_A'].max(),df['box_B'].max()))),\n",
    "                                     legend = (ax == None),\n",
    "                                     colormap = 'RdYlGn',\n",
    "                                     ax = ax\n",
    "                                  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For some analyses, it might make sense to look at the data in aggregate. Take the median for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x108868b00>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 15316 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dfmc_median = df.groupby(['timestep', 'substep']).median().reset_index()\n",
    "dfmc_median.plot('timestep', ['box_A', 'box_B'], \n",
    "                 grid=True,\n",
    "                 xticks=list(dfmc_median['timestep'].drop_duplicates()), \n",
    "                 yticks=list(range(int(1+max(dfmc_median['box_A'].max(),dfmc_median['box_B'].max())))),\n",
    "                 colormap = 'RdYlGn'\n",
    "                )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or look at edge cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 13600 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "max_final_A = df[df['timestep']==df['timestep'].max()]['box_A'].max()\n",
    "# max_final_A\n",
    "slow_runs = df[(df['timestep']==df['timestep'].max()) & \n",
    "         (df['box_A']==max_final_A)]['run']\n",
    "slow_runs = list(slow_runs)\n",
    "slow_runs\n",
    "\n",
    "ax = None\n",
    "for i in slow_runs:\n",
    "    ax = df[df['run']==i].plot('timestep', ['box_A', 'box_B'],\n",
    "                                     grid=True,\n",
    "                                     xticks=list(df['timestep'].drop_duplicates()), \n",
    "                                     yticks=list(range(1+max(df['box_A'].max(),df['box_B'].max()))),\n",
    "                                     legend = (ax == None),\n",
    "                                     colormap = 'RdYlGn',\n",
    "                                     ax = ax\n",
    "                                  )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We invite the reader to fork this code and come up with answers for their other questions that might be interesting to look at. For example:\n",
    "* How often does box B momentarily contain more marbles than box A?\n",
    "* What's the frequency distribution of the time to reach equilibrium?\n",
    "* What's the probability distribution of the waiting times of each one of the robots?"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
