[adegenet-commits] r813 - pkg/R

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

Author: jombart
Date: 2011-02-21 14:46:11 +0100 (Mon, 21 Feb 2011)
New Revision: 813

Modified:
   pkg/R/glFunctions.R
Log:
- glPca does not give the good results
- however, glDotProd gives the good matrix of dot products, with/without centring/scaling.


Modified: pkg/R/glFunctions.R
===================================================================
--- pkg/R/glFunctions.R	2011-02-18 18:15:57 UTC (rev 812)
+++ pkg/R/glFunctions.R	2011-02-21 13:46:11 UTC (rev 813)
@@ -208,6 +208,22 @@
                   useC=TRUE, multicore=require("multicore"), n.cores=NULL){
     if(!inherits(x, "genlight")) stop("x is not a genlight object")
 
+    ## COMPUTE MEANS AND VARIANCES ##
+    if(center) {
+        vecMeans <- glMean(x, alleleAsUnit=FALSE)
+        if(any(is.na(vecMeans))) stop("NAs detected in the vector of means")
+    }
+
+    if(scale){
+        vecVar <- glVar(x, alleleAsUnit=FALSE)
+        if(any(is.na(vecVar))) stop("NAs detected in the vector of variances")
+    }
+
+
+    myPloidy <- ploidy(x)
+
+
+    ## == if non-C code is used ==
     if(!useC){
         if(multicore && !require(multicore)) stop("multicore package requested but not installed")
         if(multicore && is.null(n.cores)){
@@ -215,23 +231,9 @@
         }
 
 
-        ## COMPUTE MEANS AND VARIANCES ##
-        if(center) {
-            vecMeans <- glMean(x, alleleAsUnit=FALSE)
-            if(any(is.na(vecMeans))) stop("NAs detected in the vector of means")
-        }
-
-        if(scale){
-            vecVar <- glVar(x, alleleAsUnit=FALSE)
-            if(any(is.na(vecVar))) stop("NAs detected in the vector of variances")
-        }
-
-
         ## COMPUTE DOT PRODUCTS BETWEEN GENOTYPES ##
         ## to be fast, a particular function is defined for each case of centring/scaling
 
-        myPloidy <- ploidy(x)
-
         ## NO CENTRING / NO SCALING
         if(!center & !scale){
             dotProd <- function(a,b, ploid.a, ploid.b){ # a and b are two SNPbin objects
@@ -305,7 +307,7 @@
         }
         diag(allProd) <- temp
     } else { # === use C computations ====
-        glDotProd(x, center=center, scale=scale, alleleAsUnit=alleleAsUnit)
+        allProd <- glDotProd(x, center=center, scale=scale, alleleAsUnit=alleleAsUnit)
     }
 
 
@@ -366,11 +368,19 @@
 
     ## FORMAT OUTPUT ##
     colnames(res$scores) <- paste("PC", 1:nf, sep="")
-    rownames(res$scores) <- indNames(x)
+    if(!is.null(indNames(x))){
+        rownames(res$scores) <- indNames(x)
+    } else {
+         rownames(res$scores) <- 1:nInd(x)
+    }
 
     if(!is.null(res$loadings)){
         colnames(res$loadings) <- paste("Axis", 1:nf, sep="")
-        rownames(res$loadings) <- paste(locNames(x),alleles(x), sep=".")
+        if(!is.null(locNames(x)) & !is.null(alleles(x))){
+            rownames(res$loadings) <- paste(locNames(x),alleles(x), sep=".")
+        } else {
+            rownames(res$loadings) <- 1:nLoc(x)
+        }
     }
 
     res$call <- match.call()
@@ -500,7 +510,29 @@
 
 
 
+#### TESTING DOT PRODUCTS ####
+M <- matrix(sample(c(0,1), 100*1e3, replace=TRUE), nrow=100)
+rownames(M) <- paste("ind", 1:100)
+x <- new("genlight",M)
 
+## not centred, not scaled
+res1 <- glDotProd(x,alleleAsUnit=FALSE, center=FALSE, scale=FALSE)
+res2 <- M %*% t(M)
+all.equal(res1,res2) # must be TRUE
+
+##  centred, not scaled
+res1 <- glDotProd(x,alleleAsUnit=FALSE, center=TRUE, scale=FALSE)
+M <- scalewt(M,center=TRUE,scale=FALSE)
+res2 <- M %*% t(M)
+all.equal(res1,res2) # must be TRUE
+
+##  centred, scaled
+res1 <- glDotProd(x,alleleAsUnit=FALSE, center=TRUE, scale=TRUE)
+M <- scalewt(M,center=TRUE,scale=TRUE)
+res2 <- M %*% t(M)
+all.equal(res1,res2) # must be TRUE
+
+
 #### TESTING PCA ####
 ## M <- matrix(sample(c(0,1), 1e5, replace=TRUE), nrow=100)
 ## rownames(M) <- paste("ind", 1:100)
@@ -534,10 +566,10 @@
 
 ## ## LARGE SCALE TEST ##
 ## ## perform glPca
-## M <- matrix(sample(c(0,1), 100*1e6, replace=TRUE), nrow=100)
+## M <- matrix(sample(c(0,1), 100*1e5, replace=TRUE), nrow=100)
 ## x <- new("genlight",M)
-## system.time(toto <- glPca(x, nf=4, n.core=6))
 ## system.time(titi <- dudi.pca(M,center=TRUE,scale=FALSE, scannf=FALSE, nf=4))
+## system.time(toto <- glPca(x, ,center=TRUE,scale=FALSE, useC=TRUE, nf=4))