o evolve.phylo() and plot.ancestral() have been removed.
+ o chronogram(), ratogram(), and NPRS.criterion() have been removed.
+
OTHER CHANGES
Package: ape
Version: 2.5-2
-Date: 2010-05-06
+Date: 2010-05-14
Title: Analyses of Phylogenetics and Evolution
-Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne
+Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer, Damien de Vienne
Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
Depends: R (>= 2.6.0)
Suggests: gee
-## ace.R (2010-04-23)
+## ace.R (2010-05-12)
## Ancestral Character Estimation
print(x$call)
cat("\n")
cat(" Log-likelihood:", x$loglik, "\n\n")
- cat("Rate index matrix:\n")
ratemat <- x$index.matrix
- dimnames(ratemat)[1:2] <- dimnames(x$lik.anc)[2]
- print(ratemat, na.print = ".")
- cat("\n")
- npar <- length(x$rates)
- estim <- data.frame(1:npar, round(x$rates, digits), round(x$se, digits))
- cat("Parameter estimates:\n")
- names(estim) <- c("rate index", "estimate", "std-err")
- print(estim, row.names = FALSE)
- cat("\nScaled likelihoods at the root (type 'x$lik.anc' to get them for all nodes):\n")
- print(x$lik.anc[1, ])
+ if (is.null(ratemat)) { # to be improved
+ class(x) <- NULL
+ x$loglik <- x$call <- NULL
+ print(x)
+ } else {
+ dimnames(ratemat)[1:2] <- dimnames(x$lik.anc)[2]
+ cat("Rate index matrix:\n")
+ print(ratemat, na.print = ".")
+ cat("\n")
+ npar <- length(x$rates)
+ estim <- data.frame(1:npar, round(x$rates, digits), round(x$se, digits))
+ cat("Parameter estimates:\n")
+ names(estim) <- c("rate index", "estimate", "std-err")
+ print(estim, row.names = FALSE)
+ cat("\nScaled likelihoods at the root (type 'x$lik.anc' to get them for all nodes):\n")
+ print(x$lik.anc[1, ])
+ }
}
plot.ancestral <- function(...)
.Defunct(msg = '\'plot.ancestral\' has been removed from ape,
see help("ape-defunct") for details.')
+
+chronogram <- function(...)
+ .Defunct(msg = '\'chronogram\' has been removed from ape,
+ see help("ape-defunct") for details.')
+
+ratogram <- function(...)
+ .Defunct(msg = '\'ratogram\' has been removed from ape,
+ see help("ape-defunct") for details.')
+
+NPRS.criterion <- function(...)
+ .Defunct(msg = '\'NPRS.criterion\' has been removed from ape,
+ see help("ape-defunct") for details.')
+++ /dev/null
-## nprs.R (2003-07-11)
-
-## Nonparametric Rate Smoothing Method by Sanderson
-
-## Copyright 2003 Gangolf Jobb and Korbinian Strimmer
-
-## This file is part of the R-package `ape'.
-## See the file ../COPYING for licensing issues.
-
-setTree <-
- function(lowerNodes,upperNodes,edgeLengths,minEdgeLength,tipLabels)
- .C(
- "setTree",
- as.integer(lowerNodes),
- as.integer(upperNodes),
- as.double(edgeLengths),
- as.double(minEdgeLength),
- as.integer(length(edgeLengths)),
- as.character(tipLabels),
- as.integer(length(tipLabels)),
- result=integer(1),
- PACKAGE = "ape"
- )$result
-
-getNFreeParams <-
- function()
- .C(
- "getNFreeParams",
- result=integer(1),
- PACKAGE = "ape"
- )$result
-
-getNEdges <-
- function()
- .C(
- "getNEdges",
- result=integer(1),
- PACKAGE = "ape"
- )$result
-
-getEdgeLengths <-
- function()
- .C(
- "getEdgeLengths",
- result=double(getNEdges()),
- PACKAGE = "ape"
- )$result
-
-objFuncLogScale <-
- function(params,expo)
- .C(
- "objFuncLogScale",
- as.double(params),
- as.integer(expo),
- result=double(1),
- PACKAGE = "ape"
- )$result
-
-getDurations <-
- function(params,scale)
- .C(
- "getDurations",
- as.double(params),
- as.double(scale),
- result=double(getNEdges()),
- PACKAGE = "ape"
- )$result
-
-getRates <-
- function(params,scale)
- .C(
- "getRates",
- as.double(params),
- as.double(scale),
- result=double(getNEdges()),
- PACKAGE = "ape"
- )$result
-
-getExternalParams <-
- function()
- .C(
- "getExternalParams",
- result=double(getNFreeParams()),
- PACKAGE = "ape"
- )$result
-
-### private functions
-
-prepareTree <- function(phy, minEdgeLength = 1e-06)
-{
- len <- phy$edge.length
- if (length(len) > 2048) stop("Only 2048 branches in tree allowed!")
- low <- phy$edge[, 1] # edges in the tree
- upp <- phy$edge[, 2]
- setTree(low, upp, len, minEdgeLength, phy$tip.labels)
-}
-
-optimTree <- function(phy, expo = 2) # call prepareTree first
-{
- dur <- rep(log(0.5), getNFreeParams() ) # start value
- objL <- function(d) objFuncLogScale(d, expo)
- opt <- optim(dur, objL, method = "BFGS")
- return(opt)
-}
-
-### this is just for testing purposes, to get the tree we are
-### actually using when there are many small branch lengths
-phylogram <- function(phy, ...)
-{
- if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
-
- ## added by EP for the new coding of "phylo" (2006-10-04):
- phy <- new2old.phylo(phy)
- ## End
-
- prepareTree(phy, ...)
- ##opt <- optimTree(phy, ...)
-
- newTree <- phy
- newTree$edge.length <- getEdgeLengths()
-
- ans <- newTree
- old2new.phylo(ans)
-}
-
-### public functions
-
-chronogram <- function(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
-{
- if (class(phy) != "phylo") stop("object \"phy\" is not of class \"phylo\"")
-
- ## added by EP for the new coding of "phylo" (2006-10-04):
- phy <- new2old.phylo(phy)
- ## End
-
- prepareTree(phy, minEdgeLength = minEdgeLength)
- opt <- optimTree(phy, expo = expo)
-
- newTree <- phy
- newTree$edge.length <- getDurations(opt$par, scale)
-
- ans <- newTree
- old2new.phylo(ans)
-}
-
-ratogram <- function(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
-{
- if (class(phy) != "phylo")
- stop("object \"phy\" is not of class \"phylo\"")
-
- ## added by EP for the new coding of "phylo" (2006-10-04):
- phy <- new2old.phylo(phy)
- ## End
-
- prepareTree(phy, minEdgeLength = minEdgeLength)
- opt <- optimTree(phy, expo = expo)
-
- newTree <- phy
- newTree$edge.length <- getRates(opt$par, scale)
-
- ans <- newTree
- old2new.phylo(ans)
-}
-
-NPRS.criterion <- function(phy, chrono, expo = 2, minEdgeLength = 1e-06)
-{
- if (!is.ultrametric(chrono))
- stop("tree \"chrono\" is not ultrametric (clock-like)")
-
- ## added by EP for the new coding of "phylo" (2006-10-04):
- phy <- new2old.phylo(phy)
- chrono <- new2old.phylo(chrono)
- ## End
-
- prepareTree(chrono, minEdgeLength = minEdgeLength)
- parms <- getExternalParams()
- prepareTree(phy, minEdgeLength = minEdgeLength)
- objFuncLogScale(parms, expo)
-}
+++ /dev/null
-\name{NPRS.criterion}
-\alias{NPRS.criterion}
-\title{Objective Function Employed in Nonparametric Rate Smoothing}
-\usage{
-NPRS.criterion(phy, chrono, expo = 2, minEdgeLength = 1e-06)
-}
-\arguments{
- \item{phy}{A non-clock-like phylogenetic tree (i.e. an object of class
- \code{"phylo"}), where the branch lengths are measured in
- substitutions.}
- \item{chrono}{A chronogram, i.e. a clock-like tree (i.e. an object of
- class \code{"phylo"}), where the branch lengths are measured in
- absolute time.}
- \item{expo}{Exponent in the objective function (default value: 2)}
- \item{minEdgeLength}{Minimum edge length in the phylogram (default
- value: 1e-06). If any branch lengths are smaller then they will be
- set to this value.}
-}
-\description{
- \code{NPRS.criterion} computes the objective function to be minimized
- in the NPRS (nonparametric rate smoothing) algorithm described in
- Sanderson (1997).
-}
-\details{
- Please refer to Sanderson (1997) for mathematical details. Note that
- is is not computationally efficient to optimize the branch lengths in
- a chronogram by using \code{NPRS.criterion} - please use
- \code{\link{chronogram}} instead.
-}
-\value{
- \code{NPRS.criterion} returns the value of the objective function given
- a phylogram and a chronogram.
-}
-\author{Gangolf Jobb (\url{http://www.treefinder.de}) and Korbinian
- Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})
-}
-\seealso{
- \code{\link{ratogram}}, \code{\link{chronogram}}
-}
-\references{
- Sanderson, M. J. (1997) A nonparametric approach to estimating
- divergence times in the absence of rate constancy. \emph{Molecular
- Biology and Evolution}, \bold{14}, 1218--1231.
-}
-\examples{
-# get tree
-data("landplants.newick") # example tree in NH format
-tree.landplants <- read.tree(text = landplants.newick)
-
-# plot tree
-tree.landplants
-plot(tree.landplants, label.offset = 0.001)
-
-# estimate chronogram
-chrono.plants <- chronogram(tree.landplants)
-
-# plot
-plot(chrono.plants, label.offset = 0.001)
-
-# value of NPRS function for our estimated chronogram
-NPRS.criterion(tree.landplants, chrono.plants)
-}
-\keyword{manip}
\name{ace}
\alias{ace}
+\alias{print.ace}
\alias{logLik.ace}
\alias{deviance.ace}
\alias{AIC.ace}
ace(x, phy, type = "continuous", method = "ML", CI = TRUE,
model = if (type == "continuous") "BM" else "ER",
scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1)
+\method{print}{ace}(x, digits = 4, ...)
\method{logLik}{ace}(object, ...)
\method{deviance}{ace}(object, ...)
\method{AIC}{ace}(object, ..., k = 2)
\method{anova}{ace}(object, ...)
}
\arguments{
- \item{x}{a vector or a factor.}
+ \item{x}{a vector or a factor; an object of class \code{"ace"} in the
+ case of \code{print}.}
\item{phy}{an object of class \code{"phylo"}.}
\item{type}{the variable type; either \code{"continuous"} or
\code{"discrete"} (or an abbreviation of these).}
structure to be used (this also gives the assumed model).}
\item{ip}{the initial value(s) used for the ML estimation procedure
when \code{type == "discrete"} (possibly recycled).}
+ \item{digits}{the number of digits to be printed.}
\item{object}{an object of class \code{"ace"}.}
\item{k}{a numeric value giving the penalty per estimated parameter;
the default is \code{k = 2} which is the classical Akaike
\alias{theta.s}
\alias{evolve.phylo}
\alias{plot.ancestral}
+\alias{chronogram}
+\alias{ratogram}
+\alias{NPRS.criterion}
\title{Defunct Ape Functions}
\description{
These functions have been removed from \pkg{ape} or moved to another
theta.s(s, n, variance = FALSE)
evolve.phylo(phy, value, var)
plot.ancestral(...)
+chronogram(...)
+ratogram(...)
+NPRS.criterion(...)
}
\details{
\code{klastorin} has been removed because it does not seem to be used
and \code{theta.s} have been moved to \pkg{pegas}.
\code{evolve.phylo} and \code{plot.ancestral} have been deprecated by
- the new function \link{\code{rTraitCont}}.
+ the new function \code{\link{rTraitCont}}.
+
+ \code{chronogram}, \code{ratogram}, and \code{NPRS.criterion} have
+ ceased to be maintained: consider using \code{\link{chronopl}}.
}
\keyword{internal}
\alias{mant.zstat}
\alias{lower.triang}
\alias{clado.build}
-\alias{getDurations}
-\alias{getEdgeLengths}
-\alias{getExternalParams}
-\alias{getNEdges}
-\alias{getNFreeParams}
-\alias{getRates}
-\alias{objFuncLogScale}
-\alias{optimTree}
-\alias{phylogram}
-\alias{prepareTree}
-\alias{setTree}
-\alias{nEdges}
-\alias{nNodes}
\alias{phylogram.plot}
\alias{cladogram.plot}
\alias{circular.plot}
\author{
Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard
- Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Gangolf Jobb,
- Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim
- Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian
- Strimmer, Damien de Vienne
+ Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph
+ Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon,
+ Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Korbinian Strimmer,
+ Damien de Vienne
Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
}
plot(z, show.node.label = TRUE, font = 1, root.edge = TRUE)
title("z <- bind.tree(x, y, 2, .1)")
-x <- rtree(100)
-y <- rtree(100)
+x <- rtree(50)
+y <- rtree(50)
x$root.edge <- y$root.edge <- .2
z <- x + y
plot(y, show.tip.label = FALSE, root.edge = TRUE); axisPhylo()
}
\author{Emmanuel Paradis}
\seealso{
- \code{\link{chronogram}}, \code{\link{ratogram}},
- \code{\link{NPRS.criterion}}, \code{\link{chronopl}}
+ \code{\link{chronopl}}
}
\examples{
tr <- rtree(10)
+++ /dev/null
-\name{chronogram}
-\alias{chronogram}
-
-\title{Chronogram Computed by Nonparametric Rate Smoothing}
-\usage{
-chronogram(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
-}
-\arguments{
- \item{phy}{A phylogenetic tree (i.e. an object of class
- \code{"phylo"}), where the branch lengths are measured in substitutions.}
- \item{scale}{Age of the root in the inferred chronogram (default value: 0). }
- \item{expo}{Exponent in the objective function (default value: 2)}
- \item{minEdgeLength}{Minimum edge length in the phylogram (default
- value: 1e-06). If any branch lengths are smaller then they will be
- set to this value.}
-}
-\description{
- \code{chronogram} computes a chronogram from a phylogram by applying the NPRS
- (nonparametric rate smoothing) algorithm described in Sanderson (1997).
-}
-\details{
- Please refer to Sanderson (1997) for mathematical details
-}
-\value{
-\code{chronogram} returns an object of class \code{"phylo"}. The branch lengths of this
-tree will be clock-like and scaled so that the root node has age 1 (or the value
-set by the option \code{scale}
-}
-\author{
- Gangolf Jobb (\url{http://www.treefinder.de}) and
-Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})
-}
-\seealso{
-\code{\link{ratogram}}, \code{\link{NPRS.criterion}}.
-}
-\references{
- Sanderson, M. J. (1997) A nonparametric approach to estimating
- divergence times in the absence of rate constancy. \emph{Molecular
- Biology and Evolution}, \bold{14}, 1218--1231.
-}
-\examples{
-# get tree
-data("landplants.newick") # example tree in NH format
-tree.landplants <- read.tree(text = landplants.newick)
-
-# plot tree
-tree.landplants
-plot(tree.landplants, label.offset = 0.001)
-
-# estimate chronogram
-chrono.plants <- chronogram(tree.landplants)
-
-# plot
-plot(chrono.plants, label.offset = 0.001)
-
-# value of NPRS function for our estimated chronogram
-NPRS.criterion(tree.landplants, chrono.plants)
-}
-\keyword{manip}
}
\author{Emmanuel Paradis}
\seealso{
- \code{\link{chronogram}}, \code{\link{ratogram}},
- \code{\link{NPRS.criterion}}, \code{\link{chronoMPL}}
+ \code{\link{chronoMPL}}
}
\keyword{models}
data example in the software package r8s
(\url{http://ginger.ucdavis.edu/r8s/}).
}
-\seealso{
-\code{\link{chronogram}}, \code{\link{ratogram}}, \code{\link{NPRS.criterion}}.
-}
\references{
Sanderson, M. J. (1997) A nonparametric approach to estimating
divergence times in the absence of rate constancy. \emph{Molecular
+++ /dev/null
-\name{ratogram}
-\alias{ratogram}
-
-\title{Ratogram Computed by Nonparametric Rate Smoothing}
-\usage{
-ratogram(phy, scale = 1, expo = 2, minEdgeLength = 1e-06)
-}
-\arguments{
- \item{phy}{A phylogenetic tree (i.e. an object of class \code{"phylo"}), where
- the branch lengths are measured in substitutions.}
-
- \item{scale}{Age of the root in the chronogram corresponding to the inferred ratogram(default value: 0). }
-
- \item{expo}{Exponent in the objective function (default value: 2)}
- \item{minEdgeLength}{Minimum edge length in the phylogram (default value: 1e-06). If any branch lengths are
- smaller then they will be set to this value. }
-}
-\description{
-
- \code{ratogram} computes a ratogram from a phylogram by applying the NPRS
- (nonparametric rate smoothing) algorithm described in Sanderson (1997).
-}
-\details{
- Please refer to Sanderson (1997) for mathematical details
-}
-\value{
-\code{chronogram} returns an object of class \code{"phylo"}. The branch lengths of this
-tree will be the absolute rates estimated for each branch.
-}
-\author{Gangolf Jobb (\url{http://www.treefinder.de}) and
-Korbinian Strimmer (\url{http://www.stat.uni-muenchen.de/~strimmer/})
-}
-\seealso{
-\code{\link{chronogram}}, \code{\link{NPRS.criterion}}.
-}
-\references{
- Sanderson, M. J. (1997) A nonparametric approach to estimating
- divergence times in the absence of rate constancy. \emph{Molecular
- Biology and Evolution}, \bold{14}, 1218--1231.
-}
-\examples{
-# get tree
-data("landplants.newick") # example tree in NH format
-tree.landplants <- read.tree(text = landplants.newick)
-
-# plot tree
-tree.landplants
-plot(tree.landplants, label.offset = 0.001)
-
-# estimate ratogram
-rato.plants <- ratogram(tree.landplants)
-
-# plot
-plot(rato.plants, label.offset = 0.001)
-}
-\keyword{manip}
+++ /dev/null
-/*
- * nprsfunc.c
- *
- * (c) 2003 Gangolf Jobb and Korbinian Strimmer
- *
- * Functions for nonparametric rate smoothing (NPRS)
- * (see MJ Sanderson. 1997. MBE 14:1218-1231)
- *
- * This code may be distributed under the GNU GPL
- */
-
-
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-
-/*============================================================================*/
-
-#define EPS 1e-6
-#define LARGENUM 1e99
-
-/* original scale */
-#define MINPARAM EPS
-#define MAXPARAM 1-EPS
-#define STARTPARAM 0.5
-
-
-/* log scale */
-#define MINLPARAM log(MINPARAM)
-#define MAXLPARAM log(MAXPARAM)
-#define STARTLPARAM log(STARTPARAM)
-
-#define ARRLEN 2048
-#define LABLEN 64
-
-#define index aperates_index
-
-int tree_lowerNodes[ARRLEN];
-int tree_upperNodes[ARRLEN];
-double tree_edgeLengths[ARRLEN];
-int tree_nedges;
-char tree_tipLabels[ARRLEN][LABLEN];
-int tree_ntips;
-
-int nparams,index[ARRLEN],ancestor[ARRLEN];
-
-/*============================================================================*/
-
-void setTree(
- int *lowerNodes,
- int *upperNodes,
- double *edgeLengths,
- double *minEdgeLength,
- int *nedges,
- char **tipLabels,
- int *ntips,
- int *result
-) {
- int i;
-
- tree_nedges=*nedges;
-
- for(i=0;i<tree_nedges;i++) {
- tree_lowerNodes[i]=lowerNodes[i];
- tree_upperNodes[i]=upperNodes[i];
- tree_edgeLengths[i]=(edgeLengths[i]<*minEdgeLength)?*minEdgeLength:edgeLengths[i];
- }
-
- tree_ntips=*ntips;
-
- for(i=0;i<tree_ntips;i++) {
- strcpy(tree_tipLabels[i],tipLabels[i]);
- }
-
- nparams=0;
- for(i=0;i<tree_nedges;i++) {
- if(tree_lowerNodes[i]<0&&tree_upperNodes[i]<0) {index[-tree_upperNodes[i]]=nparams++;}
- if(tree_upperNodes[i]<0) ancestor[-tree_upperNodes[i]]=tree_lowerNodes[i];
- }
-
- *result=0;
-}
-
-/*============================================================================*/
-
-void getNFreeParams(int *result) { *result=nparams; }
-
-/*============================================================================*/
-
-void getNEdges(int *result) { *result=tree_nedges; }
-
-/*============================================================================*/
-
-void getEdgeLengths(double *result) {
- int i;
-
- for(i=0;i<tree_nedges;i++) result[i]=tree_edgeLengths[i];
-
-}
-
-/*============================================================================*/
-
-double age(double *params,int node) {
- double prod;
-
- if(node>=0) return(0.0);
- if(node==-1) return(1.0);
-
- prod=0.0;
- while(node!=-1) {prod+=params[index[-node]]; node=ancestor[-node];}
-
- return(exp(prod));
-}
-
-
-void getDurations(double *params,double *scale,double *result) {
- int i,low,upp;
-
- for(i=0;i<tree_nedges;i++) {
- low=tree_lowerNodes[i]; upp=tree_upperNodes[i];
- if(low<0&&upp<0) result[i]=*scale*(age(params,low)-age(params,upp));
- else result[i]=*scale*age(params,low);
- }
-
-}
-
-/*============================================================================*/
-
-void getRates(double *params,double *scale,double *result) {
- int i,low,upp;
-
- for(i=0;i<tree_nedges;i++) {
- low=tree_lowerNodes[i]; upp=tree_upperNodes[i];
- if(low<0&&upp<0) result[i]=tree_edgeLengths[i]/(*scale*(age(params,low)-age(params,upp)));
- else result[i]=tree_edgeLengths[i]/(*scale*age(params,low));
- }
-
-}
-
-/*============================================================================*/
-
-
-/* note on the choice of parameters:
- *
- * in order to obtain a chronogram we optimize the
- * relative heights of each internal node, i.e. to
- * each internal node (minus the root) we assign a number
- * between 0 and 1 (MINPARAM - MAXPARAM).
- *
- * To obtain the actual height of a node the relative heights
- * of the nodes have to multiplied (in fact we work on the log-scale
- * so we simply sum).
- */
-
-
-/* internal objective function - parameters are on log scale */
-void nprsObjectiveFunction(
- double *params,
- int *expo,
- double *result
-) {
- int p,k,j;
- double me,rk,rj,sum,scale=1.0,durations[ARRLEN],rates[ARRLEN];
-
- p=*expo;
-
- getDurations(params,&scale,durations);
-
- for(k=0;k<tree_nedges;k++) rates[k]=tree_edgeLengths[k]/durations[k];
-
- me=0.; for(k=0;k<tree_nedges;k++) me+=rates[k]; me/=tree_nedges;
-
- sum=0.;
-
- for(k=0;k<tree_nedges;k++) { rk=rates[k];
- for(j=0;j<tree_nedges;j++) { if(j==k) continue; rj=rates[j];
- if(tree_lowerNodes[j]==-1) sum+=pow(fabs(me-rj),p); else
- if(tree_upperNodes[k]==tree_lowerNodes[j]) sum+=pow(fabs(rk-rj),p);
- }
- }
-
- *result=sum;
-
-}
-
-
-/* check parameter bounds on log scale */
-int checkLogParams(double *params)
-{
- int i;
- for(i=0; i<nparams; i++)
- {
- if(params[i] > MAXLPARAM || params[i] < MINLPARAM) return 0; /* out of bounds */
- }
- return 1; /* within bounds */
-}
-
-
-
-/*
- * public objective function
- * - parameters are on log scale
- * - if parameters are out of bounds function returns large value
- */
-void objFuncLogScale(
- double *params,
- int *expo,
- double *result
-)
-{
-
- if( checkLogParams(params) == 0 ) /* out of bounds */
- {
- *result = LARGENUM;
- }
- else
- {
- nprsObjectiveFunction(params, expo, result);
- }
-}
-
-
-/*============================================================================*/
-
-
-double ageAlongEdges(int node) { /* tree must be clock-like */
- int i;
-
- if(node>=0) return(0.);
-
- for(i=0;i<tree_nedges;i++)
- if(tree_lowerNodes[i]==node)
- return(tree_edgeLengths[i]+ageAlongEdges(tree_upperNodes[i]));
-
- return(0.);
-}
-
-/* compute set of parameters for a given clock-like tree */
-void getExternalParams(double *result) {
- int i;
-
- for(i=0;i<tree_nedges;i++) {
- if(tree_lowerNodes[i]<0&&tree_upperNodes[i]<0)
- result[index[-tree_upperNodes[i]]]
- =-log(ageAlongEdges(tree_lowerNodes[i])/ageAlongEdges(tree_upperNodes[i]));
- }
-
-}
-
-/*============================================================================*/
projects / ape.git / commitdiff