faster for large sample size (in principle)

diff --git a/bcftools/bcf.c b/bcftools/bcf.c
index 2cb98761cf570d7e6ba6260217977c536158301e..c3313881ef90fc2930c72d218c090b35e5cf9466 100644 (file)
--- a/bcftools/bcf.c
+++ b/bcftools/bcf.c
@@ -223,6 +223,7 @@ void bcf_fmt_core(const bcf_hdr_t *h, bcf1_t *b, kstring_t *s)
                fmt_str(b->fmt, s);
        }
        x = b->n_alleles * (b->n_alleles + 1) / 2;
+       if (b->n_gi == 0) return;
        for (j = 0; j < h->n_smpl; ++j) {
                kputc('\t', s);
                for (i = 0; i < b->n_gi; ++i) {

diff --git a/bcftools/prob1.c b/bcftools/prob1.c
index 1acd49b48295d4d315f52c2f3969bd6c3b66e36f..4afda56e33d3f8a8f4f9fe511e1d510aade79974 100644 (file)
--- a/bcftools/prob1.c
+++ b/bcftools/prob1.c
@@ -148,20 +148,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 +175,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;