X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2Fdist_dna.c;h=c087212dc1b7477435a1fcf6b69630c6811a3466;hb=c488b74490ee3d9d200de0e471881f002a18fe4f;hp=4e02393334a11ccd8d749294e19cff9d6a41f9f1;hpb=c827059eeafc8cbe41c812b26979543ab287803e;p=ape.git

diff --git a/src/dist_dna.c b/src/dist_dna.c
index 4e02393..c087212 100644
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,18 +1,18 @@
-/* dist_dna.c       2007-12-01 */
+/* dist_dna.c       2012-11-28 */
 
-/* Copyright 2005-2007 Emmanuel Paradis
+/* Copyright 2005-2012 Emmanuel Paradis
 
 /* This file is part of the R-package `ape'. */
 /* See the file ../COPYING for licensing issues. */
 
-#include <R.h>
 #include <R_ext/Lapack.h>
+#include "ape.h"
 
 /* from R: print(log(4), d = 22) */
 #define LN4 1.386294361119890572454
 
 /* returns 8 if the base is known surely, 0 otherwise */
-#define KnownBase(a) a & 8
+#define KnownBase(a) (a & 8)
 
 /* returns 1 if the base is adenine surely, 0 otherwise */
 #define IsAdenine(a) a == 136
@@ -64,23 +64,122 @@ double detFourByFour(double *x)
     if (KnownBase(x[s1]) && KnownBase(x[s2])) L++;\
     else continue;
 
-void distDNA_raw(unsigned char *x, int *n, int *s, double *d)
+#define COUNT_TS_TV\
+    if (SameBase(x[s1], x[s2])) continue;\
+    Nd++;\
+    if (IsPurine(x[s1]) && IsPurine(x[s2])) {\
+        Ns++;\
+        continue;\
+    }\
+    if (IsPyrimidine(x[s1]) && IsPyrimidine(x[s2])) Ns++;
+
+void distDNA_indel(unsigned char *x, int *n, int *s, double *d)
 {
-    int i1, i2, s1, s2, target, Nd;
+	int i1, i2, s1, s2, target, N;
 
-    target = 0;
-    for (i1 = 1; i1 < *n; i1++) {
-        for (i2 = i1 + 1; i2 <= *n; i2++) {
-	    Nd = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n)
-	      if (DifferentBase(x[s1], x[s2])) Nd++;
-	    d[target] = ((double) Nd / *s);
-	    target++;
+	target = 0;
+	for (i1 = 1; i1 < *n; i1++) {
+		for (i2 = i1 + 1; i2 <= *n; i2++) {
+			N = 0;
+
+			for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
+				if ((x[s1] ^ x[s2]) & 4) N++;
+
+			d[target] = ((double) N);
+			target++;
+		}
+	}
+}
+
+#define X(i, j) i - 1 + *n * (j - 1)
+
+void distDNA_indelblock(unsigned char *x, int *n, int *s, double *d)
+{
+	int i1, i2, s1, s2, target, N, start_block, end_block;
+
+	for (i1 = 1; i1 <= *n; i1++) {
+
+/* find a block of one or more '-' */
+
+		s1 = 1;
+		while (s1 < *s) {
+			if (x[X(i1, s1)] == 4) {
+				start_block = s1;
+				while (x[X(i1, s1)] == 4) s1++;
+				end_block = s1 - 1;
+
+/* then look if the same block is present in all the other sequences */
+
+				for (i2 = 1; i2 <= *n; i2++) {
+					if (i2 == i1) continue;
+
+					target = give_index(i1, i2, *n);
+
+					if (start_block > 1) {
+						if (x[X(i2, start_block - 1)] == 4) {
+							d[target]++;
+							continue;
+						}
+					}
+					if (end_block < *s) {
+						if (x[X(i2, end_block + 1)] == 4) {
+							d[target]++;
+							continue;
+						}
+					}
+					for (s2 = start_block; s2 <= end_block; s2++) {
+						if (x[X(i2, s2)] != 4) {
+							d[target]++;
+							continue;
+						}
+					}
+				}
+				s1 = end_block + 1;
+			}
+			s1++;
+		}
 	}
-    }
 }
 
