/*  An efficient software for solving the median order problem

   Irene Charon and Olivier Hudry
    ----------------------------------------------------------

    This software solves problems modelized in the following form; 
a weighted tournament T with N vertices is given; the set of vertices 
of T is denoted by X;  the weights are  non-negative integers; the 
value of a total order O is the sum of the weights of the arcs (x, y) 
of T such that  y > x in O.
The problem is to find an order of minimum value.
 
   The vertices of T are denoted by 0, 1, ..., N-1.
T is coded by an asymmetric matrix which name is 
"tournament" ; if (i,j) is an arc of T, tournament[i][j] is
equal to the weight w of this arc and tournament[j][i]=-w. 
For every i,tournament[i][i]=0. We say that "tournament" contains
the signed values of the weights.  
   The data are read on a file. This file must contain :
first, the number N of vertices
after, for 0<=i<N , for 0<=j<=i, the value of tournament[i][j]
For example, if T has 4 verices 0,1,2,3 and the arcs :
(1,0) of weight 3
(0,2) of weight 2
(3,0) of weight 4
(1,2) of weight 1
(3,1) of weight 0
(3,2) of weight 7
the input file is :
4
 0
 3  0
-2 -1 0
 4  0 7 0
When the program begins, the name of the input file is asked.
After the tournament is get, the transitivity index is computed
and given.

    After the noising method is performed, the best value and the
best order found by this method are given.
    After the first relax applied at the root of the branch and 
bound tree, the maximum value computed for the dual function
is indicated; it gives a lower bound of the best value for
the median order problem.
    If the branch and bound method takes too much time, the user
of the software already has some informations on the solution of 
his problem. 

   At the end, the output indicates the best value and the 
best order.
*/
 
#include <string.h>
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <sys/time.h>
#include <sys/resource.h>

#define eps  0.00001 

void initialization(void);
int OrderValue(int *,int);
void CopyOrder(int *,int *,int);
char IsValid(int,int);
void AddToBenefits(int,int);
void DeleteFromBenefits(int,int);
void change(int,int,int *, int * ); 
void InitBenefits(void);
double ComputeIndex(void);
void PutOrder(void);       
void noising(double);
int descent(int,int *);
void BB(int,int);
void InitBB(double);
int relax(int *, int, int, int);
int  ComputeSigma(int * , int);
void InitRelax(double);
char exist(int,int);
void InitHead(void);
void PutLowerBound(int);
int BegSecLowerBound(void);

/********************* Global variables *******************/

int  N;   /*the number of vertices of the tournament*/
int ** tournament;    /* gives the address of a matrix which contains
			 the signed  values of the weights*/ 
int *  CurrentOrder;  /* gives the address of an array which
			 contains orders which are useful
			 during the noising and the B andB*/
int  * BestOrder;     /* gives the address of an array which
			 always contains the best encountered order*/
int bound=INT_MAX;    /* contains the value of BestOrder*/
int MaxWeight=0;      /* contains the maximum weight of the arcs
			 of the tournament*/
char DoSigma=1;       /* is equal to 1 if the evaluation with sigma 
			 is performed, 0 otherwise*/
char DoRelax=1;       /* is equal to 1 if the evaluation with
			 relaxation is performed, 0 otherwise*/
char DoExists=1;      /* is equal to 1 if the lexicographic tree
			 containing encountered beginning sections
			 is built and used, 0 otherwise*/
int * score;          /* gives the address of an array which
			 contains scores of vertices in the function 
			 computing sigma and in relaxation*/
char * IsOutBegSec;   /* gives the address of an array which
			 indicates, for each vertex, if it
			 is (value 1) or not (value 0) in the 
			 current beginning section of the B 
			 and B; this array is useful for 
			 relaxations and for "exist"*/
int * place;          /* gives the address of an array which
			 is used in noising and descents*/
int drift;

/******************** End of global variables  *************/	

/***************** Functions for the descent method **********/

/*************** prototypes of local functions **************/

void GetTournament(FILE *); /* allocates place for "tournament"
			       and get the tournament signed 
			       weights  from the input file*/
void PutResults();          /* indicates, at the end, the optimal 
			       value and the optimal order*/
/*********** end of prototypes of local functions **********/

int BestPlace(int OldPlace, int *NewPlaceAddr,int LocalN, int * order)
{
  int i;
  int minimum = 0;
  int delta=0;
  int vertex=order[OldPlace];

  *NewPlaceAddr=OldPlace;
  for (i= OldPlace-1; i>=0;i--)
    {
      delta+=tournament[order[i]][vertex];
      if (delta < minimum)
	{
	  minimum=delta;
	  *NewPlaceAddr=i;
	}
    }
  delta=0;
  for (i= OldPlace+1; i< LocalN; i++)
    {
      delta+=tournament[vertex][order[i]];
      if (delta < minimum)
	{
	  minimum=delta;
	  *NewPlaceAddr=i;
	}
    }
  return minimum;
}



int descent (int LocalN, int * order)
{
  int index1,index2;
  int delta;
  int i;
  int NbUnchanged=0;
  int vertex;
  int index=0;

  for (i=0;i<LocalN;i++)
    {
      place[order[i]]=i;
    }
  index1=0;
  vertex==order[0];
  do
    {
      delta=BestPlace(index1,&index2,LocalN,order);
      if (delta<0)
	{
	  change(index1,index2,order,place); 
	  NbUnchanged=0;
	}
      else NbUnchanged++;
      index = (index + 1) % LocalN;
      vertex=order[index];
      index1=place[vertex];
  } while(NbUnchanged<LocalN);
  return(OrderValue(order,LocalN));
}  



typedef struct
{
  unsigned char vertex1,vertex2,vertex3;
  double value;
}Circuit;


double DualFunction(int);
char  UpdateDiffLambda(int, double,int,char);
void SaveOrder(int *, int);
void Decrease_ro(double);
double  ComputeGain(int);
void SortCircuits(void);
void InitOneRelax(int *, int);

int RelaxNum=0;          /* increased by one at each relaxation */
int  * SortedVertices;   /* to perform relaxation, it is necessary to 
			    have the list of the vertices of the studied
			    sub-tournament, in lexicographic order;
			    they are at the address SortedVertices */
double *** DiffLambda;   /* address of the 3-dimensionnal table which
			    contains the lambda[i][j][k]-mu[i][j][k] */
char ** X;               /* address of a matrix which gives the  
			    minimizing L(lambda,mu) for given tournament
			    lambda and mu. X[i][j] = 1 if the tournament 
			    contains  the arc 
			    (SortedVertices[i],SortedVertices[j])
			    and 0 otherwise */
			    
int PositiveTriangle;  /* constant used to compute the dual function */
int NbIterMax1, NbIterMax2, NbIterMax;  /* number of iterations respectively 
			     for the first relaxation, for the
			     other relaxation, for the current relaxation */ 
double ro1, ro2, ro;/* decreasing coefficient for the size of a step
		       of relaxation respectively for the first relaxation,
		       for the  other relaxation,s for the current relaxation*/
			    
int IterNum; /* increased by one at each step of relaxation */
double ** SaveCoeff; /* an array used to compute gain; it contains
		      the value of coeff[i][j] described in the paper*/
int NbCircuits; /* number of circuits, used in ComputeGain */
Circuit * circuits; /* list of circuits */
int MaxNbCircuits; 
char DoGain=1;  /* equal to 1 if gain is computed, 0 otherwise */
double OldL=0; /* keeps the previous value of L */


