<jupyter_start><jupyter_text>CX 4230, Spring 2016 [12]: Markov chain analysis of the SIR-CA infection modelThe goal of this notebook is to analyze theoretically the SIR-CA model of [Lab 11](http://nbviewer.jupyter.org/github/rvuduc/cx4230-ipynbs/blob/master/11--infection.ipynb). To complete it, you will need the SIR-CA machinery from Lab 9; we've packaged that up in a module that one of the code cells below will download. Recall our motivation for a different kind of analysis, alluded to in the [slides for Lab 11](https://t-square.gatech.edu/access/content/group/gtc-59b8-dc03-5a67-a5f4-88b8e4d5b69a/cx4230-sp16--11-infection-and-diffusion.pdf): there is a combinatorial explosion in the number of possible states as the recovery time or problem size increase; if you just need some very high-level summary statistics, these simulations might be overkill.  SetupRun the following code cells to get started.<jupyter_code># Our usual multidimensional array tools
import numpy as np
import scipy as sp
import scipy.sparse

# Some handy combinatorial functions, e.g.,
# for generating permutations and combinations.
# See: https://docs.python.org/3/library/itertools.html
import itertools
 # Core plotting support
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>We packaged up the SIR-CA implementation from Lab 11 into a module. The following code cell will download and import some of the functions and constants of this module.<jupyter_code>import os
os.makedirs ('cx4230', exist_ok=True)

# http://stackoverflow.com/questions/7243750/download-file-from-web-in-python-3
import urllib.request
MOD_URL = 'https://github.com/rvuduc/cx4230-ipynbs/raw/master/cx4230/sirca.py'
MOD_LOCAL = 'cx4230/sirca.py'
urllib.request.urlretrieve (MOD_URL, MOD_LOCAL)

%reload_ext autoreload
%autoreload 2

from cx4230 import sirca
from cx4230.sirca import \
    create_new_grid, set_recovery_time, infected, exposed, susceptible, \
    show_peeps, sim, isim<jupyter_output><empty_output><jupyter_text>Run this cell to verify the module is working.<jupyter_code>from ipywidgets import interact

interact (isim
          , m=(1, 10, 1)
          , n=(1, 10, 1)
          , max_steps=(0, 100, 1)
          , tau=(0.0, 1.0, 0.1)
          , dt=(0, 10, 1)
          , seed=(0, 100, 1)
         )<jupyter_output><empty_output><jupyter_text>Target system to analyze: a 1-D, 3-cell gridPer the lecture notes, let's analyze the convergence of a 1-D system consisting of just 3 cells, where the conditional probability of infection for any exposed cell is $\tau=0.2$ and the recovery time is $k=2$ days.<jupyter_code>M, N = 1, 3 # M x N grid
TAU = 0.2
K = 2

set_recovery_time (K)<jupyter_output><empty_output><jupyter_text>Defining the state spaceTo carry out the Markov chain analysis, recall that you need to do the following.1. Define the state space, $S$.2. Construct a probability transition matrix, $P \equiv \left(p_{ij}\right)$ where $p_{ij} \equiv \mathrm{Pr}[s_j\left|s_i\right.]$ is the conditional probability of moving to state $s_j$ starting from state $s_i$ for every pair of states $s_i, s_j \in S$. For the 1-D, 3-cell system, the state space $S$ consists of all possible grid states, $s \equiv (g_1, g_2, g_3)$ where $g_i \in \left\{ \mathtt{S}, \mathtt{I}_1, \ldots, \mathtt{I}_k, \mathtt{R} \right\}$.<jupyter_code>from cx4230.sirca import EMPTY, SUSCEPTIBLE, INFECTED, RECOVERED
CELL_VALUES = range (SUSCEPTIBLE, RECOVERED+1)
print ("Possible cell values:", list (CELL_VALUES))<jupyter_output><empty_output><jupyter_text>**Exercise.** Generate a Python [set()](https://docs.python.org/3/library/stdtypes.htmlset) named `S` that contains all possible grid-state triples.> Hint: Take a look at [`itertools.product()`](https://docs.python.org/3/library/itertools.htmlitertools.product).<jupyter_code># @YOUSE: Generate S
S = set (...)

# Test:
print (S)<jupyter_output><empty_output><jupyter_text>To construct the probability transition matrix, $P$, we will eventually need to map individual states to integer indices. The following code builds two lookup tables for converting grid states to integers and vice-versa.<jupyter_code>IND_TO_S = dict (enumerate (S))
S_TO_IND = {val: key for (key, val) in IND_TO_S.items ()}

print ("==> Index to state table:\n", IND_TO_S, "\n")
print ("==> State to index table:\n", S_TO_IND)<jupyter_output><empty_output><jupyter_text>The states are tuples whereas the grid in our simulation is a 2-D array (with an artificial empty boundary). Here are some handy functions to convert between a grid and a state.<jupyter_code>def convert_state_to_grid (s, m, n):
    """
    Converts a state, given as a tuple `s`, into an equivalent
    `m` x `n` grid.
    """
    G = create_new_grid (m, n)
    G_interior = G[1:-1, 1:-1]
    G_interior[:, :] = np.array (list (s)).reshape ((m, n))
    return G

def convert_grid_to_state (G):
    """
    Converts a an `m` x `n` grid `G` (+ its empty boundary)
    into a `m*n`-tuple.
    """
    G_interior = G[1:-1, 1:-1]
    return tuple (G_interior.reshape (np.prod (G_interior.shape)))
# Test
G0 = create_new_grid (M, N)
G0[int ((M+2)/2), int ((N+2)/2)] = INFECTED

print ("==> Original grid:")
print (G0, "\n")

print ("==> Equivalent state (id):")
s0 = convert_grid_to_state (G0)
i0 = S_TO_IND[s0]
print ("State", s0, "--> index", i0, "\n")

print ("==> Convert back to grid:")
s0_prime = IND_TO_S[i0]
G0_prime = convert_state_to_grid (s0_prime, M, N)
print (G0_prime, "\n")

assert (G0 == G0_prime).all ()
print ("==> Passed!")<jupyter_output><empty_output><jupyter_text>ReachabilityThe last thing you need to compute $p_{ij}$ is a way to enumerate all _reachable_ $j$ values, given $i$. The following code cells build a function to determine reachability.First, here is an auxiliary function to get a list of tuples that are the locations of any non-zero value of a given grid.<jupyter_code>def find_nonzero_locs (G):
    """
    Given a grid G, returns a list of all locations (i, j)
    where G[i, j] is not equal to zero.
    """
    (I, J) = np.nonzero (G)
    return list (zip (I, J))

# Test
print (find_nonzero_locs (exposed (G0)))<jupyter_output><empty_output><jupyter_text>**Exercise.** Let $\tau$ be the conditional probability that an exposed individual becomes infected. Further suppose the grid contains $n_e$ exposed individuals.1. How many different ways are there for $n_s$ individuals to become infected? (Assume $0 \leq n_s \leq n_e$.)2. What is the probability that $n_s$ of the $n_e$ individuals will become infected? > _Answer_:> > 1. ??> 2. ??<jupyter_code>def reachable (G, tau):
    # Find all infected cells
    I = infected (G)
    if np.sum (I) == 0: # No infected cells
        yield (1.0, G)
    else: # >= 1 infected cell
        # Find number and locations of all exposed cells
        locs_E = find_nonzero_locs (exposed (G) * susceptible (G))
        n_E = len (locs_E)

        # Try all possible spreading combinations
        for k in range (0, n_E+1):
            # Probability that a particular set of k exposed people become infected:
            prob_k = (tau**k) * ((1.0 - tau)**(n_E - k))

            # Enumerate all reachable infection states
            for spread_locs_k in itertools.combinations (locs_E, k):
                G_next = np.copy (G) + I
                for (i, j) in spread_locs_k:
                    G_next[i, j] = INFECTED
                yield (prob_k, G_next)
n_states = len (S)
nz_i = []
nz_j = []
nz_val = []
for s_i in S:
    i = S_TO_IND[s_i]
    G_i = convert_state_to_grid (s_i, M, N)
    
    for (p_ij, G_j) in reachable (G_i, TAU):
        s_j = convert_grid_to_state (G_j)
        j = S_TO_IND[s_j]
        nz_i.append (i)
        nz_j.append (j)
        nz_val.append (p_ij)
            
P = sp.sparse.coo_matrix ((nz_val, (nz_i, nz_j)), (n_states, n_states))
print ("no. of states ==", n_states)
print ("tau ==", TAU, "; P = (p_{ij}) =\n", P)
print ("nnz(P) ==", P.nnz)
print ("row sums ==", P.sum (axis=1))
# Index case:
x0 = np.zeros ((n_states, 1))
x0[S_TO_IND[(0, 1, 0)]] = 1.0

# Compute the steady-state distribution
x = x0
for t in range (5):
    x = P.transpose ().dot (x)

# Summarize the "interesting" states
x_sorted = sorted (list (enumerate (x)), key=lambda t: t[1], reverse=True)
x_interesting = list (filter (lambda t: t[1] > 0.0, x_sorted))
print ("%d state(s) have a non-zero steady-state probability:" % len (x_interesting))
for (i, x_i) in x_interesting:
    print ("%d:" % i, "Pr[%s] ==" % str(IND_TO_S[i]), x_i[0])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.126153486742203, 'fraction_numerical': 0.022894521668029435, 'mean_word_length': 3.8564946114577423, 'pattern_counts': {'":': 0, '<': 41, '<?xml version=': 0, '>': 55, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '209603', 'license': 'MIT', 'n_tokens_mistral': 3155, 'n_tokens_neox': 2988, 'n_words': 989, 'path': 'Lab/12--infection-markov.ipynb', 'repo_name': 'dem123456789/Computer-Simulation', 'size': '8561'}