File: CustomModel.py

package info (click to toggle)
epigrass 2.5.0+dfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 10,980 kB
  • sloc: python: 4,742; makefile: 83
file content (324 lines) | stat: -rw-r--r-- 16,118 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
# This is a custom model to used in place of Epigrass' built-in models. Custom
# models must always be on a file named CustomModel.py and contain at least 
# a function named Model. Both the File name and the function Names are case-sensitive,
# so be careful. Please refer to the manual for intructions on how to write your 
# own custom models.
from numpy.random import poisson,  multinomial
from numpy import *


def Model(inits, simstep, totpop, theta=0, npass=0, bi=None, bp=None, values=None):
        """
        Calculates the model SIR, and return its values.
        - inits = (E,I,S)
        - par = (Beta, alpha, E,r,delta,B, w, p) see docs.
        - theta = infectious individuals from neighbor sites (CONSIDERANDO QUE CHEGUE COMO UMA lista
            'Ip0','Ip2','Ip10','Ip15','Ip20','Ip40','Is0','Is2','Is10','Is15','Is20','Is40'
        - npass deve ser tb um vetor com populacao por faixa etaria
        """
        #Initializing

        vnames = ('S0','S2','S10','S15','S20','S40',
                'Ep0','Ep2','Ep10','Ep15','Ep20','Ep40',                  
                'Ip0','Ip2','Ip10','Ip15','Ip20','Ip40',
                'Es0','Es2','Es10','Es15','Es20','Es40',
                'Is0','Is2','Is10','Is15','Is20','Is40',
                'Rmax0','Rmax2','Rmax10','Rmax15','Rmax20','Rmax40',
                'Rmed0','Rmed2','Rmed10','Rmed15','Rmed20','Rmed40',
                'Rmin0','Rmin2','Rmin10','Rmin15','Rmin20','Rmin40',
                'Vmax0','Vmax2','Vmax10','Vmax15','Vmax20','Vmax40',
                'Vmed0','Vmed2','Vmed10','Vmed15','Vmed20','Vmed40',
                'Vmin0','Vmin2','Vmin10','Vmin15','Vmin20','Vmin40')
 
        # A distribuicao etaria e sitio-especificas
        if simstep == 1: #get initial values
            S0,S2,S10,S15,S20,S40 = (self.bi['s0'],self.bi['s2'], self.bi['s10'],self.bi['s15'],self.bi['s20'],self.bi['s40'])
            Ep0,Ep2,Ep10,Ep15,Ep20,Ep40 = (self.bi['ep0'],self.bi['ep2'],self.bi['ep10'],self.bi['ep15'],self.bi['ep20'],self.bi['ep40'])
            Ip0,Ip2,Ip10,Ip15,Ip20,Ip40 = (self.bi['ip0'],self.bi['ip2'],self.bi['ip10'],self.bi['ip15'],self.bi['ip20'],self.bi['ip40'])
            Es0,Es2,Es10,Es15,Es20,Es40 = (self.bi['es0'],self.bi['es2'],self.bi['es10'],self.bi['es15'],self.bi['es20'],self.bi['es40'])
            Is0,Is2,Is10,Is15,Is20,Is40 = (self.bi['is0'],self.bi['is2'],self.bi['is10'],self.bi['is15'],self.bi['is20'],self.bi['is40'])
            Rmax0,Rmax2,Rmax10,Rmax15,Rmax20,Rmax40=(self.bi['rmax0'],self.bi['rmax2'],self.bi['rmax10'],self.bi['rmax15'],self.bi['rmax20'],self.bi['rmax40'])
            Rmed0,Rmed2,Rmed10,Rmed15,Rmed20,Rmed40=(self.bi['rmed0'],self.bi['rmed2'],self.bi['rmed10'],self.bi['rmed15'],self.bi['rmed20'],self.bi['rmed40'])
            Rmin0,Rmin2,Rmin10,Rmin15,Rmin20,Rmin40=(self.bi['rmin0'],self.bi['rmin2'],self.bi['rmin10'],self.bi['rmin15'],self.bi['rmin20'],self.bi['rmin40'])
            Vmax0,Vmax2,Vmax10,Vmax15,Vmax20,Vmax40=(self.bi['vmax0'],self.bi['vmax2'],self.bi['vmax10'],self.bi['vmax15'],self.bi['vmax20'],self.bi['vmax40'])
            Vmed0,Vmed2,Vmed10,Vmed15,Vmed20,Vmed40=(self.bi['vmed0'],self.bi['vmed2'],self.bi['vmed10'],self.bi['vmed15'],self.bi['vmed20'],self.bi['vmed40'])
            Vmin0,Vmin2,Vmin10,Vmin15,Vmin20,Vmin40=(self.bi['vmin0'],self.bi['vmin2'],self.bi['vmin10'],self.bi['vmin15'],self.bi['vmin20'],self.bi['vmin40'])

            # sem fluxo no primeiro dia
            Mp0,Mp2,Mp10,Mp15,Mp20,Mp40,Ms0,Ms2,Ms10,Ms15,Ms20,Ms40=0,0,0,0,0,0,0,0,0,0,0,0

            # calculando o dicionario de fluxos:
            self.parentSite.pdest = setPassDest(self)
            npass0, npass2, npass10, npass15, npass20, npass40 = npass
            
            
        else:
            S0,S2,S10,S15,S20,S40,Ep0,Ep2,Ep10,Ep15,Ep20,Ep40,Ip0,Ip2,Ip10,Ip15,Ip20,Ip40,Es0,Es2,Es10,Es15,Es20,Es40,Is0,Is2,Is10,Is15,Is20,Is40,Rmax0,Rmax2,Rmax10,Rmax15,Rmax20,Rmax40,Rmed0,Rmed2,Rmed10,Rmed15,Rmed20,Rmed40,Rmin0,Rmin2,Rmin10,Rmin15,Rmin20,Rmin40,Vmax0,Vmax2,Vmax10,Vmax15,Vmax20,Vmax40,Vmed0,Vmed2,Vmed10,Vmed15,Vmed20,Vmed40,Vmin0,Vmin2,Vmin10,Vmin15,Vmin20,Vmin40 = inits
            # obtendo os visitantes
            Mp0,Mp2,Mp10,Mp15,Mp20,Mp40,Ms0,Ms2,Ms10,Ms15,Ms20,Ms40 = self.parentSite.infectedvisiting
            npass0,npass2,npass10,npass15,npass20,npass40 = npass
            # resetando os visitantes:
            self.parentSite.infectedvisiting=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
        
        totpass=npass0+npass2+npass10+npass15+npass20+npass40
        N = totpop
        # Parametros
        #beta,alpha,e,r,delta,B,w,p = par
        
        pp, ps, sp, ss = (self.bp['pp'],self.bp['ps'],self.bp['sp'],self.bp['ss'])
        pv, e, gp, gs, a, r = (self.bp['pv'],self.bp['e'],self.bp['gp'],self.bp['gs'],self.bp['a'],self.bp['r'],) 
        ep, es = (self.bp['ep'],self.bp['es'])
        c0, c1, c2, c3, c4, c5= (self.bp['c0'],self.bp['c1'],self.bp['c2'],self.bp['c3'],self.bp['c4'],self.bp['c5'])
        n0, n2, n10, n15, n20, n40 = (self.bp['n0'], self.bp['n2'],self.bp['n10'],self.bp['n15'],self.bp['n20'],self.bp['n40'])

        #Vacination event  # apenas na primeira faixa etaria
        if self.parentSite.vaccineNow:
            S0 -= self.parentSite.vaccov*S0
            Vmax0 += self.parentSite.vaccov*S0
        
        # numero de pessoas com infeccao primaria (locais e visitas)
        Ip=array([Ip0+Mp0, Ip2+Mp2, Ip10+Mp10, Ip15+Mp15, Ip20+Mp20, Ip40+Mp40])
        Is=array([Is0+Ms0, Is2+Ms2, Is10+Ms10, Is15+Ms15, Is20+Ms20, Is40+Ms40])
        Iptot = sum(Ip)
        Istot = sum(Is)
        
        # distribuicao dos contatos
        matriz = 'delvalle'
        if (matriz =='luz'):
            mc0 = array([c0, c0, c3, c5, c5, c3])
            mc2 = array([c0, c5, c0, c3, c5, c3])
            mc10 =array([c0, c0, c5, c3, c5, c3])
            mc15=array([c0, c0, c3, c5, c5, c3])
            mc20=array([c5, c3, c0, c0, c5, c3])
            mc40=array([c0, c3, c3, c0, c5, c5])
        
        if (matriz =='delvalle'):
            
            mc0 = array([c0, c2, c2, c1, c1, c0])
            mc2 = array([c2, c5, c5, c1, c1, c0])
            mc10 =array([c2, c5, c5, c1, c1, c0])
            mc15=array([c1, c1, c3, c3, c2, c1])
            mc20=array([c1, c1, c0, c2, c5, c2])
            mc40=array([c0, c0, c1, c2, c2, c1])
        
        # probabilidade de infeccao por contato (a divisao e para padronizar os vetores de contato)
        x = (1/(N+totpass))*(pp*Ip+ps*Is)
        p0 = sum(x*mc0/sum(mc0))
        p2 = sum(x*mc2/sum(mc2))
        p10 = sum(x*mc10/sum(mc10))
        p15 = sum(x*mc15/sum(mc15))
        p20 = sum(x*mc20/sum(mc20))
        p40 = sum(x*mc40/sum(mc40))
        
       
        # calculando a incidencia
        
        Lp0 = min(poisson(sp*S0*(1-(1-p0)**n0), 1)[0], S0)
        Ls0r = min(poisson(ss*Rmin0*(1-(1-p0)**n0), 1)[0], Rmin0)
        Ls0v = min(poisson(ss*Vmin0*(1-(1-p0)**n0), 1)[0], Vmin0)
        
        Lp2 = min(poisson(sp*S2*(1-(1-p2)**n0), 1)[0], S2)
        Ls2r = min(poisson(ss*Rmin2*(1-(1-p2)**n2), 1)[0], Rmin2)
        Ls2v = min(poisson(ss*Vmin2*(1-(1-p2)**n2), 1)[0], Vmin2)
        
        Lp10 = min(poisson(sp*S10*(1-(1-p10)**n10), 1)[0], S10)
        Ls10r = min(poisson(ss*Rmin10*(1-(1-p10)**n10), 1)[0], Rmin10)
        Ls10v = min(poisson(ss*Vmin10*(1-(1-p10)**n10), 1)[0], Vmin10)
        
        Lp15 = min(poisson(sp*S15*(1-(1-p15)**n15), 1)[0], S15)
        Ls15r = min(poisson(ss*Rmin15*(1-(1-p15)**n15), 1)[0], Rmin15)
        Ls15v = min(poisson(ss*Vmin15*(1-(1-p15)**n15), 1)[0], Vmin15)
        
        Lp20 = min(poisson(sp*S20*(1-(1-p20)**n20), 1)[0], S20)
        Ls20r = min(poisson(ss*Rmin20*(1-(1-p20)**n20), 1)[0], Rmin20)
        Ls20v = min(poisson(ss*Vmin20*(1-(1-p20)**n20), 1)[0], Vmin20)
        
        Lp40 = min(poisson(sp*S40*(1-(1-p40)**n40), 1)[0], S40)
        Ls40r = min(poisson(ss*Rmin40*(1-(1-p40)**n40), 1)[0], Rmin40)
        Ls40v = min(poisson(ss*Vmin40*(1-(1-p40)**n40), 1)[0], Vmin40)
            
        Lposp = Lp0+Lp2+Lp10+Lp15+Lp20+Lp40
        Lposs = Ls0r+Ls2r+Ls10r+Ls15r+Ls20r+Ls40r+Ls0v+Ls2v+Ls10v+Ls15v+Ls20v+Ls40v
        Lpos = Lposp+Lposs
        
        self.parentSite.totalcases += Lpos #update number of cases (DA PARA FAZER MAIS DETALHADO?)
        
        # Model (vacina entra por fora).  e sem progressao etaria! tirei a infeccao secundaria tb porque estava errado, e preciso repensar esta parte : a entrada do Ls no modelo, acrescentei classe com latencia
            
        S0pos = S0 -  Lp0
        Ep0pos = (1-ep) * Ep0 + Lp0
        Ip0pos = (1-gp) * Ip0 + ep*Ep0
        Rmax0pos = (1-a) * Rmax0 + gp * Ip0 + gs * Is0
        Rmed0pos = (1-a) * Rmed0 + a * Rmax0
        Rmin0pos = Rmin0 + a * Rmed0 - Ls0r
        Es0pos = (1-es)*Es0 + Ls0r + Ls0v 
        Is0pos = (1-gs)*Is0 + es*Es0
        Vmax0pos = (1- r) * Vmax0
        Vmed0pos = (1-r) * Vmed0 + r*Vmax0
        Vmin0pos =  Vmin0 + r*Vmed0 - Ls0v
        
        S2pos = S2 -  Lp2
        Ep2pos = (1-ep) * Ep2 + Lp2
        Ip2pos = (1-gp) * Ip2 + ep*Ep2
        Rmax2pos = (1-a) * Rmax2 + gp * Ip2 + gs * Is2
        Rmed2pos = (1-a) * Rmed2 + a * Rmax2
        Rmin2pos = Rmin2 + a * Rmed2 - Ls2r
        Es2pos = (1-es)*Es2 + Ls2r + Ls2v 
        Is2pos = (1-gs)*Is2 + es*Es2
        Vmax2pos = (1- r) * Vmax2
        Vmed2pos = (1-r) * Vmed2 + r*Vmax2
        Vmin2pos =  Vmin2 + r*Vmed2 - Ls2v
        
        S10pos = S10 -  Lp10
        Ep10pos = (1-ep) * Ep10 + Lp10
        Ip10pos = (1-gp) * Ip10 + ep*Ep10
        Rmax10pos = (1-a) * Rmax10 + gp * Ip10 + gs * Is10
        Rmed10pos = (1-a) * Rmed10 + a * Rmax10
        Rmin10pos = Rmin10 + a * Rmed10 - Ls10r
        Es10pos = (1-es)*Es10 + Ls10r + Ls10v 
        Is10pos = (1-gs)*Is10 + es*Es10
        Vmax10pos = (1- r) * Vmax10
        Vmed10pos = (1-r) * Vmed10 + r*Vmax10
        Vmin10pos =  Vmin10 + r*Vmed10 - Ls10v
        
        S15pos = S15 -  Lp15
        Ep15pos = (1-ep) * Ep15 + Lp15
        Ip15pos = (1-gp) * Ip15 + ep*Ep15
        Rmax15pos = (1-a) * Rmax15 + gp * Ip15 + gs * Is15
        Rmed15pos = (1-a) * Rmed15 + a * Rmax15
        Rmin15pos = Rmin15 + a * Rmed15 - Ls15r
        Es15pos = (1-es)*Es15 + Ls15r + Ls15v 
        Is15pos = (1-gs)*Is15 + es*Es15
        Vmax15pos = (1- r) * Vmax15
        Vmed15pos = (1-r) * Vmed15 + r*Vmax15
        Vmin15pos =  Vmin15 + r*Vmed15 - Ls15v
        
        S20pos = S20 -  Lp20
        Ep20pos = (1-ep) * Ep20 + Lp20
        Ip20pos = (1-gp) * Ip20 + ep*Ep20
        Rmax20pos = (1-a) * Rmax20 + gp * Ip20 + gs * Is20
        Rmed20pos = (1-a) * Rmed20 + a * Rmax20
        Rmin20pos = Rmin20 + a * Rmed20 - Ls20r
        Es20pos = (1-es)*Es20 + Ls20r + Ls20v 
        Is20pos = (1-gs)*Is20 + es*Es20
        Vmax20pos = (1- r) * Vmax20
        Vmed20pos = (1-r) * Vmed20 + r*Vmax20
        Vmin20pos =  Vmin20 + r*Vmed20 - Ls20v
        
        S40pos = S40 -  Lp40
        Ep40pos = (1-ep) * Ep40 + Lp40
        Ip40pos = (1-gp) * Ip40 + ep*Ep40
        Rmax40pos = (1-a) * Rmax40 + gp * Ip40 + gs * Is40
        Rmed40pos = (1-a) * Rmed40 + a * Rmax40
        Rmin40pos = Rmin40 + a * Rmed40 - Ls40r
        Es40pos = (1-es)*Es40 + Ls40r + Ls40v 
        Is40pos = (1-gs)*Is40 + es*Es40
        Vmax40pos = (1- r) * Vmax40
        Vmed40pos = (1-r) * Vmed40 + r*Vmax40
        Vmin40pos =  Vmin40 + r*Vmed40 - Ls40v

           
        # Updating stats
        self.parentSite.incidence.append(Lpos)
        # Raises site infected flag and adds parent site to the epidemic history list.
        if not self.parentSite.infected: 
            if Lpos > 0:
                print 'infected', self.parentSite.sitename
                #if not self.parentSite.infected:
                self.parentSite.infected = self.parentSite.parentGraph.simstep
                self.parentSite.parentGraph.epipath.append((self.parentSite.parentGraph.simstep,self.parentSite,self.parentSite.infector))
                #print self.parentSite.ftheta
        #infectious
        Ipos = Ip0pos+Is0pos+Ip2pos+Is2pos+Ip10pos+Is10pos+Ip15pos+Is15pos+Ip20pos+Is20pos+Ip40pos+Is40pos
        self.parentSite.migInf.append(Ipos)
        
        # send infected to others! 
        if (Ipos>0):
            tot0 = (S0pos+Ep0pos+Es0pos+Ip0pos+ Is0pos+Rmax0pos+Rmed0pos+Rmin0pos+Vmax0pos+Vmed0pos+Vmin0pos)
            tot2 = (S2pos+Ep2pos+Es2pos+Ip2pos+ Is2pos+Rmax2pos+Rmed2pos+Rmin2pos+Vmax2pos+Vmed2pos+Vmin2pos)
            tot10 = (S10pos+Ep10pos+Es10pos+Ip10pos+ Is10pos+Rmax10pos+Rmed10pos+Rmin10pos+Vmax10pos+Vmed10pos+Vmin10pos)
            tot15 = (S15pos+Ep15pos+Es15pos+Ip15pos+ Is15pos+Rmax15pos+Rmed15pos+Rmin15pos+Vmax15pos+Vmed15pos+Vmin15pos)
            tot20 = (S20pos+Ep20pos+Es20pos+Ip20pos+ Is20pos+Rmax20pos+Rmed20pos+Rmin20pos+Vmax20pos+Vmed20pos+Vmin20pos)
            tot40 = (S40pos+Ep40pos+Es40pos+Ip40pos+ Is40pos+Rmax40pos+Rmed40pos+Rmin40pos+Vmax40pos+Vmed40pos+Vmin40pos)
            
            infmig = [npass0*(Ip0pos/tot0),npass2*(Ip2pos/tot2),npass10*(Ip10pos/tot10),npass15*(Ip15pos/tot15),npass20*(Ip20pos/tot20),npass40*(Ip40pos/tot40),npass0*(Is0pos/tot0),npass2*(Is2pos/tot2),npass10*(Is10pos/tot10),npass15*(Is15pos/tot15),npass20*(Is20pos/tot20),npass40*(Is40pos/tot40)]
            sendInfectedPass(self, infmig, self.parentSite.pdest)
            #print self.parentSite.sitename,  self.parentSite.ts

        return [S0pos,S2pos,S10pos,S15pos,S20pos,S40pos,
                Ep0pos,Ep2pos,Ep10pos,Ep15pos,Ep20pos,Ep40pos,
                Ip0pos,Ip2pos,Ip10pos,Ip15pos,Ip20pos,Ip40pos,
                Es0pos,Es2pos,Es10pos,Es15pos,Es20pos,Es40pos,
                Is0pos,Is2pos,Is10pos,Is15pos,Is20pos,Is40pos,
                Rmax0pos,Rmax2pos,Rmax10pos,Rmax15pos,Rmax20pos,Rmax40pos,
                Rmed0pos,Rmed2pos,Rmed10pos,Rmed15pos,Rmed20pos,Rmed40pos,
                Rmin0pos,Rmin2pos,Rmin10pos,Rmin15pos,Rmin20pos,Rmin40pos,
                Vmax0pos,Vmax2pos,Vmax10pos,Vmax15pos,Vmax20pos,Vmax40pos,
                Vmed0pos,Vmed2pos,Vmed10pos,Vmed15pos,Vmed20pos,Vmed40pos,
                Vmin0pos,Vmin2pos,Vmin10pos,Vmin15pos,Vmin20pos,Vmin40pos]
                

