R/AF.Echidna.R

Defines functions esr_poisson esr_binomial esr_gaussian rho.par write.relationshipMatrix .AGH.inv AGH.inv .GenomicRel GenomicRel vm2r0 vm2r vmat met.vmat met.biplot met.corr met.plot converge.esR converge summary.esR summary Var.AR Var.esR Var trace.esR trace get.IC IC.esR IC waldT wald.esR wald coeff11 coef.esR coef predict0 predict.esR predict raneff.acc model.comp11 model.comp output sig.level pin22 pin33 pin11 pin.esR pin plot1f plot.esR plot evp.res esr.res esRT0 esRT vcms.fun b2s update.esR update init.tr run.mod echidna get.es0.file

Documented in AGH.inv b2s coef.esR converge.esR echidna esRT GenomicRel get.es0.file IC.esR met.biplot met.corr met.plot met.vmat model.comp pin pin11 pin22 pin.esR plot plot.esR predict.esR raneff.acc rho.par summary trace.esR update.esR Var Var.esR wald.esR

## version: public

# update: 2023-03-01
#' @title Summary of added functions for Echidna
#' 
#' @param dat.file	 data file to generate .es file.
#' @param es.file	   the .es file to generate .es0 file.
#' @param es0.file	 the .es0 file.
#' @param object	   Echidna result object in R.
#' @param path    the path for data files.	
#' @param softp   the path for Echidna software.
#  @param update  update for Echidna software.
#' @param trace	  show iteration procedure,FALSE(default). 
#' @param maxit   maximum number of iterations, 30(default).
#' @param Fmv     make missing values into fixed terms, FALSE(default).
#' @param mu.delete     delete term mu or Trait from model, FALSE(default).
#' @param mulT	  multi-trait model,FALSE(default).
#' @param met	    multi-environment trial model,FALSE(default).
#' @param mulN	  trait number for multi-trait analysis at one time, 2(default).
#' @param mulp	  multi-pin formula to run at one time, NULL(default). 
#' @param cycle	  Echidna result from qualifier cycle,FALSE(default).
#' @param trait	   aim trait for analysis, such as, 'h3', 'h3 h4',~h3+h4, etc, NULL(default).
#' @param family  such as esr_binomial(), esr_poisson().
#' @param selfing  the probability of selfing for parent, such as 0.1.
#' @param weights   A variable used as weights in the fit.
#' @param fixed      fixed effects, such as, c('Rep'), c('Site', 'Site.Rep') or 'Site Site.Rep', h3~1+Rep, etc.	
#' @param random	   random effects, such as,'Mum','Mum Mum.Rep',~Mum+Mum:Rep, etc. 
#' @param residual	 residual effects, such as,'units','ar1(row).ar1(col)',~ar1(row):ar1(col), etc.
#' @param batch     run batch analysis for more than 2 trait at one time, FALSE(default).
#' @param batch.G   run more than 2 G structures at one time, FALSE(default). 
#' @param batch.R   run more than 2 R structures at one time, FALSE(default). 
#' @param subF      run subF function for MET data sets,FALSE(default).
#' @param subV.org  original variable for subF.
#' @param res.no    number to show results.
#' @param foldN	    new folder name to store each run's results, only works when delf is 'FALSE'.
#' @param delf      delete all Echidna result files from the folder of .es0 file, TRUE(default).	
#' @param message	  show running procedure,FALSE(default). 
#' @param run.purrr  using purrr packages for batch analysis,FALSE(default).
#' @param predict	   prediction for model terms. 
#' @param vpredict	 run vpredict statements with Echidna soft.
#' @param jobqualf	 header line qualifiers, mainly '!view'. 
#' @param qualifier	 model qualifiers, such as '!extra 5'.
#' 
#' @details
#'   This package would supply some functions for Echidna. Details as following:
#' \tabular{ll}{
#' \strong{Function} \tab \strong{Description} \cr
#' \code{get.es0.file}  \tab generate .es0 file. \cr
#' \code{echidna}    \tab run mixed models. \cr
#' \code{wald}       \tab output wald results. \cr
#' \code{Var}        \tab output variance components. \cr
#' \code{summary}    \tab output summary results. \cr
#' \code{IC}         \tab output AIC and BIC values. \cr
#' \code{pin}        \tab run pin functions.\cr
#' \code{predict}    \tab output predict results.\cr
#' \code{plot}       \tab output model diagnose results.\cr
#' \code{coef}       \tab output fixed and random effects.\cr
#' \code{update}     \tab update mixed models.\cr
#' \code{b2s}        \tab transform batch esR results to single esR.\cr
#' \code{model.comp} \tab Model comparison for different mixed models.
#' }
#'
#' @author Yuanzhen Lin <yzhlinscau@@163.com>
#' @references
#'  Yuanzhen Lin. R & ASReml-R Statistics. China Forestry Publishing House. 2016 \cr
#'  Gilmour, A.R. (2020) Echidna Mixed Model Software www.EchidnaMMS.org
#' @name AF.Echidna
#' @examples
#' \dontrun{
#'
#'  library(AFEchidna)
#'
#'  ##  Echidna
#'  path='D:/Echidna/Jobs'
#'
#'  setwd(path)
#'  
#'  ## generate .es0 file
#'  get.es0.file(dat.file='fm.csv')
#'  get.es0.file(es.file='fm.es')
#'  # file.edit('fm.es0')
#'
#' res<-echidna(trait='h3',
#'               fixed='Rep',random='Fam',
#'               residual=NULL,predict=c('Fam'),
#'               es0.file="fm.es0")
#' 
#' ## method 2                           
#' # res<-echidna(fixed=h3~1+Rep,random=~Fam,
#' #              residual=NULL,predict=c('Fam'),
#' #              es0.file="fm.es0")
#'
#'  names(res)
#'  class(res)
#'
#'  # model diagnose
#'  plot(res) 
#'
#'  # wald result
#'  wald(res)
#'  waldT(res, term=c('mu','Rep'))
#'  
#'
#'  # variance components
#'  Var(res)
#'
#'  # summary result
#'  summary(res)
#'
#'  # AIC,BIC result
#'  IC(res)
#'
#'  # fixed and random effects
#'  coef(res)$fixed
#'  coef(res)$random
#'
#'  # predict results if using predict functions
#'  mm<-predict(res)
#'  mm$pred
#'
#'  # show vc results by using vpredict statements
#'  pin(res)
#'
#'  # run pin function to count genetic parameters
#'  pin11(res,h2~V1/(V1+V2))
#'  pin(res,mulp=c(h2~V1/(V1+V2),h2f~V1/(V1+V2/4)))
#'
#'  # model converge stage
#'  trace(res)
#'  res$Converge
#'
#' }
NULL

############### following functions to input results from Echidna

#' @rdname  AF.Echidna
#' @usage get.es0.file(dat.file=NULL,es.file=NULL,
#'                         path=NULL,message= FALSE,
#'                         softp=NULL,
#'                         faS=NULL,pedS=NULL,Rsuffix=FALSE)
#' @export
get.es0.file <- function(dat.file=NULL,es.file=NULL,path=NULL, 
                       message = FALSE,softp=NULL,#update=FALSE,
                       faS=NULL,pedS=NULL,Rsuffix=FALSE) { 
  
  if(is.null(softp)) Echsf<-AFEchidna::loadsoft()
  if(!is.null(softp)) Echsf<-softp
   
  if(!is.null(dat.file)){ #nextp==FALSE
    
    if (!is.null(path))  setwd(path)
    
    if (message == TRUE)  mess1 <- FALSE
    else mess1 <- TRUE
    
    if(!file.exists(dat.file)) stop('data file does not exist.')
     else {
      esjob <- dat.file# esjob<-'fm.csv'
      runes <- paste(Echsf, esjob, sep = " ")
      #ifelse(.Platform$OS.type == "windows",
      #   system(runes, show.output.on.console = TRUE, wait = FALSE, 
      #       invisible = mess1),
      #   system(runes))
      system(runes, show.output.on.console = TRUE, wait = FALSE,invisible = mess1)
      #system2(Echsf, args=esjob, wait = FALSE, invisible = mess1)
      
      #flst <- dir()
      #temp <- flst[grep("\\.es$", flst)]
      #temp<-sub('.csv','.es',dat.file)
      #if(file.exists(temp))
        cat("Generating .es temple for",dat.file,": ", "--done!\n") # temp,
      #else stop('Generating .es temple fails.')
    } 
    
  }
  if(!is.null(es.file)){
    
    if(!file.exists(es.file)) 
      stop('es file does not exist.\n Error reason: Echidna may not work.')
    
    if (!is.null(path)) path0<-path 
    else path0<-getwd()
    setwd(path0)
    
    #flst<-dir()
    #if(is.null(es.file)) es.file <- flst[grep("\\.es$", flst)]
    #else es.file<-es.file#"barley.es" 
    
    new_dir<-'tempd'
    dir.create(new_dir)
    file.copy(es.file,new_dir, overwrite=TRUE)
    
    path1<-paste(path0,new_dir,sep='/')
    
    setwd(path1)
    tempf <- base::readLines(es.file) #es.file<-"dm.es"
    tempf <- sub("\\!DOPART \\$1", "#!DOPART $1", tempf)
    
    tempf <- sub(" correctly classified", "", tempf)
    
    ## data fields
    ## head 4 lines: jobq+title+header+data1
    ## variable: 4+orgS
    ## faS example: 1~6; faS=1:6
    
    #if(!is.null(pedS)) tempf[4+pedS]<-sub('  #','!P #',tempf[4+pedS])
    #if(!is.null(faS)) tempf[4+faS]<-sub('  #','!I #',tempf[4+faS])
    #if(!is.null(AfaS)) tempf[4+AfaS]<-sub('  #','!A #',tempf[4+AfaS])
    
    ## data fields
    #if(is.null(pedS) ){
      b<-NULL
      b[1]<-4
      b[2]<-which(grepl("(# Verify data.*)", tempf))
      datf<-tempf[4:b[2]]
      
      c<-NULL
      if(!is.null(faS)) c <- faS+1
        else c<-which(grepl("(^[A-Z])", datf))
      c0<-which(grepl("(.*!A.*)", datf))
      c<-c[!c %in% c0]
      tempf[3+c]<-sub('  #','!I #',tempf[3+c])
    #}
    
    a <- NULL
    a[1] <- which(grepl("(Y ~ mu Fixed.*)", tempf))
    a[2] <- which(grepl("(residual.*)", tempf))
    tempf <- tempf[-a]
    
    dat.fl <- which(grepl("(!SKIP)", tempf))
    if(!is.null(pedS)) {
      tempf[4+pedS]<-sub(' !I #',' !P #',tempf[4+pedS])
      tempf <- append(tempf,tempf[dat.fl])
      tempf[length(tempf)-1] <- paste0(tempf[length(tempf)-1], ' !PARTIALSELF 0')
    }
      
    
    flst <- dir()
    flst <- gsub("\\.es", "\\.es0", flst)
    temp <- flst[grep("\\.es0$", flst)]
    if(Rsuffix) temp <- gsub("\\.es0", "\\.es0.R", temp)
    write(tempf, file = temp)
    
    file.copy(temp,path0,overwrite = TRUE)
    #file.remove(dir())
    
    setwd(path0)
    unlink(new_dir,recursive =TRUE, force=TRUE)
    
    esf<-dir(pattern='.es$')
    file.remove(esf)
    
    if(!is.null(pedS)) 
      message('Make sure there has a correct pedigree file in the .es0 file.\n')
    
    cat("Generating .es0 file:", "-- done!\n") # temp,
    
    #invisible()
  }
  
}

########## main function echidna()
#' @rdname  AF.Echidna
#' @usage echidna(fixed,random,residual,
#'                    trait,family,weights, 
#'                    es0.file,softp,  
#'                    delf,foldN,
#'                    trace,maxit,
#'                    Fmv,mu.delete,
#'                    mulT,met,cycle,
#'                    batch,mulN,mulp,
#'                    batch.G,batch.R,
#'                    subF,subV.org,
#'                    res.no,dat.file,
#'                    run.purrr,selfing,
#'                    predict,vpredict,
#'                    qualifier,jobqualf) 
#' @export
echidna <- function(fixed=NULL,random=NULL,residual=NULL,
                  trait=NULL,family=NULL,#trait.mod=NULL,
                  weights=NULL,
                  es0.file,softp=NULL,
                  delf=TRUE,foldN=NULL,
                  trace=TRUE,maxit=30,
                  Fmv=FALSE,mu.delete=FALSE,
                  mulT=FALSE,met=FALSE,cycle=FALSE,
                  batch=FALSE,mulN=NULL,mulp=NULL, 
                  batch.G=FALSE,batch.R=FALSE,
                  subF=FALSE,subV.org=NULL,
                  res.no=NULL,dat.file=NULL,
                  run.purrr=FALSE,selfing=NULL,
                  predict=NULL,vpredict=NULL,
                  qualifier=NULL,jobqualf=NULL){
  
  require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
  
  if(!is.null(trait)){
    #trait=~h3+h4+h5
    #trait=c(~h3,~h4,~h5)
    #trait=c('h3','h4','h5')
    if(class(trait)=="formula") {
      trait1 <- as.character(trait)
      #trait1<-trait1[trait1!='~']
      trait <- strsplit(trait1[2],'\\+')[[1]]
    }
    
    if(class(trait)=="list") {
      #trait1<-list();length(trait1)<-length(trait)
      trait1 <- vector("list", length(trait))
      trait1 <- lapply(trait, as.character)
      #trait1<-as.character(trait)
      if(length(trait1[[1]])==3){
        trait0 <- unlist(lapply(trait1, function(x) x[2]))
        trait  <- unlist(lapply(trait1, function(x) x[3]))
        #trait1<-gsub('\\+',' ',trait1)
      }
      if(length(trait1[[1]])==2){
        ran0 <- NA
        trait <- unlist(lapply(trait1, function(x) x[2]))
        #trait1<-gsub('\\+',' ',trait1)
      } 
    }    
  }
                               
  # get data file, maybe problem here!!!
  es0.txt <- base::readLines(es0.file)
  datL <- es0.txt[grep('\\![Ss][Kk][Ii][Pp]',es0.txt)] # >=1
  if(grepl('\\#',datL)) datL <- datL[-grep('\\#',datL)]
  lth <- length(datL)
  dat.file0 <- sub('\\s+\\![Ss][Kk][Ii][Pp].*','',datL[lth])
  if(is.null(dat.file)) dat.file <- dat.file0 else dat.file <- dat.file                               
  
  test<-function(mode=c("batch.Y", "batch.G", "batch.R", 'subF' )){
    tt<-NULL
    mode <- switch(mode, "batch.Y" = 1, "batch.G" = 2,
                    "batch.R" = 3,'subF' = 4 )
    
    if (as.numeric(mode)==1) {
      
      if(batch==FALSE){ # AFEchidna::  ###  batch.Y, batch.G, batch.R
        ttN <- 1
        batch0 <- FALSE
      } else{ 
        # batch
        require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
        require(tidyr,warn.conflicts=FALSE,quietly=TRUE) # unite
        batch0 <- TRUE
        
        cat('\nProgram starts running batch analysis ------ \n')
        
        if(trace==FALSE & !is.null(mulp)){
          cat('\npin formula: \n');for(i in 1:length(mulp)) print(mulp[[i]])
          #cat('\n')
        }
        
        cc <- trait; ccN <-length(cc)
        if(mulT==FALSE){ ttN <- trait}
        
        if(mulT==TRUE){
          if(is.null(mulN)) mulN <- 2
          
          if((ccN/mulN)>1) { bb <- utils::combn(cc,mulN);bbn <- ncol(bb)}
          if((ccN/mulN)==1){ bb <- cc;bbn <- 1}
          if((ccN/mulN<1)){ stop("\nThe trait No is less than in the model!\n")}
          
          ttN <- 1:bbn
        }
      }
      
      run.fun1 <- function(x){
        
        if(trace==TRUE & batch==TRUE) cat('\nrun',x,'-- -- --:')
        if(batch==FALSE) {
          x <- trait
        }
        
        if(batch==TRUE){
          if(mulT==FALSE) x <- x else {
            if(bbn>1)  x <- paste(bb[,x],collapse=' ')
            if(bbn==1) x <- paste(bb,collapse=' ')
          }
        }
        
        AFEchidna::run.mod(es0.file=es0.file,softp=softp,
                           trait=x,family=family,#trait.mod=trait.mod,
                           weights=weights,
                           fixed=fixed,random=random,residual=residual,
                           mulT=mulT,met=met,selfing=selfing,
                           trace=trace,delf=delf,foldN=foldN,
                           predict=predict,vpredict=vpredict,
                           maxit=maxit,cycle=cycle,
                           Fmv=Fmv,mu.delete=mu.delete,
                           qualifier=qualifier,jobqualf=jobqualf)
       } 
        
      if(!run.purrr) ss <- lapply(ttN, run.fun1)
        else ss <- ttN %>% purrr::map( run.fun1 )
      
      #
      if(batch==FALSE) {
        ss <- ss[[1]]
        if(cycle==FALSE) NTrait <- ss$Traits
        if(cycle==TRUE) NTrait <- names(ss$esr.all)
      }
      
      if(batch==TRUE){
        if(mulT==FALSE) names(ss) <- trait
        if(mulT==TRUE) names(ss) <- unlist(lapply(ss,function(x) x$Traits))
        
        NTrait <- names(ss) 
      }
      NTrait <- gsub(' ','-',NTrait)
      NTrait <- sub('-$','',NTrait)
      
      #tt <- tt1 <- NULL
      
      if(batch==FALSE) {
        tt <- ss
      }
      if(batch==TRUE){
        tt$res.all <- ss
        Converge <- sapply(tt$res.all,function(x) x$Converge)
        maxit <- sapply(tt$res.all,function(x) x$maxit)
      }
    }
    
    # for 2 more G-structures
    if (as.numeric(mode)==2) {
      batch0 <- TRUE
      
      cat('\nProgram runs for 2 more G-structure at one time. ------ \n')
      
      random1 <- lapply(random, as.character)
      if(length(random1[[1]])==3){
        ran0 <- unlist(lapply(random1, function(x) x[2]))
        ran1 <- unlist(lapply(random1, function(x) x[3]))
        ran1 <- gsub('\\+',' ',ran1)
      }
      if(length(random1[[1]])==2){
        ran0 <- NA
        ran1 <- unlist(lapply(random1, function(x) x[2]))
        ran1 <- gsub('\\+',' ',ran1)
      }
      
      ttN <- length(ran1)
      if(is.na(ran0)) ran0 <- paste0('G',1:ttN)
        
      run.fun2 <- function(x){
        #x=1
        if(trace==TRUE) cat('\nrun',x,'-- random effects:', ran1[x])
        # 
        if(!is.na(ran1[x]))
          AFEchidna::run.mod(es0.file=es0.file,softp=softp,
                  fixed=fixed,random=ran1[x],residual=residual,
                  mulT=mulT,met=met,selfing=selfing,
                  trace=trace,delf=delf,foldN=foldN,
                  predict=predict,vpredict=vpredict,
                  maxit=maxit,cycle=cycle,
                  Fmv=Fmv,mu.delete=mu.delete,
                  qualifier=qualifier,jobqualf=jobqualf)
        else
          AFEchidna::run.mod(es0.file=es0.file,softp=softp,#trait=x,
                             fixed=fixed,random=NULL,residual=residual,
                             mulT=mulT,met=met,selfing=selfing,
                             trace=trace,delf=delf,foldN=foldN,
                             predict=predict,vpredict=vpredict,
                             maxit=maxit,cycle=cycle,
                             Fmv=Fmv,mu.delete=mu.delete,
                             #batch0=batch0,
                             qualifier=qualifier,jobqualf=jobqualf)
      
      }
      
      if(!run.purrr) ss <- lapply(1:ttN, run.fun2) 
        else ss <- 1:ttN %>% purrr::map( run.fun2 )       
      names(ss) <- ran0     
      #tt <- NULL      
      tt$res.all <- ss 
    }
    
    # for 2 more R-strucutre
    if (as.numeric(mode)==3) {
      batch0 <- TRUE
      cat('\nProgram runs for 2 more R-structure at one time. ------ \n')
      
      residual1 <- lapply(residual, as.character)
      if(length(residual1[[1]])==3){
        resid0 <- unlist(lapply(residual1, function(x) x[2]))
        resid1 <- unlist(lapply(residual1, function(x) x[3]))
        resid1 <- gsub('\\+',' ',resid1)
      }
      if(length(residual1[[1]])==2){
        resid0 <- NA
        resid1 <- unlist(lapply(residual1, function(x) x[2]))
        resid1 <- gsub('\\+',' ',resid1)
      }
      
      ttN <- length(resid1)
      if(is.na(resid0)) resid0 <- paste0('R',1:ttN)
      
        
      run.fun3 <- function(x){
        #x=1
        if(trace==TRUE) cat('\nrun',x,'-- residual effects:', resid1[x])

        AFEchidna::run.mod(es0.file=es0.file,softp=softp,
                           fixed=fixed,random=random,residual=resid1[x],
                           mulT=mulT,met=met,selfing=selfing,
                           trace=trace,delf=delf,foldN=foldN,
                           predict=predict,vpredict=vpredict,
                           maxit=maxit,cycle=cycle,
                           Fmv=Fmv,mu.delete=mu.delete,
                           qualifier=qualifier,jobqualf=jobqualf)

      }
      
      if(!run.purrr) ss <- lapply(1:ttN, run.fun3 ) 
      else ss <- 1:ttN %>% purrr::map(run.fun3 )
      
      names(ss) <- resid0      
      tt <- NULL      
      tt$res.all <- ss
    }

        # subF function
    if (as.numeric(mode)==4) {
      batch0 <- TRUE
      
      cat('Starting analysis.\n')
      
      # data file
      dat <- read.csv(file=dat.file)
      
      # copy original data file and rename to an old.file
      org.datf <- paste0('old.',dat.file)
      write.csv(dat,file=org.datf,row.names=FALSE)
      dat <- read.csv(file=org.datf)
      
      dat$nSite <- dat[,subV.org]
      cc <- unique(dat$nSite)
      
      if(is.null(mulN)) mulN <- 2 
      else mulN <- mulN
      
      bb <- utils::combn(cc,mulN)
      if(is.null(res.no)) bbn <- ncol(bb) 
      else bbn <-res.no
     
      run.fun4 <- function(x){
        cat('\nAnalysing---- ',paste(append(paste0('Site-'),bb[,x]), collapse = ":"))
        
        dat22 <- dat %>% filter(.,nSite %in% bb[,x])
        write.csv(dat22,file=dat.file,row.names=FALSE)
        
        mm <- AFEchidna::run.mod(fixed=fixed,
                      random=random,
                      residual=residual,
                      #qualifier = subsetcc,
                      trace=trace,met=met,
                      es0.file = es0.file)
      }
      
      ss<-vector("list", bbn)
      if(!run.purrr) ss <- lapply(1:bbn, run.fun4 ) 
      else ss <- 1:bbn %>% purrr::map(run.fun4 )
      
      names(ss)<- lapply(1:bbn, function(x) paste(append(paste0('Site-'),bb[,x]), collapse = ":"))
      
      dat$nSite <- NULL
      write.csv(dat,file=dat.file,row.names=FALSE)
      file.remove(org.datf) 
      unlink(org.datf,force=TRUE)                   
      cat('works done.\n')
      
      #tt <- NULL      
      tt$res.all <- ss      
    }
                         
    call <- list(fixed=fixed,random=random,residual=residual)    
    org.par <- list(es0.file=es0.file,softp=softp,
                    trait=trait,family=family,#trait.mod=trait.mod,
                    weights=weights,selfing=selfing,
                    fixed=fixed,random=random,residual=residual,
                    mulT=mulT,mulN=mulN,mulp=mulp,
                    met=met,trace=trace,delf=delf,
                    Fmv=Fmv,mu.delete=mu.delete,
                    cycle=cycle,call=call,
                    run.purrr=run.purrr,
                    batch0=batch0,batch=batch,
                    batch.G=batch.G,batch.R=batch.R,
                    subF=subF,subV.org=subV.org, 
                    #subV.Lv=subV.Lv,subV.new=subV.new,
                    res.no=res.no,dat.file=dat.file,
                    predict=predict,vpredict=vpredict,
                    qualifier=qualifier,jobqualf=jobqualf)                      
    
    tt$org.par <- org.par
    
    class(tt) <- c('esR') 
    
    return(tt)
  }
  
  res<-NULL
  if(subF==FALSE){
    if(batch.G==FALSE & batch.R==FALSE)  res <- test(mode="batch.Y")
    if(batch.G==TRUE  & batch.R==FALSE)  res <- test(mode="batch.G")
    if(batch.R==TRUE  & batch.G==FALSE)  res <- test(mode="batch.R")
  }
  if(subF==TRUE) res <- test(mode="subF")
      
  return(res)
}


