src/dist_dna.c

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