// oneBranch.C
// program to run MCMC to sample from posterior of two-taxon tree
// assuming the Jukes-Cantor model
//
// Compile using:
// g++ -O4 oneBranch.C -o oneb
//
// To execute, type oneb with possible options.
// 
#include <iostream>
#include <iomanip>
#include <vector>
#include <string>
#include <sstream>
#include <fstream>

using namespace std;

const double LOG_ONE_FOURTH = log(0.25);
const double FOUR_THIRDS = 4.0/3.0;

// *******************************************************
//
// Random number generator modified from R, version 1.8.1
// http://www.r-project.org/
// R-1.8.1/src/nmath/standalone/sunif.c
//

static unsigned int I1=1234, I2=5678;

double runif(void)
{
    I1= 36969*(I1 & 0177777) + (I1>>16);
    I2= 18000*(I2 & 0177777) + (I2>>16);
    return ((I1 << 16)^(I2 & 0177777)) * 2.328306437080797e-10; /* in [0,1) */
}

void printSeeds(const unsigned int I1,const unsigned int I2)
{
  cout << "Seeds: " << I1 << " " << I2 << endl;
}
  
// *******************************************************

double logJC(double t,int n1,int n2)
{
  double a = exp(-FOUR_THIRDS*t);
  return 2*(n1+n2)*LOG_ONE_FOURTH + n1*log(1.0 + 3.0*a) + n2*log(1.0 - a);
}

void usageError(char *name)
{
  cerr << "Usage: " << name << " [-c numCycles] [-l lambda] [-n1 numSame] [-n2 numDifferent]"
       << " [-s1 seedOne] [-s2 seedTwo] [-t edge Length] [-w windlwSize]"
       << endl;
  cerr << " Options:" << endl;
  cerr << " -c:  number of cycles to run MCMC (default is 10000)." << endl;
  cerr << " -l:  lambda (reciprocal of prior mean, default is 1.0)." << endl;
  cerr << " -n1: number of unvaried sites (default is 0)." << endl;
  cerr << " -n2: number of varied sites (default is 0)." << endl;
  cerr << " -s1: first seed for random number generator, any positive integer (default is 1234)." << endl;
  cerr << " -s2: second seed for random number generator, any positive integer (default is 5678)." << endl;
  cerr << " -t:  initial value of edge length (default is 1.0)." << endl;
  cerr << " -w:  window size for update (default is 0.1)." << endl;
  exit(1);
}

void initialize(int argc, char *argv[],int& numCycles,double& lambda,
		int& n1,int& n2,double& t,double& w)
{
  cout.setf(ios::fixed, ios::floatfield);
  cout.setf(ios::showpoint);
  char **p;
  if(argc > 1) {
    for(p=argv+1;p<argv+argc;p++) {
      if(strcmp(*p,"-s1")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> I1))
	  usageError(argv[0]);
	if(I1==0) {
	  cerr << "Error: random seeds must be positive." << endl;
	  usageError(argv[0]);
	}
      }
      else if(strcmp(*p,"-s2")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> I2))
	  usageError(argv[0]);
	if(I2==0) {
	  cerr << "Error: random seeds must be positive." << endl;
	  usageError(argv[0]);
	}
      }
      else if(strcmp(*p,"-c")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> numCycles))
	  usageError(argv[0]);
	if(numCycles<0) {
	  cerr << "Error: number of cycles must be nonnegative." << endl;
	  usageError(argv[0]);
	}
      }
      else if(strcmp(*p,"-n1")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> n1))
	  usageError(argv[0]);
	if(n1<0) {
	  cerr << "Error: number of unvaried sites must be nonnegative." << endl;
	  usageError(argv[0]);
	}
      }
      else if(strcmp(*p,"-n2")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> n2))
	  usageError(argv[0]);
	if(n2<0) {
	  cerr << "Error: number of varied sites must be nonnegative." << endl;
	  usageError(argv[0]);
	}
      }
      else if(strcmp(*p,"-l")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> lambda))
	  usageError(argv[0]);
	if(lambda < 1.0e-12) {
	  cerr << "Error: lambda must be positive." << endl;
	  exit(1);
	}
      }
      else if(strcmp(*p,"-t")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> t))
	  usageError(argv[0]);
	if(t < 1.0e-12) {
	  cerr << "Error: initial edge length must be positive." << endl;
	  exit(1);
	}
      }
      else if(strcmp(*p,"-w")==0) {
	if(++p>=argv+argc)
	  usageError(argv[0]);
	istringstream f(*p);
	if(!(f >> w))
	  usageError(argv[0]);
	if(w < 1.0e-12) {
	  cerr << "Error: window size must be positive." << endl;
	  exit(1);
	}
      }
      else usageError(argv[0]);
    }
  }
}

int main(int argc, char *argv[]) {
  
  int numCycles=10000,n1=0,n2=0;
  double lambda = 1.0;
  double t = 2.0;
  double w = 0.5;
  initialize(argc,argv,numCycles,lambda,n1,n2,t,w);

  double logh = logJC(t,n1,n2) - lambda*t;

  int acceptCount=0;
  for(int cycle=0;cycle<numCycles;cycle++) {
    double tstar = fabs(t + (2*w*runif()-w));
    double loghstar = logJC(tstar,n1,n2) -lambda*tstar;
    if(runif() < exp(loghstar - logh)) {
      t = tstar;
      logh = loghstar;
      acceptCount++;
    }
    cout << setprecision(6) << setw(8) << t << endl;
  }
  cerr << "Acceptance rate is " << acceptCount << "/" << numCycles
       << " = " << setprecision(4) << setw(6) << (double)acceptCount/(double)numCycles << endl;
  return 0;
}