richfitz / diversitree Goto Github PK

diversitree: comparative phylogenetic analyses of diversification

Home Page: http://www.zoology.ubc.ca/prog/diversitree

R 58.88% Shell 0.15% TeX 13.17% C 10.01% C++ 2.84% Fortran 9.63% Python 0.10% CSS 0.43% Makefile 0.28% HTML 4.27% M4 0.24%

diversitree's Introduction

diversitree: comparative phylogenetic analyses of diversification

This repository contains "diversitree". It is my experimental sources though, and may not compile or work for you. You may prefer the released version from CRAN:

> install.packages("diversitree")

The interesting directories are:

inst/tests: testing functions that exercise most of the package's main features. Running ./run_tests.R will run the tests. These take too long to run on CRAN (over a minute), so are not set up in the usual way.
doc: Vignettes, and their required data files.

Installing diversitree

Clone the repository with

git clone git://github.com/richfitz/diversitree.git

The package should then be installable the usual way. You'll need a C, C++ and Fortran compiler.

To install and specify the location of the fftw library in a non-standard place, a line like this is required: R CMD INSTALL diversitree --configure-args='--with-fftw=/path/to/fftw' where the path will be the path so that the files /path/to/fftw/include/fftw3.h /path/to/fftw/lib/lib/fftw3.so exist.

On Windows, set the environment variable LIB_FFTW to point to the directory that contains include/fftw.h, and install the package.

If fftw is not found, installation will continue, but the (relatively) fast C based QuaSSE integration will not be available. The R based fft integrator and the method-of-lines integrator will be available.

Branches

The "master" branch contains the bleeding edge version of diversitree. It may or may not work for you. The "release" branch contains stable releases.

diversitree's People

Contributors

Stargazers

Watchers

Forkers

eeg gburin acgerstein lslsnow yanhaifei jakejing lzhangss rdinnager jmhallas

diversitree's Issues

OU calculations on non-ultrametric trees

See email from Pennell (7/12/2013), http://www.zoology.ubc.ca/prog/diversitree/examples/ou-nonultrametric/ and

library(diversitree)
set.seed(2)
tree <- rtree(100)
x <- rTraitCont(tree)
lik.ou.v <- make.ou(tree, x, control=list(method="vcv"))
lik.ou.p  <- make.ou(tree, x,
                     control=list(method="pruning", backend="C"))

Bug in custom root state with find.mle

Surprised that this is not tested, but it's a bug.

library(diversitree)
pars <- c(.1, .1, .03, .03, .1, .2)
set.seed(3)
phy <- trees(pars, "bisse", max.taxa=25, max.t=Inf, x0=0)[[1]]
#fitting a Mk2 model without and with root constraint
lik <- make.mk2(phy, phy$tip.state)
fit.mk2 <- find.mle(lik, c(.1, .1))
fit.mk2.root <- find.mle(lik,c(.1,.1), root=ROOT.GIVEN, root.p=c(0,1))
#fitting an ER model without and with root constraint
lik2=constrain(lik,q01~q10)
fit.er <- find.mle(lik2, 0.1)
fit.er.root <- find.mle(lik2, 0.1, root=ROOT.GIVEN, root.p=c(0,1), keep.func=TRUE) # error
fit.er.root <- find.mle(lik2, 0.1, root=ROOT.GIVEN, root.p=c(0,1), keep.func=FALSE)

Fix is just a bit off faffing in do.mle.search

Reported by Sébastien Lavergne

Error : 'tree' must be ultrametric

I'm using Rsudio and when I run the script of GeoSSE I got this error:
Error in check.tree(tree) : 'tree' must be ultrametric

Although the tree is ultrametric.
Thank you!

Installed diversitree from source - Error: external pointer is not valid

Hi there,

I am using R 3.1.3 GUI 1.65 Mavericks build on OSX Yosemite 10.10.1. The past few days I have been throwing constant errors only when using diversitree similar to that described in richfitz/wood#1. The error is along the lines *** caught segfault ***
address 0x0, cause 'memory not mapped'. This has occurred on all 3 of my computers, and often got to the stage where I couldn't even load diversitree without throwing this error or R completely crashing. To fix this problem, I tried installing Rcpp and diversitree as source packages as you suggested in the previous issues post, which I had some issues with but eventually got to work.

