// SIS - Spatial Invasion Simulator
//
// WS Harpole, v3 rewritten in C 7/1/2007 from R prototype
//	
//	The Spatial Invasion Simulator (SIS) is a spatially explicit, two-species lattice model. To explore alternative restoration and invasion scenarios, the two species are designated as Exotic and Native.
//	Individuals of each species occupy discrete cells on a square lattice. At each time step individuals have some probability of mortality (100% if they are annuals). Unoccupied cells are colonized by either species, depending on their dispersal mode, their abundance, and the positive soil feedback state generated by the previous cell occupant.
//	Dispersal can be either local (nearest 8 neighbors) or global (all individuals in the lattice can send propagules to all cells). Boundaries are absorbing and the lattice is updated synchronously.
//	Positive soil feedback decreases the probability of the other species’ establishment and survival. The strength of the feedback increases with time of occupancy of the species up to the maximum value.
//	The Model assumes the following:
//	- The environment is homogenous (there is no variation in the geography, soil type, etc.) and that species are identical in their environmental requirements.
//	- Net feedback is positive (or zero, i.e., positive feedback effects are greater than or equal to negative feedback effects).
//	- There is no long-term seedbanking – all propagules come from the previous time step.
//
//	Alternative management scenario options include:
//	- Solarizing the soil (setting the initial soil feedback state).
//	- Control of the exotic species (increasing exotic mortality).
//	- Augmenting native seed production (decreasing exotic/native seed ratio).
//	- Including an external source of exotic seeds (one edge of the lattice provides a persistent source of seeds).
//
// 
// run in command-line or from web-form with 13 parameter arguments, including process ID from perl script sis.pl, e.g.:
//		sis sp1dis sp2dis sp1fb sp2fb sp1m sp2m scenario initabund seed edge dim tmax pid
//			example cmd line: ./sis 1 0 30 30 5 5 1 80 1 0 50 100 4976


#include <time.h>
#include <stdio.h>
#include <stdlib.h> 

