/*
  RREA_Complexity_y_vs_beta_1.cpp

  DESCRIPTION:
  Population dynamics simulations to calculate the complexity (at y=beta) as a
  function of temperature.
  
  This program is for the ensemble of random regular networks in which each
  vertex has $K$ attached edges. Each interaction $a$ involves $L$ vertices,
  with a coupling constant $J_a \in {-J, +J}$.

  Some of the results obtained by this program are reported in Section 4.7
  of the book:
  Hai-Jun Zhou, "Spin Glass and Message Passing" 
  (Science Press, Beijing, 2016).

  LOG:
  21.04.2013--01.05.2013
  
  PROGRAMMER:
  Haijun Zhou
  Institute of Theoretical Physics, Chinese Academy of Sciences,
  Zhong-Guan-Cun East Road 55, Beijing 100190, China
  http://power.itp.ac.cn/~zhouhj
  email: zhouhj@itp.ac.cn
*/



#include <cmath>
#include <exception>
#include <fstream>
#include <iostream>
#include <string>
#include <valarray>

//#include "/cygdrive/d/bin/Programs/zhjrandom.h"
//#include "/home/zhouhj/bin/Programs/zhjrandom.h"
#include "zhjrandom.h"

using namespace std;


struct DT //double_triple
{
  double q_i_a_plus_mean; // mean value of q_{i->a}(+)
  double q_i_a_plus_sp;   // q_{i->a}^{(+)}(+)
  double q_i_a_plus_sm;   // q_{i->a}^{(-)}(+)
  DT(double, double, double);
};

DT::DT(double a=0, double b=0, double c=0)
{
  q_i_a_plus_mean=a;
  q_i_a_plus_sp  =b;
  q_i_a_plus_sm  =c;
}



class RREA
{

public:

  RREA(ZHJRANDOMv3*, int, int, int);

  void SetBeta(double);
  void PositionReset(void);
  void PopulationInitialization(void);
  bool PopulationInitialization( string&);  // read a population file

  void PopulationDynamics0(void);
  void PopulationDynamics(  string&);

  void Thermodynamics(int,  string&);

  void EraseBuffer( string&); /*write the remaining elements of the buffer
				     to file */
  void ReportPopulation( string&);


private:
  int ClauseDegree;        //=L
  int ClauseCavityDegree;  //=L-1
  int VariableDegree;       //=K
  int VariableCavityDegree; //K-1
  double Alpha;   /*clause density=clause number/variable number 
		    = VariableDegree/ClauseDegree */ 

  double Cvalue; //constant part of free energy density

  double Temperature; //temperature
  double Beta;       //inverse temperature = 1/Temperature
  double TanhBetaJ;  //=tanh(Beta*J), J=1

  valarray<DT> CavityProbabilities; //population array
  valarray<double> ProbabilityMultiple; //auxiliary container
  valarray<int>    ProbabilityPosition; //auxiliary container

  int PopulationSize;
  int PopulationBegin1;
  int PopulationBegin2;

  valarray<double> Buffer; //buffer for convenience of output
  int BuffCapacity;        //buffer capacity;
  int BuffLength;          //number of actually stored elements in baffer

  ZHJRANDOMv3* rdptr; //pointer to random number generator

  void UpdateCavityProbabilities0(const int&);
  void UpdateCavityProbabilities(const int&, double&);
  void UpdateCavityProbabilities(const int&,double&,double&,double&,double&);

};


