// Patch Selection Model -- C++ version
// Clark, C.W. & M. Mangel. 1999. Dynamic State Variable Models in Ecology:
//    Methods and Applications. Oxford University Press. /

// **********************************************************************
//
// This is a C++ version of the code for the patch-selection model
// discussed in Chapter 1. (The program is essentially written in C, with only
// minor uses of features of C++.)
//
// This version also includes Monte-Carlo simulation of the reserves state X(t),
// as described in Section 2.3 and illustrated in Fig. 2.3.
//
// **********************************************************************

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

const int  XMAX = 30;				// Capacity
const int  XREP = 4;				// Critical reserves level for reproduction
const int  TMAX = 20;			  	// Time horizon

double 	m[3] = {.01,.05,.02};			// Predation probabilities, per period
double 	p[3] = {.2,.5, 0};				// Discovery probabilities, per period
int 	y[3] = {2,4,0};				// Food sizes
int 	cost = 1;	  				// Metabolic cost, per period
int	c3 = 4;					// Maximum per-period reproduction
double	acap = 60.0;				// Terminal fitness asymptote
double	x0 = 7.5;					// Terminal fitness coeff.
int	xinit = 2;				// Initial reserves, for simulations
int	num_sims = 30;				// No. of simulations

long		addr(int x, int t);
int		chop(int x, int a, int b);
double		V(int m, int x, int t, double* f);
void		Initialize(double* f);
void		Update(double* f, int* dec, int t);
void		Simulate(int* dec, int* state, int* tot_rep, int* fate);
void		Output1(double* f, int* dec, int t);
void		Output2(int* state, int tot_rep, int fate);
void		Output3(int* vital, int sum1, int sum2);

void printarr(double *f) {
  int x,t;
  for ( x=0; x<=XMAX; x++ ) {
    for ( t=0; t<=TMAX; t++ ) {
      printf("%4.1f ", *(f + addr(x,t)));
    }
    printf("\n");
  }
}


//  **************************  Main Program  ***************************

int main()
{
   int	ltime,utime,fate,sum1,sum2;
   int	tot_rep;
   int xx=3;

   double		f[(XMAX + 1)*(TMAX+1)] = {0};
   int		dec[(XMAX + 1)*(TMAX+1)] = {0};
   int		state[TMAX + 1] = {0};
   int		vital[3] ={0};

   FILE *file1;                     //Clear output file
   if ((file1 = fopen("Patchnew_out.txt", "wt")) == NULL)
   {
      printf("%s\n", "Error opening file");
      exit(1);
   }
   fprintf(file1, "");
   fclose(file1);

   Initialize(f);					// Terminal fitness
   for (int t=TMAX-1; t>=1; t--)    			// Solve DPE and output results
  {
	Update(f,dec,t);
        //printf("%g\n", *(f+addr(xx,t)));
   	Output1(f,dec,t);
  }
  //xx=3;
  //printf("%g %g %g\n", V(0,TMAX-1,xx,f), V(1,TMAX-1,xx,f), V(2,TMAX-1,xx,f) );
  
  //printarr(f);
  return 0;

   ltime = time(NULL);              			// Initialize Random Number Generator
   utime = (unsigned int) ltime/2;
   srand(utime);

   sum1 = 0; sum2 = 0;              			// Run simulations
   for (int it=1; it<=num_sims; it++)
   {
      	Simulate(dec,state,&tot_rep,&fate);
      	Output2(state,tot_rep,fate);
      	*(vital + fate) += 1;
      	sum1 += tot_rep; sum2 += tot_rep * tot_rep;
   }
   Output3(vital,sum1,sum2);
}

// ****************************  Functions  *****************************

long addr(int x, int t)					//	Address [x,t]
{    
return (x * (TMAX+1) + t);          			//   	e.g., f[x,t] is *(f+addr(x,t))
}

void Initialize(double f[])           				// Terminal fitness
{
    	 int x;

     	for (x=0; x<=XMAX; x++)
		*(f + addr(x,TMAX)) = acap * x / (x + x0);
}

int chop(int x, int a, int b)					// Chop function (integers)
{
     	if (a <= x && x <= b) return x;
	else if (x < a) return a;
     	else return b;
}

double V(int n,int x,int t,double f[])				 //RHS of DPE, V(n,x,t)
{
	int x1,x2,reproduction;

  	 if (n < 2)   					  // Forages
   	{
		x1 = chop(x + y[n] - cost, 0, XMAX);
		x2 = chop(x - cost, 0, XMAX);
		return (1.0 - m[n]) * (p[n] * *(f + addr(x1,t+1)) + (1.0 - p[n]) * *(f + addr(x2,t+1)));
  	 }
   	else                    	 				// Tries to reproduce
   	{
		if (x < XREP)        	  		// Can't reproduce
      		{
      			reproduction = 0;
         			x1 = chop(x-cost,0,XMAX);
      		}
     	 	else if (x < XREP + c3)  			// Limited reproduction
      		{
        			reproduction = x - XREP;
         			x1 = chop(XREP-cost,0,XMAX);
      		}
     	 	else                      				// Full reproduction
      		{
      			reproduction = c3;
         			x1 = chop(x-cost-c3,0,XMAX);
      		}
     		return reproduction + (1.0 - m[n]) * *(f + addr(x1,t+1));
  	 }
}

