diff --git a/CHANGELOG.md b/CHANGELOG.md index 9bc927a..fc7fa8c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * `12_arrival_classes.ipynb`: simulate unique processes for different classes of arrival to the model. * `distributions.py`: module containing some distributions to reduce code in notebooks. * `basic_model.py`: contains a single activity version of the call centre model to use with `07_exercise.ipynb` +* `13_warm_up.ipynb`: contains explanation of implementing warm up +* `14_initial_conditions.ipynb`: contains explanation of manually setting up processes before a model is run +* `15_resource_store.ipynb`: introduction to `Store` and `FilterStore` for advanced resource modelling +* `sim_utility.py`: added `trace`, `set_trace`, and `spawn_seeds` functions to use across notebooks. ## [v0.2.0 - 11/02/2024](https://github.com/pythonhealthdatascience/intro-open-sim/releases/tag/v0.2.0) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.14849934.svg)](https://doi.org/10.5281/zenodo.14849934) diff --git a/content/10_multiple_arrival_processes.ipynb b/content/10_multiple_arrival_processes.ipynb index 522c90e..2c9560f 100644 --- a/content/10_multiple_arrival_processes.ipynb +++ b/content/10_multiple_arrival_processes.ipynb @@ -11,7 +11,8 @@ "\n", "We will work with a hypothetical hospital that provides emergency orthopedic surgery to different classes of patient.\n", "\n", - "![model image](img/multiple_arrivals.png \"Urgent care call centre\")\n", + "\"model\n", + "\n", "\n", "| ID | Arrival Type | Distribution | Mean (mins) |\n", "|----|-----------------|--------------|-------------|\n", @@ -811,7 +812,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.11" + "version": "3.11.14" } }, "nbformat": 4, diff --git a/content/13_warm_up.ipynb b/content/13_warm_up.ipynb index ddd6971..dacf063 100644 --- a/content/13_warm_up.ipynb +++ b/content/13_warm_up.ipynb @@ -514,7 +514,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.11" + "version": "3.11.15" } }, "nbformat": 4, diff --git a/content/14_initial_conditions.ipynb b/content/14_initial_conditions.ipynb index 4cac409..b760af2 100644 --- a/content/14_initial_conditions.ipynb +++ b/content/14_initial_conditions.ipynb @@ -588,7 +588,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.11" + "version": "3.11.15" } }, "nbformat": 4, diff --git a/content/15_resource_stores.ipynb b/content/15_resource_stores.ipynb new file mode 100644 index 0000000..184aeaf --- /dev/null +++ b/content/15_resource_stores.ipynb @@ -0,0 +1,1183 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4da08d47-1de7-4634-a8fc-fc05657475f7", + "metadata": {}, + "source": [ + "# Advanced Resources\n", + "\n", + "## Limitations of `Resource`\n", + "\n", + "Up to this point the models we have used a basic `simpy.Resource`. This has been very useful to model queues, but one potential downside to `Resource` is that it *does not allow you to model individual resource attributes or any type of complex behaviour*. For many models in healthcare, this is sufficient, but there may be instances where you need to track and control individual resources. For example, ambulances, or different types of staff. In this notebook we will explore how to add more complex behaviour using the `Store` and `FilterStore` objects provided by `simpy`.\n", + "\n", + "🎓 The good news is that both `Store` and `FilterStore` are easy to use. Usually this is within the context of a complex simulation, but we will keep our models simple here and focus on how to use them.\n", + "\n", + "## When to use each SimPy resource type\n", + "\n", + "| Type | Best used for | Key idea |\n", + "|---|---|---|\n", + "| `Resource` | Identical servers | Only capacity matters. |\n", + "| `Store` | Individual objects | Retrieve the next available object. |\n", + "| `FilterStore` | Individual objects with selection rules | Retrieve an object matching a condition. |\n" + ] + }, + { + "cell_type": "markdown", + "id": "8d4bd324-bdd2-4b60-ac94-862446c29aca", + "metadata": {}, + "source": [ + "## 1. Imports" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "bcedfd8f-49ac-4e0c-bf11-2e4d244db365", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import itertools\n", + "import simpy\n", + "import math" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "f1a26662-34e6-4ac7-b186-00c013183f16", + "metadata": {}, + "outputs": [], + "source": [ + "# to reduce code these classes can be found in distribution.py\n", + "# you can also pip `install sim-tools`\n", + "from distributions import (\n", + " Exponential, \n", + " DiscreteEmpirical\n", + ")\n", + "\n", + "from sim_utility import set_trace, trace, spawn_seeds" + ] + }, + { + "cell_type": "markdown", + "id": "b03d45cf-d0c9-4370-8922-561ea89e2f02", + "metadata": {}, + "source": [ + "## 2. PART 1: Using a `simpy.Store`\n", + "\n", + "A `Store` can be thought of as a container that holds Python objects. We can `put()` objects into the store and later `get()` them back out again, usually on a First In First Out (FIFO) basis.\n", + "\n", + "This is useful when resources are not all identical. With a normal `simpy.Resource`, we know how many units are available, but we do not know which specific unit is being used. With a `Store`, we can model individual resources with their own attributes.\n", + "\n", + "We create a `Store` as follows (`env` is a `simpy.Environment`):\n", + "\n", + "```python\n", + "store = simpy.Store(env, capacity=2)\n", + "```\n", + "We have to `put` something in the store. Let's take a very simple example where people experiencing a medical emergency are assigned either a Rapid Response Vehicle or a \"normal\" Type 1 Ambulance. A process flow describing use of the `Store` is below.\n", + "\n", + "\"model\n", + "\n", + "This is a two part process. We first create the instances of ambulances and then we call `Store.put()`\n", + "\n", + "\n", + "```python\n", + "ambulances = [\n", + " Ambulance(ambulance_id=1, vehicle_type=\"rrv\"), \n", + " Ambulance(ambulance_id=2, vehicle_type=\"type_1\")\n", + "]\n", + "\n", + "# loop through list of ambulances and put them in the store object\n", + "for amb in ambulances:\n", + " store.put(amb)\n", + "\n", + "```\n", + "\n", + "To get an an ambulance we used the `get()` method along with the `yield` keyword.\n", + "\n", + "```python\n", + "ambulance = yield store.get()\n", + "```\n", + "\n", + "Like normal `Resource` objects the use of `yield` means that we can simulate queues when no `Ambulance` objects are left in the store.\n", + "\n", + "The code below implements the full model. In the model we will track\n", + "\n", + "* The utilisation and number of emergencies handled by each individual vehicle.\n", + "* Overall patient waiting time and associated statistics.\n", + " \n", + "### 2.1 Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "d5b3eb67-9c8a-4fd8-9312-d337a65c52ee", + "metadata": {}, + "outputs": [], + "source": [ + "NUM_AMBULANCES = 10\n", + "FLEET = {\n", + " \"rrv\": 4,\n", + " \"type_1\": 6,\n", + "}\n", + "\n", + "RRV_SERVICE_TIME = 50.0 # minutes\n", + "TYPE1_SERVICE_TIME = 65.0 \n", + "RUN_LENGTH = 1_000 # minutes\n", + "RANDOM_SEED = 42\n", + "\n", + "# exponential IAT \n", + "MEAN_INTERARRIVAL = 6" + ] + }, + { + "cell_type": "markdown", + "id": "03114c1d-faf6-4c95-b7d2-64d6a17003a8", + "metadata": {}, + "source": [ + "### 2.2 Entity classes\n", + "\n", + "We will model Ambulances using a simple Python class called `Ambulance`. For simplicity we will give the ambulance an attribute called `vehicle_type`. We will use that to look up the service time distribution from a Python dictionary parameter." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "7956189c-c777-4c0a-afcc-ed06ec24ee3e", + "metadata": {}, + "outputs": [], + "source": [ + "class Ambulance:\n", + " \"\"\"\n", + " An ambulance resource\n", + "\n", + " Parameters:\n", + " ----------\n", + " ambulance_id: int\n", + " Unique id of the ambulance\n", + " vehicle_type: str\n", + " Set to either \"rrv\" or \"type_1\"\n", + " \"\"\"\n", + " def __init__(self, ambulance_id: int, vehicle_type: str):\n", + " self.ambulance_id = ambulance_id\n", + "\n", + " # vehicle type can be \"rrv\" or \"type_1\"\n", + " self.vehicle_type = vehicle_type\n", + " \n", + " self.total_jobs = 0\n", + " # cumulative busy time\n", + " self.total_busy = 0.0 \n", + "\n", + " self.emojis = {\n", + " \"rrv\": \"🏎️\",\n", + " \"type_1\": \"🚑\"\n", + " }\n", + " \n", + " self.emoji = self.emojis[self.vehicle_type]\n", + " \n", + " def __repr__(self):\n", + " \"\"\"\n", + " A text representation of the ambulance to help debugging\n", + " \"\"\"\n", + " return f\"Ambulance({self.ambulance_id}, {self.emoji})\"" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "fa690897-a9e5-473b-8f01-35991b1b4213", + "metadata": {}, + "outputs": [], + "source": [ + "class Patient:\n", + " \"\"\"\n", + " A class to hold patient attributes\n", + "\n", + " Parameters\n", + " ----------\n", + " patient_id: int\n", + " Unique patient id\n", + "\n", + " arrival_time: float\n", + " Time of arrival to the simulation \n", + " \"\"\"\n", + " def __init__(self, patient_id: int, arrival_time: float):\n", + " self.patient_id = patient_id\n", + " self.arrival_time = arrival_time\n", + "\n", + " def __repr__(self):\n", + " \"\"\"\n", + " A text representation of the Patient for debug\n", + " \"\"\"\n", + " return f\"Patient({self.patient_id})\"" + ] + }, + { + "cell_type": "markdown", + "id": "bc5b1808-3f8e-43e2-b395-151de78bb718", + "metadata": {}, + "source": [ + "## 2.3 Ambulance dispatch and service process\n", + "\n", + "The code looks a little different from when using a standard `Resource` because we do not use a `with` context manager when using a `Store`. This means we need to explicitly \"release\" an `Ambulance` by putting it back into the `Store`,\n", + "\n", + "The function `dispatch_ambulance` is a `simpy` process that accepts a `dict` called `dists`. It contains a service time distribution for each type of ambulance i.e. \"rrv\" and \"type_1\". For example if we had a rapid response vehicle we could sample a service time as follows:\n", + "\n", + "For example:\n", + "\n", + "```python\n", + "ambulance = Ambulance(1, \"rrv\")\n", + "\n", + "# vehicle type is a str with value \"rrv\" It returns a Exponential(50)\n", + "service_dist = dists[ambulance.vehicle_type]\n", + "service_time = service_dist.sample()\n", + "\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "08e53cbb-66dc-4e75-bcab-5ff1c6bf2a8b", + "metadata": {}, + "outputs": [], + "source": [ + "def dispatch_ambulance(\n", + " env: simpy.Environment,\n", + " store: simpy.Store,\n", + " patient: Patient,\n", + " dists: dict,\n", + " log: dict\n", + ") -> None:\n", + " \"\"\"Simulate ambulance dispatch and service: \n", + " \n", + " 1. queues a patient up for the next free Ambulance (FIFO), \n", + " 2. simulates service (as a single distribution), \n", + " 3. returns the Ambulance to the store.\n", + "\n", + " Parameters:\n", + " ----------\n", + " env: simpy.Environment\n", + " The simpy environment for the simulation\n", + " store: simpy.Store\n", + " A store of Ambulance objects\n", + " dists: dict\n", + " Contains the \"service\" distribution\n", + " log: dict\n", + " Audit dictionary\n", + " \"\"\"\n", + "\n", + " # Wait for an available Ambulance\n", + " # note we `get()` an Ambulance from the store\n", + " ambulance: Ambulance = yield store.get()\n", + "\n", + " wait_time = env.now - patient.arrival_time\n", + " log[\"wait_times\"].append(wait_time)\n", + "\n", + " # Service varies by vehicle type (travel + on-scene + return)\n", + " service_dist = dists[ambulance.vehicle_type]\n", + " service_time = service_dist.sample()\n", + "\n", + " # debug\n", + " trace(\n", + " f\"{env.now:.1f}: {ambulance.emoji} {ambulance.ambulance_id} → {patient} \"\n", + " f\"(waited {wait_time:.1f} min, service {service_time:.1f} min)\"\n", + " )\n", + "\n", + " yield env.timeout(service_time)\n", + "\n", + " # Update ambulance stats and return (put) to store\n", + " ambulance.total_jobs += 1\n", + " ambulance.total_busy += service_time\n", + " store.put(ambulance)\n", + "\n", + " log[\"service_times\"].append(service_time)\n", + " log[\"assignments\"].append(\n", + " (patient.patient_id, ambulance.ambulance_id, ambulance.vehicle_type)\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "af2955b7-fc92-474d-8238-6b3f78b47f0c", + "metadata": {}, + "source": [ + "### 2.4 Patient arrival generator\n", + "\n", + "Our patient generator process follows the simple template approach we have used before: an infinite loop where we sample the time until the next arrival and then create and schedule new `dispatch_ambulance` process." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "92544964-085b-4309-b138-0b7c8e9a1c14", + "metadata": {}, + "outputs": [], + "source": [ + "def patient_arrivals_generator(\n", + " env: simpy.Environment,\n", + " store: simpy.Store,\n", + " dists: dict,\n", + " log: dict\n", + ") -> None:\n", + " \"\"\"\n", + " Arrival process for patients to the ambulance sim.\n", + "\n", + " Parameters:\n", + " ------\n", + " env: simpy.Environment\n", + " The simpy environment for the simulation\n", + "\n", + " store: simpy.Store\n", + " A store of Ambulance objects\n", + "\n", + " dists: dict\n", + " Contains \"arrival\" and \"service\" distributions\n", + "\n", + " log: dict\n", + " Results dictionary\n", + " \"\"\"\n", + " for patient_id in itertools.count(start=1):\n", + "\n", + " # time until next patient arrival\n", + " inter_arrival_time = dists[\"arrival\"].sample()\n", + " yield env.timeout(inter_arrival_time)\n", + "\n", + " log[\"n_arrivals\"] += 1\n", + " patient = Patient(patient_id, env.now)\n", + "\n", + " # debug info\n", + " trace(f\"{env.now:.1f}: 📞 {patient}\")\n", + "\n", + " # create ambulance dispatch + service process\n", + " env.process(dispatch_ambulance(env, store, patient, dists, log))" + ] + }, + { + "cell_type": "markdown", + "id": "9dcc1021-dbdb-4c39-a847-f263f09f6d54", + "metadata": {}, + "source": [ + "### 2.5 Single run function\n", + "\n", + "The `single_run` function creates our dictionary of distributions, ambulances, `Store` resource and runs the simulation returning a log of results and the ambulances that we will use the pass to a simple function called `results_summary` that will print out some formatted statistics." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "a5d11438-61b0-4345-9ecc-83db0db32d17", + "metadata": {}, + "outputs": [], + "source": [ + "def single_run(\n", + " mean_iat: float = MEAN_INTERARRIVAL,\n", + " mean_rrv_service: float = RRV_SERVICE_TIME,\n", + " mean_type_1_service: float = TYPE1_SERVICE_TIME,\n", + " n_ambulances: int = NUM_AMBULANCES,\n", + " fleet: dict = FLEET,\n", + " run_length: float = RUN_LENGTH, \n", + " random_seed: int = 1\n", + "):\n", + " \"\"\"\n", + " Set up and perform a single replication of the MMS model\n", + " \"\"\"\n", + "\n", + " # generate 3 rng seeds\n", + " seeds = spawn_seeds(n_streams=3, main_seed=random_seed)\n", + " \n", + " # 1. distribution objects\n", + " dists = {\n", + " \"arrival\": Exponential(mean_iat, random_seed=seeds[0]),\n", + " \"rrv\": Exponential(mean_rrv_service, random_seed=seeds[1]),\n", + " \"type_1\": Exponential(mean_type_1_service, random_seed=seeds[2]),\n", + " }\n", + "\n", + " # 2. simpy environment \n", + " env = simpy.Environment()\n", + "\n", + " # 3. Initialise Store\n", + " # 3.1 Create empty Store with sufficient slots\n", + " store = simpy.Store(env, capacity=n_ambulances)\n", + "\n", + " # 3.2 Create Ambulance objects\n", + " ambulances = []\n", + " ambulance_id = 0\n", + " \n", + " for vehicle_type, count in fleet.items():\n", + " for _ in range(count):\n", + " ambulance_id += 1\n", + " ambulances.append(Ambulance(ambulance_id, vehicle_type))\n", + " \n", + " # 3.3 `put` Ambulance objects into the store\n", + " for amb in ambulances:\n", + " store.put(amb)\n", + "\n", + " # 4. results dictionary\n", + " log = {\"n_arrivals\": 0, \"wait_times\": [], \"service_times\": [], \"assignments\": []}\n", + "\n", + " env.process(patient_arrivals_generator(env, store, dists, log))\n", + " env.run(until=run_length)\n", + "\n", + " return ambulances, log" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "2b32b8dc-59cb-4680-afa2-064ee5bd3c99", + "metadata": {}, + "outputs": [], + "source": [ + "def results_summary(\n", + " log: dict,\n", + " ambulances: list[Ambulance],\n", + " run_length: float,\n", + " mean_iat: float\n", + "):\n", + " waits = np.array(log[\"wait_times\"])\n", + " n_served = len(waits)\n", + "\n", + " print(\"\\n\" + \"═\" * 55)\n", + " print(f\" Mean inter-arrival : {mean_iat:.2f} min\")\n", + " print(\"─\" * 55)\n", + " print(f\" Patients arrived : {log['n_arrivals']}\")\n", + " print(f\" Patients served : {n_served}\")\n", + "\n", + " if n_served > 0:\n", + " print(f\" Mean wait time : {waits.mean():.2f} min\")\n", + " print(f\" P(wait > 0) : {(waits > 0).mean():.2%}\")\n", + " print(f\" 95th pct wait : {np.percentile(waits, 95):.2f} min\")\n", + " else:\n", + " print(\" No patients were served.\")\n", + "\n", + " print(\"─\" * 55)\n", + " print(f\" {'Ambulance':<15} {'Jobs':>6} {'Utilisation':>12}\")\n", + " print(\"─\" * 55)\n", + "\n", + " for amb in ambulances:\n", + " util = amb.total_busy / run_length\n", + " print(f\" {amb.emoji} {amb.ambulance_id:<5} {amb.total_jobs:>6} {util:>11.2%}\")\n", + "\n", + " print(\"═\" * 55)" + ] + }, + { + "cell_type": "markdown", + "id": "3c663f15-1df1-4652-b202-90378d4e81fb", + "metadata": {}, + "source": [ + "### 2.6 Run the model and view results \n", + "\n", + "Use `set_trace` to toggle the printing of the model event debug on and off. " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "0506ba12-8039-4d13-a3a1-4d7e08efa44f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Simulation tracing set to: False\n", + "\n", + "═══════════════════════════════════════════════════════\n", + " Mean inter-arrival : 6.00 min\n", + "───────────────────────────────────────────────────────\n", + " Patients arrived : 168\n", + " Patients served : 157\n", + " Mean wait time : 38.29 min\n", + " P(wait > 0) : 78.98%\n", + " 95th pct wait : 90.39 min\n", + "───────────────────────────────────────────────────────\n", + " Ambulance Jobs Utilisation\n", + "───────────────────────────────────────────────────────\n", + " 🏎️ 1 19 79.09%\n", + " 🏎️ 2 18 86.62%\n", + " 🏎️ 3 13 85.86%\n", + " 🏎️ 4 9 79.45%\n", + " 🚑 5 10 82.41%\n", + " 🚑 6 18 81.24%\n", + " 🚑 7 15 87.30%\n", + " 🚑 8 17 90.95%\n", + " 🚑 9 15 89.03%\n", + " 🚑 10 13 86.19%\n", + "═══════════════════════════════════════════════════════\n" + ] + } + ], + "source": [ + "set_trace(False)\n", + "\n", + "ambulances, log = single_run(run_length=RUN_LENGTH, random_seed=42)\n", + "results_summary(log, ambulances, run_length=RUN_LENGTH, mean_iat=MEAN_INTERARRIVAL)" + ] + }, + { + "cell_type": "markdown", + "id": "17ba6461-a625-42d4-b4cf-13fd546c7df5", + "metadata": {}, + "source": [ + "## 3. PART 2: Using a `FilterStore`\n", + "\n", + "In this example we will modify the ambulance dispatch simulation so that patients and ambulances are located within one of five **nodes** that are specified by coordinates.\n", + "\n", + "| ID | Label | Coordinates |\n", + "|---|---|---|\n", + "| 0 | South-West | `(0, 0)` |\n", + "| 1 | South-East | `(4, 0)` |\n", + "| 2 | North-West | `(0, 4)` |\n", + "| 3 | North-East | `(4, 4)` |\n", + "| 4 | Centre | `(2, 2)` |\n", + "| 5 | Hospital | `(2, 0)` |\n", + "\n", + "\n", + "Patients in need will be assigned the closest available ambulance. This is why we use a `FilterStore`. We iterate through items in the store to locate the closest ambulance.\n", + "\n", + "An assigned ambulance has to travel to the patient (constant speed), pick them up (Exponential), travel to the hospital, and then travel back to their home node.\n", + "\n", + "If no ambulances are available they are assigned the next available ambulance when it has returned to its home node." + ] + }, + { + "cell_type": "markdown", + "id": "468a41d0-42e9-4fd1-813c-76b7457452b0", + "metadata": {}, + "source": [ + "### 3.1 Geographic information" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "0eabc7f2-bb0e-4368-a245-2a123c0ba420", + "metadata": {}, + "outputs": [], + "source": [ + "# All locations share a single distance matrix: nodes 0-4 plus hospital (5)\n", + "HOSPITAL_ID = 5\n", + "\n", + "LOCATIONS = {\n", + " 0: (0.0, 0.0), # South-West\n", + " 1: (4.0, 0.0), # South-East\n", + " 2: (0.0, 4.0), # North-West\n", + " 3: (4.0, 4.0), # North-East\n", + " 4: (2.0, 2.0), # Centre\n", + " HOSPITAL_ID: (2.0, 0.0), # Hospital (south-centre)\n", + "}\n", + "\n", + "# nodes patients can appear in\n", + "PATIENT_NODES = [0, 1, 2, 3, 4] \n", + "NODE_PROBS = [0.15, 0.25, 0.20, 0.20, 0.20]\n", + "\n", + "# Pre-compute ALL pairwise distances before simulation\n", + "DISTANCES = {\n", + " (i, j): math.dist(LOCATIONS[i], LOCATIONS[j])\n", + " for i in LOCATIONS\n", + " for j in LOCATIONS\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "88223687-d608-457b-8464-542c020c594d", + "metadata": {}, + "source": [ + "### 3.2 Parameters" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "1079ec2e-1d51-447d-99b1-11d9f081e36e", + "metadata": {}, + "outputs": [], + "source": [ + "NUM_AMBULANCES = 10\n", + "TRAVEL_SPEED = 0.25 # distance units per minute\n", + "MEAN_SCENE_TIME = 20.0 \n", + "MEAN_INTERARRIVAL = 6 \n", + "RUN_LENGTH = 1_000\n", + "RANDOM_SEED = 42\n", + "\n", + "# Ambulance positioning parameter\n", + "# 2 ambulances stationed at each of the 5 nodes\n", + "AMBULANCE_HOME_NODES = [\n", + " node\n", + " for node in PATIENT_NODES\n", + " for _ in range(NUM_AMBULANCES // len(PATIENT_NODES))\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "16800dd7-7f83-4c42-bde6-0473fa15a1de", + "metadata": {}, + "source": [ + "### 3.3 Entity Classes\n", + "\n", + "The `Ambulance` and `Patient` classes have been slightly modified to include a `node` attribute. For simplicity all ambulances have the same type" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "723855d2-1389-46cc-9a7f-663369dbd686", + "metadata": {}, + "outputs": [], + "source": [ + "class Ambulance:\n", + " def __init__(self, ambulance_id: int, home_node: int):\n", + " self.ambulance_id = ambulance_id\n", + " # modification = a home node or 'base' for the ambulance\n", + " self.home_node = home_node\n", + " self.total_jobs = 0\n", + " self.total_busy = 0.0\n", + "\n", + " def __repr__(self):\n", + " return f\"Ambulance(id={self.ambulance_id}, node={self.home_node})\"" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "e7031697-436e-4b85-bbde-91bde6ed3806", + "metadata": {}, + "outputs": [], + "source": [ + "class Patient:\n", + " def __init__(self, patient_id: int, arrival_time: float, node: int):\n", + " self.patient_id = patient_id\n", + " self.arrival_time = arrival_time\n", + " self.node = node\n", + "\n", + " def __repr__(self):\n", + " return f\"Patient(id={self.patient_id}, node={self.node})\"" + ] + }, + { + "cell_type": "markdown", + "id": "5cd315ba-0cf9-45bd-bee0-9c6796159622", + "metadata": {}, + "source": [ + "### 3.4 Travel functions" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "7a6c9fe6-4c86-4b6c-a758-2dba9fa622de", + "metadata": {}, + "outputs": [], + "source": [ + "def travel_time(loc_a: int, loc_b: int) -> float:\n", + " \"\"\"Minutes to travel between any two location IDs.\"\"\"\n", + " return DISTANCES[loc_a, loc_b] / TRAVEL_SPEED" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "17c05a56-2f6a-4ffa-b65b-02d8ecc94671", + "metadata": {}, + "outputs": [], + "source": [ + "def closest_ambulance_id(\n", + " available: list[Ambulance],\n", + " patient_node: int\n", + ") -> int | None:\n", + " \"\"\"\n", + " Finds the closest ambulance and return the ID.\n", + " Note this code is not designed for efficiency. \n", + " \"\"\"\n", + "\n", + " # if list of ambulance is empty\n", + " if not available:\n", + " return None\n", + "\n", + " # return the closest ambulance id\n", + " return min(available, key=lambda a: DISTANCES[a.home_node, patient_node]).ambulance_id" + ] + }, + { + "cell_type": "markdown", + "id": "9d7f1394-e263-4f61-ab8b-f27c22ffc86f", + "metadata": {}, + "source": [ + "### 3.5 Modified simpy processes\n", + "\n", + "The biggest update to our ambulance dispatch function is that it will now use a `FilterStore` to take account of `Patient` and `Ambulance` node location.\n", + "\n", + "Our code also needs to be able to handle a situation where all ambulances are in service. In these instances a `Patient` would just wait for the next available ambulance on a FIFO basis. \n", + "\n", + "Examples. Let's assume ambulance 2 is the closest to the patient. We `get` that ambulance from the filter store like so\n", + "\n", + "```python\n", + "best_id = 2\n", + "ambulance = yield store.get(lambda a: a.ambulance_id == best_id)\n", + "```\n", + "\n", + "In an instance where there are no ambulances available we can force the `FilterStore` to work on a FIFO basis like so\n", + "\n", + "```python\n", + "ambulance = yield store.get(lambda a: True)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "383ad331-bba1-4571-9386-28d60d16635e", + "metadata": {}, + "outputs": [], + "source": [ + "def dispatch_ambulance(env, store, patient, dists, log):\n", + " \"\"\"Modified ambulance dispatch process\n", + " \"\"\"\n", + " \n", + " # find the closest ambulance in the store\n", + " # note we pass store.items which is a list of Ambulance objects\n", + " # it may be empty! This means all ambulances are in use.\n", + " best_id = closest_ambulance_id(store.items, patient.node)\n", + "\n", + " # filter store code\n", + " if best_id is not None:\n", + " # if an ambulance is available get the closest ambulance from the store\n", + " ambulance = yield store.get(lambda a: a.ambulance_id == best_id)\n", + " else:\n", + " # otherwise just wait for the next available ambulance. FIFO\n", + " # we are forcing a FilterStore to behave like a standard Store\n", + " ambulance = yield store.get(lambda a: True)\n", + "\n", + " wait_time = env.now - patient.arrival_time\n", + "\n", + " # Leg 1: travel from home node to patient \n", + " t_to_patient = travel_time(ambulance.home_node, patient.node)\n", + "\n", + " # Leg 2: on scene\n", + " scene_time = dists[\"on_scene\"].sample()\n", + "\n", + " # Leg 3: transport patient to hospital\n", + " t_to_hospital = travel_time(patient.node, HOSPITAL_ID)\n", + "\n", + " # Leg 4: return to home base \n", + " t_to_base = travel_time(HOSPITAL_ID, ambulance.home_node)\n", + "\n", + " # total turnaround time\n", + " total_service = t_to_patient + scene_time + t_to_hospital + t_to_base\n", + "\n", + " # debug\n", + " trace(\n", + " f\"{env.now:.1f}: 🚑 {ambulance.ambulance_id} → {patient} \"\n", + " f\"(waited {wait_time:.1f} min, service {total_service:.1f} min)\"\n", + " )\n", + "\n", + " yield env.timeout(total_service)\n", + " \n", + " ambulance.total_jobs += 1\n", + " ambulance.total_busy += total_service\n", + " store.put(ambulance)\n", + "\n", + " log[\"wait_times\"].append(wait_time)\n", + "\n", + " # create row in the assignments table\n", + " log[\"assignments\"].append(\n", + " dict(patient_id=patient.patient_id,\n", + " patient_node=patient.node,\n", + " ambulance_id=ambulance.ambulance_id,\n", + " ambulance_node=ambulance.home_node,\n", + " dispatch_dist=DISTANCES[ambulance.home_node, patient.node],\n", + " t_to_patient=t_to_patient,\n", + " scene_time=scene_time,\n", + " t_to_hospital=t_to_hospital,\n", + " t_to_base=t_to_base,\n", + " total_service=total_service,\n", + " immediate_dispatch=best_id is not None,\n", + " wait=wait_time)\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "4155f0a2-4ee7-4c69-954d-a6aee59a0a12", + "metadata": {}, + "outputs": [], + "source": [ + "def patient_arrivals_generator(\n", + " env: simpy.Environment,\n", + " store: simpy.Store,\n", + " dists: dict,\n", + " log: dict\n", + ") -> None: \n", + " \"\"\"Modified Arrival process for patients to the ambulance sim\n", + " Patient now includes arrival node\"\"\"\n", + " for patient_id in itertools.count(start=1):\n", + "\n", + " # time until next patient arrival\n", + " inter_arrival_time = dists[\"arrival\"].sample()\n", + " yield env.timeout(inter_arrival_time)\n", + "\n", + " log[\"n_arrivals\"] += 1\n", + " node = dists[\"arrival_node\"].sample()\n", + " patient = Patient(patient_id, env.now, node)\n", + "\n", + " # debug info\n", + " trace(f\"{env.now:.1f}: 📞 {patient}\")\n", + "\n", + " # create ambulance dispatch + service process\n", + " env.process(dispatch_ambulance(env, store, patient, dists, log))" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "df33aaf6-df42-4fb8-918f-fefae73fc9eb", + "metadata": {}, + "outputs": [], + "source": [ + "def single_run(\n", + " mean_iat: float = MEAN_INTERARRIVAL,\n", + " mean_on_scene: float = MEAN_SCENE_TIME,\n", + " n_ambulances: int = NUM_AMBULANCES,\n", + " run_length: float = RUN_LENGTH, \n", + " random_seed: int = 42\n", + "):\n", + " \"\"\"\n", + " Set up and perform a single replication of the MMS model\n", + " \"\"\"\n", + "\n", + " # generate 3 rng seeds\n", + " seeds = spawn_seeds(n_streams=3, main_seed=random_seed)\n", + " \n", + " # 1. distribution objects\n", + " # We have included a new arrival_node an on_scene dists\n", + " dists = {\n", + " \"arrival\": Exponential(mean_iat, random_seed=seeds[0]),\n", + " \"arrival_node\": DiscreteEmpirical(\n", + " values=[0, 1, 2, 3, 4], \n", + " freq=[p * 100 for p in NODE_PROBS], # need to freqs not probs\n", + " random_seed=seeds[1]),\n", + " \"on_scene\": Exponential(mean_on_scene, random_seed=seeds[2]),\n", + " }\n", + "\n", + " # 2. simpy environment \n", + " env = simpy.Environment()\n", + "\n", + " # 3. Initialise Store\n", + " # 3.1 Create empty Store with sufficient slots\n", + " store = simpy.FilterStore(env, capacity=n_ambulances)\n", + "\n", + " # 3.2 Create Ambulance objects now with home nodes\n", + " ambulances = [Ambulance(i + 1, home) for i, home in enumerate(AMBULANCE_HOME_NODES)]\n", + "\n", + " # 3.3 `put` Ambulance objects into the filter store\n", + " for amb in ambulances:\n", + " store.put(amb)\n", + "\n", + " # 4. results dictionary\n", + " log = {\"n_arrivals\": 0, \"wait_times\": [], \"service_times\": [], \"assignments\": []}\n", + "\n", + " env.process(patient_arrivals_generator(env, store, dists, log))\n", + " env.run(until=run_length)\n", + "\n", + " return ambulances, log" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "1aa60d11-3dfb-4d2a-9bd1-b75259f9aab5", + "metadata": {}, + "outputs": [], + "source": [ + "def single_run_summary(assignments, ambulances, run_length):\n", + " waits = assignments[\"wait\"].to_numpy()\n", + " \n", + " # immediate versus delayed dispatches of Ambulance\n", + " immediate = assignments[assignments[\"immediate_dispatch\"]]\n", + " delayed = assignments[~assignments[\"immediate_dispatch\"]]\n", + " \n", + " print(\"\\n\" + \"═\" * 65)\n", + " print(f\" Patients served : {len(assignments)}\")\n", + " print(f\" Immediate Dispatch : {len(immediate)} ({len(immediate)/len(assignments):.1%})\")\n", + " print(f\" Delayed Dispatch : {len(delayed)} ({len(delayed)/len(assignments):.1%})\")\n", + " print(f\" Overall mean wait (min) : {waits.mean():.2f}\")\n", + " print(\"─\" * 65)\n", + " print(f\" {'':<30} {'Ambulance Dispatch':>20}\")\n", + " print(f\" {'Metric':<30} {'Immediate':>10} {'Delayed':>10}\")\n", + " print(\"─\" * 65)\n", + " for metric, col in [(\"Mean wait (min)\", \"wait\"),\n", + " (\"Mean dispatch dist\", \"dispatch_dist\"),\n", + " (\"Mean travel to patient\", \"t_to_patient\"),\n", + " (\"Mean scene time\", \"scene_time\"),\n", + " (\"Mean travel to hospital\", \"t_to_hospital\"),\n", + " (\"Mean return to base\", \"t_to_base\"),\n", + " (\"Mean total service\", \"total_service\")]:\n", + " c = immediate[col].mean() if len(immediate) else float(\"nan\")\n", + " f = delayed[col].mean() if len(delayed) else float(\"nan\")\n", + " \n", + " print(f\" {metric:<30} {c:>10.2f} {f:>10.2f}\")\n", + " print(\"─\" * 65)\n", + " print(f\" {'Ambulance':<14} {'Node':>10} {'Jobs':>8} {'Util':>10}\")\n", + " print(\"─\" * 65)\n", + " for amb in ambulances:\n", + " print(f\" Ambulance {amb.ambulance_id:<4} \"\n", + " f\"{amb.home_node:>5} {amb.total_jobs:>6} \"\n", + " f\"{amb.total_busy / run_length:>8.2%}\")\n", + " print(\"═\" * 65)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "19c9b8e8-7a2b-4d51-b070-9ba71463f44f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Simulation tracing set to: False\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
patient_idpatient_nodeambulance_idambulance_nodedispatch_distt_to_patientscene_timet_to_hospitalt_to_basetotal_serviceimmediate_dispatchwait
15215201042.82842711.31370817.3922338.0000008.00000044.705942True0.000000
1531493830.0000000.00000036.45605317.88854417.88854472.233141True0.000000
1541553624.00000016.0000000.48310717.88854417.88854452.260194True0.000000
1551601204.00000016.0000000.0670888.0000008.00000032.067088False13.713547
1561511310.0000000.00000055.8669578.0000008.00000071.866957True0.000000
\n", + "
" + ], + "text/plain": [ + " patient_id patient_node ambulance_id ambulance_node dispatch_dist \\\n", + "152 152 0 10 4 2.828427 \n", + "153 149 3 8 3 0.000000 \n", + "154 155 3 6 2 4.000000 \n", + "155 160 1 2 0 4.000000 \n", + "156 151 1 3 1 0.000000 \n", + "\n", + " t_to_patient scene_time t_to_hospital t_to_base total_service \\\n", + "152 11.313708 17.392233 8.000000 8.000000 44.705942 \n", + "153 0.000000 36.456053 17.888544 17.888544 72.233141 \n", + "154 16.000000 0.483107 17.888544 17.888544 52.260194 \n", + "155 16.000000 0.067088 8.000000 8.000000 32.067088 \n", + "156 0.000000 55.866957 8.000000 8.000000 71.866957 \n", + "\n", + " immediate_dispatch wait \n", + "152 True 0.000000 \n", + "153 True 0.000000 \n", + "154 True 0.000000 \n", + "155 False 13.713547 \n", + "156 True 0.000000 " + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "set_trace(False)\n", + "\n", + "# run model\n", + "ambulances, log = single_run(run_length=RUN_LENGTH, random_seed=42)\n", + "\n", + "# get all of the patient to ambulance assignment details\n", + "assignments = pd.DataFrame(log[\"assignments\"])\n", + "assignments.tail()" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "90dc1713-9d25-4e2c-adc8-f51a123c4534", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "═════════════════════════════════════════════════════════════════\n", + " Patients served : 157\n", + " Immediate Dispatch : 73 (46.5%)\n", + " Delayed Dispatch : 84 (53.5%)\n", + " Overall mean wait (min) : 13.34\n", + "─────────────────────────────────────────────────────────────────\n", + " Ambulance Dispatch\n", + " Metric Immediate Delayed\n", + "─────────────────────────────────────────────────────────────────\n", + " Mean wait (min) 0.00 24.93\n", + " Mean dispatch dist 1.35 2.79\n", + " Mean travel to patient 5.38 11.16\n", + " Mean scene time 20.32 19.71\n", + " Mean travel to hospital 11.12 11.30\n", + " Mean return to base 11.25 11.41\n", + " Mean total service 48.07 53.58\n", + "─────────────────────────────────────────────────────────────────\n", + " Ambulance Node Jobs Util\n", + "─────────────────────────────────────────────────────────────────\n", + " Ambulance 1 0 17 76.70%\n", + " Ambulance 2 0 18 85.23%\n", + " Ambulance 3 1 18 81.13%\n", + " Ambulance 4 1 14 82.45%\n", + " Ambulance 5 2 13 83.17%\n", + " Ambulance 6 2 14 78.95%\n", + " Ambulance 7 3 13 77.54%\n", + " Ambulance 8 3 13 76.57%\n", + " Ambulance 9 4 19 79.82%\n", + " Ambulance 10 4 18 79.39%\n", + "═════════════════════════════════════════════════════════════════\n" + ] + } + ], + "source": [ + "single_run_summary(assignments, ambulances, RUN_LENGTH)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "97b66363-bbda-4a91-b24a-4a10689bee5a", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/content/img/two_resource_types.drawio b/content/img/two_resource_types.drawio new file mode 100644 index 0000000..91f54a7 --- /dev/null +++ b/content/img/two_resource_types.drawio @@ -0,0 +1,188 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/content/img/two_resource_types.png b/content/img/two_resource_types.png new file mode 100644 index 0000000..a2b6514 Binary files /dev/null and b/content/img/two_resource_types.png differ diff --git a/content/sim_utility.py b/content/sim_utility.py new file mode 100644 index 0000000..906ddc5 --- /dev/null +++ b/content/sim_utility.py @@ -0,0 +1,87 @@ +import numpy as np +from typing import Optional + + +# Internal module state for tracing (default to True) +_trace_state = {"enabled": True} + +# ---------------------------------------------------------------------------- +# HELPER FUNCTIONS & DISTRIBUTIONS +# ---------------------------------------------------------------------------- + +def set_trace(state: bool): + ''' + Toggles the simulation trace output on or off. + + Usage: + ------ + Call this before running your simulation to enable/disable printing. + + >>> set_trace(True) + + Params: + ------- + state: bool + True to enable print output, False to disable. + ''' + _trace_state["enabled"] = state + print(f"Simulation tracing set to: {state}") + +def trace(msg): + ''' + Prints an event message if tracing is enabled. + + Params: + ------- + msg: str + string to print to screen. + ''' + if _trace_state["enabled"]: + print(msg) + +def spawn_seeds(n_streams: int, main_seed: Optional[int] = None) -> list[np.random.SeedSequence]: + """ + Taken from `sim-tools` + + Generate multiple statistically independent random seeds. + + This function creates a set of SeedSequence objects that are guaranteed + to produce independent streams of random numbers. This is crucial for + ensuring that multiple random number generators don't produce correlated + outputs, which could bias simulation results. + + Parameters + ---------- + n_streams : int + The number of independent seed sequences to generate. + Must be a positive integer. + + main_seed : Optional[int], default=None + Master seed that determines all generated sequences. + If None, a random entropy source is used, making results + non-reproducible across runs. Providing a value enables + reproducible sequences. + + Returns + ------- + List[np.random.SeedSequence] + A list of n_streams SeedSequence objects that can be used + to initialize random number generators with independent streams. + + Notes + ----- + This approach is preferred over manually creating seeds because + it uses NumPy's entropy pool management to guarantee statistical + independence between streams, avoiding subtle correlations that + might occur with manually chosen seeds. + + Examples + -------- + >>> seeds = spawn_seeds(3, main_seed=12345) + >>> rng1 = np.random.default_rng(seeds[0]) + >>> rng2 = np.random.default_rng(seeds[1]) + >>> rng3 = np.random.default_rng(seeds[2]) + """ + seed_sequence = np.random.SeedSequence(main_seed) + seeds = seed_sequence.spawn(n_streams) + return seeds