Now I am able to run functions in diversitree without the segfault error (or I haven't discovered it yet), however I cannot use the find.mle() without throwing the following error:
Error: external pointer is not valid.
The other functions I am using make.bisse() and starting.point.bisse() do seem to work fine though.
Is there something in the find.mle() that could cause this error?

If you could help at all I would be so grateful, as I really need these functions to work on my computer urgently!

Thank you in advance,

Emma

MKn model in diversitree with deSolve and custom initial conditions not working

I'm trying to reproduce the old 'extending diversitree' analysis with the current diversitree. I think I have all the pieces together, but my example dies due to the output of the derivatives not being what is expected. The traceback is as follows.

 Error in checkFunc(Func2, times, y, rho) : 
  The number of derivatives returned by func() (1) must equal the length of the initial conditions vector (2) 
11 stop(paste("The number of derivatives returned by func() (", 
    length(tmp[[1]]), ") must equal the length of the initial conditions vector (", 
    length(y), ")", sep = "")) 
10 checkFunc(Func2, times, y, rho) 
9 lsoda(...) 
8 t(lsoda(...)[-1, -1, drop = FALSE]) 
7 lsoda.trim(vars, times, derivs, pars, rtol = rtol, atol = atol) 
6 ode(y, c(t0, t0 + len), pars) 
5 branches(y, len, pars, t0, idx) 
4 branches(x$y, x$t, pars, 0, idx) 
3 all.branches.matrix(pars, cache, initial.conditions, branches, 
    preset) 
2 all.branches(qmat, intermediates) 
1 ll(c(1, 0), root = ROOT.GIVEN, root.p = c(1, 0), intermediates = TRUE)

when fixing CRAN issue, future proof things a bit

wch/r-source@ea54ad8

Minor bug when node is not matched in get.descendants

When calling get.descendants() if the node is not matched, the following error happens.
Error in sprintf("Node '%s' not found in tree") : too few arguments

This is due to a small typo when calling sprintf().
stop(sprintf("Node '%s' not found in tree"), node)

It should be
stop(sprintf("Node '%s' not found in tree", node))

make.musse should bail if given a clade.tree

Reported by Daniel Caetano:

tree <- read.nexus("Thomas-tree.nex")
d <- read.csv("Lislevand-states.csv", as.is=TRUE)
states <- d$dimorph
names(states) <- d$species
states <- states[tree$tip.label]
names(states) <- tree$tip.label
states.15 <- (abs(states) > 0.15) + 0
tree.clade <- clades.from.polytomies(tree)

## Making the data for MuSSE:
states.15[which(states.15 == 1)] <- 2
states.15[which(states.15 == 0)] <- 1
states.15[sample(1:350, size = 50)] <- 3
lik <- make.musse(tree.clade, states = states.15, k = 3)
> Error in dt.tips.grouped(y, y.i, cache) : y must be same length as tips

crash in fortran code with new clang

docker pull ghcr.io/r-hub/containers/clang17
docker run -v ${PWD}:/src -it --rm ghcr.io/r-hub/containers/clang17
apt-get update && apt-get -y install gdb git libfftw3-dev libgsl-dev vim
R -e 'options(repos = "https://cloud.r-project.org"); install.packages(c("Rcpp", "ape", "deSolve", "subplex", "Rcpp"))'

# Fresh
git clone https://github.com/richfitz/diversitree
git -C diversitree checkout cran-fixes

# There are lots of concerning warnings in about the problematic part
# of the fortran code here.
R CMD INSTALL diversitree

cat <<EOF > err.R
library(diversitree)
RNGkind(sample.kind = "Rounding")
pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01)
set.seed(4)
phy <- tree.bisse(pars, max.t=30, x0=0)
unresolved <- data.frame(tip.label=I(c("sp25", "sp30", "sp40", "sp56", "sp20")),
                         Nc =c(10, 9, 6, 5, 2),
                         n0=0, n1=0)
unresolved.bd <- structure(unresolved[["Nc"]], names=unresolved[["tip.label"]])
lik.bisse3 <- make.bisse(phy, phy[["tip.state"]], unresolved)
lik.bisse3(pars)
EOF
Rscript err.R

Confirmed with images 8c7713c7033e and 4a6a9279096b

Remove deprecated autoconf option

It is no longer necessary to use these in configure scripts, and 'R CMD
config CPP' has been deprecated. R does not set CPP so R CMD config CPP is either "$CC -E" or empty, depending on the 'make' program. And
autoconf will call AC_PROG_CPP to determine CPP if it needs it.

The only use by a configure.ac script I noticed was in nloptr, which
passes it down. If needed, set CPP="${CC} -E".

QuaSSE gives confusing error message when diffusion ommitted

