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