/***********************************************************************/
/* The following code computes alt. paths for metabolic conversions via*/
/* Johnson's algorithm for enumerating cycles in a directed graph.     */
/* This program requires one input file containing information on the  */
/* directed graph we are analyzing: the first line contains the total  */
/* number of vertices separated by a space from the total number of    */
/* edges of the graph. the other lines specify                         */
/* the edges b.m.o. the starting vertex (first column), the ending     */
/* vertex (second column) and the weight of the edge (third column).   */
/* The output consists in the cycle distribution of the graph and      */
/* various (1-,2-,3-point) marginal distributions derived form this.   */
/***********************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef unsigned long long RANDOM_TYPE;

#define NUM_VERTICES_MAX 200000
#define NUM_EDGES_MAX  (NUM_VERTICES_MAX*3)
#define B_SPACE 10000
#define MAX_STACK_DEPTH B_SPACE
/* NUM_EDGES is the number of oriented edges            */
/*           the number of oriented edges is 2*NUM_EDGES  */

#define MYRAND64A  myrand64a = (6364136223846793005UL * myrand64a)
#define MASSIMO_RAND64     0XFFFFFFFFFFFFFFFFUL
RANDOM_TYPE myrand64a = 578991;
double inv_massimo_rand64 = 1.0L / (double)MASSIMO_RAND64;

typedef int MY_VERTEX;
typedef int MY_INT;

FILE *f_ver;

char filename[200];

double averageNumEdges;

long NUM_VERTICES, numVerticesRead, NUM_EDGES, i_sample;
long lengthOne = 0;
int NUM_SAMPLES, connectivity;
float fconnectivity;
unsigned long long num_circuits;

MY_VERTEX num_vertices_now, num_edges_now, max_vertex;
MY_VERTEX vertices_of_edge_list0[NUM_EDGES_MAX][2];
MY_VERTEX vertices_of_edge_list[NUM_EDGES_MAX][2];
MY_VERTEX tedge_first_of_vertex[NUM_VERTICES_MAX];
MY_VERTEX tedge_last_of_vertex[NUM_VERTICES_MAX];

MY_VERTEX vertices_of_b_list[B_SPACE];
MY_VERTEX tvertex_b_next_of_tvertex[B_SPACE];
MY_VERTEX tvertex_b_former_of_tvertex[B_SPACE];
MY_VERTEX tvertex_b_first_of_vertex[NUM_VERTICES_MAX];
MY_VERTEX tvertex_b_last_of_vertex[NUM_VERTICES_MAX];
MY_VERTEX tvertex_b_first_of_empty, tvertex_b_last_of_empty;
MY_VERTEX num_b_empty, num_b[B_SPACE], num_tot_b;

MY_VERTEX startVertexConnection[NUM_VERTICES_MAX];
MY_VERTEX endVertexConnection[NUM_VERTICES_MAX];
MY_VERTEX numberVertexConnection[NUM_VERTICES_MAX];
MY_VERTEX maxOccupied;

MY_VERTEX my_stack[MAX_STACK_DEPTH], stack_depth;

unsigned long long vertex_marginal[NUM_VERTICES_MAX], edge_marginal[NUM_EDGES_MAX];
int adjacency[100][100];
int vertices_of_original_edges[2000][2];
long double edge_weight[2000];
long double tot_weighted_cycles;
long double weighted_vertex_marginal[1000], weighted_edge_marginal[2000];
long double edge_adjacency_marginal[2000][2000];

char blocked[NUM_VERTICES_MAX];

MY_VERTEX this_circuit_vertex;

unsigned long long his_circuits[NUM_VERTICES_MAX+2];

#define EMPTY_VERTEX -1
#define MY_FIRST -2
#define MY_LAST -3
#define EMPTY_TVERTEX -4
#define NO_VERTEX -5
#define MY_FALSE 0
#define MY_TRUE 1

#define MY_DEBUG
#undef MY_DEBUG
#define MY_DEBUG2
#undef MY_DEBUG2
#define PRUNING
#undef PRUNING
#define WRITE_GRAPH
#undef WRITE_GRAPH
#define LOCAL_HYS
#undef LOCAL_HYS
#define MY_TEST_RANDOM
#undef MY_TEST_RANDOM
#define FIX_RAN
#undef FIX_RAN
#define ER_P_RAN
#undef ER_P_RAN
#define ER_N_RAN
#undef ER_N_RAN

#define PRUNING
#define ER_N_RAN

#ifdef LOCAL_HYS
unsigned long long numberCircuits[NUM_VERTICES_MAX+2];
#endif

void read_params(int argc, char *argv[]);
void my_init(void);
void my_end(void);
void my_init_inside(void);
void my_init_ran(void);
void my_printout(void);
void chop_vertices(MY_VERTEX v_min);
void read_vertices_from_edges_dat(void);
void pruneItUp(void);
int check_and_set_vertices_of_edge_list(MY_VERTEX v_min);
void init_b_list_to_empty(void);
void print_b_list(void);
void add_to_b_list(MY_VERTEX v1, MY_VERTEX v2);
void erase_from_b_list(MY_VERTEX v1, MY_VERTEX v2);
void test_add_erase_b(void);
void init_stack(void);
void print_stack(void);
void print_stack_debug(void);
void print_stack_gt_two_plus(MY_VERTEX v);
void push_stack(MY_VERTEX v);
MY_VERTEX pop_stack(void);
void find_circuits(void);

void unblock(MY_VERTEX u);
char circuit(MY_VERTEX v);

void my_init_marginals(void);
void adjust_marginal(MY_VERTEX v);
void compute_ranking(void);
/**********************************************************************/
int main(int argc, char **argv ){

  read_params(argc, argv);

  my_init();

  for(i_sample=0; i_sample<NUM_SAMPLES; i_sample++){
    my_init_inside();
    my_init_marginals();

    read_vertices_from_edges_dat();

    find_circuits();
    my_printout();
  }
  my_end();
 
  return(0);
}

/***********************************************************************/
void read_params(int argc, char *argv[])
{
  int l;

  if (argc>1)
    {
      for (l=0;argv[1][l]!='\0';l++)    
	filename[l]=argv[1][l];
      filename[l]='\0';
      fflush(stdout);
    }
  else
    {
      printf("# Usage: <name>\n");
      exit(0);
    }
}

/***********************************************************************/
void my_end(void){
  //printf("AVERAGE NUM (UNORIENTED) EDGES %lg\n", 
  //averageNumEdges/(double)NUM_SAMPLES);
}

