1 /* fast-dbviterbi.c
2 * SRE, Fri Sep 30 15:09:06 1994
3 *
4 * Search with a covariance model
5 * Fast "banded search" version of dbviterbi.c; subsequence
6 * lengths are bounded by probabilistically determined bounds.
7 *
8 **********
9 *
10 * To optimize memory access patterns, the score storage is implemented
11 * as a two-matrix version. amx is the
12 * main storage. bmx is a smaller auxiliary matrix with a different
13 * access pattern, holding scores of BEGIN state alignments; it
14 * is used when calculating BIFURC scores.
15 *
16 * amx is [j = 0..1] [y = 0..statenum] [diff = 0..j]
17 * diff == 0 is for off-diagonal boundary conditions (this is why diff is shifted +1)
18 * diff == 1 is for the diagonal, i==j
19 * We only need to keep two j rows in memory (current and previous).
20 * Note that this is yet *another* memory access pattern and it's different
21 * from dbviterbi.c!!!
22 *
23 * bmx is [y = 0..statenum] [j = 0..N] [ diff = 0..j]
24 * a j,diff matrix exists only where y is a BEGIN state
25 *
26 * The 2.0 implementation allows variable storage per node rather
27 * than storing and calculating a fixed max number of states per node,
28 * which should save up to 2x in both time and space.
29 *
30 * An optimization is made which requires END states to be explicitly
31 * added, so statenum (the number of states in the integer model)
32 * is *inclusive* of ENDs.
33 */
34
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <math.h>
39
40 #include "squid.h"
41 #include "structs.h"
42 #include "funcs.h"
43
44 #ifdef MEMDEBUG
45 #include "dbmalloc.h"
46 #endif
47
48 static int allocate_mx(struct istate_s *icm, int statenum, int window,
49 int ****ret_amx, int ****ret_bmx);
50 static int init_mx (struct istate_s *icm, int statenum, int N,
51 int ***amx, int ***bmx);
52 static int recurse_mx (struct istate_s *icm, int statenum, int *minb, int *maxb,
53 char *seq, int seqlen,
54 int window, int ***amx, int ***bmx, int ithresh,
55 int (*gotone_f)(int, int, double));
56 static void free_mx (int ***amx, int ***bmx, int statenum, int window);
57
58
59 /* Function: FastViterbiScan()
60 *
61 * Purpose: Scanning version of the Viterbi alignment algorithm,
62 * for finding matches in a long sequence.
63 *
64 * Args: icm - the model to align sequence to (int log-odds)
65 * statenum - length of model in states (inclusive of END)
66 * minb - minimum length bounds for states
67 * maxb - maximum length bounds for states
68 * seq - sequence to align model to
69 * window - scanning window size (nucleotides)
70 * thresh - scores above this are reported through gotone_f()
71 * gotone_f - function which gets told about a match
72 *
73 * Return: 1 on success, 0 on failure.
74 */
75 int
FastViterbiScan(struct istate_s * icm,int statenum,int * minb,int * maxb,char * seq,int window,double thresh,int (* gotone_f)(int,int,double))76 FastViterbiScan(struct istate_s *icm, int statenum, int *minb, int *maxb,
77 char *seq, int window, double thresh,
78 int (*gotone_f)(int, int, double))
79 {
80 int ***amx; /* the main score matrix */
81 int ***bmx; /* the BEGIN score matrix */
82 int N; /* length of sequence */
83 int ithresh; /* thresh, converted and scaled to int */
84
85 N = strlen(seq);
86 seq--; /* convert to 1..N. Ugh! */
87 ithresh = (int) (thresh * INTPRECISION);
88
89 if (! allocate_mx(icm, statenum, window, &amx, &bmx)) return 0;
90 #ifdef DEBUG
91 printf("allocated matrices\n");
92 #endif
93
94 if (! init_mx(icm, statenum, window, amx, bmx)) return 0;
95 #ifdef DEBUG
96 printf("matrices initialized\n");
97 #endif
98
99 if (! recurse_mx(icm, statenum, minb, maxb, seq, N, window, amx, bmx, ithresh, gotone_f))
100 return 0;
101 #ifdef DEBUG
102 printf("recursion finished\n");
103 #endif
104 /* terminate scanning hit reporting */
105 ReportScanHit(-1,-1, 0.0, gotone_f);
106 free_mx(amx, bmx, statenum, window);
107 return 1;
108 }
109
110
111
112 /* Function: allocate_mx()
113 *
114 * Purpose: Malloc space for the score matrices.
115 * amx and atr are indexed as j, i, y.
116 * bmx and btr are indexed as k, j, i.
117 * In the two sequence dimensions j, i they are
118 * diagonal (+1 off diagonal) matrices with
119 * rows j = 0..N, i = 1..j+1.
120 * In the node dimension k bmx and btr are k = 0..M.
121 * In the state dimension y amx and atr are y = 0..numstates.
122 *
123 * Args: icm - the int, log-odds, state-based model
124 * statenum - number of states in model
125 * window - length of scanning window
126 * ret_amx - RETURN: main score matrix
127 * ret_bmx - RETURN: BEGIN score matrix
128 *
129 * Return: Ptr to allocated scoring matrix, or
130 * dies and exits.
131 */
132 static int
allocate_mx(struct istate_s * icm,int statenum,int window,int **** ret_amx,int **** ret_bmx)133 allocate_mx(struct istate_s *icm,
134 int statenum,
135 int window,
136 int ****ret_amx,
137 int ****ret_bmx)
138
139 {
140 int ***amx;
141 int ***bmx;
142 int diag, j, y;
143
144 /* Main matrix, amx: fastest varying index is y (j,i,y)
145 * we only keep two rows for j, 0 and 1.
146 */
147 /* malloc for j = 0..1 rows */
148 if ((amx = (int ***) malloc (2 * sizeof(int **))) == NULL)
149 Die("Memory allocation error in %s line %d", __FILE__, __LINE__);
150
151 for (j = 0; j <= 1; j++) /* loop over rows j = 0..1 */
152 {
153 /* malloc for diag = 0..window cols */
154 if ((amx[j] = (int **) malloc ((window + 1) * sizeof(int *))) == NULL)
155 Die("Memory allocation error in %s line %d", __FILE__, __LINE__);
156
157 /* loop over cols diag = 0..window */
158 for (diag = 0; diag <= window; diag++)
159 /* malloc for y = 0..statenum-1 decks */
160 if ((amx[j][diag] = (int *) malloc ((statenum) * sizeof (int ))) == NULL)
161 Die("Memory allocation error in %s line %d", __FILE__, __LINE__);
162 }
163
164
165 /* B auxiliary matrix: fastest varying index is diag (y,j,diag)
166 * bmx keeps score decks for BEGIN states
167 */
168 /* 0..statenum-1 decks */
169 if ((bmx = (int ***) malloc (statenum * sizeof(int **))) == NULL)
170 Die("Memory allocation error in %s line %d", __FILE__, __LINE__);
171
172 for (y = 0; y < statenum; y++)
173 {
174 bmx[y] = NULL;
175 /* we keep score info for BEGIN states */
176 if (icm[y].statetype == uBEGIN_ST)
177 {
178 /* j= 0..window-1 rows */
179 if ((bmx[y] = (int **) malloc ((window) * sizeof(int *))) == NULL)
180 Die("Memory allocation error in %s line %d", __FILE__, __LINE__);
181 /* diff = 0..window columns */
182 for (j = 0; j < window; j++)
183 if ((bmx[y][j] = (int *) malloc ((window+1) * sizeof(int ))) == NULL)
184 Die("Memory allocation error in %s line %d", __FILE__, __LINE__);
185 }
186 }
187 *ret_amx = amx;
188 *ret_bmx = bmx;
189 return 1;
190 }
191
192
193
194 /* Function: free_mx()
195 *
196 * Purpose: Free the space allocated to the scoring and traceback matrices.
197 * Precisely mirrors the allocations above in allocate_cvmx().
198 *
199 * Return: (void)
200 */
201 static void
free_mx(int *** amx,int *** bmx,int statenum,int window)202 free_mx(int ***amx,
203 int ***bmx,
204 int statenum,
205 int window)
206 {
207 int diag, j, y;
208
209 /* Free the main matrix, amx:
210 * amx[j][i][y] = [0..1] [0..window] [0..statenum-1]
211 */
212 for (j = 0; j <= 1; j++)
213 {
214 for (diag = 0; diag <= window; diag++)
215 free(amx[j][diag]);
216 free(amx[j]);
217 }
218 free(amx);
219
220 /* Free the auxiliary matrix, bmx
221 * bmx[y][j][i] = [0..statenum-1] [0..window] [0..window]
222 */
223 for (y = 0; y < statenum; y++)
224 {
225 if (bmx[y] != NULL)
226 {
227 for (j = 0; j < window; j++)
228 free(bmx[y][j]);
229 free(bmx[y]);
230 }
231 }
232 free(bmx);
233 }
234
235
236
237 /* Function: init_mx()
238 *
239 * Purpose: Initialization of the scoring matrices. We initialize the off-diagonal,
240 * the diagonal, and the "floor" (end states) of the cube.
241 *
242 * Return: 1 on success, 0 on failure.
243 */
244 static int
init_mx(struct istate_s * icm,int statenum,int window,int *** amx,int *** bmx)245 init_mx(struct istate_s *icm, /* integer model */
246 int statenum, /* number of states in icm */
247 int window, /* size of scanning window on sequence */
248 int ***amx,
249 int ***bmx)
250 {
251 int diag, j, y; /* counters for indices over the cvmx */
252 int ynext; /* index of next state k+1 */
253 int *beam; /* z-axis vector of numbers in amx */
254
255 /* Init the whole amx to -Infinity. We do this with memcpy, trying
256 * to be fast. We fill in j=0,diag=0 by hand, then memcpy() the other
257 * columns.
258 */
259 for (y = 0; y < statenum; y++)
260 amx[0][0][y] = amx[1][0][y] = NEGINFINITY;
261 for (diag = 1; diag <= window; diag++)
262 {
263 memcpy(amx[0][diag], amx[0][0], statenum * sizeof(int));
264 memcpy(amx[1][diag], amx[0][0], statenum * sizeof(int));
265 }
266
267 /* Init the whole bmx to -Inf. We know state 0 is a begin (it's ROOT), so we
268 * start there, and memcpy rows as needed.
269 */
270 for (diag = 0; diag <= window; diag++)
271 bmx[0][0][diag] = NEGINFINITY;
272 for (j = 1; j < window; j++)
273 memcpy(bmx[0][j], bmx[0][0], (window+1) * sizeof(int));
274
275 for (y = 1; y < statenum; y++)
276 if (bmx[y] != NULL)
277 for (j = 0; j < window; j++)
278 memcpy(bmx[y][j], bmx[0][0], (window+1) * sizeof(int));
279
280 /* Init the off-diagonal (j = 0..window-1; diag == 0) with -log P scores.
281 * End state = 0;
282 * del, bifurc states are calc'ed
283 * begin states same as del's
284 * THIS IS WASTEFUL AND SHOULD BE CHANGED.
285 */
286 for (j = 0; j < window; j++)
287 for (y = statenum-1; y >= 0; y--)
288 {
289 /* Set the alignment of END states to the off-diagonal (diag = 0)
290 * to be zero, and never touch them again.
291 */
292 if (icm[y].statetype == uEND_ST)
293 amx[j%2][0][y] = 0;
294
295 else if (icm[y].statetype == uBIFURC_ST)
296 amx[j%2][0][y] = bmx[icm[y].tmx[0]][j][0] + bmx[icm[y].tmx[1]][j][0];
297
298 else if (icm[y].statetype == uDEL_ST || icm[y].statetype == uBEGIN_ST)
299 {
300 /* only calc DEL-DEL and BEGIN-DEL transitions. Since
301 * we optimized the state transition tables, removing
302 * the unused ones, we don't know where the number
303 * for "to DEL" is! But we can find it, because it'll
304 * be the connection to a non-infinite score */
305 beam = amx[j%2][0] + y + icm[y].offset;
306 for (ynext = 0; ynext < icm[y].connectnum; ynext++)
307 {
308 if (*beam != NEGINFINITY)
309 amx[j%2][0][y] = *beam + icm[y].tmx[ynext];
310 beam++;
311 }
312 }
313 /* make a copy into bmx if y is a BEGIN */
314 if (icm[y].statetype == uBEGIN_ST)
315 bmx[y][j][0] = amx[j%2][0][y];
316 }
317
318 return 1;
319 }
320
321
322
323 /* Function: recurse_mx()
324 *
325 * Purpose: Carry out the fill stage of the dynamic programming
326 * algorithm. After each j row is filled in, check the score
327 * of best full alignment ending at this row; if greater
328 * than threshold (ithresh), report it.
329 *
330 * Returns: 1 on success, 0 on failure.
331 */
332 static int
recurse_mx(struct istate_s * icm,int statenum,int * minb,int * maxb,char * seq,int seqlen,int window,int *** amx,int *** bmx,int ithresh,int (* gotone_f)(int,int,double))333 recurse_mx(struct istate_s *icm, /* integer, state-form model */
334 int statenum, /* number of states in icm */
335 int *minb,
336 int *maxb,
337 char *seq, /* sequence, 1..seqlen */
338 int seqlen, /* length of seq */
339 int window, /* length of scanning window on seq */
340 int ***amx, /* main scoring matrix */
341 int ***bmx, /* bifurc scoring matrix */
342 int ithresh, /* reporting threshold */
343 int (*gotone_f)(int, int, double))
344 {
345 int i, j, y; /* indices for 3 dimensions */
346 int aj; /* 0 or 1, index for j in A matrix */
347 int bj; /* 0..window-1, index for j in B matrix */
348 int diff; /* loop counter for difference: diff = j-i + 1 */
349 int symi, symj; /* symbol indices for seq[i], seq[j] */
350 int sc; /* tmp for a score */
351 int ynext; /* index of next state y */
352 int bestdiff, bestscore;
353
354 int *beam; /* ptr to a beam (z-axis vector) */
355 int leftdiff; /* diff coord of BEGIN_L of a bifurc */
356 int leftj; /* j coord of BEGIN_L of a bifurc */
357 int **left_p; /* pointer into whole 2D deck of BEGINL's of a bifurc */
358 int *right_p; /* ptr into row of BEGIN_R's of a bifurc */
359 int *scp; /* score pointer: ptr into beam of scores being calc'ed */
360 struct istate_s *st; /* state pointer: ptr at current state in icm */
361 int *tmx;
362 int emitsc;
363
364 for (j = 1; j <= seqlen; j++)
365 {
366 aj = j % 2; /* 0 or 1 index in amx */
367 bj = j % window; /* 0..window-1 index in bmx */
368 symj = SymbolIndex(seq[j]);
369
370 for (y = statenum-1; y >= 0; y--)
371 {
372 st = &icm[y];
373 for (diff = minb[y]; diff <= maxb[y] && diff <= j; diff++)
374 {
375 i = j - diff + 1;
376 symi = SymbolIndex(seq[i]);
377 scp = &amx[aj][diff][y];
378
379 if (st->statetype != uBIFURC_ST) /* a normal (non-BIFURC) state */
380 {
381 /* Connect the "beam" pointer to the appropriate
382 * starting place in the ynext scores we're connecting
383 * y to
384 */
385 switch (st->statetype) {
386 case uBEGIN_ST:
387 case uDEL_ST:
388 beam = amx[aj][diff];
389 emitsc = 0;
390 break;
391 case uMATP_ST: /* !aj toggles from 0 to 1 and vice versa */
392 if (diff == 1) continue;
393 beam = amx[!aj][diff-2];
394 emitsc = st->emit[symi * ALPHASIZE + symj];
395 break;
396 case uMATR_ST:
397 case uINSR_ST:
398 beam = amx[!aj][diff-1];
399 emitsc = st->emit[symj];
400 break;
401 case uMATL_ST:
402 case uINSL_ST:
403 beam = amx[aj][diff-1];
404 emitsc = st->emit[symi];
405 break;
406 case uEND_ST:
407 continue;
408 default: Die("no such state type %d", st->statetype);
409 }
410 beam += y + st->offset;
411 tmx = st->tmx;
412
413 /* Init for ynext == 0 case
414 */
415 *scp = *beam + *tmx;
416
417 /* Calculate remaining cases
418 */
419 for (ynext = 1; ynext < st->connectnum; ynext++)
420 {
421 beam++;
422 tmx++;
423 if (*beam > *scp)
424 {
425 sc = *beam + *tmx;
426 if (sc > *scp)
427 *scp = sc;
428 }
429 }
430
431 /* Add emission scores now
432 */
433 *scp += emitsc;
434
435 /* Make a copy into bmx, btr if necessary
436 */
437 if (st->statetype == uBEGIN_ST)
438 bmx[y][bj][diff] = *scp;
439 } /* end block of normal state stuff */
440
441 else /* a BIFURC state */
442 {
443 leftdiff = diff;
444 leftj = bj;
445 right_p = bmx[st->tmx[1]][leftj];
446 left_p = bmx[st->tmx[0]];
447
448 /* init w/ case that left branch emits it all */
449 *scp = left_p[leftj][leftdiff] + *right_p;
450 while (leftdiff > 0)
451 {
452 leftdiff--;
453 leftj = leftj ? leftj-1 : window-1; /* scan window wraparound */
454 right_p++;
455
456 sc = left_p[leftj][leftdiff] + *right_p;
457 if (sc > *scp)
458 *scp = sc;
459 }
460 }
461
462 } /* end loop over states */
463 } /* end loop over diff */
464
465 /* We've completed a row. Now we can examine the scores in diff,
466 * aj, ROOT_ST to decide whether to report this row. If we do,
467 * we report the 1..seqlen i, j coords of the matching subsequence
468 * in seq, as well as the score converted to double-precision bits.
469 */
470 bestdiff = 1;
471 bestscore = bmx[0][bj][1];
472 for (diff = 2; diff <= window; diff++)
473 if (bmx[0][bj][diff] > bestscore)
474 {
475 bestscore = bmx[0][bj][diff];
476 bestdiff = diff;
477 }
478 if (bestscore > ithresh)
479 if (! ReportScanHit(j - bestdiff + 1, j, (double)(bestscore / INTPRECISION), gotone_f))
480 Warn("caller ignored report of a match!");
481 } /* end loop over j */
482
483 return 1;
484 }
485
486