//*****************************************************************************
//
//  main.C
//  author: Eric Michael Ryckman
//  e-mail: eryckman@umich.edu
//
//  Code for Permanent Approximations, according to
//
//  "The Distance Approach to Approximate Combinatorial Counting"
//
//  by
//
//  A. Barvinok and A. Samorodnitsky
//
//
//  function LAP code by R. Jonker and A. Volgenant
//  e-mail: roy_jonker@majiclogic.com
//
//*****************************************************************************



#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#include <stdlib.h>
#include <math.h>
#include <string>
#include <time.h>
#include <climits>

using namespace std;

const int MAXDIM = 1000;                  //Dimension of matrix, or |V|/2
const int MAX = 10000;                    //const int MAX = DIM*ceil(log(DIM))+1,
                                          //largest plausible value of m
const int SETS = 2;                       //number of sets of vertices
 

int sumRows(int[MAXDIM][MAXDIM], int, int);
int sumCols(int[MAXDIM][MAXDIM], int, int);
int makem(int[MAXDIM], int[MAXDIM], int, int[MAXDIM][MAXDIM], int[MAXDIM]);
int assignWeights(int[MAXDIM][MAXDIM], int[MAX], int[MAXDIM][MAXDIM], 
		  int[MAXDIM], int);
int LAP(int, int[MAXDIM][MAXDIM], int [MAXDIM], int[MAXDIM], int[MAXDIM], 
	int[MAXDIM]);

float H(float);

void generateStrings(int[MAX], int);



//******************************************************************************
//******************************************************************************
//******************************************************************************

