src/nj.c

   1 /* nj.c       2009-07-17 */
   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 void nj(double *D, int *N, int *edge1, int *edge2, double *edge_length)
  62 {
  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;
  65
  66         S = &Sdist;
  67         new_dist = &Ndist;
  68         otu_label = &o_l;
  69         DI = &d_i;
  70
  71         n = *N;
  72         cur_nod = 2*n - 2;
  73
  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));
  78
  79         for (i = 0; i < n; i++) otu_label[i] = i + 1;
  80         k = 0;
  81
  82         while (n > 3) {
  83
  84                 for (i = 0; i < n; i++)
  85                         S[i] = sum_dist_to_i(n, D, i + 1);
  86
  87                 ij = 0;
  88                 smallest_S = 1e50;
  89                 B = n - 2;
  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;
  93                                 if (A < smallest_S) {
  94                                         OTU1 = i + 1;
  95                                         OTU2 = j + 1;
  96                                         smallest_S = A;
  97                                         smallest = ij;
  98                                 }
  99                                 ij++;
 100                         }
 101                 }
 102
 103                 edge2[k] = otu_label[OTU1 - 1];
 104                 edge2[k + 1] = otu_label[OTU2 - 1];
 105                 edge1[k] = edge1[k + 1] = cur_nod;
 106
 107                 /* get the distances between all OTUs but the 2 selected ones
 108                    and the latter:
 109                    a) get the sum for both
 110                    b) compute the distances for the new OTU */
 111                 A = B = ij = 0;
 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;
 117                         A += x;
 118                         B += y;
 119                         ij++;
 120                 }
 121                 /* compute the branch lengths */
 122                 A /= n - 2;
 123                 B /= n - 2;
 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];
 127
 128                 /* update before the next loop */
 129                 if (OTU1 > OTU2) { /* make sure that OTU1 < OTU2 */
 130                         i = OTU1;
 131                         OTU1 = OTU2;
 132                         OTU2 = i;
 133                 }
 134                 if (OTU1 != 1)
 135                         for (i = OTU1 - 1; i > 0; i--)
 136                                 otu_label[i] = otu_label[i - 1];
 137                 if (OTU2 != n)
 138                         for (i = OTU2; i <= n; i++)
 139                                 otu_label[i - 1] = otu_label[i];
 140                 otu_label[0] = cur_nod;
 141
 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)];
 147                                 ij++;
 148                         }
 149                 }
 150
 151                 n--;
 152                 for (i = 0; i < n*(n - 1)/2; i++) D[i] = new_dist[i];
 153
 154                 cur_nod--;
 155                 k = k + 2;
 156         }
 157
 158         for (i = 0; i < 3; i++) {
 159                 edge1[*N*2 - 4 - i] = cur_nod;
 160                 edge2[*N*2 - 4 - i] = otu_label[i];
 161         }
 162
 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;
 166
 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;
 172         }
 173 }