/* Simulation of a parasitic disease in C */
* Use “select all”, “copy” and “paste” to bring this program into your environment.
* I am President of the Libertarian Programming Association: I have the freedom to write programs of catastrophical style which I do not understand once I have written and used them. En bref: Do not ask me why I add 1 in line i and not in line j.
/* Mathematical specification: modell SNM*, F1 Poisson; output is the distribution of parasites per individual
The program consists of the following parts:
Head including Imports
Declarations
Description of variables and parameters
Initialisations and initial values
Main including Calculation of elapsed time
Calculation of the type of event
Readjusting of the states
End including Calculation of the output
Output */
/*-----*/
/*-Head-*/
/*-----*/
/*-Imports-*/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
main()
{
int random(); /* used to generate U(0,1), number between 0 and 2E31-1 */
/*-Declarations-*/
double tts,tr,dt,tnp,at,h,lambda,maxlife,maxsl,mlm,my,kappa;
double r1,r2,r3,r4,r5,r6;
double hsts,sts,sum1,sum2,t,t1,theta;
double d[101],nschisto[101],thuman[200];
int te,M,i,j,k,kadenz,msm,no,shuman[200];
/* -Description of variables and parameters-*/
/*
at Actual Time that the program has run for
d[i] proportion of individuals with i worms (Distribution)
dt time until the next event happens (math. small time unit)
h local aid for calculations (Help)
hsts local aid for calculation (Help for Successive Test Summation)
i locally used index
j locally used index
k locally used index
kappa death rate of individuals ("kappa" see section 1.4)
lambda contact rate of humans ("lambda" see section 1.4)
M total number of individuals (numbered 0 to (M-1), "M" see section 1.4)
maxlife MAXimal age (LIFE) ever achieved by a person
maxsl MAXimal load of parasites ever achieved by a person (MAXimal Schistosome Load)
mlm Maximal age that can be registered in this program (Maximal Life Measurement)
msm Maximal load of parasites that can be registered in this program (Max.Sch.Meas.)
my death rate of parasites ("my" see section 1.4)
no Number of person that gets newly infected
nschisto[p] Number of people with p parasites (Number of SCHISTOsomes)
r1 - r6 random numbers (more information where used), realisations of "random()")
shuman[t] number of parasites in person t (Schistosomes in HUMAN)
sum1 total rate at which parasites die in the whole system
sum2 total rate at which people get infected
sts Sucessive Test Summation
t local aid for calculations
t1 local aid for calculations
te variable telling which Type of Event happened
theta "multiplicator" of the infection ("theta" see section 1.4)
thuman[j] age of person j (Time of HUMAN)
tnp actual Total Number of Parasites
tr Total Rate at which an event happens
tts Total Time until the end of the Simulation */
/* -Initialisations and initial values-*/
/* The following initialisations and initial values must be consistent. There are no built in test for their consistency. "Garbage in, garbage out!" */
/* Initialisations concerning the individuals: */
M = 200;
tnp = 224;
maxlife = 0; /* is going to be corrected immediately */
maxsl = 5;
mlm = 5000;
msm = 100; /* do not enlarge because of calculation of faculty (!) */
for (i=0; i<50; ++i) {shuman[i]=1;}
for (i=50; i<60; ++i) {shuman[i]=2;}
for (i=60; i<75; ++i) {shuman[i]=3;}
for (i=75; i<91; ++i) {shuman[i]=4;}
for (i=91; i<100; ++i) {shuman[i]=5;}
for (i=100; i<M; ++i) {shuman[i]=0;}
nschisto[0]=100;nschisto[1]=50;nschisto[2]=10;nschisto[3]=15;nschisto[4]=16;
nschisto[5]=9;
for (i=(maxsl+1); i<=msm; ++i) {nschisto[i]=0;}
/* The following initialisations lead to a age distribution amongst the 200 individuals which is a "typical" realisation if there are 200 individuals each with a lifetime modelled with a exponential random variable of rate 1/40. */
thuman[0]=0;thuman[199]=0;thuman[55]=0;thuman[77]=0;thuman[109]=0;
thuman[110]=0;thuman[2]=1;thuman[80]=1;thuman[26]=1;thuman[59]=1;
thuman[76]=1;thuman[98]=1;thuman[51]=2;thuman[35]=2;thuman[160]=2;
thuman[50]=2;thuman[91]=2;thuman[177]=2;thuman[158]=3;thuman[44]=3;
thuman[122]=3;thuman[61]=3;thuman[86]=3;thuman[11]=3;thuman[194]=4;
thuman[167]=4;thuman[139]=4;thuman[100]=4;thuman[102]=4;thuman[66]=4;
thuman[106]=5;thuman[147]=5;thuman[5]=5;thuman[23]=5;thuman[174]=5;
thuman[198]=5;thuman[171]=6;thuman[170]=6;thuman[183]=6;thuman[120]=6;
thuman[46]=6;thuman[124]=6;thuman[96]=7;thuman[72]=7;thuman[71]=7;
thuman[128]=7;thuman[187]=7;thuman[176]=7;thuman[89]=8;thuman[191]=8;
thuman[73]=8;thuman[81]=8;thuman[95]=8;thuman[161]=8;thuman[92]=9;
thuman[54]=9;thuman[56]=9;thuman[12]=9;thuman[90]=9;thuman[82]=9;
thuman[175]=10;thuman[166]=10;thuman[104]=10;thuman[78]=10;
thuman[103]=10;thuman[114]=10;thuman[39]=11;thuman[68]=11;
thuman[108]=11;thuman[8]=11;thuman[131]=11;thuman[143]=11;
thuman[64]=12;thuman[105]=12;thuman[1]=12;thuman[10]=12;thuman[134]=12;
thuman[155]=12;thuman[178]=13;thuman[181]=13;thuman[189]=13;
thuman[168]=13;thuman[107]=13;thuman[197]=13;thuman[52]=14;
thuman[93]=14;thuman[163]=14;thuman[173]=14;thuman[97]=15;
thuman[179]=15;thuman[9]=15;thuman[101]=15;thuman[141]=16;
thuman[41]=16;thuman[13]=16;thuman[74]=16;thuman[19]=17;thuman[169]=17; thuman[20]=17;thuman[14]=17;thuman[18]=18;thuman[75]=18;thuman[117]=18; thuman[172]=18;thuman[6]=19;thuman[27]=19;thuman[7]=19;thuman[184]=19;
thuman[182]=20;thuman[165]=20;thuman[70]=20;thuman[40]=20;
thuman[180]=21;thuman[29]=21;thuman[132]=21;thuman[150]=21;
thuman[192]=22;thuman[112]=22;thuman[157]=22;thuman[149]=22;
thuman[162]=23;thuman[125]=23;thuman[138]=23;thuman[34]=23;
thuman[25]=24;thuman[137]=24;thuman[145]=24;thuman[43]=24;
thuman[57]=25;thuman[4]=25;thuman[146]=25;thuman[37]=25;thuman[69]=26;
thuman[156]=26;thuman[148]=26;thuman[142]=26;thuman[133]=27;
thuman[30]=27;thuman[36]=27;thuman[144]=27;thuman[152]=28;
thuman[154]=28;thuman[111]=29;thuman[127]=29;thuman[118]=30;
thuman[28]=30;thuman[186]=31;thuman[24]=31;thuman[99]=32;
thuman[121]=32;thuman[84]=33;thuman[22]=33;thuman[38]=34;thuman[58]=34; thuman[79]=35;thuman[60]=35;thuman[126]=36;thuman[3]=36;thuman[135]=37; thuman[15]=37;thuman[140]=38;thuman[65]=38;thuman[153]=39;
thuman[63]=39;thuman[159]=40;thuman[32]=40;thuman[190]=41;
thuman[21]=41;thuman[196]=42;thuman[185]=43;thuman[195]=44;
thuman[45]=45;thuman[193]=46;thuman[88]=47;thuman[67]=48;thuman[16]=49; thuman[48]=50;thuman[49]=51;thuman[47]=52;thuman[115]=53;
thuman[123]=54;thuman[164]=55;thuman[188]=57;thuman[17]=59;
thuman[31]=62;thuman[33]=66;thuman[42]=70;thuman[53]=74;thuman[62]=74;
thuman[94]=80;thuman[113]=85;thuman[136]=87;thuman[151]=90;
thuman[116]=92;thuman[87]=96;thuman[85]=99;thuman[130]=100;
thuman[129]=103;thuman[119]=104;thuman[83]=108;
/* Epidemiologically relevant initialisations: */
lambda= 0.3;
my = 0.25;
theta = 1.8;
kappa= 0.025;
/* Purely technical initialisations: */
tts=1000000;
tr = 0;
dt = 0;
te = 0;
at=0;
h=0;i=0;j=0;k=0;t=0;t1=0;
no=0;
r1 = 0; r2 = 0; r3 = 0; r4 = 0; r5 = 0; r6 = 0;
sum1 = 0;sum2 = 0;
hsts=0; sts = 0;
/* Initialisations for the output: */
for (i=0; i<=msm; ++i) {d[i]=0;}
/*-----*/
/*-Main-*/
/*-----*/
/*-Calculation of elapsed time-*/
/* begin of while loop of the entire Main program*/
while (at<tts) {
/* When does the next event occur? */
r1 = ((((double)random())+1)/2147483650.0);
/* Calculation of sum1: {\sum_{i=0}^{i=maxsl} my*i*nschisto[i]} */
sum1=0;
for (i=1; i<=maxsl; ++i) {
sum1=(sum1+(my*i*nschisto[i]));
}
/*
Calculation of sum2: {\lambda*nschisto[0]*(1-(nschisto[0]/M))} */
sum2=0;
sum2=(lambda*nschisto[0]*(1-(nschisto[0]/M)));
/* Calculation of the total rate */
tr = (sum1 + sum2 + (M*kappa));
/* Calculation of elapsed time and new actual time */
dt = (-log(r1)/tr);
at = (at + dt); /* new actual time */
/* Adjustment of time vector */
for (i=0; i<M; ++i) {
thuman[i]=(thuman[i]+dt);
}
/* Calculation of maxlife */
for (i=0; i<M; ++i) {
if (thuman[i]>maxlife) {
maxlife=thuman[i];
}
}
/* Calculation of the output, small inacurracy: must be that dt << 1 */
for(i=0;i<=msm;++i) {
d[i]=(d[i]+(nschisto[i]*dt));
}
/*-Calculation of the type of event-*/
r2 = ((((double)random())+1)/2147483650.0);
if (r2< (sum1/tr))
te=0; /* a parasite dies */
else if ((r2 >= (sum1/tr)) && (r2 < ((sum1+sum2)/tr)))
te=1; /* there is a new infection */
else
te=2; /* an individual dies */
/*-Readjusting of the states-*/
/* first case: a parasite has died (te=0) */
if (te==0) {
/* Programming strategy: put all individuals in a row and weigh them with the number of parasites they carry */
tnp=0;
for (i=1; i<=maxsl; ++i) {
tnp=(tnp+(i*nschisto[i]));
}
r3=(((((double)random())+1)/2147483650.0)*tnp);
sts=0;
for (i=0; i< M; ++i) {
if((sts <= r3) && (r3 < (sts+shuman[i]))) {
/* parasite in individual i died */
nschisto[shuman[i]]=(nschisto[shuman[i]]-1); nschisto[(shuman[i]-1)]=(nschisto[(shuman[i]- 1)]+1); sts=(sts+shuman[i]);
shuman[i]=(shuman[i]-1);
}
sts=(sts+shuman[i]);
}
}
/* second case: a new infection occurred (te=1) */
else if (te==1) {
/* which individual gets infected? */
sts=0;
r4 = (((((double)random())+1)/2147483650.0)*nschisto[0]);
no=0;
while(sts<=r4) {
if(shuman[no]==0) {
++sts;
}
++no;
}
no=(no-1); /* individual no gets infected */
/* with how many parasites will no be infected? */
sts=0;
for (i=0;i<=msm;++i) {
hsts=0;
if(i==0) {
t1=1;
}
else {
t1=t1*i; /* t1 is the faculty of i (=i!)*/
}
for (j=1;j<=maxsl;++j) {
h=j*theta;
t=i;
hsts=hsts+nschisto[j]*(exp((-h)))*pow(h,t)/t1;
}
sts=sts+hsts;
} /* sts is up to a constant the total force of infection */
r5 = (((((double)random())+1)/2147483650.0)*sts);
sts=0; i=0;
while(r5>sts) {
hsts=0;
if(i==0) {
t1=1;
}
else {
t1=t1*i; /* t1 is the faculty of i (=i!) */
}
for (j=1;j<=maxsl;++j) {
h=j*theta;
t=i;
hsts=hsts+nschisto[j]*(exp((-h)))*pow(h,t)/t1;
}
sts=sts+hsts;
++i;
}
k=i-1; /* individual no is infected with k parasites */
nschisto[k]=(nschisto[k]+1);
nschisto[0]=(nschisto[0]-1);
shuman[no]=k;
if (k>maxsl) {
maxsl=k;
}
} /* end of case else if (te==1) */
/* third case: an individual dies (te=2) */
else if (te==2) {
/* Programming strategy: put all individuals in a row */
sts=0;
r6 = (((((double)random())+1)/2147483650.0)*M);
for(i=0;i<M;++i) {
if((sts<=r6) && (r6<(sts+1))) {
/* individual i dies */
nschisto[0]=(nschisto[0]+1);
/* may be wrong, see next line */
nschisto[shuman[i]]=(nschisto[shuman[i]]-1);
/* incl. correction if uninfected dies */
thuman[i]=0;
shuman[i]=0;
}
sts=(sts+1);
}
}
} /* This ends the loop of the main program */
/*-----*/
/*-End-*/
/*-----*/
/*-Calculation of the output-*/
/* Calculation of the ratios */
for(i=0;i<=msm;++i) {
d[i]=d[i]/(tts*M);
}
/*-Output-*/
printf("The number of uninfected individuals at the end is (from a total of 200)");
printf("\n");
printf(" %8.8f",
nschisto[0]);
printf("\n");
printf("The maximal number of worms in an individual was");
printf("\n");
printf(" %8.8f", maxsl);
printf("\n");
printf("The distribution of parasites per individual is:");
printf("\n");
for (i=0; i<=msm; ++i) {
printf("The proportion of individuals with");
printf(" %3d", i);
printf("parasites is ");
printf("%8.8f", d[i]);
printf("\n");
}
}