I'm vaguely experimenting with a C++ library for doing reasonably efficient discrete time simulation whilst keeping things at a slightly-sub-MATLAB level of readability. (I want it to be reasonably efficient because I want to do lots of simulation runs so the statistics/conclusions drawn from it are reliable.) I'm currently a bit stuck about how to handle systems with "probabilistic control flow" in a SIMD way, so I'm experimenting with a simple system derived from the foxes and rabbits stuff Graham, Staffan, Blake and others have been experimenting with. (This is just to see if the framework would be well-suited for these kind of tasks.) Here's a typical example of what a model-programmer would do. One of the goals is to be able to see how the system behaves as a function of various parameter values, so uppercase like F0 denotes an hyper-grid of system values at various parameters, whilst lowercase like f0 is a "block" designed for optimal low-level machine usage.

void

simluateSystem(int noRuns,int noTimesteps,float rfLo,float rfHi,float ffLo,float ffHi,

float initialFoxes,float rabbitsToSustainFox,float maxSupportableRabbits)

{

SysDescription sysDescription;

Parameter RF(LinSpace(rfLo,rfHi,4));

Parameter FF(LinSpace(ffLo,ffHi,4));

sysDescription.cartProdParameters(2,&RF,&FF);

FullTemporary F0(&sysDescription);

FullTemporary F1(&sysDescription);

FullTemporary R0(&sysDescription);

FullTemporary R1(&sysDescription);

Constant c(rabbitsToSustainFox);

Constant capR(maxSupportableRabbits);

Constant zero(0);

URandTemporary u1;

URandTemporary u2;

int i;

for(i=0;i
F0.initialiseConstant(initialFoxes);

F1.initialiseConstant(0);

R0.initialiseConstant(initialFoxes*rabbitsToSustainFox);

R1.initialiseConstant(0);

do{

LL_SIMD f0(F0),f1(F1),ff(FF),r0(R0),r1(R1),rf(RF);

int j;

for(j=0;j
LL_SIMD f=f0+f1;//total foxes

LL_SIMD r=r0+r1;//total rabbits

LL_SIMD e=clampAbove(c*f,r);//rabbits eaten this year

//young rabbits/foxes surviing become mature rabbits/foxes

r1=clampBelow(r0-clampBelow(e-r1,zero),zero);

f1=clampAbove(e/c,f0);

//now mature rabbits/foxes breed (if not already wiped out)

r0=onlyIf(r1>zero,clampAbove(rf*r1*u1,capR-r1));

f0=onlyIf(f1>zero,ff*f1*u2);

}

F0.store(f0);

F1.store(f1);

R0.store(r0);

R1.store(r1);

}while(sysDescription.stillWork());

//TODO: do incorporate array of final array of values into accumulator

}

//TODO: do some analysis

}

void

simluateSystem(int noRuns,int noTimesteps,float rfLo,float rfHi,float ffLo,float ffHi,

float initialFoxes,float rabbitsToSustainFox,float maxSupportableRabbits)

{

SysDescription sysDescription;

Parameter RF(LinSpace(rfLo,rfHi,4));

Parameter FF(LinSpace(ffLo,ffHi,4));

sysDescription.cartProdParameters(2,&RF,&FF);

FullTemporary F0(&sysDescription);

FullTemporary F1(&sysDescription);

FullTemporary R0(&sysDescription);

FullTemporary R1(&sysDescription);

Constant c(rabbitsToSustainFox);

Constant capR(maxSupportableRabbits);

Constant zero(0);

URandTemporary u1;

URandTemporary u2;

int i;

for(i=0;i

F1.initialiseConstant(0);

R0.initialiseConstant(initialFoxes*rabbitsToSustainFox);

R1.initialiseConstant(0);

do{

LL_SIMD f0(F0),f1(F1),ff(FF),r0(R0),r1(R1),rf(RF);

int j;

for(j=0;j

LL_SIMD r=r0+r1;//total rabbits

LL_SIMD e=clampAbove(c*f,r);//rabbits eaten this year

//young rabbits/foxes surviing become mature rabbits/foxes

r1=clampBelow(r0-clampBelow(e-r1,zero),zero);

f1=clampAbove(e/c,f0);

//now mature rabbits/foxes breed (if not already wiped out)

r0=onlyIf(r1>zero,clampAbove(rf*r1*u1,capR-r1));

f0=onlyIf(f1>zero,ff*f1*u2);

}

F0.store(f0);

F1.store(f1);

R0.store(r0);

R1.store(r1);

}while(sysDescription.stillWork());

//TODO: do incorporate array of final array of values into accumulator

}

//TODO: do some analysis

}