// **********************************************************************
// author   : Alexander Yong
// email    : ayong@umich.edu
// version  : April 30, 2003
// 
// Code for approximating number of spanning forests of a graph, based on
// "Random weighting, asymptotic counting and inverse isoperimetry" by:
// Alexander Barvinok and Alex Samorodnitsky.
//
// Standard Kruskal algorithm implementation.
// **********************************************************************

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

const int MAXDIM = 300;            // maximal number of vertices of G
const int MAXEDGES = 10000;        // make at least MAXDIM choose 2
const int MAXIMAL_RAND = RAND_MAX;
const float BIG_NEGATIVE_FOR_NONEDGE = -6789;   
const float SOLVER_ACCURACY = 0.000001;      
const float PI = 3.14159254;

// ----------------------------------------------------------------------
// KRUSKAL'S ALGORITHM FUNCTIONS

int U[MAXDIM];        
int used_nonedge = 0;  

float search(float i){  
  
   int j;

   j = (int)i;
   while (U[j] != j)
      j = U[j];
   return (j);
}

float kruskal(float W[MAXDIM][MAXDIM], int N, int num_edges_forest, int M){
  
   float E[MAXEDGES][3];            // complete set of edges           
   float F[MAXEDGES][3];            // set of edges in max span. forest
   int num_edges = 0;               
   int next_edge = 0;               
   float weight = 0;                
   int i, j, k, q;                     
   float a, b, c;

   // initialize set of edges 
   k = 0;
   for (i = 0; i < N; i++){
     for (j = 0; j < N; j++){
         if (j > i){  
	   E[k][0] = i;          // first vertex of edge 
	   E[k][1] = j;          // second vertex of edge 
	   E[k][2] = W[i][j];    // weight of edge 
           k=k+1;
         }
     }
   }
   
   // bubblesort set of edges in non-increasing order by weight 
   for (i = M - 1; i > 0; i--){
      for (j = 0; j < i; j++){
	 if (E[j][2] < E[j+1][2]){
            a = E[j][0];
            b = E[j][1];
            c = E[j][2];
            E[j][0] = E[j+1][0];
            E[j][1] = E[j+1][1];
            E[j][2] = E[j+1][2];
            E[j+1][0] = a;
            E[j+1][1] = b;
            E[j+1][2] = c;   
         }
     }
   }

   // create N disjoint subsets 
   for(q=0;q<N;q++){
     U[q]=q;
   }

   // initialize set of edges in max. span. forest to empty 
   for (i = 0; i < N - 1; i++){
     for (j = 0; j < 3; j++){
	 F[i][j] = -1;            // -1 denotes empty 
     }
   }

   // find maximal spanning forest
   while (num_edges < num_edges_forest){  
      a = E[next_edge][0];
      b = E[next_edge][1];      
      i = (int)search(a);
      j = (int)search(b);
      
      if(i!=j){
	 if(i<j){
	   U[j]=i;
	 }
	 else{
	   U[i]=j;
	 }
         F[num_edges][0] = E[next_edge][0];
         F[num_edges][1] = E[next_edge][1];
         F[num_edges][2] = E[next_edge][2];
         num_edges=num_edges+1;
      }
      
      next_edge=next_edge+1;
   }
   
   // compute maximal weight and check for use of nonedges
   for (i = 0; i < num_edges_forest; i++){  
      if(F[i][2]!=BIG_NEGATIVE_FOR_NONEDGE){
	weight = weight + F[i][2];
      }
      else{
	used_nonedge=1;
	weight=0;
      }
   }
   return(weight);
}

