/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  Quadratic Monte Carlo Method                                           //
//                                                                         //
//  B. Militzer                                     Berkeley 02-15-22      //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "Random.h"
#include "ParseCommandLineArguments.h"
#include "MC.h"
#include "Timer.h"
#include "Analysis.h"
#include "Blocker.h"
#include "Autocorrelation.h"
#include "AveragesSmall.h"

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#ifdef HARMONIC

class QMCState {
 public:
  int n;
  Array1 <double> p;  
  Array1 <double> dP; // step sizes
  double y;

  QMCState():n(0){
    Reset();
  }
  QMCState(const int n_){
    Resize(n_);
    Reset();
  }
  void Resize(const int nn) {
    n = nn;
    p.Resize(nn);
    dP.Resize(nn);
  }

  int nDim() const {
    return n;
  }
  bool Valid(const bool printFlag=true) const {
    return true;
  }
  void Reset() {
    y = Evaluate();
  }
  double operator[](const int i) const {
    return p[i];
  }
  double & operator[](const int i) {
    return p[i];
  }
  double Parameter(const int i) const {
    return p[i];
  }
  double & Parameter(const int i) {
    return p[i];
  }
  double ParameterStep(const int i) const {
    return dP[i];
  }
  double & ParameterStep(const int i) {
    return dP[i];
  }

  double Evaluate() {
    return (*this)();
  }
  double operator()() {
    double s2 = 0.0;
    for(int i=0; i<n; i++) {
      s2 += sqr(p[i]);
      //      s2 += pow(abs(p[i]),2.5);
    }
    y = s2;
    return s2;
  }

  void PrintParametersAndStepSizes() const {
    for(int i=0; i<nDim(); i++) {
      cout << " QMC state parameter " << IntToStringMaxNumber(i+1,nDim()) << " = " << Parameter(i) << "    step= " << ParameterStep(i) << endl;
    }
  }


  friend ostream & operator<<(ostream & os, const QMCState & s) {
    os << s.p;
    return os;
  }

};

#endif // HARMONIC

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#ifdef RING

class QMCState { // ring potential
 public:
  int n;
  Array1 <double> p;  
  Array1 <double> dP; // step sizes
  double y;
  double R, m, m2, B;

  QMCState():n(0) {
    InitializeConstants();
    Reset();
  }
  QMCState(const int n_) {
    Resize(n_);
    InitializeConstants();
    Reset();
  }
  QMCState(const int n_, const double R_, const double m_, const double B_):R(R_),m(m_),m2(2.0*m),B(B_) {
    Resize(n_);
    Reset();
  }
  void InitializeConstants() {
    R = 1.0;
    m = 1.0;
    m2= 2.0*m;
    B = 0.01;
  }
  void Resize(const int n_) {
    n = n_;
    p.Resize(n_);
    dP.Resize(n_);
    //    Write3(n,p,dP);
  }

  int nDim() const {
    return n;
  }
  bool Valid(const bool printFlag=true) const {
    return true;
  }
  void Reset() {
    Evaluate();
  }
  double operator[](const int i) const {
    return p[i];
  }
  double & operator[](const int i) {
    return p[i];
  }
  double Parameter(const int i) const {
    return p[i];
  }
  double & Parameter(const int i) {
    return p[i];
  }
  double ParameterStep(const int i) const {
    return dP[i];
  }
  double & ParameterStep(const int i) {
    return dP[i];
  }

  double operator()() { // This is where we calculate the potential, return value and store it in 'y'
    if (n<2) return 0.0; // n=0 can happen during default constructor called with Array1<QMCState> a.Resize(10);
    double rho = sqrt(sqr(p[0])+sqr(p[1])); // rho^2 = x^2 + y^2
    //    y = pow(rho-R,m2);
    //    for(int i=2; i<n; i++) { // all parameters except x and y
    //      y += pow(p[i],m2);
    //    }
    y = pow( m2*(rho-R) , m2);
    for(int i=2; i<n; i++) { // all parameters except x and y
      y += pow( m2*p[i] , m2);
    }

    y -= B*p[0]; 
    return y;
  }
  double Evaluate() {
    return (*this)();
  }

  void PrintParametersAndStepSizes() const {
    for(int i=0; i<nDim(); i++) {
      cout << " QMC state parameter " << IntToStringMaxNumber(i+1,nDim()) << " = " << Parameter(i) << "    step= " << ParameterStep(i) << endl;
    }
  }


  friend ostream & operator<<(ostream & os, const QMCState & s) {
    os << "R=" << s.R << " m= " << s.m << " B= " << s.B << " y= " << s.y << " r= " << s.p;
    return os;
  }

};

#endif // RING

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#ifdef ROSENBROCK

class QMCState { // Rosenbrock potential
 public:
  double c,d;
  double x1,x2;
  double dx1,dx2;
  double y;

  QMCState() {
    SetCAndD(5.0,1.0/20.0);
  }
    
  void SetAAndB(const double A, const double B) {
    double c_ = A/B;   // Good and Weare: c=A/B with A=100 and B=20   
    double d_ = 1.0/B; // Good and Weare: d=1/20.0 and B=20.0
    SetCAndD(c_,d_);
  }

