3 /* Copyright 2006-2009 Emmanuel Paradis
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
10 #define DINDEX(i, j) n*(i - 1) - i*(i - 1)/2 + j - i - 1
12 int give_index(int i, int j, int n)
14 if (i > j) return(DINDEX(j, i));
15 else return(DINDEX(i, j));
18 double sum_dist_to_i(int n, double *D, int i)
19 /* returns the sum of all distances D_ij between i and j
20 with j = 1...n and j != i */
22 /* we use the fact that the distances are arranged sequentially
23 in the lower triangle, e.g. with n = 6 the 15 distances are
24 stored as (the C indices are indicated):
35 so that we sum the values of the ith column-1st loop-and those of
36 (i - 1)th row (labelled 'i')-2nd loop */
42 /* the expression below CANNOT be factorized
43 because of the integer operations (it took
44 me a while to find out...) */
45 start = n*(i - 1) - i*(i - 1)/2;
47 for (j = start; j < end; j++) sum += D[j];
52 for (j = 1; j <= i - 1; j++) {
61 void nj(double *D, int *N, int *edge1, int *edge2, double *edge_length)
63 double *S, Sdist, Ndist, *new_dist, A, B, smallest_S, *DI, d_i, x, y;
64 int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
74 S = (double*)R_alloc(n, sizeof(double));
75 new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
76 otu_label = (int*)R_alloc(n, sizeof(int));
77 DI = (double*)R_alloc(n - 2, sizeof(double));
79 for (i = 0; i < n; i++) otu_label[i] = i + 1;
84 for (i = 0; i < n; i++)
85 S[i] = sum_dist_to_i(n, D, i + 1);
90 for (i = 0; i < n - 1; i++) {
91 for (j = i + 1; j < n; j++) {
92 A = D[ij] - (S[i] + S[j])/B;
103 edge2[k] = otu_label[OTU1 - 1];
104 edge2[k + 1] = otu_label[OTU2 - 1];
105 edge1[k] = edge1[k + 1] = cur_nod;
107 /* get the distances between all OTUs but the 2 selected ones
109 a) get the sum for both
110 b) compute the distances for the new OTU */
112 for (i = 1; i <= n; i++) {
113 if (i == OTU1 || i == OTU2) continue;
114 x = D[give_index(i, OTU1, n)]; /* dist between OTU1 and i */
115 y = D[give_index(i, OTU2, n)]; /* dist between OTU2 and i */
116 new_dist[ij] = (x + y)/2;
121 /* compute the branch lengths */
124 edge_length[k] = (D[smallest] + A - B)/2;
125 edge_length[k + 1] = (D[smallest] + B - A)/2;
126 DI[cur_nod - *N - 1] = D[smallest];
128 /* update before the next loop */
129 if (OTU1 > OTU2) { /* make sure that OTU1 < OTU2 */
135 for (i = OTU1 - 1; i > 0; i--)
136 otu_label[i] = otu_label[i - 1];
138 for (i = OTU2; i < n; i++)
139 otu_label[i - 1] = otu_label[i];
140 otu_label[0] = cur_nod;
142 for (i = 1; i < n; i++) {
143 if (i == OTU1 || i == OTU2) continue;
144 for (j = i + 1; j <= n; j++) {
145 if (j == OTU1 || j == OTU2) continue;
146 new_dist[ij] = D[DINDEX(i, j)];
152 for (i = 0; i < n*(n - 1)/2; i++) D[i] = new_dist[i];
158 for (i = 0; i < 3; i++) {
159 edge1[*N*2 - 4 - i] = cur_nod;
160 edge2[*N*2 - 4 - i] = otu_label[i];
163 edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
164 edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
165 edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;
167 for (i = 0; i < *N*2 - 3; i++) {
168 if (edge2[i] <= *N) continue;
169 /* In case there are zero branch lengths: */
170 if (DI[edge2[i] - *N - 1] == 0) continue;
171 edge_length[i] -= DI[edge2[i] - *N - 1]/2;