## ----initialsetup, include=FALSE----------------------------------------------
knitr::opts_chunk$set(cache=FALSE)

## ----eval=FALSE---------------------------------------------------------------
#  library(data.table)
#  setwd("/mnt/scratch/quantgen/TENSOR_EVD/pipeline")
#  
#  PHENO <- read.csv("source/PHENO.csv")
#  GENO <- fread("source/GENO.csv", data.table=FALSE)
#  ECOV <- read.csv("source/ECOV.csv", row.names=1)
#  
#  # Select North region
#  PHENO <- PHENO[PHENO$region %in% 'North',]
#  PHENO$year_loc <- factor(as.character(PHENO$year_loc))
#  PHENO$genotype <- factor(as.character(PHENO$genotype))
#  
#  save(PHENO, file="data/pheno.RData")
#  
#  # Calculate the GRM
#  ID <- GENO[,1]
#  GENO <- as.matrix(GENO[,-1])
#  rownames(GENO) <- ID
#  X <- scale(GENO, center=TRUE, scale=FALSE)
#  KG <- tcrossprod(X)
#  KG <- KG[levels(PHENO$genotype),levels(PHENO$genotype)]
#  KG <- KG/mean(diag(KG))
#  
#  save(KG, file="data/GRM.RData")
#  
#  # Calculate the ERM
#  ECOV <- ECOV[,-grep("HI30_",colnames(ECOV))]
#  
#  KE <- tcrossprod(scale(ECOV))
#  KE <- KE[levels(PHENO$year_loc),levels(PHENO$year_loc)]
#  KE <- KE/mean(diag(KE))
#  
#  save(KE, file="data/ERM.RData")

## ----eval=FALSE---------------------------------------------------------------
#  JOBS <- expand.grid(nG = c(100,500,1000),
#                      nE = c(10,30,50),
#                      n = c(10000,20000,30000),
#                      alpha = c(0.90,0.95,0.98),
#                      replicate = 1:10)
#  dim(JOBS); head(JOBS)
#  #[1] 810   5
#  #    nG nE     n alpha replicate
#  #1  100 10 10000   0.9         1
#  #2  500 10 10000   0.9         1
#  #3 1000 10 10000   0.9         1
#  #4  100 30 10000   0.9         1
#  #5  500 30 10000   0.9         1
#  #6 1000 30 10000   0.9         1
#  
#  save(JOBS, file="/mnt/scratch/quantgen/TENSOR_EVD/pipeline/parms/JOBS1.RData")

## ----eval=TRUE----------------------------------------------------------------
out <- read.csv(url("https://raw.githubusercontent.com/MarcooLopez/tensorEVD/main/inst/extdata/results_simulation.txt"))
head(out[,1:7])
head(out[,8:12])   # results from the eigen function
head(out[,13:17])  # results from the tensorEVD function

## ----eval=T-------------------------------------------------------------------
source("https://raw.githubusercontent.com/MarcooLopez/tensorEVD/main/misc/functions.R")

## ----eval=T, fig.width = 7, fig.height = 5, fig.align='center'----------------
# Some data edits
out$alpha <- factor(100*out$alpha)

out$nG <- paste0("n[G]*' = '*",out$nG,"L")
out$nE <- paste0("n[E]*' = '*",out$nE,"L")
out$nG <- factor(out$nG, levels = unique(out$nG))
out$nE <- factor(out$nE, levels = unique(out$nE))

# Reshaping the data
measure <- c("time","Frobenius","CMD","nPC","pPC")
dat <- melt_data(out, id=c("nG","nE","n","alpha"),
                 measure=paste0(measure,"_"),
                 value.name=measure, variable.name="method")

color1 <- c('90%'="navajowhite2", '95%'="chocolate1", '98%'="red4")
color2 <- c(eigen="#E69F00", tensorEVD="#009E73", eigs="#56B4E9",
            trlan="#CC79A7", chol="#D55E00")

# Figure 1: Computation time ratio (eigen/tensorEVD)                  
dat0 <- out[out$alpha != "100",]
dat0$alpha <- factor(paste0(dat0$alpha,"%"))
dat0$ratio <- log10(dat0$time_eigen/dat0$time_tensorEVD)
dat0$n <- dat0$n/1000
breaks0 <- seq(1,4,by=1)

figure1 <- make_plot(dat0, type="line", x='n', y='ratio', group="alpha",
                     group.label=NULL, facet="nG", facet2="nE", facet.type="grid",
                     xlab="Sample size (x1000)",
                     ylab="Computation time ratio (eigen/tensorEVD)",
                     group.color=color1, nSD=0, errorbar.size=0,
                     breaks.y=breaks0, labels.y=sprintf("%.f",10^breaks0),
                     scales="fixed")
