本文整理汇总了C++中p7_bg_SetLength函数的典型用法代码示例。如果您正苦于以下问题：C++ p7_bg_SetLength函数的具体用法？C++ p7_bg_SetLength怎么用？C++ p7_bg_SetLength使用的例子？那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了p7_bg_SetLength函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: main
int
main(int argc, char **argv)
{
char *msg = "p7_gmx unit test driver failed";
ESL_GETOPTS *go = p7_CreateDefaultApp(options, 0, argc, argv, banner, usage);
ESL_RANDOMNESS *r = esl_randomness_CreateFast(esl_opt_GetInteger(go, "-s"));
ESL_ALPHABET *abc = esl_alphabet_Create(eslAMINO);
P7_BG *bg = p7_bg_Create(abc);
P7_HMM *hmm = NULL;
P7_PROFILE *gm = NULL;
int M = esl_opt_GetInteger(go, "-M");
int L = esl_opt_GetInteger(go, "-L");
float tol = esl_opt_GetReal (go, "-t");
p7_FLogsumInit();
if (p7_hmm_Sample(r, M, abc, &hmm) != eslOK) esl_fatal(msg);
if ((gm = p7_profile_Create(hmm->M, abc)) == NULL) esl_fatal(msg);
if (p7_bg_SetLength(bg, L) != eslOK) esl_fatal(msg);
if (p7_ProfileConfig(hmm, bg, gm, L, p7_UNILOCAL) != eslOK) esl_fatal(msg);
utest_GrowTo();
utest_Compare(r, gm, bg, L, tol);
esl_getopts_Destroy(go);
esl_randomness_Destroy(r);
esl_alphabet_Destroy(abc);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
p7_profile_Destroy(gm);
return eslOK;
}示例2: main
int
main(int argc, char **argv)
{
ESL_GETOPTS *go = p7_CreateDefaultApp(options, 1, argc, argv, banner, usage);
char *hmmfile = esl_opt_GetArg(go, 1);
ESL_STOPWATCH *w = esl_stopwatch_Create();
ESL_RANDOMNESS *r = esl_randomness_CreateFast(esl_opt_GetInteger(go, "-s"));
ESL_ALPHABET *abc = NULL;
P7_HMMFILE *hfp = NULL;
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_GMX *gx1 = NULL;
P7_GMX *gx2 = NULL;
int L = esl_opt_GetInteger(go, "-L");
int N = esl_opt_GetInteger(go, "-N");
ESL_DSQ *dsq = malloc(sizeof(ESL_DSQ) * (L+2));
float null2[p7_MAXCODE];
int i;
float fsc, bsc;
double Mcs;
if (p7_hmmfile_OpenE(hmmfile, NULL, &hfp, NULL) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
bg = p7_bg_Create(abc);
p7_bg_SetLength(bg, L);
gm = p7_profile_Create(hmm->M, abc);
p7_ProfileConfig(hmm, bg, gm, L, p7_LOCAL);
gx1 = p7_gmx_Create(gm->M, L);
gx2 = p7_gmx_Create(gm->M, L);
esl_rsq_xfIID(r, bg->f, abc->K, L, dsq);
p7_GForward (dsq, L, gm, gx1, &fsc);
p7_GBackward(dsq, L, gm, gx2, &bsc);
p7_GDecoding(gm, gx1, gx2, gx2);
esl_stopwatch_Start(w);
for (i = 0; i < N; i++)
p7_GNull2_ByExpectation(gm, gx2, null2);
esl_stopwatch_Stop(w);
Mcs = (double) N * (double) L * (double) gm->M * 1e-6 / w->user;
esl_stopwatch_Display(stdout, w, "# CPU time: ");
printf("# M = %d\n", gm->M);
printf("# %.1f Mc/s\n", Mcs);
free(dsq);
p7_gmx_Destroy(gx1);
p7_gmx_Destroy(gx2);
p7_profile_Destroy(gm);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
p7_hmmfile_Close(hfp);
esl_alphabet_Destroy(abc);
esl_stopwatch_Destroy(w);
esl_randomness_Destroy(r);
esl_getopts_Destroy(go);
return 0;
}示例3: main
int
main(int argc, char **argv)
{
ESL_GETOPTS *go = p7_CreateDefaultApp(options, 1, argc, argv, banner, usage);
ESL_RANDOMNESS *rng = esl_randomness_CreateFast(esl_opt_GetInteger(go, "-s"));
char *hmmfile = esl_opt_GetArg(go, 1);
int L = esl_opt_GetInteger(go, "-L");
int N = esl_opt_GetInteger(go, "-N");
ESL_ALPHABET *abc = NULL;
P7_HMMFILE *hfp = NULL;
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_TRACE *tr = p7_trace_Create();
ESL_SQ *sq = NULL;
char errbuf[eslERRBUFSIZE];
int i;
int status;
status = p7_hmmfile_OpenE(hmmfile, NULL, &hfp, errbuf);
if (status == eslENOTFOUND) p7_Fail("File existence/permissions problem in trying to open HMM file %s.\n%s\n", hmmfile, errbuf);
else if (status == eslEFORMAT) p7_Fail("File format problem in trying to open HMM file %s.\n%s\n", hmmfile, errbuf);
else if (status != eslOK) p7_Fail("Unexpected error %d in opening HMM file %s.\n%s\n", status, hmmfile, errbuf);
status = p7_hmmfile_Read(hfp, &abc, &hmm);
if (status == eslEFORMAT) p7_Fail("Bad file format in HMM file %s:\n%s\n", hfp->fname, hfp->errbuf);
else if (status == eslEINCOMPAT) p7_Fail("HMM in %s is not in the expected %s alphabet\n", hfp->fname, esl_abc_DecodeType(abc->type));
else if (status == eslEOF) p7_Fail("Empty HMM file %s? No HMM data found.\n", hfp->fname);
else if (status != eslOK) p7_Fail("Unexpected error in reading HMMs from %s\n", hfp->fname);
p7_hmmfile_Close(hfp);
bg = p7_bg_Create(abc); p7_bg_SetLength(bg, L);
gm = p7_profile_Create(hmm->M, abc); p7_ProfileConfig(hmm, bg, gm, L, p7_LOCAL);
sq = esl_sq_CreateDigital(abc);
for (i = 0; i < N; i++)
{
p7_ProfileEmit(rng, hmm, gm, bg, sq, tr);
esl_sq_FormatName(sq, "%s-sample%d", hmm->name, i);
esl_sqio_Write(stdout, sq, eslSQFILE_FASTA, FALSE);
if (p7_trace_Validate(tr, abc, sq->dsq, errbuf) != eslOK) esl_fatal(errbuf);
esl_sq_Reuse(sq);
p7_trace_Reuse(tr);
}
esl_sq_Destroy(sq);
p7_trace_Destroy(tr);
p7_profile_Destroy(gm);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
esl_alphabet_Destroy(abc);
esl_randomness_Destroy(rng);
esl_getopts_Destroy(go);
return 0;
}示例4: utest_oprofileSendRecv
static void
utest_oprofileSendRecv(int my_rank, int nproc)
{
ESL_RANDOMNESS *r = esl_randomness_CreateFast(42);
ESL_ALPHABET *abc = esl_alphabet_Create(eslAMINO);
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_OPROFILE *om = NULL;
P7_OPROFILE *om2 = NULL;
int M = 200;
int L = 400;
char *wbuf = NULL;
int wn = 0;
int i;
char errbuf[eslERRBUFSIZE];
p7_hmm_Sample(r, M, abc, &hmm); /* master and worker's sampled profiles are identical */
bg = p7_bg_Create(abc);
gm = p7_profile_Create(hmm->M, abc);
om = p7_oprofile_Create(hmm->M, abc);
p7_ProfileConfig(hmm, bg, gm, L, p7_LOCAL);
p7_oprofile_Convert(gm, om);
p7_bg_SetLength (bg, L);
if (my_rank == 0)
{
for (i = 1; i < nproc; i++)
{
ESL_DPRINTF1(("Master: receiving test profile\n"));
p7_oprofile_MPIRecv(MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &wbuf, &wn, &abc, &om2);
ESL_DPRINTF1(("Master: test profile received\n"));
if (p7_oprofile_Compare(om, om2, 0.001, errbuf) != eslOK)
p7_Die("Received profile not identical to what was sent\n%s", errbuf);
p7_oprofile_Destroy(om2);
}
}
else
{
ESL_DPRINTF1(("Worker %d: sending test profile\n", my_rank));
p7_oprofile_MPISend(om, 0, 0, MPI_COMM_WORLD, &wbuf, &wn);
ESL_DPRINTF1(("Worker %d: test profile sent\n", my_rank));
}
free(wbuf);
p7_profile_Destroy(gm);
p7_oprofile_Destroy(om);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
esl_alphabet_Destroy(abc);
esl_randomness_Destroy(r);
return;
}示例5: emit_sequences
static void
emit_sequences(ESL_GETOPTS *go, FILE *ofp, int outfmt, ESL_RANDOMNESS *r, P7_HMM *hmm)
{
ESL_SQ *sq = NULL;
P7_TRACE *tr = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
int do_profile = esl_opt_GetBoolean(go, "-p");
int N = esl_opt_GetInteger(go, "-N");
int L = esl_opt_GetInteger(go, "-L");
int mode = p7_LOCAL;
int nseq;
int status;
if (esl_opt_GetBoolean(go, "--local")) mode = p7_LOCAL;
else if (esl_opt_GetBoolean(go, "--unilocal")) mode = p7_UNILOCAL;
else if (esl_opt_GetBoolean(go, "--glocal")) mode = p7_GLOCAL;
else if (esl_opt_GetBoolean(go, "--uniglocal")) mode = p7_UNIGLOCAL;
if ((sq = esl_sq_CreateDigital(hmm->abc)) == NULL) esl_fatal("failed to allocate sequence");
if ((tr = p7_trace_Create()) == NULL) esl_fatal("failed to allocate trace");
if ((bg = p7_bg_Create(hmm->abc)) == NULL) esl_fatal("failed to create null model");
if ((gm = p7_profile_Create(hmm->M, hmm->abc)) == NULL) esl_fatal("failed to create profile");
if (p7_ProfileConfig(hmm, bg, gm, L, mode) != eslOK) esl_fatal("failed to configure profile");
if (p7_bg_SetLength(bg, L) != eslOK) esl_fatal("failed to reconfig null model length");
if (p7_hmm_Validate (hmm, NULL, 0.0001) != eslOK) esl_fatal("whoops, HMM is bad!");
if (p7_profile_Validate(gm, NULL, 0.0001) != eslOK) esl_fatal("whoops, profile is bad!");
for (nseq = 1; nseq <= N; nseq++)
{
if (do_profile) status = p7_ProfileEmit(r, hmm, gm, bg, sq, tr);
else status = p7_CoreEmit (r, hmm, sq, tr);
if (status) esl_fatal("Failed to emit sequence\n");
status = esl_sq_FormatName(sq, "%s-sample%d", hmm->name, nseq);
if (status) esl_fatal("Failed to set sequence name\n");
status = esl_sqio_Write(ofp, sq, outfmt, FALSE);
if (status != eslOK) esl_fatal("Failed to write sequence\n");
p7_trace_Reuse(tr);
esl_sq_Reuse(sq);
}
esl_sq_Destroy(sq);
p7_trace_Destroy(tr);
p7_bg_Destroy(bg);
p7_profile_Destroy(gm);
return;
}示例6: p7_ViterbiMu
/* Function: p7_ViterbiMu()
* Synopsis: Determines the local Viterbi Gumbel mu parameter for a model.
* Incept: SRE, Tue May 19 10:26:19 2009 [Janelia]
*
* Purpose: Identical to p7_MSVMu(), above, except that it fits
* Viterbi scores instead of MSV scores.
*
* The difference between the two mus is small, but can be
* up to ~1 bit or so for large, low-info models [J4/126] so
* decided to calibrate the two mus separately [J5/8].
*
* Args: r : source of random numbers
* om : score profile (length config is changed upon return!)
* bg : null model (length config is changed upon return!)
* L : length of sequences to simulate
* N : number of sequences to simulate
* lambda : known Gumbel lambda parameter
* ret_vmu : RETURN: ML estimate of location param mu
*
* Returns: <eslOK> on success, and <ret_mu> contains the ML estimate
* of $\mu$.
*
* Throws: (no abnormal error conditions)
*/
int
p7_ViterbiMu(ESL_RANDOMNESS *r, P7_OPROFILE *om, P7_BG *bg, int L, int N, double lambda, double *ret_vmu)
{
P7_OMX *ox = p7_omx_Create(om->M, 0, 0); /* DP matrix: 1 row version */
ESL_DSQ *dsq = NULL;
double *xv = NULL;
int i;
float sc, nullsc;
#ifndef p7_IMPL_DUMMY
float maxsc = (32767.0 - om->base_w) / om->scale_w; /* if score overflows, use this [J4/139] */
#endif
int status;
if (ox == NULL) { status = eslEMEM; goto ERROR; }
ESL_ALLOC(xv, sizeof(double) * N);
ESL_ALLOC(dsq, sizeof(ESL_DSQ) * (L+2));
p7_oprofile_ReconfigLength(om, L);
p7_bg_SetLength(bg, L);
for (i = 0; i < N; i++)
{
if ((status = esl_rsq_xfIID(r, bg->f, om->abc->K, L, dsq)) != eslOK) goto ERROR;
if ((status = p7_bg_NullOne(bg, dsq, L, &nullsc)) != eslOK) goto ERROR;
status = p7_ViterbiFilter(dsq, L, om, ox, &sc);
#ifndef p7_IMPL_DUMMY
if (status == eslERANGE) { sc = maxsc; status = eslOK; }
#endif
if (status != eslOK) goto ERROR;
xv[i] = (sc - nullsc) / eslCONST_LOG2;
}
if ((status = esl_gumbel_FitCompleteLoc(xv, N, lambda, ret_vmu)) != eslOK) goto ERROR;
p7_omx_Destroy(ox);
free(xv);
free(dsq);
return eslOK;
ERROR:
*ret_vmu = 0.0;
if (ox != NULL) p7_omx_Destroy(ox);
if (xv != NULL) free(xv);
if (dsq != NULL) free(dsq);
return status;
}示例7: p7_Tau
/* Function: p7_Tau()
* Synopsis: Determine Forward tau by brief simulation.
* Incept: SRE, Thu Aug 9 15:08:39 2007 [Janelia]
*
* Purpose: Determine the <tau> parameter for an exponential tail fit
* to the Forward score distribution for model <om>, on
* random sequences with the composition of the background
* model <bg>. This <tau> parameter is for an exponential
* distribution anchored from $P=1.0$, so it's not really a
* tail per se; but it's only an accurate fit in the tail
* of the Forward score distribution, from about $P=0.001$
* or so.
*
* The determination of <tau> is done by a brief simulation
* in which we fit a Gumbel distribution to a small number
* of Forward scores of random sequences, and use that to
* predict the location of the tail at probability <tailp>.
*
* The Gumbel is of course inaccurate, but we can use it
* here solely as an empirical distribution to determine
* the location of a reasonable <tau> more accurately on a
* smaller number of samples than we could do with raw
* order statistics.
*
* Typical choices are L=100, N=200, tailp=0.04, which
* typically yield estimates $\hat{\mu}$ with a precision
* (standard deviation) of $\pm$ 0.2 bits, corresponding to
* a $\pm$ 15\% error in E-values. See [J1/135].
*
* The use of Gumbel fitting to a small number of $N$
* samples and the extrapolation of $\hat{\mu}$ from the
* estimated location of the 0.04 tail mass are both
* empirical and carefully optimized against several
* tradeoffs. Most importantly, around this choice of tail
* probability, a systematic error introduced by the use of
* the Gumbel fit is being cancelled by systematic error
* introduced by the use of a higher tail probability than
* the regime in which the exponential tail is a valid
* approximation. See [J1/135] for discussion.
*
* This function changes the length configuration of both
* <om> and <bg>. The caller must remember to reconfigure
* both of their length models appropriately for any
* subsequent alignments.
*
* Args: r : source of randomness
* om : configured profile to sample sequences from
* bg : null model (for background residue frequencies)
* L : mean length model for seq emission from profile
* N : number of sequences to generate
* lambda : expected slope of the exponential tail (from p7_Lambda())
* tailp : tail mass from which we will extrapolate mu
* ret_mu : RETURN: estimate for the Forward mu (base of exponential tail)
*
* Returns: <eslOK> on success, and <*ret_fv> is the score difference
* in bits.
*
* Throws: <eslEMEM> on allocation error, and <*ret_fv> is 0.
*/
int
p7_Tau(ESL_RANDOMNESS *r, P7_OPROFILE *om, P7_BG *bg, int L, int N, double lambda, double tailp, double *ret_tau)
{
P7_OMX *ox = p7_omx_Create(om->M, 0, L); /* DP matrix: for ForwardParser, L rows */
ESL_DSQ *dsq = NULL;
double *xv = NULL;
float fsc, nullsc;
double gmu, glam;
int status;
int i;
ESL_ALLOC(xv, sizeof(double) * N);
ESL_ALLOC(dsq, sizeof(ESL_DSQ) * (L+2));
if (ox == NULL) { status = eslEMEM; goto ERROR; }
p7_oprofile_ReconfigLength(om, L);
p7_bg_SetLength(bg, L);
for (i = 0; i < N; i++)
{
if ((status = esl_rsq_xfIID(r, bg->f, om->abc->K, L, dsq)) != eslOK) goto ERROR;
if ((status = p7_ForwardParser(dsq, L, om, ox, &fsc)) != eslOK) goto ERROR;
if ((status = p7_bg_NullOne(bg, dsq, L, &nullsc)) != eslOK) goto ERROR;
xv[i] = (fsc - nullsc) / eslCONST_LOG2;
}
if ((status = esl_gumbel_FitComplete(xv, N, &gmu, &glam)) != eslOK) goto ERROR;
/* Explanation of the eqn below: first find the x at which the Gumbel tail
* mass is predicted to be equal to tailp. Then back up from that x
* by log(tailp)/lambda to set the origin of the exponential tail to 1.0
* instead of tailp.
*/
*ret_tau = esl_gumbel_invcdf(1.0-tailp, gmu, glam) + (log(tailp) / lambda);
free(xv);
free(dsq);
p7_omx_Destroy(ox);
return eslOK;
ERROR:
*ret_tau = 0.;
//.........这里部分代码省略.........
示例8: main
int
main(int argc, char **argv)
{
ESL_GETOPTS *go = esl_getopts_CreateDefaultApp(options, 1, argc, argv, banner, usage);
char *hmmfile = esl_opt_GetArg(go, 1);
ESL_STOPWATCH *w = esl_stopwatch_Create();
ESL_ALPHABET *abc = NULL;
P7_HMMFILE *hfp = NULL;
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_OPROFILE *om = NULL;
int L = esl_opt_GetInteger(go, "-L");
int N = esl_opt_GetInteger(go, "-N");
int i;
if (p7_hmmfile_Open(hmmfile, NULL, &hfp) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
bg = p7_bg_Create(abc);
p7_bg_SetLength(bg, L);
gm = p7_profile_Create(hmm->M, abc);
p7_ProfileConfig(hmm, bg, gm, L, p7_LOCAL);
om = p7_oprofile_Create(gm->M, abc);
esl_stopwatch_Start(w);
for (i = 0; i < N; i++)
p7_oprofile_Convert(gm, om);
esl_stopwatch_Stop(w);
esl_stopwatch_Display(stdout, w, "# CPU time: ");
printf("# M = %d\n", gm->M);
p7_oprofile_Destroy(om);
p7_profile_Destroy(gm);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
p7_hmmfile_Close(hfp);
esl_alphabet_Destroy(abc);
esl_stopwatch_Destroy(w);
esl_getopts_Destroy(go);
return 0;
}示例9: serial_master
static void
serial_master(ESL_GETOPTS *go, struct cfg_s *cfg)
{
P7_HMM *hmm = NULL;
double *xv = NULL; /* results: array of N scores */
int *av = NULL; /* optional results: array of N alignment lengths */
char errbuf[eslERRBUFSIZE];
int status;
if ((status = init_master_cfg(go, cfg, errbuf)) != eslOK) p7_Fail(errbuf);
if ((xv = malloc(sizeof(double) * cfg->N)) == NULL) p7_Fail("allocation failed");
if (esl_opt_GetBoolean(go, "-a") &&
(av = malloc(sizeof(int) * cfg->N)) == NULL) p7_Fail("allocation failed");
while ((status = p7_hmmfile_Read(cfg->hfp, &(cfg->abc), &hmm)) != eslEOF)
{
if (status == eslEOD) p7_Fail("read failed, HMM file %s may be truncated?", cfg->hmmfile);
else if (status == eslEFORMAT) p7_Fail("bad file format in HMM file %s", cfg->hmmfile);
else if (status == eslEINCOMPAT) p7_Fail("HMM file %s contains different alphabets", cfg->hmmfile);
else if (status != eslOK) p7_Fail("Unexpected error in reading HMMs from %s", cfg->hmmfile);
if (cfg->bg == NULL) {
if (esl_opt_GetBoolean(go, "--bgflat")) cfg->bg = p7_bg_CreateUniform(cfg->abc);
else cfg->bg = p7_bg_Create(cfg->abc);
p7_bg_SetLength(cfg->bg, esl_opt_GetInteger(go, "-L")); /* set the null model background length in both master and workers. */
}
if (esl_opt_GetBoolean(go, "--recal")) {
if (recalibrate_model(go, cfg, errbuf, hmm) != eslOK) p7_Fail(errbuf);
}
if (process_workunit(go, cfg, errbuf, hmm, xv, av) != eslOK) p7_Fail(errbuf);
if (output_result (go, cfg, errbuf, hmm, xv, av) != eslOK) p7_Fail(errbuf);
p7_hmm_Destroy(hmm);
}
free(xv);
if (av != NULL) free(av);
}示例10: main
int
main(int argc, char **argv)
{
ESL_GETOPTS *go = p7_CreateDefaultApp(options, 1, argc, argv, banner, usage);
char *hmmfile = esl_opt_GetArg(go, 1);
ESL_STOPWATCH *w = esl_stopwatch_Create();
ESL_RANDOMNESS *r = esl_randomness_CreateFast(esl_opt_GetInteger(go, "-s"));
ESL_ALPHABET *abc = NULL;
P7_HMMFILE *hfp = NULL;
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_OPROFILE *om = NULL;
P7_GMX *gx1 = NULL;
P7_GMX *gx2 = NULL;
P7_OMX *ox1 = NULL;
P7_OMX *ox2 = NULL;
P7_TRACE *tr = NULL;
int L = esl_opt_GetInteger(go, "-L");
int N = esl_opt_GetInteger(go, "-N");
ESL_DSQ *dsq = malloc(sizeof(ESL_DSQ) * (L+2));
int i;
float fsc, bsc, accscore;
float fsc_g, bsc_g, accscore_g;
double Mcs;
p7_FLogsumInit();
if (p7_hmmfile_OpenE(hmmfile, NULL, &hfp, NULL) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
bg = p7_bg_Create(abc); p7_bg_SetLength(bg, L);
gm = p7_profile_Create(hmm->M, abc); p7_ProfileConfig(hmm, bg, gm, L, p7_LOCAL);
om = p7_oprofile_Create(gm->M, abc); p7_oprofile_Convert(gm, om);
p7_oprofile_ReconfigLength(om, L);
if (esl_opt_GetBoolean(go, "-x") && p7_FLogsumError(-0.4, -0.5) > 0.0001)
p7_Fail("-x here requires p7_Logsum() recompiled in slow exact mode");
ox1 = p7_omx_Create(gm->M, L, L);
ox2 = p7_omx_Create(gm->M, L, L);
tr = p7_trace_CreateWithPP();
esl_rsq_xfIID(r, bg->f, abc->K, L, dsq);
p7_Forward (dsq, L, om, ox1, &fsc);
p7_Backward(dsq, L, om, ox1, ox2, &bsc);
p7_Decoding(om, ox1, ox2, ox2);
esl_stopwatch_Start(w);
for (i = 0; i < N; i++)
{
p7_OptimalAccuracy(om, ox2, ox1, &accscore);
if (! esl_opt_GetBoolean(go, "--notrace"))
{
p7_OATrace(om, ox2, ox1, tr);
p7_trace_Reuse(tr);
}
}
esl_stopwatch_Stop(w);
Mcs = (double) N * (double) L * (double) gm->M * 1e-6 / (double) w->user;
esl_stopwatch_Display(stdout, w, "# CPU time: ");
printf("# M = %d\n", gm->M);
printf("# %.1f Mc/s\n", Mcs);
if (esl_opt_GetBoolean(go, "-c") || esl_opt_GetBoolean(go, "-x") )
{
gx1 = p7_gmx_Create(gm->M, L);
gx2 = p7_gmx_Create(gm->M, L);
p7_GForward (dsq, L, gm, gx1, &fsc_g);
p7_GBackward(dsq, L, gm, gx2, &bsc_g);
p7_GDecoding(gm, gx1, gx2, gx2);
p7_GOptimalAccuracy(gm, gx2, gx1, &accscore_g);
printf("generic: fwd=%8.4f bck=%8.4f acc=%8.4f\n", fsc_g, bsc_g, accscore_g);
printf("VMX: fwd=%8.4f bck=%8.4f acc=%8.4f\n", fsc, bsc, accscore);
p7_gmx_Destroy(gx1);
p7_gmx_Destroy(gx2);
}
free(dsq);
p7_omx_Destroy(ox1);
p7_omx_Destroy(ox2);
p7_trace_Destroy(tr);
p7_oprofile_Destroy(om);
p7_profile_Destroy(gm);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
p7_hmmfile_Close(hfp);
esl_alphabet_Destroy(abc);
esl_stopwatch_Destroy(w);
esl_randomness_Destroy(r);
esl_getopts_Destroy(go);
return 0;
}示例11: serial_loop
static int
serial_loop(WORKER_INFO *info, P7_HMMFILE *hfp)
{
int status;
int i;
int seq_len = 0;
int prev_hit_cnt = 0;
P7_OPROFILE *om = NULL;
P7_SCOREDATA *scoredata = NULL; /* hmm-specific data used by nhmmer */
ESL_ALPHABET *abc = NULL;
P7_DOMAIN *dcl;
#ifdef eslAUGMENT_ALPHABET
ESL_SQ *sq_revcmp = NULL;
if (info->pli->strand != p7_STRAND_TOPONLY && info->qsq->abc->complement != NULL ) {
sq_revcmp = esl_sq_CreateDigital(info->qsq->abc);
esl_sq_Copy(info->qsq,sq_revcmp);
esl_sq_ReverseComplement(sq_revcmp);
info->pli->nres += info->qsq->n;
}
#endif /*eslAUGMENT_ALPHABET*/
/* Main loop: */
while ((status = p7_oprofile_ReadMSV(hfp, &abc, &om)) == eslOK)
{
seq_len = 0;
p7_pli_NewModel(info->pli, om, info->bg);
p7_bg_SetLength(info->bg, info->qsq->n);
p7_oprofile_ReconfigLength(om, info->qsq->n);
scoredata = p7_hmm_ScoreDataCreate(om, FALSE);
#ifdef eslAUGMENT_ALPHABET
//reverse complement
if (info->pli->strand != p7_STRAND_TOPONLY && info->qsq->abc->complement != NULL )
{
p7_Pipeline_LongTarget(info->pli, om, scoredata, info->bg, sq_revcmp, info->th, 0);
p7_pipeline_Reuse(info->pli); // prepare for next search
seq_len = info->qsq->n;
for (i = prev_hit_cnt; i < info->th->N ; i++)
{
dcl = info->th->unsrt[i].dcl;
// modify hit positions to account for the position of the window in the full sequence
dcl->ienv = seq_len - dcl->ienv + 1;
dcl->jenv = seq_len - dcl->jenv + 1;
dcl->iali = seq_len - dcl->iali + 1;
dcl->jali = seq_len - dcl->jali + 1;
dcl->ad->sqfrom = seq_len - dcl->ad->sqfrom + 1;
dcl->ad->sqto = seq_len - dcl->ad->sqto + 1;
}
}
#endif
if (info->pli->strand != p7_STRAND_BOTTOMONLY) {
p7_Pipeline_LongTarget(info->pli, om, scoredata, info->bg, info->qsq, info->th, 0);
p7_pipeline_Reuse(info->pli);
seq_len += info->qsq->n;
}
for (i = prev_hit_cnt; i < info->th->N ; i++)
{
info->th->unsrt[i].lnP += log((float)seq_len / (float)om->max_length);
info->th->unsrt[i].dcl[0].lnP = info->th->unsrt[i].lnP;
info->th->unsrt[i].sortkey = -1.0 * info->th->unsrt[i].lnP;
info->th->unsrt[i].dcl[0].ad->L = om->M;
}
prev_hit_cnt = info->th->N;
p7_oprofile_Destroy(om);
p7_hmm_ScoreDataDestroy(scoredata);
}
esl_alphabet_Destroy(abc);
#ifdef eslAUGMENT_ALPHABET
esl_sq_Destroy(sq_revcmp);
#endif
return status;
}示例12: main
int
main(int argc, char **argv)
{
ESL_GETOPTS *go = esl_getopts_CreateDefaultApp(options, 1, argc, argv, banner, usage);
char *hmmfile = esl_opt_GetArg(go, 1);
ESL_STOPWATCH *w = esl_stopwatch_Create();
ESL_RANDOMNESS *r = esl_randomness_CreateFast(esl_opt_GetInteger(go, "-s"));
ESL_ALPHABET *abc = NULL;
P7_HMMFILE *hfp = NULL;
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_OPROFILE *om = NULL;
P7_OMX *ox = NULL;
P7_GMX *gx = NULL;
int L = esl_opt_GetInteger(go, "-L");
int N = esl_opt_GetInteger(go, "-N");
ESL_DSQ *dsq = malloc(sizeof(ESL_DSQ) * (L+2));
int i;
float sc1, sc2;
double base_time, bench_time, Mcs;
if (p7_hmmfile_Open(hmmfile, NULL, &hfp) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
bg = p7_bg_Create(abc);
p7_bg_SetLength(bg, L);
gm = p7_profile_Create(hmm->M, abc);
p7_ProfileConfig(hmm, bg, gm, L, p7_LOCAL);
om = p7_oprofile_Create(gm->M, abc);
p7_oprofile_Convert(gm, om);
p7_oprofile_ReconfigLength(om, L);
if (esl_opt_GetBoolean(go, "-x")) p7_profile_SameAsVF(om, gm);
ox = p7_omx_Create(gm->M, 0, 0);
gx = p7_gmx_Create(gm->M, L);
/* Get a baseline time: how long it takes just to generate the sequences */
esl_stopwatch_Start(w);
for (i = 0; i < N; i++)
esl_rsq_xfIID(r, bg->f, abc->K, L, dsq);
esl_stopwatch_Stop(w);
base_time = w->user;
/* Run the benchmark */
esl_stopwatch_Start(w);
for (i = 0; i < N; i++)
{
esl_rsq_xfIID(r, bg->f, abc->K, L, dsq);
p7_ViterbiFilter(dsq, L, om, ox, &sc1);
if (esl_opt_GetBoolean(go, "-c"))
{
p7_GViterbi(dsq, L, gm, gx, &sc2);
printf("%.4f %.4f\n", sc1, sc2);
}
if (esl_opt_GetBoolean(go, "-x"))
{
p7_GViterbi(dsq, L, gm, gx, &sc2);
sc2 /= om->scale_w;
if (om->mode == p7_UNILOCAL) sc2 -= 2.0; /* that's ~ L \log \frac{L}{L+2}, for our NN,CC,JJ */
else if (om->mode == p7_LOCAL) sc2 -= 3.0; /* that's ~ L \log \frac{L}{L+3}, for our NN,CC,JJ */
printf("%.4f %.4f\n", sc1, sc2);
}
}
esl_stopwatch_Stop(w);
bench_time = w->user - base_time;
Mcs = (double) N * (double) L * (double) gm->M * 1e-6 / (double) bench_time;
esl_stopwatch_Display(stdout, w, "# CPU time: ");
printf("# M = %d\n", gm->M);
printf("# %.1f Mc/s\n", Mcs);
free(dsq);
p7_omx_Destroy(ox);
p7_gmx_Destroy(gx);
p7_oprofile_Destroy(om);
p7_profile_Destroy(gm);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
p7_hmmfile_Close(hfp);
esl_alphabet_Destroy(abc);
esl_stopwatch_Destroy(w);
esl_randomness_Destroy(r);
esl_getopts_Destroy(go);
return 0;
}示例13: main
int
main(int argc, char **argv)
{
char *hmmfile = argv[1];
char *seqfile = argv[2];
ESL_ALPHABET *abc = NULL;
P7_HMMFILE *hfp = NULL;
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_OPROFILE *om = NULL;
P7_OMX *ox = NULL;
P7_GMX *gx = NULL;
ESL_SQ *sq = NULL;
ESL_SQFILE *sqfp = NULL;
int format = eslSQFILE_UNKNOWN;
float sc;
int status;
/* Read in one HMM */
if (p7_hmmfile_Open(hmmfile, NULL, &hfp) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
/* Read in one sequence */
sq = esl_sq_CreateDigital(abc);
status = esl_sqfile_Open(seqfile, format, NULL, &sqfp);
if (status == eslENOTFOUND) p7_Fail("No such file.");
else if (status == eslEFORMAT) p7_Fail("Format unrecognized.");
else if (status == eslEINVAL) p7_Fail("Can't autodetect stdin or .gz.");
else if (status != eslOK) p7_Fail("Open failed, code %d.", status);
if (esl_sqio_Read(sqfp, sq) != eslOK) p7_Fail("Failed to read sequence");
/* create default null model, then create and optimize profile */
bg = p7_bg_Create(abc);
p7_bg_SetLength(bg, sq->n);
gm = p7_profile_Create(hmm->M, abc);
p7_ProfileConfig(hmm, bg, gm, sq->n, p7_LOCAL);
om = p7_oprofile_Create(gm->M, abc);
p7_oprofile_Convert(gm, om);
p7_oprofile_Logify(om);
/* allocate DP matrices, both a generic and an optimized one */
ox = p7_omx_Create(gm->M, 0, sq->n);
gx = p7_gmx_Create(gm->M, sq->n);
/* Useful to place and compile in for debugging:
p7_oprofile_Dump(stdout, om); dumps the optimized profile
p7_omx_SetDumpMode(ox, TRUE); makes the fast DP algorithms dump their matrices
p7_gmx_Dump(stdout, gx); dumps a generic DP matrix
*/
p7_ViterbiScore(sq->dsq, sq->n, om, ox, &sc); printf("viterbi (non-optimized): %.2f nats\n", sc);
p7_GViterbi (sq->dsq, sq->n, gm, gx, &sc); printf("viterbi (generic): %.2f nats\n", sc);
/* now in a real app, you'd need to convert raw nat scores to final bit
* scores, by subtracting the null model score and rescaling.
*/
/* cleanup */
esl_sq_Destroy(sq);
esl_sqfile_Close(sqfp);
p7_omx_Destroy(ox);
p7_gmx_Destroy(gx);
p7_oprofile_Destroy(om);
p7_profile_Destroy(gm);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
p7_hmmfile_Close(hfp);
esl_alphabet_Destroy(abc);
return 0;
}示例14: pipeline_thread
static void
pipeline_thread(void *arg)
{
int i, j;
int status;
int workeridx;
WORKER_INFO *info;
ESL_THREADS *obj;
P7_OM_BLOCK *block;
void *newBlock;
P7_OPROFILE *om = NULL;
P7_SCOREDATA *scoredata = NULL; /* hmm-specific data used by nhmmer */
P7_DOMAIN *dcl;
int seq_len = 0;
int prev_hit_cnt = 0;
#ifdef eslAUGMENT_ALPHABET
ESL_SQ *sq_revcmp = NULL;
#endif /*eslAUGMENT_ALPHABET*/
impl_ThreadInit();
obj = (ESL_THREADS *) arg;
esl_threads_Started(obj, &workeridx);
info = (WORKER_INFO *) esl_threads_GetData(obj, workeridx);
status = esl_workqueue_WorkerUpdate(info->queue, NULL, &newBlock);
if (status != eslOK) esl_fatal("Work queue worker failed");
#ifdef eslAUGMENT_ALPHABET
//reverse complement
if (info->pli->strand != p7_STRAND_TOPONLY && info->qsq->abc->complement != NULL ) {
sq_revcmp = esl_sq_CreateDigital(info->qsq->abc);
esl_sq_Copy(info->qsq,sq_revcmp);
esl_sq_ReverseComplement(sq_revcmp);
info->pli->nres += info->qsq->n;
}
#endif /*eslAUGMENT_ALPHABET*/
/* loop until all blocks have been processed */
block = (P7_OM_BLOCK *) newBlock;
while (block->count > 0)
{
/* Main loop: */
for (i = 0; i < block->count; ++i)
{
om = block->list[i];
seq_len = 0;
p7_pli_NewModel(info->pli, om, info->bg);
p7_bg_SetLength(info->bg, info->qsq->n);
p7_oprofile_ReconfigLength(om, info->qsq->n);
scoredata = p7_hmm_ScoreDataCreate(om, FALSE);
#ifdef eslAUGMENT_ALPHABET
//reverse complement
if (info->pli->strand != p7_STRAND_TOPONLY && info->qsq->abc->complement != NULL )
{
p7_Pipeline_LongTarget(info->pli, om, scoredata, info->bg, sq_revcmp, info->th, 0);
p7_pipeline_Reuse(info->pli); // prepare for next search
seq_len = info->qsq->n;
for (j = prev_hit_cnt; j < info->th->N ; j++)
{
dcl = info->th->unsrt[j].dcl;
// modify hit positions to account for the position of the window in the full sequence
dcl->ienv = seq_len - dcl->ienv + 1;
dcl->jenv = seq_len - dcl->jenv + 1;
dcl->iali = seq_len - dcl->iali + 1;
dcl->jali = seq_len - dcl->jali + 1;
dcl->ad->sqfrom = seq_len - dcl->ad->sqfrom + 1;
dcl->ad->sqto = seq_len - dcl->ad->sqto + 1;
}
}
#endif
if (info->pli->strand != p7_STRAND_BOTTOMONLY) {
p7_Pipeline_LongTarget(info->pli, om, scoredata, info->bg, info->qsq, info->th, 0);
p7_pipeline_Reuse(info->pli);
seq_len += info->qsq->n;
}
for (j = prev_hit_cnt; j < info->th->N ; j++)
{
info->th->unsrt[j].lnP += log((float)seq_len / (float)om->max_length);
info->th->unsrt[j].dcl[0].lnP = info->th->unsrt[j].lnP;
info->th->unsrt[j].sortkey = -1.0 * info->th->unsrt[j].lnP;
info->th->unsrt[j].dcl[0].ad->L = om->M;
}
prev_hit_cnt = info->th->N;
p7_hmm_ScoreDataDestroy(scoredata);
p7_oprofile_Destroy(om);
block->list[i] = NULL;
//.........这里部分代码省略.........
示例15: main
int
main(int argc, char **argv)
{
ESL_GETOPTS *go = p7_CreateDefaultApp(options, 2, argc, argv, banner, usage);
char *hmmfile = esl_opt_GetArg(go, 1);
char *seqfile = esl_opt_GetArg(go, 2);
ESL_ALPHABET *abc = NULL;
P7_HMMFILE *hfp = NULL;
P7_HMM *hmm = NULL;
P7_BG *bg = NULL;
P7_PROFILE *gm = NULL;
P7_GMX *gx1 = NULL;
P7_GMX *gx2 = NULL;
ESL_SQ *sq = NULL;
ESL_SQFILE *sqfp = NULL;
P7_TRACE *tr = NULL;
int format = eslSQFILE_UNKNOWN;
char errbuf[eslERRBUFSIZE];
float fsc, bsc, vsc;
float accscore;
int status;
/* Read in one HMM */
if (p7_hmmfile_OpenE(hmmfile, NULL, &hfp, NULL) != eslOK) p7_Fail("Failed to open HMM file %s", hmmfile);
if (p7_hmmfile_Read(hfp, &abc, &hmm) != eslOK) p7_Fail("Failed to read HMM");
p7_hmmfile_Close(hfp);
/* Read in one sequence */
sq = esl_sq_CreateDigital(abc);
status = esl_sqfile_OpenDigital(abc, seqfile, format, NULL, &sqfp);
if (status == eslENOTFOUND) p7_Fail("No such file.");
else if (status == eslEFORMAT) p7_Fail("Format unrecognized.");
else if (status == eslEINVAL) p7_Fail("Can't autodetect stdin or .gz.");
else if (status != eslOK) p7_Fail("Open failed, code %d.", status);
if (esl_sqio_Read(sqfp, sq) != eslOK) p7_Fail("Failed to read sequence");
esl_sqfile_Close(sqfp);
/* Configure a profile from the HMM */
bg = p7_bg_Create(abc);
p7_bg_SetLength(bg, sq->n);
gm = p7_profile_Create(hmm->M, abc);
p7_ProfileConfig(hmm, bg, gm, sq->n, p7_LOCAL); /* multihit local: H3 default */
/* Allocations */
gx1 = p7_gmx_Create(gm->M, sq->n);
gx2 = p7_gmx_Create(gm->M, sq->n);
tr = p7_trace_CreateWithPP();
p7_FLogsumInit();
/* Run Forward, Backward; do OA fill and trace */
p7_GForward (sq->dsq, sq->n, gm, gx1, &fsc);
p7_GBackward(sq->dsq, sq->n, gm, gx2, &bsc);
p7_GDecoding(gm, gx1, gx2, gx2); /* <gx2> is now the posterior decoding matrix */
p7_GOptimalAccuracy(gm, gx2, gx1, &accscore); /* <gx1> is now the OA matrix */
p7_GOATrace(gm, gx2, gx1, tr);
if (esl_opt_GetBoolean(go, "-d")) p7_gmx_Dump(stdout, gx2, p7_DEFAULT);
if (esl_opt_GetBoolean(go, "-m")) p7_gmx_Dump(stdout, gx1, p7_DEFAULT);
p7_trace_Dump(stdout, tr, gm, sq->dsq);
if (p7_trace_Validate(tr, abc, sq->dsq, errbuf) != eslOK) p7_Die("trace fails validation:\n%s\n", errbuf);
printf("fwd = %.4f nats\n", fsc);
printf("bck = %.4f nats\n", bsc);
printf("acc = %.4f (%.2f%%)\n", accscore, accscore * 100. / (float) sq->n);
p7_trace_Reuse(tr);
p7_GViterbi(sq->dsq, sq->n, gm, gx1, &vsc);
p7_GTrace (sq->dsq, sq->n, gm, gx1, tr);
p7_trace_SetPP(tr, gx2);
p7_trace_Dump(stdout, tr, gm, sq->dsq);
printf("vit = %.4f nats\n", vsc);
printf("acc = %.4f\n", p7_trace_GetExpectedAccuracy(tr));
/* Cleanup */
esl_sq_Destroy(sq);
p7_trace_Destroy(tr);
p7_gmx_Destroy(gx1);
p7_gmx_Destroy(gx2);
p7_profile_Destroy(gm);
p7_bg_Destroy(bg);
p7_hmm_Destroy(hmm);
esl_alphabet_Destroy(abc);
esl_getopts_Destroy(go);
return 0;
}