// Short description:  Example subroutines in C for perfect simulation of an autologictic model 
// described in "Spatial-temporal modeling of forest gaps generated by colonization 
// from below- and above-ground bark beetle species" 
// by Zhu, Rasmussen, Moeller, Aukema, and Raffa in JASA 103: 162-177 (2008).  

// File name: perfect_simulation.cpp
// Version: June 10, 2008
// Submitted to: JASA Software Archive
// Written by Jakob G. Rasmussen <jgr@math.aau.dk>
// Documented by Jun Zhu <jzhu@stat.wisc.edu>
// The code here can be freely used and redistributed for non-commercial purposes only.

// HEADER

void GibbsAuto(int NumSteps, vector<double> & c, vector<short int> & d, 
	       int NumSites, vector<short int> & NumNeigh, 
	       short int ** & Neigh, double phi, vector<short int> & y);
void PSAuto(vector<double> & c, vector<short int> & d, int NumSites, 
	    vector<short int> & NumNeigh, short int ** & Neigh, 
	    double phi, vector<short int> & y);

// FUNCTIONS
// Notes: Randoms() generates a sample from uniform(0,1)
// Notes: NumNeigh is a vector of # of neighbors at each site
// Notes: Neigh contains the site ID of neighbors at each site

// Gibbs steps for autologistic model
void GibbsAuto(int NumSteps, vector<double> & c, vector<short int> & d, int NumSites, vector<short int> & NumNeigh, short int ** & Neigh, double phi, vector<short int> & y)
{
  int i, j, step;
  double p;
  int sum;

  for (step=0; step<NumSteps; step++){
    for (i=0; i<NumSites; i++){
      if (d[i]==0){
	sum = 0;
	for (j=0; j<NumNeigh[i]; j++){
	  sum += y[Neigh[i][j]];
	}
	p = exp(c[i] + phi*sum)/(1+exp(c[i] + phi*sum));
	if (Randoms()<p){
	  y[i] = 1;
	}
	else {
	  y[i] = 0;
	}
      }
      else {
	y[i] = 0;
      }
    }
  }
}


// Perfect simulation for autologistic model
// Note: upper and lower chains are incorrectly named in the anti-monotone case
// but the program still works correctly 
void PSAuto(vector<double> & c, vector<short int> & d, int NumSites, vector<short int> & NumNeigh, short int ** & Neigh, double phi, vector<short int> & y)
{
  int step, t, i, j;
  int coalescence;
  int maxstep = 100;
  double seed0 = GetSeed();
  vector<double> seed = vector<double>(maxstep);
  vector<short int> U = vector<short int>(NumSites);
  vector<short int> L = vector<short int>(NumSites);

  for (step=0; step<maxstep; step++){
    //cout << step << endl;

    // Initialise upper and lower chain
    for (i=0; i<NumSites; i++){
      L[i] = 0;
      U[i] = 1-d[i];
    } 

    // Saving seed for when chain revisits this step
    seed[step] = seed0;

    // Upper chain, Gibbs steps using new random numbers
    SetSeed(seed0);
    for (t=(int)pow(2.0,step); t>(int)pow(2.0,step-1); t--){
      GibbsAuto(1, c, d, NumSites, NumNeigh, Neigh, phi, U);
    }

    // Lower chain, Gibbs steps using new random numbers
    SetSeed(seed0);
    for (t=(int)pow(2.0,step); t>(int)pow(2.0,step-1); t--){
      GibbsAuto(1, c, d, NumSites, NumNeigh, Neigh, phi, L);
    }

    // Saving seed - needed when returning to creating new random numbers
    seed0 = GetSeed();

    // Gibbs steps reusing old random numbers
    for (j=step-1; j>=0; j--){
      // Upper chain, reusing old random numbers
      SetSeed(seed[j]);
      for (t=(int)pow(2.0,j); t>(int)pow(2.0,j-1); t--){
	GibbsAuto(1, c, d, NumSites, NumNeigh, Neigh, phi, U);
      }
      // Lower chain, reusing old random numbers
      SetSeed(seed[j]);
      for (t=(int)pow(2.0,j); t>(int)pow(2.0,j-1); t--){
	GibbsAuto(1, c, d, NumSites, NumNeigh, Neigh, phi, L);
      }
    }

    // Checking whether chains have reached coalescence
    coalescence = 1;
    for (i=0; i<NumSites; i++){
      if (U[i] != L[i]) coalescence = 0;
    }
    if (coalescence==1){
      for (i=0; i<NumSites; i++){
	y[i] = U[i];
      }
      SetSeed(seed0);

      return;
    }
  }
  cout << "more than 2^" << maxstep << " steps";
  exit(1);
}

# endif