/***********************************************************************/
void my_init_ran(void){
  printf("# Please input num vertices\n"); fflush(stdout);
  scanf("%ld",&numVerticesRead);
  printf("# num vertices is %ld\n",numVerticesRead);
  fflush(stdout);
  NUM_VERTICES = numVerticesRead;

  printf("# Please input num connectivity (long int)\n"); fflush(stdout);
  scanf("%d",&connectivity);
  printf("# connectivity is %d\n",connectivity);

  printf("# Please input num connectivity (float)\n"); fflush(stdout);
  scanf("%f",&fconnectivity);
  printf("# connectivity is %f\n",fconnectivity);
  fflush(stdout);
}

/***********************************************************************/
void my_init(void){

  NUM_SAMPLES=1;
  averageNumEdges = 0.0;
}


/***********************************************************************/
void my_init_inside(void){
  int i;

  for(i=0; i<NUM_VERTICES_MAX+2; i++){
    his_circuits[i]   = 0;
#ifdef LOCAL_HYS
    numberCircuits[i] = 0;
#endif
  }
  num_circuits = 0ULL;
}
/***********************************************************************/
void my_init_marginals(void){
  int i,j;
  
  tot_weighted_cycles=0;
  for(i=0; i<NUM_VERTICES_MAX; i++){
    vertex_marginal[i]   = 0;
  }
  for(i=0; i<NUM_EDGES_MAX; i++){
    edge_marginal[i]   = 0;
  }
  for(i=0; i<1000; i++){
    weighted_vertex_marginal[i]   = 0;
  }
  for(i=0; i<2000; i++){
    weighted_edge_marginal[i]   = 0;
    for(j=0; j<2000; j++){
      edge_adjacency_marginal[i][j]=0;
    }
  }
}

/***********************************************************************/
void my_printout(void){
  unsigned long long tot_cyc;
  unsigned long long t2=0;
  int i,j;
  long double marginal;

  tot_cyc=0;
  printf ("# [ L 2 ]= %llu\n", his_circuits[2]);
  for(i=3; i<NUM_VERTICES_MAX+2; i++){
    if(his_circuits[i]>0){
      t2+=his_circuits[i];
      printf ("# [ L %d ]= %llu\n", i, his_circuits[i]);
    }
  }

  for(i=2; i<NUM_VERTICES_MAX+2; i++){
    if(his_circuits[i]>0){
      fflush(stdout);
      tot_cyc+=his_circuits[i];
    }
  }
  printf("# TOT %lld %lld\t%lld %3.5Le\n", (t2+his_circuits[2]), t2, tot_cyc, tot_weighted_cycles);

#ifdef LOCAL_HYS
  for(i=0; i<NUM_VERTICES_MAX+2; i++){
    if(numberCircuits[i]>0){
      printf("LOCAL S %d I %d N %llu \n",
	     i_sample, i, numberCircuits[i]);
    }
    fflush(stdout);
  }
#endif

  for(i=0; i<NUM_VERTICES; i++){
    marginal=(long double)(vertex_marginal[i])/(long double)(t2+1);
    weighted_vertex_marginal[i]/=(long double)(tot_weighted_cycles+1.);
  }
  for(i=0; i<NUM_EDGES; i++){
    marginal=(long double)(edge_marginal[i])/(long double)(t2+1);
    weighted_edge_marginal[i]/=(long double)(tot_weighted_cycles+1.);
  }

  for(i=0; i<NUM_EDGES; i++){
    for(j=0; j<NUM_EDGES; j++){
      edge_adjacency_marginal[i][j]/=(long double)(tot_weighted_cycles+1.);
    }
  }

  compute_ranking();

}


/***********************************************************************/
void find_circuits(void){
  MY_VERTEX vw;
  char dummy;
  
  init_stack();
  for(this_circuit_vertex = 0; this_circuit_vertex<NUM_VERTICES-1;
      this_circuit_vertex++){
#ifdef MY_DEBUG
    printf("DBG FC this_circuit_vertex %d\n", this_circuit_vertex);
#endif
    if(check_and_set_vertices_of_edge_list(this_circuit_vertex)==0){
      break;
    }
    if((tedge_first_of_vertex[this_circuit_vertex]==EMPTY_TVERTEX) &&
       (tedge_last_of_vertex[this_circuit_vertex]==EMPTY_TVERTEX)){
      continue;
    }
#ifdef MY_DEBUG
    printf("DBG FC tedge first %d last %d\n", 
	   tedge_first_of_vertex[this_circuit_vertex],
	   tedge_last_of_vertex[this_circuit_vertex]);
#endif
    init_b_list_to_empty();
    for(vw=this_circuit_vertex; vw<NUM_VERTICES; vw++){
      blocked[vw]=(char)MY_FALSE;
#ifdef MY_DEBUG
      printf("DBG FC setting blocked[%d] to false\n",vw); 
#endif
    }
#if ( (defined MY_DEBUG) ||  (defined MY_DEBUG2 ) )
    printf("DBG FC calling circuit(%d)\n",this_circuit_vertex); 
    fflush(stdout);
#endif
    dummy = circuit(this_circuit_vertex);
  }
}

/***********************************************************************/
char circuit(MY_VERTEX v){
  MY_VERTEX tw, w, vscratch;
  char circuit_flag;
  long i;

  circuit_flag = (char)MY_FALSE;
  push_stack(v);
  blocked[v] = (char)MY_TRUE;

#ifdef MY_DEBUG
  printf("DBG CIR ENTERING with v = %d that I PUSH in stack and BLOCK \n",v);
#endif
  if ((tedge_first_of_vertex[v]!=EMPTY_TVERTEX) && (tedge_last_of_vertex[v]!=EMPTY_TVERTEX)){
    for(tw=tedge_first_of_vertex[v]; tw<=tedge_last_of_vertex[v]; tw++){
      w = vertices_of_edge_list[tw][1];
#ifdef MY_DEBUG
      printf("DBG CIR   for w in Ak tw %d starting %d final %d\n",
	     tw,vertices_of_edge_list[tw][0],w);
#endif
      if(w==this_circuit_vertex){
#ifdef MY_DEBUG
	printf("DBG CIR   I AM BACK SO I GOT IT tw %d\n",tw);
#endif
	adjust_marginal(this_circuit_vertex);
	//print_stack_gt_two_plus(this_circuit_vertex);
	//print_stack();
	num_circuits++;
	his_circuits[stack_depth]++;
	/* QUI LOCAL */
	for(i=0; i<stack_depth; i++){
#ifdef LOCAL_HYS
	  numberCircuits[my_stack[i]]++;
#endif
	}
	circuit_flag = (char)MY_TRUE;
      }else if(blocked[w]==(char)MY_FALSE){
#ifdef MY_DEBUG
	printf("DBG CIR    since blocked[%d] is FALSE I call circuit(%d)\n",w,w);
#endif
	if(circuit(w)==(char)MY_TRUE){
	  circuit_flag=(char)MY_TRUE;
	}
      }
#ifdef MY_DEBUG
      else {
	printf("DBG CIR    since blocked[%d] is TRUE I did not call circuit(%d)\n",w,w);
      }
#endif
    }
  }
  if(circuit_flag==(char)MY_TRUE){
#ifdef MY_DEBUG
    printf("DBG CIR     CIRCUIT_FLAG IS TRUE so I unblock(%d)\n",v);
#endif
    unblock(v);
  } else{
#ifdef MY_DEBUG
    printf("DBG CIR      CIRCUIT_FLAG IS FALSE\n");
#endif
    for(tw=tedge_first_of_vertex[v]; tw<=tedge_last_of_vertex[v]; tw++){
      w = vertices_of_edge_list[tw][1];
      add_to_b_list(w,v);
#ifdef MY_DEBUG
      printf("DBG CIR        for tw %d added v %d to w %d\n",tw,v,w);
#endif
    }
  }
  vscratch = pop_stack();
  if(v!=vscratch){
    printf("abort internal error in circuit, wrong w popped %d %d\n",
	   v,vscratch);
    exit(-9);
  }

#ifdef MY_DEBUG
  printf("DBG CIR  vertex %d returning away to caller with cf = %d\n",
	 v,circuit_flag);
  print_stack_debug();
#endif
  return circuit_flag;
}