#print(figure1)

## ----eval=T, fig.width = 7, fig.height = 5, fig.align='center'----------------
dat0 <- dat[dat$method %in% c("eigen","tensorEVD") & dat$alpha!="100",]
dat0$method <- factor(as.character(dat0$method))
dat0$alpha <- factor(as.character(dat0$alpha))

# Figure 2: Approximation accuracy using Frobenious norm
figure2 <- make_plot(dat0, x='alpha', y='Frobenius',
                     group="method", by="n", facet="nG", facet2="nE", facet.type="grid",
                     xlab=bquote(alpha~"x100% of variance of K"),
                     ylab=expression("Frobenius norm ("~abs(abs(K-hat(K)))[F]~")"),
                     by.label="Sample size", breaks.y=seq(0,500,by=100),
                     group.color=color2, rect.by.height=-0.05, ylim=c(0,NA))
#print(figure2)

## ----eval=T, fig.width = 7, fig.height = 5, fig.align='center'----------------
figure3 <- make_plot(dat0, x='alpha', y='pPC',
                     group="method", by="n", facet="nG", facet2="nE", facet.type="grid",
                     xlab=bquote(alpha~"x100% of variance of K"),
                     ylab="Number of eigenvectors/rank", by.label="Sample size",
                     group.color=color2, rect.by.height=-0.05,
                     hline=1, hline.color="red2", ylim=c(0,NA))
#print(figure3)

## ----eval=FALSE---------------------------------------------------------------
#  JOBS <- expand.grid(alpha = c(0.90,0.95,0.98,1.00),
#                      replicate = 1:5)
#  dim(JOBS); head(JOBS)
#  #[1] 20  2
#  #  alpha replicate
#  #1  0.90         1
#  #2  0.95         1
#  #3  0.98         1
#  #4  1.00         1
#  #5  0.90         2
#  #6  0.95         2
#  
#  save(JOBS, file="/mnt/scratch/quantgen/TENSOR_EVD/pipeline/parms/JOBS2.RData")

## ----eval=FALSE---------------------------------------------------------------
#  JOBS <- expand.grid(trait = c("yield","anthesis","silking","ASI"),
#                      method = c("eigen","tensorEVD"),
#                      alpha = c(0.90,0.95,0.98,1.00),
#                      replicate = 1:5)
#  JOBS <- JOBS[-which(JOBS$alpha==1.00 & JOBS$method=="tensorEVD"),]
#  dim(JOBS); head(JOBS)
#  #[1] 140   4
#  #     trait    method alpha replicate
#  #1    yield     eigen   0.9         1
#  #2 anthesis     eigen   0.9         1
#  #3  silking     eigen   0.9         1
#  #4      ASI     eigen   0.9         1
#  #5    yield tensorEVD   0.9         1
#  #6 anthesis tensorEVD   0.9         1
#  
#  save(JOBS, file="/mnt/scratch/quantgen/TENSOR_EVD/pipeline/parms/JOBS3.RData")

## ----eval=FALSE---------------------------------------------------------------
#  setwd("/mnt/scratch/quantgen/TENSOR_EVD/pipeline")
#  load("parms/JOBS3.RData")
#  prefix <- "output/genomic_prediction/ANOVA"
#  
#  out <- c()
#  for(k in 1:nrow(JOBS))
#  {
#    trait <- as.vector(JOBS[k,"trait"])
#    method <- as.vector(JOBS[k,"method"])
#    alpha <- as.vector(JOBS[k,"alpha"])
#    replicate <- as.vector(JOBS[k,"replicate"])
#  
#    suffix <- paste0(trait,"/",method,"/alpha_",100*alpha,"/rep_",replicate,"/VC.RData")
#    filename <- paste0(prefix,"/",suffix)
#    if(file.exists(filename)){
#      load(filename)
#      out <- rbind(out, VC)
#    }else{
#      message("File not found: '",suffix,"'")
#    }
#  }

## ----eval=TRUE----------------------------------------------------------------
out <- read.csv(url("https://raw.githubusercontent.com/MarcooLopez/tensorEVD/main/inst/extdata/results_ANOVA.txt"))
head(out)

## ----eval=T, fig.width = 6.5, fig.height = 4.5, fig.align='center'------------
out$alpha <- factor(paste0(100*out$alpha,"%"), levels=c("100%","98%","95%","90%"))
out$source <- factor(out$source, levels=c("G","E","GE","Error"))

