src/dist_dna.c

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