/***********************************************************************/
void unblock(MY_VERTEX u){
  int si, sos_n;
  MY_VERTEX sos_v[B_SPACE];
  MY_VERTEX tw, w, tfirst, tlast;

  blocked[u] = (char)MY_FALSE;

  tfirst = tvertex_b_first_of_vertex[u];
  tlast = tvertex_b_last_of_vertex[u];
#ifdef MY_DEBUG
  printf("DBG UNB unblocking u = %d tfirst = %d tlast = %d\n",u,tfirst,tlast);
#endif  
  if(tfirst!=EMPTY_TVERTEX){
    sos_n=0;
    for(tw=tfirst;;tw=tvertex_b_next_of_tvertex[tw]){
      w = vertices_of_b_list[tw];
      if(sos_n>=B_SPACE){
	printf("abort unblock sos_n too large\n");
	exit(-9);
      }
      sos_v[sos_n]=w;
      sos_n++;
      if(tw==tlast)break;
    }
    for(si=0;si<sos_n;si++){
      w = sos_v[si];
      erase_from_b_list(u,w);
      if(blocked[w]==(char)MY_TRUE){unblock(w);}
    }
  }
#ifdef MY_DEBUG
  printf("DBG UNB returning from happy unblock\n");
#endif
}

/***********************************************************************/
MY_VERTEX pop_stack(void){
  stack_depth--;
  return my_stack[stack_depth];
}

/***********************************************************************/
void push_stack(MY_VERTEX v){
  my_stack[stack_depth] = v;
  stack_depth++;
}

/***********************************************************************/
void print_stack(void){
  long i;
  printf("  ========= ******** ======== ");
  for(i=0; ; i++){
    if(i==stack_depth){printf("  =========== ********** ======\n"); return;}
    printf("%d ",my_stack[i]);
  }
}

/***********************************************************************/
void print_stack_debug(void){
  long i;
  printf("DBG PRSD stack_depth is %d ___ ",stack_depth);
  for(i=0; ; i++){
    if(i==stack_depth){printf("\n");return;}
    printf("%d ",my_stack[i]);
  }
}

/***********************************************************************/
void print_stack_gt_two_plus(MY_VERTEX v){
  long i;
  int l, found;
  //printf("DBG PRSGTP stack_depth is %d\n",stack_depth);
  if(stack_depth<2)return;
  printf("L %d\n", stack_depth);
  for(i=0; ; i++){
    if(i==stack_depth){printf("%d \n",v); return;}
    printf("%d ",my_stack[i]);
    found=0;
    for (l=0;l<NUM_EDGES;l++){
      if ((vertices_of_edge_list0[l][0]==my_stack[i])&&(vertices_of_edge_list0[l][1]==my_stack[(i+1)%stack_depth])){
	found=1;
      }
    }
    if (found==0){
      printf("ERROR: I didn't find the edge %d-%d.\n", my_stack[i], my_stack[i+1]);
      exit(-9);
    }
  }
}
/***********************************************************************/
void adjust_marginal(MY_VERTEX v){

  long i;
  int edge, previous_vertex, previous_edge;
  long double weight_cycle;

  weight_cycle=1;

  if (stack_depth>2){
    tot_weighted_cycles+=weight_cycle;
    for(i=0; ; i++){
      if(i==stack_depth){
	edge=adjacency[previous_vertex][v];
	edge_marginal[edge]++;      
	weighted_edge_marginal[edge]+=weight_cycle;      
	edge_adjacency_marginal[previous_edge][edge]+=weight_cycle;
	return;
      }
      vertex_marginal[my_stack[i]]++;
      weighted_vertex_marginal[my_stack[i]]+=weight_cycle;
      if (i!=0){
	edge=adjacency[previous_vertex][my_stack[i]];
	if (i!=1){
	  edge_adjacency_marginal[previous_edge][edge]+=weight_cycle;
	}
	else {
	  previous_edge=adjacency[my_stack[stack_depth-1]][v];
	  edge_adjacency_marginal[previous_edge][edge]+=weight_cycle;
	}
	edge_marginal[edge]++;      
	weighted_edge_marginal[edge]+=weight_cycle;      
	previous_edge=edge;
      }
      previous_vertex=my_stack[i];
    }
  }
  else {
    return;
  }
}

