/*
  DSbpdV03.cpp
  
  DESCRIPTION:
  Belief-propagation guided decimation (BPD) as a solver for the minimal
  dominating set (MDS) problem on an undirected simple graph. The program is
  applicable on a single graph instance.
  
  Each vertex has four states (unoccupied and unobserved, 0; occupied, 1;
  unoccupied but observed and having unobserved neighbors, 2; unoccupied but
  observed and having no unobserved neighbor, 3). If a vertex is occupied, 
  there is a unit energy cost. On the other hand, if a vertex is unoccupied, 
  then at least one of its nearest neighbors must be occupied.
  
  The belief propagation iteration is only carried out within the subnetwork of
  active nodes (whose states are either 0 or 2). Only the edges between two 0
  nodes or between one 0 node and one 2 node are considered. All the other 
  edges are not considered in the belief propagation iteration.
  
  A generalized leaf-removal (GLR) process is performed before the belief
  propagation iteration. This GLR process simplifies the network.
  
  For the theoretical details behind this algorithm please consult
  Jin-Hua Zhao, Yusupjan Habibulla, Hai-Jun Zhou,
  "Statistical mechanics of the minimum dominating set problem",
  Journal of Statistical Physics 159: 1154--1174 (2015),
  DOI: 10.1007/s10955-015-1220-2

  LOG:
  02.10.2014: last modified (the problem of possible underflow for nodes with
  extremely large degree is solved).
  30.09.2014:copied from DSbpdV02b.cpp to DSbpdV03.cpp. In DSbpdV03.cpp, the
  candidate nodes to be occupied must be an active node (either unobserved or
  having an unobserved neighbor). GLR process added.
  30.09.2014: copied from DSbpdV02bMPI.cpp to DSbpdV02b.cpp and modified (for
  single graph instances).
  23.09.2014: last modified.
  DSbpdV02bMPI.cpp differs from DSbpdV02MPI.cpp only in that data analysis is
  included in the program.
  23.09.2014 copied from DSbpdV02MPI.cpp to DSbpdV02bMPI.cpp.
  last modified 16.07.2014.
  copied from DSbpdV01.cpp to DSbpdV02MPI.cpp on 16.07.2014.
  copied from DFVSbpdV03MPI.cpp into DSbpdV01.cpp on 12.07.2014.
  
  PROGRAMMER:
  Hai-Jun Zhou
  Institute of Theoretical Physics, Chinese Academy of Sciences
  Zhong-Guan-Cun East Road 55, Beijing 100190
  zhouhj@itp.ac.cn
*/


#include <cmath>
#include <exception>
#include <fstream>
#include <iostream>
#include <set>
#include <string>
#include <valarray>
//#include "/Users/zhouhj/Programs/zhjrandom.h"
#include "zhjrandom.h"
using namespace std;


/*---               random real uniformly distributed in [0,1)            ---*/
double u01prn(ZHJRANDOMv3 *rd)
{
  return rd->rdflt();
}


struct OrderedIntPair
{
  int first;
  int second;
  OrderedIntPair(void);
  OrderedIntPair(int, int);
};

OrderedIntPair::OrderedIntPair(void)
{
  first=0;
  second=0;
  return;
}

OrderedIntPair::OrderedIntPair(int a, int b)
{
  if(a<=b) {    first=a;    second=b;  }
  else     {    first=b;    second=a;  }
  return;
}

bool operator<(OrderedIntPair a, OrderedIntPair b)
{
  if(a.first < b.first)      return true;
  else if(a.first > b.first) return false;
  else {                                                /* a.first = b.first */
    if(a.second<b.second)    return true;
    else                     return false;
  }
}

bool operator==(OrderedIntPair a, OrderedIntPair b)
{
  return (a.first==b.first) && (a.second==b.second);
}


/* message passed during belief propagation iteration.  notice that 
   q_{(j,i)->i}^{(c_j=1,c_i=1)} = q_{(j,i)->i}^{(c_j=1,c_i=0)},
   namely q_11=q_10 */
struct message
{
  double q_00;                                   //q_{(j,i)->i}^{(c_j=0,c_i=0)}
  double q_01;                                   //q_{(j,i)->i}^{(c_j=0,c_i=1)}
  double q_10;                                   //q_{(j,i)->i}^{(c_j=1,c_i=0)}
  message(void);
  message(double, double, double);
};

message::message(void)
{
  q_00=0;
  q_01=0;
  q_10=0;
  return;
}

message::message(double a, double b, double c)
{
  q_00=a;
  q_01=b;
  q_10=c;
  return;
}