def getNPass(self):
    '''
    Calculando quantos passageiros chegam ao site por dia
    '''
    npass = [0, 0, 0, 0, 0, 0]
    fes = [0, 2, 10, 15, 20, 40]
    for e in self.parentSite.getInEdges():
        fe = fes.index(e.length)
        npass[fe]+=e.fmig
    
    return npass
    
def setPassDest(self):
    '''
    Define, para cada ZT, os destinos possiveis e a probabilidade de cada um deles por faixa etaria
    '''
    
    # criando uma lista de dicts, um para cada FE, contendo os outedges e of fluxos
    dests = [{}, {}, {}, {}, {}, {}]
    fes = [0, 2, 10, 15, 20, 40]
    
    for e in self.parentSite.getOutEdges():
        fe = fes.index(e.length)
        dests[fe][e]=e.fmig
    
    # transformando os fluxos em probabilidade de destino
    pdests = dests
    for i in xrange(6):
        tot=sum(dests[i].values())
        di = dests[i]
        for k, v in di.iteritems():
            pdests[i][k]=v/tot
            
    
            
    return pdests

    
def sendInfectedPass(self, infgroup, pdests):
    '''
    if there is people infected, check where they are going to...
    inf = [Ip0pos,Ip2pos,Ip10pos,Ip15pos,Ip20pos,Ip40pos,Is0pos,Is2pos,Is10pos,Is15pos,Is20pos,Is40pos]
    '''
    # Quantos vao sair?
    
    a = [0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5]
    if (pdests):
        n = 0 # indexador do inf 
        for i in infgroup:
            posdest = a[n]
            #print 'posdest=', posdest
            #print 'enviando %s para os destinos:'%i , pdests[posdest]
            if (int(i) & len(pdests[posdest])):
                destlist = pdests[posdest]
                #print 'destinos gerais de %s:'%a[n], destlist
                choosedest = multinomial(int(i),destlist.values())
                k=0
                for j in destlist: # enviando os infectados para seus destinos
                    j.dest.infectedvisiting[n] += choosedest[k]
                    k+=1
            n+=1

    return 1