int relax(int * LocalVertices, int LocalN, int LocalBound, 
	  int LocalBestOrderValue)
{
  /* LocalVertices contains the list of the vertices on which 
     we apply relaxation; LocalBestOrderValue is the best known value 
     for an order on the vertices contained by LocalVertices; 
     ifIntegerL, the smallest integer bigger than a computed 
     value of the dual function, is greater or equal to 
     LocalBound,  the function returns 0; this  returned value indicates
     that there is no optimal order beginning by the current beginning 
     section; if, after NbIterMax iterations, IntegerrL is never greater
     or equal to LocalBound, the function returns 1*/
  double L;       /* for a value of the dual function */
  int IntegerL;
  char end=0;
  int i,j,k;
  double gain;
  int BestL=0;
  int CurrentOrderValue;
  char PrepareGain;  /* 1 if, in function UpdateDiffLambda, some work
		       must be done to prepare computing of gain,
		       0 otherwise */
  double OldML = LocalBound, ML;
  double rohigh = ro1;
  double min = 0.05;

  IterNum=0;
  NbCircuits=0;
  if (DoGain) PrepareGain=1;
  else PrepareGain=0;
  if (RelaxNum == 0)
    {
      NbIterMax = NbIterMax1;
      ro = ro1 ;
      ro = ro2;
    }
  else 
    {
      NbIterMax = NbIterMax2;
      ro = ro2;
    }
  InitOneRelax(LocalVertices,LocalN);
  if (DoGain) PrepareGain=1;
  else PrepareGain=0;
  IntegerL = 0;
  
  ML = OldL = L = DualFunction(LocalN);
  if (L < 0) IntegerL = -1;
  else 
    if (L - (int)L > eps) IntegerL = (int)L+1;
    else IntegerL = (int)L;
  if (IntegerL > BestL) BestL = IntegerL;
  if (BestL >= LocalBound) end = 2;
  while (!end)
    {
      IterNum++;
      if (UpdateDiffLambda(LocalN,L,LocalBestOrderValue,PrepareGain)) {
	/* UpdateDiffLambda returns 1 if the tournament given by X
	   is transitive, 0 otherwise;
	   so, here, X gives an order; its value is computed below */
	CurrentOrderValue=PositiveTriangle;
	for (i=0;i<LocalN-1;i++)
	  for (j=i+1;j<LocalN;j++) 
	    if (X[i][j]) CurrentOrderValue+=
			   tournament[SortedVertices[j]][SortedVertices[i]];
	if (CurrentOrderValue < LocalBound) 
	  {
	    bound+=CurrentOrderValue-LocalBound;
	    LocalBound=CurrentOrderValue;
	    LocalBestOrderValue=CurrentOrderValue;
	    SaveOrder(LocalVertices,LocalN);
	    printf("An order of value %d is found.\n it is :"
		   , bound);
	    PutOrder();
	    if (IntegerL>=LocalBound)  end=2;
	  }
      }
      else if (PrepareGain)
	{
	  gain=ComputeGain(LocalN);
	  L+=gain;
	  if (L<0) IntegerL=-1;
	  else 
	    if (L-(int)L>eps) IntegerL=(int)L+1;
	    else IntegerL=(int)L;
	  if (IntegerL>BestL) BestL=IntegerL;
	  if (BestL >= LocalBound)  
	    {
	      end=2;
	      break;
	    }
	  PrepareGain=0; /*for each relaxation, after a step,
			   we stop the calculus of the gain, because
			   it takes too much time for too few chance 
			   of cutting a branch */
	  gain=0;
	}
      L = DualFunction(LocalN); 
      if (L > ML) ML = L;
      if (L < 0) IntegerL = -1;
      else 
	if (L - (int)L > eps) IntegerL = (int)L+1;
	else IntegerL = (int)L;
      if (IntegerL > BestL) BestL = IntegerL;
      if (BestL >= LocalBound) end = 2;
      if (RelaxNum == 0) {
	if (L < 0) {
	  ro *= 0.5;
	  ro2 *= 0.99;
	  for (i=0;i<N-2;i++)
	    for (j=i+1;j<N-1;j++) 
	      for (k=j+1;k<N;k++) DiffLambda[i][j][k]=0;
	  OldL = 0;
	}
	if (IterNum > NbIterMax - 500) ro *= 0.999;
	if ((L - OldL)/L > 0.1) ro *= 1.1;
	else if ((L - OldL)/L < - 0.95) ro *= 0.8;
	else ro *= 1 + (L - OldL)/L;
	if (IterNum == 400) OldML = BestL;
	else if ((IterNum >= 500 ) && (!(IterNum% 100))) {
	  if ((ML - OldML)/(LocalBound - OldML) > min) { 
	    NbIterMax += 100;
	    //printf("NbIterMax=%d\n",NbIterMax); 
	  }
	  OldML = ML;
	  if (ro > rohigh) rohigh = ro;
	}
      }
      else {
	if (L < 0) {
	  ro *= 0.5;
	  ro2 *= 0.99;
	}
	else if (L < OldL * 0.98) {
	  ro *= 0.9;
	  ro2 *= 0.99;
	}
	ro *= 0.95;
      }
      OldL = L;
      /*if ((!(IterNum%500))||((!RelaxNum)&&(IterNum<=10))) {
	printf("iter=%d, L %lf, ro %lf, p= %d,  BorneLocal = %d\n",IterNum, L, ro, N-LocalN,  LocalBound);
	fflush(stdout);
      }*/
  if (IterNum >= NbIterMax)   end = 1; 
    }
  if (RelaxNum == 0) {
    ro2 = rohigh * 0.8;
    printf("\nAt the beginning, Lagrangean relaxation "
	   "gives the lower bound %d.\n",BestL);  
    fflush(stdout);
  }
  RelaxNum++;
  if (DoExists) PutLowerBound(BestL);
  if (end==2) return 1;
  else return 0;
}
/************************ end of relax ****************************/



/************  computes a value  of the lagrangean function **************/

double DualFunction(int LocalN)
{

  int i,j,k,s,SVi,SVj;
  double coeff;
  double diff;
  double L;
  double gain;
  char * Xi;
  
  L=PositiveTriangle;
  for (i=0;i<LocalN-1;i++)
    {
      Xi=X[i];
      SVi=SortedVertices[i];
      for (j=i+1;j<LocalN;j++)
	{
	  SVj=SortedVertices[j];
	  coeff=(double)tournament[SVj][SVi];   
	  for (k=j+1;k<LocalN;k++)
	    {
	      diff=DiffLambda[SVi][SVj][SortedVertices[k]];
	      if (diff>0)  L-=diff;
	      coeff+=diff;
	    }
	  for (k=0;k<i;k++)
	    coeff+=DiffLambda[SortedVertices[k]][SVi][SVj];
	  for (k=i+1;k<j;k++)
	    coeff-=DiffLambda[SVi][SortedVertices[k]][SVj]; 
	  if (coeff<0)
	    {
	      L+=coeff;
	      Xi[j]=1;
	    }
	  else Xi[j]=0;
	  if (DoGain) SaveCoeff[i][j]=coeff;
	}
    }
  return L;
}
/****************** end of DualFunction ***************************/



/************ this function essentially computes new values of the 
  variables DiffLambda, making a step; it also returns 1 if the tournament
  corresponding to X is transitive, and  0 otherwise;
  furthermore, it keeps in memory the circuits used by the function
  ComputeGain, when it  is necessary.
  The function can sees a bit complex: ***************************/