void Update(double f[],int dec[],int t)			// Updates fitness function
{
	int	x, nopt, n;
	double	maxrhs, test;

	*(f + addr(0,t)) = 0.0;
	for (x=1; x<=XMAX; x++)
	{
		nopt = 0;
		maxrhs = V(0,x,t,f);
                //printf("%g ", maxrhs);
		for (n=1; n<=2; n++)
		{
			test = V(n,x,t,f);
                        //printf("%g ", test);
			if (test > maxrhs)
			{	nopt = n; maxrhs = test;}
		}
                //printf("\n");
		*(f + addr(x,t)) = maxrhs;
		*(dec + addr(x,t)) = nopt;
	}
}

void Simulate(int dec[], int state[], int *tot_rep, int *fate)
							// Monte Carlo simulation
                        						// fate = 0 if alive, 1 if predated, 2 if starved
{
	int	x,t,sum,i,s;

   	x = xinit; t = 1; sum = 0;
   	*(state + t) = x;
   	*fate = 0;
   	while (t <= TMAX)
   	{
   		i = *(dec + addr(x,t));
      		if (i < 2)																			// Forages
      		{
    			if (rand() < m[i] * RAND_MAX)			// Dies of predation
      			{
      				for (s=t+1; s<=TMAX; s++) *(state+s) = 0;
         				*fate = 1;
        	 			t = TMAX + 1;
         			}
        			else if (rand() < p[i] * RAND_MAX)		// Finds food
         				x = chop(x+y[i]-cost,0,XMAX);
         			else 	                           				// No food
         				x = chop(x-cost,0,XMAX);
         			if (x==0)					// Dies of starvation
         			{
         				for (s=t+1; s<=TMAX; s++) *(state+s) = 0;
            				*fate = 2;
           		 		t = TMAX + 1;
         			}
        			 else
         			{
         				t += 1;
            				*(state + t) = x;
         			}
      		}
      		else                             					// Tries to Reproduce
      		{
      			if (x < XREP)        	  			// Can't reproduce
         			x = chop(x-cost,0,XMAX);
        			else if (x < XREP + c3)   				// Limited reproduction
         			{
         				sum += x - XREP;
           		 		x = chop(XREP-cost,0,XMAX);
         			}
         			else                      					// Full reproduction
         			{
         				sum += c3;
            				x = chop(x-cost-c3,0,XMAX);
        	 		}
      			if (x == 0)                      				// Dies of starvation
      			{
      				for (s=t+1; s<=TMAX; s++) *(state+s) = 0;
         				*fate = 2;
         				t = TMAX + 1;
      			}
   			else if (rand() < m[i] * RAND_MAX)  // Dies of predation after reproducing
      			{
      				for (s=t+1; s<=TMAX; s++) *(state+s) = 0;
         				*fate = 1;
         				t = TMAX + 1;
      			}
      			else
      			{
      				t += 1;
         				*(state + t) = x;
      			}
      		}
   	} // end of while loop
   	*tot_rep = sum;
}

void Output1(double f[], int dec[], int t)    			// Output fitness and decisions
{
	int i,x;

   	FILE *file1;
   	if ((file1 = fopen("Patchnew_out.txt", "at")) == NULL)
   	{
   		printf("%s\n", "Error opening file Patchnew_out");
      		exit(1);
  	 }

	fprintf(file1, "Time period %3d\n", t);

	for (i=1; i<=22; i++) fprintf(file1,"-"); fprintf(file1,"\n");
   	fprintf(file1," x  dec(x,t)  F(x,t)\n");
	for (i=1; i<=22; i++) fprintf(file1,"-"); fprintf(file1,"\n");
	for (x = 1; x <= XMAX; x++)
		fprintf(file1, "%3d %5d   %8.3f\n", x, *(dec + addr(x,t)) + 1,
      		*(f + addr(x,t)));

	for (i=1; i<=22; i++) fprintf(file1,"-"); fprintf(file1,"\n\n");
  	fclose(file1);
}

void Output2(int state[],int tot_rep, int fate)        		//Output simulations
{
   	FILE *file1;
   	if ((file1 = fopen("Patchnew_out.txt", "at")) == NULL)
   	{
   		printf("%s\n", "Error opening file Patchnew_out");
      		exit(1);
   	}

  	fprintf(file1, "\nState trajectory:\n");
   	for (int t=1; t<=TMAX; t++)
   		fprintf(file1, "%3d,",*(state + t));
   	fprintf(file1, "\nTotal Reproduction: %6d; Fate: %3d\n", tot_rep,fate);
	fclose(file1);
}

void Output3(int vital[],int sum1, int sum2)			//Output simulation statistics
{
	double 	mean, SD;

   	FILE *file1;
   	mean = sum1 / num_sims;
   	SD = sqrt(sum2 / num_sims - mean * mean);
   	if ((file1 = fopen("Patchnew_out.txt", "at")) == NULL)
   	{
   		printf("%s\n", "Error opening file Patchnew_out");
      		exit(1);
   	}

   	fprintf(file1, "\nTotal Reproduction: mean = %6.2f; SD = %6.2f\n",
   	mean, SD);
   	double survived = *vital * 100.0 / num_sims;
   	double predated = *(vital+1) * 100.0 / num_sims;
   	double starved = *(vital+2) * 100.0 / num_sims;
   	fprintf(file1, "Survived: %5.1f, Predated: %5.1f, Starved: %5.1f\n",
   	 	survived,predated,starved);
	fclose(file1);
}