src/nj.c

   1 /* nj.c       2010-04-21 */
   2
   3 /* Copyright 2006-2010 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 /* works if i < j strictly, and i = 1, ...;
  12    see give_index() below */
  13
  14 int give_index(int i, int j, int n)
  15 {
  16         if (i > j) return(DINDEX(j, i));
  17         else return(DINDEX(i, j));
  18 }
  19
  20 double sum_dist_to_i(int n, double *D, int i)
  21 /* returns the sum of all distances D_ij between i and j
  22    with j = 1...n and j != i */
  23 {
  24 /* we use the fact that the distances are arranged sequentially
  25    in the lower triangle, e.g. with n = 6 the 15 distances are
  26    stored as (the C indices are indicated):
  27
  28            i
  29      1  2  3  4  5
  30
  31   2  0
  32   3  1  5
  33 j 4  2  6  9
  34   5  3  7 10 12
  35   6  4  8 11 13 14
  36
  37   so that we sum the values of the ith column--1st loop--and those of
  38   (i - 1)th row (labelled 'i')--2nd loop */
  39
  40         double sum = 0;
  41         int j, start, end;
  42
  43         if (i < n) {
  44                 /* the expression below CANNOT be factorized
  45                    because of the integer operations (it took
  46                    me a while to find out...) */
  47                 start = n*(i - 1) - i*(i - 1)/2;
  48                 end = start + n - i;
  49                 for (j = start; j < end; j++) sum += D[j];
  50         }
  51
  52         if (i > 1) {
  53                 start = i - 2;
  54                 for (j = 1; j <= i - 1; j++) {
  55                         sum += D[start];
  56                         start += n - j - 1;
  57                 }
  58         }
  59
  60         return(sum);
  61 }
  62
  63 void nj(double *D, int *N, int *edge1, int *edge2, double *edge_length)
  64 {
  65         double *S, Sdist, Ndist, *new_dist, A, B, smallest_S, x, y;
  66         int n, i, j, k, ij, smallest, OTU1, OTU2, cur_nod, o_l, *otu_label;
  67
  68         S = &Sdist;
  69         new_dist = &Ndist;
  70         otu_label = &o_l;
  71
  72         n = *N;
  73         cur_nod = 2*n - 2;
  74
  75         S = (double*)R_alloc(n + 1, sizeof(double));
  76         new_dist = (double*)R_alloc(n*(n - 1)/2, sizeof(double));
  77         otu_label = (int*)R_alloc(n + 1, sizeof(int));
  78
  79         for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
  80
  81         k = 0;
  82
  83         while (n > 3) {
  84
  85                 for (i = 1; i <= n; i++)
  86                         S[i] = sum_dist_to_i(n, D, i); /* S[0] is not used */
  87
  88                 ij = 0;
  89                 smallest_S = 1e50;
  90                 B = n - 2;
  91                 for (i = 1; i < n; i++) {
  92                         for (j = i + 1; j <= n; j++) {
  93                                 A = B*D[ij] - S[i] - S[j];
  94                                 if (A < smallest_S) {
  95                                         OTU1 = i;
  96                                         OTU2 = j;
  97                                         smallest_S = A;
  98                                         smallest = ij;
  99                                 }
 100                                 ij++;
 101                         }
 102                 }
 103
 104                 edge2[k] = otu_label[OTU1];
 105                 edge2[k + 1] = otu_label[OTU2];
 106                 edge1[k] = edge1[k + 1] = cur_nod;
 107
 108                 /* get the distances between all OTUs but the 2 selected ones
 109                    and the latter:
 110                    a) get the sum for both
 111                    b) compute the distances for the new OTU */
 112
 113                 A = D[smallest];
 114                 ij = 0;
 115                 for (i = 1; i <= n; i++) {
 116                         if (i == OTU1 || i == OTU2) continue;
 117                         x = D[give_index(i, OTU1, n)]; /* dist between OTU1 and i */
 118                         y = D[give_index(i, OTU2, n)]; /* dist between OTU2 and i */
 119                         new_dist[ij] = (x + y - A)/2;
 120                         ij++;
 121                 }
 122                 /* compute the branch lengths */
 123                 B = (S[OTU1] - S[OTU2])/B; /* don't need B anymore */
 124                 edge_length[k] = (A + B)/2;
 125                 edge_length[k + 1] = (A - B)/2;
 126
 127                 /* update before the next loop
 128                    (we are sure that OTU1 < OTU2) */
 129                 if (OTU1 != 1)
 130                         for (i = OTU1; i > 1; i--)
 131                                 otu_label[i] = otu_label[i - 1];
 132                 if (OTU2 != n)
 133                         for (i = OTU2; i < n; i++)
 134                                 otu_label[i] = otu_label[i + 1];
 135                 otu_label[1] = cur_nod;
 136
 137                 for (i = 1; i < n; i++) {
 138                         if (i == OTU1 || i == OTU2) continue;
 139                         for (j = i + 1; j <= n; j++) {
 140                                 if (j == OTU1 || j == OTU2) continue;
 141                                 new_dist[ij] = D[DINDEX(i, j)];
 142                                 ij++;
 143                         }
 144                 }
 145
 146                 n--;
 147                 for (i = 0; i < n*(n - 1)/2; i++) D[i] = new_dist[i];
 148
 149                 cur_nod--;
 150                 k = k + 2;
 151         }
 152
 153         for (i = 0; i < 3; i++) {
 154                 edge1[*N*2 - 4 - i] = cur_nod;
 155                 edge2[*N*2 - 4 - i] = otu_label[i + 1];
 156         }
 157
 158         edge_length[*N*2 - 4] = (D[0] + D[1] - D[2])/2;
 159         edge_length[*N*2 - 5] = (D[0] + D[2] - D[1])/2;
 160         edge_length[*N*2 - 6] = (D[2] + D[1] - D[0])/2;
 161 }