diff --git a/content/02_basic_simpy.ipynb b/content/02_basic_simpy.ipynb index b2887a8..488db42 100644 --- a/content/02_basic_simpy.ipynb +++ b/content/02_basic_simpy.ipynb @@ -230,7 +230,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.9" + "version": "3.11.11" } }, "nbformat": 4, diff --git a/content/03a_exercise1.ipynb b/content/03a_exercise1.ipynb index 754412c..615d044 100644 --- a/content/03a_exercise1.ipynb +++ b/content/03a_exercise1.ipynb @@ -142,7 +142,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.9" + "version": "3.11.11" } }, "nbformat": 4, diff --git a/content/13_warm_up.ipynb b/content/13_warm_up.ipynb new file mode 100644 index 0000000..c3ea9cc --- /dev/null +++ b/content/13_warm_up.ipynb @@ -0,0 +1,522 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0be7dabf-cb34-4faf-abb1-e2c8e735beda", + "metadata": {}, + "source": [ + "# Coding a warm-up period in SimPy\n", + "\n", + "## Why do you need a warm-up period?\n", + "\n", + "Typically when you are modelling a non-terminating system, you will need to deal with **initialisation bias**. That is the real system always has work-in-progress (e.g. patients in queues and in service), but the model starts from empty. One way to deal with this is to split the model's run length into warm-up and data collection periods. We discard all results in the warm-up period.\n", + "\n", + "> In this tutorial we will focus on coding a warm-up period in Python and SimPy rather than analysis to determine its length\n", + "\n", + "## But how do you code it?\n", + "\n", + "💪 We will implement warm-up as a **single event** that resets all of our results collection variables. \n", + "\n", + "> This is a simpler approach than including lots of if statements in `simpy` processes.\n", + "\n", + "## Illustrative example model\n", + "\n", + "We will use a very simple model for this example. This is a acute stroke pathway with a single arrival processes, a single type of resource, and a single treatment process. This is a non-terminating system. There are always patients in the system - it does not start up from empty\n", + "\n", + "\n", + "![model image](img/acute_stroke_pathway.png \"stroke pathway\")" + ] + }, + { + "cell_type": "markdown", + "id": "0d9383eb-420c-49f8-b178-f2fe9e6b2a90", + "metadata": {}, + "source": [ + "## 1. Imports " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1cee9f9-8696-4b13-94ff-bee2a2a2e5f8", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import itertools\n", + "import simpy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea3d507f-9e6d-4ff0-8b90-f9c63c8a8bdf", + "metadata": {}, + "outputs": [], + "source": [ + "# to reduce code these classes can be found in distribution.py\n", + "from distributions import (\n", + " Exponential, \n", + " Lognormal\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "c422046d-488a-4743-8ad4-97e9f3dab420", + "metadata": {}, + "source": [ + "## 2. Constants" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ecf0429-f03f-4ad2-abb4-46692a74e559", + "metadata": {}, + "outputs": [], + "source": [ + "# default mean inter-arrival times(exp)\n", + "# time is in days\n", + "IAT_STROKES = 1.0\n", + "\n", + "# resources\n", + "N_ACUTE_BEDS = 9\n", + "\n", + "# Acute LoS (Lognormal)\n", + "ACUTE_LOS_MEAN = 7.0\n", + "ACUTE_LOC_STD = 1.0\n", + "\n", + "# sampling settings\n", + "N_STREAMS = 2\n", + "DEFAULT_RND_SET = 0\n", + "\n", + "# Boolean switch to simulation results as the model runs\n", + "TRACE = False\n", + "\n", + "# run variables (units = days)\n", + "WU_PERIOD = 0.0\n", + "RC_PERIOD = 100.0" + ] + }, + { + "cell_type": "markdown", + "id": "5f2a4ad9-6d5e-480d-850f-84d4882a738b", + "metadata": {}, + "source": [ + "## 2. Helper classes and functions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52c9271f-1d05-454d-a199-8768bdf5b6e8", + "metadata": {}, + "outputs": [], + "source": [ + "def trace(msg):\n", + " \"\"\"\n", + " Turing printing of events on and off.\n", + "\n", + " Params:\n", + " -------\n", + " msg: str\n", + " string to print to screen.\n", + " \"\"\"\n", + " if TRACE:\n", + " print(msg)" + ] + }, + { + "cell_type": "markdown", + "id": "5a8c050c-4bb6-408f-a805-3a4aaab56916", + "metadata": {}, + "source": [ + "## 3. Experiment class" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "576ae9b4-b21b-4ed0-9b13-e5898d423173", + "metadata": {}, + "outputs": [], + "source": [ + "class Experiment:\n", + " \"\"\"\n", + " Encapsulates the concept of an experiment 🧪 for the stroke pathway\n", + " bed blocking simulator. Manages parameters, PRNG streams and results.\n", + " \"\"\"\n", + " def __init__(\n", + " self,\n", + " random_number_set=DEFAULT_RND_SET,\n", + " n_streams=N_STREAMS,\n", + " iat_strokes=IAT_STROKES,\n", + " acute_los_mean=ACUTE_LOS_MEAN,\n", + " acute_los_std=ACUTE_LOC_STD,\n", + " n_acute_beds=N_ACUTE_BEDS, \n", + " ):\n", + " \"\"\"\n", + " The init method sets up our defaults.\n", + " \"\"\"\n", + " # sampling\n", + " self.random_number_set = random_number_set\n", + " self.n_streams = n_streams\n", + "\n", + " # store parameters for the run of the model\n", + " self.iat_strokes = iat_strokes\n", + " self.acute_los_mean = acute_los_mean\n", + " self.acute_los_std = acute_los_std\n", + "\n", + " # place holder for resources\n", + " self.acute_ward = None\n", + " self.n_acute_beds = n_acute_beds\n", + " \n", + " # initialise results to zero\n", + " self.init_results_variables()\n", + "\n", + " # initialise sampling objects\n", + " self.init_sampling()\n", + "\n", + " def set_random_no_set(self, random_number_set):\n", + " \"\"\"\n", + " Controls the random sampling\n", + " Parameters:\n", + " ----------\n", + " random_number_set: int\n", + " Used to control the set of pseudo random numbers used by\n", + " the distributions in the simulation.\n", + " \"\"\"\n", + " self.random_number_set = random_number_set\n", + " self.init_sampling()\n", + "\n", + " def init_sampling(self):\n", + " \"\"\"\n", + " Create the distributions used by the model and initialise\n", + " the random seeds of each.\n", + " \"\"\"\n", + " # produce n non-overlapping streams\n", + " seed_sequence = np.random.SeedSequence(self.random_number_set)\n", + " self.seeds = seed_sequence.spawn(self.n_streams)\n", + "\n", + " # create distributions\n", + "\n", + " # inter-arrival time distributions\n", + " self.arrival_strokes = Exponential(\n", + " self.iat_strokes, random_seed=self.seeds[0]\n", + " )\n", + "\n", + " self.acute_los = Lognormal(\n", + " self.acute_los_mean, self.acute_los_std, random_seed=self.seeds[1]\n", + " )\n", + "\n", + " def init_results_variables(self):\n", + " \"\"\"\n", + " Initialise all of the experiment variables used in results\n", + " collection. This method is called at the start of each run\n", + " of the model\n", + " \"\"\"\n", + " # variable used to store results of experiment\n", + " self.results = {}\n", + " self.results[\"n_arrivals\"] = 0\n", + " self.results[\"waiting_acute\"] = []" + ] + }, + { + "cell_type": "markdown", + "id": "7ff9beae-89cc-419c-b584-c05b81086865", + "metadata": {}, + "source": [ + "## 🥵 Warm-up period\n", + "\n", + "The acute stroke pathway model starts from empty. As it is a non-terminating system our estimate of waiting time is biased due to the empty period at the start of the simulation. We can remove this initialisation bias using a warm-up period. \n", + "\n", + "We will implement a warm-up through an **event** that happens once in a single run of the model. The model will be run for the **warm-up period + results collection period**. At the end of the warm-up period an event will happen where all variables in the current experiment are reset (e.g. empty lists and set quantitative values to 0)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dff74a08-37fd-4b18-8bcd-97994f38369a", + "metadata": {}, + "outputs": [], + "source": [ + "def warmup_complete(warm_up_period, env, args):\n", + " \"\"\"\n", + " End of warm-up period event. Used to reset results collection variables.\n", + "\n", + " Parameters:\n", + " ----------\n", + " warm_up_period: float\n", + " Duration of warm-up period in simulation time units\n", + "\n", + " env: simpy.Environment\n", + " The simpy environment\n", + "\n", + " args: Experiment\n", + " The simulation experiment that contains the results being collected.\n", + " \"\"\"\n", + " yield env.timeout(warm_up_period)\n", + " trace(f\"{env.now:.2f}: 🥵 Warm up complete.\")\n", + " \n", + " args.init_results_variables()" + ] + }, + { + "cell_type": "markdown", + "id": "94f0f9c5-22cb-493a-9f1f-4e2a8325beaa", + "metadata": {}, + "source": [ + "## 4. Stroke pathway process logic\n", + "\n", + "The key things to recognise are \n", + "\n", + "* We include a optional parameter called `collection_results` that defaults to `True`. We may set this `False` in our functions that setup initial conditions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "911528e1-e4eb-4307-bb26-632faf7769d1", + "metadata": {}, + "outputs": [], + "source": [ + "def acute_stroke_pathway(patient_id, env, args):\n", + " \"\"\"Process a patient through the acute ward\n", + " Simpy generator function.\n", + " \n", + " Parameters:\n", + " -----------\n", + " patient_id: int\n", + " A unique id representing the patient in the process\n", + "\n", + " env: simpy.Environment\n", + " The simulation environment\n", + "\n", + " args: Experiment\n", + " Container class for the simulation parameters/results.\n", + " \"\"\"\n", + " arrival_time = env.now\n", + "\n", + " with args.acute_ward.request() as acute_bed_request:\n", + " yield acute_bed_request\n", + " \n", + " acute_admit_time = env.now\n", + " wait_for_acute = acute_admit_time - arrival_time\n", + " \n", + " args.results['waiting_acute'].append(wait_for_acute)\n", + " \n", + " trace(f\"{env.now:.2f}: Patient {patient_id} admitted to acute ward.\" \\\n", + " + f\"(waited {wait_for_acute:.2f} days)\")\n", + " \n", + " # Simulate acute care treatment\n", + " acute_care_duration = args.acute_los.sample()\n", + " yield env.timeout(acute_care_duration)\n", + " \n", + " trace(f\"{env.now:.2f}: Patient {patient_id} discharged.\")" + ] + }, + { + "cell_type": "markdown", + "id": "de8990c2-a330-4c02-ac77-26c30d3e0a41", + "metadata": {}, + "source": [ + "## 4. Arrivals generator\n", + "\n", + "This is a standard arrivals generator. We create stroke arrivals according to their distribution." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3e686ce-5371-4471-a052-b9d43309bc85", + "metadata": {}, + "outputs": [], + "source": [ + "def stroke_arrivals_generator(env, args):\n", + " \"\"\"\n", + " Arrival process for strokes.\n", + "\n", + " Parameters:\n", + " ------\n", + " env: simpy.Environment\n", + " The simpy environment for the simulation\n", + "\n", + " args: Experiment\n", + " The settings and input parameters for the simulation.\n", + " \"\"\"\n", + " # use itertools as it provides an infinite loop\n", + " # with a counter variable that we can use for unique Ids\n", + " for patient_id in itertools.count(start=1):\n", + "\n", + " # the sample distribution is defined by the experiment.\n", + " inter_arrival_time = args.arrival_strokes.sample()\n", + " yield env.timeout(inter_arrival_time)\n", + "\n", + " args.results[\"n_arrivals\"] += 1\n", + " \n", + " trace(f\"{env.now:.2f}: Stroke arrival.\")\n", + "\n", + " # patient enters pathway\n", + " env.process(acute_stroke_pathway(patient_id, env, args))" + ] + }, + { + "cell_type": "markdown", + "id": "6058571e-9fdb-4961-be27-8a3b8c2fe26e", + "metadata": {}, + "source": [ + "## 5. Single run function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d0ea6cf-7d95-4d2c-9690-fcdbdae35d84", + "metadata": {}, + "outputs": [], + "source": [ + "def single_run(\n", + " experiment, \n", + " rep=0,\n", + " wu_period=WU_PERIOD,\n", + " rc_period=RC_PERIOD\n", + "):\n", + " \"\"\"\n", + " Perform a single run of the model and return the results\n", + "\n", + " Parameters:\n", + " -----------\n", + "\n", + " experiment: Experiment\n", + " The experiment/paramaters to use with model\n", + "\n", + " rep: int\n", + " The replication number.\n", + "\n", + " wu_period: float, optional (default=WU_PERIOD)\n", + " Warm-up period\n", + "\n", + " rc_period: float, optional (default=RC_PERIOD)\n", + " The run length of the model\n", + " \"\"\"\n", + "\n", + " # reset all results variables to zero and empty\n", + " experiment.init_results_variables()\n", + "\n", + " # set random number set to the replication no.\n", + " # this controls sampling for the run.\n", + " experiment.set_random_no_set(rep)\n", + "\n", + " # environment is (re)created inside single run\n", + " env = simpy.Environment()\n", + "\n", + " # simpy resources\n", + " experiment.acute_ward = simpy.Resource(env, experiment.n_acute_beds)\n", + "\n", + " # schedule a warm_up period\n", + " env.process(warmup_complete(wu_period, env, experiment))\n", + " \n", + " # we pass all arrival generators to simpy \n", + " env.process(stroke_arrivals_generator(env, experiment))\n", + "\n", + " # run model\n", + " env.run(until=wu_period+rc_period)\n", + "\n", + " # quick stats\n", + " results = {}\n", + " results['mean_acute_wait'] = np.array(\n", + " experiment.results[\"waiting_acute\"]\n", + " ).mean()\n", + "\n", + " # return single run results\n", + " return results" + ] + }, + { + "cell_type": "markdown", + "id": "c13f5e57-723c-409b-a1ce-cdb831b4e166", + "metadata": {}, + "source": [ + "## No warm-up" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "caf52390-5455-4fa1-bb22-60b5b91ad8d0", + "metadata": {}, + "outputs": [], + "source": [ + "TRACE = True\n", + "experiment = Experiment()\n", + "results = single_run(experiment, rep=2, wu_period=10.0, rc_period=6.0)\n", + "results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ddedb4f1-207d-4295-9ae4-c49b2c7cdcaf", + "metadata": {}, + "outputs": [], + "source": [ + "# check how many patient waiting times recorded.\n", + "experiment.results" + ] + }, + { + "cell_type": "markdown", + "id": "660ea2e1-d9c2-4355-876c-43dfd9dab0fe", + "metadata": {}, + "source": [ + "## Include a warm-up" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72b5284a-1fcb-4126-b663-c0ef0002e4bf", + "metadata": {}, + "outputs": [], + "source": [ + "TRACE = True\n", + "experiment = Experiment()\n", + "results = single_run(experiment, rep=0, wu_period=5.0, rc_period=1.0)\n", + "results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f5e282b-0f41-41df-bdca-f128e7d418c1", + "metadata": {}, + "outputs": [], + "source": [ + "# check how many patient waiting times recorded.\n", + "experiment.results" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/content/14_initial_conditions.ipynb b/content/14_initial_conditions.ipynb new file mode 100644 index 0000000..4ae93ad --- /dev/null +++ b/content/14_initial_conditions.ipynb @@ -0,0 +1,644 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0be7dabf-cb34-4faf-abb1-e2c8e735beda", + "metadata": {}, + "source": [ + "# Setting initial conditions\n", + "\n", + "\n", + "**IMPORTANT: There is still an initialisation bias problem:**\n", + "\n", + "* We load patients into the model before we begin the run. \n", + "* If there are $n$ acute stroke beds then the first $n$ patients loaded into the queue will begin service immediately and have a zero queuing time.\n", + "* This applies if we have one queue in the model or if we have multiple queues and activities: the time in system metrics will be biased for initial condition patients.\n", + "* This is the same problem faced when starting the model from time zero with no patients.\n", + "\n", + "**Some Options:**\n", + "* Include a setting in a process to switch results collection on and off.\n", + "* Code a seperate process for initial conditions that does not include results collection.\n", + "* Mixed initial conditions and (a shorter) warm-up period.\n", + " * You will need to do some analysis to ensure this is working acceptably.\n", + " * The warm-up pushes patients into service and also resets results collection. (deletes an initial transient)." + ] + }, + { + "cell_type": "markdown", + "id": "0d9383eb-420c-49f8-b178-f2fe9e6b2a90", + "metadata": {}, + "source": [ + "## 1. Imports " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "c1cee9f9-8696-4b13-94ff-bee2a2a2e5f8", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import itertools\n", + "import simpy" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "ea3d507f-9e6d-4ff0-8b90-f9c63c8a8bdf", + "metadata": {}, + "outputs": [], + "source": [ + "# to reduce code these classes can be found in distribution.py\n", + "from distributions import (\n", + " Exponential, \n", + " Lognormal, \n", + " DiscreteEmpirical,\n", + " FixedDistribution\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "c422046d-488a-4743-8ad4-97e9f3dab420", + "metadata": {}, + "source": [ + "## 2. Constants" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "1ecf0429-f03f-4ad2-abb4-46692a74e559", + "metadata": {}, + "outputs": [], + "source": [ + "# default mean inter-arrival times(exp)\n", + "# time is in days\n", + "IAT_STROKES = 1.0\n", + "\n", + "# resources\n", + "N_ACUTE_BEDS = 9\n", + "\n", + "# Acute LoS (Lognormal)\n", + "ACUTE_LOS_MEAN = 7.0\n", + "ACUTE_LOC_STD = 1.0\n", + "\n", + "# initial conditions for acute queue + service \n", + "# if there are 9 beds then 10 = 1 queuing\n", + "# < 9 = 0 queuing etc.\n", + "# these can be fixed or random\n", + "# Note we are adding N_ACUTE_BEDS patients to queue lengths\n", + "\n", + "INIT_COND_PARAMS = {\n", + " \"mode\": \"fixed\",\n", + " \"fixed\": 8,\n", + " \"rnd\": {\n", + " \"values\":[8, 9, 10, 11, 12, 13],\n", + " \"freq\":[25, 25, 2, 1, 1, 1]\n", + " }\n", + "}\n", + " \n", + "# sampling settings\n", + "N_STREAMS = 3\n", + "DEFAULT_RND_SET = 0\n", + "\n", + "# Boolean switch to simulation results as the model runs\n", + "TRACE = False\n", + "\n", + "# run variables (units = days)\n", + "WU_PERIOD = 0.0\n", + "RC_PERIOD = 100" + ] + }, + { + "cell_type": "markdown", + "id": "5f2a4ad9-6d5e-480d-850f-84d4882a738b", + "metadata": {}, + "source": [ + "## 2. Helper classes and functions" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "52c9271f-1d05-454d-a199-8768bdf5b6e8", + "metadata": {}, + "outputs": [], + "source": [ + "def trace(msg):\n", + " \"\"\"\n", + " Turing printing of events on and off.\n", + "\n", + " Params:\n", + " -------\n", + " msg: str\n", + " string to print to screen.\n", + " \"\"\"\n", + " if TRACE:\n", + " print(msg)" + ] + }, + { + "cell_type": "markdown", + "id": "5a8c050c-4bb6-408f-a805-3a4aaab56916", + "metadata": {}, + "source": [ + "## 3. Experiment class" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "576ae9b4-b21b-4ed0-9b13-e5898d423173", + "metadata": {}, + "outputs": [], + "source": [ + "class Experiment:\n", + " \"\"\"\n", + " Encapsulates the concept of an experiment 🧪 for the stroke pathway\n", + " bed blocking simulator. Manages parameters, PRNG streams and results.\n", + " \"\"\"\n", + " def __init__(\n", + " self,\n", + " random_number_set=DEFAULT_RND_SET,\n", + " n_streams=N_STREAMS,\n", + " iat_strokes=IAT_STROKES,\n", + " acute_los_mean=ACUTE_LOS_MEAN,\n", + " acute_los_std=ACUTE_LOC_STD,\n", + " n_acute_beds=N_ACUTE_BEDS,\n", + " init_cond_params=INIT_COND_PARAMS, \n", + " ):\n", + " \"\"\"\n", + " The init method sets up our defaults.\n", + " \"\"\"\n", + " # sampling\n", + " self.random_number_set = random_number_set\n", + " self.n_streams = n_streams\n", + "\n", + " # store parameters for the run of the model\n", + " self.iat_strokes = iat_strokes\n", + " self.acute_los_mean = acute_los_mean\n", + " self.acute_los_std = acute_los_std\n", + "\n", + " # stored initial conditions\n", + " self.init_cond_params = init_cond_params\n", + "\n", + " # place holder for resources\n", + " self.acute_ward = None\n", + " self.n_acute_beds = n_acute_beds\n", + " \n", + " # initialise results to zero\n", + " self.init_results_variables()\n", + "\n", + " # initialise sampling objects\n", + " self.init_sampling()\n", + "\n", + " def set_random_no_set(self, random_number_set):\n", + " \"\"\"\n", + " Controls the random sampling\n", + " Parameters:\n", + " ----------\n", + " random_number_set: int\n", + " Used to control the set of pseudo random numbers used by\n", + " the distributions in the simulation.\n", + " \"\"\"\n", + " self.random_number_set = random_number_set\n", + " self.init_sampling()\n", + "\n", + " def init_sampling(self):\n", + " \"\"\"\n", + " Create the distributions used by the model and initialise\n", + " the random seeds of each.\n", + " \"\"\"\n", + " # produce n non-overlapping streams\n", + " seed_sequence = np.random.SeedSequence(self.random_number_set)\n", + " self.seeds = seed_sequence.spawn(self.n_streams)\n", + "\n", + " # create distributions\n", + "\n", + " # inter-arrival time distributions\n", + " self.arrival_strokes = Exponential(\n", + " self.iat_strokes, random_seed=self.seeds[0]\n", + " )\n", + "\n", + " self.acute_los = Lognormal(\n", + " self.acute_los_mean, self.acute_los_std, random_seed=self.seeds[1]\n", + " )\n", + "\n", + " if self.init_cond_params[\"mode\"] == \"fixed\":\n", + " self.init_conds = FixedDistribution(\n", + " self.init_cond_params[\"fixed\"]\n", + " )\n", + " else:\n", + " self.init_conds = DiscreteEmpirical(\n", + " values = self.init_cond_params[\"rnd\"][\"values\"],\n", + " freq = self.init_cond_params[\"rnd\"][\"freq\"],\n", + " random_seed=self.seeds[2]\n", + " )\n", + " \n", + "\n", + " def init_results_variables(self):\n", + " \"\"\"\n", + " Initialise all of the experiment variables used in results\n", + " collection. This method is called at the start of each run\n", + " of the model\n", + " \"\"\"\n", + " # variable used to store results of experiment\n", + " self.results = {}\n", + " self.results[\"n_arrivals\"] = 0\n", + " self.results[\"waiting_acute\"] = []" + ] + }, + { + "cell_type": "markdown", + "id": "7ff9beae-89cc-419c-b584-c05b81086865", + "metadata": {}, + "source": [ + "## 🥵 Warm-up period" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "dff74a08-37fd-4b18-8bcd-97994f38369a", + "metadata": {}, + "outputs": [], + "source": [ + "def warmup_complete(warm_up_period, env, args):\n", + " \"\"\"\n", + " End of warm-up period event. Used to reset results collection variables.\n", + "\n", + " Parameters:\n", + " ----------\n", + " warm_up_period: float\n", + " Duration of warm-up period in simultion time units\n", + "\n", + " env: simpy.Environment\n", + " The simpy environment\n", + "\n", + " args: Experiment\n", + " The simulation experiment that contains the results being collected.\n", + " \"\"\"\n", + " yield env.timeout(warm_up_period)\n", + " trace(f\"{env.now:.2f}: 🥵 Warm up complete.\")\n", + " \n", + " args.init_results_variables()" + ] + }, + { + "cell_type": "markdown", + "id": "94f0f9c5-22cb-493a-9f1f-4e2a8325beaa", + "metadata": {}, + "source": [ + "## 4. Pathway process logic\n", + "\n", + "The key things to recognise are \n", + "\n", + "* We include a optional parameter called `collection_results` that defaults to `True`. We may set this `False` in our functions that setup initial conditions" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "911528e1-e4eb-4307-bb26-632faf7769d1", + "metadata": {}, + "outputs": [], + "source": [ + "def acute_stroke_pathway(patient_id, env, args, collect_results=True):\n", + " \"\"\"Process a patient through the acute ward\n", + " Simpy generator function.\n", + " \n", + " Parameters:\n", + " -----------\n", + " patient_id: int\n", + " A unique id representing the patient in the process\n", + "\n", + " env: simpy.Environment\n", + " The simulation environment\n", + "\n", + " args: Experiment\n", + " Container class for the simulation parameters/results.\n", + " \"\"\"\n", + " arrival_time = env.now\n", + "\n", + " with args.acute_ward.request() as acute_bed_request:\n", + " yield acute_bed_request\n", + " \n", + " acute_admit_time = env.now\n", + " wait_for_acute = acute_admit_time - arrival_time\n", + "\n", + " if collect_results:\n", + " args.results['waiting_acute'].append(wait_for_acute)\n", + " \n", + " trace(f\"{env.now:.2f}: Patient {patient_id} admitted to acute ward.\" \\\n", + " + f\"(waited {wait_for_acute:.2f} days)\")\n", + " \n", + " # Simulate acute care treatment\n", + " acute_care_duration = args.acute_los.sample()\n", + " yield env.timeout(acute_care_duration)\n", + " \n", + " trace(f\"{env.now:.2f}: Patient {patient_id} discharged.\")" + ] + }, + { + "cell_type": "markdown", + "id": "de8990c2-a330-4c02-ac77-26c30d3e0a41", + "metadata": {}, + "source": [ + "## 4. Arrivals generator\n", + "\n", + "This is a standard arrivals generator. We create stroke arrivals according to their distribution." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "b3e686ce-5371-4471-a052-b9d43309bc85", + "metadata": {}, + "outputs": [], + "source": [ + "def stroke_arrivals_generator(env, args):\n", + " \"\"\"\n", + " Arrival process for strokes.\n", + "\n", + " Parameters:\n", + " ------\n", + " env: simpy.Environment\n", + " The simpy environment for the simulation\n", + "\n", + " args: Experiment\n", + " The settings and input parameters for the simulation.\n", + " \"\"\"\n", + " # use itertools as it provides an infinite loop\n", + " # with a counter variable that we can use for unique Ids\n", + " for patient_id in itertools.count(start=1):\n", + "\n", + " # the sample distribution is defined by the experiment.\n", + " inter_arrival_time = args.arrival_strokes.sample()\n", + " yield env.timeout(inter_arrival_time)\n", + "\n", + " args.results[\"n_arrivals\"] += 1\n", + " \n", + " trace(f\"{env.now:.2f}: Patient {patient_id}. Stroke arrival.\")\n", + "\n", + " # patient enters pathway\n", + " env.process(acute_stroke_pathway(patient_id, env, args))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "ca37c9cd-158c-416b-8490-b4c5c2e63412", + "metadata": {}, + "outputs": [], + "source": [ + "def setup_initial_conditions(\n", + " env: simpy.Environment, \n", + " args: Experiment,\n", + " collect_results: bool = False,\n", + "):\n", + " \"\"\"Set up initial conditions with patients already in the acute\n", + " stroke queue\n", + "\n", + " Parmaters:\n", + " ---------\n", + " env: simpy.Environment\n", + " The simpy environment for the simulation\n", + "\n", + " args: Experiment\n", + " The settings and input parameters for the simulation.\n", + "\n", + " collect_results: bool, optional (default = False)#\n", + " Should results be collected for initial conditions patients?\n", + " \n", + " \"\"\"\n", + " # sample the no. patients to load into queue\n", + " patients_to_load = args.init_conds.sample()\n", + "\n", + " for initial_id in range(1, patients_to_load+1):\n", + " # Create patients with negative IDs to distinguish them as init cond.\n", + " # we may or may not want collect results for initial conditions\n", + " env.process(acute_stroke_pathway(-initial_id, env, args, collect_results))\n", + " trace(f\"{env.now:.2f}: Patient {-initial_id} loaded into queue\")" + ] + }, + { + "cell_type": "markdown", + "id": "6058571e-9fdb-4961-be27-8a3b8c2fe26e", + "metadata": {}, + "source": [ + "## 5. Single run function" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "0d0ea6cf-7d95-4d2c-9690-fcdbdae35d84", + "metadata": {}, + "outputs": [], + "source": [ + "def single_run(\n", + " experiment, \n", + " rep=0,\n", + " wu_period=WU_PERIOD,\n", + " rc_period=RC_PERIOD\n", + "):\n", + " \"\"\"\n", + " Perform a single run of the model and return the results\n", + "\n", + " Parameters:\n", + " -----------\n", + "\n", + " experiment: Experiment\n", + " The experiment/paramaters to use with model\n", + "\n", + " rep: int\n", + " The replication number.\n", + "\n", + " wu_period: float, optional (default=WU_PERIOD)\n", + " Warm-up period\n", + "\n", + " rc_period: float, optional (default=RC_PERIOD)\n", + " The run length of the model\n", + " \"\"\"\n", + "\n", + " # reset all results variables to zero and empty\n", + " experiment.init_results_variables()\n", + "\n", + " # set random number set to the replication no.\n", + " # this controls sampling for the run.\n", + " experiment.set_random_no_set(rep)\n", + "\n", + " # environment is (re)created inside single run\n", + " env = simpy.Environment()\n", + "\n", + " # simpy resources\n", + " experiment.acute_ward = simpy.Resource(env, experiment.n_acute_beds)\n", + "\n", + " # load the acute stroke queue\n", + " setup_initial_conditions(env, experiment)\n", + "\n", + " # schedule a warm_up period\n", + " env.process(warmup_complete(wu_period, env, experiment))\n", + " \n", + " # we pass all arrival generators to simpy \n", + " env.process(stroke_arrivals_generator(env, experiment))\n", + "\n", + " # run model\n", + " env.run(until=wu_period+rc_period)\n", + "\n", + " # quick stats\n", + " results = {}\n", + " results['mean_acute_wait'] = np.array(\n", + " experiment.results[\"waiting_acute\"]\n", + " ).mean()\n", + "\n", + " # return single run results\n", + " return results" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "caf52390-5455-4fa1-bb22-60b5b91ad8d0", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.00: Patient -1 loaded into queue\n", + "0.00: Patient -2 loaded into queue\n", + "0.00: Patient -3 loaded into queue\n", + "0.00: Patient -4 loaded into queue\n", + "0.00: Patient -5 loaded into queue\n", + "0.00: Patient -6 loaded into queue\n", + "0.00: Patient -7 loaded into queue\n", + "0.00: Patient -8 loaded into queue\n", + "0.00: Patient -9 loaded into queue\n", + "0.00: Patient -10 loaded into queue\n", + "0.00: Patient -1 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -2 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -3 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -4 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -5 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -6 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -7 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -8 admitted to acute ward.(waited 0.00 days)\n", + "0.00: Patient -9 admitted to acute ward.(waited 0.00 days)\n", + "0.00: 🥵 Warm up complete.\n", + "2.74: Patient 1. Stroke arrival.\n", + "2.78: Patient 2. Stroke arrival.\n", + "3.36: Patient 3. Stroke arrival.\n", + "4.42: Patient 4. Stroke arrival.\n", + "4.68: Patient 5. Stroke arrival.\n", + "5.80: Patient -3 discharged.\n", + "5.80: Patient -10 admitted to acute ward.(waited 5.80 days)\n", + "5.80: Patient -8 discharged.\n", + "5.80: Patient 1 admitted to acute ward.(waited 3.06 days)\n", + "5.97: Patient 6. Stroke arrival.\n", + "6.08: Patient -6 discharged.\n", + "6.08: Patient 2 admitted to acute ward.(waited 3.30 days)\n", + "6.22: Patient -7 discharged.\n", + "6.22: Patient 3 admitted to acute ward.(waited 2.86 days)\n", + "6.33: Patient 7. Stroke arrival.\n", + "7.28: Patient 8. Stroke arrival.\n", + "7.41: Patient -4 discharged.\n", + "7.41: Patient 4 admitted to acute ward.(waited 2.99 days)\n", + "7.43: Patient 9. Stroke arrival.\n", + "7.55: Patient 10. Stroke arrival.\n", + "7.94: Patient -5 discharged.\n", + "7.94: Patient 5 admitted to acute ward.(waited 3.26 days)\n", + "8.11: Patient -2 discharged.\n", + "8.11: Patient 6 admitted to acute ward.(waited 2.14 days)\n", + "8.26: Patient 11. Stroke arrival.\n", + "8.42: Patient 12. Stroke arrival.\n", + "8.68: Patient 13. Stroke arrival.\n", + "8.72: Patient 14. Stroke arrival.\n", + "8.82: Patient -9 discharged.\n", + "8.82: Patient 7 admitted to acute ward.(waited 2.50 days)\n", + "8.95: Patient 15. Stroke arrival.\n", + "9.18: Patient 16. Stroke arrival.\n", + "9.87: Patient -1 discharged.\n", + "9.87: Patient 8 admitted to acute ward.(waited 2.58 days)\n", + "9.92: Patient 17. Stroke arrival.\n" + ] + }, + { + "data": { + "text/plain": [ + "{'mean_acute_wait': 2.8362718457918676}" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "TRACE = True\n", + "init_cond_params = INIT_COND_PARAMS.copy()\n", + "init_cond_params[\"mode\"] = \"fixed\"\n", + "\n", + "# uncomment to vary the fixed amount 10 = 1 in queue.\n", + "init_cond_params[\"fixed\"] = 10\n", + "\n", + "experiment = Experiment(init_cond_params=init_cond_params)\n", + "results = single_run(experiment, rep=1, wu_period=0.0, rc_period=10.0)\n", + "results" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "0aaef408-09ca-49e0-8d39-f4a088a4ef1b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'n_arrivals': 17,\n", + " 'waiting_acute': [3.0586175577655674,\n", + " 3.303449502025928,\n", + " 2.857579305401165,\n", + " 2.9908615153438225,\n", + " 3.2622390398417513,\n", + " 2.136685286636955,\n", + " 2.496306611009193,\n", + " 2.58443594831056]}" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "experiment.results" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/content/distributions.py b/content/distributions.py index 33d8ab5..fc2e943 100644 --- a/content/distributions.py +++ b/content/distributions.py @@ -3,9 +3,9 @@ """ import numpy as np -from typing import Optional, Union, Tuple +from typing import Optional, Union, Tuple, Any from numpy.random import SeedSequence -from numpy.typing import NDArray +from numpy.typing import NDArray, ArrayLike import math class Lognormal: @@ -184,4 +184,199 @@ def sample(self, size=None): ------- float or np.ndarray (if size >=1) """ - return self.rand.binomial(n=1, p=self.p, size=size) \ No newline at end of file + return self.rand.binomial(n=1, p=self.p, size=size) + +class DiscreteEmpirical: + """ + DiscreteEmpirical distribution implementation. + + A probability distribution that samples values with specified frequencies. + Useful for modeling categorical data or discrete outcomes with known + probabilities. + + Example uses: + ------------- + 1. Routing percentages + 2. Classes of entity + 3. Batch sizes of arrivals + 4. Initial conditions - no. entities in a queue. + """ + + def __init__( + self, + values: ArrayLike, + freq: ArrayLike, + random_seed: Optional[Union[int, SeedSequence]] = None, + ): + """ + Initialize a discrete distribution. + + Parameters + ---------- + values : ArrayLike + List of possible outcome values. Must be of equal length to freq. + + freq : ArrayLike + List of observed frequencies or probabilities. Must be of equal + length to values. These will be normalized to sum to 1. + + random_seed : Optional[Union[int, SeedSequence]], default=None + A random seed or SeedSequence to reproduce samples. If None, a + unique sample sequence is generated. + + Raises + ------ + TypeError + If values or freq are not positive arrays + ValueError + If values and freq have different lengths. + """ + + # convert to array first + self.values = np.asarray(values) + self.freq = np.asarray(freq) + + if len(self.values) != len(self.freq): + raise ValueError( + "values and freq arguments must be of equal length" + ) + + self.rng = np.random.default_rng(random_seed) + self.probabilities = self.freq / self.freq.sum() + + def __repr__(self): + values_repr = ( + str(self.values.tolist()) + if len(self.values) < 4 + else f"[{', '.join(str(x) for x in self.values[:3])}, ...]" + ) + freq_repr = ( + str(self.freq.tolist()) + if len(self.freq) < 4 + else f"[{', '.join(str(x) for x in self.freq[:3])}, ...]" + ) + return f"Discrete(values={values_repr}, freq={freq_repr})" + + def sample( + self, size: Optional[Union[int, Tuple[int, ...]]] = None + ) -> Union[Any, NDArray]: + """ + Generate random samples from the discrete distribution. + + Parameters + ---------- + size : Optional[Union[int, Tuple[int, ...]]], default=None + The number/shape of samples to generate: + - If None: returns a single sample + - If int: returns a 1-D array with that many samples + - If tuple of ints: returns an array with that shape + + Returns + ------- + Union[Any, NDArray] + Random samples from the discrete distribution: + - A single value (of whatever type was in the values array) when + size is None + - A numpy array of values with shape determined by size parameter + """ + sample = self.rng.choice(self.values, p=self.probabilities, size=size) + + if size is None: + return sample.item() + return sample + +class FixedDistribution: + """ + Fixed distribution implementation. + + A degenerate distribution that always returns the same fixed value. + Useful for constants or deterministic parameters in models. + + Provides a method to + sample a constant value regardless of the number of samples requested. + """ + + def __init__(self, value: float): + """ + Initialize a fixed distribution. + + Parameters + ---------- + value : float + The constant value that will be returned by sampling. + """ + self.value = value + + def __repr__(self): + return f"FixedDistribution(value={self.value})" + + def sample( + self, size: Optional[Union[int, Tuple[int, ...]]] = None + ) -> Union[float, NDArray[np.float64]]: + """ + Generate "samples" from the fixed distribution (always the same value). + + Parameters + ---------- + size : Optional[Union[int, Tuple[int, ...]]], default=None + The number/shape of samples to generate: + - If None: returns the fixed value as a float + - If int: returns a 1-D array filled with the fixed value + - If tuple of ints: returns an array with that shape filled with + the fixed value + + Returns + ------- + Union[float, NDArray[np.float64]] + The fixed value: + - A single float when size is None + - A numpy array filled with the fixed value with shape + determined by size parameter + """ + if size is not None: + return np.full(size, self.value) + return self.value + +class Bernoulli: + """ + Convenience class for the Bernoulli distribution. + packages up distribution parameters, seed and random generator. + + The Bernoulli distribution is a special case of the binomial distribution + where a single trial is conducted + + Use the Bernoulli distribution to sample success or failure. + """ + + def __init__(self, p, random_seed=None): + """ + Constructor + + Params: + ------ + p: float + probability of drawing a 1 + + random_seed: int | SeedSequence, optional (default=None) + A random seed to reproduce samples. If set to none then a unique + sample is created. + """ + self.rand = np.random.default_rng(seed=random_seed) + self.p = p + + def sample(self, size=None): + """ + Generate a sample from the exponential distribution + + Params: + ------- + size: int, optional (default=None) + the number of samples to return. If size=None then a single + sample is returned. + + Returns: + ------- + float or np.ndarray (if size >=1) + """ + return self.rand.binomial(n=1, p=self.p, size=size) + diff --git a/content/img/acute_stroke_pathway.png b/content/img/acute_stroke_pathway.png new file mode 100644 index 0000000..5a3ef95 Binary files /dev/null and b/content/img/acute_stroke_pathway.png differ