src/nj.c

   1 /* nj.c       2009-07-09 */
   2
   3 /* Copyright 2006-2009 Emmanuel Paradis
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include <R.h>
   9
  10 #define DINDEX(i, j) n*(i - 1) - i*(i - 1)/2 + j - i - 1
  11
  12 int give_index(int i, int j, int n)
  13 {
  14         if (i > j) return(DINDEX(j, i));
  15         else return(DINDEX(i, j));
  16 }
  17
  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 */
  21 {
  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):
  25
  26            i
  27      1  2  3  4  5
  28
  29   2  0
  30   3  1  5
  31 j 4  2  6  9
  32   5  3  7 10 12
  33   6  4  8 11 13 14
  34
  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 */
  37
  38         double sum = 0;
  39         int j, start, end;
  40
  41         if (i < n) {
  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;
  46                 end = start + n - i;
  47                 for (j = start; j < end; j++) sum += D[j];
  48         }
  49
  50         if (i > 1) {
  51                 start = i - 2;
  52                 for (j = 1; j <= i - 1; j++) {
  53                         sum += D[start];
  54                         start += n - j - 1;
  55                 }
  56         }
  57
  58         return(sum);
  59 }
  60
  61 #define GET_I_AND_J                                               \
  62 /* Find the 'R' indices of the two corresponding OTUs */          \
  63 /* The indices of the first element of the pair in the            \
  64    distance matrix are n-1 times 1, n-2 times 2, n-3 times 3,     \
  65    ..., once n-1. Given this, the algorithm below is quite        \
  66    straightforward.*/                                             \
  67     i = 0;                                                        \
  68     for (OTU1 = 1; OTU1 < n; OTU1++) {                            \
  69         i += n - OTU1;                                            \
  70         if (i >= smallest + 1) break;                             \
  71     }                                                             \
  72     /* Finding the second OTU is easier! */                       \
  73     OTU2 = smallest + 1 + OTU1 - n*(OTU1 - 1) + OTU1*(OTU1 - 1)/2
  74
  75 #define SET_CLADE                           \
  76 /* give the node and tip numbers to edge */ \
  77     edge2[k] = otu_label[OTU1 - 1];         \
  78     edge2[k + 1] = otu_label[OTU2 - 1];     \
  79     edge1[k] = edge1[k + 1] = cur_nod
  80
  81 void nj(double *D, int *N, int *edge1, int *edge2, double *edge_length)
  82 {
  83         double SUMD, *S, Sdist, Ndist, *new_dist, A, B, *DI, d_i, x, y;
  84         int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
  85
  86         S = &Sdist;
  87         new_dist = &Ndist;
  88         otu_label = &o_l;
  89         DI = &d_i;
  90
  91         n = *N;
  92         cur_nod = 2*n - 2;
  93
  94         S = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  95         new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  96         otu_label = (int*)R_alloc(n, sizeof(int));
  97         DI = (double*)R_alloc(n - 2, sizeof(double));
  98
  99         for (i = 0; i < n; i++) otu_label[i] = i + 1;
 100         k = 0;
 101
 102         while (n > 3) {
 103
 104                 SUMD = 0;
 105                 for (i = 0; i < n*(n - 1)/2; i++) SUMD += D[i];
 106
 107                 ij = 0;
 108                 for (i = 1; i < n; i++) {
 109                         for (j = i + 1; j <= n; j++) {
 110                                 A = sum_dist_to_i(n, D, i) - D[ij];
 111                                 B = sum_dist_to_i(n, D, j) - D[ij];
 112                                 S[ij] = (A + B)/(2*n - 4) + 0.5*D[ij]
 113                                         + (SUMD - A - B - D[ij])/(n - 2);
 114                                 ij++;
 115                         }
 116                 }
 117
 118                 /* find the 'C' index of the smallest value of S */
 119                 smallest = 0;
 120                 for (i = 1; i < n*(n - 1)/2; i++)
 121                         if (S[smallest] > S[i]) smallest = i;
 122
 123                 GET_I_AND_J;
 124                 SET_CLADE;
 125
 126                 /* get the distances between all OTUs but the 2 selected ones
 127                    and the latter:
 128                    a) get the sum for both
 129                    b) compute the distances for the new OTU */
 130                 A = B = ij = 0;
 131                 for (i = 1; i <= n; i++) {
 132                         if (i == OTU1 || i == OTU2) continue;
 133                         x = D[give_index(i, OTU1, n)]; /* dist between OTU1 and i */
 134                         y = D[give_index(i, OTU2, n)]; /* dist between OTU2 and i */
 135                         new_dist[ij] = (x + y)/2;
 136                         A += x;
 137                         B += y;
 138                         ij++;
 139                 }
 140                 /* compute the branch lengths */
 141                 A /= n - 2;
 142                 B /= n - 2;
 143                 edge_length[k] = (D[smallest] + A - B)/2;
 144                 edge_length[k + 1] = (D[smallest] + B - A)/2;
 145                 DI[cur_nod - *N - 1] = D[smallest];
 146
 147                 /* update before the next loop */
 148                 if (OTU1 > OTU2) { /* make sure that OTU1 < OTU2 */
 149                         i = OTU1;
 150                         OTU1 = OTU2;
 151                         OTU2 = i;
 152                 }
 153                 if (OTU1 != 1)
 154                         for (i = OTU1 - 1; i > 0; i--) otu_label[i] = otu_label[i - 1];
 155                 if (OTU2 != n)
 156                         for (i = OTU2; i <= n; i++) otu_label[i - 1] = otu_label[i];
 157                 otu_label[0] = cur_nod;
 158
 159                 for (i = 1; i < n; i++) {
 160                         if (i == OTU1 || i == OTU2) continue;
 161                         for (j = i + 1; j <= n; j++) {
 162                                 if (j == OTU1 || j == OTU2) continue;
 163                                 new_dist[ij] = D[DINDEX(i, j)];
 164                                 ij++;
 165                         }
 166                 }
 167
 168                 n--;
 169                 for (i = 0; i < n*(n - 1)/2; i++) D[i] = new_dist[i];
 170
 171                 cur_nod--;
 172                 k = k + 2;
 173         }
 174
 175         for (i = 0; i < 3; i++) {
 176                 edge1[*N*2 - 4 - i] = cur_nod;
 177                 edge2[*N*2 - 4 - i] = otu_label[i];
 178         }
 179
 180         edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
 181         edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
 182         edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;
 183
 184         for (i = 0; i < *N*2 - 3; i++) {
 185                 if (edge2[i] <= *N) continue;
 186                 /* In case there are zero branch lengths: */
 187                 if (DI[edge2[i] - *N - 1] == 0) continue;
 188                 edge_length[i] -= DI[edge2[i] - *N - 1]/2;
 189         }
 190 }