  void SetCAndD(const double c_, const double d_) {
    c = c_;
    d = d_;

    double a = sqrt(d); // should equal 1.0
    //    double b = c/(d*d);
    double f = 1.0/sqrt(d);
    double g = 1.0/d;
    double xx = a;   // should equal 1.0
    double yy = a*a; // should equal 1.0

    x1 = f*xx;
    x2 = g*yy;
    //    x1 = 0.0;
    //    x2 = 0.0;
    dx1 = 0.1*x1;
    dx2 = 0.1*x2;

    Reset();
    //    Write(y); Quit("XX");
  }

  int nDim() const {
    return 2;
  }
  bool Valid(const bool printFlag=true) const {
    return true;
  }
  void Reset() {
    Evaluate();
  }
  double operator[](const int i) const {
    return (i==0) ? x1 : x2;
  }
  double & operator[](const int i) {
    return (i==0) ? x1 : x2;
  }
  double Parameter(const int i) const {
    return (i==0) ? x1 : x2;
  }
  double & Parameter(const int i) {
    return (i==0) ? x1 : x2;
  }
  double ParameterStep(const int i) const {
    return (i==0) ? dx1 : dx2;
  }
  double & ParameterStep(const int i) {
    return (i==0) ? dx1 : dx2;
  }

  double operator()() { // This is where we calculate the potential, return value and store it in 'y'
    y = c*sqr(x2-x1*x1) + d*sqr(1.0-x1);
    return y;
  }
  double Evaluate() {
    return (*this)();
  }

  void PrintParametersAndStepSizes() const {
    for(int i=0; i<nDim(); i++) {
      cout << " QMC state parameter " << IntToStringMaxNumber(i+1,nDim()) << " = " << Parameter(i) << "    step= " << ParameterStep(i) << endl;
    }
  }


  friend ostream & operator<<(ostream & os, const QMCState & s) {
    os << "c=" << s.c << " d= " << s.d << " x1= " << s.x1 << " x2= " << s.x2 << " y= " << s.y;
    return os;
  }

};

#endif // ROSENBROCK

/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////

class QMCInfo {
  public:
  double temp;
  double a;
  double minDeltaYFractionAtStart;
  double maxDeltaYFractionAtStart;
  double minDeltaYAtStart;
  double maxDeltaYAtStart;
  bool flagWalk;
  bool flagAffine;
  bool flagQMC;
  bool flagGaussian;
  int nWalkers;
  int nMovers; // default 0, if >0 used to exactly determine size of subset in 'walk' moves. Was called 'nSubset' before
  double fraction;

  int nRuns;
  int nBlocks;
  int nSteps;
  int blockMin; 
  bool print; // IMPORTANT THIS TURNS OFF THE END-OF-BLOCK OUTPUT, which is what we want for long runs **7344**

  friend ostream & operator<<(ostream & os, const QMCInfo & qmc) {
    os << " QMC method: ";

    if (qmc.flagWalk) os << "walk  ";
    if (qmc.flagAffine) os << "affine";
    if (qmc.flagQMC && !qmc.flagGaussian) os << "qmcLin";
    if (qmc.flagQMC &&  qmc.flagGaussian) os << "qmcGau";

    os << " temp= " << qmc.temp << " a= " << qmc.a;
    os << endl;

    os << " nRuns= " << qmc.nRuns;
    os << " nBlocks= " << qmc.nBlocks;
    os << " blockMin= " << qmc.blockMin;
    os << " nSteps= " << qmc.nSteps;
    os << " nWalkers= " << qmc.nWalkers;
    if (qmc.fraction>0.0) os << " moverFraction= " << qmc.fraction;
    if (qmc.nMovers>0)    os << " nMovers= " << qmc.nMovers;
    os << " moves/block= " << qmc.nSteps*qmc.nWalkers; // assumes we try to move every walker once per step
    os << endl;
    
    os << " minDeltaYFractionAtStart= " << qmc.minDeltaYFractionAtStart;
    os << " maxDeltaYFractionAtStart= " << qmc.maxDeltaYFractionAtStart;
    os << endl;
    os << " minDeltaYAtStart= " << qmc.minDeltaYAtStart;
    os << " minDeltaYAtStart= " << qmc.minDeltaYAtStart;
    os << endl;
    os << " print= " << qmc.print;
    os << endl;

    return os;
  }
};

/////////////////////////////////////////////////////////////////////////////////////////////