/***********************************************************************/
void compute_ranking()
{
  int i, j, temp_rank, nb_vertex_rank;
  int temp_edge1, temp_edge2;
  long double temp_marg;
  int *vertex_ranking, *edge_ranking, *edge_adjacency_ranking;
  int nb_edge_adjacency;
  struct {
    long double marginal;
    int edge1;
    int edge2;
  } *weighted_edge_adjacency_marginal;

  vertex_ranking=calloc(NUM_VERTICES,sizeof(int));
  edge_ranking=calloc(NUM_EDGES,sizeof(int));

  nb_vertex_rank=0;
  for(i=0;i<NUM_VERTICES;i++)
    vertex_ranking[i]=i;
  for(i=0; i<(NUM_VERTICES-1); i++)
    {
      for (j=(NUM_VERTICES-1); j>i; j--)
	{
	  if (weighted_vertex_marginal[j-1]<weighted_vertex_marginal[j])
	    {
	      temp_rank=vertex_ranking[j-1];
	      vertex_ranking[j-1]=vertex_ranking[j];
	      vertex_ranking[j]=temp_rank;	
	      temp_marg=weighted_vertex_marginal[j-1];
	      weighted_vertex_marginal[j-1]=weighted_vertex_marginal[j];
	      weighted_vertex_marginal[j]=temp_marg;
	    }
	}
    }

  for(i=0;i<NUM_EDGES;i++)
    edge_ranking[i]=i;

  for(i=0; i<(NUM_EDGES-1); i++)
    {
      for (j=(NUM_EDGES-1); j>i; j--)
	{
	  if (weighted_edge_marginal[j-1]<weighted_edge_marginal[j])
	    {
	      temp_rank=edge_ranking[j-1];
	      edge_ranking[j-1]=edge_ranking[j];
	      edge_ranking[j]=temp_rank;	
	      temp_marg=weighted_edge_marginal[j-1];
	      weighted_edge_marginal[j-1]=weighted_edge_marginal[j];
	      weighted_edge_marginal[j]=temp_marg;
	    }
	}
    }

  nb_edge_adjacency=0;
  for(i=0;i<NUM_EDGES;i++){
    for(j=0;j<NUM_EDGES;j++){
      if (edge_adjacency_marginal[i][j]>0.){
	nb_edge_adjacency++;
      }
    }
    if (edge_adjacency_marginal[i][i]>0.){
      printf("ERROR: edge_adjacency_marginal[%d][%d] %3.7Le != 0\n", i, i, edge_adjacency_marginal[i][i]);
      exit(-9);
    }
  }
  weighted_edge_adjacency_marginal=calloc(nb_edge_adjacency,(sizeof(long double)+2*sizeof(int)));
  nb_edge_adjacency=0;
  for(i=0;i<NUM_EDGES;i++)
    for(j=0;j<NUM_EDGES;j++)
      if (edge_adjacency_marginal[i][j]>0.){
	weighted_edge_adjacency_marginal[nb_edge_adjacency].marginal=edge_adjacency_marginal[i][j];
	weighted_edge_adjacency_marginal[nb_edge_adjacency].edge1=i;
	weighted_edge_adjacency_marginal[nb_edge_adjacency].edge2=j;
	nb_edge_adjacency++;
      }
  edge_adjacency_ranking=calloc(nb_edge_adjacency,sizeof(int));
  for(i=0;i<nb_edge_adjacency;i++){
    edge_adjacency_ranking[i]=i;
  }
  for(i=0; i<(nb_edge_adjacency-1); i++){
    for (j=(nb_edge_adjacency-1); j>i; j--){
      if (weighted_edge_adjacency_marginal[j-1].marginal<weighted_edge_adjacency_marginal[j].marginal){
	temp_rank=edge_adjacency_ranking[j-1];
	edge_adjacency_ranking[j-1]=edge_adjacency_ranking[j];
	edge_adjacency_ranking[j]=temp_rank;	
	temp_marg=weighted_edge_adjacency_marginal[j-1].marginal;
	weighted_edge_adjacency_marginal[j-1].marginal=weighted_edge_adjacency_marginal[j].marginal;
	weighted_edge_adjacency_marginal[j].marginal=temp_marg;
	temp_edge1=weighted_edge_adjacency_marginal[j-1].edge1;
	weighted_edge_adjacency_marginal[j-1].edge1=weighted_edge_adjacency_marginal[j].edge1;
	weighted_edge_adjacency_marginal[j].edge1=temp_edge1;
	temp_edge2=weighted_edge_adjacency_marginal[j-1].edge2;
	weighted_edge_adjacency_marginal[j-1].edge2=weighted_edge_adjacency_marginal[j].edge2;
	weighted_edge_adjacency_marginal[j].edge2=temp_edge2;
      }
    }
  }

  printf("# Rank\t[V-V-V]\t\tE_adjacency_marginal\n");
  for(i=0;i<nb_edge_adjacency;i++){
    temp_edge1=weighted_edge_adjacency_marginal[i].edge1;
    temp_edge2=weighted_edge_adjacency_marginal[i].edge2;
    if (vertices_of_original_edges[temp_edge1][0]==vertices_of_original_edges[temp_edge2][0]){
      printf("%d %d %d %d", vertices_of_original_edges[temp_edge1][1], vertices_of_original_edges[temp_edge1][0], vertices_of_original_edges[temp_edge2][0], vertices_of_original_edges[temp_edge2][1]);
    }
    else if (vertices_of_original_edges[temp_edge1][0]==vertices_of_original_edges[temp_edge2][1]){
      printf("%d %d %d %d ", vertices_of_original_edges[temp_edge1][1], vertices_of_original_edges[temp_edge1][0], vertices_of_original_edges[temp_edge2][1], vertices_of_original_edges[temp_edge2][0]);
    }
    else if (vertices_of_original_edges[temp_edge1][1]==vertices_of_original_edges[temp_edge2][0]){
      printf("%d %d %d %d", vertices_of_original_edges[temp_edge1][0], vertices_of_original_edges[temp_edge1][1], vertices_of_original_edges[temp_edge2][0], vertices_of_original_edges[temp_edge2][1]);
    }
    else if (vertices_of_original_edges[temp_edge1][1]==vertices_of_original_edges[temp_edge2][1]){
      printf("%d %d %d %d", vertices_of_original_edges[temp_edge1][0], vertices_of_original_edges[temp_edge1][1], vertices_of_original_edges[temp_edge2][1], vertices_of_original_edges[temp_edge2][0]);
    }
    else {
      printf("ERROR:no match found between the vertices of edge %d and %d\n", temp_edge1, temp_edge2);
      exit(-9);
    }
    printf("\t%3.7Le\n", weighted_edge_adjacency_marginal[i].marginal);
  }
  printf("# Rank\tE_label\t\tE_marginal\t\tV_label\t\tV_marginal\n");
  for(i=0;i<NUM_VERTICES;i++)
    printf("%d\tE %d\t%3.15Le V\t%d\t%3.15Le\n", (i+1), edge_ranking[i], weighted_edge_marginal[i], vertex_ranking[i], weighted_vertex_marginal[i]);
  for(i=NUM_VERTICES;i<NUM_EDGES;i++)
    printf("%d\tE %d\t%3.15Le\n", (i+1), edge_ranking[i], weighted_edge_marginal[i]);
  free(vertex_ranking);
  free(edge_ranking);
}

/***********************************************************************/
void init_stack(void){
  long i;
  stack_depth = 0;
  for(i=0; i<MAX_STACK_DEPTH; i++){
    my_stack[i] = EMPTY_VERTEX;
  }
  //printf("DBG INS stack_depth is %d\n",stack_depth);
}

