[Adephylo-commits] r177 - in pkg: R src

[noreply at r-forge.r-project.org]

Author: jombart
Date: 2012-02-14 20:56:23 +0100 (Tue, 14 Feb 2012)
New Revision: 177

Modified:
   pkg/R/tree.group.R
   pkg/src/misc.c
Log:
First worling version of phylogenetic clustering.


Modified: pkg/R/tree.group.R
===================================================================
--- pkg/R/tree.group.R	2012-02-13 18:10:12 UTC (rev 176)
+++ pkg/R/tree.group.R	2012-02-14 19:56:23 UTC (rev 177)
@@ -1,25 +1,27 @@
 ##############
 ## treeGroup
 ##############
-treeGroup <- function(x, grp, dat, FUN, boot=FALSE, n.boot=10,
+treeGroup <- function(x, grp, dat=NULL, FUN=NULL, boot=FALSE, n.boot=10,
                       n.dens=4096, plot=TRUE, warn.lab=FALSE, ...){
     if(!require(ape)) stop("ape package is required")
     if(!inherits(x,"phylo")) stop("x is not a phylo object")
-    if(warn.lab && !identical(x$tip.label,rownames(dat))) warning("Tip labels in x and rownames of dat differ \nplease make sure the same order is used in x, grp, and dat")
+    if(boot && (is.null(dat) || is.null(FUN))) stop("dat and FUN must be provided for the bootstrap procedure")
+    if(warn.lab && !is.null(dat) && !identical(x$tip.label,rownames(dat))) warning("Tip labels in x and rownames of dat differ \nplease make sure the same order is used in x, grp, and dat")
     grp <- factor(grp)
     K <- length(LEV <- levels(grp))
-    N <- nrow(dat)
+    N <- length(x$tip.label)
     D <- cophenetic.phylo(x)
     THRES <- 1e-320 # densities < THRES will be set to THRES to avoid log(x)=-Inf
 
 
     ## RE-ORDER GRP AND DATA MATRIX AFTER TIP LABELS ##
-    if(is.null(rownames(dat))) rownames(dat) <- x$tip.label
-    if(!all(x$tip.label %in% rownames(dat))) stop("some tips do not have data matching their label")
-    grp <- grp[match(x$tip.label, rownames(dat))] # grp is assumed to be in the same order as 'dat'
-    dat <- dat[x$tip.label,,drop=FALSE]
+    if(!is.null(dat)){
+        if(is.null(rownames(dat))) rownames(dat) <- x$tip.label
+        if(!all(x$tip.label %in% rownames(dat))) stop("some tips do not have data matching their label")
+        grp <- grp[match(x$tip.label, rownames(dat))] # grp is assumed to be in the same order as 'dat'
+        dat <- dat[x$tip.label,,drop=FALSE]
+    }
 
-
     #### AUXILIARY FUNCTIONS ####
     ## FUNCTION TO ESTIMATE A DENSITY AT A SERIES OF POINTS ##
     compute.dens <- function(dens, values){
@@ -115,7 +117,7 @@
     annot <- rep(" ", N)
     annot[as.character(grp)!=as.character(temp)] <- "!"
     groups <- data.frame(observed=grp, inferred=temp, annot=annot)
-    rownames(groups) <- x$tip.labels
+    rownames(groups) <- rownames(prob)
 
 
     ## BUILD / RETURN RESULT ##
@@ -123,7 +125,7 @@
     ## propcorrect.bygroup <- tapply(annot==" ", grp, mean)
     assignability <- mean((apply(prob,1,max)-.5)/.5)
     ##res <- list(prob=prob,groups=groups, mean.correct=propcorrect, prop.correct=propcorrect.bygroup)
-    res <- list(prob=prob,groups=groups, assigndex=assignability, mean.correct=propcorrect)
+    res <- list(prob=prob, groups=groups, assigndex=assignability, mean.correct=propcorrect)
 
     return(res)
 } # end treeGroup

Modified: pkg/src/misc.c
===================================================================
--- pkg/src/misc.c	2012-02-13 18:10:12 UTC (rev 176)
+++ pkg/src/misc.c	2012-02-14 19:56:23 UTC (rev 177)
@@ -2,8 +2,8 @@
 #include <time.h>
 #include <string.h>
 #include <stdlib.h>
+#include <stdio.h>
 
-
 /* interpolate values in a density */
 /* x contais n values for which we want y=f(x) */
 void predict_density(double *densx, double *densy, int *densn, double *x, double *y, int *n){
@@ -14,10 +14,15 @@
 		while(idx < *densn && x[i]>densx[idx]){
 			idx++;
 		}
-		if(idx==0 || idx == *densn){
-			y[i] = densy[idx]/2;
+		if(idx==0){
+			y[i] = densy[idx]/2.0;
+			/* printf("\nx: %.5f, idx: %d, x.chosen: %.5f, %.5f = %.5f / 2.0", x[i], idx, densx[idx], y[i], densy[idx]); */
+		} else if(idx==*densn){
+			y[i] = densy[idx-1]/2.0;
+			/* printf("\nx: %.5f, idx: %d, x.chosen: %.5f, %.5f = %.5f / 2.0", x[i], idx, densx[idx-1], y[i], densy[idx-1]); */
 		} else {
-			y[i] = (densy[idx] + densy[idx+1])/2;
+			y[i] = (densy[idx-1] + densy[idx])/2.0;
+			/* printf("\nx: %.5f, idx: %d, x.before: %.5f, x.after: %.5f, %.5f = %.5f + %.5f / 2.0", x[i], idx, densx[idx-1], densx[idx], y[i], densy[idx-1], densy[idx]); */
 		}
 	}
 }