src/dist_dna.c

   1 /* dist_dna.c       2012-02-14 */
   2
   3 /* Copyright 2005-2012 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_ext/Lapack.h>
   9 #include "ape.h"
  10
  11 /* from R: print(log(4), d = 22) */
  12 #define LN4 1.386294361119890572454
  13
  14 /* returns 8 if the base is known surely, 0 otherwise */
  15 #define KnownBase(a) (a & 8)
  16
  17 /* returns 1 if the base is adenine surely, 0 otherwise */
  18 #define IsAdenine(a) a == 136
  19
  20 /* returns 1 if the base is guanine surely, 0 otherwise */
  21 #define IsGuanine(a) a == 72
  22
  23 /* returns 1 if the base is cytosine surely, 0 otherwise */
  24 #define IsCytosine(a) a == 40
  25
  26 /* returns 1 if the base is thymine surely, 0 otherwise */
  27 #define IsThymine(a) a == 24
  28
  29 /* returns 1 if the base is a purine surely, 0 otherwise */
  30 #define IsPurine(a) a > 63
  31
  32 /* returns 1 if the base is a pyrimidine surely, 0 otherwise */
  33 #define IsPyrimidine(a) a < 64
  34
  35 /* returns 1 if both bases are different surely, 0 otherwise */
  36 #define DifferentBase(a, b) (a & b) < 16
  37
  38 /* returns 1 if both bases are the same surely, 0 otherwise */
  39 #define SameBase(a, b) KnownBase(a) && a == b
  40
  41 /* computes directly the determinant of a 4x4 matrix */
  42 double detFourByFour(double *x)
  43 {
  44     double det, a33a44, a34a43, a34a42, a32a44, a32a43, a33a42, a34a41, a31a44, a31a43, a33a41, a31a42, a32a41;
  45
  46     a33a44 = x[10]*x[15]; a34a43 = x[14]*x[11];
  47     a34a42 = x[14]*x[7];  a32a44 = x[6]*x[15];
  48     a32a43 = x[6]*x[11];  a33a42 = x[10]*x[7];
  49     a34a41 = x[14]*x[3];  a31a44 = x[2]*x[15];
  50     a31a43 = x[2]*x[11];  a33a41 = x[10]*x[3];
  51     a31a42 = x[2]*x[7];   a32a41 = x[6]*x[3];
  52
  53     det = x[0]*x[5]*(a33a44 - a34a43) + x[0]*x[9]*(a34a42 - a32a44) +
  54       x[0]*x[13]*(a32a43 - a33a42) + x[4]*x[9]*(a31a44 - a34a41) +
  55       x[4]*x[13]*(a33a41 - a31a43) + x[4]*x[1]*(a34a43 - a33a44) +
  56       x[8]*x[13]*(a31a42 - a32a41) + x[8]*x[1]*(a32a44 - a34a42) +
  57       x[8]*x[5]*(a34a41 - a31a44) + x[12]*x[1]*(a33a42 - a32a43) +
  58       x[12]*x[5]*(a31a43 - a33a41) + x[12]*x[9]*(a32a41 - a31a42);
  59
  60     return det;
  61 }
  62
  63 #define CHECK_PAIRWISE_DELETION\
  64     if (KnownBase(x[s1]) && KnownBase(x[s2])) L++;\
  65     else continue;
  66
  67 #define COUNT_TS_TV\
  68     if (SameBase(x[s1], x[s2])) continue;\
  69     Nd++;\
  70     if (IsPurine(x[s1]) && IsPurine(x[s2])) {\
  71         Ns++;\
  72         continue;\
  73     }\
  74     if (IsPyrimidine(x[s1]) && IsPyrimidine(x[s2])) Ns++;
  75
  76 void distDNA_indel(unsigned char *x, int *n, int *s, double *d)
  77 {
  78         int i1, i2, s1, s2, target, N;
  79
  80         target = 0;
  81         for (i1 = 1; i1 < *n; i1++) {
  82                 for (i2 = i1 + 1; i2 <= *n; i2++) {
  83                         N = 0;
  84
  85                         for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
  86                                 if ((x[s1] ^ x[s2]) & 4) N++;
  87
  88                         d[target] = ((double) N);
  89                         target++;
  90                 }
  91         }
  92 }
  93
  94 #define X(i, j) i - 1 + *n * (j - 1)
  95
  96 #define DINDEX(i, j) n * (i - 1) - i*(i - 1)/2 + j - i - 1
  97
  98 int give_index(int i, int j, int n)
  99 {
 100         if (i > j) return(DINDEX(j, i));
 101         else return(DINDEX(i, j));
 102 }
 103
 104 void distDNA_indelblock(unsigned char *x, int *n, int *s, double *d)
 105 {
 106         int i1, i2, s1, s2, target, N, start_block, end_block;
 107
 108         for (i1 = 1; i1 <= *n; i1++) {
 109
 110 /* find a block of one or more '-' */
 111
 112                 s1 = 1;
 113                 while (s1 < *s) {
 114                         if (x[X(i1, s1)] == 4) {
 115                                 start_block = s1;
 116                                 while (x[X(i1, s1)] == 4) s1++;
 117                                 end_block = s1 - 1;
 118
 119 /* then look if the same block is present in all the other sequences */
 120
 121                                 for (i2 = 1; i2 <= *n; i2++) {
 122                                         if (i2 == i1) continue;
 123
 124                                         target = give_index(i1, i2, *n);
 125
 126                                         if (start_block > 1) {
 127                                                 if (x[X(i2, start_block - 1)] == 4) {
 128                                                         d[target]++;
 129                                                         continue;
 130                                                 }
 131                                         }
 132                                         if (end_block < *s) {
 133                                                 if (x[X(i2, end_block + 1)] == 4) {
 134                                                         d[target]++;
 135                                                         continue;
 136                                                 }
 137                                         }
 138                                         for (s2 = start_block; s2 <= end_block; s2++) {
 139                                                 if (x[X(i2, s2)] != 4) {
 140                                                         d[target]++;
 141                                                         continue;
 142                                                 }
 143                                         }
 144                                 }
 145                                 s1 = end_block + 1;
 146                         }
 147                         s1++;
 148                 }
 149         }
 150 }
 151
 152 #undef X
 153
 154 void distDNA_TsTv(unsigned char *x, int *n, int *s, double *d, int Ts, int pairdel)
 155 {
 156         int i1, i2, s1, s2, target, Nd, Ns;
 157
 158         target = 0;
 159         for (i1 = 1; i1 < *n; i1++) {
 160                 for (i2 = i1 + 1; i2 <= *n; i2++) {
 161                         Nd = Ns = 0;
 162                         for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 163                                 if (pairdel && !(KnownBase(x[s1]) && KnownBase(x[s2]))) continue;
 164                                 COUNT_TS_TV
 165                         }
 166                         if (Ts) d[target] = ((double) Ns); /* output number of transitions */
 167                         else d[target] = ((double) Nd - Ns); /* output number of transversions */
 168                         target++;
 169                 }
 170         }
 171 }
 172
 173 void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled)
 174 {
 175         int i1, i2, s1, s2, target, Nd;
 176
 177         target = 0;
 178         for (i1 = 1; i1 < *n; i1++) {
 179                 for (i2 = i1 + 1; i2 <= *n; i2++) {
 180                         Nd = 0;
 181                         for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
 182                                 if (DifferentBase(x[s1], x[s2])) Nd++;
 183                         if (scaled) d[target] = ((double) Nd / *s);
 184                         else d[target] = ((double) Nd);
 185                         target++;
 186                 }
 187         }
 188 }
 189
 190 void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d, int scaled)
 191 {
 192     int i1, i2, s1, s2, target, Nd, L;
 193
 194     target = 0;
 195     for (i1 = 1; i1 < *n; i1++) {
 196         for (i2 = i1 + 1; i2 <= *n; i2++) {
 197             Nd = L = 0;
 198             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 199                 CHECK_PAIRWISE_DELETION
 200                 if (DifferentBase(x[s1], x[s2])) Nd++;
 201             }
 202             if (scaled) d[target] = ((double) Nd/L);
 203             else d[target] = ((double) Nd);
 204             target++;
 205         }
 206     }
 207 }
 208
 209 #define COMPUTE_DIST_JC69\
 210     p = ((double) Nd/L);\
 211     if (*gamma)\
 212       d[target] = 0.75 * *alpha*(pow(1 - 4*p/3, -1/ *alpha) - 1);\
 213     else d[target] = -0.75*log(1 - 4*p/3);\
 214     if (*variance) {\
 215         if (*gamma) var[target] = p*(1 - p)/(pow(1 - 4*p/3, -2/(*alpha + 1)) * L);\
 216         else var[target] = p*(1 - p)/(pow(1 - 4*p/3, 2)*L);\
 217     }
 218
 219 void distDNA_JC69(unsigned char *x, int *n, int *s, double *d,
 220                   int *variance, double *var, int *gamma, double *alpha)
 221 {
 222     int i1, i2, s1, s2, target, Nd, L;
 223     double p;
 224
 225     L = *s;
 226
 227     target = 0;
 228     for (i1 = 1; i1 < *n; i1++) {
 229         for (i2 = i1 + 1; i2 <= *n; i2++) {
 230             Nd = 0;
 231             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
 232               if (DifferentBase(x[s1], x[s2])) Nd++;
 233             COMPUTE_DIST_JC69
 234             target++;
 235         }
 236     }
 237 }
 238
 239 void distDNA_JC69_pairdel(unsigned char *x, int *n, int *s, double *d,
 240                           int *variance, double *var, int *gamma, double *alpha)
 241 {
 242     int i1, i2, s1, s2, target, Nd, L;
 243     double p;
 244
 245     target = 0;
 246     for (i1 = 1; i1 < *n; i1++) {
 247         for (i2 = i1 + 1; i2 <= *n; i2++) {
 248             Nd = L = 0;
 249             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 250                 CHECK_PAIRWISE_DELETION
 251                 if (DifferentBase(x[s1], x[s2])) Nd++;
 252             }
 253             COMPUTE_DIST_JC69
 254             target++;
 255         }
 256     }
 257 }
 258
 259 #define COMPUTE_DIST_K80\
 260     P = ((double) Ns/L);\
 261     Q = ((double) (Nd - Ns)/L);\
 262     a1 = 1 - 2*P - Q;\
 263     a2 = 1 - 2*Q;\
 264     if (*gamma) {\
 265         b = -1 / *alpha;\
 266         d[target] = *alpha * (pow(a1, b) + 0.5*pow(a2, b) - 1.5)/2;\
 267     }\
 268     else d[target] = -0.5 * log(a1 * sqrt(a2));\
 269     if (*variance) {\
 270         if (*gamma) {\
 271             b = -(1 / *alpha + 1);\
 272             c1 = pow(a1, b);\
 273             c2 = pow(a2, b);\
 274             c3 = (c1 + c2)/2;\
 275         } else {\
 276           c1 = 1/a1;\
 277           c2 = 1/a2;\
 278           c3 = (c1 + c2)/2;\
 279         }\
 280         var[target] = (c1*c1*P + c3*c3*Q - pow(c1*P + c3*Q, 2))/L;\
 281     }
 282
 283 void distDNA_K80(unsigned char *x, int *n, int *s, double *d,
 284                  int *variance, double *var, int *gamma, double *alpha)
 285 {
 286     int i1, i2, s1, s2, target, Nd, Ns, L;
 287     double P, Q, a1, a2, b, c1, c2, c3;
 288
 289     L = *s;
 290
 291     target = 0;
 292     for (i1 = 1; i1 < *n; i1++) {
 293         for (i2 = i1 + 1; i2 <= *n; i2++) {
 294             Nd = Ns = 0;
 295             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 296                 COUNT_TS_TV
 297             }
 298             COMPUTE_DIST_K80
 299             target++;
 300         }
 301     }
 302 }
 303
 304 void distDNA_K80_pairdel(unsigned char *x, int *n, int *s, double *d,
 305                          int *variance, double *var, int *gamma, double *alpha)
 306 {
 307     int i1, i2, s1, s2, target, Nd, Ns, L;
 308     double P, Q, a1, a2, b, c1, c2, c3;
 309
 310     target = 0;
 311     for (i1 = 1; i1 < *n; i1++) {
 312         for (i2 = i1 + 1; i2 <= *n; i2++) {
 313             Nd = Ns = L = 0;
 314             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 315                 CHECK_PAIRWISE_DELETION
 316                 COUNT_TS_TV
 317             }
 318             COMPUTE_DIST_K80
 319             target++;
 320         }
 321     }
 322 }
 323
 324 #define COMPUTE_DIST_F81\
 325     p = ((double) Nd/L);\
 326     if (*gamma) d[target] = E * *alpha * (pow(1 - p/E, -1/ *alpha) - 1);\
 327     else d[target] = -E*log(1 - p/E);\
 328     if (*variance) {\
 329         if (*gamma) var[target] = p*(1 - p)/(pow(1 - p/E, -2/(*alpha + 1)) * L);\
 330         else var[target] = p*(1 - p)/(pow(1 - p/E, 2)*L);\
 331     }
 332
 333 void distDNA_F81(unsigned char *x, int *n, int *s, double *d, double *BF,
 334                  int *variance, double *var, int *gamma, double *alpha)
 335 {
 336     int i1, i2, s1, s2, target, Nd, L;
 337     double p, E;
 338
 339     L = *s;
 340     E = 1 - BF[0]*BF[0] - BF[1]*BF[1] - BF[2]*BF[2] - BF[3]*BF[3];
 341
 342     target = 0;
 343     for (i1 = 1; i1 < *n; i1++) {
 344         for (i2 = i1 + 1; i2 <= *n; i2++) {
 345             Nd = 0;
 346             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n)
 347               if (DifferentBase(x[s1], x[s2])) Nd++;
 348             COMPUTE_DIST_F81
 349             target++;
 350         }
 351     }
 352 }
 353
 354 void distDNA_F81_pairdel(unsigned char *x, int *n, int *s, double *d, double *BF,
 355                          int *variance, double *var, int *gamma, double *alpha)
 356 {
 357     int i1, i2, s1, s2, target, Nd, L;
 358     double p, E;
 359
 360     E = 1 - BF[0]*BF[0] - BF[1]*BF[1] - BF[2]*BF[2] - BF[3]*BF[3];
 361
 362     target = 0;
 363     for (i1 = 1; i1 < *n; i1++) {
 364         for (i2 = i1 + 1; i2 <= *n; i2++) {
 365             Nd = L = 0;
 366             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 367                 CHECK_PAIRWISE_DELETION
 368                 if (DifferentBase(x[s1], x[s2])) Nd++;
 369             }
 370             COMPUTE_DIST_F81
 371             target++;
 372         }
 373     }
 374 }
 375
 376 #define COUNT_TS_TV1_TV2\
 377     if (SameBase(x[s1], x[s2])) continue;\
 378     Nd++;\
 379     if ((x[s1] | x[s2]) == 152 || (x[s1] | x[s2]) == 104) {\
 380         Nv1++;\
 381         continue;\
 382     }\
 383     if ((x[s1] | x[s2]) == 168 || (x[s1] | x[s2]) == 88) Nv2++;
 384
 385
 386 #define COMPUTE_DIST_K81\
 387     P = ((double) (Nd - Nv1 - Nv2)/L);\
 388     Q = ((double) Nv1/L);\
 389     R = ((double) Nv2/L);\
 390     a1 = 1 - 2*P - 2*Q;\
 391     a2 = 1 - 2*P - 2*R;\
 392     a3 = 1 - 2*Q - 2*R;\
 393     d[target] = -0.25*log(a1*a2*a3);\
 394     if (*variance) {\
 395         a = (1/a1 + 1/a2)/2;\
 396         b = (1/a1 + 1/a3)/2;\
 397         c = (1/a2 + 1/a3)/2;\
 398       var[target] = (a*a*P + b*b*Q + c*c*R - pow(a*P + b*Q + c*R, 2))/2;\
 399     }
 400
 401 void distDNA_K81(unsigned char *x, int *n, int *s, double *d,
 402                  int *variance, double *var)
 403 {
 404     int i1, i2, Nd, Nv1, Nv2, L, s1, s2, target;
 405     double P, Q, R, a1, a2, a3, a, b, c;
 406
 407     L = *s;
 408
 409     target = 0;
 410     for (i1 = 1; i1 < *n; i1++) {
 411         for (i2 = i1 + 1; i2 <= *n; i2++) {
 412             Nd = Nv1 = Nv2 = 0;
 413             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 414                 COUNT_TS_TV1_TV2
 415             }
 416             COMPUTE_DIST_K81
 417             target++;
 418         }
 419     }
 420 }
 421
 422 void distDNA_K81_pairdel(unsigned char *x, int *n, int *s, double *d,
 423                          int *variance, double *var)
 424 {
 425     int i1, i2, Nd, Nv1, Nv2, L, s1, s2, target;
 426     double P, Q, R, a1, a2, a3, a, b, c;
 427
 428     target = 0;
 429     for (i1 = 1; i1 < *n; i1++) {
 430         for (i2 = i1 + 1; i2 <= *n; i2++) {
 431             Nd = Nv1 = Nv2 = L = 0;
 432             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 433                 CHECK_PAIRWISE_DELETION
 434                 COUNT_TS_TV1_TV2
 435             }
 436             COMPUTE_DIST_K81
 437             target++;
 438         }
 439     }
 440 }
 441
 442 #define PREPARE_BF_F84\
 443     A = (BF[0]*BF[2])/(BF[0] + BF[2]) + (BF[1]*BF[3])/(BF[1] + BF[3]);\
 444     B = BF[0]*BF[2] + BF[1]*BF[3];\
 445     C = (BF[0] + BF[2])*(BF[1] + BF[3]);
 446
 447 #define COMPUTE_DIST_F84\
 448    P = ((double) Ns/L);\
 449    Q = ((double) (Nd - Ns)/L);\
 450    d[target] = -2*A*log(1 - P/(2*A) - (A - B)*Q/(2*A*C)) + 2*(A - B - C)*log(1 - Q/(2*C));\
 451    if (*variance) {\
 452        t1 = A*C;\
 453        t2 = C*P/2;\
 454        t3 = (A - B)*Q/2;\
 455        a = t1/(t1 - t2 - t3);\
 456        b = A*(A - B)/(t1 - t2 - t3) - (A - B - C)/(C - Q/2);\
 457        var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/L;\
 458    }
 459
 460 void distDNA_F84(unsigned char *x, int *n, int *s, double *d,
 461                  double *BF, int *variance, double *var)
 462 {
 463     int i1, i2, Nd, Ns, L, target, s1, s2;
 464     double P, Q, A, B, C, a, b, t1, t2, t3;
 465
 466     PREPARE_BF_F84
 467     L = *s;
 468
 469     target = 0;
 470     for (i1 = 1; i1 < *n; i1++) {
 471         for (i2 = i1 + 1; i2 <= *n; i2++) {
 472             Nd = Ns = 0;
 473             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 474                 COUNT_TS_TV
 475             }
 476             COMPUTE_DIST_F84
 477             target++;
 478         }
 479     }
 480 }
 481
 482 void distDNA_F84_pairdel(unsigned char *x, int *n, int *s, double *d,
 483                          double *BF, int *variance, double *var)
 484 {
 485     int i1, i2, Nd, Ns, L, target, s1, s2;
 486     double P, Q, A, B, C, a, b, t1, t2, t3;
 487
 488     PREPARE_BF_F84
 489
 490     target = 0;
 491     for (i1 = 1; i1 < *n; i1++) {
 492         for (i2 = i1 + 1; i2 <= *n; i2++) {
 493             Nd = Ns = L = 0;
 494             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 495                 CHECK_PAIRWISE_DELETION
 496                 COUNT_TS_TV
 497             }
 498             COMPUTE_DIST_F84
 499             target++;
 500         }
 501     }
 502 }
 503
 504 #define COMPUTE_DIST_T92\
 505     P = ((double) Ns/L);\
 506     Q = ((double) (Nd - Ns)/L);\
 507     a1 = 1 - P/wg - Q;\
 508     a2 = 1 - 2*Q;\
 509     d[target] = -wg*log(a1) - 0.5*(1 - wg)*log(a2);\
 510     if (*variance) {\
 511         c1 = 1/a1;\
 512         c2 = 1/a2;\
 513         c3 = wg*(c1 - c2) + c2;\
 514         var[target] = (c1*c1*P + c3*c3*Q - pow(c1*P + c3*Q, 2))/L;\
 515     }
 516
 517 void distDNA_T92(unsigned char *x, int *n, int *s, double *d,
 518                  double *BF, int *variance, double *var)
 519 {
 520     int i1, i2, Nd, Ns, L, target, s1, s2;
 521     double P, Q, wg, a1, a2, c1, c2, c3;
 522
 523     L = *s;
 524     wg = 2 * (BF[1] + BF[2]) * (1 - (BF[1] + BF[2]));
 525
 526     target = 0;
 527     for (i1 = 1; i1 < *n; i1++) {
 528         for (i2 = i1 + 1; i2 <= *n; i2++) {
 529             Nd = Ns = 0;
 530             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 531                 COUNT_TS_TV
 532             }
 533             COMPUTE_DIST_T92
 534             target++;
 535         }
 536     }
 537 }
 538
 539 void distDNA_T92_pairdel(unsigned char *x, int *n, int *s, double *d,
 540                          double *BF, int *variance, double *var)
 541 {
 542     int i1, i2, Nd, Ns, L, target, s1, s2;
 543     double P, Q, wg, a1, a2, c1, c2, c3;
 544
 545     wg = 2 * (BF[1] + BF[2]) * (1 - (BF[1] + BF[2]));
 546
 547     target = 0;
 548     for (i1 = 1; i1 < *n; i1++) {
 549         for (i2 = i1 + 1; i2 <= *n; i2++) {
 550             Nd = Ns = L = 0;
 551             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 552                 CHECK_PAIRWISE_DELETION
 553                 COUNT_TS_TV
 554             }
 555             COMPUTE_DIST_T92
 556             target++;
 557         }
 558     }
 559 }
 560
 561 /* returns 1 if one of the base is adenine and
 562    the other one is guanine surely, 0 otherwise */
 563 #define AdenineAndGuanine(a, b) (a | b) == 200
 564
 565 /* returns 1 if one of the base is cytosine and
 566    the other one is thymine surely, 0 otherwise */
 567 #define CytosineAndThymine(a, b) (a | b) == 56
 568
 569 #define PREPARE_BF_TN93\
 570     gR = BF[0] + BF[2];\
 571     gY = BF[1] + BF[3];\
 572     k1 = 2 * BF[0] * BF[2] / gR;\
 573     k2 = 2 * BF[1] * BF[3] / gY;\
 574     k3 = 2 * (gR * gY - BF[0]*BF[2]*gY/gR - BF[1]*BF[3]*gR/gY);
 575
 576 #define COUNT_TS1_TS2_TV\
 577     if (DifferentBase(x[s1], x[s2])) {\
 578         Nd++;\
 579         if (AdenineAndGuanine(x[s1], x[s2])) {\
 580             Ns1++;\
 581         continue;\
 582         }\
 583         if (CytosineAndThymine(x[s1], x[s2])) Ns2++;\
 584     }
 585
 586 #define COMPUTE_DIST_TN93\
 587     P1 = ((double) Ns1/L);\
 588     P2 = ((double) Ns2/L);\
 589     Q = ((double) (Nd - Ns1 - Ns2)/L);\
 590     w1 = 1 - P1/k1 - Q/(2*gR);\
 591     w2 = 1 - P2/k2 - Q/(2*gY);\
 592     w3 = 1 - Q/(2*gR*gY);\
 593     if (*gamma) {\
 594         k4 = 2*(BF[0]*BF[2] + BF[1]*BF[3] + gR*gY);\
 595         b = -1 / *alpha;\
 596         c1 = pow(w1, b);\
 597         c2 = pow(w2, b);\
 598         c3 = pow(w3, b);\
 599         c4 = k1*c1/(2*gR) + k2*c2/(2*gY) + k3*c3/(2*gR*gY);\
 600         d[target] = *alpha * (k1*pow(w1, b) + k2*pow(w2, b) + k3*pow(w3, b) - k4);\
 601     } else {\
 602         k4 = 2*((BF[0]*BF[0] + BF[2]*BF[2])/(2*gR*gR) + (BF[2]*BF[2] + BF[3]*BF[3])/(2*gY*gY));\
 603         c1 = 1/w1;\
 604         c2 = 1/w2;\
 605         c3 = 1/w3;\
 606         c4 = k1 * c1/(2 * gR) + k2 * c2/(2 * gY) + k4 * c3;\
 607         d[target] = -k1*log(w1) - k2*log(w2) - k3*log(w3);\
 608     }\
 609     if (*variance)\
 610       var[target] = (c1*c1*P1 + c2*c2*P2 + c4*c4*Q - pow(c1*P1 + c2*P2 + c4*Q, 2))/L;
 611
 612 void distDNA_TN93(unsigned char *x, int *n, int *s, double *d,
 613                   double *BF, int *variance, double *var,
 614                   int *gamma, double *alpha)
 615 {
 616     int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
 617     double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 618
 619     L = *s;
 620
 621     PREPARE_BF_TN93
 622
 623     target = 0;
 624     for (i1 = 1; i1 < *n; i1++) {
 625         for (i2 = i1 + 1; i2 <= *n; i2++) {
 626             Nd = Ns1 = Ns2 = 0;
 627             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 628                 COUNT_TS1_TS2_TV
 629             }
 630             COMPUTE_DIST_TN93
 631             target++;
 632         }
 633     }
 634 }
 635
 636 void distDNA_TN93_pairdel(unsigned char *x, int *n, int *s, double *d,
 637                           double *BF, int *variance, double *var,
 638                           int *gamma, double *alpha)
 639 {
 640     int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
 641     double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 642
 643     PREPARE_BF_TN93
 644
 645     target = 0;
 646     for (i1 = 1; i1 < *n; i1++) {
 647         for (i2 = i1 + 1; i2 <= *n; i2++) {
 648             Nd = Ns1 = Ns2 = L = 0;
 649             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 650                 CHECK_PAIRWISE_DELETION
 651                 COUNT_TS1_TS2_TV
 652             }
 653             COMPUTE_DIST_TN93
 654             target++;
 655         }
 656     }
 657 }
 658
 659 void distDNA_GG95(unsigned char *x, int *n, int *s, double *d,
 660                   int *variance, double *var)
 661 {
 662     int i1, i2, s1, s2, target, GC, Nd, Ns, tl, npair;
 663     double *theta, gcprop, *P, pp, *Q, qq, *tstvr, svr, A, sum, ma /* mean alpha */, K1, K2;
 664
 665     theta = &gcprop;
 666     P = &pp;
 667     Q = &qq;
 668     tstvr = &svr;
 669
 670     npair = *n * (*n - 1) / 2;
 671
 672     theta = (double*)R_alloc(*n, sizeof(double));
 673     P = (double*)R_alloc(npair, sizeof(double));
 674     Q = (double*)R_alloc(npair, sizeof(double));
 675     tstvr = (double*)R_alloc(npair, sizeof(double));
 676
 677     /* get the proportion of GC (= theta) in each sequence */
 678     for (i1 = 1; i1 <= *n; i1++) {
 679         GC = 0;
 680         for (s1 = i1 - 1; s1 < i1 + *n*(*s - 1); s1 += *n)
 681           if (IsCytosine(x[s1]) || IsGuanine(x[s1])) GC += 1;
 682         theta[i1 - 1] = ((double) GC / *s);
 683     }
 684
 685     /* get the proportions of transitions and transversions,
 686        and the estimates of their ratio for each pair */
 687     target = 0;
 688     for (i1 = 1; i1 < *n; i1++) {
 689         for (i2 = i1 + 1; i2 <= *n; i2++) {
 690             Nd = Ns = 0;
 691             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 692                 COUNT_TS_TV
 693             }
 694             P[target] = ((double) Ns / *s);
 695             Q[target] = ((double) (Nd - Ns) / *s);
 696             A = log(1 - 2*Q[target]);
 697             tstvr[target] = 2*(log(1 - 2*P[target] - Q[target]) - 0.5*A)/A;
 698             target++;
 699         }
 700     }
 701
 702     /* compute the mean alpha (ma) = mean Ts/Tv */
 703     sum = 0;
 704     tl = 0;
 705     for (i1 = 0; i1 < npair; i1++)
 706     /* some values of tstvr are -Inf if there is no
 707        transversions observed: we exclude them */
 708       if (R_FINITE(tstvr[i1])) {
 709           sum += tstvr[i1];
 710           tl += 1;
 711       }
 712     ma = sum/tl;
 713
 714     /* compute the distance for each pair */
 715     target = 0;
 716     for (i1 = 1; i1 < *n; i1++) {
 717         for (i2 = i1 + 1; i2 <= *n; i2++) {
 718             A = 1 - 2*Q[target];
 719             K1 = 1 + ma*(theta[i1 - 1]*(1 - theta[i1 - 1]) + theta[i2 - 1]*(1 - theta[i2 - 1]));
 720             K2 = ma*pow(theta[i1 - 1] - theta[i2 - 1], 2)/(ma + 1);
 721             d[target] = -0.5*K1*log(A) + K2*(1 - pow(A, 0.25*(ma + 1)));
 722             if (*variance)
 723               var[target] = pow(K1 + K2*0.5*(ma + 1)*pow(A, 0.25*(ma + 1)), 2)*Q[target]*(1 - Q[target])/(A*A * *s);
 724             target++;
 725         }
 726     }
 727 }
 728
 729 void distDNA_GG95_pairdel(unsigned char *x, int *n, int *s, double *d,
 730                           int *variance, double *var)
 731 {
 732     int i1, i2, s1, s2, target, *L, length, GC, Nd, Ns, tl, npair;
 733     double *theta, gcprop, *P, pp, *Q, qq, *tstvr, svr, A, sum, ma /* mean alpha */, K1, K2;
 734
 735     theta = &gcprop;
 736     L = &length;
 737     P = &pp;
 738     Q = &qq;
 739     tstvr = &svr;
 740
 741     npair = *n * (*n - 1) / 2;
 742
 743     theta = (double*)R_alloc(*n, sizeof(double));
 744     L = (int*)R_alloc(npair, sizeof(int));
 745     P = (double*)R_alloc(npair, sizeof(double));
 746     Q = (double*)R_alloc(npair, sizeof(double));
 747     tstvr = (double*)R_alloc(npair, sizeof(double));
 748
 749     /* get the proportion of GC (= theta) in each sequence */
 750     for (i1 = 1; i1 <= *n; i1++) {
 751         tl = GC = 0;
 752         for (s1 = i1 - 1; s1 < i1 + *n*(*s - 1); s1 += *n) {
 753             if (KnownBase(x[s1])) tl++;
 754             else continue;
 755             if (IsCytosine(x[s1]) || IsGuanine(x[s1])) GC += 1;
 756         }
 757         theta[i1 - 1] = ((double) GC / tl);
 758     }
 759
 760     /* get the proportions of transitions and transversions,
 761        and the estimates of their ratio for each pair; we
 762        also get the sample size for each pair in L */
 763     target = 0;
 764     for (i1 = 1; i1 < *n; i1++) {
 765         for (i2 = i1 + 1; i2 <= *n; i2++) {
 766             Nd = Ns = L[target] = 0;
 767             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 768                 if (KnownBase(x[s1]) && KnownBase(x[s2])) L[target]++;
 769                 else continue;
 770                 COUNT_TS_TV
 771             }
 772             P[target] = ((double) Ns/L[target]);
 773             Q[target] = ((double) (Nd - Ns)/L[target]);
 774             A = log(1 - 2*Q[target]);
 775             tstvr[target] = 2*(log(1 - 2*P[target] - Q[target]) - 0.5*A)/A;
 776             target++;
 777         }
 778     }
 779
 780     /* compute the mean alpha (ma) = mean Ts/Tv */
 781     sum = 0;
 782     tl = 0;
 783     for (i1 = 0; i1 < npair; i1++)
 784     /* some values of tstvr are -Inf if there is no
 785        transversions observed: we exclude them */
 786       if (R_FINITE(tstvr[i1])) {
 787           sum += tstvr[i1];
 788           tl += 1;
 789       }
 790     ma = sum/tl;
 791
 792     /* compute the distance for each pair */
 793     target = 0;
 794     for (i1 = 1; i1 < *n; i1++) {
 795         for (i2 = i1 + 1; i2 <= *n; i2++) {
 796             A = 1 - 2*Q[target];
 797             K1 = 1 + ma*(theta[i1 - 1]*(1 - theta[i1 - 1]) + theta[i2 - 1]*(1 - theta[i2 - 1]));
 798             K2 = ma*pow(theta[i1 - 1] - theta[i2 - 1], 2)/(ma + 1);
 799             d[target] = -0.5*K1*log(A) + K2*(1 - pow(A, 0.25*(ma + 1)));
 800             if (*variance)
 801               var[target] = pow(K1 + K2*0.5*(ma + 1)*pow(A, 0.25*(ma + 1)), 2)*Q[target]*(1 - Q[target])/(A*A*L[target]);
 802             target++;
 803         }
 804     }
 805 }
 806
 807 #define DO_CONTINGENCY_NUCLEOTIDES\
 808     switch (x[s1]) {\
 809     case 136 : m = 0; break;\
 810     case 72 : m = 1; break;\
 811     case 40 : m = 2; break;\
 812     case 24 : m = 3; break;\
 813     }\
 814     switch (x[s2]) {\
 815     case 72 : m += 4; break;\
 816     case 40 : m += 8; break;\
 817     case 24 : m += 12; break;\
 818     }\
 819     Ntab[m]++;
 820
 821 #define COMPUTE_DIST_LogDet\
 822     for (k = 0; k < 16; k++) Ftab[k] = ((double) Ntab[k]/L);\
 823     d[target] = -log(detFourByFour(Ftab))/4 - LN4;\
 824     if (*variance) {\
 825         /* For the inversion, we first make U an identity matrix */\
 826         for (k = 1; k < 15; k++) U[k] = 0;\
 827         U[0] = U[5] = U[10] = U[15] = 1;\
 828         /* The matrix is not symmetric, so we use 'dgesv'. */\
 829         /* This subroutine puts the result in U. */\
 830         F77_CALL(dgesv)(&ndim, &ndim, Ftab, &ndim, ipiv, U, &ndim, &info);\
 831         var[target] = (U[0]*U[0]*Ftab[0] + U[1]*U[1]*Ftab[4] +\
 832                        U[2]*U[2]*Ftab[8] + U[3]*U[3]*Ftab[12] +\
 833                        U[4]*U[4]*Ftab[1] + U[5]*U[5]*Ftab[5] +\
 834                        U[6]*U[6]*Ftab[9] + U[7]*U[7]*Ftab[13] +\
 835                        U[8]*U[8]*Ftab[2] + U[9]*U[9]*Ftab[6] +\
 836                        U[10]*U[10]*Ftab[10] + U[11]*U[11]*Ftab[14] +\
 837                        U[12]*U[12]*Ftab[3] + U[13]*U[13]*Ftab[7] +\
 838                        U[14]*U[14]*Ftab[11] + U[15]*U[15]*Ftab[15] - 16)/(L*16);\
 839     }
 840
 841 void distDNA_LogDet(unsigned char *x, int *n, int *s, double *d,
 842                     int *variance, double *var)
 843 {
 844     int i1, i2, k, m, s1, s2, target, L, Ntab[16], ndim = 4, info, ipiv[16];
 845     double Ftab[16], U[16];
 846
 847     L = *s;
 848
 849     target = 0;
 850     for (i1 = 1; i1 < *n; i1++) {
 851         for (i2 = i1 + 1; i2 <= *n; i2++) {
 852             for (k = 0; k < 16; k++) Ntab[k] = 0;
 853             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 854                 DO_CONTINGENCY_NUCLEOTIDES
 855             }
 856             COMPUTE_DIST_LogDet
 857             target++;
 858         }
 859     }
 860 }
 861
 862 void distDNA_LogDet_pairdel(unsigned char *x, int *n, int *s, double *d,
 863                             int *variance, double *var)
 864 {
 865     int i1, i2, k, m, s1, s2, target, L, Ntab[16], ndim = 4, info, ipiv[16];
 866     double Ftab[16], U[16];
 867
 868     target = 0;
 869     for (i1 = 1; i1 < *n; i1++) {
 870         for (i2 = i1 + 1; i2 <= *n; i2++) {
 871             for (k = 0; k < 16; k++) Ntab[k] = 0;
 872             L = 0;
 873             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 874                 CHECK_PAIRWISE_DELETION
 875                 DO_CONTINGENCY_NUCLEOTIDES
 876             }
 877             COMPUTE_DIST_LogDet
 878             target++;
 879         }
 880     }
 881 }
 882
 883 void distDNA_BH87(unsigned char *x, int *n, int *s, double *d,
 884                   int *variance, double *var)
 885 /* <FIXME>
 886    For the moment there is no need to check for pairwise deletions
 887    since DO_CONTINGENCY_NUCLEOTIDES considers only the known nucleotides.
 888    In effect the pairwise deletion has possibly been done before.
 889    The sequence length(s) are used only to compute the variances, which is
 890    currently not available.
 891    </FIXME> */
 892 {
 893     int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4];
 894     double P12[16], P21[16];
 895
 896     for (i1 = 1; i1 < *n; i1++) {
 897         for (i2 = i1 + 1; i2 <= *n; i2++) {
 898             for (k = 0; k < 16; k++) Ntab[k] = 0;
 899             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 900                 DO_CONTINGENCY_NUCLEOTIDES
 901             }
 902
 903             /* get the rowwise sums of Ntab */
 904             ROWsums[0] = Ntab[0] + Ntab[4] + Ntab[8] + Ntab[12];
 905             ROWsums[1] = Ntab[1] + Ntab[5] + Ntab[9] + Ntab[13];
 906             ROWsums[2] = Ntab[2] + Ntab[6] + Ntab[10] + Ntab[14];
 907             ROWsums[3] = Ntab[3] + Ntab[7] + Ntab[11] + Ntab[15];
 908
 909             for (k = 0; k < 16; k++)
 910               P12[k] = ((double) Ntab[k]);
 911
 912             /* scale each element of P12 by its rowwise sum */
 913             for (k = 0; k < 4; k++)
 914               for (kb = 0; kb < 16; kb += 4)
 915                 P12[k + kb] = P12[k + kb]/ROWsums[k];
 916
 917             d[*n*(i2 - 1) + i1 - 1] = -log(detFourByFour(P12))/4;
 918
 919             /* compute the columnwise sums of Ntab: these
 920                are the rowwise sums of its transpose */
 921             ROWsums[0] = Ntab[0] + Ntab[1] + Ntab[2] + Ntab[3];
 922             ROWsums[1] = Ntab[4] + Ntab[5] + Ntab[6] + Ntab[7];
 923             ROWsums[2] = Ntab[8] + Ntab[9] + Ntab[10] + Ntab[11];
 924             ROWsums[3] = Ntab[12] + Ntab[13] + Ntab[14] + Ntab[15];
 925
 926             /* transpose Ntab and store the result in P21 */
 927             for (k = 0; k < 4; k++)
 928                for (kb = 0; kb < 4; kb++)
 929                  P21[kb + 4*k] = Ntab[k + 4*kb];
 930
 931             /* scale as above */
 932             for (k = 0; k < 4; k++)
 933               for (kb = 0; kb < 16; kb += 4)
 934                 P21[k + kb] = P21[k + kb]/ROWsums[k];
 935
 936             d[*n*(i1 - 1) + i2 - 1] = -log(detFourByFour(P21))/4;
 937         }
 938     }
 939 }
 940
 941 #define COMPUTE_DIST_ParaLin\
 942     for (k = 0; k < 16; k++) Ftab[k] = ((double) Ntab[k]/L);\
 943     d[target] = -log(detFourByFour(Ftab)/\
 944                      sqrt(find[0][i1 - 1]*find[1][i1 - 1]*find[2][i1 - 1]*find[3][i1 - 1]*\
 945                           find[0][i2 - 1]*find[1][i2 - 1]*find[2][i2 - 1]*find[3][i2 - 1]))/4;\
 946     if (*variance) {\
 947         /* For the inversion, we first make U an identity matrix */\
 948         for (k = 1; k < 15; k++) U[k] = 0;\
 949         U[0] = U[5] = U[10] = U[15] = 1;\
 950         /* The matrix is not symmetric, so we use 'dgesv'. */\
 951         /* This subroutine puts the result in U. */\
 952         F77_CALL(dgesv)(&ndim, &ndim, Ftab, &ndim, ipiv, U, &ndim, &info);\
 953         var[target] = (U[0]*U[0]*Ftab[0] + U[1]*U[1]*Ftab[4] +\
 954                        U[2]*U[2]*Ftab[8] + U[3]*U[3]*Ftab[12] +\
 955                        U[4]*U[4]*Ftab[1] + U[5]*U[5]*Ftab[5] +\
 956                        U[6]*U[6]*Ftab[9] + U[7]*U[7]*Ftab[13] +\
 957                        U[8]*U[8]*Ftab[2] + U[9]*U[9]*Ftab[6] +\
 958                        U[10]*U[10]*Ftab[10] + U[11]*U[11]*Ftab[14] +\
 959                        U[12]*U[12]*Ftab[3] + U[13]*U[13]*Ftab[7] +\
 960                        U[14]*U[14]*Ftab[11] + U[15]*U[15]*Ftab[15] -\
 961                        4*(1/sqrt(find[0][i1 - 1]*find[0][i2 - 1]) +\
 962                        1/sqrt(find[1][i1 - 1]*find[1][i2 - 1]) +\
 963                        1/sqrt(find[2][i1 - 1]*find[2][i2 - 1]) +\
 964                        1/sqrt(find[3][i1 - 1]*find[3][i2 - 1])))/(L*16);\
 965     }
 966
 967 void distDNA_ParaLin(unsigned char *x, int *n, int *s, double *d,
 968                      int *variance, double *var)
 969 {
 970     int i1, i2, k, s1, s2, m, target, L, Ntab[16], ndim = 4, info, ipiv[16];
 971     double Ftab[16], U[16], *find[4];
 972
 973     L = *s;
 974
 975     for (k = 0; k < 4; k++)
 976       find[k] = (double*)R_alloc(*n, sizeof(double));
 977
 978     for (i1 = 0; i1 < *n; i1++)
 979       for (k = 0; k < 4; k++) find[k][i1] = 0.0;
 980
 981     for (i1 = 0; i1 < *n; i1++) {
 982         for (s1 = i1; s1 < i1 + *n*(*s - 1) + 1; s1+= *n) {
 983             switch (x[s1]) {
 984             case 136 : find[0][i1]++; break;
 985             case 40 : find[1][i1]++; break;
 986             case 72 : find[2][i1]++; break;
 987             case 24 : find[3][i1]++; break;
 988             }
 989         }
 990         for (k = 0; k < 4; k++) find[k][i1] /= L;
 991     }
 992
 993     target = 0;
 994     for (i1 = 1; i1 < *n; i1++) {
 995         for (i2 = i1 + 1; i2 <= *n; i2++) {
 996             for (k = 0; k < 16; k++) Ntab[k] = 0;
 997             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 998                 DO_CONTINGENCY_NUCLEOTIDES
 999             }
1000             COMPUTE_DIST_ParaLin
1001             target++;
1002         }
1003     }
1004 }
1005
1006 void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
1007                              int *variance, double *var)
1008 {
1009     int i1, i2, k, s1, s2, m, target, L, Ntab[16], ndim = 4, info, ipiv[16];
1010     double Ftab[16], U[16], *find[4];
1011
1012     L = 0;
1013
1014     for (k = 0; k < 4; k++)
1015       find[k] = (double*)R_alloc(*n, sizeof(double));
1016
1017     for (i1 = 0; i1 < *n; i1++)
1018       for (k = 0; k < 4; k++) find[k][i1] = 0.0;
1019
1020     for (i1 = 0; i1 < *n; i1++) {
1021         L = 0;
1022         for (s1 = i1; s1 < i1 + *n*(*s - 1) + 1; s1+= *n) {
1023             if (KnownBase(x[s1])) {
1024                 L++;
1025                 switch (x[s1]) {
1026                 case 136 : find[0][i1]++; break;
1027                 case 40 : find[1][i1]++; break;
1028                 case 72 : find[2][i1]++; break;
1029                 case 24 : find[3][i1]++; break;
1030                 }
1031             }
1032         }
1033         for (k = 0; k < 4; k++) find[k][i1] /= L;
1034     }
1035
1036     target = 0;
1037     for (i1 = 1; i1 < *n; i1++) {
1038         for (i2 = i1 + 1; i2 <= *n; i2++) {
1039             L = 0;
1040             for (k = 0; k < 16; k++) Ntab[k] = 0;
1041             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
1042                 CHECK_PAIRWISE_DELETION
1043                 DO_CONTINGENCY_NUCLEOTIDES
1044             }
1045             COMPUTE_DIST_ParaLin
1046             target++;
1047         }
1048     }
1049 }
1050
1051 /* void BaseProportion(unsigned char *x, int *n, double *BF, int *freq) */
1052 /* { */
1053 /*     int i, m; */
1054
1055 /*     m = 0; */
1056 /*     for (i = 0; i < *n; i++) { */
1057 /*         if (KnownBase(x[i])) { */
1058 /*          m++; */
1059 /*          switch (x[i]) { */
1060 /*          case 136 : BF[0]++; break; */
1061 /*          case 40 : BF[1]++; break; */
1062 /*          case 72 : BF[2]++; break; */
1063 /*          case 24 : BF[3]++; break; */
1064 /*          } */
1065 /*      } */
1066 /*     } */
1067 /*     if (! *freq) for (i = 0; i < 4; i++) BF[i] /= m; */
1068 /* } */
1069
1070 /* void BaseProportion(unsigned char *x, int *n, double *BF) */
1071 /* { */
1072 /*     int i; */
1073
1074 /*     for (i = 0; i < *n; i++) { */
1075 /*          switch (x[i]) { */
1076 /*          case 136 : BF[0]++; break; */
1077 /*          case 40 : BF[1]++; break; */
1078 /*          case 72 : BF[2]++; break; */
1079 /*          case 24 : BF[3]++; break; */
1080 /*          case 192 : BF[4]++; break; */
1081 /*          case 160 : BF[5]++; break; */
1082 /*          case 144 : BF[6]++; break; */
1083 /*          case 96 : BF[7]++; break; */
1084 /*          case 80 : BF[8]++; break; */
1085 /*          case 48 : BF[9]++; break; */
1086 /*          case 224 : BF[10]++; break; */
1087 /*          case 176 : BF[11]++; break; */
1088 /*          case 208 : BF[12]++; break; */
1089 /*          case 112 : BF[13]++; break; */
1090 /*          case 240 : BF[14]++; break; */
1091 /*          case 4 : BF[15]++; break; */
1092 /*          case 2 : BF[16]++; break; */
1093 /*          } */
1094 /*     } */
1095 /* } */
1096
1097 /* a hash table is much faster than switch (2012-01-10) */
1098 void BaseProportion(unsigned char *x, int *n, double *BF)
1099 {
1100         int i;
1101         double count[256];
1102
1103         memset(count, 0, 256*sizeof(double));
1104
1105         for (i = 0; i < *n; i++) count[x[i]]++;
1106
1107         BF[0] = count[136];
1108         BF[1] = count[40];
1109         BF[2] = count[72];
1110         BF[3] = count[24];
1111         BF[4] = count[192];
1112         BF[5] = count[160];
1113         BF[6] = count[144];
1114         BF[7] = count[96];
1115         BF[8] = count[80];
1116         BF[9] = count[48];
1117         BF[10] = count[224];
1118         BF[11] = count[176];
1119         BF[12] = count[208];
1120         BF[13] = count[112];
1121         BF[14] = count[240];
1122         BF[15] = count[4];
1123         BF[16] = count[2];
1124 }
1125
1126 void SegSites(unsigned char *x, int *n, int *s, int *seg)
1127 {
1128     int i, j;
1129     unsigned char basis;
1130
1131     for (j = 0; j < *s; j++) {
1132         i = *n * j;
1133         while (!KnownBase(x[i])) i++;
1134         basis = x[i];
1135         i++;
1136         while (i < *n * (j + 1)) {
1137             if (!KnownBase(x[i]) || x[i] == basis) i++;
1138             else {
1139                 seg[j] = 1;
1140                 break;
1141             }
1142         }
1143     }
1144 }
1145
1146 void GlobalDeletionDNA(unsigned char *x, int *n, int *s, int *keep)
1147 {
1148     int i, j;
1149
1150     for (j = 0; j < *s; j++) {
1151         i = *n * j;
1152         while (i < *n * (j + 1)) {
1153             if (KnownBase(x[i])) i++;
1154             else {
1155                 keep[j] = 0;
1156                 break;
1157             }
1158         }
1159     }
1160 }
1161
1162 void dist_dna(unsigned char *x, int *n, int *s, int *model, double *d,
1163               double *BF, int *pairdel, int *variance, double *var,
1164               int *gamma, double *alpha)
1165 {
1166     switch (*model) {
1167     case 1 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 1);
1168              else distDNA_raw(x, n, s, d, 1); break;
1169     case 2 : if (pairdel) distDNA_JC69_pairdel(x, n, s, d, variance, var, gamma, alpha);
1170              else distDNA_JC69(x, n, s, d, variance, var, gamma, alpha); break;
1171     case 3 : if (pairdel) distDNA_K80_pairdel(x, n, s, d, variance, var, gamma, alpha);
1172              else distDNA_K80(x, n, s, d, variance, var, gamma, alpha); break;
1173     case 4 : if (pairdel) distDNA_F81_pairdel(x, n, s, d, BF, variance, var, gamma, alpha);
1174              else distDNA_F81(x, n, s, d, BF, variance, var, gamma, alpha); break;
1175     case 5 : if (pairdel) distDNA_K81_pairdel(x, n, s, d, variance, var);
1176              else distDNA_K81(x, n, s, d, variance, var); break;
1177     case 6 : if (pairdel) distDNA_F84_pairdel(x, n, s, d, BF, variance, var);
1178              else distDNA_F84(x, n, s, d, BF, variance, var); break;
1179     case 7 : if (pairdel) distDNA_T92_pairdel(x, n, s, d, BF, variance, var);
1180              else distDNA_T92(x, n, s, d, BF, variance, var); break;
1181     case 8 : if (pairdel) distDNA_TN93_pairdel(x, n, s, d, BF, variance, var, gamma, alpha);
1182              else distDNA_TN93(x, n, s, d, BF, variance, var, gamma, alpha); break;
1183     case 9 : if (pairdel) distDNA_GG95_pairdel(x, n, s, d, variance, var);
1184              else distDNA_GG95(x, n, s, d, variance, var); break;
1185     case 10 : if (pairdel) distDNA_LogDet_pairdel(x, n, s, d, variance, var);
1186               else distDNA_LogDet(x, n, s, d, variance, var); break;
1187     case 11 : distDNA_BH87(x, n, s, d, variance, var); break;
1188     case 12 : if (pairdel) distDNA_ParaLin_pairdel(x, n, s, d, variance, var);
1189               else distDNA_ParaLin(x, n, s, d, variance, var); break;
1190     case 13 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 0);
1191              else distDNA_raw(x, n, s, d, 0); break;
1192     case 14 : if (pairdel) distDNA_TsTv(x, n, s, d, 1, 1);
1193             else distDNA_TsTv(x, n, s, d, 1, 0); break;
1194     case 15 : if (pairdel) distDNA_TsTv(x, n, s, d, 0, 1);
1195             else distDNA_TsTv(x, n, s, d, 0, 0); break;
1196     case 16 : distDNA_indel(x, n, s, d); break;
1197     case 17 : distDNA_indelblock(x, n, s, d); break;
1198     }
1199 }