//
// PopulationDynamicsSATbp.cpp
//
// This program complements Section 6.2 of the book
// "Spin Glass and Message Passing" by Hai-Jun Zhou
// (Science Press, Beijing, 2015).
// Program last modified on 13.08.2012
//   
// DESCRIPTION:
// This program works for the random K-SAT problem. Each clause has
// degree K; the degree of the variable nodes follow the Poisson
// distribution.
//
// Population dynamics simulation of the entropic belief-propagation
// to obtain the mean values of entropy density as a function of the
// clause density alpha. The mean value of the Edwards-Anderson
// order parameter as a function of alpha is also calculated.
//
// compile this program using
// c++ -O3 -o SATpopdyn.exe PopulationDynamicsSATbp.cpp
//
// Programmer:
// Hai-Jun Zhou
// Institute of Theoretical Physics, Chinese Academy of Sciences,
// Beijing 100190, China
// Webpage http://power.itp.ac.cn/~zhouhj
// Email zhouhj@itp.ac.cn
//

#include <cmath>
#include <exception>
#include <fstream>
#include <iostream>
#include <string>
#include <valarray>
#include "zhjrandom.h"

using namespace std;


class PopDynSATbp
{
public:
  PopDynSATbp(ZHJRANDOMv3*, int, int, double);
  bool PopulationInitialization( string&); //read a population file
  bool PopulationInitialization(void); //population initialization (paramagnetic)
  bool PopulationDynamics(void); //population dynamics for reaching steady state
  bool PopulationDynamics(const int,  string&); //population dynamics and calculation
  void ReportPopulation( string&); //output population
  void ReportMagnetizationPopulation( string&); //output magnetization population
  bool Bootstrap(const int,  string&,  string&); //Bootstrap analysis of simulation data

private:
  double Alpha; //clause density
  int Clause_degree; //=K
  int Clause_cavity_degree; //=K-1
  int Population_size; //number of elements in the population
  double Variable_mean_cavity_degree_p; //=alpha*K/2, variable mean positive degree
  double Variable_mean_cavity_degree_n; //=alpha*K/2, variable mean negative degree

  ZHJRANDOMv3 *rdptr; //random number generator

  valarray<double> Loglikelihoods; //population of cavity messages

  bool PopulationUpdate(const int&);
  bool PopulationUpdate(const int&, double&, double&, double&);
};


int main(int argc, char ** argv)
{
  int rdseed;
  cout<<"rdseed= "; cout.flush();
  cin>>rdseed;

  ZHJRANDOMv3 rdgenerator(rdseed);
  for(int i=0; i<rdseed; ++i) rdgenerator.rdflt();

  double alpha;
  cout<<"Constraint density (>=0) = "; cout.flush();
  cin>>alpha;

  int clause_degree;
  cout<<"K-sat with K= "; cout.flush();
  cin>>clause_degree;

  int pop_size;
  cout<<"Population size = "; cout.flush();
  cin>>pop_size;

  PopDynSATbp bp_zhj(&rdgenerator, pop_size, clause_degree, alpha);

  string outputfilename;
  cout<<"entropy output file = "; cout.flush();
  cin>>outputfilename;

  string likelihoodfile;
  cout<<"population output file = "; cout.flush();
  cin>>likelihoodfile;

  int intype=0;
  cout<<"Population initialize style (0 (para) /1 (from file) / = ";
  cout.flush();
  cin>>intype;

  bool bpgood=true;
  if(intype==0) bpgood=bp_zhj.PopulationInitialization();
  else bpgood=bp_zhj.PopulationInitialization(likelihoodfile);

  if(bpgood==false) return -1;

  cout<<"Population equilibrium:\n"; cout.flush();
  int iteration_eq;
  cout<<"number of equilibrium steps (200 could be enough) = "; cout.flush();
  cin>>iteration_eq;

  for(int iteration=1; iteration<=iteration_eq && bpgood; ++iteration)
    {
      bpgood=bp_zhj.PopulationDynamics();
      cout<<iteration<<endl;
    }

  cout<<"Calculation:\n"; cout.flush();
  int iteration_cl;
  cout<<"number of calculation steps (2000 might be enough) = "; cout.flush();
  cin>>iteration_cl;

  for(int iteration=1; iteration<=iteration_cl && bpgood; ++iteration)
    {
      bpgood=bp_zhj.PopulationDynamics(iteration, outputfilename);
      cout<<iteration<<endl;
    }

  if(bpgood)
    {
      bp_zhj.ReportPopulation(likelihoodfile);
      //      bp_zhj.ReportMagnetizationPopulation("temp.dat");

      cout<<"Bootstrap analysis of simulation data.\n";
      cout<<"Mean value output file name = ";
      string meanvaluefile;
      cin>>meanvaluefile;
      bp_zhj.Bootstrap(iteration_cl, outputfilename, meanvaluefile);
      return 1;
    }
  else return -1;
}


