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 mvr(double *D, double* v,int *N, int *edge1, int *edge2, double *edge_length)
12 double *S, Sdist, *new_v, 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;
21 S = (double*)R_alloc(n + 1, sizeof(double));
22 new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
23 new_v = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
24 otu_label = (int*)R_alloc(n + 1, sizeof(int));
26 for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
36 Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
41 for (i = 1; i <= n; i++)
45 //Rprintf("index(%i,%i)=%i\n",i,j,give_index(i,j,n));
46 //Rprintf("D[%i,%i]=%f\n",i,j,D[give_index(i,j,n)]);
47 sum+=D[give_index(i,j,n)];
51 //Rprintf("S[%i]=%f\n",i,S[i]);
57 for (i = 1; i < n; i++) {
58 for (j = i + 1; j <= n; j++) {
60 A = B*D[ij] - S[i] - S[j];
61 /*Rprintf("D[ij]=%f\n",D[ij]);
62 Rprintf("S[%i]=%f\n",i,S[i]);
63 Rprintf("S[%i]=%f\n",j,S[j]);
65 Rprintf("A=%f\n",(B*D[ij] - S[i] - S[j]));
66 Rprintf("Q[%i,%i]=%f\n",i,j,A);*/
77 //Rprintf("agglomerating %i and %i, Q=%f \n",OTU1,OTU2,smallest_S);
83 Rprintf("d[%i,%i]=%f ",i,j,D[give_index(i,j,n)]);
92 Rprintf("v[%i,%i]=%f ",i,j,v[give_index(i,j,n)]);
97 edge2[k] = otu_label[OTU1];
98 edge2[k + 1] = otu_label[OTU2];
99 edge1[k] = edge1[k + 1] = cur_nod;
101 /* get the distances between all OTUs but the 2 selected ones
103 a) get the sum for both
104 b) compute the distances for the new OTU */
109 if(i == OTU1 || i==OTU2)continue;
110 //Rprintf("index(%i,%i)=%i index(%i,%i)=%i",i,OTU1,give_index(i,OTU1,n),i,OTU2,give_index(i,OTU2,n));
111 miuSum+=(1/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]));
119 if(i == OTU1 || i==OTU2)continue;
121 double wi=miu/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]);
122 eLenSum+=wi*(D[give_index(i,OTU1,n)]-D[give_index(i,OTU2,n)]);
125 edge_length[k]=D[give_index(OTU1,OTU2,n)]/2 + eLenSum;
130 if(i == OTU1 || i==OTU2)continue;
132 double wi=miu/(v[give_index(i,OTU1,n)]+v[give_index(i,OTU2,n)]);
133 eLenSum+=wi*(D[give_index(i,OTU2,n)]-D[give_index(i,OTU1,n)]);
136 edge_length[k+1]=D[give_index(OTU1,OTU2,n)] - edge_length[k];
140 for (i = 1; i <= n; i++) {
141 if (i == OTU1 || i == OTU2) continue;
142 double xi = D[give_index(i, OTU1, n)]; /* dist between OTU1 and i */
143 double yi = D[give_index(i, OTU2, n)]; /* dist between OTU2 and i */
144 double lamb=v[give_index(i,OTU2,n)]/(v[give_index(i,OTU2,n)]+v[give_index(i,OTU1,n)]);
145 new_dist[ij] = lamb*(xi-edge_length[k])+(1-lamb)*(yi-edge_length[k+1]);
146 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)]);
149 /* compute the branch lengths */
150 //Rprintf("l2:%f \n",edge_length[k+1]);
151 /* update before the next loop
152 (we are sure that OTU1 < OTU2) */
154 for (i = OTU1; i > 1; i--)
155 otu_label[i] = otu_label[i - 1];
157 for (i = OTU2; i < n; i++)
158 otu_label[i] = otu_label[i + 1];
159 otu_label[1] = cur_nod;
161 for (i = 1; i < n; i++) {
162 if (i == OTU1 || i == OTU2) continue;
163 for (j = i + 1; j <= n; j++) {
164 if (j == OTU1 || j == OTU2) continue;
165 new_dist[ij] = D[DINDEX(i, j)];
166 new_v[ij]=v[give_index(i,j,n)];
172 for (i = 0; i < n*(n - 1)/2; i++)
180 for (i = 0; i < 3; i++) {
181 edge1[*N*2 - 4 - i] = cur_nod;
182 edge2[*N*2 - 4 - i] = otu_label[i + 1];
185 edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
186 edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
187 edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;