/***********************************************************************/
void test_add_erase_b(void){
  /* this only works if #define MY_TEST_RANDOM */
  int num_test=0, my_op, np=0,nm=0;
  MY_VERTEX v1,v2;
  srandom(878713);
#ifndef MY_TEST_RANDOM
  printf("Abort. Cannot run test_add_erase_b if MY_TEST_RANDOM not set\n");
  exit(-9);
#endif
  for(;;){
    double fr;
    fr = (double)random()/(double)RAND_MAX;
    v1 = (MY_VERTEX)(fr*(double)num_vertices_now);
    if(v1==num_vertices_now)continue;
    fr = (double)random()/(double)RAND_MAX;
    v2 = (MY_VERTEX)(fr*(double)num_vertices_now);
    if(v2==num_vertices_now)continue;
    if(v1==v2)continue;
    fr = (double)random()/(double)RAND_MAX;
    my_op = (MY_VERTEX)(fr*2.)-1;
    if(num_tot_b>=B_SPACE-1)my_op = -1;
    //    printf("%g %d %d %d\n",fr, v1,v2,my_op);
    if(my_op==0){add_to_b_list(v1,v2);np++;}
    if(my_op==-1){erase_from_b_list(v1,v2);nm++;}
    num_test++;
    
    if(num_test>=10000000){
      print_b_list();
      printf("np %d nm %d nt %d\n",np,nm,num_test);
      break;
    }
  }
}

/***********************************************************************/
void erase_from_b_list(MY_VERTEX v1, MY_VERTEX v2){

  MY_VERTEX tvertex_to_erase, tv, tfirst, tlast;
  MY_VERTEX tv_x, tv_x1, tv_x2, tv_q;

  tfirst = tvertex_b_first_of_vertex[v1];
  tlast = tvertex_b_last_of_vertex[v1];

  /********************************************/
#ifdef MY_DEBUG
  if(v2==v1){
    printf("abort v2 must be different from v1 in add_to_b_list\n");
    exit(-9);
  }
  if((v1>=NUM_VERTICES) || (v2>=NUM_VERTICES)){
    printf("v1 = %d OR v2 = %d too large in erase_from_b_list\n",v1,v2);
    exit(-9);
  }
  
  if((tlast==EMPTY_TVERTEX) || (tfirst==EMPTY_TVERTEX) || (num_b[v1]==0)){
#ifdef MY_TEST_RANDOM
    return; /* ============================================= */
#endif
    printf("abort empty vertex in erase_from_b_list v1 %d v2 %d f %d l %d\n",
	   v1,v2,tlast,tfirst);
    exit(-9);
  }
#endif
  /********************************************/

  tvertex_to_erase = EMPTY_TVERTEX;
  for(tv=tfirst;;tv=tvertex_b_next_of_tvertex[tv]){
    if(vertices_of_b_list[tv]==v2){tvertex_to_erase=tv; break;}
    if(tv==tlast){break;}
  }
  if(tvertex_to_erase == EMPTY_TVERTEX){
#ifdef MY_TEST_RANDOM
    return; /* ============================================= */
#endif
    printf("abort erase_from_b_list v %d to be removed not in list of %d\n",
	   v2,v1);
    exit(-9);
  }

  /* now four different cases */
  if( (tvertex_b_former_of_tvertex[tvertex_to_erase] != MY_FIRST) 
     && (tvertex_b_next_of_tvertex[tvertex_to_erase] != MY_LAST) ){

    /* x1 o x2 , E: q L  ==> x1 x2 , E: q o L*/
    tv_x1 = tvertex_b_former_of_tvertex[tvertex_to_erase];
    tv_x2 = tvertex_b_next_of_tvertex[tvertex_to_erase];
    tv_q = tvertex_b_last_of_empty;

    tvertex_b_next_of_tvertex[tv_x1] = tv_x2;
    tvertex_b_next_of_tvertex[tv_q] = tvertex_to_erase;
    tvertex_b_next_of_tvertex[tvertex_to_erase] = MY_LAST;
    tvertex_b_former_of_tvertex[tv_x2] = tv_x1;
    tvertex_b_former_of_tvertex[tvertex_to_erase] = tv_q;
    /* here in v1 we are not first, not last */
    tvertex_b_last_of_empty = tvertex_to_erase;

  }else if( (tvertex_b_next_of_tvertex[tvertex_to_erase]== MY_LAST) 
&&  (tvertex_b_former_of_tvertex[tvertex_to_erase] != MY_FIRST) ){
    /* x o L , E: q L  ==> x L , E: q o L*/
    tv_x = tvertex_b_former_of_tvertex[tvertex_to_erase];
    tv_q = tvertex_b_last_of_empty;

    tvertex_b_next_of_tvertex[tv_x] = MY_LAST;
    tvertex_b_next_of_tvertex[tv_q] = tvertex_to_erase;
    tvertex_b_next_of_tvertex[tvertex_to_erase] = MY_LAST;
    tvertex_b_former_of_tvertex[tvertex_to_erase] = tv_q;
    /* here in v1 we are not last but not first*/
    tvertex_b_last_of_empty = tvertex_to_erase;
    tvertex_b_last_of_vertex[v1] = tv_x;
  }else if( (tvertex_b_next_of_tvertex[tvertex_to_erase]!= MY_LAST) 
&&  (tvertex_b_former_of_tvertex[tvertex_to_erase] == MY_FIRST) ){
    /* F o x , E: q L  ==> F x , E: q o L*/
    tv_x = tvertex_b_next_of_tvertex[tvertex_to_erase];
    tv_q = tvertex_b_last_of_empty;

    tvertex_b_next_of_tvertex[tv_q] = tvertex_to_erase;
    tvertex_b_next_of_tvertex[tvertex_to_erase] = MY_LAST;
    tvertex_b_former_of_tvertex[tv_x] = MY_FIRST;
    tvertex_b_former_of_tvertex[tvertex_to_erase] = tv_q;
    /* here in v1 we are not first but not last*/
    tvertex_b_last_of_empty = tvertex_to_erase;
    tvertex_b_first_of_vertex[v1] = tv_x;
  }else if( (tvertex_b_next_of_tvertex[tvertex_to_erase]== MY_LAST) 
&&  (tvertex_b_former_of_tvertex[tvertex_to_erase] == MY_FIRST) ){
    /* F o L , E: q L  ==> NOTHING , E: q o L*/
    tv_q = tvertex_b_last_of_empty;

    tvertex_b_next_of_tvertex[tv_q] = tvertex_to_erase;
    tvertex_b_former_of_tvertex[tvertex_to_erase] = tv_q;
    /* here in v1 we are both first and last */
    tvertex_b_last_of_empty = tvertex_to_erase;
    tvertex_b_first_of_vertex[v1] = EMPTY_TVERTEX;
    tvertex_b_last_of_vertex[v1] = EMPTY_TVERTEX;
  }else{
    printf("abort erase_from_b_list unkown tvertex case\n");
    exit(-9);
  }

  num_tot_b--;
  num_b[v1]--;
  num_b_empty++;
  vertices_of_b_list[tvertex_to_erase] = EMPTY_VERTEX;
}

