]> git.donarmstrong.com Git - ape.git/blobdiff - src/me.c
modified ace( method = ...)
[ape.git] / src / me.c
index 5b383f33b2c685f0fd8887ada72514bb94e2b09c..8a00591ab680cb6846fc9ea7765c4a9b7e893c0f 100644 (file)
--- a/src/me.c
+++ b/src/me.c
@@ -1,26 +1,39 @@
+/* me.c    2012-04-30 */
+
+/* Copyright 2007-2008 Olivier Gascuel, Rick Desper,
+   R port by Vincent Lefort and Emmanuel Paradis */
+
+/* This file is part of the R-package `ape'. */
+/* See the file ../COPYING for licensing issues. */
+
 #include "me.h"
 
 //functions from me_balanced.c
-tree *BMEaddSpecies(tree *T,node *v, double **D, double **A);
+tree *BMEaddSpecies(tree *T, node *v, double **D, double **A);
 void assignBMEWeights(tree *T, double **A);
 void makeBMEAveragesTable(tree *T, double **D, double **A);
 //functions from me_ols.c
-tree *GMEaddSpecies(tree *T,node *v, double **D, double **A);
+tree *GMEaddSpecies(tree *T, node *v, double **D, double **A);
 void assignOLSWeights(tree *T, double **A);
 void makeOLSAveragesTable(tree *T, double **D, double **A);
 //functions from bNNI.c
 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, char **treeStr, int *nni)
+void me_b(double *X, int *N, int *labels,
+         int *nni, int *spr, int *tbr,
+         int *edge1, int *edge2, double *el)
 {
   double **D, **A;
   set *species, *slooper;
   node *addNode;
   tree *T;
-  char *str;
   int n, nniCount;
 
   n = *N;
@@ -29,53 +42,37 @@ void me_b(double *X, int *N, char **labels, char **treeStr, int *nni)
   species = (set *) malloc(sizeof(set));
   species->firstNode = NULL;
   species->secondNode = NULL;
-  str = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
-  /* added by EP */
-  if (strlen(str))
-    strncpy(str, "", strlen(str));
-  /* end */
-  D = loadMatrix (X, labels, n, species);
-  A = initDoubleMatrix(2 * n - 2);
+  D = loadMatrix(X, labels, n, species);
+  A = initDoubleMatrix(2*n - 2);
 
   for(slooper = species; NULL != slooper; slooper = slooper->secondNode)
   {
     addNode = copyNode(slooper->firstNode);
-    T = BMEaddSpecies(T,addNode,D,A);
+    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);
 
-  NewickPrintTreeStr(T,str);
+  if (*spr) SPR(T, D, A, &nniCount);
+  if (*tbr) TBR(T, D, A);
 
-  if (strlen (str) < MAX_INPUT_SIZE -1)
-    {
-      *treeStr = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
-      /* added by EP */
-      if (strlen(*treeStr))
-       strncpy(*treeStr, "", strlen(*treeStr));
-      /* end */
-      strncpy (*treeStr, str, strlen(str));
-    }
+  tree2phylo(T, edge1, edge2, el, labels, n);
 
-/*   free (str); */
   freeMatrix(D,n);
   freeMatrix(A,2*n - 2);
   freeSet(species);
   freeTree(T);
   T = NULL;
-
-  /* return; */
 }
 
-void me_o(double *X, int *N, char **labels, char **treeStr, int *nni)
+void me_o(double *X, int *N, int *labels, int *nni,
+         int *edge1, int *edge2, double *el)
 {
   double **D, **A;
   set *species, *slooper;
   node *addNode;
   tree *T;
-  char *str;
   int n, nniCount;
 
   n = *N;
@@ -84,11 +81,6 @@ void me_o(double *X, int *N, char **labels, char **treeStr, int *nni)
   species = (set *) malloc(sizeof(set));
   species->firstNode = NULL;
   species->secondNode = NULL;
-  str = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
-  /* added by EP */
-  if (strlen(str))
-    strncpy(str, "", strlen(str));
-  /* end */
 
   D = loadMatrix (X, labels, n, species);
   A = initDoubleMatrix(2 * n - 2);
@@ -100,37 +92,24 @@ void me_o(double *X, int *N, char **labels, char **treeStr, int *nni)
   }
   makeOLSAveragesTable(T,D,A);
   // Compute NNI
-  if (*nni == 1)
+  if (*nni)
     NNI(T,A,&nniCount,D,n);
   assignOLSWeights(T,A);
 
