The Azimuth Project
Discrete simulation code tutorial

Discrete simulation code tutorial

Idea

This page is a work in progress giving a very brief introduction to the discrete-time simulation library code in the Azimuth Code Project. The code can be found in these directories within the Azimuth Code Project.

Overview

The library is a C++ library which consists of partly a compiled files but with some important template meta-programming elements in header files 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 strikes a balance between ease of use and very high run-time performance. Models are written in C++ but using operator overloading so it looks like the system’s mathematical equations, hiding many details. With any change to the model “structure”, as opposed to parameter ranges, the C++ code is recompiled. (This takes about a second, which given even preliminary test runs take a minute or two to run isn’t a problem.) This approach acheives very high performance.

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(θ)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 yearf 0f_0
no of 1 year old foxesf 1f_1
half of maximum offspring from pair of foxesf ff_f
no of rabbits a fox needs to eat in a year to survivecc
no of rabbits born this yearr 0r_0
no of 1 year old rabbitsr 1r_1
half of maximum offspring from pair of rabbitsr fr_f
maximum no rabbits vegetation can supportK rK_r

The populations evolve from year to year with some stochastic equations (in discrete time):

(1)e =clampAbove(c(f 0 t+f 1 t),r 0 t+r 1 t) r 1 t+1 =vwhenv2wherev=r 0 tclampAbove(er 1 t,0) f 1 t+1 =vwhenv2wherev=clampAbove(e/c,f 0 t) r 0 t+1 =clampAbove(r 1 t+1r fU,K rr 1 t+1) f 0 t+1 =f 1 t+1f fU\begin{aligned} e &= clampAbove(c (f_0^t + f_1^t),r_0^t + r_1^t)\\ r_1^{t+1} &= v \quad when \quad v \geq 2 \quad where \quad v=r_0^t - clampAbove (e - r_1^t ,0)\\ f_1^{t+1} &= v \quad when \quad v \geq 2 \quad where \quad v=clampAbove (e/c,f_0^t)\\ r_0^{t+1} &= clampAbove(r_1^{t+1} r_f U, K_r - r_1^{t+1})\\ f_0^{t+1} &= f_1^{t+1} f_f U \end{aligned}

where UU denotes a fresh random variable uniformly distributed on [0,1)[0,1) for each occurrence.

The inner loop 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 22 or something like the number of rabbits a fox 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 applied to kU[0,1)+lk U[0,1)+l gives U[l,k+l)U[l,k+l) variables) composed of Constant values.
    • initialised with a (possibly degenerate) affine transformation (eg, mapping applied to kU[0,1)+lk U[0,1)+l gives U[l,k+l)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 particular 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