[med-svn] [Git][med-team/bali-phy][upstream] New upstream version 3.1.5+dfsg
Benjamin Redelings
gitlab at salsa.debian.org
Wed Jun 13 23:50:57 BST 2018
Benjamin Redelings pushed to branch upstream at Debian Med / bali-phy
Commits:
e7366411 by Benjamin Redelings at 2018-06-13T18:42:12-04:00
New upstream version 3.1.5+dfsg
- - - - -
26 changed files:
- doc/README.itex.xml
- + doc/man/tree-tool.md
- + functions/Frequencies/uniform.json
- functions/Rates/free.json
- functions/fMutSel.json
- functions/fMutSel0.json
- functions/functions/dna.json
- meson.build
- modules/SModel.hs
- scripts/bp-analyze
- src/alignment/alignment-constraint.H
- src/bali-phy.cc
- src/builtins/SModel.cc
- src/computation/computation.cc
- src/dp/dp-matrix.H
- src/dp/dp-matrix.cc
- src/mcmc/sample-alignment.cc
- src/meson.build
- src/models/parameters.cc
- src/startup/A-T-model.H
- src/startup/A-T-model.cc
- src/startup/help.cc
- src/substitution/substitution.cc
- src/tools/alignment-draw.cc
- src/tools/alignments-diff.cc
- + src/tools/tree-tool.cc
Changes:
=====================================
doc/README.itex.xml
=====================================
--- a/doc/README.itex.xml
+++ b/doc/README.itex.xml
@@ -1,5 +1,5 @@
<!DOCTYPE book [
-<!ENTITY version "3.1" >
+<!ENTITY version "3.1.4" >
<!ENTITY source.file "bali-phy-&version;.tar.gz" >
<!ENTITY linux64.file "bali-phy-&version;-linux64.tar.gz" >
@@ -161,7 +161,8 @@ However, note that this might conflict with R installed from other places, such
<section><info><title>Install on Linux</title></info>
- <section><info><title>Install BAli-Phy using <command>apt-get</command> (recommended if using <emphasis>Debian: testing or unstable</emphasis>)</title></info>
+ <section><info><title>Install BAli-Phy using <command>apt-get</command></title></info>
+ BAli-Phy is available on Ubuntu <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://launchpad.net/ubuntu/+source/bali-phy/">("Cosmic Cuttlefish" or later)</link>, and Debian (<link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://packages.debian.org/search?keywords=bali-phy&searchon=names§ion=all">testing and unstable</link>).
<screen><prompt>%</prompt> <userinput>sudo apt-get install bali-phy</userinput></screen>
Check that the executable runs:
<screen><prompt>%</prompt> <userinput>bali-phy --version</userinput></screen>
=====================================
doc/man/tree-tool.md
=====================================
--- /dev/null
+++ b/doc/man/tree-tool.md
@@ -0,0 +1,35 @@
+% tree-tool(1)
+% Benjamin Redelings
+% Feb 2018
+
+# NAME
+
+**tree-tool** -- Perform various operations on Newick trees.
+
+# SYNOPSIS
+
+**tree-tool** _tree-file_ [OPTIONS]
+
+# DESCRIPTION
+
+Perform various operations on Newick trees.
+
+# GENERAL OPTIONS:
+**-h**, **--help**
+: produce help message
+
+**-v**, **--verbose**
+: Output more log messages on stderr.
+
+**--prune** _arg_
+: Comma-separated taxa to remove
+
+**--resolve**
+: Comma-separated taxa to remove
+
+
+# REPORTING BUGS:
+ BAli-Phy online help: <http://www.bali-phy.org/docs.php>.
+
+Please send bug reports to <bali-phy-users at googlegroups.com>.
+
=====================================
functions/Frequencies/uniform.json
=====================================
--- /dev/null
+++ b/functions/Frequencies/uniform.json
@@ -0,0 +1,13 @@
+{
+ "name": "Frequencies.uniform",
+ "result_type": "List[Pair[String,Double]]",
+ "call": "SModel.uniform_frequencies[a]",
+ "args": [
+ {
+ "arg_name": "a",
+ "arg_type": "a",
+ "default_value": "getAlphabet"
+ }
+ ],
+ "examples": ["hky85[pi=Frequencies.uniform]"]
+}
=====================================
functions/Rates/free.json
=====================================
--- a/functions/Rates/free.json
+++ b/functions/Rates/free.json
@@ -30,5 +30,13 @@
],
"title": "Free rates model",
"description":"Rate heterogeneity model where the rate and frequency of each category can be estimated from the data.",
- "examples": ["hky85+Rates.free[n=3]"]
+ "examples": ["hky85+Rates.free[n=3]"],
+ "citation":{"type": "article",
+ "title": "A space-time process model for the evolution of DNA sequences.",
+ "year": "1995",
+ "author": [{"name": "Yang, Ziheng"}],
+ "journal": {"name": "Genetics", "volume": "139", "number": "2", "pages": "993--1005"},
+ "identifier": [{"type":"pmid","id":"7713447"},
+ {"type":"pmcid","id":"PMC1206396"}]
+ }
}
=====================================
functions/fMutSel.json
=====================================
--- a/functions/fMutSel.json
+++ b/functions/fMutSel.json
@@ -6,24 +6,35 @@
{
"arg_name": "omega",
"arg_type": "Double",
- "default_value": "~uniform[0,1]"
+ "default_value": "~uniform[0,1]",
+ "description": "Excess dN/dS"
},
{
"arg_name": "w",
"arg_type": "List[Pair[String,Double]]",
- "default_value": "zip[letters[A],~iid[length[letters[A]],log_normal[0,0.5]]]"
+ "default_value": "zip[letters[A],~iid[length[letters[A]],log_normal[0,0.5]]]",
+ "description": "Fitness of codons"
},
{
"arg_name": "submodel",
"arg_type": "RevCTMC[a]",
- "default_value": "hky85_sym",
- "alphabet": "getNucleotides[A]"
+ "default_value": "hky85",
+ "alphabet": "getNucleotides[A]",
+ "description": "Model of neutral nucleotide substitution"
},
{
"arg_name": "A",
"arg_type": "Codons[a]",
- "default_value": "getAlphabet"
+ "default_value": "getAlphabet",
+ "description": "The alphabet"
}
],
- "extract": "all"
+ "extract": "all",
+ "citation":{"type": "article",
+ "title": "Mutation-Selection Models of Codon Substitution and Their Use to Estimate Selective Strengths on Codon Usage",
+ "year": "2008",
+ "author": [{"name": "Yang, Ziheng"},{"name": "Nielsen, Rasmus"}],
+ "journal": {"name": "Molecular Biology and Evolution", "volume": "25", "number": "3", "pages": "568--579"},
+ "identifier": [{"type":"doi","id":"10.1093/molbev/msm284"}]
+ }
}
=====================================
functions/fMutSel0.json
=====================================
--- a/functions/fMutSel0.json
+++ b/functions/fMutSel0.json
@@ -6,24 +6,35 @@
{
"arg_name": "omega",
"arg_type": "Double",
- "default_value": "~uniform[0,1]"
+ "default_value": "~uniform[0,1]",
+ "description": "Excess dN/dS"
},
{
"arg_name": "w",
"arg_type": "List[Pair[String,Double]]",
- "default_value": "zip[letters[getAminoAcids[A]],~iid[length[letters[getAminoAcids[A]]],log_normal[0,0.5]]]"
+ "default_value": "zip[letters[getAminoAcids[A]],~iid[length[letters[getAminoAcids[A]]],log_normal[0,0.5]]]",
+ "description": "Fitness of amino acids"
},
{
"arg_name": "submodel",
"arg_type": "RevCTMC[a]",
- "default_value": "hky85_sym",
- "alphabet": "getNucleotides[A]"
+ "default_value": "hky85",
+ "alphabet": "getNucleotides[A]",
+ "description": "Model of neutral nucleotide substitution"
},
{
"arg_name": "A",
"arg_type": "Codons[a]",
- "default_value": "getAlphabet"
+ "default_value": "getAlphabet",
+ "description": "The alphabet"
}
],
- "extract": "all"
+ "extract": "all",
+ "citation":{"type": "article",
+ "title": "Mutation-Selection Models of Codon Substitution and Their Use to Estimate Selective Strengths on Codon Usage",
+ "year": "2008",
+ "author": [{"name": "Yang, Ziheng"},{"name": "Nielsen, Rasmus"}],
+ "journal": {"name": "Molecular Biology and Evolution", "volume": "25", "number": "3", "pages": "568--579"},
+ "identifier": [{"type":"doi","id":"10.1093/molbev/msm284"}]
+ }
}
=====================================
functions/functions/dna.json
=====================================
--- a/functions/functions/dna.json
+++ b/functions/functions/dna.json
@@ -1,8 +1,8 @@
{
"name": "dna",
- "synonyms": [ "DNA" ],
+ "synonyms": ["DNA"],
"result_type": "DNA",
"call": "Alphabet.dna",
"args": [],
- "title": "The RNA alphabet"
+ "title": "The DNA alphabet"
}
=====================================
meson.build
=====================================
--- a/meson.build
+++ b/meson.build
@@ -1,5 +1,5 @@
project('bali-phy', ['cpp','c'],
- version: '3.1.4',
+ version: '3.1.5',
default_options : [
'buildtype=release',
'cpp_std=c++14'
=====================================
modules/SModel.hs
=====================================
--- a/modules/SModel.hs
+++ b/modules/SModel.hs
@@ -37,7 +37,7 @@ builtin calc_root_probability 7 "calc_root_probability" "SModel";
builtin bitmask_from_alignment 2 "bitmask_from_alignment" "Alignment";
builtin peel_leaf_branch_SEV 4 "peel_leaf_branch_SEV" "SModel";
builtin peel_internal_branch_SEV 4 "peel_internal_branch_SEV" "SModel";
-builtin calc_root_probability_SEV 4 "calc_root_probability_SEV" "SModel";
+builtin calc_root_probability_SEV 5 "calc_root_probability_SEV" "SModel";
builtin peel_likelihood_1 3 "peel_likelihood_1" "SModel";
builtin peel_likelihood_2 6 "peel_likelihood_2" "SModel";
@@ -130,6 +130,9 @@ frequency_matrix (MixtureModel d) = let {model = MixtureModel d}
frequency_matrix (MixtureModels (m:ms)) = frequency_matrix m;
--
+uniform_frequencies a = zip letters (repeat $ 1.0/(intToDouble n_letters)) where {letters = alphabet_letters a;
+ n_letters = alphabetSize a};
+
plus_f_equal_frequencies a = plus_f a (replicate n_letters (1.0/intToDouble n_letters)) where {n_letters=alphabetSize a};
jukes_cantor a = reversible_markov (equ a) (plus_f_equal_frequencies a);
@@ -227,13 +230,6 @@ get_element_freqs ((key,value):rest) x = if (key == x) then value else get_eleme
get_element_exchange [] x y = error ("No exchangeability specified for '" ++ x ++ "'");
get_element_exchange ((key,value):rest) x y = if key == x || key == y then value else get_element_exchange rest x y;
-constant_frequencies_model freqs a = sequence [get_element_freqs freqs l|l <- alphabet_letters a];
-
-perform_hash [] = return [];
-perform_hash ((x,y):xys) = do {y' <- y; xys' <- perform_hash xys; return ((x,y'):xys')};
-
-constant_frequencies_model2 freqs a = do {freqs' <- freqs; return [get_element_freqs freqs' l|l <- alphabet_letters a]};
-
select_element x ((key,value):rest) = if x == key then value else select_element x rest;
select_element x [] = error $ "Can't find element " ++ show x ++ " in dictionary!";
@@ -386,6 +382,6 @@ cached_conditional_likelihoods_SEV t seqs alpha ps f a =
}
in lc;
-peel_likelihood_SEV t cl f root = let {branches_in = map (reverseEdge t) (edgesOutOfNode t root);} in
- case branches_in of {[b1,b2,b3]-> calc_root_probability_SEV (cl!b1) (cl!b2) (cl!b3) f};
+peel_likelihood_SEV t cl f root counts = let {branches_in = map (reverseEdge t) (edgesOutOfNode t root);} in
+ case branches_in of {[b1,b2,b3]-> calc_root_probability_SEV (cl!b1) (cl!b2) (cl!b3) f counts};
}
=====================================
scripts/bp-analyze
=====================================
--- a/scripts/bp-analyze
+++ b/scripts/bp-analyze
@@ -126,6 +126,7 @@ for(my $i=1;$i<=$#commands;$i++)
print "WARNING: Commands differ!\n $commands[0]\n $commands[$i]\n";
}
}
+
@directories = get_header_attributes("directory", at out_files);
my $directories_differ=0;
for(my $i=1;$i<=$#directories;$i++)
@@ -133,6 +134,13 @@ for(my $i=1;$i<=$#directories;$i++)
$directories_differ=1 if ($directories[$i] ne $directories[0]);
}
+my @versions = get_all_versions();
+my $versions_differ=0;
+for(my $i=1;$i<=$#versions;$i++)
+{
+ $versions_differ=1 if ($versions[$i] ne $versions[0]);
+}
+
@subdirs = get_header_attributes("subdirectory", at out_files);
my $betas = get_cmdline_attribute("beta");
@@ -384,8 +392,8 @@ sub html_header
th {padding-left: 0.5em; padding-right:0.5em}
td {padding: 0.1em;}
- td {padding-left: 0.3em;}
- td {padding-right: 0.3em;}
+ .backlit2 td {padding-left: 0.5em;}
+ .backlit2 td {padding-right: 1.0em;}
#topbar {
background-color: rgb(201,217,233);
@@ -393,6 +401,8 @@ sub html_header
padding: 0.5em;
display: table;
width: 100%;
+ position: fixed;
+ top: 0;
}
#topbar #menu {
@@ -407,11 +417,12 @@ sub html_header
white-space: nowrap;
}
- .content {margin:1em}
+ .content {margin:1em; margin-top: 3em}
.backlit td {background: rgb(220,220,220);}
+ .anchor {position: relative; top: -3em}
h1 {font-size: 150%;}
- h2 {font-size: 130%; margin-top:2em; margin-bottom: 0.2em}
+ h2 {font-size: 130%; margin-top:2.5em; margin-bottom: 1em}
h3 {font-size: 110%; margin-bottom: 0.2em}// margin-top: 0.3em}
ul {margin-top: 0.4em;}
@@ -445,8 +456,8 @@ sub topbar
<span id="menu">
[<a href="#data">Data+Model</a>]
[<a href="#parameters">Parameters</a>]
- [<a href="#alignment">Alignment</a>]
[<a href="#topology">Phylogeny</a>]
+ [<a href="#alignment">Alignment</a>]
[<a href="#mixing">Mixing</a>]
[<a href="#analysis">Analysis</a>]
[<a href="#models">Models+Priors</a>]
@@ -496,7 +507,7 @@ sub print_models
{
my $model_name = $name.$i;
my ($assign, $model_html) = print_model($model);
- $section .= "<tr><td><a name=\"$model_name\" class=\"modelname\">$model_name</a></td><td>$assign</td><td>\n";
+ $section .= "<tr><td><a name=\"$model_name\" class=\"anchor\"></a><span class=\"modelname\">$model_name</span></td><td>$assign</td><td>\n";
$section .= $model_html;
$section .= "</td></tr>\n";
$i = $i + 1;
@@ -509,7 +520,7 @@ sub print_models
sub print_data_and_model
{
my $section = "";
- $section .= "<h2 style=\"clear:both\"><a name=\"data\">Data & Model</a></h2>\n";
+ $section .= "<h2 style=\"clear:both\"><a class=\"anchor\" name=\"data\"></a>Data & Model</h2>\n";
$section .= "<table class=\"backlit2\">\n";
$section .= "<tr><th>Partition</th><th>Sequences</th><th>Lengths</th><th>Alphabet</th><th>Substitution Model</th><th>Indel Model</th><th>Scale Model</th></tr>\n";
for(my $p=0;$p<=$#input_file_names;$p++)
@@ -592,22 +603,22 @@ sub print_data_and_model
sub print_model_section
{
- my $section = "<h2 style=\"clear:both\"><a name=\"models\">Model and priors</h2>\n";
+ my $section = "<h2 style=\"clear:both\"><a class=\"anchor\" name=\"models\"></a>Model and priors</h2>\n";
- $section .= "<h3 style=\"clear:both\"><a name=\"tree\">Tree (+priors)</a></h3>\n";
+ $section .= "<h3 style=\"clear:both\"><a class=\"anchor\" name=\"tree\"></a>Tree (+priors)</h3>\n";
# $section .= "<table class=\"backlit2\">\n";
$section .= "<table>\n";
$section .= "<tr><td class=\"modelname\">topology</td><td>$topology_prior</td></tr>" if (defined($topology_prior));
$section .= "<tr><td class=\"modelname\">branch lengths</td><td>$branch_prior</td></tr>" if (defined($branch_prior));
$section .= "</table>\n";
- $section .= "<h3 style=\"clear:both\"><a name=\"smodel\">Substitution model (+priors) </a></h3>\n";
+ $section .= "<h3 style=\"clear:both\"><a class=\"anchor\" name=\"smodel\"></a>Substitution model (+priors)</h3>\n";
$section .= print_models("S",\@smodels);
- $section .= "<h3 style=\"clear:both\"><a name=\"imodel\">Indel model (+priors)</a></h3>\n";
+ $section .= "<h3 style=\"clear:both\"><a class=\"anchor\" name=\"imodel\"></a>Indel model (+priors)</h3>\n";
$section .= print_models("I",\@imodels);
- $section .= "<h3 style=\"clear:both\"><a name=\"scales\">Scales (+priors)</h3>\n";
+ $section .= "<h3 style=\"clear:both\"><a class=\"anchor\" name=\"scales\"></a>Scales (+priors)</h3>\n";
$section .= print_models("scale",\@scale_models);
return $section;
@@ -617,14 +628,19 @@ sub section_analysis
{
my $section = "";
$section .= "<br/><hr/><br/>\n";
- $section .= "<h2><a name=\"analysis\">Analysis</a></h2>\n";
- $section .= "<p style=\"margin-bottom:0\"><b>command line</b>: $commands[0]</p>\n" if (!$commands_differ);
- $section .= "<p style=\"margin-top:0\"><b>directory</b>: $directories[0]</p>\n" if (!$directories_differ);
+ $section .= "<h2><a class=\"anchor\" name=\"analysis\"></a>Analysis</h2>\n";
+ $section .= "<p>" if (!$commands_differ || !$directories_differ || !$versions_differ);
+ $section .= "<b>command line</b>: $commands[0]</br>\n" if (!$commands_differ);
+ $section .= "<b>directory</b>: $directories[0]</br>\n" if (!$directories_differ);
+ $section .= "<b>version</b>: $versions[0]\n" if (!$versions_differ);
+ $section .= "</p>" if (!$commands_differ || !$directories_differ || !$versions_differ);
#$section .= '<table style="width:100%">'."\n";
#$section .= '<table class="backlit2 center" style="width:100%">'."\n";
$section .= '<table class="backlit2 center">'."\n";
- $section .= "<tr><th>chain #</th><th>burnin</th><th>subsample</th><th>Iterations (after burnin)</th>";
+ $section .= "<tr><th>chain #</th>";
+ $section .= "<th>version</th>" if ($versions_differ);
+ $section .= "<th>burnin</th><th>subsample</th><th>Iterations (after burnin)</th>";
$section .= "<th>command line</th>" if ($commands_differ);
$section .= "<th>subdirectory</th>";
$section .= "<th>directory</th>" if ($directories_differ);
@@ -634,6 +650,7 @@ sub section_analysis
$section .= "<tr>\n";
$section .= " <td>".($i+1)."</td>\n";
+$section .= " <td>$versions[$i]</td>\n" if ($versions_differ);
$section .= " <td>$burnin</td>\n";
$section .= " <td>$subsample</td>\n";
@@ -645,6 +662,7 @@ $section .= " <td>$commands[$i]</td>\n" if ($commands_differ);
$section .= " <td>$subdirs[$i]</td>\n";
$section .= " <td>$directories[$i]</td>\n" if ($directories_differ);
+
$section .= "</tr>\n";
}
$section .= "</table>\n";
@@ -714,14 +732,14 @@ sub html_svg
sub section_phylogeny_distribution
{
my $section = "";
+ $section .= "<h2><a class=\"anchor\" name=\"topology\"></a>Phylogeny Distribution</h2>\n";
$section .= html_svg("c-levels.svg","35%","",["r_floating_picture"]);
- $section .= "<h2><a name=\"topology\">Phylogeny Distribution</a></h2>\n";
$section .= html_svg("c50-tree.svg","25%","",["floating_picture"]);
$section .= '<table>'."\n";
- $section .= "<tr><td>Partition support: <a href=\"consensus\">Summary</a></td></tr>\n";
- $section .= "<tr><td><span title=\"How many partitions are supported at each level of Posterior Log Odds (LOD)?\">Partition support graph:</span> <a href=\"c-levels.svg\">SVG</a></td></tr>\n";
+ $section .= "<tr><td>Partition support: <a href=\"consensus\">Summary</a></td><td><a href=\"partitions.bs\">Across chains</a></td></tr>\n";
+# $section .= "<tr><td><span title=\"How many partitions are supported at each level of Posterior Log Odds (LOD)?\">Partition support graph:</span> <a href=\"c-levels.svg\">SVG</a></td></tr>\n";
$section .= "</table>\n";
$section .= "<table>\n";
@@ -772,7 +790,7 @@ sub section_alignment_distribution
{
return "" if ($n_partitions == 0);
- my $section .= "<h2 class=\"clear\"><a name=\"alignment\">Alignment Distribution</a></h2>\n";
+ my $section .= "<h2 class=\"clear\"><a class=\"anchor\" name=\"alignment\"></a>Alignment Distribution</h2>\n";
for(my $i=0;$i<$n_partitions;$i++)
{
@@ -837,27 +855,23 @@ sub section_mixing
{
my $section = "";
-#$section .= '<object class="r_floating_picture" data="partitions.SRQ.svg" type="image/svg+xml" height="200pt"></object>';
-$section .= '<img src="partitions.SRQ.png" class="r_floating_picture" alt="SRQ plot for support of each partition."/>';
-#$section .= '<object class="r_floating_picture" data="c50.SRQ.svg" type="image/svg+xml" height="200pt"></object>';
-#$section .= '<embed class="r_floating_picture" src="c50.SRQ.svg" type="image/svg+xml" height="200" />';
-#$section .= '<embed class="r_floating_picture" src="partitions.SRQ.svg" type="image/svg+xml" height="200" />';
-$section .= '<img src="c50.SRQ.png" class="r_floating_picture" alt="SRQ plot for supprt of 50% consensus tree."/>';
+ $section .= "<h2><a class=\"anchor\" name=\"mixing\"></a>Mixing</h2>\n";
- $section .= "<h2><a name=\"mixing\">Mixing</a></h2>\n";
+ $section .= '<table style="width:100%;clear:both"><tr>';
+ $section .= '<td style="vertical-align:top">';
$section .= "<ol>\n";
- $section .= "<li><a href=\"partitions.bs\">Partition uncertainty</a></li>\n";
- for my $srq (@SRQ) {
- $section .= "<li><a href=\"$srq.SRQ.png\">SRQ plot: $srq</a></li>\n";
- }
- $section .= '<li><a href="convergence-PP.pdf">Variation in split frequency estimates</a></li>'."\n" if (-f "Results/convergence-PP.pdf");
+# for my $srq (@SRQ) {
+# $section .= "<li><a href=\"$srq.SRQ.png\">SRQ plot: $srq</a></li>\n";
+# }
+# $section .= '<li><a href="convergence-PP.pdf">Variation in split frequency estimates</a></li>'."\n" if (-f "Results/convergence-PP.pdf");
$section .= "</ol>\n";
- $section .= "<table>";
- $section .= "<tr><th>burnin (scalar)</th><th>ESS (scalar)</th><th>ESS (partition)</th><th>ASDSF</th><th>MSDSF</th><th>PSRF-CI80%</th><th>PSRF-RCF</th></tr>";
- $section .= "<tr>";
+
+ $section .= "<p><b>Statistics:</b></p>";
+ $section .= "<table class=\"backlit2\">";
+# $section .= "<tr><th>burnin (scalar)</th><th>ESS (scalar)</th><th>ESS (partition)</th><th>ASDSF</th><th>MSDSF</th><th>PSRF-CI80%</th><th>PSRF-RCF</th></tr>";
my $burnin_before = "NA";
my $min_NE = "NA";
@@ -868,18 +882,18 @@ $section .= '<img src="c50.SRQ.png" class="r_floating_picture" alt="SRQ plot for
$min_NE = get_value_from_file('Results/Report','Ne >=');
}
- $section .= "<td>$burnin_before</td>";
- $section .= "<td>$min_NE</td>";
+ $section .= "<tr><td><b>scalar burnin</b></td><td>$burnin_before</td></tr>";
+ $section .= "<tr><td><b>scalar ESS</b></td><td>$min_NE</td></tr>";
my $min_NE_partition = get_value_from_file('Results/partitions.bs','min Ne =');
- $section .= "<td>$min_NE_partition</td>";
+ $section .= "<tr><td><b title=\"Effective Sample Size for bit-vectors of partition support (smallest)\">topological ESS</b></td><td>$min_NE_partition</td></tr>";
my $asdsf = get_value_from_file('Results/partitions.bs','ASDSF\[min=0.100\] =');
$asdsf = "NA" if (!defined($asdsf));
- $section .= "<td>$asdsf</td>";
+ $section .= "<tr><td><b title=\"Average Standard Deviation of Split Frequencies\">ASDSF</b></td><td>$asdsf</td></tr>";
my $msdsf = get_value_from_file('Results/partitions.bs','MSDSF =');
$msdsf = "NA" if (!defined($msdsf));
- $section .= "<td>$msdsf</td>";
+ $section .= "<tr><td><b title=\"Maximum Standard Deviation of Split Frequencies\">MSDSF</b></td><td>$msdsf</td></tr>";
my $psrf_80 = "NA";
my $psrf_rcf = "NA";
@@ -888,10 +902,9 @@ $section .= '<img src="c50.SRQ.png" class="r_floating_picture" alt="SRQ plot for
$psrf_80 = get_value_from_file('Results/Report','PSRF-80%CI <=');
$psrf_rcf = get_value_from_file('Results/Report','PSRF-RCF <=');
}
- $section .= "<td>$psrf_80</td>";
- $section .= "<td>$psrf_rcf</td>";
+ $section .= "<tr><td><b>PSRF CI80%</b></td><td>$psrf_80</td></tr>";
+ $section .= "<tr><td><b>PSRF RCF</b></td><td>$psrf_rcf</td></tr>";
- $section .= "</tr>";
$section .= "</table>\n";
###### What file should we show for MDS?
@@ -916,14 +929,27 @@ $section .= '<img src="c50.SRQ.png" class="r_floating_picture" alt="SRQ plot for
$MDS_figure_3d = "tree-3D1-1.points";
$title = "1 chain";
}
+ $section .= '</td>';
+#$section .= '<object class="r_floating_picture" data="partitions.SRQ.svg" type="image/svg+xml" height="200pt"></object>';
+ $section .= '<td>';
+# $section .= '</td>';
+#$section .= '<object class="r_floating_picture" data="c50.SRQ.svg" type="image/svg+xml" height="200pt"></object>';
+#$section .= '<embed class="r_floating_picture" src="c50.SRQ.svg" type="image/svg+xml" height="200" />';
+#$section .= '<embed class="r_floating_picture" src="partitions.SRQ.svg" type="image/svg+xml" height="200" />';
+# $section .= '<td>';
+ $section .= '<img src="c50.SRQ.png" alt="SRQ plot for supprt of 50% consensus tree."/>';
+ $section .= '<img src="partitions.SRQ.png" alt="SRQ plot for support of each partition."/>';
+ $section .= '</td>';
+ $section .= '</tr></table>';
+
###### Table of MDS versus
$section .= '<table style="width:100%;clear:both"><tr>';
$section .= '<td style="width:40%;vertical-align:top">';
$section .= "<h4 style=\"text-align:center\">Projection of RF distances for $title (<a href=\"https://doi.org/10.1080/10635150590946961\">Hillis et al 2005</a>)</h4>";
$section .= html_svg($MDS_figure,"90%","",[]) if (defined($MDS_figure));
- $section .= "<a href='file:${MDS_figure_3d}.html'>3D</a>";
+ $section .= "<a href='${MDS_figure_3d}.html'>3D</a>";
$section .= '</td>';
$section .= '<td style="width:40%;vertical-align:top">';
@@ -944,6 +970,7 @@ $section .= '<img src="c50.SRQ.png" class="r_floating_picture" alt="SRQ plot for
}
$section .= '</td>';
+
$section .= '</tr></table>';
my $tne_string = exec_show("pickout -n Ne < Results/partitions.bs");
@@ -968,7 +995,7 @@ sub section_scalar_variables
{
my $section = "";
if ($#var_names != -1) {
- $section .= "<h2 class=\"clear\"><a name=\"parameters\">Scalar variables</a></h2>\n";
+ $section .= "<h2 class=\"clear\"><a class=\"anchor\" name=\"parameters\"></a>Scalar variables</h2>\n";
$section .= "<table class=\"backlit2\">\n";
$section .= "<tr><th>Statistic</th><th>Median</th><th title=\"95% Bayesian Credible Interval\">95% BCI</th><th title=\"Auto-Correlation Time\">ACT</th><th title=\"Effective Sample Size\">ESS</th><th>burnin</th><th title=\"Potential Scale Reduction Factor based on width of 80% credible interval\">PSRF-CI80%</th><th>PSRF-RCF</th></tr>\n";
@@ -1613,7 +1640,7 @@ sub draw_alignments
}
if (! more_recent_than("Results/$alignment-diff.html","Results/$alignment-diff.AU")) {
- exec_show("alignment-draw Results/$alignment.fasta --scale=identity --AU Results/$alignment-diff.AU --show-ruler --color-scheme=diff[3]+contrast > Results/$alignment-diff.html");
+ exec_show("alignment-draw Results/$alignment.fasta --scale=identity --AU Results/$alignment-diff.AU --show-ruler --color-scheme=diff[1]+contrast > Results/$alignment-diff.html");
}
}
print "done.\n";
@@ -2253,6 +2280,43 @@ sub get_models_for_file
return [@models];
}
+sub get_version_for_file
+{
+ my $file = shift;
+ my $version = get_header_attribute("VERSION",[$file]);
+ $version = (split /\s+/,$version)[0];
+ return $version;
+}
+
+sub get_version_for_run_file
+{
+ my $file = shift;
+ my $j = get_json_from_file($file);
+ return ${$j}{'program'}{'version'};
+}
+
+sub get_all_versions
+{
+ return [] if ($#out_files == -1);
+ my @versions;
+
+ if (!@run_files)
+ {
+ foreach my $out_file (@out_files)
+ {
+ push @versions, get_version_for_file($out_file);
+ }
+ }
+ else
+ {
+ foreach my $run_file (@run_files)
+ {
+ push @versions, get_version_for_run_file($run_file);
+ }
+ }
+ return @versions;
+}
+
sub get_models
{
my $name1 = shift;
@@ -2882,8 +2946,8 @@ sub tree_MDS
my $outfile = "Results/tree-1-2.svg";
if (! more_recent_than($outfile, $tf1) || ! more_recent_than($outfile, $tf2))
{
- my $L1 = min([$N, get_n_lines($tf1)-$burnin]);
- my $L2 = min([$N, get_n_lines($tf2)-$burnin]);
+ my $L1 = min([$N, int((get_n_lines($tf1)-$burnin)/$subsample)]);
+ my $L2 = min([$N, int((get_n_lines($tf2)-$burnin)/$subsample)]);
my $matfile = "Results/tree-1-2.M";
# print "L1 = $L1 L2 = $L2\n";
exec_show("trees-distances matrix --max=$N --jitter=0.3 $subsample_string $skip $tf1 $tf2 > $matfile");
@@ -2904,9 +2968,9 @@ sub tree_MDS
my $tf2 = $tree_files[1];
my $tf3 = $tree_files[2];
my $N = 400;
- my $L1 = min([$N, get_n_lines($tf1)-$burnin]);
- my $L2 = min([$N, get_n_lines($tf2)-$burnin]);
- my $L3 = min([$N, get_n_lines($tf3)-$burnin]);
+ my $L1 = min([$N, int((get_n_lines($tf1)-$burnin)/$subsample)]);
+ my $L2 = min([$N, int((get_n_lines($tf2)-$burnin)/$subsample)]);
+ my $L3 = min([$N, int((get_n_lines($tf3)-$burnin)/$subsample)]);
my $matfile = "Results/tree-1-2-3.M";
my $outfile = "Results/tree-1-2-3.svg";
# print "L1 = $L1 L2 = $L2\n";
=====================================
src/alignment/alignment-constraint.H
=====================================
--- a/src/alignment/alignment-constraint.H
+++ b/src/alignment/alignment-constraint.H
@@ -1,3 +1,5 @@
+#ifndef ALIGNMENT_CONSTRAINT_H
+#define ALIGNMENT_CONSTRAINT_H
/*
Copyright (C) 2004-2007,2009 Benjamin Redelings
@@ -48,3 +50,5 @@ std::vector< std::pair<int,int> > get_y_ranges_for_band(int D,
std::vector< std::pair<int,int> > get_yboundaries_from_pins(int I, int J, const std::vector< std::vector<int> >& pins);
std::vector< std::pair<int,int> > boundaries_intersection(const std::vector< std::pair<int,int> >&,const std::vector< std::pair<int,int> >&);
+
+#endif
=====================================
src/bali-phy.cc
=====================================
--- a/src/bali-phy.cc
+++ b/src/bali-phy.cc
@@ -533,7 +533,8 @@ int main(int argc,char* argv[])
info["subdirectory"] = dir_name;
files = init_files(proc_id, dir_name, argc, argv);
loggers = construct_loggers(M, subsample, Rao_Blackwellize, proc_id, dir_name);
- write_initial_alignments(*M, proc_id, dir_name);
+
+ if (args.count("align")) write_initial_alignments(args, proc_id, dir_name);
}
else {
files.push_back(shared_ptr<ostream>(new ostream(cout.rdbuf())));
=====================================
src/builtins/SModel.cc
=====================================
--- a/src/builtins/SModel.cc
+++ b/src/builtins/SModel.cc
@@ -1110,7 +1110,8 @@ namespace substitution {
log_double_t calc_root_probability_SEV(const Likelihood_Cache_Branch* LCB1,
const Likelihood_Cache_Branch* LCB2,
const Likelihood_Cache_Branch* LCB3,
- const Matrix& F);
+ const Matrix& F,
+ const vector<int>& counts);
}
extern "C" closure builtin_function_calc_root_probability(OperationArgs& Args)
@@ -1139,11 +1140,13 @@ extern "C" closure builtin_function_calc_root_probability_SEV(OperationArgs& Arg
auto arg1 = Args.evaluate(1);
auto arg2 = Args.evaluate(2);
auto arg3 = Args.evaluate(3);
+ auto arg4 = Args.evaluate(4);
log_double_t Pr = substitution::calc_root_probability_SEV(&arg0.as_<Likelihood_Cache_Branch>(),
&arg1.as_<Likelihood_Cache_Branch>(),
&arg2.as_<Likelihood_Cache_Branch>(),
- arg3.as_<Box<Matrix>>());
+ arg3.as_<Box<Matrix>>(),
+ arg4.as_<Vector<int>>());
return {Pr};
}
=====================================
src/computation/computation.cc
=====================================
--- a/src/computation/computation.cc
+++ b/src/computation/computation.cc
@@ -12,7 +12,12 @@ int OperationArgs::reg_for_slot(int slot) const
int OperationArgs::n_args() const {return current_closure().exp.size();}
-const expression_ref& OperationArgs::reference(int slot) const {return current_closure().exp.sub()[slot];}
+const expression_ref& OperationArgs::reference(int slot) const
+{
+ assert(0 <= slot);
+ assert(slot < current_closure().exp.sub().size());
+ return current_closure().exp.sub()[slot];
+}
const closure& OperationArgs::evaluate_slot_to_closure(int slot)
{
=====================================
src/dp/dp-matrix.H
=====================================
--- a/src/dp/dp-matrix.H
+++ b/src/dp/dp-matrix.H
@@ -119,17 +119,9 @@ public:
void forward_first_cell(int,int);
virtual void forward_cell(int,int)=0;
- /// Compute the forward probabilities for a square
- void forward_square_first(int,int,int,int);
- void forward_square(int,int,int,int);
- void forward_square();
-
- /// Compute the forward probabilities for a square
+ /// Compute the forward probabilities between y1(x) and y2(x)
void forward_band(const std::vector< std::pair<int,int> >& boundaries);
- /// compute FP for entire matrix, with some points on path pinned
- void forward_constrained(const std::vector<std::vector<int> >&);
-
/// Sample a path from the HMM
std::vector<int> sample_path() const;
@@ -141,21 +133,6 @@ public:
};
-/// 2D Dynamic Programming Matrix for chains which only emit or don't emit
-class DPmatrixNoEmit: public DPmatrix {
-public:
- /// Compute the forward probabilities for a cell
- void forward_cell(int,int);
-
- log_double_t path_Q_subst(const std::vector<int>&) const {return 1;}
-
- using DPmatrix::DPmatrix;
-
- virtual ~DPmatrixNoEmit() {}
-};
-
-
-
/// 2D Dynamic Programming Matrix for chains which emit different things
class DPmatrixEmit : public DPmatrix {
protected:
@@ -202,7 +179,7 @@ public:
/// 2D Dynamic Programming matrix with no constraints on states at each cell
-class DPmatrixSimple: public DPmatrixEmit {
+class DPmatrixSimple final: public DPmatrixEmit {
public:
void forward_cell(int,int);
@@ -214,7 +191,7 @@ public:
/// Dynamic Programming matrix with constraints on the states
-class DPmatrixConstrained: public DPmatrixEmit
+class DPmatrixConstrained final: public DPmatrixEmit
{
int order_of_computation() const;
std::vector< std::vector<int> > allowed_states;
@@ -233,8 +210,6 @@ public:
void clear_cell(int,int);
void forward_cell(int,int);
- void prune();
-
DPmatrixConstrained(const HMM& M,
EmissionProbs&& d1,
EmissionProbs&& d2,
=====================================
src/dp/dp-matrix.cc
=====================================
--- a/src/dp/dp-matrix.cc
+++ b/src/dp/dp-matrix.cc
@@ -113,34 +113,6 @@ inline void DPmatrix::forward_first_cell(int i2,int j2)
}
}
-inline void DPmatrix::forward_square_first(int x1,int y1,int x2,int y2) {
- assert(0 < x1);
- assert(0 < y1);
- assert(x1 <= x2 or y1 <= y2);
- assert(x2 < size1());
- assert(y2 < size2());
-
- // Since we are using M(0,0) instead of S(0,0), we need to run only the silent states at (0,0)
- // We can only use non-silent states at (0,0) to simulate S
-
- // clear left border
- for(int y=y1;y<=y2;y++)
- clear_cell(x1-1,y);
-
- // forward first row, with exception for S(0,0)
- clear_cell(x1,y1-1);
- forward_first_cell(x1,y1);
- for(int y=y1+1;y<=y2;y++)
- forward_cell(x1,y);
-
- // forward other rows
- for(int x=x1+1;x<=x2;x++) {
- clear_cell(x,y1-1);
- for(int y=y1;y<=y2;y++)
- forward_cell(x,y);
- }
-}
-
void DPmatrix::forward_band(const vector< pair<int,int> >& yboundaries)
{
// note: (x,y) is located at (x+1,y+1) in the matrix.
@@ -207,27 +179,6 @@ void DPmatrix::forward_band(const vector< pair<int,int> >& yboundaries)
compute_Pr_sum_all_paths();
}
-inline void DPmatrix::forward_square(int x1,int y1,int x2,int y2) {
- assert(0 < x1);
- assert(0 < y1);
- assert(x1 <= x2 or y1 <= y2);
- assert(x2 < size1());
- assert(y2 < size2());
-
- // Since we are using M(0,0) instead of S(0,0), we need to run only the silent states at (0,0)
- // We can only use non-silent states at (0,0) to simulate S
-
- // clear left border
- for(int y=y1;y<=y2;y++)
- clear_cell(x1-1,y);
-
- for(int x=x1;x<=x2;x++) {
- clear_cell(x,y1-1);
- for(int y=y1;y<=y2;y++)
- forward_cell(x,y);
- }
-}
-
void DPmatrix::compute_Pr_sum_all_paths()
{
const int I = size1()-1;
@@ -244,50 +195,6 @@ void DPmatrix::compute_Pr_sum_all_paths()
assert(Pr_total <= 1.0);
}
-void DPmatrix::forward_square()
-{
- const int I = size1()-1;
- const int J = size2()-1;
-
- forward_square_first(1,1,I,J);
-
- compute_Pr_sum_all_paths();
-}
-
-// FIXME - fix up pins for new matrix coordinates
-void DPmatrix::forward_constrained(const vector< vector<int> >& pins)
-{
- using std::pair;
-
- const int I = size1()-1;
- const int J = size2()-1;
-
- vector< pair<int,int> > yboundaries = get_yboundaries_from_pins(I, J, pins);
-
- const vector<int>& x = pins[0];
- const vector<int>& y = pins[1];
-
- forward_band(yboundaries);
- return;
-
- if (x.size() == 0)
- forward_square();
- else
- {
- // Propogate from S to first pin
- forward_square_first(1,1,x[0],y[0]);
-
- // Propogate from first pin to other pins (if any)
- for(int i=0;i<(int)x.size()-1;i++)
- forward_square(x[i]+1,y[i]+1,x[i+1],y[i+1]);
-
- int p = x.size()-1;
- forward_square(x[p]+1,y[p]+1,I,J);
- }
-
- compute_Pr_sum_all_paths();
-}
-
log_double_t DPmatrix::path_P(const vector<int>& path) const
{
const int I = size1()-1;
@@ -415,60 +322,6 @@ DPmatrix::DPmatrix(int i1,
(*this)(I,J,state1) = 0;
}
-inline void DPmatrixNoEmit::forward_cell(int i2,int j2)
-{
- assert(0 < i2 and i2 < size1());
- assert(0 < j2 and j2 < size2());
-
- // determine initial scale for this cell
- scale(i2,j2) = max(scale(i2-1,j2), max( scale(i2-1,j2-1), scale(i2,j2-1) ) );
-
- double maximum = 0;
-
- for(int s2=0;s2<n_dp_states();s2++)
- {
- int S2 = dp_order(s2);
-
- //--- Get (i1,j1) from (i2,j2) and S2
- int i1 = i2;
- if (di(S2)) i1--;
-
- int j1 = j2;
- if (dj(S2)) j1--;
-
- //--- compute arrival probability ----
- int MAX = n_dp_states();
- if (not di(S2) and not dj(S2)) MAX = s2;
-
- double temp = 0;
- for(int s1=0;s1<MAX;s1++) {
- int S1 = dp_order(s1);
-
- temp += (*this)(i1,j1,S1) * GQ(S1,S2);
- }
-
- // rescale result to scale of this cell
- if (scale(i1,j1) != scale(i2,j2))
- temp *= pow2(scale(i1,j1)-scale(i2,j2));
-
- // record maximum
- if (temp > maximum) maximum = temp;
-
- // store the result
- (*this)(i2,j2,S2) = temp;
- }
-
- //------- if exponent is too high or too low, rescale ------//
- if (maximum > fp_scale::hi_cutoff or (maximum > 0 and maximum < fp_scale::lo_cutoff))
- {
- int logs = -(int)log2(maximum);
- double scale_ = pow2(logs);
- for(int S2=0;S2<n_dp_states();S2++)
- (*this)(i2,j2,S2) *= scale_;
- scale(i2,j2) -= logs;
- }
-}
-
inline double sum(const valarray<double>& v) {
return v.sum();
}
@@ -490,17 +343,8 @@ log_double_t DPmatrixEmit::path_Q_subst(const vector<int>& path) const
if (dj(state2))
j++;
- double sub;
if (di(state2) and dj(state2))
- sub = emitMM(i,j);
- else if (di(state2))
- sub = emitM_(i,j);
- else if (dj(state2))
- sub = emit_M(i,j);
- else
- sub = emit__(i,j);
-
- P_sub *= sub;
+ P_sub *= emitMM(i,j);
}
assert(i == size1()-1 and j == size2()-1);
return P_sub * Pr_extra_subst;
@@ -606,18 +450,8 @@ void DPmatrixSimple::forward_cell(int i2,int j2)
temp += (*this)(i1,j1,S1) * GQ(S1,S2);
//--- Include Emission Probability----
- double sub;
if (i1 != i2 and j1 != j2)
- sub = emitMM(i2,j2);
- else if (i1 != i2)
- sub = emitM_(i2,j2);
- else if (j1 != j2)
- sub = emit_M(i2,j2);
- else // silent state - nothing emitted
-
- sub = emit__(i2,j2);
-
- temp *= sub;
+ temp *= emitMM(i2,j2);
// rescale result to scale of this cell
if (scale(i1,j1) != scale(i2,j2))
@@ -687,17 +521,8 @@ inline void DPmatrixConstrained::forward_cell(int i2,int j2)
}
//--- Include Emission Probability----
- double sub;
if (i1 != i2 and j1 != j2)
- sub = emitMM(i2,j2);
- else if (i1 != i2)
- sub = emitM_(i2,j2);
- else if (j1 != j2)
- sub = emit_M(i2,j2);
- else // silent state - nothing emitted
- sub = emit__(i2,j2);
-
- temp *= sub;
+ temp *= emitMM(i2,j2);
// rescale result to scale of this cell
if (scale(i1,j1) != scale(i2,j2))
@@ -880,37 +705,6 @@ int DPmatrixConstrained::order_of_computation() const {
}
-void DPmatrixConstrained::prune() {
-
- std::abort();
-
- unsigned order1 = order_of_computation();
-
- // For column
- for(int c = 1;c<allowed_states.size();c++) {
-
- // and for each allowed state in that column
- for(int s2=allowed_states[c].size()-1;s2 >= 0;s2--) {
-
- // FIXME! this assumes that dj(s2) is true.
-
- // check to see whether any states in the previous column are connected to it
- bool used = false;
- for(int s1 = 0;s1<allowed_states[c-1].size() and not used;s1++) {
- if (connected(allowed_states[c-1][s1],allowed_states[c][s2]))
- used = true;
- }
-
- // if nothing is connected, then remove it
- if (not used)
- allowed_states[c].erase(allowed_states[c].begin()+s2);
- }
- }
-
- unsigned order2 = order_of_computation();
- std::cerr<<" order1 = "<<order1<<" order2 = "<<order2<<" fraction = "<<double(order2)/double(order1)<<endl;
-}
-
DPmatrixConstrained::DPmatrixConstrained(const HMM& M,
EmissionProbs&& d1,
EmissionProbs&& d2,
=====================================
src/mcmc/sample-alignment.cc
=====================================
--- a/src/mcmc/sample-alignment.cc
+++ b/src/mcmc/sample-alignment.cc
@@ -71,16 +71,18 @@ boost::shared_ptr<DPmatrixSimple> sample_alignment_forward(data_partition P, con
state_emit[2] |= (1<<0);
state_emit[3] |= 0;
+ int I = P.seqlength(t.source(b));
+ int J = P.seqlength(t.source(bb));
+
boost::shared_ptr<DPmatrixSimple>
Matrices( new DPmatrixSimple(HMM(state_emit, hmm.start_pi(), hmm, P.get_beta()),
std::move(dists1), std::move(dists2), P.WeightedFrequencyMatrix())
);
//------------------ Compute the DP matrix ---------------------//
- // vector<vector<int> > pins= get_pins(P.alignment_constraint,A,group1,~group1,seq1,seq2);
- vector<vector<int> > pins(2);
+ vector<std::pair<int,int>> yboundaries(I+1, std::pair<int,int>(0,J));
- Matrices->forward_constrained(pins);
+ Matrices->forward_band(yboundaries);
return Matrices;
}
=====================================
src/meson.build
=====================================
--- a/src/meson.build
+++ b/src/meson.build
@@ -92,19 +92,24 @@ subdir('builtins')
test('bali-phy version', baliphy, args:['--version'])
test('bali-phy help', baliphy, args:['--help'])
test('bali-phy 5d test', baliphy, args:[small_fasta,'--test', packagepath])
-test('bali-phy 5d +A 50', baliphy, args:[small_fasta,'--iter=50', packagepath])
-test('bali-phy 5d -A 200', baliphy, args:[small_fasta,'--iter=200', packagepath, '-Inone'])
+# When running on very slow autobuilders these tests could take a long time.
+test('bali-phy 5d +A 50', baliphy, args:[small_fasta,'--iter=50', packagepath], timeout: 180)
+test('bali-phy 5d -A 200', baliphy, args:[small_fasta,'--iter=200', packagepath, '-Inone'], timeout:120)
#--------- Build rules for tools that are always installed ------------#
simple_tools = ['model_P','statreport','stats-select','alignment-gild','alignment-consensus', 'alignment-max','alignment-chop-internal',
'alignment-indices','alignment-info','alignment-cat','alignment-translate','alignment-find','trees-consensus',
'tree-mean-lengths','mctree-mean-lengths','trees-to-SRQ','pickout','subsample','cut-range','trees-distances', 'alignment-thin',
- 'alignments-diff'
- ]
+ 'alignments-diff','tree-tool']
foreach tool: simple_tools
- tool_exe = executable(tool, ['tools/' + tool + '.cc'], dependencies:boost, link_with: libbaliphy, install: true)
+ tool_exe = executable(tool,
+ ['tools/' + tool + '.cc'],
+ dependencies:boost,
+ link_with: libbaliphy,
+ install_rpath: extra_rpath,
+ install: true)
test(tool + ' --help', tool_exe, args:['--help'])
endforeach
@@ -123,7 +128,12 @@ if get_option('extra-tools')
'alignment-convert', 'alignment-find-conserved','partitions-supported','draw-graph','trees-pair-distances','tree-partitions','tree-reroot',
'path-graph']
foreach tool : extra_tools
- executable( tool, ['tools/' + tool + '.cc'], dependencies: boost, link_with:libbaliphy, install:true)
+ executable( tool,
+ ['tools/' + tool + '.cc'],
+ dependencies: boost,
+ link_with:libbaliphy,
+ install_rpath: extra_rpath,
+ install:true)
endforeach
all_progs = all_progs + extra_tools
endif
=====================================
src/models/parameters.cc
=====================================
--- a/src/models/parameters.cc
+++ b/src/models/parameters.cc
@@ -542,22 +542,19 @@ data_partition_constants::data_partition_constants(Parameters* p, int i, const a
transition_p_method_indices[b] = p->add_compute_expression( {dummy("Prelude.!"), transition_ps, b} );
}
- // Add parameters for observed leaf sequence objects
- vector<vector<int>> seq_counts = alignment_letters_counts(AA, t.n_leaves(), counts);
- vector<expression_ref> counts_;
+ // Add expressions for leaf sequences
for(int i=0; i<leaf_sequence_indices.size(); i++)
- {
leaf_sequence_indices[i] = p->add_compute_expression(Vector<int>(sequences[i]));
- counts_.push_back( Vector<int>(seq_counts[i]) );
- }
+ // Add array of leaf sequences
vector<expression_ref> seqs_;
for(int index: leaf_sequence_indices)
+ {
seqs_.push_back( p->get_expression(index) );
+ }
auto seqs_array = p->get_expression( p->add_compute_expression({dummy("Prelude.listArray'"),get_list(seqs_)}) );
- auto counts_array = p->get_expression( p->add_compute_expression({dummy("Prelude.listArray'"),get_list(counts_)}) );
-
- // Add methods indices for sequence lengths
+
+ // Add array of pairwise alignments
vector<expression_ref> as_;
for(int b=0;b<2*B;b++)
{
@@ -566,6 +563,7 @@ data_partition_constants::data_partition_constants(Parameters* p, int i, const a
}
expression_ref as = p->get_expression( p->add_compute_expression({dummy("Prelude.listArray'"),get_list(as_)}) );
+ // Precompute sequence lengths
for(int n=0;n<t.n_nodes();n++)
{
expression_ref L = {dummy("Alignment.seqlength"), as, p->my_tree(), n};
@@ -580,6 +578,12 @@ data_partition_constants::data_partition_constants(Parameters* p, int i, const a
}
else if (likelihood_calculator == 0)
{
+ vector<vector<int>> seq_counts = alignment_letters_counts(AA, t.n_leaves(), counts);
+ vector<expression_ref> counts_;
+ for(int i=0; i<leaf_sequence_indices.size(); i++)
+ counts_.push_back( Vector<int>(seq_counts[i]) );
+ auto counts_array = p->get_expression( p->add_compute_expression({dummy("Prelude.listArray'"),get_list(counts_)}) );
+
auto t = p->my_tree();
auto f = p->get_expression(p->PC->SModels[smodel_index].weighted_frequency_matrix);
cl_index = p->add_compute_expression({dummy("SModel.cached_conditional_likelihoods"),t,seqs_array,counts_array,as,*a,transition_ps,f}); // Create and set conditional likelihoods for each branch
@@ -587,6 +591,7 @@ data_partition_constants::data_partition_constants(Parameters* p, int i, const a
for(int b=0;b<conditional_likelihoods_for_branch.size();b++)
conditional_likelihoods_for_branch[b] = p->add_compute_expression({dummy("Prelude.!"),cls,b});
+ // FIXME: broken for fixed alignments of 2 sequences.
if (p->t().n_nodes() == 2)
{
expression_ref seq1 = {dummy("Prelude.!"), seqs_array, 0};
@@ -612,6 +617,9 @@ data_partition_constants::data_partition_constants(Parameters* p, int i, const a
for(int b=0;b<conditional_likelihoods_for_branch.size();b++)
conditional_likelihoods_for_branch[b] = p->add_compute_expression({dummy("Prelude.!"),cls,b});
+ object_ptr<Vector<int>> Counts(new Vector<int>(counts));
+
+ // FIXME: broken for fixed alignments of 2 sequences.
if (p->t().n_nodes() == 2)
{
expression_ref seq1 = {dummy("Prelude.!"), seqs_array, 0};
@@ -625,7 +633,7 @@ data_partition_constants::data_partition_constants(Parameters* p, int i, const a
else
{
auto root = parameter("*subst_root");
- likelihood_index = p->add_compute_expression({dummy("SModel.peel_likelihood_SEV"), t, cls, f, root});
+ likelihood_index = p->add_compute_expression({dummy("SModel.peel_likelihood_SEV"), t, cls, f, root, Counts});
}
}
=====================================
src/startup/A-T-model.H
=====================================
--- a/src/startup/A-T-model.H
+++ b/src/startup/A-T-model.H
@@ -36,7 +36,7 @@ owned_ptr<Model> create_A_and_T_model(const Rules& R, boost::program_options::va
std::ostream& out_cache, std::ostream& out_screen, std::ostream& out_both, json& info,
int proc_id);
-void write_initial_alignments(const owned_ptr<Model>& M, int proc_id, std::string dir_name);
+void write_initial_alignments(boost::program_options::variables_map& args, int proc_id, const std::string& dir_name);
#endif
=====================================
src/startup/A-T-model.cc
=====================================
--- a/src/startup/A-T-model.cc
+++ b/src/startup/A-T-model.cc
@@ -11,6 +11,9 @@
#include "startup/help.hh"
#include "computation/expression/lambda.H"
#include <boost/optional/optional_io.hpp>
+#include <boost/filesystem/operations.hpp>
+
+namespace fs = boost::filesystem;
namespace po = boost::program_options;
using po::variables_map;
@@ -651,14 +654,14 @@ owned_ptr<Model> create_A_and_T_model(const Rules& R, variables_map& args, const
int c = likelihood_calculators[i];
if (c != 0 and c != 1)
throw myexception()<<"Calculator "<<c<<" not recognized!";
- if (c != 0 and imodel_mapping[i] != -1)
+ if (c != 0 and imodel_mapping[i])
throw myexception()<<"Calculator "<<c<<" does not work with a variable alignment!";
}
}
bool unalign = args.count("unalign");
for(int i=0;i<A.size();i++)
- if (unalign and imodel_mapping[i] != -1)
+ if (unalign and imodel_mapping[i])
if (likelihood_calculators[i] != 0)
throw myexception()<<"Can't unalign with calculator "<<likelihood_calculators[i]<<"!";
@@ -711,22 +714,21 @@ owned_ptr<Model> create_A_and_T_model(const Rules& R, variables_map& args, const
return P;
}
-void write_initial_alignments(const owned_ptr<Model>& M, int proc_id, string dir_name)
+void write_initial_alignments(variables_map& args, int proc_id, const string& dir_name)
{
- const Parameters* P = M.as<const Parameters>();
- if (not P) return;
+ vector<string> filenames = args["align"].as<vector<string> >();
- if (P->t().n_nodes() < 2) return;
+ fs::path dir(dir_name);
- vector<alignment> A;
- for(int i=0;i<P->n_data_partitions();i++)
- A.push_back((*P)[i].A());
+ string base = "C" + convertToString(proc_id+1);
- string base = dir_name + "/" + "C" + convertToString(proc_id+1);
- for(int i=0;i<A.size();i++)
+ int i=0;
+ for(auto& filename: filenames)
{
- checked_ofstream file(base+".P"+convertToString(i+1)+".initial.fasta");
- file<<A[i]<<endl;
+ auto source = fs::path(filename);
+ auto target = fs::path(base+".P"+convertToString(i+1)+".initial.fasta");
+ fs::copy_file(source,dir/target);
+ i++;
}
}
=====================================
src/startup/help.cc
=====================================
--- a/src/startup/help.cc
+++ b/src/startup/help.cc
@@ -308,7 +308,7 @@ optional<string> get_citation(const Rule& rule, bool show_title)
-optional<string> get_citation_doi(const Rule& rule)
+optional<string> get_citation_id(const Rule& rule, const string& idtype)
{
// 1. Check if there is a citation field.
auto citation = rule.get_child_optional("citation");
@@ -320,7 +320,7 @@ optional<string> get_citation_doi(const Rule& rule)
for(auto& identifier: *identifiers)
{
auto type = identifier.second.get_child_optional("type");
- if (not type or type->get_value<string>() != "doi") continue;
+ if (not type or type->get_value<string>() != idtype) continue;
auto id = identifier.second.get_child_optional("id");
if (id)
@@ -353,10 +353,18 @@ optional<string> get_citation_url(const Rule& rule)
}
- if (auto doi = get_citation_doi(rule))
+ if (auto doi = get_citation_id(rule,"doi"))
{
return "https://doi.org/"+*doi;
}
+ else if (auto pmcid = get_citation_id(rule,"pmcid"))
+ {
+ return "https://www.ncbi.nlm.nih.gov/pmc/articles/"+*pmcid;
+ }
+ else if (auto pmid = get_citation_id(rule,"pmid"))
+ {
+ return "https://www.ncbi.nlm.nih.gov/pubmed/"+*pmid;
+ }
return boost::none;
}
@@ -390,7 +398,8 @@ string get_help_for_rule(const Rules& rules, const Rule& rule)
}
help<<header("Usage");
help<<" "<<bold(name);
- if (args_names_types.size()) help<<"["<<join(args_names_types,", ")<<"]\n\n";
+ if (args_names_types.size()) help<<"["<<join(args_names_types,", ")<<"]";
+ help<<"\n\n";
if (auto synonyms = rule.get_child_optional("synonyms"))
{
@@ -401,7 +410,8 @@ string get_help_for_rule(const Rules& rules, const Rule& rule)
help<<indent_and_wrap(3, terminal_width(),join(syn,", "))<<"\n\n";
}
- help<<header("Arguments");
+ if (not args.empty())
+ help<<header("Arguments");
for(auto& argpair: args)
{
auto& arg = argpair.second;
=====================================
src/substitution/substitution.cc
=====================================
--- a/src/substitution/substitution.cc
+++ b/src/substitution/substitution.cc
@@ -527,7 +527,8 @@ namespace substitution {
log_double_t calc_root_probability_SEV(const Likelihood_Cache_Branch* LCB1,
const Likelihood_Cache_Branch* LCB2,
const Likelihood_Cache_Branch* LCB3,
- const Matrix& F)
+ const Matrix& F,
+ const vector<int>& counts)
{
total_calc_root_prob++;
@@ -557,6 +558,7 @@ namespace substitution {
assert(L > 0);
assert(L == bits2.size());
assert(L == bits3.size());
+ assert(L == counts.size());
total_root_clv_length += L;
@@ -630,7 +632,7 @@ namespace substitution {
if (non_gap2) i2++;
if (non_gap3) i3++;
- total.mult_with_count(p_col,1);
+ total.mult_with_count(p_col,counts[c]);
// std::clog<<" i = "<<i<<" p = "<<p_col<<" total = "<<total<<"\n";
}
=====================================
src/tools/alignment-draw.cc
=====================================
--- a/src/tools/alignment-draw.cc
+++ b/src/tools/alignment-draw.cc
@@ -184,6 +184,8 @@ public:
{
if (x == 0) return white;
+ if (n_blocks == 1) return HSV(0,1,1);
+
x = (x-1)/(n_blocks-1);
// otherwise use the rainbow
=====================================
src/tools/alignments-diff.cc
=====================================
--- a/src/tools/alignments-diff.cc
+++ b/src/tools/alignments-diff.cc
@@ -89,7 +89,7 @@ variables_map parse_cmd_line(int argc,char* argv[])
("fill",value<string>()->default_value("gap"),"blank columns filled with: gap or unknown")
("differences-file,d",value<string>(),"Filename to store differences in AU format")
("blocksize",value<int>()->default_value(20),"Width of blocks of same color")
- ("classes",value<int>()->default_value(3),"Number of groups for non-matching characters")
+ ("classes",value<int>()->default_value(1),"Number of groups for non-matching characters")
;
options_description all("All options");
=====================================
src/tools/tree-tool.cc
=====================================
--- /dev/null
+++ b/src/tools/tree-tool.cc
@@ -0,0 +1,152 @@
+/*
+ Copyright (C) 2005-2009 Benjamin Redelings
+
+ This file is part of BAli-Phy.
+
+ BAli-Phy is free software; you can redistribute it and/or modify it under
+ the terms of the GNU General Public License as published by the Free
+ Software Foundation; either version 2, or (at your option) any later
+ version.
+
+ BAli-Phy is distributed in the hope that it will be useful, but WITHOUT ANY
+ WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with BAli-Phy; see the file COPYING. If not see
+ <http://www.gnu.org/licenses/>. */
+
+#include <iostream>
+#include <list>
+#include <utility>
+#include "tree/tree.H"
+#include "tree/sequencetree.H"
+#include "tree/tree-util.H"
+#include "tree-dist.H"
+#include "myexception.H"
+#include "util/assert.hh"
+
+#include <boost/program_options.hpp>
+
+using boost::dynamic_bitset;
+
+namespace po = boost::program_options;
+using po::variables_map;
+
+using std::cout;
+using std::cerr;
+using std::endl;
+using std::string;
+using std::vector;
+using std::list;
+using std::valarray;
+using std::pair;
+
+using boost::optional;
+
+variables_map parse_cmd_line(int argc,char* argv[])
+{
+ using namespace po;
+
+ // named options
+ options_description invisible("Invisible options");
+ invisible.add_options()
+ ("tree", value<string>(),"tree to operate on");
+
+ options_description general("General options");
+ general.add_options()
+ ("help,h", "produce help message")
+ ("verbose,v","Output more log messages on stderr.")
+ ;
+
+ options_description commands("Commands");
+ commands.add_options()
+ ("prune",value<string>(),"Comma-separated taxa to remove")
+ ("resolve","Comma-separated taxa to remove")
+ ;
+ options_description visible("All options");
+ visible.add(general).add(commands);
+ options_description all("All options");
+ all.add(general).add(commands).add(invisible);
+
+ // positional options
+ positional_options_description p;
+ p.add("tree", 1);
+
+ variables_map args;
+ store(command_line_parser(argc, argv).
+ options(all).positional(p).run(), args);
+ notify(args);
+
+ if (args.count("help")) {
+ cout<<"Perform various operations on Newick trees.\n\n";
+ cout<<"Usage: tree-tool <tree-file> [OPTIONS]\n\n";
+ cout<<general<<"\n";
+ cout<<commands<<"\n";
+ exit(0);
+ }
+
+ if (args.count("verbose")) log_verbose = 1;
+
+ return args;
+}
+
+void resolve(Tree& T, int node)
+{
+ while(true)
+ {
+ auto neighbors = T.neighbors(node);
+ if (neighbors.size() <= 3) return;
+
+ int new_node = T.create_node_on_branch(T.directed_branch(neighbors[0],node)).name();
+ if (not T.directed_branch(new_node,node).has_length())
+ T.directed_branch(new_node,node).set_length(0.0);
+ if (not T.directed_branch(new_node,neighbors[0]).has_length())
+ T.directed_branch(new_node,neighbors[0]).set_length(0.0);
+
+ T.reconnect_branch(neighbors[1],node,new_node);
+ assert(T.node(node).degree() < neighbors.size());
+ }
+}
+
+vector<int> polytomies(const Tree& T)
+{
+ vector<int> p;
+ for(int n=0;n<T.n_nodes();n++)
+ if (T.node(n).degree() > 3)
+ p.push_back(n);
+ return p;
+}
+
+
+int main(int argc,char* argv[])
+{
+ try {
+ //----------- Parse command line ----------//
+ variables_map args = parse_cmd_line(argc,argv);
+
+ vector<string> prune;
+ if (args.count("prune")) {
+ string p = args["prune"].as<string>();
+ prune = split(p,',');
+ }
+
+
+ //----------- Read the topology -----------//
+ SequenceTree T = load_T(args);
+
+ if (args.count("resolve"))
+ {
+ for(int n: polytomies(T))
+ resolve(T,n);
+ assert(polytomies(T).empty());
+ }
+ std::cout<<T<<std::endl;
+ }
+ catch (std::exception& e) {
+ std::cerr<<"tree-tool: Error! "<<e.what()<<endl;
+ exit(1);
+ }
+ return 0;
+}
View it on GitLab: https://salsa.debian.org/med-team/bali-phy/commit/e73664112c1ff560964d8b60a4ba295719eebc05
--
View it on GitLab: https://salsa.debian.org/med-team/bali-phy/commit/e73664112c1ff560964d8b60a4ba295719eebc05
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/20180613/0b02da21/attachment-0001.html>
More information about the debian-med-commit
mailing list