struct outmessage
{
  struct message* m_ptr;                     //pointer to message storage place
  struct vertexstruct* v_ptr;           //pointer to nearest neighboring vertex
  outmessage(void);
};

outmessage::outmessage(void)
{
  m_ptr=0;
  v_ptr=0;
  return;
}

struct vertexstruct
{  
  int index;                                //index of vertex, positive integer
  int degree;                                                 //number of edges
  int state; /* 0: unoccupied & unobserved (active); 
		1: occupied (inactive); 
		2: unoccup but observed and have unobserved neighbor (active);
		3: unoccup but observed and have no unobserved neighbor (inact)
		In belief-propagation, only the edges between two unobserved
		nodes (whose states are 0) or between an unobserved (0) node 
		and an active observed (2) node are considered. All the other
		edges are not considered.
	     */
  
  struct message *im_ptr;      //start position of input message from neighbors
  struct outmessage *om_ptr;    //start position of output message to neighbors
  vertexstruct(void);
};

vertexstruct::vertexstruct(void)
{
  index=0;
  degree=0;
  state=0;
  im_ptr=0;
  om_ptr=0;
  return;
}


class DomSET
{
public:
  DomSET(ZHJRANDOMv3*);                                           //constructor
  bool Graph( string&, int);               //read graph connection pattern
  bool Percolation(void);        //perform a generalized leaf-removal initially
  void SetBeta(double);                               //set inverse temperature
  void SetDampingFactor(double);                           //set damping factor
  void Initialization(void);                            //initialize population
  bool BeliefPropagation(double, int);                    //population updating
  void Thermodynamics(void);                         //thermodynamic quantities
  bool Fix(void);                           //fix some vertex as being occupied
  bool CheckDS(int&, int&,  string&);                       //check result
  
private:
  int VertexNumber;                                  //total number of vertices
  int ActiveVertexNumber;                     //total number of active vertices
  int NumUnobserved;                      //total number of unobserved vertices
  int EdgeNumber;                                       //total number of edges
  int MaxDegree;                           //maximal vertex degree in the graph
  double Beta;                                            //inverse temperature
  double WeightOccupy;                                            //=exp(-Beta)
  double Rnew;                        //message updating, weight of new message
  double Rold;           //message updating, weight of old message, Rnew+Rold=1
  int MaxFixNumber;             //maximal number of fixed vertices in one round
  int MinRange;                //minimal range of occupy_prob of fixed vertices
  int MaxRange;                //maximal range of occupy_prob of fixed vertices
  float FixFraction;  /* fraction of active vertices to be fixed to be occupied
			 in each round of the decimation process */
  ZHJRANDOMv3 * PRNG;                                 //random number generator
  
  valarray<vertexstruct> Vertex;
  valarray<message> InMessage;
  valarray<outmessage> MessageAddress;
  
  valarray<double> Weight00;      //auxiliary array needed for message updating
  valarray<double> Weightx0;      //auxiliary array needed for message updating
  valarray<double> Weightx1;      //auxiliary array needed for message updating
  valarray<bool> IsActive;        //auxiliary array needed for message updating
  valarray<int> Permutation;         //array used in random sequential updating
  valarray<int> CandidateVertices;     //list of candidate vertices to be fixed
  valarray<int> CandidateSize;          /* number of candidate vertices in each
					   occupy_prob range */
  void UpdateMessage(struct vertexstruct *, double&);
};


int main(int argc, char ** argv)
{
  int rd_seed=38443855;                                   //random seed
  int VertexNumber =100000;                               //vertex number
  int EdgeNumber   =500000;                               //edge number
  string graphname="ERn100kM1.g100";   //graph input file
  string dsname="myresult.dat";        //dominating set output file

  ZHJRANDOMv3 PRNG( rd_seed );
  for(int i=0; i<2300059; ++i) PRNG.rdflt();
  
  DomSET system(&PRNG);
  
  bool succeed=system.Graph(graphname, EdgeNumber);
  if(succeed==false) {
    cerr<<"Some problem in graph reading.\n";
    return -1;
  }
  
  bool hascore=system.Percolation();
  if(hascore) {                                               //there is a core
    system.Initialization();
    double Beta=10.0e0;
    system.SetBeta(Beta);
    double DampingFactor=0.95e0;
    system.SetDampingFactor(DampingFactor);
    float epsilon = 1.0e-5;
    int iterations = 200;
    system.BeliefPropagation(epsilon, iterations);
    epsilon = 1.0e-5;
    iterations = 10;
    bool needcontinue=true;
    while(needcontinue) {
      system.BeliefPropagation(epsilon, iterations);
      system.Thermodynamics();
      needcontinue=system.Fix();
    }
  }
  int num_occup, num_obser;
  if(system.CheckDS(num_occup, num_obser, dsname)==false) {
    cerr<<"BPD failed, not a DS solution.\n";
  }
  return 1;
}


