Dear mailing list, I have a C function that gives me a wrong result when I run it under Windows 7. This is the code under Linux (RHEL5):> library(phenoTest) > data(epheno) > sign <- sample(featureNames(epheno))[1:20] > score <- getFc(epheno)[,1] > head(score)1007_s_at 1053_at 117_at 121_at 1255_g_at 1294_at -1.183019 1.113544 1.186186 -1.034779 -1.044456 -1.023471> s <- which(names(score) %in% sign) > es.c <- .Call('getEs',score,s,PACKAGE='phenoTest') > head(es.c)[1] -0.001020408 -0.002040816 -0.003061224 -0.004081633 -0.005102041 [6] -0.006122449> es.c <- .Call('getEs',score,s,PACKAGE='phenoTest') > head(es.c)[1] -0.001020408 -0.002040816 -0.003061224 -0.004081633 -0.005102041 [6] -0.006122449> sessionInfo()R version 2.13.0 (2011-04-13) Platform: x86_64-redhat-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] phenoTest_1.1.1 RSQLite_0.8-4 DBI_0.2-5 [4] Heatplus_1.22.0 annotate_1.30.1 AnnotationDbi_1.14.1 [7] Biobase_2.12.2 loaded via a namespace (and not attached): [1] affyio_1.20.0 biomaRt_2.8.1 Biostrings_2.21.6 [4] Category_2.18.0 cluster_1.13.3 gdata_2.7.1 [7] genefilter_1.34.0 gplots_2.10.1 graph_1.30.0 [10] grid_2.13.0 GSEABase_1.14.0 gtools_2.6.2 [13] hgu133a.db_2.5.0 Hmisc_3.8-3 hopach_2.12.0 [16] IRanges_1.11.11 lattice_0.19-23 limma_3.8.3 [19] Matrix_0.9996875-3 mgcv_1.7-5 nlme_3.1-100 [22] oligoClasses_1.14.0 RBGL_1.22.0 RCurl_1.6-10 [25] SNPchip_1.16.0 splines_2.13.0 survival_2.36-5 [28] tools_2.13.0 XML_3.4-0 xtable_1.5-6 As you see es.c is correct. I checked it doing the same computation with R. It also runs without problems under Mac. I run valgrind on the same piece of code and got no errors. This is the same piece of code under Windows 7:> library(phenoTest) > data(epheno) > sign <- sample(featureNames(epheno))[1:20] > score <- getFc(epheno)[,1] > head(score)1007_s_at 1053_at 117_at 121_at 1255_g_at 1294_at -1.183019 1.113544 1.186186 -1.034779 -1.044456 -1.023471> s <- which(names(score) %in% sign) > es.c <- .Call('getEs',score,s,PACKAGE='phenoTest') > head(es.c)[1] 1.447208e+215 1.447208e+215 1.447208e+215 1.447208e+215 1.447208e+215 1.447208e+215> es.c <- .Call('getEs',score,s,PACKAGE='phenoTest') > head(es.c)[1] 3.176615e+170 3.176615e+170 3.176615e+170 3.176615e+170 3.176615e+170 3.176615e+170> sessionInfo()R version 2.14.0 alpha (2011-10-13 r57240) Platform: i386-pc-mingw32/i386 (32-bit) locale: [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252 [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C [5] LC_TIME=Spanish_Spain.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] phenoTest_1.1.1 RSQLite_0.10.0 DBI_0.2-5 [4] Heatplus_1.99.0 annotate_1.31.1 AnnotationDbi_1.15.36 [7] Biobase_2.13.10 loaded via a namespace (and not attached): [1] affyio_1.21.2 biomaRt_2.9.3 Biostrings_2.21.11 [4] Category_2.19.1 cluster_1.14.0 gdata_2.8.2 [7] genefilter_1.35.0 gplots_2.10.1 graph_1.31.2 [10] grid_2.14.0 GSEABase_1.15.0 gtools_2.6.2 [13] hgu133a.db_2.6.3 Hmisc_3.8-3 hopach_2.13.1 [16] IRanges_1.11.31 lattice_0.19-33 limma_3.9.21 [19] Matrix_1.0-0 mgcv_1.7-8 nlme_3.1-102 [22] oligoClasses_1.15.56 RBGL_1.29.0 RCurl_1.6-10.1 [25] SNPchip_1.17.0 splines_2.14.0 survival_2.36-10 [28] tools_2.14.0 XML_3.4-2.2 xtable_1.6-0 [31] zlibbioc_0.1.8 es.c is not correct under Windows. It also gives a different result when i rerun the same function. This is the C code: #include "getEs.h" #include <R.h> #include <Rinternals.h> double absolute(double x) { if (x<0) return (-x); else return (x); } void cumsum(double *x, int len) { int i; for (i = 1; i < len; ++i) { *(x + i) = *(x + i) + *(x + i -1); } } double getNr(double *fchr, int *sign, int signLen) { int i; double nr; nr = 0.0; for (i = 0; i < signLen; ++i) { nr = absolute(fchr[sign[i] -1]) + nr; } return nr; } void getPhit(double *fchr, int *sign, int signLen, double nr, double *phit) { int i; for (i = 0; i < signLen; ++i) { *(phit + sign[i]-1) = absolute(*(fchr + sign[i]-1)) / nr; } } void getPmiss(int *sign, int fchrLen, int signLen, double *pmiss) { int i; double tmp = 1.0 / (fchrLen-signLen); for (i = 0; i < fchrLen; ++i) { *(pmiss + i) = tmp; } for (i = 0; i < signLen; ++i) { *(pmiss + sign[i]-1) = 0; } } SEXP getEs(SEXP fchr, SEXP sign) { int i, nfchr, nsign; double *rfchr = REAL(fchr), *res; int *rsign = INTEGER(sign); PROTECT(fchr = coerceVector(fchr, REALSXP)); PROTECT(sign = coerceVector(sign, REALSXP)); nfchr = length(fchr); nsign = length(sign); SEXP es; PROTECT(es = allocVector(REALSXP, nfchr)); res = REAL(es); double nr = getNr(rfchr, rsign, nsign); SEXP phit; PROTECT(phit = allocVector(REALSXP, nfchr)); double *rphit; rphit = REAL(phit); for(i = 0; i < nfchr; i++) rphit[i] = 0.0; getPhit(rfchr, rsign, nsign, nr, rphit); cumsum(rphit, nfchr); SEXP pmiss; PROTECT(pmiss = allocVector(REALSXP, nfchr)); double *rpmiss; rpmiss = REAL(pmiss); for(i = 0; i < nfchr; i++) rpmiss[i] = 0.0; getPmiss(rsign, nfchr, nsign, rpmiss); cumsum(rpmiss, nfchr); for (i = 0; i < nfchr; ++i) { res[i] = rphit[i] - rpmiss[i]; } UNPROTECT(5); return es; } Could you please help me to find out what I am doing wrong? Many thanks in advance, -- Evarist Planet Research officer, Bioinformatics and Biostatistics unit IRB Barcelona Tel (+34) 93 402 0553 Fax (+34) 93 402 0257 evarist.planet@irbbarcelona.org http://www.irbbarcelona.org/bioinformatics [[alternative HTML version deleted]]
On Mon, Oct 24, 2011 at 6:04 AM, Evarist Planet <evarist.planet at irbbarcelona.org> wrote:> Dear mailing list, > > I have a C function that gives me a wrong result when I run it under Windows > 7.The fact that you extract pointers to the contents of fchr and sign before coercing them to REALSXP is deeply suspicious to me, though I admit I don't see how it causes the problem, since you don't use the new versions, you don't overwrite them, and the old versions should still be protected as arguments to the function.> SEXP getEs(SEXP fchr, SEXP sign) > { > ?int i, nfchr, nsign; > ?double *rfchr = REAL(fchr), *res; > ?int *rsign = INTEGER(sign); > > ?PROTECT(fchr = coerceVector(fchr, REALSXP)); > ?PROTECT(sign = coerceVector(sign, REALSXP)); > > ?nfchr = length(fchr); > ?nsign = length(sign); > > ?SEXP es; > ?PROTECT(es = allocVector(REALSXP, nfchr)); > ?res = REAL(es); > > ?double nr = getNr(rfchr, rsign, nsign); > > ?SEXP phit; > ?PROTECT(phit = allocVector(REALSXP, nfchr)); > ?double *rphit; > ?rphit = REAL(phit); > ?for(i = 0; i < nfchr; i++) rphit[i] = 0.0; > ?getPhit(rfchr, rsign, nsign, nr, rphit); > ?cumsum(rphit, nfchr); > > ?SEXP pmiss; > ?PROTECT(pmiss = allocVector(REALSXP, nfchr)); > ?double *rpmiss; > ?rpmiss = REAL(pmiss); > ?for(i = 0; i < nfchr; i++) rpmiss[i] = 0.0; > ?getPmiss(rsign, nfchr, nsign, rpmiss); > ?cumsum(rpmiss, nfchr); > > ?for (i = 0; i < nfchr; ++i) { > ? res[i] = rphit[i] - rpmiss[i]; > ? } > > ?UNPROTECT(5); > ?return es; > } > > Could you please help me to find out what I am doing wrong? > > Many thanks in advance, > > -- > Evarist Planet > Research officer, Bioinformatics and Biostatistics unit > IRB Barcelona > Tel (+34) 93 402 0553 > Fax (+34) 93 402 0257 > > evarist.planet at irbbarcelona.org > http://www.irbbarcelona.org/bioinformatics > > ? ? ? ?[[alternative HTML version deleted]] > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel >-- Thomas Lumley Professor of Biostatistics University of Auckland
On 10/24/2011 06:04 AM, Evarist Planet wrote:> Dear mailing list, > > I have a C function that gives me a wrong result when I run it under WindowsHi Evarist -- It seems like this can be written reasonably efficiently in R? getEs <-function(fchr, sign) { nfchr <- length(fchr) nsign <- length(sign) nr <- sum(abs(fchr[sign])) phit <- numeric(nfchr) phit[sign] <- abs(fchr[sign]) / nr phit <- cumsum(phit) pmiss <- numeric(nfchr) pmiss[-sign] <- 1 / (nfchr - nsign) pmiss <- cumsum(pmiss) phit - pmiss } es.c <- .Call('getEs',score,s,PACKAGE='phenoTest') all.equal(es.c, getEs(score, s)) (for your C code, it would help to have a simple reproducible example that doesn't rely on phenoTest). Martin> 7. > > This is the code under Linux (RHEL5): >> library(phenoTest) >> data(epheno) >> sign<- sample(featureNames(epheno))[1:20] >> score<- getFc(epheno)[,1] >> head(score) > 1007_s_at 1053_at 117_at 121_at 1255_g_at 1294_at > -1.183019 1.113544 1.186186 -1.034779 -1.044456 -1.023471 >> s<- which(names(score) %in% sign) >> es.c<- .Call('getEs',score,s,PACKAGE='phenoTest') >> head(es.c) > [1] -0.001020408 -0.002040816 -0.003061224 -0.004081633 -0.005102041 > [6] -0.006122449 >> es.c<- .Call('getEs',score,s,PACKAGE='phenoTest') >> head(es.c) > [1] -0.001020408 -0.002040816 -0.003061224 -0.004081633 -0.005102041 > [6] -0.006122449 >> sessionInfo() > R version 2.13.0 (2011-04-13) > Platform: x86_64-redhat-linux-gnu (64-bit) > > locale: > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > [5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8 > [7] LC_PAPER=en_US.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] phenoTest_1.1.1 RSQLite_0.8-4 DBI_0.2-5 > [4] Heatplus_1.22.0 annotate_1.30.1 AnnotationDbi_1.14.1 > [7] Biobase_2.12.2 > > loaded via a namespace (and not attached): > [1] affyio_1.20.0 biomaRt_2.8.1 Biostrings_2.21.6 > [4] Category_2.18.0 cluster_1.13.3 gdata_2.7.1 > [7] genefilter_1.34.0 gplots_2.10.1 graph_1.30.0 > [10] grid_2.13.0 GSEABase_1.14.0 gtools_2.6.2 > [13] hgu133a.db_2.5.0 Hmisc_3.8-3 hopach_2.12.0 > [16] IRanges_1.11.11 lattice_0.19-23 limma_3.8.3 > [19] Matrix_0.9996875-3 mgcv_1.7-5 nlme_3.1-100 > [22] oligoClasses_1.14.0 RBGL_1.22.0 RCurl_1.6-10 > [25] SNPchip_1.16.0 splines_2.13.0 survival_2.36-5 > [28] tools_2.13.0 XML_3.4-0 xtable_1.5-6 > > As you see es.c is correct. I checked it doing the same computation with R. > It also runs without problems under Mac. I run valgrind on the same piece of > code and got no errors. > > This is the same piece of code under Windows 7: >> library(phenoTest) >> data(epheno) >> sign<- sample(featureNames(epheno))[1:20] >> score<- getFc(epheno)[,1] >> head(score) > 1007_s_at 1053_at 117_at 121_at 1255_g_at 1294_at > -1.183019 1.113544 1.186186 -1.034779 -1.044456 -1.023471 >> s<- which(names(score) %in% sign) >> es.c<- .Call('getEs',score,s,PACKAGE='phenoTest') >> head(es.c) > [1] 1.447208e+215 1.447208e+215 1.447208e+215 1.447208e+215 1.447208e+215 > 1.447208e+215 >> es.c<- .Call('getEs',score,s,PACKAGE='phenoTest') >> head(es.c) > [1] 3.176615e+170 3.176615e+170 3.176615e+170 3.176615e+170 3.176615e+170 > 3.176615e+170 >> sessionInfo() > > R version 2.14.0 alpha (2011-10-13 r57240) > Platform: i386-pc-mingw32/i386 (32-bit) > > locale: > [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252 > [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C > [5] LC_TIME=Spanish_Spain.1252 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] phenoTest_1.1.1 RSQLite_0.10.0 DBI_0.2-5 > [4] Heatplus_1.99.0 annotate_1.31.1 AnnotationDbi_1.15.36 > [7] Biobase_2.13.10 > > loaded via a namespace (and not attached): > [1] affyio_1.21.2 biomaRt_2.9.3 Biostrings_2.21.11 > [4] Category_2.19.1 cluster_1.14.0 gdata_2.8.2 > [7] genefilter_1.35.0 gplots_2.10.1 graph_1.31.2 > [10] grid_2.14.0 GSEABase_1.15.0 gtools_2.6.2 > [13] hgu133a.db_2.6.3 Hmisc_3.8-3 hopach_2.13.1 > [16] IRanges_1.11.31 lattice_0.19-33 limma_3.9.21 > [19] Matrix_1.0-0 mgcv_1.7-8 nlme_3.1-102 > [22] oligoClasses_1.15.56 RBGL_1.29.0 RCurl_1.6-10.1 > [25] SNPchip_1.17.0 splines_2.14.0 survival_2.36-10 > [28] tools_2.14.0 XML_3.4-2.2 xtable_1.6-0 > [31] zlibbioc_0.1.8 > > es.c is not correct under Windows. It also gives a different result when i > rerun the same function. > > This is the C code: > #include "getEs.h" > #include<R.h> > #include<Rinternals.h> > > double absolute(double x) > { > if (x<0) > return (-x); > else > return (x); > } > > void cumsum(double *x, int len) > { > int i; > for (i = 1; i< len; ++i) { > *(x + i) = *(x + i) + *(x + i -1); > } > } > > double getNr(double *fchr, int *sign, int signLen) > { > int i; > double nr; > nr = 0.0; > for (i = 0; i< signLen; ++i) { > nr = absolute(fchr[sign[i] -1]) + nr; > } > return nr; > } > > void getPhit(double *fchr, int *sign, int signLen, double nr, double *phit) > { > int i; > for (i = 0; i< signLen; ++i) { > *(phit + sign[i]-1) = absolute(*(fchr + sign[i]-1)) / nr; > } > } > > void getPmiss(int *sign, int fchrLen, int signLen, double *pmiss) > { > int i; > double tmp = 1.0 / (fchrLen-signLen); > for (i = 0; i< fchrLen; ++i) { > *(pmiss + i) = tmp; > } > for (i = 0; i< signLen; ++i) { > *(pmiss + sign[i]-1) = 0; > } > } > > SEXP getEs(SEXP fchr, SEXP sign) > { > int i, nfchr, nsign; > double *rfchr = REAL(fchr), *res; > int *rsign = INTEGER(sign); > > PROTECT(fchr = coerceVector(fchr, REALSXP)); > PROTECT(sign = coerceVector(sign, REALSXP)); > > nfchr = length(fchr); > nsign = length(sign); > > SEXP es; > PROTECT(es = allocVector(REALSXP, nfchr)); > res = REAL(es); > > double nr = getNr(rfchr, rsign, nsign); > > SEXP phit; > PROTECT(phit = allocVector(REALSXP, nfchr)); > double *rphit; > rphit = REAL(phit); > for(i = 0; i< nfchr; i++) rphit[i] = 0.0; > getPhit(rfchr, rsign, nsign, nr, rphit); > cumsum(rphit, nfchr); > > SEXP pmiss; > PROTECT(pmiss = allocVector(REALSXP, nfchr)); > double *rpmiss; > rpmiss = REAL(pmiss); > for(i = 0; i< nfchr; i++) rpmiss[i] = 0.0; > getPmiss(rsign, nfchr, nsign, rpmiss); > cumsum(rpmiss, nfchr); > > for (i = 0; i< nfchr; ++i) { > res[i] = rphit[i] - rpmiss[i]; > } > > UNPROTECT(5); > return es; > } > > Could you please help me to find out what I am doing wrong? > > Many thanks in advance, > > -- > Evarist Planet > Research officer, Bioinformatics and Biostatistics unit > IRB Barcelona > Tel (+34) 93 402 0553 > Fax (+34) 93 402 0257 > > evarist.planet at irbbarcelona.org > http://www.irbbarcelona.org/bioinformatics > > [[alternative HTML version deleted]] > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel-- Computational Biology Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: M1-B861 Telephone: 206 667-2793