//constructor for PopDynSATbp
PopDynSATbp::PopDynSATbp
(ZHJRANDOMv3* rd, int pop_size, int c_degree, double ratio)
{
  Alpha=ratio; //clause density
  Clause_degree=c_degree; //K
  Clause_cavity_degree=Clause_degree-1;
  Population_size=pop_size;
  rdptr=rd;
  Variable_mean_cavity_degree_n=Alpha*Clause_degree*0.5e0;
  Variable_mean_cavity_degree_p=Variable_mean_cavity_degree_n;
  return ;
}


bool PopDynSATbp::PopulationInitialization(void)
{
  try { Loglikelihoods.resize(Population_size); } catch(bad_alloc) { return false; }

  for(int position=0; position<Population_size; ++position)
    Loglikelihoods[position]=2.0e0; 
  /*
    loglikelihood=\log((probability to be positive)/(probability to be negative))
    value>=2  : unfrozen, loglikelihood=value-2
    value<=-2 : unfrozen, loglikelihood=value+2
    value=1   : positively frozen, loglikelihood=infinity
    value=-1  : negatively frozen, loglikelihood=-infinity
  */
  return true;
}


bool PopDynSATbp::PopulationInitialization( string& filename)
{
  try { Loglikelihoods.resize(Population_size); } catch(bad_alloc) { return false; }

  int pop;
  ifstream input(filename.c_str());
  if(!input.good())
    {
      cerr<<"population file probably non-existant.\n";
      return false;
    }
  input >> pop;
  if(pop!=Population_size)
    {
      input.close();
      cerr<<"population file does not contain the same number of items as required.\n";
      return false;
    }

  double value;
  for(int position=0; position<Population_size; ++position)
    {
      input>>value;
      Loglikelihoods[position]=value;
    }
  input.close();
  return true;
}


//output population
void PopDynSATbp::ReportPopulation( string& filename)
{
  ofstream output(filename.c_str() );
  output.precision(12);
  output << Population_size<<endl<<endl;

  double *log_ptr=&Loglikelihoods[0];
  for(int position=0; position<Population_size; ++position, ++log_ptr)
    output<< *log_ptr <<endl;
  output.close();
  return ;
}


//output magnetization population
void PopDynSATbp::ReportMagnetizationPopulation( string& filename)
{
  ofstream output(filename.c_str() );
  output.precision(10);
  for(int position=0; position<Population_size; ++position)
    {
      double value=Loglikelihoods[position];
      if(value>=1.5) //value>=2
	{
	  value=exp(-(value-2.0e0));
	  output<<(1.0e0-value)/(1.0e0+value)<<endl;
	}
      else if(value<=-1.5) //value<=-2
	{
	  vlaue=exp(value+2.0e0);
	  output<<(value-1.0e0)/(value+1.0e0)<<endl;
	}
      else if(value>0.5) output<<1<<endl; //value=1
      else output<<-1<<endl; //value=-1
    }
  output.close();
  return ;
}


