[adegenet-commits] r821 - in pkg: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 24 16:29:50 CET 2011
Author: jombart
Date: 2011-02-24 16:29:50 +0100 (Thu, 24 Feb 2011)
New Revision: 821
Modified:
pkg/R/glFunctions.R
pkg/src/GLfunctions.c
pkg/src/snpbin.c
Log:
Coded GLsum in C...
but is slower than the R version (!) - will prob. be interesting for many individuals though.
All tested, all works.
Modified: pkg/R/glFunctions.R
===================================================================
--- pkg/R/glFunctions.R 2011-02-24 13:24:11 UTC (rev 820)
+++ pkg/R/glFunctions.R 2011-02-24 15:29:50 UTC (rev 821)
@@ -5,27 +5,47 @@
## compute col sums
## removing NAs
##
-glSum <- function(x, alleleAsUnit=TRUE){
+glSum <- function(x, alleleAsUnit=TRUE, useC=FALSE){
if(!inherits(x, "genlight")) stop("x is not a genlight object")
- ## DEFAULT, VECTOR-WISE PROCEDURE ##
- ## use ploidy (sum absolute frequencies)
- if(alleleAsUnit){
- res <- integer(nLoc(x))
- for(e in x at gen){
- temp <- as.integer(e)
- temp[is.na(temp)] <- 0L
- res <- res + temp
+ if(useC){
+ ## use ploidy (sum absolute frequencies)
+ if(alleleAsUnit){
+ vecbyte <- unlist(lapply(x at gen, function(e) e$snp))
+ nbVec <- sapply(x at gen, function(e) length(e$snp))
+ nbNa <- sapply(NA.posi(x), length)
+ naPosi <- unlist(NA.posi(x))
+ res <- .C("GLsumInt", vecbyte, nbVec, length(x at gen[[1]]@snp[[1]]), nbNa, naPosi, nInd(x), nLoc(x), ploidy(x),
+ integer(nLoc(x)), PACKAGE="adegenet")[[9]]
+ } else {
+ ## sum relative frequencies
+ vecbyte <- unlist(lapply(x at gen, function(e) e$snp))
+ nbVec <- sapply(x at gen, function(e) length(e$snp))
+ nbNa <- sapply(NA.posi(x), length)
+ naPosi <- unlist(NA.posi(x))
+ res <- .C("GLsumFreq", vecbyte, nbVec, length(x at gen[[1]]@snp[[1]]), nbNa, naPosi, nInd(x), nLoc(x), ploidy(x),
+ double(nLoc(x)), PACKAGE="adegenet")[[9]]
}
} else {
- ## sum relative frequencies
- res <- numeric(nLoc(x))
- myPloidy <- ploidy(x)
- for(i in 1:nInd(x)){
- temp <- as.integer(x at gen[[i]]) / myPloidy[i]
- temp[is.na(temp)] <- 0
- res <- res + temp
+ ## use ploidy (sum absolute frequencies)
+ if(alleleAsUnit){
+ res <- integer(nLoc(x))
+ for(e in x at gen){
+ temp <- as.integer(e)
+ temp[is.na(temp)] <- 0L
+ res <- res + temp
+ }
+ } else {
+ ## sum relative frequencies
+ res <- numeric(nLoc(x))
+ myPloidy <- ploidy(x)
+ for(i in 1:nInd(x)){
+ temp <- as.integer(x at gen[[i]]) / myPloidy[i]
+ temp[is.na(temp)] <- 0
+ res <- res + temp
+ }
}
+
}
names(res) <- locNames(x)
return(res)
@@ -619,3 +639,21 @@
## names(t) <- rownames(res$by.total)
## par(mar=c(7,4,4,2))
## barplot(t,las=3, cex.names=.7)
+
+
+## test GLsum:
+## library(adegenet)
+## x <- new("genlight", lapply(1:50, function(i) sample(c(0,1,NA), 1e5, prob=c(.5, .49, .01), replace=TRUE)))
+## res1 <- glSum(x, useC=FALSE)
+## res2 <- glSum(x, useC=TRUE)
+## res3 <- apply(as.matrix(x),2,sum,na.rm=TRUE)
+## all(res1==res3) # must be TRUE
+## all(res2==res3) # must be TRUE
+
+## library(adegenet)
+## x <- new("genlight", lapply(1:50, function(i) sample(c(0,1,2,NA), 1e5, prob=c(.5, .40, .09, .01), replace=TRUE)))
+## res1 <- glSum(x, alleleAsUnit=FALSE, useC=FALSE)
+## res2 <- glSum(x, alleleAsUnit=FALSE, useC=TRUE)
+## res3 <- apply(as.matrix(x)/ploidy(x),2,sum,na.rm=TRUE)
+## all.equal(res1,res3)
+## all.equal(res2,res3)
Modified: pkg/src/GLfunctions.c
===================================================================
--- pkg/src/GLfunctions.c 2011-02-24 13:24:11 UTC (rev 820)
+++ pkg/src/GLfunctions.c 2011-02-24 15:29:50 UTC (rev 821)
@@ -83,21 +83,21 @@
int i, j;
int *vecIntTemp;
vecIntTemp = (int *) calloc(*nloc, sizeof(int));
-
+
/* set res to zeros */
- for(j=0;j< *nloc;j++){
- res[j] = 0;
- }
+ /* for(j=0;j< *nloc;j++){ */
+ /* res[j] = 0; */
+ /* } */
/* Internal C representation of the genlight object */
dat = genlightTogenlightC(gen, nbvecperind, byteveclength, nbnaperind, naposi, nind, nloc, ploidy);
/* === working on frequencies === */
/* Lower triangle - without the diagonal */
- for(i=0; i < (*nind-1); i++){ /* for all individuals*/
+ for(i=0; i < (*nind); i++){ /* for all individuals*/
/* conversion to integers of current indiv */
snpbin2intvec(&(dat.x[i]), vecIntTemp);
-
+
for(j=0; j < *nloc; j++){ /* for all loci */
if(!snpbin_isna(&(dat.x[i]), j)) res[j] += vecIntTemp[j];
}
@@ -108,6 +108,35 @@
+void GLsumFreq(unsigned char *gen, int *nbvecperind, int *byteveclength, int *nbnaperind, int *naposi, int *nind, int *nloc, int *ploidy, double *res){
+ struct genlightC dat;
+ int i, j;
+ double *vecFreqTemp;
+ vecFreqTemp = (double *) calloc(*nloc, sizeof(double));
+
+ /* set res to zeros */
+ /* for(j=0;j< *nloc;j++){ */
+ /* res[j] = 0.0; */
+ /* } */
+
+ /* Internal C representation of the genlight object */
+ dat = genlightTogenlightC(gen, nbvecperind, byteveclength, nbnaperind, naposi, nind, nloc, ploidy);
+
+ /* === working on frequencies === */
+ /* Lower triangle - without the diagonal */
+ for(i=0; i < (*nind); i++){ /* for all individuals*/
+ /* conversion to frequencies of current indiv */
+ snpbin2freq(&(dat.x[i]), vecFreqTemp);
+
+ for(j=0; j < *nloc; j++){ /* for all loci */
+ if(!snpbin_isna(&(dat.x[i]), j)) res[j] += vecFreqTemp[j];
+ }
+ }
+}
+
+
+
+
/* TESTING in R */
/*
Modified: pkg/src/snpbin.c
===================================================================
--- pkg/src/snpbin.c 2011-02-24 13:24:11 UTC (rev 820)
+++ pkg/src/snpbin.c 2011-02-24 15:29:50 UTC (rev 821)
@@ -182,7 +182,12 @@
temp = (int *) calloc(8, sizeof(int));
+ /* set result vector to 0 */
+ for(i=0; i < (*veclength) * 8; i++){
+ vecres[i]=0;
+ }
+ /* build output */
for(k=0;k<*nbvec;k++){ /* for all input vector */
idres = 0;
for(i=0;i<*veclength;i++){ /* for one input vector */
@@ -205,6 +210,10 @@
double *temp;
temp = (double *) calloc(8, sizeof(double));
+ /* set result vector to 0 */
+ for(i=0; i < (*veclength) * 8; i++){
+ vecres[i]=0.0;
+ }
for(k=0;k<*nbvec;k++){ /* for all input vector */
idres = 0;
More information about the adegenet-commits
mailing list