/*---                    constructor of DomSET cluster                    ---*/
DomSET::DomSET(ZHJRANDOMv3* rd)
{
  PRNG=rd;                                            //random number generator
  FixFraction=0.01;
  return;
}


/*-                            Read graph from file                        - */
bool DomSET::Graph( string& gname,                       //graph file name
		   int enumber)               //number of edges to be read into
{
  ifstream graphf(gname.c_str());
  if( !graphf.good() ) {
    cerr<<"Graph probably non-existant.\n";
    return false;      }
  graphf >> VertexNumber >> EdgeNumber;
  if(EdgeNumber<enumber) {
    cerr<<"No so many edges in the graph.\n";
    graphf.close();
    return false;        }
  EdgeNumber=enumber;              //only the first enumber edges are read into
  try { Vertex.resize(VertexNumber+1); }
  catch(bad_alloc) {
    cerr<<"Vertex construction failed.\n";
    graphf.close();
    return false;  }

  /* (1) Read graph file for the 1st time for checking */
  bool succeed=true;
  set<OrderedIntPair> EdgeSet;
  for(int eindex=0; eindex<EdgeNumber && succeed; ++eindex) {
    int v1, v2;
    graphf >> v1 >> v2;
    if(v1==v2 || v1==0 || v1>VertexNumber || v2==0 || v2>VertexNumber) {
      cerr<<"Graph incorrect at line " << eindex+1 <<endl;
      succeed=false;                                                   }
    else if(EdgeSet.find(OrderedIntPair(v1,v2) ) != EdgeSet.end() ) {
      cerr<<"Multiple edges.\n";
      succeed=false;                                                }
    else                                      {
      EdgeSet.insert(OrderedIntPair(v1,v2));
      ++(Vertex[v1].degree);
      ++(Vertex[v2].degree);                  }
  }
  graphf.close();
  if(succeed==false)  return false;
  EdgeSet.clear();
  try { InMessage.resize( 2*EdgeNumber ); }
  catch(bad_alloc) { 
    cerr<<"InMessage construction failed.\n";
    return false;  }
  try { MessageAddress.resize( 2*EdgeNumber ); }
  catch(bad_alloc) {
    cerr<<"MessageAddress construction failed.\n";
    return false;  }
  int position=0;
  MaxDegree=0;
  struct vertexstruct *v_ptr=&Vertex[1];
  for(int v=1; v<=VertexNumber; ++v) {
    v_ptr->index=v;
    v_ptr->state=0;
    v_ptr->im_ptr = &InMessage[position];
    v_ptr->om_ptr = &MessageAddress[position];
    position += v_ptr->degree;
    if(v_ptr->degree > MaxDegree) MaxDegree = v_ptr->degree;
    v_ptr->degree=0;
    ++v_ptr;
  }
  
  /* (2) Read graph file again for construction. */
  graphf.open(gname.c_str() );
  graphf >> VertexNumber >> EdgeNumber;
  EdgeNumber=enumber;
  struct outmessage *om_ptr=&MessageAddress[0];
  for(int eindex=0; eindex<EdgeNumber; ++eindex) {
    int v1,v2;
    graphf >> v1 >> v2;
    om_ptr=Vertex[v1].om_ptr + Vertex[v1].degree;
    om_ptr->m_ptr = Vertex[v2].im_ptr + Vertex[v2].degree;
    om_ptr->v_ptr = &Vertex[v2];
    om_ptr=Vertex[v2].om_ptr + Vertex[v2].degree;
    om_ptr->m_ptr = Vertex[v1].im_ptr + Vertex[v1].degree;
    om_ptr->v_ptr = &Vertex[v1];
    ++(Vertex[v1].degree);
    ++(Vertex[v2].degree);
  }
  graphf.close();
  cout<< VertexNumber <<'\t'
      << EdgeNumber   <<'\t'
      << (2.0e0*EdgeNumber)/(1.0e0*VertexNumber)
      <<'\t'<<"MaxDegree="<<MaxDegree<<endl;
  try { Weight00.resize(MaxDegree); }
  catch(bad_alloc) {
    cerr<<"Weight00 construction failed.\n";
    return false;  }
  try { Weightx0.resize(MaxDegree); }
  catch(bad_alloc) {
    cerr<<"Weightx0 construction failed.\n";
    return false;  }
  try { Weightx1.resize(MaxDegree); }
  catch(bad_alloc) {
    cerr<<"Weightx1 construction failed.\n";
    return false;  }
  try { IsActive.resize(MaxDegree); }
  catch(bad_alloc) {
    cerr<<"IsActive construction failed.\n";
    return false;  }
  try { Permutation.resize(VertexNumber); }
  catch(bad_alloc) {
    cerr<<"Permutation construction failed.\n";
    return false;  }
  NumUnobserved=VertexNumber;
  return true;
}


