/*****************************************************************************************
 *
 * Author : Olivier Elemento, elemento@lirmm.fr
 *
 * LICENSE : Permission is granted to copy, modify, distribute this source for any purpose. 
 * INPUT   : distance file in PHYLIP format
 * SORTIES : a single tree in Newick format
 ****************************************************************************************/
#include <stdio.h> 
#include <stdlib.h>
#include <limits.h>
#include <string.h>
#include <assert.h>
#include <time.h>
#include <sys/timeb.h>

#define max(x,y) (x>y?x:y)
#define min(x,y) (x<y?x:y)


typedef struct _tree {

  
  int maxnodes;

  int numnodes;
  
  int numintnodes;
  
  int numleaves;
  
  int maxintnodes;
  
  int maxleaves;
  
  /** precedent de chaque node **/
  int* nodes;
  
  /** nodes fils (navigation, traverse) **/
  int** children;

  struct _tree* next;
} tree;



typedef struct _SEQS {

  /** séquences **/
  short* sequences;

  /** longueur du locus **/
  short size;
  
  /** pointeur sur le prochain locus **/
  struct _SEQS* next;

} SEQS;

void    printCostMatrix(FILE* fp, float** dists, SEQS* seqs); 

void    printDistances(FILE* fp, float** dists, int nb_spec);
float** loadDistances(char* file, int* nb_spec);
void    showDistances(float** dists, int nb_spec);
SEQS*   allocSeqs(short numleaves);
void    initSeqs(SEQS* seqs, short numleaves);
float   getCost(int i, int l, SEQS* seqs, float** dists);
void    mergeSeqs(SEQS* seqs, short u, short posi, short posj);
void    mergeSequences(SEQS* seqs, short* u, short i, short l);


void    printSeqs(SEQS* loc);
tree*   allocTree(int T);
float** allocMatrix(int nb_species);
void    saveTree(tree* t, int numnode, char* str);

void    showAllNodes(tree* t);
void    duplicateByUneqRecomb(tree* t, int i, int j);
void    printTree(tree* t);
void    updateDistances(SEQS* seqs, short* u, int posi, int posl, float** dists);
short   updateTree(tree* t, SEQS* seqs, short posi, short posj);
char*   get_parameter(int argc, char** argv, char* param);
void    printOrderFile(int n, char* s_orderFile);
void    showCostMatrix(float** dists, SEQS* seqs) ;

float   getWorstPairCost(int i, int l, SEQS* seqs, float** dists) ;
float   getPairCost(
		  float** crdata,
		  float** data,
		  int     p, 
		  int     l,
		  int     lastopti, 
		  int     lastoptl, 
		  SEQS*   seqs,
		  short*  u, 
		  short*  v,
		  short*  w,
		  short*  uv);
void    getMinMoyCost(int i, int l, SEQS* seqs, float** dists, float* min, float* moy); 
void    fprintTree(tree* t, char* outfile); 

/** distances **/
float** data;

/** distances corrigées **/
float** crdata;

char** names;

int verbose;

int correction;

int outorder;

int costtype;


  
int nb_species;

char* tracefile = "";
FILE* fp_trace;
int   trace = 0;


/** nombre de paires calculées **/
int nbpairscalc = 0;

/** nombre de doubles sommes **/ 
int nbdblesums  = 0;

int nbpairs = 0;


