Greetings,

Thanks Jeff. I appreciate your 'to the point' explanation. Will read into
it more.

Best,
Heramb Gadgil


2013/8/19 Jeff Newmiller <jdnewmil@dcn.davis.ca.us>

> 1. Keeping the number of variables down encourages you to structure your
> data, which allows you to re-use code more efficiently. However, I don't
> believe that the number of variables intrinsically slows down your code
> significantly.. e.g. storing a computed value in a local variable is almost
> always better than repeating the calculation.
>
> 2. Apply functions are not primarily intended for performance
> optimization; they are effective for compactly expressing the idea of your
> algorithms. Using vectorization and indexing is much more effective for
> enhancing performance (and in my mind is even better at expressing
> algorithms than using apply functions).
>
> 3. Agree that too many functions can slow things down, but that is
> normally less the fault of the functions than of the way some programmers
> handle data. It is more productive to focus on the positive idea of using
> vectors to compute as much as possible rather than the negative idea of
> eliminating functions.
>
> 4. I don't think I agree that this statement about memory usage generally
> applies to all uses of R. However, it can be crucial to learn to avoid
> mixing input/output with data processing if you are interested in
> performance.
>
> To the OP: I will pass on trying to analyse your large code base... it is
> requested in the Posting Guide that you isolate your issues and ask
> specific questions.
> ---------------------------------------------------------------------------
> Jeff Newmiller                        The     .....       .....  Go Live...
> DCN:<jdnewmil@dcn.davis.ca.us>        Basics: ##.#.       ##.#.  Live
> Go...
>                                       Live:   OO#.. Dead: OO#..  Playing
> Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
> /Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k
> ---------------------------------------------------------------------------
> Sent from my phone. Please excuse my brevity.
>
> Heramb Gadgil <heramb.gadgil@gmail.com> wrote:
> >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@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@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@aghmed.fsnet.co.uk
> >>>>
> >http://www.aghmed.fsnet.co.uk/**home.html<
> http://www.aghmed.fsnet.co.uk/home.html>
> >>>>
> >>>> ______________________________**________________
> >>>> R-help@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@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@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.
>
>

	[[alternative HTML version deleted]]