-void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d)
+#undef X
+
+void distDNA_TsTv(unsigned char *x, int *n, int *s, double *d, int Ts, int pairdel)
+{
+	int i1, i2, s1, s2, target, Nd, Ns;
+
+	target = 0;
+	for (i1 = 1; i1 < *n; i1++) {
+		for (i2 = i1 + 1; i2 <= *n; i2++) {
+			Nd = Ns = 0;
+			for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
+				if (pairdel && !(KnownBase(x[s1]) && KnownBase(x[s2]))) continue;
+				COUNT_TS_TV
+			}
+			if (Ts) d[target] = ((double) Ns); /* output number of transitions */
+			else d[target] = ((double) Nd - Ns); /* output number of transversions */
+			target++;
+		}
+	}
+}
+
+void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled)
+{
+	int i1, i2, s1, s2, target, Nd;
+
+	target = 0;
+	for (i1 = 1; i1 < *n; i1++) {
+		for (i2 = i1 + 1; i2 <= *n; i2++) {
+			Nd = 0;
+			for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
+				if (DifferentBase(x[s1], x[s2])) Nd++;
+			if (scaled) d[target] = ((double) Nd / *s);
+			else d[target] = ((double) Nd);
+			target++;
+		}
+	}
+}
+
+void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d, int scaled)
 {
     int i1, i2, s1, s2, target, Nd, L;
 
@@ -88,11 +187,12 @@ void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d)
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
                 CHECK_PAIRWISE_DELETION
 		if (DifferentBase(x[s1], x[s2])) Nd++;
 	    }
-	    d[target] = ((double) Nd/L);
+	    if (scaled) d[target] = ((double) Nd/L);
+	    else d[target] = ((double) Nd);
 	    target++;
 	}
     }
@@ -120,7 +220,7 @@ void distDNA_JC69(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
   	    Nd = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n)
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n)
 	      if (DifferentBase(x[s1], x[s2])) Nd++;
 	    COMPUTE_DIST_JC69
 	    target++;
@@ -148,15 +248,6 @@ void distDNA_JC69_pairdel(unsigned char *x, int *n, int *s, double *d,
     }
 }
 
-#define COUNT_TS_TV\
-    if (SameBase(x[s1], x[s2])) continue;\
-    Nd++;\
-    if (IsPurine(x[s1]) && IsPurine(x[s2])) {\
-        Ns++;\
-        continue;\
-    }\
-    if (IsPyrimidine(x[s1]) && IsPyrimidine(x[s2])) Ns++;
-
 #define COMPUTE_DIST_K80\
     P = ((double) Ns/L);\
     Q = ((double) (Nd - Ns)/L);\
@@ -212,7 +303,7 @@ void distDNA_K80_pairdel(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns = L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        CHECK_PAIRWISE_DELETION
 		COUNT_TS_TV
 	    }
@@ -264,7 +355,7 @@ void distDNA_F81_pairdel(unsigned char *x, int *n, int *s, double *d, double *BF
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
   	    Nd = L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        CHECK_PAIRWISE_DELETION
 		if (DifferentBase(x[s1], x[s2])) Nd++;
 	    }
@@ -311,7 +402,7 @@ void distDNA_K81(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
   	    Nd = Nv1 = Nv2 = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        COUNT_TS_TV1_TV2
 	    }
 	    COMPUTE_DIST_K81
@@ -330,7 +421,7 @@ void distDNA_K81_pairdel(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
   	    Nd = Nv1 = Nv2 = L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        CHECK_PAIRWISE_DELETION
 	        COUNT_TS_TV1_TV2
 	    }
