improving/speeding up a very large, slow simulation

David.


> 
> The below was for five iterations of the simulation, and I imagine was the
> reason it was hanging up and stalling out for large simulations:
> 
> 
> Original:
> 
> $by.total
>                         total.time total.pct self.time self.pct
> simulation.unpaired.c       2857.54     99.97    111.90     3.91
> [<-                         2555.74     89.41      2.78     0.10
> [<-.data.frame              2552.96     89.32   2543.84    89.00
> sequential.unpaired.c        149.48      5.23      2.66     0.09
> unpaired.test.c               92.30      3.23     18.80     0.66
> data.frame                    70.72      2.47      9.68     0.34
> as.data.frame                 42.36      1.48      3.74     0.132
> 
> 
> only using dataframe as the output:
> 
> $by.total
>                      total.time total.pct self.time self.pct
> simulation.unpaired.c     109.14    100.00      0.90     0.82
> sequential.unpaired.c      92.16     84.44      3.28     3.01
> unpaired.test.c            84.48     77.41     21.08    19.31
> sd                         43.68     40.02      5.42     4.97
> 
> 
> Thanks again folks.
> 
> 
> On Tue, Feb 12, 2013 at 9:45 PM, Joshua Wiley <jwiley.psych at gmail.com>wrote:
> 
>> Hi Ben,
>> 
>> That is correct about vectorizing----to some speed tests to see the
>> effects.  They can be quite dramatic (like I have easily seen 300fold
>> speedups from doing that).  The apply functions are essentially
>> loops---they tend to be about the same speed as loops, though slightly
>> less code.
>> 
>> Compiler is nice---not it will not help already vectorized code
>> (because vectorized code really just means code that is linked to C
>> code that uses a compiled loop, but C is much faster than R.
>> 
>> Sure, check out the bootstrapping (midway down) section on this page I
>> wrote: http://www.ats.ucla.edu/stat/r/dae/melogit.htm
>> 
>> Cheers,
>> 
>> Josh
>> 
>> On Tue, Feb 12, 2013 at 12:57 PM, Benjamin Caldwell
>> <btcaldwell at berkeley.edu> wrote:
>>> Joshua,
>>> 
>>> So, to your first suggestion - try to figure out whether some operation
>>> performed on every element of the vector at once could do the same thing
>> -
>>> the answer is yes! Where can I see examples of that implemented?
>>> 
>>> e.g. would it just be
>>> 
>>> a <- 1:10000
>>> b <- 1:10000
>>> c <- 1:10000
>>> d <- data.frame(b,c)
>>> 
>>> vectorize<-function(a, d) {
>>> 
>>> f <- d[,1]+d[,2]+a
>>> f
>>> }
>>> 
>>> out<-vectorize(a=a,d=d)
>>> out
>>> 
>>> vs
>>> 
>>> loop<- function(a,d){
>>> 
>>> for(i in 1:length(d[,1])){
>>> a[i]<-d[i,1]+d[i,2]+a[i]
>>> }
>>> a
>>> }
>>> 
>>> out.l<-loop(a,d)
>>> 
>>> out.l
>>> 
>>> If so, it's vectorization that's the big advantage, not something
>> mysterious
>>> that's happening under the hood with the apply functions?
>>> 
>>> To your second : Compiler - thanks for the tip, that's great to know!
>>> 
>>> To your last, could you please give an example or point me to one that
>> was
>>> useful for you? I don't understand that well enough to use it.
>>> 
>>> Thanks!
>>> 
>>> Ben Caldwell
>>> 
>>> Graduate Fellow
>>> University of California, Berkeley
>>> 130 Mulford Hall #3114
>>> Berkeley, CA 94720
>>> Office 223 Mulford Hall
>>> (510)859-3358
>>> 
>>> 
>>> On Mon, Feb 11, 2013 at 10:10 PM, Joshua Wiley <jwiley.psych at gmail.com>
>>> wrote:
>>>> 
>>>> Hi Ben,
>>>> 
>>>> I appreciate that the example is reproducible, and I understand why
>>>> you shared the entire project.  At the same time, that is an awful lot
>>>> of code for volunteers to go through and try to optimize.
>>>> 
>>>> I would try to think of the structure of your task, see what the
>>>> individual pieces are---figure out what can be reused and optimized.
>>>> Look at loops and try to think whether some operation performed on
>>>> every element of the vector at once could do the same thing.
>>>> Sometimes the next iteration of a loop depends on the previous, so it
>>>> is difficult to vectorize, but often that is not the case.
>>>> 
>>>> Make use of the compiler package, especially cmpfun().  It can give
>>>> fairly substantial speedups, particularly with for loops in R.
>>>> 
>>>> Finally if you want 1000 reps, do you actually need to repeat all the
>>>> steps 1000 times, or could you just simulate 1000x the data?  If the
>>>> latter, do the steps once but make more data.  If the former, you
>>>> ought to be able to parallelize it fairly easily, see the parallel
>>>> package.
>>>> 
>>>> Good luck,
>>>> 
>>>> Josh
>>>> 
>>>> 
>>>> On Mon, Feb 11, 2013 at 9:22 PM, Benjamin Caldwell
>>>> <btcaldwell at berkeley.edu> wrote:
>>>>> Dear R help;
>>>>> 
>>>>> I'll preface this by saying that the example I've provided below is
>>>>> pretty
>>>>> long, turgid, and otherwise a deep dive into a series of functions I
>>>>> wrote
>>>>> for a simulation study. It is, however, reproducible and
>> self-contained.
>>>>> 
>>>>> I'm trying to do my first simulation study that's quite big, and so
>> I'll
>>>>> say that the output of this simulation as I'd like it to be is 9.375
>>>>> million rows (9375 combinations of the inputs * 1000). So, first
>> thing,
>>>>> maybe I should try to cut down on the number of dimensions and do it a
>>>>> piece at a time. If that's a conclusion people come to here, I'll look
>>>>> into
>>>>> paring down/doing it a chunk at a time.
>>>>> 
>>>>> I'm hoping, though, that this is an opportunity for me to learn some
>>>>> more
>>>>> efficient, elegant ways to a. vectorize and b. use the apply
>> functions,
>>>>> which still seem to elude me when it comes to their use with more
>>>>> complex,
>>>>> custom functions that have multiple inputs.
>>>>> 
>>>>> If having looked the below over you can give me suggestions to help
>> make
>>>>> this and future code I write run more quickly/efficiently, I will be
>>>>> most grateful, as as this is currently written it would take what
>> looks
>>>>> like days to run a thousand simulations of each possible combination
>> of
>>>>> variables of interest.
>>>>> 
>>>>> Best
>>>>> 
>>>>> Ben Caldwell
>>>>> 
>>>>> -----------------------------------------------
>>>>> #unpaired
>>>>> verification.plots<-rnorm(30, 100, 10)
>>>>> # writeClipboard(as.character(verification.plots))
>>>>> 
>>>>> unpaired.test<- function(verification.plots, project.n, project.mean,
>>>>> project.sd, allowed.deviation, project.acres, alpha=.05){
>>>>> 
>>>>> verification.plots <-as.numeric(as.character(verification.plots))
>>>>> a <- qnorm(alpha/2)
>>>>> d <- alpha*project.mean
>>>>> 
>>>>> # verifier plot number
>>>>> 
>>>>> n<-length(verification.plots)
>>>>> verifier.plot.number <- c(1:n)
>>>>> 
>>>>> 
>>>>> #all plots (verifier and project)
>>>>> 
>>>>> all.plots.n <- rep(0,n)
>>>>> for(i in 1:n){
>>>>> all.plots.n[i] <- project.n + verifier.plot.number[i]
>>>>> }
>>>>> 
>>>>> #running mean
>>>>> X_bar <- rep(1,n)
>>>>> 
>>>>> for(i in 1:n){
>>>>> X_bar[i]<- mean(verification.plots[1:i])
>>>>> }
>>>>> #running sd
>>>>> sd2 <- NULL
>>>>> 
>>>>> for(i in 2:n){
>>>>> sd2[i]<-sd(verification.plots[1:i])
>>>>> }
>>>>> #Tn
>>>>> Tn<-NULL
>>>>> 
>>>>> for(i in 2:n){
>>>>> Tn[i] <- project.mean-X_bar[i]
>>>>> }
>>>>> #n_Threshold
>>>>> n_thresh<-NULL
>>>>> 
>>>>> for(i in 2:n){
>>>>> n_thresh[i] <- (((a^2)/(d^2))*((project.sd^2)+(sd2[i]^2)))
>>>>> }
>>>>> #Upper
>>>>> upper <- NULL
>>>>> for (i in 2:n){
>>>>> upper[i] <- Tn[i]-d
>>>>> }
>>>>> #Lower
>>>>> lower<- NULL
>>>>> for(i in 2:n){
>>>>> lower[i] <- Tn[i]+d
>>>>> }
>>>>> 
>>>>> #Success.fail
>>>>> success.fail <- NULL
>>>>> 
>>>>> 
>>>>> 
>>>>> success.fail.num <- rep(0,n)
>>>>> 
>>>>> 
>> if(abs(verification.plots[1]-project.mean)<(allowed.deviation*project.mean)){
>>>>> success.fail[1] <- "Pass"
>>>>> success.fail.num[1] <- 1
>>>>> } else {
>>>>> success.fail[1] <- "Fail"
>>>>> success.fail.num[1] <-0
>>>>> }
>>>>> for(i in 2:n){
>>>>> if(all.plots.n[i]>=n_thresh[i]){
>>>>> if(upper[i] <= 0){
>>>>> if(lower[i] >= 0){
>>>>> success.fail[i]<-"Pass"
>>>>> success.fail.num[i]<-1
>>>>> } else {
>>>>> success.fail[i] <- "Fail"
>>>>> success.fail.num[i] <-0}
>>>>> } else {
>>>>> success.fail[i] <- "Fail"
>>>>> success.fail.num[i] <- 0
>>>>> }
>>>>> 
>>>>> } else {
>>>>> success.fail[i] <- "Inconclusive"
>>>>> success.fail.num[i] <- 0
>>>>> }
>>>>> }
>>>>> 
>>>>> out.unpaired.test <- data.frame(success.fail, success.fail.num)
>>>>> }
>>>>> 
>>>>> 
>>>>> 
>> out.unpaired.test<-unpaired.test(verification.plots=verification.plots,
>>>>> project.n=30, project.mean=100, project.sd=10, allowed.deviation=.05,
>>>>> project.acres=99)
>>>>> out.unpaired.test
>>>>> #
>>>>> sequential.unpaired<- function(number.strata, project.n, project.mean,
>>>>> project.sd, verification.plots,allowed.deviation, project.acres,
>>>>> min.plots="default", alpha=.05){
>>>>> 
>>>>> if(min.plots=="default"){
>>>>> min.plots<-NULL # minimum passing
>>>>> if(number.strata == 1 & project.acres <= 100){
>>>>>    min.plots = 8
>>>>> } else if(number.strata == 1 & project.acres > 100 & project.acres <
>>>>> 500){
>>>>> min.plots = 12
>>>>> } else if(number.strata == 1 & project.acres > 500 & project.acres <
>>>>> 5000){
>>>>>    min.plots = 16
>>>>>    } else if(number.strata == 1 & project.acres > 5000 &
>> project.acres
>>>>> <
>>>>> 10000){
>>>>>    min.plots = 20
>>>>> } else if(number.strata == 1 & project.acres > 10000){
>>>>>    min.plots = 24
>>>>>    } else if(number.strata == 2 & project.acres < 100){
>>>>>    min.plots = 4
>>>>> } else if(number.strata == 2 & project.acres > 100 & project.acres <
>>>>> 500){
>>>>>    min.plots = 6
>>>>> } else if(number.strata == 2 & project.acres > 501 & project.acres <
>>>>> 5000){
>>>>>    min.plots = 8
>>>>> } else if(number.strata == 2 & project.acres > 5001 & project.acres <
>>>>> 10000){
>>>>>    min.plots = 10
>>>>>    } else if(number.strata == 2 & project.acres > 10000){
>>>>>    min.plots = 12
>>>>> } else if(number.strata == 3 & project.acres < 100){
>>>>>    min.plots = 2
>>>>>    } else if(number.strata == 3 & project.acres > 100 &
>> project.acres <
>>>>> 500){
>>>>> min.plots = 3
>>>>> } else if(number.strata == 3 & project.acres > 501 & project.acres <
>>>>> 5000){
>>>>> min.plots = 4
>>>>> } else if(number.strata == 3 & project.acres > 5001 & project.acres <
>>>>> 10000){
>>>>> min.plots = 5
>>>>> } else if(number.strata == 3 & project.acres > 10000){
>>>>> min.plots = 6
>>>>> } else {min.plots=NULL}
>>>>> } else (min.plots = min.plots)
>>>>> 
>>>>> 
>>>>> if(number.strata > 3){print("max three strata")}
>>>>> if(number.strata > 3){break}
>>>>> 
>>>>> 
>>>>> p <- length(verification.plots)
>>>>> mp <- min.plots
>>>>> run <- unpaired.test(verification.plots, project.n, project.mean,
>>>>> project.sd,
>>>>> allowed.deviation, project.acres, alpha=.05)
>>>>> number.passing.plots <- function(verification.plots,
>> success.fail.num){
>>>>> n<-length(verification.plots)
>>>>> number.passing<-rep(0,n)
>>>>> number.passing[1]<-ifelse(sum(success.fail.num[1]> mp),1,0)
>>>>> for(i in 2:n){
>>>>> if(success.fail.num[i] == 1){
>>>>> v <- i-1
>>>>> number.passing[i]<-number.passing[v]+1
>>>>> } else{number.passing[i] <- 0}
>>>>> }
>>>>> number.passing
>>>>> }
>>>>> 
>>>>> 
>> number.passing<-number.passing.plots(verification.plots=verification.plots,
>>>>> success.fail.num=run$success.fail.num)
>>>>> sucessful.unsucessful<-rep("blank",p)
>>>>> one.zero <- rep(0,p)
>>>>> result <- "blank"
>>>>> resultL <- 0
>>>>> n <-length(verification.plots)
>>>>> for(i in 1:n){
>>>>> if(number.passing[i] >= mp) {
>>>>> sucessful.unsucessful[i] <- "successful"
>>>>> one.zero[i] <- 1
>>>>> } else {
>>>>> sucessful.unsucessful[i]<- "unsuccessful"
>>>>> one.zero[i] <- 0
>>>>> }
>>>>> }
>>>>> 
>>>>> if(sum(one.zero) > 0 ) {
>>>>> result <- "successful"
>>>>> resultL <- 1
>>>>> }
>>>>> 
>>>>> plots.success <- function(one.zero){
>>>>> 
>>>>> plots.to.success<-NULL
>>>>> 
>>>>> for(i in 1:n){
>>>>> if(sum(one.zero)==0){plots.to.success= NA
>>>>> } else if(one.zero[i]==1){plots.to.success<-append(plots.to.success,
>>>>> i)}
>>>>> }
>>>>> plots.to.success<-min(plots.to.success)
>>>>> plots.to.success
>>>>> }
>>>>> 
>>>>> plots.to.success <- plots.success(one.zero=one.zero)
>>>>> 
>>>>> complete <-data.frame (min.plots, number.passing,
>> sucessful.unsucessful,
>>>>> one.zero)
>>>>> collapsed <- data.frame(result, resultL, plots.to.success)
>>>>> out <- c(as.list(complete),as.list(collapsed))
>>>>> }
>>>>> 
>>>>> 
>>>>> unpaired.out<-sequential.unpaired(number.strata=2,
>>>>> verification.plots=verification.plots, project.n=15, project.mean=100,
>>>>> project.sd=10, allowed.deviation=.05, project.acres=99,
>>>>> min.plots="default")
>>>>> 
>>>>> str(unpaired.out)
>>>>> #unpaired.out[[7]][1]
>>>>> 
>>>>> ##
>>>>> 
>>>>> simulation.unpaired <- function(reps, project.acreage.range = c(99,
>>>>> 110,510,5100,11000), project.mean, project.n.min, project.n.max,
>>>>> project.sd.min, project.sd.max, verification.mean.min,
>>>>> verification.mean.max, number.verification.plots.min,
>>>>> number.verification.plots.max, allowed.deviation=.1, alpha=.05,
>>>>> length.out
>>>>> = 5, min.plots="default") {
>>>>> 
>>>>> number.strata.range<-c(1:3) # set by CAR
>>>>> 
>>>>> project.n.range <- seq(project.n.min, project.n.max,
>>>>> length.out=length.out)
>>>>> 
>>>>> project.sd.range <- seq(project.sd.min, project.sd.max,
>>>>> length.out=length.out) # assume verification sd is the same as the
>>>>> project
>>>>> sd to simplify - seems a reasonable simpification
>>>>> 
>>>>> number.verification.plots<- seq(number.verification.plots.min,
>>>>> number.verification.plots.max, length.out=length.out)
>>>>> 
>>>>> verification.range <- seq(verification.mean.min,
>> verification.mean.max,
>>>>> length.out=length.out)
>>>>> 
>>>>> permut.grid<-expand.grid(number.strata.range, project.n.range,
>>>>> project.acreage.range, project.mean, project.sd.range,
>>>>> number.verification.plots, verification.range, allowed.deviation) #
>>>>> create
>>>>> a matrix with all combinations of the supplied vectors
>>>>> 
>>>>> #assign names to the colums of the grid of combinations
>>>>> names.permut<-c("number.strata", "project.n.plots", "project.acreage",
>>>>> "project.mean", "project.sd", "number.verification.plots",
>>>>> "verification.mean", "allowed.deviation")
>>>>> 
>>>>> names(permut.grid)<-names.permut # done
>>>>> 
>>>>> combinations<-length(permut.grid[,1])
>>>>> 
>>>>> size <-reps*combinations #need to know the size of the master matrix,
>>>>> which
>>>>> is the the number of replications * each combination of the supplied
>>>>> factors
>>>>> 
>>>>> # we want a df from which to read all the data into the simulation,
>> and
>>>>> record the results
>>>>> permut.col<-ncol(permut.grid)
>>>>> col.base<-ncol(permut.grid)+2
>>>>> base <- matrix(nrow=size, ncol=col.base)
>>>>> base <-data.frame(base)
>>>>> 
>>>>> # supply the names
>>>>> names.base<-c("number.strata", "project.n.plots", "project.acreage",
>>>>> "project.mean", "project.sd", "number.verification.plots",
>>>>> "verification.mean", "allowed.deviation","success.fail",
>>>>> "plots.to.success")
>>>>> 
>>>>> names(base)<-names.base
>>>>> 
>>>>> # need to create index vectors for the base, master df
>>>>> ends <- seq(reps+1, size+1, by=reps)
>>>>> begins <- ends-reps
>>>>> index <- cbind(begins, ends-1)
>>>>> #done
>>>>> 
>>>>> # next, need to assign the first 6 columns and number of rows = to the
>>>>> number of reps in the simulation to be the given row in the
>> permut.grid
>>>>> matrix
>>>>> 
>>>>> pb <- winProgressBar(title="Create base, big loop 1 of 2", label="0%
>>>>> done",
>>>>> min=0, max=100, initial=0)
>>>>> 
>>>>> for (i in 1:combinations) {
>>>>> 
>>>>> base[index[i,1]:index[i,2],1:permut.col] <- permut.grid[i,]
>>>>> #progress bar
>>>>> info <- sprintf("%d%% done", round((i/combinations)*100))
>>>>> setWinProgressBar(pb, (i/combinations)*100, label=info)
>>>>> }
>>>>> 
>>>>> close(pb)
>>>>> 
>>>>> # now, simply feed the values replicated the number of times we want
>> to
>>>>> run
>>>>> the simulation into the sequential.unpaired function, and assign the
>>>>> values
>>>>> to the appropriate columns
>>>>> 
>>>>> out.index1<-ncol(permut.grid)+1
>>>>> out.index2<-ncol(permut.grid)+2
>>>>> 
>>>>> #progress bar
>>>>> pb <- winProgressBar(title="fill values, big loop 2 of 2", label="0%
>>>>> done",
>>>>> min=0, max=100, initial=0)
>>>>> 
>>>>> for (i in 1:size){
>>>>> 
>>>>> scalar.base <- base[i,]
>>>>> verification.plots <- rnorm(scalar.base$number.verification.plots,
>>>>> scalar.base$verification.mean, scalar.base$project.sd)
>>>>> result<- sequential.unpaired(scalar.base$number.strata,
>>>>> scalar.base$project.n.plots, scalar.base$project.mean, scalar.base$
>>>>> project.sd, verification.plots, scalar.base$allowed.deviation,
>>>>> scalar.base$project.acreage, min.plots='default', alpha)
>>>>> 
>>>>> base[i,out.index1] <- result[[6]][1]
>>>>> base[i,out.index2] <- result[[7]][1]
>>>>> info <- sprintf("%d%% done", round((i/size)*100))
>>>>> setWinProgressBar(pb, (i/size)*100, label=info)
>>>>> }
>>>>> 
>>>>> close(pb)
>>>>> #return the results
>>>>> return(base)
>>>>> 
>>>>> }
>>>>> 
>>>>> # I would like reps to = 1000, but that takes a really long time right
>>>>> now
>>>>> test5 <- simulation.unpaired(reps=5, project.acreage.range = c(99,
>>>>> 110,510,5100,11000), project.mean=100, project.n.min=10,
>>>>> project.n.max=100,
>>>>> project.sd.min=.01,  project.sd.max=.2, verification.mean.min=100,
>>>>> verification.mean.max=130, number.verification.plots.min=10,
>>>>> number.verification.plots.max=50, length.out = 5)
>>>>> 
>>>>>        [[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.
>>>> 
>>>> 
>>>> 
>>>> --
>>>> Joshua Wiley
>>>> Ph.D. Student, Health Psychology
>>>> Programmer Analyst II, Statistical Consulting Group
>>>> University of California, Los Angeles
>>>> https://joshuawiley.com/
>>> 
>>> 
>> 
>> 
>> 
>> --
>> Joshua Wiley
>> Ph.D. Student, Health Psychology
>> Programmer Analyst II, Statistical Consulting Group
>> University of California, Los Angeles
>> https://joshuawiley.com/
>> 
> 
> 	[[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.

David Winsemius
Alameda, CA, USA