int main()
{
  ifstream inFile;
 
  string fileName;                             //name of input file

  int DIM;                                     //dimension of matrix
  int i, j;                                    //for iterations
  int m, M;                                    //min(m+,m-)
  int lapcost;                                 //weight of optimal matching
  
  long K, K1;                                  //number of strings to create
  long b;                                      //for large iterations
  
  float epsilon;                               //for precision
  float alphaaverage;                          //sum of alphas
  float alphaaverage1;                         //average of alphas
  float beta;                                  
  float lowerbound;
  float upperbound;
  float center;
 
  int rowSum[MAXDIM];                          
  int colSum[MAXDIM]; 
  int mi[MAXDIM];                              //mi[i]=lengths of strings
  int randomStrings[MAX];
  int rowsol[MAXDIM];                          //col matched to row
  int colsol[MAXDIM];                          //row matched to col
  int u[MAXDIM];                               //dual variable
  int v[MAXDIM];                               //dual varuable 
  int alpha[LONG_MAX];                         //bth optimal matching weight
  
  int matrix[MAXDIM][MAXDIM];
  int weights[MAXDIM][MAXDIM];                 //weight matrix
  

cout << "\n\n**************************************************************"
     << "\nThis program embeds the set of all perfect matchings of a"
     << "\ngraph into the Boolean Cube.  It then estimates the average"
     << "\nHamming distance from a random point in the cube to the set."
     << "\nBased on this distance, it provides a lower bound (achieved"
     << "\nwhen the set of perfect matchings is sparse) and an upper"
     << "\nbound (achieved when the set of perfect matchings is dense)"
     << "\nfor the number of perfect matchings of the graph."
     << "\n**************************************************************\n";


cout << "\nEnter the dimension of the matrix: ";
cin >> DIM;
while(DIM<=0)
  {
    cout << "\nThis must be greater than zero.";
    cout << "\n\nEnter the dimension of the matrix: ";
    cin >> DIM;
  }

cout << "\nEnter the name of the file to open: ";
cin >> fileName;

inFile.open(fileName.c_str(), ios::in);
 
//if can't open
while(!inFile)
  {
    cout << "\n\nCould not open the file '" << fileName << "' ";
    cout << "please try again.";
    inFile.close();
    cout << "\n\nEnter the name of the file to open: ";
    cin >> fileName;
    inFile.open(fileName.c_str(), ios::in);
  }


//get a valid K
cout << "\nEnter the number of random points to sample from the cube: ";
cin >> K;
cout << endl;
while(K<=0)
  {
    cout << "\nThis must be greater than zero.";
    cout << "\n\nEnter the number of strings to create: ";
    cin >> K;
  }


//read in matrix from inFile
for(i=0; i<DIM; i++)
  for(j=0; j<DIM; j++)
    inFile >> matrix[i][j];

//sum rows and columns
for(i=0; i<DIM; i++)
  {
    rowSum[i] = sumRows(matrix, i, DIM);
    colSum[i] = sumCols(matrix, i, DIM);
  }


//check to see if any row/col sum is zero
for(i=0; i<DIM; i++)
  if(rowSum[i]==0)
    {
      cout << "\n\n\tper A = 0";
      cout<<"\n\n**************************************************************";
      cout << "\n\n\n";
      exit(0);
    }

for(i=0; i<DIM; i++)
  if(colSum[i]==0)
    {
      cout << "\n\n\tper A = 0";
      cout<<"\n\n**************************************************************";
      cout << "\n\n\n";
      exit(0);
    }

//calculate m
m=makem(rowSum, colSum, DIM, matrix, mi);

//check to see if m is zero
if(m==0)
  {
    cout << "\n\n\tper A = 1";
    cout<<"\n\n**************************************************************";
    cout << "\n\n\n";
    exit(0);
  }

//calculate the resulting precision (output later)
//float MK = m*K;
//epsilon = sqrt(48/(MK));

M = m;
K1 = K;
alphaaverage =0;

//seed rand()
srand((unsigned)time(NULL));

//for each instance (ie, K times)
for(b=0; b<K; b++)
  {

    //create a random string of length m
    generateStrings(randomStrings, m);

    //assign weights based on RS and phi strings
    assignWeights(matrix, randomStrings, weights, mi, DIM);
    
    //assign huge weight to places with no edge
    for(i=0; i<DIM; i++)
      for(j=0; j<DIM; j++)
	if(matrix[i][j]==0)
	  weights[i][j] = M+1;

 
    for(i=0; i<DIM; i++)
      for(j=0; j<DIM; j++)
	if(matrix[i][j]!=0 && matrix[i][j]!=1)
	  {
	    cerr << "\n\n***Invalid Matrix Input***\n\n\n";
	    exit(1);
	  }

    //lapcost is the optimal weight, computed by LAP
    lapcost = LAP(DIM, weights, rowsol, colsol, u, v);  
    
    alphaaverage = alphaaverage + lapcost;      
  }

//calculate the average
alphaaverage1 = alphaaverage/K1;

 cout << "\n\n\t    m = dimension of Boolean Cube"
      << "\n\t      = " << m;


cout << "\n\n\talpha = estimated average Hamming distance"
     << "\n\t      = " << setprecision(5) << alphaaverage1;

//calculate beta
beta = 0.5 - (alphaaverage1/M);

cout << "\n\n\t Beta = 1/2 - alpha/m" 
     << "\n\t      = " << beta;


//cout << "\n\nThe precision resulting from generating " << K1
//    << " strings is:\n\n\tepsilon = " << epsilon << ".";

if(beta<0)
  {
    if(alphaaverage1>=M+1)
      cout << "\n\n\tper A = 0"; 
    else if(alphaaverage1<M+1)
      cout << "\n\n\tper A = 1";
  }

else
  {
    //calculate lower and upper bounds on log_2(per A)
    lowerbound = (1-H(0.5-beta))*M;
    upperbound = (H(beta))*M;
    
    cout << "\n\n\nThe bounds for log_2(per A) are:";
    
    cout << "\n\n\tlower bound = " << lowerbound;
    cout << "\n\n\tupper bound = " << upperbound;
    
    center = (upperbound + lowerbound)/2;
    cout << "\n\n\t     center = " << center;
  }
 
cout<<"\n\n**************************************************************";
cout << "\n\n\n";

return 0;
}

//*************************************************************************
//*************************************************************************
//*************************************************************************

int sumRows(int matrix[MAXDIM][MAXDIM], int i, int DIM)
{
int j;
int total=0;

for(j=0; j<DIM; j++)
  total = total + matrix[i][j];

return total;
}

//************************************************************************

int sumCols(int matrix[MAXDIM][MAXDIM], int j, int DIM)
{
int i;
int total=0;

for(i=0; i<DIM; i++)
  total = total + matrix[i][j];

return total;
}

//************************************************************************