-  NewickPrintTreeStr(T,str);
-
-  if (strlen (str) < MAX_INPUT_SIZE -1)
-    {
-      *treeStr = (char *)R_alloc(MAX_INPUT_SIZE, sizeof(char));
-      /* added by EP */
-      if (strlen(*treeStr))
-       strncpy(*treeStr, "", strlen(*treeStr));
-      /* end */
-      strncpy (*treeStr, str, strlen (str));
-    }
+  tree2phylo(T, edge1, edge2, el, labels, n);
 
- /*  free (str); */
   freeMatrix(D,n);
   freeMatrix(A,2*n - 2);
   freeSet(species);
   freeTree(T);
   T = NULL;
-
-  return;
 }
 
-/*************************************************************************
+/*
 
-                           MATRIX FUNCTIONS
+  -- MATRIX FUNCTIONS --
 
-*************************************************************************/
+*/
 
 double **initDoubleMatrix(int d)
 {
@@ -146,9 +125,10 @@ double **initDoubleMatrix(int d)
   return(A);
 }
 
-double **loadMatrix (double *X, char **labels, int n, set *S)
+//double **loadMatrix (double *X, char **labels, int n, set *S)
+double **loadMatrix (double *X, int *labels, int n, set *S)
 {
-  char nextString[MAX_LABEL_LENGTH];
+//  char nextString[MAX_LABEL_LENGTH];
   node *v;
   double **table;
   int i, j, a, b;
@@ -159,9 +139,10 @@ double **loadMatrix (double *X, char **labels, int n, set *S)
 
   for(i=0; i<n; i++)
     {
-      strncpy (nextString, labels[i], MAX_LABEL_LENGTH);
+//      strncpy (nextString, labels[i], MAX_LABEL_LENGTH);
 //      ReplaceForbiddenChars (nextString, '_');
-      v = makeNewNode(nextString,-1);
+//      v = makeNewNode(nextString,-1);
+      v = makeNewNode(labels[i], -1);
       v->index2 = i;
       S = addToSet(v,S);
       for (j=i; j<n; j++) {
@@ -176,11 +157,11 @@ double **loadMatrix (double *X, char **labels, int n, set *S)
   return (table);
 }
 
-/*************************************************************************
+/*
 
-                           GRAPH FUNCTIONS
+  -- GRAPH FUNCTIONS --
 
-*************************************************************************/
+*/
 
 set *addToSet(node *v, set *X)
 {
@@ -197,16 +178,19 @@ set *addToSet(node *v, set *X)
   return(X);
 }
 
-node *makeNewNode(char *label, int i)
+//node *makeNewNode(char *label, int i)
+node *makeNewNode(int label, int i)
 {
   return(makeNode(label,NULL,i));
 }
 
-node *makeNode(char *label, edge *parentEdge, int index)
+//node *makeNode(char *label, edge *parentEdge, int index)
+node *makeNode(int label, edge *parentEdge, int index)
 {
   node *newNode;  /*points to new node added to the graph*/
   newNode = (node *) malloc(sizeof(node));
-  strncpy(newNode->label,label,NODE_LABEL_LENGTH);
+//  strncpy(newNode->label,label,NODE_LABEL_LENGTH);
+  newNode->label = label;
   newNode->index = index;
   newNode->index2 = -1;
   newNode->parentEdge = parentEdge;
@@ -295,8 +279,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 %d is poorly rooted.", v->label);
     }
   for(e = v->middleEdge, v->middleEdge = NULL; NULL != e; e = f )
     {
@@ -325,7 +308,8 @@ void compareSets(tree *T, set *S)
       for(X = S; NULL != X; X = X->secondNode)
        {
          w = X->firstNode;
-         if (0 == strcmp(v->label,w->label))
+//       if (0 == strcmp(v->label,w->label))
+         if (v->label == w->label)
            {
              v->index2 = w->index2;
            w->index2 = -1;
@@ -338,7 +322,8 @@ void compareSets(tree *T, set *S)
   for(X = S; NULL != X; X = X->secondNode)
     {
       w = X->firstNode;
-      if (0 == strcmp(v->label,w->label))
+//      if (0 == strcmp(v->label,w->label))
+      if (v->label == w->label)
        {
          v->index2 = w->index2;
          w->index2 = -1;
@@ -347,8 +332,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 %d in tree not in distance matrix.", v->label);
     }
   e = depthFirstTraverse(T,NULL);
   while (NULL != e)
@@ -356,16 +340,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 %d 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 %d in matrix but not a leaf in tree.", X->firstNode->label);
       }
   return;
 }
@@ -516,4 +498,3 @@ void freeSubTree(edge *e)
   e->head = NULL;
   free(e);
 }
-