int main(int argc, char** argv) {
  /** tete de liste de SEQS **/
  SEQS* seqs;  
  
  /** conteneur de sequences **/
  float cost, min;

  /** cout max **/
  float maxcost;

  /** conteneur de sequences **/
  float moy, maxmoy;
  
   /** taille et taille maxi d'une dupli**/
  short d_maxlength;

  /** compteurs **/
  short i, l, j, p; //, t;

  int opti, optl, lastopti, lastoptl;
  
  /** noeud inter **/
  short* u;
  short* v;
  short* w;
  short* uv;

  /** arbre final **/
  tree* tr;

  char* orderfile;

  char* filename;
  char* outfile;


  /** raccourci **/
  short* s;

  struct timeb t1,t2;

  FILE* fp;
  
  ftime(&t1);

  verbose = 0;  
  if (strcmp(get_parameter(argc, argv, "-mode"), "verbose") == 0)
    verbose = 1;
  
  tracefile = get_parameter(argc, argv, "-trace");
  if (strlen(tracefile) > 0) {
    trace = 1;
    fp_trace = fopen(tracefile, "w");
  }
 
  
  filename = (char*)calloc(30, sizeof(char));
  outfile  = (char*)calloc(30, sizeof(char));
 

  if ((argv[1] == 0) || !(fp = fopen(argv[1], "r"))) {
  
    printf("Input file : ");
    scanf("%s", filename);
    printf("Output file : ");
    scanf("%s", outfile);
 
    
  } else {
    strcat(filename, argv[1]);
    fclose (fp);
  }

 
  
  
  data = loadDistances(filename, &nb_species);
  
  
  if (verbose)
    showDistances(data, 2 * nb_species - 1);

  
  /** alloue un nouveau seqs **/
  seqs = allocSeqs(nb_species);
  initSeqs(seqs, nb_species);


  /** nouvel arbre de nb_species feuilles **/
  tr = allocTree(nb_species);
  
  
  /** alloue une premiere matrice de couts dont la taille est égale au nombre final de noeuds **/
  crdata = allocMatrix(2 * nb_species - 1);

  /** rempli avec des -1 **/
  for (i=0; i<nb_species * 2 - 1; i++)
    for (j=0; j<nb_species * 2 - 1; j++)
      crdata[i][j] = -1.0;
	
  
  while (seqs->size != 2) {


    if (verbose) {
      showDistances(crdata, nb_species * 2 - 1);
      //getchar();
    }
    
    /** initialise le nombre de paires calculées **/
    nbpairscalc = 0;
    

    /** initialise le nombre de doubles sommes **/
    nbdblesums = 0;


    nbpairs = 0;
    
    /** raccourci **/
    s = seqs->sequences;
    
     /** longueur maxi de la fenetre **/
    d_maxlength = seqs->size / 2;
    
    /** cout mini **/
    maxcost = - SHRT_MAX;
    
    /** pour chaque longueur de duplication **/
    for (l=1; l<=d_maxlength; l++) {
      if (verbose)
	printf("L=%d\n", l);
      p = 0;
      while (p+l < seqs->size) {
	//for (p=0; p+2*l<=seqs->size; p++) {
	
	crdata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ] = getPairCost(crdata, data, p, l, lastopti, lastoptl, seqs, u, v, w, uv);
	
	
	if (verbose)
	  printf("Cout de la paire (%d, %d) = %3.1f\n", s[p], s[p+l], crdata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ]);

	cost = crdata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ];

	if (cost > maxcost) {

	  /** NOW      : determiner le cout de la duplication englobant la paire p,p+l, 
	                 sachant que celui-ci peut-etre inferieur ou cout max **/

	  
	  /** BACKWARD : calcule le cout de la duplication la plus a gauche si c'est possible **/
	  if (verbose)
	    printf("BACKWARD : calcule le cout de la duplication la plus a gauche si c'est possible et si elles n'ont pas ete calculées ...\n");
	  
	  for (i=max(0,p-l+1); i<p; i++) {
	    if (crdata[ min(s[i], s[i+l]) ][ max(s[i], s[i+l]) ] < -0.5)
	      crdata[ min(s[i], s[i+l]) ][ max(s[i], s[i+l]) ] = getPairCost(crdata, data, i, l, lastopti, lastoptl, seqs, u, v, w, uv);
	    else {
	      if (verbose)
		printf("la paire (%d, %d) a déjà été calculée : cout = %d\n", s[i], s[i+l], 
		       (int)crdata[ min(s[i], s[i+l]) ][ max(s[i], s[i+l]) ]);
	    }
	  }
	  
	  
	  /** FORWARD : calcule le cout des duplications a droite, et avance tant que le cout est superieur au cout max **/
	  if (verbose)
	    printf("FORWARD : calcule le cout des duplications a droite ...\n");

	  for (i=p+1; i<min(p+l, seqs->size-l); i++) { 
	   crdata[ min(s[i], s[i+l]) ][ max(s[i], s[i+l]) ] = getPairCost(crdata, data, i, l, lastopti, lastoptl, seqs, u, v, w, uv);
	  }

	  
	  /** CALCUL : fenetres en commencant par celle située à i = max(0,p-l+1) et finissant a p **/
	  if (verbose)
	    printf("calcule toutes les fenetres de i=%d à i=%d\n", max(0,p-l+1), p);
	  for (i=max(0,p-l+1); i<=min(p, seqs->size-2*l); i++) {
	    getMinMoyCost(i, l, seqs, crdata, &min, &moy);
	    if (verbose)
	      printf("Fenetre i=%d, l=%d, cout = %d\n", i, l, (int)min);
	    if (min > maxcost) {
	      opti    = i;
	      optl    = l;
	      maxcost = min;
	      maxmoy  = moy;
	    } else if ((min == maxcost) && (moy > maxmoy)) {
	      opti    = i;
	      optl    = l;
	      maxmoy  = moy;
	    }
	  }

	} else {
	  
	  /** met les couts a -1 entre p+1 et p+l-1 **/
	  for (i=p+1; i<min(p+l, seqs->size-l); i++) {
	    crdata[ min(s[i], s[i+l]) ][ max(s[i], s[i+l]) ]  = -1.0;
	  }
	  
	  if (verbose)
	    printf("Ca ameliore pas ... saute de l\n");
	  
	}
	
	p += l;
	
      }
    }

    if (verbose)
      printf("Meilleure duplication trouvée : opti = %d (%d), optl = %d pour un cout = %d\n", opti, s[opti], optl, (int)maxcost);
    
    lastopti = opti;
    lastoptl = optl;

    if (verbose) {
      printf("Matrice de couts modifée = \n");
      showCostMatrix(crdata, seqs);
      getchar();
    }

    /** verif a tejer **/
    if (nbpairs != (nbpairscalc + nbdblesums)) {
	printf("OOUPS ...\n");
	exit(0);
    }
	    

    if (trace)
      fprintf(fp_trace, "%4d %4d\n", nbpairscalc, nbdblesums);
    //printCostMatrix(fp_trace, crdata, seqs);
    
    

    /** cree une fenetre de nodes internes w **/
    w = (short*)malloc(optl*sizeof(short));
    
    /** cree les nodes internes (u = ensemble des nouveaux noeuds) **/
    for (i=0; i<optl; i++) {
      w[i] = updateTree(tr, seqs, opti+i, opti+optl+i);
      if (verbose)
	printf("w[%d] = %d\n", i, w[i]);
    }
    /** met a jour les distances entre les sequences de u et les autres dans le tableau **/
    updateDistances(seqs, w, opti, optl, data);
    
    
    /** stocke les u et les v **/
    u = (short*)malloc(optl*sizeof(short));
    memcpy(u, seqs->sequences+opti, optl * sizeof(short));
    
    v = (short*)malloc(optl*sizeof(short));
    memcpy(v, seqs->sequences+opti+optl, optl * sizeof(short));
    
    uv = (short*)malloc(2*optl*sizeof(short));
    memcpy(uv, u, optl * sizeof(short));
    memcpy(uv+optl, v, optl * sizeof(short));
    
    
    for (i=0; i<optl; i++) {
      
      if (verbose)
	printf("u[%d] = %d\n", i, u[i]);
    }

    for (i=0; i<optl; i++) {
      
      if (verbose)
	printf("v[%d] = %d\n", i, v[i]);
    }

   
    //getchar();
    
    
    
    
    if (verbose) {
      printf("Ancienne liste de seqs = \n");
      printSeqs(seqs);
    }
    
    /** merge le vecteur de feuilles **/
    mergeSequences(seqs, w, opti, optl);
    
    //mergeSequences(seqs, u, opti, optl);
  
    if (verbose) {
      printf("New liste de seqs = \n");
      printSeqs(seqs);
    }
    
    if (verbose) {

      printf("seqs->size = %d\n", seqs->size);
    }
  }

  if (verbose)
    printf("seqs->sequences[0] = %d, seqs->sequences[1] = %d ...\n", seqs->sequences[0], seqs->sequences[1]);

  tr->children[0][0] = seqs->sequences[0];
  tr->children[0][1] = seqs->sequences[1];
  
  tr->nodes[seqs->sequences[0]] = 0;
  tr->nodes[seqs->sequences[1]] = 0;
  if (verbose) {
    getchar(); printf("tree = \n");
    showAllNodes(tr);
  }
  

  if (strlen(outfile) == 0)  
    printTree(tr);
  else
    fprintTree(tr, outfile);

  if (outorder) 
    printOrderFile(nb_species, orderfile);

  ftime(&t2);
  
  if (trace)
    fprintf(fp_trace, "temps = %d\n", (int)((t2.time - t1.time)*1000 + (t2.millitm - t1.millitm)));
  

  if (trace)
    fclose(fp_trace);

  return 0;
}