/***********************************************************************/
void add_to_b_list(MY_VERTEX v1, MY_VERTEX v2){
  MY_VERTEX tv, tfirst, tlast, tstore;

  tfirst = tvertex_b_first_of_vertex[v1];
  tlast = tvertex_b_last_of_vertex[v1];

  /********************************************/
#ifdef MY_DEBUG
  if(v2==v1){
      printf("abort v2 must be different from v1 in add_to_b_list\n");
      exit(-9);
  }
  if((v1>=NUM_VERTICES) || (v2>=NUM_VERTICES)){
      printf("v1 OR v2 too large in add_to_b_list\n");
      exit(-9);
  }

  if(tlast!=EMPTY_TVERTEX){
    if(tfirst==EMPTY_TVERTEX){
      printf("abort internal error 1A in add_to_b_list\n");
      exit(-9);
    }
  }

  if(tfirst!=EMPTY_TVERTEX){
    if(tlast==EMPTY_TVERTEX){
      printf("abort internal error 1B in add_to_b_list\n");
      exit(-9);
    }
  }
  if((num_b[v1]==0) && (tlast!=EMPTY_TVERTEX)){
    printf("abort internal error 2 in add_b_to_list v1 %d num_b %d tlask %d\n",
	   v1,num_b[v1],tlast);
    exit(-9);
  }
#endif
  /********************************************/

  if(tfirst!=EMPTY_TVERTEX){/*will check only if non-empty*/
    for(tv=tfirst;;tv=tvertex_b_next_of_tvertex[tv]){
      if(vertices_of_b_list[tv]==v2){
	return;
      }
      if(tv==tlast)break;
    }
  }
  if(num_tot_b>=B_SPACE-1){
    printf("not enough B_SPACE 1 __ %d %d %d\n",B_SPACE,num_tot_b,num_b_empty);
    exit(-9);
  }
  if(num_b_empty<2){
    printf("not enough B_SPACE 2 __ %d %d %d\n",B_SPACE,num_tot_b,num_b_empty);
    exit(-9);
  }

  /* all that checked and done... */

  tstore = tvertex_b_last_of_empty;
  if(tlast!=EMPTY_TVERTEX){
    tvertex_b_next_of_tvertex[tlast] = tvertex_b_last_of_empty;
  }
  tvertex_b_next_of_tvertex
    [tvertex_b_former_of_tvertex
     [tvertex_b_last_of_empty]] = MY_LAST;
  tvertex_b_last_of_empty =
    tvertex_b_former_of_tvertex[tvertex_b_last_of_empty];
  if(num_b[v1]>0){
    tvertex_b_former_of_tvertex[tstore] = 
      tvertex_b_last_of_vertex[v1]; 
  }
  if(num_b[v1]==0){
    tvertex_b_former_of_tvertex[tstore] = 
		     MY_FIRST; 
    tvertex_b_first_of_vertex[v1] = tstore;
  }

  tvertex_b_last_of_vertex[v1] = tstore; 

  vertices_of_b_list[tstore] = v2;

  num_tot_b++;
  num_b_empty--;
  num_b[v1]++;
}

/***********************************************************************/
void print_b_list(void){
  long i, norm, norm2, norm3;
  MY_VERTEX tv;
  
  printf("\n# *** HERE THE B-LIST ***\n");
  norm  = 0;
  norm2 = 0;
  norm3 = 0;
  for(i=0; i<NUM_VERTICES; i++){
    norm+=num_b[i];
    if(num_b[i]>0){
      if(tvertex_b_former_of_tvertex[tvertex_b_first_of_vertex[i]]!=MY_FIRST){
	printf("abort print_b former(first)!=FIRST v %ld tfirst %d former %d\n",
	       i, tvertex_b_first_of_vertex[i],
	       tvertex_b_former_of_tvertex[tvertex_b_first_of_vertex[i]]);
	exit(-9);
      }
      if(tvertex_b_next_of_tvertex[tvertex_b_last_of_vertex[i]]!=MY_LAST){
	printf("abort in print_b former(last)!=LAST v %ld \n",i);
	exit(-9);
      }
      for(tv=tvertex_b_first_of_vertex[i];;tv=tvertex_b_next_of_tvertex[tv]){
	norm2++;
	printf("# non empty FROM VERTEX %ld BOX %d TO VERTEX %d first %d last %d\n",
	       i,tv,vertices_of_b_list[tv]
	       ,tvertex_b_first_of_vertex[i], tvertex_b_last_of_vertex[i]);
	if(vertices_of_b_list[tv]==EMPTY_VERTEX){
	  printf("abort internal error in print_b_list EMPTY_VERTEX\n");
	  exit(-9);
	}
	if(tv==tvertex_b_last_of_vertex[i])break;
      }
      for(tv=tvertex_b_last_of_vertex[i];;tv=tvertex_b_former_of_tvertex[tv]){
	norm3++;
	if(tv==tvertex_b_first_of_vertex[i])break;
      }
    }
  }  
  norm+= num_b_empty;
  if(num_b_empty>0){
    for(tv=tvertex_b_first_of_empty;;tv=tvertex_b_next_of_tvertex[tv]){
      norm2++;
      printf("#     empty BOX %d TO VERTEX %d first %d last %d\n",tv,vertices_of_b_list[tv],tvertex_b_first_of_empty,tvertex_b_last_of_empty);
      if(tv==tvertex_b_last_of_empty)break;
    }
    for(tv=tvertex_b_last_of_empty;;tv=tvertex_b_former_of_tvertex[tv]){
      norm3++;
      if(tv==tvertex_b_first_of_empty)break;
    }
  }
  if(norm!=B_SPACE){
    printf("abort norm != B_SPACE %ld %d\n",norm,B_SPACE);
    exit(-9);
  }
  if(norm2!=B_SPACE){
    printf("abort norm2 != B_SPACE %ld %d\n",norm,B_SPACE);
    exit(-9);
  }
  if(norm3!=B_SPACE){
    printf("abort norm3 != B_SPACE %ld %d\n",norm,B_SPACE);
    exit(-9);
  }
}

