4 * (c) 2003 Gangolf Jobb (http://www.treefinder.de)
6 * Various data structures and methods for manipulating
8 * (e.g., to classify genes)
10 * This code may be distributed under the GNU GPL
17 /*============================================================================*/
23 struct tree *branches,*next;
33 b=malloc(sizeof(Tree)); if(!b) return(NULL);
35 *b->label='\0'; b->branches=NULL; b->next=NULL; b->length=0.0;
42 void FreeTree(Tree *t) {
45 while(t->branches) {b=t->branches; t->branches=b->next; FreeTree(b);}
51 /*============================================================================*/
53 static Tree *current=NULL;
56 static int edge_index;
57 static int node_index;
59 enum {OK=0,ERROR}; /* one may invent more descriptive error codes */
62 /*============================================================================*/
64 Tree *buildTreeFromPhylo_(
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;}
85 buildTreeFromPhylo_(j,lowerNodes,upperNodes,edgeLengths,nedges,tipLabels,ntips);
87 b->length=edgeLengths[i];
88 *bb=b; bb=&b->next; n++;
89 } if(n<2) {error=ERROR; goto err;}
98 void buildTreeFromPhylo(
110 if(current) {FreeTree(current); current=NULL;}
112 if(*nedges<2||*ntips<2) {error=ERROR; *result=error; return;}
114 current=buildTreeFromPhylo_(-1,lowerNodes,upperNodes,edgeLengths,*nedges,tipLabels,*ntips);
116 if(error&¤t) {FreeTree(current); current=NULL;}
122 /*============================================================================*/
124 void destroyTree(int *result) {
128 if(current) {FreeTree(current); current=NULL;}
134 /*============================================================================*/
136 void getError(int *result) {
144 /*============================================================================*/
146 int nTips_(Tree *t) {
150 if(!t->branches) return(1);
152 n=0; for(b=t->branches;b;b=b->next) n+=nTips_(b);
159 void nTips(int *result) {
163 if(!current) {error=ERROR; *result=0; return;}
165 *result=nTips_(current);
169 /*============================================================================*/
171 int nNodes_(Tree *t) {
175 if(!t->branches) return(1);
177 n=1; for(b=t->branches;b;b=b->next) n+=nNodes_(b);
184 void nNodes(int *result) {
188 if(!current) {error=ERROR; *result=0; return;}
190 *result=nNodes_(current);
196 void nEdges(int *result) {
200 if(!current) {error=ERROR; *result=0; return;}
202 *result=nNodes_(current)-1;
205 /*============================================================================*/
207 double markClasses_(Tree *t) {
208 Tree *b; double destinct,sum;
210 /* all tips above a marked ( == 1) node belong to the same class */
212 if(!t->branches) {t->mark=1; return(t->length);}
214 sum=0.; for(b=t->branches;b;b=b->next) sum+=markClasses_(b);
216 destinct=nTips_(t)*(t->length); /* (t->length) == 0. at root */
218 if(destinct>sum) { t->mark=1; return(destinct); }
220 t->mark=0; return(sum);
226 void getMisawaTajima__(Tree *t,int ignore,int *result) { /* maps tips to marked classes */
229 if(t->mark&&!ignore) {node_index++; ignore=1;} /* marked nodes above a marked node will be ignored */
231 if(!t->branches) {result[tip_index++]=node_index; return;}
233 for(b=t->branches;b;b=b->next) getMisawaTajima__(b,ignore,result);
239 void getMisawaTajima_(Tree *t,int *result) {
243 tip_index=0; node_index=0;
245 getMisawaTajima__(t,0,result);
251 void getMisawaTajima(int *result) {
255 if(!current) {error=ERROR; return;}
257 getMisawaTajima_(current,result);