Keymer Lab in silico contact process
We start by the creation of the 1d array that will represent the 2d lattice: <syntaxhighlight lang="c" line="1" > // the int * syntax indicates that this function will return a pointer int * createlandscape2d(int x_size, int y_size) { // this creates the pointer to the array landscape, calloc allocates the size of array on the memory for the whole programm until it is freed int *landscape = calloc(y_size*x_size,y_size*x_size); // this sets all the values in the array to zero memset(landscape, 0, sizeof(landscape)); printf("Your landscape of x = %d and y = %d was successfully created.\n It is still empty. \n", x_size, y_size); // this returns the pointer to the array return landscape; } </syntaxhighlight>
We then initialize this lattice with a single 'cell' in the middle:
<syntaxhighlight lang="c" line="1" >
int * inilandscapesingle(int x_size, int y_size, int *ptom)
{
// this makes landscape the pointer to the array pointed at by ptom
int *landscape = ptom;
// this line creates a single cell at the middle of the lattice
landscape[((y_size*x_size)/2)+ (x_size /2)] = 1;
printf("Your landscape was successfully initialize with a single cell in the middle of the landscape\n");
return landscape;
}
</syntaxhighlight>
We can also initialise the landscape at random by a certain occupancy rate. This uses a random function based on the Mersenne twister. The implementation used here can be found at this adress
<syntaxhighlight lang="c" line="1" >
int * inilandscaperandom(int x_size, int y_size,int *ptom, double occupancy, long long seed)
{
int *landscape = ptom;
// this intializes the Mersenne twister random function
init_genrand64(seed);
int j = 0;
// this is the loop that goes through the whole lattice
while (j < x_size*y_size)
{
// this line generates a random real between 0 and 1
double r = genrand64_real2();
// Here if this condition is verified the position in the array is filled with a 1
if (r < occupancy)
{
landscape[j] = 1;
j++;
}
else
{
j++;
}
} printf("Your landscape was successfully inintialized randomly with occupancy rate equal to %f \n",occupancy); return landscape; } </syntaxhighlight>
This is the randomwalk through the array to update random sites inside the array depending on their neighbours: <syntaxhighlight lang="c" line="1" > int * randomwalk2(int x_size, int y_size,int *ptom, double survival, long long seed) { int *landscape = ptom; init_genrand64(seed); int sites = 0; int size = y_size * x_size; int randomx = 0; int neighbour = 0; // a loop running for total nuber of sites for (sites = 0; sites < size; sites ++) { // creating a random number that will be the site of interest randomx = (int) floor((genrand64_real2()*size)); // switch is like an if except faster in c so here it returns the value of site randomx switch(landscape[randomx]) { case 0: // choosing one of four neighbours of site randomx at random switch((int) ceil(genrand64_real3()*4)) { //Picking the neighbour on the right case 1: neighbour = landscape[(randomx+1)%x_size + (randomx/x_size)*x_size]; break;
//Picking the neighbour on the left case 2: if (randomx-1 < 0) { neighbour = landscape[x_size-1]; break; } else { neighbour = landscape[(randomx-1)%x_size + (randomx/x_size)*x_size]; break; }
// Picking the southern neighbour case 3: neighbour = landscape[(randomx + x_size)%(x_size*y_size)]; break;
//Picking the northen neighbour case 4: if (randomx-x_size < 0) { neighbour = landscape[randomx+x_size*(y_size-1)]; break; } else { neighbour = landscape[randomx-x_size]; break; }
} if (neighbour ==1) { landscape[randomx] = 1; } break; case 1: // computing if an occupied site should die or not depending on death rate if (genrand64_real2() < survival) { landscape[randomx] = 0; } break; }
} return landscape; } </syntaxhighlight>
We then run the simulation for a certain time t : <syntaxhighlight lang="c" line="1" > int runsimulationbmp(int x_size, int y_size,int *ptom, int timeofsim, double survival, long long seed) { init_genrand64(seed); int i = 0; // creating a string to stock the name of the files we wish to create char name[11] = ""; // runs for timesteps = timeofsim while (i < timeofsim) { // updates the name string with iteration number or time t sprintf(name,"images/simu%d.bmp",(i+200)); // run the randomwalk algorithm on the array ptom = randomwalk2( x_size, y_size, ptom, survival, genrand64_int63()); // creates a bmp image this function is descibed bellow createfrommatrix( x_size, y_size, 0, 255, 0, ptom, name); i++; } printf("Your simulation ran for time %d is finished and bmp files were created\n",timeofsim ); // frees the memory taken up by the array Very Important !!! void free(void * ptom); return 0; } </syntaxhighlight>
This is the function that allows the user to transform the lattice in a bitmap image and save it in a file: <syntaxhighlight lang="c" line="1" > int createfrommatrix(int x_size, int y_size, int red, int green, int blue,int *ptom,char name[]) { int *landscape = ptom; // creating a bmp object and it's space in memory BMP* bmp = BMP_Create( x_size, y_size, 24 ); int j = 0; // loop that runs throughout the whole array fed into this function while (j < x_size*y_size) { if (landscape[j] == 1) { // allows the user to set rgb values to pixels BMP_SetPixelRGB(bmp, j%x_size, j/x_size , red , green , blue ); j++; } else if(landscape[j] == 0) { BMP_SetPixelRGB(bmp, j%x_size, j/x_size , 0, 40 , 0 ); j++; } else { BMP_SetPixelRGB(bmp, j%x_size, j/x_size , 255 , 0 , 0 ); j++; } } // creates the file at the indicated name location BMP_WriteFile( bmp, name ); // frees the memory used by the bmp BMP_Free(bmp); return 0 ;
} </syntaxhighlight>
The code to run this model can be found here. If you are using a debian based system: <syntaxhighlight lang="powershell" > tar zxvf Contact_model.tar.gz </syntaxhighlight> Then: <syntaxhighlight lang="powershell" line="1" > cd Contact_model/ make ./test 200 200 500 0.05 0.2 </syntaxhighlight> The first argument is the width, the second the height, the third the number of time steps you want to run it for, the fourth the occupancy rate, and the last the death rate.
