1 /* bionjs.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. */
10 void bionjs(double *D, int *N, int *edge1, int *edge2, double *edge_length,int* fsS)
11 { //assume missing values are denoted by -1
12 double *S,*R , *v,*new_v, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
13 int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
14 /*for(i=0;i<n*(n-1)/2;i++)
15 {if(isNA(D[i])){D[i]=-1;}
17 int *s;//s contains |Sxy|, which is all we need for agglomeration
28 R = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
29 v = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
30 new_v = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
31 S = (double*)R_alloc(n + 1, sizeof(double));
32 newR = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
33 new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
34 otu_label = (int*)R_alloc(n + 1, sizeof(int));
35 s = (int*)R_alloc(n*(n - 1)/2, sizeof(int));
36 newS = (int*)R_alloc(n*(n - 1)/2, sizeof(int));
38 for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
41 //populate the v matrix
45 v[give_index(i,j,n)]=D[give_index(i,j,n)];
48 for(i=0;i<n*(n-1)/2;i++)
57 {//algorithm assumes i,j /in Sij, so skip pair if it is not known
58 if(D[give_index(i,j,n)]==-1)
63 {//ij is the pair for which we compute
64 //skip k if we do not know the distances between it and i AND j
68 /* added 2012-04-02: */
69 if(i!=k)R[give_index(i,j,n)]+=D[give_index(i,k,n)];
70 if(j!=k)R[give_index(i,j,n)]+=D[give_index(j,k,n)];
72 s[give_index(i,j,n)]++;
75 if(D[give_index(i,k,n)]==-1 || D[give_index(j,k,n)]==-1)continue;
77 s[give_index(i,j,n)]++;
78 R[give_index(i,j,n)]+=D[give_index(i,k,n)];
79 R[give_index(i,j,n)]+=D[give_index(j,k,n)];
87 Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
96 Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
102 int sw=1;//if 1 then incomplete
108 {newR[give_index(i,j,n)]=0;
109 newS[give_index(i,j,n)]=0;
120 choosePair(D,n,R,s,&sw,&OTU1,&OTU2,fS);
122 else{ //Rprintf("distance matrix is now complete\n");
126 //Rprintf("give_index(%i,%i)=%i\n",i,j,give_index(i,j,n));
127 //Rprintf("D[%i,%i]=%f\n",i,j,D[give_index(i,j,n)]);
128 S[i]+=D[give_index(i,j,n)];
131 //Rprintf("n=%i,B=%f",n,B);
132 for (i = 1; i < n; i++) {
133 for (j = i + 1; j <= n; j++) {
134 //Rprintf("S[%i]=%f, S[%i]=%f, D[%i,%i]=%f, B=%f",i,S[i],j,S[j],i,j,D[give_index(i,j,n)],B);
135 A=S[i]+S[j]-B*D[give_index(i,j,n)];
136 //Rprintf("Q[%i,%i]=%f\n",i,j,A);
137 if (A > smallest_S) {
148 //Rprintf("agglomerating %i and %i, Q=%f \n",OTU1,OTU2,smallest_S);
154 Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
163 Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
172 Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
176 //update R and S, only if matrix still incomplete
179 {if(i==OTU1 || i==OTU2)continue;
181 {if(j==OTU1 || j==OTU2)continue;
182 if(D[give_index(i,j,n)]==-1)continue;
183 if(D[give_index(i,OTU1,n)]!=-1 && D[give_index(j,OTU1,n)]!=-1)
184 {//OTU1 was considered for Rij, so now subtract
185 R[give_index(i,j,n)]-=(D[give_index(i,OTU1,n)]+D[give_index(j,OTU1,n)]);
186 s[give_index(i,j,n)]--;
188 if(D[give_index(i,OTU2,n)]!=-1 && D[give_index(j,OTU2,n)]!=-1)
189 {//OTU2 was considered for Rij, so now subtract
190 R[give_index(i,j,n)]-=(D[give_index(i,OTU2,n)]+D[give_index(j,OTU2,n)]);
191 s[give_index(i,j,n)]--;
196 edge2[k] = otu_label[OTU1];
197 edge2[k + 1] = otu_label[OTU2];
198 edge1[k] = edge1[k + 1] = cur_nod;
201 double lamb=0;//the parameter used for matrix reduction
204 {if(i==OTU1 || i==OTU2)continue;
205 if(D[give_index(OTU1,i,n)]==-1 || D[give_index(OTU2,i,n)]==-1)continue;
206 sum+=(D[give_index(OTU1,i,n)]-D[give_index(OTU2,i,n)]);
207 lambSum+=(v[give_index(OTU2,i,n)]-v[give_index(OTU1,i,n)]);
209 //if we stil have incomplete distances
212 lamb=0.5+(1/(2*(s[give_index(OTU1,OTU2,n)]-2)*v[give_index(OTU2,OTU1,n)]))*lambSum;
214 lamb=0.5+(1/(2*(n-2)*v[give_index(OTU2,OTU1,n)]))*lambSum;
217 //although s was updated above, s[otu1,otu2] has remained unchanged
218 //so it is safe to use it here
219 //if complete distanes, use N-2, else use S
221 if(sw==1){down=s[give_index(OTU1,OTU2,n)]-2;}
223 {error("distance information insufficient to construct a tree, leaves %i and %i isolated from tree",OTU1,OTU2);
225 //Rprintf("down=%f\n",B);
226 sum*=(1.0/(2*(down)));
227 //Rprintf("sum=%f\n",sum);
228 double dxy=D[give_index(OTU1,OTU2,n)]/2;
230 //Rprintf("R[%i,%i]:%f \n",OTU1,OTU2,sum);
231 edge_length[k] = dxy+sum;//OTU1
232 //Rprintf("l1:%f \n",edge_length[k]);
233 edge_length[k + 1] = dxy-sum;//OTU2
234 //Rprintf("l2:%f \n",edge_length[k+1]);
235 //no need to change distance matrix update for complete distance
236 //case, as pairs will automatically fall in the right cathegory
238 //OTU1=x, OTU2=y from formulas
239 A = D[give_index(OTU1,OTU2,n)];
241 for (i = 1; i <= n; i++) {
242 if (i == OTU1 || i == OTU2) continue;
243 if(D[give_index(OTU1,i,n)]!=-1 && D[give_index(OTU2,i,n)]!=-1)
245 new_dist[ij]= lamb*(D[give_index(OTU1,i,n)]-edge_length[k])+(1-lamb)*(D[give_index(OTU2,i,n)]-edge_length[k+1]);
246 new_v[ij]=lamb*v[give_index(OTU1,i,n)]+(1-lamb)*v[give_index(OTU2,i,n)]-lamb*(1-lamb)*v[give_index(OTU1,OTU2,n)];
248 if(D[give_index(OTU1,i,n)]!=-1)
250 new_dist[ij]=D[give_index(OTU1,i,n)]-edge_length[k];
251 new_v[ij]=v[give_index(OTU1,i,n)];
253 if(D[give_index(OTU2,i,n)]!=-1)
255 new_dist[ij]=D[give_index(OTU2,i,n)]-edge_length[k+1];
256 new_v[ij]=v[give_index(OTU2,i,n)];
257 }else{new_dist[ij]=-1;new_v[ij]=-1;}
264 for (i = 1; i < n; i++) {
265 if (i == OTU1 || i == OTU2) continue;
266 for (j = i + 1; j <= n; j++) {
267 if (j == OTU1 || j == OTU2) continue;
268 new_dist[ij] = D[DINDEX(i, j)];
269 new_v[ij]=v[give_index(i,j,n)];
276 for(j=i+1;j<=n-1;j++)
277 {Rprintf("%f ",new_dist[give_index(i,j,n-1)]);
281 //compute Rui, only if distance matrix is still incomplete
287 if(new_dist[give_index(i,1,n-1)]==-1)continue;
293 /* added 2012-04-02 */
294 if(i!=j)newR[give_index(1,i,n-1)]+=new_dist[give_index(i,j,n-1)];
295 if(1!=j)newR[give_index(1,i,n-1)]+=new_dist[give_index(1,j,n-1)];
296 /* end of addition */
297 newS[give_index(1,i,n-1)]++;
300 if(new_dist[give_index(i,j,n-1)]!=-1 && new_dist[give_index(1,j,n-1)]!=-1)
302 newS[give_index(1,i,n-1)]++;
303 newR[give_index(1,i,n-1)]+=new_dist[give_index(i,j,n-1)];
304 newR[give_index(1,i,n-1)]+=new_dist[give_index(1,j,n-1)];
308 //fill in the rest of R and S, again only if distance matrix still
312 {if(i==OTU1 || i==OTU2)continue;
314 {if(j==OTU1 || j==OTU2)continue;
315 newR[ij]=R[give_index(i,j,n)];
316 newS[ij]=s[give_index(i,j,n)];
320 //update newR and newS with the new taxa, again only if distance
321 //matrix is still incomplete
324 {if(new_dist[give_index(1,i,n-1)]==-1)continue;
325 for(j=i+1;j<=n-1;j++)
326 {if(new_dist[give_index(1,j,n-1)]==-1)continue;
327 if(new_dist[give_index(i,j,n-1)]==-1)continue; /* added 2012-04-02 */
328 newR[give_index(i,j,n-1)]+=(new_dist[give_index(1,i,n-1)]+new_dist[give_index(1,j,n-1)]);
329 newS[give_index(i,j,n-1)]++;
332 /* compute the branch lengths */
336 /* update before the next loop
337 (we are sure that OTU1 < OTU2) */
339 for (i = OTU1; i > 1; i--)
340 otu_label[i] = otu_label[i - 1];
342 for (i = OTU2; i < n; i++)
343 otu_label[i] = otu_label[i + 1];
344 otu_label[1] = cur_nod;
349 for (i = 0; i < n*(n - 1)/2; i++)
362 int dK=0;//number of known distances in final distance matrix
363 int iUK=-1;//index of unkown distance, if we have one missing distance
364 int iK=-1;//index of only known distance, only needed if dK==1
365 for (i = 0; i < 3; i++) {
366 edge1[*N*2 - 4 - i] = cur_nod;
367 edge2[*N*2 - 4 - i] = otu_label[i + 1];
368 if(D[i]!=-1){dK++;iK=i;}else{iUK=i;}
371 {//if two distances are known: assume our leaves are x,y,z, d(x,z) unknown
372 //and edge weights of three edges are a,b,c, then any b,c>0 that
373 //satisfy c-b=d(y,z)-d(x,y) a+c=d(y,z) are good edge weights, but for
374 //simplicity we assume a=c if d(yz)<d(xy) a=b otherwise, and after some
375 //algebra we get that we can set the missing distance equal to the
376 //maximum of the already present distances
380 if(D[i]>max)max=D[i];
385 {//through similar motivation as above, if we have just one known distance
386 //we set the other two distances equal to it
393 {//no distances are known, we just set them to 1
398 edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
399 edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
400 edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;