cleaning of C files and update on simulation of OU process

[ape.git] / src / dist_dna.c

diff --git a/src/dist_dna.c b/src/dist_dna.c
index a289f159e5bafad3c353845b7a364c3c99d3c6e2..210479e37b2ba29c623e9379d7dfffe621ba24db 100644 (file)
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,6 +1,6 @@
-/* dist_dna.c       2010-07-14 */
+/* dist_dna.c       2011-06-23 */
 
-/* Copyright 2005-2010 Emmanuel Paradis
+/* Copyright 2005-2011 Emmanuel Paradis
 
 /* This file is part of the R-package `ape'. */
 /* See the file ../COPYING for licensing issues. */
@@ -535,8 +535,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;
 
@@ -559,8 +559,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
 
@@ -742,7 +742,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;\
@@ -812,8 +812,8 @@ 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++) {
@@ -970,23 +970,50 @@ void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
     }
 }
 
-void BaseProportion(unsigned char *x, int *n, double *BF, int *freq)
+/* void BaseProportion(unsigned char *x, int *n, double *BF, int *freq) */
+/* { */
+/*     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; */
+/*         } */
+/*     } */
+/*     } */
+/*     if (! *freq) for (i = 0; i < 4; i++) BF[i] /= m; */
+/* } */
+
+void BaseProportion(unsigned char *x, int *n, double *BF)
 {
-    int i, m;
+    int i;
 
-    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;
+           case 192 : BF[4]++; break;
+           case 160 : BF[5]++; break;
+           case 144 : BF[6]++; break;
+           case 96 : BF[7]++; break;
+           case 80 : BF[8]++; break;
+           case 48 : BF[9]++; break;
+           case 224 : BF[10]++; break;
+           case 176 : BF[11]++; break;
+           case 208 : BF[12]++; break;
+           case 112 : BF[13]++; break;
+           case 240 : BF[14]++; break;
+           case 4 : BF[15]++; break;
+           case 2 : BF[16]++; break;
            }
-       }
     }
-    if (! *freq) for (i = 0; i < 4; i++) BF[i] /= m;
 }
 
 void SegSites(unsigned char *x, int *n, int *s, int *seg)