Should just check number of parameters and avoid this error.

    library(diversitree)
    lambda <- function(x) sigmoid.x(x, 0.1, 0.2,  0, 2.5)
    mu <- function(x) constant.x(x, 0.03)
    char <- make.brownian.with.drift(0, 0.025)

    set.seed(1)
    phy <- tree.quasse(c(lambda, mu, char), max.taxa=15, x0=0,
                       single.lineage=FALSE, verbose=TRUE)

    lik <- make.quasse(phy, phy$tip.state, 0.5, constant.x, constant.x)
    p <- starting.point.quasse(phy, phy$tip.state)
    lik(p)

models in GeoSSE

Hi FitzJohn,
How many models do we have in GeoSSE? Could we use constraint model like dA ~dB?
Thank you!
Best,
Uzma

Multiple records with same name

When a vector contains multiple species with the same name, diversitree uses only the first observation. The correct behavior should be to throw an error. I just corrected this in geiger
mwpennell/geiger-v2@b7a7d45

MWE

library(diversitree)

## simulate a continuous dataset
set.seed(1)
t <- tree.bd(c(1,0), max.taxa=10)
d <- sim.character(t,1)

## using the simulated dataset
lik <- make.bm(t,d)
lik(1)

## adding a species with same name to data
d2 <- c(d, "sp10"=1)
lik2 <- make.bm(t,d2)

## likelihood is the same
lik2(1)

## however, swapping order of species
d3 <- d2[c(1:9,11,10)]
lik3 <- make.bm(t,d3)

## different likelihood
lik3(1)

Problem with asr.stoch and irreversible model

Dear @richfitz

Nice job on diversitree - I've just started to play with it.

I'm trying to run stochastic mapping using the example of an irreversible mk2 model given here

The following code works fine with asr.joint and asr.marginal, but throws an error with asr.stoch

pars <- c(.1, .1, .03, .03, .03, .06)
set.seed(1)
phy <- trees(pars, "bisse", max.taxa=4, max.t=Inf, x0=0)[[1]]
h <- history.from.sim.discrete(phy, 0:1)
plot(h, phy, main="True history")
lik <- make.mk2(phy, phy$tip.state)
argnames(lik)
lik2<-constrain(lik, q10~0)
argnames(lik2)
fit2 = find.mle(lik2,c(.1))
st.s2 <- asr.stoch(lik2, fit2$par)

The above gives the following error:

Error in asr(lik(pars, pars.only = TRUE), ...) : 
  Exceeding maximum allowed history size

I don't get this for the unconstrained model when doing stochastic mapping. I've also tried calling make.asr.stoch.constrained explicitly, but I got the same error.

How to plot QuaSSE output

Hi,

I've seen in some papers plots of the model fits for quasse. In my analysis a hump shaped model with drift fit best. Do you know how I can plot this?

Jamie

Help with MuSSE csv file

Hi everyone!

I am struggling with a constant problem with a csv extension while preparing data for MuSSE model analysis: have tried to do a bunch of stuff to fix a problem but no success - always the same thing ("All names must be length 1"). I would very grateful for your help! :)