void printOrderFile(int n, char* s_orderFile) {
  int i;
  FILE* fp = fopen(s_orderFile, "w");

  fprintf(fp, "%d\n", n);

  for (i=0; i<n; i++) {
    fprintf(fp, "%s\n", names[i]);
  }

  fclose(fp);
} 



void showAllNodes(tree* t) {
  int i;
  printf("node\tparent\tchildren (numintnodes=%d, numleaves=%d) \tlongueurs\n", t->numintnodes, t->numleaves);
  for (i=0; i<t->maxnodes; i++) {
    printf("%d\t%d\t%d,%d\n", i, t->nodes[i], t->children[i][0], t->children[i][1]);
  }
}

void printTree(tree* t) {
	char* str = (char*)calloc(4096, sizeof(char));
	saveTree(t, 0, str);	
	printf("%s;\n", str);
	free(str);
	str = 0;
}


/*
 * fprintTree : sauve un arbre dans un fichier
 *
 */
void fprintTree(tree* t, char* outfile) 
{
  FILE* fp;
  char* str = (char*)calloc(4096, sizeof(char));
  saveTree(t, 0, str);	
  
  fp = fopen(outfile, "w");
  

  fprintf(fp, "%s;\n", str);
  
  
  fclose(fp);

  free(str);
  str = 0;
  
}

/** avec les longueurs de branches **/
void saveTree(tree* t, int numnode, char* str) {
  //char tmp[10];
  if (t->children[numnode][0] == -1) {
    //sprintf(tmp, "rep%d", numnode - t->maxintnodes + 1);
    //strcat(str, tmp);
    strcat(str, names[numnode - t->maxintnodes]);
  }
  else {
    strcat(str, "(");
    saveTree(t, t->children[numnode][0], str);
    strcat(str, ",");
    saveTree(t, t->children[numnode][1], str);
    strcat(str, ")");
  }
    
}

/** alloue une matrice n*n **/
float** allocMatrix(int nb_species) {
  float** dists;
  int lig, col;
  
  dists = (float**)malloc(nb_species*sizeof(float*));
  
  for(lig=0; lig < nb_species; lig++) {
    dists[lig] = (float*)malloc(nb_species*sizeof(float));
    for (col=0; col<nb_species; col++)
      dists[lig][col] = 0.0;
  }
  
  return dists;
}

/** alloue un seqs **/
SEQS*  allocSeqs(short numleaves) {
  
  SEQS* seqs;

  seqs = (SEQS*)malloc(sizeof(SEQS));

  seqs->sequences = (short*)malloc(numleaves*sizeof(short));
  seqs->size      = numleaves;

  return seqs;
  
}

/** initialise un seqs **/
void initSeqs(SEQS* seqs, short numleaves) {
  short i;

  for (i=0; i<numleaves; i++)
    seqs->sequences[i] = i + numleaves - 1;
  if (verbose)
    printSeqs(seqs);
}


short updateTree(tree* t, SEQS* seqs, short posi, short posj) {
  int u;
  // nouvelle node  = t->numintnodes
  u = t->numintnodes;

  if (verbose) {
    printf("new node (u) = %d\n", u); getchar(); 
  }
  t->children[u][0] = seqs->sequences[posi];
  t->children[u][1] = seqs->sequences[posj];

  t->nodes[seqs->sequences[posi]] = u;
  t->nodes[seqs->sequences[posj]] = u;

  t->numintnodes++;

  return u;
  
}


