1 /* mvrs.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 mvrs(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_length,int* fsS)
11 { //assume missing values are denoted by -1
13 double *S,*R ,*new_v, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
14 int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
15 /*for(i=0;i<n*(n-1)/2;i++)
16 {if(isNA(D[i])){D[i]=-1;}
18 int *s;//s contains |Sxy|, which is all we need for agglomeration
29 R = (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 */
43 for(i=0;i<n*(n-1)/2;i++)
52 {//algorithm assumes i,j /in Sij, so skip pair if it is not known
53 if(D[give_index(i,j,n)]==-1)
58 {//ij is the pair for which we compute
59 //skip k if we do not know the distances between it and i AND j
63 s[give_index(i,j,n)]++;
64 //Rprintf("%i",s[give_index(i,j,n)]);
68 if(D[give_index(i,k,n)]==-1 || D[give_index(j,k,n)]==-1)continue;
70 s[give_index(i,j,n)]++;
71 R[give_index(i,j,n)]+=D[give_index(i,k,n)];
72 R[give_index(i,j,n)]+=D[give_index(j,k,n)];
73 //Rprintf("%i",s[give_index(i,j,n)]);
74 //Rprintf("%f",R[give_index(i,j,n)]);
83 Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
92 Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
97 int sw=1;//if 1 then incomplete
103 {newR[give_index(i,j,n)]=0;
104 newS[give_index(i,j,n)]=0;
116 choosePair(D,n,R,s,&sw,&OTU1,&OTU2,fS);
118 else{ //Rprintf("distance matrix is now complete\n");
122 //Rprintf("give_index(%i,%i)=%i\n",i,j,give_index(i,j,n));
123 //Rprintf("D[%i,%i]=%f\n",i,j,D[give_index(i,j,n)]);
124 S[i]+=D[give_index(i,j,n)];
127 //Rprintf("n=%i,B=%f",n,B);
128 for (i = 1; i < n; i++) {
129 for (j = i + 1; j <= n; j++) {
130 //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);
131 A=S[i]+S[j]-B*D[give_index(i,j,n)];
132 Rprintf("Q[%i,%i]=%f\n",i,j,A);
133 if (A > smallest_S) {
143 if(s[give_index(OTU1,OTU2,n)]<=2)
144 {error("distance information insufficient to construct a tree, leaves %i and %i isolated from tree",OTU1,OTU2);
146 //Rprintf("agglomerating %i and %i, Q=%f \n",OTU1,OTU2,smallest_S);
152 Rprintf("R[%i,%i]=%f ",i,j,R[give_index(i,j,n)]);
161 Rprintf("s[%i,%i]=%i ",i,j,s[give_index(i,j,n)]);
170 Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
175 //update R and S, only if matrix still incomplete
178 {if(i==OTU1 || i==OTU2)continue;
180 {if(j==OTU1 || j==OTU2)continue;
181 if(D[give_index(i,j,n)]==-1)continue;
182 if(D[give_index(i,OTU1,n)]!=-1 && D[give_index(j,OTU1,n)]!=-1)
183 {//OTU1 was considered for Rij, so now subtract
184 R[give_index(i,j,n)]-=(D[give_index(i,OTU1,n)]+D[give_index(j,OTU1,n)]);
185 s[give_index(i,j,n)]--;
187 if(D[give_index(i,OTU2,n)]!=-1 && D[give_index(j,OTU2,n)]!=-1)
188 {//OTU2 was considered for Rij, so now subtract
189 R[give_index(i,j,n)]-=(D[give_index(i,OTU2,n)]+D[give_index(j,OTU2,n)]);
190 s[give_index(i,j,n)]--;
195 edge2[k] = otu_label[OTU1];
196 edge2[k + 1] = otu_label[OTU2];
197 edge1[k] = edge1[k + 1] = cur_nod;
203 if(i == OTU1 || i==OTU2)continue;
204 if(D[give_index(i,OTU1,n)]==-1 || D[give_index(i,OTU2,n)]==-1)continue;
205 //Rprintf("index(%i,%i)=%i index(%i,%i)=%i",i,OTU1,give_index(i,OTU1,n),i,OTU2,give_index(i,OTU2,n));
206 miuSum+=(1/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]));
214 if(i == OTU1 || i==OTU2)continue;
215 if(D[give_index(i,OTU1,n)]==-1 || D[give_index(i,OTU2,n)]==-1)continue;
216 double wi=miu/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]);
217 eLenSum+=wi*(D[give_index(i,OTU1,n)]-D[give_index(i,OTU2,n)]);
220 edge_length[k]=D[give_index(OTU1,OTU2,n)]/2 + eLenSum;
225 if(i == OTU1 || i==OTU2)continue;
226 if(D[give_index(i,OTU1,n)]==-1 || D[give_index(i,OTU2,n)]==-1)continue;
227 double wi=miu/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]);
228 eLenSum+=wi*(D[give_index(i,OTU2,n)]-D[give_index(i,OTU1,n)]);
231 edge_length[k+1]=D[give_index(OTU1,OTU2,n)]/2 + eLenSum;
233 //no need to change distance matrix update for complete distance
234 //case, as pairs will automatically fall in the right cathegory
236 //OTU1=x, OTU2=y from formulas
237 A = D[give_index(OTU1,OTU2,n)];
239 for (i = 1; i <= n; i++) {
240 if (i == OTU1 || i == OTU2) continue;
241 if(D[give_index(OTU1,i,n)]!=-1 && D[give_index(OTU2,i,n)]!=-1)
242 { double lamb=v[give_index(i,OTU2,n)]/(v[give_index(i,OTU2,n)]+v[give_index(i,OTU1,n)]);
243 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]);
244 new_v[ij]=(v[give_index(i,OTU2,n)]*v[give_index(i,OTU1,n)])/(v[give_index(i,OTU2,n)]+v[give_index(i,OTU1,n)]);
246 if(D[give_index(OTU1,i,n)]!=-1)
248 new_dist[ij]=D[give_index(OTU1,i,n)]-edge_length[k];
249 new_v[ij]=v[give_index(OTU1,i,n)];
251 if(D[give_index(OTU2,i,n)]!=-1)
253 new_dist[ij]=D[give_index(OTU2,i,n)]-edge_length[k+1];
254 new_v[ij]=v[give_index(OTU2,i,n)];
255 }else{new_dist[ij]=-1;new_v[ij]=-1;}
262 for (i = 1; i < n; i++) {
263 if (i == OTU1 || i == OTU2) continue;
264 for (j = i + 1; j <= n; j++) {
265 if (j == OTU1 || j == OTU2) continue;
266 new_dist[ij] = D[DINDEX(i, j)];
267 new_v[ij]=v[give_index(i,j,n)];
274 for(j=i+1;j<=n-1;j++)
275 {Rprintf("%f ",new_dist[give_index(i,j,n-1)]);
279 //compute Rui, only if distance matrix is still incomplete
285 if(new_dist[give_index(i,1,n-1)]==-1)continue;
291 newS[give_index(1,i,n-1)]++;
294 if(new_dist[give_index(i,j,n-1)]!=-1 && new_dist[give_index(1,j,n-1)]!=-1)
296 newS[give_index(1,i,n-1)]++;
297 newR[give_index(1,i,n-1)]+=new_dist[give_index(i,j,n-1)];
298 newR[give_index(1,i,n-1)]+=new_dist[give_index(1,j,n-1)];
302 //fill in the rest of R and S, again only if distance matrix still
306 {if(i==OTU1 || i==OTU2)continue;
308 {if(j==OTU1 || j==OTU2)continue;
309 newR[ij]=R[give_index(i,j,n)];
310 newS[ij]=s[give_index(i,j,n)];
314 //update newR and newS with the new taxa, again only if distance
315 //matrix is still incomplete
318 {if(new_dist[give_index(1,i,n-1)]==-1)continue;
319 for(j=i+1;j<=n-1;j++)
320 {if(new_dist[give_index(1,j,n-1)]==-1)continue;
321 newR[give_index(i,j,n-1)]+=(new_dist[give_index(1,i,n-1)]+new_dist[give_index(1,j,n-1)]);
322 newS[give_index(i,j,n-1)]++;
325 /* compute the branch lengths */
329 /* update before the next loop
330 (we are sure that OTU1 < OTU2) */
332 for (i = OTU1; i > 1; i--)
333 otu_label[i] = otu_label[i - 1];
335 for (i = OTU2; i < n; i++)
336 otu_label[i] = otu_label[i + 1];
337 otu_label[1] = cur_nod;
342 for (i = 0; i < n*(n - 1)/2; i++)
355 int dK=0;//number of known distances in final distance matrix
356 int iUK=-1;//index of unkown distance, if we have one missing distance
357 int iK=-1;//index of only known distance, only needed if dK==1
358 for (i = 0; i < 3; i++) {
359 edge1[*N*2 - 4 - i] = cur_nod;
360 edge2[*N*2 - 4 - i] = otu_label[i + 1];
361 if(D[i]!=-1){dK++;iK=i;}else{iUK=i;}
364 {//if two distances are known: assume our leaves are x,y,z, d(x,z) unknown
365 //and edge weights of three edges are a,b,c, then any b,c>0 that
366 //satisfy c-b=d(y,z)-d(x,y) a+c=d(y,z) are good edge weights, but for
367 //simplicity we assume a=c if d(yz)<d(xy) a=b otherwise, and after some
368 //algebra we get that we can set the missing distance equal to the
369 //maximum of the already present distances
373 if(D[i]>max)max=D[i];
378 {//through similar motivation as above, if we have just one known distance
379 //we set the other two distances equal to it
386 {//no distances are known, we just set them to 1
391 edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
392 edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
393 edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;