//population dynamics to reach a steady state
bool PopDynSATbp::PopulationDynamics(void)
{
  bool bpgood=true;
  for(int position=0; position<Population_size && bpgood; ++position)
    bpgood=PopulationUpdate(position);

  return bpgood;
}


//population dynamics to calculate entropy and self-overlap
bool PopDynSATbp::PopulationDynamics(const int iteration_index,  string& filename)
{
  double
    S_i =0.0e0,
    S_a =0.0e0,
    mean_qEA =0.0e0,
    delta_S_i,
    delta_S_a,
    delta_qEA;

  bool bpgood=true;
  for(int position=0; position<Population_size && bpgood; ++position)
    {
      bpgood=PopulationUpdate(position, delta_S_i, delta_S_a, delta_qEA);
      if(bpgood)
	{
	  S_i += delta_S_i;
	  S_a += delta_S_a;
	  mean_qEA += delta_qEA;
	}
    }
  if(bpgood)
    {
      S_i /= Population_size;
      S_a /= Population_size;
      mean_qEA /= Population_size;

      ofstream output(filename.c_str(), ios_base::app);
      output
	<< iteration_index    <<'\t'
	<< S_i-Alpha*Clause_cavity_degree*S_a <<'\t'
	<< mean_qEA             <<endl;
      output.close();
    }
  return bpgood;
}