int makem(int rowSum[MAXDIM], int colSum[MAXDIM], int DIM, 
	  int matrix[MAXDIM][MAXDIM], int mi[MAXDIM])
{
int i, j;
int mplus=0;
int mminus=0;
int m=0;

float temp;

int log2row[DIM];
int log2col[DIM];

int dummy[DIM][DIM];
      
for(i=0; i<DIM; i++)
  {
    temp = log(rowSum[i])/log(2);
    temp = ceil(temp);
    log2row[i] = int(temp);

    temp = log(colSum[i])/log(2);
    temp = ceil(temp);
    log2col[i] = int(temp);
  }

for(i=0; i<DIM; i++)
  {  
    mplus = mplus + log2row[i];
    mminus = mminus + log2col[i];
  }

if(mplus<=mminus)
  {
    if(mminus<0)
      {
	cerr << "\n\n***m+ and m- are both less than zero***\n\n";
	exit(1);
      }
    else if(mplus>=0)
      {
	m = mplus;
	for(i=0; i<DIM; i++)
	  mi[i] = log2row[i];
      }

    else
      {
	m = mminus;
	for(i=0; i<DIM; i++)
	  mi[i] = log2col[i];
      
	for(i=0; i<DIM; i++)
	  for(j=0; j<DIM; j++)
	    dummy[j][i] = matrix[i][j];
	
	for(i=0; i<DIM; i++)
	  for(j=0; j<DIM; j++)
	    matrix[i][j] = dummy[i][j];
      }
  }

else
  {
    if(mplus<0)
      {
	cerr << "\n\n***m+ and m- are both less than zero***\n\n";
	exit(1);
      }

    else if(mminus>=0)
      {
	m = mminus;
	for(i=0; i<DIM; i++)
	  mi[i] = log2col[i];

	for(i=0; i<DIM; i++)
	  for(j=0; j<DIM; j++)
	    dummy[j][i] = matrix[i][j];

	for(i=0; i<DIM; i++)
	  for(j=0; j<DIM; j++)
	    matrix[i][j] = dummy[i][j];
      }

    else
      {
	m = mplus;
	for(i=0; i<DIM; i++)
	  mi[i] = log2row[i];
      }
  }

return m;
}

//************************************************************************

void generateStrings(int randomStrings[MAX], int m)
{
int i;
int temp;

for(i=0; i<m; i++)
  {
    temp = rand()%2;
    randomStrings[i] = temp;
  }

randomStrings[m] = '\0';

return;
}

//************************************************************************

int assignWeights(int matrix[MAXDIM][MAXDIM], int randomStrings[MAX], 
		  int weights[MAXDIM][MAXDIM], int mi[MAXDIM], int DIM)
{
bool done;

int i, j, q;
int edgeWeight=0;
int position=0;
int value=0;
int phiStrings[MAX];



 for(i=0; i<DIM; i++)
   {     
     //initialize phistring to all 0's at beginning of new row
     for(q=0; q<mi[i]; q++)
       phiStrings[q] = 0;

     //calculate position in random string
     position = 0;
     for(q=0; q<i; q++)
       position = position + mi[q];

     for(j=0; j<DIM; j++)
       if(matrix[i][j]==1)
	 {	 
	   //initialize bool variable and weight variable
	   done = false;
	   edgeWeight=0;
	   
	   for(q=0; q<mi[i]; q++)
	     {
	       //compare phiStrings to value by Hamming distance
	       value = randomStrings[position+q];	       
	       
	       if(phiStrings[q]!=value)
		 edgeWeight = edgeWeight + 1;
	     }
	   
	   weights[i][j] = edgeWeight;
	   
	   q = 0;

	   while(!done && q<mi[i])
	     {
	       if(phiStrings[q]==1)
		 phiStrings[q] = 0;
	       
	       else if(phiStrings[q]==0)
		 {
		   phiStrings[q] = 1;
		   done = true;
		 }
	       q++;
	     }
	 }
   }
	 
}

//************************************************************************

