{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## CSCS530 Winter 2015\n", "#### Complex Systems 530 - Computer Modeling of Complex Systems (Winter 2015)\n", "\n", " * Course ID: CMPLXSYS 530\n", " * Course Title: Computer Modeling of Complex Systems\n", " * Term: Winter 2015\n", " * Schedule: Wednesdays and Friday, 1:00-2:30PM ET\n", " * Location: 120 West Hall (http://www.lsa.umich.edu/cscs/research/computerlab)\n", " * Teachers: [Mike Bommarito](https://www.linkedin.com/in/bommarito) and [Sarah Cherng](https://www.linkedin.com/pub/sarah-cherng/35/1b7/316)\n", "\n", "#### [View this repository on NBViewer](http://nbviewer.ipython.org/github/mjbommar/cscs-530-w2015/tree/master/)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## HIV Model Outline\n", "\n", "### Goal\n", " We will examine policy interventation strategies to address and control the spread of HIV in a sexual network.\n", "\n", "### Justification\n", " Given the heterogeneity of behavior and network effects of sexually-transmitted diseases, this problem requires ABM and network methodology. Differential equation compartment models will not allow for incomplete or scale-free graph structures or time-varying agent behavior, and so cannot be used. \n", " \n", "### Processes and Mechanisms of Interest\n", "\n", " We are primarily interested in how different ranges of the condom subsidy and condom budget affect the proportion of people that are infected with HIV at each time step. Additionally, we are interested in the rate at which HIV infection could increase given increased condom costs. \n", " \n", "### Outline\n", "\n", " * Implement condom subsidy mechanics\n", " * Implement gossip network to share status\n", " * Vary condom subsidy and gossip network structure\n", " * Observe how number of steps to 50% infected\n", " * Observe the mean and standard deviation of sexual partner degree distribution\n", " \n", " To formally describe our model, let's break it down into pieces:\n", "\n", "### I. Space\n", " In this model, our space will be a two-dimensional (2D) square grid. Each grid cell will contain zero or one people. Edges of the grid will wrap around.\n", "\n", "### II. Actors\n", "\n", "#### A. People\n", " In this case, people are the carriers of HIV. We are modeling them as simply as possible, using only the following properties:\n", " \n", " ** Properties **\n", " * _is_infected_: is the person infected with HIV?\n", " * _condom_budget_: what is the person's budget for protection?\n", " * _prob_hookup_: what is the probability that a person will want to hookup given contact?\n", " \n", "For their step function, agents will perform the following:\n", " \n", " * take a unit-distance step in a random direction\n", " * evaluate neighboring cells for any potential people to hookup with\n", " * for each pair, check against _prob_hookup_; if both people sample True, then hookup\n", " * for a hookup, check against _max_budget_ and the global condom cost to see if either partner will purchase\n", "\n", "#### B. Institution\n", "\n", " A \"public health\" institution will manage the level of subsidy for condoms. For now, the subsidy level will be set at $t_0$ and constant over model run.\n", "\n", "\n", "### III. Initial Conditions\n", "\n", "#### A. People\n", " \n", " * People will be randomly distributed throughout the grid by sampling from a uniform discrete distribution with replacement. In the event that an agent has already been placed at the sampled location, we continue drawing a new nandom position until it is unoccupied.\n", " * People will have their ``prob_hookup`` randomly initialized to a value from a uniform continuous distribution.\n", " * People will have their ``condom_budget`` randomly initialized to a value from a uniform continuous distribution.\n", " * A single person will be selected at random at $t_0$ to have HIV.\n", "\n", "#### B. Insitution\n", " \n", " The public health institution will set a level of subsidy for condoms using a random uniform continuous distribution.\n", " \n", "### IV. Model Parameters\n", "\n", " Based on the description above, we need the following model parameters:\n", " \n", " * ``grid_size``: size of the two-dimensional square grid, i.e., dimension length\n", " * ``num_people``: number of persons to create; must be less than ${grid\\_size}^2$\n", " * ``min_subsidy``, ``max_subsidy``: the lower and upper bounds for the institution's subsidy\n", " * ``min_condom_budget``, ``max_condom_budget``: the lower and upper bounds for initial conditions of people's ``condom_budget``\n", " * ``condom_cost``: cost of a condom; fixed to $1.0$\n", " * ``min_prob_hookup``, ``max_prob_hookup``: the lower and upper bounds for initial conditions of people's ``prob_hookup``\n", " * ``prob_transmit``, ``prob_transmit_condom``: the probability of transmitting the disease without and with a condom" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "# Standard imports\n", "import copy\n", "import itertools\n", "\n", "# Scientific computing imports\n", "import numpy\n", "import matplotlib.pyplot as plt\n", "import networkx\n", "import pandas\n", "import seaborn; seaborn.set()\n", "\n", "# Import widget methods\n", "from IPython.html.widgets import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Person Class\n", "\n", " Below, we will define our person class. This can be broken up as follows:\n", " \n", " * __constructor__: class constructor, which \"initializes\" or \"creates\" the person when we call ``Person()``. This is in the ``__init__`` method.\n", " * ``decide_condom``: decide if the person will purchase a condom, i.e., by checking that $condom\\_budget >= condom\\_cost - subsidy$\n", " * ``decide_hookup``: decide if the person will hookup, i.e., by sampling with $p=prob\\_hookup$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "class Person(object):\n", " \"\"\"\n", " Person class, which encapsulates the entire behavior of a person.\n", " \"\"\"\n", " \n", " def __init__():\n", " \"\"\"\n", " Constructor for Person class. By default,\n", " * not infected\n", " * will always buy condoms\n", " * will hookup 50% of the time\n", " \n", " Note that we must \"link\" the Person to their \"parent\" Model object.\n", " \"\"\"\n", " # Set model link and ID\n", " self.model = model\n", " self.person_id = person_id\n", " \n", " # Set Person parameters.\n", " self.is_infected = is_infected\n", " self.condom_budget = condom_budget\n", " self.prob_hookup = prob_hookup\n", " \n", " def decide_condom(self):\n", " \"\"\"\n", " Decide if we will use a condom.\n", " \"\"\"\n", " if self.condom_budget >= (self.model.condom_cost - self.model.condom_subsidy):\n", " return True\n", " else:\n", " return False\n", " \n", " def decide_hookup(self):\n", " \"\"\"\n", " Decide if we want to hookup with a potential partner.\n", " \"\"\"\n", " if numpy.random.random() <= self.prob_hookup:\n", " return True\n", " else:\n", " return False\n", " \n", " def get_position(self):\n", " \"\"\"\n", " Return position, calling through model.\n", " \"\"\"\n", " return self.model.get_person_position(self.person_id)\n", " \n", " def get_neighbors(self):\n", " \"\"\"\n", " Return neighbors, calling through model.\n", " \"\"\"\n", " return self.model.get_person_neighbors(self.person_id)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Model Class\n", "\n", " Below, we will define our model class. This can be broken up as follows:\n", " \n", " * __constructor__: class constructor, which \"initializes\" or \"creates\" the model when we call ``Model()``. This is in the ``__init__`` method.\n", " * __``setup_space``__: method to create our \"space\"\n", " * __``setup_people``__: method to create our \"people\"\n", " * __``setup_institution``__: method to create our \"institution\"\n", " * __``get_neighbors``__: method to get neighboring agents based on position\n", " * __``get_person_neighbors``__: method to get neighboring agents based on agent ID\n", " * __``get_person_position``__: method to get position based on agent ID\n", " * __``move_person``__: method to move an agent to a new position\n", " * __``step_move``__: method to step through agent moves\n", " * __``step_interact``__: method to step through agent interaction\n", " * __``step``__: main step method to control each time step simulation" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "class Model(object):\n", " \"\"\"\n", " Model class, which encapsulates the entire behavior of a single \"run\" in our HIV ABM.\n", " \"\"\"\n", " \n", " def __init__( ):\n", " \"\"\"\n", " Class constructor.\n", " \"\"\"\n", " # Set our model parameters; this is long but simple!\n", " self.grid_size = grid_size\n", " self.num_people = num_people\n", " self.min_subsidy = min_subsidy\n", " self.max_subsidy = max_subsidy\n", " self.min_condom_budget = min_condom_budget\n", " self.max_condom_budget = max_condom_budget\n", " self.condom_cost = condom_cost\n", " self.min_prob_hookup = min_prob_hookup\n", " self.max_prob_hookup = max_prob_hookup\n", " self.prob_transmit = prob_transmit\n", " self.prob_transmit_condom = prob_transmit_condom\n", " \n", " # Set our state variables\n", " self.t = 0\n", " self.space = numpy.array((0,0))\n", " self.condom_subsidy = 0.0\n", " self.people = []\n", " self.num_interactions = 0\n", " self.num_interactions_condoms = 0\n", " self.num_infected = 0\n", " \n", " # Setup our history variables.\n", " self.history_space = []\n", " self.history_space_infected = []\n", " self.history_interactions = []\n", " self.history_num_infected = []\n", " self.history_num_interactions = []\n", " self.history_num_interactions_condoms = []\n", " \n", " # Call our setup methods to initialize space, people, and institution.\n", " self.setup_space()\n", " self.setup_people()\n", " self.setup_institution()\n", " \n", " def setup_space(self):\n", " \"\"\"\n", " Method to setup our space.\n", " \"\"\"\n", " # Initialize a space with a NaN's\n", " self.space = numpy.full((self.grid_size, self.grid_size), numpy.nan)\n", " \n", " def setup_people(self):\n", " \"\"\"\n", " Method to setup our people.\n", " \"\"\"\n", " \n", " # First, begin by creating all agents without placing them.\n", " for i in xrange(self.num_people):\n", " self.people.append(Person(model=self,\n", " person_id=i,\n", " is_infected=False,\n", " condom_budget=numpy.random.uniform(self.min_condom_budget, self.max_condom_budget),\n", " prob_hookup=numpy.random.uniform(self.min_prob_hookup, self.max_prob_hookup)))\n", " \n", " # Second, once created, place them into the space.\n", " for person in self.people:\n", " # Loop until unique\n", " is_occupied = True\n", " while is_occupied:\n", " # Sample location\n", " random_x = numpy.random.randint(0, self.grid_size)\n", " random_y = numpy.random.randint(0, self.grid_size)\n", " \n", " # Check if unique\n", " if numpy.isnan(self.space[random_x, random_y]):\n", " is_occupied = False\n", " else:\n", " is_occupied = True\n", " \n", " # Now place the person there by setting their ID.\n", " self.space[random_x, random_y] = person.person_id\n", " \n", " # Third, pick one person to be infected initially.\n", " random_infected = numpy.random.choice(range(self.num_people))\n", " self.people[random_infected].is_infected = True\n", " self.num_infected += 1\n", "\n", " def setup_institution(self):\n", " \"\"\"\n", " Method to setup our space.\n", " \"\"\"\n", " # Randomly sample a subsidy level\n", " self.condom_subsidy = numpy.random.uniform(self.min_subsidy, self.max_subsidy)\n", " \n", " def get_neighborhood(self, x, y, distance=1):\n", " \"\"\"\n", " Get a Moore neighborhood of distance from (x, y).\n", " \"\"\"\n", " neighbor_pos = [ ( x % self.grid_size, y % self.grid_size)\n", " for x, y in itertools.product(xrange(x-distance, x+distance+1),\n", " xrange(y-distance, y+distance+1))]\n", " return neighbor_pos\n", " \n", " def get_neighbors(self, x, y, distance=1):\n", " \"\"\"\n", " Get any neighboring persons within distance from (x, y).\n", " \"\"\"\n", " neighbor_pos = self.get_neighborhood(x, y, distance)\n", " neighbor_list = []\n", " for pos in neighbor_pos:\n", " # Skip identity\n", " if pos[0] == x and pos[1] == y:\n", " continue\n", " \n", " # Check if empty\n", " if not numpy.isnan(self.space[pos[0], pos[1]]):\n", " neighbor_list.append(int(self.space[pos[0], pos[1]]))\n", " \n", " return neighbor_list\n", " \n", " def get_person_position(self, person_id):\n", " \"\"\"\n", " Get the position of a person based on their ID.\n", " \"\"\"\n", " # Find the value that matches our ID in self.space, then reshape to a 2-element list.\n", " return numpy.reshape(numpy.where(m.space == person_id), (1, 2))[0].tolist()\n", "\n", " def get_person_neighbors(self, person_id, distance=1):\n", " \"\"\"\n", " Get the position of a person based on their ID.\n", " \"\"\"\n", " # Find the value that matches our ID in self.space, then reshape to a 2-element list.\n", " x, y = self.get_person_position(person_id)\n", " return self.get_neighbors(x, y, distance) \n", " \n", " def move_person(self, person_id, x, y):\n", " \"\"\"\n", " Move a person to a new (x, y) location.\n", " \"\"\"\n", " \n", " # Get original\n", " original_position = self.get_person_position(person_id)\n", " \n", " # Check target location\n", " if not numpy.isnan(self.space[x, y]):\n", " raise ValueError(\"Unable to move person {0} to ({1}, {2}) since occupied.\".format(person_id, x, y))\n", " \n", " # Otherwise, move by emptying and setting.\n", " self.space[original_position[0], original_position[1]] = numpy.nan\n", " self.space[x, y] = person_id\n", " \n", " def step_move(self):\n", " \"\"\"\n", " Model step move function, which handles moving agents randomly around.\n", " \"\"\"\n", " \n", " def step_interact(self):\n", " \"\"\"\n", " \"Interact\" the agents by seeing if they will hookup and spread.\n", " \"\"\"\n", " \n", " # Get a random order for the agents.\n", " random_order = range(self.num_people)\n", " numpy.random.shuffle(random_order)\n", " \n", " # Track which pairs we've tested. Don't want to \"interact\" them twice w/in one step.\n", " seen_pairs = []\n", " \n", " # Iterate in random order.\n", " for i in random_order:\n", " # Get neighbors\n", " neighbors = self.get_person_neighbors(i)\n", " \n", " # Iterate over neighbors\n", " for neighbor in neighbors:\n", " # Check if we've already seen.\n", " a = min(i, neighbor)\n", " b = max(i, neighbor)\n", " if (a, b) not in seen_pairs:\n", " seen_pairs.append((a, b))\n", " else:\n", " continue\n", " \n", " # Check if hookup if not seen.\n", " hookup_a = self.people[a].decide_hookup()\n", " hookup_b = self.people[b].decide_hookup()\n", " if hookup_a and hookup_b:\n", " # Hookup going to happen. \n", " self.num_interactions += 1\n", " \n", " # Check now for condoms and use resulting rate.\n", " if self.people[a].decide_condom() or self.people[b].decide_condom():\n", " # Using a condom.\n", " self.num_interactions_condoms += 1\n", " use_condom = True\n", " \n", " if self.people[a].is_infected or self.people[b].is_infected:\n", " is_transmission = numpy.random.random() <= self.prob_transmit_condom\n", " else:\n", " is_transmission = False\n", " else:\n", " # Not using a condom.\n", " use_condom = False\n", " if self.people[a].is_infected or self.people[b].is_infected:\n", " is_transmission = numpy.random.random() <= self.prob_transmit\n", " else:\n", " is_transmission = False\n", " \n", " # Now infect.\n", " self.history_interactions.append((self.t, a, b, use_condom, is_transmission))\n", " if is_transmission:\n", " self.people[a].is_infected = True\n", " self.people[b].is_infected = True\n", " \n", " def get_num_infected(self):\n", " \"\"\"\n", " Get the number of infected persons.\n", " \"\"\"\n", " # Count\n", " infected = 0\n", " for person in self.people:\n", " if person.is_infected:\n", " infected += 1\n", " \n", " return infected\n", " \n", " def step(self):\n", " \"\"\"\n", " Model step function.\n", " \"\"\"\n", " \n", " # \"Interact\" agents.\n", " self.step_interact()\n", " \n", " # Move agents\n", " self.step_move()\n", " \n", " # Increment steps and track history.\n", " self.t += 1\n", " self.history_space.append(copy.deepcopy(self.space))\n", " self.history_space_infected.append(self.get_space_infected())\n", " self.num_infected = self.get_num_infected()\n", " self.history_num_infected.append(self.num_infected)\n", " self.history_num_interactions.append(self.num_interactions)\n", " self.history_num_interactions_condoms.append(self.num_interactions_condoms)\n", "\n", " def get_space_infected(self, t=None):\n", " \"\"\"\n", " Return a projection of the space that shows which cells have an infected person.\n", " \"\"\"\n", " if t == None:\n", " # Initialize empty\n", " infected_space = numpy.zeros_like(self.space)\n", " \n", " # Iterate over persons and set.\n", " for p in self.people:\n", " x, y = self.get_person_position(p.person_id)\n", " if p.is_infected:\n", " infected_space[x, y] = +1\n", " else:\n", " infected_space[x, y] = -1\n", " \n", " # Return\n", " return infected_space\n", " else:\n", " # Return historical step\n", " return self.history_space_infected[t]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Results \n", "\n", "\n", "### Overview \n", " We hope to present plots that demonstrate how sweeping over percentage changes in subsidy rates at: [0.8,0.9,1.0,1.1,1.1,1.2,1.3] will result in changes to HIV infection rates. These plots will include: 1) HIV Infection Rate by Condom Subsidy Rate at time step 100, 2) HIV Infection Rate over time by Condom Subsidy rate for 3 selected case\n", "\n", "### Hypothetical Results\n", " While adjusting condom subsidy up, we expect that condom use will increase and subsequently, the HIV infection rates will decrease over time. We believe that there may be a limit for which HIV infection rates will decrease to, and then remain at over time. " ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.6" } }, "nbformat": 4, "nbformat_minor": 0 }