char UpdateDiffLambda(int LocalN, double L, int LocalBestOrderValue,
		      char PrepareGain)
{
  int i,j,k,Ok;
  double SimpleStep, DoubleStep, TripleStep;
  double NegSimpleStep,NegDoubleStep,NegTripleStep;
  double coij;
  int S;
  double D;
  char transitive=1;
  Circuit UnCircuit;
  double gain;
  double ** diffi, * diffij, * coi, * coj;
  char * Xi, * Xj;
 
  SimpleStep=(ro*(LocalBestOrderValue-L))/(LocalN*(LocalN-1)*(LocalN-2));
  DoubleStep=2*SimpleStep;
  TripleStep=3*SimpleStep;
  NegSimpleStep=-SimpleStep;
  NegDoubleStep=-DoubleStep;
  NegTripleStep=-TripleStep;
  NbCircuits=0;
  if (PrepareGain)
    {
      for (i=0;i<LocalN-2;i++)
	{
	  diffi=DiffLambda[SortedVertices[i]];
	  Xi=X[i];
	  coi=SaveCoeff[i];
	  for (j=i+1;j<LocalN-1;j++)
	    {
	      Xj=X[j];
	      coj=SaveCoeff[j];
	      diffij=diffi[SortedVertices[j]];
	      coij=coi[j];
	      for (k=j+1;k<LocalN;k++)
		{
		  Ok=SortedVertices[k];
		  S= Xi[j]+ Xj[k]- Xi[k];
		  if (S==-1)  
		    {
		      transitive=0;
		      UnCircuit.vertex1=i;
		      UnCircuit.vertex2=j;
		      UnCircuit.vertex3=k;
		      if (coij<coj[k]) gain=coij;
		      else gain=coj[k];
		      if (-coi[k]<gain) gain=-coi[k];
		      UnCircuit.value=gain;
		      circuits[NbCircuits++]= UnCircuit;
		    }
		  else if (S==2)
		    {
		      transitive=0;
		      UnCircuit.vertex1=i;
		      UnCircuit.vertex2=j;
		      UnCircuit.vertex3=k;
		      if (coij<coj[k]) gain=-coj[k];
		      else gain=-coij;
		      if (coi[k]<gain) gain=coi[k];
		      UnCircuit.value=gain;
		      circuits[NbCircuits++]= UnCircuit;
		    }
		  D=diffij[Ok];
		  switch(S)
		    {
		    case 2 :
		      if (D>0)
			diffij[Ok]+=SimpleStep;
		      else if (D>NegDoubleStep) 
			diffij[Ok]=SimpleStep;
		      else  diffij[Ok]+=TripleStep;
		      break;
		    case 1 :
		      if (D<NegSimpleStep)
			diffij[Ok]+=SimpleStep;
		      else  if (D<0) 
			diffij[Ok]=0;
		      break;
		    case 0 :
		      if (D>SimpleStep)
			diffij[Ok]+=NegSimpleStep;
		      else  if (D>0) 
			diffij[Ok]=0;
		      break;
		    case -1 :
		      if (D>DoubleStep)
			diffij[Ok]+=NegTripleStep;
		      else if (D>0) 
			diffij[Ok]=NegSimpleStep;
		      else  diffij[Ok]+=NegSimpleStep;
		      break;
		    }
		}
	    }
	}
    }
  else
    {
      for (i=0;i<LocalN-2;i++)
	{
	  diffi=DiffLambda[SortedVertices[i]];
	  Xi=X[i];
	  for (j=i+1;j<LocalN-1;j++)
	    {
	      Xj=X[j];
	      diffij=diffi[SortedVertices[j]];
	      for (k=j+1;k<LocalN;k++)
		{
		  Ok=SortedVertices[k];
		  S= X[i][j]+ X[j][k]- X[i][k];
		  if ((S==-1)||(S==2))  transitive=0;
		  D=diffij[Ok];
		  switch(S)
		    {
		    case 2 :
		      if (D>0)
			diffij[Ok]+=SimpleStep;
		      else if (D>NegDoubleStep) 
			diffij[Ok]=SimpleStep;
		      else  diffij[Ok]+=TripleStep;
		      break;
		    case 1 :
		      if (D<NegSimpleStep)
			diffij[Ok]+=SimpleStep;
		      else  if (D<0) 
			diffij[Ok]=0;;
		      break;
		    case 0 :
		      if (D>SimpleStep)
			diffij[Ok]+=NegSimpleStep;
		      else  if (D>0) 
			diffij[Ok]=0;
		      break;
		    case -1 :
		      if (D>DoubleStep)
			diffij[Ok]+=NegTripleStep;
		      else if (D>0) 
			diffij[Ok]=NegSimpleStep;
		      else  diffij[Ok]+=NegSimpleStep;
		      break;
		    }
		}
	    }
	}
    }
  return transitive;
}  
/*********************** end of UpdateDiffLambda ******************/



/******** used to save the order when relaxation finds a best order
         than the current BestOrder ******************************/

void SaveOrder(int * LocalVertices,int LocalN)
{
  int i,j;
  int index;

  for (i=0;i<LocalN;i++) score[i]=0;
  for (i=0;i<LocalN-1;i++)
    for (j=i+1;j<LocalN;j++)
      if (X[i][j]) score[i]++;
      else score[j]++;
  for (i=0;i<LocalN;i++)
    {
      index = 0;
      while (score[index]!=LocalN-i-1) index++;
      LocalVertices[i]=SortedVertices[index];
    }
  CopyOrder(BestOrder,LocalVertices-N+LocalN,N);
}
/********************* end of SaveOrder **************************/



/************* this function computes the gain, 
          as explained in the paper ******************************/

double ComputeGain(int LocalN)
{
  double gain;
  double A,B,C;
  Circuit ACircuit;
  int i,j;

  SortCircuits();
  gain = 0;
  for (i=0;i<NbCircuits;i++)
    {
      ACircuit=circuits[i];
      if ((X[ACircuit.vertex1][ACircuit.vertex2]<10)&&
	  (X[ACircuit.vertex2][ACircuit.vertex3]<10)&&
	  (X[ACircuit.vertex1][ACircuit.vertex3]<10))
	{
	  gain += ACircuit.value;
	  X[ACircuit.vertex1][ACircuit.vertex2]+=10;
	  X[ACircuit.vertex2][ACircuit.vertex3]+=10;
	  X[ACircuit.vertex1][ACircuit.vertex3]+=10;
	}
    }
  for(i=0;i<LocalN-1;i++)
    for(j=i+1;j<LocalN;j++)
      if (X[i][j]>=10) X[i][j]-=10;
  return gain;
}
/********************** end of ComputeGain ************************/



/********************** this function sorts the circuits 
               used by computeGain, using Quick Sort ***************/

void SortCircuits()
{
  int i,j;
  double  KeyValue;
  int IndexSmallChild;
  char end;
  Circuit ACircuit;

  circuits-=1;
  for (i=2;i<= NbCircuits;i++)
    {
      ACircuit=circuits[i];
      KeyValue=ACircuit.value;
      j=i;
      while((j>1)&&(circuits[j/2].value>KeyValue))
	{
	  circuits[j]=circuits[j/2];
	  j=j/2;
	}
      circuits[j]=ACircuit;
    }
  for (i=NbCircuits;i>=2;i--)
    {
      end=0;
      ACircuit=circuits[i];
      circuits[i]=circuits[1];
      circuits[1]=ACircuit;
      KeyValue=ACircuit.value;
      j=1;
      while((2*j<=i-1)&&(!end))
      {
	if (2*j==i-1) IndexSmallChild=i-1;
	else
	  if (circuits[2*j].value>circuits[2*j+1].value)
	    IndexSmallChild=2*j+1;
	  else IndexSmallChild=2*j;
	if (circuits[IndexSmallChild].value<KeyValue)
	  {
	    circuits[j]=circuits[IndexSmallChild];
	    j=IndexSmallChild;
	  }
	else 
	  {
	    end=1;
	  }
      }
      circuits[j]=ACircuit;
    }
  circuits+=1;
}
/************************* end of SortCircuits *******************/



/********************** allocations and initializations 
            before the first relaxation***************************/
	  