/** 
 * merge un locus 
 * i position de la bande a supprimer
**/
void mergeSequences(SEQS* seqs, short* u, short i, short l) {
  
  short* liste;
  
  liste = (short*)malloc((seqs->size - l) * sizeof(short));

  /** copie de 0 a i vers 0 **/
  memcpy(liste, seqs->sequences, i*sizeof(short));
  

  /** copie de i a i+l vers i **/
  memcpy(liste + i, u, l * sizeof(short));

  
  /** copie de i+l vers i **/
  memcpy(liste + i + l, seqs->sequences + i + 2*l, (seqs->size - i - 2*l) *sizeof(short));
  

  seqs->sequences = liste;
  seqs->size -= l;
}


/** 
 * merge un seqs 
 * i, j paire a merger
 * u    nouvelle node
**/
void mergeSeqs(SEQS* seqs, short u, short posi, short posj) {
  int a,b;
  short* liste;

  if (verbose)
    printf("u=%d, posi=%d, posj=%d \n", u, posi, posj);
  /** nouvelle liste de seqs contenant seqs->size - 1 indiv **/
  liste = (short*)malloc((seqs->size - 1) * sizeof(short));
  
  /** recopier dans la nouvelle liste **/
  b = 0;
  for (a=0; a<seqs->size; a++) {
    if (a == posi) {
      liste[b] = u; 
      b++;
    } 
    else if (a == posj) {
      // ne fait rien
    }
    else {
      liste[b] = seqs->sequences[a];
      b++;
    }
    
  }

  seqs->sequences = liste;
  seqs->size--;

  if (verbose)
    printSeqs(seqs);
}


void printSeqs(SEQS* loc) {
  int i;
  for (i=0; i<loc->size; i++) {
    printf("%d ", loc->sequences[i]);
  }
  printf("\n");

}




/* 
 * calcule le cout d'une duplication (ajouté cr)
 * @param i debut de la duplication sur le locus
 * @param l taille de la duplication (fenetre / 2)
 * @param locus 
 * @param data distances
 */
float getCost(int i, int l, SEQS* seqs, float** dists) {
  int j;
  float cost = 0.0;
  short* s = seqs->sequences;

  if (verbose)
    printf("cost = ");
  for (j=i; j<i+l; j++) {
    if (verbose)
      printf("%5.4f + ", dists[ min(s[j], s[j+l]) ][ max(s[j], s[j+l]) ]);
    cost += dists[ min(s[j], s[j+l]) ][ max(s[j], s[j+l]) ]; //dists[ j ][ j+l ];
  }
  if (verbose)
    printf("\n");
    

  return cost / (float)l;

}

/* 
 * calcule le cout d'une duplication (ajouté cr) comme étant le cout de la plus mauvaise paire ...
 * @param i debut de la duplication sur le locus
 * @param l taille de la duplication (fenetre / 2)
 * @param locus 
 * @param data distances
 */
float getWorstPairCost(int i, int l, SEQS* seqs, float** dists) 
{
  
  int j;
  float mincost = 10000000.0;
  float cost    = 0.0;
  short* s      = seqs->sequences;

  if (verbose)
    printf("cost = ");
  for (j=i; j<i+l; j++) {
    if (verbose)
      printf("%5.4f + ", dists[ min(s[j], s[j+l]) ][ max(s[j], s[j+l]) ]);
    cost = dists[ min(s[j], s[j+l]) ][ max(s[j], s[j+l]) ]; //dists[ j ][ j+l ];
    if (cost < mincost)
      mincost = cost;
  }
  if (verbose)
    printf("\n");
    

  return mincost;

}






/* 
 * calcule le cout d'une duplication (ajouté cr) comme étant le cout de la plus mauvaise paire ...
 * @param i debut de la duplication sur le locus
 * @param l taille de la duplication (fenetre / 2)
 * @param locus 
 * @param data distances
 * @param min pointeur sur la valeur min (renvoyée)
 * @param moy pointeur sur la valeur moy (renvoyée)
 */
void getMinMoyCost(int i, int l, SEQS* seqs, float** dists, float* min, float* moy) 
{
  
  int j;
  float mincost = 10000000.0;
  float cost    = 0.0;
  float sum     = 0.0;
  short* s      = seqs->sequences;


  for (j=i; j<i+l; j++) {
    
    cost = dists[ min(s[j], s[j+l]) ][ max(s[j], s[j+l]) ]; 
    sum += cost;
    
    if (cost < mincost)
      mincost = cost;
  }

  *min = mincost;
  *moy = sum / l;
}







/**
 *
 *
 * @param locus = locus avant d'etre mergé
 * @param n     = nombre de sequences initialement
 */
