allow to read the prior from the error output. EM iteration is working.

[samtools.git] / bcftools / prob1.c

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 1acd49b48295d4d315f52c2f3969bd6c3b66e36f..4611af1f995fff8832300af10f733bffa41dd29d 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -2,8 +2,12 @@
 #include <stdlib.h>
 #include <string.h>
 #include <stdio.h>
+#include <errno.h>
 #include "prob1.h"
 
+#include "kseq.h"
+KSTREAM_INIT(gzFile, gzread, 16384)
+
 #define MC_AVG_ERR 0.007
 #define MC_MAX_EM_ITER 16
 #define MC_EM_EPS 1e-4
@@ -54,6 +58,43 @@ void bcf_p1_init_prior(bcf_p1aux_t *ma, int type, double theta)
        }
 }
 
+int bcf_p1_read_prior(bcf_p1aux_t *ma, const char *fn)
+{
+       gzFile fp;
+       kstring_t s;
+       kstream_t *ks;
+       long double sum;
+       int dret, k;
+       memset(&s, 0, sizeof(kstring_t));
+       fp = strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
+       ks = ks_init(fp);
+       memset(ma->phi, 0, sizeof(double) * (ma->M + 1));
+       while (ks_getuntil(ks, '\n', &s, &dret) >= 0) {
+               if (strstr(s.s, "[afs] ") == s.s) {
+                       char *p = s.s + 6;
+                       for (k = 0; k <= ma->M; ++k) {
+                               int x;
+                               double y;
+                               x = strtol(p, &p, 10);
+                               if (x != k && (errno == EINVAL || errno == ERANGE)) return -1;
+                               ++p;
+                               y = strtod(p, &p);
+                               if (y == 0. && (errno == EINVAL || errno == ERANGE)) return -1;
+                               ma->phi[ma->M - k] += y;
+                       }
+               }
+       }
+       ks_destroy(ks);
+       gzclose(fp);
+       free(s.s);
+       for (sum = 0., k = 0; k <= ma->M; ++k) sum += ma->phi[k];
+       fprintf(stderr, "[prior]");
+       for (k = 0; k <= ma->M; ++k) ma->phi[k] /= sum;
+       for (k = 0; k <= ma->M; ++k) fprintf(stderr, " %d:%.3lg", k, ma->phi[ma->M - k]);
+       fputc('\n', stderr);
+       return 0;
+}
+
 bcf_p1aux_t *bcf_p1_init(int n) // FIXME: assuming diploid
 {
        bcf_p1aux_t *ma;
@@ -148,20 +189,22 @@ int bcf_p1_call_gt(const bcf_p1aux_t *ma, double f0, int k)
 
 static void mc_cal_y(bcf_p1aux_t *ma)
 {
-       double *z[2], *tmp, *pdg, last_min, last_max;
-       int k, j;
+       double *z[2], *tmp, *pdg;
+       int k, j, last_min, last_max;
        z[0] = ma->z;
        z[1] = ma->zswap;
        pdg = ma->pdg;
-       z[0][0] = 1.; z[0][1] = z[0][2] = 0.;
+       memset(z[0], 0, sizeof(double) * (ma->M + 1));
+       memset(z[1], 0, sizeof(double) * (ma->M + 1));
+       z[0][0] = 1.;
        last_min = last_max = 0;
        for (j = 0; j < ma->n; ++j) {
                int _min = last_min, _max = last_max;
                double p[3], sum;
                pdg = ma->pdg + j * 3;
                p[0] = pdg[0]; p[1] = 2. * pdg[1]; p[2] = pdg[2];
-//             for (; _min < _max && z[0][_min] < TINY; ++_min) z[1][_min] = 0.;
-//             for (; _max > _min && z[0][_max] < TINY; --_max) z[1][_max] = 0.;
+               for (; _min < _max && z[0][_min] < TINY; ++_min) z[0][_min] = z[1][_min] = 0.;
+               for (; _max > _min && z[0][_max] < TINY; --_max) z[0][_max] = z[1][_max] = 0.;
                _max += 2;
                if (_min == 0) 
                        k = 0, z[1][k] = (2*j+2-k)*(2*j-k+1) * p[0] * z[0][k];
@@ -173,6 +216,8 @@ static void mc_cal_y(bcf_p1aux_t *ma)
                                + k*(k-1)* p[2] * z[0][k-2];
                for (k = _min, sum = 0.; k <= _max; ++k) sum += z[1][k];
                for (k = _min; k <= _max; ++k) z[1][k] /= sum;
+               if (_min >= 1) z[1][_min-1] = 0.;
+               if (_min >= 2) z[1][_min-2] = 0.;
                if (j < ma->n - 1) z[1][_max+1] = z[1][_max+2] = 0.;
                tmp = z[0]; z[0] = z[1]; z[1] = tmp;
                last_min = _min; last_max = _max;