debian/rules

@@ -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

 <!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

+% 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

+{
+    "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

@@ -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

@@ -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

@@ -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

@@ -4,5 +4,5 @@
     "result_type": "DNA",
     "call": "Alphabet.dna",
     "args": [],
-    "title": "The RNA alphabet"
+    "title": "The DNA alphabet"
 }

meson.build

 project('bali-phy', ['cpp','c'],
-	version: '3.1.4',
+	version: '3.1.5',
 	default_options : [
 	  'buildtype=release',
 	  'cpp_std=c++14'

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

@@ -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

+#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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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

@@ -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});
 	}
     }