void updateDistances(SEQS* seqs, short* u, int i, int l, float** dists) {
  int r, s, q, t;

  if (verbose) {
    printf("updateDistances avec i=%d et l=%d\n", i, l);
    printf("new nodes = ");
    for (r=0; r<l; r++)
      printf("%d ", u[r]);
    printf("\n");
  }
  //printLocus(locus);


  /** distances entre séquences intérieures à la duplication et séquences extérieures **/
  /** r = seq du locus **/
  /** q = seq ext      **/
  for (r=i, s=0; r<i+l; r++, s++) {

    /** sequences avant la partie dupliquée **/
    for (q=0; q<i; q++) {
      /** somme sur 2 des distances de la sequence aux deux sequences intérieures **/
      //printf("modifie la distance entre %d et %d\n", locus->sequences[r], locus->sequences[q]);
      
      /** q est nécéssairement avant r (voir les for) **/
      /** d(r, q) = d(r, q) + d(r+l, q) - d(r, r+l) **/
      
      if (!correction)
	dists[ min(u[s], seqs->sequences[q]) ][ max(u[s], seqs->sequences[q])  ] = 
	  (dists[ min(seqs->sequences[q], seqs->sequences[r])   ][ max(seqs->sequences[q], seqs->sequences[r])   ] + 
	   dists[ min(seqs->sequences[q], seqs->sequences[r+l]) ][ max(seqs->sequences[q], seqs->sequences[r+l]) ] ) / 2.0; 
      else 
	
       dists[ min(u[s], seqs->sequences[q]) ][ max(u[s], seqs->sequences[q])  ] = 
	(dists[ min(seqs->sequences[q], seqs->sequences[r])   ][ max(seqs->sequences[q  ], seqs->sequences[r])   ] + 
	 dists[ min(seqs->sequences[q], seqs->sequences[r+l]) ][ max(seqs->sequences[q  ], seqs->sequences[r+l]) ] -  
	 dists[ min(seqs->sequences[r], seqs->sequences[r+l]) ][ max(seqs->sequences[r]  , seqs->sequences[r+l]) ]) / 2.0; 
      
      if (verbose)
	printf("dist[%d][%d] = %5.4f\n", min(u[s] , seqs->sequences[q]),
	       max(u[s], seqs->sequences[q]),
	       dists[ min(u[s], seqs->sequences[q]) ][ max(u[s] , seqs->sequences[q]) ]);
	 /** -  dists[ seqs->sequences[r] ][ seqs->sequences[r+l] ]) / 2; **/

   
    }


    /** sequencess après la partie dupliquée **/
    for (q=i+2*l; q<seqs->size; q++) {
      //printf("modifie la distance entre %d et %d\n", seqs->sequences[r], seqs->sequences[q]);

      /** q est nécéssairement après r (voir les for) **/
      /** d(r, q) = d(r, q) + d(r+l, q) - d(r, r+l) **/
      
      
      /**
      dists[ min(u[s] , seqs->sequences[q]) ][ max(u[s] , seqs->sequences[q]) ] = 
	(dists[ min(seqs->sequences[r],   seqs->sequences[q]) ][ max(seqs->sequences[r],   seqs->sequences[q]) ] + 
	 dists[ min(seqs->sequences[r+l], seqs->sequences[q]) ][ max(seqs->sequences[r+l], seqs->sequences[q]) ] -
	 dists[ min(seqs->sequences[r], seqs->sequences[r+l]) ][ max(seqs->sequences[r]  , seqs->sequences[r+l]) ]) / 2.0; 
      **/
      
      if (!correction)
	dists[ min(u[s] , seqs->sequences[q]) ][ max(u[s] , seqs->sequences[q]) ] = 
	  (dists[ min(seqs->sequences[r],   seqs->sequences[q]) ][ max(seqs->sequences[r],   seqs->sequences[q]) ] + 
	   dists[ min(seqs->sequences[r+l], seqs->sequences[q]) ][ max(seqs->sequences[r+l], seqs->sequences[q]) ]) / 2.0; 
      else
	dists[ min(u[s] , seqs->sequences[q]) ][ max(u[s] , seqs->sequences[q]) ] = 
	(dists[ min(seqs->sequences[r],   seqs->sequences[q]) ][ max(seqs->sequences[r],   seqs->sequences[q]) ] + 
	 dists[ min(seqs->sequences[r+l], seqs->sequences[q]) ][ max(seqs->sequences[r+l], seqs->sequences[q]) ] -
	 dists[ min(seqs->sequences[r], seqs->sequences[r+l]) ][ max(seqs->sequences[r]  , seqs->sequences[r+l]) ]) / 2.0; 


       if (verbose)
	printf("dist[%d][%d] = %5.4f\n", min(u[s] , seqs->sequences[q]),
	       max(u[s], seqs->sequences[q]),
	       dists[ min(u[s], seqs->sequences[q]) ][ max(u[s] , seqs->sequences[q]) ]);
	

      /**
      dists[ min(u[s] , seqs->sequences[q]) ][ max(u[s] , seqs->sequences[q]) ] = 
	(dists[ min(seqs->sequences[r],   seqs->sequences[q]) ][ max(seqs->sequences[r],   seqs->sequences[q]) ] + 
	 dists[ min(seqs->sequences[r+l], seqs->sequences[q]) ][ max(seqs->sequences[r+l], seqs->sequences[q]) ]) / 2.0; 
      **/

	 /** - dists[ seqs->sequences[r] ][ seqs->sequences[r+l] ]) / 2; **/

   
    }
  }
  
  /** distances entre séquences intérieures à la duplication **/
  for (r=i, s=0; r<i+l; r++, s++) {
    for (q=r+1, t=s+1; q<i+l; q++, t++) {
      //printf("modifie la distance entre %d et %d\n", seqs->sequences[r], seqs->sequences[q]);

      if (verbose) {
	printf("dist[%d][%d] = ( dist[%d][%d](%5.4f) + dist[%d][%d](%5.4f) ) / 4.0 + ( dist[%d][%d](%5.4f) + dist[%d][%d](%5.4f) ) / 4.0\n", 
	       min( u[s] , u[t] ), max(u[s] , u[t] ),
	       min(seqs->sequences[r  ], seqs->sequences[q  ]), max(seqs->sequences[r  ], seqs->sequences[q  ]),
	       dists[ min(seqs->sequences[r  ], seqs->sequences[q  ]) ][ max(seqs->sequences[r  ], seqs->sequences[q  ]) ],
	       min(seqs->sequences[r+l], seqs->sequences[q  ]), max(seqs->sequences[r+l], seqs->sequences[q  ]),
	       dists[ min(seqs->sequences[r+l], seqs->sequences[q  ]) ][ max(seqs->sequences[r+l], seqs->sequences[q  ]) ],
	       min(seqs->sequences[r  ], seqs->sequences[q+l  ]), max(seqs->sequences[r  ], seqs->sequences[q+l  ]),
	       dists[ min(seqs->sequences[r  ], seqs->sequences[q+l  ]) ][ max(seqs->sequences[r  ], seqs->sequences[q+l  ]) ],
	       min(seqs->sequences[r+l  ], seqs->sequences[q+l]), max(seqs->sequences[r+l  ], seqs->sequences[q+l]),
	       dists[ min(seqs->sequences[r+l  ], seqs->sequences[q+l]) ][ max(seqs->sequences[r+l  ], seqs->sequences[q+l]) ]
	       );
	
      }
      
      if (!correction)
	dists[ min( u[s] , u[t] ) ][ max(u[s] , u[t] ) ] = 
	  (dists[ min(seqs->sequences[r  ], seqs->sequences[q  ]) ][ max(seqs->sequences[r  ], seqs->sequences[q  ]) ] + 
	 dists[ min(seqs->sequences[r+l], seqs->sequences[q  ]) ][ max(seqs->sequences[r+l], seqs->sequences[q  ]) ] ) / 4.0 +
	  
	  (dists[ min(seqs->sequences[r  ], seqs->sequences[q+l  ]) ][ max(seqs->sequences[r  ], seqs->sequences[q+l  ]) ] + 
	   dists[ min(seqs->sequences[r+l  ], seqs->sequences[q+l]) ][ max(seqs->sequences[r+l  ], seqs->sequences[q+l]) ] ) / 4.0;
      
      else
      
      
      dists[ min( u[s] , u[t] ) ][ max(u[s] , u[t] ) ] = 
	(dists[ min(seqs->sequences[r  ], seqs->sequences[q  ]) ][ max(seqs->sequences[r  ], seqs->sequences[q  ]) ] + 
	 dists[ min(seqs->sequences[r+l], seqs->sequences[q  ]) ][ max(seqs->sequences[r+l], seqs->sequences[q  ]) ] - 
	 dists[ min(seqs->sequences[r]  , seqs->sequences[r+l]) ][ max(seqs->sequences[r  ], seqs->sequences[r+l]) ]) / 4.0 +
 
	(dists[ min(seqs->sequences[r  ], seqs->sequences[q+l  ]) ][ max(seqs->sequences[r  ], seqs->sequences[q+l  ]) ] + 
	 dists[ min(seqs->sequences[r+l  ], seqs->sequences[q+l]) ][ max(seqs->sequences[r+l  ], seqs->sequences[q+l]) ] - 
	 dists[ min(seqs->sequences[q  ], seqs->sequences[q+l]) ][ min(seqs->sequences[q  ], seqs->sequences[q+l]) ]) / 4.0;
      

      /** -
	  dists[ min(seqs->sequences[r]  , seqs->sequences[r+l]) ][ max(seqs->sequences[r  ], seqs->sequences[r+l]) ] **/
      /** -
	  dists[ min(seqs->sequences[q  ], seqs->sequences[q+l]) ][ min(seqs->sequences[q  ], seqs->sequences[q+l]) ] **/
    }
  }
}






