kaln.c

   1 /* The MIT License
   2
   3    Copyright (c) 2003-2006, 2008, 2009, by Heng Li <lh3lh3@gmail.com>
   4
   5    Permission is hereby granted, free of charge, to any person obtaining
   6    a copy of this software and associated documentation files (the
   7    "Software"), to deal in the Software without restriction, including
   8    without limitation the rights to use, copy, modify, merge, publish,
   9    distribute, sublicense, and/or sell copies of the Software, and to
  10    permit persons to whom the Software is furnished to do so, subject to
  11    the following conditions:
  12
  13    The above copyright notice and this permission notice shall be
  14    included in all copies or substantial portions of the Software.
  15
  16    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  18    MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  20    BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  21    ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
  22    CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  23    SOFTWARE.
  24 */
  25
  26 #include <stdlib.h>
  27 #include <stdio.h>
  28 #include <string.h>
  29 #include <stdint.h>
  30 #include "kaln.h"
  31
  32 #define FROM_M 0
  33 #define FROM_I 1
  34 #define FROM_D 2
  35
  36 typedef struct {
  37         int i, j;
  38         unsigned char ctype;
  39 } path_t;
  40
  41 int aln_sm_blosum62[] = {
  42 /*       A  R  N  D  C  Q  E  G  H  I  L  K  M  F  P  S  T  W  Y  V  *  X */
  43          4,-1,-2,-2, 0,-1,-1, 0,-2,-1,-1,-1,-1,-2,-1, 1, 0,-3,-2, 0,-4, 0,
  44         -1, 5, 0,-2,-3, 1, 0,-2, 0,-3,-2, 2,-1,-3,-2,-1,-1,-3,-2,-3,-4,-1,
  45         -2, 0, 6, 1,-3, 0, 0, 0, 1,-3,-3, 0,-2,-3,-2, 1, 0,-4,-2,-3,-4,-1,
  46         -2,-2, 1, 6,-3, 0, 2,-1,-1,-3,-4,-1,-3,-3,-1, 0,-1,-4,-3,-3,-4,-1,
  47          0,-3,-3,-3, 9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1,-4,-2,
  48         -1, 1, 0, 0,-3, 5, 2,-2, 0,-3,-2, 1, 0,-3,-1, 0,-1,-2,-1,-2,-4,-1,
  49         -1, 0, 0, 2,-4, 2, 5,-2, 0,-3,-3, 1,-2,-3,-1, 0,-1,-3,-2,-2,-4,-1,
  50          0,-2, 0,-1,-3,-2,-2, 6,-2,-4,-4,-2,-3,-3,-2, 0,-2,-2,-3,-3,-4,-1,
  51         -2, 0, 1,-1,-3, 0, 0,-2, 8,-3,-3,-1,-2,-1,-2,-1,-2,-2, 2,-3,-4,-1,
  52         -1,-3,-3,-3,-1,-3,-3,-4,-3, 4, 2,-3, 1, 0,-3,-2,-1,-3,-1, 3,-4,-1,
  53         -1,-2,-3,-4,-1,-2,-3,-4,-3, 2, 4,-2, 2, 0,-3,-2,-1,-2,-1, 1,-4,-1,
  54         -1, 2, 0,-1,-3, 1, 1,-2,-1,-3,-2, 5,-1,-3,-1, 0,-1,-3,-2,-2,-4,-1,
  55         -1,-1,-2,-3,-1, 0,-2,-3,-2, 1, 2,-1, 5, 0,-2,-1,-1,-1,-1, 1,-4,-1,
  56         -2,-3,-3,-3,-2,-3,-3,-3,-1, 0, 0,-3, 0, 6,-4,-2,-2, 1, 3,-1,-4,-1,
  57         -1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4, 7,-1,-1,-4,-3,-2,-4,-2,
  58          1,-1, 1, 0,-1, 0, 0, 0,-1,-2,-2, 0,-1,-2,-1, 4, 1,-3,-2,-2,-4, 0,
  59          0,-1, 0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1, 1, 5,-2,-2, 0,-4, 0,
  60         -3,-3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1, 1,-4,-3,-2,11, 2,-3,-4,-2,
  61         -2,-2,-2,-3,-2,-1,-2,-3, 2,-1,-1,-2,-1, 3,-3,-2,-2, 2, 7,-1,-4,-1,
  62          0,-3,-3,-3,-1,-2,-2,-3,-3, 3, 1,-2, 1,-1,-2,-2, 0,-3,-1, 4,-4,-1,
  63         -4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4, 1,-4,
  64          0,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2, 0, 0,-2,-1,-1,-4,-1
  65 };
  66
  67 int aln_sm_blast[] = {
  68         1, -3, -3, -3, -2,
  69         -3, 1, -3, -3, -2,
  70         -3, -3, 1, -3, -2,
  71         -3, -3, -3, 1, -2,
  72         -2, -2, -2, -2, -2
  73 };
  74
  75 ka_param_t ka_param_blast = {  5,  2,   5, 2, aln_sm_blast, 5, 50 };
  76 ka_param_t ka_param_aa2aa = { 10,  2,  10, 2, aln_sm_blosum62, 22, 50 };
  77
  78 static uint32_t *ka_path2cigar32(const path_t *path, int path_len, int *n_cigar)
  79 {
  80         int i, n;
  81         uint32_t *cigar;
  82         unsigned char last_type;
  83
  84         if (path_len == 0 || path == 0) {
  85                 *n_cigar = 0;
  86                 return 0;
  87         }
  88
  89         last_type = path->ctype;
  90         for (i = n = 1; i < path_len; ++i) {
  91                 if (last_type != path[i].ctype) ++n;
  92                 last_type = path[i].ctype;
  93         }
  94         *n_cigar = n;
  95         cigar = (uint32_t*)calloc(*n_cigar, 4);
  96
  97         cigar[0] = 1u << 4 | path[path_len-1].ctype;
  98         last_type = path[path_len-1].ctype;
  99         for (i = path_len - 2, n = 0; i >= 0; --i) {
 100                 if (path[i].ctype == last_type) cigar[n] += 1u << 4;
 101                 else {
 102                         cigar[++n] = 1u << 4 | path[i].ctype;
 103                         last_type = path[i].ctype;
 104                 }
 105         }
 106
 107         return cigar;
 108 }
 109
 110 /***************************/
 111 /* START OF common_align.c */
 112 /***************************/
 113
 114 #define SET_INF(s) (s).M = (s).I = (s).D = MINOR_INF;
 115
 116 #define set_M(MM, cur, p, sc)                                                   \
 117 {                                                                                                               \
 118         if ((p)->M >= (p)->I) {                                                         \
 119                 if ((p)->M >= (p)->D) {                                                 \
 120                         (MM) = (p)->M + (sc); (cur)->Mt = FROM_M;       \
 121                 } else {                                                                                \
 122                         (MM) = (p)->D + (sc); (cur)->Mt = FROM_D;       \
 123                 }                                                                                               \
 124         } else {                                                                                        \
 125                 if ((p)->I > (p)->D) {                                                  \
 126                         (MM) = (p)->I + (sc); (cur)->Mt = FROM_I;       \
 127                 } else {                                                                                \
 128                         (MM) = (p)->D + (sc); (cur)->Mt = FROM_D;       \
 129                 }                                                                                               \
 130         }                                                                                                       \
 131 }
 132 #define set_I(II, cur, p)                                                               \
 133 {                                                                                                               \
 134         if ((p)->M - gap_open > (p)->I) {                                       \
 135                 (cur)->It = FROM_M;                                                             \
 136                 (II) = (p)->M - gap_open - gap_ext;                             \
 137         } else {                                                                                        \
 138                 (cur)->It = FROM_I;                                                             \
 139                 (II) = (p)->I - gap_ext;                                                \
 140         }                                                                                                       \
 141 }
 142 #define set_end_I(II, cur, p)                                                   \
 143 {                                                                                                               \
 144         if (gap_end_ext >= 0) {                                                         \
 145                 if ((p)->M - gap_end_open > (p)->I) {                   \
 146                         (cur)->It = FROM_M;                                                     \
 147                         (II) = (p)->M - gap_end_open - gap_end_ext;     \
 148                 } else {                                                                                \
 149                         (cur)->It = FROM_I;                                                     \
 150                         (II) = (p)->I - gap_end_ext;                            \
 151                 }                                                                                               \
 152         } else set_I(II, cur, p);                                                       \
 153 }
 154 #define set_D(DD, cur, p)                                                               \
 155 {                                                                                                               \
 156         if ((p)->M - gap_open > (p)->D) {                                       \
 157                 (cur)->Dt = FROM_M;                                                             \
 158                 (DD) = (p)->M - gap_open - gap_ext;                             \
 159         } else {                                                                                        \
 160                 (cur)->Dt = FROM_D;                                                             \
 161                 (DD) = (p)->D - gap_ext;                                                \
 162         }                                                                                                       \
 163 }
 164 #define set_end_D(DD, cur, p)                                                   \
 165 {                                                                                                               \
 166         if (gap_end_ext >= 0) {                                                         \
 167                 if ((p)->M - gap_end_open > (p)->D) {                   \
 168                         (cur)->Dt = FROM_M;                                                     \
 169                         (DD) = (p)->M - gap_end_open - gap_end_ext;     \
 170                 } else {                                                                                \
 171                         (cur)->Dt = FROM_D;                                                     \
 172                         (DD) = (p)->D - gap_end_ext;                            \
 173                 }                                                                                               \
 174         } else set_D(DD, cur, p);                                                       \
 175 }
 176
 177 typedef struct {
 178         uint8_t Mt:3, It:2, Dt:2;
 179 } dpcell_t;
 180
 181 typedef struct {
 182         int M, I, D;
 183 } dpscore_t;
 184
 185 /***************************
 186  * banded global alignment *
 187  ***************************/
 188 uint32_t *ka_global_core(uint8_t *seq1, int len1, uint8_t *seq2, int len2, const ka_param_t *ap, int *_score, int *n_cigar)
 189 {
 190         int i, j;
 191         dpcell_t **dpcell, *q;
 192         dpscore_t *curr, *last, *s;
 193         int b1, b2, tmp_end;
 194         int *mat, end, max = 0;
 195         uint8_t type, ctype;
 196         uint32_t *cigar = 0;
 197
 198         int gap_open, gap_ext, gap_end_open, gap_end_ext, b;
 199         int *score_matrix, N_MATRIX_ROW;
 200
 201         /* initialize some align-related parameters. just for compatibility */
 202         gap_open = ap->gap_open;
 203         gap_ext = ap->gap_ext;
 204         gap_end_open = ap->gap_end_open;
 205         gap_end_ext = ap->gap_end_ext;
 206         b = ap->band_width;
 207         score_matrix = ap->matrix;
 208         N_MATRIX_ROW = ap->row;
 209
 210         *n_cigar = 0;
 211         if (len1 == 0 || len2 == 0) return 0;
 212
 213         /* calculate b1 and b2 */
 214         if (len1 > len2) {
 215                 b1 = len1 - len2 + b;
 216                 b2 = b;
 217         } else {
 218                 b1 = b;
 219                 b2 = len2 - len1 + b;
 220         }
 221         if (b1 > len1) b1 = len1;
 222         if (b2 > len2) b2 = len2;
 223         --seq1; --seq2;
 224
 225         /* allocate memory */
 226         end = (b1 + b2 <= len1)? (b1 + b2 + 1) : (len1 + 1);
 227         dpcell = (dpcell_t**)malloc(sizeof(dpcell_t*) * (len2 + 1));
 228         for (j = 0; j <= len2; ++j)
 229                 dpcell[j] = (dpcell_t*)malloc(sizeof(dpcell_t) * end);
 230         for (j = b2 + 1; j <= len2; ++j)
 231                 dpcell[j] -= j - b2;
 232         curr = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1));
 233         last = (dpscore_t*)malloc(sizeof(dpscore_t) * (len1 + 1));
 234
 235         /* set first row */
 236         SET_INF(*curr); curr->M = 0;
 237         for (i = 1, s = curr + 1; i < b1; ++i, ++s) {
 238                 SET_INF(*s);
 239                 set_end_D(s->D, dpcell[0] + i, s - 1);
 240         }
 241         s = curr; curr = last; last = s;
 242
 243         /* core dynamic programming, part 1 */
 244         tmp_end = (b2 < len2)? b2 : len2 - 1;
 245         for (j = 1; j <= tmp_end; ++j) {
 246                 q = dpcell[j]; s = curr; SET_INF(*s);
 247                 set_end_I(s->I, q, last);
 248                 end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1;
 249                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 250                 ++s; ++q;
 251                 for (i = 1; i != end; ++i, ++s, ++q) {
 252                         set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */
 253                         set_I(s->I, q, last + i);
 254                         set_D(s->D, q, s - 1);
 255                 }
 256                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 257                 set_D(s->D, q, s - 1);
 258                 if (j + b1 - 1 > len1) { /* bug fixed, 040227 */
 259                         set_end_I(s->I, q, last + i);
 260                 } else s->I = MINOR_INF;
 261                 s = curr; curr = last; last = s;
 262         }
 263         /* last row for part 1, use set_end_D() instead of set_D() */
 264         if (j == len2 && b2 != len2 - 1) {
 265                 q = dpcell[j]; s = curr; SET_INF(*s);
 266                 set_end_I(s->I, q, last);
 267                 end = (j + b1 <= len1 + 1)? (j + b1 - 1) : len1;
 268                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 269                 ++s; ++q;
 270                 for (i = 1; i != end; ++i, ++s, ++q) {
 271                         set_M(s->M, q, last + i - 1, mat[seq1[i]]); /* this will change s->M ! */
 272                         set_I(s->I, q, last + i);
 273                         set_end_D(s->D, q, s - 1);
 274                 }
 275                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 276                 set_end_D(s->D, q, s - 1);
 277                 if (j + b1 - 1 > len1) { /* bug fixed, 040227 */
 278                         set_end_I(s->I, q, last + i);
 279                 } else s->I = MINOR_INF;
 280                 s = curr; curr = last; last = s;
 281                 ++j;
 282         }
 283
 284         /* core dynamic programming, part 2 */
 285         for (; j <= len2 - b2 + 1; ++j) {
 286                 SET_INF(curr[j - b2]);
 287                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 288                 end = j + b1 - 1;
 289                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i != end; ++i, ++s, ++q) {
 290                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 291                         set_I(s->I, q, last + i);
 292                         set_D(s->D, q, s - 1);
 293                 }
 294                 set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 295                 set_D(s->D, q, s - 1);
 296                 s->I = MINOR_INF;
 297                 s = curr; curr = last; last = s;
 298         }
 299
 300         /* core dynamic programming, part 3 */
 301         for (; j < len2; ++j) {
 302                 SET_INF(curr[j - b2]);
 303                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 304                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) {
 305                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 306                         set_I(s->I, q, last + i);
 307                         set_D(s->D, q, s - 1);
 308                 }
 309                 set_M(s->M, q, last + len1 - 1, mat[seq1[i]]);
 310                 set_end_I(s->I, q, last + i);
 311                 set_D(s->D, q, s - 1);
 312                 s = curr; curr = last; last = s;
 313         }
 314         /* last row */
 315         if (j == len2) {
 316                 SET_INF(curr[j - b2]);
 317                 mat = score_matrix + seq2[j] * N_MATRIX_ROW;
 318                 for (i = j - b2 + 1, q = dpcell[j] + i, s = curr + i; i < len1; ++i, ++s, ++q) {
 319                         set_M(s->M, q, last + i - 1, mat[seq1[i]]);
 320                         set_I(s->I, q, last + i);
 321                         set_end_D(s->D, q, s - 1);
 322                 }
 323                 set_M(s->M, q, last + len1 - 1, mat[seq1[i]]);
 324                 set_end_I(s->I, q, last + i);
 325                 set_end_D(s->D, q, s - 1);
 326                 s = curr; curr = last; last = s;
 327         }
 328
 329         *_score = last[len1].M;
 330         if (n_cigar) { /* backtrace */
 331                 path_t *p, *path = (path_t*)malloc(sizeof(path_t) * (len1 + len2 + 2));
 332                 i = len1; j = len2;
 333                 q = dpcell[j] + i;
 334                 s = last + len1;
 335                 max = s->M; type = q->Mt; ctype = FROM_M;
 336                 if (s->I > max) { max = s->I; type = q->It; ctype = FROM_I; }
 337                 if (s->D > max) { max = s->D; type = q->Dt; ctype = FROM_D; }
 338
 339                 p = path;
 340                 p->ctype = ctype; p->i = i; p->j = j; /* bug fixed 040408 */
 341                 ++p;
 342                 do {
 343                         switch (ctype) {
 344                         case FROM_M: --i; --j; break;
 345                         case FROM_I: --j; break;
 346                         case FROM_D: --i; break;
 347                         }
 348                         q = dpcell[j] + i;
 349                         ctype = type;
 350                         switch (type) {
 351                         case FROM_M: type = q->Mt; break;
 352                         case FROM_I: type = q->It; break;
 353                         case FROM_D: type = q->Dt; break;
 354                         }
 355                         p->ctype = ctype; p->i = i; p->j = j;
 356                         ++p;
 357                 } while (i || j);
 358                 cigar = ka_path2cigar32(path, p - path - 1, n_cigar);
 359                 free(path);
 360         }
 361
 362         /* free memory */
 363         for (j = b2 + 1; j <= len2; ++j)
 364                 dpcell[j] += j - b2;
 365         for (j = 0; j <= len2; ++j)
 366                 free(dpcell[j]);
 367         free(dpcell);
 368         free(curr); free(last);
 369
 370         return cigar;
 371 }
 372
 373 #define get_k(b, i, kk) ((i) < (b)? (kk)/3 : (kk)/3+(i)-(b))
 374 #define set_u(u, b, i, k) { int x=(i)-(b); x=x>0?x:0; (u)=((k)-x+1)*3; }
 375
 376 ka_probpar_t ka_probpar_def = { 0.0001, 0.1, 10 };
 377
 378 int ka_prob_extend(uint8_t *_ref, int l_ref, uint8_t *_query, int l_query, float *_qual,
 379                                    const ka_probpar_t *c, int *state, uint8_t *q)
 380 {
 381         double **f, **b, m[9], f_sink, sI, sM, pf, pb;
 382         float *qual;
 383         uint8_t *ref, *query;
 384         int bw, bw2, i, k;
 385         /*** initialization ***/
 386         ref = _ref - 1;
 387         query = _query - 1;
 388         qual = _qual - 1;
 389         bw = c->bw;
 390         bw2 = c->bw * 2 + 1;
 391         f = calloc(l_query+1, sizeof(void*));
 392         b = calloc(l_query+1, sizeof(void*));
 393         for (i = 0; i <= l_query; ++i) {
 394                 f[i] = calloc(bw2 * 3 + 6, sizeof(double));
 395                 b[i] = calloc(bw2 * 3 + 6, sizeof(double));
 396         }
 397         // initialize m
 398         sM = sI = 1. / (2 * l_query + 1);
 399         m[0*3+0] = (1 - c->d - c->d) * (1 - sM);
 400         m[0*3+1] = m[0*3+2] = c->d * (1 - sM);
 401         m[1*3+0] = (1 - c->e) * (1 - sI);
 402         m[1*3+1] = c->e * (1 - sI);
 403         m[1*3+2] = 0.;
 404         m[2*3+0] = 1 - c->e;
 405         m[2*3+1] = 0.;
 406         m[2*3+2] = c->e;
 407         /*** forward ***/
 408         // f[0]
 409         set_u(k, bw, 0, 0);
 410         f[0][k] = 1.;
 411         // f[1..l_query]; core loop
 412         for (i = 1; i <= l_query; ++i) {
 413                 double *fi = f[i], *fi1 = f[i-1];
 414                 int beg = 0, end = l_ref, x;
 415                 x = i - bw; beg = beg > x? beg : x;
 416                 x = i + bw; end = end < x? end : x;
 417                 for (k = beg; k <= end; ++k) {
 418                         int u, v11, v01, v10;
 419                         double e;
 420                         e = (ref[k] == 4 || query[i] == 4)? 1. : ref[k] == query[i]? 1. - qual[i] : qual[i]; // FIXME: can be better
 421                         set_u(u, bw, i, k);
 422                         set_u(v11, bw, i-1, k-1);
 423                         set_u(v10, bw, i-1, k);
 424                         set_u(v01, bw, i, k-1);
 425                         fi[u+0] = e * (m[0] * fi1[v11+0] + m[3] * fi1[v11+1] + m[6] * fi1[v11+2]);
 426                         fi[u+1] = .25 * (m[1] * fi1[v10+0] + m[4] * fi1[v10+1]);
 427                         fi[u+2] = m[2] * fi[v01+0] + m[8] * fi[v01+2];
 428                         fprintf(stderr, "{%d,%d}@%d %.10lg,%.10lg,%.10lg\n", i, k, u, fi[u], fi[u+1], fi[u+2]);
 429                 }
 430         }
 431         // sink
 432         for (k = 0, f_sink = 0.; k <= l_ref; ++k) {
 433                 int u;
 434                 set_u(u, bw, l_query, k);
 435                 if (u < 3 || u >= bw2*3+3) continue;
 436                 f_sink += f[l_query][u+0] * sM + f[l_query][u+1] * sI;
 437         }
 438         pf = f_sink;
 439         /*** backward ***/
 440         // sink
 441         for (k = 0; k <= l_ref; ++k) {
 442                 int u;
 443                 double *bi = b[l_query];
 444                 set_u(u, bw, l_query, k);
 445                 if (u < 3 || u >= bw2*3+3) continue;
 446                 bi[u+0] = sM; bi[u+1] = sI;
 447                 fprintf(stderr, "[%d,%d]@%d %.10lg,%.10lg,%.10lg\n", l_query, k, u, bi[u], bi[u+1], bi[u+2]);
 448         }
 449         // b[l_query..1]; core loop
 450         for (i = l_query - 1; i > 0; --i) {
 451                 int beg = 0, end = l_ref, x;
 452                 double *bi = b[i], *bi1 = b[i+1];
 453                 x = i - bw; beg = beg > x? beg : x;
 454                 x = i + bw; end = end < x? end : x;
 455                 for (k = end; k >= beg; --k) {
 456                         int u, v11, v01, v10;
 457                         double e;
 458                         e = k >= l_ref? 0 : (ref[k+1] == 4 || query[i+1] == 4)? 1. : ref[k+1] == query[i+1]? 1. - qual[i+1] : qual[i+1];
 459                         set_u(u, bw, i, k);
 460                         set_u(v11, bw, i+1, k+1);
 461                         set_u(v10, bw, i+1, k);
 462                         set_u(v01, bw, i, k+1);
 463                         bi[u+0] = e * m[0] * bi1[v11+0] + .25 * m[1] * bi1[v10+1] + m[2] * bi[v01+2];
 464                         bi[u+1] = e * m[3] * bi1[v11+0] + .25 * m[4] * bi1[v10+1];
 465                         bi[u+2] = e * m[6] * bi1[v11+0] + m[8] * bi[v01+2];
 466                         fprintf(stderr, "[%d,%d]@%d %.10lg,%.10lg,%.10lg\n", i, k, u, bi[u], bi[u+1], bi[u+2]);
 467                 }
 468         }
 469         {
 470                 int beg = 0, end = l_ref, x;
 471                 double *bi = b[0], *bi1 = b[1];
 472                 x = i - bw; beg = beg > x? beg : x;
 473                 x = i + bw; end = end < x? end : x;
 474                 for (k = end; k >= beg; --k) {
 475                         int u, v11, v01, v10;
 476                         double e;
 477                         e = k >= l_ref? 0 : (ref[k+1] == 4 || query[1] == 4)? 1. : ref[k+1] == query[1]? 1. - qual[1] : qual[1];
 478                         set_u(u, bw, 0, k);
 479                         set_u(v11, bw, 1, k+1);
 480                         set_u(v10, bw, 1, k);
 481                         set_u(v01, bw, 0, k+1);
 482                         bi[u+0] = e * m[0] * bi1[v11+0] + .25 * m[1] * bi1[v10+1] + m[2] * bi[v01+2];
 483                         fprintf(stderr, "[%d,%d]@%d %.10lg,%.10lg,%.10lg\n", i, k, u, bi[u], bi[u+1], bi[u+2]);
 484                 }
 485         }
 486         set_u(k, bw, 0, 0);
 487         pb = b[0][k];
 488         fprintf(stderr, "(%.10lg,%.10lg):(%.10lg,%.10lg)\n", pf, pb, pb-pf, pb/pf-1);
 489         /*** MAP ***/
 490         return 0;
 491 }
 492
 493 #ifdef _MAIN
 494 int main()
 495 {
 496         int l_ref = 5, l_query = 4;
 497         uint8_t *ref = (uint8_t*)"\0\1\3\3\1";
 498         uint8_t *query = (uint8_t*)"\0\1\3\1";
 499 //      uint8_t *query = "\1\3\3\1";
 500         static float qual[4] = {.2, .01, .01, .01};
 501         ka_prob_extend(ref, l_ref, query, l_query, qual, &ka_probpar_def, 0, 0);
 502         return 0;
 503 }
 504 #endif