/*  Initial GLR process, return true if there is a core */
bool DomSET::Percolation(void)
{
  valarray<int> ActiveDegree;
  try { ActiveDegree.resize(VertexNumber+1); }
  catch(bad_alloc) {
    cerr<<"Failure of constructing ActiveDegree.\n";
    return false;  }
  struct vertexstruct *v_ptr =&Vertex[1];
  struct vertexstruct *v2_ptr;
  struct outmessage *om_ptr;
  struct outmessage *om2_ptr;
  set<int> Leafvertices;
  set<int> Newlyoccupied;
  set<int> Focusobserved;
  for(int v=1; v<=VertexNumber; ++v) {
    ActiveDegree[v]=v_ptr->degree;
    if(v_ptr->degree==0) {
      v_ptr->state=1;                                    //occupy isolated node
      --NumUnobserved;
    }
    else if(v_ptr->degree==1) Leafvertices.insert(v);
    ++v_ptr;
  }
  while( !Leafvertices.empty() ) {
    /* (1) treat leaf vertices */
    Newlyoccupied.clear();
    do {
      v_ptr=&Vertex[ *(Leafvertices.begin() )];
      Leafvertices.erase(v_ptr->index);
      if(v_ptr->state==0) {                            //leaf node not observed
	om_ptr = v_ptr->om_ptr;
	int degree=v_ptr->degree;
	bool found=false;
	for(int d=0; d<degree && !found; ++d) {
	  v2_ptr=om_ptr->v_ptr;
	  if( v2_ptr->state % 2 ==0) {
	    found=true;
	    /* effectively deleting the edges associated with a newly occupied
	       (directly observed) vertices. */
	    if(v2_ptr->state==0) {
	      --NumUnobserved;
	      om2_ptr=v2_ptr->om_ptr;
	      int degree2=v2_ptr->degree;
	      for(int d2=0; d2<degree2; ++d2) {
		if( om2_ptr->v_ptr->state % 2 ==0 )
		  --(ActiveDegree[ om2_ptr->v_ptr->index]);
		++om2_ptr;
	      }
	    }
	    else {                                           // v2_ptr->state=2
	      om2_ptr=v2_ptr->om_ptr;
	      int degree2=v2_ptr->degree;
	      for(int d2=0; d2<degree2; ++d2) {
		if( om2_ptr->v_ptr->state==0) 
		  --(ActiveDegree[om2_ptr->v_ptr->index]);
		++om2_ptr;
	      }
	    }
	    /* occupy the vertex */
	    v2_ptr->state=1;
	    ActiveDegree[v2_ptr->index]=0;
	    Newlyoccupied.insert(v2_ptr->index);
	  }
	  ++om_ptr;
	}
      }
    } while( !Leafvertices.empty() );
    /* (2) treat newly occupied vertices */
    Focusobserved.clear();
    while( !Newlyoccupied.empty() ) {
      v_ptr=&Vertex[ *(Newlyoccupied.begin() )];
      Newlyoccupied.erase(v_ptr->index);
      om_ptr=v_ptr->om_ptr;
      int degree=v_ptr->degree;
      for(int d=0; d<degree; ++d) {
	v2_ptr=om_ptr->v_ptr;
	if(v2_ptr->state==0) {
	  v2_ptr->state=2;
	  --NumUnobserved;
	  Focusobserved.insert(v2_ptr->index);
	  //delete the edge between two observed vertices
	  om2_ptr=v2_ptr->om_ptr;
	  int degree2=v2_ptr->degree;
	  for(int k=0; k<degree2; ++k) {
	    if(om2_ptr->v_ptr->state==2) {
	      --(ActiveDegree[om2_ptr->v_ptr->index]);
	      --(ActiveDegree[v2_ptr->index]);
	      Focusobserved.insert( om2_ptr->v_ptr->index);
	    }
	    ++om2_ptr;
	  }
	}
	else if(v2_ptr->state==2) Focusobserved.insert(v2_ptr->index);
	++om_ptr;
      }
    }
    /* (3) treat focused observed vertices */
    Leafvertices.clear();
    while( !Focusobserved.empty() ) {
      v_ptr=&Vertex[ *(Focusobserved.begin() )];
      Focusobserved.erase(v_ptr->index);
      if(ActiveDegree[v_ptr->index]==0) v_ptr->state=3;
      else if(ActiveDegree[v_ptr->index]==1) {
	v_ptr->state=3;             //unoccupied, observed, and inactive vertex
	ActiveDegree[v_ptr->index]=0;
	om_ptr=v_ptr->om_ptr;
	int degree=v_ptr->degree;
	bool found=false;
	for(int d=0; d<degree && !found; ++d) {
	  v2_ptr=om_ptr->v_ptr;
	  if(v2_ptr->state==0) {
	    found=true;
	    --(ActiveDegree[v2_ptr->index]);
	    if( ActiveDegree[v2_ptr->index] ==0 ) {
	      v2_ptr->state=1;                              //ocupy this vertex
	      --NumUnobserved;
	      Leafvertices.erase(v2_ptr->index);
	    }
	    else if( ActiveDegree[v2_ptr->index]==1 )
	      Leafvertices.insert(v2_ptr->index);
	  }
	  ++om_ptr;
	}
      }
    }
  }
  if(NumUnobserved==0) {
    cout<<"No need to perform BPD.\n";
    return false;
  }
  v_ptr=&Vertex[1];
  ActiveVertexNumber=0;
  for(int v=1; v<=VertexNumber; ++v) {
    if( v_ptr->state % 2 == 0) {
      Permutation[ActiveVertexNumber]=v_ptr->index;
      ++ActiveVertexNumber;
    }
    ++v_ptr;
  }
  MaxFixNumber=static_cast<int>(NumUnobserved*FixFraction);
  if(MaxFixNumber==0) MaxFixNumber = 1;
  try { CandidateVertices.resize( MaxFixNumber * 101 ); }
  catch( bad_alloc ) {
    cerr<<"CandidateVertices construction failed.\n";
    return false;    }
  CandidateSize.resize(101); 
  MinRange=0;
  MaxRange=0;
  cout<<"Core size: "<<NumUnobserved<<endl;
  return true;
}


