[Bioc-devel] C function is wrong under Windows 7.

Steve Lianoglou mailinglist.honeypot at gmail.com
Mon Oct 24 16:23:46 CEST 2011


I think I'd pay more attention to the architecture (32 vs. 64 bit).

These days I find myself doing all of my C(++) <--> R bridge-work via
Rcpp, so I'm a bit rusty w/ the "normal" conventions here, but out of
curiosity: for the vars that essentially never "work with R", like
pmiss for example, why do you go through the whole PROTECT/UNPROTECT
allocVector rigmarole with it instead of just creating a double array
of the appropriate length right off of the bat (then freeing it,
obviously)?

-steve

On Mon, Oct 24, 2011 at 9:40 AM, Evarist Planet
<evarist.planet at irbbarcelona.org> wrote:
> I run the same code on Windows 32-bit R version 2.13.x and got a similar
> result. I do not remember which exact version under 2.13 I used. I'll check
> it and try to run the code on the same version of R.
> Thanks for the quick reply,
> Evarist
>
> On Mon, Oct 24, 2011 at 3:20 PM, Steve Lianoglou
> <mailinglist.honeypot at gmail.com> wrote:
>>
>> Hi,
>>
>> I didn't look at your code, but:
>>
>> Session info for you linux machine shows that it is 64-bit and R version
>> 2.13.0
>> You window box is running in 32-bit and is R version 2.14.0-alpha
>>
>> Maybe you can can try to control one or both of these factors and test
>> again?
>>
>> -steve
>>
>> On Mon, Oct 24, 2011 at 8:33 AM, Evarist Planet
>> <evarist.planet at irbbarcelona.org> wrote:
>> > Ok.
>> >
>> > This is the result of sessionInfo() under linux:
>> >> 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
>> >
>> > And under Windows 7:
>> >> 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
>> >
>> > 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;
>> > }
>> >
>> > Many thanks,
>> >
>> > Evarist
>> >
>> > On Mon, Oct 24, 2011 at 2:22 PM, Vincent Carey
>> > <stvjc at channing.harvard.edu>wrote:
>> >
>> >> your attachment is automatically scrubbed so we cannot see it.  this
>> >> is not really a bioconductor question; you might
>> >> consider r-devel.  you would also do well to include the result of
>> >> sessionInfo().
>> >>
>> >> On Mon, Oct 24, 2011 at 8: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.
>> >> >
>> >> > 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
>> >> >
>> >> > As you see es.c is correct. I checked it doing the same computation
>> >> > with
>> >> > R. It 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
>> >> >
>> >> > es.c is not correct under Windows. It also gives a different result
>> >> > when
>> >> i
>> >> > rerun the same function.
>> >> >
>> >> > I attached the C code.
>> >> >
>> >> > 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
>> >> >
>> >> > _______________________________________________
>> >> > Bioc-devel at r-project.org mailing list
>> >> > https://stat.ethz.ch/mailman/listinfo/bioc-devel
>> >> >
>> >> >
>> >>
>> >
>> >
>> >
>> > --
>> > 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]]
>> >
>> > _______________________________________________
>> > Bioc-devel at r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/bioc-devel
>> >
>>
>>
>>
>> --
>> Steve Lianoglou
>> Graduate Student: Computational Systems Biology
>>  | Memorial Sloan-Kettering Cancer Center
>>  | Weill Medical College of Cornell University
>> Contact Info: http://cbio.mskcc.org/~lianos/contact
>
>
>
> --
> 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
>
>



-- 
Steve Lianoglou
Graduate Student: Computational Systems Biology
 | Memorial Sloan-Kettering Cancer Center
 | Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact



More information about the Bioc-devel mailing list