//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
// BINARY REBLOCKING BUFFER.  YLL, based on NDD.  JAVA VERSION 2005/05/07.
// THIS IS A NON-STATIC, OBJECT-ORIENTED VERSION.
//
// To compile: javac brbdemo.java     or "m brbdemo"
// To run:     java brbdemo
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
import java.io.*;
import java.util.Random;


//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
// BINARY REBLOCKING BUFFER.  YLL, based on NDD.  JAVA VERSION 2005/05/07.
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
class brb_buffer {
  int Lmax,N;
  double a[],b[];
  double brb_mean,brb_err; // used for returning values  
  FileOutputStream fout;// DEBUUGGING
  PrintStream pout; //HACK (REMOVE IN REAL VESRIN)
  //========== CONSTRUCTOR: ALLOCATE BUFFER Lmax LEVELS DEEP ==============
  brb_buffer(final int Lmax) {
    this.Lmax = Lmax;
    N = 0;
    a = new double[Lmax+1];
    b = new double[Lmax+1];
    int l;
    for (l=0; l<=Lmax; l++) b[l] = 0.0d;  // no need to init a
    try {
      fout = new FileOutputStream ("brb.out");
      pout = new PrintStream (fout);
    } catch (Exception e) {};
  }
  //========== APPEND DATA xx AND UPDATE MULTIHEAP
  public void append (double xx) {
    int l,ivar;
     double yy;
    ivar = N;
    yy = xx;
    l = 0;
    for (;;) {
      b[l] += yy*yy;
      if ((ivar & 1) == 0) {a[l] = yy; break;}
      else {yy += a[l]; ivar>>=1; l++;}
    }
    N++;
  }
  //========== REPORT MEAN AND STDEV in mean and err
  public void report () {
    int l,S,M,Lmax_actual,key;
    double SoB,VoBM,VoMoBM;
    double[] aEoMoBM,aEoEoMoBM;
    boolean flag_warning;
    aEoMoBM = new double[Lmax+1];  
    aEoEoMoBM = new double[Lmax+1];
    if (N<=1) {brb_mean = brb_err = -999.0d; return;}
    System.out.println ("brb_report with N="+N);
    //----- Find the level of the largest heap, Lmax_actual.
    Lmax_actual=-1; // blah
    for (l=0; l<=Lmax; l++) {
      if ( (N & (1<<l) ) != 0) Lmax_actual = l;
    }
    //----- Accumulate sum from largest heap to smallest.
    // Compute variances at each level.
    SoB = 0.0d;
    for (l=Lmax_actual; l>=0; l--) {
      if ( (N & (1<<l) ) != 0) SoB += a[l];
      S = 1<<l; M=N/S; if (M<=1) continue; // can't do stats with 1 block
      VoBM = (b[l] - SoB*SoB/M) / ((double)S*S * (M-1));  //refactored
      VoMoBM = VoBM / M;
      aEoMoBM[l] = Math.sqrt(VoMoBM);
      aEoEoMoBM[l] = aEoMoBM[l] / Math.sqrt((double)(2*(M-1)));
    }
    for (l=0; l<=Lmax_actual-2; l++) { // careful here
      System.out.println (l + "\t" + aEoMoBM[l] + "\t" + aEoEoMoBM[l]);
      pout.println (l + "\t" + aEoMoBM[l] + "\t" + aEoEoMoBM[l]); //out to file
    }
    //----- 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 ...
    flag_warning = false;
    for (l=0; l<=Lmax_actual-2; l++) { // careful here
      if (  (aEoMoBM[l+1] - aEoMoBM[l]) < 0.25d * aEoEoMoBM[l] ) break;
    }
    brb_mean = SoB/N;
    brb_err = aEoMoBM[l];
    if (l==Lmax_actual-2) flag_warning=true;
    if (!flag_warning) 
      System.out.println ("Optimal block size = "+l);
    else
      System.out.println ("Warning: Inter-block correlation remains! Longer MC run may be required.");
  }
}


  





//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
//             TEST DRIVE
//&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
public class brbdemo { 
  //static Random random;
  public static void main (String[] args) {
    Random random;
    brb_buffer buf;

    int Lmax=32;
    int imax=100000;
    int i;
    double x=0.0d;

    random = new Random();
    buf = new brb_buffer(Lmax);
    for (i=0; i<imax; i++) {
      x = x*0.9d + 0.1d* random.nextDouble();
      buf.append (x);
    }
    buf.report ();
    System.out.println ("Mean = "+buf.brb_mean);
    System.out.println ("Err  = "+buf.brb_err);

    try {
      Runtime.getRuntime().exec("xmgrace -settype xydy brb.out");
    } catch (Exception e) {};

  }
}