/*--- set the value of Beta ---*/
void DomSET::SetBeta(double bval)
{
  Beta=bval;
  WeightOccupy=exp(-Beta);
  return;
}


/*---                set the value of DampingFactor                       ---*/
void DomSET::SetDampingFactor(double dp)
{
  if(dp >= 1.0e0)            {    Rnew=1.0e0;    Rold=0.0e0;     }
  else if(dp <= 0.0001e0)    {    Rnew=0.0001e0; Rold=0.9999e0;  }
  else                       {    Rnew=dp;       Rold=1.0e0-dp;  }
  return;
}


/*---                      initialize population                          ---*/
void DomSET::Initialization(void)
{
  struct vertexstruct *v_ptr=&Vertex[1];
  for(int v=1; v<=VertexNumber; ++v) {
    struct message *im_ptr=v_ptr->im_ptr;
    for(int d=0; d<v_ptr->degree; ++d) {
      im_ptr->q_00=0.25e0;
      im_ptr->q_01=0.25e0;
      im_ptr->q_10=0.25e0;
      ++im_ptr;
    }
    ++v_ptr;
  }
  return;
}


/*-                Belief propagation                                     -*/
bool DomSET::BeliefPropagation(double error, int count)
{
  int iter=0;
  double max_error=0;
  do {
    max_error=0;
    for(int quant=ActiveVertexNumber; quant>0; --quant) {
      int iii = static_cast<int>(quant * u01prn(PRNG) );
      int v = Permutation[iii];
      Permutation[iii] = Permutation[quant-1];
      Permutation[quant-1] = v;
      UpdateMessage(&Vertex[v], max_error);
    }
    /*
      if((iter % 10) == 0) { 
      cerr<<' '<<max_error<<' ';
      cerr.flush();
      }
    */
    ++iter;
  } while(max_error>error && iter<count);
  if(max_error<=error) {
    //    cerr<<' '<<max_error<<"  :-)\n" ;
    return true;       }
  else             {
    //    cerr<<' '<<max_error<<"  :-(\n";
    return false;  }
}


