Skip to content
Commits on Source (4)
Show whitespace changes
Inline Side-by-side
.travis.yml 0 → 100644
View file @ bfcf260b
sudo: required
cache: apt
language: c
os:
- linux
- osx
compiler:
- gcc
script:
# build and install software into $HOME/local
- ./autogen.sh
- ./configure
- make
- sudo make install
- ./configure --enable-phyrex
- make clean
- make
- ./configure --enable-phytime
- make clean
- make
# put the installed binary into the $PATH
- export PATH=$HOME/local/bin:$PATH
# install testiphy
- cd
- git clone https://gitlab.com/testiphy/testiphy.git
# run testiphy
- cd testiphy
- ./testiphy phyml
\ No newline at end of file
README.md
View file @ bfcf260b
[![Build Status](https://travis-ci.org/stephaneguindon/phyml.svg?branch=master)](https://travis-ci.org/stephaneguindon/phyml)
[![install with bioconda](https://img.shields.io/badge/install%20with-bioconda-brightgreen.svg?style=flat)](http://bioconda.github.io/recipes/phyml/README.html)
### 1. You downloaded a stable release of PhyML from [here](https://github.com/stephaneguindon/phyml/releases)
......
configure.ac
View file @ bfcf260b
# -*- Autoconf -*-
# Process this file with autoconf to produce a configure script.
AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; echo "20190321" | tr -d '\n'"]),[guindon@lirmm.fr])
AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; echo "20190909" | tr -d '\n'"]),[guindon@lirmm.fr])
AM_INIT_AUTOMAKE([foreign])
AC_CONFIG_SRCDIR([src/simu.c],[doc/phyml-manual.tex])
......@@ -46,9 +46,11 @@ AC_FUNC_VPRINTF
AC_CHECK_FUNCS([floor pow rint sqrt strchr strstr])
dnl DEFAULT_VECTOR_FLAG=
dnl DEFAULT_VECTOR_FLAG=-mavx
dnl DEFAULT_VECTOR_FLAG="-mavx -mfma"
dnl DEFAULT_VECTOR_FLAG=-msse3
DEFAULT_VECTOR_FLAG=-march=native
DEFAULT_VECTOR_FLAG="-march=native"
AC_ARG_ENABLE([native],
[AS_HELP_STRING([--disable-native],
[Remove AVX or SSE compilation options.])])
......@@ -64,15 +66,17 @@ dnl LT_INIT
dnl AC_PROG_LIBTOOL
DEFAULT_C_FLAG="-std=c99 -O3 -fomit-frame-pointer -funroll-loops -Wall -Winline ${VECTOR_FLAG}"
dnl default_C_FLAG="-O2 ${VECTOR_FLAG}"
dnl DEFAULT_C_FLAG="-O3 -fPIC"
dnl DEFAULT_C_FLAG="-fno-inline-functions -fno-inline-functions-called-once -fno-optimize-sibling-calls -O2 -fno-omit-frame-pointer -g ${VECTOR_FLAG}"
DEFAULT_C_FLAG="-std=c99 -O3 -fomit-frame-pointer -funroll-loops -Wall -Winline -finline ${VECTOR_FLAG}"
dnl DEFAULT_C_FLAG="-O2 ${VECTOR_FLAG}"
dnl DEFAULT_C_FLAG="${VECTOR_FLAG}"
dnl DEFAULT_C_FLAG="-pthread -O3 -g -ffunction-sections -fdata-sections"
AC_ARG_ENABLE([debug],
[AS_HELP_STRING([--enable-debug],
[Remove optimization options and add debug informations.])])
AS_IF([test "x$enable_debug" = "xyes"],
[CFLAGS="-ansi ${VECTOR_FLAG} -pedantic -Wall -std=c99 -O0 -g -Wuninitialized"],
[CFLAGS="-ansi ${VECTOR_FLAG} -pedantic -Wall -std=c99 -O0 -g -Wuninitialized -Wunused -Wunreachable-code"],
[CFLAGS="${DEFAULT_C_FLAG}"])
AC_ARG_ENABLE([safemode],
......@@ -104,11 +108,11 @@ AS_IF([test "x$enable_win" = "xyes"],AC_DEFINE([WIN32],[1],[WIN32 tag on]))
AM_CONDITIONAL([WANT_WIN], [test "$enable_win" = yes])
AC_ARG_ENABLE([beagle], [AS_HELP_STRING([--enable-beagle], [Compute likelihoods using BEAGLE library.])], [beagle=yes],[beagle=no])
AS_IF([test "x$enable_beagle" = "xyes"],[PKG_CHECK_MODULES(BEAGLE, hmsbeagle-1)])
AS_IF([test "x$enable_beagle" = "xyes"],
[CFLAGS="${CFLAGS} ${BEAGLE_CFLAGS} ${BEAGLE_LIBS} -lstdc++"],
[CFLAGS=${CFLAGS}])
dnl AC_ARG_ENABLE([beagle], [AS_HELP_STRING([--enable-beagle], [Compute likelihoods using BEAGLE library.])], [beagle=yes],[beagle=no])
dnl AS_IF([test "x$enable_beagle" = "xyes"],[PKG_CHECK_MODULES([BEAGLE],[hmsbeagle-1])])
dnl AS_IF([test "x$enable_beagle" = "xyes"],
dnl [CFLAGS="${CFLAGS} ${BEAGLE_CFLAGS} ${BEAGLE_LIBS} -lstdc++"],
dnl [CFLAGS=${CFLAGS}])
AS_IF([test "x$enable_beagle" = "xyes"], AC_DEFINE([BEAGLE],[1],[BEAGLE tag on]))
AM_CONDITIONAL([WANT_BEAGLE], [test "$enable_beagle" = yes])
......
debian/changelog
View file @ bfcf260b
phyml (3:3.3.20190909-1) unstable; urgency=medium
* Team upload.
* New upstream version
* Standards-Version: 4.4.1
* Reported autogenerated files in archive to upstream
(https://github.com/stephaneguindon/phyml/issues/113)
-- Steffen Moeller <moeller@debian.org> Wed, 06 Nov 2019 21:46:04 +0100
phyml (3:3.3.20190321-2) unstable; urgency=medium
* Fix autopkgtest
......
debian/control
View file @ bfcf260b
......@@ -10,7 +10,7 @@ Build-Depends: debhelper-compat (= 12),
# libhmsbeagle-dev [any-amd64 any-i386 arm64 armhf],
libopenmpi-dev,
chrpath
Standards-Version: 4.4.0
Standards-Version: 4.4.1
Vcs-Browser: https://salsa.debian.org/med-team/phyml
Vcs-Git: https://salsa.debian.org/med-team/phyml.git
Homepage: http://www.atgc-montpellier.fr/phyml
......
debian/patches/cross.patch
View file @ bfcf260b
......@@ -3,9 +3,11 @@ Last-Update: Tue, 26 Jun 2018 17:32:08 +0200
Bug-Debian: https://bugs.debian.org/902433
Description: Use MPI_C to enable cross building
--- a/configure.ac
+++ b/configure.ac
@@ -101,7 +101,7 @@ AS_IF([test "x$enable_gprof" = "xyes"],
Index: phyml/configure.ac
===================================================================
--- phyml.orig/configure.ac
+++ phyml/configure.ac
@@ -105,7 +105,7 @@ AS_IF([test "x$enable_gprof" = "xyes"],
AC_ARG_ENABLE([mpi],
[AS_HELP_STRING([--enable-mpi],
[Compile with mpicc instead of gcc.])])
......@@ -14,8 +16,10 @@ Description: Use MPI_C to enable cross building
AS_IF([test "x$enable_mpi" = "xyes"],AC_DEFINE([MPI],[1],[MPI tag on]))
AM_CONDITIONAL([WANT_MPI], [test "$enable_mpi" = yes])
--- a/src/Makefile.am
+++ b/src/Makefile.am
Index: phyml/src/Makefile.am
===================================================================
--- phyml.orig/src/Makefile.am
+++ phyml/src/Makefile.am
@@ -315,7 +315,8 @@ sse.c sse.h\
avx.c avx.h\
ancestral.c ancestral.h\
......
debian/patches/fix_declaration_of_temp.patch
View file @ bfcf260b
......@@ -2,9 +2,11 @@ Author: Andreas Tille <tille@debian.org>
Last-Update: Sat, 24 Jun 2017 15:29:56 +0200
Description: Add missing declaration, fix name of BEAGLE constant
--- a/src/utilities.c
+++ b/src/utilities.c
@@ -2196,6 +2196,9 @@ matrix *K80_dist(calign *data, phydbl g_
Index: phyml/src/utilities.c
===================================================================
--- phyml.orig/src/utilities.c
+++ phyml/src/utilities.c
@@ -2198,6 +2198,9 @@ matrix *K80_dist(calign *data, phydbl g_
}
}
......@@ -14,7 +16,7 @@ Description: Add missing declaration, fix name of BEAGLE constant
for(i=0;i<data->n_otu-1;i++)
for(j=i+1;j<data->n_otu;j++)
@@ -6664,7 +6667,10 @@ void Swap_Partial_Lk(t_edge *a, t_edge *
@@ -6640,7 +6643,10 @@ void Swap_Partial_Lk(t_edge *a, t_edge *
int *buff_p_lk_loc, *buff_patt_id;
phydbl *buff_p_lk_tip;
int *buff_ui;
......@@ -26,8 +28,10 @@ Description: Add missing declaration, fix name of BEAGLE constant
if(side_a == LEFT && side_b == LEFT)
{
--- a/src/beagle_utils.c
+++ b/src/beagle_utils.c
Index: phyml/src/beagle_utils.c
===================================================================
--- phyml.orig/src/beagle_utils.c
+++ phyml/src/beagle_utils.c
@@ -53,7 +53,7 @@ void print_beagle_flags(long inFlags) {
if (inFlags & BEAGLE_FLAG_SCALING_ALWAYS) fprintf(stdout, " SCALING_ALWAYS");
if (inFlags & BEAGLE_FLAG_SCALING_DYNAMIC) fprintf(stdout, " SCALING_DYNAMIC");
......
debian/patches/no-configure-ac-date.patch
View file @ bfcf260b
......@@ -9,7 +9,7 @@ Description: Enable reproducible build based on version number
# -*- Autoconf -*-
# Process this file with autoconf to produce a configure script.
-AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; echo "20190321" | tr -d '\n'"]),[guindon@lirmm.fr])
-AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; echo "20190909" | tr -d '\n'"]),[guindon@lirmm.fr])
+AC_INIT([PhyML],esyscmd([sh -c "echo "3.3." | tr -d '\n' ; dpkg-parsechangelog | egrep '^Version:' | cut -f 2 -d ' ' | tr -d '\n'"]),[guindon@lirmm.fr])
AM_INIT_AUTOMAKE([foreign])
......
debian/patches/no-march-native.patch
View file @ bfcf260b
......@@ -2,14 +2,16 @@ Description: Do not emit -march=native
Bug-Debian: https://bugs.debian.org/866139
Author: Graham Inggs <ginggs@debian.org>
Last-Update: 2017-06-30
--- a/configure.ac
+++ b/configure.ac
Index: phyml/configure.ac
===================================================================
--- phyml.orig/configure.ac
+++ phyml/configure.ac
@@ -58,7 +58,7 @@ AC_CHECK_FUNCS([floor pow rint sqrt strc
dnl DEFAULT_VECTOR_FLAG=
dnl DEFAULT_VECTOR_FLAG=-mavx
dnl DEFAULT_VECTOR_FLAG="-mavx -mfma"
dnl DEFAULT_VECTOR_FLAG=-msse3
-DEFAULT_VECTOR_FLAG=-march=native
-DEFAULT_VECTOR_FLAG="-march=native"
+DEFAULT_VECTOR_FLAG=""
AC_ARG_ENABLE([native],
[AS_HELP_STRING([--disable-native],
[Remove AVX or SSE compilation options.])])
debian/patches/no-stdc++.patch
View file @ bfcf260b
......@@ -2,14 +2,16 @@ Author: Andreas Tille <tille@debian.org>
Last-Update: Thu, 27 Aug 2015 14:31:04 +0200
Description: Do not uselessly link against libstdc++
--- a/configure.ac
+++ b/configure.ac
@@ -114,7 +114,7 @@ AM_CONDITIONAL([WANT_WIN], [test "$enabl
AC_ARG_ENABLE([beagle], [AS_HELP_STRING([--enable-beagle], [Compute likelihoods using BEAGLE library.])], [beagle=yes],[beagle=no])
AS_IF([test "x$enable_beagle" = "xyes"],[PKG_CHECK_MODULES(BEAGLE, hmsbeagle-1)])
AS_IF([test "x$enable_beagle" = "xyes"],
- [CFLAGS="${CFLAGS} ${BEAGLE_CFLAGS} ${BEAGLE_LIBS} -lstdc++"],
+ [CFLAGS="${CFLAGS} ${BEAGLE_CFLAGS} ${BEAGLE_LIBS}"],
[CFLAGS=${CFLAGS}])
Index: phyml/configure.ac
===================================================================
--- phyml.orig/configure.ac
+++ phyml/configure.ac
@@ -121,7 +121,7 @@ AM_CONDITIONAL([WANT_WIN], [test "$enabl
dnl AC_ARG_ENABLE([beagle], [AS_HELP_STRING([--enable-beagle], [Compute likelihoods using BEAGLE library.])], [beagle=yes],[beagle=no])
dnl AS_IF([test "x$enable_beagle" = "xyes"],[PKG_CHECK_MODULES([BEAGLE],[hmsbeagle-1])])
dnl AS_IF([test "x$enable_beagle" = "xyes"],
-dnl [CFLAGS="${CFLAGS} ${BEAGLE_CFLAGS} ${BEAGLE_LIBS} -lstdc++"],
+dnl [CFLAGS="${CFLAGS} ${BEAGLE_CFLAGS} ${BEAGLE_LIBS}"],
dnl [CFLAGS=${CFLAGS}])
AS_IF([test "x$enable_beagle" = "xyes"], AC_DEFINE([BEAGLE],[1],[BEAGLE tag on]))
AM_CONDITIONAL([WANT_BEAGLE], [test "$enable_beagle" = yes])
debian/rules
View file @ bfcf260b
......@@ -61,6 +61,7 @@ override_dh_auto_clean:
mv phyml-manual.pdf doc
# remove temporary bin dir for phyml (see above hack)
rm -rf debian/bin
rm -rf doc/Makefile
override_dh_installdocs:
dh_installdocs
......
doc/Makefile.in
View file @ bfcf260b
......@@ -160,8 +160,6 @@ AUTOCONF = @AUTOCONF@
AUTOHEADER = @AUTOHEADER@
AUTOMAKE = @AUTOMAKE@
AWK = @AWK@
BEAGLE_CFLAGS = @BEAGLE_CFLAGS@
BEAGLE_LIBS = @BEAGLE_LIBS@
CC = @CC@
CCDEPMODE = @CCDEPMODE@
CFLAGS = @CFLAGS@
......@@ -198,9 +196,6 @@ PACKAGE_URL = @PACKAGE_URL@
PACKAGE_VERSION = @PACKAGE_VERSION@
PATH_SEPARATOR = @PATH_SEPARATOR@
PDFLATEX = @PDFLATEX@
PKG_CONFIG = @PKG_CONFIG@
PKG_CONFIG_LIBDIR = @PKG_CONFIG_LIBDIR@
PKG_CONFIG_PATH = @PKG_CONFIG_PATH@
SET_MAKE = @SET_MAKE@
SHELL = @SHELL@
STRIP = @STRIP@
......
doc/phyml-manual.tex
View file @ bfcf260b
......@@ -1635,7 +1635,7 @@ PhyML in batch mode.
\end{itemize}
\subsubsection{{\tt topology} component}\index{XML options!{\tt topology} component}
Options:
Each instance of the \x{topology} component (there should be only one...) has the following options:
\begin{itemize}
\item \x{init.tree="bionj"|"user"|"random"}. Starting tree. Default is \x{bionj}.
\item \x{n.rand.starts="X"}. Number of random starting trees. Default is 5.
......@@ -1659,7 +1659,7 @@ Options:
\end{Verbatim}
\subsubsection{{\tt ratematrices} component}\index{XML options!{\tt ratematrices} component}\label{sec:xmlratematrices}
Options:
Each instance of a \x{ratematrices} component have the following options:
\begin{itemize}
\item \x{model="JC69"|"K80"|"F81"|"F84"|"HKY85"|"TN93"|"GTR"|"custom"} for nucleotide data. The default is \x{"HKY85"}.\\
\x{model="LG"|"WAG"|"JTT"|"MtREV"|"Dayhoff"|"DCMut"|"RtREV"|"CpREV"|"VT"}\\\x{|"Blosum62"|"MtMam"|"MtArt"|"HIVw"|"HIVb"|"customaa"}
......@@ -1680,7 +1680,11 @@ for amino-acid sequences. The default is \x{"LG"}.
given value.
\end{itemize}
The {\tt ratematrices} component has the attribute {\tt optimise.weights=yes/no} (default is {\tt
An instance of a \x{ratematrices} component where a GTR or a custom model of substitutions between
nucleotides is used can have a \x{rr} component associated to it in order to specificy the relative
rates of substitutions (see example below).
Also, the {\tt ratematrices} component has the attribute {\tt optimise.weights=yes/no} (default is {\tt
no}). If {\tt optimise.weights=yes}, then the probabilities (or weights) or each matrix in the
set of matrices defined by this component (see Equation \ref{equ:weights}), will be estimated from the data.
......@@ -1690,6 +1694,7 @@ set of matrices defined by this component (see Equation \ref{equ:weights}), will
<ratematrices id="RM1" optimise.weights="yes">
<instance id="M1" model="custom" model.code="000000"/>
<rr AC="1.2" AG="1.5" AT="2.3" CT="1.2" CG="1.5"/>
<instance id="M2" model="GTR" optimise.rr="yes"/>
<instance id="M3" model="WAG"/>
</ratematrices>
......@@ -1697,7 +1702,7 @@ set of matrices defined by this component (see Equation \ref{equ:weights}), will
\end{Verbatim}
\subsubsection{{\tt equfreqs} component}\index{XML options!{\tt equfreqs} component}
Options:
Each instance of a \x{equfreqs} component has the following options:
\begin{itemize}
\item \x{base.freqs="a,b,c,d"} where \x{a-d} are nucleotide frequencies. Make sure that these
frequencies are separated by comas and no space character is inserted.
......@@ -1744,11 +1749,11 @@ Options:
\end{Verbatim}
\subsubsection{{\tt siterates} component}\index{XML options!{\tt siterates} component}
Options:
Each instance defines a class of relative rate of substitution. It has the following options:
\begin{itemize}
\item \x{value="val"}, where \x{"val"} is the relative substitution rate for the corresponding class.
\end{itemize}
A \x{siterate} component generally includes a \x{weights} element that specifies the probabilitic
A \x{siterates} component generally includes a \x{weights} element that specifies the probabilitic
distribution of the relative rates. The available options for such element are:
\begin{itemize}
\item \x{family="gamma|gamma+inv|freerates"}. \x{gamma} indicates that the distribution of the
......@@ -2420,7 +2425,7 @@ corresponding to \x{cal1}, the other to \x{cal2}). The MRCA of
\x{Gymno\_Araucaria}, \x{Gymno\_Ginko}, \x{Gymno\_Juniperus} and \x{Gymno\_Juniperus} has an age
that falls here in the $[50,60]$ interval.
\subsubsection{MCMC settings}
\subsubsection{MCMC settings}\label{sec:phytimesettings}
PhyTime estimates the joint posterior distribution of the phylogenetic model parameters using an
MCMC algorithm. This algorithm relies on sampling values of these parameters iteratively. The
......@@ -2756,10 +2761,12 @@ statistical phylogeographic framework''. 2014. Systematic Biology.
\subsection{PhyREX}\index{PhyRex} PhyREX is a program that implements the spatial
\subsection{PhyREX}\index{PhyREX} PhyREX is a program that implements the spatial
$\Lambda$-Fleming-Viot model
\cite{etheridge2008,berestycki2009,barton2010,barton2010b,veber2012,barton2013} or \sfv\ for short.
According to this model, the spatial distribution of individuals in a population is uniform and does
This model can be thought of as a structured-coalescent where space is considered as continuous
rather than organized into separate demes.
Under the \sfv\ model, the spatial distribution of individuals in a population is uniform and does
not change during the course of evolution, as opposed to other models such as the very popular
isolation by distance model proposed by Wright and Mal\'ecot. The \sfv\ model does not
suffer from the ``pain in the torus'' \cite{felsenstein1975} and clumps of individuals with
......@@ -2788,11 +2795,14 @@ make;
\subsubsection{Running PhyREX}
Example input files for PhyREX can be found in the \x{examples/phyrex\_input\_files/} directory.
PhyRex takes as input a sequence alignment, in a PhyML-compatible format as well as a text file
giving the spatial coordinates of each taxon in two dimension. A example of the coordinates file is given
Example input files for PhyREX can be found in the \x{examples/phyrex/} directory.
PhyREX takes as input a XML file (See Section \ref{sec:xmlio}), in a PhyML-compatible format. This
file defines the models and various parameter settings for conducting an analysis. We give more
information below about this file. The XML file gives the names of a sequence alignment file in the
standard PHYLIP (or NEXUS) format along with a file providing the spatial coordinates of the
corresponding sequences. An example of the coordinates file is given
below:
\begin{Verbatim}[frame=single, label=Valid PhyRex spatial location file, samepage=true, baselinestretch=0.5]
\begin{Verbatim}[frame=single, label=Valid PhyREX spatial location file, samepage=true, baselinestretch=0.5]
# state.name lon lat
|SouthWest| 0 0
|NorthEast| 10 10
......@@ -2810,27 +2820,95 @@ pfrnpk 4.117791 4.489819
elwdvr 5.649958 4.824092
lptxiv 4.563302 4.005124
\end{Verbatim}
The first row in this file gives the names of the columns. It starts with the `\#' character signaling a
comment. The second and third rows define the limits of the population's habitat. At the
moment, this habitat is assumed to be a rectangle. The position in space of that rectangle are
determined by the coordinates of the bottom-left corner (\x{|SouthWest| 0 0}) and the top-right
one (\x{|NorthEast| 10 10}). The `\x{|}' characters help identify the terms in the list of coordinate
corresponding indeed to these two particular points and should thus not be omitted. The remaining
row give the coordinates of each taxon. Every taxon in the sequence file should also be listed in
the coordinate file.
PhyRex uses the same command-line arguments as that used with PhyML. The only exception being the
mandatory argument corresponding to the coordinate file. A typical command-line will this look as
follows: \x{./phyrex -i sequence\_file --coord\_file=./spatial\_location\_file -c 3 --freerates
--il}. With this particular set of options, the sequence alignment will be analysed under a
FreeRate model \cite{soubrier12} with three rate classes and a covarion-like model of rate variation
across lineaged \cite{guindon13}. Please note that the command-line option \x{--coord\_file} is mandatory.
This first column in this file gives the sequence names, the second is the longitude and the third
column gives the latitude. The first row in this file gives the names of the columns. It starts with
the `\#' character signaling a comment. The second and third rows define the limits of the
population's habitat. At the moment, this habitat is assumed to be a rectangle. The position in
space of that rectangle are determined by the coordinates of the bottom-left corner (\x{|SouthWest|
0 0}) and the top-right one (\x{|NorthEast| 10 10}). The `\x{|}' characters help identify the terms
in the list of coordinate corresponding indeed to these two particular points and should thus not be
omitted. The remaining row give the coordinates of each taxon. Every taxon in the sequence file
should also be listed in the coordinate file.
Below is an excerpt of a XML file that can be used as input to run a PhyREX analysis. The \sfv\
model parameters are estimated using a Bayesian sampler which hyperparameters (number of
samples, frequency at which samples are collected, length of the ``burnin'' phase, etc.) are set
through the corresponding XML attributes (see \ref{sec:phytimesettings}). PhyREX is similar to
PhyTime in the sense that substitution rates and times are separate parameters. The rates can
vary across lineages (using for instance a log-normal model) or not (strict clock). In case all
sequences were collected at (approximately) the same time, it is recommended to set the average
rate of evolution to a specific value (\x{<clockrate value="1E-4"/>}). Yet, with serially-sampled
data, the average rate of evolution can be estimated directly from the data and setting it to
a particular value is thus not mandatory. Finally, it is necessary to give information about the
time at which each sequence was collected. This is done through seeting appropriate calibration
constraints. From a technical perspective, these constraints apply to clades where each clade here
is a single tip.
\begin{Verbatim}[frame=single, label=Example (excerpt) of PhyREX XML file, samepage=true, baselinestretch=0.5,
fontsize=\small, numbers=left]
<phyrex run.id="example" output.file="out"
mcmc.chain.len="1E+7"
mcmc.sample.every="500"
mcmc.print.every="100"
mcmc.burnin="10000"
mutmap="no">
<!-- ``Standard'' stuff goes here, i.e., <topology>, -->
<!-- <partitionelem>, etc -->
<!-- Model of rate variation across lineages -->
<lineagerates model="lognormal"/>
<!-- Set the average (clock) rate of substitution -->
<!-- to a fixed value (optional) -->
<clockrate value="1E-4"/>
<-- File where 2-D spatial coordinates are found -->
<coordinates id="coordinates" file.name="./usa_coord.txt"/>
<!-- Set the date for each tip. Useful in serial sample -->
<!-- analysis -->
<clade id="clad1">
<taxon value="CY130177|South_Carolina|12_13|H1N1"/>
</clade>
<!-- Clade 1 (made of a single sequence) has time set to 0 -->
<calibration id="cal1">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad1"/>
</calibration>
<clade id="clad2">
<taxon value="KF647978|Idaho|12_13|H1N1"/>
</clade>
<!-- Clade 2 has time set to 10 -->
<calibration id="cal2">
<lower>10</lower>
<upper>10</upper>
<appliesto clade.id="clad2"/>
</calibration>
<!-- More calibration info goes here. At least, one -->
<!-- calibrated clade per tip. -->
\end{Verbatim}
PhyREX generates correlated samples from the posterior distribution of the \sfv\ model
parameters. The estimated distributions and related summary statistics can be monitored during the
analysis using the MCMC vizualization software Icylog. Figure \ref{fig:phyrextrace} shows the traces
for two parameters after about an hour of calculation on the example data set (36 taxa, 50
geographic locations).
geographic locations). The names of the sampled parameters are self-explanatory. However, if you
are unsure about the precise definition of one of these parameters or have any questions regarding
how to interpret the results returned by PhyREX, please do not hesitate to contact me (\url{guindon@lirmm.fr}).
\begin{figure}
\begin{center}
......@@ -2843,6 +2921,15 @@ the dispersal intensity (bottom) are shown here.}
\label{fig:phyrextrace}
\end{figure}
\subsubsection{Citing PhyREX}
The ``default citation'' is:
\begin{itemize}
\item Guindon S., H Guo, D Welch. ``Demographic inference under the coalescent in a spatial
continuum'', {\it Theoretical Population Biology}, 2016, 111:43-50.
\end{itemize}
\section{Recommendations on program usage}\label{sec:progusage}
\subsection{PhyML}
......
examples/phyrex/flu.xml 0 → 100644
View file @ bfcf260b
<phyrex run.id="example" output.file="out" mcmc.chain.len="1E+7" mcmc.sample.every="500"
mcmc.print.every="100" mcmc.burnin="10000" mutmap="no">
<!-- Tree topology -->
<topology>
<instance id="T1" init.tree="bionj" optimise.tree="yes"/>
</topology>
<!-- Model of rate variation across lineages -->
<lineagerates model="lognormal"/>
<!-- Average (clock) rate of substitution -->
<clockrate value="1E-4"/>
<!-- Substitution model -->
<ratematrices id="RM1">
<instance id="M1" model="HKY85" optimise.tstv="no" tstv="4.0"/>
</ratematrices>
<!-- Freerate model of variation of rates across sites -->
<siterates id="SR1">
<instance id="R1" init.value="1.0"/>
<instance id="R2" init.value="1.0"/>
<instance id="R3" init.value="1.0"/>
<instance id="R4" init.value="1.0"/>
<weights id="D1" family="gamma" optimise.freerates="no">
<instance appliesto="R1" value="1.00"/>
</weights>
</siterates>
<!-- Nucleotide frequencies -->
<equfreqs id="EF1">
<instance id="F1" optimise.freqs="no"/>
</equfreqs>
<!-- Vector of edge lengths -->
<branchlengths id="BL1" >
<instance id="L1" optimise.lens="yes"/>
</branchlengths>
<!-- Model assembly -->
<partitionelem id="partition1" file.name="./h1n1.nxs" data.type="nt" interleaved="no">
<mixtureelem list="T1,T1,T1,T1"/>
<mixtureelem list="M1,M1,M1,M1"/>
<mixtureelem list="F1,F1,F1,F1"/>
<mixtureelem list="R1,R2,R3,R4"/>
<mixtureelem list="L1,L1,L1,L1"/>
</partitionelem>
<coordinates id="coordinates" file.name="./usa_coord.txt"/>
<clade id="clad1">
<taxon value="CY130177|South_Carolina|12_13|H1N1"/>
</clade>
<clade id="clad2">
<taxon value="KF647978|Idaho|12_13|H1N1"/>
</clade>
<clade id="clad3">
<taxon value="KF648268|New_Mexico|12_13|H1N1"/>
</clade>
<clade id="clad4">
<taxon value="KF648064|North_Carolina|12_13|H1N1"/>
</clade>
<clade id="clad5">
<taxon value="KF647995|New_Hampshire|12_13|H1N1"/>
</clade>
<clade id="clad6">
<taxon value="KF648000|Colorado|12_13|H1N1"/>
</clade>
<clade id="clad7">
<taxon value="KF648133|Florida|12_13|H1N1"/>
</clade>
<clade id="clad8">
<taxon value="KF648145|Georgia|12_13|H1N1"/>
</clade>
<clade id="clad9">
<taxon value="KF648155|Vermont|12_13|H1N1"/>
</clade>
<clade id="clad10">
<taxon value="KF648158|Utah|12_13|H1N1"/>
</clade>
<clade id="clad11">
<taxon value="KF648159|Louisiana|12_13|H1N1"/>
</clade>
<clade id="clad12">
<taxon value="KF648008|Rhode_Island|12_13|H1N1"/>
</clade>
<clade id="clad13">
<taxon value="KF648011|New_York|12_13|H1N1"/>
</clade>
<clade id="clad14">
<taxon value="KF647921|Maryland|12_13|H1N1"/>
</clade>
<clade id="clad15">
<taxon value="KF647926|Missouri|12_13|H1N1"/>
</clade>
<clade id="clad16">
<taxon value="KF647928|Texas|12_13|H1N1"/>
</clade>
<clade id="clad17">
<taxon value="KF647936|Wyoming|12_13|H1N1"/>
</clade>
<clade id="clad18">
<taxon value="KF648197|Indiana|12_13|H1N1"/>
</clade>
<clade id="clad19">
<taxon value="KF648228|Mississippi|12_13|H1N1"/>
</clade>
<clade id="clad20">
<taxon value="KF648239|Tennessee|12_13|H1N1"/>
</clade>
<clade id="clad21">
<taxon value="KF648034|Arizona|12_13|H1N1"/>
</clade>
<clade id="clad22">
<taxon value="KF648037|Kentucky|12_13|H1N1"/>
</clade>
<clade id="clad23">
<taxon value="KF648038|Kansas|12_13|H1N1"/>
</clade>
<clade id="clad24">
<taxon value="KF647919|Alabama|12_13|H1N1"/>
</clade>
<clade id="clad25">
<taxon value="KF647931|Wisconsin|12_13|H1N1"/>
</clade>
<clade id="clad26">
<taxon value="KF647932|Pennsylvania|12_13|H1N1"/>
</clade>
<clade id="clad27">
<taxon value="KF647958|Nebraska|12_13|H1N1"/>
</clade>
<clade id="clad28">
<taxon value="KF647971|Oklahoma|12_13|H1N1"/>
</clade>
<clade id="clad29">
<taxon value="KF647902|Iowa|12_13|H1N1"/>
</clade>
<clade id="clad30">
<taxon value="KF648002|Minnesota|12_13|H1N1"/>
</clade>
<clade id="clad31">
<taxon value="KF648040|New_Jersey|12_13|H1N1"/>
</clade>
<clade id="clad32">
<taxon value="KF648101|Ohio|12_13|H1N1"/>
</clade>
<clade id="clad33">
<taxon value="KF648260|Washington|12_13|H1N1"/>
</clade>
<clade id="clad34">
<taxon value="CY170697|California|12_13|H1N1"/>
</clade>
<clade id="clad35">
<taxon value="CY168881|Massachusetts|12_13|H1N1"/>
</clade>
<clade id="clad36">
<taxon value="CY171161|Illinois|12_13|H1N1"/>
</clade>
<calibration id="cal1">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad1"/>
</calibration>
<calibration id="cal2">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad2"/>
</calibration>
<calibration id="cal3">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad3"/>
</calibration>
<calibration id="cal4">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad4"/>
</calibration>
<calibration id="cal5">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad5"/>
</calibration>
<calibration id="cal6">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad6"/>
</calibration>
<calibration id="cal7">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad7"/>
</calibration>
<calibration id="cal8">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad8"/>
</calibration>
<calibration id="cal9">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad9"/>
</calibration>
<calibration id="cal10">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad10"/>
</calibration>
<calibration id="cal11">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad11"/>
</calibration>
<calibration id="cal12">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad12"/>
</calibration>
<calibration id="cal13">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad13"/>
</calibration>
<calibration id="cal14">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad14"/>
</calibration>
<calibration id="cal15">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad15"/>
</calibration>
<calibration id="cal16">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad16"/>
</calibration>
<calibration id="cal17">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad17"/>
</calibration>
<calibration id="cal18">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad18"/>
</calibration>
<calibration id="cal19">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad19"/>
</calibration>
<calibration id="cal20">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad20"/>
</calibration>
<calibration id="cal21">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad21"/>
</calibration>
<calibration id="cal22">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad22"/>
</calibration>
<calibration id="cal23">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad23"/>
</calibration>
<calibration id="cal24">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad24"/>
</calibration>
<calibration id="cal25">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad25"/>
</calibration>
<calibration id="cal26">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad26"/>
</calibration>
<calibration id="cal27">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad27"/>
</calibration>
<calibration id="cal28">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad28"/>
</calibration>
<calibration id="cal29">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad29"/>
</calibration>
<calibration id="cal30">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad30"/>
</calibration>
<calibration id="cal31">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad31"/>
</calibration>
<calibration id="cal32">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad32"/>
</calibration>
<calibration id="cal33">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad33"/>
</calibration>
<calibration id="cal34">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad34"/>
</calibration>
<calibration id="cal35">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad35"/>
</calibration>
<calibration id="cal36">
<lower>0</lower>
<upper>0</upper>
<appliesto clade.id="clad36"/>
</calibration>
</phyrex>
examples/phytime/dating_example.xml
View file @ bfcf260b
<phytime run.id="example" output.file="out" mcmc.chain.len="1E+7" mcmc.sample.every="100"
mcmc.print.every="50" mcmc.burnin="10000" mutmap="no">
mcmc.print.every="50" mcmc.burnin="10000" mutmap="no" r.seed="1567773784">
<!-- Tree topology -->
<topology>
......@@ -12,7 +12,7 @@
<!-- Substitution model -->
<ratematrices id="RM1">
<instance id="M1" model="HKY85" optimise.tstv="no" tstv="4.0"/>
<instance id="M1" model="HKY85" optimise.tstv="yes" tstv="4.0"/>
</ratematrices>
......@@ -21,13 +21,14 @@
<instance id="R3" init.value="0.5"/>
<instance id="R2" init.value="1.0"/>
<instance id="R1" init.value="2.0"/>
<weights id="D1" family="freerates" optimise.freerates="no">
<weights id="D1" family="freerates" optimise.freerates="yes">
<instance appliesto="R3" value="0.33"/>
<instance appliesto="R2" value="0.33"/>
<instance appliesto="R1" value="0.33"/>
</weights>
</siterates>
<!-- Nucleotide frequencies -->
<equfreqs id="EF1">
<instance id="F1" optimise.freqs="no"/>
......@@ -64,6 +65,7 @@
<calibration id="cal1">
<lower>40</lower>
<upper>60</upper>
<!-- <appliesto clade.id="Gymno2"/> -->
<appliesto clade.id="Gymno1" probability="0.8"/>
<appliesto clade.id="Gymno2" probability="0.2"/>
</calibration>
......
src/Makefile.in
View file @ bfcf260b
......@@ -836,8 +836,6 @@ AUTOCONF = @AUTOCONF@
AUTOHEADER = @AUTOHEADER@
AUTOMAKE = @AUTOMAKE@
AWK = @AWK@
BEAGLE_CFLAGS = @BEAGLE_CFLAGS@
BEAGLE_LIBS = @BEAGLE_LIBS@
CC = @CC@
CCDEPMODE = @CCDEPMODE@
CFLAGS = @CFLAGS@
......@@ -874,9 +872,6 @@ PACKAGE_URL = @PACKAGE_URL@
PACKAGE_VERSION = @PACKAGE_VERSION@
PATH_SEPARATOR = @PATH_SEPARATOR@
PDFLATEX = @PDFLATEX@
PKG_CONFIG = @PKG_CONFIG@
PKG_CONFIG_LIBDIR = @PKG_CONFIG_LIBDIR@
PKG_CONFIG_PATH = @PKG_CONFIG_PATH@
SET_MAKE = @SET_MAKE@
SHELL = @SHELL@
STRIP = @STRIP@
......
src/alrt.c
View file @ bfcf260b
......@@ -368,9 +368,9 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
v_init = Duplicate_Scalar_Dbl(b_fcus->l_var);
lk_init = tree->c_lnL;
lk_temp = UNLIKELY;
b_fcus->nni->score = .0;
lk0 = lk1 = lk2 = UNLIKELY;
v1 = v2 = v3 = v4 = NULL;
Init_NNI_Score(0.0,b_fcus,tree);
/*! vertices */
v1 = b_fcus->left->v[b_fcus->l_v1];
......@@ -380,13 +380,13 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
if(v1->num < v2->num)
{
PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
Exit("");
}
if(v3->num < v4->num)
{
PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
Exit("");
}
......@@ -417,28 +417,28 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
if(b_fcus->left->v[i] != b_fcus->rght) //Only consider left_1 and left_2
{
Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus->left->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
}
Update_Partial_Lk(tree,b_fcus,b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus,tree);
lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
for(i=0;i<3;i++)
if(b_fcus->rght->v[i] != b_fcus->left) //Only consider right_1 and right_2
{
Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
lk_temp = Br_Len_Opt(b_fcus->rght->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
}
Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
if(lk_temp < lk0 - tree->mod->s_opt->min_diff_lk_local)
{
PhyML_Fprintf(stderr,"\n== lk_temp = %f lk0 = %f\n",lk_temp,lk0);
PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n\n",__FILE__,__LINE__);
PhyML_Fprintf(stderr,"\n. lk_temp = %f lk0 = %f\n",lk_temp,lk0);
PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
Exit("\n");
}
}
while(FABS(lk_temp-lk0) > tree->mod->s_opt->min_diff_lk_global);
while(fabs(lk_temp-lk0) > tree->mod->s_opt->min_diff_lk_global);
//until no significative improvement is detected
lk0 = tree->c_lnL;
......@@ -479,8 +479,8 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
MIXT_Set_Alias_Subpatt(NO,tree);
if((lk_temp > lk_init + tree->mod->s_opt->min_diff_lk_local) || (lk_temp < lk_init - tree->mod->s_opt->min_diff_lk_local))
{
PhyML_Fprintf(stderr,"\n== lk_temp = %f lk_init = %f",lk_temp,lk_init);
PhyML_Fprintf(stderr,"\n== Err. in file %s at line %d\n",__FILE__,__LINE__);
PhyML_Fprintf(stderr,"\n. lk_temp = %f lk_init = %f",lk_temp,lk_init);
PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n",__FILE__,__LINE__);
Exit("\n");
}
......@@ -513,26 +513,26 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
if(b_fcus->left->v[i] != b_fcus->rght)
{
Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus->left->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
}
Update_Partial_Lk(tree,b_fcus,b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus,tree);
lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
for(i=0;i<3;i++)
if(b_fcus->rght->v[i] != b_fcus->left)
{
Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
lk_temp = Br_Len_Opt(b_fcus->rght->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
}
Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
if(lk_temp < lk1 - tree->mod->s_opt->min_diff_lk_local)
{
PhyML_Fprintf(stderr,"\n== lk_temp = %f lk1 = %f\n",lk_temp,lk1);
PhyML_Fprintf(stderr,"\n== Err in file %s at line %d\n\n",__FILE__,__LINE__);
Warn_And_Exit("");
PhyML_Fprintf(stderr,"\n. lk_temp = %f lk1 = %f\n",lk_temp,lk1);
PhyML_Fprintf(stderr,"\n. Err. in file %s at line %d\n\n",__FILE__,__LINE__);
Exit("");
}
}
while(FABS(lk_temp-lk1) > tree->mod->s_opt->min_diff_lk_global);
......@@ -612,7 +612,7 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
if(b_fcus->left->v[i] != b_fcus->rght)
{
Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus->left->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
if(lk_temp < lk2 - tree->mod->s_opt->min_diff_lk_local)
{
......@@ -623,7 +623,7 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
}
}
Update_Partial_Lk(tree,b_fcus,b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus,tree);
lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
if(lk_temp < lk2 - tree->mod->s_opt->min_diff_lk_local)
{
......@@ -637,7 +637,7 @@ int NNI_Neigh_BL(t_edge *b_fcus, t_tree *tree)
if(b_fcus->rght->v[i] != b_fcus->left)
{
Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
lk_temp = Br_Len_Opt(b_fcus->rght->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
if(lk_temp < lk2 - tree->mod->s_opt->min_diff_lk_local)
{
......@@ -844,18 +844,18 @@ void Make_Target_Swap(t_tree *tree, t_edge *b_fcus, int swaptodo)
if(b_fcus->left->v[i] != b_fcus->rght)
{
Update_Partial_Lk(tree,b_fcus->left->b[i],b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus->left->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->left->b[i]->l->v),b_fcus->left->b[i],tree);
}
Update_Partial_Lk(tree,b_fcus,b_fcus->left);
lk_temp = Br_Len_Opt(b_fcus,tree);
lk_temp = Br_Len_Opt(&(b_fcus->l->v),b_fcus,tree);
for(i=0;i<3;i++)
if(b_fcus->rght->v[i] != b_fcus->left)
{
Update_Partial_Lk(tree,b_fcus->rght->b[i],b_fcus->rght);
lk_temp = Br_Len_Opt(b_fcus->rght->b[i],tree);
lk_temp = Br_Len_Opt(&(b_fcus->rght->b[i]->l->v),b_fcus->rght->b[i],tree);
}
Update_Partial_Lk(tree,b_fcus,b_fcus->rght);
......
src/ancestral.c
View file @ bfcf260b
......@@ -570,7 +570,7 @@ void Ancestral_Sequences(t_tree *tree, int print)
PhyML_Fprintf(tree->io->fp_out_ancestral_seq,"\n\n");
PhyML_Fprintf(tree->io->fp_out_ancestral_seq,"Site\tNodeLabel\t");
for(i=0;i<tree->mod->ns;i++) PhyML_Fprintf(tree->io->fp_out_ancestral_seq,"%c\t",Reciproc_Assign_State(i,tree->io->datatype));
for(i=0;i<tree->mod->ns;i++) PhyML_Fprintf(tree->io->fp_out_ancestral_seq,"%10c\t",Reciproc_Assign_State(i,tree->io->datatype));
PhyML_Fprintf(tree->io->fp_out_ancestral_seq,"MPEE\t");
PhyML_Fprintf(tree->io->fp_out_ancestral_seq,"\n");
......@@ -850,7 +850,7 @@ void Ancestral_Sequences_One_Node(t_node *d, t_tree *tree, int print)
sum_probas = .0;
for(i=0;i<ns;i++)
{
PhyML_Fprintf(fp,"%.4f\t",p[i]);
PhyML_Fprintf(fp,"%10g\t",p[i]);
sum_probas += p[i];
}
if(Are_Equal(sum_probas,1.0,0.01) == NO)
......
src/avx.c
View file @ bfcf260b
......@@ -103,30 +103,6 @@ void AVX_Update_Eigen_Lr(t_edge *b, t_tree *tree)
if(b->left->tax == YES) p_lk_left += ns;
if(b->rght->tax == YES) p_lk_rght += ns;
/* p_lk_left -= b->left->tax ? ns : ns*ncatg; */
/* p_lk_rght -= b->rght->tax ? ns : ns*ncatg; */
/* dot_prod -= ns; */
/* PhyML_Printf("\n. EIGEN %d [%g %g %g %g] %p: [%g %g %g %g] %p: [%g %g %g %g] ", */
/* site, */
/* dot_prod[0], */
/* dot_prod[1], */
/* dot_prod[2], */
/* dot_prod[3], */
/* p_lk_left, */
/* p_lk_left[0], */
/* p_lk_left[1], */
/* p_lk_left[2], */
/* p_lk_left[3], */
/* p_lk_rght, */
/* p_lk_rght[0], */
/* p_lk_rght[1], */
/* p_lk_rght[2], */
/* p_lk_rght[3]); */
/* p_lk_left += b->left->tax ? ns : ns*ncatg; */
/* p_lk_rght += b->rght->tax ? ns : ns*ncatg; */
/* dot_prod += ns; */
}
else
{
......@@ -151,61 +127,112 @@ void AVX_Update_Eigen_Lr(t_edge *b, t_tree *tree)
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
phydbl AVX_Lk_Core_One_Class_No_Eigen_Lr(const phydbl *p_lk_left, const phydbl *p_lk_rght, const phydbl *Pij, const phydbl *pi, const int ns, const int ambiguity_check, const int observed_state)
phydbl AVX_Lk_Core_One_Class_No_Eigen_Lr(const phydbl *p_lk_left, const phydbl *p_lk_rght, const phydbl *Pij, const phydbl *tPij, const phydbl *pi, const int ns, const int ambiguity_check, const int observed_state)
{
const unsigned int sz = (int)BYTE_ALIGN / 8;
const unsigned nblocks = ns/sz;
if(nblocks == 1)
{
__m256d _plk_r,_plk_l;
if(ambiguity_check == NO) // tip case.
{
Pij += observed_state*ns;
_plk_r = _mm256_mul_pd(_mm256_load_pd(Pij),_mm256_load_pd(p_lk_left));
return pi[observed_state] * AVX_Vect_Norm(_plk_r);
}
else
{
unsigned int i;
__m256d _pijplk,_pij;
_plk_r = _mm256_mul_pd(_mm256_load_pd(p_lk_rght),_mm256_load_pd(pi));
_pijplk = _mm256_setzero_pd();
for(i=0;i<ns;++i)
{
_pij = _mm256_load_pd(tPij);
#if (defined(__FMA__))
_pijplk = _mm256_fmadd_pd(_pij,_mm256_set1_pd(p_lk_left[i]),_pijplk);
#else
_pijplk = _mm256_add_pd(_pijplk,_mm256_mul_pd(_pij,_mm256_set1_pd(p_lk_left[i])));
#endif
tPij += ns;
}
_plk_l = _mm256_mul_pd(_pijplk,_plk_r);
return(AVX_Vect_Norm(_plk_l));
}
return UNLIKELY;
}
else
{
__m256d _plk_l[nblocks],_plk_r[nblocks];
__m256d _plk; // parent partial likelihood
__m256d _plk;
/* [ Pi . Lkr ]' x Pij x Lkl */
if(ambiguity_check == NO) // tip case.
{
unsigned int i;
for(i=0;i<nblocks;++i) _plk_l[i] = _mm256_load_pd(p_lk_left + i*sz);
for(i=0;i<nblocks;++i) _plk_r[i] = _mm256_load_pd(Pij + observed_state*ns + i*sz);
for(i=0;i<nblocks;++i) _plk_r[i] = _mm256_mul_pd(_plk_r[i],_mm256_set1_pd(pi[observed_state]));
for(i=0;i<nblocks;++i) _plk_r[i] = _mm256_mul_pd(_plk_r[i],_plk_l[i]);
Pij += observed_state*ns;
for(i=0;i<nblocks;++i)
{
_plk_l[i] = _mm256_load_pd(p_lk_left);
_plk_r[i] = _mm256_load_pd(Pij);
_plk_r[i] = _mm256_mul_pd(_plk_r[i],_plk_l[i]);
p_lk_left += sz;
Pij += sz;
}
_plk = _mm256_setzero_pd();
for(i=0;i<nblocks;++i) _plk = _mm256_add_pd(_plk,_plk_r[i]);
return AVX_Vect_Norm(_plk);
return pi[observed_state] * AVX_Vect_Norm(_plk);
}
else
{
unsigned int i,j,k;
__m256d _pij[sz];
__m256d _pijplk[sz];
phydbl lk=.0;
for(i=0;i<nblocks;++i) _plk_r[i] = _mm256_load_pd(p_lk_rght + i*sz);
for(i=0;i<nblocks;++i) _plk_r[i] = _mm256_mul_pd(_plk_r[i],_mm256_load_pd(pi + i*sz));
for(i=0;i<nblocks;++i) _plk_l[i] = _mm256_load_pd(p_lk_left + i*sz);
unsigned int i,j;
__m256d _pij[nblocks],_pijplk[nblocks];
phydbl lk;
for(j=0;j<nblocks;++j)
for(i=0;i<nblocks;++i)
{
for(i=0;i<sz;++i) _pijplk[i] = _mm256_setzero_pd();
_plk_r[i] = _mm256_mul_pd(_mm256_load_pd(p_lk_rght),_mm256_load_pd(pi));
p_lk_rght += sz;
pi += sz;
}
for(i=0;i<nblocks;++i)
for(i=0;i<nblocks;++i) _pijplk[i] = _mm256_setzero_pd();
for(i=0;i<ns;++i)
{
for(k=0;k<sz;++k)
for(j=0;j<nblocks;++j)
{
_pij[k] = _mm256_load_pd(Pij + j*nblocks*sz*sz + i*sz + k*ns);
_pij[k] = _mm256_mul_pd(_pij[k],_plk_l[i]);
_pijplk[k] = _mm256_add_pd(_pijplk[k],_pij[k]);
_pij[j] = _mm256_load_pd(tPij);
tPij += sz;
#if (defined(__FMA__))
_pijplk[j] = _mm256_fmadd_pd(_pij[j],_mm256_set1_pd(p_lk_left[i]),_pijplk[j]);
#else
_pijplk[j] = _mm256_add_pd(_pijplk[j],_mm256_mul_pd(_pij[j],_mm256_set1_pd(p_lk_left[i])));
#endif
}
}
_plk = AVX_Horizontal_Add(_pijplk);
_plk = _mm256_mul_pd(_plk,_plk_r[j]);
lk += AVX_Vect_Norm(_plk);
}
lk = 0.0;
for(i=0;i<nblocks;++i) lk += AVX_Vect_Norm(_mm256_mul_pd(_pijplk[i],_plk_r[i]));
return lk;
}
return UNLIKELY;
}
}
//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
......@@ -217,8 +244,20 @@ phydbl AVX_Lk_Core_One_Class_Eigen_Lr(const phydbl *dot_prod, const phydbl *expl
const unsigned nblocks = ns/sz;
__m256d _prod[nblocks],_x;
for(l=0;l<nblocks;++l) _prod[l] = _mm256_load_pd(dot_prod + l*sz);
if(expl != NULL) for(l=0;l<nblocks;++l) _prod[l] = _mm256_mul_pd(_prod[l],_mm256_load_pd(expl + l*sz));
for(l=0;l<nblocks;++l)
{
_prod[l] = _mm256_load_pd(dot_prod); dot_prod += sz;
}
if(expl != NULL)
{
for(l=0;l<nblocks;++l)
{
_prod[l] = _mm256_mul_pd(_prod[l],_mm256_load_pd(expl));
expl += sz;
}
}
_x = _mm256_setzero_pd();
for(l=0;l<nblocks;++l) _x = _mm256_add_pd(_x,_prod[l]);
......@@ -244,7 +283,11 @@ void AVX_Lk_dLk_Core_One_Class_Eigen_Lr(const phydbl *dot_prod, const phydbl *ex
_y = _mm256_load_pd(expl + 4*i);
#if (defined(__FMA__))
_z = _mm256_fmadd_pd(_x,_y,_z);
#else
_z = _mm256_add_pd(_z,_mm256_mul_pd(_x,_y));
#endif
}
*lk = ((double *)&_z)[0] + ((double *)&_z)[2];
......@@ -257,11 +300,19 @@ void AVX_Lk_dLk_Core_One_Class_Eigen_Lr(const phydbl *dot_prod, const phydbl *ex
phydbl AVX_Vect_Norm(__m256d _z)
{
phydbl r;
__m256d _x = _mm256_hadd_pd(_z,_z);
__m256d _y = _mm256_permute2f128_pd(_x,_x,0x21);
_mm_store_sd(&r,_mm256_castpd256_pd128(_mm256_add_pd(_x,_y)));
return r;
__m128d vlow = _mm256_castpd256_pd128(_z);
__m128d vhigh = _mm256_extractf128_pd(_z, 1); // high 128
vlow = _mm_add_pd(vlow, vhigh); // reduce down to 128
__m128d high64 = _mm_unpackhi_pd(vlow, vlow);
return _mm_cvtsd_f64(_mm_add_sd(vlow, high64));
/* phydbl r; */
/* __m256d _x = _mm256_hadd_pd(_z,_z); */
/* __m256d _y = _mm256_permute2f128_pd(_x,_x,0x21); */
/* _mm_store_sd(&r,_mm256_castpd256_pd128(_mm256_add_pd(_x,_y))); */
/* return r; */
}
//////////////////////////////////////////////////////////////
......@@ -415,12 +466,14 @@ void AVX_Update_Partial_Lk(t_tree *tree, t_edge *b, t_node *d)
AVX_Partial_Lk_Inin(_tPij1,plk1,_pmat1plk1,
_tPij2,plk2,_pmat2plk2,
ns,_plk0);
}
for(k=0;k<nblocks;++k) _mm256_store_pd(plk0+sz*k,_plk0[k]);
_tPij1 += nsns / sz;
_tPij2 += nsns / sz;
plk0 += ns;
plk1 += (n_v1->tax) ? 0 : ns;
plk2 += (n_v2->tax) ? 0 : ns;
......@@ -580,7 +633,14 @@ void AVX_Matrix_Vect_Prod(const __m256d *_m_transpose, const phydbl *_v, const i
for(i=0;i<ns;++i)
{
_x = _mm256_set1_pd(_v[i]);
for(j=0;j<nblocks;++j) _u[j] = _mm256_add_pd(_u[j],_mm256_mul_pd(_m_transpose[j],_x));
for(j=0;j<nblocks;++j)
{
#if (defined(__FMA__))
_u[j] = _mm256_fmadd_pd(_m_transpose[j],_x,_u[j]);
#else
_u[j] = _mm256_add_pd(_u[j],_mm256_mul_pd(_m_transpose[j],_x));
#endif
}
_m_transpose = _m_transpose + nblocks;
}
}
......