various bug fixes since the release of ape 3.0

[ape.git] / src / me.c

diff --git a/src/me.c b/src/me.c
index 69d1cb210bb5d17c3a69c928870457e4032e9098..6841802ae9a62b222af8e6faeccde35a12600f2a 100644 (file)
--- a/src/me.c
+++ b/src/me.c
@@ -1,4 +1,4 @@
-/* me.c    2008-01-14 */
+/* me.c    2012-02-09 */
 
 /* Copyright 2007-2008 Olivier Gascuel, Rick Desper,
    R port by Vincent Lefort, me_*() below modified by Emmanuel Paradis */
@@ -20,9 +20,14 @@ void makeOLSAveragesTable(tree *T, double **D, double **A);
 void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies);
 //functions from NNI.c
 void NNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies);
+//functions from SPR.c
+void SPR(tree *T, double **D, double **A, int *count);
+//functions from TBR.c
+void TBR(tree *T, double **D, double **A);
 
 
-void me_b(double *X, int *N, char **labels, int *nni,
+void me_b(double *X, int *N, char **labels,
+         int *nni, int *spr, int *tbr,
          int *edge1, int *edge2, double *el, char **tl)
 {
   double **D, **A;
@@ -46,10 +51,12 @@ void me_b(double *X, int *N, char **labels, int *nni,
     T = BMEaddSpecies(T, addNode, D, A);
   }
   // Compute bNNI
-  if (*nni == 1)
-    bNNI(T, A, &nniCount, D, n);
+  if (*nni) bNNI(T, A, &nniCount, D, n);
   assignBMEWeights(T,A);
 
+  if (*spr) SPR(T, D, A, &nniCount);
+  if (*tbr) TBR(T, D, A);
+
   tree2phylo(T, edge1, edge2, el, tl, n);
 
   freeMatrix(D,n);
@@ -85,7 +92,7 @@ void me_o(double *X, int *N, char **labels, int *nni,
   }
   makeOLSAveragesTable(T,D,A);
   // Compute NNI
-  if (*nni == 1)
+  if (*nni)
     NNI(T,A,&nniCount,D,n);
   assignOLSWeights(T,A);
 
@@ -267,8 +274,7 @@ tree *detrifurcate(tree *T)
     return(T);
   if (NULL != v->parentEdge)
     {
-      Rprintf ("Error: root %s is poorly rooted.\n",v->label);
-      exit(0);
+      error("root %s is poorly rooted.", v->label);
     }
   for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
     {
@@ -319,8 +325,7 @@ void compareSets(tree *T, set *S)
     }
   if (-1 == v->index2)
     {
-      Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
-      exit(0);
+      error("leaf %s in tree not in distance matrix.", v->label);
     }
   e = depthFirstTraverse(T,NULL);
   while (NULL != e)
@@ -328,16 +333,14 @@ void compareSets(tree *T, set *S)
       v = e->head;
       if ((leaf(v)) && (-1 == v->index2))
        {
-         Rprintf("Error leaf %s in tree not in distance matrix.\n",v->label);
-         exit(0);
+         error("leaf %s in tree not in distance matrix.", v->label);
        }
       e = depthFirstTraverse(T,e);
       }
   for(X = S; NULL != X; X = X->secondNode)
     if (X->firstNode->index2 > -1)
       {
-       Rprintf("Error node %s in matrix but not a leaf in tree.\n",X->firstNode->label);
-       exit(0);
+       error("node %s in matrix but not a leaf in tree.", X->firstNode->label);
       }
   return;
 }