/*---            Update the output messages from a vertex                 ---*/
void DomSET::UpdateMessage(struct vertexstruct *v_ptr,       //the focal vertex
			   double& maxdiff)      //max diff. of output messages
{
  if(v_ptr->state==0) {                      //vertex unoccupied and unobserved
    struct outmessage *om_ptr=v_ptr->om_ptr;
    for(int j=0; j<v_ptr->degree; ++j) {
      if(om_ptr->v_ptr->state % 2 == 0) {
	IsActive[j]=true;
	Weight00[j]=1.0e0;
	Weightx0[j]=1.0e0;
	Weightx1[j]=WeightOccupy;
      }
      else IsActive[j]=false;
      ++om_ptr;
    }
    struct message *im_ptr=v_ptr->im_ptr;
    for(int j=0; j<v_ptr->degree; ++j) {
      if(IsActive[j]) {
	double q_00 = im_ptr->q_00;
	double q_x0 = im_ptr->q_00+im_ptr->q_10;
	double q_x1 = im_ptr->q_01+im_ptr->q_10;
	for(int k=0; k<v_ptr->degree; ++k) {
	  if(k != j && IsActive[k]) {
	    Weight00[k] *= q_00;
	    Weightx0[k] *= q_x0;
	    Weightx1[k] *= q_x1;
	    double maxval=max(Weightx0[k], Weightx1[k]);
	    Weight00[k] /= maxval;                             //avoid underflow
	    Weightx0[k] /= maxval;                             //avoid underflow
	    Weightx1[k] /= maxval;                             //avoid underflow
	  }
	}
      }
      ++im_ptr;
    }
    om_ptr=v_ptr->om_ptr;
    for(int j=0; j<v_ptr->degree; ++j) {
      if(IsActive[j]) {
	double q_00=Weightx0[j]-Weight00[j];
	double q_01=Weightx0[j];
	double q_10=Weightx1[j];
	double norm=q_00+q_01+2.0e0*q_10;
	q_00 /= norm;
	q_01 /= norm;
	q_10 /= norm;
	double diff=abs(q_00-om_ptr->m_ptr->q_00)+abs(q_01-om_ptr->m_ptr->q_01)
	  +abs(q_10-om_ptr->m_ptr->q_10);
	if(diff>maxdiff) maxdiff=diff;
	q_00=Rnew*q_00+Rold*om_ptr->m_ptr->q_00;
	q_01=Rnew*q_01+Rold*om_ptr->m_ptr->q_01;
	q_10=Rnew*q_10+Rold*om_ptr->m_ptr->q_10;
	norm=q_00+q_01+2.0e0*q_10;
	om_ptr->m_ptr->q_00 = q_00/norm;
	om_ptr->m_ptr->q_01 = q_01/norm;
	om_ptr->m_ptr->q_10 = q_10/norm;
      }
      ++om_ptr;
    }
  }
  else {                        //v_ptr->state=2, vertex unoccupied but observed
    struct outmessage *om_ptr=v_ptr->om_ptr;
    for(int j=0; j<v_ptr->degree; ++j) {
      if( om_ptr->v_ptr->state != 0) IsActive[j]=false;
      else {
	IsActive[j]=true;
	Weightx0[j]=1.0e0;
	Weightx1[j]=WeightOccupy;
      }
      ++om_ptr;
    }
    struct message *im_ptr=v_ptr->im_ptr;
    for(int j=0; j<v_ptr->degree; ++j) {
      if( IsActive[j]) {
	double q_x0=im_ptr->q_00+im_ptr->q_10;
	double q_x1=im_ptr->q_01+im_ptr->q_10;
	for(int k=0; k<v_ptr->degree; ++k) {
	  if((k != j) && IsActive[k]) {
	    Weightx0[k] *= q_x0;
	    Weightx1[k] *= q_x1;
	    double maxval=max(Weightx0[k], Weightx1[k]);
	    Weightx0[k] /= maxval;                             //avoid underflow
	    Weightx1[k] /= maxval;                             //avoid underflow
	  }
	}
      }
      ++im_ptr;
    }
    om_ptr = v_ptr->om_ptr;
    for(int j=0; j<v_ptr->degree; ++j) {
      if(IsActive[j]) {
	double q_00=Weightx0[j];
	double q_10=Weightx1[j];
	double norm=2.0e0*(q_00+q_10);
	q_00 /= norm;
	q_10 /= norm;
	double diff=abs(q_00-om_ptr->m_ptr->q_00)+abs(q_00-om_ptr->m_ptr->q_01)
	  +abs(q_10-om_ptr->m_ptr->q_10);
	if(diff>maxdiff) maxdiff=diff;
	q_00=Rnew*q_00+Rold*om_ptr->m_ptr->q_00;
	q_10=Rnew*q_10+Rold*om_ptr->m_ptr->q_10;
	norm=2.0e0*(q_00+q_10);
	om_ptr->m_ptr->q_00 = q_00/norm;
	om_ptr->m_ptr->q_01 = q_00/norm;
	om_ptr->m_ptr->q_10 = q_10/norm;
      }
      ++om_ptr;
    }
  }
  return;
}