void InitRelax(double TransitivityIndex)
{
  int i,j,k,h;

  ro2 = ro1 = 2;
  NbIterMax1 = 1000; 
  NbIterMax2 = 10;
  SortedVertices = (int *)malloc(N*sizeof(int));
  if (SortedVertices == NULL) 
    {
      DoRelax=0;
      return;
    }
  X = (char **) malloc((N-1)*sizeof(char *));
  if (X == NULL) 
    {
      DoRelax = 0;
      free(SortedVertices);
      return;
    }
  for (i=0;i<N-1;i++)
    {
      X[i]=(char *) malloc((N-1-i)*sizeof(char));
      if (X[i]==NULL)
	{
	  DoRelax=0;
	  for (j=0;j<i;j++) free(X[j]+j+1);
	  free(X);
	  free(SortedVertices);
	  return;
	}
      X[i]-=i+1;
    }
  DiffLambda=(double ***)malloc((N-2)*sizeof(double **));
  if (DiffLambda==NULL)
    {
      DoRelax=0;
      for (i=0;i<N-1;i++) free(X[i]+i+1);
      free(X);
      free(SortedVertices);
      return;
    }
  for (i=0;i<N-2;i++)
    {
      DiffLambda[i]=(double **)malloc((N-2-i)*sizeof(double *));
      if (DiffLambda[i]==NULL) 
	{
	  DoRelax=0;
	  for (j=0;j<i;j++) free(DiffLambda[j]+j+1);
	  free(DiffLambda);
	  for (j=0;j<N-1;j++) free(X[j]+j+1);
	  free(X);
	  free(SortedVertices);
	  return;
	}
      DiffLambda[i]-=i+1;
      for (j=i+1;j<N-1;j++) 
	{
	  DiffLambda[i][j]=(double *)malloc((N-1-j)*sizeof(double));
	  if (DiffLambda[i][j]==NULL) 
	    {
	      DoRelax=0;
	      for (k=0;k<i;k++) 
		{
		  for (h=k+1;h<N-1;h++)free(DiffLambda[k][h]+h+1);
		  free(DiffLambda[k]+k+1);
		}
	      for (h=0;h<j;h++)free(DiffLambda[i][h]+h+1);
	      free(DiffLambda[i]+i+1);
	      free(DiffLambda);
	      for (j=0;j<N-1;j++) free(X[j]+j+1);
	      free(X);
	      free(SortedVertices);
	      return;
	    }
	  DiffLambda[i][j]-=j+1; 
	  for (k=j+1;k<N;k++) DiffLambda[i][j][k]=0;
	}
    }
  SaveCoeff=(double **) malloc((N-1)*sizeof(double *));
  if (SaveCoeff==NULL)  
    {
      DoGain=0;
      return;
    }
  for (i=0;i<N-1;i++)
    {
      SaveCoeff[i] = (double *) malloc((N-1-i)*sizeof(double));
      if (SaveCoeff[i]==NULL)  
	{
	  DoGain=0;
	  for (j=0;j<i;j++) free(SaveCoeff[j]+j+1);
	  free(SaveCoeff);
	  return;
	}
      SaveCoeff[i]-=i+1;
    }
  if (N%2)  MaxNbCircuits=(N*N*N-N)/24;
  else  MaxNbCircuits=(N*N*N-4*N)/24;
  circuits=(Circuit *) malloc(MaxNbCircuits*sizeof(Circuit));
  if (circuits==NULL) 
    {
      DoGain=0;
      for (i=0;i<N-1;i++) free(SaveCoeff[i]+i+1);
      free(SaveCoeff);
      for (i=0;i<N-2;i++) 
	{
	  for (j=i+1;j<N-1;j++)free(DiffLambda[i][j]+j+1);
	  free(DiffLambda[i]+i+1);
	}
      free(DiffLambda);
      for (j=0;j<N-1;j++) free(X[j]+j+1);
      free(X);
      free(SortedVertices);
    }
}
/********************* end of InitRelax ************************/



/***********  initializations before each relaxation ************/

void InitOneRelax(int * LocalVertices, int LocalN)
{
  int i,j;

  PositiveTriangle=0;
  j = 0;
  for (i = 0; i < N; i++)
    if (IsOutBegSec[i]) SortedVertices[j++] = i;
  for(i = 0; i < LocalN-1; i++) 
    for(j = i + 1; j < LocalN; j++) 
      if (tournament[SortedVertices[i]][SortedVertices[j]] > 0)
	PositiveTriangle +=
	  tournament[SortedVertices[i]][SortedVertices[j]];
  OldL = 0;
}
/********************** end of InitOneRelax ********************/

int * OutBegSecBenefit;  /* contains the address of an array; 
			    if i is a given  vertex, OutBegSecBenefit[i]
			    gives the sum of tournament[i][j] for  
			    vertices j  which are not in the current 
			    beginning section.
			    OutBegSecBenefit[i] also provides the following 
			    difference :
			    value of an order where i is just after the set of
			    vertices which are not in the current beginning 
			    section  - value of the same order except that i 
			    is just before the set of  vertices which are 
			    not in the   current beginning section*/
int ** benefit; /*  contains the address of a matrix;
		    if 1 <= i < j <(the size of the current beginning section,
		    benefit[i][j] gives the sum of 
		    tournament[CurrentOrder[i]][CurrentOrder[k]]
		    for i < k <= j.
		    It also gives the signed gain we would have if we 
		    take the vertex which is at the place i and we put it
		    just after the vertex which is at the place j*/


/************************ allocates some useful arrays *****************/

void initialization()
{ 
  int i;

  CurrentOrder = (int *)malloc(N*sizeof(int));
  if (CurrentOrder == NULL) 
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
  for (i = 0; i < N; i++) CurrentOrder[i] = i;
  BestOrder = (int *)malloc(N*sizeof(int));
  if (BestOrder == NULL) 
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
  for (i = 0;i < N;i++) BestOrder[i] = i;
  bound = OrderValue(BestOrder,N);
  score = (int *) malloc(N*sizeof(int));
  if (score == NULL) 
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    } 
  place = (int *) malloc(N*sizeof(int));
  if (place == NULL) 
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
}
/******************** end of  initialization ***************************/



/************** computes the value of a given order *********************/

int OrderValue(int * order,  int q)
{
  int i,j;
  int weight = 0;

  for (i = 0; i < q - 1; i++)
    for (j = i + 1; j < q; j++)
      if (tournament[order[i]][order[j]] < 0)
	weight +=tournament[order[j]][order[i]];
  return weight;
}
/************************* end of OrderValue **************************/



/************* copies an array of addres Order, of length p, 
               in an array of address OrderCopy ***********************/

void  CopyOrder(int * OrderCopy,int * Order, int q)
{ 
  int i;

  for (i = 0; i < q; i++) OrderCopy[i] = Order[i];
}
/************************* end of CopyOrder **************************/
	


/************* puts the contain of BestOrder ***************************/
void PutOrder()
{
  int i;
  
  for (i = 0; i < N-1; i++) printf("%d > ", BestOrder[i]);
  printf("%d", BestOrder[N - 1]);
  printf("\n");
}
/************************* end of PutOrder ***************************/



/************** this function performs the tests named 
        lex, ham, Smoves, OSmoves in the paper **************************/

