/*
  GLRandGreedyDSV02b.cpp
  
  A hybrid GLR and greedy (occupy the vertex of highest active degree) algorithm
  for the dominating set problem.
  
  Input:  a simple graph G with N vertices and M edges.
  Output: a dominating set.
  
  GLR:
  (0) S is an empty set; each vertex is being unoccupied and unobserved.
  (1) Add all the isolated vertices (whose degree is zero) to S and delete them
  from G. These vertices are now in the occupied state.
  (2) For each leaf vertex i (with a single attached edge) of G, if it is in the
  unoccupied and unobserved state, then change its unique nearest neighboring
  vertex (j) into the occupied state, and change all the unoccupied and 
  unobserved nearest neighboring vertices of j into the observed state.
  (3) After step (2) has been finished for all the leaf vertices, then simplify
  the graph G by deleting all the occupied vertices and the associated edges,
  and delete all the edges between any two observed vertices.
  (4) If in the simplified graph, an observed vertex has only one edge attached,
  then remove this vertex and the associated edge from the graph. Continue this
  process until no such observed vertex remains in the graph.
  (5) For the simplified graph, if there is at least one isolated or leaf
  vertex, then repeat steps (1)--(5).
  (6) Output the set of occupied vertices, the set of observed vertices, and the
  core (the set of unoccupied and unobserved vertices).
  
  Greedy:
  If there is isolated or leaf vertex, perform the GLR process. If no such
  vertices and the graph is not empty, then randomly choose a vertex among the
  subset of vertices of highest impact and occupy that vertex.
  The impact of an unobserved vertex is the total number of nearest neighboring
  unobserved vertices pluse one (the focal vertex itself).
  The impact of an observed but not occupied vertex is the total number of
  nearest neighboring unobserved vertices.
 
  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:
  30.09.2014: last modified.
  30.09.2014: copied from GLRandGreedyDSmpiV02b.cpp to GLRandGreedyDSV02b.cpp
  and modied to single mode (treat single graph instances).
  19.09.2014: last modified.
  19.09.2014: copied from GLRandGreedyDSmpiV02.cpp to GLRandGreedyDSmpiV02b.cpp
  and further checked (an additional slight bug was discovered).
  20--23.08.2014: further checked (a bug was discovered while discussing with
  Chuang Wang on 21.08.2014).
  31.07.2014: copied from GLRandGreedyDSv01.cpp (modified for parallel).
  30.07.2014: last modified.
  29.07.2014: copied from GLRforDSv01.cpp and modified.
  
  PROGRAMMER:
  Hai-Jun Zhou
  Institute of Theoretical Physics, the Chinese Academy of Sciences,
  Zhong-Guan-Cun-Dong-Lu 55, Beijing 100190, China
  zhouhj@itp.ac.cn
*/

#include <cmath>
#include <exception>
#include <fstream>
#include <iostream>
#include <map>
#include <set>
#include <string>
#include <valarray>
#include "/Users/zhouhj/Programs/zhjrandom.h"

using namespace std;


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


struct OrderedPair
{
  int first;
  int second;
  OrderedPair(void);
  OrderedPair(int a, int b);
};

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

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

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

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


struct vstruct
{
  int vertex;                                        //vertex index 1, 2, 3, ...
  int degree;                                      //number of nearest neighbors
  int active_degree;   /* for an unoccupied and unobserved vertex, active_degree
			  is the number of unobserved and active observed
			  neighboring vertices; for an observed vertex, active_
			  degree is the number of unobserved neighboring 
			  vertices. */
  int *nv_ptr;                     //starting position of nearest neighbor array
  int state;               /* state = 0: unoccupied and unobserved
			      state = 1: unoccupied but observed, and active
			      state = 2: occupied
			      state = 3: unoccupied but observed and inactive */
  vstruct(void);
};

vstruct::vstruct(void)
{
  vertex=0;
  degree=0;
  active_degree=0;
  nv_ptr=0;
  state=0;
}


