Running 'mmec_examples.R' [156s]
  Running 'testthat.R' [2s]
Running the tests in 'tests/mmec_examples.R' failed.
Complete output:
  > # plotMonitor <- function(x,...){
  > #   mmin <- min(x)
  > #   mmax <- max(x)
  > #   plot(x[1,], type="l", ylim=c(mmin,mmax))
  > #   if(nrow(x) > 1){
  > #     for(i in 2:nrow(x)){
  > #       par(new=TRUE)
  > #       plot(x[i,], col=i, ylim=c(mmin,mmax), ylab="",xlab="", type="l",...)
  > #     }
  > #   }
  > #   legend("topright",legend = rownames(x), col=1:(nrow(x)), bty="n", lty=1)
  > # }
  > # 
  > # library(sommer)
  > # data("DT_cpdata")
  > # DT <- DT_cpdata
  > # GT <- GT_cpdata
  > # MP <- MP_cpdata
  > # DT$Yield <- imputev(DT$Yield)
  > # #### create the variance-covariance matrix
  > # K <- A.mat(GT) # additive relationship matrix
  > # K <- K + diag(1e-3,nrow(K))
  > # #### look at the data and fit the model
  > # head(DT)
  > # mix1 <- mmer(Yield~1,
  > #              random=~vsr(id,Gu=K)
  > #              + Rowf,
  > #              rcov=~units,
  > #              data=DT)
  > # summary(mix1)$varcomp
  > # 
  > # Ki <- as(solve(K), Class="sparseMatrix")
  > # mix2 <- mmec(fixed=Yield~1,
  > #              random=~vsc(ids(id),Gu=Ki) + Rowf,
  > #              rcov=~units, return.param = F,
  > #              data=DT)
  > # 
  > # mix2$monitor[,15]
  > # plotMonitor(mix2$monitor)
  > # pp <- predict.mmec(object=mix2,classify = c("id","(Intercept)"))#, ignore="Rowf")
  > # head(pp$pvals)
  > # 
  > # Z <- list(model.matrix(~id-1, data=DT),model.matrix(~Rowf-1,data=DT)
  > # )
  > # Zind <- 1:2
  > # 
  > # A <- list(
  > #   K, 
  > #   diag(13)
  > # ) #
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # thetaI <- list(
  > #   matrix(1000,1,1),
  > #   matrix(1000,1,1),
  > #   matrix(1000,1,1)
  > # );thetaI
  > # 
  > # thetaC <- list(
  > #   matrix(1,1,1),
  > #   matrix(1,1,1),
  > #   matrix(1,1,1)
  > # );thetaC
  > # 
  > # X <- model.matrix(~1, data=DT)
  > # 
  > # y <- as.matrix(DT$Yield)
  > # 
  > # S <- list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # thetaF
  > # 
  > # ## apply the function
  > # weightInf=rep(1,40) # weights for the information matrix
  > # weightEmInf=c(seq(.9,.1,-.2),rep(0,36));weightEmInf # weights for the EM information matrix
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=TRUE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # # compare results
  > # res3$monitor
  > # summary(mix1)$varcomp
  > # 
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # # CS + DIAGONAL MODEL
  > # 
  > # library(sommer)
  > # data("DT_example")
  > # DT <- DT_example
  > # K <- A_example
  > # #### look at the data and fit the model
  > # head(DT)
  > # DT$Yield <- scale(DT$Yield)
  > # mix1 <- mmer(Yield~Env,
  > #              random= ~Name + vsr(dsr(Env),Name),
  > #              rcov= ~ vsr(dsr(Env),units),
  > #              data=DT)
  > # summary(mix1)$varcomp
  > # 
  > # mix2 <- mmec(Yield~Env,
  > #              random= ~ Name + vsc(dsr(Env),ids(Name)),
  > #              rcov= ~ vsc(dsr(Env),ids(units)),
  > #              return.param = F,
  > #              data=DT)
  > # plotMonitor(mix2$monitor)
  > # summary(mix1)$varcomp
  > # mix2$sigma
  > # 
  > # zz <- with(DT, vsr(dsr(Env),Name))
  > # 
  > # Z <- c(list(model.matrix(~Name-1, data=DT)),zz$Z)
  > # 
  > # Zind <- c(1,2,2,2)
  > # 
  > # A <- list(diag(41), diag(41))#rep(list(diag(41)),4)
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # thetaI <- list(
  > #   matrix(10,1,1),
  > #   diag(10,3,3),
  > #   diag(10,3,3)
  > # );thetaI
  > # 
  > # thetaC <- list(
  > #   matrix(1,1,1),
  > #   diag(1,3,3),
  > #   diag(1,3,3)
  > # );thetaC
  > # 
  > # X <- model.matrix(~Env, data=DT)
  > # 
  > # y <- as.matrix(DT$Yield)
  > # 
  > # DTx <- DT; DTx$units <- as.factor(1:nrow(DTx))
  > # ss <- with(DTx, vsr(dsr(Env),units) )
  > # 
  > # S <- ss$Z #list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # 
  > # ## apply the function
  > # weightInf=rep(1,40); # weights for the information matrix
  > # weightEmInf=c(seq(.9,.1,-.1),rep(0,36)); # weights for the EM information matrix
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=TRUE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # # compare results
  > # res3$monitor
  > # summary(mix1)$varcomp
  > # 
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # # UNSTRUCTURED MODEL
  > # library(sommer)
  > # data("DT_example")
  > # DT <- DT_example
  > # K <- A_example
  > # #### look at the data and fit the model
  > # head(DT)
  > # DT$Yield <- scale(DT$Yield)
  > # mix1 <- mmer(Yield~Env,
  > #              random= ~ vsr(usr(Env),Name),
  > #              rcov= ~ vsr(dsr(Env),units),
  > #              data=DT)
  > # summary(mix1)$varcomp[,1]
  > # mix1$Beta
  > # 
  > # # > summary(mix1)$varcomp[,1]
  > # # [1] 0.76660279 0.29901471 0.22166682 0.31243923 0.01923234 0.42072739 0.24321083 0.27761802
  > # # [9] 0.12513727
  > # 
  > # mix2 <- mmec(fixed=Yield~Env-1,
  > #               random= ~ vsc(usr(Env),ids(Name)),
  > #               rcov= ~ vsc(dsr(Env),ids(units)),
  > #               return.param = F,
  > #              # emweight = rep(1,30),
  > #              nIters = 30,
  > #               data=DT)
  > # mix2$theta
  > # plotMonitor(mix2$monitor, cex=1)
  > # mix2$sigma
  > # mix2$monitor[,12]
  > # plot(summary(mix1)$varcomp[,1],mix2$sigma)
  > # 
  > # zz <- with(DT, vsr(dsr(Env),Name))
  > # 
  > # Z <- zz$Z
  > # 
  > # Zind <- rep(1,length(Z))
  > # 
  > # A <- list( diag(41) )#rep(list(diag(41)),4)
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # tt = ((unsm(3)/5) + diag(.8,3,3)) * 10 ;# tt[lower.tri(tt)]=0;tt
  > # thetaI <- list(
  > #   tt,
  > #   diag(10,3,3)
  > # );thetaI
  > # 
  > # ttc= unsm(3);ttc[lower.tri(ttc)]=0;ttc
  > # thetaC <- list(
  > #   ttc,
  > #   diag(1,3,3)
  > # );thetaC
  > # 
  > # X <- model.matrix(~Env, data=DT)
  > # 
  > # y <- as.matrix(DT$Yield)
  > # 
  > # DTx <- DT; DTx$units <- as.factor(1:nrow(DTx))
  > # ss <- with(DTx, vsr(dsr(Env),units) )
  > # 
  > # S <- ss$Z #list(diag(length(y)))
  > # 
  > # Sind <- rep(1,3)
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # 
  > # ## apply the function
  > # weightInf=rep(1,50);weightInf # weights for the information matrix
  > # weightEmInf=rep(1,50)
  > # weightEmInf=c(seq(1,.1,-.1),rep(0,50));weightEmInf # weights for the EM information matrix
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=TRUE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # summary(mix1)$varcomp
  > # res3$monitor
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # 
  > # 
  > # library(sommer)
  > # data(DT_technow)
  > # DT <- DT_technow
  > # Md <- Md_technow
  > # Mf <- Mf_technow
  > # Md <- (Md*2) - 1
  > # Mf <- (Mf*2) - 1
  > # Ad <- A.mat(Md)
  > # Af <- A.mat(Mf)
  > # DT$GY <- scale(DT$GY)
  > # ####=========================================####
  > # ####=========================================####
  > # mix1 <- mmer(GY~1,
  > #              random=~vsr(dent,Gu=Ad) + vsr(flint,Gu=Af),
  > #              rcov=~units,
  > #              data=DT)
  > # summary(mix1)$varcomp[,1]
  > # 
  > # Adi <- as(solve(Ad + diag(ncol(Ad))*1e-6 ),Class="sparseMatrix")
  > # Afi <- as(solve(Af + diag(ncol(Af))*1e-6 ),Class="sparseMatrix")
  > # 
  > # mix2 <- mmec(GY~1,
  > #              random=~vsc(ids(dent),Gu=Adi) + vsc(ids(flint),Gu=Afi),
  > #              rcov=~units, return.param = F,
  > #              data=DT)
  > # mix2$monitor
  > # mix2$sigma
  > # 
  > # z1=model.matrix(~dent-1, data=DT); colnames(z1) <- gsub("dent","",colnames(z1))
  > # z2=model.matrix(~flint-1, data=DT); colnames(z2) <- gsub("flint","",colnames(z2))
  > # 
  > # Z <- list(z1,z2)
  > # 
  > # Zind <- 1:2
  > # 
  > # A <- list(
  > #   Ad[colnames(z1),colnames(z1)], 
  > #   Af[colnames(z2),colnames(z2)]
  > # ) 
  > # A <- lapply(A, function(x){x + diag(1e-3,ncol(x),ncol(x))})
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # thetaI <- list(
  > #   matrix(10,1,1),
  > #   matrix(10,1,1),
  > #   matrix(10,1,1)
  > # );thetaI
  > # 
  > # thetaC <- list(
  > #   matrix(1,1,1),
  > #   matrix(1,1,1),
  > #   matrix(1,1,1)
  > # );thetaC
  > # 
  > # X <- model.matrix(~1, data=DT)
  > # 
  > # y <- as.matrix(DT$GY)
  > # 
  > # S <- list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # thetaF
  > # ## apply the function
  > # weightInf=rep(1,40);weightInf # weights for the information matrix
  > # weightEmInf=c(seq(.9,.1,-.2),rep(0,36));weightEmInf # weights for the EM information matrix
  > # # v2=rep(1,40)
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=FALSE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # summary(mix1)$varcomp
  > # res3$monitor
  > # 
  > # dim(res3$Ci)
  > # 
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # 
  > # library(sommer)
  > # data(DT_ige)
  > # DT <- DT_ige
  > # Af <- A_ige
  > # An <- A_ige
  > # DT$trait <- scale(imputev(DT$trait))
  > # ### Direct genetic effects model
  > # mix1 <- mmer(trait ~ block,
  > #              random = ~vsr(focal,Gu=Af) + vsr(neighbour,Gu=An),
  > #              rcov = ~ units,
  > #              data = DT)
  > # summary(mix1)$varcomp
  > # 
  > # Afi <- as(solve(Af + diag(ncol(Af))*1e-6 ),Class="sparseMatrix")
  > # Ani <- as(solve(An + diag(ncol(An))*1e-6 ),Class="sparseMatrix")
  > # 
  > # mix2 <- mmec(trait ~ block,
  > #              random = ~vsc(ids(focal),Gu=Afi) + vsc(ids(neighbour),Gu=Ani),
  > #              rcov = ~ units, return.param = F,
  > #              data = DT)
  > # mix2$monitor
  > # plotMonitor(mix2$monitor)
  > # 
  > # g=(unsm(2)*-.15) + (diag(2)*.4)
  > # g
  > # mix3 <- mmec(trait ~ block,
  > #              random = ~vsc(ids(focal),ids(neighbour),Gu=Afi, meN=2, meThetaC = unsm(2)),# + vsc(ids(neighbour),Gu=Ani),
  > #              rcov = ~ units, return.param = F, emweight = rep(1,20),
  > #              tolParConv = 1e-6, nIters=40,
  > #              data = DT)
  > # mix3$monitor
  > # mix3$sigma
  > # mix3$rTermsNames
  > # plotMonitor(mix3$monitor)
  > # 
  > # 
  > # z1=model.matrix(~focal-1,data=DT); colnames(z1) <- gsub("focal","",colnames(z1))
  > # z2=model.matrix(~neighbour-1,data=DT);  colnames(z2) <- gsub("neighbour","",colnames(z2))
  > # 
  > # Z <- list(z1,z2)
  > # 
  > # Zind <- 1:2
  > # 
  > # A <- list( 
  > #   Af[colnames(z1),colnames(z1)],
  > #   An[colnames(z2),colnames(z2)]
  > # )
  > # 
  > # A <- lapply(A, function(x){x + diag(1e-3,ncol(x),ncol(x))})
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # tt = diag(1)*10000
  > # thetaI <- list(
  > #   tt,
  > #   tt,
  > #   diag(10000,1,1)
  > # );thetaI
  > # 
  > # ttc=diag(1)
  > # thetaC <- list(
  > #   ttc,
  > #   ttc,
  > #   diag(1,1,1)
  > # );thetaC
  > # 
  > # X <- model.matrix(~block, data=DT)
  > # 
  > # y <- as.matrix(DT$trait)
  > # 
  > # S <- list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # thetaF
  > # 
  > # ## apply the function
  > # weightInf=rep(1,40);weightInf # weights for the information matrix
  > # weightEmInf=c(seq(.9,.1,-.1),rep(0,36));weightEmInf # weights for the EM information matrix
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=TRUE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # # compare results
  > # summary(mix1)$varcomp
  > # res3$monitor
  > # 
  > # 
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # # UNSTRUCTURED MODEL
  > # 
  > # library(sommer)
  > # data(DT_ige)
  > # DT <- DT_ige
  > # Af <- A_ige
  > # An <- A_ige
  > # DT$trait <- imputev(DT$trait)
  > # ### Direct genetic effects model
  > # xx <- c(rep(1,4),rep(0,40))
  > # mix1 <- mmer(trait ~ block,
  > #              random = ~ gvs(focal, neighbour, Gu=list(Af,Af)), 
  > #              rcov = ~ units, iters=10,tolpar = 1e-4,
  > #              data = DT, emupdate = xx)
  > # summary(mix1)$varcomp
  > # 
  > # z1=model.matrix(~focal-1,data=DT); colnames(z1) <- gsub("focal","",colnames(z1))
  > # z2=model.matrix(~neighbour-1,data=DT);  colnames(z2) <- gsub("neighbour","",colnames(z2))
  > # 
  > # Z <- list(z1,z2)
  > # 
  > # Zind <- rep(1,2)
  > # A <- list( 
  > #   Af
  > # )
  > # 
  > # A <- lapply(A, function(x){x + diag(1e-3,ncol(x),ncol(x))})
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # tt = ((-unsm(2)/5) + diag(.8,2,2)) * 1000 ; #tt[lower.tri(tt)]=0;tt
  > # # tt = diag(2)*10000
  > # thetaI <- list(
  > #   tt,
  > #   diag(10000,1,1)
  > # );thetaI
  > # 
  > # ttc= unsm(2);ttc[lower.tri(ttc)]=0;ttc
  > # # ttc=diag(2)
  > # thetaC <- list(
  > #   ttc,
  > #   diag(1,1,1)
  > # );thetaC
  > # 
  > # X <- model.matrix(~block, data=DT)
  > # 
  > # y <- as.matrix(DT$trait)
  > # 
  > # S <- list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # thetaF
  > # 
  > # ## apply the function
  > # 
  > # weightInf=c(rep(1,100));weightInf
  > # weightEmInf=c(rep(1,4),seq(1,.1,-.2),rep(0.05,3),rep(0.05,30));weightEmInf # weights for the EM information matrix
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=TRUE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # # compare results
  > # summary(mix1)$varcomp
  > # res3$monitor
  > # 
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # # UNSTRUCTURED MODEL
  > # 
  > # library(orthopolynom)
  > # data(DT_legendre)
  > # DT <- DT_legendre
  > # DT$Y <- scale(imputev(DT$Y))
  > # mix1<-mmer(Y~ 1 + Xf
  > #            , random=~ vsr(usr(leg(X,1)),SUBJECT)
  > #            , rcov=~vsr(units), tolpar = 1e-6
  > #            , data=DT)
  > # summary(mix1)$varcomp
  > # 
  > # mix2<-mmec(Y~ 1 + Xf
  > #            , random=~ vsc(usr(leg(X,1)),ids(SUBJECT))
  > #            , rcov=~units,return.param = F
  > #            , data=DT)
  > # mix2$monitor
  > # mix2$sigma
  > # 
  > # xx=with(DT,vsr(usr(leg(X,1)),SUBJECT))
  > # 
  > # Z <- xx$Z[c(1,3)]
  > # 
  > # Zind <- rep(1,2)
  > # 
  > # A <- list(
  > #   diag(100)
  > # ) #
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # tt = ((unsm(2)/5) + diag(.8,2,2))  ; #tt[lower.tri(tt)]=0;tt
  > # thetaI <- list(
  > #   tt,
  > #   diag(1,1,1)
  > # );thetaI
  > # 
  > # ttc= unsm(2);ttc[lower.tri(ttc)]=0;ttc
  > # # ttc=diag(2)
  > # thetaC <- list(
  > #   ttc,
  > #   diag(1,1,1)
  > # );thetaC
  > # 
  > # X <- model.matrix(~1 + Xf, data=DT)
  > # 
  > # y <- as.matrix(DT$Y)
  > # 
  > # S <- list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # thetaF
  > # 
  > # ## apply the function
  > # weightInf=rep(1,40) # weights for the information matrix
  > # weightEmInf=rep(1,40)
  > # weightEmInf=c(seq(1,.1,-.2),rep(0,36));weightEmInf # weights for the EM information matrix
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=TRUE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # # compare results
  > # summary(mix1)$varcomp
  > # res3$monitor
  > # 
  > # 
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # 
  > # data(DT_cornhybrids)
  > # DT <- DT_cornhybrids
  > # DT<-DT[which(!is.na(DT$Yield)),]
  > # nrow(DT)
  > # DTi <- DTi_cornhybrids
  > # GT <- GT_cornhybrids
  > # hybrid2 <- DT # extract cross data
  > # A <- GT
  > # K1 <- A[levels(hybrid2$GCA1), levels(hybrid2$GCA1)]; dim(K1)
  > # K2 <- A[levels(hybrid2$GCA2), levels(hybrid2$GCA2)]; dim(K2)
  > # S <- kronecker(K1, K2) ; dim(S)
  > # rownames(S) <- colnames(S) <- levels(hybrid2$SCA)
  > # 
  > # hybrid2$Yield <- scale(hybrid2$Yield)
  > # 
  > # mix1 <- mmer(Yield ~ Location,
  > #              random = ~ vsr(GCA1,Gu=K1) + vsr(GCA2,Gu=K2),# + vsr(SCA,Gu=S),
  > #              rcov=~units,
  > #              data=hybrid2)
  > # summary(mix1)$varcomp
  > # 
  > # K1i <- as(solve(K1 + diag(ncol(K1))*1e-6 ),Class="sparseMatrix")
  > # K2i <- as(solve(K2 + diag(ncol(K2))*1e-6 ),Class="sparseMatrix")
  > # Si <- as(solve(S + diag(ncol(S))*1e-6 ),Class="sparseMatrix")
  > # 
  > # mix2 <- mmec(Yield ~ Location,
  > #              random = ~ vsc(ids(GCA1),Gu=K1i) + vsc(ids(GCA2),Gu=K2i),# + vsX(ids(SCA),Gu=Si),
  > #              rcov=~units, return.param = F,
  > #              data=hybrid2)
  > # mix2$monitor
  > # mix2$sigma
  > # 
  > # z1 <- model.matrix(~GCA1-1,data=DT);colnames(z1) <- gsub("GCA1","",colnames(z1))
  > # z2 <- model.matrix(~GCA2-1,data=DT);colnames(z2) <- gsub("GCA2","",colnames(z2))
  > # z3 <- model.matrix(~SCA-1,data=DT);colnames(z3) <- gsub("SCA","",colnames(z3))
  > # 
  > # Z <- list(
  > #   z1,z2,z3
  > # )
  > # 
  > # Zind <- 1:3
  > # 
  > # A <- list(
  > #   K1[colnames(z1),colnames(z1)], K2[colnames(z2),colnames(z2)], S[colnames(z3),colnames(z3)]
  > # ) #
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # thetaI <- list(
  > #   diag(1)*10,
  > #   diag(1)*10,
  > #   diag(1)*10,
  > #   matrix(50,1,1)
  > # );thetaI
  > # 
  > # thetaC <- list(
  > #   diag(1),
  > #   diag(1),
  > #   diag(1),
  > #   matrix(1,1,1)
  > # );thetaC
  > # 
  > # X <- model.matrix(~Location, data=DT)
  > # 
  > # y <- as.matrix(DT$Yield)
  > # 
  > # S <- list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # thetaF
  > # 
  > # ## apply the function
  > # weightInf=rep(1,40) # weights for the information matrix
  > # weightEmInf=rep(1,40)
  > # weightEmInf=c(seq(.9,.5,-.1),rep(0,36));weightEmInf # weights for the EM information matrix
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp2(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=TRUE,
  > #                       nIters=20, tolParConv=1e-5,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # # compare results
  > # summary(mix1)$varcomp
  > # res3$monitor
  > # 
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # ###############################################################################
  > # 
  > # library(sommer)
  > # 
  > # data(DT_h2)
  > # DT <- DT_h2
  > # 
  > # head(DT)
  > # length(table(DT$Env))
  > # nrow(DT)
  > # DT$y <- scale(DT$y)
  > # 
  > # ans1 <- mmer(y~Env,
  > #              random=~vsr(dsr(Env),Name) + vsr(dsr(Env),Block),
  > #              rcov=~vsr(dsr(Env),units),
  > #              data=DT)
  > # summary(ans1)$varcomp
  > # 
  > # ans2 <- mmec(y~Env,
  > #              random=~vsc(dsr(Env),ids(Name)) + vsc(dsr(Env),ids(Block)),
  > #              rcov=~vsc(dsr(Env),ids(units)), return.param = F,
  > #              # emweight = c(1,rep(0,20)),
  > #              # stepweight = c(1,1,1,rep(0,30)),
  > #              nIters=30,
  > #              data=DT)
  > # 
  > # plotMonitor(ans2$monitor, cex=.1)
  > # ans2$sigma
  > # 
  > # d=ans2$monitor
  > # plotMonitor(d)
  > # plot(d[,(ncol(d)-1)], summary(ans1)$varcomp[,1])
  > # 
  > # DT2 <- DT[with(DT, order(Env)), ]
  > # library(asreml)
  > # ans3 <- asreml(y~Env,
  > #               random=~diag(Env):Name + diag(Env):Block,
  > #               residual=~dsum(~units | Env),
  > #               
  > #               data=DT2)
  > # 
  > # plot(d[,(ncol(d)-1)], summary(ans3)$varcomp[,1])
  > # 
  > # 
  > # xx=with(DT,vsr(dsr(Env),Name))
  > # 
  > # Z <- xx$Z
  > # 
  > # Zind <- rep(1,length(Z))
  > # 
  > # A <- list(
  > #   
  > #   diag(41)
  > # ) #
  > # 
  > # Ai <- lapply(A, function(x){solve(x)})
  > # 
  > # thetaI <- list(
  > #   diag(15)*5,#+diag(rnorm(15)),
  > #   matrix(10,1,1)
  > # );thetaI
  > # 
  > # thetaC <- list(
  > #   diag(15),
  > #   matrix(1,1,1)
  > # );thetaC
  > # 
  > # X <- model.matrix(~Env, data=DT)
  > # 
  > # y <- as.matrix(DT$y)
  > # 
  > # S <- list(diag(length(y)))
  > # 
  > # H <- diag(length(y))
  > # 
  > # addScaleParam <- 0
  > # nn <- unlist(lapply(thetaC, function(x){length(which(x > 0))}))
  > # nn2 <- sum(nn[1:max(Zind)])
  > # ff <- diag(nn2)
  > # thetaF <- cbind(ff,matrix(0,nn2,1))
  > # thetaF
  > # 
  > # ## apply the function
  > # weightInf=rep(1,40) # weights for the information matrix
  > # 
  > # weightEmInf =c(seq(1,.1,-.1),rep(.01,30))
  > # # weightEmInf =rep(.01,30)
  > # # weightInfEMv=c(rep(1,8),rep(0,36));weightInfEMv # weights for the EM information matrix
  > # 
  > # 
  > # Z <- lapply(Z, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # S <- lapply(S, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # Ai <- lapply(Ai, function(x){
  > #   return(as(x, Class = "sparseMatrix"))
  > # })
  > # 
  > # H <- as(H, Class = "sparseMatrix")
  > # 
  > # X <- as(X, Class = "sparseMatrix")
  > # 
  > # y <- as(y, Class = "sparseMatrix")
  > # 
  > # tt=system.time( 
  > #   expr = res3<-ai_mme_sp(X=X,Z=Z, Zind=Zind, 
  > #                       Ai=Ai,y=y,
  > #                       S=S, 
  > #                       H=H, useH=FALSE,
  > #                       nIters=15, tolParConv=1e-4,
  > #                       tolParInv=1e-6,thetaI=thetaI, 
  > #                       thetaC=thetaC,thetaF=thetaF, 
  > #                       addScaleParam=addScaleParam, weightEmInf = weightEmInf, 
  > #                       weightInf = weightInf,
  > #                       verbose=TRUE
  > #                       
  > #   )
  > # )
  > # 
  > # plot(res3$monitor[,ncol(res3$monitor)], summary(ans1)$varcomp[,1] )
  > # 
  > # plot(res3$llik[1,])
  > # 
  > # #######################################################
  > # 
  > # data(DT_yatesoats)
  > # DT <- DT_yatesoats
  > # head(DT)
  > # DT$Y <- scale(DT$Y)
  > # m3 <- mmer(fixed=Y ~ V + N + V:N,
  > #            # random = ~ B + B:MP,
  > #            rcov=~units,
  > #            data = DT)
  > # summary(m3)$varcomp
  > # 
  > # m4 <- mmec(fixed=Y ~ V + N + V:N,
  > #            # random = ~ B + B:MP,
  > #            rcov=~units, return.param = F,
  > #            data = DT)
  > # m4$sigma
  > # plotMonitor(m4$monitor)
  > # str(m4)
  > # 
  > # 
  > # ####################################################
  > # 
  > # library(sommer)
  > # data("DT_cpdata")
  > # DT <- DT_cpdata
  > # GT <- GT_cpdata
  > # MP <- MP_cpdata
  > # DT$Yield <- imputev(DT$Yield)
  > # #### create the variance-covariance matrix
  > # K <- A.mat(GT) # additive relationship matrix
  > # K <- K + diag(1e-3,nrow(K))
  > # #### look at the data and fit the model
  > # head(DT)
  > # mix1 <- mmer(Yield~1,
  > #              random=~vsr(id,Gu=K)
  > #              + Rowf,
  > #              rcov=~units,
  > #              data=DT)
  > # summary(mix1)$varcomp
  > # 
  > # # Ki <- as(solve(K), Class="sparseMatrix")
  > # m <- GT[,1:300]
  > # mix2 <- mmec(fixed=Yield~1,
  > #               random=~vsc(ids(m)) + Rowf,
  > #               rcov=~units, return.param = F,
  > #               data=DT)
  > # mix2$monitor
  > # summary(mix1)$varcomp
  > # 
  > # 
  > # 
  > # ?overlay
  > # 
  > # data("DT_halfdiallel")
  > # DT <- DT_halfdiallel
  > # head(DT)
  > # DT$femalef <- as.factor(DT$female)
  > # DT$malef <- as.factor(DT$male)
  > # DT$genof <- as.factor(DT$geno)
  > # 
  > # A <- diag(7); colnames(A) <- rownames(A) <- 1:7;A # if you want to provide a covariance matrix
  > # #### model using overlay
  > # modh <- mmer(sugar~1, 
  > #              random=~vsr(overlay(femalef,malef), Gu=A) 
  > #              + genof,
  > #              data=DT)
  > # 
  > # Ai <- as(solve(A), Class="sparseMatrix")
  > # modh <- mmec(sugar~1, 
  > #              random=~vsc(ids(overlay(femalef,malef)), Gu=Ai) 
  > #              + genof, #return.param = T,
  > #              data=DT)
  > # modh$monitor
  > # 
  > # #####################################################
  > # 
  > # ?DT_cpdata
  > # 
  > # data(DT_cpdata)
  > # DT <- DT_cpdata
  > # GT <- GT_cpdata
  > # MP <- MP_cpdata
  > # #### create the variance-covariance matrix
  > # A <- A.mat(GT) # additive relationship matrix
  > # Ai <- as( solve(A+diag(1e-5,ncol(A),ncol(A))), Class = "sparseMatrix")
  > # 
  > # head(DT)
  > # DT <- DT[,-c(7:8)]
  > # DT$color <- as.vector(scale(DT$color))
  > # DT$Yield <- as.vector(scale(DT$Yield))
  > # head(DT)
  > # DT2 <- reshape(DT, idvar = "id", varying = list(5:6),
  > #                v.names = "y", direction = "long", timevar = "trait", times =colnames(DT)[5:6] )
  > # DT2$trait <- as.factor(DT2$trait)
  > # head(DT2)
  > # 
  > # g=diag(2)*.05 + matrix(.1,2,2);g
  > # g2 <- diag(2)*.75;g2
  > # mix1 <- mmec(y~trait-1,
  > #              random=~vsc(usr(trait, theta = g),ids(id),Gu=Ai),
  > #              stepweight = ss,
  > #              emweight =rep(1,30),
  > #              return.param = F, tolParConv = 1e-6,
  > #              rcov=~vsc(dsr(trait),ids(units)),
  > #              data=DT2)
  > # 
  > # mix1$theta
  > # mix1$thetaC
  > # mix1$sigma
  > # mix1$b
  > # plot(mix1$u)
  > # plot(mix1$llik[1,])
  > # 
  > # # > summary(mix2)$varcomp
  > # # VarComp  VarCompSE    Zratio Constraint
  > # # u:id.color-color    0.7557346 0.15068202  5.015426   Positive
  > # # u:id.color-Yield    0.0848963 0.07171699  1.183768   Unconstr
  > # # u:id.Yield-Yield    0.1394332 0.07012037  1.988484   Positive
  > # # u:units.color-color 0.3779553 0.04195655  9.008254   Positive
  > # # u:units.Yield-Yield 0.8808832 0.07522796 11.709519   Positive
  > # # # 
  > # 
  > 
  > proc.time()
     user  system elapsed 
     0.12    0.07    0.18