Hi R-users,
### reproducible example:
library(gmp)
library(matlab)
library(Matrix)
library(foreach)
library(MASS)
library(mvtnorm)
#####################################
## A function to create a field experiment
## Returns a data frame, called matdf
######################################
setup<-function(b,g,rb,cb,r,c)
{
# where
# b = number of blocks
# g = number of treatments per block
# rb = number of rows per block
# cb = number of columns per block
# r = total number of rows for the layout
# c = total number of columns for the layout
### Check points
stopifnot(is.numeric(b),is.whole(b),is.numeric(g),g>1)
## Compatibility checks
genot<<-seq(1,g,1)
stopifnot(rb*cb==length(genot),r/rb * c/cb == b)
## Generate the design
genotypes<-times(b) %do% sample(genot,g)
block<-rep(1:b,each=length(genot))
genotypes<-factor(genotypes)
block<-factor(block)
### generate the base design
k<-c/cb # number of blocks on the x-axis
y<-rep(rep(1:r,each=cb),k) # X-coordinate
#w<-rb
l<-cb
p<-r/rb
m<-l+1
d<-l*b/p
x<-c(rep(1:l,r),rep(m:d,r)) # Y-coordinate
## compact
matdf<-data.frame(x,y,block,genotypes)
matdf[order(matdf$x),]
}
#########################################################################
## a function to calculate trace from the original data frame
## Returns trace of the variance-covariance matrix, named C22
######################################################################
mainF<-function(matdf,h2=h2,rhox=rhox,rhoy=rhoy,ped="F",s2e=1,c=c,r=r)
{
N<-nrow(matdf)
## Identity matrices
X<<-model.matrix(~matdf$block)
s2g<-varG(s2e,h2)
## calculate G and Z
IG<<-(1/s2g)*eye(length(genot))
z<-model.matrix(~matdf$genotypes-1) # changes everytime there is a
swap
## calculate R and IR
# rhox=rhoy=0.3;s2e=1
# # calculate R and IR
sigx <- diag(c)
sigx <- rhox^ abs(row(sigx) - col(sigx))
sigy <- diag(r)
sigy <- rhoy ^ abs(row(sigy) - col(sigy))
R<- s2e * kronecker(sigx, sigy) # takes 0.01 second
################
# find inverse of R by choleski decomposition
IR<<-chol2inv(chol(R)) # takes about 20 seconds
####
#### brute force matrix multiplication
C11<-t(X)%*%IR%*%X
C11inv<-chol2inv(chol(C11))
k1<<-IR%*%X # 0.2 seconds
k2<-C11inv%*%t(X) # 0 seconds
k3<-k2%*%IR # 0.2 seconds
K<<-k1%*%k3 # 0.16 seconds
### Variance covariance matrices
temp<-t(z)%*%IR%*%z+IG - t(z)%*%K%*%z
C22<-chol2inv(chol(temp))
##########################
## Optimality Criteria
#########################
traceI=sum(diag(C22)) # A-Optimality
traceI
}
## ################################################
### My function to randomly swap two elements in the same block of a
dataframe
### returns a dataframe, called newmatdf
####################################################
swapsimple<-function(matdf)
{
## now, new design after swapping is
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)
names(newmatdf)[names(newmatdf)=="genotypes"] <- "old_G"
names(newmatdf)[names(newmatdf)=="NewG"] <- "genotypes"
#newmatdf <- remove.vars(newmatdf, "old_G")
newmatdf$old_G <- newmatdf$old_G <- NULL
newmatdf[order(newmatdf$x),]
}
#############################################################
### My function to re-calculate trace after swaping the pairs of
elements
################################################
swapmainF<-function(newmatdf)
{
Z<-model.matrix(~newmatdf$genotypes-1)
### Variance covariance matrices
temp<-t(Z)%*%IR%*%Z+IG - t(Z)%*%K%*%Z
C22<-chol2inv(chol(temp))
##########################
## Optimality Criteria
#########################
traceI=sum(diag(C22)) # A-Optimality
traceI
}
#######################################
#I need help in the function below
## I need to avoid the for loops and if loops here :
########################################################
########################################################
#Take an original matrix, called matrix0, calculate its trace. Generate
#a new matrix, called newmatdf after swapping two elements of the old
#one and calculate the trace. If the trace of the newmatrix is smaller
#than that of the previous matrix, store both the current trace
together with
#the older trace and their iteration values. If the newer matrix has a
#trace larger than the previous trace, drop this trace and drop this
#matrix too (but count its iteration).
#Re-swap the old matrix that you stored previously and recalculate the
#trace. Repeat the process many times, say 10,000. The final results
should be a list
#with the original initial matrix and its trace, the final best
#matrix that had the smallest trace after the 10000 simulations and a
#dataframe showing the values of the accepted traces that
#were smaller than the previous and their respective iterations.
####################################################################
funF<- function(newmatdf,n,traceI)
{
matrix0<-newmatdf
trace<-traceI
res <- list(mat = NULL, Design_best = newmatdf, Original_design =
matrix0)
res$mat <- rbind(res$mat, c(value = trace, iterations = 0))
Des<-list()
for(i in seq_len(n)){
ifelse(i==1,
newmatdf<-swapsimple(matrix0),newmatdf<-swapsimple(newmatdf))
Des[[i]]<-newmatdf
if(swapmainF(newmatdf) < trace){
newmatdf<-Des[[i]]
Des[[i]]<-newmatdf
trace<- swapmainF(newmatdf)
res$mat <- rbind(res$mat, c(trace = trace, iterations = i))
res$Design_best <- newmatdf
}
if(swapmainF(newmatdf) > trace & nrow(res$mat)<=1){
newmatdf<-matrix0
Des[[i]]<-matrix0
res$Design_best<-matrix0
}
if(swapmainF(newmatdf)> trace & nrow(res$mat)>1){
newmatdf<-Des[[length(Des)-1]]
Des[[i]]<-newmatdf
res$Design_best<-newmatdf
}
}
res
}
######################################
## Call the function setup, generate 100 designs,
## calculate their traces using the function mainF,
## choose the dataframe with the smallest trace
## store only that dataframe and its trace for later use
######################################
M2F<-function(D,ped="F",k=1,b,g,rb,cb,r,c,
h2,rhox,rhoy,s2e=1)
{
matrix0<-list()
start0 <- c()
value0 <- c()
for (i in 1:D)
{
print(sprintf("generating initial design: %d", i, "complete\n",
sep=""))
flush.console()
matrix0[[i]]<-setup(b=b,g=g,rb=rb,cb=cb,r=r,c=c)
start0[i]<-mainF(matdf=matrix0[[i]],h2=h2,rhox=rhox,rhoy=rhoy,s2e=s2e,c=c,r=r)
s<-which.min(start0)
newmatdf<-matrix0[s][[1]]
trace0<-start0[s][[1]]
}
list(newmatdf=newmatdf,start0=start0,trace0=trace0,index=s)
}
################################################
#### Test my functions ### works perfectly but takes too long
###########################################
b=16;g=196;rb=14;cb=14;r=56;c=56;h2=0.1;rhox=0.3;rhoy=0.3
h2=0.1;rhox=0.3;rhoy=0.3;s2e=1
#
tic() # takes 42.020000 seconds for D==2. but for D==100 , takes
about 30 minutes !!!
res1<-M2F(D=2,ped="F",k=1,b=b,g=g,rb=rb,cb=cb,r=r,c=c,
h2=0.1,rhox=0.3,rhoy=0.3,s2e=1)
toc()
tic() # takes 37.720000 seconds for n==5 but I need for n==4000 or more
takes >7hours
ans1<-funF(res1$newmatdf,traceI=res1$trace0,n=5)
toc()
ans1$mat
regards,
Laz
______________________________________________
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.