class DomSetCore
{
public:
  DomSetCore( string&, int, bool&, ZHJRANDOMv3* );
  bool Percolation(void);
  void Greedy(void);
  bool Results( string&);
  
private:
  int VertexNumber;                                         //number of vertices
  int EdgeNumber;                                              //number of edges
  int NumUnobserved;                             //number of unobserved vertices
  int NumOccupied;                        //number of directly observed vertices
  int NumCandidate;             //number of candidate vertices in greedy process
  int NumActive;                                     //number of active vertices
  int MaxImpact;               //maximal vertex impact (used for greedy process)
  valarray<vstruct> Vertex;                              //vertices of the graph
  valarray<int> Neighbors;                                  //edges of the graph
  valarray<int> ActiveArray;                           //list of active vertices
  valarray<int> CandidateArray;    //list of candidate vertices (maximal impact)
  valarray<int> Impact;                                  //impact of each vertex
  set<int> Leafvertices;                                  //set of leaf vertices
  set<int> Newlyoccupied;                       //set of newly occupied vertices
  set<int> Focusobserved;     //set of observed vertices that need to be treated
  ZHJRANDOMv3 * PRNG;                                  //random number generator
};


int main(int argc, char ** argv)
{

  int rd_seed=27443854;
  ZHJRANDOMv3 PRNG( rd_seed );
  for(int i=0; i<2200058; ++i) PRNG.rdflt();
  
  int VertexNumber =100000;                                     //vertex number
  int EdgeNumber   =500000;                                      //edge number
  
  bool succeed=true;
  DomSetCore sample("ERn100kM1m.g100", EdgeNumber, succeed, &PRNG);
  if(succeed) {
    succeed=sample.Percolation();
    if(succeed) {
      sample.Greedy();
      succeed=sample.Results("myresult2.dat");
    }
  }
  
  return 1;
}

/*                             Read graph                                     */
DomSetCore::DomSetCore( string& graphname, int e_num, bool& succeed,
		       ZHJRANDOMv3* rd)
{
  PRNG=rd;
  succeed=true;
  
  ifstream graph(graphname.c_str() );
  if(!graph.good() ) {
    cerr<<"Graph probably non-existant.\n";
    succeed=false;
    return;
  }
  graph >>VertexNumber >>EdgeNumber;
  if(EdgeNumber<e_num) {
    cerr<<"Not so many edges in the graph.\n";
    graph.close();
    succeed=false;
    return;
  }
  try { Vertex.resize(VertexNumber+1); }
  catch(bad_alloc) {
    cerr<<"Error in constructing Vertex.\n";
    succeed=false;
    return;
  }
  try { Impact.resize(VertexNumber+1); }
  catch(bad_alloc) {
    cerr<<"Error in constructing Impact.\n";
    succeed=false;
    return;
  }
  
  EdgeNumber=e_num;
  try { Neighbors.resize(2*EdgeNumber); }
  catch(bad_alloc) {
    cerr<<"Error in constructing Neighbors.\n";
    succeed=false;
    return;
  }
  
  //first reading for counting
  for(int eindex=0; eindex<EdgeNumber; ++eindex) {
    int v1,v2;
    graph >>v1 >>v2;
    if(v1==v2 || v1 <= 0 || v1>VertexNumber || v2 <=0 || v2 > VertexNumber) {
      cerr<<"Input graph not correct at line "<<eindex+1<<endl;
      graph.close();
      succeed=false;
      return;
    }
    ++(Vertex[v1].degree);
    ++(Vertex[v2].degree);
  }
  graph.close();
  int position=0;
  Leafvertices.clear();
  NumUnobserved=VertexNumber;
  NumOccupied=0;
  struct vstruct *v_ptr=&Vertex[1];
  for(int v=1; v<=VertexNumber; ++v) {
    v_ptr->active_degree=0;
    v_ptr->state=0;
    v_ptr->vertex=v;
    v_ptr->nv_ptr=&Neighbors[position];
    position += v_ptr->degree;
    if(v_ptr->degree==0) {
      v_ptr->state=2;                                   //occupy isolated vertex
      --NumUnobserved;
      ++NumOccupied;
    }
    else if(v_ptr->degree ==1) Leafvertices.insert(v_ptr->vertex);
    ++v_ptr;
  }
  
  //second reading for setting up neighborhood
  graph.open(graphname.c_str());
  graph >>VertexNumber >>EdgeNumber;
  EdgeNumber=e_num;
  for(int eindex=0; eindex<EdgeNumber; ++eindex) {
    int v1, v2;
    graph  >>v1 >>v2;
    int *nv_ptr = Vertex[v1].nv_ptr + Vertex[v1].active_degree;
    *nv_ptr = v2;
    ++(Vertex[v1].active_degree);
    nv_ptr=Vertex[v2].nv_ptr + Vertex[v2].active_degree;
    *nv_ptr = v1;
    ++(Vertex[v2].active_degree);
  }
  graph.close();
  
  //check for edge multiplicity
  set<int> nearest_neighbors;
  for(int v=1; v<=VertexNumber; ++v) {
    nearest_neighbors.clear();
    int* nv_ptr=Vertex[v].nv_ptr;
    int degree=Vertex[v].degree;
    if(degree != Vertex[v].active_degree) {
      cout<<"Error in graph construction.\n";
      succeed=false;
    }
    for(int d=0; d<degree; ++d) {
      if(nearest_neighbors.find( *nv_ptr ) != nearest_neighbors.end()) {
	cout<<"multiple edges.\n";
	succeed=false;
      }
      else nearest_neighbors.insert( *nv_ptr);
      ++nv_ptr;
    }
  }
  // if(succeed) cout<<"Graph simple. No multiple edge or self connection.\n";
  if(NumUnobserved==0) {
    succeed=false;
    // cout<<"Graph trivial, all vertices should be directly observed.\n";
  }

  NumActive=VertexNumber-NumOccupied;     //isolated vertices have been occupied
  
  return;
}


