//
// BeliefPropDecimationSAT.cpp
//
// This program complements Section 6.3 of the book
// "Spin Glass and Message Passing" by Hai-Jun Zhou
// (Science Press, Beijing, 2015).
// Program last modified on 20.08.2012
//
// DESCRIPTION:
// Entropic belief-propagation-guided decimation for a general SAT formula.
//
// compile this program using
// c++ -O3 -o SATbpd.exe BeliefPropDecimationSAT.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 <set>
#include <string>
#include <valarray>
#include "zhjrandom.h"

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
  int active_degree;
  clausestruct(void);
};
clausestruct::clausestruct(void)
{
  actionctov=0;
  degree=0;
  active_degree=0;
}


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


struct variablestruct //variable
{
  int degree; //number of edges attached
  int spin;
  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;
  spin=0;
}



class BPDforSAT
{
public:
  BPDforSAT(ZHJRANDOMv3*);
  int Readformula              ( string&);
  int Sequential_converge      (float, int, float&);
  bool Compute_entropy          (float&);
  bool Decimation               (int);
  bool Decimation               (void);
  int Report_number_freespin   (void);
  void Report_partialsolution  ( string&);
  void Restart                 (void);

private:
  int                  Variable_number;
  int                  Clause_number;
  int                  Number_edges;
  int                  Max_clause_degree;
  int                  Number_freespin;
  int                  Iterations;
  int                  Is_frozen;
  int                  Number_candidate_variables;
  int                  Number_isolated_variables;
  double               Epsilon;
  float                Maximal_field;

  ZHJRANDOMv3*         rdptr;

  valarray<variablestruct>               Variables;
  valarray<clausestruct>                 Clauses;
  valarray<actionvtocstruct>             Alledgesvtoc;
  valarray<actionctovstruct>             Alledgesctov;
  valarray<etastruct>                    Etaset;
  valarray<int>                          Nclausetype;
  valarray<int>                          Permutation;
  valarray<int>                          Set_candidate_variables;
  valarray<uarraystruct>                 Uarray;

  set<int>                               Special_clauses;
  set<int>                               Isolated_variables;

  void   Compute_eta  (void);
  void   Compute_eta  (struct variablestruct* );
  void   Fix_isolated_variables(void);
  int    Compute_eta  (const int&, double&);
  int    Fix          (int, int);
  int    Simplify     (int);
  float  Update_u     (struct clausestruct*);
  float  Iterate      (void);
};


int main(int argc, char ** argv)
{
  cout<<"Belief-propagation-guided-decimation for a SAT formula\n"
      <<"Haijun Zhou (last modified 20.08.2012)\n\n";
  cout.flush();

  unsigned int rdseed;  cout<<"rdseed= "; 
  cout.flush(); 
  cin>>rdseed;
  if(rdseed>=999999998 || rdsee<0)  rdseed=323456789;

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

  BPDforSAT sp_zhj(&rdgenerator);

  string inputfilename;
  cout<<"cnf formula file --:  ";  cout.flush();
  cin>>inputfilename;
  sp_zhj.Readformula(inputfilename);

  string solutionfilename;
  cout<<"output file for solution --:  ";  cout.flush();
  cin>>solutionfilename;

  string trajectoryfilename;
  cout<<"output file for iteration logs --:  ";  cout.flush();
  cin>>trajectoryfilename;

  int    iterations       = 20;
  int    total_trial      = 100;
  float  epsilon          = 5e-6;
  float  entropy_density  = -1;
  float  epsval;

  cout<<"changing default for iterations (y/n)? ";  cout.flush();
  char sig;
  cin>>sig;
  if(sig=='y' || sig=='Y')
    {
      cout<<"reset epsilon (old value was"<<epsilon<<") ";  cout.flush();
      cin>>epsilon;
      if(epsilon>=1 || epsilon<=0)
	epsilon=5e-6;

      cout<<"reset number of iterations (old value was "<<iterations<<") ";   cout.flush();
      cin>>iterations;
      if(iterations<=0)
	iterations=20;
      else if(iterations>1000)
	iterations=1000;

      cout<<"reset maximal number of trials (old value was"<<total_trial<<") "; cout.flush();
      cin>>total_trial;
      if(total_trial<=0)
	total_trial=10;
      else if(total_trial>100)
	total_trial=100;
    }

  // a long iteration to drive to fixed point
  sp_zhj.Sequential_converge(epsilon, iterations*200, epsval);

  int trial_index=0;
  int failure=1;
  ofstream output(trajectoryfilename.c_str());
  while((failure==1 || entropy_density<0) && trial_index<total_trial)
    {
      failure=0;
      entropy_density=1;

      int repeat_time   =0;
      while(sp_zhj.Report_number_freespin()>0 && failure==0 && entropy_density>=0)
	{
	  sp_zhj.Sequential_converge(epsilon, iterations, epsval);
	  failure=1-sp_zhj.Compute_entropy(entropy_density);
	  output<<trial_index<<' '<<repeat_time<<'\t'<<epsval<<'\t'<<entropy_density<<endl;
	  output.flush();
	  cout<<trial_index<<" "<<repeat_time<<" "<<entropy_density<<" "
	      <<sp_zhj.Report_number_freespin()<<" ";
	  cout.flush();
	  if(failure==0 && entropy_density>=0)
	    {
	      failure=1-sp_zhj.Decimation(); /* each step only fix one */
	      cout<<sp_zhj.Report_number_freespin()<<endl;
	      cout.flush();
	    }
          ++repeat_time;
	  //	  if(repeat_time % 1000 ==0)
	  //	    sp_zhj.Report_partialsolution(solutionfilename.c_str());
	}
      if(failure==1 || entropy_density<0)
	{
	  cout<<"Trial --"<<trial_index<<" failed. Restarting!\n";
	  cout.flush();
	  sp_zhj.Restart();
	  sp_zhj.Sequential_converge(epsilon, iterations*200, epsval);
	}
      ++trial_index;
    }
  output.close();

  if(sp_zhj.Report_number_freespin()==0)
    {
      cout<<"Solution found!\n";
      cout.flush();
      sp_zhj.Report_partialsolution(solutionfilename.c_str());
      return 1;
    }
  else
    {
      cout<<"All trials failed.\n";
      cout.flush();
      return 0;
    }
}


