a fix in cophyloplot()

[ape.git] / src / dist_dna.c

diff --git a/src/dist_dna.c b/src/dist_dna.c
index 7174af6d013d961c80dba33b453c78ad5e3c8bfe..61f32d29ca3a55e5d526fe3b3de8a09d7e3f1c95 100644 (file)
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,4 +1,4 @@
-/* dist_dna.c       2012-01-10 */
+/* dist_dna.c       2012-02-13 */
 
 /* Copyright 2005-2012 Emmanuel Paradis
 
@@ -73,6 +73,9 @@ double detFourByFour(double *x)
     }\
     if (IsPyrimidine(x[s1]) && IsPyrimidine(x[s2])) Ns++;
 
+
+
+
 void distDNA_TsTv(unsigned char *x, int *n, int *s, double *d, int Ts, int pairdel)
 {
        int i1, i2, s1, s2, target, Nd, Ns;
@@ -81,7 +84,7 @@ void distDNA_TsTv(unsigned char *x, int *n, int *s, double *d, int Ts, int paird
        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) {
                                if (pairdel && !(KnownBase(x[s1]) && KnownBase(x[s2]))) continue;
                                COUNT_TS_TV
                        }
@@ -100,7 +103,7 @@ void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled)
        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++;
                        if (scaled) d[target] = ((double) Nd / *s);
                        else d[target] = ((double) Nd);
@@ -117,7 +120,7 @@ void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d, int scaled
     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++;
            }
@@ -150,7 +153,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++;
@@ -233,7 +236,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
            }
@@ -285,7 +288,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++;
            }
@@ -332,7 +335,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
@@ -351,7 +354,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
            }
@@ -392,7 +395,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
@@ -413,7 +416,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
            }
@@ -449,7 +452,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
@@ -470,7 +473,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
            }
@@ -546,7 +549,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
@@ -568,7 +571,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
            }
@@ -610,7 +613,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);
@@ -686,7 +689,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
@@ -772,7 +775,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
@@ -792,7 +795,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
            }
@@ -818,7 +821,7 @@ void distDNA_BH87(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
            }
 
@@ -916,7 +919,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
@@ -960,7 +963,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
            }
@@ -1115,6 +1118,7 @@ void dist_dna(unsigned char *x, int *n, int *s, int *model, double *d,
            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;
     }
 }