/*  Initial GLR process, return true if greedy process needed  */
bool DomSetCore::Percolation(void)
{
  bool succeed=true;
  struct vstruct *v_ptr =&Vertex[1];
  struct vstruct *v2_ptr=&Vertex[1];
  
  while( !Leafvertices.empty() ) {                 /* (1) treat leaf vertices */
    Newlyoccupied.clear();
    do {
      v_ptr=&Vertex[ *(Leafvertices.begin() )];
      Leafvertices.erase(v_ptr->vertex);
      if(v_ptr->state==0) { //if state !=0, then state=2 (neighboring to a leaf)
	int *nv_ptr = v_ptr->nv_ptr;
	int degree=v_ptr->degree;
	bool found=false;
	for(int d=0; d<degree && !found; ++d) {
	  v2_ptr=&Vertex[ *nv_ptr ];
	  if(v2_ptr->state <=1) {                                     //=0 or =1
	    found=true;
	    ++NumOccupied;
	    --NumActive;
	   
	    /* effectively deleting the edges associated with a newly occupied
	       (directly observed) vertices. */
	    if(v2_ptr->state==0) {
	      --NumUnobserved;
	      int *nv2_ptr=v2_ptr->nv_ptr;
	      int degree2=v2_ptr->degree;
	      for(int d2=0; d2<degree2; ++d2) {
		if(Vertex[ *nv2_ptr].state <= 1)
		  --(Vertex[ *nv2_ptr ].active_degree);
		++nv2_ptr;
	      }
	    }
	    else {                                            // v2_ptr->state=1
	      int *nv2_ptr=v2_ptr->nv_ptr;
	      int degree2=v2_ptr->degree;
	      for(int d2=0; d2<degree2; ++d2) {
		if(Vertex[ *nv2_ptr ].state==0) 
		  --(Vertex[ *nv2_ptr ].active_degree);
		++nv2_ptr;
	      }
	    }

	    /* occupy the vertex */
	    v2_ptr->state=2;
	    v2_ptr->active_degree=0;
	    Newlyoccupied.insert(v2_ptr->vertex);
	  }
	  ++nv_ptr;
	}
      }
    } while( !Leafvertices.empty() );
    
    Focusobserved.clear();
    while( !Newlyoccupied.empty() ) {    /* (2) treat newly occupied vertices */
      v_ptr=&Vertex[ *(Newlyoccupied.begin() )];
      Newlyoccupied.erase(v_ptr->vertex);
      int *nv_ptr=v_ptr->nv_ptr;
      int degree=v_ptr->degree;
      for(int d=0; d<degree; ++d) {
	v2_ptr=&Vertex[ *nv_ptr ];
	if(v2_ptr->state==0) {
	  v2_ptr->state=1;
	  --NumUnobserved;
	  Focusobserved.insert(v2_ptr->vertex);
	  //delete the edge between two observed vertices
	  int *nv2_ptr=v2_ptr->nv_ptr;
	  int degree2=v2_ptr->degree;
	  for(int k=0; k<degree2; ++k) {
	    if(Vertex[*nv2_ptr].state==1) {
	      --(Vertex[*nv2_ptr].active_degree);
	      --(v2_ptr->active_degree);
	      Focusobserved.insert( *nv2_ptr);
	    }
	    ++nv2_ptr;
	  }
	}
	else if(v2_ptr->state==1) Focusobserved.insert(v2_ptr->vertex);
	++nv_ptr;
      }
    }
    
    Leafvertices.clear();
    while( !Focusobserved.empty() ) {  /* (3) treat focused observed vertices */
      v_ptr=&Vertex[ *(Focusobserved.begin() )];
      Focusobserved.erase(v_ptr->vertex);
      if(v_ptr->active_degree==0) {
	v_ptr->state=3;
	--NumActive;
      }
      else if(v_ptr->active_degree==1) {
	v_ptr->state=3;              //unoccupied, observed, and inactive vertex
	--NumActive;
	v_ptr->active_degree=0;
	int *nv_ptr=v_ptr->nv_ptr;
	int degree=v_ptr->degree;
	bool found=false;
	for(int d=0; d<degree && !found; ++d) {
	  v2_ptr=&Vertex[ *nv_ptr ];
	  if(v2_ptr->state==0) {
	    found=true;
	    --(v2_ptr->active_degree);
	    if(v2_ptr->active_degree==0) {
	      v2_ptr->state=2;                               //ocupy this vertex
	      --NumUnobserved;
	      ++NumOccupied;
	      --NumActive;
	      Leafvertices.erase(v2_ptr->vertex);
	    }
	    else if(v2_ptr->active_degree==1)
	      Leafvertices.insert(v2_ptr->vertex);
	  }
	  ++nv_ptr;
	}
      }
    }
  }
  if(NumUnobserved==0) {
    //cout<<"No need to perform greedy process.\n";
    return false;
  }
  
  // cout<<"Number of unobserved vertices after the intial GLR process: "
  // <<NumUnobserved<<endl;
  
  try { ActiveArray.resize(NumActive); }
  catch(bad_alloc) {
    cerr<<"Error in ActiveArray construction.\n";
    return false;
  }
  try { CandidateArray.resize(NumActive); }
  catch(bad_alloc) {
    cerr<<"Error in CandidateArray construction.\n";
    return false;
  }
  
  int position=0;
  succeed=true;
  for(int v=1; v<=VertexNumber; ++v) {
    if(Vertex[v].state<=1) {
      if(position<NumActive) {
	ActiveArray[position]=v;
	++position;
      }
      else if(position>=NumActive) {
	cerr<<"Number of active vertices does not consistent.\n";
	succeed=false;
      }
    }
  }
  return succeed;
}