int main(int argc, char ** argv)
{
  int rdseed =323456789;
  int rdseed2=197654321;
  cout <<"rdseed= ";  cout.flush();
  cin >>rdseed;
  cout <<"rdvalue equilibriation =";  cout.flush();
  cin >>rdseed2;
  ZHJRANDOMv3 rdgenerator(rdseed);
  for(int i=0; i<rdseed2; ++i) rdgenerator.rdflt();

  int ClauseDegree=2;
  int VariableDegree=4;
  int PopulationSize = 10000;
  RREA bp_zhj(&rdgenerator, VariableDegree, ClauseDegree, PopulationSize);

  double Beta;
  do {
    cout<<"Beta = ";
    cin>>Beta;
    if(Beta<0) continue;

    bp_zhj.SetBeta(Beta);

    string sigfilename;
    cout  <<"td output file "; cout.flush();
    cin  >>sigfilename;
    string popfilename;
    cout  <<"population output file ";  cout.flush();
    cin  >>popfilename;
    string valfilename;
    cout  <<"measurement output file "; cout.flush();
    cin  >>valfilename;

    int equil_time;
    cout<<"equilibrium time ";
    cin>>equil_time;
    int measure_time;
    cout<<"trajectory length ";
    cin>>measure_time;

    int itype=0;
    do {
      cout<<"New initial population (0-yes/1-from file/2-no)?";
      cin>>itype;
    }
    while(itype != 0 && itype != 1 && itype != 2);
    if(itype==0)
      bp_zhj.PopulationInitialization();
    else if(itype==1)
      {
	if(bp_zhj.PopulationInitialization(popfilename) == false)
	  {
	    cout<<"Population input failed.\n";
	    return -1;
	  }
      }
    bp_zhj.PositionReset();

    //    for(int iteration=1; iteration<=equil_time; ++iteration)
    //      bp_zhj.PopulationDynamics0();
    for(int iteration=1; iteration<=equil_time; ++iteration)
      bp_zhj.PopulationDynamics(sigfilename);
    bp_zhj.EraseBuffer(sigfilename);

    bp_zhj.Thermodynamics(measure_time, valfilename);

    bp_zhj.ReportPopulation(popfilename);
  }
  while(Beta>=0);

  return 1;
}



RREA::RREA
(ZHJRANDOMv3* rd, int v_degree, int c_degree, int pop_size)
{
  VariableDegree       = v_degree;
  VariableCavityDegree = VariableDegree-1;
  ClauseDegree         = c_degree;
  ClauseCavityDegree   = ClauseDegree-1;

  Alpha = 
    static_cast<double>(VariableDegree)/static_cast<double>(ClauseDegree);

  PopulationSize=pop_size;
  PopulationBegin1 = 0;
  PopulationBegin2 = PopulationSize;

  rdptr = rd;  

  try { CavityProbabilities.resize(2*PopulationSize); }
  catch(bad_alloc)
    { cerr<<"CavityProbabilities construction failed.\n";
      exit(-1);
    }

  try { ProbabilityMultiple.resize(ClauseDegree+1); }
  catch(bad_alloc) 
    { cerr<<"ProbabilityMultiple construction failed.\n";
      exit(-1);
    }

  try { ProbabilityPosition.resize(ClauseDegree); }
  catch(bad_alloc)
    {
      cerr<<"ProbabilityPosition construction failed.\n";
      exit(-1);
    }

  BuffCapacity=1000;
  try { Buffer.resize(BuffCapacity); }
  catch(bad_alloc)
    {
      cerr<<"Buffer construction failed.\n";
      exit(-1);
    }
  BuffLength=0;

  return ;
}



void RREA::SetBeta(double b)
{
  Beta = b;
  Temperature = 1.0e0/Beta;
  TanhBetaJ=tanh(Beta);
  BuffLength=0;
  Cvalue = -Alpha*(Beta+log(0.5e0*(1+exp(-2.0e0*Beta))));
  return ;
}



void RREA::PositionReset(void)
{
  PopulationBegin1=0;
  PopulationBegin2=PopulationSize;
  return ;
}



void RREA::PopulationInitialization(void)
{
  for(int position=0; position<PopulationSize; ++position)
    CavityProbabilities[position]=DT(rdptr->rdflt(), 1.0e0, 0.0e0);
  return ;
}



bool RREA::PopulationInitialization( string& filename)
{
  int pop;
  ifstream input(filename.c_str());
  if(!input.good())
    {
      cerr<<"population file probably non-existant.\n";
      return false;
    }

  input >> pop ;
  if(pop!=PopulationSize)
    {
      input.close();
      cerr<<"number of cavity probabilitiess incorrect.\n";
      return false;
    }

  double p_1, p_2, p_3;
  for(int position=0; position<PopulationSize; ++position)
    {
      input >> p_1 >> p_2 >> p_3;
      CavityProbabilities[position]=DT( p_1, p_2, p_3);
    }
  input.close();
  return true;
}



