From 6d8376bfd915df64ff4fe804bd5e9152ebd4d5d2 Mon Sep 17 00:00:00 2001 From: Travis Date: Sat, 5 Sep 2026 00:02:06 -0700 Subject: [PATCH] Remove unused contrib/, libdivsufsort/, and .travis.yml contrib/ held Lee Newberg's p7_thermo code, parked there in HMMER 3.0b1 and never compiled by any Makefile. libdivsufsort/ was nhmmer's suffix-array library; BATH removed nhmmer, and nothing in src/ includes or calls it, but every binary still linked against it. Removing it also drops its wiring from configure.ac and the Makefiles (rerun autoconf after pulling this). .travis.yml described HMMER's Travis setup; CI now lives in .github/workflows/release.yml. Also replace the HMMER- and frahmmer-era binary names in .gitignore with the bath* programs, and ignore the configure-generated titlepage_daemon.tex. Implemented with assistance from Claude (Fable 5) --- .gitignore | 38 +- .travis.yml | 14 - Makefile.in | 7 - configure.ac | 4 - contrib/p7_thermo.c | 930 ----------------- contrib/p7_thermo.h | 37 - libdivsufsort/AUTHORS | 3 - libdivsufsort/COPYING | 27 - libdivsufsort/Makefile.in | 54 - libdivsufsort/README | 185 ---- libdivsufsort/VERSION | 1 - libdivsufsort/divsufsort.c | 1794 --------------------------------- libdivsufsort/divsufsort.h.in | 63 -- src/Makefile.in | 6 +- 14 files changed, 9 insertions(+), 3154 deletions(-) delete mode 100644 .travis.yml delete mode 100644 contrib/p7_thermo.c delete mode 100644 contrib/p7_thermo.h delete mode 100644 libdivsufsort/AUTHORS delete mode 100644 libdivsufsort/COPYING delete mode 100644 libdivsufsort/Makefile.in delete mode 100644 libdivsufsort/README delete mode 100644 libdivsufsort/VERSION delete mode 100644 libdivsufsort/divsufsort.c delete mode 100644 libdivsufsort/divsufsort.h.in diff --git a/.gitignore b/.gitignore index 0bb9cf89..6c7f7ec5 100644 --- a/.gitignore +++ b/.gitignore @@ -18,38 +18,14 @@ src/impl*/*_utest src/impl*/*_benchmark documentation/userguide/copyright.tex documentation/userguide/titlepage.tex +documentation/userguide/titlepage_daemon.tex documentation/man/*.man -libdivsufsort/divsufsort.h -libdivsufsort/libdivsufsort.a src/libhmmer.a src/p7_config.h -profmark/create-profmark -profmark/rocplot -src/alimask -src/hmmalign -src/hmmbuild -src/hmmc2 -src/hmmconvert -src/hmmemit -src/hmmerfm-exactmatch -src/hmmfetch -src/hmmlogo -src/hmmpgmd -src/hmmpgmd_shard -src/hmmpress -src/hmmscan -src/hmmsearch -src/hmmsim -src/hmmstat src/itest_brute -src/jackhmmer -src/makehmmerdb -src/nhmmer -src/nhmmscan -src/phmmer -src/frahmmbuild -src/frahmmconvert -src/frahmmemit -src/frahmmer -src/frahmmfetch -src/frahmmstat +src/bathalign +src/bathbuild +src/bathconvert +src/bathfetch +src/bathsearch +src/bathstat diff --git a/.travis.yml b/.travis.yml deleted file mode 100644 index f84f0f04..00000000 --- a/.travis.yml +++ /dev/null @@ -1,14 +0,0 @@ -language: c - -env: - - SQC_NONZERO_EXIT=1 - -script: - - git clone -b develop https://github.com/${TRAVIS_REPO_SLUG/hmmer/easel}.git - - ln -s easel/aclocal.m4 aclocal.m4 - - autoconf - - ./configure - - make - - make dev - - make check - diff --git a/Makefile.in b/Makefile.in index dd2c8abc..49a95509 100644 --- a/Makefile.in +++ b/Makefile.in @@ -26,8 +26,6 @@ SHELL = /bin/sh # location of easel ESLDIR = @HMMER_ESLDIR@ -# location of libdivsufsort for suffix array creation -SADIR = @HMMER_SADIR@ # Package information @@ -88,7 +86,6 @@ endif # all: ${QUIET_SUBDIR0}${ESLDIR} ${QUIET_SUBDIR1} all - ${QUIET_SUBDIR0}${SADIR} ${QUIET_SUBDIR1} all ${QUIET_SUBDIR0}src ${QUIET_SUBDIR1} all # ${QUIET_SUBDIR0}profmark ${QUIET_SUBDIR1} all @@ -96,7 +93,6 @@ all: # dev: ${QUIET_SUBDIR0}${ESLDIR} ${QUIET_SUBDIR1} dev - ${QUIET_SUBDIR0}${SADIR} ${QUIET_SUBDIR1} all ${QUIET_SUBDIR0}src ${QUIET_SUBDIR1} dev # ${QUIET_SUBDIR0}profmark ${QUIET_SUBDIR1} dev @@ -111,7 +107,6 @@ tests: check: @command -v python3 >/dev/null 2>&1 || { echo >&2 "python3 is required for 'make check', but is not in your PATH. Aborting."; exit 1; } ${QUIET_SUBDIR0}${ESLDIR} ${QUIET_SUBDIR1} tests - ${QUIET_SUBDIR0}${SADIR} ${QUIET_SUBDIR1} all ${QUIET_SUBDIR0}src ${QUIET_SUBDIR1} tests ${QUIET_SUBDIR0}${ESLDIR} ${QUIET_SUBDIR1} check ${QUIET_SUBDIR0}testsuite ${QUIET_SUBDIR1} check @@ -157,7 +152,6 @@ clean: ${QUIET_SUBDIR0}testsuite ${QUIET_SUBDIR1} clean ${QUIET_SUBDIR0}documentation ${QUIET_SUBDIR1} clean ${QUIET_SUBDIR0}${ESLDIR} ${QUIET_SUBDIR1} clean - ${QUIET_SUBDIR0}${SADIR} ${QUIET_SUBDIR1} clean ${QUIET}-rm -f *.o *~ Makefile.bak core TAGS TAGS.part gmon.out ifndef V @echo ' ' CLEAN hmmer @@ -171,7 +165,6 @@ distclean: ${QUIET_SUBDIR0}testsuite ${QUIET_SUBDIR1} distclean ${QUIET_SUBDIR0}documentation ${QUIET_SUBDIR1} distclean ${QUIET_SUBDIR0}${ESLDIR} ${QUIET_SUBDIR1} distclean - ${QUIET_SUBDIR0}${SADIR} ${QUIET_SUBDIR1} distclean ${QUIET}-rm -f config.log config.status ${QUIET}-rm -rf autom4te.cache ${QUIET}-rm -f *.o *~ Makefile.bak core TAGS TAGS.part gmon.out diff --git a/configure.ac b/configure.ac index 829e5601..23861804 100644 --- a/configure.ac +++ b/configure.ac @@ -120,7 +120,6 @@ BATH_VERSION="2.0" BATH_URL="https://github.com/TravisWheelerLab/BATH" HMMER_ESLDIR="easel" -HMMER_SADIR="libdivsufsort" EASEL_DATE="Nov 2020" EASEL_COPYRIGHT="Copyright (C) 2020 Howard Hughes Medical Institute." @@ -139,7 +138,6 @@ AC_SUBST(BATH_VERSION) AC_SUBST(BATH_URL) AC_SUBST(HMMER_ESLDIR) -AC_SUBST(HMMER_SADIR) AC_SUBST(EASEL_DATE) AC_SUBST(EASEL_COPYRIGHT) @@ -669,7 +667,6 @@ AC_CONFIG_FILES([ \ documentation/man/Makefile \ documentation/userguide/Makefile \ documentation/userguide/inclusions/Makefile \ - libdivsufsort/Makefile \ ]) # Easel Makefiles. @@ -729,7 +726,6 @@ AC_CONFIG_FILES([ \ AC_CONFIG_HEADERS([easel/decoy_config.h]) # Put this first to keep gnu 'autoheader' from overwriting a real config file. AC_CONFIG_HEADERS([src/p7_config.h]) AC_CONFIG_HEADERS([easel/esl_config.h]) -AC_CONFIG_HEADERS([libdivsufsort/divsufsort.h]) # the following incantation establishes a symlink of # src/impl_{whatever} to src/impl in the *build* directory. diff --git a/contrib/p7_thermo.c b/contrib/p7_thermo.c deleted file mode 100644 index 8806dfc4..00000000 --- a/contrib/p7_thermo.c +++ /dev/null @@ -1,930 +0,0 @@ -#include "hmmer.h" -#include - -/* Check that esl exceptions, esl_fatal, -eslINFINITY, etc. are used - in the correct contexts!! */ - -/* Note that all calculations assume L is fixed. That is, for these - p-value and sensitivity calculations, a score for a sequence is - compared to the scores of sequences of the same length only. Is - this okay?!! */ - -/* Function: utest_thermo() - * Synopsis: Test error rate calculations - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - */ - -static void -processEstimates(const float numSamplesM1, float v1Up, float v1Dn, float *v1Bt, float v2Up, float v2Dn, float *v2Bt); - -/* printMantissa(x): For printing the m.mm part of a m.mmEx.xx - floating point number, exp(x)*/ -static float -printMantissa(const float x); - -/* printExponent(x): For printing the xxx part of a m.mmExxx floating - point number, exp(x)*/ -static int -printExponent(const float x); - -void -utest_thermo(ESL_GETOPTS *go, ESL_RANDOMNESS *r, const ESL_ALPHABET *abc, const P7_BG *bg, const P7_PROFILE *gm, int nseq, int L) -{ - ESL_DSQ *dsq = NULL; - P7_GMX *gx = NULL; - P7_TRACE *tr = NULL; - float *scores[p7M_NSCORETYPES]; - P7_THERMO *thermos[p7M_NSCORETYPES]; /* all the thermos */ - P7_THERMO *thermo; /* the thermo of current interest */ - int idx; - int support; /* number of nonzero terms in importance sampling sum */ - float pv; /* p-value */ - float pvsp; /* p-value or specificity, whichever is smaller */ - float pvstd; /* standard deviation of pv (and also of sp) */ - float sn; /* sensitivity */ - float snfn; /* sensitivity or false negative rate, whichever is smaller */ - float snstd; /* standard deviation of sn (and also of fnr) */ - float stretch; /* For generating interesting scores to test */ - float sc = 0; - enum p7m_scoretypes_e p7m_score; - const int numSamples = 100; /* Should be a parameter?!! */ - - for (p7m_score = 0; p7m_score < p7M_NSCORETYPES; p7m_score++) { - scores[p7m_score] = NULL; - thermos[p7m_score] = NULL; - } - - if (esl_opt_GetBoolean(go, "--vv")) { - printf("utest_thermo: M = %d, L = %d, nseq = %d, numSamples = %d. Generating typical scores for this profile-HMM...\n", gm->M, L, nseq, numSamples); - } - - if ((dsq = malloc(sizeof(ESL_DSQ) *(L+2))) == NULL) esl_fatal("dsq malloc failed"); - if ((gx = p7_gmx_Create(gm->M, L)) == NULL) esl_fatal("matrix creation failed"); - for (p7m_score = 0; p7m_score < p7M_NSCORETYPES; p7m_score++) { - if ((scores[p7m_score] = malloc(sizeof(float) * nseq)) == NULL) esl_fatal("scores malloc failed"); - } - - /* Generate some typical scores */ - for (idx = 0; idx < nseq; idx++) { - if (esl_rnd_xfIID(r, bg->f, abc->K, L, dsq) != eslOK) esl_fatal("seq generation failed"); - if (p7_GForward(dsq, L, gm, gx, &scores[p7M_FORWARD][idx]) != eslOK) esl_fatal("forward failed"); - if (p7_GViterbi(dsq, L, gm, gx, &scores[p7M_VITERBI][idx]) != eslOK) esl_fatal("viterbi failed"); - } - for (p7m_score = 0; p7m_score < p7M_NSCORETYPES; p7m_score++) { - esl_vec_FSortIncreasing(scores[p7m_score], nseq); - } - - /* Loop through the score types */ - for (p7m_score = 0; p7m_score < p7M_NSCORETYPES; p7m_score++) { - if (esl_opt_GetBoolean(go, "--vv")) { - printf("utest_thermo: Calibrating temperatures ...\n"); - } - /* Build and calibrate a P7_THERMO structure */ - if ((thermo = p7_thermo_Create()) == NULL) esl_fatal("thermo creation failed"); - thermos[p7m_score] = thermo; - /* p7_thermoCalibrate() is worth it only if we're going to - evaluate several different score thresholds, which we are - going to do. */ - if (p7_thermoCalibrate(thermo, r, p7m_score, bg, gm, gx) != eslOK) esl_fatal("thermo calibration failed"); - - if (esl_opt_GetBoolean(go, "--vv")) { - printf("utest_thermo: p7_thermo[%s] (temperature, score) calibration curve:", (p7m_score == p7M_FORWARD ? "forward" : "viterbi")); - for (idx = 0; idx < thermo->numTemperatures; idx++) { - printf(" (%8.4f, %8.4f)", thermo->temperatures[idx], thermo->scores[idx]); - } - printf("\n"); - } - - /* Use the current P7_THERMO structure to guide evaluation - of the sequences we generated. */ - if (esl_opt_GetBoolean(go, "--vv")) { - printf("utest_thermo: ROC curve ...\n"); - } - for (idx = 0; idx < nseq; idx++) { - /* We'll stretch sc's distance from the minimum score - we've seen, up to the maximum score we've seen, just to - be interesting */ - stretch = (float) idx / (float) (nseq - 1); - stretch *= stretch; - stretch *= stretch; - sc = stretch * (thermo->scores[0] - scores[p7m_score][0]) + scores[p7m_score][0]; - - if (p7_thermoEstimate(r, thermo, bg, gm, gx, numSamples, sc, &support, &pv, &pvstd, &sn, &snstd) != eslOK) esl_fatal("p7_thermoEstimate failed"); - if (esl_opt_GetBoolean(go, "--vv")) { - /* If expf(pv) <= 0.5 display it; otherwise display - specificity =1-pv. Likewise, if expf(sn) <= 0.5 - display it; otherwise display false negative rate - =1-sn. Note that sp=logf(1.0f-expf(pv)) is - log(-pv) + pv/2 +pv^2/24 + ... when pv is near - zero. */ - if (pv <= logf(0.5f)) pvsp = pv; /* pv */ - else if (pv <= -1e-10) pvsp = logf(1.0f - expf(pv)); /* sp */ - else pvsp = logf(-pv) + 0.5f*pv; /* sp */ - if (sn <= logf(0.5f)) snfn = sn; /* sn */ - else if (sn <= -1e-10) snfn = logf(1.0f - expf(sn)); /* fnr */ - else snfn = logf(-sn) + 0.5*sn; /* fnr */ - printf("utest_thermo: %s threshold %9.4f: %s = %5.3fe%+04d +- %3.0f%%, %s = %5.3fe%+04d +- %3.0f%%, support = %d/%d\n", - (p7m_score == p7M_FORWARD ? "Forward" : "Viterbi"), - sc, - pv > log(0.5f) ? "specificity" : " p-value", - printMantissa(pvsp), - printExponent(pvsp), - 100.0f * expf(pvstd-pvsp), - sn > log(0.5f) ? "f-neg. rate" : "sensitivity", - printMantissa(snfn), - printExponent(snfn), - 100.0f * expf(snstd-snfn), - support, - numSamples); - } - } - } - - if (dsq) free(dsq); - p7_gmx_Destroy(gx); - p7_trace_Destroy(tr); - for (p7m_score = 0; p7m_score < p7M_NSCORETYPES; p7m_score++) { - if (scores[p7m_score]) free(scores[p7m_score]); - p7_thermo_Destroy(thermos[p7m_score]); - } - return; -} - -/* Function: p7_thermoEstimate() - * Synopsis: Estimates both the p-value and sensitivity of a score threshold. - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: - * For any forward score threshold (or likewise for viterbi - * scores) we may be interested in what fraction of - * sequences of length L, drawn from the background model - * would have that score or higher; this is called - * p-value. Likewise we may be interested in what fraction - * of sequences of length L, drawn from the foreground - * model would have that score or higher; this is called - * sensitivity. Here we are calculating p-value and - * sensitivity. - * - * We must first create and configure , , , - * , and . In particular, includes - * whether the score is forward or viterbi. - * is the number of samples drawn from the - * importance sampling distribution, and confidence limits - * are inversely proportional to its square root. - * - * Returns: is the logf of the computed p-value. is - * the logf of the standard deviation of the computed . - * is the logf of the computed sensitivity. - * is the logf of the standard deviation of the computed - * . is the number of the that - * contributed a non-zero value to the importance sampling - * sum. As a general rule for p-values under 0.01, at - * least 20 samples should *fail* to be non-zero, otherwise - * we worry that the region just above the threshold may - * not be adequately represented in the importance sampling - * sum. - */ - -int -p7_thermoEstimate(ESL_RANDOMNESS *r, const P7_THERMO *thermo, const P7_BG *bg, const P7_PROFILE *gm, P7_GMX *gx, int numSamples, float threshold, int *support, float *pv, float *pvstd, float *sn, float *snstd) -{ - P7_PROFILE *gmT = NULL; /* For sampling a sequence at a specified temperature. */ - P7_PROFILE *gmDT = NULL; /* For evaluating at a specified temperature. */ - P7_GMX *gxT = NULL; - P7_GMX *gxDT = NULL; - ESL_DSQ *dsqX = NULL; /* Sequence of all X/N */ - ESL_DSQ *dsq = NULL; /* Sampled sequence */ - P7_TRACE *tr = NULL; /* Sampled path through states */ - float *pvUp = NULL; /* pv contributions at threshold or higher */ - float *pvDn = NULL; /* pv contributions below threshold */ - float *snUp = NULL; /* sn contributions at threshold or higher */ - float *snDn = NULL; /* sn contributions below threshold */ - float temperature = 0.0f; - float sc = 0.0f; - float ZT = 0.0f; /* Forward sum of gmT(dsqX) */ - float Z1 = 0.0f; /* Forward sum of gmT(dsqX) for T = 1 */ - float ZDT = 0.0f; /* Forward sum of gmDT(dsq) */ - float pv1Up = 0.0f; /* "mean" of pv samples at or above threshold */ - float pv1Dn = 0.0f; /* "mean" of pv samples below threshold */ - float pv1Bt = 0.0f; /* mean of all pv samples */ - float pv2Up = 0.0f; /* "variance" of pv samples at or above threshold */ - float pv2Dn = 0.0f; /* "variance" of pv samples below threshold */ - float pv2Bt = 0.0f; /* variance of all pv samples */ - float sn1Up = 0.0f; /* "mean" of sn samples at or above threshold */ - float sn1Dn = 0.0f; /* "mean" of sn samples below threshold */ - float sn1Bt = 0.0f; /* mean of all sn samples */ - float sn2Up = 0.0f; /* "variance" of sn samples at or above threshold */ - float sn2Dn = 0.0f; /* "variance" of sn samples below threshold */ - float sn2Bt = 0.0f; /* variance of all sn samples */ - const int L = gm->L; - const int deg = gm->abc->Kp-2; /* maximally degenerate character (X/N) */ - const int p7m_score = thermo->p7m_score; - const float numSamplesLog = logf((float) numSamples); - int z; - int status; - - if ((gmT = p7_profile_Clone(gm)) == NULL) esl_fatal("failed to create gmT"); - if ((gmDT = p7_profile_Clone(gm)) == NULL) esl_fatal("failed to create gmDT"); - if ((gxT = p7_gmx_Create(gm->M, L)) == NULL) esl_fatal("failed to create gxT"); - if ((gxDT = p7_gmx_Create(gm->M, L)) == NULL) esl_fatal("failed to create gxDT"); - if ((dsqX = malloc(sizeof(ESL_DSQ) * (L+2))) == NULL) esl_fatal("failed to create dsqX"); - if ((dsq = malloc(sizeof(ESL_DSQ) * (L+2))) == NULL) esl_fatal("failed to create dsq"); - if ((tr = p7_trace_Create()) == NULL) esl_fatal("failed to create tr"); - if ((pvUp = malloc(sizeof(float) * numSamples)) == NULL) esl_fatal("failed to create samples"); - if ((pvDn = malloc(sizeof(float) * numSamples)) == NULL) esl_fatal("failed to create samples"); - if ((snUp = malloc(sizeof(float) * numSamples)) == NULL) esl_fatal("failed to create samples"); - if ((snDn = malloc(sizeof(float) * numSamples)) == NULL) esl_fatal("failed to create samples"); - - /* Create string of all degenerate characters */ - for (z = 1; z <= L; z++) dsqX[z] = deg; /* The X/N character */ - - /* Need to normalize sensitivity value by the foreground model - probability of length . Note that is the same even for - subsequent calles to p7_thermoEstimate(), even if the score - threshold is significanctly different; should be - cached, perhaps in !! */ - p7_profileAdjustClones(1.0, bg, gm, gmT, gmDT); - p7_GForward(dsqX, L, gmT, gxT, &Z1); - - /* Choose a temperature for efficient estimation. */ - p7_thermoSuggestTemp(thermo, threshold, &temperature); - - /* Make two profiles at this temperature. will be used to - evaluate . will be used to evaluate . */ - p7_profileAdjustClones(temperature, bg, gm, gmT, gmDT); - /* Fill up gxT matrix with useful values for a later trace. */ - p7_GForward(dsqX, L, gmT, gxT, &ZT); - - /* Note that the calculated 1/ZDT scores should be cached (perhaps - in ) for later calls to p7_thermoEstimate, because they - can be re-used if a subsequent call is for score threshold - nearby to the current one. Under these circumstances, this - would save nearly 100% of the run time for subsequent calls to - p7_thermoEstimate!! */ - int nonZeros = 0; - for (z = 0; z < numSamples; z++) { - /* Trace back a set of states and emissions */ - p7_StochasticDsqTrace(r, dsqX, L, bg, gmT, gxT, tr, dsq); - /* Evaluate sampled sequence of emissions */ - switch(p7m_score) { - case p7M_FORWARD: - p7_GForward(dsq, L, gm, gx, &sc); - break; - case p7M_VITERBI: - p7_GViterbi(dsq, L, gm, gx, &sc); - break; - } - p7_GForward(dsq, L, gmDT, gxDT, &ZDT); - if (sc >= threshold) { - if (p7m_score != p7M_FORWARD) p7_GForward(dsq, L, gm, gx, &sc); /* needed for sensitivities only */ - nonZeros++; - pvUp[z] = ZT - ZDT; /* Contribution to pv importance sampling sum */ - pvDn[z] = -eslINFINITY; - snUp[z] = ZT - ZDT + sc - Z1; /* Contribution to sn importance sampling sum */ - snDn[z] = -eslINFINITY; - } else { - if (p7m_score != p7M_FORWARD) p7_GForward(dsq, L, gm, gx, &sc); /* needed for sensitivities only */ - pvDn[z] = ZT - ZDT; /* Contribution to pv Importance Sampling sum */ - pvUp[z] = -eslINFINITY; - snDn[z] = ZT - ZDT + sc - Z1; /* Contribution to sn Importance Sampling sum */ - snUp[z] = -eslINFINITY; - } - } - - /* Compute means. */ - pv1Up = esl_vec_FLogSum(pvUp, numSamples) - numSamplesLog; - pv1Dn = esl_vec_FLogSum(pvDn, numSamples) - numSamplesLog; - sn1Up = esl_vec_FLogSum(snUp, numSamples) - numSamplesLog; - sn1Dn = esl_vec_FLogSum(snDn, numSamples) - numSamplesLog; - - /* Compute expected squares. */ - esl_vec_FScale(pvUp, numSamples, 2.0f); - esl_vec_FScale(pvDn, numSamples, 2.0f); - esl_vec_FScale(snUp, numSamples, 2.0f); - esl_vec_FScale(snDn, numSamples, 2.0f); - pv2Up = esl_vec_FLogSum(pvUp, numSamples) - numSamplesLog; - pv2Dn = esl_vec_FLogSum(pvDn, numSamples) - numSamplesLog; - sn2Up = esl_vec_FLogSum(snUp, numSamples) - numSamplesLog; - sn2Dn = esl_vec_FLogSum(snDn, numSamples) - numSamplesLog; - - /* Process from the perspective of the Up data or the Down data, - which ever gives the smaller mean. However, don't use the Down - data if nonZeros equals numSamples. At the other extreme, when - nonZeros equals zero we are toast; we should try a lower - temperature, but instead we'll just "compute" a p-value and - sensitivity of exactly zero. */ - if (pv1Up < pv1Dn || nonZeros == numSamples) { - processEstimates(numSamples - 1.0f, pv1Up, pv1Dn, &pv1Bt, pv2Up, pv2Dn, &pv2Bt); - *pvstd = 0.5f * pv2Bt; - *pv = pv1Bt; - } else { - processEstimates(numSamples - 1.0f, pv1Dn, pv1Up, &pv1Bt, pv2Dn, pv2Up, &pv2Bt); - *pvstd = 0.5f * pv2Bt; - /* restore original perspective */ - if (pv1Bt < -1e-10) *pv = logf(1.0f - expf(pv1Bt)); - else *pv = logf(-pv1Bt) + 0.5 * pv1Bt; - } - if (sn1Up < sn1Dn || nonZeros == numSamples) { - processEstimates(numSamples - 1.0f, sn1Up, sn1Dn, &sn1Bt, sn2Up, sn2Dn, &sn2Bt); - *snstd = 0.5f * sn2Bt; - *sn = sn1Bt; - } else { - processEstimates(numSamples - 1.0f, sn1Dn, sn1Up, &sn1Bt, sn2Dn, sn2Up, &sn2Bt); - *snstd = 0.5f * sn2Bt; - /* restore original perspective */ - if (sn1Bt < -1e-10) *sn = logf(1.0f - expf(sn1Bt)); - else *sn = logf(-sn1Bt) + 0.5 * sn1Bt; - } - - /* Report to the user */ - *support = nonZeros; - status = eslOK; - - /* Clean up */ - p7_profile_Destroy(gmT); - p7_profile_Destroy(gmDT); - p7_gmx_Destroy(gxT); - p7_gmx_Destroy(gxDT); - if (dsqX) free (dsqX); - if (dsq) free (dsq); - p7_trace_Destroy(tr); - if (pvUp) free (pvUp); - if (pvDn) free (pvDn); - if (snUp) free (snUp); - if (snDn) free (snDn); - - return status; -} - - -/* Function: processEstimates() - - * Synopsis: File static function that estimates p-value (or - * sensitivity) from summary statistics from importance - * samples. - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: expf(v1Up) is the importance sampling sum. expf(v1Dn) is - * the sum over the rejected samples, which should come to - * approximately 1.0-. expf(v2Up) and expf(v2Dn) are - * the sum of sampled value squares for the accepted and - * rejected samples, respectively; they are useful for - * computing the variances of exp(v1Up) and exp(v1Dn), - * respectively. - * - * Returns: expf(*v1Bt) is the estimate for p-value (or sensitivity). - * expf(*v2Bt) is the variance of that estimate. - */ - -void -processEstimates(const float numSamplesM1, float v1Up, float v1Dn, float *v1Bt, float v2Up, float v2Dn, float *v2Bt) -{ - /* This could be a fancy function that gleans precision from a - proper combination of v1Up and v1Dn, but instead we just use - v1Up, knowing only that it is smaller than v1Dn. */ - - /* Convert sums of squares to variances of means */ - if ((v2Up > 2.0f * v1Up) && (numSamplesM1 > 0.0f)) - v2Up += logf((1.0f - expf(2.0f * v1Up - v2Up)) / numSamplesM1); - else v2Up = -eslINFINITY; - - *v1Bt = v1Up; - *v2Bt = v2Up; - -#if 0 - /* If we wanted to try to get fancy, instead of: - - expf(v1Up) - - we could use - - 1-expf(v1Dn), or - expf(v1Up)/(expf(v1Up)+expf(v1Dn), or - - similar, or some linear combination of the above. For added - glamour we could have a linear combination of, say, expf(v1Up) - and 1-expf(v1Dn) where the relative weights for combining these - two values depend upon v2Up and v2Dn ... or v1Up and v1Dn - themselves. - - If we are to estimate the variance of a function of both - expf(v1Up) and expf(v1Dn), then we may need to know their - covariance. Fortunately we can estimate this via the variance - of the whole set of sampled values, as follows: - */ - - /* Combine vUp and vDn data sets */ - *v1Bt = v1Dn + logf(1.0f + expf(v1Up - v1Dn)); /* Sum of sums */ - *v2Bt = v2Dn + logf(1.0f + expf(v2Up - v2Dn)); /* Sum of sums of squares */ - - /* Compute variances of the data sets separately (v2Up, v2Dn) and - together (v2Bt). (v2Up is already done, above.) */ - if ((v2Dn > 2.0f * v1Dn) && (numSamplesM1 > 0.0f)) - v2Dn += logf((1.0f - expf(2.0f * v1Dn - v2Dn)) / numSamplesM1); - else v2Dn = -eslINFINITY; - if ((*v2Bt > 2.0f * *v1Bt) && (numSamplesM1 > 0.0f)) - *v2Bt += logf((1.0f - expf(2.0f * *v1Bt - *v2Bt)) / numSamplesM1); - else *v2Bt = -eslINFINITY; - - /* Covariance between expf(v1Up) and expf(v1Dn) is negative - because each sample sequence that contributes to vUp will - be one that does not contribute to vDn, and vice versa. It - turns out that the logarithm of its negative can be estimated - by: */ - - float negativeCovariance = (v2Up > v2Dn ? v2Up + logf(1.0f + expf(v2Dn - v2Up)) : v2Dn + logf(1.0f + expf(v2Up - v2Dn))); - negativeCovariance += logf(0.5f * (1.0f - expf(*v2Bt - negativeCovariance))); -#endif - -} - -/* Function: p7_profileAdjustClones() - * Synopsis: Configure two profiles for use with the - * p-value calculations. - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: A profile , recently cloned from a profile , - * is modified to reflect the temperature. In concrete - * terms: all scores of are divided by - * and put in . A profile , also recently - * cloned from a profile , is then built from the - * scores of ; the expf(scores) for each emission - * are averaged according to the background model prior - * distribution. The average is stored in the location for - * the totally degenerate character X/N. The weighted - * contributions to the average are cached in where - * the unweighted values were in , for use by - * p7_StochasticDsqTrace(). - * - * Returns: on success; the profile now contains - * scores and is ready for searching target sequences. - * - * Throws: on allocation error. - */ - -int -p7_profileAdjustClones(float temperature, const P7_BG *bg, const P7_PROFILE *src, P7_PROFILE *dstT, P7_PROFILE *dstDT) -{ - int x, y, n; - float dot; - const int K = src->abc->K; /* Alphabet size */ - const int Kp = src->abc->Kp; /* Alphabet size including degenerates, gaps, etc. */ - const int deg = Kp-2; /* maximally degenerate character (X/N) */ - float *bg_flog = NULL; /* For caching logf(bg->f[x]) values */ - int status; - - /* Make some cursory checks that and were once - clones of , or are otherwise ready to receive this - adjustment. Perhaps we should be more rigorous!! */ - if (src->M > dstDT->allocM) ESL_XEXCEPTION(eslEINVAL, " profile is too small to hold a copy of profile"); - if (Kp != dstDT->abc->Kp) ESL_XEXCEPTION(eslEINVAL, " and profiles have different alphabets"); - if (src->M > dstT->allocM) ESL_XEXCEPTION(eslEINVAL, " profile is too small to hold a copy of profile"); - if (Kp != dstT->abc->Kp) ESL_XEXCEPTION(eslEINVAL, " and profiles have different alphabets"); - - /* - * Update - */ - - /* Divide *all* scores from the profile by the temperature, - and record them in the profile. Did we miss any - scores?!! */ - esl_vec_FCopy (src->tsc, src->M*p7P_NTRANS, dstDT->tsc); - esl_vec_FScale(dstDT->tsc, src->M*p7P_NTRANS, 1.0f/temperature); - for (x = 0; x < Kp; x++) { - esl_vec_FCopy (src->rsc[x], (src->M+1)*p7P_NR, dstDT->rsc[x]); - esl_vec_FScale(dstDT->rsc[x], (src->M+1)*p7P_NR, 1.0f / temperature); - } - for (x = 0; x < p7P_NXSTATES; x++) { - esl_vec_FCopy (src->xsc[x], p7P_NXTRANS, dstDT->xsc[x]); - esl_vec_FScale(dstDT->xsc[x], p7P_NXTRANS, 1.0f / temperature); - } - - /* - * Update - */ - - /* Copy all scores from to . Did we miss any - scores?!! */ - esl_vec_FCopy(dstDT->tsc, src->M*p7P_NTRANS, dstT->tsc); - for (x = 0; x < Kp; x++) esl_vec_FCopy(dstDT->rsc[x], (src->M+1)*p7P_NR, dstT->rsc[x]); - for (x = 0; x < p7P_NXSTATES; x++) esl_vec_FCopy(dstDT->xsc[x], p7P_NXTRANS, dstT->xsc[x]); - - /* Compute degenerate character's score for match and insert state - emissions. Other emissions (i.e., N, C, and J loop - transitions) have score zero for all emissions, yes?!!, so - nothing needs to be done here for them. */ - - /* Cache some useful logarithms */ - ESL_ALLOC(bg_flog, sizeof(float) * K); - esl_vec_FCopy(bg->f, K, bg_flog); - esl_vec_FLog(bg_flog, K); - - /* Compute degenerate character's score for each match and - insert state using the background model as a prior - distribution for sequences. That is, this will be used for - p-value calculations. */ - n = (src->M+1)*p7P_NR; - for (y = 0; y < n; y++) /* Are loop bounds correct?!! */ { - dot = -eslINFINITY; - for (x = 0; x < K; x++) { - /* Note "+=" assignment within */ - dot = p7_FLogsum(dot, dstT->rsc[x][y] += bg_flog[x]); - } - dstT->rsc[deg][y] = dot; - } - - status = eslOK; /* Fall through to ERROR */ - - ERROR: - if (bg_flog) free (bg_flog); - return status; -} - -/* Function: p7_StochasticDsqTrace() - * Synopsis: Stochastic traceback, producing both a state trace and - * an emission trace. - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: Stochastic traceback of a grand forward sum. First - * p7_StochasticDsqTrace() samples a state trace for a - * run on a degenerate sequence . Then from the - * state trace it samples a digital sequence (of - * length ) from the profile , using the background - * model for distribution for N, C, and J loop - * emissions. For match and insert states, the probability - * that a letter is chosen is proportional to its - * contribution to the dot product computed in - * p7_profileAdjustClones(); the contribution was computed - * from both the temperature and the background prior - * model. - * - * Returns: on success. - */ - -int -p7_StochasticDsqTrace(ESL_RANDOMNESS *r, const ESL_DSQ *dsqX, int L, const P7_BG *bg, const P7_PROFILE *gmT, const P7_GMX *gxT, P7_TRACE *tr, ESL_DSQ *dsq) -{ - int x, z; - float *sc = NULL; /* array to send to esl_rnd_FChoose */ - const int K = gmT->abc->K; /* alphabet size */ - int kz; /* index of match or insert state */ - int status; - - ESL_ALLOC(sc, sizeof(float) * K); - - /* Compute a state trace from and */ - tr->N = 0; - p7_StochasticTrace(r, dsqX, L, gmT, gxT, tr); - - /* - * Use the state trace to get an emission trace - */ - - /* Because p7_StochasticTrace() is called many times with a given - (gmT, gxT) pair, it would be more efficient to somehow cache - the FLogNorm'ed vectors below!! (For instance, perhaps this - routine should compute these FlogNorms once, but return a - user-requested number of sample sequences. Alternatively, we - could store the FlogNorm'ed values in gmT inside of - p7_profileAdjustClones().) Note that what the caller does with - each we return is usually significanctly slower then our - generation of it, so this inefficiency isn't critical path. */ - for (z = 0; z < tr->N; ++z) { - if (tr->i[z] == 0) continue; /* no emission */ - switch (tr->st[z]) { - case p7T_M: - kz = tr->k[z]; /* index of match state */ - for (x = 0; x < K; x++) sc[x] = p7P_MSC(gmT, kz, x); - esl_vec_FLogNorm(sc, K); /* now sc is a prob vector */ - dsq[tr->i[z]] = esl_rnd_FChoose(r, sc, K); /* choose the emission */ - break; - case p7T_I: - kz = tr->k[z]; /* index of insert state */ - for (x = 0; x < K; x++) sc[x] = p7P_ISC(gmT, kz, x); - esl_vec_FLogNorm(sc, K); /* now sc is a prob vector */ - dsq[tr->i[z]] = esl_rnd_FChoose(r, sc, K); /* choose the emission */ - break; - case p7T_N: - case p7T_C: - case p7T_J: - /* Use background model probabilities */ - dsq[tr->i[z]] = esl_rnd_FChoose(r, bg->f, K); /* choose the emission */ - break; - } - } - - status = eslOK; /* fall through to ERROR */ - - ERROR: - if (sc != NULL) free (sc); - return status; -} - -/* Function: p7_thermo_Create() - * Synopsis: Create a P7_THERMO structure - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: For storing the relationship between temperature and - * score. Score can be forward or viterbi. - * - * Returns: The pointer to the allocated structure on success; NULL on - * failure. - * - * Throws: on allocation error. - */ - -P7_THERMO * -p7_thermo_Create(void) -{ - P7_THERMO *thermo = NULL; - int status; - - ESL_ALLOC(thermo, sizeof(P7_THERMO)); - thermo->p7m_score = p7M_NSCORETYPES; /* i.e., not set */ - thermo->numTemperatures = 0; - thermo->temperatures = NULL; - thermo->scores = NULL; - - return thermo; - - ERROR: - p7_thermo_Destroy(thermo); - return NULL; -} - -/* Function: p7_thermo_Destroy() - * Synopsis: Frees a P7_THERMO - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: Frees a P7_THERMO . - * - * Returns: nothing - */ - -void -p7_thermo_Destroy(P7_THERMO *thermo) -{ - if (thermo) { - if (thermo->temperatures) free (thermo->temperatures); - if (thermo->scores) free (thermo->scores); - free(thermo); - } -} - -/* Function: p7_thermoCalibrate() - * Synopsis: Populate a P7_THERMO structure - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: Computes the relationship between temperature and score. - * Score can be forward or viterbi. - * - * Returns: - * - * Throws: allocation errors - */ - -int -p7_thermoCalibrate(P7_THERMO *thermo, ESL_RANDOMNESS *r, int p7m_score, const P7_BG *bg, const P7_PROFILE *gm, P7_GMX *gx) -{ - P7_PROFILE *gmT = NULL; /* For computing Z(T) */ - P7_PROFILE *gmDT = NULL; /* For computing Z(D,T) */ - P7_GMX *gxT = NULL; /* For backtrace through Z(T) calculation */ - ESL_DSQ *dsqX = NULL; /* Sequence of all X/N */ - ESL_DSQ *dsq = NULL; /* Sampled sequence */ - P7_TRACE *tr = NULL; /* For backtrace through Z(T) calculation */ - float *sc = NULL; /* Array of scores at one temperature */ - float ZT = 0.0f; /* Computed Z(T) value */ - const int L = gm->L; /* Length of sequence to be scanned */ - const int deg = gm->abc->Kp-2; /* maximally degenerate character (X/N) */ - int numTemperatures = 50; /* Number of (x,y) points for the curve */ - const float firstTemperature = 0.01f; /* Lowest x=temperature value */ - const float lastTemperature = 100.f; /* Highest x=temperature value */ - int t, j, g, maxG, prMaxG, z; - float maxGap; /* Difference between adjacent y values */ - float *temperatures; /* Array of temperatures */ - float *scores; /* Array of scores by temperature */ - const int numScores = 2; /* Number of scores to try at each temperature */ - int status; - - /* Clear memory */ - if (thermo->temperatures) free (thermo->temperatures); - thermo->temperatures = NULL; - if (thermo->scores) free (thermo->scores); - thermo->scores = NULL; - - /* Allocate memory */ - if (!(numTemperatures > 0 && firstTemperature > 0.0f && lastTemperature > firstTemperature)) - ESL_XEXCEPTION(eslEINVAL, "Bad p7_thermoCalibrate parameter(s)"); - ESL_ALLOC(thermo->temperatures, sizeof(float) * numTemperatures); - ESL_ALLOC(thermo->scores, sizeof(float) * numTemperatures); - if ((gmT = p7_profile_Clone(gm)) == NULL) esl_fatal("failed to create gmT"); - if ((gmDT = p7_profile_Clone(gm)) == NULL) esl_fatal("failed to create gmDT"); - if ((gxT = p7_gmx_Create(gm->M, L)) == NULL) esl_fatal("failed to create gxT"); - if ((dsqX = malloc(sizeof(ESL_DSQ) * (L+2))) == NULL) esl_fatal("failed to create dsqX"); - if ((dsq = malloc(sizeof(ESL_DSQ) * (L+2))) == NULL) esl_fatal("failed to create dsq"); - if ((tr = p7_trace_Create()) == NULL) esl_fatal("failed to create tr"); - ESL_ALLOC(sc, sizeof(float) * numScores); - - /* Create string of all degenerate characters */ - for (z = 1; z <= L; z++) dsqX[z] = deg; /* The X/N character */ - - /* Start filling in */ - thermo->p7m_score = p7m_score; - thermo->numTemperatures = numTemperatures; - - /* local variable makes things a little faster */ - temperatures = thermo->temperatures; - scores = thermo->scores; - - prMaxG = -10; - for (t = 0; t < numTemperatures ; ++t) { - /* - * Set x coordinates, - */ - switch(t) { - case 0: - temperatures[t] = firstTemperature; - break; - case 1: - temperatures[t] = lastTemperature; - break; - default: - /* Look for biggest gap in scores. Recall that scores are - decreasing with increasing temperature. */ - maxG = 1; /* The right endpoint of the best interval found */ - maxGap = scores[maxG-1] - scores[maxG]; - for (g = 2; g < t; g++) { - if (maxGap < scores[g-1] - scores[g]) { - maxG = g; - maxGap = scores[maxG-1] - scores[maxG]; - } - } - /* There could be a big gap for a number of reasons. Two - are: (1) we don't have enough x values, and we need to - put one within this gap or (2), one of the endpoints of - this interval has a score that is atypical of the - temperature it is supposed to represent and we need to - readjust the endpoint. In this case, because of the - sorting we do, the needed new x value is on the other - side of this atypical endpoint. */ - - if (t >= 5 && maxG == prMaxG && maxG > 1) { - /* Case 2: Same left endpoint as before. Choose a new - x on the other side of that endpoint. */ - maxG--; - prMaxG = -10; /* Force case 1 next time */ - } else if (t >= 5 && maxG == prMaxG + 1 && maxG + 1 < t) { - /* Case 2: Same right endpoint as before. Choose a - new x on the other side of that endpoint. */ - maxG++; - prMaxG = -10; /* Force case 1 next time */ - } else { - /* Case 1: Put the new point at the geometric mean of - the largest interval. That is, leave - unchanged. */ - prMaxG = maxG; /* Allow either case next time. */ - } - temperatures[t] = sqrtf(temperatures[maxG-1] * temperatures[maxG]); - } - - /* - * Compute y coordinates - */ - - /* Make and profiles work for this temperature */ - p7_profileAdjustClones(temperatures[t], bg, gm, gmT, gmDT); - /* Fill up matrix with useful values for the subsequent - stochastic trace. */ - p7_GForward(dsqX, L, gmT, gxT, &ZT); - - for (j = 0; j < numScores; j++) { - /* Trace back a set of states and emissions */ - p7_StochasticDsqTrace(r, dsqX, L, bg, gmT, gxT, tr, dsq); - /* Find the score of and save it as sc[j]. */ - switch(thermo->p7m_score) { - case p7M_FORWARD: - p7_GForward(dsq, L, gm, gx, &sc[j]); - break; - case p7M_VITERBI: - p7_GViterbi(dsq, L, gm, gx, &sc[j]); - break; - } - } - - /* In an actual p-value calculation, we want most of the - sampled scores to fail to meet the threshold. (When too - many scores meet the threshold, this is an indication that - some of the contributions to the importance sampling sum - are relatively huge. That is, the ratio of background - model probability to importance sampling model probability - is high in the region just exceeding the threshold.) So we - will plug in the larger/largest of sampled scores. */ - esl_vec_FSortDecreasing(sc, numScores); - scores[t] = sc[0]; /* Largest of numScores scores */ - - /* We want the (x,y) points sorted by increasing x. But whoa, - * psychedelic, instead of simply sorting (x,y) pairs by x, we - * will separate the x and y of each pair and independently - * sort the y values, in decreasing order! Reasoning: The - * curve should be monotonically decreasing, but the inherent - * randomness in the use of a sampled , can yield scores - * that do not decrease with increasing temperature. - * Therefore we want to smooth the curve. We smooth by - * sorting the scores independently of the temperatures. This - * is equivalent to not sorting the scores, but when seeking a - * temperature for a supplied score, defining the best - * temperature as the one that has the same number of too low - * scores at lower temperatures as it has too high scores at - * higher temperatures. The sorting may cause some extra - * slope between close values of x if the true slope is small, - * but so be it. - */ - - esl_vec_FSortIncreasing(temperatures, t+1); - esl_vec_FSortDecreasing(scores, t+1); - } - - status = eslOK; - - CLEAN: - /* Free up local variables */ - p7_profile_Destroy(gmT); - p7_profile_Destroy(gmDT); - p7_gmx_Destroy(gxT); - if (dsqX) free (dsqX); - if (dsq) free (dsq); - p7_trace_Destroy(tr); - if (sc) free (sc); - - return status; - - ERROR: - /* Reset to an unused state */ - thermo->p7m_score = p7M_NSCORETYPES; - thermo->numTemperatures = 0; - if (thermo->temperatures) free (thermo->temperatures); - thermo->temperatures = NULL; - if (thermo->scores) free (thermo->scores); - thermo->scores = NULL; - goto CLEAN; -} - -/* Function: p7_thermoSuggestTemp() - * Synopsis: Suggest a temperature for use with a supplied score - * Incept: LAN, Mon Aug 18 12:31:53 EDT 2008 [Wadsworth] - * - * Purpose: Given a score threshold (forward or viterbi) for which - * we wish to compute a p-value, suggest a temperature for - * the importance sampling distribution. - * - * Returns: - */ - -int -p7_thermoSuggestTemp(const P7_THERMO *thermo, float score, float *temperature) -{ - /* Binary search on the monotonically decreasing thermo->scores - array */ - const float *temperatures = thermo->temperatures; - const float *scores = thermo->scores; - const int numTemperatures = thermo->numTemperatures; - int low; /* scores[low] <= score <= scores[hi] */ - int high; - int mid; - - low = -1; /* just out of array bounds! */ - high = numTemperatures; /* just out of array bounds! */ - while (high > low + 1) { - mid = (low + high) / 2; /* within array bounds */ - if (score >= scores[mid]) high = mid; - if (score <= scores[mid]) low = mid; - } - if (low == -1) - *temperature = temperatures[0]; - else if (high == numTemperatures) - *temperature = temperatures[numTemperatures-1]; - else if (scores[low] == scores[high]) - *temperature = 0.5f * (temperatures[low] + temperatures[high]); - else *temperature = temperatures[low] + - (scores[low] - score) / (scores[low] - scores[high]) - * (temperatures[high] - temperatures[low]); - - return eslOK; -} - -/* printMantissa(x): For printing the m.mm part of a m.mmEx.xx - floating point number, exp(x)*/ -float -printMantissa(const float x) -{ - return expf(x - floorf(x / logf(10.0f)) * logf(10.0f)); -} - -/* printExponent(x): For printing the xxx part of a m.mmExxx floating - point number, exp(x)*/ -int -printExponent(const float x) -{ - return (int) floorf(x / logf(10.0f)); -} - diff --git a/contrib/p7_thermo.h b/contrib/p7_thermo.h deleted file mode 100644 index 65dd0a92..00000000 --- a/contrib/p7_thermo.h +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef P7_THERMOH_INCLUDED -#define P7_THERMOH_INCLUDED - -#include "esl_vectorops.h" -#include "esl_getopts.h" - -/* - * Structure: P7_THERMO - * - * A curve that estimates the relationship between temperature - * (x-coordinate) and a score (forward or viterbi) it is good for. - */ - -enum p7m_scoretypes_e { - p7M_FORWARD = 0, - p7M_VITERBI = 1 -}; -#define p7M_NSCORETYPES 2 - -typedef struct p7_thermo_s { - enum p7m_scoretypes_e p7m_score; /* Whether curve describes forward or viterbi scores */ - int numTemperatures; /* array length for and */ - float *temperatures; /* x-coordinate of curve */ - float *scores; /* y-coordinates of curve */ -} P7_THERMO; - -extern void utest_thermo (ESL_GETOPTS *go, ESL_RANDOMNESS *r, const ESL_ALPHABET *abc, const P7_BG *bg, const P7_PROFILE *gm, int nseq, int L); -extern int p7_thermoEstimate (ESL_RANDOMNESS *r, const P7_THERMO *thermo, const P7_BG *bg, const P7_PROFILE *gm, P7_GMX *gx, int numSamples, float threshold, int *support, float *pv, float *pvstd, float *sn, float *snstd); -extern int p7_profileAdjustClones(float temperature, const P7_BG *bg, const P7_PROFILE *src, P7_PROFILE *dstT, P7_PROFILE *dstDT); -extern int p7_StochasticDsqTrace (ESL_RANDOMNESS *r, const ESL_DSQ *dsqX, int L, const P7_BG *bg, const P7_PROFILE *gmT, const P7_GMX *gxT, P7_TRACE *tr, ESL_DSQ *dsq); -extern P7_THERMO *p7_thermo_Create(void); -extern void p7_thermo_Destroy (P7_THERMO *thermo); -extern int p7_thermoCalibrate (P7_THERMO *thermo, ESL_RANDOMNESS *r, int p7m_score, const P7_BG *bg, const P7_PROFILE *gm, P7_GMX *gx); -extern int p7_thermoSuggestTemp (const P7_THERMO *thermo, float score, float *temperature); - -#endif /*P7_THERMOH_INCLUDED*/ - diff --git a/libdivsufsort/AUTHORS b/libdivsufsort/AUTHORS deleted file mode 100644 index 429a6a51..00000000 --- a/libdivsufsort/AUTHORS +++ /dev/null @@ -1,3 +0,0 @@ --- AUTHORS for libdivsufsort-lite - -Yuta Mori diff --git a/libdivsufsort/COPYING b/libdivsufsort/COPYING deleted file mode 100644 index 13023f75..00000000 --- a/libdivsufsort/COPYING +++ /dev/null @@ -1,27 +0,0 @@ -The libdivsufsort-lite copyright is as follows: - -Copyright (c) 2003-2008 Yuta Mori All Rights Reserved. - -Permission is hereby granted, free of charge, to any person -obtaining a copy of this software and associated documentation -files (the "Software"), to deal in the Software without -restriction, including without limitation the rights to use, -copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the -Software is furnished to do so, subject to the following -conditions: - -The above copyright notice and this permission notice shall be -included in all copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES -OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND -NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT -HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, -WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING -FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR -OTHER DEALINGS IN THE SOFTWARE. - -See also the libdivsufsort web site: - http://libdivsufsort.googlecode.com/ for more information. diff --git a/libdivsufsort/Makefile.in b/libdivsufsort/Makefile.in deleted file mode 100644 index 37eeb7c1..00000000 --- a/libdivsufsort/Makefile.in +++ /dev/null @@ -1,54 +0,0 @@ -# VPATH and shell configuration -top_srcdir = @top_srcdir@ -srcdir = @srcdir@ -VPATH = @srcdir@ -SHELL = /bin/sh - -# sources -OBJS = divsufsort.o -TARGET = libdivsufsort.a -MAKEFILE = Makefile - -# options -CC = @CC@ -CFLAGS = @CFLAGS@ -SSE_CFLAGS = @SSE_CFLAGS@ -VMX_CFLAGS = @VMX_CFLAGS@ -PTHREAD_CFLAGS = @PTHREAD_CFLAGS@ -CPPFLAGS = @CPPFLAGS@ -AR = @AR@ rc -RANLIB = @RANLIB@ -INSTALL = @INSTALL@ - -# beautification magic stolen from git (added within hmmer source) -ifndef V - QUIET_CC = @echo ' ' CC $@; - QUIET_AR = @echo ' ' AR $@; -endif - - -# targets -all: $(TARGET) -$(TARGET): $(OBJS) - -.c.o: - ${QUIET_CC}${CC} -I. ${CFLAGS} ${SSE_CFLAGS} ${VMX_CFLAGS} ${PTHREAD_CFLAGS} ${CPPFLAGS} -o $@ -c $< - - -libdivsufsort.a: $(OBJS) - ${QUIET_AR}${AR} libdivsufsort.a $(OBJS) - @${RANLIB} libdivsufsort.a - @chmod 644 libdivsufsort.a - - -clean: - $(RM) $(TARGET) $(OBJS) -ifndef V - @echo ' ' CLEAN libdivsufsort -endif - - -distclean: - $(RM) $(TARGET) $(OBJS) $(MAKEFILE) divsufsort.h - - diff --git a/libdivsufsort/README b/libdivsufsort/README deleted file mode 100644 index 63c6655d..00000000 --- a/libdivsufsort/README +++ /dev/null @@ -1,185 +0,0 @@ -libdivsufsort - A lightweight suffix-sorting library. ------------------------------------------------------ - -Introduction: -------------- - -The libdivsufsort project provides a fast, lightweight, and robust -C API library to construct the suffix array and the Burrows-Wheeler -transformed string for any input string of a constant-size alphabet. - -The suffix-sorting algorithm runs in O(n log n) worst-case time -using only 5n+O(1) bytes of memory space, where n is the length of -the input string. - -The latest version of libdivsufsort is available at: - http://libdivsufsort.googlecode.com/ - - -License: --------- - -libdivsufsort is released under the MIT/X11 license. See the file -COPYING for more details. - - -APIs: ------ - - * Data types - typedef int32_t saint_t; - typedef int32_t saidx_t; - typedef uint8_t sauchar_t; - - * Constructs the suffix array of a given string. - * @param T[0..n-1] The input string. - * @param SA[0..n-1] The output array or suffixes. - * @param n The length of the given string. - * @return 0 if no error occurred, -1 or -2 otherwise. - saint_t - divsufsort(const sauchar_t *T, saidx_t *SA, saidx_t n); - - * Constructs the burrows-wheeler transformed string of a given string. - * @param T[0..n-1] The input string. - * @param U[0..n-1] The output string. (can be T) - * @param A[0..n-1] The temporary array. (can be NULL) - * @param n The length of the given string. - * @return The primary index if no error occurred, -1 or -2 otherwise. - saidx_t - divbwt(const sauchar_t *T, sauchar_t *U, saidx_t *A, saidx_t n); - - -Benchmark: ------------------- - -= Specifications = -Processor: 2.66 GHz Intel Core 2 Duo E6750 -L1 Cache: (32 Kb + 32 Kb) x 2 -L2 Cache: 4 Mb -RAM: 2 Gb main memory -Operating system: Windows XP Home SP 3 (with Cygwin) -Compiler: GCC version 4.3.1 - -= Programs = -Archon4r0 kvark's sorting algorithm http://forum.compression.ru/viewtopic.php?t=352 -BPR Bucket-Pointer Refinement algorithm http://bibiserv.techfak.uni-bielefeld.de/bpr/ -DC Difference-Cover algorithm (v = 32) http://www.cs.helsinki.fi/juha.karkkainen/publications/cpm03.tar.gz -DS Deep-Shallow sorting algorithm http://www.mfn.unipmn.it/~manzini/lightweight/ -divsufsort1 libdivsufsort version 1.2.3 http://libdivsufsort.googlecode.com/ -divsufsort2 libdivsufsort version 2.0.0 http://libdivsufsort.googlecode.com/ -KA Ko-Aluru algorithm http://ko.pang.cn.googlepages.com/software2 -KS Kärkkäinen-Sanders algorithm http://www.mpi-inf.mpg.de/~sanders/programs/suffix/ -MSufSort3 MSufSort version 3.1.1 beta http://www.michael-maniscalco.com/msufsort.htm -qsufsort Larsson-Sadakane algorithm http://www.larsson.dogma.net/research.html -sais Induced Sorting algorithm http://yuta.256.googlepages.com/sais - -All programs were compiled with gcc/g++ using '-O3 -fomit-frame-pointer -DNDEBUG' -optimization options. The times are the average of five runs, in seconds, and were -measured using the standard Unix/Cygwin 'time' command. (user + system) The spaces -were measured using the 'memusage' command. - -= Testfiles = -Manzini's Large Corpus http://www.mfn.unipmn.it/~manzini/lightweight/corpus/ -The Gauntlet http://www.michael-maniscalco.com/testset/gauntlet/ - -= Running times = - -== Manzini's Corpus == -Files Size Archon4r0 BPR DC DS divsufsort1 divsufsort2 KA KS MSufSort3 qsufsort sais -chr22.dna 34553758 6.030 6.196 22.694 7.514 5.404 5.362 16.980 50.006 7.132 10.642 10.796 -etext99 105277340 22.160 32.582 79.872 34.264 18.758 18.064 73.236 202.684 24.106 56.612 38.748 -gcc-3.0.tar 86630400 13.856 20.692 61.690 35.822 10.382 10.084 40.908 135.174 14.952 40.766 20.990 -howto 39422105 5.806 8.326 25.432 8.288 5.472 5.320 20.694 64.834 5.672 16.366 11.388 -jdk13c 69728899 18.106 22.252 61.234 32.182 9.260 9.010 34.172 101.096 11.314 39.792 16.396 -linux-2.4.5.tar 116254720 18.174 26.226 82.830 25.912 14.672 14.290 58.586 194.412 19.890 54.054 29.614 -rctail96 114711151 32.490 55.826 119.026 62.502 18.500 17.914 70.072 190.562 21.060 70.456 33.248 -rfc 116421901 20.736 35.404 91.284 29.666 16.116 15.658 64.390 196.500 17.936 61.436 32.224 -sprot34.dat 109617186 22.832 36.720 93.122 32.096 17.894 17.404 68.084 187.594 23.352 56.946 34.092 -w3c2 104201579 27.264 29.384 89.352 54.682 13.866 13.486 52.660 162.582 17.090 77.804 25.498 -totals 896819039 187.454 273.608 726.536 322.928 130.324 126.592 499.782 1485.444 162.504 484.874 252.994 - -== The Gauntlet == -Files Size Archon4r0 BPR DC DS divsufsort1 divsufsort2 KA KS MSufSort3 qsufsort sais -abac 200000 0.044 0.064 0.104 27.914 0.042 0.036 0.058 0.048 0.050 0.062 0.044 -abba 10500600 3.270 5.124 10.766 30.702 1.714 1.602 2.570 7.952 3.514 15.272 1.460 -book1x20 15375420 4.392 3.530 13.872 97.468 2.312 2.154 7.442 15.756 3.542 22.376 3.912 -fib_s14930352 14930352 12.728 10.830 18.524 179.040 3.638 3.588 3.544 10.232 6.700 18.224 2.542 -fss10 12078908 11.390 8.974 15.130 85.328 2.828 2.824 3.344 8.646 4.618 14.754 2.076 -fss9 2851443 1.002 1.210 1.644 5.256 0.410 0.416 0.618 1.290 0.554 2.836 0.336 -houston 3840000 0.344 0.708 2.226 118.960 0.118 0.128 0.520 0.744 0.242 1.230 0.238 -paper5x80 981924 0.110 0.154 0.454 0.806 0.092 0.090 0.210 0.256 0.144 0.448 0.110 -test1 2097152 0.332 2.132 1.108 8.680 0.268 0.280 0.376 1.066 1.302 2.762 0.202 -test2 2097152 0.710 0.616 1.110 8.682 0.180 0.176 0.374 1.076 3.354 2.768 0.206 -test3 2097152 0.488 213.154 1.164 1.772 0.220 0.226 0.388 1.082 0.922 3.246 0.212 -totals 67050103 34.810 246.496 66.102 564.608 11.822 11.520 19.444 48.148 24.942 83.978 11.338 - -= Space (in MiBytes) = - -== Manzini's Corpus == -Files Size Archon4r0 BPR DC DS divsufsort1 divsufsort2 KA KS MSufSort3 qsufsort sais -chr22.dna 34553758 174.66 296.88 193.60 165.18 165.02 165.02 289.97 428.39 199.72 263.62 164.77 -etext99 105277340 531.13 915.48 589.85 503.23 502.25 502.25 907.34 1305.20 604.45 803.20 502.00 -gcc-3.0.tar 86630400 437.14 756.43 485.38 415.87 413.34 413.34 709.50 1074.01 497.79 660.94 413.09 -howto 39422105 199.20 367.53 220.88 188.45 188.23 188.23 331.54 488.75 227.67 300.77 187.98 -jdk13c 69728899 351.96 603.99 390.68 333.40 332.74 332.74 609.71 864.48 401.04 531.99 332.49 -linux-2.4.5.tar 116254720 586.46 1061.83 651.36 555.76 554.60 554.60 977.81 1441.30 667.39 886.95 554.35 -rctail96 114711151 578.68 987.64 642.71 548.32 547.24 547.24 1004.98 1422.16 658.43 875.18 546.99 -rfc 116421901 587.30 1005.85 652.29 556.53 555.39 555.39 956.52 1443.37 668.26 888.23 555.14 -sprot34.dat 109617186 553.01 941.95 614.17 524.03 522.95 522.95 930.06 1359.01 629.26 836.31 522.70 -w3c2 104201579 525.71 958.37 583.82 498.09 497.12 497.12 912.00 1291.87 598.82 795.00 496.87 -totals 896819039 4525.25 7895.95 5024.74 4288.86 4278.88 4278.88 7629.43 11118.54 5152.83 6842.19 4276.38 -mean - 5.29 9.23 5.88 5.01 5.00 5.00 8.92 13.00 6.02 8.00 5.00 - -== The Gauntlet == -Files Size Archon4r0 BPR DC DS divsufsort1 divsufsort2 KA KS MSufSort3 qsufsort sais -abac 200000 1.51 1.73 1.12 0.98 1.21 1.20 1.75 2.48 3.15 1.53 0.95 -abba 10500600 53.43 90.19 58.83 50.21 50.32 50.32 86.20 130.18 62.09 80.11 50.07 -book1x20 15375420 78.00 134.00 86.15 73.52 73.57 73.57 132.42 190.62 89.99 117.31 73.32 -fib_s14930352 14930352 75.75 128.15 83.65 71.71 71.44 71.44 117.16 185.10 87.43 113.91 71.19 -fss10 12078908 61.38 103.68 67.68 58.05 57.85 57.85 107.05 149.75 71.12 92.16 57.60 -fss9 2851443 14.87 24.48 15.98 13.71 13.85 13.85 25.27 35.35 18.32 21.76 13.60 -houston 3840000 19.85 36.96 21.52 18.46 18.56 18.56 28.79 47.58 23.98 29.30 18.31 -paper5x80 981924 5.45 11.40 5.50 4.72 4.93 4.93 8.59 12.17 7.63 7.49 4.68 -test1 2097152 11.07 82.00 11.75 10.10 10.25 10.25 18.34 25.99 14.01 16.00 10.00 -test2 2097152 11.07 82.00 11.75 10.10 10.25 10.25 18.34 25.99 14.01 16.00 10.00 -test3 2097152 11.07 82.00 11.75 10.05 10.25 10.25 18.34 26.00 14.63 16.00 10.12 -totals 67050103 343.45 776.59 375.68 321.61 322.48 322.47 562.25 831.21 406.36 511.57 319.84 -mean - 5.37 12.14 5.88 5.03 5.04 5.04 8.79 13.00 6.35 8.00 5.00 - - -Algorithm: ----------- - -libdivsufsort uses the following algorithms for suffix sorting. - - The improved version of Itho-Tanaka two-stage sorting algorithm. [2][6] - - A substring sorting/encoding technique. [1][3] - - Maniscalco's tandem repeat sorting algorithm. [5] - - Larsson-Sadakane sorting algorithm. [4] - - -References: ------------ - - 1. Stefan Burkhardt and Juha K"arkk"ainen. Fast lightweight suffix - array construction and checking. Proceedings of the 14th Annual - Symposium on Combinatorial Pattern Matching, LNCS 2676, - Springer, pp. 55-69, 2003. - - 2. Hideo Itoh and Hozumi Tanaka, An Efficient Method for in Memory - Construction of Suffix Arrays, Proceedings of the IEEE String - Processing and Information Retrieval Symposium, pp. 81-88, 1999. - - 3. Pang Ko and Srinivas Aluru, Space-efficient linear time - construction of suffix arrays, Proceedings of the 14th Annual - Symposium on Combinatorial Pattern Matching, pp. 200-210, 2003. - - 4. Jesper Larsson and Kunihiko Sadakane, Faster suffix sorting. - Technical report LU-CS-TR:99-214, Department of Computer - Science, Lund University, Sweden, 1999. - - 5. Michael Maniscalco, MSufSort. - http://www.michael-maniscalco.com/msufsort.htm - - 6. Yuta Mori, Short description of improved two-stage suffix sorting - algorithm, 2005. - http://homepage3.nifty.com/wpage/software/itssort.txt diff --git a/libdivsufsort/VERSION b/libdivsufsort/VERSION deleted file mode 100644 index 227cea21..00000000 --- a/libdivsufsort/VERSION +++ /dev/null @@ -1 +0,0 @@ -2.0.0 diff --git a/libdivsufsort/divsufsort.c b/libdivsufsort/divsufsort.c deleted file mode 100644 index 4499b665..00000000 --- a/libdivsufsort/divsufsort.c +++ /dev/null @@ -1,1794 +0,0 @@ -/* - * divsufsort.c for libdivsufsort-lite - * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved. - * - * Permission is hereby granted, free of charge, to any person - * obtaining a copy of this software and associated documentation - * files (the "Software"), to deal in the Software without - * restriction, including without limitation the rights to use, - * copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following - * conditions: - * - * The above copyright notice and this permission notice shall be - * included in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES - * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT - * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, - * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING - * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR - * OTHER DEALINGS IN THE SOFTWARE. - */ - -#include -#include -#include -#ifdef _OPENMP -# include -#endif -#include "divsufsort.h" - - -/*- Constants -*/ -#define INLINE __inline -#if defined(ALPHABET_SIZE) && (ALPHABET_SIZE < 1) -# undef ALPHABET_SIZE -#endif -#if !defined(ALPHABET_SIZE) -# define ALPHABET_SIZE (256) -#endif -#define BUCKET_A_SIZE (ALPHABET_SIZE) -#define BUCKET_B_SIZE (ALPHABET_SIZE * ALPHABET_SIZE) -#if defined(SS_INSERTIONSORT_THRESHOLD) -# if SS_INSERTIONSORT_THRESHOLD < 1 -# undef SS_INSERTIONSORT_THRESHOLD -# define SS_INSERTIONSORT_THRESHOLD (1) -# endif -#else -# define SS_INSERTIONSORT_THRESHOLD (8) -#endif -#if defined(SS_BLOCKSIZE) -# if SS_BLOCKSIZE < 0 -# undef SS_BLOCKSIZE -# define SS_BLOCKSIZE (0) -# elif 32768 <= SS_BLOCKSIZE -# undef SS_BLOCKSIZE -# define SS_BLOCKSIZE (32767) -# endif -#else -# define SS_BLOCKSIZE (1024) -#endif -/* minstacksize = log(SS_BLOCKSIZE) / log(3) * 2 */ -#if SS_BLOCKSIZE == 0 -# define SS_MISORT_STACKSIZE (96) -#elif SS_BLOCKSIZE <= 4096 -# define SS_MISORT_STACKSIZE (16) -#else -# define SS_MISORT_STACKSIZE (24) -#endif -#define SS_SMERGE_STACKSIZE (32) -#define TR_INSERTIONSORT_THRESHOLD (8) -#define TR_STACKSIZE (64) - - -/*- Macros -*/ -#ifndef SWAP -# define SWAP(_a, _b) do { t = (_a); (_a) = (_b); (_b) = t; } while(0) -#endif /* SWAP */ -#ifndef MIN -# define MIN(_a, _b) (((_a) < (_b)) ? (_a) : (_b)) -#endif /* MIN */ -#ifndef MAX -# define MAX(_a, _b) (((_a) > (_b)) ? (_a) : (_b)) -#endif /* MAX */ -#define STACK_PUSH(_a, _b, _c, _d)\ - do {\ - assert(ssize < STACK_SIZE);\ - stack[ssize].a = (_a), stack[ssize].b = (_b),\ - stack[ssize].c = (_c), stack[ssize++].d = (_d);\ - } while(0) -#define STACK_PUSH5(_a, _b, _c, _d, _e)\ - do {\ - assert(ssize < STACK_SIZE);\ - stack[ssize].a = (_a), stack[ssize].b = (_b),\ - stack[ssize].c = (_c), stack[ssize].d = (_d), stack[ssize++].e = (_e);\ - } while(0) -#define STACK_POP(_a, _b, _c, _d)\ - do {\ - assert(0 <= ssize);\ - if(ssize == 0) { return; }\ - (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\ - (_c) = stack[ssize].c, (_d) = stack[ssize].d;\ - } while(0) -#define STACK_POP5(_a, _b, _c, _d, _e)\ - do {\ - assert(0 <= ssize);\ - if(ssize == 0) { return; }\ - (_a) = stack[--ssize].a, (_b) = stack[ssize].b,\ - (_c) = stack[ssize].c, (_d) = stack[ssize].d, (_e) = stack[ssize].e;\ - } while(0) -#define BUCKET_A(_c0) bucket_A[(_c0)] -#if ALPHABET_SIZE == 256 -#define BUCKET_B(_c0, _c1) (bucket_B[((_c1) << 8) | (_c0)]) -#define BUCKET_BSTAR(_c0, _c1) (bucket_B[((_c0) << 8) | (_c1)]) -#else -#define BUCKET_B(_c0, _c1) (bucket_B[(_c1) * ALPHABET_SIZE + (_c0)]) -#define BUCKET_BSTAR(_c0, _c1) (bucket_B[(_c0) * ALPHABET_SIZE + (_c1)]) -#endif - - -/*- Private Functions -*/ - -static const int lg_table[256]= { - -1,0,1,1,2,2,2,2,3,3,3,3,3,3,3,3,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4, - 5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5,5, - 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, - 6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6,6, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7, - 7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7,7 -}; - -#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) - -static INLINE -int -ss_ilg(int n) { -#if SS_BLOCKSIZE == 0 - return (n & 0xffff0000) ? - ((n & 0xff000000) ? - 24 + lg_table[(n >> 24) & 0xff] : - 16 + lg_table[(n >> 16) & 0xff]) : - ((n & 0x0000ff00) ? - 8 + lg_table[(n >> 8) & 0xff] : - 0 + lg_table[(n >> 0) & 0xff]); -#elif SS_BLOCKSIZE < 256 - return lg_table[n]; -#else - return (n & 0xff00) ? - 8 + lg_table[(n >> 8) & 0xff] : - 0 + lg_table[(n >> 0) & 0xff]; -#endif -} - -#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */ - -#if SS_BLOCKSIZE != 0 - -static const int sqq_table[256] = { - 0, 16, 22, 27, 32, 35, 39, 42, 45, 48, 50, 53, 55, 57, 59, 61, - 64, 65, 67, 69, 71, 73, 75, 76, 78, 80, 81, 83, 84, 86, 87, 89, - 90, 91, 93, 94, 96, 97, 98, 99, 101, 102, 103, 104, 106, 107, 108, 109, -110, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, -128, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, -143, 144, 144, 145, 146, 147, 148, 149, 150, 150, 151, 152, 153, 154, 155, 155, -156, 157, 158, 159, 160, 160, 161, 162, 163, 163, 164, 165, 166, 167, 167, 168, -169, 170, 170, 171, 172, 173, 173, 174, 175, 176, 176, 177, 178, 178, 179, 180, -181, 181, 182, 183, 183, 184, 185, 185, 186, 187, 187, 188, 189, 189, 190, 191, -192, 192, 193, 193, 194, 195, 195, 196, 197, 197, 198, 199, 199, 200, 201, 201, -202, 203, 203, 204, 204, 205, 206, 206, 207, 208, 208, 209, 209, 210, 211, 211, -212, 212, 213, 214, 214, 215, 215, 216, 217, 217, 218, 218, 219, 219, 220, 221, -221, 222, 222, 223, 224, 224, 225, 225, 226, 226, 227, 227, 228, 229, 229, 230, -230, 231, 231, 232, 232, 233, 234, 234, 235, 235, 236, 236, 237, 237, 238, 238, -239, 240, 240, 241, 241, 242, 242, 243, 243, 244, 244, 245, 245, 246, 246, 247, -247, 248, 248, 249, 249, 250, 250, 251, 251, 252, 252, 253, 253, 254, 254, 255 -}; - -static INLINE -int -ss_isqrt(int x) { - int y, e; - - if(x >= (SS_BLOCKSIZE * SS_BLOCKSIZE)) { return SS_BLOCKSIZE; } - e = (x & 0xffff0000) ? - ((x & 0xff000000) ? - 24 + lg_table[(x >> 24) & 0xff] : - 16 + lg_table[(x >> 16) & 0xff]) : - ((x & 0x0000ff00) ? - 8 + lg_table[(x >> 8) & 0xff] : - 0 + lg_table[(x >> 0) & 0xff]); - - if(e >= 16) { - y = sqq_table[x >> ((e - 6) - (e & 1))] << ((e >> 1) - 7); - if(e >= 24) { y = (y + 1 + x / y) >> 1; } - y = (y + 1 + x / y) >> 1; - } else if(e >= 8) { - y = (sqq_table[x >> ((e - 6) - (e & 1))] >> (7 - (e >> 1))) + 1; - } else { - return sqq_table[x] >> 4; - } - - return (x < (y * y)) ? y - 1 : y; -} - -#endif /* SS_BLOCKSIZE != 0 */ - - -/*---------------------------------------------------------------------------*/ - -/* Compares two suffixes. */ -static INLINE -int -ss_compare(const unsigned char *T, - const int *p1, const int *p2, - int depth) { - const unsigned char *U1, *U2, *U1n, *U2n; - - for(U1 = T + depth + *p1, - U2 = T + depth + *p2, - U1n = T + *(p1 + 1) + 2, - U2n = T + *(p2 + 1) + 2; - (U1 < U1n) && (U2 < U2n) && (*U1 == *U2); - ++U1, ++U2) { - } - - return U1 < U1n ? - (U2 < U2n ? *U1 - *U2 : 1) : - (U2 < U2n ? -1 : 0); -} - - -/*---------------------------------------------------------------------------*/ - -#if (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) - -/* Insertionsort for small size groups */ -static -void -ss_insertionsort(const unsigned char *T, const int *PA, - int *first, int *last, int depth) { - int *i, *j; - int t; - int r; - - for(i = last - 2; first <= i; --i) { - for(t = *i, j = i + 1; 0 < (r = ss_compare(T, PA + t, PA + *j, depth));) { - do { *(j - 1) = *j; } while((++j < last) && (*j < 0)); - if(last <= j) { break; } - } - if(r == 0) { *j = ~*j; } - *(j - 1) = t; - } -} - -#endif /* (SS_BLOCKSIZE != 1) && (SS_INSERTIONSORT_THRESHOLD != 1) */ - - -/*---------------------------------------------------------------------------*/ - -#if (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) - -static INLINE -void -ss_fixdown(const unsigned char *Td, const int *PA, - int *SA, int i, int size) { - int j, k; - int v; - int c, d, e; - - for(v = SA[i], c = Td[PA[v]]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) { - d = Td[PA[SA[k = j++]]]; - if(d < (e = Td[PA[SA[j]]])) { k = j; d = e; } - if(d <= c) { break; } - } - SA[i] = v; -} - -/* Simple top-down heapsort. */ -static -void -ss_heapsort(const unsigned char *Td, const int *PA, int *SA, int size) { - int i, m; - int t; - - m = size; - if((size % 2) == 0) { - m--; - if(Td[PA[SA[m / 2]]] < Td[PA[SA[m]]]) { SWAP(SA[m], SA[m / 2]); } - } - - for(i = m / 2 - 1; 0 <= i; --i) { ss_fixdown(Td, PA, SA, i, m); } - if((size % 2) == 0) { SWAP(SA[0], SA[m]); ss_fixdown(Td, PA, SA, 0, m); } - for(i = m - 1; 0 < i; --i) { - t = SA[0], SA[0] = SA[i]; - ss_fixdown(Td, PA, SA, 0, i); - SA[i] = t; - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Returns the median of three elements. */ -static INLINE -int * -ss_median3(const unsigned char *Td, const int *PA, - int *v1, int *v2, int *v3) { - int *t; - if(Td[PA[*v1]] > Td[PA[*v2]]) { SWAP(v1, v2); } - if(Td[PA[*v2]] > Td[PA[*v3]]) { - if(Td[PA[*v1]] > Td[PA[*v3]]) { return v1; } - else { return v3; } - } - return v2; -} - -/* Returns the median of five elements. */ -static INLINE -int * -ss_median5(const unsigned char *Td, const int *PA, - int *v1, int *v2, int *v3, int *v4, int *v5) { - int *t; - if(Td[PA[*v2]] > Td[PA[*v3]]) { SWAP(v2, v3); } - if(Td[PA[*v4]] > Td[PA[*v5]]) { SWAP(v4, v5); } - if(Td[PA[*v2]] > Td[PA[*v4]]) { SWAP(v2, v4); SWAP(v3, v5); } - if(Td[PA[*v1]] > Td[PA[*v3]]) { SWAP(v1, v3); } - if(Td[PA[*v1]] > Td[PA[*v4]]) { SWAP(v1, v4); SWAP(v3, v5); } - if(Td[PA[*v3]] > Td[PA[*v4]]) { return v4; } - return v3; -} - -/* Returns the pivot element. */ -static INLINE -int * -ss_pivot(const unsigned char *Td, const int *PA, int *first, int *last) { - int *middle; - int t; - - t = last - first; - middle = first + t / 2; - - if(t <= 512) { - if(t <= 32) { - return ss_median3(Td, PA, first, middle, last - 1); - } else { - t >>= 2; - return ss_median5(Td, PA, first, first + t, middle, last - 1 - t, last - 1); - } - } - t >>= 3; - first = ss_median3(Td, PA, first, first + t, first + (t << 1)); - middle = ss_median3(Td, PA, middle - t, middle, middle + t); - last = ss_median3(Td, PA, last - 1 - (t << 1), last - 1 - t, last - 1); - return ss_median3(Td, PA, first, middle, last); -} - - -/*---------------------------------------------------------------------------*/ - -/* Binary partition for substrings. */ -static INLINE -int * -ss_partition(const int *PA, - int *first, int *last, int depth) { - int *a, *b; - int t; - for(a = first - 1, b = last;;) { - for(; (++a < b) && ((PA[*a] + depth) >= (PA[*a + 1] + 1));) { *a = ~*a; } - for(; (a < --b) && ((PA[*b] + depth) < (PA[*b + 1] + 1));) { } - if(b <= a) { break; } - t = ~*b; - *b = *a; - *a = t; - } - if(first < a) { *first = ~*first; } - return a; -} - -/* Multikey introsort for medium size groups. */ -static -void -ss_mintrosort(const unsigned char *T, const int *PA, - int *first, int *last, - int depth) { -#define STACK_SIZE SS_MISORT_STACKSIZE - struct { int *a, *b, c; int d; } stack[STACK_SIZE]; - const unsigned char *Td; - int *a, *b, *c, *d, *e, *f; - int s, t; - int ssize; - int limit; - int v, x = 0; - - for(ssize = 0, limit = ss_ilg(last - first);;) { - - if((last - first) <= SS_INSERTIONSORT_THRESHOLD) { -#if 1 < SS_INSERTIONSORT_THRESHOLD - if(1 < (last - first)) { ss_insertionsort(T, PA, first, last, depth); } -#endif - STACK_POP(first, last, depth, limit); - continue; - } - - Td = T + depth; - if(limit-- == 0) { ss_heapsort(Td, PA, first, last - first); } - if(limit < 0) { - for(a = first + 1, v = Td[PA[*first]]; a < last; ++a) { - if((x = Td[PA[*a]]) != v) { - if(1 < (a - first)) { break; } - v = x; - first = a; - } - } - if(Td[PA[*first] - 1] < v) { - first = ss_partition(PA, first, a, depth); - } - if((a - first) <= (last - a)) { - if(1 < (a - first)) { - STACK_PUSH(a, last, depth, -1); - last = a, depth += 1, limit = ss_ilg(a - first); - } else { - first = a, limit = -1; - } - } else { - if(1 < (last - a)) { - STACK_PUSH(first, a, depth + 1, ss_ilg(a - first)); - first = a, limit = -1; - } else { - last = a, depth += 1, limit = ss_ilg(a - first); - } - } - continue; - } - - /* choose pivot */ - a = ss_pivot(Td, PA, first, last); - v = Td[PA[*a]]; - SWAP(*first, *a); - - /* partition */ - for(b = first; (++b < last) && ((x = Td[PA[*b]]) == v);) { } - if(((a = b) < last) && (x < v)) { - for(; (++b < last) && ((x = Td[PA[*b]]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - } - for(c = last; (b < --c) && ((x = Td[PA[*c]]) == v);) { } - if((b < (d = c)) && (x > v)) { - for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - for(; b < c;) { - SWAP(*b, *c); - for(; (++b < c) && ((x = Td[PA[*b]]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - for(; (b < --c) && ((x = Td[PA[*c]]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - - if(a <= d) { - c = b - 1; - - if((s = a - first) > (t = b - a)) { s = t; } - for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - if((s = d - c) > (t = last - d - 1)) { s = t; } - for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - - a = first + (b - a), c = last - (d - c); - b = (v <= Td[PA[*a] - 1]) ? a : ss_partition(PA, a, c, depth); - - if((a - first) <= (last - c)) { - if((last - c) <= (c - b)) { - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - STACK_PUSH(c, last, depth, limit); - last = a; - } else if((a - first) <= (c - b)) { - STACK_PUSH(c, last, depth, limit); - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - last = a; - } else { - STACK_PUSH(c, last, depth, limit); - STACK_PUSH(first, a, depth, limit); - first = b, last = c, depth += 1, limit = ss_ilg(c - b); - } - } else { - if((a - first) <= (c - b)) { - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - STACK_PUSH(first, a, depth, limit); - first = c; - } else if((last - c) <= (c - b)) { - STACK_PUSH(first, a, depth, limit); - STACK_PUSH(b, c, depth + 1, ss_ilg(c - b)); - first = c; - } else { - STACK_PUSH(first, a, depth, limit); - STACK_PUSH(c, last, depth, limit); - first = b, last = c, depth += 1, limit = ss_ilg(c - b); - } - } - } else { - limit += 1; - if(Td[PA[*first] - 1] < v) { - first = ss_partition(PA, first, last, depth); - limit = ss_ilg(last - first); - } - depth += 1; - } - } -#undef STACK_SIZE -} - -#endif /* (SS_BLOCKSIZE == 0) || (SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE) */ - - -/*---------------------------------------------------------------------------*/ - -#if SS_BLOCKSIZE != 0 - -static INLINE -void -ss_blockswap(int *a, int *b, int n) { - int t; - for(; 0 < n; --n, ++a, ++b) { - t = *a, *a = *b, *b = t; - } -} - -static INLINE -void -ss_rotate(int *first, int *middle, int *last) { - int *a, *b, t; - int l, r; - l = middle - first, r = last - middle; - for(; (0 < l) && (0 < r);) { - if(l == r) { ss_blockswap(first, middle, l); break; } - if(l < r) { - a = last - 1, b = middle - 1; - t = *a; - do { - *a-- = *b, *b-- = *a; - if(b < first) { - *a = t; - last = a; - if((r -= l + 1) <= l) { break; } - a -= 1, b = middle - 1; - t = *a; - } - } while(1); - } else { - a = first, b = middle; - t = *a; - do { - *a++ = *b, *b++ = *a; - if(last <= b) { - *a = t; - first = a + 1; - if((l -= r + 1) <= r) { break; } - a += 1, b = middle; - t = *a; - } - } while(1); - } - } -} - - -/*---------------------------------------------------------------------------*/ - -static -void -ss_inplacemerge(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int depth) { - const int *p; - int *a, *b; - int len, half; - int q, r; - int x; - - for(;;) { - if(*(last - 1) < 0) { x = 1; p = PA + ~*(last - 1); } - else { x = 0; p = PA + *(last - 1); } - for(a = first, len = middle - first, half = len >> 1, r = -1; - 0 < len; - len = half, half >>= 1) { - b = a + half; - q = ss_compare(T, PA + ((0 <= *b) ? *b : ~*b), p, depth); - if(q < 0) { - a = b + 1; - half -= (len & 1) ^ 1; - } else { - r = q; - } - } - if(a < middle) { - if(r == 0) { *a = ~*a; } - ss_rotate(a, middle, last); - last -= middle - a; - middle = a; - if(first == middle) { break; } - } - --last; - if(x != 0) { while(*--last < 0) { } } - if(middle == last) { break; } - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Merge-forward with internal buffer. */ -static -void -ss_mergeforward(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int *buf, int depth) { - int *a, *b, *c, *bufend; - int t; - int r; - - bufend = buf + (middle - first) - 1; - ss_blockswap(buf, first, middle - first); - - for(t = *(a = first), b = buf, c = middle;;) { - r = ss_compare(T, PA + *b, PA + *c, depth); - if(r < 0) { - do { - *a++ = *b; - if(bufend <= b) { *bufend = t; return; } - *b++ = *a; - } while(*b < 0); - } else if(r > 0) { - do { - *a++ = *c, *c++ = *a; - if(last <= c) { - while(b < bufend) { *a++ = *b, *b++ = *a; } - *a = *b, *b = t; - return; - } - } while(*c < 0); - } else { - *c = ~*c; - do { - *a++ = *b; - if(bufend <= b) { *bufend = t; return; } - *b++ = *a; - } while(*b < 0); - - do { - *a++ = *c, *c++ = *a; - if(last <= c) { - while(b < bufend) { *a++ = *b, *b++ = *a; } - *a = *b, *b = t; - return; - } - } while(*c < 0); - } - } -} - -/* Merge-backward with internal buffer. */ -static -void -ss_mergebackward(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int *buf, int depth) { - const int *p1, *p2; - int *a, *b, *c, *bufend; - int t; - int r; - int x; - - bufend = buf + (last - middle) - 1; - ss_blockswap(buf, middle, last - middle); - - x = 0; - if(*bufend < 0) { p1 = PA + ~*bufend; x |= 1; } - else { p1 = PA + *bufend; } - if(*(middle - 1) < 0) { p2 = PA + ~*(middle - 1); x |= 2; } - else { p2 = PA + *(middle - 1); } - for(t = *(a = last - 1), b = bufend, c = middle - 1;;) { - r = ss_compare(T, p1, p2, depth); - if(0 < r) { - if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; } - *a-- = *b; - if(b <= buf) { *buf = t; break; } - *b-- = *a; - if(*b < 0) { p1 = PA + ~*b; x |= 1; } - else { p1 = PA + *b; } - } else if(r < 0) { - if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; } - *a-- = *c, *c-- = *a; - if(c < first) { - while(buf < b) { *a-- = *b, *b-- = *a; } - *a = *b, *b = t; - break; - } - if(*c < 0) { p2 = PA + ~*c; x |= 2; } - else { p2 = PA + *c; } - } else { - if(x & 1) { do { *a-- = *b, *b-- = *a; } while(*b < 0); x ^= 1; } - *a-- = ~*b; - if(b <= buf) { *buf = t; break; } - *b-- = *a; - if(x & 2) { do { *a-- = *c, *c-- = *a; } while(*c < 0); x ^= 2; } - *a-- = *c, *c-- = *a; - if(c < first) { - while(buf < b) { *a-- = *b, *b-- = *a; } - *a = *b, *b = t; - break; - } - if(*b < 0) { p1 = PA + ~*b; x |= 1; } - else { p1 = PA + *b; } - if(*c < 0) { p2 = PA + ~*c; x |= 2; } - else { p2 = PA + *c; } - } - } -} - -/* D&C based merge. */ -static -void -ss_swapmerge(const unsigned char *T, const int *PA, - int *first, int *middle, int *last, - int *buf, int bufsize, int depth) { -#define STACK_SIZE SS_SMERGE_STACKSIZE -#define GETIDX(a) ((0 <= (a)) ? (a) : (~(a))) -#define MERGE_CHECK(a, b, c)\ - do {\ - if(((c) & 1) ||\ - (((c) & 2) && (ss_compare(T, PA + GETIDX(*((a) - 1)), PA + *(a), depth) == 0))) {\ - *(a) = ~*(a);\ - }\ - if(((c) & 4) && ((ss_compare(T, PA + GETIDX(*((b) - 1)), PA + *(b), depth) == 0))) {\ - *(b) = ~*(b);\ - }\ - } while(0) - struct { int *a, *b, *c; int d; } stack[STACK_SIZE]; - int *l, *r, *lm, *rm; - int m, len, half; - int ssize; - int check, next; - - for(check = 0, ssize = 0;;) { - if((last - middle) <= bufsize) { - if((first < middle) && (middle < last)) { - ss_mergebackward(T, PA, first, middle, last, buf, depth); - } - MERGE_CHECK(first, last, check); - STACK_POP(first, middle, last, check); - continue; - } - - if((middle - first) <= bufsize) { - if(first < middle) { - ss_mergeforward(T, PA, first, middle, last, buf, depth); - } - MERGE_CHECK(first, last, check); - STACK_POP(first, middle, last, check); - continue; - } - - for(m = 0, len = MIN(middle - first, last - middle), half = len >> 1; - 0 < len; - len = half, half >>= 1) { - if(ss_compare(T, PA + GETIDX(*(middle + m + half)), - PA + GETIDX(*(middle - m - half - 1)), depth) < 0) { - m += half + 1; - half -= (len & 1) ^ 1; - } - } - - if(0 < m) { - lm = middle - m, rm = middle + m; - ss_blockswap(lm, middle, m); - l = r = middle, next = 0; - if(rm < last) { - if(*rm < 0) { - *rm = ~*rm; - if(first < lm) { for(; *--l < 0;) { } next |= 4; } - next |= 1; - } else if(first < lm) { - for(; *r < 0; ++r) { } - next |= 2; - } - } - - if((l - first) <= (last - r)) { - STACK_PUSH(r, rm, last, (next & 3) | (check & 4)); - middle = lm, last = l, check = (check & 3) | (next & 4); - } else { - if((next & 2) && (r == middle)) { next ^= 6; } - STACK_PUSH(first, lm, l, (check & 3) | (next & 4)); - first = r, middle = rm, check = (next & 3) | (check & 4); - } - } else { - if(ss_compare(T, PA + GETIDX(*(middle - 1)), PA + *middle, depth) == 0) { - *middle = ~*middle; - } - MERGE_CHECK(first, last, check); - STACK_POP(first, middle, last, check); - } - } -#undef STACK_SIZE -} - -#endif /* SS_BLOCKSIZE != 0 */ - - -/*---------------------------------------------------------------------------*/ - -/* Substring sort */ -static -void -sssort(const unsigned char *T, const int *PA, - int *first, int *last, - int *buf, int bufsize, - int depth, int n, int lastsuffix) { - int *a; -#if SS_BLOCKSIZE != 0 - int *b, *middle, *curbuf; - int j, k, curbufsize, limit; -#endif - int i; - - if(lastsuffix != 0) { ++first; } - -#if SS_BLOCKSIZE == 0 - ss_mintrosort(T, PA, first, last, depth); -#else - if((bufsize < SS_BLOCKSIZE) && - (bufsize < (last - first)) && - (bufsize < (limit = ss_isqrt(last - first)))) { - if(SS_BLOCKSIZE < limit) { limit = SS_BLOCKSIZE; } - buf = middle = last - limit, bufsize = limit; - } else { - middle = last, limit = 0; - } - for(a = first, i = 0; SS_BLOCKSIZE < (middle - a); a += SS_BLOCKSIZE, ++i) { -#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE - ss_mintrosort(T, PA, a, a + SS_BLOCKSIZE, depth); -#elif 1 < SS_BLOCKSIZE - ss_insertionsort(T, PA, a, a + SS_BLOCKSIZE, depth); -#endif - curbufsize = last - (a + SS_BLOCKSIZE); - curbuf = a + SS_BLOCKSIZE; - if(curbufsize <= bufsize) { curbufsize = bufsize, curbuf = buf; } - for(b = a, k = SS_BLOCKSIZE, j = i; j & 1; b -= k, k <<= 1, j >>= 1) { - ss_swapmerge(T, PA, b - k, b, b + k, curbuf, curbufsize, depth); - } - } -#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE - ss_mintrosort(T, PA, a, middle, depth); -#elif 1 < SS_BLOCKSIZE - ss_insertionsort(T, PA, a, middle, depth); -#endif - for(k = SS_BLOCKSIZE; i != 0; k <<= 1, i >>= 1) { - if(i & 1) { - ss_swapmerge(T, PA, a - k, a, middle, buf, bufsize, depth); - a -= k; - } - } - if(limit != 0) { -#if SS_INSERTIONSORT_THRESHOLD < SS_BLOCKSIZE - ss_mintrosort(T, PA, middle, last, depth); -#elif 1 < SS_BLOCKSIZE - ss_insertionsort(T, PA, middle, last, depth); -#endif - ss_inplacemerge(T, PA, first, middle, last, depth); - } -#endif - - if(lastsuffix != 0) { - /* Insert last type B* suffix. */ - int PAi[2]; PAi[0] = PA[*(first - 1)], PAi[1] = n - 2; - for(a = first, i = *(first - 1); - (a < last) && ((*a < 0) || (0 < ss_compare(T, &(PAi[0]), PA + *a, depth))); - ++a) { - *(a - 1) = *a; - } - *(a - 1) = i; - } -} - - -/*---------------------------------------------------------------------------*/ - -static INLINE -int -tr_ilg(int n) { - return (n & 0xffff0000) ? - ((n & 0xff000000) ? - 24 + lg_table[(n >> 24) & 0xff] : - 16 + lg_table[(n >> 16) & 0xff]) : - ((n & 0x0000ff00) ? - 8 + lg_table[(n >> 8) & 0xff] : - 0 + lg_table[(n >> 0) & 0xff]); -} - - -/*---------------------------------------------------------------------------*/ - -/* Simple insertionsort for small size groups. */ -static -void -tr_insertionsort(const int *ISAd, int *first, int *last) { - int *a, *b; - int t, r; - - for(a = first + 1; a < last; ++a) { - for(t = *a, b = a - 1; 0 > (r = ISAd[t] - ISAd[*b]);) { - do { *(b + 1) = *b; } while((first <= --b) && (*b < 0)); - if(b < first) { break; } - } - if(r == 0) { *b = ~*b; } - *(b + 1) = t; - } -} - - -/*---------------------------------------------------------------------------*/ - -static INLINE -void -tr_fixdown(const int *ISAd, int *SA, int i, int size) { - int j, k; - int v; - int c, d, e; - - for(v = SA[i], c = ISAd[v]; (j = 2 * i + 1) < size; SA[i] = SA[k], i = k) { - d = ISAd[SA[k = j++]]; - if(d < (e = ISAd[SA[j]])) { k = j; d = e; } - if(d <= c) { break; } - } - SA[i] = v; -} - -/* Simple top-down heapsort. */ -static -void -tr_heapsort(const int *ISAd, int *SA, int size) { - int i, m; - int t; - - m = size; - if((size % 2) == 0) { - m--; - if(ISAd[SA[m / 2]] < ISAd[SA[m]]) { SWAP(SA[m], SA[m / 2]); } - } - - for(i = m / 2 - 1; 0 <= i; --i) { tr_fixdown(ISAd, SA, i, m); } - if((size % 2) == 0) { SWAP(SA[0], SA[m]); tr_fixdown(ISAd, SA, 0, m); } - for(i = m - 1; 0 < i; --i) { - t = SA[0], SA[0] = SA[i]; - tr_fixdown(ISAd, SA, 0, i); - SA[i] = t; - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Returns the median of three elements. */ -static INLINE -int * -tr_median3(const int *ISAd, int *v1, int *v2, int *v3) { - int *t; - if(ISAd[*v1] > ISAd[*v2]) { SWAP(v1, v2); } - if(ISAd[*v2] > ISAd[*v3]) { - if(ISAd[*v1] > ISAd[*v3]) { return v1; } - else { return v3; } - } - return v2; -} - -/* Returns the median of five elements. */ -static INLINE -int * -tr_median5(const int *ISAd, - int *v1, int *v2, int *v3, int *v4, int *v5) { - int *t; - if(ISAd[*v2] > ISAd[*v3]) { SWAP(v2, v3); } - if(ISAd[*v4] > ISAd[*v5]) { SWAP(v4, v5); } - if(ISAd[*v2] > ISAd[*v4]) { SWAP(v2, v4); SWAP(v3, v5); } - if(ISAd[*v1] > ISAd[*v3]) { SWAP(v1, v3); } - if(ISAd[*v1] > ISAd[*v4]) { SWAP(v1, v4); SWAP(v3, v5); } - if(ISAd[*v3] > ISAd[*v4]) { return v4; } - return v3; -} - -/* Returns the pivot element. */ -static INLINE -int * -tr_pivot(const int *ISAd, int *first, int *last) { - int *middle; - int t; - - t = last - first; - middle = first + t / 2; - - if(t <= 512) { - if(t <= 32) { - return tr_median3(ISAd, first, middle, last - 1); - } else { - t >>= 2; - return tr_median5(ISAd, first, first + t, middle, last - 1 - t, last - 1); - } - } - t >>= 3; - first = tr_median3(ISAd, first, first + t, first + (t << 1)); - middle = tr_median3(ISAd, middle - t, middle, middle + t); - last = tr_median3(ISAd, last - 1 - (t << 1), last - 1 - t, last - 1); - return tr_median3(ISAd, first, middle, last); -} - - -/*---------------------------------------------------------------------------*/ - -typedef struct _trbudget_t trbudget_t; -struct _trbudget_t { - int chance; - int remain; - int incval; - int count; -}; - -static INLINE -void -trbudget_init(trbudget_t *budget, int chance, int incval) { - budget->chance = chance; - budget->remain = budget->incval = incval; -} - -static INLINE -int -trbudget_check(trbudget_t *budget, int size) { - if(size <= budget->remain) { budget->remain -= size; return 1; } - if(budget->chance == 0) { budget->count += size; return 0; } - budget->remain += budget->incval - size; - budget->chance -= 1; - return 1; -} - - -/*---------------------------------------------------------------------------*/ - -static INLINE -void -tr_partition(const int *ISAd, - int *first, int *middle, int *last, - int **pa, int **pb, int v) { - int *a, *b, *c, *d, *e, *f; - int t, s; - int x = 0; - - for(b = middle - 1; (++b < last) && ((x = ISAd[*b]) == v);) { } - if(((a = b) < last) && (x < v)) { - for(; (++b < last) && ((x = ISAd[*b]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - } - for(c = last; (b < --c) && ((x = ISAd[*c]) == v);) { } - if((b < (d = c)) && (x > v)) { - for(; (b < --c) && ((x = ISAd[*c]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - for(; b < c;) { - SWAP(*b, *c); - for(; (++b < c) && ((x = ISAd[*b]) <= v);) { - if(x == v) { SWAP(*b, *a); ++a; } - } - for(; (b < --c) && ((x = ISAd[*c]) >= v);) { - if(x == v) { SWAP(*c, *d); --d; } - } - } - - if(a <= d) { - c = b - 1; - if((s = a - first) > (t = b - a)) { s = t; } - for(e = first, f = b - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - if((s = d - c) > (t = last - d - 1)) { s = t; } - for(e = b, f = last - s; 0 < s; --s, ++e, ++f) { SWAP(*e, *f); } - first += (b - a), last -= (d - c); - } - *pa = first, *pb = last; -} - -static -void -tr_copy(int *ISA, const int *SA, - int *first, int *a, int *b, int *last, - int depth) { - /* sort suffixes of middle partition - by using sorted order of suffixes of left and right partition. */ - int *c, *d, *e; - int s, v; - - v = b - SA - 1; - for(c = first, d = a - 1; c <= d; ++c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *++d = s; - ISA[s] = d - SA; - } - } - for(c = last - 1, e = d + 1, d = b; e < d; --c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *--d = s; - ISA[s] = d - SA; - } - } -} - -static -void -tr_partialcopy(int *ISA, const int *SA, - int *first, int *a, int *b, int *last, - int depth) { - int *c, *d, *e; - int s, v; - int rank, lastrank, newrank = -1; - - v = b - SA - 1; - lastrank = -1; - for(c = first, d = a - 1; c <= d; ++c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *++d = s; - rank = ISA[s + depth]; - if(lastrank != rank) { lastrank = rank; newrank = d - SA; } - ISA[s] = newrank; - } - } - - lastrank = -1; - for(e = d; first <= e; --e) { - rank = ISA[*e]; - if(lastrank != rank) { lastrank = rank; newrank = e - SA; } - if(newrank != rank) { ISA[*e] = newrank; } - } - - lastrank = -1; - for(c = last - 1, e = d + 1, d = b; e < d; --c) { - if((0 <= (s = *c - depth)) && (ISA[s] == v)) { - *--d = s; - rank = ISA[s + depth]; - if(lastrank != rank) { lastrank = rank; newrank = d - SA; } - ISA[s] = newrank; - } - } -} - -static -void -tr_introsort(int *ISA, const int *ISAd, - int *SA, int *first, int *last, - trbudget_t *budget) { -#define STACK_SIZE TR_STACKSIZE - struct { const int *a; int *b, *c; int d, e; }stack[STACK_SIZE]; - int *a, *b, *c; - int t; - int v, x = 0; - int incr = ISAd - ISA; - int limit, next; - int ssize, trlink = -1; - - for(ssize = 0, limit = tr_ilg(last - first);;) { - - if(limit < 0) { - if(limit == -1) { - /* tandem repeat partition */ - tr_partition(ISAd - incr, first, first, last, &a, &b, last - SA - 1); - - /* update ranks */ - if(a < last) { - for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; } - } - if(b < last) { - for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } - } - - /* push */ - if(1 < (b - a)) { - STACK_PUSH5(NULL, a, b, 0, 0); - STACK_PUSH5(ISAd - incr, first, last, -2, trlink); - trlink = ssize - 2; - } - if((a - first) <= (last - b)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd, b, last, tr_ilg(last - b), trlink); - last = a, limit = tr_ilg(a - first); - } else if(1 < (last - b)) { - first = b, limit = tr_ilg(last - b); - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } else { - if(1 < (last - b)) { - STACK_PUSH5(ISAd, first, a, tr_ilg(a - first), trlink); - first = b, limit = tr_ilg(last - b); - } else if(1 < (a - first)) { - last = a, limit = tr_ilg(a - first); - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } else if(limit == -2) { - /* tandem repeat copy */ - a = stack[--ssize].b, b = stack[ssize].c; - if(stack[ssize].d == 0) { - tr_copy(ISA, SA, first, a, b, last, ISAd - ISA); - } else { - if(0 <= trlink) { stack[trlink].d = -1; } - tr_partialcopy(ISA, SA, first, a, b, last, ISAd - ISA); - } - STACK_POP5(ISAd, first, last, limit, trlink); - } else { - /* sorted partition */ - if(0 <= *first) { - a = first; - do { ISA[*a] = a - SA; } while((++a < last) && (0 <= *a)); - first = a; - } - if(first < last) { - a = first; do { *a = ~*a; } while(*++a < 0); - next = (ISA[*a] != ISAd[*a]) ? tr_ilg(a - first + 1) : -1; - if(++a < last) { for(b = first, v = a - SA - 1; b < a; ++b) { ISA[*b] = v; } } - - /* push */ - if(trbudget_check(budget, a - first)) { - if((a - first) <= (last - a)) { - STACK_PUSH5(ISAd, a, last, -3, trlink); - ISAd += incr, last = a, limit = next; - } else { - if(1 < (last - a)) { - STACK_PUSH5(ISAd + incr, first, a, next, trlink); - first = a, limit = -3; - } else { - ISAd += incr, last = a, limit = next; - } - } - } else { - if(0 <= trlink) { stack[trlink].d = -1; } - if(1 < (last - a)) { - first = a, limit = -3; - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - continue; - } - - if((last - first) <= TR_INSERTIONSORT_THRESHOLD) { - tr_insertionsort(ISAd, first, last); - limit = -3; - continue; - } - - if(limit-- == 0) { - tr_heapsort(ISAd, first, last - first); - for(a = last - 1; first < a; a = b) { - for(x = ISAd[*a], b = a - 1; (first <= b) && (ISAd[*b] == x); --b) { *b = ~*b; } - } - limit = -3; - continue; - } - - /* choose pivot */ - a = tr_pivot(ISAd, first, last); - SWAP(*first, *a); - v = ISAd[*first]; - - /* partition */ - tr_partition(ISAd, first, first + 1, last, &a, &b, v); - if((last - first) != (b - a)) { - next = (ISA[*a] != v) ? tr_ilg(b - a) : -1; - - /* update ranks */ - for(c = first, v = a - SA - 1; c < a; ++c) { ISA[*c] = v; } - if(b < last) { for(c = a, v = b - SA - 1; c < b; ++c) { ISA[*c] = v; } } - - /* push */ - if((1 < (b - a)) && (trbudget_check(budget, b - a))) { - if((a - first) <= (last - b)) { - if((last - b) <= (b - a)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - STACK_PUSH5(ISAd, b, last, limit, trlink); - last = a; - } else if(1 < (last - b)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - first = b; - } else { - ISAd += incr, first = a, last = b, limit = next; - } - } else if((a - first) <= (b - a)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd, b, last, limit, trlink); - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - last = a; - } else { - STACK_PUSH5(ISAd, b, last, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } else { - STACK_PUSH5(ISAd, b, last, limit, trlink); - STACK_PUSH5(ISAd, first, a, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } else { - if((a - first) <= (b - a)) { - if(1 < (last - b)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - STACK_PUSH5(ISAd, first, a, limit, trlink); - first = b; - } else if(1 < (a - first)) { - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - last = a; - } else { - ISAd += incr, first = a, last = b, limit = next; - } - } else if((last - b) <= (b - a)) { - if(1 < (last - b)) { - STACK_PUSH5(ISAd, first, a, limit, trlink); - STACK_PUSH5(ISAd + incr, a, b, next, trlink); - first = b; - } else { - STACK_PUSH5(ISAd, first, a, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } else { - STACK_PUSH5(ISAd, first, a, limit, trlink); - STACK_PUSH5(ISAd, b, last, limit, trlink); - ISAd += incr, first = a, last = b, limit = next; - } - } - } else { - if((1 < (b - a)) && (0 <= trlink)) { stack[trlink].d = -1; } - if((a - first) <= (last - b)) { - if(1 < (a - first)) { - STACK_PUSH5(ISAd, b, last, limit, trlink); - last = a; - } else if(1 < (last - b)) { - first = b; - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } else { - if(1 < (last - b)) { - STACK_PUSH5(ISAd, first, a, limit, trlink); - first = b; - } else if(1 < (a - first)) { - last = a; - } else { - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } - } else { - if(trbudget_check(budget, last - first)) { - limit = tr_ilg(last - first), ISAd += incr; - } else { - if(0 <= trlink) { stack[trlink].d = -1; } - STACK_POP5(ISAd, first, last, limit, trlink); - } - } - } -#undef STACK_SIZE -} - - - -/*---------------------------------------------------------------------------*/ - -/* Tandem repeat sort */ -static -void -trsort(int *ISA, int *SA, int n, int depth) { - int *ISAd; - int *first, *last; - trbudget_t budget; - int t, skip, unsorted; - - trbudget_init(&budget, tr_ilg(n) * 2 / 3, n); -/* trbudget_init(&budget, tr_ilg(n) * 3 / 4, n); */ - for(ISAd = ISA + depth; -n < *SA; ISAd += ISAd - ISA) { - first = SA; - skip = 0; - unsorted = 0; - do { - if((t = *first) < 0) { first -= t; skip += t; } - else { - if(skip != 0) { *(first + skip) = skip; skip = 0; } - last = SA + ISA[t] + 1; - if(1 < (last - first)) { - budget.count = 0; - tr_introsort(ISA, ISAd, SA, first, last, &budget); - if(budget.count != 0) { unsorted += budget.count; } - else { skip = first - last; } - } else if((last - first) == 1) { - skip = -1; - } - first = last; - } - } while(first < (SA + n)); - if(skip != 0) { *(first + skip) = skip; } - if(unsorted == 0) { break; } - } -} - - -/*---------------------------------------------------------------------------*/ - -/* Sorts suffixes of type B*. */ -static -int -sort_typeBstar(const unsigned char *T, int *SA, - int *bucket_A, int *bucket_B, - int n) { - int *PAb, *ISAb, *buf; -#ifdef _OPENMP - int *curbuf; - int l; -#endif - int i, j, k, t, m, bufsize; - int c0, c1; -#ifdef _OPENMP - int d0, d1; - int tmp; -#endif - - /* Initialize bucket arrays. */ - for(i = 0; i < BUCKET_A_SIZE; ++i) { bucket_A[i] = 0; } - for(i = 0; i < BUCKET_B_SIZE; ++i) { bucket_B[i] = 0; } - - /* Count the number of occurrences of the first one or two characters of each - type A, B and B* suffix. Moreover, store the beginning position of all - type B* suffixes into the array SA. */ - i = n - 1; - m = n; - c0 = T[n - 1]; - while ( 0 <= i) { - /* type A suffix. */ - do { - ++BUCKET_A(c1 = c0); - } while((0 <= --i) && ((c0 = T[i]) >= c1)); - if(0 <= i) { - /* type B* suffix. */ - ++BUCKET_BSTAR(c0, c1); - SA[--m] = i; - /* type B suffix. */ - for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { - ++BUCKET_B(c0, c1); - } - } - } - m = n - m; -/* -note: - A type B* suffix is lexicographically smaller than a type B suffix that - begins with the same first two characters. -*/ - - /* Calculate the index of start/end point of each bucket. */ - for(c0 = 0, i = 0, j = 0; c0 < ALPHABET_SIZE; ++c0) { - t = i + BUCKET_A(c0); - BUCKET_A(c0) = i + j; /* start point */ - i = t + BUCKET_B(c0, c0); - for(c1 = c0 + 1; c1 < ALPHABET_SIZE; ++c1) { - j += BUCKET_BSTAR(c0, c1); - BUCKET_BSTAR(c0, c1) = j; /* end point */ - i += BUCKET_B(c0, c1); - } - } - - if(0 < m) { - /* Sort the type B* suffixes by their first two characters. */ - PAb = SA + n - m; ISAb = SA + m; - for(i = m - 2; 0 <= i; --i) { - t = PAb[i], c0 = T[t], c1 = T[t + 1]; - SA[--BUCKET_BSTAR(c0, c1)] = i; - } - t = PAb[m - 1], c0 = T[t], c1 = T[t + 1]; - SA[--BUCKET_BSTAR(c0, c1)] = m - 1; - - /* Sort the type B* substrings using sssort. */ -#ifdef _OPENMP - tmp = omp_get_max_threads(); - buf = SA + m, bufsize = (n - (2 * m)) / tmp; - c0 = ALPHABET_SIZE - 2, c1 = ALPHABET_SIZE - 1, j = m; -#pragma omp parallel default(shared) private(curbuf, k, l, d0, d1, tmp) - { - tmp = omp_get_thread_num(); - curbuf = buf + tmp * bufsize; - k = 0; - for(;;) { - #pragma omp critical(sssort_lock) - { - if(0 < (l = j)) { - d0 = c0, d1 = c1; - do { - k = BUCKET_BSTAR(d0, d1); - if(--d1 <= d0) { - d1 = ALPHABET_SIZE - 1; - if(--d0 < 0) { break; } - } - } while(((l - k) <= 1) && (0 < (l = k))); - c0 = d0, c1 = d1, j = k; - } - } - if(l == 0) { break; } - sssort(T, PAb, SA + k, SA + l, - curbuf, bufsize, 2, n, *(SA + k) == (m - 1)); - } - } -#else - buf = SA + m, bufsize = n - (2 * m); - for(c0 = ALPHABET_SIZE - 2, j = m; 0 < j; --c0) { - for(c1 = ALPHABET_SIZE - 1; c0 < c1; j = i, --c1) { - i = BUCKET_BSTAR(c0, c1); - if(1 < (j - i)) { - sssort(T, PAb, SA + i, SA + j, - buf, bufsize, 2, n, *(SA + i) == (m - 1)); - } - } - } -#endif - - /* Compute ranks of type B* substrings. */ - for(i = m - 1; 0 <= i; --i) { - if(0 <= SA[i]) { - j = i; - do { ISAb[SA[i]] = i; } while((0 <= --i) && (0 <= SA[i])); - SA[i + 1] = i - j; - if(i <= 0) { break; } - } - j = i; - do { ISAb[SA[i] = ~SA[i]] = j; } while(SA[--i] < 0); - ISAb[SA[i]] = j; - } - - /* Construct the inverse suffix array of type B* suffixes using trsort. */ - trsort(ISAb, SA, m, 1); - - /* Set the sorted order of tyoe B* suffixes. */ - for(i = n - 1, j = m, c0 = T[n - 1]; 0 <= i;) { - for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) >= c1); --i, c1 = c0) { } - if(0 <= i) { - t = i; - for(--i, c1 = c0; (0 <= i) && ((c0 = T[i]) <= c1); --i, c1 = c0) { } - SA[ISAb[--j]] = ((t == 0) || (1 < (t - i))) ? t : ~t; - } - } - - /* Calculate the index of start/end point of each bucket. */ - BUCKET_B(ALPHABET_SIZE - 1, ALPHABET_SIZE - 1) = n; /* end point */ - for(c0 = ALPHABET_SIZE - 2, k = m - 1; 0 <= c0; --c0) { - i = BUCKET_A(c0 + 1) - 1; - for(c1 = ALPHABET_SIZE - 1; c0 < c1; --c1) { - t = i - BUCKET_B(c0, c1); - BUCKET_B(c0, c1) = i; /* end point */ - - /* Move all type B* suffixes to the correct position. */ - for(i = t, j = BUCKET_BSTAR(c0, c1); - j <= k; - --i, --k) { SA[i] = SA[k]; } - } - BUCKET_BSTAR(c0, c0 + 1) = i - BUCKET_B(c0, c0) + 1; /* start point */ - BUCKET_B(c0, c0) = i; /* end point */ - } - } - - return m; -} - -/* Constructs the suffix array by using the sorted order of type B* suffixes. */ -static -void -construct_SA(const unsigned char *T, int *SA, - int *bucket_A, int *bucket_B, - int n, int m) { - int *i, *j, *k; - int s; - int c0, c1, c2; - - if(0 < m) { - /* Construct the sorted order of type B suffixes by using - the sorted order of type B* suffixes. */ - for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) { - /* Scan the suffix array from right to left. */ - for(i = SA + BUCKET_BSTAR(c1, c1 + 1), - j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; - i <= j; - --j) { - if(0 < (s = *j)) { - assert(T[s] == c1); - assert(((s + 1) < n) && (T[s] <= T[s + 1])); - assert(T[s - 1] <= T[s]); - *j = ~s; - c0 = T[--s]; - if((0 < s) && (T[s - 1] > c0)) { s = ~s; } - if(c0 != c2) { - if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; } - k = SA + BUCKET_B(c2 = c0, c1); - } - assert(k < j); - *k-- = s; - } else { - assert(((s == 0) && (T[s] == c1)) || (s < 0)); - *j = ~s; - } - } - } - } - - /* Construct the suffix array by using - the sorted order of type B suffixes. */ - k = SA + BUCKET_A(c2 = T[n - 1]); - *k++ = (T[n - 2] < c2) ? ~(n - 1) : (n - 1); - /* Scan the suffix array from left to right. */ - for(i = SA, j = SA + n; i < j; ++i) { - if(0 < (s = *i)) { - assert(T[s - 1] >= T[s]); - c0 = T[--s]; - if((s == 0) || (T[s - 1] < c0)) { s = ~s; } - if(c0 != c2) { - BUCKET_A(c2) = k - SA; - k = SA + BUCKET_A(c2 = c0); - } - assert(i < k); - *k++ = s; - } else { - assert(s < 0); - *i = ~s; - } - } -} - -/* Constructs the burrows-wheeler transformed string directly - by using the sorted order of type B* suffixes. */ -static -int -construct_BWT(const unsigned char *T, int *SA, - int *bucket_A, int *bucket_B, - int n, int m) { - int *i, *j, *k, *orig; - int s; - int c0, c1, c2; - - if(0 < m) { - /* Construct the sorted order of type B suffixes by using - the sorted order of type B* suffixes. */ - for(c1 = ALPHABET_SIZE - 2; 0 <= c1; --c1) { - /* Scan the suffix array from right to left. */ - for(i = SA + BUCKET_BSTAR(c1, c1 + 1), - j = SA + BUCKET_A(c1 + 1) - 1, k = NULL, c2 = -1; - i <= j; - --j) { - if(0 < (s = *j)) { - assert(T[s] == c1); - assert(((s + 1) < n) && (T[s] <= T[s + 1])); - assert(T[s - 1] <= T[s]); - c0 = T[--s]; - *j = ~((int)c0); - if((0 < s) && (T[s - 1] > c0)) { s = ~s; } - if(c0 != c2) { - if(0 <= c2) { BUCKET_B(c2, c1) = k - SA; } - k = SA + BUCKET_B(c2 = c0, c1); - } - assert(k < j); - *k-- = s; - } else if(s != 0) { - *j = ~s; -#ifndef NDEBUG - } else { - assert(T[s] == c1); -#endif - } - } - } - } - - /* Construct the BWTed string by using - the sorted order of type B suffixes. */ - k = SA + BUCKET_A(c2 = T[n - 1]); - *k++ = (T[n - 2] < c2) ? ~((int)T[n - 2]) : (n - 1); - /* Scan the suffix array from left to right. */ - for(i = SA, j = SA + n, orig = SA; i < j; ++i) { - if(0 < (s = *i)) { - assert(T[s - 1] >= T[s]); - c0 = T[--s]; - *i = c0; - if((0 < s) && (T[s - 1] < c0)) { s = ~((int)T[s - 1]); } - if(c0 != c2) { - BUCKET_A(c2) = k - SA; - k = SA + BUCKET_A(c2 = c0); - } - assert(i < k); - *k++ = s; - } else if(s != 0) { - *i = ~s; - } else { - orig = i; - } - } - - return orig - SA; -} - - -/*---------------------------------------------------------------------------*/ - -/*- Function -*/ - -int -divsufsort(const unsigned char *T, int *SA, int n) { - int *bucket_A, *bucket_B; - int m; - int err = 0; - - /* Check arguments. */ - if((T == NULL) || (SA == NULL) || (n < 0)) { return -1; } - else if(n == 0) { return 0; } - else if(n == 1) { SA[0] = 0; return 0; } - else if(n == 2) { m = (T[0] < T[1]); SA[m ^ 1] = 0, SA[m] = 1; return 0; } - - bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int)); - bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int)); - - /* Suffixsort. */ - if((bucket_A != NULL) && (bucket_B != NULL)) { - m = sort_typeBstar(T, SA, bucket_A, bucket_B, n); - construct_SA(T, SA, bucket_A, bucket_B, n, m); - } else { - err = -2; - } - - free(bucket_B); - free(bucket_A); - - return err; -} - -/* Constructs the burrows-wheeler transformed string of a given string. -* @param T[0..n-1] The input string. -* @param U[0..n-1] The output string. (can be T) -* @param A[0..n-1] The temporary array. (can be NULL) -* @param n The length of the given string. -* @return The primary index if no error occurred, -1 or -2 otherwise. -*/ -int -divbwt(const unsigned char *T, unsigned char *U, int *A, int n) { - int *B; - int *bucket_A, *bucket_B; - int m, pidx, i; - - /* Check arguments. */ - if((T == NULL) || (U == NULL) || (n < 0)) { return -1; } - else if(n <= 1) { if(n == 1) { U[0] = T[0]; } return n; } - - if((B = A) == NULL) { B = (int *)malloc((size_t)(n + 1) * sizeof(int)); } - bucket_A = (int *)malloc(BUCKET_A_SIZE * sizeof(int)); - bucket_B = (int *)malloc(BUCKET_B_SIZE * sizeof(int)); - - /* Burrows-Wheeler Transform. */ - if((B != NULL) && (bucket_A != NULL) && (bucket_B != NULL)) { - m = sort_typeBstar(T, B, bucket_A, bucket_B, n); - pidx = construct_BWT(T, B, bucket_A, bucket_B, n, m); - - /* Copy to output string. */ - U[0] = T[n - 1]; - for(i = 0; i < pidx; ++i) { U[i + 1] = (unsigned char)B[i]; } - for(i += 1; i < n; ++i) { U[i] = (unsigned char)B[i]; } - pidx += 1; - } else { - pidx = -2; - } - - free(bucket_B); - free(bucket_A); - if(A == NULL) { free(B); } - - return pidx; -} diff --git a/libdivsufsort/divsufsort.h.in b/libdivsufsort/divsufsort.h.in deleted file mode 100644 index 8d8952e9..00000000 --- a/libdivsufsort/divsufsort.h.in +++ /dev/null @@ -1,63 +0,0 @@ -/* - * divsufsort.h for libdivsufsort-lite - * Copyright (c) 2003-2008 Yuta Mori All Rights Reserved. - * - * Permission is hereby granted, free of charge, to any person - * obtaining a copy of this software and associated documentation - * files (the "Software"), to deal in the Software without - * restriction, including without limitation the rights to use, - * copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following - * conditions: - * - * The above copyright notice and this permission notice shall be - * included in all copies or substantial portions of the Software. - * - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES - * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT - * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, - * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING - * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR - * OTHER DEALINGS IN THE SOFTWARE. - */ - -#ifndef _DIVSUFSORT_H -#define _DIVSUFSORT_H 1 - -#ifdef __cplusplus -extern "C" { -#endif /* __cplusplus */ - - -/*- Prototypes -*/ - -/** - * Constructs the suffix array of a given string. - * @param T[0..n-1] The input string. - * @param SA[0..n-1] The output array of suffixes. - * @param n The length of the given string. - * @return 0 if no error occurred, -1 or -2 otherwise. - */ -int -divsufsort(const unsigned char *T, int *SA, int n); - -/** - * Constructs the burrows-wheeler transformed string of a given string. - * @param T[0..n-1] The input string. - * @param U[0..n-1] The output string. (can be T) - * @param A[0..n-1] The temporary array. (can be NULL) - * @param n The length of the given string. - * @return The primary index if no error occurred, -1 or -2 otherwise. - */ -int -divbwt(const unsigned char *T, unsigned char *U, int *A, int n); - - -#ifdef __cplusplus -} /* extern "C" */ -#endif /* __cplusplus */ - -#endif /* _DIVSUFSORT_H */ diff --git a/src/Makefile.in b/src/Makefile.in index 6beceb19..bd501b24 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -29,7 +29,7 @@ AVX_CFLAGS = @AVX_CFLAGS@ CPPFLAGS = @CPPFLAGS@ LDFLAGS = @LDFLAGS@ DEFS = @DEFS@ -LIBS = -lhmmer -leasel -ldivsufsort @LIBS@ @LIBGSL@ @PTHREAD_LIBS@ -lm +LIBS = -lhmmer -leasel @LIBS@ @LIBGSL@ @PTHREAD_LIBS@ -lm AR = @AR@ RANLIB = @RANLIB@ @@ -38,11 +38,9 @@ INSTALL = @INSTALL@ IMPLDIR = impl_@IMPL_CHOICE@ ESLDIR = @HMMER_ESLDIR@ -SADIR = @HMMER_SADIR@ -MYLIBDIRS = -L../${ESLDIR} -L./${IMPLDIR} -L../${SADIR} -L. +MYLIBDIRS = -L../${ESLDIR} -L./${IMPLDIR} -L. MYINCDIRS = -I../${ESLDIR} \ - -I../${SADIR} \ -I${top_srcdir}/easel \ -I. \ -I${srcdir}