//update one element of the population array Loglikelihoods
bool PopDynSATbp::PopulationUpdate(const int& position0)
{
  double
    psum=0, //=log(1-prod_{k in \partial c\i} [1-J_c^k m_{k->c})/2]), with J_c^i=+1
    msum=0; //=log(1-prod_{j in \partial b\i} [1-J_b^j m_{j->b})/2]), with J_b^i=-1

  int
    punsat=0, //number of warning messages from positive edges
    munsat=0, //number of warning messages from negative edges
    variable_cavity_degree_p=static_cast<int>(rdptr->poidev(Variable_mean_cavity_degree_p)),
    variable_cavity_degree_n=static_cast<int>(rdptr->poidev(Variable_mean_cavity_degree_n));

  if(variable_cavity_degree_p>=1) //variable node is is attached with one or more positive edges
    for(int clause_c=0; clause_c<variable_cavity_degree_p; ++clause_c)
      {
	double
	  expv=0, //expv= -\sum_{k\in \partial c\i} log[(1-I_kc)/2+e^{2-eta_kc} (1+I_kc)/2]
	  minv=0,
	  resv=0; //prod_{k\parital c\i}(1+e^{2-eta_kc})=1+exp(-minv)*resv

	for(int variable_k=0; variable_k<Clause_cavity_degree; ++variable_k)
	  {
	    //consider all the other neighboring variables k of clause c besides i
	    int I_kc=(rdptr->rdflt()<0.5 ? +1 : -1);
	    double eta_kc=Loglikelihoods[static_cast<int>(Population_size*rdptr->rdflt())];
	    if(eta_kc<0)
	      {
		eta_kc *= -1;
		I_kc *= -1;
	      } //eta_kc=2*(cavity field) on variable k. If it is negative, a gauge transform is applied

	    if(eta_kc>1.5) //variable k not frozen in absence of c, (eta_kc>=2)
	      {
		double val=1.0e0+exp(-(eta_kc-2.0e0)); //m_{k->c}=[1-exp(2-eta_kc)]/[1+exp(2-eta_kc)]
		if(minv>1.5) //minv>=2
		  {
		    if(I_kc>0) expv += eta_kc-2.0e0; 
		    if(minv<=eta_kc) resv=val*resv+exp(minv-eta_kc);
		    else {
		      resv=1.0+exp(eta_kc-minv)*val*resv;
		      minv=eta_kc;
		    }
		  }
		else if(minv < 0.5) // not yet initialized (minv=0)
		  {
		    minv=eta_kc;
		    expv=(I_kc>0? eta_kc-2.0 : 0);
		    resv=1.0e0;
		  }
	      }
	    else // eta_kc=1, frozen (the direction is defined as positive)
	      {
		if(I_kc>0) minv=1; //clause for sure satisfied by k, expv=+infinity
		//else I_kc<0, clause can not satisfied by variable k
	      }
	  }

	if(minv>=1.5e0) //minv>=2
	  psum += (expv>0? (log(1.0e0-exp(-expv)/(1.0e0+exp(2.0e0-minv)*resv))) :
		   (2.0e0-minv+log(resv/(1.0e0+exp(2.0e0-minv)*resv))) );
	else if(minv < 0.5) //==0, clause can not be satisfied
	  ++punsat;
	// else minv=1, clause satisfied
      }

  if(variable_cavity_degree_n>=1) //variable node i is attached with one or more negative edges
    for(int clause_b=0; clause_b<variable_cavity_degree_n; ++clause_b)
      {
	double
	  expv=0, //expv= -\sum_{j\in \partial b\i} log[(1-I_jb)/2+e^{2-eta_jb} (1+I_jb)/2]
	  minv=0,
	  resv=0; //prod_{j\parital b\i}(1+e^{2-eta_jb})=1+exp(-minv)*resv

	for(int variable_j=0; variable_j<Clause_cavity_degree; ++variable_j)
	  {
	    //consider all the other neighboring variables j of clause b
	    int I_jb=(rdptr->rdflt()<0.5 ? +1 : -1);
	    double eta_jb=Loglikelihoods[static_cast<int>(Population_size*rdptr->rdflt())];
	    if(eta_jb<0)
	      {
		eta_jb *= -1;
		I_jb *= -1;
	      } //eta_jb=2*(cavity field) on variable j. A gauge transform is applied if it is negative

	    if(eta_jb>1.5) //variable j not frozen in absence of b, (eta_jb>=2)
	      {
		double val=1.0e0+exp(-(eta_jb-2.0e0)); //m_{j->b}=[1-exp(2-eta_jb)]/[1+exp(2-eta_jb)]
		if(minv>1.5) //minv>=2
		  {
		    if(I_jb>0) expv += eta_jb-2.0e0; 
		    if(minv<=eta_jb) resv=val*resv+exp(minv-eta_jb);
		    else {
		      resv=1.0+exp(eta_jb-minv)*val*resv;
		      minv=eta_jb;
		    }
		  }
		else if(minv < 0.5) // not yet initialized (minv=0)
		  {
		    minv=eta_jb;
		    expv=(I_jb>0? eta_jb-2.0 : 0);
		    resv=1.0e0;
		  }
	      }
	    else // eta_jb=1, frozen (the direction is defined as positive)
	      {
		if(I_jb>0) minv=1; //clause for sure satisfied by j, expv=+infinity
		//else I_jb<0, clause can not satisfied by variable j
	      }
	  }

	if(minv>=1.5e0) //minv>=2
	  msum += (expv>0? (log(1.0e0-exp(-expv)/(1.0e0+exp(2.0e0-minv)*resv))) :
		   (2.0e0-minv+log(resv/(1.0e0+exp(2.0e0-minv)*resv))) );
	else if(minv < 0.5) //==0, clause can not be satisfied
	  ++munsat;
	// else minv=1, clause satisfied
      }

  double newlog;
  if(punsat>0 && munsat>0)
    {
      cout<<"Contradiction.\n";
      return false;
    }
  if(punsat>0) newlog=1; //cavity field +infinity
  else if(munsat>0) newlog=-1; //cavity field -infinity
  else //punsat=munsat=0
    {
      newlog=msum-psum;
      if(newlog>=0) newlog += 2.0e0;
      else newlog -= 2.0e0;
    }
  Loglikelihoods[position0]=newlog;

  return true;
}


