Remove redundancy in vcv and contrasts methods for lambda and eb.

diversitree/R/continuous.R

@@ -76,3 +76,71 @@ toC.cache.continuous <- function(cache) {
 
   cache
 }
+
+make.all.branches.rescale.vcv <- function(cache, control) {
+  n.tip <- cache$n.tip
+  states <- cache$states
+  states.sd <- cache$states.sd
+
+  rescale <- make.rescale.phylo(cache$info$phy, cache$info$name)
+
+  function(pars, intermediates) {
+    s2 <- pars[[1]]
+
+    vcv <- vcv.phylo(rescale(pars[-1]))
+    vv <- s2 * vcv
+
+    diag(vv) <- diag(vv) + states.sd^2
+    mu <- phylogMean(vv, states)
+    dmvnorm2(states, rep(mu, n.tip), vv, solve(vv), log=TRUE)
+  }
+}
+
+make.all.branches.rescale.contrasts <- function(cache, control) {
+  if (any(cache$states.sd > 0))
+    stop("Cannot (yet) do contrasts based bm models with state error")
+
+  ## This is annoying; we'll have to do this in make.bm.  If we don't,
+  ## then we will reorder every time and that's going to hurt,
+  ## timewise.  If I rewrite the pic() calculations to use my
+  ## ordering, we can skip this step.  Given that eventually I think I
+  ## want access to more of the pic components, that seems like a good
+  ## idea.
+  tree <- reorder(cache$info$phy, "pruningwise")
+  states <- cache$states[tree$tip.label] # does reorder change this?
+  rescale <- make.rescale.phylo(tree, cache$info$name)
+  n <- length(tree$tip.label)
+
+  function(pars, intermediates, preset=NULL) {
+    s2 <- pars[[1]]
+    tree <- rescale(pars[-1])
+
+    ## Copied from model-bm.R:make.cache.bm
+    ## There is duplication here with make.cache.pgls; perhaps merge
+    ## these?  That might help with some of the root treatment
+    ## things.
+    pics <- pic(cache$states, tree, var.contrasts=TRUE)
+    u <- pics[,"contrasts"]
+    V <- pics[,"variance"]
+    V0 <- pgls.root.var.bm(tree)
+    ## This step is brutal.  We could rewrite the branch lengths there
+    ## once the lookups are done properly.  This is particularly bad
+    ## because it will involve doing another tree reordering.
+    ## Costly!  What's more is we can actually get it from the PIC
+    ## calculation.  So this part needs to be done separately, and for
+    ## all the pic methods.  This entire section of code is horribly
+    ## duplicated!
+    root.x <- pgls.root.mean.bm(tree, cache$states)
+
+    ## This is the bit that is shared with all.branches.bm
+    ll <- -(n * log(2 * pi * s2) +
+            sum(log(V)) +
+            log(V0) +
+            sum(u * u) / s2) * 0.5
+    list(loglik=ll,
+         root.x=root.x,
+         root.v=s2 * V0,
+         # not sure if these are needed...
+         V=V, V0=V0)
+  }
+}

diversitree/R/model-eb.R

