Step 10. Make things move spatially
To make cells talk to each other is a little more complex than what we have been doing before. Whereas so far we were simply using some predefined commands, and the model we built in Stella, from now on if we are to define some meaningful spatial interaction we will need to do some programming.
There are some modules that we can use in the Library of Hydro-Ecological Modules, however there are not too many things we can do with those predesigned modules. If you really want to be able to build complex spatial models you probably will need to be capable of some level of c++ programming. SME supports so called User Code and offers full access to its classes and methods, which can significantly help you design your own code for spatial dynamics.
Suppose for our Rabbits & Grass model we wish to allow rabbits to move between cells in search of better grazing conditions. We will assume that whenever rabbits find that there is more grass in the neighbor cell a certain proportion of rabbits from the current cell will move to the cell with more grass. Let us write the code that will describe this behavior of the predator.
/**************************************************************************/
#include "Rabbit.h"
/**************************************************************************/
void MoveRabbits( CVariable& Rabbits, CVariable& Grass, CVariable& Rate )
// moves rabbits toward more grass, if there are less rabbits there
// arguments come from MML.config file, first arg is always variable being configured.
{
Grid_Direction il;
float fr, R_moved = 0.;
DistributedGrid& grid = Rabbits.Grid();
grid.SetPointOrdering(0);
// sets grid ordering to default ordering (row-col) (ordering #0)
Rabbits.LinkEdges();
Grass.LinkEdges();
static CVariable* R_Flux = NULL;
if(R_Flux == NULL )
R_Flux = Grass.GetSimilarVariable("R_Flux");
// intermediate increment to Rabbits
R_Flux->Set(0.0);
for( Pix p = grid.first(); p; grid.next(p) )
{
const OrderedPoint& pt = grid.GetPoint(p);
// sets currentPoint
if( !grid.onGrid(pt) ) continue;
// (onGrid == False) -> Ghost Point
float g_max = Grass(pt);
Pix p_max = p;
// for each point calculate where is the max Grass
// in the vicinity
for( il = firstGD(); moreGD(il); incrGD(il) )
{
// enum Grid_Direction {NE=2, EE, SE, SS, SW, WW, NW, NN};
Pix rp = grid.NeighborPix( p, il );
// relative to pt, takes enum Grid_Direction as arg
if( rp )
{
const OrderedPoint& rpt = grid.GetPoint(rp);
if ( Grass(rpt) > g_max )
{ g_max = Grass(rpt);
p_max = rp;
}
}
}
const OrderedPoint& pt_max = grid.GetPoint(p_max);
// sets currentPoint
if ( g_max > Grass(pt) )
fr = ( Rabbits(pt) > Rabbits(pt_max) ) ?
(Rabbits(pt) - Rabbits(pt_max)) * Rate(pt) : 0;
(*R_Flux)(pt_max) += fr;
(*R_Flux)(pt) -= fr;
R_moved += fr;
}// end area loop
Rabbits.AddData(*R_Flux);
printf ("\ninfo: Rabbits moved = %f", R_moved);
}
/**************************************************************************/
So here we have only Rabbits moving horizontally from one cell to another in search of a better life. How do we tell SME that there is something new that the model wants to take into account?
First we go all the way back to the MML.config file that we can find in the Config directory. In this file we add a command for Rabbits:
* RABBITS UF( Rabbit,MoveRabbits,GRASS,RATE)
Here Rabbit is the name of the file that contains the above c++ code. Actually its name is Rabbit.cc and it lives in the UserCode directory. MoveRabbits is the name of the function in this file that we use. GRASS and RATE are two variables that are passed to this function. While GRASS was always there, RATE is actually new. The way I got into this config file is by modifying the Stella model and adding another variable into there. Then once again I had to export the equations, and then I had to do the SME import command. Alternatively I could have modified the equation file at the very beginning of this manual. I simply needed to add one line:
rate = 0.5
And then I could also do this by hand in the R_G1.MML.config file. I would caution you about this however, since it is very easy to forget about some of these small modification to your equation file. There is no way you can import these modifications from the equations to the Stella model. As a result once you finally decide that you wish to modify the Stella model for some other reason later on, most likely you will forget then about this modification of the equations file that you have just done. You will then take the eqations from Stella, create a new equtions file and loose all these previous changes. When you start running the new model you will find out that it performs quite contrary to your expectations, and will take you quite a while to figure out why and to do all the little updates once again. So while every now and then it seems very simple to modify just the equation file, actually you will be much better off if all the modifications are done directly to the Stella model.
As we remember, whenever the equations or the MML.config file is changed we need to do the SME import command. Now we can do the SME build command, and update the config file to add the RATE parameter to it as well. Remember - either you do it by hand, or rename the file to use the R_G1.xxx.conf.out instead. As a result you get:
* RATE pm(0.5)
Are we ready to run? We can try, but I bet that the driver will choke, most likely with a Buss Error. This will leave you quite perplexed and frustrated because it is pretty hard to debug the code at this stage. However, I will give you a hint.
The variables that you have been passing to the newly designed function to move Rabbits were all assumed to be spatial.
MoveRabbits( CVariable& Rabbits, CVariable& Grass, CVariable& Rate )
However the RATE parameter as we defined above is a scalar. There is an easy fix. Just add the override command:
* RATE pm(0.5) oi(0,1)
and you will be back in game. Alternatively you could also define this parameter as a map:
* RATE d(A,${RMAP},${RMAP}) S(.5e+00,0.0)
Here I used the Study Area map to initialize this parmater, which with this scaling factor is identical to what we did above. However, this could be any map, which would probably be the only reasonable way to define this parameter if we wanted it to be spatially heterogeneous.
Or if you do not want this parameter to be spatial, do not refer to it as if it were spatial in the code. Replace Rate(pt) for Rate.Value(). Rate.Value() is a scalar, it will not need to be initialized by a map or a spatial variable. It will take pm(0.5). That's it.
Finally we are ready to hit the SME run button and watch something moving across the landscape, Rabbits hoping from one place to another, Grass dying and regrowing back when the predators leave, and so on. Lots of fun!
E-mail to Alexey Voinov