//update one element of the population array Loglikelyhoods
bool PopDynSATbp::PopulationUpdate
(const int& position0, double& deltaSi, double& deltaSa, double& deltaqEA)
{
  double
    psum=0, //=log(1-prod_{k in \partial c\i} [1-J_c^k m_{k->c})/2]), with J_c^i=+1
    msum=0; //=log(1-prod_{j in \partial b\i} [1-J_b^j m_{j->b})/2]), with J_b^i=-1

  int
    punsat=0, //number of warning messages from positive edges
    munsat=0, //number of warning messages from negative edges
    variable_cavity_degree_p=static_cast<int>(rdptr->poidev(Variable_mean_cavity_degree_p)),
    variable_cavity_degree_n=static_cast<int>(rdptr->poidev(Variable_mean_cavity_degree_n));

  if(variable_cavity_degree_p>=1) //variable node i has one or more positive links
    for(int clause_c=0; clause_c<variable_cavity_degree_p; ++clause_c)
      {
	double
	  expv=0, //expv= -\sum_{k\in \partial c\i} log[(1-I_kc)/2+e^{2-eta_kc} (1+I_kc)/2]
	  minv=0,
	  resv=0; //prod_{k\parital c\i}(1+e^{2-eta_kc})=1+exp(-minv)*resv

	for(int variable_k=0; variable_k<Clause_cavity_degree; ++variable_k)
	  {
	    //consider all the other neighboring variables k of clause c except i
	    int I_kc=(rdptr->rdflt()<0.5 ? +1 : -1);
	    double eta_kc=Loglikelihoods[static_cast<int>(Population_size*rdptr->rdflt())];
	    if(eta_kc<0)
	      {
		eta_kc *= -1;
		I_kc *= -1;
	      } //eta_kc=2*(cavity field) on variable k, it is gauge-transformed to be positive if needed

	    if(eta_kc>1.5) //variable k not frozen in absence of c, (eta_kc>=2)
	      {
		double val=1.0e0+exp(-(eta_kc-2.0e0)); //m_{k->c}=[1-exp(2-eta_kc)]/[1+exp(2-eta_kc)]
		if(minv>1.5) //minv>=2
		  {
		    if(I_kc>0) expv += eta_kc-2.0e0; 
		    if(minv<=eta_kc) resv=val*resv+exp(minv-eta_kc);
		    else {
		      resv=1.0+exp(eta_kc-minv)*val*resv;
		      minv=eta_kc;
		    }
		  }
		else if(minv < 0.5) // not yet initialized (minv=0)
		  {
		    minv=eta_kc;
		    expv=(I_kc>0? eta_kc-2.0 : 0);
		    resv=1.0e0;
		  }
	      }
	    else // eta_kc=1, frozen (the direction is defined as positive)
	      {
		if(I_kc>0) minv=1; //clause for sure satisfied by k, expv=+infinity
		//else I_kc<0, clause can not satisfied by variable k
	      }
	  }

	if(minv>=1.5e0) //minv>=2
	  psum += (expv>0? (log(1.0e0-exp(-expv)/(1.0e0+exp(2.0e0-minv)*resv))) :
		   (2.0e0-minv+log(resv/(1.0e0+exp(2.0e0-minv)*resv))) );
	else if(minv < 0.5) //==0, clause can not be satisfied
	  ++punsat;
	// else minv=1, clause satisfied
      }

  if(variable_cavity_degree_n>=1) //variable node i has one or more negative links
    for(int clause_b=0; clause_b<variable_cavity_degree_n; ++clause_b)
      {
	double
	  expv=0, //expv= -\sum_{j\in \partial b\i} log[(1-I_jb)/2+e^{2-eta_jb} (1+I_jb)/2]
	  minv=0,
	  resv=0; //prod_{j\parital b\i}(1+e^{2-eta_jb})=1+exp(-minv)*resv

	for(int variable_j=0; variable_j<Clause_cavity_degree; ++variable_j)
	  {
	    //consider all the other neighboring variables j of clause b
	    int I_jb=(rdptr->rdflt()<0.5 ? +1 : -1);
	    double eta_jb=Loglikelihoods[static_cast<int>(Population_size*rdptr->rdflt())];
	    if(eta_jb<0)
	      {
		eta_jb *= -1;
		I_jb *= -1;
	      } //eta_jb=2*(cavity field) on variable j. A gauge transform is performed

	    if(eta_jb>1.5) //variable j not frozen in absence of b, (eta_jb>=2)
	      {
		double val=1.0e0+exp(-(eta_jb-2.0e0)); //m_{j->b}=[1-exp(2-eta_jb)]/[1+exp(2-eta_jb)]
		if(minv>1.5) //minv>=2
		  {
		    if(I_jb>0) expv += eta_jb-2.0e0; 
		    if(minv<=eta_jb) resv=val*resv+exp(minv-eta_jb);
		    else {
		      resv=1.0+exp(eta_jb-minv)*val*resv;
		      minv=eta_jb;
		    }
		  }
		else if(minv < 0.5) // not yet initialized (minv=0)
		  {
		    minv=eta_jb;
		    expv=(I_jb>0? eta_jb-2.0 : 0);
		    resv=1.0e0;
		  }
	      }
	    else // eta_jb=1, frozen (the direction is defined as positive)
	      {
		if(I_jb>0) minv=1; //clause for sure satisfied by j, expv=+infinity
		//else I_jb<0, clause can not satisfied by variable j
	      }
	  }

	if(minv>=1.5e0) //minv>=2
	  msum += (expv>0? (log(1.0e0-exp(-expv)/(1.0e0+exp(2.0e0-minv)*resv))) :
		   (2.0e0-minv+log(resv/(1.0e0+exp(2.0e0-minv)*resv))) );
	else if(minv < 0.5) //==0, clause can not be satisfied
	  ++munsat;
	// else minv=1, clause satisfied
      }

  double newlog;
  if(punsat>0 && munsat>0)
    {
      cout<<"Contradiction.\n"; cout.flush();
      return false;
    }
  if(punsat>0)
    {
      newlog=1; //cavity field +infinity
      deltaSi=msum;
      deltaqEA=1;
    }
  else if(munsat>0)
    {
      newlog=-1; //cavity field -infinity
      deltaSi=psum;
      deltaqEA=1;
    }
  else //punsat=munsat=0
    {
      newlog=msum-psum;
      if(newlog>=0)
	{
	  double value=exp(-newlog);
	  deltaSi=msum+log(1.0e0+value);
	  double meanspin=(1.0e0-value)/(1.0e0+value);
	  deltaqEA=meanspin*meanspin;
	  newlog += 2.0e0;
	}
      else
	{
	  double value=exp(newlog);
	  deltaSi=psum+log(1.0e0+value);
	  double meanspin=-(1.0e0-value)/(1.0e0+value);
	  deltaqEA=meanspin*meanspin;
	  newlog -= 2.0e0;
	}
    }
  Loglikelihoods[position0]=newlog;

  //calculate the entropy contribution of clause a
  double
    expva=0, //expva= -\sum_{j\in \partial a} log[(1-I_ja)/2+e^{2-eta_ja} (1+I_ja)/2]
    minva=0,
    resva=0; //prod_{j\parital a}(1+e^{2-eta_ja})=1+exp(-minva)*resva

  for(int variable_j=0; variable_j<Clause_degree; ++variable_j)
    {
      //consider all the neighboring variables j of clause a
      int I_ja=(rdptr->rdflt()<0.5 ? +1 : -1);
      double eta_ja=Loglikelihoods[static_cast<int>(Population_size*rdptr->rdflt())];
      if(eta_ja<0)
	{
	  eta_ja *= -1;
	  I_ja *= -1;
	} //eta_ja=2*(cavity field) on variable j, if needed a gauge transform is performed

      if(eta_ja>1.5) //variable j not frozen in absence of a, (eta_ja>=2)
	{
	  double val=1.0e0+exp(-(eta_ja-2.0e0)); //m_{j->a}=[1-exp(2-eta_ja)]/[1+exp(2-eta_ja)]
	  if(minva>1.5) //minv>=2
	    {
	      if(I_ja>0) expva += eta_ja-2.0e0; 
	      if(minva<=eta_ja) resva=val*resva+exp(minva-eta_ja);
	      else {
		resva=1.0+exp(eta_ja-minva)*val*resva;
		minva=eta_ja;
	      }
	    }
	  else if(minva < 0.5) // not yet initialized (minva=0)
	    {
	      minva=eta_ja;
	      expva=(I_ja>0? eta_ja-2.0 : 0);
	      resva=1.0e0;
	    }
	}
      else // eta_ja=1, frozen (the direction is defined as positive)
	{
	  if(I_ja>0) minva=1; //clause for sure satisfied by j, expva=+infinity
	  //else I_ja<0, clause can not satisfied by variable j
	}
    }

  if(minva>=1.5e0) //minva>=2
    deltaSa=(expva>0? (log(1.0e0-exp(-expva)/(1.0e0+exp(2.0e0-minva)*resva))) :
	     (2.0e0-minva+log(resva/(1.0e0+exp(2.0e0-minva)*resva))) );
  else if(minva>0.5) //minva=0, clause satisfied
    deltaSa=0;
  else //minva==0, clause can not be satisfied
    {
      cout<<"Contradiction.\n"; cout.flush();
      return false;
    }

  return true;
}