@@ -9,7 +9,7 @@ make.eb <- function(tree, states, states.sd=0, control=list()) {
   cache <- make.cache.eb(tree, states, states.sd, control)
 
   if (control$method == "vcv") {
-    all.branches <- make.all.branches.eb.vcv(cache, control)
+    all.branches <- make.all.branches.rescale.vcv(cache, control)
     rootfunc <- rootfunc.bm.vcv
   } else {
     all.branches <- make.all.branches.eb.pruning(cache, control)
@@ -70,30 +70,6 @@ check.pars.eb <- function(pars) {
   TRUE
 }
 
-## VCV approach:
-##
-## Notice here that this now works for all rescaling functions - so
-## long as they take parameters except for the first one (s2).
-make.all.branches.eb.vcv <- function(cache, control) {
-  n.tip <- cache$n.tip
-  states <- cache$states
-  states.sd <- cache$states.sd
-
-  rescale <- make.rescale.phylo.eb(cache$info$phy)
-
-  function(pars, intermediates) {
-    s2 <- pars[[1]]
-    a  <- pars[[2]]
-
-    vcv <- vcv.phylo(rescale(a))
-    vv <- s2 * vcv
-    # Below here is identical to make.all.branches.ou.vcv
-    diag(vv) <- diag(vv) + states.sd^2
-    mu <- phylogMean(vv, states)
-    dmvnorm2(states, rep(mu, n.tip), vv, solve(vv), log=TRUE)
-  }
-}
-
 make.all.branches.eb.pruning <- function(cache, control) {
   ## NOTE: This is a hack, but allow here for the extra paramter:
   cache$info$np <- 3L
@@ -113,7 +89,6 @@ make.all.branches.eb.pruning <- function(cache, control) {
     all.branches(c(pars, pars.extra), ...)
 }
 
-
 ## The issue that I have here is that time is computed against the
 ## root, so we'll need to know that.  t0 is the *depth* of the
 ## branch *tip* so t0+len is the *depth* of the branch base:

diversitree/R/model-lambda.R

@@ -4,13 +4,13 @@ make.lambda <- function(tree, states, states.sd=0, control=list()) {
   cache <- make.cache.lambda(tree, states, states.sd, control)
 
   if (control$method == "vcv") {
-    all.branches <- make.all.branches.lambda.vcv(cache, control)
+    all.branches <- make.all.branches.rescale.vcv(cache, control)
     rootfunc <- rootfunc.bm.vcv
   } else if (control$method == "pruning") {
     all.branches <- make.all.branches.lambda.pruning(cache, control)
     rootfunc <- rootfunc.bm.pruning
   } else if (control$method == "contrasts") {
-    all.branches <- make.all.branches.lambda.contrasts(cache, control)
+    all.branches <- make.all.branches.rescale.contrasts(cache, control)
     rootfunc <- rootfunc.bm.contrasts
   } else {
     stop("Unknown method", control$method)
@@ -72,27 +72,6 @@ check.pars.lambda <- function(pars) {
   TRUE
 }
 
-## VCV approach:
-make.all.branches.lambda.vcv <- function(cache, control) {
-  n.tip <- cache$n.tip
-  states <- cache$states
-  states.sd <- cache$states.sd
-
-  rescale <- make.rescale.phylo.lambda(cache$info$phy)
-
-  function(pars, intermediates) {
-    s2     <- pars[[1]]
-    lambda <- pars[[2]]
-
-    vcv <- vcv.phylo(rescale(lambda))
-    vv <- s2 * vcv
-    # Below here is identical to make.all.branches.ou.vcv
-    diag(vv) <- diag(vv) + states.sd^2
-    mu <- phylogMean(vv, states)
-    dmvnorm2(states, rep(mu, n.tip), vv, solve(vv), log=TRUE)
-  }
-}
-
 make.all.branches.lambda.pruning <- function(cache, control) {
   ## NOTE: This is a hack, but allow here for the extra parameters
   cache$info$np <- 4L

diversitree/R/rescale.R

@@ -90,3 +90,12 @@ rescale.phylo.se <- function(phy, se) {
   phy$edge.length[tips] <- phy$edge.length[tips] + se
   phy
 }
+
+## A little wrapping helper function.
+make.rescale.phylo <- function(phy, model) {
+  switch(model,
+         ou=make.rescale.phylo.eb,
+         eb=make.rescale.phylo.eb,
+         lambda=make.rescale.phylo.lambda,
+         stop("Unknown model ", dQuote(model)))(phy)
+}

diversitree/inst/tests/test-lambda.R

@@ -22,17 +22,17 @@ lik.pru.R <- make.lambda(phy, states,
                      control=list(method="pruning", backend="R"))
 lik.pru.C <- make.lambda(phy, states,
                      control=list(method="pruning", backend="C"))
-## lik.con <- make.lambda(phy, states, control=list(method="contrasts"))
+lik.con <- make.lambda(phy, states, control=list(method="contrasts"))
 
 lik.vcv.se <- make.lambda(phy, states, se, control=list(method="vcv"))
 lik.pru.R.se <- make.lambda(phy, states, se,
                      control=list(method="pruning", backend="R"))
 lik.pru.C.se <- make.lambda(phy, states, se,
                      control=list(method="pruning", backend="C"))
-## ## Not yet supported
-## expect_that(make.lambda(phy, states, se,
-##                         control=list(method="contrasts")),
-##             throws_error())
+## Not yet supported
+expect_that(make.lambda(phy, states, se,
+                        control=list(method="contrasts")),
+            throws_error())
 
 test_that("Likelihood calculations agree on known case", {
   ## Start with BM; set lambda to one: