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
}