o The functions heterozygosity, nuc.div, theta.h, theta.k and
theta.s have been moved from ape to pegas.
+ o The functions mlphylo, DNAmodel and sh.test have been removed.
+
CHANGES IN APE VERSION 2.3-3
Package: ape
Version: 2.4
-Date: 2009-09-30
+Date: 2009-10-03
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
Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
.Defunct(msg = '\'klastorin\' has been removed from ape,
see help("ape-defunct") for details.')
+mlphylo <-
+ function(x, phy, model = DNAmodel(), search.tree = FALSE,
+ quiet = FALSE, value = NULL, fixed = FALSE)
+ .Defunct(msg = '\'mlphylo\' has been removed from ape,
+ see help("ape-defunct") for details.')
+
+DNAmodel <- function(model = "K80", partition = 1,
+ ncat.isv = 1, invar = FALSE,
+ equal.isv = TRUE, equal.invar = 1)
+ .Defunct(msg = '\'DNAmodel\' has been removed from ape,
+ see help("ape-defunct") for details.')
+
+sh.test <- function(..., x, model = DNAmodel(), B = 100)
+ .Defunct(msg = '\'sh.test\' has been removed from ape,
+ see help("ape-defunct") for details.')
+
heterozygosity <- function (x, variance = FALSE)
.Defunct(msg = '\'heterozygosity\' has been moved from ape to pegas,
see help("ape-defunct") for details.')
+++ /dev/null
-## mlphylo.R (2008-07-15)
-
-## Estimating Phylogenies by Maximum Likelihood
-
-## Copyright 2006-2008 Emmanuel Paradis
-
-## This file is part of the R-package `ape'.
-## See the file ../COPYING for licensing issues.
-
-logLik.phylo <- function(object, ...) attr(object, "loglik")
-
-deviance.phylo <- function(object, ...) -2*attr(object, "loglik")
-
-AIC.phylo <- function(object, ..., k = 2)
-{
- np <- length(object$edge.length) +
- length(attr(object, "rates")) +
- length(attr(object, "alpha")) +
- length(attr(object, "invar")) +
- length(attr(object, "xi"))
- if (!attr(object, "model") %in% c("JC69", "F81"))
- np <- np + 3
- -2*attr(object, "loglik") + k*np
-}
-
-.subst.model <- structure(c(0, 1, 0, 1, 1, 1, 2, 5),
- names = c("JC69", "K80", "F81", "F84",
- "HKY85", "T92", "TN93", "GTR"))
-
-mlphylo <-
- function(x, phy, model = DNAmodel(), search.tree = FALSE,
- quiet = FALSE, value = NULL, fixed = FALSE)
-{
- ## not yet generic....
- if (class(x) != "DNAbin") stop("DNA sequences not in binary format")
- if (!is.binary.tree(phy))
- stop("the initial tree must be dichotomous.")
- if (!quiet && is.rooted(phy)) {
- warning("the initial tree is rooted: it will be unrooted.")
- phy <- unroot(phy)
- }
- if (is.null(phy$edge.length))
- stop("the initial tree must have branch lengths.")
- if (any(phy$edge.length > 1))
- stop("some branch lengths are greater than one.")
- phy <- reorder(phy, "pruningwise")
- if (!quiet) cat("Preparing the sequences...\n")
- if (is.list(x)) x <- as.matrix(x)
- if (is.null(rownames(x)))
- stop("DNA sequences have no names") # safe...
- if (!all(names(x) %in% phy$tip.label))
- stop("the names of the DNA sequences and the tip labels
-of the tree do not match") # safe here also
- x <- x[phy$tip.label, ]
- Y <- prepareDNA(x, model)
- BF <- if (Y$model %in% 1:2) rep(0.25, 4) else base.freq(x)
- S <- length(Y$weight)
- npart <- dim(Y$partition)[2] # the number of overall partitions
- ## in case of negative branch lengths:
- phy$edge.length <- abs(phy$edge.length)
- nb.tip <- length(phy$tip.label)
- para <- if (Y$npara) rep(1, Y$npara) else 0
- alpha <- if (Y$nalpha) rep(.5, Y$nalpha) else 0
- invar <- if (Y$ninvar) rep(0.5, Y$ninvar) else 0
-
- if (!is.null(value)) {
- if (para && !is.null(value$rates))
- para <- value$rates[1:Y$npara]
- if (alpha && !is.null(value$alpha))
- alpha <- value$alpha[1:Y$nalpha]
- if (invar && !is.null(value$invar))
- invar <- value$invar[1:Y$ninvar]
- }
-
- loglik <- 0
- if (!quiet) cat("Fitting in progress... ")
- res <- .C("mlphylo_DNAmodel", as.integer(nb.tip), as.integer(S),
- as.raw(Y$SEQ), as.double(Y$ANC), as.double(Y$w),
- as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
- as.double(phy$edge.length), as.integer(npart),
- as.integer(Y$partition), as.integer(Y$model),
- as.double(Y$xi), as.double(para), as.integer(Y$npara),
- as.double(alpha), as.integer(Y$nalpha),
- as.integer(Y$ncat), as.double(invar), as.integer(Y$ninvar),
- as.double(BF), as.integer(search.tree), as.integer(fixed),
- as.double(loglik), NAOK = TRUE, PACKAGE = "ape")
- if (!quiet) cat("DONE!\n")
- phy$edge.length = res[[8]]
- attr(phy, "loglik") <- res[[23]]
- attr(phy, "npart") <- npart
- attr(phy, "model") <- names(Y$npara)
- if (para) attr(phy, "rates") <- res[[13]]
- if (alpha) attr(phy, "alpha") <- res[[15]]
- if (invar) attr(phy, "invar") <- res[[18]]
- if (npart > 1) attr(phy, "xi") <- res[[12]]
- phy
-}
-
-DNAmodel <- function(model = "K80", partition = 1,
- ncat.isv = 1, invar = FALSE,
- equal.isv = TRUE, equal.invar = 1)
-{
- if (ncat.isv > 10)
- stop("number of categories for inter-site variation cannot exceed 10")
- structure(list(model = model, partition = partition,
- ncat.isv = ncat.isv, invar = invar,
- equal.isv = equal.isv, equal.invar = equal.invar),
- class = "DNAmodel")
-}
-
-prepareDNA <- function(X, DNAmodel)
-{
- L <- dim(X)[2] # already converted as a matrix in mlphylo()
-
- npart <- length(unique(DNAmodel$partition))
-
- ## find which substitution model:
- mo <- which(names(.subst.model) == DNAmodel$model)
- npara <- .subst.model[mo] # keeps the 'names'
-
- ## inter-sites variation:
- nalpha <- as.numeric(DNAmodel$ncat.isv > 1)
- if (!DNAmodel$equal.isv) nalpha <- npart * nalpha
-
- ## proportion of invariants:
- ninvar <- as.numeric(DNAmodel$invar)
- if (!DNAmodel$equal.invar) ninvar <- npart * ninvar
-
- SEQ <- weight <- part <- NULL
-
- ## For each partition...
- for (i in 1:npart) {
- ## extracts the sites in this partition:
- M <- X[, DNAmodel$partition == i, drop = FALSE]
- ## convert each column as a character string:
- M <- apply(M, 2, rawToChar)
- ## get their frequencies:
- w <- table(M)
- ## convert back to raw the unique(M):
- M <- sapply(dimnames(w)[[1]], charToRaw)
- ## remove useless attributes:
- colnames(M) <- dimnames(w) <- NULL
- w <- unclass(w)
- ## bind everything:
- SEQ <- cbind(SEQ, M)
- weight <- c(weight, w)
- part <- c(part, length(w)) # the length of each partition
- }
-
- class(SEQ) <- "DNAbin"
- ANC <- array(0, c(nrow(SEQ) - 2, ncol(SEQ), 4))
-
- ## 'partition' gives the start and end of each partition:
- partition <- matrix(1, 2, npart)
- partition[2, ] <- cumsum(part)
- if (npart > 1) {
- partition[1, 2:npart] <- partition[2, 1:(npart - 1)] + 1
- partition[2, npart] <- length(weight)
- xi <- rep(1, npart - 1)
- } else xi <- 0
- list(SEQ = SEQ, ANC = ANC, weight = weight, partition = partition,
- model = mo, xi = xi, npara = npara, nalpha = nalpha,
- ncat = DNAmodel$ncat.isv, ninvar = ninvar)
-}
-## scales.R (2009-07-23)
+## scales.R (2009-10-02)
## Add a Scale Bar or Axis to a Phylogeny Plot
+++ /dev/null
-## sh.test.R (2008-07-11)
-
-## Shimodaira-Hasegawa Test
-
-## Copyright 2006-2008 Emmanuel Paradis
-
-## This file is part of the R-package `ape'.
-## See the file ../COPYING for licensing issues.
-
-sh.test <- function(..., x, model = DNAmodel(), B = 100)
-{
- ## Prepare the list of trees:
- phy <- list(...)
- if (length(phy) == 1 && class(phy[[1]]) != "phylo")
- phy <- unlist(phy, recursive = FALSE)
- ntree <- length(phy)
-
- ## Arrange the sequences as a matrix:
- if (is.list(x)) {
- nm <- names(x)
- n <- length(x)
- x <- unlist(x)
- nL <- length(x)
- x <- matrix(x, n, nL/n, byrow = TRUE)
- rownames(x) <- nm
- }
-
- ## Step 1:
- foo <- function(PHY)
- attr(mlphylo(x, PHY, model, search.tree = FALSE, quiet = TRUE), "loglik")
- Talpha <- sapply(phy, foo)
- Talpha <- max(Talpha) - Talpha
-
- ## Do the bootstrap resampling (Step 2):
- M <- matrix(NA, ntree, B)
- for (i in 1:B) {
- boot.samp <- x[, sample(ncol(x), replace = TRUE)]
- for (j in 1:ntree)
- M[j, i] <- attr(mlphylo(boot.samp, phy[[j]], model,
- search.tree = FALSE, quiet = TRUE),
- "loglik")
- }
- M <- M - rowMeans(M) # Step 3
- ## Step 4: <FIXME> This can greatly simplified </FIXME>
- for (i in 1:B)
- for (j in 1:ntree)
- M[j, i] <- max(M[j, i] - M[, i])
- ## Step 5:
- count <- numeric(ntree)
- for (j in 1:ntree)
- count[j] <- sum(M[j, ] > Talpha[j])
- count <- count/B
- names(count) <- names(phy)
- count
-}
+++ /dev/null
-\name{DNAmodel}
-\alias{DNAmodel}
-\title{Defines Models of DNA Evolution}
-\usage{
-DNAmodel(model = "K80", partition = 1,
- ncat.isv = 1, invar = FALSE,
- equal.isv = TRUE, equal.invar = 1)
-}
-\arguments{
- \item{model}{a vector of mode character giving the substition model.}
- \item{partition}{a vector of integers defining the partitions for the
- substition models (eventually recycled).}
- \item{ncat.isv}{the number of categories in each partition.}
- \item{invar}{a logical value specifying whether there are invariants.}
- \item{equal.isv}{a logical value specifying whether the `alpha'
- parameter is the same in all partitions; has no effet if \code{ncat
- = 1} or if \code{partition = 1}.}
- \item{equal.invar}{similar to the argument above but for the
- proportion of invariants.}
-}
-\description{
- This function defines a model of evolution for a set of DNA sequences
- with possible partitions.
-}
-\details{
- \code{partition} is recycled along the sequence: thus by default there
- is a single partition. For instance, to partition a sequence of 1000
- sites into two partitions of equal length, one will use
- \code{partition = c(rep(1, 500), rep(2, 500))}. The partitions must be
- numbered with a series of integers (1, 2, 3, ...). To partition the
- codon positions, one could do \code{partition = c(1, 1, 2)}.
-
- The substition models are the same in all partitions. Branch lengths
- are the same in all partitions up to a multiplying coefficient (the
- contrast parameter, denoted 'xi').
-
- The substitution models must be among the followings: \code{"JC69"}
- \code{"K80"}, \code{"F81"}, \code{"F84"}, \code{"HKY85"},
- \code{"T92"}, \code{"TN93"}, and \code{"GTR"}. These models (except
- HKY85 and GTR) are described in the help page of
- \code{\link{dist.dna}}.
-
- Inter-sites variation in substitution rates (ISV) is allowed by
- specifying \code{ncat.isv} greater than one.
-}
-\note{
- The result of this function is not intended to be used by the user,
- but rather to be passed to \code{\link{mlphylo}}.
-}
-\value{
- an object of class \code{"DNAmodel"} with components defined by the
- arguments of the function call.
-}
-\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
-\seealso{
- \code{\link{mlphylo}}, \code{\link{dist.dna}}
-}
-\examples{
-### the K80 model:
-mod <- DNAmodel()
-### the simplest substitution model:
-mod <- DNAmodel("JC69")
-### the classical GTR + G4 + I:
-mod <- DNAmodel("GTR", ncat.isv = 4, invar = TRUE)
-### codon-partitioning (with K80):
-mod <- DNAmodel(partition = c(1, 1, 2))
-### the same but adding inter-sites variation (the alpha parameter
-### is the same for both partitions):
-mod <- DNAmodel(partition = c(1, 1, 2), ncat.isv = 4)
-### ... and with different `alpha' for each partition:
-mod <- DNAmodel(partition = c(1, 1, 2), ncat.isv = 4, equal.isv = FALSE)
-}
-\keyword{models}
\name{ape-defunct}
\alias{klastorin}
+\alias{mlphylo}
+\alias{DNAmodel}
+\alias{sh.test}
\alias{heterozygosity}
\alias{H}
\alias{nuc.div}
}
\usage{
klastorin(phy)
+mlphylo(x, phy, model = DNAmodel(), search.tree = FALSE,
+ quiet = FALSE, value = NULL, fixed = FALSE)
+DNAmodel(model = "K80", partition = 1,
+ ncat.isv = 1, invar = FALSE,
+ equal.isv = TRUE, equal.invar = 1)
+sh.test(..., x, model = DNAmodel(), B = 100)
heterozygosity(x, variance = FALSE)
H(x, variance = FALSE)
nuc.div(x, variance = FALSE, pairwise.deletion = FALSE)
and this helped to clear some internal C code (this function may be
put back with a different coding).
+ \code{mlphylo}, \code{DNAmodel} and \code{sh.test} have been removed:
+ see the package \pkg{phangorn} for a much better implementation of
+ these methods (and others).
+
\code{heterozygosity}, \code{nuc.div}, \code{theta.h}, \code{theta.k}
and \code{theta.s} have been moved to \pkg{pegas}.
}
\alias{pos.move}
\alias{BOTHlabels}
\alias{chronopl.cv}
-\alias{prepareDNA}
\alias{floating.pie.asp}
\alias{checkLabel}
\alias{getMRCA}
\seealso{
\code{\link{nj}}, \code{\link{fastme}},
\code{\link{write.tree}}, \code{\link{read.tree}},
- \code{\link{dist.dna}}, \code{\link{mlphylo}}
+ \code{\link{dist.dna}}
}
\examples{
### From Saitou and Nei (1987, Table 1):
\seealso{
\code{\link{nj}}, \code{\link{bionj}},
\code{\link{write.tree}}, \code{\link{read.tree}},
- \code{\link{dist.dna}}, \code{\link{mlphylo}}
+ \code{\link{dist.dna}}
}
\examples{
### From Saitou and Nei (1987, Table 1):
+++ /dev/null
-\name{mlphylo}
-\alias{mlphylo}
-\alias{logLik.phylo}
-\alias{deviance.phylo}
-\alias{AIC.phylo}
-\title{Estimating Phylogenies by Maximum Likelihood}
-\usage{
-mlphylo(x, phy, model = DNAmodel(), search.tree = FALSE,
- quiet = FALSE, value = NULL, fixed = FALSE)
-\method{logLik}{phylo}(object, ...)
-\method{deviance}{phylo}(object, ...)
-\method{AIC}{phylo}(object, ..., k = 2)
-}
-\arguments{
- \item{x}{an object of class \code{"DNAbin"} giving the (aligned) DNA
- sequence data.}
- \item{phy}{an object of class \code{"phylo"} giving the tree.}
- \item{model}{an object of class \code{"DNAmodel"} giving the model to
- be fitted.}
- \item{search.tree}{a logical specifying whether to search for the best
- tree (defaults to FALSE) (not functional for the moment).}
- \item{quiet}{a logical specifying whether to display the progress of
- the analysis.}
- \item{value}{a list with elements named \code{rates}, \code{alpha},
- and \code{invar}, or at least one of these, giving the initial
- values of the parameters of the model. If \code{NULL}, some initial
- values are given internally.}
- \item{fixed}{a logical specifying whether to optimize parameters given
- in \code{value}.}
- \item{object}{an object of class \code{"phylo"}.}
- \item{k}{a numeric value giving the penalty per estimated parameter;
- the default is \code{k = 2} which is the classical Akaike
- information criterion.}
- \item{...}{further arguments passed to or from other methods.}
-}
-\description{
- \code{mlphylo} estimates a phylogenetic tree by maximum likelihood
- given a set of DNA sequences. The model of evolution is specified with
- the function \code{\link{DNAmodel}}.
-
- \code{logLik}, \code{deviance}, and \code{AIC} are generic functions
- used to extract the log-likelihood, the deviance (-2*logLik), or the
- Akaike information criterion of a tree. If no such values are
- available, \code{NULL} is returned.
-}
-\details{
- The model specified by \code{\link{DNAmodel}} is fitted using the
- standard ``pruning'' algorithm of Felsenstein (1981).
-
- The implementation of the inter-sites variation in substitution rates
- follows the methodology developed by Yang (1994).
-
- The difference among partitions is parametrized with a contrast
- parameter (denoted \eqn{\xi}{xi}) that specifies the contrast in mean
- susbtitution rate among the partitions. This methodology is inspired
- from one introduced by Yang (1996).
-
- The substitution rates are indexed column-wise in the rate matrix: the
- first rate is set to one.
-}
-\note{
- For the moment, it is not possible to estimate neither branch lengths,
- nor the topology with \code{mlphylo}. The function may estimate all other
- parameters: substitution rates, shape (\eqn{\alpha}{alpha}) of the
- inter-sites variation in substitution rates, the proportion of
- invariants, and the ``contrast'' parameter (\eqn{\xi}{xi}) among
- partitions.
-
- Alternative topologies can also be compared using likelihood-ratio
- tests (LRTs) or AICs.
-}
-\value{
- an object of class \code{"phylo"}. There are possible additional
- attributes:
-
- \item{loglik}{the maximum log-likelihood.}
- \item{npart}{the number of partitions.}
- \item{model}{the substitution model.}
- \item{rates}{the estimated substitution rates.}
- \item{invar}{the estimated proportion of invariants.}
- \item{alpha}{the estimated shape parameter of the inter-sites
- variation in substitution rates.}
-}
-\references{
- Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a
- maximum likelihood approach. \emph{Journal of Molecular Evolution},
- \bold{17}, 368--376.
-
- Yang, Z. (1994) Maximum likelihood phylogenetic estimation from DNA
- sequences with variable rates over sites: approximate methods.
- \emph{Journal of Molecular Evolution}, \bold{39}, 306--314.
-
- Yang, Z. (1996) Maximum-likelihood models for combined analyses of
- multiple sequence data. \emph{Journal of Molecular Evolution},
- \bold{42}, 587--596.
-}
-\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
-\seealso{
- \code{\link{DNAmodel}}, \code{\link{nj}}, \code{\link{read.dna}},
- \code{\link{summary.phylo}}, \code{\link{bionj}}, \code{\link{fastme}}
-}
-\keyword{models}
\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
\seealso{
\code{\link{write.tree}}, \code{\link{read.tree}},
- \code{\link{dist.dna}}, \code{\link{mlphylo}}, \code{\link{bionj}},
+ \code{\link{dist.dna}}, \code{\link{bionj}},
\code{\link{fastme}}
}
\examples{
+++ /dev/null
-\name{sh.test}
-\alias{sh.test}
-\title{Shimodaira-Hasegawa Test}
-\usage{
-sh.test(..., x, model = DNAmodel(), B = 100)
-}
-\arguments{
- \item{...}{either a series of objects of class \code{"phylo"}
- separated by commas, or a list containing such objects.}
- \item{x}{a list or a matrix containing the (aligned) DNA sequences.}
- \item{model}{the model to be fitted to each tree (as an object of
- \code{"DNAmodel"}).}
- \item{B}{the number of bootstrap replicates.}
-}
-\description{
- This function computes the Shimodaira--Hasegawa test for a set of
- trees.
-}
-\details{
- The present implementation follows the original formulation of
- Shimodaira and Hasegawa (1999) with the difference that the bootstrap
- resampling is done on the original sequence data rather than the RELL
- method suggested by Shimodaira and Hasegawa.
-}
-\value{
- a numeric vector with the P-value associated with each tree given in
- \code{...}.
-}
-\references{
- Shimodaira, H. and Hasegawa, M. (1999) Multiple comparisons of
- log-likelihoods with applications to phylogenetic
- inference. \emph{Molecular Biology and Evolution}, \bold{16},
- 1114--1116.
-}
-\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
-\seealso{
- \code{\link{mlphylo}}, \code{\link{DNAmodel}}
-}
-\examples{
-data(woodmouse)
-t1 <- nj(dist.dna(woodmouse))
-t2 <- rtree(15, tip.label = t1$tip.label)
-t3 <- rtree(15, tip.label = t1$tip.label)
-### Are the NJ tree and two random tress significantly different?
-\dontrun{sh.test(t1, t2, t3, x = woodmouse, B = 100)}
-}
-\keyword{models}
optional elements (branch lengths, node labels, and root edge) are not
found in the tree.
- If the tree was estimated by maximum likelihood with
- \code{\link{mlphylo}}, a summary of the model fit and the parameter
- estimated is printed.
-
\code{summary} simply prints its results on the standard output and is
not meant for programming.
}
BIRTH.logis <- function(t) log(exp(-a*t) + exp(b))/a + t
BIRTH.step <- function(t)
{
- if (t <= Tcl) return(t*l1)
- Tcl*l1 + (t - Tcl)*l2
+ out <- numeric(length(t))
+ sel <- t <= Tcl
+ if (any(sel)) out[sel] <- t[sel] * l1
+ if (any(!sel)) out[!sel] <- Tcl * l1 + (t[!sel] - Tcl) * l2
+ out
}
data(bird.families)
### fit both models:
yule.time(bird.families, birth.logis)
yule.time(bird.families, birth.logis, BIRTH.logis) # same but faster
\dontrun{yule.time(bird.families, birth.step)} # fails
-yule.time(bird.families, birth.step, BIRTH.step)
+yule.time(bird.families, birth.step, BIRTH.step,
+ opti = "nlminb", start = c(.01, .01, 100))
}
\keyword{models}
+++ /dev/null
-/* mlphylo.c 2008-06-18 */
-
-/* Copyright 2006-2008 Emmanuel Paradis */
-
-/* This file is part of the R-package `ape'. */
-/* See the file ../COPYING for licensing issues. */
-
-#include <R.h>
-#include <Rmath.h>
-#include <R_ext/Applic.h>
-#include <R_ext/Lapack.h>
-
-typedef struct {
- int *n;
- int *s;
- double *w;
- unsigned char *seq;
- double *anc;
-} dna_matrix;
-
-typedef struct {
- int *edge1;
- int *edge2;
- double *el;
-} phylo;
-
-typedef struct {
- int *npart;
- int *partition;
- int *model;
- double *xi;
- double *para;
- int *npara;
- double *alpha;
- int *nalpha;
- int *ncat;
- double *invar;
- int *ninvar;
-} DNAmodel;
-
-typedef struct {
- dna_matrix X;
- phylo PHY;
- DNAmodel MOD;
- double *BF;
-} DNAdata;
-
-typedef struct {
- DNAdata *D; int i;
-} info;
-
-
-void tQ_unbalBF(double *BF, double *P)
-/* This function computes the rate matrix Q multiplied by
- time t in the case of unbalanced base frequencies.
- The arguments are:
- BF: the base frequencies
- P: (input) the matrix of substitution rates
- (output) tQ
- NOTE: P must already be multiplied by t */
-{
- P[1] *= BF[0]; P[2] *= BF[0]; P[3] *= BF[0];
- P[4] *= BF[1]; P[6] *= BF[1]; P[7] *= BF[1];
- P[8] *= BF[2]; P[9] *= BF[2]; P[11] *= BF[2];
- P[12] *= BF[3]; P[13] *= BF[3]; P[14] *= BF[3];
-
- P[0] = -P[4] - P[8] - P[12];
- P[5] = -P[1] - P[9] - P[13];
- P[10] = -P[2] - P[6] - P[14];
- P[15] = -P[3] - P[7] - P[11];
-}
-
-void mat_expo4x4(double *P)
-/* This function computes the exponential of a 4x4 matrix */
-{
- double U[16], vl[4], WR[4], Uinv[16], WI[4], work[32];
- int i, info, ipiv[16], n = 4, lw = 32, ord[4];
- char yes = 'V', no = 'N';
-
- /* The matrix is not symmetric, so we use 'dgeev'. */
- /* We take the real part of the eigenvalues -> WR */
- /* and the right eigenvectors (vr) -> U */
- F77_CALL(dgeev)(&no, &yes, &n, P, &n, WR, WI, vl, &n,
- U, &n, work, &lw, &info);
-
- /* It is not necessary to sort the eigenvalues... */
- /* Copy U -> P */
- for (i = 0; i < 16; i++) P[i] = U[i];
-
- /* For the inversion, we first make Uinv an identity matrix */
- for (i = 1; i < 15; i++) Uinv[i] = 0;
- Uinv[0] = Uinv[5] = Uinv[10] = Uinv[15] = 1;
-
- /* The matrix is not symmetric, so we use 'dgesv'. */
- /* This subroutine puts the result in Uinv (B) */
- /* (P [= U] is erased) */
- F77_CALL(dgesv)(&n, &n, P, &n, ipiv, Uinv, &n, &info);
-
- /* The matrix product of U with the eigenvalues diagonal matrix: */
- for (i = 0; i < 4; i++) U[i] *= exp(WR[0]);
- for (i = 4; i < 8; i++) U[i] *= exp(WR[1]);
- for (i = 8; i < 12; i++) U[i] *= exp(WR[2]);
- for (i = 12; i < 16; i++) U[i] *= exp(WR[3]);
-
- /* The second matrix product with U^-1 */
- P[1] = U[1]*Uinv[0] + U[5]*Uinv[1] + U[9]*Uinv[2] + U[13]*Uinv[3];
- P[2] = U[2]*Uinv[0] + U[6]*Uinv[1] + U[10]*Uinv[2] + U[14]*Uinv[3];
- P[3] = U[3]*Uinv[0] + U[7]*Uinv[1] + U[11]*Uinv[2] + U[15]*Uinv[3];
- P[4] = U[0]*Uinv[4] + U[4]*Uinv[5] + U[8]*Uinv[6] + U[12]*Uinv[7];
- P[6] = U[2]*Uinv[4] + U[6]*Uinv[5] + U[10]*Uinv[6] + U[14]*Uinv[7];
- P[7] = U[3]*Uinv[4] + U[7]*Uinv[5] + U[11]*Uinv[6] + U[15]*Uinv[7];
- P[8] = U[0]*Uinv[8] + U[4]*Uinv[9] + U[8]*Uinv[10] + U[12]*Uinv[11];
- P[9] = U[1]*Uinv[8] + U[5]*Uinv[9] + U[9]*Uinv[10] + U[13]*Uinv[11];
- P[11] = U[3]*Uinv[8] + U[7]*Uinv[9] + U[11]*Uinv[10] + U[15]*Uinv[11];
- P[12] = U[0]*Uinv[12] + U[4]*Uinv[13] + U[8]*Uinv[14] + U[12]*Uinv[15];
- P[13] = U[1]*Uinv[12] + U[5]*Uinv[13] + U[9]*Uinv[14] + U[13]*Uinv[15];
- P[14] = U[2]*Uinv[12] + U[6]*Uinv[13] + U[10]*Uinv[14] + U[14]*Uinv[15];
- P[0] = 1 - P[4] - P[8] - P[12];
- P[5] = 1 - P[1] - P[9] - P[13];
- P[10] = 1 - P[2] - P[6] - P[14];
- P[15] = 1 - P[3] - P[7] - P[11];
-}
-
-void PMAT_JC69(double t, double u, double *P)
-{
- P[1]=P[2]=P[3]=P[4]=P[6]=P[7]=P[8]=P[9]=P[11]=P[12]=P[13]=P[14]=(1 - exp(-4*u*t))/4;
- P[0] = P[5] = P[10] = P[15] = 1 - 3*P[1];
-}
-
-void PMAT_K80(double t, double b, double a, double *P)
-{
- double R, p;
-
- R = a/(2*b);
- p = exp(-2*t/(R + 1));
-
- P[1] = 0.5*(1 - p); /* A -> C */
- P[2] = 0.25 - 0.5*exp(-t*(2*R + 1)/(R + 1)) + 0.25*p; /* A -> G */
- P[0] = P[5] = P[10] = P[15] = 1 - 2*P[1] - P[2];
- P[3] = P[4] = P[6] = P[11] = P[9] = P[12] = P[14] = P[1];
- P[7] = P[8] = P[13] = P[2];
-}
-
-void PMAT_F81(double t, double u, double *BF, double *P)
-{
- double p;
- p = exp(-t*u);
-
- P[0] = p + (1 - p) * BF[0]; /* A->A */
- P[1] = P[9] = P[13] = (1 - p)*BF[1]; /* A->C, G->C, T->C */
- P[2] = P[6] = P[14] = (1 - p)*BF[2]; /* A->G, C->G, T->G */
- P[3] = P[7] = P[11] = (1 - p)*BF[3]; /* A->T, C->T, G->T */
- P[4] = P[8] = P[12] = (1 - p)*BF[0]; /* C->A, G->A, T->A */
- P[5] = p + P[1]; /* C->C */
- P[10] = p + P[2]; /* G->G */
- P[15] = p + P[3]; /* T->T */
-}
-
-void PMAT_F84(double t, double a, double b, double *BF, double *P)
-{
- double pI, pII, B, x, y;
-
- B = exp(-b*t);
- pI = B * (1 - exp(-a*t)); /* probability of at least one event of type I */
- pII = 1 - B; /* probability of at least one event of type II */
- x = pI * (BF[0] + BF[2]);
- y = pI * (BF[1] + BF[3]);
-
- P[12] = P[4] = pII * BF[0]; /* C->A, T->A */
- P[14] = P[6] = pII * BF[2]; /* C->G, T->G */
- P[9] = P[1] = pII * BF[1]; /* A->C, G->C */
- P[2] = x + P[6]; /* A->G */
- P[11] = P[3] = pII * BF[3]; /* A->T, G->T*/
- P[0] = 1 - P[1] - P[2] - P[3]; /* A->A */
- P[7] = y + P[3]; /* C->T */
- P[5] = 1 - P[4] - P[6] - P[7]; /* C->C */
- P[8] = x + P[4]; /* G->A */
- P[10] = 1 - P[8] - P[9] - P[11]; /* G->G */
- P[13] = y + P[1]; /* T->C */
- P[15] = 1 - P[12] - P[13] - P[14]; /* T->T */
-}
-
-void PMAT_HKY85(double t, double a, double b, double *BF, double *P)
-{
- P[2] = P[7] = P[8] = P[13] = t*a;
- P[1] = P[3] = P[4] = P[6] = P[9] = P[11] = P[12] = P[14] = t*b;
- tQ_unbalBF(BF, P);
- mat_expo4x4(P);
-}
-
-void PMAT_T92(double t, double a, double b, double *BF, double *P)
-{
- double theta, A, B1, B2, C, x, y;
-
- theta = BF[1] + BF[2];
- A = (1 - theta)/2;
- B1 = (1 + exp(-t));
- B2 = (1 - exp(-t));
- C = exp(-t * ((a/b + 1)/2));
- x = 0.5 * theta * B1;
- y = (1 - theta) * C;
-
- P[1] = P[6] = P[9] = P[14] = 0.5 * theta * B2; /* A->C, C->G, T->G, G->C */
- P[2] = P[13] = x - theta * C; /* A->G, T->C */
- P[3] = P[4] = P[11] = P[12] = A * B2; /* A->T, C->A, G->T, T->A */
- P[0] = P[15] = 1 - P[1] - P[2] - P[3]; /* A->A, T->T */
- P[5] = P[10] = x + y; /* C->C, G->G */
- P[7] = P[8] = A * B1 - y; /* C->T, G->A */
-}
-
-void PMAT_TN93(double t, double a1, double a2, double b,
- double *BF, double *P)
-{
-
- double A1, A2, B;
-
- A1 = (1 - exp(-a1*t)); /* transitions between purines (A <-> G) */
- A2 = (1 - exp(-a2*t)); /* transitions between pyrimidines (C <-> T) */
- B = (1 - exp(-b*t));
-
- P[1] = B * BF[1]; /* A->C */
- P[2] = A1 * BF[2]; /* A->G */
- P[3] = B * BF[3]; /* A->T */
- P[0] = 1 - P[1] - P[2] - P[3]; /* A->A */
- P[4] = B * BF[0]; /* C->A */
- P[6] = B * BF[2]; /* C->G */
- P[7] = A2 * BF[3] ; /* C->T */
- P[5] = 1 - P[4] - P[6] - P[7]; /* C->C */
- P[8] = A1 * BF[0]; /* G->A */
- P[9] = B * BF[1]; /* G->C */
- P[11] = B * BF[3]; /* G->T */
- P[10] = 1 - P[8] - P[9] - P[11]; /* G->G */
- P[12] = B * BF[0]; /* T->A */
- P[13] = A2 * BF[1]; /* T->C */
- P[14] = B * BF[2]; /* T->G */
- P[15] = 1 - P[12] - P[13] - P[14]; /* T->T */
-}
-
-void PMAT_GTR(double t, double a, double b, double c, double d, double e,
- double f, double *BF, double *P)
-{
- P[1] = P[4] = t*a;
- P[2] = P[8] = t*b;
- P[3] = P[12] = t*c;
- P[6] = P[9] = t*d;
- P[7] = P[13] = t*e;
- P[11] = P[14] = t*f;
- tQ_unbalBF(BF, P);
- mat_expo4x4(P);
-}
-
-#define GET_DNA_PARAMETERS \
- /* get the substitution parameters */ \
- model = *(D->MOD.model); \
- /* If the model is not JC69 or F81: */ \
- if (model != 1 && model != 3) { \
- for(i = 0; i < *(D->MOD.npara); i++) \
- u[i] = D->MOD.para[i]; \
- } \
- /* get the shape parameter and calculate the coefficients */ \
- ncat = *(D->MOD.ncat); \
- if (ncat > 1) { \
- /* use `tmp[0]' to store the mean of the coefficients */ \
- /* in order to rescale them */ \
- tmp[0] = 0.; \
- if (*(D->MOD.nalpha) > 1) alpha = *(D->MOD.alpha); \
- else alpha = D->MOD.alpha[k]; \
- for (j = 0; j < ncat; j++) { \
- coef_gamma[j] = qgamma((0.5 + j)/ncat, alpha, \
- 1/alpha, 1, 0); \
- tmp[0] += coef_gamma[j]; \
- } \
- tmp[0] /= ncat; \
- for (j = 0; j < ncat; j++) \
- coef_gamma[j] /= tmp[0]; \
- } else coef_gamma[0] = 1.; \
- /* get the proportion of invariants */ \
- if (*(D->MOD.ninvar)) { \
- if (*(D->MOD.ninvar) > 1) I = *(D->MOD.invar); \
- else I = D->MOD.invar[k]; \
- } else I = 0.; \
-
-void getSiteLik(int n, int d, int j, int nr, DNAdata *D, double *L)
-{
- int i;
-
- if (d <= n - 1) {
- i = d + j*n;
- memset(L, 0, 4*sizeof(double));
- if (D->X.seq[i] & 128) L[0] = 1; /* if A or else */
- if (D->X.seq[i] & 32) L[1] = 1; /* if C or else */
- if (D->X.seq[i] & 64) L[2] = 1; /* if G or else */
- if (D->X.seq[i] & 16) L[3] = 1; /* if T or else */
- } else {
- i = (d - n) + j*(n - 2);
- L[0] = D->X.anc[i];
- L[1] = D->X.anc[i + nr];
- L[2] = D->X.anc[i + 2*nr];
- L[3] = D->X.anc[i + 3*nr];
- }
-}
-
-#define LOOP_THROUGH_SITES \
- for(j = start; j <= end; j++) { \
- memset(tmp, 0, 4*sizeof(double)); \
- getSiteLik(n, d1, j, nr, D, L1); \
- getSiteLik(n, d2, j, nr, D, L2); \
- for(i = 0; i < ncat; i++) { \
- switch(model) { \
- case 1 : PMAT_JC69(l1, coef_gamma[i], P1); \
- PMAT_JC69(l2, coef_gamma[i], P2); break; \
- case 2 : PMAT_K80(l1, coef_gamma[i], u[0], P1); \
- PMAT_K80(l2, coef_gamma[i], u[0], P2); break; \
- case 3 : PMAT_F81(l1, coef_gamma[i], BF, P1); \
- PMAT_F81(l2, coef_gamma[i], BF, P2); break; \
- case 4 : PMAT_F84(l1, coef_gamma[i], u[0], BF, P1); \
- PMAT_F84(l2, coef_gamma[i], u[0], BF, P2); break; \
- case 5 : PMAT_HKY85(l1, coef_gamma[i], u[0], BF, P1); \
- PMAT_HKY85(l2, coef_gamma[i], u[0], BF, P2); break; \
- case 6 : PMAT_T92(l1, coef_gamma[i], u[0], BF, P1); \
- PMAT_T92(l2, coef_gamma[i], u[0], BF, P2); break; \
- case 7 : PMAT_TN93(l1, coef_gamma[i], u[0], u[1], BF, P1); \
- PMAT_TN93(l2, coef_gamma[i], u[0], u[1], BF, P2); break; \
- case 8 : PMAT_GTR(l1, coef_gamma[i], u[0], u[1], u[2], u[3], u[4], BF, P1); \
- PMAT_GTR(l2, coef_gamma[i], u[0], u[1], u[2], u[3], u[4], BF, P2); break; \
- } \
- tmp[0] += (L1[0]*P1[0] + L1[1]*P1[1] + L1[2]*P1[2] + L1[3]*P1[3]) * \
- (L2[0]*P2[0] + L2[1]*P2[1] + L2[2]*P2[2] + L2[3]*P2[3]); \
- tmp[1] += (L1[0]*P1[4] + L1[1]*P1[5] + L1[2]*P1[6] + L1[3]*P1[7]) * \
- (L2[0]*P2[4] + L2[1]*P2[5] + L2[2]*P2[6] + L2[3]*P2[7]); \
- tmp[2] += (L1[0]*P1[8] + L1[1]*P1[9] + L1[2]*P1[10] + L1[3]*P1[11]) * \
- (L2[0]*P2[8] + L2[1]*P2[9] + L2[2]*P2[10] + L2[3]*P2[11]); \
- tmp[3] += (L1[0]*P1[12] + L1[1]*P1[13] + L1[2]*P1[14] + L1[3]*P1[15]) * \
- (L2[0]*P2[12] + L2[1]*P2[13] + L2[2]*P2[14] + L2[3]*P2[15]); \
- } \
- if (ncat > 1) { \
- tmp[0] /= ncat; \
- tmp[1] /= ncat; \
- tmp[2] /= ncat; \
- tmp[3] /= ncat; \
- } \
- if (D->MOD.ninvar) { \
- V = 1. - I; \
- tmp[0] = V*tmp[0] + I*L1[0]*L2[0]; \
- tmp[1] = V*tmp[1] + I*L1[1]*L2[1]; \
- tmp[2] = V*tmp[2] + I*L1[2]*L2[2]; \
- tmp[3] = V*tmp[3] + I*L1[3]*L2[3]; \
- } \
- ind = anc - n + j*(n - 2); \
- D->X.anc[ind] = tmp[0]; \
- D->X.anc[ind + nr] = tmp[1]; \
- D->X.anc[ind + 2*nr] = tmp[2]; \
- D->X.anc[ind + 3*nr] = tmp[3]; \
- }
-
-void lik_dna_node(DNAdata *D, int ie)
-/*
-This function computes the likelihoods at a node for all
-nucleotides.
-*/
-{
- int d1, d2, anc, n, nr;
- int i, j, k, start, end, ind, i1, i2, ncat, model;
- double tmp[4], L1[4], L2[4], l1, l2, P1[16], P2[16], V, coef_gamma[10], u[6], I, alpha, *BF;
-
- n = *(D->X.n);
- nr = *(D->X.s) * (n - 2);
- BF = D->BF;
-
- d1 = D->PHY.edge2[ie];
- d2 = D->PHY.edge2[ie + 1];
- anc = D->PHY.edge1[ie];
-
- /* change these to use them as indices */
- d1--; d2--; anc--;
-
- l1 = D->PHY.el[ie];
- l2 = D->PHY.el[ie + 1];
-
- for(k = 0; k < *(D->MOD.npart); k++) {
- start = D->MOD.partition[k*2] - 1;
- end = D->MOD.partition[k*2 + 1] - 1;
-
- GET_DNA_PARAMETERS
-
- if (k > 0) {
- l1 *= D->MOD.xi[k - 1];
- l2 *= D->MOD.xi[k - 1];
- }
-
- LOOP_THROUGH_SITES
- }
-} /* lik_dna_node */
-
-void lik_dna_root(DNAdata *D)
-/*
-This function computes the likelihoods at the root for all
-nucleotides.
-*/
-{
- int i, j, k, start, end, ind, ncat, model, d1, d2, n, N, nr;
- double tmp[4], L1[4], L2[4], el, P[16], V, coef_gamma[10], u[6], I, alpha, *BF;
-
- n = *(D->X.n); /* number of tips */
- N = 2*n - 3; /* number of edges */
- nr = *(D->X.s) * (n - 2);
- BF = D->BF;
-
- d1 = D->PHY.edge1[N - 1]; /* it's the root */
- d2 = D->PHY.edge2[N - 1];
-
- /* change these to use them as indices */
- d1--; d2--;
-
- el = D->PHY.el[N - 1];
-
- for(k = 0; k < *(D->MOD.npart); k++) {
- start = D->MOD.partition[k*2] - 1;
- end = D->MOD.partition[k*2 + 1] - 1;
-
- GET_DNA_PARAMETERS
-
- if (k > 0) el *= D->MOD.xi[k - 1];
-
- for(j = start; j <= end; j++) {
- getSiteLik(n, d1, j, nr, D, L1);
- getSiteLik(n, d2, j, nr, D, L2);
- memset(tmp, 0, 4*sizeof(double));
- for(i = 0; i < ncat; i++) {
- switch(model) {
- case 1 : PMAT_JC69(el, coef_gamma[i], P); break;
- case 2 : PMAT_K80(el, coef_gamma[i], u[0], P); break;
- case 3 : PMAT_F81(el, coef_gamma[i], BF, P); break;
- case 4 : PMAT_F84(el, coef_gamma[i], u[0], BF, P); break;
- case 5 : PMAT_HKY85(el, coef_gamma[i], u[0], BF, P); break;
- case 6 : PMAT_T92(el, coef_gamma[i], u[0], BF, P); break;
- case 7 : PMAT_TN93(el, coef_gamma[i], u[0], u[1], BF, P); break;
- case 8 : PMAT_GTR(el, coef_gamma[i], u[0], u[1], u[2], u[3], u[4], BF, P); break;
- }
- tmp[0] += (L2[0]*P[0] + L2[1]*P[1] + L2[2]*P[2] + L2[3]*P[3]);
- tmp[1] += (L2[0]*P[4] + L2[1]*P[5] + L2[2]*P[6] + L2[3]*P[7]);
- tmp[2] += (L2[0]*P[8] + L2[1]*P[9] + L2[2]*P[10] + L2[3]*P[11]);
- tmp[3] += (L2[0]*P[12] + L2[1]*P[13] + L2[2]*P[14] + L2[3]*P[15]);
- }
- if (ncat > 1) {
- tmp[0] /= ncat;
- tmp[1] /= ncat;
- tmp[2] /= ncat;
- tmp[3] /= ncat;
- }
- if (*(D->MOD.ninvar)) {
- V = 1. - I;
- tmp[0] = V*tmp[0] + I*L2[0]*L1[0];
- tmp[1] = V*tmp[1] + I*L2[1]*L1[1];
- tmp[2] = V*tmp[2] + I*L2[2]*L1[2];
- tmp[3] = V*tmp[3] + I*L2[3]*L1[3];
- }
- ind = j*(n - 2);
- D->X.anc[ind] = tmp[0]*L1[0];
- D->X.anc[ind + nr] = tmp[1]*L1[1];
- D->X.anc[ind + 2*nr] = tmp[2]*L1[2];
- D->X.anc[ind + 3*nr] = tmp[3]*L1[3];
- }
- }
-} /* lik_dna_root */
-
-void lik_dna_tree(DNAdata *D, double *loglik)
-{
- int i, j, n, nnode, nsite, nr;
- double tmp;
-
- n = *(D->X.n);
- nnode = n - 2;
- nsite = *(D->X.s);
- nr = nsite*nnode;
-
- /* loop through the tree
- We don't do the root node here, so i goes between 0 and 2n - 5 */
- for(i = 0; i < 2*n - 4; i += 2) lik_dna_node(D, i);
-
- /* We now do the root */
- lik_dna_root(D);
- *loglik = 0.;
- for(j = 0; j < nsite; j++) {
- tmp = 0.;
- for (i = 0; i < 4; i++)
- tmp += D->BF[i] * D->X.anc[j*nnode + i*nr];
- *loglik += D->X.w[j]*log(tmp);
- }
-} /* lik_dna_tree */
-
-double fcn_mlphylo_invar(double I, info *INFO)
-{
- double loglik;
-
- INFO->D->MOD.invar[INFO->i] = I;
- lik_dna_tree(INFO->D, &loglik);
-
- return -loglik;
-}
-
-void mlphylo_invar(int N, DNAdata *D, double *loglik)
-/*
-optimize proportion of invariants
-*/
-{
- int i;
- info INFO, *infptr;
- double I;
-
- infptr = &INFO;
- INFO.D = D;
-
- for(i = 0; i < N; i++) {
- infptr->i = i;
- I = Brent_fmin(0.0, 1.0,
- (double (*)(double, void*)) fcn_mlphylo_invar,
- infptr, 1.e-9);
- D->MOD.invar[i] = I;
- }
-}
-
-double fcn_mlphylo_gamma(double a, info *INFO)
-{
- double loglik;
-
- INFO->D->MOD.alpha[INFO->i] = a;
- lik_dna_tree(INFO->D, &loglik);
-
- return -loglik;
-}
-
-void mlphylo_gamma(int N, DNAdata *D, double *loglik)
-/*
-optimize gamma (ISV) parameters
-*/
-{
- int i;
- info INFO, *infptr;
- double a;
-
- infptr = &INFO;
- INFO.D = D;
-
- for(i = 0; i < N; i++) {
- infptr->i = i;
- a = Brent_fmin(0.0, 1.e4,
- (double (*)(double, void*)) fcn_mlphylo_gamma,
- infptr, 1.e-6);
- D->MOD.alpha[i] = a;
- }
-}
-
-double fcn_mlphylo_para(double p, info *INFO)
-{
- double loglik;
-
- INFO->D->MOD.para[INFO->i] = p;
- lik_dna_tree(INFO->D, &loglik);
-
- return -loglik;
-}
-
-void mlphylo_para(int N, DNAdata *D, double *loglik)
-/*
-optimize the contrast parameter(s) xi
-*/
-{
- int i;
- info INFO, *infptr;
- double p;
-
- infptr = &INFO;
- INFO.D = D;
-
- for(i = 0; i < N; i++) {
- infptr->i = i;
- p = Brent_fmin(0, 1.e3,
- (double (*)(double, void*)) fcn_mlphylo_para,
- infptr, 1.e-6);
- D->MOD.para[i] = p;
- }
-}
-
-double fcn_mlphylo_xi(double XI, info *INFO)
-{
- double loglik;
-
- INFO->D->MOD.xi[INFO->i] = XI;
- lik_dna_tree(INFO->D, &loglik);
-
- return -loglik;
-}
-
-void mlphylo_xi(int N, DNAdata *D, double *loglik)
-/*
-optimize the contrast parameter(s) xi
-*/
-{
- int i;
- info INFO, *infptr;
- double XI;
-
- infptr = &INFO;
- INFO.D = D;
-
- /* In the following, the range of the search algo was changed from */
- /* 0-1000 to 0-100 to avoid infinite looping. (2006-04-15) */
- /* This was changed again to 0-20. (2006-07-17) */
-
- for(i = 0; i < N; i++) {
- infptr->i = i;
- XI = Brent_fmin(0.0, 2.e1,
- (double (*)(double, void*)) fcn_mlphylo_xi,
- infptr, 1.e-4);
- D->MOD.xi[i] = XI;
- }
-}
-
-double fcn_mlphylo_edgelength(double l, info *INFO)
-{
- double loglik;
-
- INFO->D->PHY.el[INFO->i] = l;
- lik_dna_tree(INFO->D, &loglik);
-
- return -loglik;
-}
-
-void mlphylo_edgelength(int N, DNAdata *D, double *loglik)
-/*
-optimize branch lengths
-*/
-{
- int i;
- info INFO, *infptr;
- double l;
-
- infptr = &INFO;
- INFO.D = D;
-
- for(i = 0; i < N; i++) {
- infptr->i = i;
- l = Brent_fmin(0.0, 0.1,
- (double (*)(double, void*)) fcn_mlphylo_edgelength,
- infptr, 1.e-6);
- D->PHY.el[i] = l;
- }
-}
-
-void mlphylo_DNAmodel(int *n, int *s, unsigned char *SEQ, double *ANC,
- double *w, int *edge1, int *edge2,
- double *edge_length, int *npart, int *partition,
- int *model, double *xi, double *para, int *npara,
- double *alpha, int *nalpha, int *ncat,
- double *invar, int *ninvar, double *BF,
- int *search_tree, int *fixed, double *loglik)
-/*
-This function iterates to find the MLEs of the substitution
-paramaters and of the branch lengths for a given tree.
-*/
-{
- DNAdata *D, data;
-
- D = &data;
-
- D->X.n = n;
- D->X.s = s;
- D->X.w = w;
- D->X.seq = SEQ;
- D->X.anc = ANC;
-
- D->PHY.edge1 = edge1;
- D->PHY.edge2 = edge2;
- D->PHY.el = edge_length;
-
- D->MOD.npart = npart;
- D->MOD.partition = partition;
- D->MOD.model = model;
- D->MOD.xi = xi;
- D->MOD.para = para;
- D->MOD.npara = npara;
- D->MOD.alpha = alpha;
- D->MOD.nalpha = nalpha;
- D->MOD.ncat = ncat;
- D->MOD.invar = invar;
- D->MOD.ninvar = ninvar;
-
- D->BF = BF;
-
- lik_dna_tree(D, loglik);
- if (! *fixed) {
- if (*npart > 1) mlphylo_xi(*npart - 1, D, loglik);
- if (*npara) mlphylo_para(*npara, D, loglik);
- if (*nalpha) mlphylo_gamma(*nalpha, D, loglik);
- if (*ninvar) mlphylo_invar(*ninvar, D, loglik);
- lik_dna_tree(D, loglik);
- }
-
-} /* mlphylo_DNAmodel */
-/*
-void jc69(double *P, double *t, double *u)
-{
- PMAT_JC69(*t, *u, P);
-}
-
-void k80(double *P, double *t, double *u)
-{
- PMAT_K80(*t, u[0], u[1], P);
-}
-*/
-/* plot_phylo.c (2006-10-13) */
+/* plot_phylo.c (2009-10-03) */
-/* Copyright 2004-2006 Emmanuel Paradis
+/* Copyright 2004-2009 Emmanuel Paradis
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
}
}
}
+
+void get_single_index_integer(int *x, int *val, int *index)
+{
+ while (x[*index] != *val) (*index)++;
+ *index += 1;
+}
+
+void get_two_index_integer(int *x, int *val, int *index)
+{
+ while (x[index[0]] != *val) index[0]++;
+ index[1] = index[0] + 1;
+ while (x[index[1]] != *val) index[1]++;
+ index[0] += 1;
+ index[1] += 1;
+}
projects / ape.git / commitdiff