template <class S> 
int SetUpSimplex(S & s0, Array1 <S> & s, const QMCInfo & qmc, const bool print) {
  //  const bool print = false;
  //  const bool print = true;
  int nDim     = s0.nDim();
  if (print) cout << "Setting up simplex in " << nDim << " dimensions from single state." << endl;

  if (!s0.Valid()) {
    cout << s0 << endl;
    s0.Valid(true);
    error("Cannot set up simplex starting with an invalid parameter set.");
    //    warning("Should not start MC with an invalid parameter set. May be a benign S2<S1 issue.");
  }

  int nWalkers = nDim+1;
  s.Resize(nWalkers); 
  
  s[0] = s0;       // State class needs copy properly. 
  s[0].Evaluate(); // evaluate s0, this takes time
  int nFunc = 1;
  double y0 = s[0].y;
  if (print) Write(y0);
  if (print) s[0].PrintParametersAndStepSizes();

  bool minDYFlag = (qmc.minDeltaYAtStart>0.0) || (qmc.minDeltaYFractionAtStart>0.0); // true if we want to try to generate an ensemble with similar y values
  double minDY1 = qmc.minDeltaYFractionAtStart*abs(s[0].y); // y0 may be less than 0 but it always provides a scale. 
  double minDY2 = qmc.minDeltaYAtStart;                     
  double minDY  = max(minDY1,minDY2); // take the larger of the two because do the same below

  bool maxDYFlag = (qmc.maxDeltaYAtStart>0.0) || (qmc.maxDeltaYFractionAtStart>0.0); // true if we want to try to generate an ensemble with similar y values
  double maxDY1 = qmc.maxDeltaYFractionAtStart*abs(s[0].y); // y0 may be less than 0 but it always provides a scale. 
  double maxDY2 = qmc.maxDeltaYAtStart;                     
  double maxDY  = max(maxDY1,maxDY2); // take the larger of the two. This solve the issue of y~1e-5 and if maxDeltaYFractionAtStart or maxDeltaYAtStart are negative. If both are <0, we need to use maxYAtStartFlag 

  for(int i=0; i<nDim; i++) { // set up initial tetrahedron of nDim points in addition to original one
    s[i+1] = s0; // here we make copies. Make sure the state copies properly.

    double step         = s0.ParameterStep(i);
    double previousStep = step;
    double previousY    = y0;
    int k=0;
    bool reduceDY = true;
    while (true) {
      k++; // increases counter k
      //  pp(i+1,i) = s[i] + step; // do not use s[i] because it may have changed by the call 's(pp[i+1])'
      s[i+1].Parameter(i) = s0.Parameter(i) + step;
      s[i+1].Evaluate(); 
      nFunc++;
      // y[i+1]    = s(pp[i+1]); // warning, depending how State is written, this call may copy pp[i+1] over into state 's' !!!!!!!!!!!
      double dy = s[i+1].y - y0;
      if (print) { cout << " ::::::"; write6(i+1,k,step,y0,s[i+1].y,dy); cout << " range= " << minDY << " ... " << maxDY << endl; }
      //      cout << endl; // print s[i+1] last because the line is sooo long

      if (dy<0.0) { // Surprise, we found a new minimum. Keep it and make it the new state to compare with. (May be lost otherwise.)
	s[0]   = s[i+1];   // needs to be copied into s[0], not into s0, to be used in simplex method
	s[0].y = s[i+1].y;
	y0     = s[i+1].y; // since we write 'y0' for convenience, overwrite this value also
	step *=  +2.0;
	k--; // reduce so that we do one more step after this one. Do not reset k to 0 because the this routine may never terminate otherwise. 
	continue;
      }

      bool dYTooSmall = minDYFlag && (abs(dy)<minDY);
      bool dYTooLarge = maxDYFlag && (dy     >maxDY); // large negative dY are ok
      if (k==1) reduceDY = dYTooLarge; // decide once (in the first step) whether we need to increase or decrease the step size

      if (!dYTooSmall && !dYTooLarge)      break; // success, not too large, not too small
      if (k>1 &&  reduceDY && !dYTooLarge) break; // success, was too large before but not anymore (do not check if it became too small now)
      if (k>1 && !reduceDY && !dYTooSmall) {      // success, was too small before but not anymore
	if (dYTooLarge) { // if it became too large, go back one step (requires k>1)
	  s[i+1].Parameter(i) = s0.Parameter(i) + previousStep;
	  s[i+1].y            = previousY;
	}
	break; 
      }

      const int kMax = 5;
      if (k>=kMax) {
	cout << " Could not adjust step size that dY became of the right magnitude. Breaking out anyway." << endl;
	if (reduceDY==false && dYTooLarge) { // if we tried increase dY but the last step went too far, go back to previous step
	  s[i+1].Parameter(i) = s0.Parameter(i) + previousStep;
	  s[i+1].y            = previousY;
	}
	break;
      }

      previousY    = s[i+1].y;
      previousStep = step;
      if (reduceDY) step *=  -0.1; // change direction and reduce step size quite a bit and try again	
      else          step *=  +2.0;
    }
  }
  return nFunc;
}

// see http://extremelearning.com.au/how-to-generate-uniformly-random-points-on-n-spheres-and-n-balls/
void SampleNSphere(Array1 <double> & x, const int n) {
  double d = 0.0;
  for(int i=0; i<n; i++) {
    x[i] = SampleGaussian(1.0); // Box-Mueller formula in MC.h
    d += sqr(x[i]);
  }
  d = 1.0/sqrt(d);
  for(int i=0; i<n; i++) {
    x[i] *= d;
  }
}