BPDforSAT::BPDforSAT(ZHJRANDOMv3* rd)
{
  rdptr                          = rd;
  Epsilon                        = 5.0e-6;
  Iterations                     = 50;

  Maximal_field                  = 0;
  Is_frozen                      = 0;
  Number_candidate_variables     = 0;
}


int BPDforSAT::Report_number_freespin(void)
{
  return Number_freespin;
}


int BPDforSAT::Readformula( string& filename)
{
  ifstream source(filename.c_str());
  if(!source.good() )
    {
      cerr<<"Error in cnf formula input (file possibly non-existent).\n";
      cerr.flush();
      exit(-1);
    }
  string abcd1, abcd2;
  source>>abcd1 >> abcd2 >> Variable_number>> Clause_number;

  try { Variables.resize(Variable_number+1); } catch(bad_alloc) { exit(-1); }
  try { Clauses.resize(Clause_number+1); } catch(bad_alloc) { exit(-1); }
  try { Etaset.resize(Variable_number); } catch(bad_alloc) { exit(-1); }
  try { Set_candidate_variables.resize(Variable_number); } catch(bad_alloc) { exit(-1); }

  //first pass for counting
  int edge_index=0;
  Max_clause_degree=0;
  Special_clauses.clear();
  for(int m=1; m<=Clause_number; ++m)
    {
      int var;
      int size=0;
      bool satisfied=false;
      source>>var;
      set<int> neighbors;
      while(var!=0)
	{
	  if(abs(var)>Variable_number)
	    {
	      cerr<<"Too many variables.\n";
	      cerr.flush();
	      exit(-1);
	    }
	  if(neighbors.find(var)==neighbors.end())
	    {
	      neighbors.insert(var);
	      ++(Variables[abs(var)].degree);
	      ++size;
	      ++edge_index;
	      if(neighbors.find(-var)!=neighbors.end())
		satisfied=true;
	    }
	  source>>var;
	}
      Clauses[m].degree=size;
      Clauses[m].active_degree=(satisfied?0:size);
      if(satisfied)
	Special_clauses.insert(m);
      if(Max_clause_degree<size)
	Max_clause_degree=size;
    }
  source.close();

  Number_edges=edge_index;

  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 *eta_ptr=&Etaset[0];
  Number_isolated_variables=0;
  for(int var=1; var<=Variable_number; ++var)
    {
      Variables[var].spin=0;
      if(Variables[var].degree)
	{
	  Variables[var].actionvtoc=clause_ptr;
	  clause_ptr += Variables[var].degree;
	  Variables[var].etaptr=eta_ptr;
	  ++eta_ptr;
	  Variables[var].degree=0;
	}
      else
	{
	  Isolated_variables.insert(var);
	  ++Number_isolated_variables;
	}
    }
  struct actionctovstruct *variable_ptr=&Alledgesctov[0];
  for(int m=1; m<=Clause_number; ++m)
    if(Clauses[m].active_degree)
      {
	Clauses[m].actionctov=variable_ptr;
	variable_ptr += Clauses[m].degree;
      }
  
  //second pass to do the actual reading
  try { Nclausetype.resize(Max_clause_degree+1); } catch(bad_alloc) { exit(-1); }
  source.open(filename.c_str());
  source>>abcd1 >> abcd2 >> Variable_number >> Clause_number;
  for(int m=1; m<=Clause_number; ++m)
    {
      int size=0;
      int var;
      source>>var;
      set<int> neighbors;
      while(var!=0)
	{
	  if(neighbors.find(var)==neighbors.end())
	    {
	      neighbors.insert(var);
	      Clauses[m].actionctov[size].bar=(var>0?1:0);
	      var=abs(var);
	      Clauses[m].actionctov[size].var=var;
	      Variables[var].actionvtoc[Variables[var].degree].clause = &Clauses[m];
	      Variables[var].actionvtoc[Variables[var].degree++].position=size++;
	    }
	  source>>var;
	}
      ++Nclausetype[Clauses[m].active_degree];
    }
  source.close();
  Number_freespin=Variable_number;
  
  try { Permutation.resize(Clause_number) ; } catch(bad_alloc) { exit(-1); }
  try { Uarray.resize(Max_clause_degree); } catch(bad_alloc) { exit(-1); }
  
  cout<<"Clause number "<<Clause_number
      <<"  Variable number "<<Variable_number
      <<"  Isolated variable "<<Number_isolated_variables<<endl;
  cout.flush();

  if(Number_isolated_variables>0)
    {
      Fix_isolated_variables();
      Number_freespin -= Number_isolated_variables;
    }
  
  return 0;
}


