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

/**********************************************************/
/*                                                        */
/* by Michael McCarthy, Australian Research Centre for    */
/* Urban Ecology c/- School of Botany, University of      */
/* Melbourne, VIC 3010, Australia. mamcca@unimelb.edu.au  */
/*                                                        */
/* This simualtes stochastic density dependent population */
/* growth, using a Ricker function:                       */
/*                                                        */
/* N(t+1) = N(t) * lambda * exp[- beta * N(t)]            */
/*                                                        */
/* beta beta controls the strength of the density         */
/* dependence - the slope of the relationship.            */
/*                                                        */
/* Stochasticity enters into the process through annual   */
/* random variation in lambda, which is drawn from a      */
/* lognormal distribution with mean mlambda, and variance */
/* varlambda.                                             */
/*                                                        */
/* Transforming gives:                                    */
/*                                                        */
/* ln[N(t+1) / N(t)] = -beta*N(t)/equil + ln(lambda)      */
/*                                                        */
/* The term ln(lambda) will now be distributed as a       */
/* Gaussian (normal) random number. The term beta can     */
/* now be estimated using standard linear regression,     */
/* with a correction for the bias.                        */
/*                                                        */
/**********************************************************/


double Normal(void);
void regress(double *, double *, int, int, double *, double *, double *);
void biascorr(double *start, int n, double *a, double *b, double *sd);


/* gloabal variables */
int DEBUG=0;  /* print debugging info? */
double heldz; /* used in Normal() */
int hold=0;   /* used in Normal() */


void main(int argc, char *argv[])
{
  double N, Nold = 30, gr;
  double Nmanage, NK;
  double randomr;
  double equil = -999.9, MinPop, EMinPop, MinPopM, EMinPopM, MinPopK, EMinPopK;
  double yy[2000], xx[2000];
  int i, j, nsurveys=10, k, ExtCount, HalfCount, ExtCountM, HalfCountM, ExtCountK, HalfCountK;
  double mu=0.3, var, sd=-999.9, beta=0.05, cv=0.2;
  double a, b, ser;
  char c;
  int NITS=1000;

  while (--argc > 0 && (*++argv)[0] == '-')
  {
    switch(c = *++argv[0])
    {
      case 'c':
        cv = atof(*++argv);
        argc--;
        break;

      case 's':
        sd = atof(*++argv);
        argc--;
        break;

      case 'e':
        equil = atof(*++argv);
        argc--;
        break;

      case 'm':
        mu = atof(*++argv);
        argc--;
        break;

      case 'b':
        beta = atof(*++argv);
        argc--;
        break;

      case 't':
        nsurveys = atoi(*++argv);
        argc--;
        break;

      case 'i':
        NITS = atoi(*++argv);
        argc--;
        break;

      default:
        fprintf(stderr, "Unknown option\n");
    }
  }

  srand48(time(NULL));

  if (sd<0)
    sd = cv*mu;
  var = sd*sd;

  if (equil < 0)
    equil = mu / beta;
  else
    beta = mu / equil;

  a = mu;
  b = -beta;
  ser = sd;

  for (j=0; j<NITS; j++)
  {
    Nold = equil;

    for (i=0; i<20; i++)
    {
      gr = - beta*Nold + mu + sd*Normal();
      Nold = exp(gr) * Nold;
    }

    for (i=0; i<nsurveys; i++)
    {
      gr = - beta*Nold + mu + sd*Normal();
      yy[i] = gr;
      xx[i] = Nold;

      Nold = exp(gr) * Nold;
    }

    if (nsurveys)
    {
      regress(xx, yy, i, 0, &a, &b, &ser);
      biascorr(xx, nsurveys, &a, &b, &ser);
    }

    ExtCount = HalfCount = 0;
    ExtCountM = HalfCountM = 0;
    ExtCountK = HalfCountK = 0;
    EMinPop = 0.0;
    EMinPopM = 0.0;
    EMinPopK = 0.0;

    for (k=0; k<1000; k++)
    {
      Nmanage = Nold = NK = equil;
      MinPop = equil;
      MinPopM = equil;
      MinPopK = equil;

      for (i=0; i<100; i++)
      {
/*
        printf("%f\n", Nold);
*/
        randomr = Normal();

        gr = b*Nold + a + ser*randomr;
        Nold = exp(gr) * Nold;

        gr = b*Nmanage + a + ser*randomr;
        Nmanage = 1.05 * exp(gr) * Nmanage;

        gr = b*NK/1.5 + a + ser*randomr;
        NK = exp(gr) * NK;

        if (Nold < MinPop)
          MinPop = Nold;
        if (Nmanage < MinPopM)
          MinPopM = Nmanage;
        if (NK < MinPopK)
          MinPopK = NK;
      }

      EMinPop = EMinPop*k/(k+1.0) + MinPop/(k+1.0);
      EMinPopM = EMinPopM*k/(k+1.0) + MinPopM/(k+1.0);
      EMinPopK = EMinPopK*k/(k+1.0) + MinPopK/(k+1.0);

      if (MinPop < 1)
        ExtCount++;

      if (MinPop < 0.5*equil)
        HalfCount++;

      if (MinPopM < 1)
        ExtCountM++;

      if (MinPopM < 0.5*equil)
        HalfCountM++;

      if (MinPopK < 1)
        ExtCountK++;

      if (MinPopK < 0.5*equil)
        HalfCountK++;
    }

    printf("a = %f , b = %0.10f , sd = %f , EMinPop = %f , ERisk = %f , HalfRisk = %f , ", a, b, ser, EMinPop, ExtCount/1000.0, HalfCount/1000.0);
    printf("ChEMinPopF = %f , ChERiskF = %f , ChHalfRiskF = %f , ", EMinPopM-EMinPop, (ExtCount-ExtCountM)/1000.0, (HalfCount-HalfCountM)/1000.0);
    printf("ChEMinPopK = %f , ChERiskK = %f , ChHalfRiskK = %f\n", EMinPopK-EMinPop, (ExtCount-ExtCountK)/1000.0, (HalfCount-HalfCountK)/1000.0);

  }
}


