[R] How can I make my functions run faster
William Dunlap
wdunlap at tibco.com
Mon Aug 19 18:24:31 CEST 2013
> I am a newbie too. I will share what I do normally for speeding up the code.
>
> 1. Restrict defining too many variables (Global/ Local)
> 2. Use apply functions (apply,sapply,lapply,tapply, etc.) whenever feasible
> 3. Having multiple user defined functions doesn't help. Try to compact
> everything in minimum number of functions
> 4. The in-memory of R is just 10% of your total RAM (Correct me if wrong).
> Make sure most of it is used for processing and not storing
(2) The apply functions do not run a whole lot faster than the equivalent loop
and sometimes they run slower. Their real advantage is that they make the
code more understandable - you don't have to wade through the boilerplate
of setting up a vector or matrix to return or worry about whether all the variables
defined in the loop are temporaries or need to be retained for the rest of the
function. You can get big speedups by avoiding loops and apply functions and
trying to vectorize the code instead (although sometimes vectorization leads
to an unacceptable increase in memory usage).
(3) I think that splitting your problem up into several functions is helpful. You
can test each function on its own to check for correctness and speed. You can
can reuse and share your tested functions. There is some overhead involved
in calling a function but you have to make a lot of calls to notice.
Make your functions "functional" - have almost all of the data come in via
the argument list and have all of the output go out via its return value. Without
this functions are hard to understand and to test. The OP's functions had a lot
of <<-'s in them - get rid of them. They also passed key variables like 'N' through
the global environment (or some other environment, depending on how they
were called). Put that in the argument list so you can easily see how function using it
(rspat) behaves as a function of N (it will be quadratic, bad if N will be getting big).
Know the type, shape, and size of you will want to work on. Some algorithms are
fastest on matrices, some on datasets with many more rows then columns and
some on other types or shapes of data. Run your function on a variety of smallish
datasets to see how it behaves as the number of rows, columns, etc. grows. It
is a waste of time to wait days for a result when you could know that its time is growing
as 0.0001 * nrow(X)^3.
Then there are all the standard recommendations on speeding things up - don't
do an expensive computation if you don't use its result, don't recompute things
every time you pass through a loop, etc.
Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf
> Of Heramb Gadgil
> Sent: Monday, August 19, 2013 8:13 AM
> To: Laz
> Cc: r-help at r-project.org
> Subject: Re: [R] How can I make my functions run faster
>
> Greetings,
>
> I am a newbie too. I will share what I do normally for speeding up the code.
>
> 1. Restrict defining too many variables (Global/ Local)
> 2. Use apply functions (apply,sapply,lapply,tapply, etc.) whenever feasible
> 3. Having multiple user defined functions doesn't help. Try to compact
> everything in minimum number of functions
> 4. The in-memory of R is just 10% of your total RAM (Correct me if wrong).
> Make sure most of it is used for processing and not storing
>
> Hope this will help. Kindly suggest if I have misunderstood anything.
>
> Thanks and Regards,
>
> Heramb Gadgil
>
>
> 2013/8/19 Laz <lmramba at ufl.edu>
>
> > Yes Bert, I am a beginner in writing R functions. I just don't know what
> > to avoid or what to use in order to make the R functions faster.
> >
> > When I run the individual functions, they run quite well.
> > However, calling all of them using the final function it becomes too slow.
> >
> > So I don't know how to make it faster.
> > I used system.time()
> >
> > Regards,
> > Laz
> >
> >
> > On 8/19/2013 10:13 AM, Bert Gunter wrote:
> >
> >> ... and read the "R Language Definition" manual. I noticed unnecessary
> >> constructs
> >> (e.g., z <- f(something); return(z)) that suggest you have more basics
> >> to learn to write efficient, well-structured R code.
> >>
> >> -- Bert
> >>
> >> On Mon, Aug 19, 2013 at 3:55 AM, Michael Dewey <info at aghmed.fsnet.co.uk>
> >> wrote:
> >>
> >>> At 10:28 19/08/2013, Laz wrote:
> >>>
> >>>> Dear R users,
> >>>>
> >>>> I have written a couple of R functions, some are through the help of
> >>>> the R
> >>>> group members. However, running them takes days instead of minutes or a
> >>>> few
> >>>> hours. I am wondering whether there is a quick way of doing that.
> >>>>
> >>>
> >>> Your example code is rather long for humans to profile. Have you thought
> >>> of
> >>> getting R to tell where it is spending most time? The R extensions manual
> >>> tells you how to do this.
> >>>
> >>>
> >>> Here are all my R functions. The last one calls almost all of the
> >>>> previous
> >>>> functions. It is the one I am interested in most. It gives me the
> >>>> correct
> >>>> output but it takes several days to run only 1000 or 2000 simulations!
> >>>> e.g. system.time(test1<-finalF(**designs=5,swaps=20));test1
> >>>> will take about 20 minutes to run but
> >>>> system.time(test1<-finalF(**designs=5,swaps=50));test1 takes about 10
> >>>> hours
> >>>> and system.time(test1<-finalF(**designs=25,swaps=2000));test1 takes
> >>>> about 3
> >>>> days to run
> >>>>
> >>>> Here are my functions
> >>>>
> >>>>
> >>>> ##############################**##############################**
> >>>> #########
> >>>>
> >>>> ls() # list all existing objects
> >>>> rm(list = ls()) # remove them all
> >>>> rm(list = ls()[!grepl("global.var.A", ls())])
> >>>> # refresh memory
> >>>> gc()
> >>>> ls()
> >>>>
> >>>> ### Define a function that requires useful input from the user
> >>>> #b=4;g=seq(1,20,1);rb=5;cb=4;**s2e=1; r=10;c=8
> >>>>
> >>>> ##############################**#######
> >>>> ##############################**######
> >>>> # function to calculate heritability
> >>>> herit<-function(varG,varR=1)
> >>>> {
> >>>> h<-4*varG/(varG+varR)
> >>>> return(c(heritability=h))
> >>>> }
> >>>>
> >>>> ##############################**#####
> >>>> # function to calculate random error
> >>>> varR<-function(varG,h2)
> >>>> {
> >>>> varR<- varG*(4-h2)/h2
> >>>> return(c(random_error=varR))
> >>>> }
> >>>>
> >>>> ##############################**############
> >>>> # function to calculate treatment variance
> >>>> varG<-function(varR=1,h2)
> >>>> {
> >>>> varG<-varR*h2/(4-h2)
> >>>> return(c(treatment_variance=**varG))
> >>>> }
> >>>>
> >>>>
> >>>> ##############################**#
> >>>>
> >>>> # calculating R inverse from spatial data
> >>>> rspat<-function(rhox=0.6,rhoy=**0.6)
> >>>> {
> >>>> s2e<-1
> >>>> R<-s2e*eye(N)
> >>>> for(i in 1:N) {
> >>>> for (j in i:N){
> >>>> y1<-y[i]
> >>>> y2<-y[j]
> >>>> x1<-x[i]
> >>>> x2<-x[j]
> >>>> R[i,j]<-s2e*(rhox^abs(x2-x1))***(rhoy^abs(y2-y1)) # Core
> >>>> AR(1)*AR(1)
> >>>> R[j,i]<-R[i,j]
> >>>> }
> >>>> }
> >>>> IR<-solve(R)
> >>>> IR
> >>>> }
> >>>>
> >>>> ped<<-read.table("ped2new.txt"**,header=FALSE)
> >>>> # Now work on the pedigree
> >>>> ## A function to return Zinverse from pedigree
> >>>>
> >>>> ZGped<-function(ped)
> >>>> {
> >>>> ped2<-data.frame(ped)
> >>>> lenp2<-length(unique(ped2$V1))**;lenp2 # how many Genotypes in
> >>>> total in
> >>>> the pedigree =40
> >>>> ln2<-length(g);ln2#ln2=nrow(**matdf)=30
> >>>> # calculate the new Z
> >>>> Zped<-model.matrix(~ matdf$genotypes -1)# has order N*t = 180 by 30
> >>>> dif<-(lenp2-ln2);dif # 40-30=10
> >>>> #print(c(lenp2,ln2,dif))
> >>>> zeromatrix<-zeros(nrow(matdf),**dif);zeromatrix # 180 by 10
> >>>> Z<-cbind(zeromatrix,Zped) # Design Matrix for random effect
> >>>> (Genotypes):
> >>>> 180 by 40
> >>>> # calculate the new G
> >>>> M<-matrix(0,lenp2,lenp2) # 40 by 40
> >>>> for (i in 1:nrow(ped2)) { M[ped2[i, 1], ped2[i, 2]] <- ped2[i, 3] }
> >>>> G<-s2g*M # Genetic Variance covariance matrix for pedigree 2: 40 by
> >>>> 40
> >>>> IG<-solve(G)
> >>>> return(list(IG=IG, Z=Z))
> >>>> }
> >>>>
> >>>> ##########################
> >>>> ## Required packages #
> >>>> ############################
> >>>> library(gmp)
> >>>> library(knitr) # load this packages for publishing results
> >>>> library(matlab)
> >>>> library(Matrix)
> >>>> library(psych)
> >>>> library(foreach)
> >>>> library(epicalc)
> >>>> library(ggplot2)
> >>>> library(xtable)
> >>>> library(gdata)
> >>>> library(gplots)
> >>>>
> >>>> #b=6;g=seq(1,30,1);rb=5;cb=6;**r=15;c=12;h2=0.3;rhox=0.6;**
> >>>> rhoy=0.6;ped=0
> >>>>
> >>>> setup<-function(b,g,rb,cb,r,c,**h2,rhox=0.6,rhoy=0.6,ped="F")
> >>>> {
> >>>> # where
> >>>> # b = number of blocks
> >>>> # t = number of treatments per block
> >>>> # rb = number of rows per block
> >>>> # cb = number of columns per block
> >>>> # s2g = variance within genotypes
> >>>> # h2 = heritability
> >>>> # r = total number of rows for the layout
> >>>> # c = total number of columns for the layout
> >>>>
> >>>> ### Check points
> >>>> if(b==" ")
> >>>> stop(paste(sQuote("block")," cannot be missing"))
> >>>> if(!is.vector(g) | length(g)<3)
> >>>> stop(paste(sQuote("treatments"**)," should be a vector and
> >>>> more than
> >>>> 2"))
> >>>> if(!is.numeric(b))
> >>>> stop(paste(sQuote("block"),"is not of class",
> >>>> sQuote("numeric")))
> >>>> if(length(b)>1)
> >>>> stop(paste(sQuote("block"),"**has to be only 1 numeric
> >>>> value"))
> >>>> if(!is.whole(b))
> >>>> stop(paste(sQuote("block"),"**has to be an",
> >>>> sQuote("integer")))
> >>>>
> >>>> ## Compatibility checks
> >>>> if(rb*cb !=length(g))
> >>>> stop(paste(sQuote("rb x cb")," should be equal to number of
> >>>> treatment", sQuote("g")))
> >>>> if(length(g) != rb*cb)
> >>>> stop(paste(sQuote("the number of treatments"), "is not equal to",
> >>>> sQuote("rb*cb")))
> >>>>
> >>>> ## Generate the design
> >>>> g<<-g
> >>>> genotypes<-times(b) %do% sample(g,length(g))
> >>>> #genotypes<-rep(g,b)
> >>>> block<-rep(1:b,each=length(g))
> >>>> genotypes<-factor(genotypes)
> >>>> block<-factor(block)
> >>>>
> >>>> ### generate the base design
> >>>> k<-c/cb # number of blocks on the x-axis
> >>>> x<<-rep(rep(1:r,each=cb),k) # X-coordinate
> >>>>
> >>>> #w<-rb
> >>>> l<-cb
> >>>> p<-r/rb
> >>>> m<-l+1
> >>>> d<-l*b/p
> >>>> y<<-c(rep(1:l,r),rep(m:d,r)) # Y-coordinate
> >>>>
> >>>> ## compact
> >>>> matdf<<-data.frame(x,y,block,**genotypes)
> >>>> N<<-nrow(matdf)
> >>>> mm<-summ(matdf)
> >>>> ss<-des(matdf)
> >>>>
> >>>> ## Identity matrices
> >>>> X<<-model.matrix(~block-1)
> >>>> h2<<-h2;rhox<<-rhox;rhoy<<-**rhoy
> >>>> s2g<<-varG(varR=1,h2)
> >>>> ## calculate G and Z
> >>>> ifelse(ped == "F",
> >>>> c(IG<<-(1/s2g)*eye(length(g)),**Z<<-model.matrix(~matdf$**
> >>>> genotypes-1)),
> >>>> c(IG<<- ZGped(ped)[[1]],Z<<-ZGped(ped)**[[2]]))
> >>>> ## calculate R and IR
> >>>> s2e<-1
> >>>> ifelse(rhox==0 | rhoy==0, IR<<-(1/s2e)*eye(N),
> >>>> IR<<-rspat(rhox=rhox,rhoy=**rhoy))
> >>>> C11<-t(X)%*%IR%*%X
> >>>> C11inv<-solve(C11)
> >>>> K<<-IR%*%X%*%C11inv%*%t(X)%*%**IR
> >>>> return(list(matdf=matdf,**summary=mm,description=ss))
> >>>>
> >>>> }
> >>>>
> >>>>
> >>>> #setup(b=6,g=seq(1,30,1),rb=5,**cb=6,r=15,c=12,h2=0.3,rhox=0.**
> >>>> 6,rhoy=0.6,ped="F")[1]
> >>>>
> >>>> #system.time(out3<-setup(b=6,**g=seq(1,30,1),rb=5,cb=6,r=15,**
> >>>> c=12,h2=0.3,rhox=0.6,rhoy=0.6,**ped="F"));out3
> >>>>
> >>>> #system.time(out4<-setup(b=16,**g=seq(1,196,1),rb=14,cb=14,r=**
> >>>> 56,c=56,h2=0.3,rhox=0.6,rhoy=**0.6,ped="F"));out4
> >>>>
> >>>>
> >>>> ##############################**######################
> >>>> # The function below uses shortcuts from textbook by Harville 1997
> >>>> # uses inverse of a partitioned matrix technique
> >>>> ##############################**######################
> >>>>
> >>>> mainF<-function(criteria=c("A"**,"D"))
> >>>> {
> >>>> ### Variance covariance matrices
> >>>> temp<-t(Z)%*%IR%*%Z+IG - t(Z)%*%K%*%Z
> >>>> C22<-solve(temp)
> >>>> ##########################
> >>>> ## Optimality Criteria
> >>>> #########################
> >>>> traceI<<-sum(diag(C22)) ## A-Optimality
> >>>> doptimI<<-log(det(C22)) # D-Optimality: minimize the det of the
> >>>> inverse
> >>>> of Inform Matrix
> >>>> #return(c(traceI,doptimI))
> >>>> if(criteria=="A") return(traceI)
> >>>> if(criteria=="D") return(doptimI)
> >>>> else{return(c(traceI,doptimI))**}
> >>>> }
> >>>>
> >>>> # system.time(res1<-mainF(**criteria="A"));res1
> >>>> # system.time(res2<-mainF(**criteria="D"));res2
> >>>> #system.time(res3<-mainF(**criteria="both"));res3
> >>>>
> >>>>
> >>>> ##############################**################
> >>>> ### Swap function that takes matdf and returns
> >>>> ## global values newnatdf and design matrices
> >>>> ### Z and IG
> >>>> ##############################**################
> >>>>
> >>>> swapsimple<-function(matdf,**ped="F")
> >>>> {
> >>>> # dataset D =mat1 generated from the above function
> >>>> ## now, new design after swapping is
> >>>> matdf<-as.data.frame(matdf)
> >>>> attach(matdf,warn.conflict=**FALSE)
> >>>> b1<-sample(matdf$block,1,**replace=TRUE);b1
> >>>> gg1<-matdf$genotypes[block==**b1];gg1
> >>>> g1<-sample(gg1,2);g1
> >>>> samp<-Matrix(c(g1=g1,block=b1)**,nrow=1,ncol=3,
> >>>> dimnames=list(NULL,c("gen1","**gen2","block")));samp
> >>>> newGen<-matdf$genotypes
> >>>> newG<-ifelse(matdf$genotypes==**samp[,1] &
> >>>> block==samp[,3],samp[,2],**matdf$genotypes)
> >>>> NewG<-ifelse(matdf$genotypes==**samp[,2] &
> >>>> block==samp[,3],samp[,1],newG)
> >>>> NewG<-factor(NewG)
> >>>>
> >>>> ## now, new design after swapping is
> >>>> newmatdf<-cbind(matdf,NewG)
> >>>> newmatdf<<-as.data.frame(**newmatdf)
> >>>> mm<-summ(newmatdf)
> >>>> ss<-des(newmatdf)
> >>>>
> >>>> ## Identity matrices
> >>>> ifelse(ped == "F",
> >>>> c(IG<<-(1/s2g)*eye(length(g)),**Z<<-model.matrix(~newmatdf$**NewG-1)),
> >>>> c(IG<<-
> >>>> ZGped(ped)[[1]],Z<<-ZGped(ped)**[[2]]))
> >>>> ## calculate R and IR
> >>>> C11<-t(X)%*%IR%*%X
> >>>> C11inv<-solve(C11)
> >>>> K<<-IR%*%X%*%C11inv%*%t(X)%*%**IR
> >>>> return(list(newmatdf=newmatdf,**summary=mm,description=ss))
> >>>> }
> >>>> #swapsimple(matdf,ped="F")[c(**2,3)]
> >>>> #which(newmatdf$genotypes != newmatdf$NewG)
> >>>> ##############################**#############
> >>>> # for one design, swap pairs of treatments
> >>>> # several times and store the traces
> >>>> # of the successive swaps
> >>>> ##############################**############
> >>>>
> >>>> optmF<-function(iterations=2,**verbose=FALSE)
> >>>> {
> >>>> trace<-c()
> >>>>
> >>>> for (k in 1:iterations){
> >>>>
> >>>> setup(b=6,g=seq(1,30,1),rb=5,**cb=6,r=15,c=12,h2=0.3,rhox=0.**
> >>>> 6,rhoy=0.6,ped="F")
> >>>> swapsimple(matdf,ped="F")
> >>>> trace[k]<-mainF(criteria="A")
> >>>> iterations[k]<-k
> >>>> mat<-cbind(trace, iterations= seq(iterations))
> >>>> }
> >>>>
> >>>> if (verbose){
> >>>> cat("***starting matrix\n")
> >>>> print(mat)
> >>>> }
> >>>> # iterate till done
> >>>> while(nrow(mat) > 1){
> >>>> high <- diff(mat[, 'trace']) > 0
> >>>> if (!any(high)) break # done
> >>>> # find which one to delete
> >>>> delete <- which.max(high) + 1L
> >>>> #mat <- mat[-delete, ]
> >>>> mat <- mat[-delete,, drop=FALSE]
> >>>> }
> >>>> mat
> >>>> }
> >>>>
> >>>> #system.time(test1<-optmF(**iterations=10));test1
> >>>>
> >>>> ##############################**##################
> >>>> ##############################**#################
> >>>>
> >>>> swap<-function(matdf,ped="F",**criteria=c("A","D"))
> >>>> {
> >>>> # dataset D =mat1 generated from the above function
> >>>> ## now, new design after swapping is
> >>>> matdf<-as.data.frame(matdf)
> >>>> attach(matdf,warn.conflict=**FALSE)
> >>>> b1<-sample(matdf$block,1,**replace=TRUE);b1
> >>>> gg1<-matdf$genotypes[block==**b1];gg1
> >>>> g1<-sample(gg1,2);g1
> >>>> samp<-Matrix(c(g1=g1,block=b1)**,nrow=1,ncol=3,
> >>>> dimnames=list(NULL,c("gen1","**gen2","block")));samp
> >>>> newGen<-matdf$genotypes
> >>>> newG<-ifelse(matdf$genotypes==**samp[,1] &
> >>>> block==samp[,3],samp[,2],**matdf$genotypes)
> >>>> NewG<-ifelse(matdf$genotypes==**samp[,2] &
> >>>> block==samp[,3],samp[,1],newG)
> >>>> NewG<-factor(NewG)
> >>>>
> >>>> ## now, new design after swapping is
> >>>> newmatdf<-cbind(matdf,NewG)
> >>>> newmatdf<<-as.data.frame(**newmatdf)
> >>>> mm<-summ(newmatdf)
> >>>> ss<-des(newmatdf)
> >>>>
> >>>> ## Identity matrices
> >>>> #X<<-model.matrix(~block-1)
> >>>> #s2g<<-varG(varR=1,h2)
> >>>> ## calculate G and Z
> >>>> ifelse(ped == "F",
> >>>> c(IG<<-(1/s2g)*eye(length(g)),**Z<<-model.matrix(~newmatdf$**NewG-1)),
> >>>> c(IG<<-
> >>>> ZGped(ped)[[1]],Z<<-ZGped(ped)**[[2]]))
> >>>> ## calculate R and IR
> >>>> C11<-t(X)%*%IR%*%X
> >>>> C11inv<-solve(C11)
> >>>> K<-IR%*%X%*%C11inv%*%t(X)%*%IR
> >>>> temp<-t(Z)%*%IR%*%Z+IG - t(Z)%*%K%*%Z
> >>>> C22<-solve(temp)
> >>>> ##########################
> >>>> ## Optimality Criteria
> >>>> #########################
> >>>> traceI<-sum(diag(C22)) ## A-Optimality
> >>>> doptimI<-log(det(C22)) #
> >>>> #return(c(traceI,doptimI))
> >>>> if(criteria=="A") return(traceI)
> >>>> if(criteria=="D") return(doptimI)
> >>>> else{return(c(traceI,doptimI))**}
> >>>> }
> >>>>
> >>>> #swap(matdf,ped="F",criteria="**both")
> >>>>
> >>>> ##############################**#############
> >>>> ### Generate 25 initial designs
> >>>> ##############################**#############
> >>>> #rspatf<-function(design){
> >>>> # arr = array(1, dim=c(nrow(matdf),ncol(matdf)+**1,design))
> >>>> # l<-list(length=dim(arr)[3])
> >>>> # for (i in 1:dim(arr)[3]){
> >>>> # l[[i]]<-swapsimple(matdf,ped="**F")[[1]][,,i]
> >>>> # }
> >>>> # l
> >>>> #}
> >>>> #matd<-rspatf(design=5)
> >>>> #matd
> >>>>
> >>>> #which(matd[[1]]$genotypes != matd[[1]]$NewG)
> >>>> #which(matd[[2]]$genotypes != matd[[2]]$NewG)
> >>>>
> >>>>
> >>>> ##############################**#################
> >>>> ##############################**#################
> >>>>
> >>>> optm<-function(iterations=2,**verbose=FALSE)
> >>>> {
> >>>> trace<-c()
> >>>>
> >>>> for (k in 1:iterations){
> >>>>
> >>>> setup(b=6,g=seq(1,30,1),rb=5,**cb=6,r=15,c=12,h2=0.3,rhox=0.**
> >>>> 6,rhoy=0.6,ped="F")
> >>>> trace[k]<-swap(matdf,ped="F",**criteria="A")
> >>>> iterations[k]<-k
> >>>> mat<-cbind(trace, iterations= seq(iterations))
> >>>> }
> >>>>
> >>>> if (verbose){
> >>>> cat("***starting matrix\n")
> >>>> print(mat)
> >>>> }
> >>>> # iterate till done
> >>>> while(nrow(mat) > 1){
> >>>> high <- diff(mat[, 'trace']) > 0
> >>>> if (!any(high)) break # done
> >>>> # find which one to delete
> >>>> delete <- which.max(high) + 1L
> >>>> #mat <- mat[-delete, ]
> >>>> mat <- mat[-delete,, drop=FALSE]
> >>>> }
> >>>> mat
> >>>> }
> >>>>
> >>>> #system.time(res<-optm(**iterations=10));res
> >>>> ##############################**###################
> >>>> ##############################**##################
> >>>> finalF<-function(designs,**swaps)
> >>>> {
> >>>> Nmatdf<-list()
> >>>> OP<-list()
> >>>> Miny<-NULL
> >>>> Maxy<-NULL
> >>>> Minx<-NULL
> >>>> Maxx<-NULL
> >>>> for (i in 1:designs)
> >>>> {
> >>>>
> >>>> setup(b=4,g=seq(1,20,1),rb=5,**cb=4,r=10,c=8,h2=0.3,rhox=0.6,**
> >>>> rhoy=0.6,ped="F")[1]
> >>>> mainF(criteria="A")
> >>>> for (j in 1:swaps)
> >>>> {
> >>>> OP[[i]]<- optmF(iterations=swaps)
> >>>> Nmatdf[[i]]<-newmatdf[,5]
> >>>> Miny[i]<-min(OP[[i]][,1])
> >>>> Maxy[i]<-max(OP[[i]][,1])
> >>>> Minx[i]<-min(OP[[i]][,2])
> >>>> Maxx[i]<-max(OP[[i]][,2])
> >>>> }
> >>>> }
> >>>> return(list(OP=OP,Miny=Miny,**Maxy=Maxy,Minx=Minx,Maxx=Maxx,**
> >>>> Nmatdf=Nmatdf))
> >>>> # gives us both the Optimal conditions and designs
> >>>> }
> >>>>
> >>>> ##############################**###################
> >>>> sink(file= paste(format(Sys.time(),
> >>>> "Final_%a_%b_%d_%Y_%H_%M_%S"),**"txt",sep="."),split=TRUE)
> >>>> system.time(test1<-finalF(**designs=25,swaps=2000));test1
> >>>> sink()
> >>>>
> >>>>
> >>>> I expect results like this below
> >>>>
> >>>> sink()
> >>>>> finalF<-function(designs,**swaps)
> >>>>>
> >>>> +{
> >>>> + Nmatdf<-list()
> >>>> + OP<-list()
> >>>> + Miny<-NULL
> >>>> + Maxy<-NULL
> >>>> + Minx<-NULL
> >>>> + Maxx<-NULL
> >>>> + for (i in 1:designs)
> >>>> + {
> >>>> +
> >>>> setup(b=4,g=seq(1,20,1),rb=5,**cb=4,r=10,c=8,h2=0.3,rhox=0.6,**
> >>>> rhoy=0.6,ped="F")[1]
> >>>> + mainF(criteria="A")
> >>>> + for (j in 1:swaps)
> >>>> + {
> >>>> + OP[[i]]<- optmF(iterations=swaps)
> >>>> + Nmatdf[[i]]<-newmatdf[,5]
> >>>> + Miny[i]<-min(OP[[i]][,1])
> >>>> + Maxy[i]<-max(OP[[i]][,1])
> >>>> + Minx[i]<-min(OP[[i]][,2])
> >>>> + Maxx[i]<-max(OP[[i]][,2])
> >>>> + }
> >>>> + }
> >>>> +
> >>>>
> return(list(OP=OP,Miny=Miny,**Maxy=Maxy,Minx=Minx,Maxx=Maxx,**Nmatdf=Nmatdf
> ))
> >>>> #
> >>>> gives us both the Optimal conditions and designs
> >>>> +}
> >>>>
> >>>>> sink(file= paste(format(Sys.time(),
> >>>>> "Final_%a_%b_%d_%Y_%H_%M_%S"),**"txt",sep="."),split=TRUE)
> >>>>> system.time(test1<-finalF(**designs=5,swaps=5));test1
> >>>>>
> >>>> user system elapsed
> >>>> 37.88 0.00 38.04
> >>>> $OP
> >>>> $OP[[1]]
> >>>> trace iterations
> >>>> [1,] 0.8961335 1
> >>>> [2,] 0.8952822 3
> >>>> [3,] 0.8934649 4
> >>>>
> >>>> $OP[[2]]
> >>>> trace iterations
> >>>> [1,] 0.893955 1
> >>>>
> >>>> $OP[[3]]
> >>>> trace iterations
> >>>> [1,] 0.9007225 1
> >>>> [2,] 0.8971837 4
> >>>> [3,] 0.8902474 5
> >>>>
> >>>> $OP[[4]]
> >>>> trace iterations
> >>>> [1,] 0.8964726 1
> >>>> [2,] 0.8951722 4
> >>>>
> >>>> $OP[[5]]
> >>>> trace iterations
> >>>> [1,] 0.8973285 1
> >>>> [2,] 0.8922594 4
> >>>>
> >>>>
> >>>> $Miny
> >>>> [1] 0.8934649 0.8939550 0.8902474 0.8951722 0.8922594
> >>>>
> >>>> $Maxy
> >>>> [1] 0.8961335 0.8939550 0.9007225 0.8964726 0.8973285
> >>>>
> >>>> $Minx
> >>>> [1] 1 1 1 1 1
> >>>>
> >>>> $Maxx
> >>>> [1] 4 1 5 4 4
> >>>>
> >>>> $Nmatdf
> >>>> $Nmatdf[[1]]
> >>>> [1] 30 8 5 28 27 29 1 26 24 22 13 6 17 18 2 19 14 11 3 23 10
> >>>> 15 21
> >>>> 9 25 4 7 20 12 16 14 17 15 5 8 6 19
> >>>> [38] 4 1 10 11 3 24 20 13 2 27 12 16 28 21 23 30 25 29 7 26 18 9
> >>>> 22
> >>>> 24 21 26 2 13 30 5 28 20 11 3 7 18 25
> >>>> [75] 22 16 4 17 19 27 29 10 23 6 12 15 14 1 9 8 12 11 3 8 5
> >>>> 20 23
> >>>> 22 7 15 19 29 24 27 13 2 6 1 21 26 25
> >>>> [112] 10 16 14 18 4 30 17 9 28 29 9 7 27 11 2 30 18 8 14 19 20 15
> >>>> 21
> >>>> 4 3 16 24 13 28 26 10 12 6 5 25 1 17
> >>>> [149] 23 22 21 2 23 16 4 10 9 22 30 24 1 27 3 20 12 5 26 17 28 11
> >>>> 7
> >>>> 14 8 25 19 13 18 29 15 6
> >>>> Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
> >>>> 25
> >>>> 26 27 28 29 30
> >>>>
> >>>> $Nmatdf[[2]]
> >>>> [1] 5 13 30 2 21 23 6 27 16 19 8 26 18 4 20 9 22 28 7 3 15
> >>>> 10 11
> >>>> 17 25 24 29 1 14 12 28 18 23 19 21 16 17
> >>>> [38] 29 13 7 15 27 25 22 10 1 2 5 30 9 20 3 14 24 26 4 6 12
> >>>> 11 8
> >>>> 8 18 25 12 5 23 21 4 9 17 20 1 2 6
> >>>> [75] 22 7 16 26 30 29 3 15 19 14 13 11 24 28 27 10 16 21 26 23 25 4
> >>>> 9
> >>>> 24 15 14 22 1 20 27 2 7 17 18 13 8 12
> >>>> [112] 5 6 19 28 3 10 30 11 29 11 30 14 9 26 5 1 10 29 28 4 18 8
> >>>> 24
> >>>> 20 13 3 23 27 6 15 16 21 2 17 7 25 12
> >>>> [149] 19 22 7 28 8 11 26 24 12 29 9 16 21 27 22 23 18 19 13 6 15 3
> >>>> 1
> >>>> 30 2 17 14 5 25 20 4 10
> >>>> Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
> >>>> 25
> >>>> 26 27 28 29 30
> >>>>
> >>>> $Nmatdf[[3]]
> >>>> [1] 7 25 4 30 12 11 14 13 26 1 10 21 15 22 29 19 27 16 2 24 28
> >>>> 20 3
> >>>> 5 23 8 18 6 17 9 6 21 9 15 11 17 13
> >>>> [38] 29 24 4 20 7 23 14 2 16 18 26 19 25 8 1 12 10 28 27 22 30 5
> >>>> 3
> >>>> 20 12 8 2 11 18 24 19 9 22 15 7 30 27
> >>>> [75] 17 29 6 3 5 1 21 25 28 14 23 4 16 26 13 10 20 29 26 25 15
> >>>> 22 9
> >>>> 10 28 17 18 21 6 16 7 1 3 24 11 2 4
> >>>> [112] 14 8 5 13 27 23 30 19 12 6 30 1 2 7 28 18 8 20 10 4 25 14
> >>>> 19
> >>>> 27 11 13 29 12 9 3 26 22 21 16 15 17 24
> >>>> [149] 5 23 17 6 25 11 21 29 5 26 13 7 15 2 9 4 18 30 3 8 20 24
> >>>> 27
> >>>> 22 19 16 28 12 1 23 14 10
> >>>> Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
> >>>> 25
> >>>> 26 27 28 29 30
> >>>>
> >>>> $Nmatdf[[4]]
> >>>> [1] 24 8 17 30 10 20 4 28 25 16 14 13 7 12 26 29 21 19 1 22 11 6
> >>>> 23
> >>>> 18 15 5 27 2 3 9 1 24 27 15 26 14 28
> >>>> [38] 20 8 5 4 29 2 25 9 13 6 21 7 22 30 17 3 10 12 19 11 18
> >>>> 16 23
> >>>> 25 18 3 29 1 4 8 6 9 30 2 14 11 16
> >>>> [75] 23 13 10 12 7 19 17 5 21 28 24 20 15 27 26 22 14 5 7 6 17 3
> >>>> 1
> >>>> 29 25 23 19 11 21 18 4 30 20 8 2 12 9
> >>>> [112] 16 10 15 27 26 13 24 28 22 19 7 17 1 12 8 18 16 14 22 3 28 27
> >>>> 25
> >>>> 10 6 4 15 30 9 11 5 20 26 24 29 21 2
> >>>> [149] 23 13 2 16 10 25 18 15 26 22 12 19 30 17 23 8 3 7 20 14 13 28
> >>>> 9
> >>>> 21 11 29 6 5 4 24 27 1
> >>>> Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
> >>>> 25
> >>>> 26 27 28 29 30
> >>>>
> >>>> $Nmatdf[[5]]
> >>>> [1] 12 18 8 22 9 21 2 1 29 13 30 25 17 6 16 5 26 7 3 14 23
> >>>> 15 28
> >>>> 27 10 24 20 11 19 4 20 30 14 27 25 4 6
> >>>> [38] 28 23 8 9 29 26 19 24 7 5 1 11 22 21 2 10 18 12 15 3 17
> >>>> 13 16
> >>>> 16 22 6 9 21 5 14 2 30 10 3 25 27 15
> >>>> [75] 28 7 17 20 11 8 19 29 12 26 24 13 1 4 18 23 4 16 10 25 5
> >>>> 13 18
> >>>> 19 22 7 28 30 23 21 11 2 14 9 20 24 8
> >>>> [112] 17 1 15 29 6 12 27 3 26 14 8 26 6 20 9 15 23 3 22 7 30 25
> >>>> 24
> >>>> 1 10 19 21 4 11 2 18 17 13 28 29 27 16
> >>>> [149] 12 5 19 2 4 5 15 21 17 7 25 8 6 16 20 29 10 18 1 12 26 28
> >>>> 27
> >>>> 11 14 23 22 9 3 13 30 24
> >>>> Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
> >>>> 25
> >>>> 26 27 28 29 30
> >>>>
> >>>>
> >>>> Michael Dewey
> >>> info at aghmed.fsnet.co.uk
> >>>
> http://www.aghmed.fsnet.co.uk/**home.html<http://www.aghmed.fsnet.co.uk/home.
> html>
> >>>
> >>> ______________________________**________________
> >>> R-help at r-project.org mailing list
> >>> https://stat.ethz.ch/mailman/**listinfo/r-
> help<https://stat.ethz.ch/mailman/listinfo/r-help>
> >>> PLEASE do read the posting guide http://www.R-project.org/**
> >>> posting-guide.html <http://www.R-project.org/posting-guide.html>
> >>> and provide commented, minimal, self-contained, reproducible code.
> >>>
> >>
> >>
> >>
> > ______________________________**________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/**listinfo/r-help<https://stat.ethz.ch/mailman/listinfo/r-
> help>
> > PLEASE do read the posting guide http://www.R-project.org/**
> > posting-guide.html <http://www.R-project.org/posting-guide.html>
> > and provide commented, minimal, self-contained, reproducible code.
> >
>
> [[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.
More information about the R-help
mailing list