src/treefunc.c

   1 /*
   2  *  treefunc.c
   3  *
   4  * (c) 2003  Gangolf Jobb (http://www.treefinder.de)
   5  *
   6  *  Various data structures and methods for manipulating
   7  *  and traversing trees
   8  *  (e.g., to classify genes)
   9  *
  10  *  This code may be distributed under the GNU GPL
  11  */
  12
  13 #include <stdlib.h>
  14 #include <string.h>
  15 #include <math.h>
  16
  17 /*============================================================================*/
  18
  19 #define LABLEN 128
  20
  21 typedef struct tree {
  22  char label[LABLEN];
  23  struct tree *branches,*next;
  24  double length;
  25  int mark;
  26 } Tree;
  27
  28 /*....*/
  29
  30 Tree *NewTree(void) {
  31  Tree *b;
  32
  33  b=malloc(sizeof(Tree)); if(!b) return(NULL);
  34
  35  *b->label='\0'; b->branches=NULL; b->next=NULL; b->length=0.0;
  36
  37  return(b);
  38 }
  39
  40 /*....*/
  41
  42 void FreeTree(Tree *t) {
  43  Tree *b;
  44
  45  while(t->branches) {b=t->branches; t->branches=b->next; FreeTree(b);}
  46
  47  free(t);
  48
  49 }
  50
  51 /*============================================================================*/
  52
  53 static Tree *current=NULL;
  54
  55 static int tip_index;
  56 static int edge_index;
  57 static int node_index;
  58
  59 enum {OK=0,ERROR}; /* one may invent more descriptive error codes */
  60 static int error=OK;
  61
  62 /*============================================================================*/
  63
  64 Tree *buildTreeFromPhylo_(
  65  int node,
  66  int *lowerNodes,
  67  int *upperNodes,
  68  double *edgeLengths,
  69  int nedges,
  70  char **tipLabels,
  71  int ntips
  72 ) {
  73  Tree *t,*b,**bb;
  74  int i,j,n;
  75
  76  t=NewTree();
  77
  78  bb=&t->branches; n=0;
  79  for(i=0;i<nedges;i++) { if(lowerNodes[i]!=node) continue;
  80   j=upperNodes[i];                                             if(j==0) {error=ERROR; goto err;}
  81   if(j>0) {                                                    if(j>ntips) {error=ERROR; goto err;}
  82    b=NewTree(); strcpy(b->label,tipLabels[j-1]);
  83   } else {                                                     if(-j>nedges) {error=ERROR; goto err;}
  84    b=
  85     buildTreeFromPhylo_(j,lowerNodes,upperNodes,edgeLengths,nedges,tipLabels,ntips);
  86   }
  87   b->length=edgeLengths[i];
  88   *bb=b; bb=&b->next; n++;
  89  }                                                             if(n<2) {error=ERROR; goto err;}
  90                                                                err:
  91  *bb=NULL;
  92
  93  return(t);
  94 }
  95
  96 /*....*/
  97
  98 void buildTreeFromPhylo(
  99  int *lowerNodes,
 100  int *upperNodes,
 101  double *edgeLengths,
 102  int *nedges,
 103  char **tipLabels,
 104  int *ntips,
 105  int *result
 106 ) {
 107
 108  error=OK;
 109
 110  if(current) {FreeTree(current); current=NULL;}
 111
 112  if(*nedges<2||*ntips<2) {error=ERROR; *result=error; return;}
 113
 114  current=buildTreeFromPhylo_(-1,lowerNodes,upperNodes,edgeLengths,*nedges,tipLabels,*ntips);
 115
 116  if(error&&current) {FreeTree(current); current=NULL;}
 117
 118  *result=error;
 119
 120 }
 121
 122 /*============================================================================*/
 123
 124 void destroyTree(int *result) {
 125
 126  error=OK;
 127
 128  if(current) {FreeTree(current); current=NULL;}
 129
 130  *result=error;
 131
 132 }
 133
 134 /*============================================================================*/
 135
 136 void getError(int *result) {
 137
 138  *result=error;
 139
 140  error=OK;
 141
 142 }
 143
 144 /*============================================================================*/
 145
 146 int nTips_(Tree *t) {
 147  Tree *b;
 148  int n;
 149
 150  if(!t->branches) return(1);
 151
 152  n=0; for(b=t->branches;b;b=b->next) n+=nTips_(b);
 153
 154  return(n);
 155 }
 156
 157 /*....*/
 158
 159 void nTips(int *result) {
 160
 161  error=OK;
 162
 163  if(!current) {error=ERROR; *result=0; return;}
 164
 165  *result=nTips_(current);
 166
 167 }
 168
 169 /*============================================================================*/
 170
 171 int nNodes_(Tree *t) {
 172  Tree *b;
 173  int n;
 174
 175  if(!t->branches) return(1);
 176
 177  n=1; for(b=t->branches;b;b=b->next) n+=nNodes_(b);
 178
 179  return(n);
 180 }
 181
 182 /*....*/
 183
 184 void nNodes(int *result) {
 185
 186  error=OK;
 187
 188  if(!current) {error=ERROR; *result=0; return;}
 189
 190  *result=nNodes_(current);
 191
 192 }
 193
 194 /*....*/
 195
 196 void nEdges(int *result) {
 197
 198  error=OK;
 199
 200  if(!current) {error=ERROR; *result=0; return;}
 201
 202  *result=nNodes_(current)-1;
 203 }
 204
 205 /*============================================================================*/
 206
 207 double markClasses_(Tree *t) {
 208  Tree *b; double destinct,sum;
 209
 210  /* all tips above a marked ( == 1) node belong to the same class */
 211
 212  if(!t->branches) {t->mark=1; return(t->length);}
 213
 214  sum=0.; for(b=t->branches;b;b=b->next) sum+=markClasses_(b);
 215
 216  destinct=nTips_(t)*(t->length); /* (t->length) == 0. at root */
 217
 218  if(destinct>sum) { t->mark=1;  return(destinct); }
 219
 220                     t->mark=0;  return(sum);
 221
 222 }
 223
 224 /*....*/
 225
 226 void getMisawaTajima__(Tree *t,int ignore,int *result) { /* maps tips to marked classes */
 227  Tree *b;
 228
 229  if(t->mark&&!ignore) {node_index++; ignore=1;} /* marked nodes above a marked node will be ignored */
 230
 231  if(!t->branches) {result[tip_index++]=node_index; return;}
 232
 233  for(b=t->branches;b;b=b->next) getMisawaTajima__(b,ignore,result);
 234
 235 }
 236
 237 /*....*/
 238
 239 void getMisawaTajima_(Tree *t,int *result) {
 240
 241  markClasses_(t);
 242
 243  tip_index=0; node_index=0;
 244
 245  getMisawaTajima__(t,0,result);
 246
 247 }
 248
 249 /*....*/
 250
 251 void getMisawaTajima(int *result) {
 252
 253  error=OK;
 254
 255  if(!current) {error=ERROR; return;}
 256
 257  getMisawaTajima_(current,result);
 258
 259 }