trait <- c("yield", "anthesis", "silking", "ASI")[1]

myfun <- function(x) sprintf('%.3f', x)

# Figure 4: Phenotypic variance of yield
dat <- out[out$trait==trait,]
figure4 <- make_plot(dat, x='alpha', y='mean', SD="SD",
                     group="method", facet="source",
                     xlab=bquote(alpha~"x100% of variance of K"),
                     ylab=paste0("Proportion of variance of ",trait),
                     group.color=color2, scales="free_y",
                     ylabels=myfun, text=myfun, ylim=c(0,NA))
#print(figure4)

## ----eval=FALSE---------------------------------------------------------------
#  setwd("/mnt/scratch/quantgen/TENSOR_EVD/pipeline")
#  load("data/pheno.RData")
#  table(PHENO$CV_10fold)
#  #   1    2    3    4    5    6    7    8    9   10
#  #6180 6277 6246 5785 6160 5858 5492 5660 5436 5975

## ----eval=FALSE---------------------------------------------------------------
#  JOBS <- expand.grid(trait = c("yield","anthesis","silking","ASI"),
#                      method = c("eigen","tensorEVD"),
#                      alpha = c(0.90,0.95,0.98),
#                      fold = 1:10)
#  dim(JOBS); head(JOBS)
#  #[1] 240   4
#  #     trait    method alpha fold
#  #1    yield     eigen   0.9    1
#  #2 anthesis     eigen   0.9    1
#  #3  silking     eigen   0.9    1
#  #4      ASI     eigen   0.9    1
#  #5    yield tensorEVD   0.9    1
#  #6 anthesis tensorEVD   0.9    1
#  
#  save(JOBS, file="/mnt/scratch/quantgen/TENSOR_EVD/pipeline/parms/JOBS4.RData")

## ----eval=FALSE---------------------------------------------------------------
#  setwd("/mnt/scratch/quantgen/TENSOR_EVD/pipeline")
#  source("https://raw.githubusercontent.com/MarcooLopez/tensorEVD/main/misc/functions.R")
#  load("parms/JOBS4.RData")
#  prefix <- "output/genomic_prediction/10F_CV"
#  
#  dat <- c()
#  for(trait in levels(JOBS$trait)){
#    for(method in levels(JOBS$method)){
#      for(alpha in unique(JOBS$alpha)){
#        out0 <- c()
#        for(fold in unique(JOBS$fold)){
#          suffix <- paste0(trait,"/",method,"/alpha_",100*alpha,"/results_fold_",fold,".RData")
#          filename <- paste0(prefix,"/",suffix)
#          if(file.exists(filename)){
#            load(filename)
#            out0 <- rbind(out0, out)
#          }else{
#            message("File not found: '",suffix,"'")
#          }
#        }
#        tmp <- get_corr(out0, by="year_loc")
#        dat <- rbind(dat, data.frame(trait,method,alpha,tmp))
#      }
#    }
#  }
#  
#  # Reshaping the data
#  dat$trait <- factor(dat$trait, levels=levels(JOBS$trait))
#  out <- reshape2::dcast(dat, trait+alpha+year_loc+nRecords~method, value.var="correlation")
#  tmp <- reshape2::dcast(dat, trait+alpha+year_loc+nRecords~method, value.var="SE")[,levels(JOBS$method)]
#  colnames(tmp) <- paste0(colnames(tmp),".SE")
#  out <- data.frame(out, tmp)

## ----eval=TRUE----------------------------------------------------------------
out <- read.csv(url("https://raw.githubusercontent.com/MarcooLopez/tensorEVD/main/inst/extdata/results_10F_CV.txt"))
head(out)

## ----eval=T, fig.width = 6.0, fig.height = 6.5, fig.align='center'------------
out$trait <- factor(out$trait, levels=unique(out$trait))
out$alpha <- factor(paste0(100*out$alpha,"%"), levels=c("98%","95%","90%"))

# Figure 5: Within environment prediction correlation
rg <- range(c(out$eigen,out$tensorEVD))
if(requireNamespace("ggplot2", quietly=TRUE)){
  figure5 <-  ggplot2::ggplot(out, ggplot2::aes(tensorEVD, eigen)) +
              ggplot2::geom_abline(color="gray70", linetype="dashed") +
              ggplot2::geom_point(fill="#56B4E9", shape=21, size=1.4) +
              ggplot2::facet_grid(trait ~ alpha) +
              ggplot2::theme_bw() + ggplot2::xlim(rg) + ggplot2::ylim(rg)
}

#print(figure5)