template <class S> 
int SetUpNWalkersGeneral(S & s0, Array1 <S> & s, const QMCInfo & qmc, const bool print) {
  //  const bool print = false;
  //  const bool print = true;
  int nDim     = s0.nDim();
  if (print) cout << "Setting up " << qmc.nWalkers << " walkers in " << nDim << " dimensions from single state." << endl;

  if (!s0.Valid()) {
    cout << s0 << endl;
    s0.Valid(true);
    error("Cannot set up simplex starting with an invalid parameter set.");
    //    warning("Should not start MC with an invalid parameter set. May be a benign S2<S1 issue.");
  }

  //  int nWalkers = nDim+1;
  s.Resize(qmc.nWalkers); 
  
  s[0] = s0;       // State class needs copy properly. 
  s[0].Evaluate(); // evaluate s0, this takes time
  int nFunc = 1;
  double y0 = s[0].y;
  if (print) Write(y0);
  if (print) s[0].PrintParametersAndStepSizes();

  bool minDYFlag = (qmc.minDeltaYAtStart>0.0) || (qmc.minDeltaYFractionAtStart>0.0); // true if we want to try to generate an ensemble with similar y values
  double minDY1 = qmc.minDeltaYFractionAtStart*abs(s[0].y); // y0 may be less than 0 but it always provides a scale. 
  double minDY2 = qmc.minDeltaYAtStart;                     
  double minDY  = max(minDY1,minDY2); // take the larger of the two because do the same below

  bool maxDYFlag = (qmc.maxDeltaYAtStart>0.0) || (qmc.maxDeltaYFractionAtStart>0.0); // true if we want to try to generate an ensemble with similar y values
  double maxDY1 = qmc.maxDeltaYFractionAtStart*abs(s[0].y); // y0 may be less than 0 but it always provides a scale. 
  double maxDY2 = qmc.maxDeltaYAtStart;                     
  double maxDY  = max(maxDY1,maxDY2); // take the larger of the two. This solve the issue of y~1e-5 and if maxDeltaYFractionAtStart or maxDeltaYAtStart are negative. If both are <0, we need to use maxYAtStartFlag 

  Array1 <double> sphere(nDim);
  for(int i=0; i<qmc.nWalkers-1; i++) { // set up initial tetrahedron of nDim points in addition to original one
    s[i+1] = s0; // here we make copies. Make sure the state copies properly.

    SampleNSphere(sphere,nDim);
    for(int j=0; j<nDim; j++) { 
      sphere[j] *=  s0.ParameterStep(j); // scale every direction
    }

    double step         = 1.0;
    double previousStep = step;
    double previousY    = y0;
    int k=0;
    bool reduceDY = true;
    while (true) {
      k++; // increases counter k
      //  pp(i+1,i) = s[i] + step; // do not use s[i] because it may have changed by the call 's(pp[i+1])'
      for(int j=0; j<nDim; j++) { 
	s[i+1].Parameter(j) = s0.Parameter(j) + sphere[j]*step;
      }
      s[i+1].Evaluate(); 
      nFunc++;
      // y[i+1]    = s(pp[i+1]); // warning, depending how State is written, this call may copy pp[i+1] over into state 's' !!!!!!!!!!!
      double dy = s[i+1].y - y0;
      if (print) { cout << " ::::::"; write6(i+1,k,step,y0,s[i+1].y,dy); cout << " range= " << minDY << " ... " << maxDY << endl; }
      //      cout << endl; // print s[i+1] last because the line is sooo long

      if (dy<0.0) { // Surprise, we found a new minimum. Keep it and make it the new state to compare with. (May be lost otherwise.)
	s[0]   = s[i+1];   // needs to be copied into s[0], not into s0, to be used in simplex method
	s[0].y = s[i+1].y;
	y0     = s[i+1].y; // since we write 'y0' for convenience, overwrite this value also
	step *=  +2.0;
	k--; // reduce so that we do one more step after this one. Do not reset k to 0 because the this routine may never terminate otherwise. 
	continue;
      }

      bool dYTooSmall = minDYFlag && (abs(dy)<minDY);
      bool dYTooLarge = maxDYFlag && (dy     >maxDY); // large negative dY are ok
      if (k==1) reduceDY = dYTooLarge; // decide once (in the first step) whether we need to increase or decrease the step size

      if (!dYTooSmall && !dYTooLarge)      break; // success, not too large, not too small
      if (k>1 &&  reduceDY && !dYTooLarge) break; // success, was too large before but not anymore (do not check if it became too small now)
      if (k>1 && !reduceDY && !dYTooSmall) {      // success, was too small before but not anymore
	if (dYTooLarge) { // if it became too large, go back one step (requires k>1)
	  for(int j=0; j<nDim; j++) { 
	    s[i+1].Parameter(j) = s0.Parameter(j) + sphere[j]*previousStep;
	  }
	  s[i+1].y            = previousY;
	}
	break; 
      }

      const int kMax = 5;
      if (k>=kMax) {
	cout << " Could not adjust step size that dY became of the right magnitude. Breaking out anyway." << endl;
	if (reduceDY==false && dYTooLarge) { // if we tried increase dY but the last step went too far, go back to previous step
	  for(int j=0; j<nDim; j++) { 
	    s[i+1].Parameter(j) = s0.Parameter(j) + sphere[j]*previousStep;
	  }
	  s[i+1].y            = previousY;
	}
	break;
      }

      previousY    = s[i+1].y;
      previousStep = step;
      if (reduceDY) step *=  -0.1; // change direction and reduce step size quite a bit and try again	
      else          step *=  +2.0;
    }
  }
  return nFunc;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////

template <class S> 
// int RunMC(S & s0, Array1 <S> s, const QMCInfo & qmc, Array1 <ScalarAverage> & energyAverage,  Array1 < Array1 <ScalarAverage> > & parameterAverage, const bool estimateTravelTimeFlag) {
int RunMC(S & s0, Array1 <S> s, const QMCInfo & qmc, Array1 <ScalarAverageSmall> & energyAverage,  Array1 < Array1 <ScalarAverageSmall> > & parameterAverage, const bool estimateTravelTimeFlag) {
  Analysis analysis(qmc.blockMin);

  //  int nWalkers = s0.nDim() + 1;
  //  if (nWalkers != s.Size()) error("RunQuadraticMC: nWalkers != s.Size()");
  if (qmc.nWalkers<3) error("nWalkers too small",qmc.nWalkers);

  const string line = " -----------------------------------------------------------------------------------";
  if (qmc.print) cout << endl << " ======= Starting Monte Carlo run with " << s0.nDim() << " dimensions and " << qmc.nWalkers << " walkers. =======" << endl << endl;

  for(int j=0; j<s.Size(); j++) { // loop over all walkers 'j'
    analysis.AddTotalEnergy(s[j].y);
    energyAverage[0].Add(s[j].y);
    //    for(int i=0; i<s[j].nDim(); i++) { // loop over all 'n' parameters of walker s[j]
    for(int i=0; i<parameterAverage.Size(); i++) { // loop over all 'n' parameters of walker s[j] but stop at size of 'parameterAverage' b/c we may only track the first parameter
      parameterAverage[i][0].Add(s[j].Parameter(i));
    }
  }

  int blockStep = 1;
  int k = 0;
  for(int iBlock=0; iBlock<qmc.nBlocks; iBlock++) { 
    //    if (qmc.print) cout << line << endl; 

    analysis.StartBlock(iBlock);
    AccRatio accRatio;
    for (int iStep=0; iStep<qmc.nSteps; iStep++) {

      if (qmc.flagWalk)   PerformWalkMoves                  (s,qmc.nWalkers,qmc.temp,accRatio,qmc.a,qmc.fraction,qmc.nMovers); 
      if (qmc.flagAffine) PerformOneAffineMonteCarloSweep   (s,qmc.nWalkers,qmc.temp,accRatio,qmc.a); 
      if (qmc.flagQMC)    PerformOneQuadraticMonteCarloSweep(s,qmc.nWalkers,qmc.temp,accRatio,qmc.a,qmc.flagGaussian); 

      for(int j=0; j<s.Size(); j++) {
	analysis.AddTotalEnergy(s[j].y);
	energyAverage[iBlock+1].Add(s[j].y);
	//	for(int i=0; i<s[j].nDim(); i++) { // loop over all 'n' parameters in walker s[j]
	for(int i=0; i<parameterAverage.Size(); i++) { // loop over all 'n' parameters of walker s[j] but stop at size of 'parameterAverage' b/c we may only track the first parameter
	  parameterAverage[i][iBlock+1].Add(s[j].Parameter(i)); // add once per step
	}
      }
    } // iStep
    
    // add only once per block, here at the end 
    //    for(int j=0; j<s.Size(); j++) { // loop over all walkers
    //      for(int i=0; i<s[j].nDim(); i++) { // loop over all 'n' parameters in walker s[j]
    //	parameterAverage[j][i][iBlock+1].Add(s[j].Parameter(i));
    //      }
    //    }

    if (estimateTravelTimeFlag && parameterAverage[0][iBlock+1].GetAverage()>0.0) { // assume we travel from -R with high energy to +R with low energy.
      return iBlock;
    }

    if (qmc.print && iBlock%blockStep==0) {
      if ((++k)%10==0) blockStep *= 2; // double step size once we have printed 10 times

      cout << " " << IntToStringMaxNumber(iBlock,qmc.nBlocks) << " : ";
      //      cout << accRatio;
      cout << "ratio= " << fix(6,4,accRatio.Ratio()); // remember that we do not print every block (blockStep>0)
      //      cout << "ratio= " << fix(8,6,accRatio.Ratio()); // adds no new information because we only have qmc.nSteps per block

      //      cout << IntToStringMaxNumber(iBlock,nBlocks) << " : ";
      //      analysis.EndBlock(iBlock,qmc.print);
      analysis.EndBlock(iBlock,false);
      analysis.PrintBlockEnergy();
      cout << endl;

      //    Quit("after first MC sweep");
    }

  } //iBlock

  if (qmc.print) {
    //  if (true) {
    cout << line << endl; 
    cout << line << endl; 
    analysis.EndSimulation();
    if (qmc.nRuns==1) { // at the moment, we do not look at this
      cout << endl;
      cout << " Last step:" << endl;
      cout << endl;
      //  cout << s;
      for(int i=0; i<s.Size(); i++) {
	cout << " state " <<  IntToStringMaxNumber(i,s.Size()) << " : " << s[i] << endl;
      }
    }
    cout << endl << " ======= End of MC chain ======= " << endl << endl;
  }

  return -1;
}

/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////

// pass in an array because we have a separate average for every block
inline Array1 <double> ExtractArray(const Array1 <ScalarAverageSmall> & sa) {
  const int n = sa.Size();
  Array1 <double> x(n);
  for(int i=0; i<n; i++) {
    x[i] = sa[i].GetAverage();
  }
  return x;
}

void BlockingAnalysis(const Array1 <ScalarAverageSmall> & sa, const string s="av") {
  cout << endl;

  Array1 <double> x = ExtractArray(sa);
  const int n2 = x.Size();
  
  int nF = 10;
  for(int iF=nF; iF>0; iF--) { // iF==1 means using the last 10%
    int n1 = n2 - (n2*iF) / nF; // integer division
    int n  = n2 - n1;
    //    Write4(iF,nF,n1,n);
    double av,err;
    //    Form f(fix,14,8);
    Form f(fix,18,12);
    //    double err1; Blocker(x,n1,n2,av,err,err1,true); // debug statement to compare with 'blocker04.c'
    BlockerFast(x,n1,n2,av,err);
    cout << " Fraction= " << fix(3,1,double(iF)/double(nF)) 
	 << " n= " << IntToStringMaxNumber(n,n2)
	 << " " << s << "= " << f(av) << " +- " << f(err) 
	 << endl;
  }
}

void AutocorrelationAnalysis(const Array1 <ScalarAverageSmall> & sa, const string s="") {
  cout << endl;

  Array1 <double> x = ExtractArray(sa);
  const int n2 = x.Size();
  
  //  const int nF = 10;
  const int nF = 5;
  for(int iF=nF; iF>0; iF--) {        // study different fraction starting from end
    const int n1 = n2 - (n2*iF) / nF; // integer division
    const int nn = n2 - n1;           // number of data points to be considered here

    const int nCorrMax = nn / 10; // integer division
    double sum = 0.0;
    double c0 = Autocorrelation(x,0,n1,n2);
    int step = 1;
    int k = 0;
    for(int iCorr=0; iCorr<nCorrMax; iCorr+=step) {
      double c = Autocorrelation(x,iCorr,n1,n2) / c0;
      //      sum += c;
      sum += c*step; // assume 'c' represents autocorrelation values during the last 'step' steps of which step-1 have been skipped over.
      if (s.length()>0) cout << " " << s << " :";
      cout << " f= " << fix(3,1,double(iF)/double(nF)) 
	   << " step= " << IntToStringMaxNumber(iCorr,nCorrMax) 
	   << " ac= " << fix(11,8,c)
	   << " sum= " << fix(12,8,sum) << endl;
      //      if ((++k)%(10*step)==0) step *= 2; // double step size every 10 steps
      if ((++k)%10==0) step *= 2; // double step size every 10 steps
    }
  }
}

/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////

int main(int argc, char **argv) {
  Timer timerTotal(Timer::total, true);
  Timer timerUser(Timer::user, true);

  cout << endl;
  cout << "Starting a QMC calculation." << endl;
  cout << endl;
  cout << "Process ID: " << getpid() << endl;
  cout << endl;

  Array1 <string> arg = CopyCommandLineArgumentsToArray(argc,argv);
  cout.precision(10);

  if (ParseCommandLineFlagWithOrWithoutDash(arg,"noRandom")) {
    //  InitRandom(false); warning("Using standard random numbers");
    InitRandom(false);
    cout << "Using standard random numbers stream: ";
  } else {
    InitRandom(true);
    cout << "Using randomized random numbers stream: ";
  }
  Write(Random());

  //  int n = 10;
  int n = 3;
  ParseCommandLineArgumentNoSpace(arg,"n=",n);
  Write(n);

#ifdef RING
  double R = 1.0;
  ParseCommandLineArgumentNoSpace(arg,"R=",R);
  double m = 8.0;
  ParseCommandLineArgumentNoSpace(arg,"m=",m);
  double B = 0.01;
  ParseCommandLineArgumentNoSpace(arg,"B=",B);
  Write3(R,m,B);
#endif // RING

#ifdef ROSENBROCK
  double A = 100.0;
  ParseCommandLineArgumentNoSpace(arg,"A=",A);
  double B = 20.0;
  ParseCommandLineArgumentNoSpace(arg,"B=",B);
  Write2(A,B);
#endif

  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  QMCInfo qmc;
  //  qmc.print = false;
  qmc.print = true;

  qmc.flagWalk = false;
  qmc.flagWalk = ParseCommandLineFlag(arg,"-walk");
  qmc.flagAffine = false;
  qmc.flagAffine = ParseCommandLineFlag(arg,"-affine");
  qmc.flagQMC = false;
  qmc.flagQMC = ParseCommandLineFlag(arg,"-qmc");
  qmc.flagGaussian = false;
  qmc.flagGaussian = ParseCommandLineFlag(arg,"-Gaussian");

  if (ParseCommandLineFlag(arg,"-qmcLin")) {
    if (qmc.flagQMC) error("-qmc and -qmcLin given");
    qmc.flagQMC = true;
    qmc.flagGaussian = false;
  }
  if (ParseCommandLineFlag(arg,"-qmcGau")) {
    if (qmc.flagQMC) error("-qmc and -qmcGau given");
    qmc.flagQMC = true;
    qmc.flagGaussian = true;
  }

  qmc.a = (qmc.flagWalk) ? 1.0 : 1.5;
  ParseCommandLineArgumentNoSpace(arg,"a=",qmc.a);

  //  double temp = 0.1;
  //  double temp = 0.01; // so that the default of maxDeltaYAtStart becomes 0.01
  //  double temp = 1.0;
  qmc.temp = 0.001;
  ParseCommandLineArgumentNoSpace(arg,"T=",qmc.temp);
  //  if (a<=1.0) error("'a' too small",a);
  Write2(qmc.temp,qmc.a);

  qmc.maxDeltaYFractionAtStart = qmc.temp;     
  qmc.minDeltaYFractionAtStart = qmc.maxDeltaYFractionAtStart/10.0; 

  qmc.maxDeltaYAtStart         = qmc.temp;      
  qmc.minDeltaYAtStart         = qmc.temp/10.0; 

  ParseCommandLineArgumentNoSpace(arg,"maxDeltaYFractionAtStart=",qmc.maxDeltaYFractionAtStart);
  ParseCommandLineArgumentNoSpace(arg,"maxDeltaYAtStart=",qmc.maxDeltaYAtStart);
  Write2(qmc.maxDeltaYFractionAtStart,qmc.maxDeltaYAtStart);
  ParseCommandLineArgumentNoSpace(arg,"minDeltaYFractionAtStart=",qmc.minDeltaYFractionAtStart);
  ParseCommandLineArgumentNoSpace(arg,"minDeltaYAtStart=",qmc.minDeltaYAtStart);
  Write2(qmc.minDeltaYFractionAtStart,qmc.minDeltaYAtStart);

  qmc.nRuns   = 1;
  //  int nRuns   = 100;
  ParseCommandLineArgumentNoSpace(arg,"nRuns=",qmc.nRuns);

  //  qmc.nBlocks = 250;
  qmc.nBlocks = 1000;
  ParseCommandLineArgumentNoSpace(arg,"nBlocks=",qmc.nBlocks);
  qmc.blockMin = qmc.nBlocks / 20;

#ifdef RING
  bool lowFlag = false; // start from high-energy point by default
  lowFlag = ParseCommandLineFlag(arg,"-low");
#endif

  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  //  qmc.nWalkers = -1;
  qmc.nWalkers = n+2; // nWalkers=n+1 does not seem to work
  bool nWalkersFlag = ParseCommandLineArgumentNoSpace(arg,"nWalkers=",qmc.nWalkers);
  Write2(n,qmc.nWalkers);

  // just for 'walk' moves
  qmc.fraction=0.5;
  bool fractionsFlag = ParseCommandLineArgumentNoSpace(arg,"fraction=",qmc.fraction);

  qmc.nMovers = 0;
  bool nMoversFlag = ParseCommandLineArgumentNoSpace(arg,"nMovers=",qmc.nMovers); // Was called 'nSubset' before

  int nStepsPerBlock = 1000;
  ParseCommandLineArgumentNoSpace(arg,"nStepsPerBlock=",nStepsPerBlock); // 's' was missing from string 

  double fraction2 = 0.0;
  bool fractionsFlag2 = ParseCommandLineArgumentNoSpace(arg,"fraction2=",fraction2);
  if (nMoversFlag && fractionsFlag2) error("nMoversFlag && fractionsFlag2");
  if (fractionsFlag2) {
    nMoversFlag = true;
    qmc.nMovers = 2 + int(round(fraction2*(qmc.nWalkers-3))); // fraction2==0.0 means nMovers=2=min. // fraction2==1.0 means nMovers=nWalkers-1=max. 
    Write2(fraction2,qmc.nMovers);
  }

  if (fractionsFlag && nMoversFlag) error("fractionsFlag && nMoversFlag");
  if (nMoversFlag) qmc.fraction=0.0;
  Write2(qmc.fraction,qmc.nMovers);

  bool estimateTravelTimeFlag = ParseCommandLineFlag(arg,"-estimateTravelTime");

  CheckForUnusedCommandLineParameters(arg);

  ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

  Array1 <ScalarAverageSmall> energyAverage(qmc.nBlocks+1);
  Array1< Array1 <ScalarAverageSmall> > parameterAverage(n); // [parameter index][block index]
  for(int j=0; j<n; j++) { 
    parameterAverage[j].Resize(qmc.nBlocks+1); // will be initialized to zero automically
  }
  /*
  Array1 < Array1< Array1 <ScalarAverage> > > parameterAverage(qmc.nWalkers);// [walker index][parameter index][block index]
  for(int i=0; i<qmc.nWalkers; i++) { 
    parameterAverage[i].Resize(n);
    for(int j=0; j<n; j++) { 
      parameterAverage[i][j].Resize(qmc.nBlocks+1); // will be initialized to zero automically
    }
  }
  */

  //////////////////////////////////////////////////////////////////////////////////////

  if (estimateTravelTimeFlag) {
#ifdef HARMONIC
    error("estimateTravelTimeFlag and HARMONIC");
#endif
#ifdef ROSENBROCK
    error("estimateTravelTimeFlag and ROSENBROCK");
#endif
#ifdef RING
    if (lowFlag) error("estimateTravelTimeFlag and lowFlag should not both be true.");

    cout << " ------ test run to estimate travel time ------" << endl;

    QMCState s0(n,R,m,B);
    for(int i=0; i<n; i++) { // Unfortunately code is replicated below
      s0.Parameter(i)     = (i>0) ? 0.0 : -R; // start from high-energy point
      s0.ParameterStep(i) = 0.1;
    }

    Array1 <QMCState> s; // Unfortunately code is replicated below
    if (nWalkersFlag==false) {
      SetUpSimplex(s0,s,qmc,true);
      qmc.nWalkers = s.Size();
    } else {
      SetUpNWalkersGeneral(s0,s,qmc,true);
    }

    // set this so late because currently SetUpSimplex() will set nWalkers=nDim+1
    qmc.nSteps = nStepsPerBlock / double(qmc.nWalkers); // steps per block

    int iBlockZero = RunMC(s0,s,qmc,energyAverage,parameterAverage,true); // true to step when parameter[0][:]=0
    //    int nBlocksNew = 2*iBlockZero;
    int nBlocksNew = 4*iBlockZero;
    if (iBlockZero>0 && nBlocksNew<qmc.nBlocks) {
      qmc.nBlocks = nBlocksNew;
      cout << " Revising block number: "; 
      Write2(iBlockZero,qmc.nBlocks);
      int nb = qmc.nBlocks+1;
      energyAverage.Resize(nb);
      for(int i=0; i<nb; i++) 
	energyAverage[i].SetToZero();
      for(int i=0; i<n; i++) {
	parameterAverage[i].Resize(nb);
	for(int j=0; j<nb; j++) 
	  parameterAverage[i][j].SetToZero();
      }
      qmc.blockMin = min(qmc.nBlocks/2,qmc.blockMin); // avoid case qmc.blockMin>qmc.nBlocks/2
    } else {
      cout << " Not revising block number: "; 
      Write2(iBlockZero,qmc.nBlocks);
    }
    
#endif
  }

  //////////////////////////////////////////////////////////////////////////////////////

  for(int iRun=0; iRun<qmc.nRuns; iRun++) {

#ifdef HARMONIC

    QMCState s0(n);
    for(int i=0; i<n; i++) {
      s0.Parameter(i)     = Random();
      s0.ParameterStep(i) = 0.1;
    }

#endif // HARMONIC

#ifdef RING

    QMCState s0(n,R,m,B);
    for(int i=0; i<n; i++) {
      if (lowFlag==false) {
	s0.Parameter(i)     = (i>0) ? 0.0 : -R; // start from high-energy point
      } else {
	s0.Parameter(i)     = (i>0) ? 0.0 : +R; // start from low-energy point
      }
      s0.ParameterStep(i) = 0.1;
    }

#endif // RING

#ifdef ROSENBROCK

    QMCState s0;
    n = s0.nDim(); // 'n' should be 2
    s0.SetAAndB(A,B);

#endif // ROSENBROCK

    const bool printDuringSetup = (iRun<10); // only for the first 10 times, print more info during setup of walkers.
    Array1 <QMCState> s;
    if (nWalkersFlag==false) {
      SetUpSimplex(s0,s,qmc,printDuringSetup);
      qmc.nWalkers = s.Size();
    } else {
      SetUpNWalkersGeneral(s0,s,qmc,printDuringSetup);
    }

    // set this so late because currently SetUpSimplex() will set nWalkers=nDim+1
    qmc.nSteps = nStepsPerBlock / double(qmc.nWalkers); // steps per block

    if (iRun==0) {
      cout << endl << qmc;
    }

    RunMC(s0,s,qmc,energyAverage,parameterAverage,false);

    if (qmc.nRuns==1) { // We typically do not look at this for runs with qmc.nRuns>1
      AutocorrelationAnalysis(energyAverage,"E");
      for(int i=0; i<n; i++) 
	AutocorrelationAnalysis(parameterAverage[i],"p["+IntToString(i+1)+"]");
      //    for(int j=0; j<qmc.nWalkers; j++) 
      //      for(int i=0; i<n; i++) 
      //	AutocorrelationAnalysis(parameterAverage[j][i],"p["+IntToString(j+1)+"]["+IntToString(i+1)+"]");
      
      BlockingAnalysis(energyAverage,"E");
      for(int i=0; i<n; i++) 
	BlockingAnalysis(parameterAverage[i],"p["+IntToString(i+1)+"]");
      //    for(int j=0; j<qmc.nWalkers; j++) 
      //      for(int i=0; i<n; i++) 
      //	BlockingAnalysis(parameterAverage[j][i],"p["+IntToString(j+1)+"]["+IntToString(i+1)+"]");
    }

  }

  //  cout << endl;

  if (qmc.nRuns>1) { // needed for travel time calculations **7341** and **7344**
    for(int iBlock=0; iBlock<=qmc.nBlocks; iBlock++) {
      Form f(fix,12,8);
      cout << " Block= " << IntToStringMaxNumber(iBlock,qmc.nBlocks);
      cout << " E= "    << f(energyAverage[iBlock].GetAverage())       << " var= " << f(energyAverage[iBlock].GetVariance());
      cout << " p[0]= " << f(parameterAverage[0][iBlock].GetAverage()) << " var= " << f(parameterAverage[0][iBlock].GetVariance()) << endl; // only print x=p[0] for travel time calculation
    }
    cout << endl;
  }

  cout << endl << " Program finished properly after " << fix(6,2,timerTotal.Read()) << " ( " << fix(6,2,timerUser.Read()) << " ) seconds. Happy QMC computing!" << endl;
  cerr << endl << " Program finished properly after " << fix(6,2,timerTotal.Read()) << " ( " << fix(6,2,timerUser.Read()) << " ) seconds. Happy QMC computing!" << endl;
  return 0;
}