[med-svn] [Git][med-team/bali-phy][upstream] New upstream version 3.6.1+dfsg
Andreas Tille (@tille)
gitlab at salsa.debian.org
Mon Aug 30 19:29:03 BST 2021
Andreas Tille pushed to branch upstream at Debian Med / bali-phy
Commits:
60599e45 by Andreas Tille at 2021-08-30T17:50:12+02:00
New upstream version 3.6.1+dfsg
- - - - -
27 changed files:
- + .vscode/c_cpp_properties.json
- + .vscode/launch.json
- README.md
- bindings/models/m2a.json
- bindings/models/m2a_test.json
- bindings/models/m3_test.json
- bindings/models/m8.json
- bindings/models/m8a.json
- bindings/models/m8a_test.json
- doc/README.itex.xml
- haskell/Probability/Distribution/Mixture.hs
- meson.build
- src/alignment/load.H
- src/alignment/load.cc
- src/bali-phy/A-T-model.cc
- src/bali-phy/bali-phy.cc
- src/builtins/PopGen.cc
- src/computation/expression/expression_ref.H
- src/computation/optimization/float-out.cc
- src/dp/alignment-sums.H
- src/math/exponential.cc
- src/math/logprod.H
- src/probability/choose.cc
- src/sequence/sequence-format.H
- src/sequence/sequence-format.cc
- src/tools/pickout.cc
- src/util/include/util/math/log-double.H
Changes:
=====================================
.vscode/c_cpp_properties.json
=====================================
@@ -0,0 +1,14 @@
+{
+ "configurations": [
+ {
+ "includePath": [
+ "${workspaceFolder}/**"
+ ],
+ "defines": [],
+ "cStandard": "c17",
+ "cppStandard": "c++17",
+ "compileCommands": "${workspaceFolder}/builddir/compile_commands.json"
+ }
+ ],
+ "version": 4
+}
\ No newline at end of file
=====================================
.vscode/launch.json
=====================================
@@ -0,0 +1,27 @@
+{
+ // Use IntelliSense to learn about possible attributes.
+ // Hover to view descriptions of existing attributes.
+ // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+ "version": "0.2.0",
+ "configurations": [
+ {
+ "name": "(gdb) Launch",
+ "type": "cppdbg",
+ "request": "launch",
+ "program": "${workspaceFolder}/builddir/src/bali-phy/bali-phy",
+ "args": ["${workspaceFolder}/examples/sequences/5S-rRNA/5d-muscle.fasta","--iter=10","--package-path=${workspaceFolder}/builddir/src/builtins:${workspaceFolder}"],
+ "stopAtEntry": false,
+ "cwd": "${workspaceFolder}/run/",
+ "environment": [],
+ "externalConsole": false,
+ "MIMode": "gdb",
+ "setupCommands": [
+ {
+ "description": "Enable pretty-printing for gdb",
+ "text": "-enable-pretty-printing",
+ "ignoreFailures": true
+ }
+ ]
+ }
+ ]
+}
\ No newline at end of file
=====================================
README.md
=====================================
@@ -4,19 +4,7 @@
Install
-------
-Please visit the [releases page](http://www.bali-phy.org/download.php) to download official binaries.
-
-You can also install via homebrew on a Mac, and using `apt-get` on recent version of Debian or Ubuntu.
-
-Documentation
-------------
-
-* [http://bali-phy.org/](http://bali-phy.org/)
-* [Manual](http://bali-phy.org/README.xhtml)
-* [Tutorial](http://bali-phy.org/Tutorial3.html)
-* [Developer's Guide](http://bali-phy.org/developer.html)
-
-The Manual describes [how to install](http://bali-phy.org/README.xhtml#installation) bali-phy in detail. Simplified instructions are below.
+If you just want to install bali-phy, please visit the [release page](http://www.bali-phy.org/download.php). If you want to compile BAli-phy from source, the quick start instructions are below.
Compiling
---------
@@ -26,26 +14,20 @@ You will need a C++ compiler that understands C++17.
* clang 8 (or higher) works
* XCode 11 (or higher) works
-You will also need to install
- * cairo graphics library (optional, but required to build the `draw-tree` program)
-
-You will also need
- * python3
- * ninja
- * meson >= 0.53
-
Install Prerequisites
---------------------
+On Ubuntu, you can use apt-get:
```bash
-sudo apt-get install g++ libcairo2-dev ninja-build python3
+sudo apt-get install g++ libcairo2-dev meson
```
-You also need to install meson. First try:
+
+On Mac (or Linux, actually) you can use homebrew:
```bash
-sudo apt-get install meson
+brew install cairo meson
```
-If the version of meson is not at least 0.49, then you need to install
-meson through pip:
+If the version of meson is not at least 0.53, then you need to install
+meson through the python package manager "pip":
```bash
python3 -m venv meson
source meson/bin/activate
@@ -57,7 +39,7 @@ Build BAli-Phy
```
git clone https://github.com/bredelings/BAli-Phy.git
cd BAli-Phy
-meson build --prefix=$HOME/Applications/bali-phy --buildtype=release
+meson build --prefix=$HOME/Applications/bali-phy
ninja -C build install
ninja -C build test
```
@@ -79,3 +61,12 @@ If you installed in `$HOME/Applications/bali-phy/` as recommended above, then fi
| ~/Applications/bali-phy/share/doc/bali-phy/ | Documentation. |
+Further Documentation
+---------------------
+
+* [http://bali-phy.org/](http://bali-phy.org/)
+* [Manual](http://bali-phy.org/README.xhtml)
+* [Tutorial](http://bali-phy.org/Tutorial4.html)
+
+The Manual describes [how to install](http://bali-phy.org/README.xhtml#installation) bali-phy in greater detail.
+
=====================================
bindings/models/m2a.json
=====================================
@@ -21,7 +21,7 @@
"arg_name": "p0",
"arg_type": "Double",
"default_value": "~uniform[0,1]",
- "description": "The fraction of conserved sites among non-positively selected sites."
+ "description": "The fraction of conserved sites among non-positively selected sites"
},
{
"arg_name": "posP",
@@ -36,7 +36,7 @@
"description": "The dN\/dS value for positively selected sites"
}
],
- "description":"A mixture of conserved, neutral, and positively-selected sites. The conserved sites dN/dS=omega0 and the positively-selected sites have dN/dS=posW. In this parameterization of the M2a model, p0 is the probability of being conserved conditional on not being positively selected.\n\nThe M2a model modifies the M2 model to avoid forcing the conserved dN/dS to 0.",
+ "description":"A mixture of conserved, neutral, and positively-selected sites. The conserved sites have dN/dS=omega0 and the positively-selected sites have dN/dS=posW. In this parameterization of the M2a model, p0 is the probability of being conserved conditional on not being positively selected.\n\nThe M2a model modifies the M2 model to avoid forcing the conserved dN/dS to 0.",
"citation":{"type": "article",
"title": "Accuracy and Power of Statistical Methods for Detecting Adaptive Evolution in Protein Coding Sequences and for Identifying Positively Selected Sites",
"year": "2004",
=====================================
bindings/models/m2a_test.json
=====================================
@@ -58,5 +58,6 @@
"id":"10.1534/genetics.104.031153 "}]
},
"examples": ["function[w,gy94[omega=w,pi=f1x4]]+m2a_test","function[w,mg94[omega=w]]+m2a_test","function[w,fMutSel0[omega=w]]+m2a_test"],
+ "see": ["m1a","m2a"],
"extract": "all"
}
=====================================
bindings/models/m3_test.json
=====================================
@@ -39,7 +39,7 @@
"arg_name": "posSelection",
"arg_type": "Int",
"default_value": "~bernoulli[0.5]",
- "description": "The model selector: 1 if positive selection, 0 if not."
+ "description": "The model selector: 1 if positive selection, 0 if not"
},
{
"arg_name": "n",
=====================================
bindings/models/m8.json
=====================================
@@ -33,7 +33,7 @@
"arg_name": "posP",
"arg_type": "Double",
"default_value": "~beta[1,10]",
- "description": "The fraction of neutral sites"
+ "description": "The fraction of positively selected sites"
},
{
"arg_name": "posW",
@@ -42,7 +42,8 @@
"description": "The dN\/dS value for positively selected sites"
}
],
- "description": "A site-heterogeneous model with a discrete Beta distribution on dN/dS values, plus a category of positively-selected sites. See the m7 model for the definition of mu and gamma.",
+ "see": ["m7", "m8a","m8a_test"],
+ "description": "The m8 model is an m7 model with an additional category of positively-selected sites. The additional category of sites has proportion=posP and dN/dS=posW.",
"citation": {"type": "article",
"title": "Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites",
"year": "2000",
=====================================
bindings/models/m8a.json
=====================================
@@ -33,10 +33,11 @@
"arg_name": "posP",
"arg_type": "Double",
"default_value": "~beta[1,10]",
- "description": "The fraction of neutral sites"
+ "description": "The fraction of *neutral* sites"
}
],
- "description": "A site-heterogeneous model with a discrete Beta distribution on dN/dS values, plus a category of neutral sites. See the m7 model for the definition of mu and gamma.",
+ "description": "The m8a model is an m7 model with an additional category of neutral sites. The additional category of sites has proportion posP.\n\nThe m8a model is also equivalent to m8[posW=1]. That is, you get m8a by constraining posW to 1 in the m8 model. Therefore, the m8a model is nested within the m8 model.\n\nThe m8a model was created as an improved null hypothesis for use in a likelihood ratio test. (See m8a_test)",
+ "see": ["m7","m8","m8a_test"],
"citation":{"type": "article",
"title": "Pervasive Adaptive Evolution in Mammalian Fertilization Proteins",
"year": "2003",
=====================================
bindings/models/m8a_test.json
=====================================
@@ -1,5 +1,6 @@
{
"name": "m8a_test",
+ "title": "Bayesian version of the M8a test for positive selection",
"synonyms": ["M8a_Test"],
"result_type": "MixtureModel[Codons[a]]",
"call": "SModel.m8a_test[@mu, at gamma, at n, at posP, at posW, at posSelection, at submodel]",
@@ -8,7 +9,8 @@
{
"arg_name": "n",
"arg_type": "Int",
- "default_value": "4"
+ "default_value": "4",
+ "description": "The number of bins for discretizing the Beta distribution"
},
{
"arg_name": "submodel",
@@ -19,28 +21,44 @@
{
"arg_name": "mu",
"arg_type": "Double",
- "default_value": "~uniform[0,1]"
+ "default_value": "~uniform[0,1]",
+ "description": "The mean of the Beta distribution"
},
{
"arg_name": "gamma",
"arg_type": "Double",
- "default_value": "~beta[1,10]"
+ "default_value": "~beta[1,10]",
+ "description": "The fraction of possible Beta variance given mu"
},
{
"arg_name": "posP",
"arg_type": "Double",
- "default_value": "~beta[1,10]"
+ "default_value": "~beta[1,10]",
+ "description": "The fraction of positively-selected or neutral sites"
},
{
"arg_name": "posW",
"arg_type": "Double",
- "default_value": "~log_gamma[4,0.25]"
+ "default_value": "~log_gamma[4,0.25]",
+ "description": "The dN\/dS value for positively-selected or neutral sites"
},
{
"arg_name": "posSelection",
"arg_type": "Int",
- "default_value": "~bernoulli[0.5]"
+ "default_value": "~bernoulli[0.5]",
+ "description": "The model indicator -- H0 or H1"
}
],
+ "description": "This is a Bayesian version of the M8a test for positive selection. The original M8a test is a likelihood-ratio test. The Bayesian version of the test chooses a prior with 50% mass on the null hypothesis H0 and 50% on the alternative hypothesis H1:\n\n H0: no positive selection\n H1: some positive selection.\n\nWhen posSelection = 0, the value of posW is ignored and dN/dS=1. When posSelection = 1, the value of posW is used and dN/dS=posW. The posterior probability of posSelection gives the posterior probability of H1. The Bayes Factor is given by Pr(posSelection=1)/Pr(posSelection=0).\n\nThe M8a test (Swanson et al, 2003) has null and alternative hypotheses given by:\n\n H0: m8[posW=1] = m8a\n H1: m8\n\nThis is an improvement over the test originally proposed by Yang (2000):\n\n H0: m8[posP=1] = m7\n H1: m8\n\nIn the originally proposed likelihood-ratio test, a mixture of conserved and neutral sites would not fit the null hypothesis (m7), leading to the erroneous inference of positive selection. In the M8a test, the null hypothesis (m8a) has been improved to handle both neutral and conserved sites, so that positive selection is not incorrectly inferred.",
+ "see": ["m7","m8a","m8a"],
+ "citation":{"type": "article",
+ "title": "Pervasive Adaptive Evolution in Mammalian Fertilization Proteins",
+ "year": "2003",
+ "author": [{"name": "Swanson, Willie J."},
+ {"name": "Nielsen, Rasmus"},
+ {"name": "Yang, Qiaofeng"}],
+ "journal": {"name": "Molecular Biology and Evolution", "volume": "20", "number": "1", "pages": "18--20"},
+ "identifier": [{"type":"doi","id":"10.1093/oxfordjournals.molbev.a004233"}]
+ },
"extract": "all"
}
=====================================
doc/README.itex.xml
=====================================
@@ -106,7 +106,7 @@ C:\cygwin64\home\<emphasis>username</emphasis>\Applications\
<section><info><title>Install on Mac OS X</title></info>
<section><info><title>Install BAli-Phy using homebrew (recommended) </title></info>
- <para>First install the <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://developer.apple.com/xcode/">XCode</link> (version 6 or higher) command line tools:
+ <para>First install the <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://developer.apple.com/xcode/">XCode</link> (version 11 or higher) command line tools:
<screen><prompt>%</prompt> <userinput>xcode-select --install</userinput></screen>
Then install <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://brew.sh/">homebrew</link> and use homebrew to compile and install <command>bali-phy</command>:
@@ -1148,7 +1148,7 @@ The majority consensus tree is usually not fully resolved, so we recommend using
<info><title>Frequencies</title></info>
<para>Frequencies are estimated by default. Frequencies can be fixed by setting the <userinput>pi</userinput> parameter to a constant value, if the model allows unequal frequencies.</para>
<para>Constant frequencies are specified as a list of pairs that associates each letter with its frequency:</para>
- <programlisting language="java">gtr[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]</programlisting>
+ <programlisting language="java">gtr[pi=List[Tuple["A",0.1],Tuple["C",0.2],Tuple["T",0.3],Tuple["G",0.4]]</programlisting>
<para>Frequencies can also be specified using functions:</para>
<para>
@@ -1283,10 +1283,10 @@ The majority consensus tree is usually not fully resolved, so we recommend using
<info><title>Frequencies</title></info>
<para>Frequencies are estimated by default. Frequencies can be fixed by setting the <userinput>pi</userinput> parameter to a constant value, if the model allows unequal frequencies.</para>
<para>Constant frequencies are specified as a list of pairs that associates each letter with its frequency:</para>
- <programlisting language="java">wag+f[pi=List[Pair["A",0.047],Pair["R",0.19],...]]</programlisting>
+ <programlisting language="java">wag+f[pi=List[Tuple["A",0.047],Tuple["R",0.19],...]]</programlisting>
<para>Frequencies can also be specified using functions:</para>
<para>
- <programlisting language="java">wag+f[pi=Frequencies.uniform]]</programlisting>
+ <programlisting language="java">wag+f[pi=Frequencies.uniform]</programlisting>
</para>
<informaltable>
@@ -1421,9 +1421,9 @@ The majority consensus tree is usually not fully resolved, so we recommend using
<info><title>Frequencies</title></info>
<para>Frequencies are estimated by default. Frequencies can be fixed by setting the <userinput>pi</userinput> parameter to a constant value, if the model allows unequal frequencies.</para>
<para>Constant frequencies are specified as a list of pairs that associates each letter with its frequency.</para>
- <programlisting language="java">hky85[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]]+x2
+ <programlisting language="java">hky85[pi=List[Tuple["A",0.1],Tuple["C",0.2],Tuple["T",0.3],Tuple["G",0.4]]]+x2
-hky85_sym+x2_sym+f[pi=List[Pair["AA",0.01],Pair["AC",0.01],Pair["AG",0.01],Pair["AU",0.22],Pair["CA",0.01],Pair["CC",0.01],Pair["CG",0.22],Pair["CU",0.01],Pair["GA",0.01],Pair["GC",0.22],Pair["GG",0.01],Pair["GU",0.01],Pair["UA",0.22],Pair["UC",0.01],Pair["UG",0.01],Pair["UU",0.01]]]</programlisting>
+hky85_sym+x2_sym+f[pi=List[Tuple["AA",0.01],Tuple["AC",0.01],Tuple["AG",0.01],Tuple["AU",0.22],Tuple["CA",0.01],Tuple["CC",0.01],Tuple["CG",0.22],Tuple["CU",0.01],Tuple["GA",0.01],Tuple["GC",0.22],Tuple["GG",0.01],Tuple["GU",0.01],Tuple["UA",0.22],Tuple["UC",0.01],Tuple["UG",0.01],Tuple["UU",0.01]]]</programlisting>
<para>Frequencies can also be specified using functions:</para>
<!-- para>
@@ -1530,7 +1530,7 @@ hky85_sym+x2_sym+f[pi=List[Pair["AA",0.01],Pair["AC",0.01],Pair["AG",0.01],Pair[
<info><title>Frequencies</title></info>
<para>Frequencies are estimated by default. Frequencies can be fixed by setting the <userinput>pi</userinput> parameter to a constant value, if the model allows unequal frequencies.</para>
<para>Constant frequencies are specified as a list of pairs that associates each letter with its frequency.</para>
- <programlisting language="java">hky85[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]]+x3</programlisting>
+ <programlisting language="java">hky85[pi=List[Tuple["A",0.1],Tuple["C",0.2],Tuple["T",0.3],Tuple["G",0.4]]]+x3</programlisting>
<para>Frequencies can also be specified using functions:</para>
<para>
<programlisting language="java">hky85_sym+x3_sym+f[pi=f1x4] // nucleotide frequencies are estimated</programlisting>
@@ -1760,8 +1760,8 @@ o <row>
<info><title>Frequencies</title></info>
<para>Frequencies are estimated by default. Frequencies can be fixed by setting the <userinput>pi</userinput> parameter to a constant value, if the model allows unequal frequencies.</para>
<para>Constant frequencies are specified as a list of pairs that associates each letter with its frequency.</para>
- <programlisting language="java">gy94[pi=List[Pair["AAA",0.01],Pair["C",0.02],...]]
-mg94[pi=List[Pair["A",0.1],Pair["C",0.2],Pair["T",0.3],Pair["G",0.4]]]
+ <programlisting language="java">gy94[pi=List[Tuple["AAA",0.01],Tuple["C",0.02],...]]
+mg94[pi=List[Tuple["A",0.1],Tuple["C",0.2],Tuple["T",0.3],Tuple["G",0.4]]]
</programlisting>
<para>Frequencies can also be specified using functions:</para>
<para>
@@ -2277,7 +2277,7 @@ The <userinput>help</userinput> command can be used to determine the default val
<para>Some types contain parameters. For example <userinput>List[Int]</userinput> indicates a list of integers and <userinput>List[Double]</userinput> indicates a list of real numbers. In order to indicate a list of unknown type, we use a <emphasis>type variable</emphasis> <userinput>a</userinput> and write <userinput>List[a]</userinput>. Type variables always begin with a lower-case letter. They are able to match any specific type, and their value is found by pattern-matching. For example, the function <userinput>add[x,y]</userinput> takes two arguments of type <userinput>a</userinput> and has a result of type <userinput>a</userinput>. Thus:
<programlisting>add[1,2] // arguments are a=Int, so result is of type Int
add[1.0,2.0] // arguments are a=Double, so result is of type Double</programlisting>
-<userinput>Pair[a,b]</userinput> is a parameterized type that can be specialized to (for example) <userinput>Pair[String,Double]</userinput> and <userinput>Pair[Int,Int]</userinput>.
+<userinput>Tuple[a,b]</userinput> is a parameterized type that can be specialized to (for example) <userinput>Tuple[String,Double]</userinput> and <userinput>Tuple[Int,Int]</userinput>.
</para>
<para>Types for components of substitution models are often parameterized by type of the alphabet. For example, hky85 has a result type of <userinput>RevCTMC[a]</userinput>, where <userinput>a</userinput> could be <userinput>DNA</userinput> or <userinput>RNA</userinput>. The use of alphabet types in substitution models prevents combining substitution models with mismatched alphabets.
@@ -2670,7 +2670,7 @@ adequately mix between modes.</para>
<section><info><title>ESS for split frequencies</title></info>
- <para>To calculate the split ESS values, run:
+ <para>As desribed in <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://doi.org/10.3852/10-120">Gaya et al (2011)</link>, we can also compute ESS values for splits on the tree:
<screen><prompt>%</prompt> trees-bootstrap <replaceable>dir-1</replaceable>/C1.trees <replaceable>dir-2</replaceable>/C1.trees ... <replaceable>dir-n</replaceable>/C1.trees > partitions.bs</screen>
To compute the ESS for a split, we consider the presence or absence
of a split in each iteration as a series of binary values. We
@@ -2901,20 +2901,20 @@ The stabilization criterion is the same one described above for numerical values
<para>In order to compile BAli-Phy, you need
<itemizedlist>
- <listitem>a <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://en.wikipedia.org/wiki/C%2B%2B14">C++14</link> compiler</listitem>
- <listitem><link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://mesonbuild.com">meson</link> (version >= 0.45)</listitem>
+ <listitem>a <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://en.wikipedia.org/wiki/C%2B%2B14">C++17</link> compiler</listitem>
+ <listitem><link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://mesonbuild.com">meson</link> (version >= 0.53)</listitem>
</itemizedlist>
- We recommend the GNU C++ Compiler (<link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://gcc.gnu.org">GCC</link>) version 5.0 (or higher) or the <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://clang.llvm.org">Clang</link> compiler version 3.5.0 or higher. The <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://www.cairographics.org/">Cairo</link> graphics library is optional, but if it is missing, the <command>drawtree</command> tool that is used to draw consensus trees won't be built. See also <xref linkend="software_req"/>. </para>
+ We recommend the GNU C++ Compiler (<link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://gcc.gnu.org">GCC</link>) version 9.0 (or higher) or the <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://clang.llvm.org">Clang</link> compiler version 8 or higher. The <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://www.cairographics.org/">Cairo</link> graphics library is optional, but if it is missing, the <command>drawtree</command> tool that is used to draw consensus trees won't be built. See also <xref linkend="software_req"/>. </para>
<section><info><title>Linux</title></info>
<para>On Debian and Ubuntu, you can type:
<screen><prompt>%</prompt> <userinput>sudo apt-get install g++ git libcairo2-dev pandoc</userinput></screen>
</para>
-If your version of Debian or Ubuntu is recent enough to contain meson version 0.45 or higher, you can install meson with apt-get:
+If your version of Debian or Ubuntu is recent enough to contain meson version 0.53 or higher, you can install meson with apt-get:
<screen><prompt>%</prompt> <userinput>sudo apt-get install meson</userinput>
<prompt>%</prompt> <userinput>dpkg -s meson | grep Version</userinput>
-Version: 0.45.1-2
+Version: 0.53
</screen>
Otherwise you can install meson through pip3:
<screen><prompt>%</prompt> <userinput>sudo apt-get install python3 python3-pip ninja</userinput>
@@ -2926,7 +2926,7 @@ Otherwise you can install meson through pip3:
</section>
<section><info><title>Mac</title></info>
- <para>On Mac OS X, the simplest way to get a compiler is to install <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://developer.apple.com/xcode/">XCode</link> (version 6 or newer) command line tools, which come with <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://clang.llvm.org">clang</link>.
+ <para>On Mac OS X, the simplest way to get a compiler is to install <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="https://developer.apple.com/xcode/">XCode</link> (version 11 or newer) command line tools, which come with <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://clang.llvm.org">clang</link>.
<screen><prompt>%</prompt> <userinput>xcode-select --install</userinput></screen> To get the other tools, first install <link xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="http://brew.sh/">homebrew</link>, and then type:
<screen><prompt>%</prompt> <userinput>brew install git meson cairo pandoc</userinput></screen>
</para>
@@ -2981,8 +2981,8 @@ The command <command>bali-phy</command> and its associated tools should then be
</section>
<section><info><title>Options: compiler and linker flags</title></info>
- <para>You can select the C++ compiler by setting the CXX variable. A useful example of this is to use <command>g++-8</command> on systems where <command>g++</command> invokes a compiler that is too old:
- <screen><prompt>%</prompt> <userinput>CXX=g++-8 meson build --prefix=$HOME/Applications/bali-phy-&version; --buildtype=release</userinput></screen>
+ <para>You can select the C++ compiler by setting the CXX variable. A useful example of this is to use <command>g++-9</command> on systems where <command>g++</command> invokes a compiler that is too old:
+ <screen><prompt>%</prompt> <userinput>CXX=g++-9 meson build --prefix=$HOME/Applications/bali-phy-&version; --buildtype=release</userinput></screen>
You may also set compiler and linker options using the CPPFLAGS, CXXFLAGS, and LDFLAGS variables. For example, you can instruct the compiler to use all the features of your chip, instead of producing generic code that will run anywhere:
<screen><prompt>%</prompt> <userinput>CXXFLAGS="-mtune=native -march=native" meson --prefix=$HOME/Applications/bali-phy-&version;</userinput></screen>
=====================================
haskell/Probability/Distribution/Mixture.hs
=====================================
@@ -5,7 +5,7 @@ import Probability.Distribution.Categorical
mixtureRange _ (d:ds) = distRange
mixture_density (p:ps) (dist:dists) x = (doubleToLogDouble p)*(density dist x) + (mixture_density ps dists x)
-mixture_density [] _ x = (doubleToLogDouble 0.0)
+mixture_density [] _ x = (doubleToLogDouble 1.0)
sample_mixture weights dists = do cat <- categorical weights
dists!!cat
=====================================
meson.build
=====================================
@@ -1,8 +1,10 @@
project('bali-phy', ['cpp','c'],
- version: '3.6.0',
+ version: '3.6.1',
default_options : [
+ 'cpp_std=c++17',
'warning_level=3',
- 'cpp_std=c++17'
+ 'buildtype=release',
+ 'b_ndebug=if-release'
],
meson_version: '>= 0.53',
license: 'GPLv2')
@@ -29,15 +31,32 @@ configure_file(output : 'config.h', configuration : conf_data)
root_inc = include_directories('.','src')
add_project_arguments('-DHAVE_CONFIG_H', language : 'cpp')
+
+# 3. Set compiler flags
+
+# 3.1 Set warning flags
warn_flags = ['-Wno-sign-compare','-Wno-maybe-uninitialized','-Woverloaded-virtual','-Wstrict-aliasing','-Wno-unknown-pragmas', '-fdiagnostics-show-template-tree']
add_project_arguments(cpp.get_supported_arguments(warn_flags), language : 'cpp')
+# 3.2 Set optimization flags
+
if get_option('optimization') == '2' or get_option('optimization') == '3'
opt_flags = ['-funroll-loops','-fno-math-errno','-fno-signed-zeros']
- add_project_arguments(['-DNDEBUG','-DNDEBUG_DP']+cpp.get_supported_arguments(opt_flags), language : 'cpp')
+ add_project_arguments(cpp.get_supported_arguments(opt_flags), language : 'cpp')
+endif
+
+# 3.3 Enable/Disable extra debugging for dynamic programming
+assertions_off = get_option('b_ndebug') == 'true' or (get_option('b_ndebug') == 'if-release' and
+ (get_option('buildtype') == 'release' or
+ get_option('buildtype') == 'plain')
+ )
+if assertions_off
+ add_project_arguments('-DNDEBUG_DP', language : 'cpp')
endif
+# 4. Dependencies
+
threads = dependency('threads')
if host_machine.system() == 'darwin'
@@ -52,18 +71,6 @@ else
used_cairo = 'system'
endif
-install_data('scripts/bp-analyze', install_mode: 'rwxr-xr-x', install_dir: 'bin')
-install_data('scripts/bali-phy-pkg', install_mode: 'rwxr-xr-x', install_dir: 'bin')
-install_data(['scripts/compare-runs.R', 'scripts/compare-runs2.R',
- 'scripts/tree-plot.R', 'scripts/tree-plot-3D.R'],
- install_dir:'lib/bali-phy/libexec')
-
-## we need to build libdl if we are on windows (cygwin or mingw64, and then find the library
-install_subdir('haskell', install_dir: 'lib/bali-phy')
-install_subdir('bindings', install_dir: 'lib/bali-phy')
-install_subdir('help', install_dir: 'lib/bali-phy')
-install_subdir('examples', install_dir: get_option('datadir')/'doc/bali-phy')
-
if get_option('with-mpi')
mpi = dependency('mpi', language='cpp')
endif
@@ -174,6 +181,21 @@ test('bali-phy testsuite',
workdir: meson.source_root()/'tests',
args:[baliphy.full_path(), packagepath])
+# Install non-executable files
+
+install_data('scripts/bp-analyze', install_mode: 'rwxr-xr-x', install_dir: 'bin')
+install_data('scripts/bali-phy-pkg', install_mode: 'rwxr-xr-x', install_dir: 'bin')
+install_data(['scripts/compare-runs.R', 'scripts/compare-runs2.R',
+ 'scripts/tree-plot.R', 'scripts/tree-plot-3D.R'],
+ install_dir:'lib/bali-phy/libexec')
+
+install_subdir('haskell', install_dir: 'lib/bali-phy')
+install_subdir('bindings', install_dir: 'lib/bali-phy')
+install_subdir('help', install_dir: 'lib/bali-phy')
+install_subdir('examples', install_dir: get_option('datadir')/'doc/bali-phy')
+
+# Print a summary of the configuration
+
summary({'host': host_machine.system()
}, section: 'Architecture')
@@ -187,7 +209,12 @@ summary({'BOOST': used_boost,
summary({'prefix': get_option('prefix'),
},section: 'Directories')
+assertions_enabled = 'enabled'
+if assertions_off
+ assertions_enabled = 'disabled'
+endif
+
summary({'optimization': get_option('optimization'),
- 'b_ndebug': get_option('b_ndebug'),
'debug': get_option('debug'),
+ 'assertions': assertions_enabled,
},section: 'Configuration')
=====================================
src/alignment/load.H
=====================================
@@ -54,7 +54,7 @@ class alignment_reader: public file_reader<std::string>
{
std::string current;
public:
- const std::string& read() const {return current;}
+ const std::string& read() const {assert(not done()); return current;}
void next();
alignment_reader(std::istream&);
=====================================
src/alignment/load.cc
=====================================
@@ -424,6 +424,8 @@ alignment find_last_alignment(std::istream& ifile, const string& alph_name)
return A;
}
+// After reading the alignment, we will have
+// `not file` if there is not a blank line.
string read_next_alignment(std::istream& file)
{
string alignment_string;
@@ -434,8 +436,13 @@ string read_next_alignment(std::istream& file)
alignment_string += "\n";
}
- assert(alignment_string.size() > 0 and alignment_string[0] == '>');
- return alignment_string;
+ if (file)
+ {
+ assert(alignment_string.size() > 0 and alignment_string[0] == '>');
+ return alignment_string;
+ }
+ else
+ return {};
}
=====================================
src/bali-phy/A-T-model.cc
=====================================
@@ -871,8 +871,8 @@ owned_ptr<Model> create_A_and_T_model(const Rules& R, variables_map& args, const
P.probability();
//-------- Set the alignments for variable partitions ---------//
for(int i=0;i<P.n_data_partitions();i++)
- if (not unalign and P.get_data_partition(i).has_IModel())
- P.get_data_partition(i).set_alignment(A[i]);
+ if (unalign and P.get_data_partition(i).has_IModel())
+ P.get_data_partition(i).set_alignment(unalign_A(A[i]));
P.probability();
// If the tree has any foreground branch attributes, then set the corresponding branch to foreground, here.
=====================================
src/bali-phy/bali-phy.cc
=====================================
@@ -348,7 +348,7 @@ std::string generate_print_program(const model_t& print, const expression_ref& a
}
if (print.code.is_action())
- program_file<<" result <- run_lazy (print_model alphabet)\n";
+ program_file<<" result <- run_lazy ("<<E<<")\n";
else
program_file<<" let result = "<<E<<"\n";
@@ -446,12 +446,11 @@ int main(int argc,char* argv[])
owned_ptr<Model> M; // Declare early so we can reference from catch block.
try {
-
-#if defined(HAVE_FEENABLEEXCEPT) && !defined(NDEBUG)
+// This should probably be a run-time switch.
+#if defined(HAVE_FEENABLEEXCEPT) && defined(DEBUG_FPE)
feenableexcept(FE_DIVBYZERO|FE_OVERFLOW|FE_INVALID);
- // feclearexcept(FE_DIVBYZERO|FE_OVERFLOW|FE_INVALID);
#endif
-#if defined(HAVE_CLEAREXCEPT) && defined(NDEBUG)
+#if defined(HAVE_CLEAREXCEPT) && !defined(DEBUG_FPE)
feclearexcept(FE_DIVBYZERO|FE_OVERFLOW|FE_INVALID);
#endif
fp_scale::initialize();
=====================================
src/builtins/PopGen.cc
=====================================
@@ -9,6 +9,7 @@
#include "util/string/split.H"
#include "util/string/strip.H"
#include "util/log-level.H"
+#include "math/logprod.H"
using std::string;
using std::optional;
@@ -534,17 +535,10 @@ extern "C" closure builtin_function_ewens_diploid_probability(OperationArgs& Arg
n /= 2;
assert(n == I.size());
- double Pr = 1.0;
- log_double_t Pr2 = 1.0;
+ log_prod_underoverflow Pr;
int n_theta_pow = 0;
for(int i=0,total=0;i<n;i++)
{
- if (Pr < 1.0e-30)
- {
- Pr2 *= Pr;
- Pr = 1.0;
- }
-
int a1 = alleles[2*i].as_int();
int a2 = alleles[2*i+1].as_int();
@@ -568,7 +562,7 @@ extern "C" closure builtin_function_ewens_diploid_probability(OperationArgs& Arg
// Heterozygotes coalesce before outbreeding with probability 0.
if (heterozygote and coalesced)
{
- Pr2 *= 0.0; // This accumulates zeros, in order to penalize them.
+ Pr *= 0.0; // This accumulates zeros, in order to penalize them.
continue;
}
@@ -581,12 +575,9 @@ extern "C" closure builtin_function_ewens_diploid_probability(OperationArgs& Arg
}
}
- Pr2 *= Pr * pow(log_double_t(theta), n_theta_pow);
-
- assert(Pr > 0.0);
- // assert(Pr2 > 0.0);
+ Pr *= pow(log_double_t(theta), n_theta_pow);
- return {Pr2};
+ return expression_ref{Pr};
}
// Pr(I|s) = \sum_t=0^\infty s^t (1-s) (1/2^t)^(L-n) (1-(1/2^t))^n
=====================================
src/computation/expression/expression_ref.H
=====================================
@@ -55,6 +55,8 @@ public:
return (i == E2.as_index_var());
break;
default:
+ if (ptr() == E2.ptr()) return true;
+
return (*ptr() == *E2.ptr());
}
}
=====================================
src/computation/optimization/float-out.cc
=====================================
@@ -198,7 +198,12 @@ float_lets(expression_ref& E, int level)
{
// 1. Var
if (is_var(E))
+ {
+ auto x = E.as_<var>();
+ x = strip_level(x);
+ E = x;
return {};
+ }
// 2. Constant
else if (not E.size())
=====================================
src/dp/alignment-sums.H
=====================================
@@ -29,8 +29,9 @@
#include <vector>
#include <memory>
-#include <memory> // for shared_ptr
-#include <vector> // for vector
+#include <memory> // for std::shared_ptr
+#include <vector> // for std::vector
+#include <optional> // for std::optional
#include "dp-engine.H" // for DPengine
#include "util/math/log-double.H" // for log_double_t
#ifndef NDEBUG
=====================================
src/math/exponential.cc
=====================================
@@ -86,11 +86,11 @@ Matrix expm1(const EigenValues& solution,double t)
E(i,j) = temp;
}
-#ifndef NDEBUG
- for(int i=0;i<E.size1();i++)
- for(int j=0;j<E.size2();j++)
- assert(E(i,j) + ((i==j)?1.0:0.0) >= -1.0e-10);
-#endif
+//#ifndef NDEBUG
+// for(int i=0;i<E.size1();i++)
+// for(int j=0;j<E.size2();j++)
+// assert(E(i,j) + ((i==j)?1.0:0.0) >= -1.0e-10);
+//#endif
return E;
}
@@ -117,11 +117,11 @@ Matrix exp(const EigenValues& solution,double t)
E(i,j) = temp;
}
-#ifndef NDEBUG
- for(int i=0;i<E.size1();i++)
- for(int j=0;j<E.size2();j++)
- assert(E(i,j) >= -1.0e-10);
-#endif
+//#ifndef NDEBUG
+// for(int i=0;i<E.size1();i++)
+// for(int j=0;j<E.size2();j++)
+// assert(E(i,j) >= -1.0e-10);
+//#endif
return E;
}
@@ -185,13 +185,16 @@ Matrix exp(const EigenValues& eigensystem, const vector<double>& pi, const doubl
{
double sum = 0;
for(int j=0;j<E.size2();j++) {
- assert(E(i,j) >= -1.0e-10);
+ // assert(E(i,j) >= -1.0e-10);
+
// Force positive
E(i,j) = std::max(0.0, E(i,j));
sum += E(i,j);
}
+
+ // assert(t > 10 or std::abs(sum - 1.0) < 1.0e-10*E.size1());
+
// Renormalize rows
- assert(t > 10 or std::abs(sum - 1.0) < 1.0e-10*E.size1());
double factor = 1.0/sum;
for(int j=0;j<E.size2();j++)
E(i,j) *= factor;
=====================================
src/math/logprod.H
=====================================
@@ -22,6 +22,7 @@
#include <cmath>
#include "util/assert.hh"
+#include "util/math/log-double.H"
#include <iostream>
struct log_prod_underflow
@@ -95,6 +96,110 @@ struct log_prod_underflow
constexpr log_prod_underflow() {}
};
+struct log_prod_underoverflow
+{
+ log_double_t accum = 1.0;
+ double prod = 1.0;
+
+ log_prod_underoverflow& operator*=(log_double_t x)
+ {
+ accum *= x;
+ return *this;
+ }
+
+ // NOTE: by not using the smallest/largest possible values, we avoid FPE exceptions on underflow / overflow (mostly).
+ static constexpr double max_value = 115792089237316195423570985008687907853269984665640564039457584007913129639936e0;
+ // static constexpr double min_value = std::numeric_limits<double>::min();
+ static constexpr double min_value = 1.0/max_value;
+
+ log_prod_underoverflow& operator*=(double x)
+ {
+ bool was_zero = is_zero();
+
+ double mult = prod * x;
+ if (mult < min_value)
+ {
+ // keep the larger of {prod,x} in the "prod" slot.
+ if (prod < x)
+ {
+ accum *= prod;
+ prod = x;
+ }
+ else
+ accum *= x;
+ }
+ else if (mult > max_value)
+ {
+ // keep the smaller of {prod,x} in the "prod" slot.
+ if (prod > x)
+ {
+ accum *= prod;
+ prod = x;
+ }
+ else
+ accum *= x;
+ }
+ else
+ prod = mult;
+
+
+ assert(not this->is_zero() or x == 0 or was_zero);
+ return *this;
+ }
+
+ log_prod_underoverflow& mult_with_count(double x, int c)
+ {
+ assert(c >= 0);
+ assert(prod >= 0);
+
+ constexpr int count_limit = 4;
+
+ bool was_zero = is_zero();
+
+ if (c >= count_limit)
+ accum *= pow(log_double_t(x),c);
+ else
+ {
+ double f = x;
+ for(int i=1;i<c;i++)
+ f *= x;
+
+ double mult = prod * f;
+ if (mult < min_value)
+ {
+ // keep the larger of {prod,x} in the "prod" slot.
+ if (prod < x)
+ {
+ accum *= prod;
+ prod = f;
+ }
+ else
+ accum *= pow(log_double_t(x),c);
+ }
+ else if (mult > max_value)
+ {
+ if (prod > x)
+ {
+ accum *= prod;
+ prod = f;
+ }
+ else
+ accum *= pow(log_double_t(x),c);
+ }
+ else
+ prod = mult;
+ }
+
+ assert(not this->is_zero() or x == 0 or was_zero);
+ return *this;
+ }
+
+ log_double_t value() const {return accum * prod;}
+ operator log_double_t () const {return value();}
+ bool is_zero() const {return value() <= 0.0;}
+ constexpr log_prod_underoverflow() {}
+};
+
struct log_prod_rescale
{
static constexpr double scale_factor = 115792089237316195423570985008687907853269984665640564039457584007913129639936e0;
@@ -115,6 +220,11 @@ struct log_prod_rescale
prod *= scale_factor;
scale++;
}
+ else if (prod > scale_factor)
+ {
+ prod *= scale_min;
+ scale--;
+ }
assert(not this->is_zero() or x == 0 or was_zero);
return *this;
=====================================
src/probability/choose.cc
=====================================
@@ -26,11 +26,20 @@ using std::vector;
int choose2(log_double_t x, log_double_t y)
{
- std::vector<log_double_t> Pr(2);
- Pr[0] = x;
- Pr[1] = y;
+ auto sum0 = x;
+ auto sum1 = x + y;
- return choose_scratch(Pr);
+ log_double_t r = log_double_t(uniform()) * sum1;
+
+ if (r < sum0)
+ return 0;
+ if (r < sum1)
+ return 1;
+
+ choose_exception<log_double_t> c(0, {x,y});
+ c.prepend(":\n");
+ c.prepend(__PRETTY_FUNCTION__);
+ throw c;
}
template <> choose_exception<log_double_t>::choose_exception(int i, const std::vector<log_double_t>& V)
=====================================
src/sequence/sequence-format.H
=====================================
@@ -1,21 +1,21 @@
/*
- Copyright (C) 2004-2005,2007,2009 Benjamin Redelings
+ Copyright (C) 2004-2005,2007,2009,2012 Benjamin Redelings
-This file is part of BAli-Phy.
+ 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 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.
+ 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/>. */
+ 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/>. */
/**
* @file sequence-format.H
@@ -33,39 +33,39 @@ along with BAli-Phy; see the file COPYING. If not see
/// A namespace containing functions for parsing and writing sequence formats
namespace sequence_format
{
- /// A typedef for functions that reads sequences from a file
- typedef std::vector<sequence> (loader_t)(std::istream& file);
+ /// A typedef for functions that reads sequences from a file
+ typedef std::vector<sequence> (loader_t)(std::istream& file);
- /// Read an alignments letters and names from a file in phylip format
- std::vector<sequence> read_guess(std::istream& file);
+ /// Read an alignments letters and names from a file in phylip format
+ std::vector<sequence> read_guess(std::istream& file);
- /// Read an alignments letters and names from a file in phylip format
- std::vector<sequence> read_phylip(std::istream& file);
+ /// Read an alignments letters and names from a file in phylip format
+ std::vector<sequence> read_phylip(std::istream& file);
- /// Read an alignments letters and names from a file in fasta format
- std::vector<sequence> read_fasta(std::istream& file);
+ /// Read an alignments letters and names from a file in fasta format
+ std::vector<sequence> read_fasta(std::istream& file);
- /// Read an alignments letters and names from a file in fasta format
- std::vector<sequence> read_fasta_entire_file(std::istream& file);
+ /// Read an alignments letters and names from a file in fasta format
+ std::vector<sequence> read_fasta_entire_file(std::istream& file);
- /// A typedef for functions that write sequences to a file
- typedef void (dumper_t)(std::ostream&, const std::vector<sequence>&);
+ /// A typedef for functions that write sequences to a file
+ typedef void (dumper_t)(std::ostream&, const std::vector<sequence>&);
- /// Read an alignments letters and names from a file in phylip format
- void write_phylip(std::ostream& file,const std::vector<sequence>& sequences);
+ /// Read an alignments letters and names from a file in phylip format
+ void write_phylip(std::ostream& file,const std::vector<sequence>& sequences);
- /// Read an alignments letters and names from a file in fasta format
- void write_fasta(std::ostream& file, const std::vector<sequence>& sequences);
+ /// Read an alignments letters and names from a file in fasta format
+ void write_fasta(std::ostream& file, const std::vector<sequence>& sequences);
- /// load a format from a file with name filename
- std::vector<sequence> load_from_file(loader_t loader,const std::string& filename);
+ /// load a format from a file with name filename
+ std::vector<sequence> load_from_file(loader_t loader,const std::string& filename);
- /// load a format from a file with name filename
- std::vector<sequence> load_from_file(const std::string& filename);
+ /// load a format from a file with name filename
+ std::vector<sequence> load_from_file(const std::string& filename);
- /// write a format to a file.
- void write_to_file(loader_t loader,const std::vector<sequence>&,const std::string&);
+ /// write a format to a file.
+ void write_to_file(loader_t loader,const std::vector<sequence>&,const std::string&);
}
#endif
=====================================
src/sequence/sequence-format.cc
=====================================
@@ -110,7 +110,7 @@ namespace sequence_format {
else
continue;
}
-
+
// add the letters in
letters += line;
}
=====================================
src/tools/pickout.cc
=====================================
@@ -17,6 +17,7 @@
along with BAli-Phy; see the file COPYING. If not see
<http://www.gnu.org/licenses/>. */
+#include <optional>
#include <iostream>
#include <string>
#include <vector>
=====================================
src/util/include/util/math/log-double.H
=====================================
@@ -43,7 +43,7 @@ public:
operator double() const {return exp(value);}
constexpr log_double_t& operator=(double y) {
- assert(y >= 0);
+ assert(std::isnan(y) or y >= 0);
if (y == 0)
value = log_0;
else if (y == 1)
View it on GitLab: https://salsa.debian.org/med-team/bali-phy/-/commit/60599e45231d39a47fb4932e98f888290bc93689
--
View it on GitLab: https://salsa.debian.org/med-team/bali-phy/-/commit/60599e45231d39a47fb4932e98f888290bc93689
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/20210830/370c6f51/attachment-0001.htm>
More information about the debian-med-commit
mailing list