#' @export
run.mod <- function(es0.file,softp=NULL,
                    trait=NULL,family=NULL,#trait.mod=NULL,
                    weights=NULL,
                    fixed=NULL,random=NULL,residual=NULL,
                    delf=TRUE,foldN=NULL,selfing=NULL,
                    trace=TRUE,maxit=30,
                    Fmv=TRUE,mu.delete=FALSE,
                    mulT=FALSE,met=FALSE,cycle=FALSE,
                    batch=FALSE, 
                    predict=NULL,vpredict=NULL,
                    qualifier=NULL,jobqualf=NULL) {
  
  if(is.null(softp))  Echsf <- AFEchidna::loadsoft()
  if(!is.null(softp)) Echsf <- softp
  
  #setwd(es0.path)
  flst0 <- dir()

  if(is.null(es0.file)) stop('Needs an .es0 file.\n')
  # if(met==TRUE) mulT <- TRUE
  
  #flst0 <- flst
  
  if(class(fixed)=="formula") {
    fixed1 <- as.character(fixed)
    if(length(fixed1)==2) # fixed=~Rep
      fixed <- gsub('\\+','',fixed1[2])
    if(length(fixed1)==3){ 
      # trait0<-sub('cbind\\(','',fixed1[2])
      # trait0<-sub('\\)$','',trait0)
      
      #fixed=cbind(h3,h4)~Trait+Trait:Rep
      patt <- '(?<=\\().+?(?=\\))' # match text in the (...)
      trait0 <- stringr::str_extract(string =fixed1[2], pattern = patt)
      if(is.na(trait0)) trait0 <- fixed1[2] # fixed=h3~1+Rep+plot
      if(is.null(family)) trait0 <- gsub(',',' ',trait0)
      #family=c(esr_binomial(),esr_binomial())
      
      if(!is.null(family)){
        trait0 <- strsplit(trait0,', ')[[1]]
        trait00<-NULL
        for(i in 1:length(trait0))
          trait00[i]<-paste0(trait0[i],family[i])
        trait0 <- paste0(trait00,collapse = ' ')
      }
      
      fixed <- gsub('\\+','',fixed1[3])
      fixed <- sub('^1\\s','',fixed)
      fixed <- sub('^Trait\\s','',fixed)
    }
  }
  
  if(is.null(trait)) trait <- trait0
   else trait <- trait
  
  #i=1;
  if(trace==TRUE) cat('\nRunning Echidna for analysis: ',trait,'\n')
  
  qualF<-NULL;qualF[1]<-' '
  qualF[2]<-' !SLN !YHT '
  qualF[3]<-paste0(' !maxit ',maxit, ' ')
  if(!is.null(qualifier)) qualF[4]<-qualifier
  qualF<-paste(qualF,sep=' ')
  
  estxt<-NULL;
  #if(!is.null(qualifier)) estxt[1]<-qualifier else estxt[1]<-''
  if(cycle==TRUE) {
    estxt[1]<-paste('\n!cycle', paste(trait,collapse = " "), sep=' ')
    #estxt[1]<-paste0(estxt[1], cytxt)
  } else estxt[1]<-''
  estxt[2]=''
  estxt<-paste(estxt,sep=' ')
  
  # linear model
  lmtxt<-NULL
  
  if(mulT==FALSE) {
    if(cycle==FALSE){
      if(is.null(fixed)) {lmtxt[1]<-paste0(trait,'~ mu ,') 
      } else lmtxt[1]<-paste0(paste(trait,paste(fixed,collapse = " "),sep='~ mu '),' ,') # fix
    } else {
      if(is.null(fixed)){ lmtxt[1]<-paste0(paste('$I',' ~ mu ,'))
      } else lmtxt[1]<-paste0(paste('$I',paste(fixed,collapse = " "),sep=' ~ mu '),' ,')                                   
    } 
  }else{
    if(is.null(fixed)){ lmtxt[1]<-paste0(trait,' ~ Trait ,')
    } else lmtxt[1]<-paste0(paste(trait,paste(fixed,collapse = " "),sep=' ~ Trait '),' ,')
  }
  
  # if(met==TRUE) lmtxt<-gsub('Trait',' mu ',lmtxt)
  lmtxt<-gsub('~ mu 1','~ mu ',lmtxt)
  
  ## weights
  if(!is.null(weights))
    lmtxt<-gsub('~ mu',paste0(' !WT ',weights,' ~ mu '),lmtxt)

  ## family
  # if(!is.null(family))
  #  lmtxt<-gsub('~ mu',paste0(' ',family,' ~ mu '),lmtxt)
  
  
  if(class(random)=="formula"){
    random1 <- as.character(random)
    
    if(!grepl('init',random1[2])){
      random <- gsub('\\+','',random1[2])
      if(grepl('\\*',random)) random <- gsub(' ','',random)
    } else { # random=~us(Trait,init=c(.1,.1,.1)):Fam
      random <- AFEchidna::init.tr(random1)
    }
  } 
  if(is.null(random)) {
    lmtxt[2] <- '!r '
  }else lmtxt[2] <- paste0('!r ',paste(random,collapse = " ")) # ran
  
  if(Fmv==TRUE) lmtxt[2] <- paste0(lmtxt[2], '  !f mv')
  
  if(class(residual)=="formula"){
    residual1 <- as.character(residual)
    
    if(!grepl('init',residual1[2])){
      residual <- gsub('\\+','',residual1[2])
    } else { # residual=~units:us(Trait,init=c(.1,.1,.1))
      residual <- AFEchidna::init.tr(residual1)
    }
  } 
  
  
  if(is.null(residual)) {# resid
    lmtxt[3]<- ' ' #'residual units'
  }  else lmtxt[3]<-paste0('residual ',residual)
  
  #### !!!!! spatial
  lmtxt<-gsub('ar1v','ar1',lmtxt) 
  lmtxt<-gsub('idv','id',lmtxt)


  #### !!!!! delete mu or Trait
  if(mu.delete==TRUE) {
      if(mulT==TRUE) lmtxt<-gsub('~ Trait','~ ',lmtxt)
      if(mulT==FALSE) lmtxt<-gsub('~ mu',' ~ ',lmtxt)
  }
  
  # predict
  predtxt<-NULL
  if(!is.null(predict)){

    predtxt<-sapply(1:length(predict), function(x) paste('predict ',predict[x],' '))
  } 
  predtxt<-c('',predtxt)
  
  # vpredict
  vptxt<-vptxt2<-NULL
  vptxt[1]<-''
  vptxt[2]<-'vpredict'
  vptxt[3]<-'W components'
  if(!is.null(vpredict)){
    vpredict <- gsub('\\:xfa','.xfa',vpredict)
    vpredict <- gsub('\\:fa','.fa',vpredict)
    vptxt2<-sapply(1:length(vpredict), function(x) vpredict[x])
  }
  vptxt<-c(vptxt,vptxt2)
  
  file.copy(es0.file,'temp.es',overwrite=TRUE)
  
  estxt<-c(qualF,estxt,lmtxt,predtxt)#,vptxt)
  estxt0<-gsub(':','.',estxt)
  estxt<-c(estxt0,vptxt)
                   
  if(!is.null(selfing)) {
    #es0.file='pine_provenance.es0'
    tempf <- base::readLines(es0.file)
    #selfing=0.5
    selfq<-paste0('!PARTIALSELF ',selfing)
    tempf<-gsub('!PARTIALSELF 0',selfq,tempf)
    utils::write.table(tempf,file='temp.es',quote=FALSE,
                       row.names=FALSE,col.names=FALSE,append=FALSE)
  }else file.copy(es0.file,'temp.es',overwrite=TRUE)                 
                   
  utils::write.table(estxt,file='temp.es',quote=FALSE,
              row.names=FALSE,col.names=FALSE,append=TRUE)
  
  if(!is.null(jobqualf)){
    write(jobqualf,file='temp')
    file.append('temp','temp.es')
    file.copy('temp','temp.es',overwrite=TRUE)
    file.remove('temp')
    unlink('temp',force=TRUE)
  }
  
  #dir()
  esjob<-'temp.es'
  runes<-paste(Echsf,esjob,sep=' ')
  
  ifelse(.Platform$OS.type == "windows",
         system(runes,show.output.on.console=FALSE), # run program
         system(runes))
  
  es0.path<-getwd()
  df<-NULL
  df<-AFEchidna::esRT(es0.path,trace=FALSE,mulT=mulT,met=met,
                      cycle=cycle,vpredict=vpredict)
  
  ## Iteration procedure
  if(trace==TRUE) {
    if(cycle==FALSE){
      cat('\n',df$StartTime,'\n')
      if(met==TRUE | !is.null(family)) cat(df$Iterations00,'\n')
        else print.data.frame(df$Iterations)
      cat(df$FinishAt,'\n\n')
    }
    if(cycle==TRUE){
      nn<-length(df$esr.all)
      trt<-unlist(lapply(df$esr.all,function(x) x$Traits))
      
      for(i in 1:nn){
        cat('\n\nIteration procedure for trait: ',trt[i],'\n')
        cat('\n',df$esr.all[[i]]$StartTime,'\n')
        print.data.frame(df$esr.all[[i]]$Iterations)
        cat(df$esr.all[[i]]$FinishAt,'\n')
      }
    }
    if(batch==TRUE){
      
      trt<-unlist(lapply(df$res.all,function(x) x$Traits))
      
      trt<-gsub(' ','-',trt)
      trt<-sub('-$','',trt)
      
      
      for(i in 1:length(trt)) {
        
        cat('\n\nIteration procedure for trait: ',trt[i],'\n')
        cat('\n',df$res.all[[i]]$StartTime,'\n')
        print.data.frame(df$res.all[[i]]$Iterations)
        cat(df$res.all[[i]]$FinishAt,'\n')
      }
    }
  }
  
  flst<-dir()
  dlst<-base::setdiff(flst,flst0)
  esv<-flst[grep('\\.esv$',flst)]
  dlst2<-dlst[dlst!=esv]
  #dlst<-dlst[!grep('\\.esv$',dlst)]
  
  if(delf==TRUE) file.remove(dlst2)
  if(delf==FALSE) {
    if(!is.null(foldN))  new_dir<-foldN
    if(is.null(foldN)) {
      if(cycle==FALSE) new_dir <-paste(trait,'result',sep='.') else new_dir<-'cyl.result'
    } 
    if(!dir.exists(new_dir)) dir.create(new_dir)
    for(file in dlst) file.copy(file,new_dir, overwrite=TRUE)
    
    file.remove(dlst2)
    #unlink(dlst2,force=TRUE)
  }
  
  dlst0<-dir(".", pattern="^temp")
  dlst0<-dlst0[dlst0!=esv]
  file.remove(dlst0)
  base::on.exit('temp.ess')
  file.remove('temp.ess')
  #unlink(dlst0,force=TRUE)
  #file.remove('fort.13')
  
  return(df)
}

## transform for initial values in G and R
#' @export
init.tr <- function(random1) {
  random2<-strsplit(random1[2],'\\+')[[1]]
  random2<-gsub(',',' ',random2)
  random2<-gsub('init',' ',random2)
  random2<-gsub('\\)\\)','\\)',random2)
  random2<-gsub('c\\(','',random2)
  random2<-gsub('=','',random2)
  random2<-paste(random2,collapse = ' ')
  
  return(random2)
}

## ==========================
#' @export
update <- function(object,trait=NULL,fixed=NULL,
                   random=NULL,residual=NULL,
                   predict=NULL,vpredict=NULL,
                   qualifier=NULL,jobqualf=NULL,
                   trace=NULL,maxit=30,
                   selfing=NULL,mu.delete=FALSE,
                   mulT=NULL,met=NULL,
                   batch=NULL,mulN=NULL, 
                   batch.G=NULL,batch.R=NULL,
                   subF=FALSE,subV.org=NULL,
                   res.no=NULL, dat.file=NULL,                  
                   delf=NULL,foldN=NULL,...){
  UseMethod("update",object)
}
#' @rdname  AF.Echidna
#' @method  update esR
#' @export  update.esR
#' @export
update.esR<-function(object,trait=NULL,fixed=NULL,
                     random=NULL,residual=NULL,
                     predict=NULL,vpredict=NULL,
                     qualifier=NULL,jobqualf=NULL,
                     trace=NULL,maxit=30,
                     selfing=NULL,mu.delete=NULL,
                     mulN=NULL,mulT=NULL,met=NULL,
                     cycle=NULL,softp=NULL,
                     batch=NULL, 
                     batch.G=NULL,batch.R=NULL,
                     subF=FALSE,subV.org=NULL,
                     res.no=NULL,dat.file=NULL,                     
                     delf=NULL,foldN=NULL,...){
  #object<-res21
  org.par<-object$org.par
  
  es0.file<-org.par$es0.file
  
  if(is.null(trait))    trait<-org.par$trait
  if(is.null(fixed))    fixed<-org.par$fixed
  if(is.null(random))   random<-org.par$random
  if(is.null(residual)) residual<-org.par$residual
  
  if(is.null(predict))   predict<-org.par$predict
  if(is.null(vpredict))  vpredict<-org.par$vpredict
  if(is.null(qualifier)) qualifier<-org.par$qualifier
  if(is.null(jobqualf))  jobqualf<-org.par$jobqualf
  if(is.null(selfing))   selfing<-org.par$selfing
  if(is.null(mu.delete))  mu.delete<-org.par$mu.delete

  if(is.null(delf))    delf<-org.par$delf
  if(is.null(trace))   trace<-org.par$trace
  if(is.null(foldN))   foldN<-org.par$foldN
  if(is.null(batch))   batch<-org.par$batch
  if(is.null(mulN))    mulN<-org.par$mulN
  if(is.null(batch.G)) batch.G<-org.par$batch.G
  if(is.null(batch.R)) batch.R<-org.par$batch.R

  if(is.null(subF))     subF<-org.par$subF
  if(is.null(subV.org)) subV.org<-org.par$subV.org
  if(is.null(dat.file)) dat.file<-org.par$dat.file
  if(is.null(res.no))   res.no<-org.par$res.no
  
  if(is.null(mulT))     mulT<-org.par$mulT
  if(is.null(met))      met<-org.par$met
  if(is.null(cycle))    cycle<-org.par$cycle
  if(is.null(softp))    softp<-org.par$softp
  
  AFEchidna::echidna(es0.file=es0.file,
                     softp=softp,trait=trait,
              fixed=fixed,random=random,residual=residual,
              mulT=mulT,met=met,cycle=cycle,
              predict=predict,vpredict=vpredict,
              qualifier=qualifier,jobqualf=jobqualf,
              trace=trace,maxit=maxit,selfing=selfing,
              mu.delete=mu.delete,
              batch=batch,mulN=mulN,
              batch.G=batch.G,batch.R=batch.R,
              subF=subF,subV.org=subV.org, 
              res.no=res.no,dat.file=dat.file,              
              delf=delf,foldN=foldN)
  
}
                                                                     

## batch results to each-single result
#' @usage b2s(object)
#' @rdname  AF.Echidna
#' @export
b2s<-function(object){
  #object<-res22
  
  org.par<-object$org.par
  org.par$batch<-FALSE
  
  data<-object$res.all
  trt<-names(data)
  trtN<-length(trt)
  
  #res<-list();length(res)<-trtN
  res<-vector("list", trtN)
  
  res<-lapply(1:trtN, function(x){
    xx<-data[[x]]
    xx$org.par<-org.par
    class(xx)<-'esR'
    xx
  })
  names(res)<-trt
  #class(res[[1]])
  return(res)
}



#' @export
vcms.fun <- function(aa1,mulN) {
  
  require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
  
  mat1<-diag(mulN)
  df.mat1<-reshape::melt(lower.tri(mat1,diag=TRUE))
  
  for(i in 1:2) {
    #df.mat0[,i]<-paste0('t',df.mat0[,i])
    df.mat1[,i]<-paste0('t',df.mat1[,i])
  }

  
  df.mat11<-df.mat1 %>% dplyr::filter(value==TRUE) %>% 
    dplyr::arrange(X1) %>% tidyr::unite('tt',X1:X2,sep='.',remove=F)
  ttn<-df.mat11$tt
  
  
  ## diag
  ttn0<-paste0('t',1:mulN)
  aa1[grepl('units\\.diag\\(Trait\\)',aa1)]<-paste0('units_',ttn0)# resid
  
  vcn.gen<-aa1[grepl('diag\\(Trait\\)\\.',aa1)]  ## genetic
  vgn<-sub('diag\\(Trait\\)\\.','',vcn.gen[1])
  vgn<-sub(';diag\\(Trait\\)','',vgn)
  aa1[grepl(';diag\\(Trait\\)',aa1)]<-paste(vgn,ttn0,sep='_')## genetic
  
  ## corgh
  aa1[grepl('units\\.corgh\\(Trait\\)',aa1)]<-paste0('units_',ttn0)# resid
  
  vcn.gen<-aa1[grepl('corgh\\(Trait\\)\\.',aa1)]  ## genetic
  vgn<-sub('corgh\\(Trait\\)\\.','',vcn.gen[1])
  vgn<-sub(';corgh\\(Trait\\)','',vgn)
  aa1[grepl(';corgh\\(Trait\\)',aa1)]<-paste(vgn,c('corr',ttn0),sep='_')## genetic
  
  # unstructure: us
  aa1[grepl('units\\.us\\(Trait\\)',aa1)]<-paste0('units_',ttn)# resid
  
  vcn.gen<-aa1[grepl('us\\(Trait\\)\\.',aa1)]  ## genetic
  vgn<-sub('us\\(Trait\\)\\.','',vcn.gen[1])
  vgn<-sub(';us\\(Trait\\)','',vgn)
  aa1[grepl(';us\\(Trait\\)',aa1)]<-paste(vgn,ttn,sep='_')## genetic
  
  return(aa1)
}