char IsValid(int vertex,int index)
     /* when this function is called, the current beginning section
	is in CurrentOrder between the indexes 0 and "index"-1 and 
	the possibility to put "vertex" at the index "index" is considered;
	this function return 0 if it finds that is is impossible 
	to have an optimal order beginning with the current beginning
	function and having the vertex "vertex" at the index "index".
	Otherwise, it returns 1*/
{
  int j;
  int i,k;
  int CurrentBenefit;
  
  /*test named "ham"*/
  if ((index >= 1) && (tournament[CurrentOrder[index-1]][vertex] < 0))  
    return 0;

 /* special case of OSmoves, when we consider the move of "vertex" at 
    the end of the total order*/
  if (OutBegSecBenefit[vertex] < 0)  return 0;

  CurrentBenefit = 0;
  for (j = index-1; j >= 0; j--)
    {  
      /* special case of Smoves:  we consider the move of "vertex"
         before the vertex which is now at the index j*/
      CurrentBenefit += tournament[vertex][CurrentOrder[j]];
      if (CurrentBenefit > 0) return 0;
      else if (!CurrentBenefit)
	if (CurrentOrder[j]>vertex)  return 0;
    }
 
  for (j = 0; j < index-1; j++)
    { 
      /* special case of OSmoves: we consider the move of the vertex 
	 which is now at the index j after the beginning section obtained
	 by adding "vertex" after the current beginning section*/
      CurrentBenefit = benefit[j][index - 1] 
	+ tournament[CurrentOrder[j]][vertex];
      if (CurrentBenefit < 0)  return 0;
      else if (!CurrentBenefit)
	  if (CurrentOrder[j] > CurrentOrder[j + 1])  return 0;
    }

  /* special case of lex: if the weight between the last vertex 
     of the current beginning section and "vertex" is equal to 0,
     and if the "name" of the last vertex of the current beginning
     section is greater (for the lexicographic order) than the "name"
     of "vertex", we do not study the beginning section obtained  by 
     adding "vertex" after the current beginning section */
  if ((index>=1)&&(!tournament[CurrentOrder[index-1]][vertex])
      &&(vertex < CurrentOrder[index-1]))  return 0;

  /*test OSmoves, only performed for intervals having no more than
   3 vertices; for more vertices, we observed that it nearly never cut*/
  CurrentBenefit = OutBegSecBenefit[vertex];
  for(i = index - 1; (i >= 0) && (i >= index-3); i--)
    {
      CurrentBenefit+=OutBegSecBenefit[CurrentOrder[i]] -
	tournament[CurrentOrder[i]][vertex];
      if (CurrentBenefit < 0)  return 0;
    }

  /*test Smoves, only performed for intervals having no more than
   3 vertices; for more vertices, we observed that it nearly never cut*/
  for (j = 1;j < index-1;j++)
    {
      CurrentBenefit = benefit[j][index-1]+tournament[CurrentOrder[j]][vertex];
      for(i = j-1;i>=0;i--)
	{
	  if ((j - i > 3) && (index - j > 4)) break;
	  CurrentBenefit += benefit[i][index - 1] +
	    tournament[CurrentOrder[i]][vertex] - benefit[i][j];
	  if (CurrentBenefit < 0)   return 0;
	  else if (!CurrentBenefit)
	    /* lexicographic test */
	      if (CurrentOrder[i]>CurrentOrder[j+1]) return 0;
	}
    }
  return 1;
}
/************************ end of IsValid *****************************/




/*********** updates the values of benefit and OutBegSecBenefit when 
we add "vertex" at the end of the current beginning section (at the index p),
to obtain the new current beginning section******************************/

void  AddToBenefits(int vertex,int p)
{
  int i;

  if (p > 0) benefit[p-1][p] = tournament[CurrentOrder[p-1]][vertex];
  for (i = 0; i < p-1; i++) 
    benefit[i][p] = benefit[i][p-1]+tournament[CurrentOrder[i]][vertex];
  for (i = 0;i < N;i++) 
    OutBegSecBenefit[i] -= tournament[i][vertex];
}
/********************** end of AddToBenefits **************************/




/****************updates the values of OutBegSecBenefit when we 
retry "vertex" at the end of the current beginning section to 
obtain the new current beginning section*******************************/

void DeleteFromBenefits(int vertex,int p)
{
  int j;

  for (j = 0; j < N; j++) 
     OutBegSecBenefit[j] += tournament[j][vertex];
}
/********************* end of DeleteFromBenefits ***********************/



/*********** allocates the matrix for benefit, the array for 
  OutBegSecBenefit and initialize this array ***************************/

void InitBenefits()
{
  int i,j;

  benefit = (int **)malloc(N*sizeof(int *));
  if (benefit==NULL) 
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
  for (i = 0;i < N;i++)
    {
      benefit[i] = (int *)malloc((N-i-1)*sizeof(int));
      if (benefit[i] ==NULL) 
	{
	  printf("sorry : the tournament is too big\n");
	  exit(1);
	}
      benefit[i] -= i + 1;
    }
  OutBegSecBenefit = (int *)malloc(sizeof(int)*N);	
  if (OutBegSecBenefit==NULL)
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
  for (i = 0; i < N; i++) 
    {
      OutBegSecBenefit[i] = 0;
      for (j = 0;j < N; j++)
	OutBegSecBenefit[i] += tournament[i][j];
    }
}
/************************** end of InitBenefits *************************/



/************* this function is used in noising and descent; 
it moves the vertex which is at the place OldPlace in order, 
to the place NewPlace; it updates the array "place" ********************/

void change(int OldPlace,int NewPlace,
	    int * order, int * place)
{ 
  
  int MovingVertex;
  int i;

  if (NewPlace == OldPlace) return;
  MovingVertex = order[OldPlace];
  if (NewPlace  <  OldPlace)
    for (i = OldPlace; i > NewPlace; i--)  
      {
	order[i] = order[i-1];
	place[order[i]] = i;
      }
	
  else
    for (i = OldPlace; i < NewPlace; i++)  
      {
	order[i] = order[i+1];
	place[order[i]] = i;
      }
  order[NewPlace] = MovingVertex;
  place[MovingVertex] = NewPlace;
}    
/***************************** end of change ****************************/



/***************** computes the transitivity index **********************/
double ComputeIndex()
{
  int i,j;
  int sum;
  int max;

  sum = (N*(N-1)*(N-2))/6;
  for (i = 0;i < N;i++) score[i] = 0;
  for (i = 0;i < N-1;i++)
      for (j = i+1;j < N;j++)
	{
	  if (tournament[i][j]>=0) score[i]++;
	  else score[j]++;
    }
  for (i = 0; i < N; i++) 
      sum -= (score[i]*(score[i]-1))/2;
  if (N%2) max = (N*N*N-N)/24;
  else max = (N*N*N-4*N)/24;
  return 1 - (double)sum/(double)max;
}
/*********************** end of ComputeIndex ****************************/
  
int **  ToExam;

void exchange(int,int);

int compteur=0;

void BB(int p,int BegSecValue) 
/* EHS=somme des poids comptes negativement des arcs de HS*/
/* p est la taille de la section commencante correspondante
ES est ce que coute la section commencante*/
{
  int V,i,j;
  int NewBegSecValue;
  int vertex;
  int q;
  int sigma;
  char ExistResult;

  if(p==N-1)
    {
      if (BegSecValue < bound)
	{
	  bound=BegSecValue;
	  CopyOrder(BestOrder, CurrentOrder, N);
	  printf("\nAn order of value %d is found.\nIt is :\n"
		 ,bound);
	  PutOrder();
	}
      return;
    }
  if (DoExists) ExistResult = exist(p, BegSecValue); /* BegSec Test */
  else ExistResult = 0;
  if (ExistResult == 1)  return;
  q=N-p;
  if (DoSigma)
    {
      sigma = ComputeSigma(CurrentOrder + p, q); /* sigma test */
      if ((sigma + BegSecValue >= bound) || (!sigma))  return;
    }
  V = descent(q, CurrentOrder + p);
  if (BegSecValue + V < bound)
    {
      bound = BegSecValue + V;
      CopyOrder(BestOrder,CurrentOrder,N);
      printf("\nAn order of value %d is found.\nIt is :\n"
	     , bound);
      PutOrder(); 
    }
  if ((DoRelax) && (q > 2))      /*relax test */
    if ((ExistResult==2)&&(bound-BegSecValue<=LowerBoundBegSec())) return;
    else if (relax(CurrentOrder+p,q,bound-BegSecValue,V)) 
    	return;
  CopyOrder(ToExam[p],CurrentOrder+p,q);
  for (i = 0; i < q; i++)
    {
      vertex = ToExam[p][i];
      if (!IsValid(vertex, p))  continue; /* ham, Smoves, OSmoves, lex tests */
      NewBegSecValue = BegSecValue;
      for (j = p; j < N; j++)
	if (tournament[vertex][CurrentOrder[j]] < 0)
	  NewBegSecValue += tournament[CurrentOrder[j]][vertex];
      exchange(vertex, p);
      AddToBenefits(vertex, p);
      IsOutBegSec[vertex] = 0;
      BB(p+1,NewBegSecValue);
      IsOutBegSec[vertex] = 1;
      DeleteFromBenefits(vertex, p);
    }
}
/*************************** end of BB *****************************/



/************ initialization of the branch and bound method *********/