/***********************************************************************/
void init_b_list_to_empty(void){
  long i;

  tvertex_b_first_of_empty = 0;
  tvertex_b_last_of_empty = B_SPACE-1;
  num_b_empty = B_SPACE;
  for(i=0; i<NUM_VERTICES; i++){
    tvertex_b_first_of_vertex[i] = EMPTY_TVERTEX;
    tvertex_b_last_of_vertex[i]  = EMPTY_TVERTEX;
    num_b[i] = 0;
  }
  for(i=0; i<B_SPACE; i++){
    tvertex_b_next_of_tvertex[i]=i+1;
    tvertex_b_former_of_tvertex[i]=i-1;
    vertices_of_b_list[i] = EMPTY_VERTEX;
  }
  tvertex_b_next_of_tvertex[B_SPACE-1] = MY_LAST;
  tvertex_b_former_of_tvertex[0] = MY_FIRST;
  num_tot_b=0;
}

/***********************************************************************/
void chop_vertices(MY_VERTEX v_min){
  long i,j; 

  num_edges_now = 0;
  num_vertices_now = NUM_VERTICES - v_min;
  max_vertex = 0;
  j=0;
  for(i=0;i<NUM_EDGES;i++){
    if((vertices_of_edge_list0[i][0]>=v_min) &&
       (vertices_of_edge_list0[i][1]>=v_min)){ 
      num_edges_now++;
      vertices_of_edge_list[j][0] =     vertices_of_edge_list0[i][0];
      vertices_of_edge_list[j][1] =     vertices_of_edge_list0[i][1];
      if(vertices_of_edge_list[j][0]>=max_vertex){
	max_vertex = vertices_of_edge_list[j][0];
      }
      j++;
    }
  }
#ifdef MY_DEBUG
  printf(" ** DBG CHOPCHOP %d ** max_vertex %d \n",
	 v_min, max_vertex);
#endif
  for(i=num_edges_now;i<NUM_EDGES;i++){
    vertices_of_edge_list[i][0] = EMPTY_VERTEX;
    vertices_of_edge_list[i][1] = EMPTY_VERTEX;
  }
#if ( (defined MY_DEBUG) ||  (defined MY_DEBUG2 ) )
  for(i=0;i<num_edges_now;i++){
    printf(" ** DBG CHOPCHOP %d ** i %ld start %d end %d \n",
	   v_min, i,vertices_of_edge_list[i][0],vertices_of_edge_list[i][1]);
  }
  printf(" ** DBG CHOPCHOP %d ** num_edges_now %d \n",
	 v_min, num_edges_now);
  fflush(stdout);
#endif
}

/***********************************************************************/
void pruneItUp(void){  

  long i, j, level=0, k;

  lengthOne = 0;

  for(;;){
    int flag0, flag1, mySecond;

    flag0=1;
    flag1=1;
    level++;
    for(i=0; i<NUM_VERTICES; i++){
      if(endVertexConnection[i]==NO_VERTEX-1){
	flag0 = 0;
	for(j=0; j<=maxOccupied; j++){
	  if(vertices_of_edge_list0[j][0]>i){
	    vertices_of_edge_list0[j][0]--;
	  } 
	  if(vertices_of_edge_list0[j][1]>i){
	    vertices_of_edge_list0[j][1]--;
	  } 
	}
	NUM_VERTICES--;
	break;
      }
    }
    if(flag0 == 0){
      for(j=i; j<NUM_VERTICES; j++){
	numberVertexConnection[j] = numberVertexConnection[j+1];
	endVertexConnection[j] = endVertexConnection[j+1];
	startVertexConnection[j] = startVertexConnection[j+1];
      }
    }
     
    for(i=0; i<NUM_VERTICES; i++){
      if(numberVertexConnection[i]==1){
	flag1 = 0;
	numberVertexConnection[i] = 0;
	startVertexConnection[i] = NO_VERTEX;
	endVertexConnection[i] = NO_VERTEX-1;
	lengthOne++;
	break;
      }
    }	

    if(flag1==0){
      for(j=0; j<=maxOccupied; j++){
	if(vertices_of_edge_list0[j][0]==i){
	  mySecond = vertices_of_edge_list0[j][1];
	  numberVertexConnection[mySecond]--;
	  if(numberVertexConnection[mySecond]==0){
	    startVertexConnection[mySecond] = NO_VERTEX;	    
	    endVertexConnection[mySecond] = NO_VERTEX-1;	    
	  }
	  break;
	} 
      }
      for(k=j; k<maxOccupied; k++){
	vertices_of_edge_list0[k][0] = vertices_of_edge_list0[k+1][0];
	vertices_of_edge_list0[k][1] = vertices_of_edge_list0[k+1][1];
      }
      for(k=0; k<NUM_VERTICES; k++){
	if(startVertexConnection[k]>j){
	  startVertexConnection[k] --;
	}
	if(endVertexConnection[k]>j){
	  endVertexConnection[k] --;
	}
      }
      
      maxOccupied--;
      
      for(j=0; j<=maxOccupied; j++){
	if((vertices_of_edge_list0[j][0]==mySecond)
	   && (vertices_of_edge_list0[j][1]==i) ){
	  break;
	} 
      }
      
      for(k=j; k<maxOccupied; k++){
	vertices_of_edge_list0[k][0] = vertices_of_edge_list0[k+1][0];
	vertices_of_edge_list0[k][1] = vertices_of_edge_list0[k+1][1];
      }
      for(k=0; k<NUM_VERTICES; k++){
	if(startVertexConnection[k]>j){
	  startVertexConnection[k] --;
	  endVertexConnection[k] --;
	}
      }
      maxOccupied--;
      
      NUM_EDGES--;
    }
    if((flag0==1) && (flag1==1))break;
  }

#ifdef WRITE_GRAPH
#ifdef PRUNING
  for (i=0; i<NUM_EDGES;i++){
    printf("PRUNED  I %d V1 %d V2 %d \n",
	   i,
	   vertices_of_edge_list0[i][0],
	   vertices_of_edge_list0[i][1]);
  }
#endif
#endif
}