/** charge une matrice de distances **/
float** loadDistances(char* file, int* nb_spec) {
  float** dists; 
  int lig, col;

  int matsize;
  
  FILE* fp = fopen(file, "r");
  if (fp == 0) { printf("impossible d'ouvrir %s\n", file); exit(0); }
  
  fscanf(fp, " %d", nb_spec);
  //printf("Matrice de %d distances...\n", *nb_spec);

  matsize = 2 * (*nb_spec) - 1;

  names = (char**)malloc(*nb_spec*sizeof(char*));
  for (lig=0; lig<*nb_spec; lig++)
    names[lig] = (char*)calloc(10, sizeof(char));
  
  dists = (float**)malloc(matsize*sizeof(float*));

  for(lig=0; lig < matsize; lig++) {

    dists[lig] = (float*)malloc(matsize*sizeof(float));
    
  }
  
  for(lig=0; lig < *nb_spec; lig++) {
    fscanf(fp,"%s",names[lig]);  

    //dists[lig] = (float*)malloc(matsize*sizeof(float));
    
    for(col= 0; col < *nb_spec; col++) {
      fscanf(fp,"%f", &(dists[lig + *nb_spec - 1][col + *nb_spec - 1]));             /* read the distance  */
    }
  }
  
  return dists;
}


/**
 *  affiche la matrice de couts
 *  
 */
void printCostMatrix(FILE* fp, float** dists, SEQS* seqs) 
{
  int i, j;
  short* s = seqs->sequences;

  for (i=0; i<seqs->size; i++) {
    for (j=0; j<seqs->size; j++) {
      if (j>i)
	fprintf(fp, "%3d ", (int)dists[ min(s[i],s[j]) ][ max(s[i],s[j]) ]);  
      else
	fprintf(fp, "%3d ", (int)0.000);  
    }
    fprintf(fp, "\n");
  }
  fprintf(fp, "\n");
}


/**
 *  affiche la matrice de couts
 *  
 */
void showCostMatrix(float** dists, SEQS* seqs) 
{
  int i, j;
  short* s = seqs->sequences;

  for (i=0; i<seqs->size; i++) {
    for (j=0; j<seqs->size; j++) {
      printf("%3d ", (int)dists[ min(s[i],s[j]) ][ max(s[i],s[j]) ]);    
    }
    printf("\n");
  }

}