/* thermodynamic quantities of the active subgraph */
void DomSET::Thermodynamics(void)
{
  if(NumUnobserved==0) return;              //all vertices occupied or observed
  cout<<NumUnobserved <<'\t' << ActiveVertexNumber <<endl; cout.flush();
  MaxFixNumber=static_cast<int>(NumUnobserved*FixFraction);
  if(MaxFixNumber==0) MaxFixNumber=1;
  CandidateSize=0;                         //CandidateSize[r] = 0 for 0<=r<=100
  MinRange=0;        //maximal r0: sum_{r>=r0} CandidateSize[r] >= MaxFixNumber
  MaxRange=0;
  /*
    double bf_vtx=0;                  //vertex contribution to beta*free_energy
    double bf_edge=0;                   //edge contribution to beta*free_energy
    double rho_vtx=0;                               //vertex occupation density
  */
  double q_1=0;
  struct vertexstruct *v_ptr;
  struct message *im_ptr;
  struct outmessage *om_ptr;
  for(int v_index=0; v_index < ActiveVertexNumber; ++v_index) {
    v_ptr= &Vertex[ Permutation[v_index] ];
    if(v_ptr->state==0) {                    //vertex unoccupied and unobserved
      double w00=1.0e0;
      double wx0=1.0e0;
      double wx1=WeightOccupy;
      im_ptr = v_ptr->im_ptr;
      om_ptr = v_ptr->om_ptr;
      for(int j=0; j<v_ptr->degree; ++j) {
	if(om_ptr->v_ptr->state != 3) {
	  w00 *= im_ptr->q_00;
	  wx0 *= (im_ptr->q_00+im_ptr->q_10);
	  wx1 *= (im_ptr->q_01+im_ptr->q_10);
	  double maxval=max(wx0, wx1);
	  w00 /= maxval;                                       //avoid underflow
	  wx0 /= maxval;                                       //avoid underflow
	  wx1 /= maxval;                                       //avoid underflow
	  /*
	    if(v_ptr->index < om_ptr->v_ptr->index)
	    bf_edge += -log(im_ptr->q_00 * om_ptr->m_ptr->q_00 + 
	    im_ptr->q_10 * om_ptr->m_ptr->q_01 +
	    im_ptr->q_01 * om_ptr->m_ptr->q_10 +
	    im_ptr->q_10 * om_ptr->m_ptr->q_10);
	  */
	}
	++im_ptr;
	++om_ptr;
      }
      double norm=(wx0-w00)+wx1;
      q_1 = wx1/norm;
      // bf_vtx += -log(norm);
      // rho_vtx += q_1;
    }
    else {              //v_ptr->state=2, vertex unoccupied, observed and active
      double wx0=1.0e0;
      double wx1=WeightOccupy;
      im_ptr = v_ptr->im_ptr;
      om_ptr = v_ptr->om_ptr;
      for(int j=0; j<v_ptr->degree; ++j) {
	if( om_ptr->v_ptr->state  == 0) {                     //neighbor active
	  wx0 *= (im_ptr->q_00+im_ptr->q_10);
	  wx1 *= (im_ptr->q_01+im_ptr->q_10);
	  double maxval=max(wx0, wx1);
	  wx0 /= maxval;                                       //avoid underflow
	  wx1 /= maxval;                                       //avoid underflow
	  /*
	    if(v_ptr->index < om_ptr->v_ptr->index)
	    bf_edge += -log(im_ptr->q_00 * om_ptr->m_ptr->q_00 + 
	    im_ptr->q_10 * om_ptr->m_ptr->q_01 +
	    im_ptr->q_01 * om_ptr->m_ptr->q_10 +
	    im_ptr->q_10 * om_ptr->m_ptr->q_10);
	  */
	}
	++im_ptr;
	++om_ptr;
      }
      double norm=wx0+wx1;
      q_1 =wx1/norm;
      //      bf_vtx += -log(norm);
      //      rho_vtx += q_1;
    }
    int rrr = static_cast<int>(q_1 * 100);
    if(rrr>=MinRange) {
      if(rrr>MaxRange) MaxRange=rrr;
      if(CandidateSize[rrr]<MaxFixNumber) {
	CandidateVertices[rrr*MaxFixNumber+CandidateSize[rrr] ] =v_ptr->index;
	CandidateSize[rrr] += 1;
      }
      else MinRange=rrr;
    }
  }
  /*
    bf_vtx  /= ActiveVertexNumber;
    bf_edge /= ActiveVertexNumber;
    rho_vtx /= ActiveVertexNumber;
    double phi = bf_vtx-bf_edge;
    output<< rho_vtx <<'\t'
    << (Beta>0? phi/Beta : 0)<<'\t'
    << Beta *rho_vtx - phi << endl;
  */
  return;
}


