simpy-getting-started.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Over the last 30 or so years, I've built numerous healthcare discrete event simulation (DES) models. I built models as a student, as an analyst for several healthcare systems and consulting firms, and as a teacher and researcher at Oakland University. I've built models of inpatient units, pharmacies, pneumatic tube systems, outpatient clinics, surgical and recovery suites, emergency rooms and a whole bunch more. I've built models in several different commercial discrete event simulation languages/packages such as SIMAN, Arena, MedModel, and Simio. While each of these are very fine products, one thing that's always bothered me was how difficult to it was to share, review, and improve upon the models built by the healthcare simulation modeling community. These packages are really expensive often have pretty restrictive licensing. So, academics and practitioners publish articles describing how a simulation model of some healthcare related system (e.g. inpatient beds, pharmacy, lab, surgical suites...) was built and used in some project. Great. Inspiring. Not super helpful for the person looking to build a similar model, or better yet, build upon that same model. Models were rarely if ever made available and even if they were, might be written in a package that you don't have, don't know, or can't afford.\n",
    "\n",
    "As the years went on, I found myself more and more attracted to the world of open source software. I fell in love with Linux and R and Python. However, open source discrete event simulation options were not plentiful. I did have success building a pretty big pneumatic tube simulation model in Java using the Simkit library developed by [Arnie Buss and others at the Naval Postgraduate School](https://www.movesinstitute.org/research/simulation-modeling-for-analysis/). While I really like Simkit, I'm not a huge fan of programming in Java. Most of my research work and teaching is focused on R and Python and ideally would love to have a Python based DES framework. So, I was pretty happy when I stumbled onto the SimPy project and followed it for a bit as it moved through its early phases of growth and change. With SimPy 3.x now available, I decided to take the plunge and try to build some DES models for an ongoing research project involving simulating a simplified inpatient obstetrical patient flow network. I originally built the model in [Simio](http://hselab.org/comparing-predictive-models-for-obstetrical-unit-occupancy-using-caret-part-1.html). I describe the project in a bit of detail in a [previous post](http://hselab.org/comparing-predictive-models-for-obstetrical-unit-occupancy-using-caret-part-1.html) on using the R package `caret` for building metamodels of simulation output.\n",
    "\n",
    "I decided to start by building SimPy models of overly simplified versions of this patient flow network. As my SimPy knowledge grew, I'd add more features to the model and overhaul its architecture. My goal is be able to build models of a similar complexity to the Java based pneumatic tube model. I'll also explore the ecosystem of Python tools that might be used alongside SimPy for things like model animation, logging and monitoring of simulation runs, simulation experiment management and management of simulation input and output data.\n",
    "\n",
    "If you aren't familiar with [SimPy](https://simpy.readthedocs.io/en/latest/), visit the main site and check out the documentation. It's got an overview, installation instructions, a basic tutorial, topic guide, example models and the API. Another site with some very nice network models and tutorials is [available from Grotto Networking](https://www.grotto-networking.com/DiscreteEventPython.html#Intro). Models in SimPy can range from purely process based with no object oriented (OO) features to almost fully OO based models.\n",
    "\n",
    "**Note:** This series of blog posts assumes some knowledge of DES modeling, Python, and SimPy. I'm not trying to teach you how to use SimPy in general, but instead to share my journey of building progressively move complex DES patient flow models in SimPy. In the end, I hope to end up with a nice set of models and accompanying explanatory text that others can use and build on. I highly recommend the [SimPy in 10 Minutes](https://simpy.readthedocs.io/en/latest/simpy_intro/basic_concepts.html) introduction in the official docs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OB patient flow system to model\n",
    "Here's a picture of the simple system I'll be modeling.\n",
    "\n",
    "![](images/obflow.png)\n",
    "\n",
    "For this first set of simple models, features/assumptions of the system include:\n",
    "\n",
    "* Patients in labor arrive according to a poisson process\n",
    "* Assume c-section rate is zero and that all patients flow from OBS --> LDR --> PP units.\n",
    "* Length of stay at each unit is exponentially distributed with a given unit specific mean\n",
    "* Each unit has capacity level (i.e. the number of beds)\n",
    "* Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.\n",
    "\n",
    "I'm going to start by building the simplest model possible and then start layering on complexity and features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## obflow_1: A Hello World level SimPy patient flow model\n",
    "\n",
    "For my first ever SimPy model I decided to model the system as a sequence of delays. In other words, even though SimPy has features for modeling capacitated resources, I'm going to ignore them for now. This model is entirely *process oriented*. SimPy models processes with generator functions that can yield for events such as time durations or some system state condition. No classes are created in this first model. You'll also note that all the functions take\n",
    "a *simulation environment* instance as their first argurment. This is what you do in SimPy. It handles the clock and next event list logic as well as providing the means for running or single stepping through your model.\n",
    "\n",
    "The model prints out a message as each patient enters each delay stage. No summary stats are computed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Patient1 entering OBS at 0.0000\n",
      "Patient1 entering LDR at 0.3475\n",
      "Patient2 entering OBS at 3.7668\n",
      "Patient3 entering OBS at 5.6439\n",
      "Patient1 entering PP at 6.1449\n",
      "Patient4 entering OBS at 7.2625\n",
      "Patient4 entering LDR at 7.3175\n",
      "Patient4 entering PP at 7.8836\n",
      "Patient5 entering OBS at 8.9829\n",
      "Patient3 entering LDR at 10.8140\n",
      "Patient2 entering LDR at 12.9060\n",
      "Patient6 entering OBS at 16.8153\n",
      "Patient5 entering LDR at 16.9444\n",
      "Patient7 entering OBS at 17.8737\n",
      "Patient7 entering LDR at 18.0987\n",
      "Patient7 entering PP at 18.1332\n",
      "Patient8 entering OBS at 18.7157\n",
      "Patient6 entering LDR at 19.4866\n",
      "Patient9 entering OBS at 20.4913\n",
      "Patient8 entering LDR at 22.8456\n",
      "Patient8 entering PP at 23.1882\n",
      "Patient5 entering PP at 24.5646\n"
     ]
    }
   ],
   "source": [
    "import simpy\n",
    "import numpy as np\n",
    "from numpy.random import RandomState\n",
    "\n",
    "\"\"\"\n",
    "Simple OB patient flow model - NOT OO\n",
    "\n",
    "Details:\n",
    "\n",
    "- Generate arrivals via Poisson process\n",
    "- Simple delays modeling stay in OBS, LDR, and PP units. No capacitated resource contention modeled.\n",
    "- Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "# Arrival rate and length of stay inputs.\n",
    "ARR_RATE = 0.4\n",
    "MEAN_LOS_OBS = 3\n",
    "MEAN_LOS_LDR = 12\n",
    "MEAN_LOS_PP = 48\n",
    "\n",
    "RNG_SEED = 6353\n",
    "\n",
    "def source(env, arr_rate, prng=RandomState(0)):\n",
    "    \"\"\"Source generates patients according to a simple Poisson process\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        env : simpy.Environment\n",
    "            the simulation environment\n",
    "        arr_rate : float\n",
    "            exponential arrival rate\n",
    "        prng : RandomState object\n",
    "            Seeded RandomState object for generating pseudo-random numbers.\n",
    "            See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    patients_created = 0\n",
    "\n",
    "    # Infinite loop for generatirng patients according to a poisson process.\n",
    "    while True:\n",
    "\n",
    "        # Generate next interarrival time\n",
    "        iat = prng.exponential(1.0 / arr_rate)\n",
    "\n",
    "        # Generate length of stay in each unit for this patient\n",
    "        los_obs = prng.exponential(MEAN_LOS_OBS)\n",
    "        los_ldr = prng.exponential(MEAN_LOS_LDR)\n",
    "        los_pp = prng.exponential(MEAN_LOS_PP)\n",
    "\n",
    "        # Update counter of patients\n",
    "        patients_created += 1\n",
    "\n",
    "        # Create a new patient flow process.\n",
    "        obp = obpatient_flow(env, 'Patient{}'.format(patients_created),\n",
    "                             los_obs=los_obs, los_ldr=los_ldr, los_pp=los_pp)\n",
    "\n",
    "        # Register the process with the simulation environment\n",
    "        env.process(obp)\n",
    "\n",
    "        # This process will now yield to a 'timeout' event. This process will resume after iat time units.\n",
    "        yield env.timeout(iat)\n",
    "\n",
    "# Define an obpatient_flow \"process\"\n",
    "def obpatient_flow(env, name, los_obs, los_ldr, los_pp):\n",
    "    \"\"\"Process function modeling how a patient flows through system.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        env : simpy.Environment\n",
    "            the simulation environment\n",
    "        name : str\n",
    "            process instance id\n",
    "        los_obs : float\n",
    "            length of stay in OBS unit\n",
    "        los_ldr : float\n",
    "            length of stay in LDR unit\n",
    "        los_pp : float\n",
    "            length of stay in PP unit\n",
    "    \"\"\"\n",
    "\n",
    "    # Note how we are simply modeling each stay as a delay. There\n",
    "    # is NO contention for any resources.\n",
    "    print(\"{} entering OBS at {:.4f}\".format(name, env.now))\n",
    "    yield env.timeout(los_obs)\n",
    "\n",
    "    print(\"{} entering LDR at {:.4f}\".format(name, env.now))\n",
    "    yield env.timeout(los_ldr)\n",
    "\n",
    "    print(\"{} entering PP at {:.4f}\".format(name, env.now))\n",
    "    yield env.timeout(los_pp)\n",
    "\n",
    "\n",
    "# Initialize a simulation environment\n",
    "env = simpy.Environment()\n",
    "\n",
    "# Initialize a random number generator.\n",
    "# See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html\n",
    "prng = RandomState(RNG_SEED)\n",
    "\n",
    "# Create a process generator and start it and add it to the env\n",
    "# Calling obpatient(env) creates the generator.\n",
    "# env.process() starts and adds it to env\n",
    "runtime = 25\n",
    "env.process(source(env, ARR_RATE, prng))\n",
    "\n",
    "# Run the simulation\n",
    "env.run(until=runtime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This first simple model illustrates the basic structure of a SimPy model with a single process function. We see how the simulation environment is needed for orchestrating the process as it yields and resumes. Now let's move on to a slightly more complex model in which the OBS, LDR, and PP units are modeled as ``Resource`` objects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## obflow_2: OBS unit modeled as a Resource\n",
    "\n",
    "To start, let's model just the OBS unit as a ``Resource`` - LDR and PP units will still be modeled with delays. By doing this we can make sure we get the ``request`` and ``release`` logic working for one unit before duplicating it for the other two units. Like ``obflow_1``, this model is process oriented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Patient1 trying to get OBS at 0.0000\n",
      "Patient1 entering OBS at 0.0000\n",
      "Patient1 leaving OBS at 0.3475\n",
      "Patient1 entering LDR at 0.3475\n",
      "Patient2 trying to get OBS at 3.7668\n",
      "Patient2 entering OBS at 3.7668\n",
      "Patient3 trying to get OBS at 5.6439\n",
      "Patient3 entering OBS at 5.6439\n",
      "Patient1 entering PP at 6.1449\n",
      "Patient4 trying to get OBS at 7.2625\n",
      "Patient5 trying to get OBS at 8.9829\n",
      "Patient3 leaving OBS at 10.8140\n",
      "Patient3 entering LDR at 10.8140\n",
      "Patient4 entering OBS at 10.8140\n",
      "Patient4 waited 3.5515 time units for OBS bed\n",
      "Patient4 leaving OBS at 10.8690\n",
      "Patient4 entering LDR at 10.8690\n",
      "Patient5 entering OBS at 10.8690\n",
      "Patient5 waited 1.8861 time units for OBS bed\n",
      "Patient4 entering PP at 11.4351\n",
      "Patient2 leaving OBS at 12.9060\n",
      "Patient2 entering LDR at 12.9060\n",
      "Patient6 trying to get OBS at 16.8153\n",
      "Patient6 entering OBS at 16.8153\n",
      "Patient7 trying to get OBS at 17.8737\n",
      "Patient8 trying to get OBS at 18.7157\n",
      "Patient5 leaving OBS at 18.8306\n",
      "Patient5 entering LDR at 18.8306\n",
      "Patient7 entering OBS at 18.8306\n",
      "Patient7 waited 0.9569 time units for OBS bed\n",
      "Patient7 leaving OBS at 19.0556\n",
      "Patient7 entering LDR at 19.0556\n",
      "Patient8 entering OBS at 19.0556\n",
      "Patient8 waited 0.3400 time units for OBS bed\n",
      "Patient7 entering PP at 19.0901\n",
      "Patient6 leaving OBS at 19.4866\n",
      "Patient6 entering LDR at 19.4866\n",
      "Patient9 trying to get OBS at 20.4913\n",
      "Patient9 entering OBS at 20.4913\n",
      "Patient8 leaving OBS at 23.1856\n",
      "Patient8 entering LDR at 23.1856\n",
      "Patient8 entering PP at 23.5282\n",
      "Patient2 entering PP at 25.7847\n",
      "Patient5 entering PP at 26.4507\n",
      "Patient10 trying to get OBS at 26.8687\n",
      "Patient10 entering OBS at 26.8687\n",
      "Patient10 leaving OBS at 27.6307\n",
      "Patient10 entering LDR at 27.6307\n",
      "Patient11 trying to get OBS at 27.6876\n",
      "Patient11 entering OBS at 27.6876\n",
      "Patient6 entering PP at 28.3114\n",
      "Patient9 leaving OBS at 30.4208\n",
      "Patient9 entering LDR at 30.4208\n",
      "Patient11 leaving OBS at 30.4652\n",
      "Patient11 entering LDR at 30.4652\n",
      "Patient11 entering PP at 30.5128\n",
      "Patient12 trying to get OBS at 31.4872\n",
      "Patient12 entering OBS at 31.4872\n",
      "Patient13 trying to get OBS at 31.9981\n",
      "Patient13 entering OBS at 31.9981\n",
      "Patient12 leaving OBS at 32.2156\n",
      "Patient12 entering LDR at 32.2156\n",
      "Patient13 leaving OBS at 32.5002\n",
      "Patient13 entering LDR at 32.5002\n",
      "Patient14 trying to get OBS at 33.5398\n",
      "Patient14 entering OBS at 33.5398\n",
      "Patient12 entering PP at 35.1766\n",
      "Patient15 trying to get OBS at 35.6094\n",
      "Patient15 entering OBS at 35.6094\n",
      "Patient13 entering PP at 35.8941\n",
      "Patient16 trying to get OBS at 36.5130\n",
      "Patient17 trying to get OBS at 36.5369\n",
      "Patient15 leaving OBS at 36.6000\n",
      "Patient15 entering LDR at 36.6000\n",
      "Patient16 entering OBS at 36.6000\n",
      "Patient16 waited 0.0870 time units for OBS bed\n",
      "Patient16 leaving OBS at 37.4484\n",
      "Patient16 entering LDR at 37.4484\n",
      "Patient17 entering OBS at 37.4484\n",
      "Patient17 waited 0.9115 time units for OBS bed\n",
      "Patient18 trying to get OBS at 39.6783\n",
      "Patient16 entering PP at 39.7335\n",
      "Patient17 leaving OBS at 39.9760\n",
      "Patient17 entering LDR at 39.9760\n",
      "Patient18 entering OBS at 39.9760\n",
      "Patient18 waited 0.2977 time units for OBS bed\n",
      "Patient9 entering PP at 40.7127\n",
      "Patient14 leaving OBS at 41.0861\n",
      "Patient14 entering LDR at 41.0861\n",
      "Patient19 trying to get OBS at 42.1910\n",
      "Patient19 entering OBS at 42.1910\n",
      "Patient10 entering PP at 42.5130\n",
      "Patient18 leaving OBS at 42.7172\n",
      "Patient18 entering LDR at 42.7172\n",
      "Patient17 entering PP at 44.1436\n",
      "Patient18 entering PP at 44.4498\n",
      "Patient15 entering PP at 44.6092\n",
      "Patient20 trying to get OBS at 44.6790\n",
      "Patient20 entering OBS at 44.6790\n",
      "Patient3 entering PP at 44.9184\n",
      "Patient20 leaving OBS at 45.2709\n",
      "Patient20 entering LDR at 45.2709\n",
      "Patient21 trying to get OBS at 46.6338\n",
      "Patient21 entering OBS at 46.6338\n",
      "Patient21 leaving OBS at 47.8777\n",
      "Patient21 entering LDR at 47.8777\n",
      "Patient19 leaving OBS at 48.7520\n",
      "Patient19 entering LDR at 48.7520\n",
      "Patient14 entering PP at 49.3856\n",
      "Patient19 entering PP at 49.7737\n",
      "Patient22 trying to get OBS at 49.9682\n",
      "Patient22 entering OBS at 49.9682\n",
      "Patient22 leaving OBS at 50.6934\n",
      "Patient22 entering LDR at 50.6934\n",
      "Patient21 entering PP at 50.6969\n",
      "Patient23 trying to get OBS at 51.9160\n",
      "Patient23 entering OBS at 51.9160\n",
      "Patient24 trying to get OBS at 52.8406\n",
      "Patient24 entering OBS at 52.8406\n",
      "Patient24 leaving OBS at 53.3548\n",
      "Patient24 entering LDR at 53.3548\n",
      "Patient23 leaving OBS at 54.9689\n",
      "Patient23 entering LDR at 54.9689\n",
      "Patient25 trying to get OBS at 54.9893\n",
      "Patient25 entering OBS at 54.9893\n",
      "Patient24 entering PP at 55.3316\n",
      "Patient20 entering PP at 55.7203\n",
      "Patient23 entering PP at 57.8120\n",
      "Patient22 entering PP at 58.5375\n",
      "Patient26 trying to get OBS at 61.0166\n",
      "Patient26 entering OBS at 61.0166\n",
      "Patient26 leaving OBS at 61.1124\n",
      "Patient26 entering LDR at 61.1124\n",
      "Patient27 trying to get OBS at 62.3948\n",
      "Patient27 entering OBS at 62.3948\n",
      "Patient25 leaving OBS at 63.4644\n",
      "Patient25 entering LDR at 63.4644\n",
      "Patient27 leaving OBS at 63.5210\n",
      "Patient27 entering LDR at 63.5210\n",
      "Patient26 entering PP at 65.5079\n",
      "Patient28 trying to get OBS at 68.1463\n",
      "Patient28 entering OBS at 68.1463\n",
      "Patient29 trying to get OBS at 68.4204\n",
      "Patient29 entering OBS at 68.4204\n",
      "Patient25 entering PP at 68.5573\n",
      "Patient28 leaving OBS at 68.5577\n",
      "Patient28 entering LDR at 68.5577\n",
      "Patient30 trying to get OBS at 72.7648\n",
      "Patient30 entering OBS at 72.7648\n",
      "Patient30 leaving OBS at 72.7858\n",
      "Patient30 entering LDR at 72.7858\n",
      "Patient31 trying to get OBS at 72.8227\n",
      "Patient31 entering OBS at 72.8227\n",
      "Patient27 entering PP at 73.2677\n",
      "Patient29 leaving OBS at 73.3031\n",
      "Patient29 entering LDR at 73.3031\n",
      "Patient31 leaving OBS at 73.5402\n",
      "Patient31 entering LDR at 73.5402\n",
      "Patient32 trying to get OBS at 73.7931\n",
      "Patient32 entering OBS at 73.7931\n",
      "Patient30 entering PP at 74.6737\n"
     ]
    }
   ],
   "source": [
    "import simpy\n",
    "import numpy as np\n",
    "from numpy.random import RandomState\n",
    "\n",
    "\"\"\"\n",
    "Simple OB patient flow model 2 - NOT OO\n",
    "\n",
    "Details:\n",
    "\n",
    "- Generate arrivals via Poisson process\n",
    "- Uses one Resource objects to model OBS, the other units are modeled as simple delays.\n",
    "- Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.\n",
    "\n",
    "\"\"\"\n",
    "# Arrival rate and length of stay inputs.\n",
    "ARR_RATE = 0.4\n",
    "MEAN_LOS_OBS = 3\n",
    "MEAN_LOS_LDR = 12\n",
    "MEAN_LOS_PP = 48\n",
    "\n",
    "CAPACITY_OBS = 2\n",
    "\n",
    "RNG_SEED = 6353\n",
    "\n",
    "def patient_generator(env, arr_rate, prng=RandomState(0)):\n",
    "    \"\"\"Generates patients according to a simple Poisson process\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        env : simpy.Environment\n",
    "            the simulation environment\n",
    "        arr_rate : float\n",
    "            exponential arrival rate\n",
    "        prng : RandomState object\n",
    "            Seeded RandomState object for generating pseudo-random numbers.\n",
    "            See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    patients_created = 0\n",
    "\n",
    "    # Infinite loop for generatirng patients according to a poisson process.\n",
    "    while True:\n",
    "\n",
    "        # Generate next interarrival time\n",
    "        iat = prng.exponential(1.0 / arr_rate)\n",
    "\n",
    "        # Generate length of stay in each unit for this patient\n",
    "        los_obs = prng.exponential(MEAN_LOS_OBS)\n",
    "        los_ldr = prng.exponential(MEAN_LOS_LDR)\n",
    "        los_pp = prng.exponential(MEAN_LOS_PP)\n",
    "\n",
    "        # Update counter of patients\n",
    "        patients_created += 1\n",
    "\n",
    "        # Create a new patient flow process.\n",
    "        obp = obpatient_flow(env, 'Patient{}'.format(patients_created),\n",
    "                             los_obs=los_obs, los_ldr=los_ldr, los_pp=los_pp)\n",
    "\n",
    "        # Register the process with the simulation environment\n",
    "        env.process(obp)\n",
    "\n",
    "        # This process will now yield to a 'timeout' event. This process will resume after iat time units.\n",
    "        yield env.timeout(iat)\n",
    "\n",
    "\n",
    "def obpatient_flow(env, name, los_obs, los_ldr, los_pp):\n",
    "    \"\"\"Process function modeling how a patient flows through system.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        env : simpy.Environment\n",
    "            the simulation environment\n",
    "        name : str\n",
    "            process instance id\n",
    "        los_obs : float\n",
    "            length of stay in OBS unit\n",
    "        los_ldr : float\n",
    "            length of stay in LDR unit\n",
    "        los_pp : float\n",
    "            length of stay in PP unit\n",
    "    \"\"\"\n",
    "\n",
    "    print(\"{} trying to get OBS at {:.4f}\".format(name, env.now))\n",
    "\n",
    "    # Timestamp when patient tried to get OBS bed\n",
    "    bed_request_ts = env.now\n",
    "    # Request an obs bed\n",
    "    bed_request = obs_unit.request()\n",
    "    # Yield this process until a bed is available\n",
    "    yield bed_request\n",
    "\n",
    "    # We got an OBS bed\n",
    "    print(\"{} entering OBS at {:.4f}\".format(name, env.now))\n",
    "    # Let's see if we had to wait to get the bed.\n",
    "    if env.now > bed_request_ts:\n",
    "        print(\"{} waited {:.4f} time units for OBS bed\".format(name, env.now - bed_request_ts))\n",
    "\n",
    "    # Yield this process again. Now wait until our length of stay elapses.\n",
    "    # This is the actual stay in the bed\n",
    "    yield env.timeout(los_obs)\n",
    "\n",
    "    # All done with OBS, release the bed. Note that we pass the bed_request object\n",
    "    # to the release() function so that the correct unit of the resource is released.\n",
    "    obs_unit.release(bed_request)\n",
    "    print(\"{} leaving OBS at {:.4f}\".format(name, env.now))\n",
    "\n",
    "    # Continue on through LDR and PP; modeled as simple delays for now.\n",
    "\n",
    "    print(\"{} entering LDR at {:.4f}\".format(name, env.now))\n",
    "    yield env.timeout(los_ldr)\n",
    "\n",
    "    print(\"{} entering PP at {:.4f}\".format(name, env.now))\n",
    "    yield env.timeout(los_pp)\n",
    "\n",
    "\n",
    "# Initialize a simulation environment\n",
    "env = simpy.Environment()\n",
    "\n",
    "# Initialize a random number generator.\n",
    "# See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html\n",
    "prng = RandomState(RNG_SEED)\n",
    "\n",
    "# Declare a Resource to model OBS unit. Default capacity is 1, we pass in desired capacity.\n",
    "obs_unit = simpy.Resource(env, CAPACITY_OBS)\n",
    "\n",
    "# Run the simulation for a while\n",
    "runtime = 75\n",
    "env.process(patient_generator(env, ARR_RATE, prng))\n",
    "env.run(until=runtime)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seems to be working. Note that with an arrival rate of 0.4 and a runtime of 75, we expect about..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Expected number of arrivals = 30.00\n"
     ]
    }
   ],
   "source": [
    "print(\"Expected number of arrivals = {:.2f}\".format(ARR_RATE * runtime))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We got 32 with this random number seed. No cause for alarm. Of course, we should run this thing longer and compute some stats so that we can make sure things are working fine. We'll do that shortly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## obflow_3: All units modeled as a Resources\n",
    "\n",
    "Let's build on `obflow_2` by modeling all three units as resources. Since these are hospital beds, there are no queues per se. If a patient finishes their stay in OBS and there are no beds available in the LDR, they wait in OBS. Likewise, when trying to get into PP after delivering the baby in the LDR, they may get blocked in the LDR waiting for an available PP bed. In more complex models we may actually do length of stay adjustments for those patients who are blocked. For now, let's just keep it simple and not do any such adjustments.\n",
    "\n",
    "In SimPy terms, we have to be careful that we don't release a bed we are in until we have successfully reserved a bed in the next unit we are to visit. As you'll see, this leads to quite a bit of repetitive code in our process oriented approach. This will motivate us to start exploring object oriented models. \n",
    "\n",
    "Finally, we'll improve our patient generator by adding the ability to model an initial delay before patient generation should start as well as to specify a stop time for the patient generation process (which can be different than the simulation model should stop)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Patient1 trying to get OBS at 0.0000\n",
      "Patient1 entering OBS at 0.0000\n",
      "Patient1 trying to get LDR at 0.3475\n",
      "Patient1 leaving OBS at 0.3474891089551544\n",
      "Patient1 entering LDR at 0.3475\n",
      "Patient2 trying to get OBS at 3.7668\n",
      "Patient2 entering OBS at 3.7668\n",
      "Patient3 trying to get OBS at 5.6439\n",
      "Patient3 entering OBS at 5.6439\n",
      "Patient1 trying to get PP at 6.1449\n",
      "Patient1 leaving LDR at 6.1449\n",
      "Patient1 entering PP at 6.1449\n",
      "Patient4 trying to get OBS at 7.2625\n",
      "Patient5 trying to get OBS at 8.9829\n",
      "Patient3 trying to get LDR at 10.8140\n",
      "Patient3 leaving OBS at 10.813985162144862\n",
      "Patient3 entering LDR at 10.8140\n",
      "Patient4 entering OBS at 10.8140\n",
      "Patient4 waited 3.5515 time units for OBS bed\n",
      "Patient4 trying to get LDR at 10.8690\n",
      "Patient4 leaving OBS at 10.869006908518232\n",
      "Patient4 entering LDR at 10.8690\n",
      "Patient5 entering OBS at 10.8690\n",
      "Patient5 waited 1.8861 time units for OBS bed\n",
      "Patient4 trying to get PP at 11.4351\n",
      "Patient4 leaving LDR at 11.4351\n",
      "Patient4 entering PP at 11.4351\n",
      "Patient2 trying to get LDR at 12.9060\n",
      "Patient2 leaving OBS at 12.90600490911474\n",
      "Patient2 entering LDR at 12.9060\n",
      "Patient6 trying to get OBS at 16.8153\n",
      "Patient6 entering OBS at 16.8153\n",
      "Patient7 trying to get OBS at 17.8737\n",
      "Patient8 trying to get OBS at 18.7157\n",
      "Patient5 trying to get LDR at 18.8306\n",
      "Patient5 leaving OBS at 18.830564537083035\n",
      "Patient5 entering LDR at 18.8306\n",
      "Patient7 entering OBS at 18.8306\n",
      "Patient7 waited 0.9569 time units for OBS bed\n",
      "Patient7 trying to get LDR at 19.0556\n",
      "Patient7 leaving OBS at 19.05561983104953\n",
      "Patient7 entering LDR at 19.0556\n",
      "Patient8 entering OBS at 19.0556\n",
      "Patient8 waited 0.3400 time units for OBS bed\n",
      "Patient7 trying to get PP at 19.0901\n",
      "Patient7 leaving LDR at 19.0901\n",
      "Patient7 entering PP at 19.0901\n",
      "Patient6 trying to get LDR at 19.4866\n",
      "Patient6 leaving OBS at 19.48663625825889\n",
      "Patient6 entering LDR at 19.4866\n",
      "Patient9 trying to get OBS at 20.4913\n",
      "Patient9 entering OBS at 20.4913\n",
      "Patient8 trying to get LDR at 23.1856\n",
      "Patient8 leaving OBS at 23.185576349161767\n",
      "Patient8 entering LDR at 23.1856\n",
      "Patient8 trying to get PP at 23.5282\n",
      "Patient8 leaving LDR at 23.5282\n",
      "Patient8 entering PP at 23.5282\n",
      "Patient2 trying to get PP at 25.7847\n",
      "Patient2 leaving LDR at 25.7847\n",
      "Patient2 entering PP at 25.7847\n",
      "Patient5 trying to get PP at 26.4507\n",
      "Patient5 leaving LDR at 26.4507\n",
      "Patient5 entering PP at 26.4507\n",
      "Patient10 trying to get OBS at 26.8687\n",
      "Patient10 entering OBS at 26.8687\n",
      "Patient10 trying to get LDR at 27.6307\n",
      "Patient10 leaving OBS at 27.63070937211375\n",
      "Patient10 entering LDR at 27.6307\n",
      "Patient11 trying to get OBS at 27.6876\n",
      "Patient11 entering OBS at 27.6876\n",
      "Patient6 trying to get PP at 28.3114\n",
      "Patient6 leaving LDR at 28.3114\n",
      "Patient6 entering PP at 28.3114\n",
      "Patient9 trying to get LDR at 30.4208\n",
      "Patient9 leaving OBS at 30.420840322655998\n",
      "Patient9 entering LDR at 30.4208\n",
      "Patient11 trying to get LDR at 30.4652\n",
      "Patient11 leaving OBS at 30.465221449496795\n",
      "Patient11 entering LDR at 30.4652\n",
      "Patient11 trying to get PP at 30.5128\n",
      "Patient11 leaving LDR at 30.5128\n",
      "Patient11 entering PP at 30.5128\n",
      "Patient12 trying to get OBS at 31.4872\n",
      "Patient12 entering OBS at 31.4872\n",
      "Patient13 trying to get OBS at 31.9981\n",
      "Patient13 entering OBS at 31.9981\n",
      "Patient12 trying to get LDR at 32.2156\n",
      "Patient12 leaving OBS at 32.215634051591636\n",
      "Patient12 entering LDR at 32.2156\n",
      "Patient13 trying to get LDR at 32.5002\n",
      "Patient13 leaving OBS at 32.500238016146994\n",
      "Patient13 entering LDR at 32.5002\n",
      "Patient14 trying to get OBS at 33.5398\n",
      "Patient14 entering OBS at 33.5398\n",
      "Patient7 leaving PP and system at 34.5316\n",
      "Patient12 trying to get PP at 35.1766\n",
      "Patient12 leaving LDR at 35.1766\n",
      "Patient12 entering PP at 35.1766\n",
      "Patient15 trying to get OBS at 35.6094\n",
      "Patient15 entering OBS at 35.6094\n",
      "Patient13 trying to get PP at 35.8941\n",
      "Patient13 leaving LDR at 35.8941\n",
      "Patient13 entering PP at 35.8941\n",
      "Patient16 trying to get OBS at 36.5130\n",
      "Patient17 trying to get OBS at 36.5369\n",
      "Patient15 trying to get LDR at 36.6000\n",
      "Patient15 leaving OBS at 36.59999969293714\n",
      "Patient15 entering LDR at 36.6000\n",
      "Patient16 entering OBS at 36.6000\n",
      "Patient16 waited 0.0870 time units for OBS bed\n",
      "Patient16 trying to get LDR at 37.4484\n",
      "Patient16 leaving OBS at 37.44839668463982\n",
      "Patient16 entering LDR at 37.4484\n",
      "Patient17 entering OBS at 37.4484\n",
      "Patient17 waited 0.9115 time units for OBS bed\n",
      "Patient11 leaving PP and system at 38.1145\n",
      "Patient18 trying to get OBS at 39.6783\n",
      "Patient16 trying to get PP at 39.7335\n",
      "Patient16 leaving LDR at 39.7335\n",
      "Patient16 entering PP at 39.7335\n",
      "Patient17 trying to get LDR at 39.9760\n",
      "Patient17 leaving OBS at 39.97599817647182\n",
      "Patient17 entering LDR at 39.9760\n",
      "Patient18 entering OBS at 39.9760\n",
      "Patient18 waited 0.2977 time units for OBS bed\n",
      "Patient9 trying to get PP at 40.7127\n",
      "Patient9 leaving LDR at 40.7127\n",
      "Patient9 entering PP at 40.7127\n",
      "Patient14 trying to get LDR at 41.0861\n",
      "Patient14 leaving OBS at 41.086116528209494\n",
      "Patient14 entering LDR at 41.0861\n",
      "Patient19 trying to get OBS at 42.1910\n",
      "Patient19 entering OBS at 42.1910\n",
      "Patient10 trying to get PP at 42.5130\n",
      "Patient10 leaving LDR at 42.5130\n",
      "Patient10 entering PP at 42.5130\n",
      "Patient18 trying to get LDR at 42.7172\n",
      "Patient18 leaving OBS at 42.71719377930215\n",
      "Patient18 entering LDR at 42.7172\n",
      "Patient1 leaving PP and system at 42.8634\n",
      "Patient17 trying to get PP at 44.1436\n",
      "Patient17 leaving LDR at 44.1436\n",
      "Patient17 entering PP at 44.1436\n",
      "Patient18 trying to get PP at 44.4498\n",
      "Patient18 leaving LDR at 44.4498\n",
      "Patient18 entering PP at 44.4498\n",
      "Patient15 trying to get PP at 44.6092\n",
      "Patient15 leaving LDR at 44.6092\n",
      "Patient15 entering PP at 44.6092\n",
      "Patient20 trying to get OBS at 44.6790\n",
      "Patient20 entering OBS at 44.6790\n",
      "Patient3 trying to get PP at 44.9184\n",
      "Patient3 leaving LDR at 44.9184\n",
      "Patient3 entering PP at 44.9184\n",
      "Patient20 trying to get LDR at 45.2709\n",
      "Patient20 leaving OBS at 45.270857202928546\n",
      "Patient20 entering LDR at 45.2709\n",
      "Patient21 trying to get OBS at 46.6338\n",
      "Patient21 entering OBS at 46.6338\n",
      "Patient21 trying to get LDR at 47.8777\n",
      "Patient21 leaving OBS at 47.8776934769442\n",
      "Patient21 entering LDR at 47.8777\n",
      "Patient19 trying to get LDR at 48.7520\n",
      "Patient19 leaving OBS at 48.75197164898342\n",
      "Patient19 entering LDR at 48.7520\n",
      "Patient14 trying to get PP at 49.3856\n",
      "Patient14 leaving LDR at 49.3856\n",
      "Patient14 entering PP at 49.3856\n",
      "Patient19 trying to get PP at 49.7737\n",
      "Patient19 leaving LDR at 49.7737\n",
      "Patient19 entering PP at 49.7737\n",
      "Patient22 trying to get OBS at 49.9682\n",
      "Patient22 entering OBS at 49.9682\n",
      "Patient22 trying to get LDR at 50.6934\n",
      "Patient22 leaving OBS at 50.69343889384862\n",
      "Patient22 entering LDR at 50.6934\n",
      "Patient21 trying to get PP at 50.6969\n",
      "Patient21 leaving LDR at 50.6969\n",
      "Patient21 entering PP at 50.6969\n",
      "Patient20 trying to get PP at 55.7203\n",
      "Patient20 leaving LDR at 55.7203\n",
      "Patient20 entering PP at 55.7203\n",
      "Patient16 leaving PP and system at 56.0439\n",
      "Patient18 leaving PP and system at 57.0161\n",
      "Patient12 leaving PP and system at 58.3965\n",
      "Patient22 trying to get PP at 58.5375\n",
      "Patient22 leaving LDR at 58.5375\n",
      "Patient22 entering PP at 58.5375\n",
      "Patient10 leaving PP and system at 59.8131\n",
      "Patient9 leaving PP and system at 60.4868\n",
      "Patient5 leaving PP and system at 63.4959\n",
      "Patient17 leaving PP and system at 69.8364\n",
      "Patient14 leaving PP and system at 69.9692\n",
      "Patient15 leaving PP and system at 70.5547\n",
      "Patient22 leaving PP and system at 72.2739\n",
      "Patient20 leaving PP and system at 73.4076\n"
     ]
    }
   ],
   "source": [
    "import simpy\n",
    "import numpy as np\n",
    "from numpy.random import RandomState\n",
    "\n",
    "\"\"\"\n",
    "Simple OB patient flow model 3 - NOT OO\n",
    "\n",
    "Details:\n",
    "\n",
    "- Generate arrivals via Poisson process\n",
    "- Uses one Resource objects to model OBS, LDR, and PP.\n",
    "- Arrival rates and mean lengths of stay hard coded as constants. Later versions will read these from input files.\n",
    "- Additional functionality added to arrival generator (initial delay and arrival stop time).\n",
    "\n",
    "\"\"\"\n",
    "ARR_RATE = 0.4\n",
    "MEAN_LOS_OBS = 3\n",
    "MEAN_LOS_LDR = 12\n",
    "MEAN_LOS_PP = 48\n",
    "\n",
    "CAPACITY_OBS = 2\n",
    "CAPACITY_LDR = 6\n",
    "CAPACITY_PP = 24\n",
    "\n",
    "RNG_SEED = 6353\n",
    "\n",
    "def patient_generator(env, arr_stream, arr_rate, initial_delay=0,\n",
    "                      stoptime=simpy.core.Infinity, prng=RandomState(0)):\n",
    "    \"\"\"Generates patients according to a simple Poisson process\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        env : simpy.Environment\n",
    "            the simulation environment\n",
    "        arr_rate : float\n",
    "            exponential arrival rate\n",
    "        initial_delay: float (default 0)\n",
    "            time before arrival generation should begin\n",
    "        stoptime: float (default Infinity)\n",
    "            time after which no arrivals are generated\n",
    "        prng : RandomState object\n",
    "            Seeded RandomState object for generating pseudo-random numbers.\n",
    "            See https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.RandomState.html\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    patients_created = 0\n",
    "\n",
    "    # Yield for the initial delay\n",
    "    yield env.timeout(initial_delay)\n",
    "\n",
    "    # Generate arrivals as long as simulation time is before stoptime\n",
    "    while env.now < stoptime:\n",
    "\n",
    "        iat = prng.exponential(1.0 / arr_rate)\n",
    "\n",
    "        # Sample los distributions\n",
    "        los_obs = prng.exponential(MEAN_LOS_OBS)\n",
    "        los_ldr = prng.exponential(MEAN_LOS_LDR)\n",
    "        los_pp = prng.exponential(MEAN_LOS_PP)\n",
    "\n",
    "\n",
    "        # Create new patient process instance\n",
    "        patients_created += 1\n",
    "        obp = obpatient_flow(env, 'Patient{}'.format(patients_created),\n",
    "                             los_obs=los_obs, los_ldr=los_ldr, los_pp=los_pp)\n",
    "\n",
    "        env.process(obp)\n",
    "\n",
    "        # Compute next interarrival time\n",
    "\n",
    "        yield env.timeout(iat)\n",
    "\n",
    "\n",
    "def obpatient_flow(env, name, los_obs, los_ldr, los_pp):\n",
    "    \"\"\"Process function modeling how a patient flows through system.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        env : simpy.Environment\n",
    "            the simulation environment\n",
    "        name : str\n",
    "            process instance id\n",
    "        los_obs : float\n",
    "            length of stay in OBS unit\n",
    "        los_ldr : float\n",
    "            length of stay in LDR unit\n",
    "        los_pp : float\n",
    "            length of stay in PP unit\n",
    "    \"\"\"\n",
    "\n",
    "    # Note the repetitive code and the use of separate request objects for each\n",
    "    # stay in the different units.\n",
    "\n",
    "    # OBS\n",
    "    print(\"{} trying to get OBS at {:.4f}\".format(name, env.now))\n",
    "    bed_request_ts = env.now\n",
    "    bed_request1 = obs_unit.request() # Request an OBS bed\n",
    "    yield bed_request1\n",
    "    print(\"{} entering OBS at {:.4f}\".format(name, env.now))\n",
    "    if env.now > bed_request_ts:\n",
    "        print(\"{} waited {:.4f} time units for OBS bed\".format(name, env.now-  bed_request_ts))\n",
    "    yield env.timeout(los_obs) # Stay in obs bed\n",
    "\n",
    "    print(\"{} trying to get LDR at {:.4f}\".format(name, env.now))\n",
    "    bed_request_ts = env.now\n",
    "    bed_request2 = ldr_unit.request()  # Request an LDR bed\n",
    "    yield bed_request2\n",
    "\n",
    "    # Got LDR bed, release OBS bed\n",
    "    obs_unit.release(bed_request1)  # Release the OBS bed\n",
    "    print(\"{} leaving OBS at {}\".format(name, env.now))\n",
    "\n",
    "    # LDR stay\n",
    "    print(\"{} entering LDR at {:.4f}\".format(name, env.now))\n",
    "    if env.now > bed_request_ts:\n",
    "        print(\"{} waited {:.4f} time units for LDR bed\".format(name, env.now - bed_request_ts))\n",
    "    yield env.timeout(los_ldr) # Stay in LDR bed\n",
    "\n",
    "    print(\"{} trying to get PP at {:.4f}\".format(name, env.now))\n",
    "    bed_request_ts = env.now\n",
    "    bed_request3 = pp_unit.request()  # Request a PP bed\n",
    "    yield bed_request3\n",
    "\n",
    "    # Got PP bed, release LDR bed\n",
    "    ldr_unit.release(bed_request2)  # Release the obs bed\n",
    "    print(\"{} leaving LDR at {:.4f}\".format(name, env.now))\n",
    "\n",
    "    # PP stay\n",
    "    print(\"{} entering PP at {:.4f}\".format(name, env.now))\n",
    "    if env.now > bed_request_ts:\n",
    "        print(\"{} waited {:.4f} time units for PP bed\".format(name, env.now - bed_request_ts))\n",
    "    yield env.timeout(los_pp) # Stay in PP bed\n",
    "    pp_unit.release(bed_request3)  # Release the PP bed\n",
    "\n",
    "    print(\"{} leaving PP and system at {:.4f}\".format(name, env.now))\n",
    "\n",
    "# Initialize a simulation environment\n",
    "env = simpy.Environment()\n",
    "\n",
    "prng = RandomState(RNG_SEED)\n",
    "\n",
    "rho_obs = ARR_RATE * MEAN_LOS_OBS / CAPACITY_OBS\n",
    "rho_ldr = ARR_RATE * MEAN_LOS_LDR / CAPACITY_LDR\n",
    "rho_pp = ARR_RATE * MEAN_LOS_PP / CAPACITY_PP\n",
    "\n",
    "# Declare Resources to model all units\n",
    "obs_unit = simpy.Resource(env, CAPACITY_OBS)\n",
    "ldr_unit = simpy.Resource(env, CAPACITY_LDR)\n",
    "pp_unit = simpy.Resource(env, CAPACITY_PP)\n",
    "\n",
    "# Run the simulation for a while. Let's shut arrivals off after 50 time units.\n",
    "runtime = 75\n",
    "stop_arrivals = 50\n",
    "env.process(patient_generator(env, \"Type1\", ARR_RATE, 0, stop_arrivals, prng))\n",
    "env.run(until=runtime)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice we get far fewer patients because we turned off the arrival spigot at time 50.\n",
    "\n",
    "As a veteran simulation modeler, I'm not really liking this process oriented approach. It's not going to scale well as we add more units (e.g. if we were modeling all the inpatient beds) and more routing complexity. Having built a pretty big object oriented DES model (in Java with Simkit) for analyzing pneumatic tube systems, I know that that's the way to go. In the next post, I'll share my first object oriented SimPy model. Not only does an object oriented approach let us avoid so much code repitition, it also makes collecting stats and logging easier."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:hselab]",
   "language": "python",
   "name": "conda-env-hselab-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}