3 /* Copyright 2006-2011 Emmanuel Paradis
5 /* This file is part of the R-package `ape'. */
6 /* See the file ../COPYING for licensing issues. */
10 double sum_dist_to_i(int n, double *D, int i)
11 /* returns the sum of all distances D_ij between i and j
12 with j = 1...n and j != i */
14 /* we use the fact that the distances are arranged sequentially
15 in the lower triangle, e.g. with n = 6 the 15 distances are
16 stored as (the C indices are indicated):
27 so that we sum the values of the ith column--1st loop--and those of
28 (i - 1)th row (labelled 'i')--2nd loop */
34 /* the expression below CANNOT be factorized
35 because of the integer operations (it took
36 me a while to find out...) */
37 start = n*(i - 1) - i*(i - 1)/2;
39 for (j = start; j < end; j++) sum += D[j];
44 for (j = 1; j <= i - 1; j++) {
53 void nj(double *D, int *N, int *edge1, int *edge2, double *edge_length)
55 double *S, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
56 int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
65 S = (double*)R_alloc(n + 1, sizeof(double));
66 new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
67 otu_label = (int*)R_alloc(n + 1, sizeof(int));
69 for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
75 for (i = 1; i <= n; i++)
76 S[i] = sum_dist_to_i(n, D, i); /* S[0] is not used */
81 for (i = 1; i < n; i++) {
82 for (j = i + 1; j <= n; j++) {
83 A = B*D[ij] - S[i] - S[j];
94 edge2[k] = otu_label[OTU1];
95 edge2[k + 1] = otu_label[OTU2];
96 edge1[k] = edge1[k + 1] = cur_nod;
98 /* get the distances between all OTUs but the 2 selected ones
100 a) get the sum for both
101 b) compute the distances for the new OTU */
105 for (i = 1; i <= n; i++) {
106 if (i == OTU1 || i == OTU2) continue;
107 x = D[give_index(i, OTU1, n)]; /* dist between OTU1 and i */
108 y = D[give_index(i, OTU2, n)]; /* dist between OTU2 and i */
109 new_dist[ij] = (x + y - A)/2;
112 /* compute the branch lengths */
113 B = (S[OTU1] - S[OTU2])/B; /* don't need B anymore */
114 edge_length[k] = (A + B)/2;
115 edge_length[k + 1] = (A - B)/2;
117 /* update before the next loop
118 (we are sure that OTU1 < OTU2) */
120 for (i = OTU1; i > 1; i--)
121 otu_label[i] = otu_label[i - 1];
123 for (i = OTU2; i < n; i++)
124 otu_label[i] = otu_label[i + 1];
125 otu_label[1] = cur_nod;
127 for (i = 1; i < n; i++) {
128 if (i == OTU1 || i == OTU2) continue;
129 for (j = i + 1; j <= n; j++) {
130 if (j == OTU1 || j == OTU2) continue;
131 new_dist[ij] = D[DINDEX(i, j)];
137 for (i = 0; i < n*(n - 1)/2; i++) D[i] = new_dist[i];
143 for (i = 0; i < 3; i++) {
144 edge1[*N*2 - 4 - i] = cur_nod;
145 edge2[*N*2 - 4 - i] = otu_label[i + 1];
148 edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
149 edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
150 edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;