# Discrete simulation code tutorial

## Idea

This page is a work in progress giving a very brief introduction to the discrete simulation library code in the azimuth code project. The code can be found in here within the Azimuth Code Project.

## Overview

The library is a C++ library which is partly a compiled files but with some important elements that need to be included in your code to allow full compiler optimisation. It was written in order to be able to run enough simulations of non-analytic stochastic systems to be able to answer questions like the probability of a predator species surviving. Thus it balances between base of use and run time speed. Models are written in C++ but using operator overloading so it looks like mathematical equations, hiding many details. With any change to the model “structure”, as opposed to parameter ranges, the code is recompiled. (This takes a second or two, which given even preliminary test runs take a minute or two isn’t a problem.)

Typical usage can be seen by looking at the source code for the simple model in testModel0.cxx and the associated makefile. It can be used to produce the results behind the visualisations on Experiments in discrete stochastic simulation using simplified predator-prey.

The idea is to consider a stochastic system $S(\theta)$ which is depends on some parameters $\theta$, where the interest is in understanding the kinds of behaviour the system has for some set of parameters. Consider the following system:

quantitysymbol
no of foxes born this year$f_0$
no of 1 year old foxes$f_1$
half of maximum offspring from pair of foxes$f_f$
no of rabbits a fog needs to eat in a year to survive$c$
no of rabbits born this year$r_0$
no of 1 year old rabbits$r_1$
half of maximum offspring from pair of rabbits$r_f$
maximum no rabbits vegetation can support$K_r$

Next years values are primed versions of existing years values:

(1)\begin{aligned} e &= clampAbove(c (f_0+f_1),r_0+r_1)\\ r_1' &= t \quad when \quad t \geq 2 \quad where \quad t=r_0 - clampAbove (e - r_1 ,0)\\ f_1' &= t \quad when \quad t \geq 2 \quad where \quad t=clampAbove (e/c,f_0)\\ r_0' &= clampAbove(r_f r_1' U, K_r - r_1')\\ f_0' &= f_f f_1' U \end{aligned}

The inner goop of an implementation is

    SysDescription sysDescription;
Parameter FF(ffLo,ffHi,paramStepsF);
Parameter RF(rfLo,rfHi,paramStepsR);
sysDescription.cartProdParameters(2,&FF,&RF);
FullTemporary F0(sysDescription,initialFoxes);
FullTemporary F1(sysDescription,0);
FullTemporary R0(sysDescription,2*initialFoxes*rabbitsToSustainFox);
FullTemporary R1(sysDescription,0);
Constant c(rabbitsToSustainFox);
Constant capR(maxSupportableRabbits);
Constant zero(0);
Constant two(2);

UniformSIMDStream ustream1(randBufferUnifromU32,noTimesteps);
UniformSIMDStream ustream2(randBufferUnifromU32,noTimesteps);

RandomTemporary rabbitReprodU("rabbitReprodU",sysDescription,&ustream1,RF);
RandomTemporary foxReprodU("foxReprodU",sysDescription,&ustream2,FF);
int run;
for(run=1;run<=noRuns;++run){
sysDescription.initialiseVarsForRun();
do{
LL_SIMD f0=F0,f1=F1,r0=R0,r1=R1;
int step=0;
do{
LL_SIMD e=clampAbove(c*(f0+f1),r0+r1);//rabbits eaten this year. note e<=r0+r1
//young rabbits/foxes that survive become mature rabbits/foxes
r1=r0-clampBelow(e-r1,zero);
f1=clampAbove(e/c,f0);
//ensure once the population drops below 2.0 it definitely gets zeroed
f1=onlyIf(f1>two,f1);
//now any surviving mature rabbits/foxes breed
r0=clampAbove(r1*rabbitReprodU,capR-r1);
r0=onlyIf(r0>two,r0);
f0=f1*foxReprodU;
}while(++step<noTimesteps);
F0.store(f0);
}while(sysDescription.stillWork());
}

Systems have five kinds of components in this implementation:

1. Constant: This is a value which is constant for the entire simulation, either a natural constant such as $2$ or something like the number of rabbits a fog needs to eat in a year.

2. Parameter: This is a value which is to be varied to explore behaviour over that range (in parallel with other Parameters). These are specified as linear ranges giving low extreme, high extreme and number of equidistant sample points (similar to MATLABs linspace). Under the hood the system controller (see later) will arrange to March through the combination of these parameter ranges.

3. FullTemporary: This denotes what is traditionally referred to as a state value which is being stored for later analysis. Under the hood the library will keep an array of these values corresponding to each combination of parameters in the search space. Extraction from a FullTemporary happens automatically, but to store an expression (which typically has type LL_SIMD) it is necessarily to use the store() method.

4. LL_SIMD: This is a low level type which should be used for any intermediate values which aren’t being saved for analysis (which would be a FullTemporary). Writing code using only the other types will work correctly, but will be very slow, so it is worth using

LL_SIMD for any value which isn’t being kept.

1. RandomTemporary: This is a random variable which produces a fresh value each time it is referred to. It can be
• initialised with nothing, in which case it produces a basic variable of its underlying type (eg, a uniform random variable),
• initialised with a (possibly degenerate) affine transformation (eg, mapping to $k U[0,1)+l$ which gives $U[l,k+l)$ variables) composed of Constant values.
• initialised with a (possibly degenerate) affine transformation (eg, mapping to $k U[0,1)+l$ which gives $U[l,k+l)$ variables) composed of Parameter values.

These last two ways of setting up a RandomTemporary are semantically equivalent, but under the hood the library moves the transformation application to the best point for performance.

A collection of various values with these types are all controlled by the SysDescription object. Parameters are “registered” using the cartProdParameters() method, which needs the number of Parameter‘s as the first argument.

for a given run through the system, everything is set up using the initialiseVarsForRun() method, followed by a do-loop where the end condition is

while(sysDescription.stillWork());

It is then useful to set any initial values from FullTemporary values into LL_SIMD values. Then any process code can be written, with store() calls wherever desired.

Another useful class is HistogramU32, which can be used to keep track of the number of times a dssfdaf condition holds (see testModel0.cxx for example usage).

Finally, sets of values can be written out in plain text format readable by, for example, Octave using the writeASCII() method of any array contained in an object, eg, FullTemporary or HistogramU32. (For technical reasons, its necessary to call flipFormat(true) method of any histogram before writing it out).

### Todo

• Expand tutorial.

• Mention bug reporting.

category: software