#' @rdname  AF.Echidna
#' @usage esRT(path,trace=FALSE,mulT=FALSE,met=FALSE,cycle=FALSE) 
#' @export
esRT<-function(path,trace=FALSE,mulT=FALSE,met=FALSE,cycle=FALSE,vpredict=NULL){
  suppressMessages(esRT0(path=path,trace=trace,mulT=mulT,met=met,
                         cycle=cycle,vpredict=vpredict))
} 
  

#' @export
esRT0 <- function(path,trace=FALSE,mulT=FALSE,met=FALSE,
                  cycle=FALSE,vpredict=NULL) {
  
  if(!require(readr,warn.conflicts=FALSE,quietly=TRUE))
    stop('Need package: readr.\n')
  if(!require(stringr,warn.conflicts=FALSE,quietly=TRUE))
    stop('Need package: stringr.\n')

  require(readr,warn.conflicts=FALSE,quietly=TRUE)

  setwd(path)

  flst<-dir()

# if(any(grepl('_e.R$',flst))){
#   Rresf<-flst[grep('_e.R$',flst)]
#   
#   ff<-readr::read_file(file=Rresf)
#   ff<-paste0('mm=function(){',ff,'}')
#   if(grepl('Vpredict',ff))  ff<-sub('Vpredict','if(0==1) Vpredict',ff)
#   readr::write_file(ff,'ff1.R')
#   tt<-source('ff1.R')$value()
#   file.remove('ff1.R')
#   
#   aa<-strsplit(tt$FinishAt,'LogL ')[[1]][2]
#   if(aa=="Converged") tt$Converge<- TRUE else tt$Converge<- FALSE
# } else 
  tt<-list()

  #tt<-tt[!sapply(tt, is.null)]

  ## input .esr results
  if(length(grep('\\.esr$',flst))==1) {
    esrf<-flst[grep('\\.esr$',flst)]
    esr<-readr::read_file(file=esrf)
    
    if(cycle==FALSE) tt<-AFEchidna::esr.res(esr, mulT=mulT, met=met)
      
    tt$esr<-esr
  }

  ## fixed and random effects
  if(length(grep('\\.ess$',flst))==1) {
    essf<-flst[grep('\\.ess$',flst)]
    if(cycle==TRUE) tt$coef<-base::readLines(essf)
     else {
       
       coef<-base::readLines(essf)
       skipn<-grep('at\\(',coef)
       
       ## new here
       if(length(skipn)!=0) 
        df <- readr::read_lines(file=essf,skip=skipn)
       else df <- readr::read_lines(file=essf)
       
       dfL<-vector('list',length(df))
       
       for(i in 1:length(df))
         dfL[[i]]<-strsplit(df[i],'\\,\\s+')[[1]]
       
       df0<-do.call('rbind',dfL)
       df1<-as.data.frame(df0[-1,])
       names(df1)<-df0[1,]
       for(i in 3:4) df1[,i]<-as.numeric(df1[,i])
       
       tt$coef <- df1
       rm(df,df0,df1,dfL)
       #file.remove(essf)
       #unlink(essf,force=TRUE)
     }
  }


  ## model diagnosis for residuals
  if(length(grep('\\.esy$',flst))==1) {
    esyf<-flst[grep('\\.esy$',flst)]
    if(cycle==TRUE) tt$yht<-base::readLines(esyf) 
     else tt$yht<-utils::read.csv(file=esyf,header=TRUE)
  }else{
    if(mulT==TRUE) 
      warning('.esy file not exists.Please use !YHT to generate it.')
  }

  ## input predict results
  if(length(grep('\\.epv$',flst))==1) {
    epvf<-flst[grep('\\.epv$',flst)] 
    if(cycle==TRUE) tt$pred<-readr::read_file(file=epvf)
     else tt$pred<-base::readLines(epvf)
  }

  ## input vpredict results
  if(length(grep('\\.evp$',flst))==1) {
    evpf<-flst[grep('\\.evp$',flst)]
    evp<-readr::read_file(file=evpf)
    evp0<-base::readLines(evpf)
    tt$evp<-evp
    tt$evp0<-evp0
  }

  ## input var.comp results
  if(length(grep('\\.vpc$',flst))==1) {
    vpcf<-flst[grep('\\.vpc$',flst)]
    vpc<-utils::read.table(file=vpcf,header=F)
    names(vpc)<-c('Vc','Term')
    tt$vpc<-vpc
  }
  
  #flst<-dir()
  ## input variances of var.comp
  if(length(grep('\\.vpv$',flst))==1) {
    vpvf<-flst[grep('\\.vpv$',flst)]
    vpv<-scan(file=vpvf,quiet=TRUE)
    #tt$vpv<-vpv

    if(cycle==FALSE){
      vpv.mat<-diag(1,nrow=nrow(vpc))
      vpv.mat[upper.tri(vpv.mat,diag=TRUE)]<-vpv
      a<-vpv.mat
      a[lower.tri(a)]<-t(a)[lower.tri(a)]
      rownames(a)<-colnames(a)<-paste0('V',1:nrow(vpc))
      vpv.mat<-a
      tt$vpv.mat<-vpv.mat
    } else tt$vpv<-vpv

  }
  
  if(cycle==TRUE){### !cycle
    res<-tt#<-res13
    
    temp<-list()
    
    ##1 esr
    esrr<-strsplit(res$esr,'\\s+?Echidna\\s+')[[1]]
    esrr<-esrr[esrr!=""]
    esrr<-paste('Echidna',esrr,sep=' ')
    
    # trait number
    trtn<-length(esrr)
    
    for(i in 1:trtn) temp[[i]]<-esrr[i]
    
    esr.all <- lapply(temp,esr.res)
    
    # trait names
    trt<-lapply(esr.all,function(x) x$Traits)
    trt<-unlist(trt)
    
    ##2 evp
    temp<-NULL
    
    evpr<-strsplit(res$evp,'\\s+?Warning:\\s+')[[1]]
    evpr<-evpr[evpr!=""]
    
    for(i in 1:trtn) temp[[i]]<-evpr[i]
    evp.all <- lapply(temp,evp.res)
    
    rm(temp)
    
    ##3 vpc
    vpc<-res$vpc
    
    termn<-nrow(vpc)/trtn
    vpc.all<-split(vpc, ceiling(seq_along(vpc$Vc)/termn))
    
    
    ##4 vpv
    vpv<-res$vpv
    matn<-termn*(termn+1)/2
    
    vpv.all<-split(vpv, ceiling(seq_along(vpv)/matn))
    
    names(esr.all)<-trt
    names(evp.all)<-trt
    names(vpc.all)<-trt
    names(vpv.all)<-trt
    
    res$esr.all<-esr.all
    res$evp.all<-evp.all
    res$vpc.all<-vpc.all
    res$vpv.all<-vpv.all
    
    tt<-res
  }

  return(tt)
}


#' @export
esr.res <- function(esr, mulT=FALSE,met=FALSE) {
    require(stringr,warn.conflicts=FALSE,quietly=TRUE)
  
    tt<-list()
    
    #esr<-res$esr
    
    vpatt<-'(Echidna\\s+.*Windows)'
    tt$Version<-stringr::str_extract(string =esr, pattern = vpatt)[[1]]
    
    jpatt<-'(?<=TITLE:\\s).+(.*)'
    tt$job<-stringr::str_extract(string =esr, pattern = jpatt)[[1]]
    
    #cat(heads)
    st.patt <- "(?<=at\\s).+(\\d+:\\d+:\\d+\\s\\d+)"
    tt$StartTime <- stringr::str_extract(string =esr, pattern = st.patt)[[1]]
    
    It.patt <- "(\\s+.*LogL=.*\\sDF.*)"
    Iter <- stringr::str_extract_all(string =esr, pattern = It.patt)[[1]]
    #cat(Iter)
    tt$Iterations00<-Iter
    Iterv<-as.numeric(unlist(regmatches(Iter,
                                        gregexpr("[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?",
                                                 Iter, perl=TRUE))))
    #if(family ==" !poisson !log !disp 1")
    
    if(mulT==FALSE){
      Iterations <- data.frame(matrix(Iterv,ncol=4,byrow=TRUE))
      names(Iterations)<-c("Iteration","LogL","eSigma","NEDF")
    }
    if(mulT==TRUE|met==TRUE){
      Iterations <- data.frame(matrix(Iterv,ncol=3,byrow=TRUE))
      names(Iterations)<-c("Iteration","LogL","NEDF")
    }
    
    tt$Iterations<-Iterations
    
    tt$LogLikelihood<-utils::tail(Iterations$LogL,1)
    tt$Residual.DF <-utils::tail(Iterations$NEDF,1)
    tt$maxit<-nrow(Iterations)
    
    IC.patt<-'(Akaike\\s+.*)|(Bayesian\\s+.*)'
    IC<-stringr::str_extract_all(string =esr, pattern = IC.patt)[[1]]
    dd<-as.numeric(unlist(regmatches(IC,
                                     gregexpr("[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?",
                                              IC, perl=TRUE))))
    tt$AkaikeIC<-dd[1];tt$BayesianIC<-dd[3]
    tt$ICparameter.Count<-dd[2]
    
    Tr.patt<-'(?<=Analysis of\\s).+(.*)'
    tt$Traits<-stringr::str_extract_all(string =esr, pattern = Tr.patt)[[1]]
    
    #if(crayon) cat(blue(heads)) else cat(heads)
    
    heads<-strsplit(esr,'\n Akaike Information Criterion')[[1]][1]
    
    mids<-strsplit(esr,'Akaike Information Criterion')[[1]][2]
    mids1<-paste0(' Akaike Information Criterion',mids)
    #cat(mids1)
    
    mids2<-strsplit(mids1,'Finished: ')[[1]][1]
    mids2a<-strsplit(mids2,'\r\nCovariance')[[1]][1]
    #cat(mids2a)
    tt$keyres<-mids2a
    
    # wald parts
    mids2b<-strsplit(mids2a,'\n Model_Term')[[1]][1]
    tt$waldT<-paste0('\t\t\tWald',strsplit(mids2b,'Wald')[[1]][2])
    #cat(tt$waldT)
    
    # varcomp parts
    mids3<-strsplit(mids2a,'\n Model_Term')[[1]][2]
    mids3a<-paste0(' Model_Term',mids3)
    #cat(mids3a)
    tt$components<-mids3a
    
    vc.patt<-"[-+]?\\d+\\.\\d+.*"
    vct<-stringr::str_extract_all(string =mids3a, pattern = vc.patt)[[1]]

    
    ft.patt <- '(?<=Finished:\\s).+(.*Converged)'
    tt$FinishAt<-stringr::str_extract(string =esr, pattern = ft.patt)[[1]]
    
    #tt='Finished: Wed Dec 30 13:05:02 2020 LogL NOT Converged  teaching'
    cg.patt<-"(LogL.*Converged)"
    if(!is.na(tt$FinishAt)){
       aa<-stringr::str_extract(string =tt$FinishAt, pattern = cg.patt)[[1]]
       if(aa=="LogL Converged") tt$Converge<- TRUE else tt$Converge<- FALSE
    } # else tt$Converge<-NA
    
  return(tt)
}  

#' @export
evp.res <- function(evp) {
  require(stringr,warn.conflicts=FALSE,quietly=TRUE)
  #evp<-gsub('Vmat 2','Vmat.2',evp)
  #evp<-mm$evp
  evp1<-evp
  
  evp1<-strsplit(evp1,'W components')[[1]][1]
  #cat(evp1)
  
  pattern<-"(\\s+.\\d+.\\w+.*\\d+.)"
  xx<-stringr::str_extract_all(string =evp1, pattern = pattern)[[1]]
  df<-utils::read.table(text=xx,header=FALSE)

  names(df)<-c('NO','Term','Sigma','SE')
  
  return(df)
}
  
#======================================================
#' @title output model diagnose results
#'
#' @description
#' \code{plot} This function output model diagnose results.
#'
#' @details
#' Test trait's norm for Echidna object,similar to asreml.
#' @aliases plot
#' @param object an object of Echidna-R result.
#' @param idx  trait order(1,2,...) when use !cycle in Echidna.
#  @param mulT	  multi-trait model,FALSE(default). 
#' @param meanN  make hat(y) higher than 4*mean(y) to NA.
#'
# @export
# plot <- function(object,...) {
#   UseMethod("plot")
# }
#' @return the result is returned directly.
#' @author Yuanzhen Lin <yzhlinscau@@163.com>
#' @references
#' Yuanzhen Lin. R & ASReml-R Statistics. China Forestry Publishing House. 2016 
#' @name  plot
#' @examples
#'
#' \dontrun{
#' library(AFEchidna)
#'
#'  ##  Echidna
#'  path='D:/Echidna/Jobs/METb11'
#'
#'  #  mainly works for '_e.R'
#'  res<-esRT(path=path,trace=T) # for single trait
#'
#'  names(res)
#'  class(res)
#'
#'  # check .esy exist or not
#'  plot(res)
#' }
#'
#' @export
plot <- function(object,idx=NULL,meanN=2.5){ # mulT=FALSE,
  UseMethod("plot",object)
}
#' @method  plot esR
#' @export  plot.esR
#' @rdname  plot
#' @export
#'
plot.esR <- function(object,idx=NULL,meanN=2.5) { # mulT=FALSE,
  #object<-res
  
  if(object$org.par$batch==FALSE){
    mulT<-object$org.par$mulT
    if(is.null(object$yht)) stop('Please use !view in Echidna.')
    
    if(mulT==TRUE){# mulT=T, each trait y: y-resid, x: y-hat
      plot1f(object,mulT=TRUE,meanN=2.5)
    }else {
      
      if(is.null(object$esr.all)){ # for single-trait
        if(!is.null(object$Traits)){ 
          title<-paste0('Plots for trait: ',object$Traits)
          
          plot1f(object$yht,title=title,meanN=meanN)
        }
      }else{ # for !cycle single-trait
        df<-list()
        
        trt<-unlist(lapply(object$esr.all,function(x) x$Traits))
        
        yhtt<-object$yht
        aa<-which(grepl('(Record,\\s+.*)', yhtt))
        
        for(i in 1:(length(aa)-1)){
          df[[i]]<-utils::read.table(text=yhtt[aa[i]:(aa[i+1]-1)],header=TRUE,sep=',')
          attr(df[[i]],'terms')<-trt[i]
        }
        df[[length(aa)]]<-utils::read.table(text=yhtt[aa[length(aa)]:length(yhtt)],
                                     header=TRUE,sep=',')
        attr(df[[length(aa)]],'terms')<-trt[length(aa)]
        
        names(df)<-trt
        
        if(is.null(idx)) lapply(df,plot1f) 
         else lapply(df[idx],plot1f)
        
      }
    }
  }
  
  if(object$org.par$batch==TRUE){
    mulT<-object$org.par$mulT
    
    dat<-lapply(object$res.all,function(x) x$yht)
    trt<-unlist(lapply(object$res.all,function(x) x$Traits))
    names(dat)<-trt
    trtN<-length(trt)
    ## single trait batch
    if(mulT==FALSE){
      
      title<-paste0('Plots for trait: ',trt)
      if(is.null(idx)) lapply(1:trtN,function(x) plot1f(dat[[x]],title=title[x]))
      else plot1f(dat[[idx]],title=title[idx])
    }else{
      
      #dat0<-list();length(dat0)<-trtN
      dat0 <- vector("list", trtN)
      for(i in 1:trtN){
        dat0[[i]]<-list(yht=dat[[i]],Traits=trt[i])
      }
      
      if(is.null(idx)) lapply(1:trtN,function(x) plot1f(dat0[[x]],mulT=TRUE,meanN=2.5))
      else plot1f(dat0[[idx]],mulT=TRUE,meanN=2.5)
    }
  }
  
}



#' @export
plot1f <- function(object,title=NULL,mulT=FALSE,meanN=2.5) {
  
  dat<-object
  
  if(mulT==FALSE){ # single trait
    if(!is.null(attr(dat,'terms'))) 
      title<-paste0('Plots for trait: ',attr(dat,'terms'))
    else title<-title
    
    graphics::par(mfrow=c(2,2))
    
    dat$HatValue[dat$HatValue>meanN*mean(dat$HatValue,na.rm=T)]<-NA
    
    graphics::hist(dat$Residual,main='',xlab='Residuals',col='blue')
    
    stats::qqnorm(dat$Residual,main='',col='blue',ylab='Residuals')
    
    graphics::plot(Residual~HatValue,data=dat,xlab='Fitted',ylab='Residuals',col='blue')
    graphics::abline(h=0)
    
    graphics::plot(Residual~Record,data=dat,xlab='Unit Number',ylab='Residuals',col='blue')
    graphics::abline(h=0)
    
    if(!is.null(title)) graphics::mtext(title,side=3,line=-02,outer=TRUE)
    
    graphics::par(mfrow=c(1,1))
  } else{ # multT==T

    df<-dat$yht
    
    trt<-dat$Traits
    trt<-strsplit(trt,' ')[[1]]
    trtN<-length(trt)
    trtL<-trt#paste0('T',1:trtN)
    
    df$trait<-trtL

    library(dplyr,warn.conflicts=FALSE,quietly=TRUE)
    
    #df0<-list();length(df0)<-trtN
    df0 <- vector("list", trtN)
    for(i in 1:trtN){
      df0[[i]]<-df %>% dplyr::filter(trait==trtL[i]) %>% 
        dplyr::mutate(HatValue = ifelse(HatValue>meanN*mean(HatValue,na.rm=TRUE), NA, HatValue))
    }
    df<-do.call(rbind,df0)
    
    require(ggplot2,warn.conflicts=FALSE,quietly=TRUE)
    
    print(ggplot2::ggplot(df,ggplot2::aes(x=HatValue,y=Residual,colour=trait))+
            ggplot2::geom_point()+ggplot2::xlab('Fitted values')+
            ggplot2::facet_grid(.~trait,scales='free'))
    
  }

}
## =====================================

