1 /* triangMtd.c 2012-04-02 */
3 /* Copyright 2011-2012 Andrei-Alin Popescu */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
9 * leafs labelled 1 to n. root labelled n+1. other nodes labelled n+1 to m
14 int give_indexx(int i, int j, int n)
15 { if (i == j) return 0;
16 if (i > j) return(n*(j - 1) - j*(j - 1)/2 + i - j - 1);
17 else return(n*(i - 1) - i*(i - 1)/2 + j - i - 1);
20 int pred(int k, int* ed1, int* ed2, int numEdges)
21 /* find the predecesor of vertex k */
24 for (i = 0; i <= numEdges; i++)
25 if (ed2[i] == k) return ed1[i];
29 int* getPathBetween(int x, int y, int n, int* ed1, int* ed2, int numEdges)
30 //get the path between vertices x and y in an array ord.
31 //ord[i]=j means that we go between i and j on the path between x and y
34 int ch[2*n-1];//ch[i]==1 implies {k,pred(k)} on path between x and y
42 k=pred(k,ed1,ed2,numEdges);
49 k=pred(k,ed1,ed2,numEdges);
51 int *ord=(int*)malloc(sizeof(int)*(2*n-1));
52 //starting from x, fill ord
60 p=pred(p,ed1,ed2,numEdges);
68 p=pred(p,ed1,ed2,numEdges);
69 ord[p]=aux;//other way
75 int getLength(int x, int y, int* ed1, int* ed2, int numEdges, int* edLen)
76 /* get length of edge {x,y}, -1 if edge does not exist */
79 for (i = 0; i <= numEdges; i++)
80 if ((ed1[i] == x && ed2[i] == y) || (ed1[i] == y && ed2[i] == x))
85 void triangMtd(double* d, int* np, int* ed1,int* ed2, double* edLen)
96 //error("%f ,giveindex= %f",d[ij],d[give_indexx(i,j,n)]);
112 {error("%f ",d[i*n+j]);
118 int x=-1,y=-1,numEdges=0;
119 int st[n+1];//array holding at position i the starting point of attachment path of leaf i
120 int ed[n+1];//array holding at position i the ending point of attachment path of leaf i
121 double l[n+1];//array holding at position i the distance of leaf i from the partially constructed tree
122 int w[n+1];//presence array:w[i]=1 if leaf i is added to the tree, 0 otherwise
123 int wSize=0;//number of leaves added so far
125 for(i=1;i<=n;i++){w[i]=0;}
126 //find the minimum distance in our matrix
129 //choose two leaves x,y such that d[x][y] is minimal
133 {if(d[give_indexx(i,j,n)]<=0)continue;
134 minW=d[give_indexx(i,j,n)];
145 if(d[give_indexx(i,j,n)]<minW)
147 minW=d[give_indexx(i,j,n)];
153 w[x]=1; w[y]=1;//mark x and y as added
155 //calculate the distance between leaves not in w and leaves in w
160 l[i]=0.5*(d[give_indexx(i,x,n)]+d[give_indexx(i,y,n)]-d[give_indexx(x,y,n)]);//distance of leaf i from the
165 wSize+=3;//since we construct a star tree on three leaves
166 int nv=n+1;//since first n numbers are reserved for leaves
169 //search for x3 that is closest to the tree
171 {if(w[i])continue;//skip if leaf already added
178 //construct initial star tree on three leaves:x,y,x3. nv is the interior
179 //vertex and also the root of the tree
182 edLen[numEdges]=0.5*(d[give_indexx(x3,x,n)]+d[give_indexx(y,x,n)]-d[give_indexx(y,x3,n)]);
186 edLen[numEdges]=0.5*(d[give_indexx(y,x3,n)]+d[give_indexx(y,x,n)]-d[give_indexx(x,x3,n)]);
190 edLen[numEdges]=0.5*(d[give_indexx(x3,x,n)]+d[give_indexx(y,x3,n)]-d[give_indexx(y,x,n)]);
193 /*Rprintf("new tree\n");
196 Rprintf("%i->%i of length:%f \n",ed1[i],ed2[i],edLen[i]);
198 Rprintf("end new tree\n");*/
200 //calculate distance of leaves not yet added to the star tree
207 double newL=0.5*(d[give_indexx(i,s,n)]+d[give_indexx(s,x3,n)]-d[give_indexx(i,x3,n)]);
217 //partial tree construction begins here
220 { double minDist=1e50;
223 //search for leaf z which is closest to partial tree
225 {if(w[i])continue;//skip if leaf already added
232 //x and y are the leaves on the path between which z is attached
236 int* ord=getPathBetween(x,y,n,ed1,ed2,numEdges);
238 for(i=1;i<=2*n-2;i++)
239 {Rprintf("ord[%i]=%i ",i,ord[i]);
242 //order the path from x to y, in an array ord, ord[i]=j means vertex i comes before j
243 //first count number of edges on path (i.e count all i s.t ch[i]==1)
246 //look for the edge on the path x to y to subdivide
251 int subdiv=-1;//index of edge to subdivide
252 //error("d[y,x]=%f,d[z,x]=%f,d[z,y]=%f\n",d[give_indexx(y,x,n)],d[give_indexx(z,x,n)],d[give_indexx(z,y,n)]);
253 double lx=0.5*(d[give_indexx(y,x,n)]+d[give_indexx(z,x,n)]-d[give_indexx(z,y,n)]);//distance of attachment
254 // Rprintf("adding %i on the path between %i and %i at a distance from x of %f and a distance of %f from tree",z,x,y,lx,minDist);
256 //Rprintf("Adding leaf %i, between %i and %i\n",z,x,y);
257 // Rprintf("%i situated at a distance of %d from tree",z,minDist);
259 //Rprintf("path between %i and %i\n",x,y);
261 while(p!=y && sum<lx)
264 //error("%i -> %i ",p,ord[p]);
267 for(i=0;i<=numEdges;i++)
269 if((ed1[i]==aux && ed2[i]==p)||(ed2[i]==aux && ed1[i]==p))
271 if(ed1[i]==aux && ed2[i]==p){sw=1;}
279 //subdivide subdiv with a node labelled nv
281 //multifurcating vertices
285 edLen[subdiv]= (sw==1)?(lx-prevSum):(sum-lx);//check which 'half' of the
286 //path the leaf belongs to
287 //and updates accordingly
288 //error("sum=%f, prevsum=%f\n",sum,prevSum);
289 //error("lx-prevSum=%f, sum-lx=%f, minDist=%f",lx-prevSum,sum-lx,minDist);
293 edLen[numEdges]= (sw==1)?(sum-lx):(lx-prevSum);
295 edLen[numEdges]=minDist;
302 /*update distance matrix, only needed for incomplete distances
306 error("path between %i and %i\n",ii,z);
307 int* ord=getPathBetween(ii,z,n,ed1,ed2,numEdges);
313 {error("ord[%i]=%i ",i,ord[i]);
318 { //error("p=%i\n",p);
321 newDist+=getLength(ii,z,ed1,ed2,numEdges,edLen);
326 d[ii][z]=d[z][ii]=newDist;
328 //update l[s] for all s not yet added
334 if(i!=x && i!=y)continue;//we only consider x and y as being other leaf
335 //and take the minimum of them as being new distance
336 double newL=0.5*(d[give_indexx(i,s,n)]+d[give_indexx(z,s,n)]-d[give_indexx(i,z,n)]);//one of leaves is
338 //all pairs not cotaining z
339 //will remain unchanged
349 /*Rprintf("new tree\n");
352 Rprintf("%i->%i of length:%f \n",ed1[i],ed2[i],edLen[i]);
354 Rprintf("end new tree\n");*/
356 //for(i=0;i<=numEdges;i++){ed1[i]++;ed2[i]++;}