[Adephylo-commits] r67 - in pkg: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Dec 2 15:50:42 CET 2008
Author: sdray
Date: 2008-12-02 15:50:42 +0100 (Tue, 02 Dec 2008)
New Revision: 67
Added:
pkg/R/orthogram.R
pkg/man/orthogram.Rd
pkg/src/
pkg/src/adesub.c
pkg/src/adesub.h
pkg/src/phylog.c
Log:
add orthogram and related files
Added: pkg/R/orthogram.R
===================================================================
--- pkg/R/orthogram.R (rev 0)
+++ pkg/R/orthogram.R 2008-12-02 14:50:42 UTC (rev 67)
@@ -0,0 +1,186 @@
+"orthogram"<- function (x, orthobas = NULL, neig = NULL, phylog = NULL,
+ nrepet = 999, posinega = 0, tol = 1e-07, cdot = 1.5, cfont.main = 1.5, lwd = 2, nclass, high.scores = 0,alter=c("greater", "less", "two-sided"))
+{
+ if(is.numeric(x)&is.vector(x)){
+ type <- "numeric"
+ } else if(is.factor(x)){
+ type <- "factor"
+ } else if (inherits(x, "dudi")){
+ type <- "dudi"
+ } else {
+ stop("x must be a numeric vector, a factor or a dudi object")
+ }
+ if(type == "dudi") {
+ nobs <- nrow(x$tab)
+ } else {
+ nobs <- length(x)
+ }
+ if (!is.null(neig)) {
+ orthobas <- scores.neig(neig)
+ } else if (!is.null(phylog)) {
+ if (!inherits(phylog, "phylog")) stop ("'phylog' expected with class 'phylog'")
+ orthobas <- phylog$Bscores
+ }
+
+ if (is.null(orthobas)){
+ stop ("'orthobas','neig','phylog' all NULL")
+ }
+
+ if (!inherits(orthobas, "data.frame")) stop ("'orthobas' is not a data.frame")
+ if (nrow(orthobas) != nobs) stop ("non convenient dimensions")
+ if (ncol(orthobas) != (nobs-1)) stop (paste("'orthobas' has",ncol(orthobas),"columns, expected:",nobs-1))
+ vecpro <- as.matrix(orthobas)
+ npro <- ncol(vecpro)
+
+ w <- t(vecpro/nobs)%*%vecpro
+ if (any(abs(diag(w)-1)>tol)) {
+
+ stop("'orthobas' is not orthonormal for uniform weighting")
+ }
+ diag(w) <- 0
+ if ( any( abs(as.numeric(w))>tol) )
+ stop("'orthobas' is not orthogonal for uniform weighting")
+ if (nrepet < 99) nrepet <- 99
+ if (posinega !=0) {
+ if (posinega >= nobs-1) stop ("Non convenient value in 'posinega'")
+ if (posinega <0) stop ("Non convenient value in 'posinega'")
+ }
+ if(type!="dudi"){
+ if (any(is.na(x)))
+ stop("missing value in 'x'")
+ }
+ if(type == "factor"){
+ dudi1 <- dudi.acm(data.frame(x), scannf = F, nf = min(nobs, nlevels(x)))
+ }
+ if(type == "dudi") {
+ if (!all.equal(x$lw, rep(1/nobs, nobs)))
+ stop("not implemented for non-uniform row weights")
+ dudi1 <- redo.dudi(x, newnf = x$rank)
+ if(any(colMeans(dudi1$li)>tol))
+ stop("not implemented for non-centered analysis")
+ }
+
+ if(type == "numeric") {
+ z <- x - mean(x)
+ et <- sqrt(mean(z * z))
+ if ( et <= tol*(max(z)-min(z))) stop ("No variance")
+ z <- z/et
+ w <- .C("VarianceDecompInOrthoBasis",
+ param = as.integer(c(nobs,npro,nrepet,posinega)),
+ observed = as.double(z),
+ vecpro = as.double(vecpro),
+ phylogram = double(npro),
+ phylo95 = double(npro),
+ sig025 = double(npro),
+ sig975 = double(npro),
+ R2Max = double(nrepet+1),
+ SkR2k = double(nrepet+1),
+ Dmax = double(nrepet+1),
+ SCE = double(nrepet+1),
+ ratio = double(nrepet+1),
+ PACKAGE="ade4"
+ )
+ } else {
+ w <- .C("MVarianceDecompInOrthoBasis",
+ param = as.integer(c(nobs,npro,nrepet,posinega)),
+ observed = as.double(as.matrix(dudi1$li)),
+ nvar = as.integer(ncol(dudi1$li)),
+ inertot = as.double(sum(dudi1$eig)),
+ vecpro = as.double(vecpro),
+ phylogram = double(npro),
+ phylo95 = double(npro),
+ sig025 = double(npro),
+ sig975 = double(npro),
+ R2Max = double(nrepet+1),
+ SkR2k = double(nrepet+1),
+ Dmax = double(nrepet+1),
+ SCE = double(nrepet+1),
+ ratio = double(nrepet+1),
+ PACKAGE="adephylo"
+ )
+ }
+ ##return(w$phylogram)
+ ## multiple graphical window (6 graphs)
+ ## 1 pgram
+ ## 2 cumulated pgram
+ ## 3-6 Randomization tests
+
+ def.par <- par(no.readonly = TRUE)
+ on.exit(par(def.par))
+ layout (matrix(c(1,1,2,2,1,1,2,2,3,4,5,6),4,3))
+ par(mar = c(0.1, 0.1, 0.1, 0.1))
+ par(usr = c(0,1,-0.05,1))
+
+
+ ylim <- max(c(w$phylogram, w$phylo95))
+ names(w$phylogram) <- as.character(1:npro)
+ phylocum <- cumsum(w$phylogram)
+ lwd0=2
+ fun <- function (y, last=FALSE) {
+ delta <- (mp[2]-mp[1])/3
+ sel <- 1:(npro - 1)
+ segments(mp[sel]-delta,y[sel],mp[sel]+delta, y[sel],lwd=lwd0)
+ if(last) segments(mp[npro]-delta,y[npro],mp[npro]+delta, y[npro],lwd=lwd0)
+ }
+ sig50 <- (1:npro)/npro
+ y0 <- phylocum - sig50
+ h.obs <- max(y0)
+ x0 <- min(which(y0 == h.obs))
+ par(mar = c(3.1, 2.5, 2.1, 2.1))
+ if(type == "numeric"){
+ z0 <- apply(vecpro, 2, function(x) sum(z * x))
+ mp <- barplot(w$phylogram, col = grey(1 - 0.3 * (sign(z0) > 0)), ylim = c(0, ylim * 1.05))
+ } else {
+ mp <- barplot(w$phylogram, ylim = c(0, ylim * 1.05))
+ }
+ scores.order <- (1:length(w$phylogram))[order(w$phylogram, decreasing=TRUE)[1:high.scores]]
+ fun(w$phylo95,TRUE)
+ abline(h = 1/npro)
+ if (posinega!=0) {
+ verti = (mp[posinega]+mp[posinega+1])/2
+ abline (v=verti, col="red",lwd=1.5)
+ }
+ title(main = "Variance decomposition",font.main=1, cex.main=cfont.main)
+ box()
+ obs0 <- rep(0, npro)
+ names(obs0) <- as.character(1:npro)
+ barplot(obs0, ylim = c(-0.05, 1.05))
+ abline(h=0,col="white")
+ if (posinega!=0) {
+ verti = (mp[posinega]+mp[posinega+1])/2
+ abline (v=verti, col="red",lwd=1.5)
+ }
+
+ title(main = "Cumulative decomposition",font.main=1, cex.main=cfont.main)
+ points(mp, phylocum, pch = 21, cex = cdot, type = "b")
+ segments(mp[1], 1/npro, mp[npro], 1, lty = 1)
+ fun(w$sig975)
+ fun(w$sig025)
+ arrows(mp[x0], sig50[x0], mp[x0], phylocum[x0], ang = 15, le = 0.15,
+ lwd = 2)
+ box()
+ if (missing(nclass)) {
+ nclass <- as.integer (nrepet/25)
+ nclass <- min(c(nclass,40))
+ }
+ plot.randtest (as.randtest (w$R2Max[-1],w$R2Max[1],call=match.call()),main = "R2Max",nclass=nclass)
+ if (posinega !=0) {
+ plot.randtest (as.randtest (w$ratio[-1],w$ratio[1],call=match.call()),main = "Ratio",nclass=nclass)
+ } else {
+ plot.randtest (as.randtest (w$SkR2k[-1],w$SkR2k[1],call=match.call()),main = "SkR2k",nclass=nclass)
+ }
+ plot.randtest (as.randtest (w$Dmax[-1],w$Dmax[1],call=match.call()),main = "DMax",nclass=nclass)
+ plot.randtest (as.randtest (w$SCE[-1],w$SCE[1],call=match.call()),main = "SCE",nclass=nclass)
+
+ w$param <- w$observed <- w$vecpro <- NULL
+ w$phylo95 <- w$sig025 <- w$sig975 <- NULL
+ if (posinega==0) {
+ w <- as.krandtest(obs=c(w$R2Max[1],w$SkR2k[1],w$Dmax[1],w$SCE[1]),sim=cbind(w$R2Max[-1],w$SkR2k[-1],w$Dmax[-1],w$SCE[-1]),names=c("R2Max","SkR2k","Dmax","SCE"),alter=alter,call=match.call())
+ } else {
+ w <- as.krandtest(obs=c(w$R2Max[1],w$SkR2k[1],w$Dmax[1],w$SCE[1],w$ratio[1]),sim=cbind(w$R2Max[-1],w$SkR2k[-1],w$Dmax[-1],w$SCE[-1],w$ratio[-1]),names=c("R2Max","SkR2k","Dmax","SCE","ratio"),alter=alter,call=match.call())
+ }
+
+ if (high.scores != 0)
+ w$scores.order <- scores.order
+ return(w)
+}
Added: pkg/man/orthogram.Rd
===================================================================
--- pkg/man/orthogram.Rd (rev 0)
+++ pkg/man/orthogram.Rd 2008-12-02 14:50:42 UTC (rev 67)
@@ -0,0 +1,100 @@
+\encoding{latin1}
+\name{orthogram}
+\alias{orthogram}
+\title{Orthonormal decomposition of variance}
+\description{
+This function performs the orthonormal decomposition of variance of a quantitative variable on an orthonormal basis. It also returns the results of five non parametric tests associated to the variance decomposition.
+It thus provides tools (graphical displays and test) for analysing phylogenetic, spatial and temporal pattern of one quantitative variable.
+}
+\usage{
+orthogram(x, orthobas = NULL, neig = NULL, phylog = NULL,
+ nrepet = 999, posinega = 0, tol = 1e-07, na.action = c("fail",
+ "mean"), cdot = 1.5, cfont.main = 1.5, lwd = 2, nclass,
+ high.scores = 0,alter=c("greater", "less", "two-sided"))
+}
+\arguments{
+ \item{x}{a numeric vector corresponding to the quantitative variable}
+ \item{orthobas}{an object of class \code{'orthobasis'}}
+ \item{neig}{an object of class \code{'neig'}}
+ \item{phylog}{an object of class \code{'phylog'}}
+ \item{nrepet}{an integer giving the number of permutations}
+ \item{posinega}{a parameter for the ratio test. If posinega > 0, the function computes the ratio test.}
+ \item{tol}{a tolerance threshold for orthonormality condition}
+ \item{na.action}{if 'fail' stops the execution of the current expression when \code{z} contains any missing value. If 'mean' replaces any missing values by mean(\code{z})}
+ \item{cdot}{a character size for points on the cumulative decomposition display}
+ \item{cfont.main}{a character size for titles}
+ \item{lwd}{a character size for dash lines}
+ \item{nclass}{a single number giving the number of cells for the histogram}
+ \item{high.scores}{a single number giving the number of vectors to
+ return. If > 0, the function returns labels of vectors that explains
+ the larger part of variance.}
+ \item{alter}{a character string specifying the alternative hypothesis,
+ must be one of "greater" (default), "less" or "two-sided"}
+}
+\details{
+The function computes the variance decomposition of a quantitative vector x on an orthonormal basis B. The variable is normalized given the uniform weight to eliminate problem of scales.
+It plots the squared correlations \eqn{R^{2}}{R^2} between x and vectors of B (variance decomposition) and the cumulated squared correlations \eqn{SR^{2}}{SR^2} (cumulative decomposition).
+The function also provides five non parametric tests to test the existence of autocorrelation. The tests derive from the five following statistics :
+ \item{R2Max}{=\eqn{\max(R^{2})}{max(R^2)}. It takes high value when a high part of the variability is explained by one score.}
+ \item{SkR2k}{=\eqn{\sum_{i=1}^{n-1}(iR^{2}_i)}{sum_i^(n-1) i*(R^2)_i}. It compares the part of variance explained by internal nodes to the one explained by end nodes.}
+ \item{Dmax}{=\eqn{\max_{m=1,...,n-1}(\sum_{j=1}^{m}R^{2}_j - \frac{m}{n-1})}{max_(m=1,...,n-1)(sum_(j=1)^m(R^2_j) - (m/n-1))}. It examines the accumulation of variance for a sequence of scores.}
+ \item{SCE}{=\eqn{\sum_{m=1}^{n-1} (\sum_{j=1}^{m}R^{2}_j - \frac{m}{n-1})^{2}}{sum_(m=1)^(n-1)(sum_(j=1)^m(R^2_j) - (m/n-1))^2}. It examines also the accumulation of variance for a sequence of scores.}
+ \item{ratio}{depends of the parameter posinega. If posinega > 0, the statistic ratio exists and equals \eqn{\sum_{i=1}^{posinega}R^{2}_i}{sum_i (R^2)_i with i < posinega + 1}. It compares the part of variance explained by internal nodes to the one explained by end nodes when we can define how many vectors correspond to internal nodes.}
+}
+\value{
+If (high.scores = 0), returns an object of class \code{'krandtest'} (randomization tests) corresponding to the five non parametric tests. \cr \cr
+If (high.scores > 0), returns a list containg :
+ \item{w}{: an object of class \code{'krandtest'} (randomization tests)}
+ \item{scores.order}{: a vector which terms give labels of vectors that explain the larger part of variance}
+}
+\references{
+Ollier, S., Chessel, D. and Couteron, P. (2005) Orthonormal Transform to Decompose the Variance of a Life-History Trait across a Phylogenetic Tree. \emph{Biometrics}, \bold{62}, 471--477.
+}
+\author{Sébastien Ollier \email{ollier at biomserv.univ-lyon1.fr} \cr
+Daniel Chessel
+}
+\seealso{\code{\link{gridrowcol}}, \code{\link{orthobasis}}, \code{\link{mld}}}
+\examples{
+# a phylogenetic example
+data(ungulates)
+ung.phy <- newick2phylog(ungulates$tre)
+FemBodyMass <- log(ungulates$tab[,1])
+NeonatBodyMass <- log((ungulates$tab[,2]+ungulates$tab[,3])/2)
+plot(FemBodyMass,NeonatBodyMass, pch = 20, cex = 2)
+abline(lm(NeonatBodyMass~FemBodyMass))
+z <- residuals(lm(NeonatBodyMass~FemBodyMass))
+dotchart.phylog(ung.phy,val = z, clabel.n = 1,
+ labels.n = ung.phy$Blabels, cle = 1.5, cdot = 2)
+table.phylog(ung.phy$Bscores, ung.phy,clabel.n = 1,
+ labels.n = ung.phy$Blabels)
+orthogram(z, ung.phy$Bscores)
+orthogram(z, phyl=ung.phy) # the same thing
+
+# a spatial example
+data(irishdata)
+neig1 <- neig(mat01 = 1*(irishdata$link > 0))
+sco1 <- scores.neig(neig1)
+z <- scalewt(irishdata$tab$cow)
+orthogram(z, sco1)
+
+# a temporal example
+data(arrival)
+w <- orthobasis.circ(24)
+orthogram(arrival$hours, w)
+par(mfrow = c(1,2))
+dotcircle(arrival$hours)
+dotcircle(w[,2])
+par(mfrow = c(1,1))
+
+data(lynx)
+ortho <- orthobasis.line(114)
+orthogram(lynx,ortho)
+attributes(lynx)$tsp
+par(mfrow = c(2,1))
+par(mar = c(4,4,2,2))
+plot.ts(lynx)
+plot(ts(ortho[,23], start = 1821, end = 1934, freq = 1), ylab = "score 23")
+par(mfrow = c(1,1))
+}
+\keyword{spatial}
+\keyword{ts}
Added: pkg/src/adesub.c
===================================================================
--- pkg/src/adesub.c (rev 0)
+++ pkg/src/adesub.c 2008-12-02 14:50:42 UTC (rev 67)
@@ -0,0 +1,1214 @@
+#include <math.h>
+#include <time.h>
+#include <string.h>
+#include <stdlib.h>
+#include "adesub.h"
+
+/***********************************************************************/
+double traceXtdLXq (double **X, double **L, double *d, double *q)
+/* Produit matriciel XtDLXQ avec LX comme lag.matrix */
+{
+ /* Declarations de variables C locales */
+ int j, i, lig, col;
+ double **auxi, **A, trace;
+
+
+
+ /* Allocation memoire pour les variables C locales */
+ lig = X[0][0];
+ col = X[1][0];
+ taballoc(&auxi, lig, col);
+ taballoc(&A, col, col);
+
+
+ /* Calcul de LX */
+ prodmatABC(L, X, auxi);
+
+ /* Calcul de DLX */
+ for (i=1;i<=lig;i++) {
+ for (j=1;j<=col;j++) {
+ auxi[i][j] = auxi[i][j] * d[i];
+ }
+ }
+
+ /* Calcul de XtDLX */
+ prodmatAtBC(X,auxi,A);
+
+ /* Calcul de trace(XtDLXQ) */
+ trace=0;
+ for (i=1;i<=col;i++) {
+ trace = trace + A[i][i] * q[i];
+ }
+
+ /* Libération des réservations locales */
+ freetab (auxi);
+ freetab (A);
+ return(trace);
+}
+
+/***********************************************************************/
+void tabintalloc (int ***tab, int l1, int c1)
+/*--------------------------------------------------
+* Allocation de memoire dynamique pour un tableau
+* d'entiers (l1, c1)
+--------------------------------------------------*/
+{
+ int i, j;
+
+ *tab = (int **) calloc(l1+1, sizeof(int *));
+
+ if ( *tab != NULL) {
+ for (i=0;i<=l1;i++) {
+
+ *(*tab+i)=(int *) calloc(c1+1, sizeof(int));
+ if ( *(*tab+i) == NULL ) {
+ for (j=0;j<i;j++) {
+ free(*(*tab+j));
+ }
+ return;
+ }
+ }
+ } else return;
+ **(*tab) = l1;
+ **(*tab+1) = c1;
+ for (i=1;i<=l1;i++) {
+ for (j=1;j<=c1;j++) {
+ (*tab)[i][j] = 0;
+ }
+ }
+}
+
+/***********************************************************************/
+void freeinttab (int **tab)
+/*--------------------------------------------------
+* Allocation de memoire dynamique pour un tableau
+--------------------------------------------------*/
+{
+ int i, n;
+
+ n = *(*(tab));
+
+ for (i=0;i<=n;i++) {
+ free((char *) *(tab+i) );
+ }
+
+ free((char *) tab);
+}
+
+
+/*********************/
+int dtodelta (double **data, double *pl)
+{
+ /* la matrice de distances d2ij dans data est associee aux poids pl
+ Elle est transformee par aij - ai. -a.j + a..
+ aij = -d2ij/2);*/
+
+ int lig, i, j;
+ double *moy, a0, moytot;
+
+ lig=data[0][0];
+ vecalloc(&moy, lig);
+
+ for (i=1; i<=lig; i++) {
+ for (j=1; j<=lig; j++) data[i][j] = 0.0 - data[i][j] * data[i][j] / 2.0;
+ }
+
+ for (i=1; i<=lig; i++) {
+ a0=0;
+ for (j=1; j<=lig; j++) a0 = a0 + pl[j]*data[i][j];
+ moy[i] = a0;
+ }
+ moytot=0;
+ for (i=1; i<=lig; i++) {
+ moytot = moytot+pl[i]*moy[i];
+ }
+ for (i=1; i<=lig; i++) {
+ for (j=1; j<=lig; j++) data[i][j] = data[i][j] - moy[i] - moy[j] + moytot;
+ }
+ freevec (moy);
+ return 1;
+}
+/***************************/
+void initvec (double *v1, double r)
+/*--------------------------------------------------
+* Initialisation des elements d'un vecteur
+--------------------------------------------------*/
+{
+ int i, c1;
+
+ c1 = v1[0];
+ for (i=1;i<=c1;i++) {
+ v1[i] = r;
+ }
+}
+/**************************/
+double alea (void)
+{
+ double w;
+ w = ((double) rand())/ (double)RAND_MAX;
+ return (w);
+}
+/*************************/
+void aleapermutmat (double **a)
+{
+ /* permute au hasard les lignes du tableau a
+ Manly p. 42 le tableau est modifié */
+ int lig, i,j, col, n, k;
+ double z;
+
+ lig = a[0][0];
+ col = a[1][0];
+ for (i=1; i<=lig-1; i++) {
+ j=lig-i+1;
+ k = (int) (j*alea ()+1);
+ /*k = (int) (j*genrand()+1);*/
+ if (k>j) k=j;
+ for (n=1; n<=col; n++) {
+ z = a[j][n];
+ a[j][n]=a[k][n];
+ a[k][n] = z;
+ }
+ }
+}
+/*************************/
+void aleapermutvec (double *a)
+{
+ /* permute au hasard les ÚlÚments du vecteur a
+ Manly p. 42 Le vecteur est modifié
+ from Knuth 1981 p. 139*/
+ int lig, i,j, k;
+ double z;
+
+ lig = a[0];
+ for (i=1; i<=lig-1; i++) {
+ j=lig-i+1;
+ k = (int) (j*alea()+1);
+ /*k = (int) (j*genrand()+1);*/
+ if (k>j) k=j;
+ z = a[j];
+ a[j]=a[k];
+ a[k] = z;
+ }
+}
+/***********************************************************************/
+void DiagobgComp (int n0, double **w, double *d, int *rang)
+/*--------------------------------------------------
+* Diagonalisation
+* T. FOUCART Analyse factorielle de tableaux multiples,
+* Masson, Paris 1984,185p., p. 62. D'apr?s VPROP et TRIDI,
+* de LEBART et coll.
+--------------------------------------------------*/
+{
+ double *s, epsilon;
+ double a, b, c, x, xp, q, bp, ab, ep, h, t, u , v;
+ double dble;
+ int ni, i, i2, j, k, jk, ijk, ij, l, ix, m, m1, isnou;
+
+ vecalloc(&s, n0);
+ a = 0.000000001;
+ epsilon = 0.0000001;
+ ni = 100;
+ if (n0 == 1) {
+ d[1] = w[1][1];
+ w[1][1] = 1.0;
+ *rang = 1;
+ freevec (s);
+ return;
+ }
+
+ for (i2=2;i2<=n0;i2++) {
+
+ b=0.0;
+ c=0.0;
+ i=n0-i2+2;
+ k=i-1;
+ if (k < 2) goto Et1;
+ for (l=1;l<=k;l++) {
+ c = c + fabs((double) w[i][l]);
+ }
+ if (c != 0.0) goto Et2;
+
+Et1: s[i] = w[i][k];
+ goto Etc;
+
+Et2: for (l=1;l<=k;l++) {
+ x = w[i][l] / c;
+ w[i][l] = x;
+ b = b + x * x;
+ }
+ xp = w[i][k];
+ ix = 1;
+ if (xp < 0.0) ix = -1;
+
+/* q = -sqrt(b) * ix; */
+ dble = b;
+ dble = -sqrt(dble);
+ q = dble * ix;
+
+ s[i] = c * q;
+ b = b - xp * q;
+ w[i][k] = xp - q;
+ xp = 0;
+ for (m=1;m<=k;m++) {
+ w[m][i] = w[i][m] / b / c;
+ q = 0;
+ for (l=1;l<=m;l++) {
+ q = q + w[m][l] * w[i][l];
+ }
+ m1 = m + 1;
+ if (k < m1) goto Et3;
+ for (l=m1;l<=k;l++) {
+ q = q + w[l][m] * w[i][l];
+ }
+
+Et3: s[m] = q / b;
+ xp = xp + s[m] * w[i][m];
+ }
+ bp = xp * 0.5 / b;
+ for (m=1;m<=k;m++) {
+ xp = w[i][m];
+ q = s[m] - bp * xp;
+ s[m] = q;
+ for (l=1;l<=m;l++) {
+ w[m][l] = w[m][l] - xp * s[l] - q * w[i][l];
+ }
+ }
+ for (l=1;l<=k;l++) {
+ w[i][l] = c * w[i][l];
+ }
+
+Etc: d[i] = b;
+ } /* for (i2=2;i2<n0;i2++) */
+
+ s[1] = 0.0;
+ d[1] = 0.0;
+
+ for (i=1;i<=n0;i++) {
+
+ k = i - 1;
+ if (d[i] == 0.0) goto Et4;
+ for (m=1;m<=k;m++) {
+ q = 0.0;
+ for (l=1;l<=k;l++) {
+ q = q + w[i][l] * w[l][m];
+ }
+ for (l=1;l<=k;l++) {
+ w[l][m] = w[l][m] - q * w[l][i];
+ }
+ }
+
+Et4: d[i] = w[i][i];
+ w[i][i] = 1.0;
+ if (k < 1) goto Et5;
+ for (m=1;m<=k;m++) {
+ w[i][m] = 0.0;
+ w[m][i] = 0.0;
+ }
+
+Et5:;
+ }
+
+ for (i=2;i<=n0;i++) {
+ s[i-1] = s[i];
+ }
+ s[n0] = 0.0;
+
+ for (k=1;k<=n0;k++) {
+
+ m = 0;
+
+Et6: for (j=k;j<=n0;j++) {
+ if (j == n0) goto Et7;
+ ab = fabs((double) s[j]);
+ ep = a * (fabs((double) d[j]) + fabs((double) d[j+1]));
+ if (ab < ep) goto Et7;
+ }
+
+Et7: isnou = 1;
+ h = d[k];
+ if (j == k) goto Eta;
+ if (m < ni) goto Etd;
+
+ /*err_message("Error: can't compute matrix eigenvalues");*/
+
+Etd: m = m + 1;
+ q = (d[k+1]-h) * 0.5 / s[k];
+
+/* t = sqrt(q * q + 1.0); */
+ dble = q * q + 1.0;
+ dble = sqrt(dble);
+ t = dble;
+
+ if (q < 0.0) isnou = -1;
+ q = d[j] - h + s[k] / (q + t * isnou);
+ u = 1.0;
+ v = 1.0;
+ h = 0.0;
+ jk = j-k;
+ for (ijk=1;ijk<=jk;ijk++) {
+ i = j - ijk;
+ xp = u * s[i];
+ b = v * s[i];
+ if (fabs((double) xp) < fabs((double) q)) goto Et8;
+ u = xp / q;
+
+/* t = sqrt(u * u + 1); */
+ dble = u * u + 1.0;
+ dble = sqrt(dble);
+ t = dble;
+
+ s[i+1] = q * t;
+ v = 1 / t;
+ u = u * v;
+ goto Et9;
+
+Et8: v = q / xp;
+
+/* t = sqrt(1 + v * v); */
+ dble = 1.0 + v * v;
+ dble = sqrt(dble);
+ t = dble;
+
+ s[i+1] = t * xp;
+ u = 1 / t;
+ v = v * u;
+
+Et9:
+ q = d[i+1] - h;
+ t = (d[i] - q) * u + 2.0 * v * b;
+ h = u * t;
+ d[i+1] = q + h;
+ q = v * t - b;
+ for (l=1;l<=n0;l++) {
+ xp = w[l][i+1];
+ w[l][i+1] = u * w[l][i] + v * xp;
+ w[l][i] = v * w[l][i] - u * xp;
+ }
+ }
+ d[k] = d[k] - h;
+ s[k] = q;
+ s[j] = 0.0;
+
+ goto Et6;
+
+Eta:;
+ } /* for (k=1;k<=n0;k++) */
+
+ for (ij=2;ij<=n0;ij++) {
+
+ i = ij - 1;
+ l = i;
+ h = d[i];
+ for (m=ij;m<=n0;m++) {
+ if (d[m] >= h) {
+ l = m;
+ h = d[m];
+ }
+ }
+ if (l == i) {
+ goto Etb;
+ } else {
+ d[l] = d[i];
+ d[i] = h;
+ }
+ for (m=1;m<=n0;m++) {
+ h = w[m][i];
+ w[m][i] = w[m][l];
+ w[m][l] = h;
+ }
+
+Etb:;
+ } /* for (ij=2;ij<=n0;ij++) */
+
+ *rang = 0;
+ for (i=1;i<=n0;i++) {
+ if (d[i] / d[1] < epsilon) d[i] = 0.0;
+ if (d[i] != 0.0) *rang = *rang + 1;
+ }
+ freevec(s);
+} /* DiagoCompbg */
+/***********************************************************************/
+void freeintvec (int *vec)
+/*--------------------------------------------------
+* liberation de memoire pour un vecteur
+--------------------------------------------------*/
+{
+
+ free((char *) vec);
+
+}
+/***********************************************************************/
+void freetab (double **tab)
+/*--------------------------------------------------
+* Allocation de memoire dynamique pour un tableau (l1, c1)
+--------------------------------------------------*/
+{
+ int i, n;
+
+ n = *(*(tab));
+ for (i=0;i<=n;i++) {
+ free((char *) *(tab+i) );
+ }
+ free((char *) tab);
+}
+/***********************************************************************/
+void freevec (double *vec)
+/*--------------------------------------------------
+* liberation de memoire pour un vecteur
+--------------------------------------------------*/
+{
+ free((char *) vec);
+}
+/***********************************************************************/
+void getpermutation (int *numero, int repet)
+/*----------------------
+* affectation d'une permutation alÚatoire des n premiers entiers
+* dans dans un vecteur d'entiers de dimension n
+* vecintalloc prÚalable exigÚ
+* *numero est un vecteur d'entier
+* repet est un entier qui peut prendre une valeur arbitraire
+* utilise dans le germe du generateur de nb pseudo-aleatoires
+* si on l'incremente dans des appels repetes (e.g. simulation) garantit
+* que deux appels donnent deux resultats distincts (seed=clock+repet)
+------------------------*/
+{
+ int i, n, seed;
+ int *alea;
+
+ n=numero[0];
+ vecintalloc (&alea,n);
+
+ /*-------------
+ * numerotation dans numero
+ -----------*/
+ for (i=1;i<=n;i++) {
+ numero[i]=i;
+ }
+
+ /*-------------
+ * affectation de nombres aleatoires dans alea
+ ----------------*/
+ seed = clock();
+ seed = seed + repet;
+ srand(seed);
+ for (i=1;i<=n;i++) {
+ alea[i]=rand();
+ }
+
+ trirapideint (alea , numero, 1, n);
+ freeintvec (alea);
+}
+/***********************************************************************/
+void matcentrage (double **A, double *poili, char *typ)
+{
+
+ if (strcmp (typ,"nc") == 0) {
+ return;
+ } else if (strcmp (typ,"cm") == 0) {
+ matmodifcm (A, poili);
+ return;
+ } else if (strcmp (typ,"cn") == 0) {
+ matmodifcn (A, poili);
+ return;
+ } else if (strcmp (typ,"cp") == 0) {
+ matmodifcp (A, poili);
+ return;
+ } else if (strcmp (typ,"cs") == 0) {
+ matmodifcs (A, poili);
+ return;
+ } else if (strcmp (typ,"fc") == 0) {
+ matmodiffc (A, poili);
+ return;
+ } else if (strcmp (typ,"fl") == 0) {
+ matmodifcm (A, poili);
+ return;
+ }
+}
+/***********************************************************************/
+ void matcentragehi (double **tab, double *poili, int *index, int *assign)
+{
+/*centrage d'un tableau de hill smith
+tab tableau avec quantitatives et qualitatives disjonctifs complets
+poili vecteur poids lignes
+index indique si chaque variables est quali (1) ou quanti (2)
+assign vecteur entier qui donne l'index de la variable pour chaque colonne
+*/
+
+ int l1,c1,i,j,nquant=0,nqual=0,jqual=1,jquant=1;
+ double **tabqual, **tabquant;
+ l1 = tab[0][0];
+ c1 = tab[1][0];
+ for(j=1;j<=c1;j++){
+ if(index[assign[j]]==1){
+ nqual=nqual+1;
+ }
+ else if (index[assign[j]]==2){
+ nquant=nquant+1;
+ }
+ }
+
+ taballoc(&tabqual,l1,nqual);
+ taballoc(&tabquant,l1,nquant);
+
+ for (j=1;j<=c1;j++){
+ if (index[assign[j]]==1) {
+ for (i=1; i<=l1;i++) {
+ tabqual[i][jqual]=tab[i][j];
+
+ }
+ jqual=jqual+1;
+ } else if (index[assign[j]]==2){
+ for (i=1; i<=l1;i++) {
+ tabquant[i][jquant]=tab[i][j];
+
+ }
+ jquant=jquant+1;
+ }
+ }
+
+
+ matmodifcm (tabqual, poili);
+ matmodifcn (tabquant, poili);
+ jqual=1;
+ jquant=1;
+
+ for (j=1;j<=c1;j++) {
+ if (index[assign[j]]==1) {
+ for (i=1;i<=l1;i++) {
+ tab[i][j] = tabqual[i][jqual];
+ }
+ jqual=jqual+1;
+ }
+ else if (index[assign[j]]==2) {
+ for (i=1;i<=l1;i++) {
+ tab[i][j] = tabquant[i][jquant];
+ }
+ jquant=jquant+1;
+ }
+ }
+ freetab(tabqual);
+ freetab(tabquant);
+ return;
+}
+
+/***********************************************************************/
+void matmodifcm (double **tab, double *poili)
+/*--------------------------------------------------
+* tab est un tableau n lignes, m colonnes
+* disjonctif complet
+* poili est un vecteur n composantes
+* la procedure retourne tab centre par colonne
+* pour la ponderation poili (somme=1)
+* centrage type correspondances multiples
+--------------------------------------------------*/
+{
+ double poid;
+ int i, j, l1, m1;
+ double *poimoda;
+ double x, z;
+
+ l1 = tab[0][0];
+ m1 = tab[1][0];
+ vecalloc(&poimoda, m1);
+
+
+ for (i=1;i<=l1;i++) {
+ poid = poili[i];
+ for (j=1;j<=m1;j++) {
+ poimoda[j] = poimoda[j] + tab[i][j] * poid;
+ }
+ }
+
+ for (j=1;j<=m1;j++) {
+ x = poimoda[j];
+ if (x==0) {
+ for (i=1;i<=l1;i++) tab[i][j] = 0;
+ } else {
+
+ for (i=1;i<=l1;i++) {
+ z = tab[i][j]/x - 1.0;
+ tab[i][j] = z;
+ }
+ }
+ }
+ freevec (poimoda);
+}
+/***********************************************************************/
+void matmodifcn (double **tab, double *poili)
+/*--------------------------------------------------
+* tab est un tableau n lignes, p colonnes
+* poili est un vecteur n composantes
+* la procedure retourne tab norme par colonne
+* pour la ponderation poili (somme=1)
+--------------------------------------------------*/
+{
+ double poid, x, z, y, v2;
+ int i, j, l1, c1;
+ double *moy, *var;
+
+ l1 = tab[0][0];
+ c1 = tab[1][0];
+
+ vecalloc(&moy, c1);
+ vecalloc(&var, c1);
+
+
+/*--------------------------------------------------
+* calcul du tableau centre/norme
+--------------------------------------------------*/
+
+ for (i=1;i<=l1;i++) {
+ poid = poili[i];
+ for (j=1;j<=c1;j++) {
+ moy[j] = moy[j] + tab[i][j] * poid;
+ }
+ }
+
+ for (i=1;i<=l1;i++) {
+ poid=poili[i];
+ for (j=1;j<=c1;j++) {
+ x = tab[i][j] - moy[j];
+ var[j] = var[j] + poid * x * x;
+ }
+ }
+
+ for (j=1;j<=c1;j++) {
+ v2 = var[j];
+ if (v2<=0) v2 = 1;
+ v2 = sqrt(v2);
+ var[j] = v2;
+ }
+
+ for (i=1;i<=c1;i++) {
+ x = moy[i];
+ y = var[i];
+ for (j=1;j<=l1;j++) {
+ z = tab[j][i] - x;
+ z = z / y;
+ tab[j][i] = z;
+ }
+ }
+
+ freevec(moy);
+ freevec(var);
+
+}
+/***********************************************************************/
+void matmodifcs (double **tab, double *poili)
+/*--------------------------------------------------
+* tab est un tableau n lignes, p colonnes
+* poili est un vecteur n composantes
+* la procedure retourne tab standardise par colonne
+* pour la ponderation poili (somme=1)
+--------------------------------------------------*/
+{
+ double poid, x, z, y, v2;
+ int i, j, l1, c1;
+ double *var;
+
+ l1 = tab[0][0];
+ c1 = tab[1][0];
+
+ vecalloc(&var, c1);
+
+
+/*--------------------------------------------------
+* calcul du tableau standardise
+--------------------------------------------------*/
+
+ for (i=1;i<=l1;i++) {
+ poid=poili[i];
+ for (j=1;j<=c1;j++) {
+ x = tab[i][j];
+ var[j] = var[j] + poid * x * x;
+ }
+ }
+
+ for (j=1;j<=c1;j++) {
+ v2 = var[j];
+ if (v2<=0) v2 = 1;
+ v2 = sqrt(v2);
+ var[j] = v2;
+ }
+
+ for (i=1;i<=c1;i++) {
+ y = var[i];
+ for (j=1;j<=l1;j++) {
+ z = tab[j][i];
+ z = z / y;
+ tab[j][i] = z;
+ }
+ }
+ freevec(var);
+}
+/***********************************************************************/
+void matmodifcp (double **tab, double *poili)
+/*--------------------------------------------------
+* tab est un tableau n lignes, p colonnes
+* poili est un vecteur n composantes
+* la procedure retourne tab centre par colonne
+* pour la ponderation poili (somme=1)
+--------------------------------------------------*/
+{
+ double poid;
+ int i, j, l1, c1;
+ double *moy, x, z;
+
+ l1 = tab[0][0];
+ c1 = tab[1][0];
+ vecalloc(&moy, c1);
+
+
+/*--------------------------------------------------
+* calcul du tableau centre
+--------------------------------------------------*/
+
+ for (i=1;i<=l1;i++) {
+ poid = poili[i];
+ for (j=1;j<=c1;j++) {
+ moy[j] = moy[j] + tab[i][j] * poid;
+ }
+ }
+
+
+ for (i=1;i<=c1;i++) {
+ x = moy[i];
+ for (j=1;j<=l1;j++) {
+ z = tab[j][i] - x;
+ tab[j][i] = z;
+ }
+ }
+ freevec(moy);
+}
+/***********************************************************************/
+void matmodiffc (double **tab, double *poili)
+/*--------------------------------------------------
+* tab est un tableau n lignes, m colonnes
+* de nombres positifs ou nuls
+* poili est un vecteur n composantes
+* la procedure retourne tab centre doublement
+* pour la ponderation poili (somme=1)
+* centrage type correspondances simples
+--------------------------------------------------*/
+{
+ double poid;
+ int i, j, l1, m1;
+ double *poimoda;
+ double x, z;
+
+ l1 = tab[0][0];
+ m1 = tab[1][0];
+ vecalloc(&poimoda, m1);
+
+
+ for (i=1;i<=l1;i++) {
+ x = 0;
+ for (j=1;j<=m1;j++) {
+ x = x + tab[i][j];
+ }
+ if (x!=0) {
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/adephylo -r 67
More information about the Adephylo-commits
mailing list