#' @title Count error for h2 and corr.
#'
#' @description
#' \code{pin} This function counts standard error(se) for heritability(h2) and
#' corr value and also outputs significent level for corr value in asreml and breedR package.
#'
#' @details
#' Count error for h2 and corr value, also outputs significent level.
#' @aliases  pin
#' @param object	 Echidna result object in R.
#' @param formula	 formula for h2 or corr.
#' @param idx	    trait order to run, 1,2...
#' @param mulp	  multi-formula for h2, corr, etc.
#' @param signif	 Index to output signif levels, FALSE(default) for non-signif.
#' @param digit	 Index for decimal number, 3(default).
#' @param Rres  Index(TRUE) to restore results, FALSE(default).
#' @param all	 Show variance component and genetic parmameter together, FALSE(default).
#  @param cycle  run models with qualifier !cycle.
#  @param org	 print original results by vpredict statements, FALSE(default).
#' @return the result is returned directly.
#' @author Yuanzhen Lin <yzhlinscau@@163.com>
#' @rdname  pin
#' @examples
#' \dontrun{
#'  library(AFEchidna)
#'
#'  ##  Echidna
#'  path='D:/Echidna/Jobs/METb11'
#'
#'  rm(ls())
#'
#'  #  mainly works for '_e.R'
#'  res<-esRT(path=path,trace=T)# for single trait
#'  #res<-esRT(path=path,trace=T,mulT=T) # for multi-traits
#'  
#'  # pin results if using vpredict function
#'  pin(res)
#'  
#'  # run pin function to count genetic parameters
#'  pin11(res,h2~V1/(V1+V2))
#'  pin(res,mulp=c(h2~V1/(V1+V2),h2f~V1/(V1+V2/4))
#' }
#' @export
pin <- function(object,mulp=NULL,idx=1,digit=3,
                all=FALSE,signif=FALSE,Rres=FALSE){ 
  UseMethod("pin",object)
}
#' @rdname  pin
#' @method  pin esR
#' @export  pin.esR
#' @export
pin.esR <- function(object,mulp=NULL,idx=1,digit=3,
                    all=FALSE,signif=FALSE,Rres=FALSE){
  
  
  if(is.null(mulp)){
    
    if(object$org.par$batch==FALSE){
      if(!is.null(object$vpc.all)){
        #object<-res13
        
        vpc0<-object$vpc.all[[idx]]
        vpc<-do.call(rbind,object$vpc.all)#[[idx]]
        vpc$vcS<-paste0('V',1:nrow(vpc0))
        
        trt<-lapply(object$esr.all,function(x) x$Traits)
        trt<-unlist(trt)
        vpc$trait<-rep(trt,each=nrow(vpc0))
        
        row.names(vpc)<-NULL
        vpc$Term<-gsub('\\.',':',vpc$Term)
        
        print(vpc)
      } else {
        #object<-res11
        
        dat<-object$vpc #object$evp
        if(is.data.frame(dat)) {
          dat$vcS<-paste0('V',1:nrow(dat))
          dat$Term<-gsub('\\.',':',dat$Term)
          print.data.frame(dat)
        } else cat(dat)
      }
    }
    if(object$org.par$batch==TRUE){
      #object<-res21
      
      vpc0<-object$res.all[[idx]]$vpc
      vpc<-lapply(object$res.all,function(x) x$vpc)
      vpc<-do.call(rbind,vpc)#[[idx]]
      vpc$vcS<-paste0('V',1:nrow(vpc0))
      
      trt<-lapply(object$res.all,function(x) x$Traits)
      trt<-unlist(trt)
      trt<-gsub(' ','-',trt)
      trt<-sub('-$','',trt)
      vpc$trait<-rep(trt,each=nrow(vpc0))
      
      row.names(vpc)<-NULL
      vpc$Term<-gsub('\\.',':',vpc$Term)
      
      print(vpc)
    } 
  } 
    
  if(!is.null(mulp)){
    if(object$org.par$batch==FALSE){
      if(is.null(object$vpc.all)){
        aa <- lapply(mulp,function(x){
          aa1 <- AFEchidna::pin11(object,formula=x,digit=digit,Rres=TRUE) 
        })
        res<-do.call(rbind,aa)
        res<-as.data.frame(res)
        
        res<-AFEchidna::fdata(res,faS=list(0,c(2:3)),FtN=TRUE)
        
        vc<-AFEchidna::Var(object)[,c('Term','Sigma','SE')]
        vc1<-vc
        vc1$vcS<-paste0('V',1:nrow(vc1))
        
        if(Rres==FALSE){
          cat("variance components are as following:\n")
          print.data.frame(vc1)
          cat('\npin formula: \n');for(i in 1:length(mulp)) print(mulp[[i]])
          cat('\n')
        }
        
        names(vc)<-c('Term','Estimate','SE')
        res<-res[,c('Term','Estimate','SE')]
        
        if(all==TRUE) res<-rbind(vc,res)
        
        if(signif==TRUE) res$Siglevel<-AFEchidna::sig.level(res$Estimate,res$SE)
        
        df<-AFEchidna::output(res,signif=signif,digit=digit,Rres=TRUE)
        
        if(Rres==FALSE){
          print.data.frame(df)
          
          if(signif==TRUE) {
            cat("---------------")
            cat("\nSig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1\n\n")
          }
        } else return(df)
        
      }else{ # for !cycle
        #object<-res13
        if(!is.null(object$org.par)) cycle<-object$org.par$cycle
        else cycle<-FALSE
        
        trt<-lapply(object$esr.all,function(x) x$Traits)
        trt<-unlist(trt)
        trtN<-length(trt)
        
        xxx<-lapply(1:trtN, function(x){
          xx<-lapply(mulp,function(y){
            x<-AFEchidna::pin11(object,idx=x,formula=y,digit=digit,all=all,Rres=TRUE) 
          })
          res<-do.call(rbind,xx)
          res<-as.data.frame(res)
          row.names(res)<-NULL
          res<-dplyr::distinct(res,Term,.keep_all = TRUE)
          #res[,-1]<-round(res[,-1],3) # digit=3
          #nres<-ncol(res)
          res<-AFEchidna::fdata(res,faS=list(0,c(2:3)),FtN=TRUE)
          if(signif==TRUE) res$Siglevel<-AFEchidna::sig.level(res$Estimate,res$SE)
          res
        })
        
        if(Rres==FALSE){
          cat('\nresults as following: \n')
          cat('\npin formula: \n');for(i in 1:length(mulp)) print(mulp[[i]])
          cat('\n')
        }
        
        names(xxx)<-trt
        
        vcms<-xxx[[1]][,1] # pin sign
        #aa1<-paste(paste0('V',1:length(aa1)),aa1,sep='-')
        if(cycle==TRUE)
          vcms<-paste(paste0('V',1:length(vcms)),vcms,sep='-',collapse='; ')
        else vcms<-paste(vcms,collapse='; ')
        
        if(Rres==FALSE) cat('terms: ',vcms);cat('\n\n')
        
        if(signif==FALSE)
          df<-do.call(rbind,lapply(xxx,function(x){unlist(x[,2:3])}))
        if(signif==TRUE) 
          df<-do.call(rbind,lapply(xxx,function(x){unlist(x[,2:4])}))
        
        df<-as.data.frame(df)
        if(length(mulp)==1) names(df)[1]<-'Estimate1'
        
        if(cycle==TRUE) names(df)<-gsub('Estimate','V',names(df))
        
        trt<-gsub(' ','-',trt)
        trt<-sub('-$','',trt)
        row.names(df)<-trt
        
        if(Rres==FALSE) {
          print.data.frame(df)
          if(signif==TRUE) {
            cat("---------------")
            cat("\nSig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1\n\n")
          }
        } else return(list(vcms=vcms,res=df))  
      } 
    }
      
  if(object$org.par$batch==TRUE){
      #object<-res22
    mulT<-object$org.par$mulT
    mulN<-object$org.par$mulN
      
      trt<-names(object$res.all)#,function(x) x$Traits)
      trt<-unlist(trt)
      trtN<-length(trt)
      
      ## problem here!!!!
      xxx<-lapply(1:trtN, function(x){
        #x=1
        xx<-lapply(mulp,function(y){ 
          x<-AFEchidna::pin33(fm=object$res.all[[x]],formula=y)
        })
        
        xx<-do.call(rbind,xx)
        
        fm<-object$res.all[[x]]
        #names(fm)
        vc<-AFEchidna::evp.res(fm$evp)
        vc1<-vc[,-1]
        names(vc1)[2]<-'Estimate'
        
        if(all==TRUE) res<-rbind(vc1,xx)
         else res<-xx
        res<-dplyr::distinct(res,Term,Estimate,.keep_all = TRUE)
        res<-AFEchidna::fdata(res,faS=list(0,c(2:3)),FtN=TRUE)
        res[,2:3]<-round(res[,2:3],digits=digit) # first: Term
        if(signif==TRUE) res$Siglevel<-AFEchidna::sig.level(res$Estimate,res$SE)
        res
      })
      
      if(Rres==FALSE){
        cat('\nresults as following: \n')
        cat('\npin formula: \n');for(i in 1:length(mulp)) print(mulp[[i]])
        cat('\n')
      }
      
      names(xxx)<-trt
      
      vcms<-lapply(xxx,function(x) unlist(x[,1]))[[1]]
      if(mulT==TRUE) vcms<-AFEchidna::vcms.fun(vcms,mulN)
      vcms<-paste(paste0('V',1:length(vcms)),vcms,sep='-',collapse='; ')
      
      if(Rres==FALSE) cat('terms: ',vcms);cat('\n\n')
      
      if(signif==FALSE)
        df<-do.call(rbind,lapply(xxx,function(x){unlist(x[,2:3])}))
      else df<-do.call(rbind,lapply(xxx,function(x){unlist(x[,2:4])}))
      
      df<-as.data.frame(df)
      if(length(mulp)==1) names(df)[1]<-'Estimate1'
      
      names(df)<-gsub('Estimate','V',names(df))
      
      trt<-gsub(' ','-',trt)
      trt<-sub('-$','',trt)
      if(mulT==FALSE) trt<-gsub('-','',trt)
      row.names(df)<-trt
      
      if(Rres==FALSE) {
        print.data.frame(df)
        if(signif==TRUE) {
          cat("---------------")
          cat("\nSig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1\n\n")
        }
      } else return(list(vcms=vcms,res=df)) 
    }
  }
}


#' @usage 
#'       pin11(object,formula=NULL,idx=1,signif=FALSE,
#'                 digit=3,all=FALSE,Rres=FALSE)
# @method  pin11
#' @rdname  pin
#' @export  pin11
#' @export
pin11 <-function(object,formula=NULL,idx=1,signif=FALSE,
                       digit=3,all=FALSE,Rres=FALSE){ 

  #object<-res11 

  if(!require(msm,warn.conflicts=FALSE,quietly=TRUE)){stop('Need package: msm.\n')}
  require(msm,warn.conflicts=FALSE,quietly=TRUE)
  
  if(is.null(formula)) {
    if(object$org.par$batch==FALSE){
      # original variance components
      #if(cycle==TRUE){
      if(!is.null(object$vpc.all)){  
        vpc<-object$vpc.all[[idx]]
        vpc$vcS<-paste0('V',1:nrow(vpc))
        print(vpc)
      } else {
        dat<-object$vpc #object$evp
        if(is.data.frame(dat)) {
          dat$vcS<-paste0('V',1:nrow(dat))
          print(dat)
        } else cat(dat)
      }
    }
    
    if(object$org.par$batch==TRUE){
      var1<-AFEchidna::evp.res(data[[1]]$evp)
      var1$vcS<-paste0('V',1:nrow(var1))
      print(var1)
    }
  }
  
  if(!is.null(formula)) {
    #object<-res11
    
    if(object$org.par$batch==FALSE){
      if(!is.null(object$vpc.all)){
        #idx=1
        #object<-res21
        if(!is.null(object$org.par)) cycle<-object$org.par$cycle
        else cycle<-FALSE
        
        vpc<-object$vpc.all[[idx]] # vc, data frame
        vpv<-object$vpv.all[[idx]] # vpv, vector
        if(cycle==TRUE) var<-object$evp.all[[idx]] # Vc with SE !!cycle
        else var<-AFEchidna::Var(object$esr.all[[idx]])
        #class(object$esr.all[[idx]])
        
        trt<-unlist(object$esr.all[[idx]]$Traits)
        if(Rres==FALSE) cat('Analysis trait is:',trt,'\n\n')  ##??
        
        # vc
        vc<-vpc[,1]
        
        # vvc.mat
        if(is.matrix(vpv)) vpv.mat<-vpv ###!!!
        else{
          vpv.mat<-diag(1,nrow=length(vc))
          vpv.mat[upper.tri(vpv.mat,diag=T)]<-vpv
          a<-vpv.mat
          a[lower.tri(a)]<-t(a)[lower.tri(a)]
          vpv<-a
        }
        
      } else{
        
        if(is.null(object$vpv.mat)) 
          stop('Please use "W components" under "vpredict".')
        
        vpc<-object$vpc
        vc <- vpc[,'Vc']
        
        vpv.mat<-object$vpv.mat
      }
      
      result<-AFEchidna::pin33(vpc=vpc,vc=vc,vpv.mat=vpv.mat,formula=formula)
      
      if(is.null(object$vpc.all)){
        vc1<-AFEchidna::Var(object)[,c('Term','Sigma','SE')]
      } else vc1<-var[,c('Term','Sigma','SE')] ##??
      
      vc0<-vc1
      names(vc1)<-c('Term','Estimate','SE')
      vc0$vcS<-paste0('V',1:nrow(vc0))
      
      if(Rres==FALSE){
        cat("variance components are as following:\n")
        print.data.frame(vc0)
        cat('\npin formula: ');print(formula)
        cat('\n\n')
      }
      if(all==TRUE) result<-rbind(vc1,result)
      
    }
    
    if(object$org.par$batch==TRUE){
      #object<-res21
      
        data<-object[[1]]
        
        #nn1<-nn2<-data.frame()
        NTrait<-names(data)
        NTrait<-gsub(' ','-',NTrait)
        NTrait<-sub('-$','',NTrait)
        
        res<-lapply(1:length(data), function(x){
          #x=1
          xx<-AFEchidna::pin33(fm=data[[x]],formula=formula)
          vc<-AFEchidna::evp.res(data[[x]]$evp)
          #names(vc)<-NTrait
          #vc1<-do.call(rbind,vc)
          vc1<-vc[,-1]
          names(vc1)[2]<-'Estimate'
          if(all==TRUE) xx<-rbind(vc1,xx)
          xx
        })
        names(res)<-NTrait
        result<-do.call(rbind,res)
        
      }
    result<-AFEchidna::fdata(result,faS=list(0,c(2:3)),FtN=TRUE)
    #result1<-result
    if(signif==TRUE) result$Siglevel<-AFEchidna::sig.level(result$Estimate,result$SE)
    
    AFEchidna::output(result,Rres=Rres,signif=signif,digit=digit) 
  }
}


#' @export
pin33 <- function(fm=NULL,vpc=NULL,vc=NULL,vpv.mat=NULL,
                  formula) {
  if(!is.null(fm)){
    vpc<-fm$vpc
    vc <- vpc[,'Vc']
    vpv.mat<-fm$vpv.mat
  }
 
  #formula<-h2~V2*4/(V2+V5)
  dd<-gsub('V','x',formula) # 'x'
  dd[2]<-as.character(formula[[2]])
  formula1<-stats::as.formula(paste(dd[2],dd[3],sep=' ~ '))
  
  transform<-formula1
  
  pframe <- as.list(vc)
  names(pframe) <- paste("x", seq(1, length(pframe)), sep = "") # 'x'
  tvalue<-eval(stats::deriv(transform[[length(transform)]], names(pframe)),pframe)
  tname <- if(length(transform)==3){transform[[2]]}else ""
  
  invAI <- vpv.mat
  se <- msm::deltamethod(transform,vc,invAI)
  
  result<-data.frame(row.names =tname, Estimate=tvalue, SE=se)
  result$Term<-rownames(result)
  result<-result[,c('Term','Estimate','SE')]
  rownames(result)<-NULL
  
  return(result)
}


#' @usage 
#'       pin22(...,mulp,signif=FALSE,digit=3)
#' @rdname  pin
#' @export  pin22
#' @export
pin22<-function(...,mulp,signif=FALSE,digit=3){
  args <- list(...)
  objs <- all.vars(match.call())
  
  for(i in 1:length(args)){
    cat("\npin results for: ", objs[i],'\n')

    AFEchidna::pin(args[[i]],mulp=mulp,signif=signif,digit=digit)

  }
  
}


# ------------------------------------
# sig.level functions
#' @export
sig.level <- function(tvalue,se,...){
  n <- length(se)
  siglevel <- 1:n
  for(i in 1:n){
    sig.se <- c(se[i]*1.645,se[i]*2.326,se[i]*3.090)
    # 1.450?

    if(abs(tvalue[i])>sig.se[3]) {siglevel[i] <- "***"}
    else if(abs(tvalue[i])>sig.se[2]) {siglevel[i] <- "**"}
    else if(abs(tvalue[i])>sig.se[1]) {siglevel[i] <- "*"}
    else {siglevel[i] <- "Not signif"}
  }
  return(siglevel)
}


#### output format
#' @export
output<-function(res,signif=TRUE,Rres=FALSE,digit=3){#
  if(Rres==FALSE){
    cat("pin results are as following:\n\n")
    if(signif==TRUE){
      if(is.data.frame(res)){
        #res$Signif<-AFfR::sig.level(res$Esmate,res$SE)
        print.data.frame(format(res, digits=digit,nsmall=digit))
      } else print(res)
      cat("---------------")
      cat("\nSig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1\n")
    }else{
      if(is.data.frame(res)) print.data.frame(format(res, digits=digit,nsmall=digit))
      else print(res)
    }
    cat("\n")
  }  
  
  if(Rres==TRUE){
    if(is.data.frame(res)) return(format(res, digits=digit,nsmall=digit))
    else return(res)
  }
  
}

#### 
#' @title Model comparison for Echidna.
#' 
#' @description 
#' \code{model.comp} This function would compare models with different random structure 
#' under the same fixed factors.
#'  
#' @usage model.comp(...,LRT=FALSE,boundary = TRUE) 
#' 
#' @param ...	 A list with more than two Echidna-R results, such as "m1,m2,m3,m4".
#' @param LRT	 Value TRUE for Likelihood ratio test (LRT), default (FALSE) for no LRT. 
#' @param boundary	 boundary If \code{TRUE} (the default) hypothesized parameter
#' values being tested lie on the boundary of the parameter space.
#' @author Yuanzhen Lin <yzhlinscau@@163.com>
#' @references
#' Yuanzhen Lin. R & ASReml-R Statistics. China Forestry Publishing House. 2016 
#' 
#' @examples
#' 
#' \dontrun{
#' library(AFEchidna)
#' 
#' es0.path="D:/teaching/"
#' 
#'               
#' #####   model comparison    #####
#' model.comp(m1,m2)
#' model.comp(m1,m2,LRT=TRUE)
#' }
#' 

#' @export
model.comp <- function(...,LRT=FALSE, boundary = TRUE)
{
  #args <- list(m1,m2)
  
    args <- list(...)
    n <- length(args)
    
    ## fixed models must be the same
    fixed.labels <- lapply(args, function(x) {
      
      if(class(x$org.par$fixed)=='formula'){
        tt<-as.character(x$org.par$fixed)[3]
        tt<-unlist(strsplit(tt,'\\+'))
        tt<-gsub('\\s+','',tt)
      } else tt<-x$org.par$fixed
      tt
    })
    
    batch.state <- sapply(args, function(x) {
      x$org.par$batch
    })
    
    if(stats::var(batch.state)==0&batch.state[1]==FALSE){
      
      ## check fixed
      trms <- length(unique(unlist(lapply(fixed.labels, function(x)x))))
      nt <- unlist(lapply(fixed.labels, function(x)length(x)))
      if(!all(nt == trms))
        stop("fixed models differ")
      
      
      pchisq.mixture <- function(x, n=2) {
        df <- 0:n
        mixprobs <- stats::dbinom(df, size=n, prob=0.5)
        p <- c()
        for (i in 1:length(x)){
          p[i] <- sum(mixprobs*pchisq(x[i], df))
        }
        p
      }
      
      objs <- all.vars(match.call())
      LL <- sapply(args, function(x)x$LogLikelihood)
      np <- sapply(args, function(x){
        x$ICparameter.Count})
      
      AIC<-sapply(args, function(x)x$AkaikeIC)
      BIC<-sapply(args, function(x)x$BayesianIC)
      
      tab0 <- matrix(nrow=n,ncol=5)
      tab0[,1] <- np
      tab0[,2] <- LL
      tab0[,3] <- AIC
      tab0[,4] <- BIC
      #tab0[,5] 
      tab0 <- data.frame(tab0)
      
      #tab0.ncol=5
      b1<- which.min(AIC)
      tab0[,5] <- ""
      tab0[b1,5] <- "better"
      b2 <- which.min(BIC)
      tab0[,6] <- ""
      tab0[b2,6] <- "better"
      
      attr(tab0, "names") <- c("parNO","LogL","AIC",'BIC','AIC.State','BIC.State')
      attr(tab0, "row.names") <- objs
      class(tab0) <- c("data.frame") #"anova",
      tab0<-tab0[base::order(tab0$parNO),]
      cat('\nModel comparison results as following:\n\n')
      print.data.frame(tab0)
      
      if(LRT==TRUE){
        cat("\n\n=====================================")
        cat("\nLikelihood ratio test (LRT) results:")
        cat('\nnote:left model before "/" is full model,right is reduced.\n\n')
        
        ## Increasing order of number or parameters
        idx <- base::order(np)
        D <- df <- prob <- numeric(n-1)
        models <- character(n-1)
        ## do successive pairs
        for(i in seq(1,n-1)) {
          ii <- idx[i]   # reduced
          jj <- idx[i+1] # full
          D[i] <- -2*(LL[ii]-LL[jj])
          df[i] <- np[jj]-np[ii]
          # constraint <- guzpfx(args[[jj]]$vparameters.con)
          # df[i] <- df[i] - sum(as.numeric(constraint == "F"))
          # df[i] <- df[i] - sum(as.numeric(constraint == "C"))
          # df[i] <- df[i] - sum(as.numeric(constraint == "S"))
          # constraint <- guzpfx(args[[ii]]$vparameters.con)
          # df[i] <- df[i] + sum(as.numeric(constraint == "F"))
          # df[i] <- df[i] + sum(as.numeric(constraint == "C"))
          # df[i] <- df[i] + sum(as.numeric(constraint == "S"))
          
          if(df[i] == 0)
            warning(ii," vs ",jj,": same number of parameters.")
          
          if(boundary) {
            if(df[i] == 1)
              prob[i] <- 0.5 * (1 - stats::pchisq(D[i],df[i]))
            else
              prob[i] <- 1 - pchisq.mixture(D[i],df[i])
          }
          else
            prob[i] <- 1 - stats::pchisq(D[i],df[i])
          
          models[i] <- paste(objs[jj],"/",objs[ii],sep="")
        }
        tab <- matrix(nrow=n-1,ncol=3)
        tab[,1] <- df
        tab[,2] <- D
        tab[,3] <- prob
        tab <- data.frame(tab)
        heading <- "Likelihood ratio test(s) assuming nested random models.\n"
        if(boundary)
          heading <- paste(heading,
                           "(See Self & Liang, 1987)\n",sep="")
        attr(tab, "heading") <- heading
        attr(tab, "names") <- c("df","LR-statistic","Pr(Chisq)")
        attr(tab, "row.names") <- models
        class(tab) <- c("anova","data.frame")
        return(tab)
        
      }

  } else{
    cat('not work for batch analysis!') ##!!!
  }
  #cat('\n')
}