void printDistances(FILE* fp, float** dists, int nb_spec) 
{
	int lig, col;
	for(lig=0; lig < nb_spec; lig++) {
		for(col= 0; col < nb_spec; col++) {
			fprintf(fp, "%3.4f ", dists[lig][col]);
		}	
		fprintf(fp, "\n");
	}
	fprintf(fp, "\n");
}


/** affiche une matrice de distances **/
void showDistances(float** dists, int nb_spec) {
  int lig, col;
  for(lig=0; lig < nb_spec; lig++) {
    for(col= 0; col < nb_spec; col++) {
      printf("%3.4f ", dists[lig][col]);     
    }
    printf("\n");
  }
}

/** alloue un nouvel arbre 
 *
 *  @param T int nombre de feuilles
 */
tree* allocTree(int T) {
  
  int i;
  
  int maxnodes = 2*T-1;
  
  tree* t     = (tree*)malloc(sizeof(tree));
  
  t->nodes    = (int*) malloc(maxnodes * sizeof(int));
  t->children = (int**)malloc(maxnodes * sizeof(int*));
  
  t->maxnodes    = maxnodes;


  for (i=0; i<maxnodes; i++)
    t->children[i] = (int*)malloc(2*sizeof(int));
  
  for (i=0; i<t->maxnodes; i++)
    t->nodes[i] = -1;

  for (i=0; i<t->maxnodes; i++) {
    t->children[i][0] = -1;
    t->children[i][1] = -1;
  }

  t->maxleaves   = T;
  t->maxintnodes = T-1;
  
  t->numleaves   = T;
  t->numintnodes = 1;

  //t->nodes[t->maxintnodes]   = 0;
  //t->nodes[t->maxintnodes+1] = 0;
  //t->children[0][0]          = t->maxintnodes;
  //t->children[0][1]          = t->maxintnodes + 1;

  
 
  return t;
}




/*
 * retourne le cout d'une paire p, p+l
 *
 */