/***********************************************************************/
void read_vertices_from_edges_dat(void){  

  FILE* ident;
  int i, j;
  MY_VERTEX a,b, temp_first, temp_second;
  double var_weight;

  ident=fopen(filename,"r");
  fscanf(ident,"%ld %ld\n",&NUM_VERTICES,&NUM_EDGES); 
  if(NUM_EDGES >= NUM_EDGES_MAX){
    printf("select vertices abort too large NUM_EDGES %ld NUM_EDGES_MAX %d\n",
	   NUM_EDGES,NUM_EDGES_MAX);
    exit(-9);
  }
  if(NUM_VERTICES >= NUM_VERTICES_MAX){
    printf("select vertices abort too large NUM_VERTICES %ld NUM_VERTICES_MAX %d\n",
	   NUM_VERTICES,NUM_VERTICES_MAX);
    exit(-9);
  }
  for(i=0;i<NUM_VERTICES;i++)
    for(j=0;j<NUM_VERTICES;j++)
      adjacency[i][j]=-1;
  for(i=0; i<NUM_EDGES; i++)
    edge_weight[i]=0.;
  for(i=0; i<NUM_EDGES; i++)
    {
      if(fscanf(ident,"%i\t%i\t%lf",&a, &b, &var_weight)==EOF){
	printf("abort f_ver a,b\n");
	exit(-9);
      }
      vertices_of_edge_list0[i][0] = a;  
      vertices_of_edge_list0[i][1] = b; 
      vertices_of_original_edges[i][0] = a;  
      vertices_of_original_edges[i][1] = b;   
      adjacency[a][b]=i;
      edge_weight[i]= (long double)var_weight;
    }
  fclose(ident);

  /*Here we need to do some bubble sort such that the vertices of verticeS_of_edge_list0 are sorted first on the first([0]) vertex and secondly on the second([1]) vertex. */
  for(i=0; i<(NUM_EDGES-1); i++)
    {
      for (j=(NUM_EDGES-1); j>i; j--)
	{
	  if (vertices_of_edge_list0[j-1][0]>vertices_of_edge_list0[j][0])
	    {
	      temp_first=vertices_of_edge_list0[j-1][0];
	      vertices_of_edge_list0[j-1][0]=vertices_of_edge_list0[j][0];
	      vertices_of_edge_list0[j][0]=temp_first;
	      temp_second=vertices_of_edge_list0[j-1][1];
	      vertices_of_edge_list0[j-1][1]=vertices_of_edge_list0[j][1];
	      vertices_of_edge_list0[j][1]=temp_second;
	    }
	}
    }
  for(i=0; i<(NUM_EDGES-1); i++)
    {
      for (j=(NUM_EDGES-1); j>i; j--)
	{
	  if ((vertices_of_edge_list0[j-1][0]==vertices_of_edge_list0[j][0]) && (vertices_of_edge_list0[j-1][1]>vertices_of_edge_list0[j][1]))
	    {
	      temp_second=vertices_of_edge_list0[j-1][1];
	      vertices_of_edge_list0[j-1][1]=vertices_of_edge_list0[j][1];
	      vertices_of_edge_list0[j][1]=temp_second;
	    }
	}
    }
}

/***********************************************************************/
int check_and_set_vertices_of_edge_list(MY_VERTEX v_min){
  long i, j;

  chop_vertices(v_min);

  for(i=0; i<num_edges_now; i++){
    if(vertices_of_edge_list[i][0]==vertices_of_edge_list[i][1]){
      printf("ABORT check_and_set_vertices_of_edge_list"
	     " not self edges allowed %ld %d %d\n",
	     i,vertices_of_edge_list[i][0],vertices_of_edge_list[i][1]);
      exit(-9);
    }
  }
  for(i=0; i<num_edges_now; i++){
    for(j=0; j<num_edges_now; j++){
      if(i!=j){
	if( (vertices_of_edge_list[i][0]==vertices_of_edge_list[j][0])
	    && (vertices_of_edge_list[i][1]==vertices_of_edge_list[j][1]) ){
	  printf("ABORT check_and_set_vertices_of_edge_list"
		 " not double edges allowed %ld %ld %d %d\n",
		 i,j,vertices_of_edge_list[i][0],vertices_of_edge_list[i][1]);
	  exit(-9);
	}
      }
    }
  }

  for(i=1; i<num_edges_now; i++){
    if(vertices_of_edge_list[i][0]<vertices_of_edge_list[i-1][0]){
    printf("ABORT check_and_set_vertices_of_edge_list"
	   " is not ordered in first vertex %ld %d %d\n",
	   i, vertices_of_edge_list[i-1][0], vertices_of_edge_list[i][0]);
    exit(-9);
    }
  }

  /* we set tedge_first and tedge_last */
  for(i=0; i<NUM_VERTICES; i++){
    tedge_first_of_vertex[i] = EMPTY_TVERTEX;
    tedge_last_of_vertex[i]  = EMPTY_TVERTEX;
  }

  if(num_edges_now>0){
    tedge_first_of_vertex[vertices_of_edge_list[0][0]] = 0;
    tedge_last_of_vertex[max_vertex] = num_edges_now-1;
  }
  for(i=1; i<num_edges_now; i++){
    if(vertices_of_edge_list[i][0]>vertices_of_edge_list[i-1][0]){
      if((vertices_of_edge_list[i-1][0]>=NUM_EDGES) || 
	 (vertices_of_edge_list[i-1][0]>=NUM_EDGES)){
      }
      tedge_last_of_vertex[vertices_of_edge_list[i-1][0]] = i-1;
      tedge_first_of_vertex[vertices_of_edge_list[i][0]] = i;
    }
  }
#if ( (defined MY_DEBUG) ||  (defined MY_DEBUG2 ) )
  for(i=0; i<NUM_VERTICES; i++){
    if ((tedge_first_of_vertex[i]!=EMPTY_TVERTEX) && (tedge_last_of_vertex[i]!=EMPTY_TVERTEX)){
      printf("# I %ld VF %d VL %d E_F %d E_L %d\n",i,
	     vertices_of_edge_list[tedge_first_of_vertex[i]][0],
	     vertices_of_edge_list[tedge_last_of_vertex[i]][0],
	     tedge_first_of_vertex[i],tedge_last_of_vertex[i]);
    }
    else {
      printf("# I %ld E_F %d E_L %d\n",i,
	     tedge_first_of_vertex[i],tedge_last_of_vertex[i]);
    } 
  }
  printf("#edges for %d =%d\n", v_min, num_edges_now);
  fflush(stdout);
#endif
  /* now we check that list in second vertex is ordered */
  for(i=0; i<num_vertices_now; i++){
    long j;
    if(tedge_last_of_vertex[i]>tedge_first_of_vertex[i]){
      for(j=tedge_first_of_vertex[i]+1; j<=tedge_last_of_vertex[i]; j++){
	if(vertices_of_edge_list[j][0]==vertices_of_edge_list[j-1][0]){
	  if(vertices_of_edge_list[j][1]<=vertices_of_edge_list[j-1][1]){
	    printf("ABORT check_and_set_vertices_of_edge_list"
		   " problem in [1] i %ld j %ld tf %d tl %d vj-1 %d vj %d\n",
		   i, j, tedge_first_of_vertex[i], 
		   tedge_last_of_vertex[i],
		   vertices_of_edge_list[j-1][1], 
		   vertices_of_edge_list[j][1]);
	    printf("E_F %d E_L %d\n",
		   tedge_first_of_vertex[i],tedge_last_of_vertex[i]);
	    exit(-9);
	  }
	}
      }
    }
  }
  return num_edges_now;
}