#' @export
model.comp11<-function(mulM,LRT=FALSE,rdDF=FALSE,obL=23){
  # obL, each echdina results length
  
  #if(!require(dplyr,warn=F,quiet=T)) stop('Need package: dplyr.\n')
  #require(poorman,warn=F,quiet=T)
  
  # if(object$org.par$batch==TRUE){
  #   cat('not work for batch analysis!') ##!!!
  # }else{
    cat("Attension:\n")
    cat("Fixed factors should be the same!\n\n\n")
    
    Mnames <- vector()#<-NULL
    LogL <- parNo <- DF <- AIC<-BIC<-vector()
    Npm <- 0
    Nml<-mulM
    
    Nmls <- ceiling(length(Nml)/obL) 
    #LogL=Pm=Nedf=vector()
    for(i in 1:Nmls){
      DF[i] <- Nml[[6+(i-1)*obL]]  
      parNo[i] <- Nml[[9+(i-1)*obL]]
      LogL[i] <- Nml[[5+(i-1)*obL]] 
      AIC[i]  <- Nml[[7+(i-1)*obL]]
      BIC[i]  <- Nml[[8+(i-1)*obL]]
      
      Mnames <- all.vars(match.call())
    }
    #print(Mnames)
    
    df<- data.frame(DF=DF,parNO=parNo,LogL=LogL,AIC=AIC,BIC=BIC)
    
    #ifelse(mulM==TRUE,df$Model<-paste("m",1:Nmls,sep="")#, df$Model<-Mnames) 
    df$Model<-Mnames[1:Nmls]
    df<-df[,c(6,1:5)]
    df.ncol<-ncol(df)
    
    b1<- which.min(AIC)
    df$AIC.State <- ""
    df[b1,(df.ncol+1)] <- "better"
    b2 <- which.min(BIC)
    df$BIC.State <- ""
    df[b2,(df.ncol+2)] <- "better"
    df <- df[base::order(df$parNO),] #dplyr::arrange(df,parNo)
    
    #invisible(df)
    
    print.data.frame(df)
    cat("-----------------------------\n")
    cat("Lower AIC and BIC is better model.\n\n")
    
    A <- utils::combn(1:Nmls,2)
    B <- Nmls*(Nmls-1)/2
    
    if(LRT==TRUE){
      #cat("\n")
      cat("\nAttension: Please check every result's length is the same (default,23);\n")
      cat("if the length is less, put the object at the end of mulM.\n")
      cat("In the present, just allow one object's length less.")
      cat("\n=====================================")
      cat("\nLikelihood ratio test (LRT) results:\n\n")
      
      for(i in 1:B){
        if(B>1) df1 <- df[A[,i],1:4] else df1 <- df[1:2,1:4]
        df1 <- df1[base::order(df1$parNO),]
        #df1 <- dplyr::arrange(df1,parNo)
        DlogL <- df1$LogL[2]-df1$LogL[1]
        Ddf <- df1$parNO[2]-df1$parNO[1]
        
        pv<-ifelse(rdDF==TRUE,round(1-stats::pchisq(2*DlogL,Ddf-0.5),3),
                   round(1-stats::pchisq(2*DlogL,Ddf),3))
        
        # pv<-ifelse(rdDF==TRUE,.5*round(1-pchisq(2*DlogL,Ddf-0.5),3),
        #                .5*round(1-pchisq(2*DlogL,Ddf),3))
        
        df1$Ddf<-c(NA,Ddf)
        df1$DlogL<-c(NA,DlogL)
        df1$pv <- c(NA,pv)
        names(df1)[7] <- "Pr(Chisq)"
        
        #cat("\nModel compared between ",df1$Model[1],"--",df1$Model[2],":\n")
        # heading <-paste0("\nModel compared between ",df1$Model[1],"--",df1$Model[2],":\n")
        attr(df1, "row.names") <- df1$Model
        class(df1) <- c("anova","data.frame")
        # return(df1)
        cat('\n')
        print(df1[,-1])
        #cat("---------------")
        #cat("\nSig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1\n\n")
      }
      cat("=====================================")
      if(rdDF==TRUE){
        cat("\nAttension: Ddf=Ddf-0.5. \n")
        cat("When for corr model, against +/-1. \n\n")
      }else {
        cat("\nAttension: Ddf did not minus 0.5. \n")
        cat("When for corr model, against 0. \n\n")
      }
    }
  # }

}

## calculate accuracy for random effects
# ped is pedigree
# ran.eff is random effects
# Var is related variance

#' @usage raneff.acc(object,ran.eff,Var, ped=NULL)
#' @rdname  AF.Echidna
#' @export
raneff.acc <- function(object, ran.eff, Var, ped=NULL) {
  
  if(!require(dplyr,warn.conflicts=FALSE,quietly=TRUE)) 
    stop('Need package: dplyr.\n')
  if(!require(pedigree,warn.conflicts=FALSE,quietly=TRUE)) 
    stop('Need package: pedigree.\n')
  
  require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
  
  #ran.eff<-ran.eff %>% filter(Term !='mv')
  
  if(object$org.par$batch==FALSE){
    if(is.null(ped)) ran.eff$Fi<- 0
    if(!is.null(ped)) {
      ped0<-ped
      names(ped)[1]<-'ID'
      ped1<-pedigree::add.Inds(ped)
      ped1$Fi<-pedigree::calcInbreeding(ped1)
      ped2<-ped1 %>% dplyr::filter(ID %in% ped$ID)
      ped0$Fi<-ped2$Fi
      
      ran.eff$Fi<- ped2 %>% 
        dplyr::filter(ID %in% ran.eff$Level) %>% 
        dplyr::select(Fi)
    }
    
    cat('Var is:',Var,'; Fi is inbreeding coefficient.')
    cat('\naccurancy formula: sqrt(1-SE^2/((1+Fi)*Var))\n')
    ran.eff$SE<-as.numeric(ran.eff$SE)
    ran.eff<- transform(ran.eff,accuracy=sqrt(1-SE^2/((1+Fi)*Var)))
    
    names(ran.eff)[ncol(ran.eff)]<-'accuracy'
    
    cat('The first six results:\n\n')
    print.data.frame(utils::head(ran.eff))
    cat('...\n\n')
    
    return(ran.eff)
  }
  

}

#=====================================
#' @export
predict <- function(object,...){
  UseMethod("predict",object)
}
#' @method  predict esR
#' @export  predict.esR
#' @rdname  AF.Echidna
#' @export

predict.esR <- function(object) {
  #object<-res11
  
  if(object$org.par$batch==FALSE){
    if(object$org.par$cycle==FALSE){
      
      preds<-AFEchidna::predict0(object)
       
    }else{
      
      trt<-names(object$esr.all)
      #preds<-list();length(preds)<-length(trt)
      preds <- vector("list", length(trt))
      
      pred<-object$pred
      
      #pred<-preds[[1]]
      predt<-unlist(strsplit(pred,'TITLE:'))
      #predt[1]
      predt<-predt[-1]# predt[predt!="  "]
      #length(predt)
      
      for(i in 1:length(trt)){
        #i=1
        readr::write_file(predt[i],'tmpf') #predt[i]
        #preds[[i]]<-mult.pred('tmpf')
        txt<-base::readLines('tmpf')
        
        a0<-which(grepl('(Predicted values\\s+.*)', txt))
        heads0<-txt[a0]
        
        aa<-which(grepl('(Prediction\\s+.*)', txt))
        bb<-which(grepl('(Average\\s+.*)', txt))
        #txt[bb]
        
        ased0<-as.numeric(unlist(regmatches(txt[bb],
                                            gregexpr("[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?",
                                                     txt[bb], perl=TRUE))))
        
        names(ased0)<-paste0('ased',1:length(bb))
        ased0<-as.list(ased0)
        
        #preds0<-list();length(preds0)<-length(aa)
        preds0 <- vector("list", length(aa))
        for(j in 1:length(aa)){
          preds0[[j]]<-utils::read.table(text=txt[aa[j]:(bb[j]-1)],header=T)
          preds0[[j]]<-list(pred=preds0[[j]],ased=ased0[[j]])
        }
        names(preds0)<-paste0('pred',1:length(aa))
        preds[[i]]<-list(heads=heads0,pred=preds0)
      }
      
      trt<-gsub(' ','-',trt)
      trt<-sub('-$','',trt)
      names(preds)<-trt

    }
  }
  
  if(object$org.par$batch==TRUE){
    trt<-unlist(lapply(object$res.all,function(x) x$Traits))
    #preds<-list();length(preds)<-length(trt)
    preds <- vector("list", length(trt))
    preds0<-preds
    
    #preds$pred<-lapply(object$res.all,function(x) x$pred)
    for(i in 1:length(trt)) {
      preds0[[i]]$pred<-object$res.all[[i]]$pred
      class(preds0[[i]])<-'esR'
    }
    #names(preds)<-trt
    
    preds<-lapply(preds0,predict0) ##???
    
    trt<-gsub(' ','-',trt)
    trt<-sub('-$','',trt)
    names(preds)<-trt

  }
  return(preds)
}

#' @export
predict0<-function(object){
  if(is.null(object$pred)) stop('Please use predict statements.')
  
  txt<-object$pred
  
  head0<-which(grepl('(Predicted values\\s+.*)', txt))
  heads<-txt[head0]
  
  aa<-which(grepl('(Prediction\\s+.*)', txt))
  bb<-which(grepl('(Average\\s+.*)', txt))
  
  preds<-list()
  for(j in 1:length(aa)){
    preds[[j]]<-utils::read.table(text=txt[aa[j]:(bb[j]-1)],header=T)
  }
  names(preds)<-paste0('pred',1:length(aa))
  
  ased<-as.numeric(unlist(regmatches(txt[bb],
                                     gregexpr("[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?",
                                              txt[bb], perl=TRUE))))
  
  names(ased)<-paste0('ased',1:length(bb))
  ased<-as.list(ased)
  
  tt<-list(heads=heads,pred=preds,ased=ased)
  
  preds<-tt
  return(preds)
}


#' @export
coef <- function(object){
  UseMethod("coef",object)
}
#' @rdname  AF.Echidna
#' @method  coef esR
#' @export  coef.esR
#' @export

coef.esR <- function(object){
  #object<-HT
  
  if(class(object$org.par$random)=='formula'){
    random <- as.character(object$org.par$random)[2]
    if(grepl('\\*',random)) tt <- Var(object)$Term[Var(object)$Term!="Residual"]
      else tt <- random
    tt<-unlist(strsplit(tt,'\\+'))
    tt<-gsub('\\s+','',tt)
  } else tt <- object$org.par$random
  ranf <- tt
  
  if(object$org.par$batch==FALSE){
    if(is.null(object$esr.all)){
      if(is.null(object$coef)) 
        stop('Please use qualifier !SLN.')
      
      trt <- object$Traits # AFEchidna::
      sol <- AFEchidna::coeff11(object$coef,ranf=ranf)
      
    }else{ # !cycle
      #object<-res13
      
      trt <- names(object$esr.all)#unlist(lapply(object$esr.all,function(x) x$Traits))
      #sol0 <- list();length(sol0) <- length(trt)
      sol0 <- vector("list", length(trt))
      
      coeft<-object$coef
      aa<-which(grepl('(Model_Term,\\s+.*)', coeft))
      
      for(i in 1:(length(aa)-1)){
        sol0[[i]]<-utils::read.table(text=coeft[aa[i]:(aa[i+1]-1)],
                              header=TRUE,sep=',')
      }
      sol0[[length(aa)]]<-read.table(text=coeft[aa[length(aa)]:length(coeft)],
                                    header=TRUE,sep=',')

      sol<-lapply(sol0,function(x) AFEchidna::coeff11(x,ranf=ranf))
      
      trt<-gsub(' ','-',trt)
      trt<-sub('-$','',trt)
      names(sol)<-trt
    }
  }
  
  if(object$org.par$batch==TRUE){
    
    trt<-unlist(lapply(object$res.all,function(x) x$Traits))
    trtN<-length(trt)
    #sol<-list();length(sol)<-trtN
    sol <- vector("list", trtN)
    sol0<-sol
    
    for(i in 1:trtN) {
      sol0[[i]]$coef<-object$res.all[[i]]$coef
      class(sol0[[i]])<-'esR'
    }

    sol<-lapply(1:trtN,function(x) AFEchidna::coeff11(sol0[[x]]$coef,ranf=ranf))
    
    trt<-gsub(' ','-',trt)
    trt<-sub('-$','',trt)
    names(sol)<-trt

  }
 
  return(sol)
}


#' @export
coeff11 <- function(object,ranf) {
  
  require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
  
  coeft<-object#$coef
  
  names(coeft)[1:4]<-c('Term','Level',"Effect",'SE')
  coeft$Term<-gsub('\\.',':',coeft$Term)
  coeft$Level<-gsub('\\t','',coeft$Level)
  coeft <-coeft %>% dplyr::filter(Term!='mv')
  
  ranf<-gsub('\\.',':',ranf)
  Term<-unique(coeft$Term)
  fixf<-base::setdiff(Term,ranf)
  
  fixeff<-AFEchidna::filterD1(coeft,Term %in% fixf) #AFEchidna::
  raneff<-AFEchidna::filterD1(coeft,Term %in% ranf)
  #head(raneff)
  
  attr(fixeff,'terms')<-'fixed'
  attr(raneff,'terms')<-'random'
  
  effects<-list(fixed=fixeff,random=raneff)
  
  return(effects)
}

#coef(res)$fixed
#coef(res)$random


#' @export
wald <- function(object,...){
  UseMethod("wald",object)
}
#' @method  wald esR
#' @export  wald.esR
#' @rdname  AF.Echidna
#' @export

wald.esR<-function(object) {
  
  #object<-mm
  
  if(object$org.par$batch==FALSE){
    if(is.null(object$esr.all)){
      cat("\nWald tests for fixed effects.\n")
      if(!is.null(object$waldT)){
        cat(object$waldT)
       }
    }else{
      nn<-length(object$esr.all)
      trt<-unlist(lapply(object$esr.all,function(x) x$Traits))
      for(i in 1:nn){
        cat("\nWald tests of fixed effects for trait: ",trt[i],'\n')
        cat(object$esr.all[[i]]$waldT,'\n')#names(object$esr.all[[1]])
        cat('\n')
      }
    }
  }
  
  if(object$org.par$batch==TRUE){
    trt<-unlist(lapply(object$res.all,function(x) x$Traits))
    
    trt<-gsub(' ','-',trt)
    trt<-sub('-$','',trt)
    
    for(i in 1:length(trt)) {
      
      cat('\nWald tests of fixed effects for trait: ',trt[i],'\n')
      cat(object$res.all[[i]]$waldT,'\n')
      cat('\n')
    }
  }

}

##

#' @export
waldT <- function(object,...){
  UseMethod("waldT",object)
}
#' @method  waldT esR
#' @export  waldT.esR
#' @rdname  AF.Echidna
#' @export
waldT.esR<-function (object, term=NULL,ncol=NULL) 
{  
  if(is.null(term)){
    waldT2<-strsplit(object$waldT,'P-inc \r\n')[[1]][2]
    term<-unlist(regmatches(waldT2, 
                      gregexpr("[A-Za-z]+\\.?\\w+", 
                               waldT2, perl = TRUE)))
  }
  term<-gsub('\\.','\\:',term) # change . to : in R
  n <- length(term)
  
  waldTv <- as.numeric(unlist(regmatches(object$waldT, 
                                         gregexpr("[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?", 
                                                  object$waldT, perl = TRUE))))
  if (object$org.par$mulT) { # TRUE
    if(is.null(ncol)) ncol <- 2
    wv <- matrix(waldTv, ncol = ncol, byrow=TRUE) # ncol=2
    wv <- data.frame(wv)
    
    denDF <- rep(object$Residual.DF,n)
    numDF <- as.numeric(as.character(wv[,1]))
    Fv <- as.numeric(as.character(wv[,2]))
  } else { #F: single trait
    if(is.null(ncol)) ncol <- 5
    wv <- matrix(waldTv, ncol = ncol, byrow=TRUE) # ncol=5
    wv <- data.frame(wv) 
    
    numDF <- as.numeric(as.character(wv[,1]))
    denDF <- as.numeric(as.character(wv[,2]))
    Fv <- as.numeric(as.character(wv[,3]))
  }
  
  pv <- NULL
  for (i in seq(1, n)) #pv[i] <- 1 - pchisq(Fv[i], numDF[i])
    pv[i] <- pf(Fv[i], numDF[i],denDF[i],lower.tail=FALSE)
  tab <- matrix(nrow = n, ncol = 3)
  tab[, 1] <- numDF[1:n]
  tab[, 2] <- Fv[1:n]
  tab[, 3] <- pv[1:n]
  tab <- data.frame(tab)
  tab[, 2] <- tab[, 1] * tab[, 2] 
  heading <- "Wald tests for fixed effects:\n"
  attr(tab, "heading") <- heading
  attr(tab, "names") <- c("DF", "Wald-F", 
                          "Pr(Chisq)")
  attr(tab, "row.names") <- term
  class(tab) <- c("anova", "data.frame")
  return(tab)
}


#' @export
IC <- function(object,...){
  UseMethod("IC",object)
}
#' @method  IC esR
#' @export  IC.esR
#' @rdname  AF.Echidna
#' @export

IC.esR <-function(object){
  #object<-res13
  
  if(object$org.par$batch==FALSE){
    if(is.null(object$esr.all)){
      df<-get.IC(object)
      #print(df)
      
    }else{ #!cycle
      
      #nn<-length(object$esr.all)
      trt<-unlist(lapply(object$esr.all,function(x) x$Traits))
      #preds<-list();length(preds)<-length(trt)
      preds<-vector("list", length(trt))
      
      for(i in 1:length(trt)) {
        preds[[i]]<-object$esr.all[[i]]
      }
      
      IC0<-lapply(preds,get.IC)   
      df<-do.call(rbind,IC0)
      
      trt<-gsub(' ','-',trt)
      trt<-sub('-$','',trt)
      rownames(df)<-trt
      
      #print(df)
    }
  }
  
  if(object$org.par$batch==TRUE){
    trt<-unlist(lapply(object$res.all,function(x) x$Traits))
    #preds<-list();length(preds)<-length(trt)
    preds<-vector("list", length(trt))
    
    for(i in 1:length(trt)) {
      preds[[i]]<-object$res.all[[i]]
    }
    #names(preds)<-trt
    
    IC0<-lapply(preds,get.IC)   
    df<-do.call(rbind,IC0)
    
    trt<-gsub(' ','-',trt)
    trt<-sub('-$','',trt)
    row.names(df)<-trt
    #print(df)
  }
  return(df)
}

get.IC<-function(object){
  data.frame(DF=object$Residual.DF,parNO=object$ICparameter.Count,
             LogL=object$LogLikelihood,
             AIC=object$AkaikeIC,BIC=object$BayesianIC)
}