library(diversitree)
dat="MuSSE_hosts.csv"
dat<- read.table("MuSSE_hosts.csv", header=TRUE, dec=".", sep=",", row.names=1)
mat <- dat[,2:ncol(dat)]
lik.0 <- make.musse.multitrait(tree, mat, depth=0)
Error in check.states.musse.multitrait(tree, states, strict = strict, :
All names must be length 1

I can't install diversitree in linux (R version 4.1.0 (2021-05-18) -- "Camp Pontanezen")

I can Install diversitree, I am always having the error

/usr/bin/ld: /usr/local/lib/libfftw3.a(assert.o): relocation R_X86_64_PC32 against symbol `stdout@@GLIBC_2.2.5' can not be used when making a shared object; recompile with -fPIC
/usr/bin/ld: final link failed: bad value
collect2: error: ld returned 1 exit status
make: *** [/usr/share/R/share/make/shlib.mk:10: diversitree.so] Error 1
ERROR: compilation failed for package ‘diversitree’

Even if I clone the github repository and install from the console, I am always getting the same error....

I am on linux mint 20.1, and yes I am sure I have the libfftw3 installed

libfftw3-dev is already the newest version (3.3.8-2ubuntu1).

Any help may be welcome

cheers

protection errors

Package diversitree version 0.9-10
Package built using 76078/R 3.6.0; x86_64-pc-linux-gnu; 2019-02-10 10:27:31 UTC; unix   
Checked with rchk version c132302c5d39a55c436a06df73e539a1bace7486
More information at https://github.com/kalibera/cran-checks/blob/master/rchk/PROTECT.md

Function GslOdeR::target(double, SEXPREC*)
  [UP] calling allocating function Rf_eval with argument allocated using Rf_lang4(?,V,?,?) diversitree/src/GslOdeR.cpp:20

CRAN LTO checks

A number of things here that need resolving, mostly fairly easy

Bug in info of musse

Edit of the argnames of musse objects fails like below.

argnames(lik)[1]<-c("l")
Error in argnames<-.dtlik(*tmp*, value = c("l", "lambda2", "lambda3", :
Wrong argnames length: expected 15, got 12

Then I checked the source and found in "model-musse.R"

make.info.musse <- function(k, phy) {
list(name="musse", # canonical name (as in make.)
name.pretty="MuSSE", # pretty name for printing
## Parameters:
np=as.integer(k * (k + 2)), # number of ODE parameters

The parameters should be composed of k lambda, k mu and k*(k-1) transition rates.
Then np=as.integer(k * (k + 1)) seems to be correct.

Problem with profiles.plot

This won't work for alpha or theta...

yule25.trees<-structure(list(edge = structure(c(26L, 27L, 28L, 29L, 29L, 30L, 
30L, 31L, 31L, 28L, 32L, 33L, 34L, 35L, 35L, 34L, 33L, 36L, 37L, 
37L, 36L, 38L, 39L, 40L, 40L, 39L, 41L, 41L, 42L, 42L, 38L, 32L, 
43L, 43L, 44L, 44L, 27L, 45L, 46L, 46L, 45L, 26L, 47L, 48L, 49L, 
49L, 48L, 47L, 27L, 28L, 29L, 1L, 30L, 2L, 31L, 3L, 4L, 32L, 
33L, 34L, 35L, 5L, 6L, 7L, 36L, 37L, 8L, 9L, 38L, 39L, 40L, 10L, 
11L, 41L, 12L, 42L, 13L, 14L, 15L, 43L, 16L, 44L, 17L, 18L, 45L, 
46L, 19L, 20L, 21L, 47L, 48L, 49L, 22L, 23L, 24L, 25L), .Dim = c(48L, 
2L)), Nnode = 24L, tip.label = c("t1", "t2", "t3", "t4", "t5", 
"t6", "t7", "t8", "t9", "t10", "t11", "t12", "t13", "t14", "t15", 
"t16", "t17", "t18", "t19", "t20", "t21", "t22", "t23", "t24", 
"t25"), edge.length = c(0.5904258688, 0.8249415727, 2.470189261, 
0.932987671, 0.7898975338, 0.1430901373, 0.01988442392, 0.1232057133, 
0.1232057133, 0.4996568724, 0.138501434, 1.739476059, 0.3999097194, 
0.6256328474, 0.6256328474, 1.025542567, 1.180119854, 1.198822986, 
0.386075786, 0.386075786, 0.3431831631, 0.2316611904, 0.3650179037, 
0.6450365148, 0.6450365148, 0.1565813483, 0.8534730703, 0.7328424265, 
0.1206306438, 0.1206306438, 1.241715609, 2.049187285, 0.8543327749, 
0.2658454834, 0.5884872915, 0.5884872915, 0.9436889333, 2.704861689, 
0.5795678831, 0.5795678831, 3.284429572, 3.050576216, 0.1083979538, 
0.1808128489, 1.478757356, 1.478757356, 1.659570205, 1.767968158
)), .Names = c("edge", "Nnode", "tip.label", "edge.length"), class = "phylo", order = "cladewise")

yule25.data<-data.frame(t1 = 2.24, t10= 2.86, t11= 0.576, t12=0.244, t13= 1.82, t14= 1.33,t15= 2.78, t16= 0.985, t17=1.77, t18= 1.65, t19 = 3.8, t2 =-0.531, t20 = 2.7, t21 = 2.06, t22= -3.73, t23= -1.19, t24= -2.29, t25= -0.507, t3 = -0.065, t4 = -0.568, t5 =1.59, t6 =1.48,  t7 =1.48, t8 = 2.45, t9 =0.688)

#Make likelihood function for OU
lik.ou <- make.ou(yule25.trees, yule25.data, control = list(method="pruning", backend = "C"))

#Get MLE
fit.ou.mle<-find.mle(lik.ou, c(0.1,0.1,0.1))

#Set start parameters (s2, alpha, theta) for MCMC from MLE
pars<-fit.ou.mle$par

#Set prior (uniformative exponential)
prior.exp<-make.prior.exponential(0.1)

#Find values for w (tuning parameter) using short MCMC run
# = 95% confidence interval of parameters from short MCMC run
#Needs abs or theta < 0

tmp <- mcmc(lik.ou, abs(pars), nsteps = 100, w = 0.1, prior = prior.exp)

    w <- diff(sapply(tmp[2:4], quantile, c(.05, .95)))

#Run MCMC
out <- mcmc(lik.ou, abs(pars), nsteps=10000, w=w, prior = prior.exp)

# Look at the marginal distribution of parameters
profiles.plot(out["alpha"], col.line="red") #Error in xy.coords(x, y) : 'x' and 'y' lengths differ
profiles.plot(out["theta"], col.line="blue") #Error in xy.coords(x, y) : 'x' and 'y' lengths differ
profiles.plot(out["s2"], col.line="green") #Runs

Missing references in make.ou

Perhaps Hansen 1997? And Felsenstein? Blomberg 200X?

Error in drift hump model

What does this mean, and how can I get around it?

Error in branches(y.in, len[i], pars, depth[i], i) :
Actual negative D value detected -- calculation failure

mcmc run

I got this error when I ran mcmc for best fitted model which is 'full.4c'
Error in mcmc(model.4c.full, p.4c.mcmc.full, nsteps = 100, prior = prior, :
object 'model.4c.full' not found

Thank you!

Upside down bell-curve in best-fitting modal model in QuaSSE

The best fitting model in QuaSSE for my data is a modal model. However, this is an upside down bell curve where speciation rates are high at all values except for a dip of the bell curve.

Is this a statistical thing? I know that some ML optimisation methods (e.g. RPANDA), resulting values should be considered as absolute values if extinction is negative. This would make sense with my data when plotting trait values against tip-rates as estimated by the DR statistic and BAMM.

constrain behaves weirdly if doubly constraining a parameter:

Reported by @dwbapst:

library(diversitree)
pars <- c(0.1, 0.2, 0.03, 0.03, 0.01, 0.01)
set.seed(2)
phy <- tree.bisse(pars, max.t=60, x0=0)
lik <- make.bisse(phy, phy$tip.state)
lik.2 <- constrain(lik, lambda0 ~ lambda1, lambda0 ~ mu0)

last line generates warning:

Warning message:
In rels[i] <- res[[2]] : (symbol) object cannot be coerced to type 'list'

but the likelihood function doubly constrains lambda0:

BiSSE likelihood function:
  * Parameter vector takes 5 elements:
     - lambda1, mu0, mu1, q01, q10
  * Function constrained (original took 6 elements):
     - lambda0 ~ lambda1
     - lambda0 ~ mu0
  * Function takes arguments (with defaults)
     - pars: Parameter vector
     - ...: Additional arguments to underlying function
     - pars.only [FALSE]: Return full parameter vector?
  * Phylogeny with 203 tips and 202 nodes
     - Taxa: sp1, sp27, sp29, sp31, sp41, sp43, sp46, sp47, sp48, ...
  * References:
     - Maddison et al. (2007) doi:10.1080/10635150701607033
     - FitzJohn et al. (2009) doi:10.1093/sysbio/syp067
R definition:
function (pars, ..., pars.only = FALSE)

lik.2(2:6, pars.only=TRUE)

gives

lambda0 lambda1     mu0     mu1     q01     q10 
      3       2       3       4       5       6

Update for CRAN

Fails to configure with compiler flags that Apple will start using:

See the configure logs at

https://www.stats.ox.ac.uk/pub/bdr/NoImplicitDecl/

These are the config.log from installing the package on a Linux machine
with the C compiler flag -Werror=implicit-function-declaration and
indicate problems with your configure script.  You should see a similar
warning even without that flag.

The reason for looking at these now is that Apple has (essentially)
added that flag by default in the toolchain (12.x) it released earlier
this month, and so configure may not work as you intended on Macs with
that toolchain.  There are comments in the logs which indicate the
corrections needed.

Other toolchains may follow suit: there has been a bug report on gcc for
quite some time and careful users may already be adding this flag.

Please correct before 2020-12-08 to safely retain your package on CRAN.

richfitz / diversitree Goto Github PK

diversitree's Introduction

diversitree: comparative phylogenetic analyses of diversification

Installing diversitree

Branches

diversitree's People

Contributors

Stargazers

Watchers

Forkers

diversitree's Issues

Recommend Projects

Recommend Topics

Recommend Org