/**********************************************/
/*                                            */
/* Normal()                                   */
/*                                            */
/* returns a standard normal deviate          */
/*                                            */
/**********************************************/
double Normal(void)
{
  double v1, v2, w, c;

  if (hold)
  {
    hold = 0;
    return(heldz);
  }
  else
  {
    do
    {
      v1 = (double) 2.0*drand48() - 1.0;
      v2 = (double) 2.0*drand48() - 1.0;
      w = v1*v1 + v2*v2;
    }
    while (w >= 1);

    c = sqrt(-2.0 * log(w) / w);
    heldz = c*v2;
    hold = 1;
    return(c*v1);
  }
}


void getstats(double *y, int n,
         double *mean, double *sd)
{
  int i;
  double sum = 0, sumsq = 0;

  for (i=0; i<n; i++)
  {
    sum += y[i];
    sumsq += y[i]*y[i];
  }

  *mean = sum/n;
  *sd = sqrt((sumsq - sum*sum/n)/(n-1));
}


void regress(double *y1, double *y2, int n, int logy, double *a, double *b, double *sd)
{
  int i;
  double sumy1, sumy2, sumy1sq, sumy2sq, sumy1y2;
  double ssy, ssx, sp;
  double correl, slope, seb, sr, srest;
  double minx, maxx, dummyd, di;


/* log transform the y's? */
  if (logy)
    for (i=0; i<n; i++)
      y2[i] = log(y2[i]);


  sumy1 = sumy1sq = sumy2 = sumy2sq = sumy1y2 = 0;

  for (i=0; i<n; i++)
  {
    sumy1 = sumy1 + y1[i];
    sumy2 = sumy2 + y2[i];
    sumy1sq = sumy1sq + y1[i]*y1[i];
    sumy2sq = sumy2sq + y2[i]*y2[i];
    sumy1y2 = sumy1y2 + y1[i]*y2[i];
  }

  sp = sumy1y2 - (sumy1*sumy2/n);
  ssx = sumy1sq - (sumy1*sumy1/n);
  ssy = sumy2sq - (sumy2*sumy2/n);

  correl = sp / sqrt(ssx*ssy);
  slope = sp / ssx;
  sr = (ssy-sp*sp/ssx) / (n-2);
  seb = sqrt(sr/ssx);

/*
  printf("n = %d\n", n);
  printf("Slope = %f\n", slope);
  printf("S.E. of slope = %f\n", seb);
  printf("residual S.E. = %f\n", sqrt(sr));
  printf("Intercept = %f\n", (sumy2-slope*sumy1) / n);
  printf("Correlation = %f\n", correl);
  printf("r - squared = %f\n", correl*correl);
  printf("std dev. of y = %f\n", sqrt(ssy/(n-1)));
*/

  *a = (sumy2-slope*sumy1) / n;
  *b = slope;
  *sd = sqrt(sr);
}


void biascorr(double *start, int n, double *a, double *b, double *sd)
{
  double goala, goalb, meana, meanb, meanx;
  double esta, estb, estsd;
  double x[1000], y[1000];
  int i, j;
  double prevta=0, prevea=0, prevtb=0, preveb=0, slope;

  goala = *a;
  goalb = *b;
  (*sd) *= sqrt(((double) n) / (n-1.0));

  for (i=0, meanx = 0.0; i<n; i++)
    meanx += start[i];
  meanx /= (double) n;
  meana = goala + goalb*meanx;

/*
  printf("meanx = %f\n", meanx);
*/

  if (*b > 0)
  {
    *a = meana;
    *b = 0;

    return;
  }

  do
  {
    for (i=0, meanb=0.0; i<1000; i++)
    {
      x[0] = start[0];

      for (j=0; j<n; j++)
      {
        y[j] = (*b)*x[j] + (*a) + (*sd)*Normal();
        x[j+1] = exp(y[j])*x[j];
      }

      regress(x, y, n, 0, &esta, &estb, &estsd);
      meanb += estb;
    }

    meanb /= 1000.0;

    slope = (prevtb-(*b)) / (preveb-meanb);

    if (slope > 100)
      slope = 100;

    if (slope < -100)
      slope = -100;

    prevtb = *b;
    preveb = meanb;
    *b += slope * (goalb - meanb);

    if (*b > 0)
      *b = 0.0;

    *a = meana - meanx*(*b);

    if ((*b) == 0.0 && (prevtb==0.0))
      return;

/*
    printf("goala = %f meana = %f;  goalb = %f meanb = %f\n", goala, meana, goalb, meanb);
    printf("goala = %f a = %f;  goalb = %f b = %f\n", goala, *a, goalb, *b);
*/

  } while (fabs((goalb - meanb)/goalb) > 0.01);
}