//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
// reblock.cc                                                      2005-10-05
//  
// In-situ binary reblocking algorithm for error estimation (Flyvbjerg) 
// demonstrated for mock Markov chain Monte Carlo simulation.
// USAGE:
//   m reblock; reblock            (m=make??)
//   x energy.out
//   x -settype xydy error.out     (x=xmgrace)
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

// Really, I should be using templates here and the vector class in the C++ std library. ...
#include <iostream> // for cout
#include <iomanip>  // for setf
#include <fstream>  // for ofstream
using namespace std;
ofstream ofs_error("error.out");

const double pi=3.141592653589793238462643383279502884197169399375105820974944;

class brb {
public:
  int Lmax,N,imax;
  double *a,*b;
  brb (int Lmax, int imax) { // Lmax=entries in buffer, imax=size of each entry
    this->imax = imax;
    this->Lmax = Lmax;
    this->N = 0;
    a = new double[imax*(Lmax+1)];  // even C++ doesn't do dynamic a[i][j]
    b = new double[imax*(Lmax+1)];
    int i;
    for (i=0; i<imax*(Lmax+1); i++) b[i]=0.0;
  }
  void append (double *xx) {
    double *yy;
    int i,l,ivar;  // ivar has nothing to do with i
    ivar = N;

    yy = new double[imax];
    for (i=0; i<imax; i++) yy[i]=xx[i];

    l = 0;
    for (;;) {
      for (i=0; i<imax; i++) b[i + imax*l] += yy[i]*yy[i];
      if ((ivar&1) == 0) {
	for (i=0; i<imax; i++) a[i + imax*l] = yy[i];
	break; // exit "for" loop
      } else {
	for (i=0; i<imax; i++) yy[i] += a[i + imax*l];
	ivar >>= 1;
	l++;
      }
    }
    N++;
  }
  void report (double *xxmean, double *xxerr) {
    // PLEASE ENSURE THAT xxmean AND xxerr ARE PROPERLY ALLOCATED BEFORE CALL!
    int i,l,S,M,Lmax_actual,key; // N is this->N 
    double fractchg;
    double *SoB;
    double *aEoMoBM,*aEoEoMoBM;
    SoB = new double[imax];
    aEoMoBM = new double[imax*(Lmax+1)];
    aEoEoMoBM = new double[imax*(Lmax+1)];

    if (N<=1) {
      for (i=0; i<imax; i++) xxmean[i] = xxerr[i] = -999.0;
      cerr << "Empty brb buffer!\n";
      return;
    }

    //----- Find the level of the largest heap, Lmax_actual.
    for (l=0; l<=Lmax; l++)  if ( ((N & (1<<l))) != 0) Lmax_actual=l;
    
    cout << Lmax_actual << endl;
    //----- Accumulate sum from largest heap to smallest.
    // Compute variances at each level.
    for (i=0; i<imax; i++) SoB[i] = 0.0;

    for (l=Lmax_actual; l>=0; l--) {
      cout << l << endl;

      if ( ((N & (1<<l))) != 0) 
	for (i=0; i<imax; i++) 
	  SoB[i] += a[i + imax*l];
      
      S=(1<<l); M=N/S; if(M<=1) continue; // cannot do stats with only 1 block
      
      for (i=0; i<imax; i++) {
	aEoMoBM[i + imax*l] = 
	  sqrt(
	       (b[i + imax*l] - SoB[i]*SoB[i]/M) 
	       / ((double)S*S*(M-1.0)  * M)
	       );  //dodgy
	
	aEoEoMoBM[i + imax*l] = aEoMoBM[i + imax*l] / sqrt(2.0*(M-1));
      }						   
      ofs_error << setw(4) << l;
      ofs_error << setw(24) << aEoMoBM[0 + imax*l];
      ofs_error << setw(24) << aEoEoMoBM[0 + imax*l];
      ofs_error << endl;
    }



    //----- Now for the AI: Computer guesses true error
    // As the block size is increased,
    // compute the CHANGE in the error as a fraction of the ERROR in the error
    // of the block means.  When this is negligible (less than 0.2, say),
    // it means that the error estimate is unaffected by
    // further increases of the block size; i.e., blocks are bigger than
    // the correlation length.  Then we can return the error ...
    cout << "Number of data values N = " << N << endl;
    for (i=0; i<imax; i++) { //i=key in F90 version
      for (l=0; l<=Lmax_actual-2; l++) { //note this carefully
	fractchg =
	  (aEoMoBM[i + imax*(l+1)] - aEoMoBM[i + imax*l])
	  / aEoEoMoBM[i + imax*l];
	if (fractchg < 0.25) break;
      }
      xxmean[i] = SoB[i]/N;
      xxerr[i] = aEoMoBM[i + imax*l];
      cout << "Optimal block size =" << l << endl;
      cout << "Mean =" << xxmean[i] << endl;
      cout << "Error =" << xxerr[i] << endl;
      if (l==Lmax_actual-2) 
	cerr << "Warning: Inter-block correlation remains!\nLonger MC run may be required.\n";
      
      // do this on hindsight ...
      ofs_error << endl;
      for (double blah=0; blah<=Lmax; blah+=0.2) {
	ofs_error << setw(4) << blah;
	ofs_error << setw(24) << aEoMoBM[0 + imax*l];
	ofs_error << setw(24) << aEoEoMoBM[0 + imax*l];
	ofs_error << endl;
      }
    }
  }
};




inline double drand() {return ((double)random()/(RAND_MAX+1.0));}
inline int irand(int n) {return (int)(drand()*n);}
double drand_normal() { //Box-Mueller
  double x,y,u,t;
  do {
    x=2*drand() - 1;
    y=2*drand() - 1;
    u=x*x + y*y;
  } while (u > 1);
  t=sqrt(-2*log(u)/u);
  return t*x; //also t*y
}
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
// Main program: Example of using the brb class
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
int main() {
  int Lmax,Nmax,imax,l,n;
  double theta,costheta,sintheta,corrlen;
  double *xx,*xxmean,*xxerr;
  ofstream ofs_energy("energy.out");

  srandom(time(NULL)); // seed random number generator
  Lmax = 16; // big and safe
  Nmax = (1<<Lmax);

  imax = 1;
  xx = new double[imax];
  xxmean = new double[imax];
  xxerr = new double[imax];
  brb brb1(Lmax, imax);  

  //%%%%% Mock MC simulation: Fabricate correlated data %%%%%%%
  corrlen=1.0;
  costheta=exp(-1/corrlen);
  sintheta=sqrt(1-costheta*costheta);
  xx[0]=drand_normal()   + 1.0;
  for (n=0; n<Nmax; n++) {
    xx[0] = costheta*xx[0] + sintheta*drand_normal() + 1.0;
    ofs_energy << setw(10) << n;
    ofs_energy << setw(24) << xx[0];
    ofs_energy << endl;
    brb1.append(xx);
  }

  //%%%%% Error analysis of artificial data! %%%%%
  brb1.report(xxmean, xxerr);
  int i;
  for (i=0; i<imax; i++) {
    cout << setw(24) << xxmean[i];
    cout << setw(24) << xxerr[i];
    cout << endl;
  }
  system ("xmgrace energy.out");
  system ("xmgrace -settype xydy error.out");

  return 0;
}