int main ( int argc, char *argv[] ) 
{


// get process ID to set output directory
	char *pid;			
		pid = argv[13] ; 

// Output Files	
	FILE *outtime ;
	FILE *outinitial ;
	FILE *outfinal ;
	FILE *outlog ;
	// tmp filenames
	char timetxt[80];
		sprintf( timetxt,	"/home/sis/public_html/tmp/sis%s/time.txt", pid);
	char initialtxt[80];
		sprintf( initialtxt,"/home/sis/public_html/tmp/sis%s/initial.txt", pid);
	char finaltxt[80];
		sprintf( finaltxt,	"/home/sis/public_html/tmp/sis%s/final.txt", pid);
	char logtxt[80];
		sprintf( logtxt,	"/home/sis/public_html/tmp/sis%s/log.txt", pid);		

	// Open files
	outtime =		fopen(timetxt, "w");
	outinitial =	fopen(initialtxt, "w");
	outfinal =		fopen(finaltxt, "w");
	outlog =		fopen(logtxt, "w");
	// print headers
	fprintf(outtime, "time spp1 spp2\n");
	fprintf( outinitial, "i,j,x\n");
	fprintf( outfinal, "i,j,x\n");

// local time
	struct tm *timeptr;
	time_t tm;
	tm = time(NULL);
	timeptr = localtime(&tm);
	
// set random number seed based on time
	srand( tm );

// set user defined parameters from command-line arguements
// Note: current error handling only checks for out-of-bounds values - doesn't change flow control
//		species 1=exotic, species 2=native
//			if( argc<14) fprintf( outlog, "ERROR: too few parameters %d\n", argc );
	int sp1dis, sp2dis;		// species dispersal: 0=local, 1=global
		sp1dis = atoi(argv[1]) ;
			if( !( sp1dis == 0 || sp1dis == 1) ) fprintf( outlog, "ERROR: sp1dis must be 0 or 1, value entered=%d\n", sp1dis );
		sp2dis = atoi(argv[2]) ;
			if( !( sp2dis == 0 || sp2dis == 1) ) fprintf( outlog, "ERROR: sp2dis must be 0 or 1, value entered=%d\n", sp2dis );
	int sp1fb, sp2fb;		// species feedback strength: 0=none, 100=overwhelmingly strong
		sp1fb = atoi(argv[3]) ;
			if( ( sp1fb < 0 || sp1fb > 100) ) fprintf( outlog, "ERROR: sp1fb must be >= 0 and <=100, value entered=%d\n", sp1fb );
		sp2fb = atoi(argv[4]) ;
			if( ( sp2fb < 0 || sp2fb > 100) ) fprintf( outlog, "ERROR: sp2fb must be >= 0 and <=100, value entered=%d\n", sp2fb );
	int sp1m, sp2m;			// species percent mortality: 0=immortal, <100=perennial, 100=annual
		sp1m = atoi(argv[5]) ;
			if( ( sp1m < 0 || sp1m > 100) ) fprintf( outlog, "ERROR: sp1m must be >= 0 and <=100, value entered=%d\n", sp1m );
		sp2m = atoi(argv[6]) ;
			if( ( sp2m < 0 || sp2m > 100) ) fprintf( outlog, "ERROR: sp2m must be >= 0 and <=100, value entered=%d\n", sp2m );
	int scenario;			// initial soil feedback state: 0=neutral, 1=exotic (spp1), 2=native (spp2)
		scenario = atoi(argv[7]) ; 
			if( !( scenario == 0 || scenario == 1 || scenario == 2 ) ) fprintf( outlog, "ERROR: scenario must be 0, 1, or 2, value entered=%d\n", scenario );
	int initabund;			// initial percent natives
		initabund = atoi(argv[8]) ; 
			if( ( initabund <= 0 || initabund > 100) ) fprintf( outlog, "ERROR: initabund must be >0 and <=100, value entered=%d\n", initabund );
	float seed;				// ratio of exotic to native seed production: 1= default (values: 1000, 100, 10, 1, 0.1, 0.01, 0.001)
		seed = atof(argv[9]) ; 
			if( ( seed < 0.001 || seed > 1000) ) fprintf( outlog, "ERROR: seed must be between 0.001 and 1000, value entered=%f\n", seed );
	int edge;				// 0=absorbing boundaries, 1=one edge all exotics (spp1)
		edge= atoi(argv[10]) ;
			if( !( edge == 0 || edge == 1) ) fprintf( outlog, "ERROR: edge must be 0 or 1, value entered=%d\n", edge );
	int dim;				// dimension of matrix (not counting the edges)
		dim = atoi(argv[11]) ; 
			if( ( dim <= 0 || dim > 1000) ) fprintf( outlog, "ERROR: dim must be > 0 and <=1000, value entered=%d\n", dim );
	int tmax;				// time steps
		tmax = atoi(argv[12]) ; 
			if( ( tmax <= 0 || tmax > 10000) ) fprintf( outlog, "ERROR: tmax must be > 0 and <=10000, value entered=%d\n", tmax );

// save log file 
	fprintf( outlog, "Spatial Invasion Simulator, parameter settings\n\n" );
	fprintf( outlog, "time, %s\n", asctime(timeptr) );
	fprintf( outlog, "sp1dis, %d\n", sp1dis );
	fprintf( outlog, "sp2dis, %d\n", sp2dis );
	fprintf( outlog, "sp1fb, %d\n", sp1fb );
	fprintf( outlog, "sp2fb, %d\n", sp2fb );
	fprintf( outlog, "sp1m, %d\n", sp1m );
	fprintf( outlog, "sp2m, %d\n", sp2m );
	fprintf( outlog, "scenario, %d\n", scenario );
	fprintf( outlog, "initabund, %d\n", initabund );
	fprintf( outlog, "seed, %f\n", seed );
	fprintf( outlog, "edge, %d\n", edge );
	fprintf( outlog, "dim, %d\n", dim );
	fprintf( outlog, "tmax, %d\n", tmax );
	fprintf( outlog, "ProcessID, %s\n", pid );

// declarations of internal variables
	int i, j, k, t;				// matrix and time indexes
	float s1, s2;				// species 1 and 2 survival probabilities
	float c1, c2;				// species 1 and 2 colonization probabilities
	float p1, p2;				// species 1 and 2 establishment probabilities
	int xt0[dim+2][dim+2];		// time t generation, a 10x10 lattice with 1-cell edge boundaries
	int xt1[dim+2][dim+2];		// time t+1 generation, a 10x10 lattice with 1-cell edge boundaries
	int near[3];				// counts of 8 nearest neighbors+cell[i][j]: empty, spp1, and spp2
	int all[tmax][3];			// counts of all cells: empty, spp1, and spp2
	int state[dim+2][dim+2];	// cell state, indicating effect of plant-soil feedback


	// initialize species matrix
	for (i =1; i< dim+1 ; ++i) {
		for (j=1; j< dim+1 ; ++j) {
		xt1[i][j] = 1 ; // set to species 1
		if ( ((double)rand() / ((double)(RAND_MAX)+(double)(1))) < initabund/100. )
			xt1[i][j] = 2; // set to species 2 with P=initial abundance
			switch (scenario) {
				case 0:  // initially neutral soils
					state[i][j] = 0; 
					break;
				case 1:  // initially exotic soils (species 1 values are NEGATIVE)
					state[i][j] = -sp1fb; 
					break;
				case 2:  // initially native soils
					state[i][j] = sp2fb; 
					break;
			}
		}
	}
		
	// absorbing edges
	for (i =0; i< dim+2 ; ++i) {
		xt1[0][i]	  =0;
		if( edge==1 && i<dim+1 ) xt1[0][i]=1; // set one edge to all exotics (spp1)
		xt1[dim+1][i] =0;
		xt1[i][0]	  =0;
		xt1[i][dim+1] =0;
	}
	
	// save initial matrix
	for (i =0; i<dim+2 ; ++i) {
		for (j=0; j<dim+2 ; ++j){
			fprintf( outinitial, "%d,%d,%d\n", i, j, xt1[i][j]);
		}
	}


	// initialize global abundance vs time array
	for (t =0; t<tmax; ++t){
		all[t][0]=0;
		all[t][1]=0;
		all[t][2]=0;
	}


	// run model for tmax time steps	
	for( t=0; t< tmax; ++t){
		
		// determine global species abundance and set previous time step = to current
		for (i =0; i< dim+2 ; ++i) {
			for (j=0; j< dim+2 ; ++j) {
				// set previous t0 matrix to current t1 matrix 
				xt0[i][j] = xt1[i][j] ;
				// global abundance			
				if (xt0[i][j]==0){ all[t][0]++ ; }
				if (xt0[i][j]==1){ all[t][1]++ ; }
				if (xt0[i][j]==2){ all[t][2]++ ; }
			}
		}

		// update current matrix (mortality then colonization)
		for (i =1; i<dim+1 ; ++i) {
			for (j=1; j<dim+1 ; ++j){
			
				// count 8 nearest neighbors plus current cell using previous time t0 matrix (empty, spp1 or spp2)
				near[0]=0;
				near[1]=0;
				near[2]=0;
				for (k =0; k<9; ++k){
					switch ( xt0[i+(k - 3*(k/3) -1)][j+(k/3-1)] ) { 
						case 0: near[0]++ ; 
						break;
						case 1: near[1]++ ; 
						break;
						case 2: near[2]++ ; 
						break;
					}
				}
				
				// feedback state effect on survival probabilities
				s1= s2= 1. ;
				if( state[i][j]>0 ) // positive state value decreases s1
					s1= (1-abs(state[i][j])/100.) ;  
				if ( state[i][j]<0 ) // negative state value decreases s2
					s2= (1-abs(state[i][j])/100.) ;
					
				// Species survival
				switch( xt1[i][j] ) {
					case 1:      // species 1 occupies cell, m% mortality + additional feedback effect
						if(((double)rand() / ((double)(RAND_MAX)+(double)(1))) < (sp1m/100. + (1-s1)) )
						{
							xt1[i][j]= 0;
						}
						if( state[i][j] > -sp1fb ) // decrement cell state in favor of species 1, up to max state -10
							state[i][j] = state[i][j]-10  ; // species 1 has NEGATIVE values
					break; 
					
					case 2: 	// species 2
						if(((double)rand() / ((double)(RAND_MAX)+(double)(1))) < (sp2m/100. + (1-s2)) )
							xt1[i][j]= 0;
						if( state[i][j] < sp2fb ) // increment cell state in favor of species 1, up to max state 10
							state[i][j] = state[i][j]+10  ;
					break;
				}
				
				if( xt1[i][j]==0){ 	// empty cell: colonization probability proportional to local or global neighborhood
								//		weighted by ratio of exotic to native seed production
						if(sp1dis== 0)  //local
							c1= seed * near[1]/9. ;
						else  //global
							c1=  seed * (1.* all[t][1]) / (all[t][1] + all[t][2]) ; 
						
						if(sp2dis==0)  //local
							c2= near[2]/9. ;
						else //global
							c2= ( 1.* all[t][2]) / (all[t][1] + all[t][2]) ; 
							
						// probabilities of establishment: Psurvival x Pcolonization
						p1= c1*s1/(c1*s1+c2*s2) ;
						p2= 1-p1 ;
						if ( ((double)rand() / ((double)(RAND_MAX)+(double)(1))) < p1) 
							xt1[i][j] = 1; // set to species 1
						else
							xt1[i][j] = 2;
				}
			}
		}
	}	

	// save abundance vs. time
	for(t=0; t<tmax; ++t){
		fprintf(outtime, "%d %d %d\n", t, all[t][1], all[t][2] );
	}

	// save final matrix 
	for (i =0; i<dim+2 ; ++i) {
		for (j=0; j<dim+2 ; ++j){
			fprintf( outfinal, "%d,%d,%d\n", i, j, xt1[i][j] );
		}
	}


// close files
	fclose(outtime);
	fclose(outinitial);
	fclose(outfinal);
	fclose(outlog);
 
return 0;
}