1 /* Routines for the P7_PROFILE structure - Plan 7's search profile
2 *
3 * 1. The P7_PROFILE object: allocation, initialization, destruction.
4 * 2. Access methods.
5 * 3. Debugging and development code.
6 * 4. Unit tests.
7 * 5. Test driver.
8 *
9 * See also:
10 * modelconfig.c : routines that configure a profile given an HMM
11 */
12
13 #include "p7_config.h"
14
15 #include <string.h>
16 #ifdef HMMER_MPI
17 #include <mpi.h>
18 #endif
19
20 #include "easel.h"
21 #include "esl_vectorops.h"
22
23 #include "hmmer.h"
24
25
26 /*****************************************************************
27 * 1. The P7_PROFILE object: allocation, initialization, destruction.
28 *****************************************************************/
29
30 /* Function: p7_profile_Create()
31 * Synopsis: Allocates a profile.
32 *
33 * Purpose: Allocates for a profile of up to <M> nodes, for digital
34 * alphabet <abc>.
35 *
36 * Because this function might be in the critical path (in
37 * hmmscan, for example), we leave much of the model
38 * unintialized, including scores and length model
39 * probabilities. The <p7_ProfileConfig()> call is what
40 * sets these.
41 *
42 * The alignment mode is set to <p7_NO_MODE>. The
43 * reference pointer <gm->abc> is set to <abc>.
44 *
45 * Returns: a pointer to the new profile.
46 *
47 * Throws: <NULL> on allocation error.
48 *
49 * Xref: STL11/125.
50 */
51 P7_PROFILE *
p7_profile_Create(int allocM,const ESL_ALPHABET * abc)52 p7_profile_Create(int allocM, const ESL_ALPHABET *abc)
53 {
54 P7_PROFILE *gm = NULL;
55 int x;
56 int status;
57
58 /* level 0 */
59 ESL_ALLOC(gm, sizeof(P7_PROFILE));
60 gm->tsc = NULL;
61 gm->rsc = NULL;
62 gm->rf = NULL;
63 gm->mm = NULL;
64 gm->cs = NULL;
65 gm->consensus = NULL;
66
67 /* level 1 */
68 ESL_ALLOC(gm->tsc, sizeof(float) * allocM * p7P_NTRANS);
69 ESL_ALLOC(gm->rsc, sizeof(float *) * abc->Kp);
70 ESL_ALLOC(gm->rf, sizeof(char) * (allocM+2)); /* yes, +2: each is (0)1..M, +trailing \0 */
71 ESL_ALLOC(gm->mm, sizeof(char) * (allocM+2));
72 ESL_ALLOC(gm->cs, sizeof(char) * (allocM+2));
73 ESL_ALLOC(gm->consensus, sizeof(char) * (allocM+2));
74 gm->rsc[0] = NULL;
75
76 /* level 2 */
77 ESL_ALLOC(gm->rsc[0], sizeof(float) * abc->Kp * (allocM+1) * p7P_NR);
78 for (x = 1; x < abc->Kp; x++)
79 gm->rsc[x] = gm->rsc[0] + x * (allocM+1) * p7P_NR;
80
81 /* Initialize some edge pieces of memory that are never used,
82 * and are only present for indexing convenience.
83 */
84 esl_vec_FSet(gm->tsc, p7P_NTRANS, -eslINFINITY); /* node 0 nonexistent, has no transitions */
85 if (allocM > 1) {
86 p7P_TSC(gm, 1, p7P_DM) = -eslINFINITY; /* delete state D_1 is wing-retracted */
87 p7P_TSC(gm, 1, p7P_DD) = -eslINFINITY;
88 }
89 for (x = 0; x < abc->Kp; x++) {
90 p7P_MSC(gm, 0, x) = -eslINFINITY; /* no emissions from nonexistent M_0... */
91 p7P_ISC(gm, 0, x) = -eslINFINITY; /* or I_0... */
92 /* I_M is initialized in profile config, when we know actual M, not just allocated max M */
93 }
94 x = esl_abc_XGetGap(abc); /* no emission can emit/score gap characters */
95 esl_vec_FSet(gm->rsc[x], (allocM+1)*p7P_NR, -eslINFINITY);
96 x = esl_abc_XGetMissing(abc); /* no emission can emit/score missing data characters */
97 esl_vec_FSet(gm->rsc[x], (allocM+1)*p7P_NR, -eslINFINITY);
98
99 /* Set remaining info */
100 gm->mode = p7_NO_MODE;
101 gm->L = 0;
102 gm->allocM = allocM;
103 gm->M = 0;
104 gm->max_length = -1;
105 gm->nj = 0.0f;
106
107 gm->roff = -1;
108 gm->eoff = -1;
109 gm->offs[p7_MOFFSET] = -1;
110 gm->offs[p7_FOFFSET] = -1;
111 gm->offs[p7_POFFSET] = -1;
112
113 gm->name = NULL;
114 gm->acc = NULL;
115 gm->desc = NULL;
116 gm->rf[0] = 0; /* RF line is optional annotation; this flags that it's not set yet */
117 gm->mm[0] = 0; /* likewise for MM annotation line */
118 gm->cs[0] = 0; /* likewise for CS annotation line */
119 gm->consensus[0] = 0;
120
121 for (x = 0; x < p7_NEVPARAM; x++) gm->evparam[x] = p7_EVPARAM_UNSET;
122 for (x = 0; x < p7_NCUTOFFS; x++) gm->cutoff[x] = p7_CUTOFF_UNSET;
123 for (x = 0; x < p7_MAXABET; x++) gm->compo[x] = p7_COMPO_UNSET;
124
125 gm->abc = abc;
126 return gm;
127
128 ERROR:
129 p7_profile_Destroy(gm);
130 return NULL;
131 }
132
133
134 /* Function: p7_profile_Copy()
135 * Synopsis: Copy a profile.
136 *
137 * Purpose: Copies profile <src> to profile <dst>, where <dst>
138 * has already been allocated to be of sufficient size.
139 *
140 * Returns: <eslOK> on success.
141 *
142 * Throws: <eslEMEM> on allocation error; <eslEINVAL> if <dst> is too small
143 * to fit <src>.
144 */
145 int
p7_profile_Copy(const P7_PROFILE * src,P7_PROFILE * dst)146 p7_profile_Copy(const P7_PROFILE *src, P7_PROFILE *dst)
147 {
148 int x,z;
149 int status;
150
151 if (src->M > dst->allocM) ESL_EXCEPTION(eslEINVAL, "destination profile is too small to hold a copy of source profile");
152
153 esl_vec_FCopy(src->tsc, src->M*p7P_NTRANS, dst->tsc);
154 for (x = 0; x < src->abc->Kp; x++) esl_vec_FCopy(src->rsc[x], (src->M+1)*p7P_NR, dst->rsc[x]);
155 for (x = 0; x < p7P_NXSTATES; x++) esl_vec_FCopy(src->xsc[x], p7P_NXTRANS, dst->xsc[x]);
156
157 dst->mode = src->mode;
158 dst->L = src->L;
159 dst->allocM = src->allocM;
160 dst->M = src->M;
161 dst->max_length = src->max_length;
162 dst->nj = src->nj;
163
164 dst->roff = src->roff;
165 dst->eoff = src->eoff;
166 for (x = 0; x < p7_NOFFSETS; ++x) dst->offs[x] = src->offs[x];
167
168 if (dst->name != NULL) free(dst->name);
169 if (dst->acc != NULL) free(dst->acc);
170 if (dst->desc != NULL) free(dst->desc);
171
172 if ((status = esl_strdup(src->name, -1, &(dst->name))) != eslOK) return status;
173 if ((status = esl_strdup(src->acc, -1, &(dst->acc))) != eslOK) return status;
174 if ((status = esl_strdup(src->desc, -1, &(dst->desc))) != eslOK) return status;
175
176 strcpy(dst->rf, src->rf); /* RF is optional: if it's not set, *rf=0, and strcpy still works fine */
177 strcpy(dst->mm, src->mm); /* MM is also optional annotation */
178 strcpy(dst->cs, src->cs); /* CS is also optional annotation */
179 strcpy(dst->consensus, src->consensus); /* consensus though is always present on a valid profile */
180
181 for (z = 0; z < p7_NEVPARAM; z++) dst->evparam[z] = src->evparam[z];
182 for (z = 0; z < p7_NCUTOFFS; z++) dst->cutoff[z] = src->cutoff[z];
183 for (z = 0; z < p7_MAXABET; z++) dst->compo[z] = src->compo[z];
184 return eslOK;
185 }
186
187
188 /* Function: p7_profile_Clone()
189 * Synopsis: Duplicates a profile.
190 *
191 * Purpose: Duplicate profile <gm>; return a pointer
192 * to the newly allocated copy.
193 */
194 P7_PROFILE *
p7_profile_Clone(const P7_PROFILE * gm)195 p7_profile_Clone(const P7_PROFILE *gm)
196 {
197 P7_PROFILE *g2 = NULL;
198 int status;
199
200 if ((g2 = p7_profile_Create(gm->allocM, gm->abc)) == NULL) return NULL;
201 if ((status = p7_profile_Copy(gm, g2)) != eslOK) goto ERROR;
202 return g2;
203
204 ERROR:
205 p7_profile_Destroy(g2);
206 return NULL;
207 }
208
209
210
211 /* Function: p7_profile_SetNullEmissions()
212 * Synopsis: Set all emission scores to zero (experimental).
213 *
214 * Purpose: Set all emission scores in profile <gm> to zero.
215 * This makes the profile a null model, with all the same
216 * length distributions as the original model, but
217 * the emission probabilities of the background.
218 *
219 * Written to test the idea that score statistics will be
220 * even better behaved when using a null model with the
221 * same length distribution as the search model.
222 *
223 * Returns: <eslOK> on success.
224 */
225 int
p7_profile_SetNullEmissions(P7_PROFILE * gm)226 p7_profile_SetNullEmissions(P7_PROFILE *gm)
227 {
228 int x;
229 for (x = 0; x <= gm->abc->K; x++) esl_vec_FSet(gm->rsc[x], (gm->M+1)*p7P_NR, 0.0); /* canonicals */
230 for (x = gm->abc->K+1; x <= gm->abc->Kp-3; x++) esl_vec_FSet(gm->rsc[x], (gm->M+1)*p7P_NR, 0.0); /* noncanonicals */
231 return eslOK;
232 }
233
234
235 /* Function: p7_profile_Reuse()
236 * Synopsis: Prepare profile to be re-used for a new HMM.
237 *
238 * Purpose: Prepare profile <gm>'s memory to be re-used
239 * for a new HMM.
240 */
241 int
p7_profile_Reuse(P7_PROFILE * gm)242 p7_profile_Reuse(P7_PROFILE *gm)
243 {
244 /* name, acc, desc annotation is dynamically allocated for each HMM */
245 if (gm->name != NULL) { free(gm->name); gm->name = NULL; }
246 if (gm->acc != NULL) { free(gm->acc); gm->acc = NULL; }
247 if (gm->desc != NULL) { free(gm->desc); gm->desc = NULL; }
248
249 /* set annotations to empty strings */
250 gm->rf[0] = 0;
251 gm->mm[0] = 0;
252 gm->cs[0] = 0;
253 gm->consensus[0] = 0;
254
255 /* reset some other things, but leave the rest alone. */
256 gm->mode = p7_NO_MODE;
257 gm->L = 0;
258 gm->M = 0;
259 gm->nj = 0.0f;
260
261 gm->roff = -1;
262 gm->eoff = -1;
263 gm->offs[p7_MOFFSET] = -1;
264 gm->offs[p7_FOFFSET] = -1;
265 gm->offs[p7_POFFSET] = -1;
266
267 return eslOK;
268 }
269
270
271 /* Function: p7_profile_Sizeof()
272 * Synopsis: Return the allocated size of a P7_PROFILE.
273 *
274 * Purpose: Return the allocated size of a <P7_PROFILE>, in bytes.
275 */
276 size_t
p7_profile_Sizeof(P7_PROFILE * gm)277 p7_profile_Sizeof(P7_PROFILE *gm)
278 {
279 size_t n = 0;
280
281 /* these mirror malloc()'s in p7_profile_Create(); maintain one:one correspondence for maintainability */
282 n += sizeof(P7_PROFILE);
283 n += sizeof(float) * gm->allocM * p7P_NTRANS; /* gm->tsc */
284 n += sizeof(float *) * gm->abc->Kp; /* gm->rsc */
285 n += sizeof(char) * (gm->allocM+2); /* gm->rf */
286 n += sizeof(char) * (gm->allocM+2); /* gm->mm */
287 n += sizeof(char) * (gm->allocM+2); /* gm->cs */
288 n += sizeof(char) * (gm->allocM+2); /* gm->consensus */
289
290 n += sizeof(float) * gm->abc->Kp * (gm->allocM+1) * p7P_NR; /* gm->rsc[0] */
291
292 return n;
293 }
294
295
296 /* Function: p7_profile_Destroy()
297 * Synopsis: Frees a profile.
298 *
299 * Purpose: Frees a profile <gm>.
300 *
301 * Returns: (void).
302 *
303 * Xref: STL11/125.
304 */
305 void
p7_profile_Destroy(P7_PROFILE * gm)306 p7_profile_Destroy(P7_PROFILE *gm)
307 {
308 if (gm != NULL) {
309 if (gm->rsc != NULL && gm->rsc[0] != NULL) free(gm->rsc[0]);
310 if (gm->tsc != NULL) free(gm->tsc);
311 if (gm->rsc != NULL) free(gm->rsc);
312 if (gm->name != NULL) free(gm->name);
313 if (gm->acc != NULL) free(gm->acc);
314 if (gm->desc != NULL) free(gm->desc);
315 if (gm->rf != NULL) free(gm->rf);
316 if (gm->mm != NULL) free(gm->mm);
317 if (gm->cs != NULL) free(gm->cs);
318 if (gm->consensus != NULL) free(gm->consensus);
319 free(gm);
320 }
321 return;
322 }
323
324
325 /*****************************************************************
326 * 2. Access methods.
327 *****************************************************************/
328
329 /* Function: p7_profile_IsLocal()
330 * Synopsis: Return TRUE if profile is in a local alignment mode.
331 *
332 * Purpose: Return <TRUE> if profile is in a local alignment mode.
333 */
334 int
p7_profile_IsLocal(const P7_PROFILE * gm)335 p7_profile_IsLocal(const P7_PROFILE *gm)
336 {
337 if (gm->mode == p7_UNILOCAL || gm->mode == p7_LOCAL) return TRUE;
338 return FALSE;
339 }
340
341 /* Function: p7_profile_IsMultihit()
342 * Synopsis: Return TRUE if profile is in a multihit alignment mode.
343 *
344 * Purpose: Return <TRUE> if profile is in a multihit alignment mode.
345 */
346 int
p7_profile_IsMultihit(const P7_PROFILE * gm)347 p7_profile_IsMultihit(const P7_PROFILE *gm)
348 {
349 if (gm->mode == p7_LOCAL || gm->mode == p7_GLOCAL) return TRUE;
350 return FALSE;
351 }
352
353
354
355
356 /* Function: p7_profile_GetT()
357 *
358 * Purpose: Convenience function that looks up a transition score in
359 * profile <gm> for a transition from state type <st1> in
360 * node <k1> to state type <st2> in node <k2>. For unique
361 * state types that aren't in nodes (<p7T_S>, for example), the
362 * <k> value is ignored, though it would be customarily passed as 0.
363 * Return the transition score in <ret_tsc>.
364 *
365 * This function would almost always be called on profile
366 * traces, of course, but it's possible to call it
367 * on core traces (for example, if you were to try to
368 * trace_Dump() during HMM construction, and you wanted
369 * to see detailed profile scores for that trace). Core traces
370 * can contain <p7T_X> "states" used solely to signal
371 * a sequence fragment, treated as missing data. Transitions
372 * involving <p7T_X> states are assigned zero score here.
373 * Other transitions that occur only in core traces
374 * (B->I0, B->D1, I_M->E) also silently get a zero score.
375 * This is safe, because we would only ever use this number
376 * for display, not as a log probability somewhere.
377 *
378 * Returns: <eslOK> on success, and <*ret_tsc> contains the requested
379 * transition score.
380 *
381 * Throws: <eslEINVAL> if a nonexistent transition is requested. Now
382 * <*ret_tsc> is set to $-\infty$.
383 *
384 */
385 int
p7_profile_GetT(const P7_PROFILE * gm,char st1,int k1,char st2,int k2,float * ret_tsc)386 p7_profile_GetT(const P7_PROFILE *gm, char st1, int k1, char st2, int k2, float *ret_tsc)
387 {
388 float tsc = 0.0f;
389 int status;
390
391 /* Detect transitions that can only come from core traces;
392 * return 0.0 as a special case (this is only done for displaying
393 * "scores" in trace dumps, during debugging.)
394 */
395 if (st1 == p7T_X || st2 == p7T_X) return eslOK;
396 if (st1 == p7T_B && st2 == p7T_I) return eslOK;
397 if (st1 == p7T_B && st2 == p7T_D) return eslOK;
398 if (st1 == p7T_I && st2 == p7T_E) return eslOK;
399
400 /* Now we're sure this is a profile trace, as it should usually be. */
401 switch (st1) {
402 case p7T_S: break;
403 case p7T_T: break;
404 case p7T_N:
405 switch (st2) {
406 case p7T_B: tsc = gm->xsc[p7P_N][p7P_MOVE]; break;
407 case p7T_N: tsc = gm->xsc[p7P_N][p7P_LOOP]; break;
408 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s->%s", p7_hmm_DecodeStatetype(st1), p7_hmm_DecodeStatetype(st2));
409 }
410 break;
411
412 case p7T_B:
413 switch (st2) {
414 case p7T_M: tsc = p7P_TSC(gm, k2-1, p7P_BM); break; /* remember, B->Mk is stored in [k-1][p7P_BM] */
415 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s->%s", p7_hmm_DecodeStatetype(st1), p7_hmm_DecodeStatetype(st2));
416 }
417 break;
418
419 case p7T_M:
420 switch (st2) {
421 case p7T_M: tsc = p7P_TSC(gm, k1, p7P_MM); break;
422 case p7T_I: tsc = p7P_TSC(gm, k1, p7P_MI); break;
423 case p7T_D: tsc = p7P_TSC(gm, k1, p7P_MD); break;
424 case p7T_E:
425 if (k1 != gm->M && ! p7_profile_IsLocal(gm)) ESL_EXCEPTION(eslEINVAL, "local end transition (M%d of %d) in non-local model", k1, gm->M);
426 tsc = 0.0f; /* by def'n in H3 local alignment */
427 break;
428 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s_%d->%s", p7_hmm_DecodeStatetype(st1), k1, p7_hmm_DecodeStatetype(st2));
429 }
430 break;
431
432 case p7T_D:
433 switch (st2) {
434 case p7T_M: tsc = p7P_TSC(gm, k1, p7P_DM); break;
435 case p7T_D: tsc = p7P_TSC(gm, k1, p7P_DD); break;
436 case p7T_E:
437 if (k1 != gm->M && ! p7_profile_IsLocal(gm)) ESL_EXCEPTION(eslEINVAL, "local end transition (D%d of %d) in non-local model", k1, gm->M);
438 tsc = 0.0f; /* by def'n in H3 local alignment */
439 break;
440 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s_%d->%s", p7_hmm_DecodeStatetype(st1), k1, p7_hmm_DecodeStatetype(st2));
441 }
442 break;
443
444 case p7T_I:
445 switch (st2) {
446 case p7T_M: tsc = p7P_TSC(gm, k1, p7P_IM); break;
447 case p7T_I: tsc = p7P_TSC(gm, k1, p7P_II); break;
448 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s_%d->%s", p7_hmm_DecodeStatetype(st1), k1, p7_hmm_DecodeStatetype(st2));
449 }
450 break;
451
452 case p7T_E:
453 switch (st2) {
454 case p7T_C: tsc = gm->xsc[p7P_E][p7P_MOVE]; break;
455 case p7T_J: tsc = gm->xsc[p7P_E][p7P_LOOP]; break;
456 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s->%s", p7_hmm_DecodeStatetype(st1), p7_hmm_DecodeStatetype(st2));
457 }
458 break;
459
460 case p7T_J:
461 switch (st2) {
462 case p7T_B: tsc = gm->xsc[p7P_J][p7P_MOVE]; break;
463 case p7T_J: tsc = gm->xsc[p7P_J][p7P_LOOP]; break;
464 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s->%s", p7_hmm_DecodeStatetype(st1), p7_hmm_DecodeStatetype(st2));
465 }
466 break;
467
468 case p7T_C:
469 switch (st2) {
470 case p7T_T: tsc = gm->xsc[p7P_C][p7P_MOVE]; break;
471 case p7T_C: tsc = gm->xsc[p7P_C][p7P_LOOP]; break;
472 default: ESL_XEXCEPTION(eslEINVAL, "bad transition %s->%s", p7_hmm_DecodeStatetype(st1), p7_hmm_DecodeStatetype(st2));
473 }
474 break;
475
476 default: ESL_XEXCEPTION(eslEINVAL, "bad state type %d in traceback", st1);
477 }
478
479 *ret_tsc = tsc;
480 return eslOK;
481
482 ERROR:
483 *ret_tsc = -eslINFINITY;
484 return status;
485 }
486
487
488 /*****************************************************************
489 * 3. Debugging and development code.
490 *****************************************************************/
491
492 /* Function: p7_profile_Validate()
493 *
494 * Purpose: Validates the internals of the generic profile structure
495 * <gm>.
496 *
497 * TODO: currently this function is incomplete, and only
498 * validates the entry distribution.
499 *
500 * Returns: <eslOK> if <gm> internals look fine. Returns <eslFAIL>
501 * if something is wrong, and leaves an error message in
502 * <errbuf> if caller passed it non-<NULL>.
503 */
504 int
p7_profile_Validate(const P7_PROFILE * gm,char * errbuf,float tol)505 p7_profile_Validate(const P7_PROFILE *gm, char *errbuf, float tol)
506 {
507 int status;
508 int k;
509 double *pstart = NULL;
510
511 ESL_ALLOC(pstart, sizeof(double) * (gm->M+1));
512 pstart[0] = 0.0;
513
514 /* Validate the entry distribution.
515 * In a glocal model, this is an explicit probability distribution,
516 * corresponding to left wing retraction.
517 * In a local model, this is an implicit probability distribution,
518 * corresponding to the implicit local alignment model, and we have
519 * to calculate the M(M+1)/2 fragment probabilities accordingly.
520 */
521 if (p7_profile_IsLocal(gm))
522 { /* the code block below is also in emit.c:sample_endpoints */
523 for (k = 1; k <= gm->M; k++)
524 pstart[k] = exp(p7P_TSC(gm, k-1, p7P_BM)) * (gm->M - k + 1); /* multiply p_ij by the number of exits j */
525 }
526 else
527 {
528 for (k = 1; k <= gm->M; k++)
529 pstart[k] = exp(p7P_TSC(gm, k-1, p7P_BM));
530 }
531
532 if (esl_vec_DValidate(pstart, gm->M+1, tol, NULL) != eslOK) ESL_XFAIL(eslFAIL, errbuf, "profile entry distribution is not normalized properly");
533 free(pstart);
534 return eslOK;
535
536 ERROR:
537 if (pstart != NULL) free(pstart);
538 return eslFAIL;
539 }
540
541 /* Function: p7_profile_Compare()
542 * Synopsis: Compare two profiles for equality.
543 *
544 * Purpose: Compare two profiles <gm1> and <gm2> to each other.
545 * Return <eslOK> if they're identical, and <eslFAIL> if
546 * they differ. Floating-point probabilities are
547 * compared for equality within a fractional tolerance
548 * <tol>. Only compares the scores, not any annotation
549 * on the profiles.
550 */
551 int
p7_profile_Compare(P7_PROFILE * gm1,P7_PROFILE * gm2,float tol)552 p7_profile_Compare(P7_PROFILE *gm1, P7_PROFILE *gm2, float tol)
553 {
554 int x;
555
556 if (gm1->mode != gm2->mode) return eslFAIL;
557 if (gm1->M != gm2->M) return eslFAIL;
558
559 if (esl_vec_FCompare(gm1->tsc, gm2->tsc, gm1->M*p7P_NTRANS, tol) != eslOK) return eslFAIL;
560 for (x = 0; x < gm1->abc->Kp; x++)
561 if (esl_vec_FCompare(gm1->rsc[x], gm2->rsc[x], (gm1->M+1)*p7P_NR, tol) != eslOK) return eslFAIL;
562
563 for (x = 0; x < p7P_NXSTATES; x++)
564 if (esl_vec_FCompare(gm1->xsc[x], gm2->xsc[x], p7P_NXTRANS, tol) != eslOK) return eslFAIL;
565
566 return eslOK;
567 }
568
569
570
571 /*****************************************************************
572 * 4. Unit tests
573 *****************************************************************/
574 #ifdef p7PROFILE_TESTDRIVE
575 #include "esl_alphabet.h"
576 #include "esl_random.h"
577
578 static void
utest_Compare(void)579 utest_Compare(void)
580 {
581 ESL_RANDOMNESS *r = esl_randomness_CreateFast(42);
582 ESL_ALPHABET *abc = esl_alphabet_Create(eslAMINO);
583 P7_HMM *hmm = NULL;
584 P7_BG *bg = NULL;
585 P7_PROFILE *gm = NULL;
586 P7_PROFILE *gm2 = NULL;
587 int M = 200;
588 int L = 400;
589
590 p7_hmm_Sample(r, M, abc, &hmm); /* master and worker's sampled profiles are identical */
591 bg = p7_bg_Create(abc);
592 gm = p7_profile_Create(hmm->M, abc);
593 gm2 = p7_profile_Create(hmm->M, abc);
594 p7_ProfileConfig(hmm, bg, gm, 400, p7_LOCAL);
595 p7_ProfileConfig(hmm, bg, gm2, 400, p7_LOCAL);
596 p7_ReconfigLength(gm, L);
597 p7_ReconfigLength(gm2, L);
598
599 if (p7_profile_Compare(gm, gm2, 0.001) != eslOK) p7_Die("identical profile comparison failed");
600
601 p7_profile_Destroy(gm);
602 p7_profile_Destroy(gm2);
603 p7_bg_Destroy(bg);
604 p7_hmm_Destroy(hmm);
605 esl_alphabet_Destroy(abc);
606 esl_randomness_Destroy(r);
607 return;
608 }
609
610
611 #endif /*p7PROFILE_TESTDRIVE*/
612
613 /*****************************************************************
614 * 5. Test driver
615 *****************************************************************/
616 #ifdef p7PROFILE_TESTDRIVE
617
618 /* gcc -o profile_utest -g -Wall -I../easel -L../easel -I. -L. -Dp7PROFILE_TESTDRIVE p7_profile.c -lhmmer -leasel -lm
619 * ./profile_utest
620 */
621 #include "esl_getopts.h"
622
623 static ESL_OPTIONS options[] = {
624 /* name type default env range toggles reqs incomp help docgroup*/
625 { "-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show brief help on version and usage", 0 },
626 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
627 };
628 static char usage[] = "[-options]";
629 static char banner[] = "test driver for p7_profile.c";
630
631 int
main(int argc,char ** argv)632 main(int argc, char **argv)
633 {
634 ESL_GETOPTS *go = p7_CreateDefaultApp(options, 0, argc, argv, banner, usage);
635
636 utest_Compare();
637
638 esl_getopts_Destroy(go);
639 return 0;
640 }
641 #endif /*p7PROFILE_TESTDRIVE*/
642
643