@@ -348,14 +439,14 @@ void distDNA_K81_pairdel(unsigned char *x, int *n, int *s, double *d,
 #define COMPUTE_DIST_F84\
    P = ((double) Ns/L);\
    Q = ((double) (Nd - Ns)/L);\
-   d[target] = -2*A*log(1 - (P/(2*A) - (A - B)*Q/(2*A*C))) + 2*(A - B - C)*log(1 - Q/(2*C));\
+   d[target] = -2*A*log(1 - P/(2*A) - (A - B)*Q/(2*A*C)) + 2*(A - B - C)*log(1 - Q/(2*C));\
    if (*variance) {\
        t1 = A*C;\
        t2 = C*P/2;\
        t3 = (A - B)*Q/2;\
        a = t1/(t1 - t2 - t3);\
        b = A*(A - B)/(t1 - t2 - t3) - (A - B - C)/(C - Q/2);\
-       var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/2;\
+       var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/L;\
    }
 
 void distDNA_F84(unsigned char *x, int *n, int *s, double *d,
@@ -371,7 +462,7 @@ void distDNA_F84(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        COUNT_TS_TV
 	    }
 	    COMPUTE_DIST_F84
@@ -392,7 +483,7 @@ void distDNA_F84_pairdel(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns = L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 		CHECK_PAIRWISE_DELETION
 		COUNT_TS_TV
 	    }
@@ -428,7 +519,7 @@ void distDNA_T92(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        COUNT_TS_TV
 	    }
 	    COMPUTE_DIST_T92
@@ -449,7 +540,7 @@ void distDNA_T92_pairdel(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns = L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        CHECK_PAIRWISE_DELETION
 	        COUNT_TS_TV
 	    }
@@ -514,8 +605,8 @@ void distDNA_TN93(unsigned char *x, int *n, int *s, double *d,
 		  double *BF, int *variance, double *var,
 		  int *gamma, double *alpha)
 {
-    int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
-    double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
+    int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
+    double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 
     L = *s;
 
@@ -525,7 +616,7 @@ void distDNA_TN93(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns1 = Ns2 = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 		COUNT_TS1_TS2_TV
 	    }
 	    COMPUTE_DIST_TN93
@@ -538,8 +629,8 @@ void distDNA_TN93_pairdel(unsigned char *x, int *n, int *s, double *d,
 			  double *BF, int *variance, double *var,
 			  int *gamma, double *alpha)
 {
-    int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
-    double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
+    int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
+    double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 
     PREPARE_BF_TN93
 
@@ -547,7 +638,7 @@ void distDNA_TN93_pairdel(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns1 = Ns2 = L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 		CHECK_PAIRWISE_DELETION
 		COUNT_TS1_TS2_TV
 	    }
@@ -589,7 +680,7 @@ void distDNA_GG95(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        COUNT_TS_TV
 	    }
 	    P[target] = ((double) Ns / *s);
@@ -665,7 +756,7 @@ void distDNA_GG95_pairdel(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    Nd = Ns = L[target] = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        if (KnownBase(x[s1]) && KnownBase(x[s2])) L[target]++;
 		else continue;
 	        COUNT_TS_TV
@@ -721,7 +812,7 @@ void distDNA_GG95_pairdel(unsigned char *x, int *n, int *s, double *d,
 
 #define COMPUTE_DIST_LogDet\
     for (k = 0; k < 16; k++) Ftab[k] = ((double) Ntab[k]/L);\
-    d[target] = (-log(detFourByFour(Ftab))/4 - LN4)/4;\
+    d[target] = -log(detFourByFour(Ftab))/4 - LN4;\
     if (*variance) {\
         /* For the inversion, we first make U an identity matrix */\
         for (k = 1; k < 15; k++) U[k] = 0;\
@@ -751,7 +842,7 @@ void distDNA_LogDet(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    for (k = 0; k < 16; k++) Ntab[k] = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 		DO_CONTINGENCY_NUCLEOTIDES
 	    }
 	    COMPUTE_DIST_LogDet
@@ -771,7 +862,7 @@ void distDNA_LogDet_pairdel(unsigned char *x, int *n, int *s, double *d,
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    for (k = 0; k < 16; k++) Ntab[k] = 0;
 	    L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 		CHECK_PAIRWISE_DELETION
 		DO_CONTINGENCY_NUCLEOTIDES
 	    }
@@ -791,13 +882,13 @@ void distDNA_BH87(unsigned char *x, int *n, int *s, double *d,
    currently not available.
    </FIXME> */
 {
-    int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4], ndim = 4, info, ipiv[16];
-    double P12[16], P21[16], U[16];
+    int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4];
+    double P12[16], P21[16];
 
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    for (k = 0; k < 16; k++) Ntab[k] = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 		DO_CONTINGENCY_NUCLEOTIDES
 	    }
 
@@ -895,7 +986,7 @@ void distDNA_ParaLin(unsigned char *x, int *n, int *s, double *d,
     for (i1 = 1; i1 < *n; i1++) {
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    for (k = 0; k < 16; k++) Ntab[k] = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 		DO_CONTINGENCY_NUCLEOTIDES
 	    }
 	    COMPUTE_DIST_ParaLin
@@ -939,7 +1030,7 @@ void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
         for (i2 = i1 + 1; i2 <= *n; i2++) {
 	    L = 0;
 	    for (k = 0; k < 16; k++) Ntab[k] = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
+	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1 += *n, s2 += *n) {
 	        CHECK_PAIRWISE_DELETION
 		DO_CONTINGENCY_NUCLEOTIDES
 	    }
@@ -949,23 +1040,33 @@ void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
     }
 }
 
+/* a hash table is much faster than switch (2012-01-10) */
 void BaseProportion(unsigned char *x, int *n, double *BF)
 {
-    int i, m;
-
-    m = 0;
-    for (i = 0; i < *n; i++) {
-        if (KnownBase(x[i])) {
-	    m++;
-	    switch (x[i]) {
-	    case 136 : BF[0]++; break;
-	    case 40 : BF[1]++; break;
-	    case 72 : BF[2]++; break;
-	    case 24 : BF[3]++; break;
-	    }
-	}
-    }
-    for (i = 0; i < 4; i++) BF[i] /= m;
+	int i;
+	double count[256];
+
+	memset(count, 0, 256*sizeof(double));
+
+	for (i = 0; i < *n; i++) count[x[i]]++;
+
+	BF[0] = count[136];
+	BF[1] = count[40];
+	BF[2] = count[72];
+	BF[3] = count[24];
+	BF[4] = count[192];
+	BF[5] = count[160];
+	BF[6] = count[144];
+	BF[7] = count[96];
+	BF[8] = count[80];
+	BF[9] = count[48];
+	BF[10] = count[224];
+	BF[11] = count[176];
+	BF[12] = count[208];
+	BF[13] = count[112];
+	BF[14] = count[240];
+	BF[15] = count[4];
+	BF[16] = count[2];
 }
 
 void SegSites(unsigned char *x, int *n, int *s, int *seg)
@@ -979,7 +1080,7 @@ void SegSites(unsigned char *x, int *n, int *s, int *seg)
 	basis = x[i];
 	i++;
 	while (i < *n * (j + 1)) {
-	    if (x[i] == basis) i++;
+	    if (!KnownBase(x[i]) || x[i] == basis) i++;
 	    else {
 	        seg[j] = 1;
 		break;
@@ -988,27 +1089,6 @@ void SegSites(unsigned char *x, int *n, int *s, int *seg)
     }
 }
 
-void NucleotideDiversity(unsigned char *x, int *n, int *s,
-			 int *pairdel, double *ans)
-{
-    int i1, i2, s1, s2, Nd, L;
-
-    if (!*pairdel) L = *s;
-
-    for (i1 = 1; i1 < *n; i1++) {
-        for (i2 = i1 + 1; i2 <= *n; i2++) {
-	    Nd = 0;
-	    if (*pairdel) L = 0;
-	    for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
-                CHECK_PAIRWISE_DELETION
-		if (DifferentBase(x[s1], x[s2])) Nd++;
-	    }
-	    *ans += ((double) Nd/L);
-	}
-    }
-    *ans /= (*n * (*n - 1)/2);
-}
-
 void GlobalDeletionDNA(unsigned char *x, int *n, int *s, int *keep)
 {
     int i, j;
@@ -1030,39 +1110,57 @@ void dist_dna(unsigned char *x, int *n, int *s, int *model, double *d,
 	      int *gamma, double *alpha)
 {
     switch (*model) {
-    case 1 : if (pairdel) distDNA_raw_pairdel(x, n, s, d);
-             else distDNA_raw(x, n, s, d); break;
-
+    case 1 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 1);
+             else distDNA_raw(x, n, s, d, 1); break;
     case 2 : if (pairdel) distDNA_JC69_pairdel(x, n, s, d, variance, var, gamma, alpha);
              else distDNA_JC69(x, n, s, d, variance, var, gamma, alpha); break;
-
     case 3 : if (pairdel) distDNA_K80_pairdel(x, n, s, d, variance, var, gamma, alpha);
              else distDNA_K80(x, n, s, d, variance, var, gamma, alpha); break;
-
     case 4 : if (pairdel) distDNA_F81_pairdel(x, n, s, d, BF, variance, var, gamma, alpha);
              else distDNA_F81(x, n, s, d, BF, variance, var, gamma, alpha); break;
-
     case 5 : if (pairdel) distDNA_K81_pairdel(x, n, s, d, variance, var);
              else distDNA_K81(x, n, s, d, variance, var); break;
-
     case 6 : if (pairdel) distDNA_F84_pairdel(x, n, s, d, BF, variance, var);
              else distDNA_F84(x, n, s, d, BF, variance, var); break;
-
     case 7 : if (pairdel) distDNA_T92_pairdel(x, n, s, d, BF, variance, var);
              else distDNA_T92(x, n, s, d, BF, variance, var); break;
-
     case 8 : if (pairdel) distDNA_TN93_pairdel(x, n, s, d, BF, variance, var, gamma, alpha);
              else distDNA_TN93(x, n, s, d, BF, variance, var, gamma, alpha); break;
-
     case 9 : if (pairdel) distDNA_GG95_pairdel(x, n, s, d, variance, var);
              else distDNA_GG95(x, n, s, d, variance, var); break;
-
     case 10 : if (pairdel) distDNA_LogDet_pairdel(x, n, s, d, variance, var);
               else distDNA_LogDet(x, n, s, d, variance, var); break;
-
     case 11 : distDNA_BH87(x, n, s, d, variance, var); break;
-
     case 12 : if (pairdel) distDNA_ParaLin_pairdel(x, n, s, d, variance, var);
               else distDNA_ParaLin(x, n, s, d, variance, var); break;
+    case 13 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 0);
+             else distDNA_raw(x, n, s, d, 0); break;
+    case 14 : if (pairdel) distDNA_TsTv(x, n, s, d, 1, 1);
+	    else distDNA_TsTv(x, n, s, d, 1, 0); break;
+    case 15 : if (pairdel) distDNA_TsTv(x, n, s, d, 0, 1);
+	    else distDNA_TsTv(x, n, s, d, 0, 0); break;
+    case 16 : distDNA_indel(x, n, s, d); break;
+    case 17 : distDNA_indelblock(x, n, s, d); break;
     }
 }
+
+void where(unsigned char *x, unsigned char *pat, int *s, int *p,
+	   int *ans, int *n)
+{
+	int i, j, k, ln;
+
+	ln = 0; /* local n */
+
+	for (i = 0; i <= *s - *p; i++) {
+		k = i; j = 0;
+		while (1) {
+			if (x[k] != pat[j]) break;
+			j++; k++;
+			if (j == *p) {
+				ans[ln++] = k - 1;
+				break;
+			}
+		}
+	}
+	*n = ln;
+}