float getPairCost(
		  float** mycrdata,
		  float** mydata,
		  int     p, 
		  int     l,
		  int     lastopti, 
		  int     lastoptl, 
		  SEQS*   seqs,
		  short*  u, 
		  short*  v,
		  short*  w,
		  short*  uv) {
  
  int     i, j, t;
  float   s1, s2, s3;
  float   cost;
  short*  s;

  nbpairs++;

  s    = seqs->sequences;
  cost = mycrdata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ];
  

  /** si la paire p,p+l est dans le cas A, cad ne chevauche pas l'ancienne duplication et qu'elle n'est pas à -1 **/
  if ((cost > -0.5) && 
      (!((p+l >= lastopti) && (p+l < lastopti+lastoptl)) || ((p >= lastopti) && (p < lastopti+lastoptl)))) {
    
    
    if (verbose)
      printf("Cout initial de la paire (%d, %d) = %3.1f\n", s[p], s[p+l], cost);
    
    nbdblesums++;
    
    /** calcul du cout de la paire p, p+l **/
    
    
    for (i=0; i<2*lastoptl-1; i++) {
      for (j=i+1; j<2*lastoptl; j++) {
	if (verbose)
	  printf("étudie le quadruplet (%d, %d) (%d, %d)", s[p], s[p+l], uv[i], uv[j]);
	
	s1 = mydata[ min(s[p], s[p+l])  ][ max(s[p], s[p+l]) ]  + mydata[ min(uv[i], uv[j])  ][ max(uv[i], uv[j])   ];
	s2 = mydata[ min(s[p], uv[j]  ) ][ max(s[p], uv[j]  ) ] + mydata[ min(uv[i], s[p+l]) ][ max(uv[i], s[p+l]) ];
	s3 = mydata[ min(s[p], uv[i]  ) ][ max(s[p], uv[i]  ) ] + mydata[ min(uv[j], s[p+l]) ][ max(uv[j], s[p+l]) ];

	//nbdblesums++;
	
	if ((s1 < s2) && (s1 < s3)) {
	  cost -= 1.0;
	  if (verbose) {
		  printf(" -1\n");
	  }
	} else {
	  if (verbose) {
	    printf("   \n");
	  }
	}
      }
    }
    
    for (i=0; i<lastoptl; i++) {
      for (t=0; t<seqs->size; t++) {
	if ((t == p) || (t == p+l) || (u[i] == s[t]))
	  continue;
	
	/** cout des paire (p,p+l)(u_i, t) **/
	if ((t < lastopti) || (t >= lastopti+lastoptl)) {  
	  
	  if (verbose)
	    printf("etudie le quadruplet (%d, %d) (%d, %d) (u)", min(s[p], s[p+l]), max(s[p], s[p+l]), min(u[i], s[t]), max(u[i], s[t]));
	  
	  s1 = mydata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ] + mydata[ min(u[i], s[t])   ][ max(u[i], s[t])   ];
	  s2 = mydata[ min(s[p], s[t]  ) ][ max(s[p], s[t]  ) ] + mydata[ min(u[i], s[p+l]) ][ max(u[i], s[p+l]) ];
	  s3 = mydata[ min(s[p], u[i]  ) ][ max(s[p], u[i]  ) ] + mydata[ min(s[t], s[p+l]) ][ max(s[t], s[p+l]) ];

	  //nbdblesums++;

	  if ((s1 < s2) && (s1 < s3)) {
	    cost -= 1.0;
	    if (verbose) {
	      printf(" -1\n");
	    }
	  } else {
	    if (verbose) {
	      printf("   \n");
	    }
	  }
	}
	
	
	
	if ((t < lastopti) || (t >= lastopti+lastoptl)) {
	  
	  if (verbose)
	    printf("etudie le quadruplet (%d, %d) (%d, %d) (v)", min(s[p], s[p+l]), max(s[p], s[p+l]), min(v[i], s[t]), max(v[i], s[t]));
	  
	  /** cout des paire (p,p+l)(v_i, t) **/
	  s1 = mydata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ] + mydata[ min(v[i], s[t])   ][ max(v[i], s[t])   ];
	  s2 = mydata[ min(s[p], s[t]  ) ][ max(s[p], s[t]  ) ] + mydata[ min(v[i], s[p+l]) ][ max(v[i], s[p+l]) ];
	  s3 = mydata[ min(s[p], v[i]  ) ][ max(s[p], v[i]  ) ] + mydata[ min(s[t], s[p+l]) ][ max(s[t], s[p+l]) ];

	  //nbdblesums++;
	  
	  if ((s1 < s2) && (s1 < s3)) {
	    cost -= 1.0;
	    if (verbose) {
	      printf(" -1\n");
	    }
	    
	  } else {
	    if (verbose) {
	      printf("   \n");
	    }
	  }
	}
	
	if ((t < lastopti) || (t >= lastopti+lastoptl)) {
	  if (verbose)
	    printf("etudie le quadruplet (%d, %d) (%d, %d) (w)", min(s[p], s[p+l]), max(s[p], s[p+l]), min(w[i], s[t]), max(w[i], s[t]));
	  
	  /** cout des paire (p,p+l)(w_i, t), valable quelque soit t **/
	  s1 = mydata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ] + mydata[ min(w[i], s[t])   ][ max(w[i], s[t])   ];
	  s2 = mydata[ min(s[p], s[t]  ) ][ max(s[p], s[t]  ) ] + mydata[ min(w[i], s[p+l]) ][ max(w[i], s[p+l]) ];
	  s3 = mydata[ min(s[p], w[i]  ) ][ max(s[p], w[i]  ) ] + mydata[ min(s[t], s[p+l]) ][ max(s[t], s[p+l]) ];
	  
	  //nbdblesums++;

	  if ((s1 < s2) && (s1 < s3)) {
	    cost += 1.0;
	    if (verbose) {
	      printf(" +1\n");
	    }
	  } else {
	    if (verbose) {
	      printf("   \n");
	    }
	  }
	  
	}
	
	
      }


      for (j=i+1; j<lastoptl; j++) {
	if (verbose)
	  printf("etudie le quadruplet (%d, %d) (%d, %d) (w)", min(s[p], s[p+l]), max(s[p], s[p+l]), min(w[i], w[j]), max(w[i], w[j]));
	/** cout des paire (p,p+l)(w_i, t), valable quelque soit t **/
	s1 = mydata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ] + mydata[ min(w[i], w[j])   ][ max(w[i], w[j])   ];
	s2 = mydata[ min(s[p], w[j]  ) ][ max(s[p], w[j]  ) ] + mydata[ min(w[i], s[p+l]) ][ max(w[i], s[p+l]) ];
	s3 = mydata[ min(s[p], w[i]  ) ][ max(s[p], w[i]  ) ] + mydata[ min(w[j], s[p+l]) ][ max(w[j], s[p+l]) ];

	//nbdblesums++;
	
	if ((s1 < s2) && (s1 < s3)) {
	  cost += 1.0;
	  if (verbose) {
	    printf(" +1\n");
	  }
	} else {
	  if (verbose) {
	    printf("   \n");
	  }
	}
      }
    }
    
  } else {

    nbpairscalc++;	  
    /** CAS B et C : au moins un taxon parmi le quadruplet appartient a w, donc on recalcule completement la paire **/
    if (verbose)
      printf("Calcul de la paire (%d, %d) from scratch... \n", s[p], s[p+l]);
    cost = 0.0;
    for (i=0; i<seqs->size; i++) {
      
      if ((i == p) || (i == p+l))
	continue;
      for (t=i+1; t<seqs->size; t++) {
	if ((t == p) || (t == p+l))
	  continue;
	
	if (verbose)
	  printf("etudie le quadruplet (%d, %d) (%d, %d) ", min(s[p], s[p+l]), max(s[p], s[p+l]), min(s[i], s[t]), max(s[i], s[t]));
	
	
	
	/** cout des paire (p,p+l)(w_i, t), valable quelque soit t **/
	s1 = mydata[ min(s[p], s[p+l]) ][ max(s[p], s[p+l]) ] + mydata[ min(s[i], s[t])   ][ max(s[i], s[t])   ];
	s2 = mydata[ min(s[p], s[t]  ) ][ max(s[p], s[t]  ) ] + mydata[ min(s[i], s[p+l]) ][ max(s[i], s[p+l]) ];
	s3 = mydata[ min(s[p], s[i]  ) ][ max(s[p], s[i]  ) ] + mydata[ min(s[t], s[p+l]) ][ max(s[t], s[p+l]) ];
	
	//nbdblesums++;

	if ((s1 < s2) && (s1 < s3)) {
	  cost += 1.0;
	  if (verbose)
	    printf("+1\n");
	  
	} else 
	  if (verbose)
	    printf("\n");
	
      }
    }
    
    
  }
  
  
  
  return cost;
}


char* get_parameter(int argc, char** argv, char* param)
{
  int i = 0;
  while ((i < argc) && (strcmp(param, argv[i]))) 
    i++;
  if (i<argc)
    return (argv[i+1]);
  else
    return "";
}