1 /* triangMtds.c 2011-10-11 */
3 /* Copyright 2011 Andrei-Alin Popescu */
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
10 void triangMtds(double* d, int* np, int* ed1,int* ed2, double* edLen)
27 if(i==j){m[i]++;continue;}
28 if(d[give_index(i,j,n)]==-1)continue;
47 while(m[ij]==n){ij++;}
49 {if(s[c[i]]==1)continue;
52 if(s[c[j]]==1)continue;
53 if(d[give_index(c[i],c[j],n)]==-1)
67 double* sub_d=(double*)R_alloc(k*(k - 1)/2, sizeof(double));
74 sub_d[ij++]=d[give_index(i,j,n)];
78 triangMtd(sub_d,&k,ed1,ed2,edLen);
99 if(s[i]==0)continue;//take only leaves not in Y
103 if(s[j]==1)continue;//take only leaves already in Y
104 if(d[give_index(i,j,n)]==-1)continue;//igonore if distance unknown
108 int numEdges=2*k-4;//0-based, so subtract 1
109 //Rprintf("numEdge=%i",numEdges);
113 //now find element in X\Y such that it has most known distances to already
114 //built tree until all leaves are added or we can not find a place to attach
116 //s[i]=1 => i not added to tree
128 s[maxPos]=0;//mark maxPos as added
129 //calculate new m values for leaves not added, i.e we just increment any
130 //already present value by 1 if we know the distance between i and maxPos
134 if(d[give_index(i,maxPos,n)]==-1)continue;
138 //find path to attach maxPos to, grow tree
143 {if(s[i]==1 || d[give_index(i,z,n)]==-1 || i==z)continue;
146 if(s[j]==1 || d[give_index(j,z,n)]==-1 || j==z)continue;
147 double tDist=(d[give_index(i,z,n)]+d[give_index(j,z,n)]-d[give_index(i,j,n)])/2;
157 {error("could not build tree from given distance matrix");
159 int* ord=getPathBetween(x,y,n,ed1,ed2,numEdges);
160 /*for(i=1;i<=2*n-2;i++)
161 {Rprintf("ord[%i]=%i ",i,ord[i]);
164 //look for the edge on the path x to y to subdivide
171 int subdiv=-1;//index of edge to subdivide
172 //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)]);
173 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
174 // 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);
176 //Rprintf("Adding leaf %i, between %i and %i\n",z,x,y);
177 // Rprintf("%i situated at a distance of %d from tree",z,minDist);
179 //Rprintf("path between %i and %i\n",x,y);
181 while(p!=y && sum<lx)
184 //error("%i -> %i ",p,ord[p]);
187 for(i=0;i<=numEdges;i++)
189 if((ed1[i]==aux && ed2[i]==p)||(ed2[i]==aux && ed1[i]==p))
191 if(ed2[i]==aux && ed1[i]==p){sw=1;}
196 //if(cc>1000)error("failed to follow path between x=%i y=%i\n",x,y);
201 //subdivide subdiv with a node labelled nv
206 edLen[subdiv]= (sw==1)?(lx-prevSum):(sum-lx);//check which 'half' of the
207 //path the leaf belongs to
208 //and updates accordingly
209 //error("sum=%f, prevsum=%f\n",sum,prevSum);
210 //error("lx-prevSum=%f, sum-lx=%f, minDist=%f",lx-prevSum,sum-lx,minDist);
211 //Rprintf("adding %i on path %i %i, at distance %f from %i, and %f from tree\n",z,x,y,lx,x,minDist);
212 // Rprintf("subdividing edge %i\n",subdiv);
216 edLen[numEdges]= (sw==1)?(sum-lx):(lx-prevSum);
218 edLen[numEdges]=minDist;