#' @export
trace <- function(object,...){
  UseMethod("trace",object)
}
#' @method  trace esR
#' @export  trace.esR
#' @rdname  AF.Echidna
#' @export
## iteration trace of runs
trace.esR<-function(object){
  #object<-spm.esr
  
  # if(object$org.par$family==" !poisson !log !disp 1")
  # cat(object$Iterations00)
  
  if(object$org.par$batch==FALSE){
    if(is.null(object$esr.all)){
      #cat('\nEchinda version: ',object$Version,'\n')
      cat('\n',object$StartTime,'\n')
      print.data.frame(object$Iterations)
      cat(object$FinishAt,'\n\n')
    }else{
      nn<-length(object$esr.all)
      trt<-unlist(lapply(object$esr.all,function(x) x$Traits))
      #cat('\nEchinda version: ',object$esr.all[[1]]$Version,'\n')
      
      for(i in 1:nn){
        cat('\n\nIteration procedure for trait: ',trt[i],'\n')
        cat('\n',object$esr.all[[i]]$StartTime,'\n')
        print.data.frame(object$esr.all[[i]]$Iterations)
        cat(object$esr.all[[i]]$FinishAt,'\n')
      }
    }
  }
  
  if(object$org.par$batch==TRUE){
    trt<-unlist(lapply(object$res.all,function(x) x$Traits))
    
    trt<-gsub(' ','-',trt)
    trt<-sub('-$','',trt)
    
    #preds$pred<-lapply(object$res.all,function(x) x$pred)
    for(i in 1:length(trt)) {
      #preds[[i]]<-object$res.all[[i]]$Iterations
      cat('\n\nIteration procedure for trait: ',trt[i],'\n')
      cat('\n',object$res.all[[i]]$StartTime,'\n')
      print.data.frame(object$res.all[[i]]$Iterations)
      cat(object$res.all[[i]]$FinishAt,'\n')
    }
  }

}



#' @export
Var <- function(object,...){
  UseMethod("Var",object)
}
#' @method  Var esR
#' @export  Var.esR 
#' @rdname  Var
#' @export
Var.esR<-function(object){
  
  batch<-object$org.par$batch
  cycle<-object$org.par$cycle
  
  mulT<-object$org.par$mulT
  mulN<-object$org.par$mulN
  mulT<-object$org.par$mulT
  trace<-object$org.par$trace
  
  if(batch==FALSE){
    if(cycle==FALSE){
      #object<-r1.res
      
      if(is.null(object$evp)) 
        stop('Please use "W components" under "vpredict".')
      tt<-AFEchidna::evp.res(object$evp)
      #cat(mm$evp)
      for(i in 3:4) tt[,i]<-as.numeric(as.character(tt[,i]))
      tt$Z.ratio<-tt$Sigma/tt$SE
      
      tt$Term<-gsub('\\.',':',tt$Term)
      
      ## for MET
      var1<-tt[,-1]
      
      #return(var1)
    }
    if(cycle==TRUE){ # !cycle
      #object<-res13b
      
      nn<-length(object$esr.all)
      trt<-unlist(lapply(object$esr.all,function(x) x$Traits))
      
      Converge<-sapply(object$esr.all,function(x) x$Converge)
      maxit<-sapply(object$esr.all,function(x) x$maxit)
      
      vpc0<-object$vpc.all[[1]]
      terms<-vpc0$Term # res,fam,...
      
      vpc<-do.call(rbind,object$evp.all)#[[idx]]
      
      varv<-vpc$Sigma
      varm<-matrix(varv,ncol=length(terms),byrow=TRUE)
      vardf<-as.data.frame(varm)
      names(vardf)<-paste0('V',1:ncol(vardf))
      
      sev<-vpc$SE
      sem<-matrix(sev,ncol=length(terms),byrow=TRUE)
      sedf<-as.data.frame(sem)
      names(sedf)<-paste0(names(vardf),'.se')
      
      var1<-cbind(vardf,sedf)
      var1$Converge<-Converge
      var1$maxit<-maxit
      
      trt<-gsub(' ','-',trt)
      trt<-sub('-$','',trt)
      row.names(var1)<-trt
      
      vcms<-terms
      vcms<-paste(paste0('V',1:length(vcms)),vcms,sep='-',collapse='; ')
      
      attr(var1, "heading") <- paste0('\n',vcms,'\n',
                                      'Converge: 1 means True; 0 means FALSE.\n')
      attr(var1, "row.names") <- trt
      class(var1) <- c("anova","data.frame")
      #return(var1)
    }
  }
  
  if(batch==TRUE){
    #object<-r2g#res21b
    
    trt<-names(object$res.all)
    trt<-gsub(' ','-',trt)
    trt<-sub('-$','',trt)
    trt<-sub('^-','',trt)
    trtN<-length(trt)
    
    Converge<-sapply(object$res.all,function(x) x$Converge)
    maxit<-sapply(object$res.all,function(x) x$maxit)
    
    #var<-list();length(var)<-trtN
    var<-vector("list", trtN)
    
    var<-lapply(1:trtN, function(x){
      var0<-AFEchidna::evp.res(object$res.all[[x]]$evp)
      #terms<-var0[,'Term']
      mat.var<-matrix(c(var0[,'Sigma'],var0[,'SE']),nrow=1)
      mat.var<-as.data.frame(mat.var)
      names(mat.var)<-c(paste0('V',1:nrow(var0)),paste0('V',1:nrow(var0),'.se'))
      mat.var
    })
    
    var0<-AFEchidna::evp.res(object$res.all[[1]]$evp)
    terms<-var0[,'Term']
    
    var1<-do.call(rbind,var)
    row.names(var1)<-trt
    var1$Converge<-Converge
    var1$maxit<-maxit
    
    if(mulT==TRUE) terms<-AFEchidna::vcms.fun(terms,mulN)
    
    terms<-paste(paste0('V',1:length(terms)),terms,sep='-',collapse='; ')
    
    head0<- paste0('\n',terms,'\n',
                   'Converge: 1 means True; 0 means FALSE.\n')
    attr(var1, "heading") <-head0 
    attr(var1, "row.names") <- trt
    class(var1) <- c("anova","data.frame")
  }
  return(var1)
}


#' @export
Var.AR<-function(object,delN=4){
  #object<-m3b
  patt<-"[-+]?[0-9]*\\.[0-9]+([eE][-+]?[0-9]+)?"
  numv<-as.numeric(unlist(regmatches(object$keyres,
                                     gregexpr(patt, object$keyres, perl=TRUE))))
  
  numm<-matrix(numv[-1:-delN],ncol=3,byrow=T)#
  
  varcomp<-as.data.frame(numm)
  names(varcomp)<-c('Gamma','Sigma','Zratio')
  varcomp$SE<-varcomp$Sigma/varcomp$Zratio
  
  varcomp <- round(varcomp,3)
  varcomp <- varcomp[,c('Sigma','SE','Zratio')] 
  varcomp$Siglevel<-sig.level(varcomp$Sigma,varcomp$SE)
  
  Term <- Var(object)$Term
  ssr <- which(grepl("^id(.*)", Term))
  if(length(ssr)!=0) Term <- Term[-ssr] 
  ss <- which(grepl("(ar1(.*))", Term))
  Term[ss] <- paste0(Term[ss],'.cor')
  varcomp$Term <- Term[c(2:length(Term),1)]
  varcomp <- varcomp[,c(length(varcomp),1:(length(varcomp)-1))]
  
  return(varcomp)
}


#' @title Summarize an esR object
#' 
#' @description
#' A \code{summary} method for objects inheriting from class
#' \code{esR}.
#' 
#' @param object
#'
#' An \code{esR} object. 
#' 
#' @return
#'
#' A list of class \code{summary.esR} with the following components:
#'
#' \describe{
#'
#' \item{org.res}{Original results from .esr file in Echidna.} 
#'   
#' \item{varcomp}{A dataframe summarising the random variance component.}
#' 
#' \item{IC}{nedf, loglik, Akaike information criterion and Bayesian information criterion.}
#' 
#' \item{coef.fixed}{A dataframe of coefficients and their standard errors
#' for fixed effects.}
#'
#' \item{coef.random}{A dataframe of coefficients and their standard errors
#' for random effects.}
#' 
#' }
#' @export
summary <- function(object,...){
   UseMethod("summary",object)
 }
#' 
#' 
#' @export  summary.esR
# @rdname  AF.Echidna
#' @method  summary esR
#' @aliases summary
#' @export
#'
summary.esR <- function(object){
  
  sum.object <- vector(mode="list")
  if(object$org.par$batch0==FALSE & is.null(object$esr.all)){
      keyres <- object$keyres
      sum.object$org.res <- keyres  # cat(keyres)
      sum.object$varcomp <- AFEchidna::Var(object)
      sum.object$IC <- AFEchidna::IC(object)
      sum.object$coef.fixed  <- AFEchidna::coef(object)$fixed
      sum.object$coef.random <- AFEchidna::coef(object)$random
    }else{  
      if(!is.null(object$esr.all)) xxxx <- object$esr.all ## !cycle
      if(!is.null(object$res.all)) xxxx <- object$res.all ## batch
      nn <- length(xxxx) #length(object$esr.all)
      trt <- unlist(lapply(xxxx,function(x) x$Traits))
     
      for(i in 1:nn){
        sum.object$org.res[[i]] <- xxxx[[i]]$keyres #object$esr.all[[i]]$keyres
        sum.object$coef.fixed[[i]]  <- AFEchidna::coef(object)[[i]]$fixed
        sum.object$coef.random[[i]] <- AFEchidna::coef(object)[[i]]$random
      }
      names(sum.object$coef.fixed) <-names(AFEchidna::coef(object))
      names(sum.object$coef.random)<-names(AFEchidna::coef(object))
      sum.object$IC <- AFEchidna::IC(object)
    }

  class(sum.object) <- "summary.esR"
  return(sum.object)
}

#' @export
converge <- function(object,...){
  UseMethod("converge",object)
}
#' @method  converge esR
#' @export  converge.esR
#' @rdname  AF.Echidna
#' @export
#'
converge.esR<-function(object){
  #object<-res
  
  if(object$org.par$batch==FALSE){
    if(is.null(object$esr.all)){
      cat(object$Converge)
    }else{
      nn<-length(object$esr.all)
      trt<-unlist(lapply(object$esr.all,function(x) x$Traits))
      for(i in 1:nn){ # i=1
        cat('\nConverge state for trait: ',trt[i],' ')
        cat(object$esr.all[[i]]$Converge,'\n')
        cat('\n')
      }
    }
  }
  
  if(object$org.par$batch==TRUE){
    trt<-unlist(lapply(object$res.all,function(x) x$Traits))
    trt<-gsub(' ','-',trt)
    trt<-sub('-$','',trt)

    for(i in 1:length(trt)) {
      cat('\nConverge state for trait: ',trt[i],'\n')
      cat('\n',object$res.all[[i]]$Converge,'\n')
      cat('\n')
    }
  }
}

#### original site codes connection with results

#' @title MET analysis
#' 
#' @description 
#' \code{met.plot} plots MET data.\code{met.corr} calculates var/cov/corr from echidna MET factor analytic 
#' results to further research the relation of trial sites.
#' \code{met.biplot} This function biplots MET factor analytic results from echidna 
#' to find the relation of trial sites and the best variety suitable to trial sites.  
#' 
#' @param data	 MET data.
#' @param plot.title	 MET plot title.
#' @param object	 echidna factor analytic results for MET, such as mm.
#  @param aimS  Specify the aim location parts of echidna object to count corr matrix.
#' @param rotate  Rotate the factor loadings, FALSE(default).
#' @param siteV	  trial site values, a vector or data frame.
#' @param kn	 Site cluster group numbers, 3(default).
#' @param horiz  output cluster site result format, horiz(default).
#' @param biplot	 output biplots, FALSE(default).
#  @param plg	 Adding labels before site, "S"(default). 
#  @param dmethod	 The distance measured method for site cluster, "manhattan"(default), more details see amap::hcluster.
#' @param dSco.u	 Least score of Variety breeding value. 
#' @param dLam.u	 Least distance from center.
#'
#' @rdname  met
#' @export met.plot
#' @author Yuanzhen Lin <yzhlinscau@@163.com>
#' @references
#' Yuanzhen Lin. R & ASReml-R Statistics. China Forestry Publishing House. 2016 
#' @examples 
#' \dontrun{
#' library(AFEchidna)
#' 
#' path="C:/Users/echi/exam" #home
#' setwd(path)
#' 
#' MET<-read.csv('MET.csv')
#' names(MET)
#' 
#' # example 1
#' # variable order: yield,genotype,site,row,col
#' MET2<-MET[,c(9,1,2,4:5)] 
#' str(MET2)
#' met.plot(MET2)
#' 
#' # example 2
#' MET3<-MET[,c(9,1,2,4:7)] # add variable order on MET2: Rep, Block
#' str(MET3)
#' met.plot(MET3,"My met trials")
#' 
#' ## running met analysis with FA model
#' mm<-echidna(es0.file="MET.es0",trait='yield',fixed='Loc',
#'       random='Genotype.xfa2(Loc)',
#'       residual='sat(Loc).units', #sat(Loc).ar1(Col).ar1(Row)
#'       #predict=c('Genotype'),
#'       vpredict=c('V Vmat Genotype.xfa1(Loc)','R cor 20:40'),
#'       qualifier='!maxit 50 !SLN',
#'       foldN='mm',
#'       met=T)
#'       
#' Var(mm)
#' 
#' siteV<-unique(MET['Loc']) # should be a data.frame or vector
#' 
#' met.corr(mm,siteV=siteV)
#' 
#' met.biplot(mm,siteV=siteV)
#' met.biplot(mm,siteV=siteV,biplot=T)
#' met.biplot(mm,siteV=siteV,biplot=T,dSco=1.0,dLam=0.8)
#' 
#' res2<-met.vmat(mm,siteV=siteV,VmatN='Vmat',corN='cor')
#' res2$res
#' res2$var
#' 
#' }
#' 

#' @usage met.plot(data, plot.title = NULL,...) 
#' @export
met.plot <-function(data,plot.title=NULL,...){

    if(!require(desplot))
      stop('Need package: desplot.\n')
  
    require(desplot,warn.conflicts=FALSE,quietly=TRUE)
    
    if(is.null(plot.title)) plot.title <- "MET data plot"
    
    dat <- data
    levels(dat[,3]) <- paste("S",1:nlevels(dat[,3]),sep="")
    names(dat)[1:5] <- c("yield","genotype","site","row","col")
    for(i in 4:5) dat[,i] <- as.numeric(dat[,i])
    
    #windows(10,8)
    # desplot(yield~ col*row|site, dat, main=plot.title)
    if(length(dat)==5){  
      desplot::desplot(yield~ col*row|site, dat, main=plot.title)
    }else{    
      names(dat)[6:7] <- c("Rep","Blk")  
      desplot::desplot(yield ~ col*row|site, dat, main=plot.title,
                       out1=Rep, out2=Blk,strip.cex=1.5,
                       out1.gpar=list(col="blue", lwd=4),
                       out2.gpar=list(col="red", lwd=1, lty=1),
                       par.settings = list(layout.heights=list(strip=2)))
    } 
}


#' @usage met.corr(object,siteV,kN=NULL,horiz=TRUE,rotate=FALSE) 
#' @rdname  met
#' @export
met.corr <- function(object,siteV,kN=NULL,horiz=TRUE,rotate=FALSE) { 
  
  if(!require(dplyr))
    stop('Need package: dplyr.\n')
  if(!require(amap))
    stop('Need package: amap.\n')
  
  require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
  require(amap,warn.conflicts=FALSE,quietly=TRUE)
  
  #mm<-object
  
  if(is.null(kN)) kN <- 3
  
  if(is.data.frame(siteV)) {
    siteV<-factor(siteV[,1])
    siteN <- n<- nlevels(siteV) 
  } else siteN <- n<- length(siteV)
  
  
  varcomp <- AFEchidna::Var(object)[c('Term','Sigma')]
  varcomp <- varcomp %>% dplyr::filter(grepl('xfa', Term)) %>% dplyr::select(Sigma)
  
  vect1 <- varcomp[1:n,]
  w.var <- diag(vect1)
  vect2 <- varcomp[(n+1):nrow(varcomp),]
  t.var <- matrix(vect2,nrow=siteN)
  #crossprod(t.var)
  
  if(rotate==TRUE){
    cat('Using rotating factor loadings.\n\n')
    
    L <- t.var
    svd.L <- svd(L)
    L.star <- L %*% svd.L$v
    t.var <- L.star
  }
  
  wt.var <- t.var%*%t(t.var)+w.var
  df <- wt.var
  
  df2<-stats::cov2cor(wt.var)
  rownames(df2)<-colnames(df2)<-paste0('S-',siteV)
  
  df[upper.tri(df)]<-df2[upper.tri(df2)]
  df<-as.data.frame(df)
  
  rownames(df)<-colnames(df)<-paste0('S-',siteV)
  
  cat("\nVar\\Cov\\corr matrix:\n")
  print.data.frame(round(df,3))
  cat("---------------")
  cat("\ndiag is Variance, lower is covariance, upper is correlation.\n")
  
  
  df2<-stats::na.omit(df2)
  chcluster <- amap::hclusterpar(df2, method="manhattan")
  graphics::plot(chcluster,main='Cluster of different sites', 
              hang=-1,ylab='',xlab='')  
  stats::rect.hclust(chcluster, k=kN)
  cat("\nSite cluster results:\n")
  id <- stats::cutree(chcluster,k=kN)
  cl.res<- data.frame(cl.No=id,row.names=rownames(df))
  if(horiz) print(t(cl.res)) else print.data.frame(cl.res)
  cat('\n')
  
  if(.Platform$OS.type == "windows") invisible(df)
}

#' @usage met.biplot(object,siteV,biplot=FALSE, dSco.u=NULL,dLam.u=NULL)
#' @rdname  met
#' @export
met.biplot <- function(object,siteV,biplot=FALSE,dSco.u=NULL,dLam.u=NULL) {
  #object<-mm
  require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
  #is.data.frame(Var(mm))
  
  component<-AFEchidna::Var(object)[c('Term','Sigma')] 
  arr<-component %>% dplyr::filter(grepl('xfa', Term)) %>% dplyr::select(Sigma)
  
  if(is.data.frame(siteV)) {
    siteV<-factor(siteV[,1])
    siteN <- n<- nlevels(siteV) 
  } else siteN <- length(siteV)
  
  sname<-paste("S-",siteV,sep="")
  
  Xfam<-matrix(arr[,1],nrow=siteN)#,(1+faN))
  faN<-ncol(Xfam)-1
  
  fa.name<-paste("FA",1:faN,sep="")
  dimnames(Xfam)<-list(sname,c("Psi",fa.name))
  #windows(8,8)
  graphics::pairs(Xfam,main="Fig.1-- pairs of Psi with FAs")
  
  ss<-svd(Xfam[,-1])
  Lam<-Xfam[,-1] %*% ss$v
  colnames(Lam)<-c(paste("FA",1:faN,sep="")) # c("Fa1","Fa2")
  Gvar<-Lam %*% t(Lam)+diag(Xfam[,1])
  cLam<-diag(1/sqrt(diag(Gvar))) %*% Lam  ##??
  varp<-round(mean(diag(Lam %*% t(Lam))/diag(Gvar))*100,2) # %variance explained
  cat("\nFA number is:",faN,",\t%Var.explained is: ",varp,".\n")
  
  
  # get each factor for each site
  dt.Lam<-as.data.frame(Lam)
  row.names(dt.Lam)<-sname
  for(i in 1:faN) dt.Lam[(faN+i)]<-dt.Lam[i]^2
  colnames(dt.Lam)[(faN+1):(2*faN)]<-c(paste("sq.FA",1:faN,sep=""))
  dt.Lam$Ps.Var<-Xfam[,1]
  dt.Lam$T.Var<-rowSums(dt.Lam[(faN+1):(2*faN+1)],na.rm=T)
  nnn<-ncol(dt.Lam)
  for(i in 1:faN) dt.Lam[(nnn+i)]<-100*dt.Lam[(faN+i)]/dt.Lam$T.Var
  names(dt.Lam)[(nnn+1):(nnn+faN)]<-c(paste("per.FA",1:faN,sep=""))
  #nnn1<-ncol(dt.Lam)
  dt.Lam[(nnn+faN+1)]<-rowSums(dt.Lam[(nnn+1):(nnn+faN)],na.rm=T)
  names(dt.Lam)[(nnn+faN+1)]<-"tper.FA"
  mmm<-nrow(dt.Lam)
  dt.Lam[mmm+1,]<-colMeans(dt.Lam)
  row.names(dt.Lam)<-c(row.names(dt.Lam)[1:mmm],'Mean')
  
  #print(format(dt.Lam[(nnn+1):(nnn+faN+1)],digits=3,nsmall=3))
  print.data.frame(format(dt.Lam,digits=3,nsmall=3))
  
  if(biplot){ ### 
    if(!require(tidyr,warn.conflicts=FALSE,quietly=TRUE)) 
      stop('Need package: tidyr.\n')
    require(tidyr,warn.conflicts=FALSE,quietly=TRUE)
    
    if(faN==1){cat("\nAttension: biplot worked when more than 2 FAs!\n\n")}

    if(faN>1){

      #object<-m7
      bv<-AFEchidna::coef(object)$random # here would be complexed!!
      #head(bv)

      Xfasln<-bv %>% dplyr::filter(grepl('\\.F', Level)) #%>% select(Effect)
      #head(Xfasln1)
      Xfasln$Level<-gsub('\t','',Xfasln$Level)
      
      # here careful!!
      Xfasln1<-tidyr::separate(data = Xfasln, col = Level,
                               into = c("Genotype", "Loc"), sep = "\\.")

      #VarietyN=70
      VarietyN<-nrow(Xfasln1)/faN
      scores<-matrix(Xfasln1$Effect,nrow=VarietyN)
      dimnames(scores)<-list(unique(Xfasln1$Genotype),paste("Fa",1:faN,sep=""))

      acb<-utils::combn(1:faN,2)
      bl<-faN*(faN-1)/2

      mLam<-rep(1/siteN,siteN) %*% Lam # get loading means
      sLam<-Lam-rep(mLam,rep(siteN,faN)) # center loadings
      dLam<-sqrt((sLam*sLam) %*% rep(1,faN)) # distance from center
      dSco<-sqrt((scores*scores) %*% rep(1,faN))

      dSco.a<-0.65*max(dSco,rm=TRUE)
      dLam.a<-max(dLam,rm=TRUE)

      if(is.null(dSco.u)) dSco.u<-round(dSco.a,1) # 2
      if(is.null(dLam.u)) dLam.u<-round(dLam.a,1) #0.1

      if(faN>2){
        for(i in 1:bl){
          #windows(18,8)
          #par(mfrow=c(1,2))
          stats::biplot(scores[,acb[,i]],Lam[,acb[,i]],cex=0.75,
                 main=paste("Fig 2-",i, " biplot with all variety",sep=""))
          graphics::abline(h=0,lty=3)
          graphics::abline(v=0,lty=3)
          if(nrow(Lam[dLam>dLam.u,1:2])>1){
            stats::biplot(scores[dSco>dSco.u,acb[,i]],Lam[dLam>dLam.u,acb[,i]],cex=0.75,
                   main=paste("Fig 3-",i, " biplot when dSco>",dSco.u,sep="")) # dSco>2
            graphics::abline(h=0,lty=3)
            graphics::abline(v=0,lty=3)
          }else {
            cat("\ndSco is:\n")
            print(tail(sort(dSco),6))
            cat("\ndLam is:\n")
            print(round(dLam,3))
            cat('\nthe second figure failure, we should set up dLam.\n')
          }
        }
      }else {
        #windows(18,8)
        #par(mfrow=c(1,2))
        stats::biplot(scores[,1:2],Lam[,1:2],cex=0.75,
               main="Fig 2 biplot with all variety")
        graphics::abline(h=0,lty=3)
        graphics::abline(v=0,lty=3)

        if(nrow(Lam[dLam>dLam.u,1:2])>1){
          stats::biplot(scores[dSco>dSco.u,1:2],Lam[dLam>dLam.u,1:2],cex=0.75,
                 main=paste("Fig 3 biplot when dSco>",dSco.u,' and dLam>',dLam.u,sep="")) # dSco>2
          graphics::abline(h=0,lty=3)
          graphics::abline(v=0,lty=3)
        }else {
          cat("\ndSco is:\n")
          print(utils::tail(sort(dSco),6))
          cat("\ndLam is:\n")
          print(round(dLam,3))
          cat('\nthe second figure failure, we should set up dLam.\n')
        }

      }

      dscores<-data.frame(scores[dSco>dSco.u,],Scores=dSco[dSco>dSco.u]) #2
      ddLam<-data.frame(Lam[dLam>dLam.u,],distFC=dLam[dLam>dLam.u]) # 0.1

      cat("\nScores.u is:",dSco.u,"\n")
      print.data.frame(round(dscores,3))
      cat("\ndistFC.u is:",dLam.u,"\n")
      print.data.frame(round(ddLam,3))
      cat("\n")
    }
  }
  
}

