1 /* Optimal accuracy alignment; VMX version.
2 *
3 * Contents:
4 * 1. Optimal accuracy alignment, DP fill
5 * 2. OA traceback
6 * 3. Benchmark driver
7 * 4. Unit tests
8 * 5. Test driver
9 * 6. Example
10 *
11 * SRE, Mon Aug 18 20:01:01 2008 [Casa de Gatos]
12 */
13 #include "p7_config.h"
14
15 #include <float.h>
16
17 #ifndef __APPLE_ALTIVEC__
18 #include <altivec.h>
19 #endif
20
21 #include "easel.h"
22 #include "esl_vmx.h"
23 #include "esl_vectorops.h"
24
25 #include "hmmer.h"
26
27
28 /*****************************************************************
29 * 1. Optimal accuracy alignment, DP fill
30 *****************************************************************/
31
32 /* Function: p7_OptimalAccuracy()
33 * Synopsis: DP fill of an optimal accuracy alignment calculation.
34 * Incept: SRE, Mon Aug 18 11:04:48 2008 [Janelia]
35 *
36 * Purpose: Calculates the fill step of the optimal accuracy decoding
37 * algorithm \citep{Kall05}.
38 *
39 * Caller provides the posterior decoding matrix <pp>,
40 * which was calculated by Forward/Backward on a target sequence
41 * of length <pp->L> using the query model <om>.
42 *
43 * Caller also provides a DP matrix <ox>, allocated for a full
44 * <om->M> by <L> comparison. The routine fills this in
45 * with OA scores.
46 *
47 * Args: gm - query profile
48 * pp - posterior decoding matrix created by <p7_GPosteriorDecoding()>
49 * gx - RESULT: caller provided DP matrix for <gm->M> by <L>
50 * ret_e - RETURN: expected number of correctly decoded positions
51 *
52 * Returns: <eslOK> on success, and <*ret_e> contains the final OA
53 * score, which is the expected number of correctly decoded
54 * positions in the target sequence (up to <L>).
55 *
56 * Throws: (no abnormal error conditions)
57 */
58 int
p7_OptimalAccuracy(const P7_OPROFILE * om,const P7_OMX * pp,P7_OMX * ox,float * ret_e)59 p7_OptimalAccuracy(const P7_OPROFILE *om, const P7_OMX *pp, P7_OMX *ox, float *ret_e)
60 {
61 vector float mpv, dpv, ipv; /* previous row values */
62 vector float sv; /* temp storage of 1 curr row value in progress */
63 vector float xEv; /* E state: keeps max for Mk->E as we go */
64 vector float xBv; /* B state: splatted vector of B[i-1] for B->Mk calculations */
65 vector float dcv;
66 float *xmx = ox->xmx;
67 vector float *dpc = ox->dpf[0]; /* current row, for use in {MDI}MO(dpp,q) access macro */
68 vector float *dpp; /* previous row, for use in {MDI}MO(dpp,q) access macro */
69 vector float *ppp; /* quads in the <pp> posterior probability matrix */
70 vector float *tp; /* quads in the <om->tfv> transition scores */
71 vector float zerov;
72 vector float infv;
73 int M = om->M;
74 int Q = p7O_NQF(M);
75 int q;
76 int j;
77 int i;
78 float t1, t2;
79
80 zerov = (vector float) vec_splat_u32(0);
81 infv = esl_vmx_set_float(-eslINFINITY);
82
83 ox->M = om->M;
84 ox->L = pp->L;
85 for (q = 0; q < Q; q++) MMO(dpc, q) = IMO(dpc,q) = DMO(dpc,q) = infv;
86 XMXo(0, p7X_E) = -eslINFINITY;
87 XMXo(0, p7X_N) = 0.;
88 XMXo(0, p7X_J) = -eslINFINITY;
89 XMXo(0, p7X_B) = 0.;
90 XMXo(0, p7X_C) = -eslINFINITY;
91
92 for (i = 1; i <= pp->L; i++)
93 {
94 dpp = dpc; /* previous DP row in OA matrix */
95 dpc = ox->dpf[i]; /* current DP row in OA matrix */
96 ppp = pp->dpf[i]; /* current row in the posterior probabilities per position */
97 tp = om->tfv; /* transition probabilities */
98 dcv = infv;
99 xEv = infv;
100 xBv = esl_vmx_set_float(XMXo(i-1, p7X_B));
101
102 mpv = vec_sld(infv, MMO(dpp,Q-1), 12); /* Right shifts by 4 bytes. 4,8,12,x becomes x,4,8,12. */
103 dpv = vec_sld(infv, DMO(dpp,Q-1), 12);
104 ipv = vec_sld(infv, IMO(dpp,Q-1), 12);
105 for (q = 0; q < Q; q++)
106 {
107 sv = vec_and(vec_cmpgt(*tp, zerov), xBv); tp++;
108 sv = vec_max(sv, vec_and(vec_cmpgt(*tp, zerov), mpv)); tp++;
109 sv = vec_max(sv, vec_and(vec_cmpgt(*tp, zerov), ipv)); tp++;
110 sv = vec_max(sv, vec_and(vec_cmpgt(*tp, zerov), dpv)); tp++;
111 sv = vec_add(sv, *ppp); ppp += 2;
112 xEv = vec_max(xEv, sv);
113
114 mpv = MMO(dpp,q);
115 dpv = DMO(dpp,q);
116 ipv = IMO(dpp,q);
117
118 MMO(dpc,q) = sv;
119 DMO(dpc,q) = dcv;
120
121 dcv = vec_and(vec_cmpgt(*tp, zerov), sv); tp++;
122
123 sv = vec_and(vec_cmpgt(*tp, zerov), mpv); tp++;
124 sv = vec_max(sv, vec_and(vec_cmpgt(*tp, zerov), ipv)); tp++;
125 IMO(dpc,q) = vec_add(sv, *ppp); ppp++;
126 }
127
128 /* dcv has carried through from end of q loop above; store it
129 * in first pass, we add M->D and D->D path into DMX
130 */
131 dcv = vec_sld(infv, dcv, 12);
132 tp = om->tfv + 7*Q; /* set tp to start of the DD's */
133 for (q = 0; q < Q; q++)
134 {
135 DMO(dpc, q) = vec_max(dcv, DMO(dpc, q));
136 dcv = vec_and(vec_cmpgt(*tp, zerov), DMO(dpc,q)); tp++;
137 }
138
139 /* fully serialized D->D; can optimize later */
140 for (j = 1; j < 4; j++)
141 {
142 dcv = vec_sld(infv, dcv, 12);
143 tp = om->tfv + 7*Q;
144 for (q = 0; q < Q; q++)
145 {
146 DMO(dpc, q) = vec_max(dcv, DMO(dpc, q));
147 dcv = vec_and(vec_cmpgt(*tp, zerov), dcv); tp++;
148 }
149 }
150
151 /* D->E paths */
152 for (q = 0; q < Q; q++) xEv = vec_max(xEv, DMO(dpc,q));
153
154 /* Specials */
155 XMXo(i,p7X_E) = esl_vmx_hmax_float(xEv);
156
157 t1 = ( (om->xf[p7O_J][p7O_LOOP] == 0.0) ? 0.0 : ox->xmx[(i-1)*p7X_NXCELLS+p7X_J] + pp->xmx[i*p7X_NXCELLS+p7X_J]);
158 t2 = ( (om->xf[p7O_E][p7O_LOOP] == 0.0) ? 0.0 : ox->xmx[ i *p7X_NXCELLS+p7X_E]);
159 ox->xmx[i*p7X_NXCELLS+p7X_J] = ESL_MAX(t1, t2);
160
161 t1 = ( (om->xf[p7O_C][p7O_LOOP] == 0.0) ? 0.0 : ox->xmx[(i-1)*p7X_NXCELLS+p7X_C] + pp->xmx[i*p7X_NXCELLS+p7X_C]);
162 t2 = ( (om->xf[p7O_E][p7O_MOVE] == 0.0) ? 0.0 : ox->xmx[ i *p7X_NXCELLS+p7X_E]);
163 ox->xmx[i*p7X_NXCELLS+p7X_C] = ESL_MAX(t1, t2);
164
165 ox->xmx[i*p7X_NXCELLS+p7X_N] = ((om->xf[p7O_N][p7O_LOOP] == 0.0) ? 0.0 : ox->xmx[(i-1)*p7X_NXCELLS+p7X_N] + pp->xmx[i*p7X_NXCELLS+p7X_N]);
166
167 t1 = ( (om->xf[p7O_N][p7O_MOVE] == 0.0) ? 0.0 : ox->xmx[i*p7X_NXCELLS+p7X_N]);
168 t2 = ( (om->xf[p7O_J][p7O_MOVE] == 0.0) ? 0.0 : ox->xmx[i*p7X_NXCELLS+p7X_J]);
169 ox->xmx[i*p7X_NXCELLS+p7X_B] = ESL_MAX(t1, t2);
170 }
171
172 *ret_e = ox->xmx[pp->L*p7X_NXCELLS+p7X_C];
173 return eslOK;
174 }
175 /*------------------- end, OA DP fill ---------------------------*/
176
177
178
179
180
181 /*****************************************************************
182 * 2. OA traceback
183 *****************************************************************/
184
185 static inline float get_postprob(const P7_OMX *pp, int scur, int sprv, int k, int i);
186
187 static inline int select_m(const P7_OPROFILE *om, const P7_OMX *ox, int i, int k);
188 static inline int select_d(const P7_OPROFILE *om, const P7_OMX *ox, int i, int k);
189 static inline int select_i(const P7_OPROFILE *om, const P7_OMX *ox, int i, int k);
190 static inline int select_n(int i);
191 static inline int select_c(const P7_OPROFILE *om, const P7_OMX *pp, const P7_OMX *ox, int i);
192 static inline int select_j(const P7_OPROFILE *om, const P7_OMX *pp, const P7_OMX *ox, int i);
193 static inline int select_e(const P7_OPROFILE *om, const P7_OMX *ox, int i, int *ret_k);
194 static inline int select_b(const P7_OPROFILE *om, const P7_OMX *ox, int i);
195
196
197 /* Function: p7_OATrace()
198 * Synopsis: Optimal accuracy decoding: traceback.
199 * Incept: SRE, Mon Aug 18 13:53:33 2008 [Janelia]
200 *
201 * Purpose: The traceback stage of the optimal accuracy decoding algorithm
202 * \citep{Kall05}.
203 *
204 * Caller provides the OA DP matrix <ox> that was just
205 * calculated by <p7_OptimalAccuracyDP()>, as well as the
206 * posterior decoding matrix <pp>, which was calculated by
207 * Forward/Backward on a target sequence using the query
208 * model <gm>. Because the calculation depends only on
209 * <pp>, the target sequence itself need not be provided.
210 *
211 * The resulting optimal accuracy decoding traceback is put
212 * in a caller-provided traceback structure <tr>, which the
213 * caller has allocated for optional posterior probability
214 * annotation on residues (with <p7_trace_CreateWithPP()>,
215 * generally). This structure will be reallocated
216 * internally if necessary.
217 *
218 * Args: om - profile
219 * pp - posterior probability matrix
220 * ox - OA matrix to trace, LxM
221 * tr - storage for the recovered traceback
222 *
223 * Returns: <eslOK> on success.
224 *
225 * Throws: <eslEMEM> on allocation error.
226 * <eslEINVAL> if the trace <tr> isn't empty (needs to be Reuse()'d).
227 */
228 int
p7_OATrace(const P7_OPROFILE * om,const P7_OMX * pp,const P7_OMX * ox,P7_TRACE * tr)229 p7_OATrace(const P7_OPROFILE *om, const P7_OMX *pp, const P7_OMX *ox, P7_TRACE *tr)
230 {
231 int i = ox->L; /* position in sequence 1..L */
232 int k = 0; /* position in model 1..M */
233 int s0, s1; /* choice of a state */
234 float postprob;
235 int status;
236
237 if (tr->N != 0) ESL_EXCEPTION(eslEINVAL, "trace not empty; needs to be Reuse()'d?");
238
239 if ((status = p7_trace_AppendWithPP(tr, p7T_T, k, i, 0.0)) != eslOK) return status;
240 if ((status = p7_trace_AppendWithPP(tr, p7T_C, k, i, 0.0)) != eslOK) return status;
241
242 s0 = tr->st[tr->N-1];
243 while (s0 != p7T_S)
244 {
245 switch (s0) {
246 case p7T_M: s1 = select_m(om, ox, i, k); k--; i--; break;
247 case p7T_D: s1 = select_d(om, ox, i, k); k--; break;
248 case p7T_I: s1 = select_i(om, ox, i, k); i--; break;
249 case p7T_N: s1 = select_n(i); break;
250 case p7T_C: s1 = select_c(om, pp, ox, i); break;
251 case p7T_J: s1 = select_j(om, pp, ox, i); break;
252 case p7T_E: s1 = select_e(om, ox, i, &k); break;
253 case p7T_B: s1 = select_b(om, ox, i); break;
254 default: ESL_EXCEPTION(eslEINVAL, "bogus state in traceback");
255 }
256 if (s1 == -1) ESL_EXCEPTION(eslEINVAL, "OA traceback choice failed");
257
258 postprob = get_postprob(pp, s1, s0, k, i);
259 if ((status = p7_trace_AppendWithPP(tr, s1, k, i, postprob)) != eslOK) return status;
260
261 if ( (s1 == p7T_N || s1 == p7T_J || s1 == p7T_C) && s1 == s0) i--;
262 s0 = s1;
263 } /* end traceback, at S state */
264 tr->M = om->M;
265 tr->L = ox->L;
266 return p7_trace_Reverse(tr);
267 }
268
269 static inline float
get_postprob(const P7_OMX * pp,int scur,int sprv,int k,int i)270 get_postprob(const P7_OMX *pp, int scur, int sprv, int k, int i)
271 {
272 int Q = p7O_NQF(pp->M);
273 int q = (k-1) % Q; /* (q,r) is position of the current DP cell M(i,k) */
274 int r = (k-1) / Q;
275 union { vector float v; float p[4]; } u;
276
277 switch (scur) {
278 case p7T_M: u.v = MMO(pp->dpf[i], q); return u.p[r];
279 case p7T_I: u.v = IMO(pp->dpf[i], q); return u.p[r];
280 case p7T_N: if (sprv == scur) return pp->xmx[i*p7X_NXCELLS+p7X_N];
281 case p7T_C: if (sprv == scur) return pp->xmx[i*p7X_NXCELLS+p7X_C];
282 case p7T_J: if (sprv == scur) return pp->xmx[i*p7X_NXCELLS+p7X_J];
283 default: return 0.0;
284 }
285 }
286
287 /* M(i,k) is reached from B(i-1), M(i-1,k-1), D(i-1,k-1), or I(i-1,k-1). */
288 static inline int
select_m(const P7_OPROFILE * om,const P7_OMX * ox,int i,int k)289 select_m(const P7_OPROFILE *om, const P7_OMX *ox, int i, int k)
290 {
291 int Q = p7O_NQF(ox->M);
292 int q = (k-1) % Q; /* (q,r) is position of the current DP cell M(i,k) */
293 int r = (k-1) / Q;
294 vector float *tp = om->tfv + 7*q; /* *tp now at start of transitions to cur cell M(i,k) */
295 vector float xBv;
296 vector float zerov;
297 vector float mpv, dpv, ipv;
298 union { vector float v; float p[4]; } u, tv;
299 float path[4];
300 int state[4] = { p7T_M, p7T_I, p7T_D, p7T_B };
301
302 xBv = esl_vmx_set_float(ox->xmx[(i-1)*p7X_NXCELLS+p7X_B]);
303 zerov = (vector float) vec_splat_u32(0);
304
305 if (q > 0) {
306 mpv = ox->dpf[i-1][(q-1)*3 + p7X_M];
307 dpv = ox->dpf[i-1][(q-1)*3 + p7X_D];
308 ipv = ox->dpf[i-1][(q-1)*3 + p7X_I];
309 } else {
310 mpv = vec_sld(zerov, ox->dpf[i-1][(Q-1)*3 + p7X_M], 12);
311 dpv = vec_sld(zerov, ox->dpf[i-1][(Q-1)*3 + p7X_D], 12);
312 ipv = vec_sld(zerov, ox->dpf[i-1][(Q-1)*3 + p7X_I], 12);
313 }
314
315 /* paths are numbered so that most desirable choice in case of tie is first. */
316 u.v = xBv; tv.v = *tp; path[3] = ((tv.p[r] == 0.0) ? -eslINFINITY : u.p[r]); tp++;
317 u.v = mpv; tv.v = *tp; path[0] = ((tv.p[r] == 0.0) ? -eslINFINITY : u.p[r]); tp++;
318 u.v = ipv; tv.v = *tp; path[1] = ((tv.p[r] == 0.0) ? -eslINFINITY : u.p[r]); tp++;
319 u.v = dpv; tv.v = *tp; path[2] = ((tv.p[r] == 0.0) ? -eslINFINITY : u.p[r]);
320 return state[esl_vec_FArgMax(path, 4)];
321 }
322
323
324 /* D(i,k) is reached from M(i, k-1) or D(i,k-1). */
325 static inline int
select_d(const P7_OPROFILE * om,const P7_OMX * ox,int i,int k)326 select_d(const P7_OPROFILE *om, const P7_OMX *ox, int i, int k)
327 {
328 int Q = p7O_NQF(ox->M);
329 int q = (k-1) % Q; /* (q,r) is position of the current DP cell D(i,k) */
330 int r = (k-1) / Q;
331 vector float zerov;
332 union { vector float v; float p[4]; } mpv, dpv, tmdv, tddv;
333 float path[2];
334
335 zerov = (vector float) vec_splat_u32(0);
336
337 if (q > 0) {
338 mpv.v = ox->dpf[i][(q-1)*3 + p7X_M];
339 dpv.v = ox->dpf[i][(q-1)*3 + p7X_D];
340 tmdv.v = om->tfv[7*(q-1) + p7O_MD];
341 tddv.v = om->tfv[7*Q + (q-1)];
342 } else {
343 mpv.v = vec_sld(zerov, ox->dpf[i][(Q-1)*3 + p7X_M], 12);
344 dpv.v = vec_sld(zerov, ox->dpf[i][(Q-1)*3 + p7X_D], 12);
345 tmdv.v = vec_sld(zerov, om->tfv[7*(Q-1) + p7O_MD], 12);
346 tddv.v = vec_sld(zerov, om->tfv[8*Q-1], 12);
347 }
348
349 path[0] = ((tmdv.p[r] == 0.0) ? -eslINFINITY : mpv.p[r]);
350 path[1] = ((tddv.p[r] == 0.0) ? -eslINFINITY : dpv.p[r]);
351 return ((path[0] >= path[1]) ? p7T_M : p7T_D);
352 }
353
354 /* I(i,k) is reached from M(i-1, k) or I(i-1,k). */
355 static inline int
select_i(const P7_OPROFILE * om,const P7_OMX * ox,int i,int k)356 select_i(const P7_OPROFILE *om, const P7_OMX *ox, int i, int k)
357 {
358 int Q = p7O_NQF(ox->M);
359 int q = (k-1) % Q; /* (q,r) is position of the current DP cell D(i,k) */
360 int r = (k-1) / Q;
361 vector float *tp = om->tfv + 7*q + p7O_MI;
362 union { vector float v; float p[4]; } tv, mpv, ipv;
363 float path[2];
364
365 mpv.v = ox->dpf[i-1][q*3 + p7X_M]; tv.v = *tp; path[0] = ((tv.p[r] == 0.0) ? -eslINFINITY : mpv.p[r]); tp++;
366 ipv.v = ox->dpf[i-1][q*3 + p7X_I]; tv.v = *tp; path[1] = ((tv.p[r] == 0.0) ? -eslINFINITY : ipv.p[r]);
367 return ((path[0] >= path[1]) ? p7T_M : p7T_I);
368 }
369
370 /* N(i) must come from N(i-1) for i>0; else it comes from S */
371 static inline int
select_n(int i)372 select_n(int i)
373 {
374 return ((i==0) ? p7T_S : p7T_N);
375 }
376
377 /* C(i) is reached from E(i) or C(i-1). */
378 static inline int
select_c(const P7_OPROFILE * om,const P7_OMX * pp,const P7_OMX * ox,int i)379 select_c(const P7_OPROFILE *om, const P7_OMX *pp, const P7_OMX *ox, int i)
380 {
381 float path[2];
382 path[0] = ( (om->xf[p7O_C][p7O_LOOP] == 0.0) ? -eslINFINITY : ox->xmx[(i-1)*p7X_NXCELLS+p7X_C] + pp->xmx[i*p7X_NXCELLS+p7X_C]);
383 path[1] = ( (om->xf[p7O_E][p7O_MOVE] == 0.0) ? -eslINFINITY : ox->xmx[ i *p7X_NXCELLS+p7X_E]);
384 return ((path[0] > path[1]) ? p7T_C : p7T_E);
385 }
386
387 /* J(i) is reached from E(i) or J(i-1). */
388 static inline int
select_j(const P7_OPROFILE * om,const P7_OMX * pp,const P7_OMX * ox,int i)389 select_j(const P7_OPROFILE *om, const P7_OMX *pp, const P7_OMX *ox, int i)
390 {
391 float path[2];
392
393 path[0] = ( (om->xf[p7O_J][p7O_LOOP] == 0.0) ? -eslINFINITY : ox->xmx[(i-1)*p7X_NXCELLS+p7X_J] + pp->xmx[i*p7X_NXCELLS+p7X_J]);
394 path[1] = ( (om->xf[p7O_E][p7O_LOOP] == 0.0) ? -eslINFINITY : ox->xmx[ i *p7X_NXCELLS+p7X_E]);
395 return ((path[0] > path[1]) ? p7T_J : p7T_E);
396 }
397
398 /* E(i) is reached from any M(i, k=1..M) or D(i, k=2..M). */
399 /* This assumes all M_k->E, D_k->E are 1.0 */
400 static inline int
select_e(const P7_OPROFILE * om,const P7_OMX * ox,int i,int * ret_k)401 select_e(const P7_OPROFILE *om, const P7_OMX *ox, int i, int *ret_k)
402 {
403 int Q = p7O_NQF(ox->M);
404 vector float *dp = ox->dpf[i];
405 union { vector float v; float p[4]; } u;
406 float max = -eslINFINITY;
407 int smax, kmax;
408 int q,r;
409
410 /* precedence rules in case of ties here are a little tricky: M beats D: note the >= max! */
411 for (q = 0; q < Q; q++)
412 {
413 u.v = *dp; dp++; for (r = 0; r < 4; r++) if (u.p[r] >= max) { max = u.p[r]; smax = p7T_M; kmax = r*Q + q + 1; }
414 u.v = *dp; dp+=2; for (r = 0; r < 4; r++) if (u.p[r] > max) { max = u.p[r]; smax = p7T_D; kmax = r*Q + q + 1; }
415 }
416 *ret_k = kmax;
417 return smax;
418 }
419
420
421 /* B(i) is reached from N(i) or J(i). */
422 static inline int
select_b(const P7_OPROFILE * om,const P7_OMX * ox,int i)423 select_b(const P7_OPROFILE *om, const P7_OMX *ox, int i)
424 {
425 float path[2];
426
427 path[0] = ( (om->xf[p7O_N][p7O_MOVE] == 0.0) ? -eslINFINITY : ox->xmx[i*p7X_NXCELLS+p7X_N]);
428 path[1] = ( (om->xf[p7O_J][p7O_MOVE] == 0.0) ? -eslINFINITY : ox->xmx[i*p7X_NXCELLS+p7X_J]);
429 return ((path[0] > path[1]) ? p7T_N : p7T_J);
430 }
431 /*---------------------- end, OA traceback ----------------------*/
432
433
434
435
436 /*****************************************************************
437 * 3. Benchmark driver
438 *****************************************************************/
439 #ifdef p7OPTACC_BENCHMARK
440 /*
441 icc -O3 -static -o optacc_benchmark -I.. -L.. -I../../easel -L../../easel -Dp7OPTACC_BENCHMARK optacc.c -lhmmer -leasel -lm
442
443 ./optacc_benchmark <hmmfile> runs benchmark on optimal accuracy fill and trace
444 ./optacc_benchmark -c -N1 <hmmfile> compare scores of VMX version to generic impl
445 ./optacc_benchmark -x -N1 <hmmfile> test that scores match trusted implementation.
446
447 RRM_1 (M=72) Caudal_act (M=136) SMC_N (M=1151)
448 ----------------- ------------------ ---------------
449 20 Aug 08: 13.11u (110 Mc/s) 23.39u (116 Mc/s) 332.62u (69 Mc/s)
450
451 */
452 #include "p7_config.h"
453
454 #include "easel.h"
455 #include "esl_alphabet.h"
456 #include "esl_getopts.h"
457 #include "esl_random.h"
458 #include "esl_randomseq.h"
459 #include "esl_stopwatch.h"
460
461 #include "hmmer.h"
462 #include "impl_vmx.h"
463
464 static ESL_OPTIONS options[] = {
465 /* name type default env range toggles reqs incomp help docgroup*/
466 { "-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show brief help on version and usage", 0 },
467 { "-c", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, "-x", "compare scores to generic implementation (debug)", 0 },
468 { "-s", eslARG_INT, "42", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>", 0 },
469 { "-x", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, "-c", "equate scores to trusted implementation (debug)", 0 },
470 { "-L", eslARG_INT, "400", NULL, "n>0", NULL, NULL, NULL, "length of random target seqs", 0 },
471 { "-N", eslARG_INT, "50000", NULL, "n>0", NULL, NULL, NULL, "number of random target seqs", 0 },
472 { "--notrace", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "only benchmark the DP fill stage", 0 },
473 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
474 };
475 static char usage[] = "[-options] <hmmfile>";
476 static char banner[] = "benchmark driver for optimal accuracy alignment, VMX version";
477
478 int
main(int argc,char ** argv)479 main(int argc, char **argv)
480 {
481 ESL_GETOPTS *go = p7_CreateDefaultApp(options, 1, argc, argv, banner, usage);
482 char *hmmfile = esl_opt_GetArg(go, 1);
483 ESL_STOPWATCH *w = esl_stopwatch_Create();
484 ESL_RANDOMNESS *r = esl_randomness_CreateFast(esl_opt_GetInteger(go, "-s"));
485 ESL_ALPHABET *abc = NULL;
486 P7_HMMFILE *hfp = NULL;
487 P7_HMM *hmm = NULL;
488 P7_BG *bg = NULL;
489 P7_PROFILE *gm = NULL;
490 P7_OPROFILE *om = NULL;
491 P7_GMX *gx1 = NULL;
492 P7_GMX *gx2 = NULL;
493 P7_OMX *ox1 = NULL;
494 P7_OMX *ox2 = NULL;
495 P7_TRACE *tr = NULL;
496 int L = esl_opt_GetInteger(go, "-L");
497 int N = esl_opt_GetInteger(go, "-N");
498 ESL_DSQ *dsq = malloc(sizeof(ESL_DSQ) * (L+2));
499 int i;
500 float fsc, bsc, accscore;
501 float fsc_g, bsc_g, accscore_g;
502 double Mcs;
503
504 p7_FLogsumInit();
505
506 if (p7_hmmfile_OpenE(hmmfile, NULL, &hfp, NULL) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
507 if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
508
509 bg = p7_bg_Create(abc); p7_bg_SetLength(bg, L);
510 gm = p7_profile_Create(hmm->M, abc); p7_ProfileConfig(hmm, bg, gm, L, p7_LOCAL);
511 om = p7_oprofile_Create(gm->M, abc); p7_oprofile_Convert(gm, om);
512 p7_oprofile_ReconfigLength(om, L);
513
514 if (esl_opt_GetBoolean(go, "-x") && p7_FLogsumError(-0.4, -0.5) > 0.0001)
515 p7_Fail("-x here requires p7_Logsum() recompiled in slow exact mode");
516
517 ox1 = p7_omx_Create(gm->M, L, L);
518 ox2 = p7_omx_Create(gm->M, L, L);
519 tr = p7_trace_CreateWithPP();
520
521 esl_rsq_xfIID(r, bg->f, abc->K, L, dsq);
522 p7_Forward (dsq, L, om, ox1, &fsc);
523 p7_Backward(dsq, L, om, ox1, ox2, &bsc);
524 p7_Decoding(om, ox1, ox2, ox2);
525
526 esl_stopwatch_Start(w);
527 for (i = 0; i < N; i++)
528 {
529 p7_OptimalAccuracy(om, ox2, ox1, &accscore);
530
531 if (! esl_opt_GetBoolean(go, "--notrace"))
532 {
533 p7_OATrace(om, ox2, ox1, tr);
534 p7_trace_Reuse(tr);
535 }
536 }
537 esl_stopwatch_Stop(w);
538
539 Mcs = (double) N * (double) L * (double) gm->M * 1e-6 / (double) w->user;
540 esl_stopwatch_Display(stdout, w, "# CPU time: ");
541 printf("# M = %d\n", gm->M);
542 printf("# %.1f Mc/s\n", Mcs);
543
544 if (esl_opt_GetBoolean(go, "-c") || esl_opt_GetBoolean(go, "-x") )
545 {
546 gx1 = p7_gmx_Create(gm->M, L);
547 gx2 = p7_gmx_Create(gm->M, L);
548
549 p7_GForward (dsq, L, gm, gx1, &fsc_g);
550 p7_GBackward(dsq, L, gm, gx2, &bsc_g);
551 p7_GDecoding(gm, gx1, gx2, gx2);
552 p7_GOptimalAccuracy(gm, gx2, gx1, &accscore_g);
553
554 printf("generic: fwd=%8.4f bck=%8.4f acc=%8.4f\n", fsc_g, bsc_g, accscore_g);
555 printf("VMX: fwd=%8.4f bck=%8.4f acc=%8.4f\n", fsc, bsc, accscore);
556
557 p7_gmx_Destroy(gx1);
558 p7_gmx_Destroy(gx2);
559 }
560
561 free(dsq);
562 p7_omx_Destroy(ox1);
563 p7_omx_Destroy(ox2);
564 p7_trace_Destroy(tr);
565 p7_oprofile_Destroy(om);
566 p7_profile_Destroy(gm);
567 p7_bg_Destroy(bg);
568 p7_hmm_Destroy(hmm);
569 p7_hmmfile_Close(hfp);
570 esl_alphabet_Destroy(abc);
571 esl_stopwatch_Destroy(w);
572 esl_randomness_Destroy(r);
573 esl_getopts_Destroy(go);
574 return 0;
575 }
576 #endif /*p7OPTACC_BENCHMARK*/
577 /*---------------- end, benchmark driver ------------------------*/
578
579
580
581
582 /*****************************************************************
583 * 4. Unit tests
584 *****************************************************************/
585 #ifdef p7OPTACC_TESTDRIVE
586 #include "esl_alphabet.h"
587 #include "esl_getopts.h"
588 #include "esl_random.h"
589 /*
590 * 1. Compare accscore to GOptimalAccuracy().
591 * 2. Compare trace to GOATrace().
592 *
593 * Note: This test is subject to some expected noise and can fail
594 * for entirely innocent reasons. Generic Forward/Backward calculations with
595 * p7_GForward(), p7_GBackward() use coarse-grain table lookups to sum
596 * log probabilities, and sufficient roundoff error can accumulate to
597 * change the optimal accuracy traceback, causing this test to fail.
598 * So, if optacc_utest fails, before you go looking for bugs, first
599 * go to ../logsum.c, change the #ifdef to activate the slow/accurate
600 * version, recompile and rerun optacc_utest. If the failure goes away,
601 * you can ignore it. - SRE, Wed Dec 17 09:45:31 2008
602 */
603 static void
utest_optacc(ESL_GETOPTS * go,ESL_RANDOMNESS * r,ESL_ALPHABET * abc,P7_BG * bg,int M,int L,int N)604 utest_optacc(ESL_GETOPTS *go, ESL_RANDOMNESS *r, ESL_ALPHABET *abc, P7_BG *bg, int M, int L, int N)
605 {
606 char *msg = "optimal accuracy unit test failed";
607 P7_HMM *hmm = NULL;
608 P7_PROFILE *gm = NULL;
609 P7_OPROFILE *om = NULL;
610 ESL_SQ *sq = esl_sq_CreateDigital(abc);
611 P7_OMX *ox1 = p7_omx_Create(M, L, L);
612 P7_OMX *ox2 = p7_omx_Create(M, L, L);
613 P7_GMX *gx1 = p7_gmx_Create(M, L);
614 P7_GMX *gx2 = p7_gmx_Create(M, L);
615 P7_TRACE *tr = p7_trace_CreateWithPP();
616 P7_TRACE *trg = p7_trace_CreateWithPP();
617 P7_TRACE *tro = p7_trace_CreateWithPP();
618 float accscore_o;
619 float fsc, bsc, accscore;
620 float fsc_g, bsc_g, accscore_g, accscore_g2;
621 float pptol = 0.01;
622 float sctol = 0.001;
623 float gtol;
624
625 p7_FLogsumInit();
626 gtol = ( (p7_FLogsumError(-0.4, -0.5) > 0.0001) ? 0.1 : 0.001);
627
628 if (p7_oprofile_Sample(r, abc, bg, M, L, &hmm, &gm, &om)!= eslOK) esl_fatal(msg);
629 while (N--)
630 {
631 if (p7_ProfileEmit(r, hmm, gm, bg, sq, tro) != eslOK) esl_fatal(msg);
632
633 if (p7_omx_GrowTo(ox1, M, sq->n, sq->n) != eslOK) esl_fatal(msg);
634 if (p7_omx_GrowTo(ox2, M, sq->n, sq->n) != eslOK) esl_fatal(msg);
635 if (p7_gmx_GrowTo(gx1, M, sq->n) != eslOK) esl_fatal(msg);
636 if (p7_gmx_GrowTo(gx2, M, sq->n) != eslOK) esl_fatal(msg);
637
638 if (p7_Forward (sq->dsq, sq->n, om, ox1, &fsc) != eslOK) esl_fatal(msg);
639 if (p7_Backward(sq->dsq, sq->n, om, ox1, ox2, &bsc) != eslOK) esl_fatal(msg);
640 if (p7_Decoding(om, ox1, ox2, ox2) != eslOK) esl_fatal(msg);
641 if (p7_OptimalAccuracy(om, ox2, ox1, &accscore) != eslOK) esl_fatal(msg);
642
643 #if 0
644 p7_omx_FDeconvert(ox1, gx1);
645 p7_gmx_Dump(stdout, gx1, p7_DEFAULT);
646 p7_omx_FDeconvert(ox2, gx1);
647 p7_gmx_Dump(stdout, gx1, p7_DEFAULT);
648 #endif
649 if (p7_OATrace(om, ox2, ox1, tr) != eslOK) esl_fatal(msg);
650
651 if (p7_GForward (sq->dsq, sq->n, gm, gx1, &fsc_g) != eslOK) esl_fatal(msg);
652 if (p7_GBackward(sq->dsq, sq->n, gm, gx2, &bsc_g) != eslOK) esl_fatal(msg);
653
654 #if 0
655 p7_gmx_Dump(stdout, gx1, p7_DEFAULT); /* fwd */
656 p7_gmx_Dump(stdout, gx2, p7_DEFAULT); /* bck */
657 #endif
658
659 if (p7_GDecoding(gm, gx1, gx2, gx2) != eslOK) esl_fatal(msg);
660 if (p7_GOptimalAccuracy(gm, gx2, gx1, &accscore_g) != eslOK) esl_fatal(msg);
661
662 #if 0
663 p7_gmx_Dump(stdout, gx1, p7_DEFAULT); /* oa */
664 p7_gmx_Dump(stdout, gx2, p7_DEFAULT); /* pp */
665 #endif
666 if (p7_GOATrace(gm, gx2, gx1, trg) != eslOK) esl_fatal(msg);
667
668 if (p7_trace_SetPP(tro, gx2) != eslOK) esl_fatal(msg);
669
670 if (esl_opt_GetBoolean(go, "--traces"))
671 {
672 p7_trace_Dump(stdout, tro, gm, sq->dsq);
673 p7_trace_Dump(stdout, tr, gm, sq->dsq);
674 p7_trace_Dump(stdout, trg, gm, sq->dsq);
675 }
676
677 if (p7_trace_Validate(tr, abc, sq->dsq, NULL) != eslOK) esl_fatal(msg);
678 if (p7_trace_Validate(trg, abc, sq->dsq, NULL) != eslOK) esl_fatal(msg);
679 if (p7_trace_Compare(tr, trg, pptol) != eslOK) esl_fatal(msg);
680
681 accscore_o = p7_trace_GetExpectedAccuracy(tro); /* according to gx2; see p7_trace_SetPP() call above */
682 accscore_g2 = p7_trace_GetExpectedAccuracy(trg);
683
684 #if 0
685 printf("%f %f %f %f\n", accscore, accscore_g, accscore_g2, accscore_o);
686 #endif
687
688 if (esl_FCompare(fsc, bsc, sctol) != eslOK) esl_fatal(msg);
689 if (esl_FCompare(fsc_g, bsc_g, gtol) != eslOK) esl_fatal(msg);
690 if (esl_FCompare(fsc, fsc_g, gtol) != eslOK) esl_fatal(msg);
691 if (esl_FCompare(accscore, accscore_g, gtol) != eslOK) esl_fatal(msg);
692 if (esl_FCompare(accscore_g, accscore_g2, gtol) != eslOK) esl_fatal(msg);
693 if (accscore_g2 < accscore_o) esl_fatal(msg);
694 /* the above deserves explanation:
695 * - accscore_o is the accuracy of the originally emitted trace, according
696 * to the generic posterior decoding matrix <gx2>. This is a lower bound
697 * on the expected # of accurately aligned residues found by a DP
698 * optimization.
699 * - accscore is the accuracy found by the fast (vector) code DP implementation.
700 * - accscore_g is the accuracy found by the generic DP implementation.
701 * accscore and accscore_g should be nearly identical,
702 * within tolerance of roundoff error accumulation and
703 * the imprecision of Logsum() tables.
704 * - accscore_g2 is the accuracy of the traceback identified by the generic
705 * DP implementation. It should be identical (within order-of-evaluation
706 * roundoff error) to accscore_g.
707 *
708 * the "accscore_g2 < accscore_o" test is carefully contrived.
709 * accscore_o is a theoretical lower bound but because of fp error,
710 * accscore and (much more rarely) even accscore_g can exceed accscore_o.
711 * accscore_g2, however, is calculated with identical order of evaluation
712 * as accscore_o if the optimal trace does turn out to be identical to
713 * the originally emitted trace. It should be extremely unlikely (though
714 * not impossible) for accscore_o to exceed accscore_g2. (The DP algorithm
715 * would have to identify a trace that was different than the original trace,
716 * which the DP algorithm, by order-of-evaluation, assigned higher accuracy,
717 * but order-of-evaluation in traceback dependent code assigned lower accuracy.
718 * [xref J5/29]
719 */
720
721 esl_sq_Reuse(sq);
722 p7_trace_Reuse(tr);
723 p7_trace_Reuse(trg);
724 p7_trace_Reuse(tro);
725 }
726
727 p7_trace_Destroy(tro);
728 p7_trace_Destroy(trg);
729 p7_trace_Destroy(tr);
730 p7_gmx_Destroy(gx2);
731 p7_gmx_Destroy(gx1);
732 p7_omx_Destroy(ox2);
733 p7_omx_Destroy(ox1);
734 esl_sq_Destroy(sq);
735 p7_oprofile_Destroy(om);
736 p7_profile_Destroy(gm);
737 p7_hmm_Destroy(hmm);
738 }
739
740 #endif /*p7OPTACC_TESTDRIVE*/
741 /*------------------- end, unit tests ---------------------------*/
742
743
744
745
746 /*****************************************************************
747 * 5. Test driver
748 *****************************************************************/
749 #ifdef p7OPTACC_TESTDRIVE
750
751 /* Failures in this test are to be expected, if you change the defaults.
752 * The default RNG seed of 41 is very carefully chosen! Most seeds will
753 * cause this test to fail. (Only 13 and 41 work for seeds 1..50.)
754 *
755 * The generic fwd/bck algorithms work in log space and suffer from a
756 * small amount of imprecision, not only from the use of FLogsum()'s
757 * table-driven approximation, but even (apparently) inherent in log()
758 * and exp(). To minimize this, the generic decoding algorithm burns
759 * time renormalizing each row. Still, it's a problem. See notes at
760 * the header of utest_optacc() for more info.
761 *
762 * Another expected failure mode is when a fwd, bck nat score are close to
763 * 0.0; FCompare() can evaluate two close-to-zero numbers as "different"
764 * even if their absolute diff is negligible. (I suppose I could fix this.)
765 */
766
767 /*
768 gcc -g -Wall -maltivec -std=gnu99 -o optacc_utest -I.. -L.. -I../../easel -L../../easel -Dp7OPTACC_TESTDRIVE optacc.c -lhmmer -leasel -lm
769 ./optacc_utest
770 */
771 #include "p7_config.h"
772
773 #include "easel.h"
774 #include "esl_alphabet.h"
775 #include "esl_getopts.h"
776
777 #include "hmmer.h"
778 #include "impl_vmx.h"
779
780 static ESL_OPTIONS options[] = {
781 /* name type default env range toggles reqs incomp help docgroup*/
782 { "-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show brief help on version and usage", 0 },
783 { "-s", eslARG_INT, "41", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>", 0 },
784 { "-L", eslARG_INT, "50", NULL, NULL, NULL, NULL, NULL, "size of random sequences to sample", 0 },
785 { "-M", eslARG_INT, "45", NULL, NULL, NULL, NULL, NULL, "size of random models to sample", 0 },
786 { "-N", eslARG_INT, "20", NULL, NULL, NULL, NULL, NULL, "number of random sequences to sample", 0 },
787 { "--traces", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump all tracebacks", 0 },
788 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
789 };
790 static char usage[] = "[-options]";
791 static char banner[] = "test driver for VMX Forward, Backward implementations";
792
793 int
main(int argc,char ** argv)794 main(int argc, char **argv)
795 {
796 ESL_GETOPTS *go = p7_CreateDefaultApp(options, 0, argc, argv, banner, usage);
797 ESL_RANDOMNESS *r = esl_randomness_CreateFast(esl_opt_GetInteger(go, "-s"));
798 ESL_ALPHABET *abc = NULL;
799 P7_BG *bg = NULL;
800 int M = esl_opt_GetInteger(go, "-M");
801 int L = esl_opt_GetInteger(go, "-L");
802 int N = esl_opt_GetInteger(go, "-N");
803
804 /* first round of tests for DNA alphabets. */
805 if ((abc = esl_alphabet_Create(eslDNA)) == NULL) esl_fatal("failed to create alphabet");
806 if ((bg = p7_bg_Create(abc)) == NULL) esl_fatal("failed to create null model");
807
808 utest_optacc(go, r, abc, bg, M, L, N); /* normal sized models */
809 utest_optacc(go, r, abc, bg, 1, L, 10); /* size 1 models */
810 utest_optacc(go, r, abc, bg, M, 1, 10); /* size 1 sequences */
811
812 esl_alphabet_Destroy(abc);
813 p7_bg_Destroy(bg);
814
815 /* Second round of tests for amino alphabets. */
816 if ((abc = esl_alphabet_Create(eslAMINO)) == NULL) esl_fatal("failed to create alphabet");
817 if ((bg = p7_bg_Create(abc)) == NULL) esl_fatal("failed to create null model");
818
819 utest_optacc(go, r, abc, bg, M, L, N);
820 utest_optacc(go, r, abc, bg, 1, L, 10);
821 utest_optacc(go, r, abc, bg, M, 1, 10);
822
823 esl_alphabet_Destroy(abc);
824 p7_bg_Destroy(bg);
825
826 esl_getopts_Destroy(go);
827 esl_randomness_Destroy(r);
828 return eslOK;
829 }
830 #endif /*p7OPTACC_TESTDRIVE*/
831 /*------------------ end, test driver ---------------------------*/
832
833
834
835
836 /*****************************************************************
837 * 6. Example
838 *****************************************************************/
839 #ifdef p7OPTACC_EXAMPLE
840 /*
841 gcc -g -Wall -o optacc_example -Dp7OPTACC_EXAMPLE -I.. -I../../easel -L.. -L../../easel optacc.c -lhmmer -leasel -lm
842 ./optacc_example <hmmfile> <seqfile>
843 */
844 #include "p7_config.h"
845
846 #include "easel.h"
847 #include "esl_alphabet.h"
848 #include "esl_getopts.h"
849 #include "esl_sq.h"
850 #include "esl_sqio.h"
851
852 #include "hmmer.h"
853 #include "impl_vmx.h"
854
855 static ESL_OPTIONS options[] = {
856 /* name type default env range toggles reqs incomp help docgroup*/
857 { "-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show brief help on version and usage", 0 },
858 { "-d", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump posterior residue decoding matrix", 0 },
859 { "-m", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump OA matrix", 0 },
860 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
861 };
862 static char usage[] = "[-options] <hmmfile> <seqfile>";
863 static char banner[] = "example of optimal accuracy alignment, VMX implementation";
864
865 int
main(int argc,char ** argv)866 main(int argc, char **argv)
867 {
868 ESL_GETOPTS *go = p7_CreateDefaultApp(options, 2, argc, argv, banner, usage);
869 char *hmmfile = esl_opt_GetArg(go, 1);
870 char *seqfile = esl_opt_GetArg(go, 2);
871 ESL_ALPHABET *abc = NULL;
872 P7_HMMFILE *hfp = NULL;
873 P7_HMM *hmm = NULL;
874 P7_BG *bg = NULL;
875 P7_PROFILE *gm = NULL;
876 P7_OPROFILE *om = NULL;
877 P7_OMX *ox1 = NULL;
878 P7_OMX *ox2 = NULL;
879 P7_GMX *gx = NULL;
880 ESL_SQ *sq = NULL;
881 ESL_SQFILE *sqfp = NULL;
882 P7_TRACE *tr = NULL;
883 int format = eslSQFILE_UNKNOWN;
884 char errbuf[eslERRBUFSIZE];
885 float fsc, bsc;
886 float accscore;
887 int status;
888
889 /* Read in one HMM */
890 if (p7_hmmfile_OpenE(hmmfile, NULL, &hfp, NULL) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
891 if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
892 p7_hmmfile_Close(hfp);
893
894 /* Read in one sequence */
895 sq = esl_sq_CreateDigital(abc);
896 status = esl_sqfile_OpenDigital(abc, seqfile, format, NULL, &sqfp);
897 if (status == eslENOTFOUND) p7_Fail("No such file.");
898 else if (status == eslEFORMAT) p7_Fail("Format unrecognized.");
899 else if (status == eslEINVAL) p7_Fail("Can't autodetect stdin or .gz.");
900 else if (status != eslOK) p7_Fail("Open failed, code %d.", status);
901 if (esl_sqio_Read(sqfp, sq) != eslOK) p7_Fail("Failed to read sequence");
902 esl_sqfile_Close(sqfp);
903
904 /* Configure a profile from the HMM */
905 bg = p7_bg_Create(abc); p7_bg_SetLength(bg, sq->n);
906 gm = p7_profile_Create(hmm->M, abc); p7_ProfileConfig(hmm, bg, gm, sq->n, p7_LOCAL); /* multihit local: H3 default */
907 om = p7_oprofile_Create(gm->M, abc); p7_oprofile_Convert(gm, om);
908
909 /* Allocations */
910 ox1 = p7_omx_Create(gm->M, sq->n, sq->n);
911 ox2 = p7_omx_Create(gm->M, sq->n, sq->n);
912 gx = p7_gmx_Create(gm->M, sq->n);
913 tr = p7_trace_CreateWithPP();
914 p7_FLogsumInit();
915
916 /* Run Forward, Backward; do OA fill and trace */
917 p7_Forward (sq->dsq, sq->n, om, ox1, &fsc);
918 p7_Backward(sq->dsq, sq->n, om, ox1, ox2, &bsc);
919 p7_Decoding(om, ox1, ox2, ox2); /* <gx2> is now the posterior decoding matrix */
920 p7_OptimalAccuracy(om, ox2, ox1, &accscore); /* <gx1> is now the OA matrix */
921 p7_OATrace(om, ox2, ox1, tr);
922
923 if (esl_opt_GetBoolean(go, "-d")) { p7_omx_FDeconvert(ox2, gx); p7_gmx_Dump(stdout, gx, p7_DEFAULT); }
924 if (esl_opt_GetBoolean(go, "-m")) { p7_omx_FDeconvert(ox1, gx); p7_gmx_Dump(stdout, gx, p7_DEFAULT); }
925
926 p7_trace_Dump(stdout, tr, gm, sq->dsq);
927
928 if (p7_trace_Validate(tr, abc, sq->dsq, errbuf) != eslOK) p7_Die("trace fails validation:\n%s\n", errbuf);
929
930 printf("fwd = %.4f nats\n", fsc);
931 printf("bck = %.4f nats\n", bsc);
932 printf("acc = %.4f (%.2f%%)\n", accscore, accscore * 100. / (float) sq->n);
933
934 /* Cleanup */
935 esl_sq_Destroy(sq);
936 p7_trace_Destroy(tr);
937 p7_omx_Destroy(ox1);
938 p7_omx_Destroy(ox2);
939 p7_gmx_Destroy(gx);
940 p7_oprofile_Destroy(om);
941 p7_profile_Destroy(gm);
942 p7_bg_Destroy(bg);
943 p7_hmm_Destroy(hmm);
944 esl_alphabet_Destroy(abc);
945 esl_getopts_Destroy(go);
946 return 0;
947 }
948 #endif /*p7OPTACC_EXAMPLE*/
949 /*-------------------- end, example -----------------------------*/
950
951
952