/* Greedy and GLR combined to obtain a dominating set */
void DomSetCore::Greedy(void)
{
  struct vstruct *v_ptr;
  struct vstruct *v2_ptr;
  
  /* determine initial impact of each vertex before greedy process starts */  
  for(int v=1; v<=VertexNumber; ++v) {
    if(Vertex[v].state>=2) Impact[v]=0;                        //vertex inactive
    else if(Vertex[v].state==1)
      Impact[v]=Vertex[v].active_degree;       // # unobserved nearest neighbors
    else {                                                             //state=0
      int value=1;                          //# unobserved nearest neighbors + 1
      int *nv_ptr=Vertex[v].nv_ptr;
      int degree=Vertex[v].degree;
      for(int d=0; d<degree; ++d) {
	if(Vertex[ *nv_ptr ].state==0) ++value;
	++nv_ptr;
      }
      Impact[v]=value;
    }
  }
  
  MaxImpact=0;
  NumCandidate=0;
  while(NumUnobserved>0) {                    //still having unobserved vertices
    if(NumCandidate>0)
      for(int index=NumCandidate-1; index>=0; --index) {
	int vtx=CandidateArray[index];
	if(Vertex[vtx].state>=2) {                            //no longer active
	  --NumCandidate;
	  CandidateArray[index]=CandidateArray[NumCandidate];
	}
	else {                                                    //still active
	  if(Impact[vtx]<MaxImpact) {
	    --NumCandidate;
	    CandidateArray[index]=CandidateArray[NumCandidate];
	  }
	}
      }
    if(NumCandidate==0) {
      MaxImpact=0;
      for(int index=NumActive-1; index>=0; --index) {
	int vtx=ActiveArray[index];
	if(Vertex[vtx].state<=1) {                                //still active
	  int impact=Impact[vtx];
	  if(impact>MaxImpact) {
	    MaxImpact=impact;
	    NumCandidate=0;
	    CandidateArray[NumCandidate]=vtx;
	    ++NumCandidate;
	  }
	  else if(impact==MaxImpact) {
	    CandidateArray[NumCandidate]=vtx;
	    ++NumCandidate;
	  }
	}
	else {                                                //no longer active
	  --NumActive;
	  ActiveArray[index]=ActiveArray[NumActive];
	}
      }
    }
    
    /* greedy process: occupy a node of maximal impact */
    int myplace=static_cast<int>(NumCandidate*u01prn(PRNG) );
    v_ptr=&Vertex[ CandidateArray[myplace] ];
    --NumCandidate;
    CandidateArray[myplace]=CandidateArray[NumCandidate];
    if(v_ptr->state==0) {                                     //not yet observed
      --NumUnobserved;                          //vertex being directly observed
      int *nv_ptr=v_ptr->nv_ptr;
      int degree=v_ptr->degree;
      for(int d=0; d<degree; ++d) {
	if(Vertex[ *nv_ptr ].state<=1) {
	  --(Vertex[ *nv_ptr ].active_degree);
	  --(Impact[ *nv_ptr ]);
	}
	++nv_ptr;
      }
    }
    else {                                                      //v_ptr->state=1
      int *nv_ptr=v_ptr->nv_ptr;
      int degree=v_ptr->degree;
      for(int d=0; d<degree; ++d) {
	if(Vertex[ *nv_ptr ].state==0) --(Vertex[ *nv_ptr].active_degree);
	++nv_ptr;
      }
    }
    v_ptr->state=2;
    ++NumOccupied;
    v_ptr->active_degree = 0;
    Impact[ v_ptr->vertex] = 0;
    Newlyoccupied.clear();
    Newlyoccupied.insert(v_ptr->vertex);
    
    /* generalized leaf removal */
    do {
      Focusobserved.clear();     //observed vertices that might be leaf vertices
      while( !Newlyoccupied.empty() ) {  /* (1) treat newly occupied vertices */
	v_ptr=&Vertex[ *(Newlyoccupied.begin() )];
	Newlyoccupied.erase(v_ptr->vertex);
	int *nv_ptr=v_ptr->nv_ptr;
	int degree=v_ptr->degree;
	for(int d=0; d<degree; ++d) {         //nearest neighbors being observed
	  v2_ptr=&Vertex[ *nv_ptr ];
	  if(v2_ptr->state==0) {
	    v2_ptr->state=1;
	    --(Impact[v2_ptr->vertex]);               //line added on 19.09.2014
	    --NumUnobserved;
	    Focusobserved.insert(v2_ptr->vertex);
	    int *nv2_ptr=v2_ptr->nv_ptr;
	    int degree2=v2_ptr->degree;
	    for(int k=0; k<degree2; ++k) { //2nd nearest neighbors should update
	      if(Vertex[*nv2_ptr].state==1) {
		--(v2_ptr->active_degree);
		--(Vertex[*nv2_ptr].active_degree);
		--(Impact[ *nv2_ptr ]);
		Focusobserved.insert( *nv2_ptr);
	      }
	      else if(Vertex[ *nv2_ptr].state==0) --(Impact[ *nv2_ptr ]);
	      ++nv2_ptr;
	    }
	  }
	  else if(v2_ptr->state==1) Focusobserved.insert(v2_ptr->vertex);
	  ++nv_ptr;
	}
      }
      
      Leafvertices.clear();
      while( !Focusobserved.empty()) { /* (2) treat focused observed vertices */
	v_ptr=&Vertex[ *(Focusobserved.begin() )];
	Focusobserved.erase(v_ptr->vertex);
	if(v_ptr->active_degree==0) {
	  v_ptr->state=3;
	  Impact[v_ptr->vertex]=0;
	}
	else if(v_ptr->active_degree==1) {
	  v_ptr->state=3;            //unoccupied, observed, and inactive vertex
	  v_ptr->active_degree=0;
	  Impact[v_ptr->vertex]=0;
	  
	  int *nv_ptr=v_ptr->nv_ptr;
	  int degree=v_ptr->degree;
	  bool found=false;
	  for(int d=0; d<degree && !found; ++d) {
	    v2_ptr=&Vertex[ *nv_ptr ];
	    if(v2_ptr->state==0) {
	      found=true;
	      --(v2_ptr->active_degree);
	      if(v2_ptr->active_degree==0) {
		v2_ptr->state=2;                             //ocupy this vertex
		Impact[v2_ptr->vertex]=0;
		--NumUnobserved;
		++NumOccupied;
		Leafvertices.erase(v2_ptr->vertex);
	      }
	      else if(v2_ptr->active_degree==1) 
		Leafvertices.insert(v2_ptr->vertex);
	    }
	    ++nv_ptr;
	  }
	}
      }
      
      Newlyoccupied.clear();
      while( !Leafvertices.empty() ) {             /* (3) treat leaf vertices */
	v_ptr=&Vertex[ *(Leafvertices.begin() )];
	Leafvertices.erase(v_ptr->vertex);
	if(v_ptr->state==0) { //if state !=0, then state=2 (neighboring to leaf)
	  int *nv_ptr = v_ptr->nv_ptr;
	  int degree=v_ptr->degree;
	  bool found=false;
	  for(int d=0; d<degree && !found; ++d) {
	    v2_ptr=&Vertex[ *nv_ptr ];
	    if(v2_ptr->state==0) {
	      found=true;
	      --NumUnobserved;
	      int *nv2_ptr=v2_ptr->nv_ptr;
	      int degree2=v2_ptr->degree;
	      for(int d2=0; d2<degree2; ++d2) {
		if(Vertex[ *nv2_ptr].state<=1) {
		  --(Vertex[ *nv2_ptr ].active_degree);
		  --(Impact[ *nv2_ptr ]);
		}
		++nv2_ptr;
	      }
	    }
	    else if(v2_ptr->state==1) {
	      found=true;
	      int *nv2_ptr=v2_ptr->nv_ptr;
	      int degree2=v2_ptr->degree;
	      for(int d2=0; d2<degree2; ++d2) {
		if(Vertex[ *nv2_ptr ].state==0) 
		  --(Vertex[ *nv2_ptr ].active_degree);
		++nv2_ptr;
	      }
	    }
	    if(found) {
	      Impact[ v2_ptr->vertex]=0;
	      v2_ptr->state=2;
	      ++NumOccupied;
	      v2_ptr->active_degree=0;
	      Newlyoccupied.insert(v2_ptr->vertex);
	    }
	    ++nv_ptr;
	  }
	}
      }
    } while( !Newlyoccupied.empty() );
  }
  return;
}


