//
// BeliefPropagationSAT.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 08.08.2012
//   
// DESCRIPTION:
// Entropic belief-propagation for a general SAT formula. The program
// start from the all-zero initial condition.
//
// compile this program using
// c++ -O3 -o SATbp.exe BeliefPropagationSAT.cpp
//
// Programmer:
// Haijun Zhou
// Institute of Theoretical Physics, Chinese Academy of Sciences,
// Beijing 100190, China
// Webpage http://www.itp.ac.cn/~zhouhj
// Email zhouhj@itp.ac.cn
//

#include <cmath>
#include <exception>
#include <fstream>
#include <iostream>
#include <set>
#include <string>
#include <valarray>
#include "zhjrandom.h"  //random number generator

using namespace std;


struct uarraystruct //data structure used in updating clause<->variable messages
{
  double  minv;
  double  expv;
  double  resv;
  uarraystruct(double, double, double);
};
uarraystruct::uarraystruct(double m=0, double e=0, double r=0)
{
  minv=m;
  expv=e;
  resv=r;
}


struct etastruct //message on each variable node
{
  double  p; // log(weight) of being positive=\sum_{b\in\partial i^{+}} u_{b->i}
  double  m; // log(weight of being negative =\sum_{c\in\partial i^{-}} u_{c->i}
  int     punsat; //number of warnings from positive edges
  int     munsat; //number of warnings from negative edges
  etastruct(double, double, int, int);
};
etastruct::etastruct(double a=0, double b=0, int c=0, int d=0)
{
  p     =a;
  m     =b;
  punsat=c;
  munsat=d;
}


struct actionctovstruct //clause to variable message
{
  int var; //index of variable node
  unsigned char bar; //edge coupling bar=1 (positive), =0 (negative)
  double u; 
  /* u: clause to variable message.
     u=log[1-\prod_{j\in \partial b\backslash i} (1-J_b^j m_{j->b})/2]
  */
  double expu; //=exp(u)
  actionctovstruct(void);
};
actionctovstruct::actionctovstruct(void)
{
  var =0;
  bar =0;
  u   =0;
  expu=1.0e0;
}


struct clausestruct //clause
{
  struct actionctovstruct *actionctov; //beginning address of output messages
  int degree; //number of edges attached
  clausestruct(void);
};
clausestruct::clausestruct(void)
{
  actionctov=0;
  degree=0;
}


struct actionvtocstruct //variable to clause message
{
  struct clausestruct * clause; //address of destination clause
  int position; //position of starting variable in neighborhood list of clause
  actionvtocstruct(void);
};
actionvtocstruct::actionvtocstruct(void)
{
  clause=0;
  position=0;
}


struct variablestruct //variable
{
  int degree; //number of edges attached
  struct actionvtocstruct *actionvtoc; //beginning address of output messages
  struct etastruct *etaptr; //address of message vector
  variablestruct(void);
};
variablestruct::variablestruct(void)
{
  actionvtoc=0;
  etaptr=0;
  degree=0;
}


class BPforSAT
{
public:
  BPforSAT(ZHJRANDOMv3* );
  void Initial_bpmessages(void);
  bool Readformula( string&, int);
  bool Sequential_converge(double, int);
  bool Compute_entropy(double&);
  void Report_magnetization( string&);

private:
  int Variable_number; //number of variable nodes
  int Clause_number; //number of clause nodes
  int Max_clause_degree; //maximal clause degree
  double Epsilon; //convergence condition
  int Iterations; //number of Belief-propagation iterations
  int Num_trivial_sat; //number of trivially satisfied clauses
  ZHJRANDOMv3* rdptr; //pointer to random number generator
  valarray<variablestruct> Variables; //variable array
  valarray<clausestruct> Clauses; //clause array
  valarray<actionvtocstruct> Alledgesvtoc; //variable to clause messages
  valarray<actionctovstruct> Alledgesctov; //clause to variable messages
  valarray<etastruct> Etaset; //variable probability array
  valarray<int> Permutation; //random sequential updating, order of updating
  valarray<uarraystruct> Uarray; //internal array to ensure no data overflow or underflow

  void Compute_eta(void);
  void Compute_eta(struct variablestruct* );
  double Update_u(struct clausestruct*);
  double Iterate(void);
};