void InitBB(double TransitivityIndex)
{
  int i;

  IsOutBegSec = (char *)malloc(N*sizeof(char));
  if (IsOutBegSec==NULL) 
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
  for(i=0;i<N;i++)  IsOutBegSec[i]=1;
  ToExam=(int **) malloc(N*sizeof(int *));
  if (ToExam==NULL) 
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
  for (i=0;i<N;i++)
    {
      ToExam[i]=(int *) malloc((N-i)*sizeof(int));
      if (ToExam[i]==NULL)  
	{
	  printf("sorry : the tournament is too big\n");
	  exit(1);
	}
    }
  InitBenefits();
  InitRelax(TransitivityIndex);
  InitHead();
}
/************************* end of InitBB **************************/




/******************* exchange the vertex "vertex" with the vertex 
                     which is at the place index *******************/

void exchange(int vertex, int index)
{
  int i=index;

  while (CurrentOrder[i]!=vertex) i++;
  CurrentOrder[i]=CurrentOrder[index];
  CurrentOrder[index]=vertex;
}
/************************ end of exchange ***************************/

/*********  For the noisy method ***********/
void  firstPreparation();
double date(); 
void recap();
void ajustRate();
void putBest();
void abstractRecap();
void changeAuto(int ,int);
int aDescent();
void fourRounds();
void noisedPlace(int, int *);
void permutation (int *);
int bestPlace(int, int *);
void FillNoiseTable(void);
double averageWeight();
void initCurrentOrder();

int currentVal;
double rate;	
int nbAcept, nbAceptNul;
int round;

int middleMin;
double theTime;	
double beginningDate;
int sup;
double increasing;
double bestRate;
double maxRate;
double sumRates = 0;
double sumImprovements = 0;
double beginningRate, lastRate;
int localBest;
double sumRates;
int nbRates = 0;
long  MaxRandomTable;
double * RandomTable;
long IndexRandomTable=0;
double purcent, maxPurcent = 0.35;

int * permut;
void noisedPlace(int , int *);
int begin;
int nbAcept;

void noising(double t)
{
  double presentDate, oldDate, noisingDuration;
  double cycleDuration, elapsedTime = 0;
  int k;
  double ratio;
  int nb4rounds = 10;
  double finishingDate;
  int doEnd = 0;
  int sumCoeffNoising = 0;
  int coeffNoising = 1;
  double sumAjustedRate = 0;

  permut = score; //we use score to save some place
  FillNoiseTable();
  round = N*(N-1);
  theTime = t;
  firstPreparation();
  finishingDate  =  beginningDate + theTime;
  oldDate = date();
  elapsedTime = oldDate - beginningDate;
  localBest = sup;
  while (elapsedTime < theTime)
    {
      initCurrentOrder();
      currentVal = aDescent();
      if (currentVal < sup) sup = currentVal;
      sumImprovements = sumRates = 0; 
      for (k = 0; k < nb4rounds; k++)    
	{
	  rate = k * (lastRate - beginningRate) / 
	    nb4rounds + beginningRate;
	  fourRounds();
	  putBest();
	  recap();
	  if (doEnd)
	    {
	      presentDate = date();
	      if (presentDate > finishingDate)  break;
	      cycleDuration = (presentDate - oldDate) 
	 	/ (k + 1);
	      nb4rounds = (finishingDate - oldDate) 
		/ cycleDuration;
	    }
	}
      presentDate = date();
      elapsedTime = presentDate - beginningDate;
      noisingDuration = presentDate - oldDate;
      if (doEnd) break;
      if (elapsedTime > theTime) break;
      abstractRecap();
      if (begin)
	{
	  ajustRate();
	  sumCoeffNoising += coeffNoising;
	  sumAjustedRate += rate * coeffNoising;
	  coeffNoising *= 2;
	  beginningRate = sumAjustedRate / 
	    sumCoeffNoising;
	}
      else beginningRate = bestRate;
	    
      if (4*noisingDuration >  theTime - elapsedTime)
	{
	  nb4rounds *= (theTime - elapsedTime) /  
	    noisingDuration;
	  doEnd = 1;
	}
      else nb4rounds *= 2;
      if (nb4rounds <= 0) break;
       oldDate = presentDate;
    }
  CopyOrder(CurrentOrder, BestOrder, N);
}

void firstPreparation()
{
  int nbDesc = 0;
  double ratio = 0;
  int valDescente;

  
  beginningDate = date();
  lastRate = 0.0;
  begin = 0;
  initCurrentOrder();
  rate = averageWeight(N, CurrentOrder);
  sup = bound = currentVal = aDescent();
  CopyOrder(BestOrder, CurrentOrder, N);
  nbDesc++;
 while ((nbDesc < 10) && (date() - beginningDate <  theTime/10))
    {
      initCurrentOrder();
      currentVal = aDescent();
      nbDesc++;
      if (currentVal < bound)
	{
	  sup = bound;
	  CopyOrder(BestOrder, CurrentOrder, N);
	}
    }
  ratio = (date() - beginningDate) /  (nbDesc * theTime);
  increasing = 1 + 10 * ratio;
  if (increasing < 1.1) increasing = 1.1;
  else if (increasing > 2) increasing = 2;
  fourRounds();
  putBest();
  if (currentVal < bound)
    {
      sup = bound;
      CopyOrder(BestOrder, CurrentOrder, N);
    }
  if (purcent < maxPurcent)
    do
      { 
	rate *= increasing;
	fourRounds();
	putBest();
      }while (purcent < maxPurcent);
  else 
    {
      do
	{ 
	  rate /= increasing;
	  fourRounds();
	  putBest();
	}while (purcent > maxPurcent);
      rate *= increasing;
    }
  beginningRate = rate;
  initCurrentOrder();
  currentVal = aDescent();
  if (currentVal < sup)  sup = currentVal;
  bestRate = rate;
}

void ajustRate()
{
  rate = bestRate;
  middleMin = (sup + bound) / 2;
  CopyOrder(CurrentOrder, BestOrder, N);
  do 
    {
      rate *= increasing;
      fourRounds();
      putBest();
      if (date() - beginningDate > theTime) return;
    }while ((currentVal < middleMin) && 
	    (purcent < maxPurcent));
  rate /= increasing;
}

double date()
{
  long  sec, misec ;
  struct rusage buftime;
  
  getrusage(RUSAGE_SELF, &buftime);
  sec = buftime.ru_utime.tv_sec;
  misec = buftime.ru_utime.tv_usec;
  return  (double)sec + (double)misec/1000000;
}


void recap()
{
  int improvement;

  if  (currentVal < localBest)
    {
      localBest = currentVal;
      sumImprovements += 1;
      sumRates +=  rate;
    }
}


void abstractRecap()
{
  
  if (sumImprovements > 0)
    {
      bestRate = sumRates / sumImprovements;
      begin = 1;
    }
}

void putBest()
{
  if  (currentVal < bound)
    {
      bound = currentVal;
      CopyOrder(BestOrder, CurrentOrder, N);
    }
}
 

int aDescent()
{
  int ind1,ind2;
  int delta;
  int i;
  int vertex;
  int notFound = 0;
  int index = 0;
 
  permutation(permut);  
  vertex = permut[0];
  ind1 = place[vertex];
  do
    { 
      delta= bestPlace(ind1,&ind2);
      if (delta<0)
	{
	  changeAuto(ind1,ind2); 
	  notFound=0;
	}
      else notFound++;
      index++;
      if (index == N)
	{
	  permutation(permut);
	  index = 0;
	}
      vertex = permut[index];
      ind1=place[vertex];
    } while(notFound < N); 
  return  OrderValue(CurrentOrder, N);
}


