/* ***********************************************************************//* author   : Alexander Yong                                             *//* email    : ayong@umich.edu                                            *//* version  : May 5, 2003                                                *//* Code for approximating the hafnian of a nonnegative integer matrix,   *//* based on                                                              *//* "Random weighting, asymptotic counting and inverse isoperimetry" by:  *//* Alexander Barvinok and Alex Samorodnitsky.                            *//* Uses the implementation of H. Gabow's N-cubed weighted matching       *//* algorithm ("Edmonds' algorithm") by Ed Rothberg                       *//* Call this file main.c and "make" it with Rothberg's package (thus     *//* replacing the main.c included there).                                 *//* The package may be found at:                                          *//* ftp://ftp.zib.de/pub/Packages/mathprog/matching/weighted/index.html   *//* ********************************************************************* */#include <stdio.h>#include <stdlib.h>#include <math.h>#include <string.h>#include <time.h>#include "graphtypes.h"#define MAXDIM 300            const int MAXIMAL_RAND=RAND_MAX;const float SOLVER_ACCURACY = 0.000001;     const float PI = 3.14159254;const int INTEGER_SCALE_FACTOR = 100;float nonedge_weight;              int used_nonedge = 0;       /* --------------------------------------------------------------------- *//* LOGISTIC DISTRIBUTION FUNCTIONS *//* Modified rand(), since returning 0 or MAXIMAL_RAND can cause problems */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);  }  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;      /* ensures at least one interation of the following loop: */      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);}main(){    long int seed = time(0);             /* random seed */    int ii,jj,kk,mm,pp;    int k,m,max_entry;    int random_data[MAXDIM][MAXDIM];    float weights[MAXDIM][MAXDIM];                int matrix[MAXDIM][MAXDIM];          float GAMMA,avg_GAMMA;                             float max_weight_matching;    float upper_bound,lower_bound;    float min_weight,max_weight,prob;    float tempa;    int tempb;    int no_matching = 0;    Graph graph;    int *Mate;    int temp_weight;    char file[25];    FILE *fp;    int *input_number;    printf("\n\n\n\n");    printf("--------------------------------------------------------------\n");    printf("This program estimates the hafnian of a given 2k x 2k\n");     printf("symmetric nonnegative integer matrix. This is done by\n");    printf("associating to the matrix an undirected, weight graph\n");    printf("G, assigning random weights to each edge of G, and\n");    printf("applying Edmonds' algorithm to find a maximal weight perfect\n");    printf("matching. This is done m times and the resulting average estimates\n");    printf("the function GAMMA. Then upper and lower bounds for log\n");    printf("of the hafnian are calculated.\n");    printf("--------------------------------------------------------------\n\n\n");    printf("Enter the number of times m>0 to sample: \n");    scanf("%i",&m);    printf("\n");    while(m<=0){       printf("Enter the number of times m>0 to sample: \n");       scanf("%i",&m);       printf("\n");    }    printf("Enter k, half the dimension of the matrix:\n");    scanf("%i",&k);    printf("\n");    printf("Enter name of input file:\n");    scanf("%24s",file);     fp = fopen(file,"r");    for(ii=0;ii<2*k;ii++){       for(jj=0;jj<2*k;jj++){	    fscanf(fp,"%i",&input_number);	    matrix[ii][jj]=input_number;       }    }    fclose(fp);    /* Make matrix symmetric and entries are zero on main diagonal */    for(ii=0; ii<2*k; ii++){      for(jj=0; jj<2*k; jj++){	if(ii==jj){	  matrix[ii][jj]=0;	}  	else if(jj<ii){	  matrix[ii][jj]=matrix[jj][ii];	}      }    }             /* Randomly weight edges and run Edmonds' algorithm m times */    srand((unsigned)time(NULL));    GAMMA=0;    for(kk=0;kk<m;kk++){        /* initialize */        max_weight_matching=0;	max_weight=0;	min_weight=0;        for(ii=0;ii<2*k;ii++){           for(jj=0;jj<2*k;jj++){  	        weights[ii][jj]=0;	   }	}	/* Randomly weight edges */	for(ii=0;ii<2*k;ii++){	  for(jj=0;jj<2*k;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];	      }	      if(weights[ii][jj]<min_weight){		min_weight=weights[ii][jj];	      }            }  	  }	}	/* assign small weight to non-edges */	min_weight=min_weight-1;	nonedge_weight=min_weight;	for(ii=0;ii<2*k;ii++){	  for(jj=0;jj<2*k;jj++){	    if(matrix[ii][jj]==0){	      weights[ii][jj]=nonedge_weight;	    }	  }	}		/* Compose graph data and apply Edmonds' algorithm */	graph = NewGraph(2*k);	for(mm=0; mm<2*k;mm++){	  NLabel(graph,mm+1)=0;	  Xcoord(graph,mm+1)=0;	  Ycoord(graph,mm+1)=0;	  for(pp=mm+1;pp<2*k;pp++){	    if(matrix[mm][pp]>0){	      /* The external Edmonds' algorithm routines          */  	      /* find a MAXIMAL weight perfect matching assuming nonnegative */	      /* integer weights, so we adjust weights accordingly */	      temp_weight=INTEGER_SCALE_FACTOR*(1-min_weight+weights[mm][pp]);	      AddEdge(graph,mm+1,pp+1,temp_weight);	    }  	  }	}         	Mate = Weighted_Match(graph,1);	for(pp=1;pp<=2*k;pp++){	  if(Mate[pp]==0){	    no_matching=1;	  }	  if(Mate[pp]>pp){	    max_weight_matching = max_weight_matching+weights[pp-1][Mate[pp]-1]; 	    if(weights[pp-1][Mate[pp]-1]==nonedge_weight){	      used_nonedge=1;	    }	  }	}        GAMMA=GAMMA+max_weight_matching;    }    avg_GAMMA=(float)GAMMA/m;    /* Compute the upper and lower estimates */    upper_bound=upper_logistic(avg_GAMMA, k);     lower_bound=lower_logistic(avg_GAMMA, k);    /* Output results */    printf("\n");    printf("\n");    printf("-----------------------------------------------------\n");    if(no_matching==1){      printf("No perfect matching was found: no results returned.\n");    }    else{      printf("Log(X) should be between the following two bounds:\n");      printf("Upper bound:\n");      printf("%f\n",upper_bound);      printf("Lower bound:\n");      printf("%f\n",lower_bound);      printf("GAMMA(X):\n");      printf("%f\n",avg_GAMMA);    }}