int LAP(int dim, int assigncost[MAXDIM][MAXDIM], int rowsol[MAXDIM], 
	int colsol[MAXDIM], int u[MAXDIM], int v[MAXDIM])
{
bool unassignedfound;

int i, imin, numfree=0, prvnumfree, f, i0, k, freerow;
int pred[MAXDIM], free[MAXDIM];
int j, j1, j2, endofpath, last, low, up;
int collist[MAXDIM], matches[MAXDIM];
int min, h, umin, usubmin, v2;
int d[MAXDIM];

 
for(i=0; i<dim; i++)
  matches[i] = 0;

//column reduction
for(j=dim-1; j>=0; j--)//reverse order gives better results
  {
    min = assigncost[0][j];
    imin = 0;
    for(i=1; i<dim; i++)
      if(assigncost[i][j]<min)
	{
	  min = assigncost[i][j];
	  imin = i;
	}
    v[j] = min;

    if(++matches[imin]==1)
      {
	rowsol[imin] = j;
	colsol[j] = imin;
      }
    else
      colsol[j] = -1;
  }

//REDUCTION TRANSFER
for(i=0; i<dim; i++)
  if(matches[i]==0)
    free[numfree++] = i;
  else 
    if(matches[i]==1)
      {
	j1 = rowsol[i];
	min = INT_MAX;
	for(j=0; j<dim; j++)
	  if(j!=j1)
	    if(assigncost[i][j]-v[j] < min)
	      min = assigncost[i][j] - v[j];
	v[j1] = v[j1] - min;
      }
  
//AUGMENGING ROW REDUCTION
int loopcnt = 0;
do
  {
    loopcnt++;
    k=0;
    prvnumfree = numfree;
    numfree = 0;
    while(k<prvnumfree)
      {
	i = free[k];
	k++;
	
	//find min and second min reduced cost over cols
	umin = assigncost[i][0] - v[0];
	j1 = 0;
	usubmin = INT_MAX;
	for(j=1; j<dim; j++)
	  {
	    h = assigncost[i][j] - v[j];
	    if(h<usubmin)
	      if(h>=umin)
		{
		  usubmin = h;
		  j2 = j;
		}
	      else
		{
		  usubmin = umin;
		  umin = h;
		  j2 = j1;
		  j1 = j;
		}
	  }
	
	i0 = colsol[j1];
	if(umin < usubmin)
	  //change the reduction of the min col to increase the min
	  //reduced cost in the row to the subminimum
	  v[j1] = v[j1] - (usubmin - umin);
	else
	  if(i0>=0)
	    {
	      j1 = j2;
	      i0 = colsol[j2];
	    }
	
	//(re)assign i to j1, possibly un-assigning an i0
	rowsol[i] = j1;
	colsol[j1] = i;
	
	if(i0 >= 0)
	  if(umin < usubmin)
	    //put in current k, and go back to that k
	    //continue augmenting path i - j1 with i0
	    free[--k] = i0;
	  else
	    //no further augmenting reduction possible
	    //store i0 in list of free rows for next phase
	    free[numfree++] = i0;
      }
  }
 while(loopcnt < 2);

//AUGMENT SOLUTION for each free row
for(f=0; f<numfree; f++)
  {
    freerow = free[f];

    //Dijkstra shortest path algorithm
    //runs until unassigned column added to shortest path tree
    for(j=0; j<dim; j++)
      {
	d[j] = assigncost[freerow][j] - v[j];
	pred[j] = freerow;
	collist[j] = j;
      }
    
    low = 0;
    up = 0;
    
    unassignedfound = false;
    
    do
      {
	if(up==low)
	  {
	    last = low - 1;
	    
	    //scan cols for up...dim-1 to find indices for which new min occurs
	    //store these indices between low...up-1 (increasing)
	    min = d[collist[up++]];
	    for(k=up; k<dim; k++)
	      {
		j = collist[k];
		h = d[j];
		if(h<=min)
		  {
		    if(h<min)//new min
		      {
			up = low;//restart list at index low
			min = h;
		      }
		    //new index with same min, put on index up, and extend list
		    collist[k] = collist[up];
		    collist[up++] = j;
		  }
	      }
	    //check if any of the min cols happen to be unassigned
	    //if so, we have augmenting path
	    for(k=low; k<up; k++)
	      if(colsol[collist[k]] < 0)
		{
		  endofpath = collist[k];
		  unassignedfound = true;
		  break;
		}
	  }
	
	if(!unassignedfound)
	  {
	    //update 'distances' between freerow and all unscanned cols
	    j1 = collist[low];
	    low++;
	    i = colsol[j1];
	    h = assigncost[i][j1] - v[j1] - min;
	    
	    for(k=up; k<dim; k++)
	      {
		j = collist[k];
		v2 = assigncost[i][j] - v[j] -h;
		if(v2<d[j])
		  {
		    pred[j] = i;
		    if(v2==min)
		      if(colsol[j]<0)
			{
			  endofpath = j;
			  unassignedfound = true;
			  break;
			}
		      else
			{
			  collist[k] = collist[up];
			  collist[up++] = j;
			}
		    d[j] = v2;
		  }
	      }
	  }
      }
    while(!unassignedfound);

    //update column prices
    for(k=0; k<=last; k++)
      {
	j1 = collist[k];
	v[j1] = v[j1] + d[j1] - min;
      }
    
    //reset row and col assignments along alternating path
    do
      {
	i = pred[endofpath];
	colsol[endofpath] = i;
	j1 = endofpath;
	endofpath = rowsol[i];
	rowsol[i] = j1;
      }
    while(i!=freerow);
  }
 
//calculate optimal cost
int lapcost = 0;
for(i=0; i<dim; i++)
  {
    j = rowsol[i];
    u[i] = assigncost[i][j] - v[j];
    lapcost = lapcost + assigncost[i][j];
  }
 
return lapcost;
}
    
//************************************************************************

float H(float x)                          //entropy function
{
float HofX;

if(x!=0)
  HofX = x*((log(1/x))/(log(2))) + (1-x)*((log(1/(1-x)))/(log(2)));

else
  HofX = 0;

return HofX;
}