int bestPlace(int oldPlace, int *addressNewPlace)
{
  int i;
  int minimum = 0;
  int delta = 0;

  *addressNewPlace = oldPlace;
  for (i = oldPlace - 1; i>=0;i--)
    {
      delta +=tournament[CurrentOrder[i]][CurrentOrder[oldPlace]];
      if (delta < minimum)
	{
	  minimum = delta;
	  *addressNewPlace = i;
	}
    }
  delta = 0;
  for (i = oldPlace+1; i < N; i++)
    {
      delta += tournament[CurrentOrder[oldPlace]][CurrentOrder[i]];
      if (delta < minimum)
	{
	  minimum = delta;
	  *addressNewPlace = i;
	}
    }
  return minimum;
}

void fourRounds()
{
  int entreDesc, k;
  int ind1, ind2;
  int vertex;
  int nbTries;

  nbTries = nbAcept = 0;
  permutation(permut);
  nbAcept = nbAceptNul = 0;
  for (entreDesc = 0; entreDesc < 4; entreDesc++)
    {
      for (k = 0; k < N; k++)
	{
	  vertex = permut[k];
	  ind1=place[vertex];
	  noisedPlace(ind1, &ind2);
	  if (ind1 != ind2)  changeAuto(ind1, ind2);
	}
    }
  currentVal = aDescent();
  nbTries = 4 * N * (N - 1);
  purcent = (double)nbAcept / (double)nbTries;
}

  
void noisedPlace(int oldPlace, int *addressNewPlace)
{
  int i, j;
  double minimum = 0;
  double deltaNoised = 0, noise;
  int * tournamentLine;
  int delta = 0;

  *addressNewPlace = oldPlace;
  tournamentLine = tournament[CurrentOrder[oldPlace]];
  for (i = oldPlace - 1; i >= 0 ;i--)
    {
      delta -= tournamentLine[CurrentOrder[i]];
      noise = RandomTable[IndexRandomTable++]*rate;
      if (IndexRandomTable == MaxRandomTable) IndexRandomTable  = 0;
      deltaNoised -= tournamentLine[CurrentOrder[i]] + noise;
      if (deltaNoised < 0)  nbAcept++;
      if (deltaNoised < minimum)
	{
	  minimum=deltaNoised;
	  *addressNewPlace=i;
	}
    }

  deltaNoised=0; 
  delta = 0;

  for (i= oldPlace+1; i < N; i++)
    {
      delta += tournamentLine[CurrentOrder[i]];
      noise = RandomTable[IndexRandomTable++]*rate;
      if (IndexRandomTable == MaxRandomTable) IndexRandomTable  = 0;
      deltaNoised += tournamentLine[CurrentOrder[i]] + noise;
      if (deltaNoised < 0)  nbAcept++;
      if (deltaNoised < minimum)
	{
	  minimum = deltaNoised;
	  *addressNewPlace = i;
	}
    }
}

void changeAuto(int oldPlace,int newPlace)
{ 
  
  int movingVertex;
  int i;

  if (newPlace == oldPlace) return;
  movingVertex = CurrentOrder[oldPlace];
  if (newPlace < oldPlace)
    for (i = oldPlace; i > newPlace; i--)  
      {
	CurrentOrder[i] = CurrentOrder[i - 1];
	place[CurrentOrder[i]] = i;
      }
	
  else
    for (i = oldPlace; i < newPlace; i++)  
      {
	CurrentOrder[i] = CurrentOrder[i+1];
	place[CurrentOrder[i]] = i;
      }
  CurrentOrder[newPlace] = movingVertex;
  place[movingVertex] = newPlace;
}


void permutation (int * table)
{
  int  i, vertex, num;
 
  for (i = 0; i < N; i++) table[i]=i;
  for (i=0;i<N-1;i++) 
  {
    num = random() % (N - i);
    vertex = table[i];
    table[i]=table[i+num];
    table[i+num]=vertex;
  }
}

#include <time.h>

void FillNoiseTable()
{
  long i;

  /*the noising method uses the table RandomTable which contains 
    "real" random numbers chosen between 0 and 1*/
  MaxRandomTable=321596;
    srandom(time((void *)&i)); /*to initialize the random generator*/
  while(1)
    {
      RandomTable=(double *)malloc(MaxRandomTable*sizeof(double));
      if (RandomTable!=NULL) break;
      else 
	{
	  printf("the length of the table of random numbers "
		 "is divided by 2\n");
	  MaxRandomTable /= 2;
	}
    }
  for (i = 0; i < MaxRandomTable; i++)
    RandomTable[i] = (double)random() / (double)INT_MAX-0.5;
}


double averageWeight()
{
  int i, j;
  int sum = 0;

  for (i = 0; i< N-1 ; i++)
    for (j = i+1; j < N; j++) sum += abs(tournament[i][j]);
  if (sum == 0) 
    {
      printf("all the weights are equal to zero, "
	     "every order is optimal\n");
      exit(1);
    }
  return (double) 2*sum / (double)(N*(N-1)) ;
}

void initCurrentOrder()
{
  int  i, vertex, num;
 
  for (i = 0; i < N; i++) CurrentOrder[i] = i;
  for (i = 0; i < N-1; i++) 
  {
    num = random()%(N-i);
    vertex = CurrentOrder[i];
    CurrentOrder[i] = CurrentOrder[i+num];
    CurrentOrder[i + num] = vertex;
  }
  for (i = 0; i < N; i++) place[CurrentOrder[i]] = i;
}
/* end of noisy method*/


/* For the beginning sections tree */

/* a branch of the tree, from the root to a leaf, gives the set of vertices
   of a beginning encountered section, written in lexicographic order.
   To find the children of a node, one have to first take the first child with
   the address "child" and the other children as brothers of the first child,
   using adresses contain in "brother"*/ 
typedef struct section
{
  int vertex; /* a vertex of a beginning section */
  int value; /* when it is the end of a set of vertices corresponding 
		to an encountered  beginning section, it is equal to 
		the smallest value of encountered beginning sections on 
		these vertices. 
		Otherwise, it is equal to -1 */
  int LowerBound; /* when it is the end of a set S of vertices corresponding 
		to an encountered  beginning section, it is equal to the
		greatest lower bound (obtained by relaxation) for value of 
		orders on X-S. 
		Otherwise, it is equal to -1 */
  struct section * child; /* gives the first child */
  struct section * brother; /* gives the first brother */
}Section;


Section * memory; /* keeps the address of the node corresponding 
		     to the current beginning section */

Section *  create(Section *, char , int);
void FreeTree(Section *);

Section * head=NULL; /* the address of the root of the tree */

/* this function searches if a set S of vertices is already encountered as 
the set of vertices of a beginning section ; if not, this set is put in the 
tree and the function returns 0; it it is, the value of the current beginning 
section, BegSecValue, is compared to the stored value (in the fieldvalue);
if BegSecValue is the biggest, the function returns 1; if BegSecValue is
the smallest, the field value is replaced by BegSecValue and  the function
returns 0 if no relaxation has already been done on X-S or 2 otherwise*/

char exist(int p, int BegSecValue)
{
  int i,j;
  char IsInTree;
  Section * current;
  char  Ok;

  j=-1;
  IsInTree=1;
  current=head;
  for (i=0;i<p;i++)
    {
      while (IsOutBegSec[++j]);
      if (IsInTree)
	{
	  Ok=0;
	  if (current->child!=NULL)
	    {
	      if (current->child->vertex<=j)
		{
		  current=current->child;
		  if (current->vertex==j) Ok=1;
		  else
		    {
		      while (current->brother!=NULL)
			if (current->brother->vertex<=j)
			  {
			    current=current->brother;
			    if (current->vertex==j) Ok=1;
			  }
			else break;
		      if (!Ok)
			{
			  current=create(current,1,j); 
			  if (current==NULL) return 0;
			  IsInTree=0;
			}
		    }
		}
	      else
		{
		  current=create(current,2,j);
		  if (current==NULL) return 0;
		  IsInTree=0;
		}
	    }
	  else
	    {
	      current=create(current,2,j);
	      if (current==NULL)   return 0;
	      IsInTree=0;
	     }
	  if (i==p-1)
	    {
	      if (Ok)
		{
		  if  (current->value>=0) 
		    if (current->value<=BegSecValue)   
		      return 1;
		    else 
		      {
			current->value=BegSecValue;
			memory=current;
			if (current->LowerBound>=0)
			  return 2;
			else  return 0;
		      }
		  else
		    {
		      current->value=BegSecValue;
		      memory=current;
		      return 0;
		    }
		}
	      else current->value=BegSecValue;
	    }
	}
      else 
	{
	  current=create(current,2,j);
	  if (current==NULL) return 0;
	  if (i==p-1)   current->value=BegSecValue;
	}
    }
  memory=current;
  return 0;
}
/********************** end of exist ***************************/



