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}