#!/usr/bin/env python # coding: utf-8 # # Table of Contents #
# After comparing Keeling's Matlab code and output, and using multiple print statements throughout my code, I am still not sure why the size of the nearest infected farm is not significant. I found a couple of bugs, but they did not change the results. I also remembered that the reason I had previously decided not to use Keeling's kernel and use a Cauchy kernel was because Keeling's kernel makes the disease progress very slowly. So I have used the Cauchy kernel in the code below. The changes that I have made result in output that is consistent with Keeling's code, however I have noticed that in this new version the number of new infected is dramatically lower than it used to be.
# In[37]:
import numpy as np
import pandas as pd
from numpy.random import randint
# In[33]:
region_size = 20
N = 75
np.random.seed(53080)
x = region_size*np.random.rand(N)
y = region_size*np.random.rand(N)
# In[12]:
cows = np.concatenate([randint(25, 51, size=15),
randint(51, 76, size=30),
randint(76, 95, size=20),
randint(95, 250, size=10)])
sheep = np.concatenate([randint(25, 51, size=15),
randint(51, 76, size=30),
randint(76, 95, size=20),
randint(95, 250, size=10)])
# In[13]:
def which_grid(x, y, x_range, y_range, x_num, y_num):
'''
Calculates which grid square a particular location is in (turn a 2-d coordinate into a scalar)
Essentially: floor(Unif[0,1)griddim)griddim+floor(Unif[0,1)griddim)+1
Returns a number from 1 to griddim^2
'''
return(np.floor(x*(x_num/x_range))*y_num + np.floor(y*(y_num/y_range)) + 1)
# In[18]:
dist_squared = 3
dist_squared = np.asarray(dist_squared)
# In[8]:
def kernel(dist_squared):
dist_squared = np.asarray(dist_squared)
is_scalar = not dist_squared.ndim
# NOT CLEAR WHAT YOU ARE DOING BELOW
dist_squared.shape = (1,)*(1-dist_squared.ndim) + dist_squared.shape
K = 1 / (np.pi * (1 + dist_squared**2))
K[(dist_squared < 0.0138)] = 0.3093
K[(dist_squared > 60*60)] = 0
if is_scalar:
return K
return K[0]
# In[9]:
def iterate(status, x, y, suscept, transmiss, grid, first_in_grid, last_in_grid, num, max_rate):
event = 0*status
infected = np.where(status>5)[0]
number_infected = infected.size # Note reported farms still infectious
infected_grids = grid[infected]-1
for i in range(number_infected):
infected_i = infected[i]
#transmissibility of infected farm to all other grid squares
trans = -Transmiss[infected_i] * num
#max number of animals to be infected in infected grid square
max_r = max_rate[infected_grids[i], :]
# Elementwise multiplication
# max number of animals to be infected in each grid square based on infected grid square
rate = trans * max_r
#Max probability that infected farm infected noninfected farm
max_prob = 1 - np.exp(rate)
u = np.random.rand(len(max_prob))
#these grid squares need further consideration
m = np.where((max_prob - u)>0)[0]
for n in range(len(m)):
s = 1
M = m[n]
#Max probability that infected farm infects noninfected farms under consideration
prob_infect = 1 - np.exp(-transmiss[infected_i] * max_rate[infected_grids[i], M])
if (prob_infect == 1):
# Calculate the infection probability for each farm in the susceptible grid
leng = last_in_grid[M] - first_in_grid[M] + 1
R = np.random.rand(leng)
for j in range(leng):
ind1 = first_in_grid[M]+j-1
Q = 1 - np.exp(-transmiss[infected_i]
*suscept[ind1]
*kernel((x[infected_i]-x[ind1])**2 + (y[infected_i]-y[ind1])**2))
if ((R[j] < Q) & (status[ind1] == 0)):
event[ind1] = 1
else:
R = np.random.rand(num[M])
# Loop through all susceptible farms in the grids where an infection event occurred.
for j in range(num[M]):
P = 1 - s*(1 - prob_infect)**(num[M] - j)
if (R[j] < (prob_infect / P)):
s = 0
ind1 = first_in_grid[M]+j-1
# HAVE YOU TRIED ANALYZING JUST THIS FUNCTION FOR SENSITIVITY TO
# FARM SIZE?
Q = 1 - np.exp(-transmiss[infected_i]
*suscept[ind1]
*kernel((x[infected_i]-x[ind1])**2 + (y[infected_i]-y[ind1])**2))
if ((R[j] < Q/P) & (status[ind1] == 0)):
event[ind1] = 1
# Evolve the infection process of those farms which have been exposed and already infectious ones.
status[status > 0] += 1
status += event
return {'Status':status, 'NI':number_infected}
# THE FUNCTION BELOW IS CUMBERSOME. BREAK IT DOWN INTO FUNCTIONS, AND BE SURE TO WRITE TESTS TO ENSURE THAT EACH PIECE IS OPERATING AS EXPECTED
# In[10]:
def Outbreaks(Size,N,Y0,farms,end,end2,x,y,Cows,Sheep,Maxtime=1000):
#This is an attempt of converting the Matlab Program 7.6 Code into Python
Status = np.array([0]*N) #Initial Status of each farm
init_ind = np.random.randint(0,N)
for i in range(Y0):
Status[init_ind] = 6 #one farm is initially infected
#Cows are 10.5 times more susceptible to disease than sheep
Suscept = Sheep + 10.5*Cows
Transmiss = 5.1e-7*Sheep + 7.7e-7*Cows
#Set up the grid
grid = WhichGrid(x,y,Size,Size,10.0,10.0)
tmp = sorted(grid) #Sort grid values
# DO WE REALLY NEED ALL OF THE BELOW? IF SO, BETTER TO INITIALIZE A PANDAS
# DATA FRAME THAN ALL THESE LISTS
#i = np.argsort(grid) #get indexed values after sort
i = [i[0] for i in sorted(enumerate(grid), key=lambda x:x[1])]
x = x[i]
y = y[i]
Status = Status[i]
grid = grid[i]
Transmiss = Transmiss[i]
Suscept = Suscept[i]
Cows = Cows[i]
Sheep = Sheep[i]
Xgrid = []
Ygrid = []
Num = []
first_in_grid = []
last_in_grid = []
Max_Sus_grid = []
index_inf = np.where(Status==6)[0].astype(int)
m2 = np.array(np.where(grid==1))
for i in range(1,int(max(grid))+1):
#turn the grid square number into an x-coordinate and y-coordinate (should not exceed XNum)
Xgrid.append(np.floor((i-1)/10))
Ygrid.append((i-1)%10)
m = np.array(np.where(grid==i))
Num.append(m.shape[1])
if Num[i-1] > 0:
first_in_grid.append(m.min()) #Add the "+1" here so the indicies match those in the Keeling code
last_in_grid.append(m.max())
Max_Sus_grid.append(Suscept[m].max())
else:
first_in_grid.append(0)
last_in_grid.append(-1)
Max_Sus_grid.append(0)
#Work out grid to maximum grid transmission probabilities
from numpy import ndarray
MaxRate = ndarray((max(grid),max(grid)))
#Determine maximum number of animals to be infected in each grid square
# YOUR INDENTATION IS MESSED UP HERE (only 3 spaces) WHICH SUGGESTS THIS NEVER GETS CALLED
for i in range (1,int(max(grid))):
for j in range(1,int(max(grid))):
if ((i-1)==(j-1)) | (Num[i-1]==0) | (Num[j-1] == 0):
MaxRate[i-1,j-1] = np.inf
else:
Dist2 = (Size*max([0,(abs(Xgrid[i-1]-Xgrid[j-1])-1)])/10)**2+(Size*max([0,(abs(Ygrid[i-1]-Ygrid[j-1])-1)])/10)**2
MaxRate[i-1,j-1] = Max_Sus_grid[j-1]*Kernel(Dist2)
#Susceptible, Exposed, Infectious, Reported.==> latent period is 4 days
i=1
S=len(np.where(Status==0))
E=len(np.where(np.logical_and(Status>0, Status<=5)))
I=len(np.where(np.logical_and(Status>5, Status<=9)))
R=len(np.where(Status==10))
R2=len(np.where(Status>9))
CullSheep=0
CullCattle=0
i=i+1
IterateFlag=1
# THIS IS A BAD APPROACH; GROWING LISTS IS EXTREMELY INEFFICIENT
S=[]
E=[]
I=[]
R=[]
R2=[]
CullSheep=[]
CullCattle=[]
t=[]
t.append(0)
results = np.c_[np.array([1]*N),np.arange(1,N+1),np.array([0]*N)]
# NEVER SEPARATE STATEENTS WITH SEMICOLONS; VERY HARD TO READ AND DEBUG
# JUST USE & INSTEAD OF logical_and
while(np.logical_and(t[-1]