richfitz / diversitree

diversitree: comparative phylogenetic analyses of diversification
http://www.zoology.ubc.ca/prog/diversitree
30 stars 9 forks source link

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

Closed sdwfrost closed 7 years ago

sdwfrost commented 8 years ago

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) 
richfitz commented 8 years ago

Terribly sorry about this. diversitree was a casualty of a major shift in research focus between PhD and postdoc, which combined with emergency patch updates for CRAN has left it in a less than ideal state and I'm not spending any time on it all anymore...

It looks like the extending diversitree manual is incorrect and you should omit the list(ret) line (just returning ret). So:

derivs <- function(t, y, pars) {
  k <- length(y)
  Q <- matrix(pars, k, k)
  ret <- Q %*% y
  dim(ret) <- NULL
  ret
}

appears to work for me.

sdwfrost commented 8 years ago

Thanks - this works for me too - not sure how though, as ode and lsoda (in the traceback) expect that the model function returns a list.

No apologies necessary! The codebase and docs are quite a lot to maintain for one person.

sdwfrost commented 7 years ago

Looks like this can be closed. A quick question; is there any hard-coded reason (such as the way that the tree is traversed) why models are restricted to ultrametric trees?

richfitz commented 7 years ago

I can't remember - most of the models in the package are restricted to ultrametric trees because they're diversification models. It's possible that the code that the mkn models uses uses code that assumes an ultrametric tree but I'd have to check. I don't think that's the case though