X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=src%2Fdist_dna.c;h=2aa6da6be9531dd1c24ba7d5c90e66c9afe27cc1;hb=24fc6c03893f85a3f9ab3d088201b3731f3035b4;hp=f10bc4debeacd4b8acf7d80da75af5913b3d0b63;hpb=6dfbab243973c0c3fa2e6d02b190aefbe5a67280;p=ape.git

diff --git a/src/dist_dna.c b/src/dist_dna.c
index f10bc4d..2aa6da6 100644
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,4 +1,4 @@
-/* dist_dna.c       2008-01-19 */
+/* dist_dna.c       2009-09-16 */
 
 /* Copyright 2005-2008 Emmanuel Paradis
 
@@ -12,7 +12,7 @@
 #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,7 +64,7 @@ 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)
+void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled)
 {
     int i1, i2, s1, s2, target, Nd;
 
@@ -74,13 +74,14 @@ void distDNA_raw(unsigned char *x, int *n, int *s, double *d)
 	    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);
+	    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)
+void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d, int scaled)
 {
     int i1, i2, s1, s2, target, Nd, L;
 
@@ -92,7 +93,8 @@ void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d)
                 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++;
 	}
     }
@@ -949,7 +951,7 @@ void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
     }
 }
 
-void BaseProportion(unsigned char *x, int *n, double *BF)
+void BaseProportion(unsigned char *x, int *n, double *BF, int *freq)
 {
     int i, m;
 
@@ -965,7 +967,7 @@ void BaseProportion(unsigned char *x, int *n, double *BF)
 	    }
 	}
     }
-    for (i = 0; i < 4; i++) BF[i] /= m;
+    if (! *freq) for (i = 0; i < 4; i++) BF[i] /= m;
 }
 
 void SegSites(unsigned char *x, int *n, int *s, int *seg)
@@ -979,7 +981,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;
@@ -1009,8 +1011,8 @@ 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;
@@ -1043,5 +1045,7 @@ void dist_dna(unsigned char *x, int *n, int *s, int *model, double *d,
 
     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;
     }
 }