int main(int argc, char ** argv)
{
  cout
    <<"SATbp by Haijun Zhou (last modified 08.08.2012)."
    <<endl;
  cout.flush();

  int rdseed = 271828183;
  cout
    <<"random number seed (positive number) = ";
  cout.flush();
  cin
    >>rdseed;

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

  BPforSAT sp_zhj(&rdgenerator);

  cout <<"Name of SAT formula file = ";
  cout.flush();
  string inputfilename;
  cin >>inputfilename;
  cout<<"How many clauses to be read in ? M=";
  cout.flush();
  int Clause_number;
  cin>>Clause_number;

  if(sp_zhj.Readformula(inputfilename, Clause_number) ==false)
    return -1;

  sp_zhj.Initial_bpmessages();

  bool tocontinue=true;
  double epsilon;
  int iterations;
  double entropy_density;
  while(tocontinue)
    {
      cout <<"BP convergence criterion epsilon = ";
      cout.flush();
      cin >>epsilon;
      cout <<"Number of iterations = ";
      cout.flush();
      cin >>iterations;

      if(sp_zhj.Sequential_converge(epsilon, iterations))
	tocontinue=false;
    }

  if(sp_zhj.Compute_entropy(entropy_density)==false)
    {
      cout <<"Contradiction (formula may not have a solution).\n";
      cout.flush();
      return -1;
    }

  string outputfilename;
  cout <<"Entropy output file ";
  cout.flush();
  cin >>outputfilename;
  ofstream outputfile(outputfilename.c_str());
  outputfile <<inputfilename
	     <<"  has entropy density  "
	     <<entropy_density
	     <<endl;
  outputfile.close();

  string magnet_name;
  cout <<"Magnetization file name = ";
  cout.flush();
  cin >>magnet_name;
  sp_zhj.Report_magnetization(magnet_name);

  return 0;
}


BPforSAT::BPforSAT( ZHJRANDOMv3* rd )
{
  rdptr=rd;
  Epsilon=1.0e-10;
  Iterations=200;
}


