// **********************************************************************// 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 Gconst 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 FUNCTIONSfloat 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 reductionfor(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 TRANSFERfor(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 REDUCTIONint 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 rowfor(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 usedfloat 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 problemsint 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 distributionfloat 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 distributionfloat 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 embedsit\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()