/* Output result */
bool DomSetCore::Results( string& filename)
{ 
  NumOccupied=0;
  NumUnobserved=0;
  struct vstruct *v_ptr=&Vertex[1];
  struct vstruct *v2_ptr;
  bool succeed=true;
  for(int vindex=1; vindex<=VertexNumber && succeed; ++vindex) {
    if(v_ptr->state==2) {                                    //directly observed
      ++NumOccupied;
      int degree=v_ptr->degree;
      int *nv_ptr=v_ptr->nv_ptr;
      for(int d=0; d<degree; ++d) {
	v2_ptr=&Vertex[ *nv_ptr ];
	if(v2_ptr->state == 0) {
	  succeed= false;
	  cerr<<"unobserved vertex being neighbor of an occupied one.\n";
	}
	++nv_ptr;
      }
    }
    else if(v_ptr->state==0) {                                      //unobserved
      ++NumUnobserved;
      int degree=v_ptr->degree;
      int *nv_ptr=v_ptr->nv_ptr;
      int num=0;                                    //number of active neighbors
      for(int d=0; d<degree; ++d) {
	v2_ptr=&Vertex[ *nv_ptr];
	if(v2_ptr->state <= 1) ++num;
	else if(v2_ptr->state==2) {
	  succeed=false;
	  cerr<<"unobserved vertex being neighbor of an occupie one.\n";
	}
	++nv_ptr;
      }
      if(num<=1) {
	succeed=false;
	cerr<<"unobserved vertex not consistent.\n";          //leaf or isolated
      }
    }
    else {                                                 //indirectly observed
      int degree=v_ptr->degree;
      int *nv_ptr=v_ptr->nv_ptr;
      int num=0;                                  //number of occupied neighbors
      for(int d=0; d<degree; ++d) {
	v2_ptr=&Vertex[*nv_ptr];
	if(v2_ptr->state==2) ++num;
	++nv_ptr;
      }
      if(num==0) {
	succeed=false;
	cerr<<"observed vertex has no occupied neighbor.\n";
      }
    }
    ++v_ptr;
  }
  if(succeed && NumUnobserved==0) {
    ofstream output(filename.c_str());
    output <<NumOccupied<<endl<<endl;
    for(int v=1; v<=VertexNumber; ++v)
      if(Vertex[v].state==2) output<<v<<endl;
    output.close();
  }
  else {
    succeed=false;
    cerr<<"Something must be wrong.\n";
  }
  return succeed;
}