//Read the sat formula
bool BPforSAT::Readformula( string& filename, int clval)
{
  Clause_number=clval;

  ifstream source(filename.c_str());
  if(!source.good() )
    {
      cerr<<"Error in SAT formula input (file possibly non-existent).\n";
      cerr.flush();
      return false;
    }

  string abcd1, abcd2;
  int Mval;
  source
    >> abcd1
    >> abcd2
    >> Variable_number
    >> Mval;
  if(Mval<Clause_number)
    {
      source.close();
      return false;
    }

  try { Variables.resize(Variable_number+1); } catch(bad_alloc) { return false; }
  try { Clauses.resize(Clause_number+1); } catch(bad_alloc) { return false; }
  try { Etaset.resize(Variable_number); } catch(bad_alloc) { return false; }

  //first pass for counting
  bool goodfile=true;
  int edge_index=0;
  int var,degree;
  bool satisfied;
  Max_clause_degree=0;
  Num_trivial_sat=0; //number of trivially satisfied clauses
  set<int> neighbors;
  for(int m=1; m<=Clause_number & goodfile; ++m)
    {
      degree=0;
      satisfied=false;
      source>>var;
      neighbors.clear();
      while(var!=0 && goodfile)
	{
	  if( (var>Variable_number) || (var<(-Variable_number)) )
	    //variable index too large
	    {
	      cerr<<"Too many variables.\n";
	      cerr.flush();
	      goodfile=false;
	    }
	  else
	    {
	      if(neighbors.find(var)==neighbors.end())
		{
		  neighbors.insert(var);
		  if(neighbors.find(-var)!=neighbors.end())
		    satisfied=true;
		}
	    }
	  source>>var;
	}
      if(!goodfile) continue;

      if(!satisfied)
	{
	  for(set<int>::const_iterator sci=neighbors.begin();
	      sci!=neighbors.end(); ++sci)
	    {
	      var= *sci;
	      if(var>0) ++(Variables[ var ].degree);
	      else ++(Variables[ - var ].degree);
	      ++degree;
	      ++edge_index;
	    }
	  Clauses[m-Num_trivial_sat].degree=degree;
	  if(Max_clause_degree<degree) Max_clause_degree=degree;
	}
      else ++Num_trivial_sat;
    }
  source.close();

  if(!goodfile) return false;

  try { Alledgesvtoc.resize(edge_index); } catch(bad_alloc) { exit(-1);}
  try { Alledgesctov.resize(edge_index); } catch(bad_alloc) { exit(-1);}
  struct actionvtocstruct *clause_ptr=&Alledgesvtoc[0];
  struct etastruct *etaptr=&Etaset[0];
  int number_isolated_variables=0;
  for(var=1; var<=Variable_number; ++var)
    {
      if(Variables[var].degree>0)
	{
	  Variables[var].actionvtoc=clause_ptr;
	  clause_ptr += Variables[var].degree;
	  Variables[var].etaptr=etaptr;
	  ++etaptr;
	  Variables[var].degree=0;
	}
      else ++number_isolated_variables;
    }
  struct actionctovstruct *variable_ptr=&Alledgesctov[0];
  for(int m=1; m<=Clause_number-Num_trivial_sat; ++m)
    {
      Clauses[m].actionctov=variable_ptr;
      variable_ptr += Clauses[m].degree;
    }

  //second pass to do the actual reading
  source.open(filename.c_str());
  source >> abcd1
	 >> abcd2
	 >> Variable_number
	 >> Mval;
  Num_trivial_sat=0;
  for(int m=1; m<=Clause_number; ++m)
    {
      degree=0;
      satisfied=false;
      source>>var;
      neighbors.clear();
      while(var!=0)
	{
	  if(neighbors.find(var)==neighbors.end())
	    {
	      neighbors.insert(var);
	      if(neighbors.find(-var)!=neighbors.end())
		satisfied=true;
	    }
	  source>>var;
	}
      if(satisfied) ++Num_trivial_sat;
      else
	{
	  for(set<int>::const_iterator sci=neighbors.begin(); sci!=neighbors.end(); ++sci)
	    {
	      var = *sci;
	      Clauses[m-Num_trivial_sat].actionctov[degree].bar=( var>0? 1:0);
	      if(var<0) var *= -1;
	      Clauses[m-Num_trivial_sat].actionctov[degree].var=var;
	      Variables[var].actionvtoc[Variables[var].degree].clause
		= &Clauses[m-Num_trivial_sat];
	      Variables[var].actionvtoc[ Variables[var].degree++ ].position
		=degree++;
	    }
	}
    }
  source.close();
  
  try { Permutation.resize(Clause_number-Num_trivial_sat) ; }
  catch(bad_alloc) { return false; }
  try { Uarray.resize(Max_clause_degree); } catch(bad_alloc) { return false; }

  cout
    <<"Clause number "
    <<Clause_number
    <<"  Variable number "
    <<Variable_number
    <<"  Isolated variable "
    <<number_isolated_variables
    <<"  Number of trivially satisfied clauses "
    <<Num_trivial_sat
    <<endl;
  cout.flush();

  return true;
}


//Message initialization
void BPforSAT::Initial_bpmessages(void)
{
  struct clausestruct *clause=&Clauses[1];
  struct actionctovstruct *actionctov;
  for(int cl=1; cl<=Clause_number-Num_trivial_sat; ++cl, ++clause)
    {
      actionctov=clause->actionctov;
      for(int dd=0; dd<clause->degree; ++dd,++actionctov)
	{
	  actionctov->u= -0.1e0*rdptr->rdflt();
	  actionctov->expu=exp(actionctov->u);
	}
    }
  return ;
}


//update the message on a variable node
void BPforSAT::Compute_eta(struct variablestruct *varptr)
{
  varptr->etaptr->p=0;
  varptr->etaptr->m=0;
  varptr->etaptr->punsat=0;
  varptr->etaptr->munsat=0;

  struct actionvtocstruct *actionvtoc=varptr->actionvtoc;
  struct actionctovstruct *edgeptr;
  for(int d=0; d<varptr->degree; ++d, ++actionvtoc)
    {
      edgeptr=(actionvtoc->clause->actionctov)+(actionvtoc->position);
      if(edgeptr->u<=0) //normal message
	{
	  if(edgeptr->bar) //positive coupling
	    varptr->etaptr->p += edgeptr->u; 
	  else
	    varptr->etaptr->m += edgeptr->u;
	}
      else //warning
	{
	  if(edgeptr->bar)  ++(varptr->etaptr->punsat);
	  else  ++(varptr->etaptr->munsat);
	}
    }
  return ;
}


//update the messages on all the variable nodes
void BPforSAT::Compute_eta(void)
{
  for(int var=1; var<=Variable_number; ++var)
    if(Variables[var].degree>0)
      Compute_eta(&Variables[var]);
  return ;
}


