{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# T11 - Advanced features\n", "\n", "This tutorial covers advanced features of Covasim, including custom population options and changing the internal computational methods." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " \n", "Click [here](https://mybinder.org/v2/gh/institutefordiseasemodeling/covasim/HEAD?urlpath=lab%2Ftree%2Fdocs%2Ftutorials%2Ftut_advanced.ipynb) to open an interactive version of this notebook.\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining populations with SynthPops\n", "\n", "For complex populations, we suggest using [SynthPops](http://synthpops.org), a Python library designed specifically for this purpose. In contrast the population methods built-in to Covasim, SynthPops uses data to produce synthetic populations that are statistically indistinguishable from real ones. For a relatively complex example of how SynthPops was used to create a complex school network for the Seattle region, see [here](https://github.com/institutefordiseasemodeling/testing-the-waters/blob/main/covasim_schools/school_pop.py).\n", "\n", "## Defining contact layers\n", "\n", "As mentioned in Tutorial 4, contact layers are the graph connecting the people in the simulation. Each person is a node, and each contact is an edge. While enormous complexity can be used to define realistic contact networks, a reasonable approximation in many cases is random connectivity, often with some age assortativity. Here is an example for generating a new contact layer, nominally representing public transportation, and adding it to a simulation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import covasim as cv\n", "cv.options(jupyter=True, verbose=0)\n", "\n", "# Create the first sim\n", "orig_sim = cv.Sim(pop_type='hybrid', n_days=120, label='Default hybrid population')\n", "orig_sim.initialize() # Initialize the population\n", "\n", "# Create the second sim\n", "sim = orig_sim.copy()\n", "\n", "# Define the new layer, 'transport'\n", "n_people = len(sim.people)\n", "n_contacts_per_person = 0.5\n", "n_contacts = int(n_contacts_per_person*n_people)\n", "contacts_p1 = cv.choose(max_n=n_people, n=n_contacts)\n", "contacts_p2 = cv.choose(max_n=n_people, n=n_contacts)\n", "beta = np.ones(n_contacts)\n", "layer = cv.Layer(p1=contacts_p1, p2=contacts_p2, beta=beta) # Create the new layer\n", "\n", "# Add this layer in and re-initialize the sim\n", "sim.people.contacts.add_layer(transport=layer)\n", "sim.reset_layer_pars() # Automatically add layer 'transport' to the parameters using default values\n", "sim.initialize() # Reinitialize\n", "sim.label = f'Transport layer with {n_contacts_per_person} contacts/person'\n", "\n", "# Run and compare\n", "msim = cv.parallel(orig_sim, sim)\n", "msim.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining dynamic layers\n", "\n", "You can also define custom layers that update dynamically, e.g. based on a supplied number of contacts per day. To do this, create a new `Layer` class and define the `update()` method. For example:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import covasim as cv\n", "import numpy as np\n", "import sciris as sc\n", "\n", "class CustomLayer(cv.Layer):\n", " ''' Create a custom layer that updates daily based on supplied contacts '''\n", "\n", " def __init__(self, layer, contact_data):\n", " ''' Convert an existing layer to a custom layer and store contact data '''\n", " for k,v in layer.items():\n", " self[k] = v\n", " self.contact_data = contact_data\n", " return\n", "\n", " def update(self, people):\n", " ''' Update the contacts '''\n", " pop_size = len(people)\n", " n_new = self.contact_data[people.t] # Pull out today's contacts\n", " self['p1'] = np.array(cv.choose_r(max_n=pop_size, n=n_new), dtype=cv.default_int) # Choose with replacement\n", " self['p2'] = np.array(cv.choose_r(max_n=pop_size, n=n_new), dtype=cv.default_int) # Paired contact\n", " self['beta'] = np.ones(n_new, dtype=cv.default_float) # Per-contact transmission (just 1.0)\n", " return\n", "\n", "\n", "# Define some simple parameters\n", "pars = sc.objdict(\n", " pop_size = 1000,\n", " n_days = 90,\n", ")\n", "\n", "# Set up and run the simulation\n", "base_sim = cv.Sim(pars, label='Default simulation')\n", "sim = cv.Sim(pars, dynam_layer=True, label='Dynamic layers')\n", "sim.initialize()\n", "\n", "# Update to custom layer\n", "for key in sim.layer_keys():\n", " contact_data = np.random.randint(pars.pop_size*10, pars.pop_size*20, size=pars.n_days+1) # Generate a number of contacts for today\n", " sim.people.contacts[key] = CustomLayer(sim.people.contacts[key], contact_data=contact_data)\n", "\n", "# Run and plot\n", "msim = cv.parallel(base_sim, sim)\n", "msim.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining custom population properties\n", "\n", "Another useful feature is adding additional features to people, for use in subtargeting. For example, this example shows how to define a subpopulation with higher baseline mortality rates. This is a simple example illustrating how you would identify and target people based on whether or not the have a prime-number index, based on the protecting the elderly example from Tutorial 1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import sciris as sc\n", "import covasim as cv\n", "\n", "def protect_the_prime(sim):\n", " if sim.t == sim.day('2020-04-01'):\n", " are_prime = sim.people.prime\n", " sim.people.rel_sus[are_prime] = 0.0\n", "\n", "pars = dict(\n", " pop_type = 'hybrid',\n", " pop_infected = 100,\n", " n_days = 90,\n", " verbose = 0,\n", ")\n", "\n", "# Default simulation\n", "orig_sim = cv.Sim(pars, label='Default')\n", "\n", "# Create the simulation\n", "sim = cv.Sim(pars, label='Protect the prime', interventions=protect_the_prime)\n", "sim.initialize() # Initialize to create the people array\n", "sim.people.prime = np.array([sc.isprime(i) for i in range(len(sim.people))]) # Define whom to target\n", "\n", "# Run and plot\n", "msim = cv.parallel(orig_sim, sim)\n", "msim.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Changing Numba options\n", "\n", "Finally, this example shows how you can change the default Numba calculation options. It's not recommended – especially running with multithreading, which is faster but gives stochastically unreproducible results – but it's there if you want it." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import covasim as cv\n", "\n", "# Create a standard 32-bit simulation\n", "sim32 = cv.Sim(label='32-bit, single-threaded (default)', verbose='brief')\n", "sim32.run()\n", "\n", "# Use 64-bit instead of 32\n", "cv.options.set(precision=64)\n", "sim64 = cv.Sim(label='64-bit, single-threaded', verbose='brief')\n", "sim64.run()\n", "\n", "# Use parallel threading\n", "cv.options.set(numba_parallel=True)\n", "sim_par = cv.Sim(label='64-bit, multi-threaded', verbose='brief')\n", "sim_par.run()\n", "\n", "# Reset to defaults\n", "cv.options.set('defaults')\n", "sim32b = cv.Sim(label='32-bit, single-threaded (restored)', verbose='brief')\n", "sim32b.run()\n", "\n", "# Plot\n", "msim = cv.MultiSim([sim32, sim64, sim_par, sim32b])\n", "msim.plot()" ] } ], "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.9.10" }, "pycharm": { "stem_cell": { "cell_type": "raw", "metadata": { "collapsed": false }, "source": [] } }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }