// **********************************************************************
// author   : Alexander Yong
// email    : ayong@umich.edu
// version  : May 1, 2003
// Approximates the permanent of a nonnegative integer matrix based on
// "Random weighting, asymptotic counting and inverse isoperimetry" by:
// Alexander Barvinok and Alex Samorodnitsky
//
// The LAP code by R. Jonker and A. Volgenant
// email: roy_jonker@majiclogic.com
// with some minor modifications to handle floating point data, nonedges.
// **********************************************************************

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

const int MAXDIM = 300;            // maximal number of vertices of G
const int MAXIMAL_RAND = RAND_MAX;
const float SOLVER_ACCURACY = 0.000001;     
const float PI = 3.14159254;
float nonedge_weight;              
int used_nonedge = 0;              // flag for using a nonedge

//----------------------------------------------------------------------
// LAP FUNCTIONS

float LAP(int dim, float assigncost[MAXDIM][MAXDIM], int rowsol[MAXDIM], 
	int colsol[MAXDIM], float u[MAXDIM], float 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];
float min, umin, h, usubmin, v2;
float  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 and check if nonedge used
float lapcost = 0;
for(i=0; i<dim; i++)
  {
    j = rowsol[i];
    u[i] = assigncost[i][j] - v[j];
    if(assigncost[i][j]==nonedge_weight){
      used_nonedge=1;
    }
    lapcost = lapcost + assigncost[i][j];
  }
 
return lapcost;
}
    
// ---------------------------------------------------------------------
// LOGISTIC DISTRIBUTION FUNCTIONS

// Modified rand(), since returning 0 or MAXIMAL_RAND can cause problems

int modified_rand(){

  int temp;

  temp=rand();
  if(temp==0){
    temp=temp+1;
  }
  else if(temp==MAXIMAL_RAND){
    temp=MAXIMAL_RAND-1;
  }
  return(temp);
}

// Redistributes a point from the uniform distribution to 
// the logistic distribution
float random_weight_logistic(float data_point, int aij){
 
  float log_rand,temp_var,temp_var2;

  temp_var2=-log(data_point)/aij;
  if(temp_var2>0.1){
        temp_var=exp(-log(data_point)/aij);
        temp_var=temp_var-1;
	log_rand=log(temp_var);
	log_rand=-log_rand;
  }
  else{
	log_rand=log(aij)-log(log(1/data_point)) - temp_var2/2-(temp_var2)*(temp_var2)/24;
  }

  if(aij==0){
	return(nonedge_weight);
  }
  else{
   return(-log_rand);
  }
}

// h(t) for the logistic distribution
float H_fn_logistic(float t){
  float Hoft;
  float interval_range;
  float left_point;
  float right_point;
  float mid_point;
  float d1, value1, d2, value2, max_value;

  left_point=0;
  right_point=1-SOLVER_ACCURACY;
  interval_range=2;      
  while(interval_range>SOLVER_ACCURACY){
      mid_point=(left_point+right_point)/2;
      d1=(left_point+mid_point)/2;
      d2=(right_point+mid_point)/2;
      value1=t*d1-log(PI*d1/(sin(PI*d1)));
      value2=t*d2-log(PI*d2/(sin(PI*d2)));
      
		    if(value1>value2){
		      right_point=mid_point;
		      max_value=value1;
		    }
		    else{
		      left_point=mid_point;
		      max_value=value2;
		    }
		    interval_range=right_point-left_point;
  }
  Hoft= max_value;
  return(Hoft);
}

float upper_logistic(float avg_GAMMA, int k){

  float upper_ans;
  
  upper_ans=avg_GAMMA;
  return(upper_ans);
}

float lower_logistic(float avg_GAMMA, int k){

  float lower_ans;
 
  lower_ans=(k-1)*H_fn_logistic((float)avg_GAMMA/(k-1));
  return(lower_ans);
}
 