//updates all eta's and etavtoc's in clause cl
double BPforSAT::Update_u(struct clausestruct *clause)
{
  double max_eps=0;

  for(int d=0; d<clause->degree; ++d)
    Uarray[d].minv=0;

  /*
    Uarray[i] is an interval data used for precise computation of clause
    to variable message
    u_{a->i}=log[1-\prod_{j\in \partial a\backslash i}(1-J_a^j m_{j->a})/2 ].
    m_{j->a}=(Nja-Pja)/(Nja+Pja), with
    log(Nja)=\sum\limits_{b\in \partial j^{-} \backslash a} exp(u_{b->j})
    log(Pja)=\sum\limits_{b\in \partial j^{+} \backslash a} exp(u_{b->j})
  */

  struct actionctovstruct *actionctov=clause->actionctov;
  struct etastruct *etaptr;
  for(int dd=0; dd<clause->degree; ++dd, ++actionctov)
    {
      int I_ja=actionctov->bar ? 1: -1;
      etaptr=Variables[actionctov->var].etaptr;
      int por=I_ja*(etaptr->punsat-etaptr->munsat);
      double eta_ja=etaptr->m-etaptr->p;
      if(actionctov->u<=0) //no warning
	{
	  if(por==0) //variable unfrozen
	    {
	      eta_ja += I_ja * actionctov->u;
	      if(eta_ja>=0) eta_ja += 2.0; //a trick: eta_ja>=2 means positive eta_ja
	      else eta_ja -= 2.0; //trick: eta_ja<=-2 means negative eta_ja
	    }
	  else //variable frozen
	    eta_ja = por>0? I_ja : -I_ja; //trick: eta_ja=1 (positive frozen), =-1 (negative frozen)
	}
      else // actionctov->u=1 (warning)
	{
	  if(por==1)
	    {
	      if(eta_ja>=0) eta_ja += 2.0;
	      else eta_ja -= 2.0;
	    }
	  else eta_ja=por>1 ? I_ja : -I_ja;
	}

      if(eta_ja<0)
	{
	  eta_ja *= -1;
	  I_ja *= -1;
	} //eta_ja=2*(cavity field) on variable j, it is positive with respect to I_ja

      if(eta_ja>1.5) //variable j not frozen in absence of a, (eta_ja>=2)
	{
	  double val=1.0e0+exp(2.0e0-eta_ja); //m_{j->a}=[1-exp(2-eta_ja)]/[1+exp(2-eta_ja)]
	  for(int bb=0; bb<clause->degree; ++bb)
	    if(bb!=dd)
	      {
		double minv=Uarray[bb].minv;
		if(minv>1.5) //minv>=2
		  {
		    if(I_ja>0) Uarray[bb].expv += eta_ja-2.0e0; 
		    /* Uarray[i].expv= -\sum_{j\neq i} log[(1-I_aj)/2+e^{2-eta_ja} (1+I_aj)/2] */

		    /* \prod_{j\neq i}(1+e^{2-eta_ja})=1+exp(-Uarray[i].minv)*Uarray[i].resv */
		    if(minv<=eta_ja)
		      Uarray[bb].resv=val*Uarray[bb].resv+exp(minv-eta_ja);
		    else
		      {
			Uarray[bb].resv=1.0+exp(eta_ja-minv)*val*Uarray[bb].resv;
			Uarray[bb].minv=eta_ja;
		      }
		  }
		else if(minv < 0.5) // not yet initialized (minv=0)
		  Uarray[bb]=uarraystruct(eta_ja,(I_ja>0? eta_ja-2.0 : 0),1);
		//min,expv,resv
		//else  min=1, no need to proceed (because clause satisfied)
	      }
	}
      else // eta_ja=1, frozen (the direction is defined as positive)
	{
	  if(I_ja>0)
	    {
	      for(int bb=0; bb<clause->degree; ++bb)
		if(bb!=dd)
		  Uarray[bb].minv=1; //clause satisfied u_{b->i}=0
	    }
	  //else I_ja<0, clause can not satisfied by variable
	}
    }

  actionctov=clause->actionctov;
  for(int dd=0; dd<clause->degree; ++dd, ++actionctov)
    {
      double new_u=0; //=log(1.0e0-\prod_{j\neq i} [(1-I_aj m_{j->a})/2])
      etaptr=Variables[actionctov->var].etaptr;
      if(Uarray[dd].minv>=1.5e0) //minv>=2
	{
	  if(Uarray[dd].expv>0)
	    new_u=log(1.0e0-exp(-Uarray[dd].expv)/
		      (1.0e0+exp(2.0e0-Uarray[dd].minv)*Uarray[dd].resv));
	  else new_u=2.0-Uarray[dd].minv
	    +log(Uarray[dd].resv/
		 (1.0e0+exp(2.0e0-Uarray[dd].minv)*Uarray[dd].resv));

	  if(actionctov->u<=0) //normal message
	    {
	      if(actionctov->bar) etaptr->p += new_u-actionctov->u;
	      else                etaptr->m += new_u-actionctov->u;
	    }
	  else //warning
	    {
	      if(actionctov->bar)
		{
		  --(etaptr->punsat);
		  etaptr->p += new_u;
		}
	      else
		{
		  --(etaptr->munsat);
		  etaptr->m += new_u;
		}
	    }
	}
      else if(Uarray[dd].minv>=0.5) //=1, clause satisfied
	{
	  new_u=0;
	  if(actionctov->u<=0)
	    {
	      if(actionctov->bar)
		etaptr->p -= actionctov->u;
	      else
		etaptr->m -= actionctov->u;
	    }
	  else
	    {
	      if(actionctov->bar) --(etaptr->punsat);
	      else                --(etaptr->munsat);
	    }
	}
      else // Uarray[edge_index].minv=0, warning
	{
	  new_u=1.0;
	  if(actionctov->u <=0)
	    {
	      if(actionctov->bar)
		{
		  ++(etaptr->punsat);
		  etaptr->p -= actionctov->u;
		}
	      else
		{
		  ++(etaptr->munsat);
		  etaptr->m -= actionctov->u;
		}
	    }
	}
      double new_expu=(new_u<=0?exp(new_u):0);
      double eps=abs(new_expu-actionctov->expu);
      if(max_eps<eps) max_eps=eps;
      actionctov->u = new_u;
      actionctov->expu=new_expu;
    }

  return max_eps;
}