// Modification of rand(): since problems occur when
// rand() returns 0 or MAXIMAL_RAND
int modified_rand(){

  int temp;

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

// ---------------------------------------------------------------------
// EXPONENTIAL DISTRIBUTION FUNCTIONS

// Redistribute a point from the uniform distribution
// to the exponential distribution
float random_weight_exp(float x){
 
  float exp_rand;

  // now get the exponential distribution
  if(x<=.5){
    exp_rand=(float)log(2*x);
  }
  else{
    exp_rand=-(float)log(2-2*x);
  }
  return(exp_rand);
}

// h(t) for the exponential distribution
float H_fn_exp(float t){
  float Hoft;

  Hoft= sqrt(1+t*t)+log(sqrt(1+t*t)-1)-2*log(t)+log(2)-1;
  return(Hoft);
}

float upper_exp(float avg_GAMMA, int k){

  float upper_ans;
  
  upper_ans=avg_GAMMA+log(2)*k;
  return(upper_ans);
}

float lower_exp(float avg_GAMMA, int k){

  float lower_ans;
 
  lower_ans=(k-1)*H_fn_exp((float)avg_GAMMA/(k-1));
  return(lower_ans);
}

// -----------------------------------------------------------------------
// LOGISTIC DISTRIBUTION FUNCTIONS

// Redistribute a point from the uniform distribution
// to the logistic distribution
float random_weight_logistic(float x){
 
  float log_rand;

  log_rand=log(x)-log(1-x);
  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, num_edges_forest;
    string fileName;
    ifstream inFile;
    int matrix[MAXDIM][MAXDIM];          // vertex 0-1 incidence matrix
    float weights[MAXDIM][MAXDIM];      
    float GAMMA, avg_GAMMA;                         
    float max_weight_forest;
    float upper_bound,lower_bound;
    float prob;
    int m;                               // number of times to sample       
    int complete_num_edges;              // equals (num_vertices choose 2)
    int dist_type;                       // which distribution to use 

    cout << "--------------------------------------------------------------\n";
    cout << "This program takes a given graph G and embeds it into the\n";
    cout << "hypersimplex. It then estimates the number of spanning forests\n";
    cout << "of G by assigning random weights to each edge of G, and\n";
    cout << "applying Kruskal's algorithm to find a maximal weight spanning\n";
    cout << "forest. 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 spanning forests 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 number of edges in the forest.\n";
    cout << "This must be less than the number of vertices.\n";
    cin >> num_edges_forest;

    cout << "Enter the name of the input file:\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 input file:\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];
	}
    }

    cout << "Which distribution?:\n";
    cout << "1. Exponential\n";
    cout << "2. Logistic\n";
    cout << "Enter 1 or 2\n";
    cin >> dist_type;

    while(dist_type!=1 && dist_type!=2){
	cout << "Please choose distribution 1 or 2.\n";
	cin >> dist_type;
    }

    complete_num_edges=(int)(num_vertices)*(num_vertices-1)/2;

    // now sample weights and run Kruskal m times
    GAMMA=0;
    srand((unsigned)time(NULL));
    for(kk=0;kk<m;kk++){
        max_weight_forest=0;

	// randomly weight the edges
        for(ii=0;ii<num_vertices;ii++){
           for(jj=0;jj<num_vertices;jj++){
  	        weights[ii][jj]=0;
	   }
	}
	for(ii=0;ii<num_vertices;ii++){
	  for(jj=0;jj<num_vertices;jj++){
	    if(matrix[ii][jj]==1){
              if(dist_type==1){
	            weights[ii][jj]=random_weight_exp((float)modified_rand()/MAXIMAL_RAND);
              }
              else if(dist_type==2){
                    weights[ii][jj]=random_weight_logistic((float)modified_rand()/MAXIMAL_RAND);
              }
	    }
	    else{
	      weights[ii][jj]=BIG_NEGATIVE_FOR_NONEDGE;  // nonedges get a 
	                                                 // very negative weight
	    }
	  }
	}
    	  
	// Apply Kruskal's algorithm to find the maximal weight 
	// of a spanning forest. Check that nonedges are never used.
     
        max_weight_forest=kruskal(weights,num_vertices,num_edges_forest,complete_num_edges);
        GAMMA=GAMMA+max_weight_forest;
    }

    // Compute the upper and lower estimates
    avg_GAMMA=(float)GAMMA/m;
    if(dist_type==1){
      upper_bound=upper_exp(avg_GAMMA, num_vertices); 
      lower_bound=lower_exp(avg_GAMMA, num_vertices);
    }
    else if(dist_type==2){
      upper_bound=upper_logistic(avg_GAMMA, num_vertices); 
      lower_bound=lower_logistic(avg_GAMMA, num_vertices);
    }

    // Output results
    cout<< "\n";
    cout<< "\n";
    cout<< "-----------------------------------------------------\n";
    if(used_nonedge==1){
      cout<< "The graph was not connected, so no results returned\n";
    }
    else{
      cout << "Bounds for Log(X) are:\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);
    }

}