/* ---       externally fixing some variables to be occupied         ---*/
bool DomSET::Fix(void)
{
  int rank = MaxRange;         /* occupy_prob distributed to 101 bins [0,0.01),
				  [0.01,0.02), ..., [0.99,1), [1,1] */
  int num_examined = 0;                 //number of examined candidate vertices
  while(num_examined < MaxFixNumber && NumUnobserved > 0) {
    int size=CandidateSize[rank];
    if(size>0) {
      int *i_ptr = &CandidateVertices[rank * MaxFixNumber];
      if(num_examined+size <= MaxFixNumber) {
	num_examined +=size;
	for(int s=0; s<size && NumUnobserved>0; ++s) {
	  struct vertexstruct *v_ptr = &Vertex[ *i_ptr ];
	  if(v_ptr->state==0) --NumUnobserved;
	  v_ptr->state=1;
	  struct outmessage *om_ptr=v_ptr->om_ptr;
	  for(int k=0; k<v_ptr->degree; ++k) {
	    if(om_ptr->v_ptr->state==0) {
	      om_ptr->v_ptr->state=2;
	      --NumUnobserved;
	    }
	    ++om_ptr;
	  }
	  ++i_ptr;
	}
      }
      else {
	while(num_examined < MaxFixNumber && NumUnobserved>0) {
	  ++num_examined;
	  struct vertexstruct *v_ptr = &Vertex[ *i_ptr ];
	  if(v_ptr->state==0) --NumUnobserved;
	  v_ptr->state=1;
	  struct outmessage *om_ptr=v_ptr->om_ptr;
	  for(int k=0; k<v_ptr->degree; ++k) {
	    if(om_ptr->v_ptr->state==0) {
	      om_ptr->v_ptr->state=2;
	      --NumUnobserved;
	    }
	    ++om_ptr;
	  }
	  ++i_ptr;
	}
      }
    }
    --rank;
  }
  if(NumUnobserved>0) {
    for(int i=ActiveVertexNumber-1; i>=0; --i) {
      struct vertexstruct *v_ptr = &Vertex[ Permutation[i] ];
      if( v_ptr->state % 2 ==1) {
	--ActiveVertexNumber;
	Permutation[i]=Permutation[ActiveVertexNumber];
      }
      else if(v_ptr->state==2) {
	int ac_n=0;
	struct outmessage *om_ptr=v_ptr->om_ptr;
	for(int k=0; k<v_ptr->degree; ++k) {
	  if(om_ptr->v_ptr->state==0) ++ac_n;
	  ++om_ptr;
	}
	if(ac_n<=1) {
	  v_ptr->state=3;
	  --ActiveVertexNumber;
	  Permutation[i]=Permutation[ActiveVertexNumber];
	}
      }
    }
    return true;
  }
  else
    return false;
}


/* check whether the final occupation pattern corresponds to a Dom-Set */
bool DomSET::CheckDS(int& num_occup,              //number of occupied vertices
		     int& num_obser,    //# of unoccupied but observed vertices
		      string& filename)               //the dominating set
{
  num_occup=0;
  num_obser=0;
  struct vertexstruct *v_ptr = &Vertex[1];
  struct outmessage *om_ptr;
  set<int> occupiedvertices;
  bool succeed=true;
  for(int v=1; v<=VertexNumber && succeed; ++v) {
    if(v_ptr->state==1) {                                     //vertex occupied
      ++num_occup;
      occupiedvertices.insert(v_ptr->index);
    }
    else if(v_ptr->state>=2) {                 //vertex unoccupied but observed
      ++num_obser;
      om_ptr=v_ptr->om_ptr;
      bool found=false;
      for(int k=0; k<v_ptr->degree && !found; ++k) {
	if(om_ptr->v_ptr->state==1) found=true;
	++om_ptr;
      }
      if(!found) succeed=false;
    }
    else                                     //vertex unoccupied and unobserved
      succeed=false;
    ++v_ptr;
  }
  if(succeed) {
    ofstream output(filename.c_str() );
    output<<num_occup <<'\t' <<num_obser <<endl<<endl;
    for(set<int>::const_iterator sci=occupiedvertices.begin();
	sci !=occupiedvertices.end(); ++sci)
      output<< *sci <<endl;
    output.close();
  }
  return succeed;
}