/***************** this function create a new node 
              for the lexicographic tree ***********************/

Section *  create(Section * current, char message, int LeSommet)
{
  Section * pointer;

  pointer=(Section *)malloc(sizeof(Section));
  if (pointer==NULL)
    {
      DoExists=0;
      printf("On arrete les sections commencantes\n");
      return NULL;
    }
  pointer->vertex=LeSommet;
  pointer->value=-1;
  pointer->LowerBound=-1;
  pointer->child=NULL;
  if (message==1)  
    {
      pointer->brother=current->brother;
      current->brother=pointer;
    }
  else 
    {
      pointer->brother=current->child;
      current->child=pointer;
    }
  return pointer; 
}
/************************ end of create *************************/



/*************** initializes the root of the tree ****************/

void InitHead()
{
  head=(Section *)malloc(sizeof(Section));
  if (head==NULL) 
    {
      DoExists=0;
      printf("We stop beginning sections\n");
      return;
    }
  head->vertex=-1;
  head->value=-1;
  head->child=NULL;
  head->brother=NULL;
}
/************************ end of InitHead ************************/



/************* this function is called at the end 
       of relaxation in order to store the best value found
       for the dual problem **************************************/

void PutLowerBound(int value)
{
  if (value>memory->LowerBound) memory->LowerBound=value;
}
/****************** end of PutLowerBound *************************/



/*************** this function is called by BB in order
          to know the greatest  value already computed for
          the dual problem on the vertices which
          are not in the beginning section ***********************/

int LowerBoundBegSec()
{
  return memory->LowerBound;
}
/******************** end of LowerBoundBegSec ********************/



/****** frees the place allocated for the lexicographic tree *****/

void FreeTree(Section * current)
{
  if (current!=NULL)
    {
      FreeTree(current->child);
      FreeTree(current->brother);
      free(current);
    }
}
/************************ end of FreeTree ************************/

/*  For sigma */
void SortScores(int);


/*************** computes the parameter sigma ********************/
int  ComputeSigma(int * OutBegSec, int q)
{
  int i,j;
  int sigma = 0;

  for (i = 0; i < q; i++) score[i] = 0;
  for (i = 0; i < q - 1; i++)
    for (j = i+1; j<q; j++)
      if (tournament[OutBegSec[i]][OutBegSec[j]] > 0) score[i]++;
      else score[j]++;
  SortScores(q);
  for (i = 0; i < q; i++)
    if (score[i] > i) sigma += score[i] - i;
  return sigma;
}
/*********************** end of ComputeSigma *********************/



/*********** sort vertices on their  scores with Quick Sort********/
void SortScores(int q)
{
  int i,j;
  int key;
  int IndexBigChild;
  char end;

  score -= 1;
  for (i = 2;i <= q; i++)
    {
      key = score[i];
      j = i;
      while((j > 1) && (score[j/2] < key))
	{
	  score[j] = score[j/2];
	  j = j/2;
	}
      score[j] = key;
    }
  for (i = q; i >= 2; i--)
    {
      end = 0;
      key = score[i];
      score[i] = score[1];
      j = 1;
      while((2*j <= i - 1) && (!end))
      {
	if (2*j ==i-1) IndexBigChild = i - 1;
	else
	  if (score[2*j] < score[2*j + 1])
	    IndexBigChild=2*j + 1;
	  else IndexBigChild = 2*j;
	if (score[IndexBigChild] > key)
	  {
	    score[j]=score[IndexBigChild];
	    j=IndexBigChild;
	  }
	else  end=1;
      }
      score[j]=key;
    }
  score+=1;
}
/*********************** end of SortScores **********************/

struct rusage buftime;
long start_sec=1,start_misec,stop_sec,stop_misec;
double duration;
double start, stop;
/************** beginning of the main function **************/

int main()					      
{
  FILE * InFile; 
  char InputFile[30];
  double TransitivityIndex;
  double timeForNoising;
  double answer;
  int car;
  int i;
	
  getrusage(RUSAGE_SELF, &buftime);
  start_sec = buftime.ru_utime.tv_sec;
  start_misec = buftime.ru_utime.tv_usec;
  start =(double)start_sec + (double)start_misec/1000000;	
  printf("Please, give the name of the input file : ");
  gets(InputFile);
  InFile=fopen(InputFile,"r");			  
  if (InFile == NULL) 
    {
      printf("There is a problem with this file.\n");
      exit(1);
    }
  fscanf(InFile,"%d",&N);
  GetTournament(InFile);
  fclose(InFile);
  initialization();
  TransitivityIndex = ComputeIndex();
  printf("The transitivity index (number of 3-circuits "
	 "divided by the maximum number\non 3-circuits in a "
	 "tournament of same order) is %lf.\n",TransitivityIndex);
  timeForNoising = 0.0001*N*N*N*(1-TransitivityIndex + eps);
  if (timeForNoising < 0) timeForNoising = 0;
  printf("We choose to begin with a noising method which will last "
	 "%.2lf seconds\n", 
	 timeForNoising);
  printf("BUT, if you want more or less time, please, give it, in seconds, "
	 "\nand, if not, please type return ");
  fflush(stdout);
  while ((!isdigit(car = getc(stdin))) && (car != '\n'));
  if (car !='\n')
    { 
      ungetc(car, stdin);
      scanf("%lf", &answer);
      timeForNoising = answer;
    }
  noising(timeForNoising);
  printf("\nThe noising method obtains the value %d\n"
	 "for the order :\n",bound);
  PutOrder();
  for (i = 0; i < N; i++) CurrentOrder[i] = i;
  InitBB(TransitivityIndex);
  if (bound > 0) BB(0,0);
  PutResults();
  getrusage(RUSAGE_SELF, &buftime);
  stop_sec = buftime.ru_utime.tv_sec;
  stop_misec = buftime.ru_utime.tv_usec;
  stop =(double)stop_sec + (double)stop_misec/1000000;
  duration = stop -start;
  printf("total duration = %lf\n", duration);
  return 0;
}	
/****************** end  of the main function ********************/



/**************allocates and fills the matrix tournament *********/

void GetTournament(FILE * InFile)	
{ 			
  int i,j;   
	
  tournament=(int **)malloc(sizeof(int *)*N);	
  if (tournament==NULL)
    {
      printf("sorry : the tournament is too big\n");
      exit(1);
    }
  for (i=0;i<N;i++) 
    {
      tournament[i]=(int *)malloc(sizeof(int)*N);
      if (tournament[i]==NULL)
	{
	  printf("sorry : the tournament is too big\n");
	  exit(1);
	}
    }			
  for (i=0;i<N;i++)		
    for (j=0;j<=i;j++) 
      {
	fscanf(InFile,"%d",&tournament[i][j]);
	if ((i!=j)&&(tournament[i][j]!=1)&&(tournament[i][j]!=-1))
	    DoSigma=0;
	if (tournament[i][j]> MaxWeight) MaxWeight=tournament[i][j];
	if (-tournament[i][j]> MaxWeight) MaxWeight=-tournament[i][j];
	tournament[j][i]=-tournament[i][j];
      }
}		
/********************* end of GetTournament ************************/	

 

/******* this function indicates the results of the software ********/

void PutResults()
{
  int i,j,k;
  
  printf("\nThe best value is %d\nwhich is the value of the order\n",
	 OrderValue(BestOrder,N));
  PutOrder();
}	
/********************* end of PutsResuls ************************/