/
polioCaseLight.m
291 lines (239 loc) · 13.1 KB
/
polioCaseLight.m
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
function [newCases, activeCases, cumulativeCases ] = polioCaseLight(Indv,nHhold, CaseIntroPopulationType, caseIntroSites, nbCases, nbDays)
%this is a light version of polioCase code with fewer outputs
% newCases is a vector of length nbDays+1. newCases(i) indicates how many new infected cases we have on day i
% activeCases vector of length nbDays+1, where activeCases(i) indicates how many active infecting cases we have on day i
% cumulativeCases is a vector indicating in each day the cumulative number
% of new infected cases
% Indv represents studied country people. it is generated using the
% following files: poliomodel_Leb.m, poliomodel_Syr.m, and poliomodel_Pal.m
% -- see README.txt for the generation of Indv and nHhold parameter.
% CaseIntroPopulationType, caseIntroSites, nbCases, nbDays
%CaseIntroPopulationType is a vector representing the population type of the introduced infected case. the values in the vector are the integers 1,2,or 3 indicating respectively
%lebanese, syrian, palestinian
%caseIntroSites is a vector of integers 1,2,3,4 indicating respectively
%Beirut, Bekaa, North, and South, same length as number of cases. this
%parameter indicate where each infected case entered.
% nbCases is an integer indicating how many infected cases are introduced
% nbDays indicate on how many days we are studying the propagation of
% infection after the introduction of the infected cases.
symptomDist=[zeros(50,1);ones(5,1);zeros(45,1);]; % Symtpoms probability 0.05
paralysisDist=[zeros(595,1); ones(5,1);zeros(400,1)]; % paralysis probability 0.005
mu=[2 10 35 7 0.1 2 2 1 1];
RangL=[0.1 0.1 20 3 0 1 1 0 0];
RangH=[7 20 50 20 0.5 5 4 3 3];
F{1}=[];
for i=1:length(mu)
Fa= mu(i)+(mu(i)-RangL(i))/2 *randn(1000,1);
Fb=mu(i)+(RangH(i)-mu(i))/2 *randn(1000,1);
inda=find(Fa<mu(i) & Fa>=0);
indb=find(Fb>=mu(i));
F{i}=[Fa(inda);Fb(indb)]; % generate skew normal distributions of mode mu, and lower boundary
% RangL, and higher boundary RangH for each of the parameters indicated at
% the end of the code.
end
tempTrackIndv1=[];% 2-D matrix of 7 columns indicating for each new infected case the following criteria:
% popType, Site, householdType, which household , which person, isAdult, has access
tempTrackIndv2=[];
% 2-D matrix of 7 columns indicating for each new infected case the following
% criteria: days of infection, latency period, incubation , disease time, symptomsFlag, paralysis flag, detected
for i=1:nbCases
sympFlag=randsample(symptomDist,1) ;
ParalyzedFlag=randsample(paralysisDist,1);
if ParalyzedFlag==1
sympFlag=1;
end
tempTrackIndv1=[tempTrackIndv1; CaseIntroPopulationType(i) caseIntroSites(i) -1 0 1 0 0]; % the introduced cases are usually children
tempTrackIndv2=[tempTrackIndv2; 0 round(randsample(F{1},1)) round(randsample(F{2},1)) round(randsample(F{3},1)) sympFlag ParalyzedFlag 0];
end
communicationPool{1}=[ones(1,4) 2 3]; %each population is more likely to communicate within itself than with others by a ratio of 2/3 internally, and 1/3 with other populations
communicationPool{2}=[ones(1,4)*2 1 3];
communicationPool{3}=[ones(1,4)*3 2 1];
sitesPool{1}=[2 3 4];%each infected individual is more likely to infect people within the same site, but if he infects from other sites the infection probability is 1/3 for each of the three other sites
sitesPool{2}=[1 3 4];
sitesPool{3}=[1 2 4];
sitesPool{4}=[1 2 3];
trackIndv=[];
trackDays=[];
tempIm=0;
tempAdult=0;
tempHasAccess=0;
rateInfection=0;
randInf=0;
trackIndv=[trackIndv;tempTrackIndv1];
trackDays=[trackDays;tempTrackIndv2(1,1)*ones(length(tempTrackIndv2(:,1)),1) tempTrackIndv2];
fprintf('Day %d ***\n',tempTrackIndv2(1,1));
fprintf('Cumulative number of cases %d ***\n',length(tempTrackIndv2(:,1)));
activeCases(1)=0;
newCases(1)=nbCases;
cumulativeCases(1)=sum(newCases);
while (nbDays>tempTrackIndv2(1,1))
tempTrackIndv2(:,1)=tempTrackIndv2(:,1)+1;
N1=length(tempTrackIndv1(:,1));
j1=tempTrackIndv2(1,1);
Sites(j1+1,:)=zeros(1,4);
populations(j1+1,:)=zeros(1,3);
nbAdults(j1+1)=0;
nbChildren(j1+1)=0;
fprintf('Day %d ***\n',tempTrackIndv2(1,1));
for i=1:N1
%if the patient is within the incubation time, before the end of the
%disease , and if not paralysed he continues to infect other people.
flagInfect=0;
if (tempTrackIndv2(i,1)<tempTrackIndv2(i,2) )
continue; %if number of infected days< latency period he doesnt infect anyone
elseif ((tempTrackIndv2(i,1)==tempTrackIndv2(i,2)) &&(i>nbCases)) % after latency period he infects his household if they are not immuned
for i1=1:tempTrackIndv1(i,3)
if i1~=tempTrackIndv1(i,5)%if not same individuel
tempIm= Indv(tempTrackIndv1(i,2),tempTrackIndv1(i,1),tempTrackIndv1(i,3),tempTrackIndv1(i,4),i1,2);
if tempIm<=0
if (tempIm==0)%if partially immuned
rateInfection=randsample(F{5},1);
partialPool=[ones(floor(rateInfection*100),1); zeros(floor((1-rateInfection)*100),1)];
if (~isempty(partialPool))
randInf=randsample(partialPool,1);
end
end
end
if (randInf ||tempIm<0)% if not immuned or if partially immuned and was picked up to be infected
tempAdult=Indv(tempTrackIndv1(i,2),tempTrackIndv1(i,1),tempTrackIndv1(i,3),tempTrackIndv1(i,4),i1,1);
tempHasAccess=Indv(tempTrackIndv1(i,2),tempTrackIndv1(i,1),tempTrackIndv1(i,3),tempTrackIndv1(i,4),i1,3);
if isempty(tempHasAccess)
tempHasAccess=0;
end
tempPop=tempTrackIndv1(i,1);
tempSite=tempTrackIndv1(i,2);
tempHholdType=tempTrackIndv1(i,3);
tempHhold=tempTrackIndv1(i,4);
tempInd=tempTrackIndv1(i,5);
tempRow=[tempPop tempSite tempHholdType tempHhold tempInd tempAdult tempHasAccess];
[ ind1,loc1]=ismember(tempRow,tempTrackIndv1,'rows');%if he was not already infected add him to the infected list
if (ind1==0)
tempTrackIndv1=[tempTrackIndv1;tempRow];
sympFlag=randsample(symptomDist,1) ;
ParalyzedFlag=randsample(paralysisDist,1);
if ParalyzedFlag==1
sympFlag=1;
end
tempTrackIndv2=[tempTrackIndv2; 0 round(randsample(F{1},1)) round(randsample(F{2},1)) round(randsample(F{3},1)) sympFlag ParalyzedFlag 0];
end
end
else
continue;
end
end
elseif (tempTrackIndv2(i,7)==1)%detected
continue;%if detected he cannot infect other people
elseif (tempTrackIndv2(i,1)>tempTrackIndv2(i,3) && tempTrackIndv2(i,1)<tempTrackIndv2(i,4))
if (tempTrackIndv2(i,6)==1)%paralizedFlag
continue;%if paralyzed he cannot infect other people
elseif (tempTrackIndv2(i,5)==1 && tempTrackIndv1(i,7)==1)%symptomatic and has access to medical care
continue;%if symptomatic and has access to medcare he cannot infect other people
else
flagInfect=1;
end
%1- if symptomatic and has access to medical care detected and stops infecting
%2- if symptomatic and paralysis stops infecting
%3- else he contiues to infect
elseif (flagInfect==1 || (tempTrackIndv2(i,1)>tempTrackIndv2(i,2) && tempTrackIndv2(i,1)<tempTrackIndv2(i,3)))
%here the person has infectivity=1 and symptoms=0 and detection=0
flagInfect=0;
if tempTrackIndv1(i,6)==1 % check if adult
nbOnSite= round(randsample(F{7},1));%same province number of people to meet if u r adult
nbOutSite =round(randsample(F{9},1));%outside province number of people to meet if u r adult
else
nbOnSite= round(randsample(F{6},1));%same province number of people to meet if u r child
nbOutSite =round(randsample(F{8},1));%outside province number of people to meet if u r child
end
tempSite=tempTrackIndv1(i,2);
for k=1:nbOnSite+nbOutSite
tempPop= randsample(communicationPool{tempTrackIndv1(i,1)},1);
if (k>nbOnSite)
tempSite= randsample(sitesPool{tempTrackIndv1(i,2)},1);
end
tempHholdType = randi([1 8],1,1);
N=0;
N= nHhold(tempSite,tempPop,tempHholdType);
tempHhold = randi([1 N],1,1);
tempInd = randi([1 tempHholdType],1,1);
flagImmune=0;
tempAdult=0;
flagImmune=Indv(tempSite,tempPop,tempHholdType,tempHhold,tempInd,2);
tempAdult=Indv(tempSite,tempPop,tempHholdType,tempHhold,tempInd,1);
tempHasAccess=Indv(tempSite,tempPop,tempHholdType,tempHhold,tempInd,3);
if flagImmune>0
continue;
end
randInf=0;
if flagImmune==0
rateInfection=randsample(F{5},1);
partialPool=[ones(floor(rateInfection*100),1); zeros(floor((1-rateInfection)*100),1)];
if (~isempty(partialPool))
randInf=randsample(partialPool,1);
end
end
if randInf==1 || flagImmune<0
if isempty(tempHasAccess)
tempHasAccess=0;
end
tempRow=[tempPop tempSite tempHholdType tempHhold tempInd tempAdult tempHasAccess];
[ ind1,loc1]=ismember(tempRow,tempTrackIndv1,'rows');
if (ind1==0)
tempTrackIndv1=[tempTrackIndv1;tempRow];
sympFlag=randsample(symptomDist,1) ;
ParalyzedFlag=randsample(paralysisDist,1);
if ParalyzedFlag==1
sympFlag=1;
end
tempTrackIndv2=[tempTrackIndv2; 0 round(randsample(F{1},1)) round(randsample(F{2},1)) round(randsample(F{3},1)) sympFlag ParalyzedFlag 0];
end
end
end
elseif ((tempTrackIndv2(i,1)>=tempTrackIndv2(i,4))&& (i>nbCases))
%the person becomes immuned after finishing the disease time
Indv(tempSite,tempPop,tempHholdType,tempHhold,tempInd,2)=1;
continue;
end
end
trackIndv=tempTrackIndv1;
trackDays=[tempTrackIndv2(1,1)*ones(length(tempTrackIndv2(:,1)),1) tempTrackIndv2];
fprintf('Cumulative number of cases %d ***\n',length(tempTrackIndv2(:,1)));
index1=find(trackDays(:,1)==tempTrackIndv2(1,1),1,'last');
tempActive=[];
for j=1:length(tempTrackIndv2(:,1))
if ((tempTrackIndv2(j,1)>=tempTrackIndv2(j,2)) &&(tempTrackIndv2(j,1)<=tempTrackIndv2(j,4)))
tempActive=[tempActive;j];
end
end
activeCases(j1+1)=length(tempActive);
[index2]=find(trackDays(1:index1,2)==0);
newCases(j1+1)=numel(index2);
cumulativeCases(j1+1)=sum(newCases);
end
end
%%
%%
% Feature Mean RangL RangH
% 1-Latency period 2 0.1 7 # infected with time
% 2-Incubation 10 0.1 20 # detected with time Symp=1 Access=1
% 3-Disease (Not Immune) 35 20 50 # paralyzed with time
% 4-Disease (Partial Immune) 7 3 20 Day of alarm
% Household infection rate
% Immune 0 0 0
% 5-Partial 0.1 0 0.5
% Non 1 1 1
% Contacts in province
% 6-child 2 1 10 67% same country, 33% other country..
% 7-Adult 2 1 4
% Contacts outside province
% 8-child 1 0 3
% 9-Adult 1 0 3
% Symtpoms probability 0.05
% paralysis 0.005
% Access to Medcare
% Syrian 75%
% Pal 75%
% Leb 98%
%%
%if paralyzed he stop infecting people. and paralysis may happen after the
%incubation period, not before it
%%