/* ***********************************************************************/
/* 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);
    }
}