#' @usage met.vmat(object,siteV,VmatN,corN)
#' @rdname  met
#' @export
met.vmat <- function(object,siteV,VmatN='Vmat',corN='cor') {
  
  if(!require(reshape,warn.conflicts=FALSE,quietly=TRUE)) 
    stop('Need package: reshape.\n')
  if(!require(dplyr,warn.conflicts=FALSE,quietly=TRUE)) 
    stop('Need package: dplyr.\n')
  
  require(reshape,warn.conflicts=FALSE,quietly=TRUE)
  require(dplyr,warn.conflicts=FALSE,quietly=TRUE)
  
  odd <- function(x) x %% 2 != 0
  
  if(is.data.frame(siteV)) {
    siteV<-factor(siteV[,1])
    siteN <- n<- nlevels(siteV) 
  } else siteN <- length(siteV)
  
  sname<-paste("S-",siteV,sep="")
  
  #path1<-paste(path,foldN,sep='/')
  #evpf<-paste(path1,'temp.evp',sep='/')
  res.evp<-object$evp0 #base::readLines(evpf)
  
  patt1<-paste0(VmatN,' 2')
  vmat.str<-res.evp[grep(patt1,res.evp)]
  #print(vmat.str)
  vmat.v<-as.numeric(unlist(regmatches(vmat.str,
                                       gregexpr("[-+]?[0-9]*\\.[0-9]+([eE][-+]?[0-9]+)?",
                                                vmat.str, perl=TRUE))))
  vn<-1:length(vmat.v)
  
  v.mat<-vmat(siteN=siteN, vec=vmat.v[odd(vn)])
  se.mat<-vmat(siteN=siteN, vec=vmat.v[!odd(vn)])
  
  
  rownames(v.mat)<-colnames(v.mat)<-sname
  rownames(se.mat)<-colnames(se.mat)<-sname
  
  patt2<-paste0(corN,' 2')
  cor.str<-res.evp[grep(patt2,res.evp)]
  #print(cor.str)
  cor.v<-as.numeric(unlist(regmatches(cor.str,
                                      gregexpr("[-+]?[0-9]*\\.[0-9]+([eE][-+]?[0-9]+)?",
                                               cor.str, perl=TRUE))))
  
  
  corn<-1:length(cor.v)
  cor.mat<-AFEchidna::vmat(siteN=siteN, vec=cor.v[odd(corn)],cor=TRUE)
  corse.mat<-AFEchidna::vmat(siteN=siteN, vec=cor.v[!odd(corn)],cor=TRUE)
  
  #cat('Var\\cov\\corr matrix:\n')
  #corn<-1:length(cor.v)
  v.mat1<-v.mat
  v.mat1[upper.tri(v.mat1,diag=F)]<-cor.v[odd(corn)]
  print(v.mat1)
  
  #cat('\nThe se for Var\\cov\\corr matrix:\n')
  se.mat1<-se.mat
  se.mat1[upper.tri(se.mat1,diag=F)]<-cor.v[!odd(corn)]
  print(se.mat1)
  
  #cat('\nThe sig.level for Var\\cov\\corr matrix:\n')
  sig.l1<-AFEchidna::siglevel(vmat.v[odd(vn)],vmat.v[!odd(vn)])
  sig.mat<-AFEchidna::vmat(siteN=siteN, vec=sig.l1)
  
  sig.l2<-AFEchidna::siglevel(cor.v[odd(corn)],cor.v[!odd(corn)])
  sig.mat[upper.tri(sig.mat,diag=F)]<-sig.l2
  
  rownames(sig.mat)<-colnames(sig.mat)<-sname
  print(noquote(sig.mat))
  
  cat('\nThe se for corr matrix:\n')
  corse.mat[upper.tri(corse.mat,diag=F)]<-cor.v[odd(corn)]
  rownames(corse.mat)<-colnames(corse.mat)<-sname
  print(corse.mat)
  
  cat('\nThe sig.level for corr matrix:\n')
  corsig.mat<-cor.mat
  corsig.mat[upper.tri(corsig.mat,diag=F)]<-sig.l2
  rownames(corsig.mat)<-colnames(corsig.mat)<-sname
  print(noquote(t(corsig.mat)))
  
  cat("=================\n")
  cat("upper is corr and lower is error (or sig.level) for corr matrix.\n")
  cat("Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1\n\n")
  
  df<-list(vmat=v.mat1,se.vmat=se.mat1,sig.vmat=sig.mat,
           cor.se.mat=corse.mat,cor.sig.mat=corsig.mat)
  
  vmat0<-v.mat1
  rownames(v.mat1)<-colnames(v.mat1)<-1:nrow(v.mat1)
  rownames(se.mat1)<-colnames(se.mat1)<-1:nrow(se.mat1)
  rownames(sig.mat)<-colnames(sig.mat)<-1:nrow(sig.mat)
  
  rr0<-reshape::melt(vmat0)
  rr1<-reshape::melt(v.mat1)
  rr2<-reshape::melt(se.mat1)
  rr3<-reshape::melt(sig.mat)
  
  rr4<-merge(rr1,rr2,by=c('X1','X2'),sort=F)
  rr5<-merge(rr4,rr3,by=c('X1','X2'),sort=F)
  
  rr5$Type<-ifelse(rr5$X1<rr5$X2,'Cor','Var')
  
  rr5$X1<-rr0$X1
  rr5$X2<-rr0$X2
  
  rr5<-dplyr::arrange(rr5,dplyr::desc(Type),X1,X2)
  names(rr5)[1:5]<-c('site1','site2','estimate','se','Sig.level')
  
  df<-list(res=df,var=rr5)
  
  return(df)
  
}

#' @export
vmat <- function(siteN, vec,cor=FALSE) {
  vmat<-diag(siteN) # siteN
  if(cor==TRUE) vmat[upper.tri(vmat,diag=F)]<-vec
  else vmat[upper.tri(vmat,diag=T)]<-vec
  a<-vmat
  a[lower.tri(a)]<-t(a)[lower.tri(a)]
  
  return(a)
}

## for batch
#' @export
vm2r <- function(object,corN=NULL) {
  if(!is.null(object$res.all)){
    tt<-names(object$res.all)
    ss<-lapply(1:length(tt),function(x) AFEchidna::vm2r0(object,idx=x,corN=corN))
    names(ss)<-tt
  }else ss<-AFEchidna::vm2r0(object,corN=corN)
  
  return(ss)
}

#' @export
vm2r0 <- function(object,idx=1,corN=NULL) {
  
  #object<-res22
  if(!is.null(object$res.all))
    evp.str<-strsplit(object$res.all[[idx]]$evp,'\r\n')[[1]] # !batch
  
  if(!is.null(object$evp))
    evp.str<-strsplit(object$evp,'\r\n')[[1]] 
  
  if(!is.null(corN)){
    #corN<-'cor'
    patt2<-paste0(corN,' 2')
    cor.str<-evp.str[grep(patt2,evp.str)]
    num.patt<-"[-+]?[0-9]*\\.[0-9]+([eE][-+]?[0-9]+)?"
    cor.v<-as.numeric(unlist(regmatches(cor.str,
                                        gregexpr(num.patt,
                                                 cor.str, perl=TRUE))))
    
    rdf<-as.data.frame.matrix(matrix(cor.v,ncol=2,byrow=T)) # corr
    rdf$Type<-'cor'
    
    cor.str1<-cor.str[grep('cor 2\\s+.*=',cor.str)]
    
    
    cor.str1<-unlist(regmatches(cor.str1,
                                gregexpr('(?<=cor 2\\s)\\s+\\d.*=\\s+us',
                                         cor.str1, perl=TRUE)))
    cor.str1<-gsub('\\s+','.t',cor.str1)
    cor.str1<-gsub('.t=.tus','',cor.str1)
    row.names(rdf)<-paste0('r',cor.str1)
    
    
    var.str<-evp.str[evp.str!=cor.str]
  } else var.str<-evp.str
  
  var.v<-as.numeric(unlist(regmatches(var.str,
                                      gregexpr("[-+]?[0-9]*\\.[0-9]+([eE][-+]?[0-9]+)?",
                                               var.str, perl=TRUE))))
  vardf<-as.data.frame.matrix(matrix(var.v,ncol=2,byrow=TRUE)) # var
  vardf$Type<-'var'
  #print(cor.str)
  
  if(!is.null(object$res.all)){
    terms<-attr(object$varcomp,'heading')
    terms<-strsplit(terms,'\n')[[1]][2]
    terms<-gsub('V\\d+-','',terms)
    terms<-unlist(strsplit(terms,'; '))
    row.names(vardf)<-terms
  }

  if(!is.null(corN)) resdf<-rbind(vardf,rdf)
  else resdf<-vardf
  
  names(resdf)[1:2]<-c('Estimate','SE')
  
  return(resdf)

}


#======================================================
#' @title Generate genomic relationship matrix.
#' 
#' @description 
#' \code{GenomicRel} This function generates 5 genomic relationship matrixs.
#' 
#' \tabular{ll}{
#' \strong{option} \tab \strong{Description} \cr
#' 1    \tab observed allele frequencies (GOF, VanRaden, 2008). \cr
#' 2    \tab weighted markers by recipricals of expected variance (GD, Forni et al., 2011). \cr
#' 3    \tab allele frequencies fixed at 0.5 (G05, Forni et al., 2011). \cr
#' 4    \tab allele frequencies fixed at mean for each locus (GMF, Forni et al., 2011).\cr
#' 5    \tab regression of MM' on A sort (Greg, VanRaden, 2008).
#' }  
#'  
#' @usage GenomicRel(marker,option,ped=NULL,
#'                    Infv=1000,Gm=NULL,G.adj=FALSE,Gres=TRUE) 
#' 
#' @param marker	 markers data.
#' @param option	 option (1~5) for different G matrixs.
#' @param ped    ped data.
#' @param Infv   A value for Inf in G matrix. 
#' @param Gm     G matrix from marker.
#' @param G.adj  Adjust G matrix with A matrix from ped data. 
#' @param Gres   return G matrix directly(True, default) or in data frame(FALSE).  
#' @author Isik Fikret
#' @references
#' Isik Fikret. Genetic data analysis for plant and animal breeding. 2017 
#' @examples 
#' \dontrun{
#' library(AFEchidna)
#' 
#' read.example(package = "AFEchidna", setpath = TRUE)
#' Markers<-read.file(file="sim_markers.txt",sep=' ' )
#' ped<-read.table( "sim_pedigree.txt", sep=' ')
#' 
#' 
#'  GOF1=GenomicRel( Markers,1)
#'  GD1=GenomicRel(  Markers,2)
#'  G051=GenomicRel( Markers,3)
#'  GMF1=GenomicRel( Markers,4)
#'  
#'  # the same result but adjust G matrix with ped data:
#'  GOF2=GenomicRel(  Markers,1, ped,G.adj=T)
#'  GD2=GenomicRel(   Markers,2, ped,G.adj=T)
#'  G052=GenomicRel(  Markers,3, ped,G.adj=T)
#'  GMF2=GenomicRel(  Markers,4, ped,G.adj=T)
#'  Greg=GenomicRel(  Markers,5, ped,G.adj=T)
#' }
#' 
#' @export .GenomicRel
#' @export
GenomicRel <- function(marker,option=NULL,ped=NULL,Infv=10000,
                       Gm=NULL,G.adj=FALSE,Gres=TRUE){
  return(.GenomicRel(marker,option,ped,Infv,Gm,G.adj,Gres))
}
## functions to generate G matrix from SNP markers.

.GenomicRel <- function(marker, option=NULL,ped=NULL,Infv=10000,
                       Gm=NULL,G.adj=FALSE, Gres=TRUE){ # ,invG=FALSE
  
  
  #data
  #	column 1 - id
  #	column 2:m markers 
  #	markers must be 0, 1, 2 for homozygous, heterozygous and other homozygous
  
  #pedigree
  #	only needed for option 5
  #	ID, Sire, Dam in columns 1,2 and 3, respectively
  #	must be sorted with animal in first column before they are in 2nd or 3rd column
  
  #NO HEADERS on input tables
  #Returns G matrix in the following form: col1 = row, col2 = col, col3=Genomic relationship, col4=pedigree relationship
  
  if(!require(GeneticsPed)){stop('Need package: GeneticsPed.\n')}
  
  require(GeneticsPed,warn.conflicts=FALSE,quietly=TRUE)
  
  cat('Generating G matrix is under going.\n')
  
  options(warn=-1)
  data<-marker
  
  if(!is.null(ped)){
    names(ped)<-c("id","father","mother")
    ped<-GeneticsPed::as.Pedigree(ped)
    fullA<-GeneticsPed::relationshipAdditive(ped)
    rowName<-as.numeric(rownames(fullA))
    A<-fullA[match(data[,1],rowName),match(data[,1],rowName)]
  }
  
  #library(MASS)
  #library(GeneticsPed)
  if(!is.null(Gm)){
    G<-Gm
  }else{
    
    M1<-data
    
    if(option!=5){
      # M matrix
      if(option %in% c(1:2)) M<- M1[,2:ncol(M1)]-1 
        else M<- M1[,2:ncol(M1)]
      
      # p-value
      if(option==3){
        p1<-array(0.5,ncol(M))
        p<-p1
      }else{ # option=1,2,4
        
        if(option==4) p1<-round((apply(M,2,mean)),3) 
        else p1<-round((apply(M,2,sum)+nrow(M))/(nrow(M)*2),3)
        
        p<-2*(p1-.5)
      }
      
      P <- matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M))
      Z <- as.matrix(M-P)
      
      if(option==2){
        D<-1/(ncol(M)*(2*p1*(1-p1)))
        D[!is.finite(D)]<-Infv
        
        G <- Z%*%(D*t(Z))
      }else{ # option=1,2,4
        b<-1-p1
        c<-p1*b
        d<-2*(sum(c))
        
        ZZt <- Z %*% t(Z)
        G <- (ZZt/d)
      }
    }else{ # option=5
      M<- M1[,2:ncol(M1)]-1
      
      M<-as.matrix(M)
      MtM<-M%*%t(M)
      
      n<-(nrow(A))^2
      sumA<-sum(A)
      SqA<-A*A
      sumSqA<-sum(SqA)
      
      rhs1<-sum(MtM)
      rhs2<-sum(MtM*A)
      
      lhs<-cbind(rbind(n,sumA),rbind(sumA,sumSqA))
      rhs<-rbind(rhs1,rhs2)
      g<-solve(lhs,rhs)
      
      one<-matrix(1,nrow=nrow(data))
      
      G<-(MtM-(g[1,1]*(one%*%t(one))))/g[2,1]
    }
  }    
  
  
  if(G.adj==TRUE){
    GW<-round(0.95*G+(1-0.95)*A,3)
    #if(invG==TRUE) A=solve(A)
  }else GW<-G
  
  col1<-col2<-col3<-col4<-NA
  # col2<-NA
  # col3<-NA
  # col4<-NA
  for (i in 1:nrow(GW)){
    for (j in 1:i){
      
      col1<-cbind(col1,i)
      col2<-cbind(col2,j)
      col3<-cbind(col3,GW[i,j])
      if(!is.null(ped))	col4<-cbind(col4,round(A[i,j],3))
    }	
    
  }
  Gmat<-cbind(t(col1),t(col2),t(col3))
  row.names(Gmat)<-c(0:(nrow(Gmat)-1))
  
  if(!is.null(ped))	{
    Amat<-cbind(t(col1),t(col2),t(col4))
    row.names(Amat)<-c(0:(nrow(Amat)-1))
  }
  
  #write.table(Ginv[-1,],"Ginv.giv",sep=" ", row.names=F, col.names=F,quote = F)
  if(!is.null(ped)){
    df<-data.frame("col"=Gmat[-1,1],"row"=Gmat[-1,2],"G"=Gmat[-1,3],"A"=Amat[-1,3])
  } else df<-data.frame("col"=Gmat[-1,1],"row"=Gmat[-1,2],"G"=Gmat[-1,3])
  
  if(Gres==TRUE) return(GW) else return(df)
  
}