//asign a spin value completely at random for each isolated variable
void BPDforSAT::Fix_isolated_variables
(void)
{
  if(Number_isolated_variables>0)
    {
      for(set<int>::const_iterator sci=Isolated_variables.begin();
	  sci!=Isolated_variables.end(); ++sci)
	Variables[*sci ].spin=(rdptr->rdflt()<=0.5? 1: -1);

      Number_freespin -= Number_isolated_variables;
    }
  return ;
}


//completely discard previous spin asignment and re-initialize the messages
void BPDforSAT::Restart
(void)
{
  for(int v=1; v<=Variable_number; ++v)
    Variables[v].spin=0;
  Number_freespin=Variable_number;
  Fix_isolated_variables();

  for(int d=0; d<=Max_clause_degree; ++d)
    Nclausetype[d]=0;

  for(int cl=1; cl<=Clause_number; ++cl)
    {
      Clauses[cl].active_degree
	=(Special_clauses.find(cl)==Special_clauses.end()? Clauses[cl].degree : 0);
      ++Nclausetype[Clauses[cl].active_degree];
    }

  for(int edge=0; edge<Number_edges; ++edge)
    {
      Alledgesctov[edge].u=0;
      Alledgesctov[edge].expu=1;
    }

  return ;
}


//update the message on a variable node
void BPDforSAT::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)
    if(actionvtoc->clause->active_degree)
      {
	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 message on a variable node and calculate its entropy contribution
int BPDforSAT::Compute_eta
(const int& var, double& entropy_variables)
{
  struct variablestruct *varptr=&Variables[var];
  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)
    if(actionvtoc->clause->active_degree)
      {
	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);
	  }
      }

  if((varptr->etaptr->punsat==0) && (varptr->etaptr->munsat==0))
    {
      double value=varptr->etaptr->m - varptr->etaptr->p;
      if(value<=0)
	entropy_variables += varptr->etaptr->p +log(1.0e0+exp(value));
      else
	entropy_variables += varptr->etaptr->m +log(1.0e0+exp(-value));

      if(Is_frozen==0)
	{
	  if(abs(value)>=Maximal_field)
	    {
	      if(abs(value)>Maximal_field)
		{
		  Maximal_field=abs(value);
		  Number_candidate_variables=0;
		}
	      Set_candidate_variables[Number_candidate_variables++ ]=var;
	    }
	}
      return 1;
    }
  else if(etaptr->punsat==0)
    {
      entropy_variables += varptr->etaptr->p;
      if(Is_frozen==0)
	{
	  Is_frozen=1;
	  Maximal_field=varptr->etaptr->munsat;
	  Set_candidate_variables[0]=var;
	  Number_candidate_variables=1;
	}
      else if(varptr->etaptr->munsat >= Maximal_field)
	{
	  if(varptr->etaptr->munsat>Maximal_field)
	    {
	      Maximal_field=varptr->etaptr->munsat;
	      Number_candidate_variables=0;
	    }
	  Set_candidate_variables[Number_candidate_variables++ ]=var;
	}
      return 1;
    }
  else if(varptr->etaptr->munsat==0)
    {
      entropy_variables += varptr->etaptr->m;
      if(Is_frozen==0)
	{
	  Is_frozen=1;
	  Maximal_field=varptr->etaptr->punsat;
	  Set_candidate_variables[0]=var;
	  Number_candidate_variables=1;
	}
      else if(varptr->etaptr->punsat>=Maximal_field)
	{
	  if(varptr->etaptr->punsat>Maximal_field)
	    {
	      Maximal_field=varptr->etaptr->punsat;
	      Number_candidate_variables=0;
	    }
	  Set_candidate_variables[Number_candidate_variables++ ]=var;
	}
      return 1;
    }
  else
    return 0;   // contradiction
}


