Skip to content
Commits on Source (3)
Show whitespace changes
Inline Side-by-side
debian/rules
View file @ a83386b6
......@@ -7,6 +7,14 @@ export DEB_BUILD_MAINT_OPTIONS = hardening=+all
export DEB_CFLAGS_MAINT_APPEND = -O3
export DEB_CXXFLAGS_MAINT_APPEND = -O3
BUILDARCH := $(shell dpkg-architecture -qDEB_BUILD_ARCH)
# Compilation with debugging can use > 4G RAM, so disable it on most architectures.
ifneq ($(BUILDARCH),amd64)
DEB_CFLAGS_MAINT_APPEND += -g0
DEB_CXXFLAGS_MAINT_APPEND += -g0
endif
# Use debhelper.
%:
dh $@ --buildsystem=meson
......
doc/README.itex.xml
View file @ a83386b6
<!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&amp;searchon=names&amp;section=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 0 → 100644
View file @ a83386b6
% 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@googlegroups.com>.
functions/Frequencies/uniform.json 0 → 100644
View file @ a83386b6
{
"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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -4,5 +4,5 @@
"result_type": "DNA",
"call": "Alphabet.dna",
"args": [],
"title": "The RNA alphabet"
"title": "The DNA alphabet"
}
meson.build
View file @ a83386b6
project('bali-phy', ['cpp','c'],
version: '3.1.4',
version: '3.1.5',
default_options : [
'buildtype=release',
'cpp_std=c++14'
......
modules/SModel.hs
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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",@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",@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 &amp; Model</a></h2>\n";
$section .= "<h2 style=\"clear:both\"><a class=\"anchor\" name=\"data\"></a>Data &amp; Model</h2>\n";
$section .= "<table class=\"backlit2\">\n";
$section .= "<tr><th>Partition</th><th>Sequences</th><th>Lengths</th><th>Alphabet</th><th>Substitution&nbsp;Model</th><th>Indel&nbsp;Model</th><th>Scale&nbsp;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,6 +929,19 @@ $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>';
......@@ -923,7 +949,7 @@ $section .= '<img src="c50.SRQ.png" class="r_floating_picture" alt="SRQ plot for
$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
View file @ a83386b6
#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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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
View file @ a83386b6
......@@ -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});
}
}
......