SIRAsynchronousDynamicsRewireDegree.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from GraphWithAsynchronousDynamics import *\n",
    "import operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class SIRAsynchronousDynamicsRewireDegree(GraphWithAsynchronousDynamics):\n",
    "    '''A graph with a particular SIR dynamics. We use probabilities\n",
    "    to express infection and recovery per timestep, and run the system\n",
    "    using synchronous dynamics.'''\n",
    "    \n",
    "    # the possible dynamics states of a node for SIR dynamics\n",
    "    SUSCEPTIBLE = 'susceptible'\n",
    "    INFECTED = 'infected'\n",
    "    RECOVERED = 'recovered'\n",
    "    \n",
    "    # list of infected nodes, the sites of all the dynamics\n",
    "    _infected = []\n",
    "    _susceptible = []\n",
    "    _recovered = []\n",
    "    \n",
    "    DEGREE_DICT = dict()\n",
    "    \n",
    "    \n",
    "    def __init__( self, time_limit = 10000, p_infect = 0.0, p_recover = 1.0, p_infected = 0.0, p_rewire = 0.0, graph = None ):\n",
    "        '''Generate a graph with dynamics for the given parameters.\n",
    "        \n",
    "        pInfect: infection probability (defaults to 0.0)\n",
    "        pRecover: probability of recovery (defaults to 1.0)\n",
    "        pInfected: initial infection probability (defaults to 0.0)\n",
    "        g: the graph to copy from (optional)'''\n",
    "        states = {self.SUSCEPTIBLE,self.INFECTED,self.RECOVERED}\n",
    "        rates = dict()\n",
    "        rates['p_infect'] = p_infect\n",
    "        rates['p_recover'] = p_recover\n",
    "        rates['p_infected'] = p_infected\n",
    "        rates['p_rewire'] = p_rewire\n",
    "        GraphWithAsynchronousDynamics.__init__(self, time_limit = time_limit, graph = graph, states = states, rates = rates)\n",
    "        self._p_infect = p_infect\n",
    "        self._p_recover = p_recover\n",
    "        self._p_infected = p_infected\n",
    "        self._p_rewire = p_rewire\n",
    "        \n",
    "        for i in self.nodes():\n",
    "            self.DEGREE_DICT[i] = self.degree(i)\n",
    "            \n",
    "    def _before( self ):\n",
    "        '''Seed the network with infected nodes, and mark all edges\n",
    "        as unoccupied by the dynamics.'''\n",
    "        self._infected = []       # in case we re-run from a dirty intermediate state\n",
    "        for n in self.node.keys():\n",
    "            if numpy.random.random() <= self._p_infected:\n",
    "                self.node[n][self.DYNAMICAL_STATE] = self.INFECTED\n",
    "            else:\n",
    "                self.node[n][self.DYNAMICAL_STATE] = self.SUSCEPTIBLE\n",
    "        for (n, m, data) in self.edges_iter(data = True):\n",
    "            data[self.OCCUPIED] = False\n",
    "    \n",
    "    def _after(self):\n",
    "        pass\n",
    "    \n",
    "    def set_timestep_rate(self):\n",
    "        '''Calculate the minimum time-jump. Choose the maximum from the recovery rate\n",
    "        or the max node degree * the infection rate. Timestep will be 1 over this value'''\n",
    "        max_degree = max(self.degree().values())\n",
    "        return max(max_degree*self._p_infect, max_degree*self._p_rewire, self._p_recover)\n",
    "     \n",
    "    def at_equilibrium( self ):\n",
    "        '''SIR dynamics is at equilibrium if there are no more\n",
    "        infected nodes left in the network or if we've exceeded\n",
    "        the default simulation length.\n",
    "        \n",
    "        returns: True if the model has stopped'''\n",
    "        \n",
    "        '''if(t%1000 == 0.0):\n",
    "            print(\"Processing timestep \",t)'''\n",
    "        \n",
    "        if len(self.POPULATION[self.INFECTED]) < (0.8 * max(self._pop_dist[self.INFECTED].values())):\n",
    "            self.STATISTICS['peak_infection'] = max(self._pop_dist[self.INFECTED].iteritems(), key=operator.itemgetter(1))\n",
    "            return True\n",
    "        \n",
    "        \n",
    "        if self.CURRENT_TIMESTEP >= self._time_limit:\n",
    "            self.STATISTICS['peak_infection'] = max(self._pop_dist[self.INFECTED].iteritems(), key=operator.itemgetter(1))\n",
    "            return True\n",
    "        else:\n",
    "            return (len(self.POPULATION[self.INFECTED]) == 0)\n",
    "            \n",
    "    \n",
    "    def asyn_node_action(self, node = 0, dt = 0.0, r = 0.0):\n",
    "        \n",
    "        # If chosen node is infected\n",
    "        if(self.node[node][self.DYNAMICAL_STATE] == self.INFECTED):\n",
    "            # If random number lower than Number of nodes x rate recovery * timestep\n",
    "            if(r < self.order()*self._p_recover*dt):\n",
    "                # Node recovers\n",
    "                self.update_node(node,self.INFECTED,self.RECOVERED)\n",
    "                return True\n",
    "        # If chosen node is susceptible\n",
    "        elif(self.node[node][self.DYNAMICAL_STATE] == self.SUSCEPTIBLE):\n",
    "            neighbours = self.neighbors(node)\n",
    "            infected_neighbours = [ ne for ne in neighbours if self.node[ne][self.DYNAMICAL_STATE] == self.INFECTED ]\n",
    "            num_inf_neighb = len(infected_neighbours)\n",
    "            \n",
    "            infect_chance = self._p_infect*num_inf_neighb\n",
    "            rewire_chance = self._p_rewire*num_inf_neighb\n",
    "            if(r < self.order()*dt*infect_chance):\n",
    "                # Node gets infected\n",
    "                self.update_node(node,self.SUSCEPTIBLE,self.INFECTED)\n",
    "                return True\n",
    "            elif(r < self.order()*dt*infect_chance + self.order()*dt*rewire_chance):\n",
    "                # Node rewires\n",
    "                self.rewire(node, neighbours, infected_neighbours)\n",
    "        \n",
    "        return False\n",
    "    \n",
    "    def rewire(self, n = 0 , neighbours = [], infected_neighbours = []):\n",
    "        \n",
    "        # Remove a link to an infected neighbour\n",
    "        neighb_index = numpy.random.randint(len(infected_neighbours))\n",
    "        \n",
    "        # Record previous degree\n",
    "        previous_degree = self.degree(infected_neighbours[neighb_index])\n",
    "        \n",
    "        # Remove edge\n",
    "        self.remove_edges_from([(n, infected_neighbours[neighb_index])])\n",
    "        \n",
    "        # Update degree dictionary\n",
    "        self.DEGREE_DICT[n] -= 1\n",
    "        self.DEGREE_DICT[infected_neighbours[neighb_index]] -= 1\n",
    "        \n",
    "        # Add a link to a new node\n",
    "        \n",
    "        # Clone the list of susceptible and recovered nodes - these will be the candidates for the new link\n",
    "        candidates = list(self.POPULATION[self.SUSCEPTIBLE])\n",
    "        candidates.extend(list(self.POPULATION[self.RECOVERED]))                \n",
    "                        \n",
    "        # Remove the current node (don't want self loops)\n",
    "        if n in candidates:\n",
    "            candidates.remove(n)\n",
    "        \n",
    "        # Remove neighbours from candidates        \n",
    "        candidates = [c for c in candidates if c not in neighbours]\n",
    "       \n",
    "        # Create a dict to store how far away each candidate is in terms of degree from the disconnected node\n",
    "        pref_cand = dict()\n",
    "        for c in candidates:\n",
    "            pref_cand[c] = abs(self.DEGREE_DICT[c] - previous_degree)\n",
    "    \n",
    "        # Reduce the candidates down to only those who are nearest to the disconncected node in terms of degree\n",
    "        min_distance = min(pref_cand.values()) \n",
    "        candidates = [c for c in pref_cand if pref_cand[c] == min_distance]\n",
    "    \n",
    "        # Pick a candidate\n",
    "        if(len(candidates) >= 1):\n",
    "            new_neighb_index = numpy.random.randint(len(candidates))\n",
    "            self.add_edge(n, candidates[new_neighb_index])\n",
    "            self.DEGREE_DICT[n] += 1\n",
    "            self.DEGREE_DICT[candidates[new_neighb_index]] += 1\n",
    "        \n",
    "        # Network topology has changed so alter the timestep rate\n",
    "        max_trans_rate = self.set_timestep_rate()\n",
    "        \n",
    "        # Timestep = 1 over number of nodes by maximum transition rate\n",
    "        self.DT = 1.0/(self.order()*max_trans_rate);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}