//update all the variable nodes
void BPDforSAT::Compute_eta
(void)
{
  for(int var=1; var<=Variable_number; ++var)
    if(Variables[var].spin==0)
      Compute_eta(&Variables[var]);

  return ;
}


//updates all eta's and etavtoc's in clause cl
double BPDforSAT::Update_u
(struct clausestruct *clause)
{
  /*
    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})
  */
  double max_eps=0;

  struct actionctovstruct *actionctov=clause->actionctov;
  int edge_index=0;
  for(int dd=0; dd<clause->degree; ++dd,++actionctov)
    if(Variables[actionctov->var].spin==0)
      Uarray[edge_index++].minv=0;

  edge_index=0;
  actionctov=clause->actionctov;
  struct etastruct *etaptr;
  for(int dd=0; dd<clause->degree; ++dd,++actionctov)
    if(Variables[actionctov->var].spin==0)
      {
	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->active_degree; ++bb)
	      if(bb!=edge_index)
		{
		  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 has been satisfied)
		}
	  }
	else // eta_ja=1, frozen (the direction is defined as positive)
	  {
	    if(I_ja>0)
	      {
		for(int bb=0; bb<clause->active_degree; ++bb)
		  if(bb!=edge_index)
		    Uarray[bb].minv=1; //clause satisfied u_{b->i}=0
	      }
	    //else I_ja<0, clause can not satisfied by variable
	  }
	++edge_index;
      }

  edge_index=0;
  actionctov=clause->actionctov;
  for(int dd=0; dd<clause->degree; ++dd,++actionctov)
    if(Variables[actionctov->var].spin==0)
      {
	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[edge_index].minv>=1.5e0) //minv>=2
	  {
	    if(Uarray[edge_index].expv>0)
	      new_u=log(1.0e0-exp(-Uarray[edge_index].expv)/
			(1.0e0+exp(2.0e0-Uarray[edge_index].minv)*Uarray[edge_index].resv));
	    else new_u=2.0-Uarray[edge_index].minv
	      +log(Uarray[edge_index].resv/
		   (1.0e0+exp(2.0e0-Uarray[edge_index].minv)*Uarray[edge_index].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[edge_index].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;
	++edge_index;
      }

  return max_eps;
}


//update u's of clauses in a random permutation order
double BPDforSAT::Iterate
(void)
{
  double max_eps=0;
  for(int quant=Clause_number-Nclausetype[0]; 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 BPDforSAT::Sequential_converge
(double newepsilon, int newiterations, double& max_eps)
{
  Epsilon=newepsilon;
  Iterations=newiterations;

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

  Compute_eta();

  for(int cl=1; cl<=Clause_number; ++cl)
    if(Clauses[cl].active_degree)
      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 && satisfied; ++cl, ++clause)
    if(clause->active_degree>1)
      {
	double minv=0;
	double expv=0;
	double resv=0;

	actionctov=clause->actionctov;
	for(int dd=0; dd<clause->degree; ++dd,++actionctov)
	  if(Variables[actionctov->var].spin==0)
	    {
	      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->active_degree-1)*log(1.0e0-exp(-expv)/(1.0e0+exp(2.0e0-minv)*resv));
	    else
	      entropy_clauses += (clause->active_degree-1)*
		(2.0-minv+log(resv/(1.0e0+exp(2.0e0-minv)*resv)));
	  }
	else if(minv==0)
	  satisfied=false;
      }
  if(!satisfied)
    return satisfied;

  // calculate the entropy contribution from variables
  Is_frozen=0;// 0 unfrozen; 1 being frozen
  Maximal_field=0;
  Number_candidate_variables=0;
  for(int var=1; var<=Variable_number && satisfied; ++var)
    if(Variables[var].spin==0)
      satisfied=Compute_eta(var, entropy_variables);
  if(!satisfied)
    return satisfied;

  entropy_density=(entropy_variables - entropy_clauses)/Number_freespin;

  return true;
}


//performing decimation
bool BPDforSAT::Decimation
(void)
{
  int v_index=
    (Number_candidate_variables>1? static_cast<int>(Number_candidate_variables*rdptr->rdflt()) : 0);
  int var=Set_candidate_variables[v_index];
  if(Is_frozen==0)
    {
      float rdval=rdptr->rdflt();
      float eta_i=(Variables[var].etaptr)->m - (Variables[var].etaptr)->p;
      if(eta_i>=0)
	return Fix(var,1);
      else
	return Fix(var,-1);
    }
  else if((Variables[var].etaptr)->munsat==0)
    return Fix(var,1);
  else if((Variables[var].etaptr)->punsat==0)
    return Fix(var,-1);
  else 
    return false;
}

  
//fix var to value spin and possibly simplify the resulting formula
bool BPDforSAT::Fix(int var, int spin)
{
  --Number_freespin;
  if(Variables[var].spin!=0)     //something must be wrong!
    return false;
  Variables[var].spin=spin;
  return Simplify(var);
}


//simplify the SAT formula after fixing var
bool BPDforSAT::Simplify(int var)
{
  struct actionvtocstruct *actionvtoc=Variables[var].actionvtoc;
  for(int cl=0; cl<Variables[var].degree; ++cl, ++actionvtoc)
    {
      struct clausestruct *clause=actionvtoc->clause;
      if(clause->active_degree==0) //clause has been satisfied
	continue;
      int position=actionvtoc->position;
      --Nclausetype[clause->active_degree];

      //check if var renders SAT the clause
      if(clause->actionctov[position].bar==(Variables[var].spin==1))
	{
	  ++Nclausetype[0];
	  clause->active_degree=0;
	  for(int i=0; i<clause->degree; ++i)
	    {
	      int var2=clause->actionctov[i].var;
	      if(Variables[var2].spin==0)
		Compute_eta(&Variables[var2]);
	    }
	  continue;
	}

      //otherwise, check for further simplifications
      //type 0, contradiction?:
      ++Nclausetype[--(clause->active_degree)];
      if(clause->active_degree==0)
	{
	  cerr<<"contradition in simplification.\n";
	  cerr.flush();
	  //	  ofstream contrafile("contradiction.tmp.cnf");
	  //	  writeformula(contrafile);
	  //	  exit(-1);
	  return false;
	}

      //no contradiction
      //type 1: unit clause propagation
      if(clause->active_degree==1)
	{
	  //find the unfixed literal
	  int aux=0,i;
	  for(i=0; i<clause->degree; ++i)
	    {
	      if(Variables[aux=clause->actionctov[i].var].spin==0)
		break;
	    }
	  if(i==clause->degree) //a clause could be unit-clause-fixed by two different paths
	    continue;
	  return Fix(aux,clause->actionctov[i].bar?1:-1);
	}
    }
  return true;
}


//Report the partial spin asignment
void BPDforSAT::Report_partialsolution( string& outputfile)
{
  ofstream output(outputfile.c_str());
  output<<Variable_number-Number_freespin<<endl<<endl; //total number of assigned variables

  for(int v=1; v<=Variable_number; ++v)
    if(Variables[v].spin!=0)
      output<< v*Variables[v].spin <<endl;
  output.close();

  return ;
}