#======================================================
#' @title Generate H-inverse matrix for SS-GBLUP.
#' 
#' @description 
#' \code{AGH.inv} This function calculate genomic relationship matrix(G),
#' full additative matrix(A) and blended relationship matrix(H) from
#' genotyped marker, genotyped pedigree and ungenotyped pedigree.
#'  
#' @usage AGH.inv(gmarker,option=1, ugped, gped) 
#' 
#' @param option	 option (1~5) for different G matrixs.
#' @param ugped	 ungenotyped pedigree, or total pedigree.
#' @param gped	 genotyped pedigree.
#' @param gmarker	 genotyped marker,column 1 should be sample ID.
#' @param asrV	 asreml version, 3(default) or 4.
#' @details 
#'   This function would return a list containing 3 elements. The types of
#' option (1~5) as following:
#' 
#' \tabular{ll}{
#' \strong{option} \tab \strong{Description} \cr
#' 1    \tab observed allele frequencies (GOF, VanRaden, 2008). \cr
#' 2    \tab weighted markers by recipricals of expected variance (GD, Forni et al., 2011). \cr
#' 3    \tab allele frequencies fixed at 0.5 (G05, Forni et al., 2011). \cr
#' 4    \tab allele frequencies fixed at mean for each locus (GMF, Forni et al., 2011).\cr
#' 5    \tab regression of MM' on A sort (Greg, VanRaden, 2008).
#' } 
#' @return 
#' \describe{ 
#' \item{Ainv}{inverse of full additative matrix(A).}
#' \item{Ginv}{inverse of genomic relationship matrix(G).}	
#' \item{Hinv}{inverse of blended relationship matrix(H).}	
#' }
#' 
#' @author Yuanzhen Lin <yzhlinscau@@163.com>
#' @references
#' Yuanzhen Lin. R & ASReml-R Statistics. China Forestry Publishing House. 2016    
# AFfR website:https://github.com/yzhlinscau/AFfR
#' @examples 
#' \dontrun{
#' library(AFEchidna)
#' 
#' data("ugped")
#' data("gped")
#' data("gmarker")
#' 
#' # get A-matrix, G-matrix and H-matrix
#' AGH1<-AGH.inv(option=1,ugped,gped,gmarker)
#' 
#'   
#' data(MET)
#' MET$yield<-0.01*MET$yield
#' levels(MET$Genotype)<-gped$ID
#' MET1<-filterD1(MET, Loc %in% c(3))
#' 
#' ## for ASReml-R V3.0 
#' library(asreml)
#' 
#' # base model
#' sm1.asr<-asreml(yield~Rep, random=~ Genotype+units, 
#'                 rcov=~ ar1(Col):ar1(Row), 
#'                 data=MET1, maxiter=50)
#' 
#' Var(sm1.asr)
#'
#' # A-BLUP
#' Ainv <- AGH1$Ainv
#' sm2.asr<-update(sm1.asr, random=~ ped(Genotype)+units, 
#'                 ginverse=list(Genotype=Ainv))
#' 
#' Var(sm2.asr)
#' 
#' # G-BLUP
#' Ginv <- AGH1$Ginv
#' sm3.asr<-update(sm1.asr, random=~ ped(Genotype)+units, 
#'                 ginverse=list(Genotype=Ginv))
#' 
#' Var(sm3.asr)
#' 
#' # H-BLUP
#' Hinv <- AGH1$Hinv
#' sm4.asr<-update(sm1.asr, random=~ ped(Genotype)+units, 
#'                 ginverse=list(Genotype=Hinv))
#' 
#' Var(sm4.asr)
#' 
#' 
#' ## for ASReml-R V4 
#' library(asreml)
#' 
#' 
#' sm1.asr<-asreml(yield~Rep, random=~ Genotype+units, 
#'                 residual=~ ar1(Col):ar1(Row), 
#'                 data=MET1, maxiter=50)
#' 
#' Var(sm1.asr)
#'
#' # A-BLUP
#' Ainv <- AGH1$Ainv
#' sm2.asr<-update(sm1.asr, 
#'              random=~ vm(Genotype,Ainv)+units)
#' 
#' Var(sm2.asr)
#' 
#' # G-BLUP
#' Ginv <- AGH1$Ginv
#' sm3.asr<-update(sm1.asr, 
#'              random=~ vm(Genotype,Ginv)+units)
#' 
#' Var(sm3.asr)
#' 
#' # H-BLUP
#' Hinv <- AGH1$Hinv
#' sm4.asr<-update(sm1.asr, 
#'              random=~ vm(Genotype,Hinv)+units)
#' 
#' Var(sm4.asr)
#' 
#' ########## if any other genotyped id without ped
#' ########## we can put their parent code to 0 or NA
#' ########## to make pedigree, then use H-matrix.
#' ## gmarker2 without pedigree
#' #
#' ## make their pedigree
#' # gid2<-gmarker2[,1]
#' # gped2<-data.frame(ID=gid2,Female=0,Male=0)
#' #
#' ## combine genotyped id's pedigree
#' # gped1<-rbind(gped,gped2)
#' #
#' ## combine all genotyped marker data
#' # gmarker1<-rbind(gmarker,gmarker2)
#' #
#' # AGH1a<-AGH.inv(option=1,tped1,gped1,gmarker1)
#' #
#' }
#' @export .AGH.inv
#' @export

AGH.inv <- function(option=1,ugped,gped,gmarker,asrV=3,tidn=NULL,gidn=NULL){
  
  return(.AGH.inv(option,ugped,gped,gmarker,asrV,tidn,gidn))
}

.AGH.inv <- function(option=1,ugped,gped,gmarker,asrV=3,
                     tidn=NULL,gidn=NULL){
  # tidn: vector of total id number
  # gidn: vector of genotyped id number
  
  # names(ped)<- c("id","father","mother")
  #asrV<-getASRemlVersionLoaded(Rsver=TRUE)  
  
  if(!require(nadiv)){stop('Need package: nadiv.\n')}
  #if(!require(synbreed)){stop('Need package: synbreed.\n')}
  if(!require(GeneticsPed)){stop('Need package: GeneticsPed.\n')}
  
  # genotyped ped
  gped<-gped[!duplicated(gped),]
  gid<-as.character(gped[,1])
  
  # total ped
  tped<-rbind(ugped,gped)
  #row.names(tped)<-tped[,1]
  tped<-tped[!duplicated(tped),]
  
  #tped1<-tped
  tped1<-nadiv::prepPed(tped)
  fullA<-as.matrix(nadiv::makeA(tped1))
  tid1<-rownames(fullA)#<-colnames(fullA)
  
  ugid<-base::setdiff(tid1,gid)
  tid1<-c(ugid,gid)
  
  rowName<-rownames(fullA)
  fullA<-fullA[match(tid1,rowName),match(tid1,rowName)]
  
  gNO<-length(gid)
  ugNO<-length(ugid)
  
  A11<-fullA[1:ugNO,1:ugNO] # ungenotyped A
  A22<-fullA[(1+ugNO):(ugNO+gNO),(1+ugNO):(ugNO+gNO)] # genotyped A
  A12<-fullA[1:ugNO,(1+ugNO):(ugNO+gNO)]
  A21<-t(A12)
  #row.names(A22)
  
  G<-AFEchidna::GenomicRel( gmarker, option, Gres=TRUE)
  #G<-AFEchidna::GenomicRel( gmarker, option,gped,G.adj=T, Gres=TRUE)
  
  H11<-A11 + A12 %*% solve(A22) %*% (G-A22) %*% solve(A22) %*% A21
  H12<-G %*% solve(A22) %*% A21
  H21<-t(H12)
  H22<-G
  
  H1<-rbind(H11,H12)
  H2<-rbind(H21,H22)
  
  H<-cbind(H1,H2)
  
  if(!is.null(tidn)){
    class(fullA)<-c("relationshipMatrix", "matrix")
    rowName<-as.numeric(rownames(fullA))
    fullA<-fullA[match(tidn,rowName),match(tidn,rowName)]
    
    class(H)<-c("relationshipMatrix", "matrix")
    rowName<-as.numeric(rownames(H))
    H<-H[match(tidn,rowName),match(tidn,rowName)]
  }
  if(!is.null(gidn)){
    class(G)<-c("relationshipMatrix", "matrix")
    rowName<-as.numeric(rownames(G))
    G<-G[match(gidn,rowName),match(gidn,rowName)]
  }
  
  ## H-inverse
  
  class(H)<-c("relationshipMatrix", "matrix")
  #summary(eigen(H)$values)
  
  Hinv <- AFEchidna::write.relationshipMatrix(H, 
                                              file =NULL,
                                              sorting=c("ASReml"), 
                                              type=c("inv"), digits=10)
  
  #head(attr(Hinv, "rowNames"),10)
  names(Hinv) <- c("row", "column", "coefficient")
  
  if(asrV=='4'){
    Hinv1<-as.matrix(Hinv)
    attr(Hinv1, "rowNames")<-attr(Hinv, "rowNames")
    class(Hinv1)<-c("matrix","ginv")
    Hinv<-Hinv1
  }
  
  ## A-inverse
  A<-fullA
  class(A)<-c("relationshipMatrix", "matrix")
  #summary(eigen(A)$values)
  
  Ainv <- AFEchidna::write.relationshipMatrix(A, 
                                              file =NULL,
                                              sorting=c("ASReml"), 
                                              type=c("inv"), digits=10)
  
  #head(attr(Ainv, "rowNames"),10)
  names(Ainv) <- c("row", "column", "coefficient")
  
  if(asrV=='4'){
    Ainv1<-as.matrix(Ainv)
    attr(Ainv1, "rowNames")<-attr(Ainv, "rowNames")
    class(Ainv1)<-c("matrix","ginv")
    Ainv<-Ainv1
  }
  
  ## G-inverse
  class(G)<-c("relationshipMatrix", "matrix")
  #summary(eigen(A)$values)
  
  Ginv <- AFEchidna::write.relationshipMatrix(G, 
                                              file =NULL,
                                              sorting=c("ASReml"), 
                                              type=c("inv"), digits=10)
  
  #head(attr(Ainv, "rowNames"),10)
  names(Ginv) <- c("row", "column", "coefficient")
  
  if(asrV=='4'){
    Ginv1<-as.matrix(Ginv)
    attr(Ginv1, "rowNames")<-attr(Ginv, "rowNames")
    class(Ginv1)<-c("matrix","ginv")
    Ginv<-Ginv1
  }
  
  return(list(Ainv=Ainv,Ginv=Ginv,Hinv=Hinv))
}


## functions to output G matrix to giv or grm file for 
## further running in ASReml or Echidna.
#' @export
write.relationshipMatrix <- function(x,file=NULL,sorting=c("WOMBAT","ASReml"),
                                     type=c("ginv","inv","none"),digits=10){
  
  type <- match.arg(type)
  sorting <- match.arg(sorting)
  
  if(sorting=="WOMBAT" & type!="ginv") stop("'type' must be 'ginv' for WOMBAT")
  
  # pass (inverse) relationship matrix
  if(type=="ginv") rMinv <- MASS::ginv(x)
  if(type=="inv")  rMinv <- solve(x)
  if(type=="none") rMinv <- x
  
  rMinv <- round(rMinv,digits)
  
  # add rownames and colnames
  res <- data.frame(Row = rep(1:nrow(rMinv), nrow(rMinv)),
                    Column = rep(1:nrow(rMinv), each = nrow(rMinv)),
                    coeff = as.numeric(rMinv),
                    lower = as.logical(lower.tri(rMinv, diag = TRUE)))
  rm(rMinv)
  
  
  
  # only use lower triangle
  res <- res[res$lower == TRUE, c("Row", "Column", "coeff")]
  
  if (sorting=="ASReml"){    
    res <-  res[order( res$Row, res$Column), ] 
  }
  if (sorting=="WOMBAT"){
    res <- res[, c(2,1,3)]
    res <-  res[order(res$Column, res$Row), ]  
  }
  res <- res[res$coeff != 0, ]
  
  # write to given file
  if (!is.null(file)){
    write.table(res, file, sep = " ", quote = FALSE, 
                col.names = FALSE, row.names = FALSE)
    rm(x)
  } else {
    attr(res, "rowNames") <- rownames(x)
    rm(x)
    return(res)
  }
  
}


#########################
# function to generate specific correlation with another variable
#' @usage rho.par(rho,x0,a=1,b=0)
#' @rdname  AF.base
#' @export
rho.par <- function(rho,x0,a=1,b=0) {
  x0[is.na(x0)]<-0.001
  n     <- length(x0)#20         # length of vector
  rho   <- rho                   # desired correlation = cos(angle)
  theta <- acos(rho)             # corresponding angle
  x1    <- x0#rnorm(n, 1, 1)     # fixed given data
  x2    <- stats::rnorm(n, 2, 0.5)      # new random data
  X     <- cbind(x1, x2)         # matrix
  Xctr  <- scale(X, center=TRUE, scale=FALSE)   # centered columns (mean 0)
  
  Id   <- diag(n)                               # identity matrix
  Q    <- qr.Q(qr(Xctr[, 1, drop=FALSE]))      # QR-decomposition, just matrix Q
  P    <- tcrossprod(Q)          # = Q Q'       # projection onto space defined by x1
  x2o  <- (Id-P) %*% Xctr[, 2]                 # x2ctr made orthogonal to x1ctr
  Xc2  <- cbind(Xctr[, 1], x2o)                # bind to matrix
  Y    <- Xc2 %*% diag(1/sqrt(colSums(Xc2^2)))  # scale columns to length 1
  
  x <- Y[, 2] + (1 / tan(theta)) * Y[, 1]     # final new vector
  
  x<-x*a+b # make values positive 
  
  return(x)
}

##========== weights ==============



##========== family ==============
#' @export
esr_gaussian <- function(link = "identity", dispersion=NA)
{
  ## gaussian family
  
  # link <- as.character(substitute(link))
  # misnames <- c("inverse", "log", "identity", "reciprocal", "1/mu",
  #               "Inverse", "Reciprocal", "Log", "Identity")
  # corresp <- c(1, 2, 3, 1, 1, 1, 1, 2, 3)
  # lmatch <- pmatch(link, misnames, FALSE)
  # if(!lmatch)
  #   stop('Gaussian links are "log", "inverse" or "identity"\n')
  # link <- misnames[corresp[lmatch]]
  # fam <- esr_makeFamily("gaussian",link=link)
  # fam$dispersion <- dispersion
  
  if(!is.na(dispersion)) fam <- paste0(' !disp ',dispersion)
   else fam <- ' '
  
   return(fam)
}


#' @export
esr_binomial <- function(link = "probit", dispersion=1, total=NULL)
{
  if(link=="logit") fam<-paste0(' !binomial !logit !disp ',dispersion)#,'!totals ')
  if(link=="probit") fam<-paste0(' !binomial !probit ')
  if(link=="comploglog") fam<-paste0(' !binomial !comploglog ')
  if(link=="marginal") fam<-paste0(' !binomial !MARG ')
  return(fam)
}

#' @export
esr_poisson <- function(link = "log", dispersion=1.0)
{
  if(link=="log") fam<-paste0(' !poisson !log !disp ',dispersion)
  if(link=="sqrt") fam<-paste0(' !poisson !sqrt !disp ',dispersion)
  #if(link=="identity") fam=paste0(' !poisson !identity !disp ',dispersion)
  return(fam)
}  


#======================================================
# update: 2020-03-13
#' @title Model Qualifiers for Echidna
#'
#' @details
#'   The Qualifiers control various aspects of the model fitting process and reporting of results:
#'
#' \cr
#' \tabular{ll}{
#' \strong{jobqualf } \tab \strong{Description} \cr
#' \code{!WORKSPACE n}       \tab where n is an integer between 1 and 32 gigabytes. \cr
#' \code{!CONTINUE [f] or !FINAL [f] }       \tab instructs Echidna to retrieve variance parameters from an .esv file. \cr
#' \code{!VIEW }    \tab  displays any Winteracter graphics directly to the screen, such as '!view !PNG' \cr
#' \code{!DEBUG}          \tab requests additional debugging output to be written. \cr
#' \code{!LOGFILE}          \tab receive the DEBUG information. \cr
#' \code{!EQN [q]}          \tab sets the !EQN qualifier which selects the equation ordering option, q could be 1~7 or 77. \cr
#' \code{others}           \tab !RENAME,!ARGS,!OUTFOLDER, no needs in R. 
#' }
#' 
#' \cr
#' \tabular{ll}{
#' \strong{qualifier} \tab \strong{Description} \cr
# \code{!MAXIT m}     \tab sets the maximum number of iterations to m,13(default). \cr
#' \code{!EXTRA n}     \tab forces n more iterations after convergence . \cr
#' \code{!SINGLE}      \tab forces Echidna not to use Parallel Processing. \cr
#' \code{!SLOW}      \tab reduces stepsize when updating variance parameters so convergence is slower but
#' possible more reliable. \cr
#' \code{!SKIP i}      \tab indicating the file contains i heading lines to be ignored. \cr
#' \code{!READ f}      \tab specifying that f data fields are to be read. \cr
# \code{!YHT}         \tab requests the residuals and fitted values. \cr
#' \code{!FILTER Variable !SELECT value !EXCLUDE value}         \tab subset the data sets. \cr
#' \code{!MVINCLUDE}   \tab design variables which are missing are assumed to be zero.\cr
#' \code{!MVREMOVE}    \tab data records are ignored if any design variables have missing values.
#' }
#' 
#' \cr
#' \tabular{ll}{
#' \strong{GRM qualifiers } \tab \strong{Description} \cr
#' \code{!SKIP i}       \tab to skip the first i lines of the grm/giv file. \cr
#' \code{!LDET d}       \tab lets you set the logDet value for the supplied giv matrix. \cr
#' \code{!ADD value}    \tab adds value to the diagonal of a GRM matrix before inversion.value can only be 1, 2 or 3. \cr
#' \code{!PSD}          \tab allows a GRM matrix to be positive semidefinite. \cr
#' \code{!NSD}          \tab allows a GRM matrix to be negative semidefinite. \cr
#' \code{!ND}           \tab allows a GRM matrix to be negative definite. 
#' }
#' 
#' \cr Data fields: if VARIABLE is a CLASS variable, the NAME must be followed by a qualifier indicating it
#' is a class variable.
#'
#' \tabular{ll}{
#' \strong{class qualifier} \tab \strong{Description} \cr
#' \code{*}               \tab variable is coded 1:n. \cr
#' \code{!I n}            \tab variable is coded with INTEGER labels, not 1:n. \cr
#' \code{!A n}            \tab variable is coded with ALPHANUMERIC labels. \cr
#' \code{!L <labels>}     \tab data is coded 1:n and <labels> is the list of n class names. \cr
#' \code{!AS <factor> }   \tab when two or more alphanumeric variables have a common list of class names. \cr
#' \code{!LL c }          \tab reset the maximum length to c characters. \cr
#' \code{!PRUNE }         \tab reset the number of levels (n) in a factor to the real maximum levels. 
#' }
#'
#' \cr PEDIGREE FILE: ifile contains 3 or 4 fields being the identity (tag, id, name) of
#'        an individual, its Sire and its Dam..
#'
#' \tabular{ll}{
#' \strong{qualifier} \tab \strong{Description} \cr
#' \code{!SKIP i}        \tab the file contains i heading lines to be ignored. \cr
#' \code{!CSV}           \tab file is strictly COMMA delimited, without .csv suffix. \cr
#' \code{!GROUPS g }     \tab first g lines are genetic groups. . \cr
#' \code{!MGS}           \tab  the third field is a maternal grandsire (not DAM). \cr
#' \code{!SAVE f }       \tab requests the inverse be written as a .giv (f=1, ascii) or .bgiv (f=3, real binary) file. 
#' }
#'
#' 
#' \cr VPREDICT statements have the syntax as: <Key Letter> <Label> <arguments>.
#'    <Key Letter> is one of the letters F, H, R, V, etc.
#'
#' \tabular{ll}{
#' \strong{Key Letter } \tab \strong{Description} \cr
#' \code{F}       \tab specifies a linear function like animal + units. \cr
#' \code{H}       \tab specifies a ratio (heritability) like animal / Total. \cr
#' \code{R}       \tab specifies correlations from a matrix or formula. \cr
#' \code{V}       \tab converts a FA structure to a US structure. \cr 
#' \code{X}       \tab is a multiply function. \cr
#' \code{S}       \tab is a square root function. \cr
#' \code{K}       \tab sets short vectors for constants (e.g. Legendre coefficients) to be used for M. \cr
#' \code{M}       \tab uses the K coefficients to convert US (us(leg(dim,3))) matrix for a specific dim. \cr
#' \code{C}       \tab moves components back e.g. C label I[:II]=J[:JJ] where I < J. \cr
#' \code{D}       \tab discards components eg D from 393. \cr
#' \code{W}       \tab writes components to f.vpc and their variance to f.vpv.
#' }
#' 
#' \cr Variance structures
#'
#' \tabular{ll}{
#' \strong{function} \tab \strong{Description} \cr
#' \code{id(),idv()}      \tab identity or scaled identity. \cr
#' \code{ar1(), ar1v()}   \tab autoregressive correlation or covariance matrix. \cr
#' \code{coru(), coruv(), coruh()}   \tab uniform correlation matrix, v for common variance,h for heterogeneous variance. \cr
#' \code{diag()}          \tab diagonal variance matrix. \cr
#' \code{grmk()}          \tab the kth GRM matrix. \cr
#' \code{mrmk()}          \tab multiple relationship matrices. \cr
#' \code{us()}            \tab unstructured variance matrix.\cr 
#' \code{facvk()}         \tab factor analytic (basic form). \cr
#' \code{xfak()}          \tab factor analytic (eXtended form). \cr
#' \code{rrk()}           \tab factor analytic: No specific variance. 
#' }
#' 
#' @author Yuanzhen Lin <yzhlinscau@@163.com>
#' @references
#' Yuanzhen Lin. R & ASReml-R Statistics. China Forestry Publishing House. 2016

#' @name Echidna.options

NULL
yzhlinscau/AAfun0s documentation built on April 18, 2023, 4:11 p.m.