From: paradis Date: Sat, 25 Jun 2011 01:01:05 +0000 (+0000) Subject: cleaning of C files and update on simulation of OU process X-Git-Url: https://git.donarmstrong.com/?p=ape.git;a=commitdiff_plain;h=6ee9e0a4e1e6bbc09187382bfdef57fafe3844c7 cleaning of C files and update on simulation of OU process git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@162 6e262413-ae40-0410-9e79-b911bd7a66b7 --- diff --git a/DESCRIPTION b/DESCRIPTION index 269ff2d..a45a55b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: ape Version: 2.7-3 -Date: 2011-06-21 +Date: 2011-06-22 Title: Analyses of Phylogenetics and Evolution 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, Klaus Schliep, Korbinian Strimmer, Damien de Vienne Maintainer: Emmanuel Paradis diff --git a/NEWS b/NEWS index a5c311c..2a2601e 100644 --- a/NEWS +++ b/NEWS @@ -1,11 +1,22 @@ CHANGES IN APE VERSION 2.7-3 +NEW FEATURES + + o mantel.test() has a new argument 'alternative' which is + "two-sided" by default. Previously, this test was one-tailed + with no possibility to change. + + BUG FIXES o Branch lengths were wrongly updated with bind.tree(, where = , position = 0). (Thanks to Liam Revell for digging this bug out.) + o Simulation of OU process with rTraitCont() did not work correctly. + This now uses formula from Gillespie (1996) reduced to a BM + process when alpha = 0 avoiding division by zero. + CHANGES IN APE VERSION 2.7-2 diff --git a/R/mantel.test.R b/R/mantel.test.R index 3185be3..9704152 100644 --- a/R/mantel.test.R +++ b/R/mantel.test.R @@ -1,4 +1,4 @@ -## mantel.test.R (2006-07-28) +## mantel.test.R (2011-06-22) ## Mantel Test for Similarity of Two Matrices @@ -23,15 +23,21 @@ lower.triang <- function(m) m[col(m) <= row(m)] } -mantel.test <- function (m1, m2, nperm = 1000, graph = FALSE, ...) +mantel.test <- function (m1, m2, nperm = 1000, graph = FALSE, + alternative = "two.sided", ...) { + alternative <- match.arg(alternative, c("two.sided", "less", "greater")) n <- nrow(m1) realz <- mant.zstat(m1, m2) nullstats <- replicate(nperm, mant.zstat(m1, perm.rowscols(m2, n))) - pval <- sum(nullstats > realz)/nperm + pval <- switch(alternative, + "two.sided" = 2 * sum(abs(nullstats) > abs(realz)), + "less" = sum(nullstats < realz), + "greater" = sum(nullstats > realz)) + pval <- pval / nperm if (graph) { plot(density(nullstats), type = "l", ...) abline(v = realz) } - list(z.stat = realz, p = pval) + list(z.stat = realz, p = pval, alternative = alternative) } diff --git a/R/rTrait.R b/R/rTrait.R index 6731c5c..12dc5bf 100644 --- a/R/rTrait.R +++ b/R/rTrait.R @@ -1,4 +1,4 @@ -## rTrait.R (2011-06-16) +## rTrait.R (2011-06-15) ## Trait Evolution @@ -82,7 +82,7 @@ rTraitDisc <- rTraitCont <- function(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0, - ancestor = FALSE, root.value = 0, linear = TRUE, ...) + ancestor = FALSE, root.value = 0, ...) { if (is.null(phy$edge.length)) stop("tree has no branch length") @@ -115,7 +115,6 @@ rTraitCont <- if (length(theta) == 1) theta <- rep(theta, N) else if (length(theta) != N) stop("'theta' must have one or Nedge(phy) elements") - if (!linear) model <- model + 1L } .C("rTraitCont", as.integer(model), as.integer(N), as.integer(anc - 1L), as.integer(des - 1L), el, diff --git a/man/mantel.test.Rd b/man/mantel.test.Rd index 5364952..babb82c 100644 --- a/man/mantel.test.Rd +++ b/man/mantel.test.Rd @@ -2,7 +2,8 @@ \alias{mantel.test} \title{Mantel Test for Similarity of Two Matrices} \usage{ -mantel.test(m1, m2, nperm = 1000, graph = FALSE, ...) +mantel.test(m1, m2, nperm = 1000, graph = FALSE, + alternative = "two.sided", ...) } \arguments{ \item{m1}{a numeric matrix giving a measure of pairwise distances, @@ -12,34 +13,38 @@ mantel.test(m1, m2, nperm = 1000, graph = FALSE, ...) \item{nperm}{the number of times to permute the data.} \item{graph}{a logical indicating whether to produce a summary graph (by default the graph is not plotted).} + \item{alternative}{a character string defining the alternative + hypothesis: \code{"two.sided"} (default), \code{"less"}, + \code{"greater"}, or any unambiguous abbreviation of these.} \item{\dots}{further arguments to be passed to \code{plot()} (to add a title, change the axis labels, and so on).} } \description{ This function computes Mantel's permutation test for similarity of two matrices. It permutes the rows and columns of the two matrices - randomly and calculates a Z-statistic. + randomly and calculates a \eqn{Z}-statistic. } \details{ - The function calculates a Z-statistic for the Mantel test, equal to + The function calculates a \eqn{Z}-statistic for the Mantel test, equal to the sum of the pairwise product of the lower triangles of the permuted matrices, for each permutation of rows and columns. It - compares the permuted distribution with the Z-statistic observed for - the actual data. + compares the permuted distribution with the \eqn{Z}-statistic observed + for the actual data. If \code{graph = TRUE}, the functions plots the density estimate of - the permutation distribution along with the observed Z-statistic as a - vertical line. + the permutation distribution along with the observed \eqn{Z}-statistic + as a vertical line. - The \code{...} argument allows the user to give further options to + The \code{\dots} argument allows the user to give further options to the \code{plot} function: the title main be changed with \code{main=}, the axis labels with \code{xlab =}, and \code{ylab =}, and so on. } \value{ - \item{z.stat}{the Z-statistic (sum of rows*columns of lower triangle) - of the data matrices.} - \item{p}{P-value (quantile of the observed Z-statistic in the - permutation distribution).} + \item{z.stat}{the \eqn{Z}-statistic (sum of rows*columns of lower + triangle) of the data matrices.} + \item{p}{\eqn{P}-value (quantile of the observed \eqn{Z}-statistic in + the permutation distribution).} + \item{alternative}{the alternative hypothesis.} } \references{ Mantel, N. (1967) The detection of disease clustering and a @@ -49,7 +54,8 @@ mantel.test(m1, m2, nperm = 1000, graph = FALSE, ...) Manly, B. F. J. (1986) \emph{Multivariate statistical methods: a primer.} London: Chapman & Hall. } -\author{Original code in S by Ben Bolker \email{bolker@zoo.ufl.edu}, ported +\author{ + Original code in S by Ben Bolker \email{bolker@zoo.ufl.edu}, ported to \R by Julien Claude \email{claude@isem.univ-montp2.fr} } \examples{ diff --git a/man/rTraitCont.Rd b/man/rTraitCont.Rd index 8a72487..f9ce291 100644 --- a/man/rTraitCont.Rd +++ b/man/rTraitCont.Rd @@ -3,7 +3,7 @@ \title{Continuous Character Simulation} \usage{ rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0, - ancestor = FALSE, root.value = 0, linear = TRUE, ...) + ancestor = FALSE, root.value = 0, ...) } \arguments{ \item{phy}{an object of class \code{"phylo"}.} @@ -20,8 +20,6 @@ rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0, values at the nodes as well (by default, only the values at the tips are returned).} \item{root.value}{a numeric giving the value at the root.} - \item{linear}{a logical indicating which parameterisation of the OU - model to use (see details).} \item{\dots}{further arguments passed to \code{model} if it is a function.} } @@ -44,23 +42,10 @@ rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0, \item{\code{"OU"}:}{an Ornstein-Uhlenbeck model is used. The above indexing rule is used for the three parameters \code{sigma}, - \code{alpha}, and \code{theta}. This may be more interesting for the - last one to model varying phenotypic optima. - - By default the following formula is used: - - \deqn{x_{t''} = x_{t'} - \alpha l (x_{t'} - \theta) + \sigma - l \epsilon}{x(t'') = x(t') - alpha l (x(t') - theta) + sigma - l epsilon} - - where \eqn{l (= t'' - t')} is the branch length, and \eqn{\epsilon - \sim N(0, 1)}{\epsilon ~ N(0, 1)}. If \eqn{\alpha > 1}{alpha > 1}, - this may lead to chaotic oscillations. Thus an alternative - parameterisation is used if \code{linear = FALSE}: - - \deqn{x_{t''} = x_{t'} - (1 - exp(-\alpha l)) * (x_{t'} - \theta) + - \sigma l \epsilon}{x(t'') = x(t') - (1 - exp(-alpha l)) * (x(t') - - theta) + sigma l epsilon}} + \code{alpha}, and \code{theta}. This may be interesting for the last + one to model varying phenotypic optima. The exact updating formula + from Gillespie (1996) are used which are reduced to BM formula if + \code{alpha = 0}. \item{A function:}{it must be of the form \code{foo(x, l)} where \code{x} is the trait of the ancestor and \code{l} is the branch @@ -73,6 +58,10 @@ rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0, present, otherwise, ``Node1'', ``Node2'', etc. } \references{ + Gillespie, D. T. (1996) Exact numerical simulation of the + Ornstein-Uhlenbeck process and its integral. \emph{Physical Review E}, + \bold{54}, 2084--2091. + Paradis, E. (2006) \emph{Analyses of Phylogenetics and Evolution with R.} New York: Springer. } diff --git a/src/SPR.c b/src/SPR.c index 82f0a67..e5e204f 100644 --- a/src/SPR.c +++ b/src/SPR.c @@ -249,7 +249,7 @@ void assignUpWeights(edge *etest, node *vtest, node *va, edge *back, node *cprev /*sib is tree C being passed by B*/ /*D is tree below etest*/ double D_AB, D_CD, D_AC, D_BD; - double thisWeight; + /* double thisWeight; */ sib = siblingEdge(etest); left = etest->head->leftEdge; right = etest->head->rightEdge; diff --git a/src/bNNI.c b/src/bNNI.c index 9363ea5..b477b67 100644 --- a/src/bNNI.c +++ b/src/bNNI.c @@ -57,7 +57,7 @@ void weighTree(tree *T) //void bNNI(tree *T, double **avgDistArray, int *count) void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies) { - edge *e, *centerEdge; + edge *e; /* , *centerEdge; */ edge **edgeArray; int *p, *location, *q; int i,j; diff --git a/src/bipartition.c b/src/bipartition.c index aeb6977..a2b34d3 100644 --- a/src/bipartition.c +++ b/src/bipartition.c @@ -1,6 +1,6 @@ -/* bipartition.c 2007-06-29 */ +/* bipartition.c 2011-06-23 */ -/* Copyright 2005-2007 Emmanuel Paradis, and 2007 R Development Core Team */ +/* Copyright 2005-2011 Emmanuel Paradis, and 2007 R Development Core Team */ /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -143,7 +143,7 @@ int SameClade(SEXP clade1, SEXP clade2) SEXP prop_part(SEXP TREES, SEXP nbtree, SEXP keep_partitions) { - int i, j, k, l, KeepPartition, Ntree, Ntip, Nnode, Npart, NpartCurrent, *no; + int i, j, k, KeepPartition, Ntree, Ntip, Nnode, Npart, NpartCurrent, *no; SEXP bp, ans, nbtip, nbnode, number; PROTECT(nbtree = coerceVector(nbtree, INTSXP)); diff --git a/src/delta_plot.c b/src/delta_plot.c index 66fadac..73afbd4 100644 --- a/src/delta_plot.c +++ b/src/delta_plot.c @@ -1,6 +1,6 @@ -/* delta_plot.c 2010-01-12 */ +/* delta_plot.c 2011-06-23 */ -/* Copyright 2010 Emmanuel Paradis +/* Copyright 2010-2011 Emmanuel Paradis /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -12,7 +12,7 @@ void delta_plot(double *D, int *size, int *nbins, int *counts, double *deltabar) int x, y, u, v; /* same notation than in Holland et al. 2002 */ int n = *size, nb = *nbins; double dxy, dxu, dxv, dyu, dyv, duv, A, B, C, delta; - int xy, xu, xv, yu, yv, uv, i, k; + int xy, xu, xv, yu, yv, uv, i; for (x = 0; x < n - 3; x++) { xy = x*n - x*(x + 1)/2; /* do NOT factorize */ diff --git a/src/dist_dna.c b/src/dist_dna.c index 2eb88a6..210479e 100644 --- a/src/dist_dna.c +++ b/src/dist_dna.c @@ -1,4 +1,4 @@ -/* dist_dna.c 2011-02-18 */ +/* dist_dna.c 2011-06-23 */ /* Copyright 2005-2011 Emmanuel Paradis @@ -535,8 +535,8 @@ void distDNA_TN93(unsigned char *x, int *n, int *s, double *d, double *BF, int *variance, double *var, int *gamma, double *alpha) { - int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2; - double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b; + int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2; + double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b; L = *s; @@ -559,8 +559,8 @@ void distDNA_TN93_pairdel(unsigned char *x, int *n, int *s, double *d, double *BF, int *variance, double *var, int *gamma, double *alpha) { - int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2; - double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b; + int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2; + double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b; PREPARE_BF_TN93 @@ -812,8 +812,8 @@ void distDNA_BH87(unsigned char *x, int *n, int *s, double *d, currently not available. */ { - int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4], ndim = 4, info, ipiv[16]; - double P12[16], P21[16], U[16]; + int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4]; + double P12[16], P21[16]; for (i1 = 1; i1 < *n; i1++) { for (i2 = i1 + 1; i2 <= *n; i2++) { diff --git a/src/mat_expo.c b/src/mat_expo.c index a5ae2db..a2038b2 100644 --- a/src/mat_expo.c +++ b/src/mat_expo.c @@ -1,6 +1,6 @@ -/* matexpo.c 2007-10-08 */ +/* matexpo.c 2011-06-23 */ -/* Copyright 2007 Emmanuel Paradis +/* Copyright 2007-2011 Emmanuel Paradis /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -12,7 +12,7 @@ void mat_expo(double *P, int *nr) /* This function computes the exponential of a nr x nr matrix */ { double *U, *vl, *WR, *Uinv, *WI, *work; - int i, j, k, l, info, *ipiv, n = *nr, nc = n*n, lw = nc << 1, *ord; + int i, j, k, l, info, *ipiv, n = *nr, nc = n*n, lw = nc << 1; char yes = 'V', no = 'N'; U = (double *)R_alloc(nc, sizeof(double)); @@ -23,7 +23,6 @@ void mat_expo(double *P, int *nr) work = (double *)R_alloc(lw, sizeof(double)); ipiv = (int *)R_alloc(nc, sizeof(int)); - ord = (int *)R_alloc(n, sizeof(int)); /* The matrix is not symmetric, so we use 'dgeev'. We take the real part of the eigenvalues -> WR diff --git a/src/plot_phylo.c b/src/plot_phylo.c index fec2601..d8666a2 100644 --- a/src/plot_phylo.c +++ b/src/plot_phylo.c @@ -1,6 +1,6 @@ -/* plot_phylo.c (2009-10-03) */ +/* plot_phylo.c (2011-06-23) */ -/* Copyright 2004-2009 Emmanuel Paradis +/* Copyright 2004-2011 Emmanuel Paradis /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -37,12 +37,11 @@ void node_depth(int *ntip, int *nnode, int *edge1, int *edge2, void node_height(int *ntip, int *nnode, int *edge1, int *edge2, int *nedge, double *yy) { - int i, k, n; + int i, n; double S; /* The coordinates of the tips have been already computed */ - k = 1; S = 0; n = 0; for (i = 0; i < *nedge; i++) { @@ -59,14 +58,13 @@ void node_height(int *ntip, int *nnode, int *edge1, int *edge2, void node_height_clado(int *ntip, int *nnode, int *edge1, int *edge2, int *nedge, double *xx, double *yy) { - int i, k, n; + int i, n; double S; node_depth(ntip, nnode, edge1, edge2, nedge, xx); /* The coordinates of the tips have been already computed */ - k = 1; S = 0; n = 0; for (i = 0; i < *nedge; i++) { diff --git a/src/rTrait.c b/src/rTrait.c index 11cc63e..0e543a6 100644 --- a/src/rTrait.c +++ b/src/rTrait.c @@ -1,6 +1,6 @@ -/* rTrait.c 2010-07-26 */ +/* rTrait.c 2011-06-25 */ -/* Copyright 2010 Emmanuel Paradis */ +/* Copyright 2010-2011 Emmanuel Paradis */ /* This file is part of the R-package `ape'. */ /* See the file ../COPYING for licensing issues. */ @@ -21,14 +21,16 @@ void rTraitCont(int *model, int *Nedge, int *edge1, int *edge2, double *el, } break; case 2 : for (i = *Nedge - 1; i >= 0; i--) { + if (alpha[i]) { + alphaT = alpha[i] * el[i]; + M = exp(-alphaT); + S = sigma[i] * sqrt((1 - exp(-2 * alphaT))/(2 * alpha[i])); + } else { + M = 1; + S = sqrt(el[i]) * sigma[i]; + } GetRNGstate(); - x[edge2[i]] = x[edge1[i]] - alpha[i]*el[i]*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand(); - PutRNGstate(); - } - break; - case 3 : for (i = *Nedge - 1; i >= 0; i--) { - GetRNGstate(); - x[edge2[i]] = x[edge1[i]] - (1 - exp(alpha[i]*el[i]))*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand(); + x[edge2[i]] = x[edge1[i]] * M + theta[i] * (1 - M) + S * norm_rand(); PutRNGstate(); } break; diff --git a/src/tree_build.c b/src/tree_build.c index eb47337..5ac3f86 100644 --- a/src/tree_build.c +++ b/src/tree_build.c @@ -1,4 +1,4 @@ -/* tree_build.c 2011-02-28 */ +/* tree_build.c 2011-06-23 */ /* Copyright 2008-2011 Emmanuel Paradis */ @@ -42,7 +42,7 @@ void decode_terminal_edge_token(const char *x, int a, int b, int *node, double * void decode_internal_edge(const char *x, int a, int b, char *lab, double *w) { - int i, k, co = a; + int co = a; char *endstr, str[100]; while (x[co] != ':') co++;