various bug fixes since the release of ape 3.0

[ape.git] / src / dist_dna.c

diff --git a/src/dist_dna.c b/src/dist_dna.c
index 61f32d29ca3a55e5d526fe3b3de8a09d7e3f1c95..289be088b37878de150195605ce39241be6bb539 100644 (file)
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,12 +1,12 @@
-/* dist_dna.c       2012-02-13 */
+/* dist_dna.c       2012-02-14 */
 
 /* 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
@@ -73,8 +73,83 @@ double detFourByFour(double *x)
     }\
     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, N;
+
+       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)
+
+#define DINDEX(i, j) n * (i - 1) - i*(i - 1)/2 + j - i - 1
+
+int give_index(int i, int j, int n)
+{
+       if (i > j) return(DINDEX(j, i));
+       else return(DINDEX(i, j));
+}
 
+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++;
+               }
+       }
+}
 
+#undef X
 
 void distDNA_TsTv(unsigned char *x, int *n, int *s, double *d, int Ts, int pairdel)
 {