- /*** Compute A ***/
- /* // compute the posterior of a gap, but I do not know how to use it...
- if (1) {
- double *a;
- a = calloc(3 * (l_ref + 2), sizeof(double));
- for (i = 1; i < l_query; ++i) {
- double sum = 0., *fi = f[i], *bi1 = b[i+1], qli1 = qual[i+1];
- int beg = 1, end = l_ref, x;
- uint8_t qyi1 = query[i+1];
- x = i - bw; beg = beg > x? beg : x;
- x = i + bw; end = end < x? end : x;
- for (k = beg; k <= end; ++k) {
- double *ak = a + 3 * k;
- int u, v11, v01, v10;
- double e;
- set_u(u, bw, i, k); set_u(v11, bw, i+1, k+1); set_u(v10, bw, i+1, k); set_u(v01, bw, i, k+1);
- e = k >= l_ref? 0 : (ref[k+1] > 3 || qyi1 > 3)? 1. : ref[k+1] == qyi1? 1. - qli1 : qli1 * EM;
- ak[0] += fi[u] * bi1[v11] * m[0] * e;
- ak[1] += fi[u] * bi1[v10+1] * m[1] * EI;
- ak[2] += fi[u] * bi1[v01+2] * m[2];
- }
- }
- for (k = 1; k < l_ref; ++k) {
- double sum = 0., *ak = a + 3 * k;
- sum += 1. / (ak[0] + ak[1] + ak[2]);
- ak[0] *= sum; ak[1] *= sum; ak[2] *= sum;
- fprintf(stderr, "%d: %lf, %lf, %lf\n", k, ak[0], ak[1], ak[2]);
- }
- free(a);
- }
- */