//
// this program implements the site percolation algorithm from
// Newman and Ziff, PRE64, 016706
// 

#include <fstream>
#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <random>


//------------------------------------------------------

void setCellNeigbhours(int length, int latticeSize, vector<int> &cellNeighbours)
{
  for(int i=0;i<latticeSize;i++) {
    cellNeighbours[i] = (i+1)%latticeSize;
    cellNeighbours[i+latticeSize] = (i+latticeSize-1)%latticeSize;
    cellNeighbours[i+2*latticeSize] = (i+length)%latticeSize;
    cellNeighbours[i+3*latticeSize] = (i+latticeSize-length)%latticeSize;
    // wrap horizontally
    if(i%length==0) cellNeighbours[i+latticeSize] = i+length-1; 
    if((i+1)%length==0) cellNeighbours[i] = i-length+1;     
  }
}

//------------------------------------------------------

int findroot(int i, vector<int> &cellPointer)
// implements path compression
{
  if (cellPointer[i]<0) return i;
  return cellPointer[i] = findroot(cellPointer[i],cellPointer);
}


 /*
int findroot(int i, vector<int> &cellPointer)
// implements path halving
{
  int r, s;
  r=s=i;
  while(cellPointer[r]>=0) {
    cellPointer[s]=cellPointer[r];
    s=r;
    r=cellPointer[r];
  }
  return r;
}
 */
//------------------------------------------------------

void sitePercolation(int length)
// length = size of one side of a 2D square lattice
{
  // initialise some internal variables
  int latticeSize = length*length; // lattice size
  int emptyCell = -(latticeSize+1); // mark empty cells
  std::vector<int> cellPointer(latticeSize,emptyCell); // assign cells to clusters

  // initialise observables
  int biggestCluster = 0;
  
  // initialise order of occupying the cells in the lattice
  std::vector<int> cellOrder(latticeSize); // create vector
  std::iota(cellOrder.begin(), cellOrder.end(), 0); // fill it form 0 to latticeSize-1
  std::shuffle(cellOrder.begin(), cellOrder.end(), std::random_device()); // suffle order

  // initialise neighbors of cell
  std::vector<int> cellNeighbours(latticeSize*4); // each cell has 4 neighbours
  setCellNeigbhours(length,latticeSize,cellNeighbours);

  // percolate
  int s1, s2, r1, r2; // internal indices for the pointers vector
  for(int i=0; i<latticeSize; i++) { // loop over cells
    r1 = s1 = cellOrder[i]; // new cell
    cellPointer[s1] = -1; // current size of the cluster
    for(int j=0;j<4;j++) {// loop over neighbours
      s2 = cellNeighbours[s1+j*latticeSize]; // index of neighbouring cell
      if(cellPointer[s2] != emptyCell) { // cell not empty, form a cluster
	r2 = findroot(s2,cellPointer); // find representative of the cluster of this cell
	if (r2 != r1) {// merge clusters: smaller cluster is absorved
	  if (cellPointer[r1] > cellPointer[r2]) { // cluster size is negative for root nodes!
	    cellPointer[r2] += cellPointer[r1];
	    cellPointer[r1] = r2;
	    r1 = r2;
	  } else {
	    cellPointer[r1] += cellPointer[r2];
	    cellPointer[r2] = r1;
	  }
	  // fill the observable
	  if (-cellPointer[r1]>biggestCluster) biggestCluster = -cellPointer[r1];
	} // end of merging
      } // end cell not empty
    } // end loop over neighbours
    // print out the observable
    std::cout << i << "  " << biggestCluster  << endl;
  } // end loop over cells
}