//Bootstrap analysis of simulation data
bool PopDynSATbp::Bootstrap(const int totalnumber,  string& inputname,  string& outputname)
{
  int totalrepeat=100000; //total number of bootstrap repeat
  //  int totalnumber=2000; //total number of measurements
  int beginsample=0; //number of measurements skipped

  int samplenumber=totalnumber-beginsample;

  valarray<double> v1_profile(0.0,samplenumber);
  valarray<double> v2_profile(0.0,samplenumber);

  ifstream inputf(inputname.c_str() );
  if(!inputf.good())
    {
      cerr<<"File "<<inputname<<" probably non-existent.\n";
      return false;
    }

  for(int i=0; i<beginsample; ++i) 
    {
      int index;
      double a;
      double b;
      inputf>> index >> a >> b ;
    }

  for(int i=0; i<samplenumber; ++i) 
    {    
      int index;
      double a;
      double b;
      inputf>>index>>a>>b;
      v1_profile[i]=a;
      v2_profile[i]=b;
    }
  inputf.close();

  double v1    =0;
  double v1_sd =0;
  double v2    =0;
  double v2_sd =0;
  double value =0;

  for(int repeat=0; repeat < totalrepeat; ++repeat)
    {
      double v1_mean=0;
      double v2_mean=0;

      for(int index=0; index<samplenumber; ++index)
	{
	  int ii=static_cast<int>(samplenumber*rdptr->rdflt());
	  v1_mean += v1_profile[ii];
	  v2_mean += v2_profile[ii];
	}
      v1_mean /= static_cast<double>(samplenumber);
      v2_mean /= static_cast<double>(samplenumber);

      v1    += v1_mean;
      v1_sd += pow(v1_mean,2);
      v2    += v2_mean;
      v2_sd += pow(v2_mean,2);
      value  += 1.0e0;
    }

  v1    /= value;
  v1_sd /= value;
  v2    /= value;
  v2_sd /= value;

  double rescale=sqrt(1.0e0+1.0e0/(value-1.0e0));
  v1_sd = sqrt(abs(v1_sd-v1*v1))*rescale;
  v2_sd = sqrt(abs(v2_sd-v2*v2))*rescale;

  ofstream outputf(outputname.c_str() );
  outputf << inputname<<'\t'<<v1 <<'\t'<< v1_sd <<'\t'
	  << v2 <<'\t'<< v2_sd <<endl;
  outputf.close();  

  return true;
}