R/nprs.R

   1 ## nprs.R (2003-07-11)
   2
   3 ##   Nonparametric Rate Smoothing Method by Sanderson
   4
   5 ## Copyright 2003 Gangolf Jobb and Korbinian Strimmer
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 setTree <-
  11   function(lowerNodes,upperNodes,edgeLengths,minEdgeLength,tipLabels)
  12    .C(
  13     "setTree",
  14     as.integer(lowerNodes),
  15     as.integer(upperNodes),
  16     as.double(edgeLengths),
  17     as.double(minEdgeLength),
  18     as.integer(length(edgeLengths)),
  19     as.character(tipLabels),
  20     as.integer(length(tipLabels)),
  21     result=integer(1),
  22     PACKAGE = "ape"
  23    )$result
  24
  25 getNFreeParams <-
  26   function()
  27    .C(
  28     "getNFreeParams",
  29     result=integer(1),
  30     PACKAGE = "ape"
  31    )$result
  32
  33 getNEdges <-
  34   function()
  35    .C(
  36     "getNEdges",
  37     result=integer(1),
  38     PACKAGE = "ape"
  39    )$result
  40
  41 getEdgeLengths <-
  42   function()
  43    .C(
  44     "getEdgeLengths",
  45     result=double(getNEdges()),
  46     PACKAGE = "ape"
  47    )$result
  48
  49 objFuncLogScale <-
  50   function(params,expo)
  51    .C(
  52     "objFuncLogScale",
  53     as.double(params),
  54     as.integer(expo),
  55     result=double(1),
  56     PACKAGE = "ape"
  57    )$result
  58
  59 getDurations <-
  60   function(params,scale)
  61    .C(
  62     "getDurations",
  63     as.double(params),
  64     as.double(scale),
  65     result=double(getNEdges()),
  66     PACKAGE = "ape"
  67    )$result
  68
  69 getRates <-
  70   function(params,scale)
  71    .C(
  72     "getRates",
  73     as.double(params),
  74     as.double(scale),
  75     result=double(getNEdges()),
  76     PACKAGE = "ape"
  77    )$result
  78
  79 getExternalParams <-
  80   function()
  81    .C(
  82     "getExternalParams",
  83     result=double(getNFreeParams()),
  84     PACKAGE = "ape"
  85    )$result
  86
  87 ### private functions
  88
  89 prepareTree <- function(phy, minEdgeLength = 1e-06)
  90 {
  91     len <- phy$edge.length
  92     if (length(len) > 2048) stop("Only 2048 branches in tree allowed!")
  93     low <- phy$edge[, 1] # edges in the tree
  94     upp <- phy$edge[, 2]
  95     setTree(low, upp, len, minEdgeLength, phy$tip.labels)
  96 }
  97
  98 optimTree <- function(phy, expo = 2) # call prepareTree first
  99 {
 100     dur <- rep(log(0.5), getNFreeParams() ) # start value
 101     objL <- function(d) objFuncLogScale(d, expo)
 102     opt <- optim(dur, objL, method = "BFGS")
 103     return(opt)
 104 }
 105
 106 ### this is just for testing purposes, to get the tree we are
 107 ### actually using when there are many small branch lengths
 108 phylogram <- function(phy, ...)
 109 {
 110     if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
 111
 112     ## added by EP for the new coding of "phylo" (2006-10-04):
 113     phy <- new2old.phylo(phy)
 114     ## End
 115
 116     prepareTree(phy, ...)
 117     ##opt <- optimTree(phy, ...)
 118
 119     newTree <- phy
 120     newTree$edge.length <- getEdgeLengths()
 121
 122     ans <- newTree
 123     old2new.phylo(ans)
 124 }
 125
 126 ### public functions
 127
 128 chronogram <- function(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
 129 {
 130     if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
 131
 132     ## added by EP for the new coding of "phylo" (2006-10-04):
 133     phy <- new2old.phylo(phy)
 134     ## End
 135
 136     prepareTree(phy, minEdgeLength = minEdgeLength)
 137     opt <- optimTree(phy, expo = expo)
 138
 139     newTree <- phy
 140     newTree$edge.length <- getDurations(opt$par, scale)
 141
 142     ans <- newTree
 143     old2new.phylo(ans)
 144 }
 145
 146 ratogram <- function(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
 147 {
 148     if (class(phy) != "phylo")
 149       stop("object \"phy\" is not of class \"phylo\"")
 150
 151     ## added by EP for the new coding of "phylo" (2006-10-04):
 152     phy <- new2old.phylo(phy)
 153     ## End
 154
 155     prepareTree(phy, minEdgeLength = minEdgeLength)
 156     opt <- optimTree(phy, expo = expo)
 157
 158     newTree <- phy
 159     newTree$edge.length <- getRates(opt$par, scale)
 160
 161     ans <- newTree
 162     old2new.phylo(ans)
 163 }
 164
 165 NPRS.criterion <- function(phy, chrono, expo = 2, minEdgeLength = 1e-06)
 166 {
 167     if (!is.ultrametric(chrono))
 168       stop("tree \"chrono\" is not ultrametric (clock-like)")
 169
 170     ## added by EP for the new coding of "phylo" (2006-10-04):
 171     phy <- new2old.phylo(phy)
 172     chrono <- new2old.phylo(chrono)
 173     ## End
 174
 175     prepareTree(chrono, minEdgeLength = minEdgeLength)
 176     parms <- getExternalParams()
 177     prepareTree(phy, minEdgeLength = minEdgeLength)
 178     objFuncLogScale(parms, expo)
 179 }