void RREA::ReportPopulation( string& filename)
{
  ofstream output(filename.c_str() );
  output << PopulationSize <<endl<<endl;
  for(int position=0; position<PopulationSize; ++position)
    output << CavityProbabilities[position].q_i_a_plus_mean << '\t'
	   << CavityProbabilities[position].q_i_a_plus_sp   << '\t'
	   << CavityProbabilities[position].q_i_a_plus_sm   << endl;
  output.close();
  return ;
}



/* belief propagation population dynamics */
void RREA::PopulationDynamics0(void)
{
  for(int position=0; position<PopulationSize; ++position)
    UpdateCavityProbabilities0(position+PopulationBegin2);

  if(PopulationBegin1==0)
    {
      PopulationBegin1=PopulationSize;
      PopulationBegin2=0;
    }
  else
    {
      PopulationBegin1=0;
      PopulationBegin2=PopulationSize;
    }
  return ;
}



/* survey propagation population dynamics at y=beta */
void RREA::PopulationDynamics(  string& filename )
{
  double mean_delta_pm=0.0e0; //difference between bp value
  double delta_pm;

  for(int position=0; position<PopulationSize; ++position)
    {
      UpdateCavityProbabilities(position+PopulationBegin2, delta_pm);
      mean_delta_pm    += delta_pm;
    }
  mean_delta_pm    /= PopulationSize;

  Buffer[BuffLength++ ]=mean_delta_pm;
  if(BuffLength==BuffCapacity)
    {
      ofstream output(filename.c_str(), ios_base::app);
      for(int i=0; i<BuffCapacity; ++i) output << Buffer[i] <<endl;
      output.close();
      BuffLength=0;
    }

  if(PopulationBegin1==0)
    {
      PopulationBegin1=PopulationSize;
      PopulationBegin2=0;
    }
  else
    {
      PopulationBegin1=0;
      PopulationBegin2=PopulationSize;
    }

  return ;
}





/* write to file the remaining elements in Buffer */
void RREA::EraseBuffer( string& filename)
{
  if(BuffLength>0)
    {
      ofstream output(filename.c_str(), ios_base::app);
      for(int i=0; i<BuffLength; ++i)
	output << Buffer[i] <<endl;
      output.close();
      BuffLength=0;
    }
  return ;
}





/* grand free energy density, mean free energy density, complexity */
void RREA::Thermodynamics(int iterations,  string& filename)
{
  int buff=1000;
  double gfe_density[buff];   //grand free energy density
  double mfe_density[buff];   //mean free energy density
  double Sigma[buff];         //complexity
  int blength=0;
  for(int repeat=0; repeat<iterations; ++repeat)
    {
      double bgi_density=0;
      double bga_density=0;
      double bfi_density=0;
      double bfa_density=0;

      double bgi, bga, bfi, bfa;

      for(int position=0; position<PopulationSize; ++position)
	{
	  UpdateCavityProbabilities(position+PopulationBegin2,bgi,bga,bfi,bfa);
	  bgi_density += bgi;
	  bga_density += bga;
	  bfi_density += bfi;
	  bfa_density += bfa;
	}
      bgi_density /= PopulationSize; 
      bga_density /= PopulationSize; 
      bfi_density /= PopulationSize; 
      bfa_density /= PopulationSize; 

      double val_bg=Cvalue+bgi_density-Alpha*(ClauseDegree-1)*bga_density;
      double val_bf=Cvalue+bfi_density-Alpha*(ClauseDegree-1)*bfa_density;

      Sigma[blength]=val_bf-val_bg;
      gfe_density[blength]=val_bg/Beta;
      mfe_density[blength]=val_bf/Beta;

      ++blength;
      if(blength==buff)
	{
	  ofstream output(filename.c_str(), ios_base::app);
	  for(int i=0; i<buff; ++i)
	    output << gfe_density[i]<<'\t'
		   << mfe_density[i]<<'\t'
		   << Sigma[i]<<endl;
	  output.close();
	  blength=0;
	}

      if(PopulationBegin1==0)
	{
	  PopulationBegin1=PopulationSize;
	  PopulationBegin2=0;
	}
      else
	{
	  PopulationBegin1=0;
	  PopulationBegin2=PopulationSize;
	}
    }

  if(blength>0)
    {
      ofstream output(filename.c_str(), ios_base::app);
      for(int i=0; i<blength; ++i)
	output << gfe_density[i]<<'\t'
	       << mfe_density[i]<<'\t'
	       << Sigma[i]<<endl;
      output.close();
      blength=0;
    }
  return ;
}

  




