src/nj.c

   1 /* nj.c       2011-10-20 */
   2
   3 /* Copyright 2006-2011 Emmanuel Paradis
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include "ape.h"
   9
  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 */
  13 {
  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):
  17
  18            i
  19      1  2  3  4  5
  20
  21   2  0
  22   3  1  5
  23 j 4  2  6  9
  24   5  3  7 10 12
  25   6  4  8 11 13 14
  26
  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 */
  29
  30         double sum = 0;
  31         int j, start, end;
  32
  33         if (i < n) {
  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;
  38                 end = start + n - i;
  39                 for (j = start; j < end; j++) sum += D[j];
  40         }
  41
  42         if (i > 1) {
  43                 start = i - 2;
  44                 for (j = 1; j <= i - 1; j++) {
  45                         sum += D[start];
  46                         start += n - j - 1;
  47                 }
  48         }
  49
  50         return(sum);
  51 }
  52
  53 void nj(double *D, int *N, int *edge1, int *edge2, double *edge_length)
  54 {
  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;
  57
  58         S = &Sdist;
  59         new_dist = &Ndist;
  60         otu_label = &o_l;
  61
  62         n = *N;
  63         cur_nod = 2*n - 2;
  64
  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));
  68
  69         for (i = 1; i <= n; i++) otu_label[i] = i; /* otu_label[0] is not used */
  70
  71         k = 0;
  72
  73         while (n > 3) {
  74
  75                 for (i = 1; i <= n; i++)
  76                         S[i] = sum_dist_to_i(n, D, i); /* S[0] is not used */
  77
  78                 ij = 0;
  79                 smallest_S = 1e50;
  80                 B = n - 2;
  81                 for (i = 1; i < n; i++) {
  82                         for (j = i + 1; j <= n; j++) {
  83                                 A = B*D[ij] - S[i] - S[j];
  84                                 if (A < smallest_S) {
  85                                         OTU1 = i;
  86                                         OTU2 = j;
  87                                         smallest_S = A;
  88                                         smallest = ij;
  89                                 }
  90                                 ij++;
  91                         }
  92                 }
  93
  94                 edge2[k] = otu_label[OTU1];
  95                 edge2[k + 1] = otu_label[OTU2];
  96                 edge1[k] = edge1[k + 1] = cur_nod;
  97
  98                 /* get the distances between all OTUs but the 2 selected ones
  99                    and the latter:
 100                    a) get the sum for both
 101                    b) compute the distances for the new OTU */
 102
 103                 A = D[smallest];
 104                 ij = 0;
 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;
 110                         ij++;
 111                 }
 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;
 116
 117                 /* update before the next loop
 118                    (we are sure that OTU1 < OTU2) */
 119                 if (OTU1 != 1)
 120                         for (i = OTU1; i > 1; i--)
 121                                 otu_label[i] = otu_label[i - 1];
 122                 if (OTU2 != n)
 123                         for (i = OTU2; i < n; i++)
 124                                 otu_label[i] = otu_label[i + 1];
 125                 otu_label[1] = cur_nod;
 126
 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)];
 132                                 ij++;
 133                         }
 134                 }
 135
 136                 n--;
 137                 for (i = 0; i < n*(n - 1)/2; i++) D[i] = new_dist[i];
 138
 139                 cur_nod--;
 140                 k = k + 2;
 141         }
 142
 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];
 146         }
 147
 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;
 151 }