//update u's of clauses in a random permutation order
double BPforSAT::Iterate(void)
{
  double max_eps=0;
  for(int quant=Clause_number-Num_trivial_sat; quant>0; --quant)
    {
      int i=static_cast<int>(quant*rdptr->rdflt());
      int cl=Permutation[i];
      if(i<(quant-1))
	{
	  Permutation[i]=Permutation[quant-1];
	  Permutation[quant-1]=cl;
	}
      double eps=Update_u(&Clauses[cl]);
      if(eps>max_eps)
	max_eps=eps;
    }

  return max_eps;
}


//Construct a fixed point for the belief-propagation equation
bool BPforSAT::Sequential_converge(double newepsilon, int newiterations)
{
  Epsilon=newepsilon;
  Iterations=newiterations;

  double max_eps=0;
  int iter=1;
  int quant=0;

  Compute_eta();

  for(int cl=1; cl<=Clause_number-Num_trivial_sat; ++cl) Permutation[quant++]=cl;
  do
    {
      max_eps=Iterate();
      Compute_eta();

      if((iter%50)==0)
	{
	  cerr<<max_eps<<' ';
	  cerr.flush();
	}
      ++iter;
    }
  while (max_eps>Epsilon && iter<=Iterations);

  if(max_eps<=Epsilon)
    {
      cerr<<'\t'<<max_eps<<" :-)\n";
      return true;
    }
  else
    {
      cerr<<'\t'<<max_eps<<" :-(\n";
      return false;
    }
}