/* updating the mean probability $\overline{q}_{i->a}(s_i)$ */
void RREA::UpdateCavityProbabilities0
(const int& position)
{
  double a_mean   = 1.0e0; // weight for s_i=+1
  double b_mean   = 1.0e0; // weight for s_i=-1

  for(int b=0; b < VariableCavityDegree;  ++b)
    {
      int J_b=(rdptr->rdflt()<0.5? 1:-1);
      ProbabilityMultiple[0]=1.0e0; 
      // meaning: prod_{j in \partial b\i} (q_{j->b}^+ - q_{j->b}^-)

      for(int j=0; j < ClauseCavityDegree; ++j)
	ProbabilityMultiple[0] *= -1 +
	  2*CavityProbabilities[static_cast<int>(PopulationSize*rdptr->rdflt())
				+PopulationBegin1].q_i_a_plus_mean;
      double val2= J_b*ProbabilityMultiple[0]*TanhBetaJ;
      a_mean *= 1.0e0+val2;
      b_mean *= 1.0e0-val2;
    }
  CavityProbabilities[position].q_i_a_plus_mean = a_mean/(a_mean+b_mean);
  return ;
}



/* updating DT */
void RREA::UpdateCavityProbabilities
(const int& position, double& delta_q_ia_pm)
{
  double a_mean   = 1.0e0; //weight for s_i=+1
  double b_mean   = 1.0e0; //weight for s_i=-1

  double a_p_mean = 1.0e0; //weight for s_i=+1 given that s_i=+1
  double b_p_mean = 1.0e0; //weight for s_i=-1 given that s_i=+1

  double a_m_mean = 1.0e0; //weight for s_i=+1 given that s_i=-1
  double b_m_mean = 1.0e0; //weight for s_i=-1 given that s_i=-1

  int position_rand ;

  for(int b=0; b < VariableCavityDegree; ++b)
    {
      int J_b=(rdptr->rdflt()<0.5? 1:-1);
      ProbabilityMultiple[0]=1.0e0;

      for(int j=ClauseCavityDegree-1; j>0; --j)
	{
	  position_rand = static_cast<int>(PopulationSize*rdptr->rdflt())
	    +PopulationBegin1;
	  ProbabilityPosition[j]=position_rand;

	  ProbabilityMultiple[0] *=
	    2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

	  ProbabilityMultiple[j]=ProbabilityMultiple[0];
	  /* ProbabilityPosition[i]= 
	     prod_{j=i+1,ClauseCavityDegree} (q_{j->b}^+ - q_{j->b}^-) */
	}
      ProbabilityMultiple[ClauseCavityDegree]=1; //introduced for convenience
      position_rand =static_cast<int>( PopulationSize*rdptr->rdflt() )
	+PopulationBegin1;
      ProbabilityPosition[0] = position_rand;
      ProbabilityMultiple[0] *=
	2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

      double val2=J_b*ProbabilityMultiple[0]*TanhBetaJ;
      a_mean *= 1.0e0+val2;
      b_mean *= 1.0e0-val2;

      int sig_p= J_b; //sigma_i =  1
      int sig_m=-J_b; //sigma_i = -1

      double a_p=1.0e0; /* = prod_{j=1,ClauseCavityDegree}
			   (q_{j->b}^+ - q_{j->b}^-)  given sigma_i=+1 */
      double a_m=1.0e0; /* = prod_{j=1,ClauseCavityDegree}
			   (q_{j->b}^+ - q_{j->b}^-)  given sigma_i=-1 */
      for(int j=1; j <=ClauseCavityDegree; ++j)
	{
	  val2=ProbabilityMultiple[j-1]*TanhBetaJ;
	  double denominator_p=1.0e0+sig_p*val2;
	  double denominator_m=1.0e0+sig_m*val2;
	  position_rand=ProbabilityPosition[j-1];

	  double ratio=
	    (CavityProbabilities[position_rand].q_i_a_plus_mean *
	     (1+sig_p*ProbabilityMultiple[j]*TanhBetaJ))/denominator_p;
	  if(rdptr->rdflt()<ratio)
	    a_p *= 2*CavityProbabilities[position_rand].q_i_a_plus_sp-1;
	  else
	    {
	      sig_p *= -1;
	      a_p *= 2*CavityProbabilities[position_rand].q_i_a_plus_sm-1;
	    }

	  ratio=(CavityProbabilities[position_rand].q_i_a_plus_mean*
		 (1+sig_m*ProbabilityMultiple[j]*TanhBetaJ))/denominator_m;
	  if(rdptr->rdflt()<ratio)
	    a_m *= 2*CavityProbabilities[position_rand].q_i_a_plus_sp-1;
	  else
	    {
	      sig_m *= -1;
	      a_m *=2*CavityProbabilities[position_rand].q_i_a_plus_sm-1;
	    }
	}
      val2=J_b*a_p*TanhBetaJ;
      a_p_mean *= 1+val2;
      b_p_mean *= 1-val2;

      val2=J_b*a_m*TanhBetaJ;
      a_m_mean *= 1+val2;
      b_m_mean *= 1-val2;
    }
  CavityProbabilities[position]=
    DT(a_mean/(a_mean+b_mean), a_p_mean/(a_p_mean+b_p_mean),
       a_m_mean/(a_m_mean+b_m_mean));

  delta_q_ia_pm = 0.5e0*(abs(CavityProbabilities[position].q_i_a_plus_mean
			     -CavityProbabilities[position].q_i_a_plus_sp)+
			 abs(CavityProbabilities[position].q_i_a_plus_mean
			     -CavityProbabilities[position].q_i_a_plus_sm));
  return ;
}




