Commit 97ed4e53 by Chris Jewell

### Change to kernel function. Now makes sense mathematically(!)

parent ee51c48d
 ... ... @@ -8,7 +8,8 @@ import matplotlib.pyplot as plt def hhNetwork(numHouseHolds, numPerHouse): """Fake household network. Returns a block diagonal matrix representing people clustered into numHouseHolds each of size numPerHouse""" nTotal = numHouseHolds * numPerHouse net = np.zeros(shape=[nTotal,nTotal]) for i in range(numHouseHolds): ... ... @@ -22,10 +23,11 @@ def hhNetwork(numHouseHolds, numPerHouse): class K: def __init__(self,net): """Represents a network""" self.net = net def __call__(self, beta0, beta1, idx): return beta0 + beta1*self.net[idx] # Background plus pairwise infection on network def __call__(self, beta1, idx): return beta1*self.net[idx] @property def shape(self): ... ... @@ -34,6 +36,12 @@ class K: class Simulation: def __init__(self, simmatrix, T, par, kernel): """Simulation represents a simulation of an epidemic. :param simmatrix: the nTxnS matrix with nT time steps and nS states :param T: a nT length vector of times :param par: the parameters used in the simulation :param kernel: the kernel matrix used.""" self.history = simmatrix self.T = T self.par = par ... ... @@ -72,6 +80,7 @@ class Simulation: ax.set_xlabel('Time/days') ax.set_ylabel('Household id') class Simulator: def __init__(self, kernel, timestep=1): self.kernel = kernel ... ... @@ -79,10 +88,10 @@ class Simulator: self.n = kernel.shape[0] self.par = {'alpha': 0., 'beta': 0., 'gamma': 0.} def calcKernel(self,par, idx): # Memoisation -- only recalculate if parameters change def calcKernel(self, par, idx): # Memoisation -- only recalculate if parameters change if (par['alpha'] != self.par['alpha']) | (par['beta'] != self.par['beta']): self.k = self.kernel(par['alpha'], par['beta'], idx) return self.k self.lambdaj = par['alpha'] + self.kernel(par['beta'], idx).sum(axis=0) return self.lambdaj # Infection rate equation def propagate(self, par, state0, nTimeSteps=1): ... ... @@ -93,8 +102,7 @@ class Simulator: I = np.where(state0 == 1)[0] # S->I k = self.calcKernel(par, np.ix_(I,S)) # pull relevant rows/cols from kernel lambdaj = k.sum(axis=0) # sum across rows lambdaj = self.calcKernel(par, np.ix_(I,S)) # pull relevant rows/cols from kernel StoI = None if S.size > 0: # draw who gets infected pSI = 1-np.exp(-lambdaj*self.timestep) # Rate -> probability ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!