[R] .Call()

[Martin Morgan]

HI Wudd Wudd -- not really answering your question, but it might pay
to refactor some of your R code first. For instance

  for(j in 2:nrow(ranklist)){
      phit <- sum(rep(1/Ns, sum(ranklist[1:j,2]==1)))
      pmiss <- sum(rep(1/(N-Ns), sum(ranklist[1:j,2]==0)))
      if((phit-pmiss)>score[i]) score[i] <- phit - pmiss
  }

might be replaced (you'd better check, I'm not sure I have your logic
correct) by

  phit <- (cumsum(ranklist[,2]==1) / Ns)[-1]
  pmiss <- (cumsum(ranklist[,2]==0) / (N - Ns))[-1]
  pmax <- max(phit - pmiss)
  if (pmax > score[i]) score[i] <- pmax

These should be much faster than your loop. Your earlier 'apply' might
also be replaced by rowMeans, etc. I think, and once simplified to
that level further higher-level reorganizations might be possible.

Hope that helps

Martin

wudd wudd <wuddster at gmail.com> writes:

> Hi guys,
> I want to transfer the following code from R into .Call compatible form. How
> can i do that?
>
> Thanks!!!
>
> INT sim;
>      for(i in 1:sim){
>      if(i>2) genemat <- genemat[,sample(1:ncol(genemat))]
>              ranklist[,1] <- apply(genemat, 1, function(x){
>                                      (mean(x[cols]) -
> mean(x[-cols]))/sd(x)})
>      ranklist <- ranklist[order(ranklist[,1]),]
>      if(ranklist[1,2]==1) score[i] <- 1/Ns
>      if(ranklist[1,2]==0) score[i] <- -(1/(N-Ns))
>   for(j in 2:nrow(ranklist)){
>                          phit <- sum(rep(1/Ns, sum(ranklist[1:j,2]==1)))
>  pmiss <- sum(rep(1/(N-Ns), sum(ranklist[1:j,2]==0)))
>  if((phit-pmiss)>score[i]) score[i] <- phit - pmiss
>          }
>      }
>
> I tried a little bit, but not enough knowledge in C.
> #include <stdio.h>
> #include <R.h>
> #include <Rdefines.h>
> #include <math.h>
> SEXP ESscore(SEXP Rgeneset, SEXP Rgenemat, SEXP Rranklist, SEXP sim)
> {
>      int nc = ncols(Rgenemat);
>  double *geneset = NUMERIC_DATA(Rgeneset);
>  double *genemat = NUMERIC_DATA(Rgenemat);
>
>  SEXP Rscore;
>  PROTECT(Rscore=NEW_NUMERIC(sim));
>  double *score = NUMERIC_DATA(Rscore);
>
>  for(i=1; i<=sim; i++){
>  if(i>2) {genemat <- genemat[,sample(1:nc)];}
>      for(k=1; k<nc )
>
>  }
>
> }
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

-- 
Martin Morgan
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793