src/dist_dna.c

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