[med-svn] [Git][med-team/phast][master] 6 commits: New upstream version 1.7+dfsg
Andreas Tille (@tille)
gitlab at salsa.debian.org
Fri Dec 6 20:35:49 GMT 2024
Andreas Tille pushed to branch master at Debian Med / phast
Commits:
204a70ae by Andreas Tille at 2024-12-06T21:26:04+01:00
New upstream version 1.7+dfsg
- - - - -
5db9636f by Andreas Tille at 2024-12-06T21:26:04+01:00
New upstream version
- - - - -
2033d1fc by Andreas Tille at 2024-12-06T21:26:07+01:00
Update upstream source from tag 'upstream/1.7+dfsg'
Update to upstream version '1.7+dfsg'
with Debian dir 0f980ac9fa47253a95471767a69236466f29c676
- - - - -
ed071d6d by Andreas Tille at 2024-12-06T21:26:11+01:00
Remove trailing whitespace in debian/changelog (routine-update)
- - - - -
52cdc94c by Andreas Tille at 2024-12-06T21:33:11+01:00
Refresh patches
- - - - -
7ba14524 by Andreas Tille at 2024-12-06T21:34:56+01:00
Upload to unstable
- - - - -
13 changed files:
- debian/changelog
- − debian/patches/fix_test_makefile.patch
- debian/patches/gcc-14.patch
- debian/patches/pcre2.patch
- debian/patches/series
- include/phast/fit_column.h
- src/dless/Makefile
- src/lib/base/phast_misc.c
- src/lib/phylo/phast_fit_column.c
- src/lib/phylo/phast_fit_feature.c
- src/prequel/Makefile
- test/Makefile
- test/hmrc_summary_correct
Changes:
=====================================
debian/changelog
=====================================
@@ -1,3 +1,10 @@
+phast (1.7+dfsg-1) unstable; urgency=medium
+
+ * New upstream version
+ * Remove trailing whitespace in debian/changelog (routine-update)
+
+ -- Andreas Tille <tille at debian.org> Fri, 06 Dec 2024 21:33:23 +0100
+
phast (1.6+dfsg-6) unstable; urgency=medium
* gcc-14.patch: fix build failure with gcc 14. (Closes: #1075383)
@@ -24,7 +31,7 @@ phast (1.6+dfsg-4) unstable; urgency=medium
* Provide alternative Build-Depends to liblapack-dev and libblas-dev
* Standards-Version: 4.6.2 (routine-update)
* Upstream email removed from d/copyright since its bouncing
-
+
[ Sébastien Villemot ]
* Complete migration away from ATLAS
Closes: #1056681
=====================================
debian/patches/fix_test_makefile.patch deleted
=====================================
@@ -1,143 +0,0 @@
-Author: Andreas Tille <tille at debian.org>
-Last-Update: Thu, 18 Oct 2018 18:16:40 +0200
-Description: Test suite makefile has syntax errors and there
- was at least one buggy data file which was fixed.
- .
- With these changes `make msa_view` works but other targets are
- throwing errors. This needs to be investigated
-Forwarded: https://github.com/CshlSiepelLab/phast/issues/47
-
---- a/test/Makefile
-+++ b/test/Makefile
-@@ -7,30 +7,30 @@ all: msa_view phyloFit phastCons
- msa_view:
- @echo "*** Testing msa_view ***"
- msa_view hmrc.ss -i SS --end 10000 > hmrc.fa
-- @if [[ -n `diff --brief hmrc.fa hmrc_correct.fa` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hmrc.fa hmrc_correct.fa ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.fa -o PHYLIP > hmrc.ph
- msa_view hmrc.ph -i PHYLIP > hmrc.fa
-- @if [[ -n `diff --brief hmrc.fa hmrc_correct.fa` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hmrc.fa hmrc_correct.fa ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.ph --out-format MPM --in-format PHYLIP > hmrc.mpm
- msa_view hmrc.mpm --in-format MPM > hmrc.fa
-- @if [[ -n `diff --brief hmrc.fa hmrc_correct.fa` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hmrc.fa hmrc_correct.fa ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.fa -o SS > hmrc_short_a.ss
- msa_view hmrc.ss --end 10000 -i SS -o SS > hmrc_short_b.ss
-- @if [[ -n `diff --brief hmrc_short_[ab].ss` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hmrc_short_[ab].ss ; then echo "ERROR" ; exit 1 ; fi
- msa_view --seqs human,cow hmrc_short_a.ss -i SS -o SS | msa_view - -i SS > hm_a.fa
- msa_view --seqs human,cow hmrc.fa > hm_b.fa
-- @if [[ -n `diff --brief hm_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hm_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.ss --gap-strip ANY -i SS -o SS | msa_view - -i SS > hmrc_nogaps_a.fa
- msa_view hmrc.ss -i SS | msa_view - --gap-strip ANY > hmrc_nogaps_b.fa
-- @if [[ -n `diff --brief hmrc_nogaps_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hmrc_nogaps_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.fa --seqs human,cow --gap-strip ALL > hm_a.fa
- msa_view hmrc_short_a.ss --seqs human,cow -i SS --gap-strip ALL > hm_b.fa
-- @if [[ -n `diff --brief hm_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hm_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.ss --start 10000 --end 20000 -i SS > hmrc_sub_a.fa
- msa_view hmrc.ss -i SS | msa_view - --start 10000 --end 20000 > hmrc_sub_b.fa
-- @if [[ -n `diff --brief hmrc_sub_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief hmrc_sub_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.ss --summary -i SS > hmrc_summary
-- @if [[ -n `diff --brief hmrc_summary hmrc_summary_correct` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief -b hmrc_summary hmrc_summary_correct ; then echo "ERROR" ; exit 1 ; fi
- @echo -e "Passed all tests.\n"
- @rm -f hmrc.ph hmrc.fa hmrc.mpm hmrc_short_[ab].ss hm_[ab].fa hmrc_nogaps_[ab].fa hmrc_sub_[ab].fa hmrc_summary
-
-@@ -43,40 +43,40 @@ msa_view:
- phyloFit:
- @echo "*** Testing phyloFit ***"
- phyloFit hmrc.ss --subst-mod JC69 --tree "(human, (mouse,rat), cow)" -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod jc.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod jc.mod ; then echo "Rounding errors found" ; diff -u phyloFit.mod jc.mod ; fi
- phyloFit hmrc.ss --subst-mod JC69 --tree "((((human,chimp), (mouse,rat)), cow), chicken)" -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod jc.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod jc.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod F81 --tree "(human, (mouse,rat), cow)" -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod f81.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod f81.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod HKY85 --tree "(human, (mouse,rat), cow)" -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod hky.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod hky.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat), cow)" -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod rev.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod UNREST --tree "(human, (mouse,rat), cow)" -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod unrest.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod unrest.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod HKY85 --tree "(human, (mouse,rat), cow)" -i SS -k 4 --quiet
-- @if [[ -n `diff --brief phyloFit.mod hky-dg.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod hky-dg.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat), cow)" -i SS -k 4 --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev-dg.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod rev-dg.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod HKY85 --tree "(human, (mouse,rat), cow)" -i SS --EM --quiet
-- @if [[ -n `diff --brief phyloFit.mod hky-em.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod hky-em.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat), cow)" -i SS --EM --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev-em.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod rev-em.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hpmrc.ss --subst-mod REV --tree "(hg16, (mm3,rn3), galGal2)" -i SS --gaps-as-bases --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev-gaps.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod rev-gaps.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod REV -i SS --init-model rev.mod --post-probs --lnl --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev-lnl.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-- @if [[ -n `diff --brief phyloFit.postprob rev.postprob` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod rev-lnl.mod ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.postprob rev.postprob ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat))" -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev-hmr.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod rev-hmr.mod ; then echo "ERROR" ; exit 1 ; fi
- msa_view hmrc.ss -i SS --seqs human,mouse,rat --unordered -o SS > hmr.ss
- phyloFit hmr.ss -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev-hmr2.mod` ]] ; then echo "ERROR" ; exit 1 ; fi # FIXME: rev-hmr.mod and rev-hmr2.mod should be equal (eq freqs being obtained from all seqs)
-+ if ! diff --brief phyloFit.mod rev-hmr2.mod ; then echo "ERROR" ; exit 1 ; fi # FIXME: rev-hmr.mod and rev-hmr2.mod should be equal (eq freqs being obtained from all seqs)
- msa_view hmrc.ss -i SS --seqs human,mouse --unordered -o SS > hm.ss
- phyloFit hm.ss -i SS --quiet
-- @if [[ -n `diff --brief phyloFit.mod rev-hm.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod rev-hm.mod ; then echo "ERROR" ; exit 1 ; fi
- phyloFit hmrc.ss --subst-mod UNREST --tree "((human, mouse), cow)" -i SS --ancestor cow --quiet
-- @if [[ -n `diff --brief phyloFit.mod unrest-cow-anc.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief phyloFit.mod unrest-cow-anc.mod ; then echo "ERROR" ; exit 1 ; fi
- @echo -e "Passed all tests.\n"
- @rm -f phyloFit.mod phyloFit.postprob hmr.ss hm.ss
-
-@@ -87,12 +87,12 @@ phyloFit:
- phastCons:
- @echo "*** Testing phastCons ***"
- phastCons hpmrc.ss hpmrc-rev-dg-global.mod --nrates 20 --transitions .08,.008 --quiet --viterbi elements.bed --seqname chr22 > cons.dat
-- @if [[ -n `diff --brief cons.dat cons_correct.dat` ]] ; then echo "ERROR" ; exit 1 ; fi
-- @if [[ -n `diff --brief elements.bed elements_correct.bed` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief cons.dat cons_correct.dat ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief elements.bed elements_correct.bed ; then echo "ERROR" ; exit 1 ; fi
- tree_doctor hpmrc-rev-dg-global.mod --prune galGal2 > hpmr.mod
- phastCons hpmrc.ss hpmr.mod --nrates 20 --transitions .08,.008 --quiet --viterbi elements-4way.bed --seqname chr22 > cons-4way.dat
-- @if [[ -n `diff --brief cons-4way.dat cons-4way_correct.dat` ]] ; then echo "ERROR" ; exit 1 ; fi
-- @if [[ -n `diff --brief elements-4way.bed elements-4way_correct.bed` ]] ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief cons-4way.dat cons-4way_correct.dat ; then echo "ERROR" ; exit 1 ; fi
-+ if ! diff --brief elements-4way.bed elements-4way_correct.bed ; then echo "ERROR" ; exit 1 ; fi
- @echo -e "Passed all tests.\n"
- @rm -f cons.dat cons-4way.dat elements.bed elements-4way.bed hpmr.mod
-
-@@ -110,7 +110,7 @@ phyloP:
- phyloP --null 10 phyloFit.mod > phyloP_null_test.txt
- phyloP -i SS phyloFit.mod hmrc.ss > phyloP_sph_test.txt
- phyloP -i SS --method LRT phyloFit.mod hmrc.ss > phyloP_lrt_test.txt
-- phyloP -i SS --method LRT --mode CONACC phyloFit.mod hmrc.ss > phyloP_lrt_conacc_test.txt
-+ phyloP -i SS --method LRT --mode CONACC phyloFit.mod hmrc.ss > phyloP_lrt_conacc_test.txt
- phyloP -i SS --method GERP phyloFit.mod hmrc.ss > phyloP_gerp_test.txt
- phyloP -i SS --method SCORE phyloFit.mod hmrc.ss > phyloP_score_test.txt
- phyloP -i SS --method LRT --wig-scores phyloFit.mod hmrc.ss > phyloP_wig_test.wig
---- a/test/hmrc_summary_correct
-+++ b/test/hmrc_summary_correct
-@@ -1,5 +1,2 @@
--descrip. A C G T G+C length all_gaps some_gaps
- descrip. A C G T G+C length all_gaps some_gaps
--
--hmrc.ss 0.3258 0.1913 0.1827 0.3001 0.3740 95927 0 82185
- hmrc.ss 0.3258 0.1913 0.1827 0.3001 0.3740 95927 0 82185
=====================================
debian/patches/gcc-14.patch
=====================================
@@ -12,9 +12,9 @@ Forwarded: no
Last-Update: 2024-07-31
---
This patch header follows DEP-3: http://dep.debian.net/deps/dep3/
---- phast.orig/include/phast/stringsplus.h
-+++ phast/include/phast/stringsplus.h
-@@ -383,7 +383,7 @@
+--- a/include/phast/stringsplus.h
++++ b/include/phast/stringsplus.h
+@@ -383,7 +383,7 @@ int str_split(String *s, const char* del
expression syntax.
@result Newly allocated and compiled Regex object.
*/
@@ -23,9 +23,9 @@ This patch header follows DEP-3: http://dep.debian.net/deps/dep3/
/** Free resources associated with regular expression object.
@param re Regex object to free
---- phast.orig/src/lib/base/phast_stringsplus.c
-+++ phast/src/lib/base/phast_stringsplus.c
-@@ -463,12 +463,12 @@
+--- a/src/lib/base/phast_stringsplus.c
++++ b/src/lib/base/phast_stringsplus.c
+@@ -463,12 +463,12 @@ int str_ends_with_charstr(String *s, con
return (strncmp(&s->chars[s->length - len], substr, len) == 0);
}
@@ -41,7 +41,7 @@ This patch header follows DEP-3: http://dep.debian.net/deps/dep3/
if (re == NULL) {
die("ERROR: cannot compile regular expression '%s' (%d): %d\n",
re_str, erroffset, errorcode);
-@@ -529,13 +529,13 @@
+@@ -529,13 +529,13 @@ int str_re_match_sub(String *s, pcre2_co
int str_re_match(String *s, pcre2_compile_context *re, List *l, int nsubexp) {
=====================================
debian/patches/pcre2.patch
=====================================
@@ -198,7 +198,7 @@ Forwarded: https://github.com/CshlSiepelLab/phast/issues/49
--- a/src/lib/base/phast_misc.c
+++ b/src/lib/base/phast_misc.c
-@@ -669,7 +669,7 @@ int draw_index(double *p, int size) {
+@@ -670,7 +670,7 @@ int draw_index(double *p, int size) {
character as well as "->" to indicate mapping. */
struct hash_table *make_name_hash(char *mapstr) {
Hashtable *retval = hsh_new(20);
=====================================
debian/patches/series
=====================================
@@ -3,6 +3,5 @@ do_not_install_to_opt.patch
use_debian_packaged_libpcre.patch
use_debian_packaged_help2man.patch
hardening.patch
-fix_test_makefile.patch
pcre2.patch
gcc-14.patch
=====================================
include/phast/fit_column.h
=====================================
@@ -143,9 +143,8 @@ double col_likelihood_wrapper_1d(double x, void *data);
@note Uses log rather than log2
*/
-double col_compute_log_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
- double **scratch);
-
+double col_compute_scaled_log_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
+ double **scratch);
/** Compute and return the likelihood of a tree model with respect
to a single column tuple in an alignment.
@@ -229,7 +228,7 @@ double col_scale_derivs(ColFitData *d, double *first_deriv,
/** Compute the first and (optionally) second derivatives with respect
to the scale parameters for the single-column log likelihood
- function (col_compute_log_likelihood).
+ function (col_compute_scaled_log_likelihood).
@param[in] d Column Data to analyze
@param[out] gradient gradient from first partial derivative
@param[out] hessian (optional) second order partial derivative
=====================================
src/dless/Makefile
=====================================
@@ -3,7 +3,7 @@ PHAST := ${PHAST}/..
#EXEC = $(addprefix ${BIN}/,$(notdir ${PWD}))
PROGS = dless dlessP
-MODULES = bd_phylo_hmm
+MODULES = phast_bd_phylo_hmm
EXEC = $(addprefix ${BIN}/,${PROGS})
# assume all *.c files are source
=====================================
src/lib/base/phast_misc.c
=====================================
@@ -495,6 +495,7 @@ double normalize_probs(double *p, int size) {
int i;
double sum = 0;
for (i = 0; i < size; i++) sum += p[i];
+ if(sum == 0) return sum;
for (i = 0; i < size; i++) p[i] /= sum;
return sum;
}
=====================================
src/lib/phylo/phast_fit_column.c
=====================================
@@ -14,6 +14,7 @@
perform single-base LRTs, score tests, phyloP, etc. */
#include <stdlib.h>
+#include <float.h>
#include <phast/fit_column.h>
#include <phast/sufficient_stats.h>
#include <phast/tree_likelihoods.h>
@@ -37,7 +38,7 @@
sufficient stats already available. Note that this function uses
natural log rather than log2. This function does allow for rate
variation. */
-double col_compute_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
+double col_compute_scaled_log_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
double **scratch) {
int i, j, k, nodeidx, rcat;
@@ -46,15 +47,17 @@ double col_compute_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
double total_prob = 0;
List *traversal = tr_postorder(mod->tree);
double **pL = NULL;
+ double log_scale = 0;
+ double scaling_threshold = DBL_MIN;
if (msa->ss->tuple_size != 1)
- die("ERROR col_compute_likelihood: need tuple size 1, got %i\n",
+ die("ERROR col_compute_scaled_log_likelihood: need tuple size 1, got %i\n",
msa->ss->tuple_size);
if (mod->order != 0)
- die("ERROR col_compute_likelihood: got mod->order of %i, expected 0\n",
+ die("ERROR col_compute_scaled_log_likelihood: got mod->order of %i, expected 0\n",
mod->order);
if (!mod->allow_gaps)
- die("ERROR col_compute_likelihood: need mod->allow_gaps to be TRUE\n");
+ die("ERROR col_compute_scaled_log_likelihood: need mod->allow_gaps to be TRUE\n");
/* allocate memory or use scratch if avail */
if (scratch != NULL)
@@ -93,8 +96,14 @@ double col_compute_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
for (k = 0; k < nstates; k++)
totr += pL[k][n->rchild->id] *
mm_get(rsubst_mat, i, k);
-
- pL[i][n->id] = totl * totr;
+
+ if (totl * totr < scaling_threshold) {
+ pL[i][n->id] = (totl / scaling_threshold) * totr;
+ log_scale -= log(scaling_threshold);
+ }
+ else {
+ pL[i][n->id] = totl * totr;
+ }
}
}
}
@@ -110,20 +119,9 @@ double col_compute_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
sfree(pL);
}
- return(total_prob);
-}
-
-
-
-/* See col_compute_likelihood above for notes.
- Note that this function uses natural log rather than log2
- */
-double col_compute_log_likelihood(TreeModel *mod, MSA *msa, int tupleidx,
- double **scratch) {
- return log(col_compute_likelihood(mod, msa, tupleidx, scratch));
+ return log(total_prob) - log_scale;
}
-
/* version of col_scale_derivs_subst that allows for the general case
of complex eigenvalues and eigenvectors */
void col_scale_derivs_subst_complex(ColFitData *d) {
@@ -298,7 +296,7 @@ void col_scale_derivs_subst(ColFitData *d) {
/* Compute the first and (optionally) second derivatives with respect
to the scale parameter for the single-column log likelihood
- function (col_compute_log_likelihood). This version assumes a
+ function (col_compute_scaled_log_likelihood). This version assumes a
single scale parameter; see below for the subtree case. Return
value is log likelihood, which is computed as a by-product. Derivs
will be stored in *first_deriv and *second_deriv. If second_deriv
@@ -449,7 +447,7 @@ double col_scale_derivs(ColFitData *d, double *first_deriv,
/* Compute the first and (optionally) second derivatives with respect
to the scale parameters for the single-column log likelihood
- function (col_compute_log_likelihood). This version assumes scale
+ function (col_compute_scaled_log_likelihood). This version assumes scale
parameters for the whole tree and for the subtree. Return value is
log likelihood, which is computed as a by-product. Derivs will be
stored in *gradient and *hessian. If hessian == NULL,
@@ -679,8 +677,8 @@ double col_likelihood_wrapper(Vector *params, void *data) {
/* reestimate subst models on edges */
tm_set_subst_matrices(d->mod);
- return -1 * col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ return -1 * col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
}
/* Wrapper for likelihood function for use in parameter estimation;
@@ -695,8 +693,8 @@ double col_likelihood_wrapper_1d(double x, void *data) {
/* reestimate subst models on edges */
tm_set_subst_matrices(d->mod);
- return -1 * col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ return -1 * col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
}
/* Wrapper for gradient function for use in parameter estimation */
@@ -771,7 +769,7 @@ void col_lrts(TreeModel *mod, MSA *msa, mode_type mode, double *tuple_pvals,
tm_set_subst_matrices(mod);
/* compute log likelihoods under null and alt hypotheses */
- null_lnl = col_compute_log_likelihood(mod, msa, i, d->fels_scratch[0]);
+ null_lnl = col_compute_scaled_log_likelihood(mod, msa, i, d->fels_scratch[0]);
vec_set(d->params, 0, d->init_scale);
d->tupleidx = i;
@@ -1368,16 +1366,16 @@ void col_scale_derivs_num(ColFitData *d, double *first_deriv,
double orig_scale = d->mod->scale;
d->mod->scale += (2*DERIV_EPSILON);
tm_set_subst_matrices(d->mod);
- lnl1 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl1 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
d->mod->scale = orig_scale + DERIV_EPSILON;
tm_set_subst_matrices(d->mod);
- lnl2 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl2 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
d->mod->scale = orig_scale;
tm_set_subst_matrices(d->mod);
- lnl3 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl3 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
*first_deriv = (lnl2 - lnl3) / DERIV_EPSILON;
*second_deriv = (lnl1 - 2*lnl2 + lnl3) / (DERIV_EPSILON * DERIV_EPSILON);
}
@@ -1391,37 +1389,37 @@ void col_scale_derivs_subtree_num(ColFitData *d, Vector *gradient,
d->mod->scale += (2*DERIV_EPSILON);
tm_set_subst_matrices(d->mod);
- lnl00 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl00 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
d->mod->scale = orig_scale + DERIV_EPSILON;
tm_set_subst_matrices(d->mod);
- lnl0 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl0 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
d->mod->scale = orig_scale + DERIV_EPSILON;
d->mod->scale_sub = orig_scale_sub + DERIV_EPSILON;
tm_set_subst_matrices(d->mod);
- lnl01 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl01 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
d->mod->scale = orig_scale;
d->mod->scale_sub = orig_scale_sub + DERIV_EPSILON;
tm_set_subst_matrices(d->mod);
- lnl1 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl1 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
d->mod->scale = orig_scale;
d->mod->scale_sub = orig_scale_sub + (2*DERIV_EPSILON);
tm_set_subst_matrices(d->mod);
- lnl11 = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl11 = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
d->mod->scale = orig_scale;
d->mod->scale_sub = orig_scale_sub;
tm_set_subst_matrices(d->mod);
- lnl = col_compute_log_likelihood(d->mod, d->msa, d->tupleidx,
- d->fels_scratch[0]);
+ lnl = col_compute_scaled_log_likelihood(d->mod, d->msa, d->tupleidx,
+ d->fels_scratch[0]);
gradient->data[0] = (lnl0 - lnl) / DERIV_EPSILON;
gradient->data[1] = (lnl1 - lnl) / DERIV_EPSILON;
@@ -1628,4 +1626,4 @@ int col_has_data_sub(TreeModel *mod, MSA *msa, int tupleidx, List *inside,
return TRUE;
return FALSE;
-}
+}
\ No newline at end of file
=====================================
src/lib/phylo/phast_fit_feature.c
=====================================
@@ -25,15 +25,15 @@
parameters (currently affects 1d parameter estimation only) */
/* Compute and return the log likelihood of a tree model with respect
- to a given feature. Calls col_compute_log_likelihood for each
+ to a given feature. Calls col_compute_scaled_log_likelihood for each
column */
double ff_compute_log_likelihood(TreeModel *mod, MSA *msa, GFF_Feature *feat,
double **scratch) {
double retval = 0;
int i;
for (i = feat->start-1; i < feat->end; i++) /* offset of one */
- retval += col_compute_log_likelihood(mod, msa, msa->ss->tuple_idx[i],
- scratch);
+ retval += col_compute_scaled_log_likelihood(mod, msa, msa->ss->tuple_idx[i],
+ scratch);
return retval;
}
=====================================
src/prequel/Makefile
=====================================
@@ -3,7 +3,7 @@ PHAST := ${PHAST}/..
# assume executable name is given by directory name
PROGS = pbsDecode pbsEncode pbsScoreMatrix pbsTrain prequel
-MODULES = simplex_grid pbs_code
+MODULES = phast_simplex_grid phast_pbs_code
EXEC = $(addprefix ${BIN}/,${PROGS})
# assume all *.c files are source
=====================================
test/Makefile
=====================================
@@ -7,30 +7,30 @@ all: msa_view phyloFit phastCons
msa_view:
@echo "*** Testing msa_view ***"
msa_view hmrc.ss -i SS --end 10000 > hmrc.fa
- @if [[ -n `diff --brief hmrc.fa hmrc_correct.fa` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hmrc.fa hmrc_correct.fa ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.fa -o PHYLIP > hmrc.ph
msa_view hmrc.ph -i PHYLIP > hmrc.fa
- @if [[ -n `diff --brief hmrc.fa hmrc_correct.fa` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hmrc.fa hmrc_correct.fa ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.ph --out-format MPM --in-format PHYLIP > hmrc.mpm
msa_view hmrc.mpm --in-format MPM > hmrc.fa
- @if [[ -n `diff --brief hmrc.fa hmrc_correct.fa` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hmrc.fa hmrc_correct.fa ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.fa -o SS > hmrc_short_a.ss
msa_view hmrc.ss --end 10000 -i SS -o SS > hmrc_short_b.ss
- @if [[ -n `diff --brief hmrc_short_[ab].ss` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hmrc_short_[ab].ss ; then echo "ERROR" ; exit 1 ; fi
msa_view --seqs human,cow hmrc_short_a.ss -i SS -o SS | msa_view - -i SS > hm_a.fa
msa_view --seqs human,cow hmrc.fa > hm_b.fa
- @if [[ -n `diff --brief hm_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hm_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.ss --gap-strip ANY -i SS -o SS | msa_view - -i SS > hmrc_nogaps_a.fa
msa_view hmrc.ss -i SS | msa_view - --gap-strip ANY > hmrc_nogaps_b.fa
- @if [[ -n `diff --brief hmrc_nogaps_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hmrc_nogaps_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.fa --seqs human,cow --gap-strip ALL > hm_a.fa
msa_view hmrc_short_a.ss --seqs human,cow -i SS --gap-strip ALL > hm_b.fa
- @if [[ -n `diff --brief hm_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hm_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.ss --start 10000 --end 20000 -i SS > hmrc_sub_a.fa
msa_view hmrc.ss -i SS | msa_view - --start 10000 --end 20000 > hmrc_sub_b.fa
- @if [[ -n `diff --brief hmrc_sub_[ab].fa` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief hmrc_sub_[ab].fa ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.ss --summary -i SS > hmrc_summary
- @if [[ -n `diff --brief hmrc_summary hmrc_summary_correct` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief -b hmrc_summary hmrc_summary_correct ; then echo "ERROR" ; exit 1 ; fi
@echo -e "Passed all tests.\n"
@rm -f hmrc.ph hmrc.fa hmrc.mpm hmrc_short_[ab].ss hm_[ab].fa hmrc_nogaps_[ab].fa hmrc_sub_[ab].fa hmrc_summary
@@ -43,40 +43,40 @@ msa_view:
phyloFit:
@echo "*** Testing phyloFit ***"
phyloFit hmrc.ss --subst-mod JC69 --tree "(human, (mouse,rat), cow)" -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod jc.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod jc.mod ; then echo "Rounding errors found" ; diff -u phyloFit.mod jc.mod ; fi
phyloFit hmrc.ss --subst-mod JC69 --tree "((((human,chimp), (mouse,rat)), cow), chicken)" -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod jc.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod jc.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod F81 --tree "(human, (mouse,rat), cow)" -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod f81.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod f81.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod HKY85 --tree "(human, (mouse,rat), cow)" -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod hky.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod hky.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat), cow)" -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod rev.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod rev.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod UNREST --tree "(human, (mouse,rat), cow)" -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod unrest.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod unrest.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod HKY85 --tree "(human, (mouse,rat), cow)" -i SS -k 4 --quiet
- @if [[ -n `diff --brief phyloFit.mod hky-dg.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod hky-dg.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat), cow)" -i SS -k 4 --quiet
- @if [[ -n `diff --brief phyloFit.mod rev-dg.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod rev-dg.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod HKY85 --tree "(human, (mouse,rat), cow)" -i SS --EM --quiet
- @if [[ -n `diff --brief phyloFit.mod hky-em.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod hky-em.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat), cow)" -i SS --EM --quiet
- @if [[ -n `diff --brief phyloFit.mod rev-em.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod rev-em.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hpmrc.ss --subst-mod REV --tree "(hg16, (mm3,rn3), galGal2)" -i SS --gaps-as-bases --quiet
- @if [[ -n `diff --brief phyloFit.mod rev-gaps.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod rev-gaps.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod REV -i SS --init-model rev.mod --post-probs --lnl --quiet
- @if [[ -n `diff --brief phyloFit.mod rev-lnl.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
- @if [[ -n `diff --brief phyloFit.postprob rev.postprob` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod rev-lnl.mod ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.postprob rev.postprob ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod REV --tree "(human, (mouse,rat))" -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod rev-hmr.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod rev-hmr.mod ; then echo "ERROR" ; exit 1 ; fi
msa_view hmrc.ss -i SS --seqs human,mouse,rat --unordered -o SS > hmr.ss
phyloFit hmr.ss -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod rev-hmr2.mod` ]] ; then echo "ERROR" ; exit 1 ; fi # FIXME: rev-hmr.mod and rev-hmr2.mod should be equal (eq freqs being obtained from all seqs)
+ if ! diff --brief phyloFit.mod rev-hmr2.mod ; then echo "ERROR" ; exit 1 ; fi # FIXME: rev-hmr.mod and rev-hmr2.mod should be equal (eq freqs being obtained from all seqs)
msa_view hmrc.ss -i SS --seqs human,mouse --unordered -o SS > hm.ss
phyloFit hm.ss -i SS --quiet
- @if [[ -n `diff --brief phyloFit.mod rev-hm.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod rev-hm.mod ; then echo "ERROR" ; exit 1 ; fi
phyloFit hmrc.ss --subst-mod UNREST --tree "((human, mouse), cow)" -i SS --ancestor cow --quiet
- @if [[ -n `diff --brief phyloFit.mod unrest-cow-anc.mod` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief phyloFit.mod unrest-cow-anc.mod ; then echo "ERROR" ; exit 1 ; fi
@echo -e "Passed all tests.\n"
@rm -f phyloFit.mod phyloFit.postprob hmr.ss hm.ss
@@ -87,12 +87,12 @@ phyloFit:
phastCons:
@echo "*** Testing phastCons ***"
phastCons hpmrc.ss hpmrc-rev-dg-global.mod --nrates 20 --transitions .08,.008 --quiet --viterbi elements.bed --seqname chr22 > cons.dat
- @if [[ -n `diff --brief cons.dat cons_correct.dat` ]] ; then echo "ERROR" ; exit 1 ; fi
- @if [[ -n `diff --brief elements.bed elements_correct.bed` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief cons.dat cons_correct.dat ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief elements.bed elements_correct.bed ; then echo "ERROR" ; exit 1 ; fi
tree_doctor hpmrc-rev-dg-global.mod --prune galGal2 > hpmr.mod
phastCons hpmrc.ss hpmr.mod --nrates 20 --transitions .08,.008 --quiet --viterbi elements-4way.bed --seqname chr22 > cons-4way.dat
- @if [[ -n `diff --brief cons-4way.dat cons-4way_correct.dat` ]] ; then echo "ERROR" ; exit 1 ; fi
- @if [[ -n `diff --brief elements-4way.bed elements-4way_correct.bed` ]] ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief cons-4way.dat cons-4way_correct.dat ; then echo "ERROR" ; exit 1 ; fi
+ if ! diff --brief elements-4way.bed elements-4way_correct.bed ; then echo "ERROR" ; exit 1 ; fi
@echo -e "Passed all tests.\n"
@rm -f cons.dat cons-4way.dat elements.bed elements-4way.bed hpmr.mod
@@ -110,7 +110,7 @@ phyloP:
phyloP --null 10 phyloFit.mod > phyloP_null_test.txt
phyloP -i SS phyloFit.mod hmrc.ss > phyloP_sph_test.txt
phyloP -i SS --method LRT phyloFit.mod hmrc.ss > phyloP_lrt_test.txt
- phyloP -i SS --method LRT --mode CONACC phyloFit.mod hmrc.ss > phyloP_lrt_conacc_test.txt
+ phyloP -i SS --method LRT --mode CONACC phyloFit.mod hmrc.ss > phyloP_lrt_conacc_test.txt
phyloP -i SS --method GERP phyloFit.mod hmrc.ss > phyloP_gerp_test.txt
phyloP -i SS --method SCORE phyloFit.mod hmrc.ss > phyloP_score_test.txt
phyloP -i SS --method LRT --wig-scores phyloFit.mod hmrc.ss > phyloP_wig_test.wig
=====================================
test/hmrc_summary_correct
=====================================
@@ -1,5 +1,2 @@
-descrip. A C G T G+C length all_gaps some_gaps
descrip. A C G T G+C length all_gaps some_gaps
-
-hmrc.ss 0.3258 0.1913 0.1827 0.3001 0.3740 95927 0 82185
hmrc.ss 0.3258 0.1913 0.1827 0.3001 0.3740 95927 0 82185
View it on GitLab: https://salsa.debian.org/med-team/phast/-/compare/90c7243d2182cb7010151855f07c0c5709809063...7ba14524e3549180e977ce78a4f533c0fdf748f7
--
View it on GitLab: https://salsa.debian.org/med-team/phast/-/compare/90c7243d2182cb7010151855f07c0c5709809063...7ba14524e3549180e977ce78a4f533c0fdf748f7
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20241206/bbc7241e/attachment-0001.htm>
More information about the debian-med-commit
mailing list