/* updating population array,
   calculate grand free energy and free energy contributions */

void RREA::UpdateCavityProbabilities
(const int& position, double& bg_i, double& bg_a, double& bf_i, double& bf_a)
{
  double a_mean   = 1.0e0; //weight for s_i=+1
  double b_mean   = 1.0e0; //weight for s_i=-1

  double a_p_mean = 1.0e0; //weight for s_i=+1 given that s_i=+1
  double b_p_mean = 1.0e0; //weight for s_i=-1 given that s_i=+1

  double a_m_mean = 1.0e0; //weight for s_i=+1 given that s_i=-1
  double b_m_mean = 1.0e0; //weight for s_i=-1 given that s_i=-1

  int position_rand, sig_p, sig_m;
  double val2, ratio, a_p, a_m;

  //now updating one element of the population

  for(int b=0; b < VariableCavityDegree; ++b)
    {
      int J_b=(rdptr->rdflt()<0.5? 1:-1);
      ProbabilityMultiple[0]=1.0e0;

      for(int j=ClauseCavityDegree-1; j>0; --j)
	{
	  position_rand = static_cast<int>(PopulationSize*rdptr->rdflt())
	    +PopulationBegin1;
	  ProbabilityPosition[j]=position_rand;

	  ProbabilityMultiple[0] *=
	    2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

	  ProbabilityMultiple[j]=ProbabilityMultiple[0];
	  /* ProbabilityPosition[i]= 
	     prod_{j=i+1,ClauseCavityDegree} (q_{j->b}^+ - q_{j->b}^-) */
	}
      ProbabilityMultiple[ClauseCavityDegree]=1; //introduced for convenience
      position_rand =static_cast<int>( PopulationSize*rdptr->rdflt() )
	+PopulationBegin1;
      ProbabilityPosition[0] = position_rand;
      ProbabilityMultiple[0] *=
	2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

      val2=J_b*ProbabilityMultiple[0]*TanhBetaJ;
      a_mean *= 1.0e0+val2;
      b_mean *= 1.0e0-val2;

      sig_p= J_b; //sigma_i =  1
      sig_m=-J_b; //sigma_i = -1

      a_p=1.0e0; /* = prod_{j=1,ClauseCavityDegree}
		    (q_{j->b}^+ - q_{j->b}^-)  given sigma_i=+1 */
      a_m=1.0e0; /* = prod_{j=1,ClauseCavityDegree}
		    (q_{j->b}^+ - q_{j->b}^-)  given sigma_i=-1 */

      for(int j=1; j <=ClauseCavityDegree; ++j)
	{
	  val2=ProbabilityMultiple[j-1]*TanhBetaJ;
	  double denominator_p=1.0e0+sig_p*val2;
	  double denominator_m=1.0e0+sig_m*val2;
	  position_rand=ProbabilityPosition[j-1];

	  ratio=
	    (CavityProbabilities[position_rand].q_i_a_plus_mean *
	     (1+sig_p*ProbabilityMultiple[j]*TanhBetaJ))/denominator_p;
	  if(rdptr->rdflt()<ratio)
	    a_p *= 2*CavityProbabilities[position_rand].q_i_a_plus_sp-1;
	  else
	    {
	      sig_p *= -1;
	      a_p *= 2*CavityProbabilities[position_rand].q_i_a_plus_sm-1;
	    }

	  ratio=(CavityProbabilities[position_rand].q_i_a_plus_mean*
		 (1+sig_m*ProbabilityMultiple[j]*TanhBetaJ))/denominator_m;
	  if(rdptr->rdflt()<ratio)
	    a_m *= 2*CavityProbabilities[position_rand].q_i_a_plus_sp-1;
	  else
	    {
	      sig_m *= -1;
	      a_m *=2*CavityProbabilities[position_rand].q_i_a_plus_sm-1;
	    }
	}
      val2=J_b*a_p*TanhBetaJ;
      a_p_mean *= 1+val2;
      b_p_mean *= 1-val2;

      val2=J_b*a_m*TanhBetaJ;
      a_m_mean *= 1+val2;
      b_m_mean *= 1-val2;
    }
  CavityProbabilities[position]=
    DT(a_mean/(a_mean+b_mean), a_p_mean/(a_p_mean+b_p_mean),
       a_m_mean/(a_m_mean+b_m_mean));

  //now calculate beta*g_{i+partial i}

  int J_a=(rdptr->rdflt()<0.5? 1:-1);
  ProbabilityMultiple[0]=1.0e0;

  for(int j=ClauseCavityDegree-1; j>0; --j)
    {
      position_rand = static_cast<int>(PopulationSize*rdptr->rdflt())
	+PopulationBegin1;
      ProbabilityPosition[j]=position_rand;

      ProbabilityMultiple[0] *=
	2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

      ProbabilityMultiple[j]=ProbabilityMultiple[0];
      /* ProbabilityPosition[i]= 
	 prod_{j=i+1,ClauseCavityDegree} (q_{j->a}^+ - q_{j->a}^-) */
    }
  ProbabilityMultiple[ClauseCavityDegree]=1;
  position_rand =static_cast<int>( PopulationSize*rdptr->rdflt() )
    +PopulationBegin1;
  ProbabilityPosition[0] = position_rand;
  ProbabilityMultiple[0] *=
    2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

  val2=J_a*ProbabilityMultiple[0]*TanhBetaJ;
  a_mean *= 1.0e0+val2;
  b_mean *= 1.0e0-val2;

  bg_i=-log(a_mean+b_mean); //constant part not yet considered

  //now calculate beta* <f_{i+partial i}>

  ratio=a_mean/(a_mean+b_mean);
  if(rdptr->rdflt()<ratio) //sigma_i=1
    {
      sig_p= J_a; //sigma_i =  1
      a_p=1.0e0; /* = prod_{j=1,ClauseCavityDegree}
		    (q_{j->a}^+ - q_{j->a}^-)  given sigma_i=+1 */
      for(int j=1; j <=ClauseCavityDegree; ++j)
	{
	  val2=ProbabilityMultiple[j-1]*TanhBetaJ;
	  double denominator_p=1.0e0+sig_p*val2;
	  position_rand=ProbabilityPosition[j-1];

	  ratio=
	    (CavityProbabilities[position_rand].q_i_a_plus_mean *
	     (1+sig_p*ProbabilityMultiple[j]*TanhBetaJ))/denominator_p;
	  if(rdptr->rdflt()<ratio)
	    a_p *= 2*CavityProbabilities[position_rand].q_i_a_plus_sp-1;
	  else
	    {
	      sig_p *= -1;
	      a_p *= 2*CavityProbabilities[position_rand].q_i_a_plus_sm-1;
	    }
	}
      val2=J_a*a_p*TanhBetaJ;
      a_p_mean *= 1+val2;
      b_p_mean *= 1-val2;

      bf_i = -log(a_p_mean+b_p_mean); //constant part not yet considered
    }
  else //sigma_i=-1
    {
      sig_m=-J_a; //sigma_i = -1
      a_m=1.0e0; /* = prod_{j=1,ClauseCavityDegree}
		    (q_{j->a}^+ - q_{j->a}^-)  given sigma_i=-1 */
      for(int j=1; j <=ClauseCavityDegree; ++j)
	{
	  val2=ProbabilityMultiple[j-1]*TanhBetaJ;
	  double denominator_m=1.0e0+sig_m*val2;
	  position_rand=ProbabilityPosition[j-1];

	  ratio=(CavityProbabilities[position_rand].q_i_a_plus_mean*
		 (1+sig_m*ProbabilityMultiple[j]*TanhBetaJ))/denominator_m;
	  if(rdptr->rdflt()<ratio)
	    a_m *= 2*CavityProbabilities[position_rand].q_i_a_plus_sp-1;
	  else
	    {
	      sig_m *= -1;
	      a_m *=2*CavityProbabilities[position_rand].q_i_a_plus_sm-1;
	    }
	}
      val2=J_a*a_m*TanhBetaJ;
      a_m_mean *= 1+val2;
      b_m_mean *= 1-val2;

      bf_i = -log(a_m_mean+b_m_mean); //constant part not yet considered
    }

  //now calculate beta*g_a

  J_a = (rdptr->rdflt()<0.5? 1:-1);
  ProbabilityMultiple[0]=1.0e0;

  for(int j=ClauseDegree-1; j>0; --j)
    {
      position_rand = static_cast<int>(PopulationSize*rdptr->rdflt())
	+PopulationBegin1;
      ProbabilityPosition[j]=position_rand;

      ProbabilityMultiple[0] *=
	2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

      ProbabilityMultiple[j]=ProbabilityMultiple[0];
      /* ProbabilityPosition[i]= 
	 prod_{j=i+1,ClauseDegree} (q_{j->a}^+ - q_{j->a}^-) */
    }
  ProbabilityMultiple[ClauseDegree]=1;
  position_rand =static_cast<int>( PopulationSize*rdptr->rdflt() )
    +PopulationBegin1;
  ProbabilityPosition[0] = position_rand;
  ProbabilityMultiple[0] *=
    2*CavityProbabilities[position_rand].q_i_a_plus_mean-1;

  bg_a = -log(1.0e0+J_a*TanhBetaJ*ProbabilityMultiple[0]);
  //constant part not yet considered

  //now calculate beta <f_a>

  int sig_pm= J_a;
  double a_pm=1.0e0; /* =prod_{j=1,ClauseDegree} (q_{j->a}^+ - q_{j->a}^-) */
  for(int j=1; j <=ClauseDegree; ++j)
    {
      val2=ProbabilityMultiple[j-1]*TanhBetaJ;
      double denominator_pm=1.0e0+sig_pm*val2;
      position_rand=ProbabilityPosition[j-1];

      ratio=
	(CavityProbabilities[position_rand].q_i_a_plus_mean *
	 (1+sig_pm*ProbabilityMultiple[j]*TanhBetaJ))/denominator_pm;
      if(rdptr->rdflt()<ratio)
	a_pm *= 2*CavityProbabilities[position_rand].q_i_a_plus_sp-1;
      else
	{
	  sig_pm *= -1;
	  a_pm *= 2*CavityProbabilities[position_rand].q_i_a_plus_sm-1;
	}
    }
  bf_a = -log(1.0e0+J_a*TanhBetaJ*a_pm); //constant part not yet considered

  return ;
}