/* SERVICE ROUTINES
Collected here are general service routines common to different versions of
the main program. These routines are attached with a "#include" statement
rather than being compiled as a separate module, which allows some of them
access to data structures in the main program.
/*----------------------------------------------------------------------------*
CONTAGION KERNEL
Two contagion kernels are in this version. The method is to generate a
probability, then run it backwards, through the inverse cumulative contagion
kernel, to obtain a spatial displacement from the source of the infection. The
individual at that displacement is the target.
This program uses continuous time but not continuous space. It is not yet clear
how to make continuous space efficient, but a regular spatial lattice is
extremely efficient. This version makes space into a torus.
For actual applications, arbitrary two-dimensional kernels of any complexity
can be added without loss of speed, using the cumulative inverse method
demonstrated here. That requires a direct-access, high-speed table lookup, to
be added in actual versions.
ENTRY: 'n' indexes the individual about to infect another.
'kernel' selects the contagion kernel.
0 Mean field - all individuals in the population are equally accessible.
1 Cauchy kernel - nearby individuals are more accessible.
2 Gaussian kernel - similar to above but tighter (not implemented).
'indiv' contains the number of individuals in the simulation.
'trho' and 'nrho' contain the total dispersal distance and the number
of dispersal events, for later calculation of the mean dispersal
distance.
'tinfections' and 'linfections' contain the total number of infections
targetted and the number that fell within the geographic area.
EXIT: 'Kernel' contains the number of the individual targeted for infection.
'trho', 'nrho', 'tinfections', and 'linfections' are updated as
necessary.
*/
#define PI 3.141592653589793238462643
//dec trho, nrho, tinfections, linfections; //Statistics for local dispersal.
/*
int Kernel(int n)
{ int i, j, k, di, dj;
dec dx, dy, rho, theta;
switch((int)kernel)
{ case 0: //MEAN-FIELD KERNEL.
tinfections += 1; linfections += 1; //Advance the number of infections.
do k = 1+Rand()*indiv; while(k==n); //Select equally from all
return k; //individuals in the system.
//-tinfections could be incremented here? Which would mean infections attempted
//-rather than accomplished.
case 1: //CAUCHY KERNEL.
while(1)
{ tinfections += 1; //Advance the number of infections.
rho = abs(Cauchy(0., sigma)); //Generate a Cauchy radius and
trho += rho; nrho += 1; //maintain dispersal statistics.
theta = 2*PI*Rand(); //Generate a random angle and
dx = rho*cos(theta); di = round(dx); //convert angle and radius to a
dy = rho*sin(theta); dj = round(dy); //lattice displacement.
i = (di+A[n].i); //Add the displacement to the
j = (dj+A[n].j); //location of this individual and
if(i<0 || i>=gi) continue; //make sure it is not out of the
if(j<0 || j>=gj) continue; //area.
linfections += 1; //Advance the number of infections
k = G[i][j]; if(k!=n) break; //that fall within the geographic
} //area and return the individual
return k; //occupying that place.
default: Error1(714., "",kernel); //(Invalid kernel number.)
}
}
*/
/*----------------------------------------------------------------------------*
UNIFORM DISTRIBUTION
This routine generates random numbers within a specified interval from a
uniform distribution --- that is, were all values are equally likely.
ENTRY: 'a' defines the lower bound of a uniform distribution.
'b' defines the upper bound.
EXIT: 'Uniform' contains a random number uniformly distributed in the
range 'a' to 'b'.
*/
dec Uniform(dec a, dec b)
{
return Rand()*(b-a)+a;
}
/*----------------------------------------------------------------------------*
EXPONENTIAL DISTRIBUTION
This function generates random time intervals between independent events such
that the events obey a Poisson distribution with mean and variance 'lambda*dt'
in any time unit 'dt'. The resulting intervals obey an exponential
distribution with mean '1/lambda' and variance '1/lambda**2'.
ENTRY: 'lambda' contains the average number of events per unit time, which
must be greater than zero.
'LIMITG' defines a multiple of the average time interval. No time step
will be greater than 'LIMITG/lambda'.
EXIT: 'Expon' contains the next Poisson time interval.
*/
#define LIMITG 10
dec Expon(dec lambda)
{ dec r, expdt;
while(1) // Generate a uniformly distributed
{ r = Rand(); if(r==0) continue; // random number.
expdt = -log(r); // Convert it to an exponential
if(expdt>LIMITG) continue; // distribution and control its
if(expdt==0) continue; // range.
return(expdt/lambda); // Scale and return the next time
} // interval.
}
/* NOTE: This routine must make practical concessions to the reality of
fixed-point arithmetic. First, the random number generator can return a value
that is precisely equal to zero. In infinite-precision arithmetic, that would
occur with probability zero, but here it occurs with finite probability, so it
must be ignored to avoid singularities that would otherwise halt execution.
Second, and related, the finite granularity of the random number generator
near zero seems able to generate excessively large time intervals, so those
are ignored also. The entry parameter 'LIMITG' helps accomplish that. Finally,
intervals of zero are prevented to avoid violating the principle that only one
event can occur at a time. An interval of zero would only happen if 'Rand()'
returned a value of unity. It should not, but this routine provides an
additional safeguard against floating-point rounding anomalies.
*/
/*----------------------------------------------------------------------------*
NORMAL DISTRIBUTION
This routine generates random numbers from the Gaussian distribution,
f(x) = (1/sqrt(2*PI))*exp(-x*x/2), then adjusts them to a specified mean and
standard deviation.
\ENTRY: 'mu' contains the mean of the distribution.
'sigma' contains the standard deviation.
The random sequence of 'Rand' is ready to use.
\EXIT: 'Gauss' contains a Gaussian-distributed random number of the
specified mean and standard deviation.
*/
dec Gauss(dec mu, dec sigma)
{ dec v1, v2, w;
do { v1 = 2.*Rand() - 1.; //Generate a point within the
v2 = 2.*Rand() - 1.; //unit circle.
w = v1*v1 + v2*v2; }
while(w>1 || w==0);
w = v2 * sqrt(-2.*log(w)/w); //Return the normal deviate.
return mu + w*sigma;
}
/* NOTE: The routine uses the polar method due to G.E.P.Box, M.E.Muller, and
G.Marsaglia, described in D.E.Knuth, The Art of Computer Programming, Volume
2 (Addison Wesley, 1969). This method involves, roughly, picking a point
uniformly distributed within the unit circle and then projecting it.
NOTE: It would be possible to save random number 'v2' in a static variable and
use it as 'v1' on the next call. That would make the routine slightly faster at
the cost of a little complexity.
*/
/*----------------------------------------------------------------------------*
LOGNORMAL DISTRIBUTION
This routine generates random numbers from the lognormal distribution, given
the mean and standard deviation of the underlying normal distribution.
\ENTRY: 'mu' contains the mean of the underlying normal distribution.
'sigma' contains the standard deviation.
The random sequence of 'Rand' is ready to use.
\EXIT: 'LogNormal' contains a lognormally-distributed random number of the
specified parameters.
*/
dec LogNormal(dec mu, dec sigma)
{
return exp(mu + sigma*Gauss(0., 1.));
}
/*----------------------------------------------------------------------------*
CAUCHY DISTRIBUTION
This routine generates random numbers from the Cauchy distribution,
f(x) = 1/(1+x*x)/PI, adjusted to specified median and width parameters.
\ENTRY: 'mu' contains the median of the distribution. Note that this is not
the mean, since the Cauchy distribution has no mean.
'sigma' contains the width of the distribution, at which the height
falls to half its maximum value. Note that this is not a standard
deviation, since the Cauchy distribution has no standard deviation.
The random sequence of 'Rand' is ready to use.
\EXIT: 'Cauchy' contains a Cauchy-distributed random number of the
specified median and width.
NOTE: Values are generated from the inverse cumulative Cauchy distribution.
I am not sure I have the width parameter is precisely right. To be checked out
further.
*/
dec Cauchy(dec mu, dec sigma)
{
return mu + sigma*tan(PI*(Rand()-0.5));
}
/*----------------------------------------------------------------------------*
PARAMETERS.
Parameters have default values that can be overridden by command-line arguments.
Each entry on the command line begins with the parameter name, followed by an
equal sign (=), followed by the parameter value. For example, to set a
parameter named 'ralpha0' to 0.21 and another named 'beta1' to 0.718, and also
to clear three parameters named 'mu0', 'mu1', and 'mu2', the following
command-line string could be used:
'ralpha0=.21 beta1=.718 mu0=mu1=mu2=0'
ENTRY: 'argc' contains the number of arguments plus one.
'argv' points to a list of arguments in its second and higher positions.
'pntab' contains a list of valid parameter names.
'pstab' contains pointers to the corresponding parameter values.
EXIT: All parameters correctly specified by 'argv' on entry are set.
Error messages are issued for any others.
*/
#define STR 100
char msg001[] = "Parameter: %s=%s\n";
char msg101[] = "E101. Parameter %d (%s) does not have the correct format"
" (name=value).\n";
char msg102[] = "E102. Parameter %d does not contain a simple numeric"
" value (contains \"%s\").\n";
char msg103[] = "E103. Parameter %d (%s) is not a recognized name.\n";
gparam(int argc, char *argv[])
{ static int lines;
int i, j, k, m, n;
char str[STR+1], *cval;
dec val, atof();
for(i=1; i=0; k--) //Scan for the beginning of the
if(argv[i][k]=='=') break; //parameter value and issue an
if(k<0) //error message if no equal signs
{ printf(msg101, i, argv[i]); //are to be found.
continue; }
k += 1;
for(m=n=0, j=k; argv[i][j]; j++) //Make sure the parameter value is
{ if(argv[i][j]=='-' && j==k) continue; //a simple decimal number.
if(argv[i][j]>='0' && argv[i][j]<='9') { m += 1; continue; }
if(argv[i][j]=='.') { n += 1; continue; }
break; }
if(argv[i][j] || m<1 || n>1) //If it is not, issue an error
{ printf(msg102, i, &argv[i][k]); //message and skip processing it.
continue; }
cval = &argv[i][k]; //Otherwise convert it to internal
val = atof(cval); //form.
for(k=0; ; k+=j+1) //Begin looping through all names.
{
for(j=0; j=0; i++)
{ w = tab[subset[i]]; if(x>w) { x = w; m = i; } }
return subset[m];
}