//calculate the Bethe-Peierls entropy
bool BPforSAT::Compute_entropy(double& entropy_density)
{
  double entropy_clauses=0;
  double entropy_variables=0;

  bool satisfied=true;

  struct clausestruct *clause=&Clauses[1];
  struct actionctovstruct *actionctov;
  struct etastruct *etaptr;
  for(int cl=1; cl<=Clause_number-Num_trivial_sat && satisfied; ++cl, ++clause)
    if(clause->degree>1)
      {
	double minv=0;
	double expv=0;
	double resv=0;

	actionctov=clause->actionctov;
	for(int dd=0; dd<clause->degree; ++dd,++actionctov)
	  {
	    int I_ja=actionctov->bar? 1:-1;
	    etaptr=Variables[actionctov->var].etaptr;
	    int por=I_ja*(etaptr->punsat-etaptr->munsat);
	    double eta_ja=etaptr->m-etaptr->p;
	    if(actionctov->u<=0)
	      {
		if(por==0)
		  {
		    eta_ja += I_ja*actionctov->u;
		    if(eta_ja>=0) eta_ja += 2.0;
		    else eta_ja -= 2.0;
		  }
		else eta_ja = por>0? I_ja : -I_ja;
	      }
	    else // actionctov->u=1 (warning)
	      {
		if(por==1)
		  {
		    if(eta_ja>=0) eta_ja += 2.0;
		    else eta_ja -= 2.0;
		  }
		else eta_ja=por>1 ? I_ja : -I_ja;
	      }

	    if(eta_ja<0)
	      {
		eta_ja *= -1;
		I_ja   *= -1;
	      }

	    if(eta_ja>1.5)
	      {
		double val=1.0+exp(2.0e0-eta_ja);
		if(minv>1.5) // not warning
		  {
		    if(I_ja>0) expv += eta_ja-2.0;
		    if(minv<=eta_ja) resv=val*resv+exp(minv-eta_ja);
		    else
		      {
			resv=1.0+exp(eta_ja-minv)*val*resv;
			minv=eta_ja;
		      }
		  }
		else if(minv < 0.5) // not yet initialized
		  {
		    minv=eta_ja;
		    expv=(I_ja>0? eta_ja-2.0:0);
		    resv=1;
		  }
		// else  min=1,  no need to proceed
	      }
	    else 	// eta_ja=1
	      {
		if(I_ja>0) minv=1;
	      }
	  }

	if(minv>=1.5)
	  {
	    if(expv>0)
	      entropy_clauses +=
		(clause->degree-1)*log(1.0e0-exp(-expv)/
				       (1.0e0+exp(2.0e0-minv)*resv));
	    else
	      entropy_clauses += (clause->degree-1)*
		(2.0-minv+log(resv/(1.0e0+exp(2.0e0-minv)*resv)));
	  }
	else if(minv==0)
	  satisfied=false;
      }
  if(!satisfied)
    return satisfied;

  for(int var=1; var<=Variable_number && satisfied; ++var)
    {
      if(Variables[var].degree==0)
	{
	  entropy_variables += log(2.0e0);
	  continue;
	}
      Compute_eta(&Variables[var]);
      etaptr=Variables[var].etaptr;
      if((etaptr->punsat==0) && (etaptr->munsat==0))
	{
	  if(etaptr->p >= etaptr->m)
	    entropy_variables += etaptr->p +log(1.0e0+exp(etaptr->m-etaptr->p));
	  else
	    entropy_variables += etaptr->m +log(1.0e0+exp(etaptr->p-etaptr->m));
	}
      else if(etaptr->punsat==0)
	entropy_variables += etaptr->p;
      else if(etaptr->munsat==0)
	entropy_variables += etaptr->m;
      else
	satisfied=false;
    }
  if(!satisfied)
    return satisfied;

  entropy_density=(entropy_variables - entropy_clauses)/Variable_number;
  cerr <<Variable_number
       <<" variables,  entropy density "
       <<entropy_density<<endl;

  return satisfied;
}


//Report the average spin value of each variable node
void BPforSAT::Report_magnetization( string& filename)
{
  ofstream output(filename.c_str() );
  struct etastruct *etaptr;
  double qEA=0;
  double meanval;
 for(int var=1; var<=Variable_number; ++var)
    {
      if(Variables[var].degree==0) output<<var<<'\t'<<0<<endl;
      else
	{
	  Compute_eta(&Variables[var]);
	  etaptr=Variables[var].etaptr;
	  int hardfield=etaptr->punsat-etaptr->munsat;
	  if(hardfield==0)
	    {
	      double value=etaptr->m-etaptr->p;
	      if(value>=0.0e0)
		{
		  value=exp(-value);
		  meanval=(1.0e0-value)/(1.0e0+value);
		}
	      else
		{
		  value=exp(value);
		  meanval=(value-1.0e0)/(value+1.0e0);
		}
	    }
	  else if(hardfield>0)
	    meanval=1;
	  else
	    meanval=-1;

	  output<<var<<'\t'<<meanval<<endl;
	  qEA += meanval*meanval;
	}
    }
  output.close();
  cout<<"qEA="<<qEA/Variable_number<<endl;

  return ;
}