// ----------------------------------------------------------------------
int main(){
using std::string;
    long int seed = time(0);             // random seed
    int ii,jj,kk;
    int num_vertices;
    string fileName;
    ifstream inFile;
    int matrix[MAXDIM][MAXDIM];                     // [a_ij]
    float weights[MAXDIM][MAXDIM];                  // weight array matrix
    float GAMMA,avg_GAMMA;                         
    float max_weight_matching;
    float upper_bound,lower_bound;
    int m;                                          // num. times to sample
    float max_weight;
    int rowsol[MAXDIM];                          //col matched to row
    int colsol[MAXDIM];                          //row matched to col
    float u[MAXDIM];                               //dual variable
    float v[MAXDIM];                               //dual varuable 
   
    cout << "\n\n\n\n";
    cout << "--------------------------------------------------------------\n";
    cout << "This program takes a given bipartite graph G and embeds
it\n"; 
    cout << "into the\n";
    cout << "hypersimplex. It then estimates the permanent corresponding \n";
    cout << "to G by assigning random weights to each edge of G, and\n";
    cout << "applying the Hungarian algorithm to find a maximal weight perfect\n";
    cout << "matching. This is done m times and the resulting average estimates\n";
    cout << "the function GAMMA. Then upper and lower bounds for the number\n";
    cout << "of the permanent  are calculated.\n";
    cout << "--------------------------------------------------------------\n\n\n";

    cout << "Enter the number of times m to sample:\n";
    cin >> m;
    cout << "\n";
    while(m<=0){
       cout << "This number must be greater than zero.\n";
       cout << "Enter the number of times m to sample: \n";
       cin >> m;
       cout << "\n";
    }

    cout << "Enter the number of vertices in the graph G\n";
    cin >> num_vertices;
    cout << "\n";

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

    inFile.open(fileName.c_str(), ios::in);

    while(!inFile){
       cout << "File open error: '" << fileName << "' ";
       cout <<"Try again.\n";
       inFile.close();
       cout << "Enter the name of the file to open:\n";
       cin >> fileName;
       inFile.open(fileName.c_str(), ios::in);
    }
 
    // read in the matrix from inFile
    for(ii=0; ii<num_vertices; ii++){
      for(jj=0; jj<num_vertices; jj++){
  	   inFile >> matrix[ii][jj];
         }
    }

    // now do the sampling and run LAP m times
    GAMMA=0;
    srand((unsigned)time(NULL));

    for(kk=0;kk<m;kk++){
        max_weight_matching=0;

        // initialize
        for(ii=0;ii<num_vertices;ii++){
           for(jj=0;jj<num_vertices;jj++){
  	        weights[ii][jj]=0;
	   }
	}

	max_weight=0;
	// do the random weighting
	for(ii=0;ii<num_vertices;ii++){
	  for(jj=0;jj<num_vertices;jj++){
	    if(matrix[ii][jj]>0){
              weights[ii][jj]=random_weight_logistic((float)modified_rand()/MAXIMAL_RAND,matrix[ii][jj]);
	      if(weights[ii][jj]>max_weight){
		max_weight=weights[ii][jj];
	      }
            }  
	  }
	}

	// assign large weight to nonedges
	max_weight=max_weight+1;
	nonedge_weight=2*num_vertices*max_weight;
	for(ii=0;ii<num_vertices;ii++){
	  for(jj=0;jj<num_vertices;jj++){
	    if(matrix[ii][jj]==0){
	      weights[ii][jj]=nonedge_weight;
	    }
	  }
	}
    
	// Apply the Hungarian algorithm using the LAP code
        max_weight_matching=-LAP(num_vertices,weights,rowsol,colsol,u,v);
        GAMMA=GAMMA+max_weight_matching;
    }

    avg_GAMMA=(float)GAMMA/m;

    // Compute the upper and lower estimates
    upper_bound=upper_logistic(avg_GAMMA, num_vertices); 
    lower_bound=lower_logistic(avg_GAMMA, num_vertices);

    // Output results
    printf("\n");
    printf("\n");
    printf("-----------------------------------------------------\n");
    if(used_nonedge==1){
      cout<< "No perfect matching was found so the permanent is 0.\n";
    }
    else{
      cout << "Log(X) should be between the following two bounds:\n";
      cout << "Upper bound:\n";
      printf("%f\n",upper_bound);
      cout << "Lower bound:\n";
      printf("%f\n",lower_bound);
      cout << "GAMMA(X):\n";
      printf("%f\n",avg_GAMMA);
    }
}// end of main()