not_in_pkg/summarize_sim_runs.R

# Nathan dot lazar at gmail dot com

args <- commandArgs(TRUE)

run_prefix <- args[[1]]

library(dplyr)
library(methods)
library(R6)
library(BaTFLED3D)

# Cut offs used when determining sparsity
m1.rem.cut <- 1e-5
m2.rem.cut <- 1e-5
m3.rem.cut <- 1e-5
core.rem.cut <- 1e-3

# Functions
#############################################################
rmse <- function(obs, pred) sqrt(mean((obs-pred)^2,na.rm=T))

# Read in a file to get the number of iterations
get_iters <- function(f1) {
  load(f1)
  ret <- list()
  return(trained$iter)
}

# Function to load in data from runs
loadData <- function(f, results){
  #loads an RData file, and returns a list with the 
  #trained model and warm & cold RMSE vectors.
  
  load(f)
  
  if(!exists('all.resp.train'))
    all.resp.train <- all.train.resp
  
  m1.means <- apply(train.data$resp, c(2,3), mean, na.rm=T)
  m2.means <- apply(train.data$resp, c(1,3), mean, na.rm=T)
  m3.means <- apply(train.data$resp, c(1,2), mean, na.rm=T)
  
  # Get std. devs of responses
  stddev <- sd(train.data$resp, na.rm=T)
  
  # Make a tensor of the means for each value of mode2 & mode3 for comparisons
  m1.mean.tens <- train.data$resp
  m2.mean.tens <- train.data$resp
  m3.mean.tens <- train.data$resp
  m1m2.mean.tens <- train.data$resp
  m1m3.mean.tens <- train.data$resp
  m2m3.mean.tens <- train.data$resp
  m1m2m3.mean.tens <- train.data$resp
  for(i in 1:dim(m1.mean.tens)[1]) m1.mean.tens[i,,] <- m1.means
  for(j in 1:dim(m2.mean.tens)[2]) m2.mean.tens[,j,] <- m2.means
  for(k in 1:dim(m2.mean.tens)[3]) m3.mean.tens[,,k] <- m3.means
  for(k in 1:dim(m1m2.mean.tens)[3]) m1m2.mean.tens[,,k] <- m3.means
  for(j in 1:dim(m1m3.mean.tens)[2]) m1m3.mean.tens[,j,] <- m2.means
  for(i in 1:dim(m2m3.mean.tens)[1]) m2m3.mean.tens[i,,] <- m1.means
  m1m2m3.mean.tens[,,] <- mean(train.data$resp, na.rm=T)

  results$mean['RMSE', 'train.m1', fold] <- rmse(train.data$resp, m1.mean.tens)/stddev
  results$mean['RMSE', 'train.m2', fold] <- rmse(train.data$resp, m2.mean.tens)/stddev
  results$mean['RMSE', 'train.m3', fold] <- rmse(train.data$resp, m3.mean.tens)/stddev
  results$mean['RMSE', 'train.m1m2', fold] <- rmse(train.data$resp, m1m2.mean.tens)/stddev
  results$mean['RMSE', 'train.m1m3', fold] <- rmse(train.data$resp, m1m3.mean.tens)/stddev
  results$mean['RMSE', 'train.m2m3', fold] <- rmse(train.data$resp, m2m3.mean.tens)/stddev
  results$mean['RMSE', 'train.m1m2m3', fold] <- rmse(train.data$resp, m1m2m3.mean.tens)/stddev
  
  results$mean['RMSE', 'warm.m1', fold] <- rmse(all.resp.train[is.na(train.data$resp)],
                                                m1.mean.tens[is.na(train.data$resp)])/stddev
  results$mean['RMSE', 'warm.m2', fold] <- rmse(all.resp.train[is.na(train.data$resp)],
                                                m2.mean.tens[is.na(train.data$resp)])/stddev
  results$mean['RMSE', 'warm.m3', fold] <- rmse(all.resp.train[is.na(train.data$resp)],
                                                m3.mean.tens[is.na(train.data$resp)])/stddev
  results$mean['RMSE', 'warm.m1m2', fold] <- rmse(all.resp.train[is.na(train.data$resp)], 
                                                  m1m2.mean.tens[is.na(train.data$resp)])/stddev
  results$mean['RMSE', 'warm.m1m3', fold] <- rmse(all.resp.train[is.na(train.data$resp)], 
                                                  m1m3.mean.tens[is.na(train.data$resp)])/stddev
  results$mean['RMSE', 'warm.m2m3', fold] <- rmse(all.resp.train[is.na(train.data$resp)], 
                                                  m2m3.mean.tens[is.na(train.data$resp)])/stddev
  results$mean['RMSE', 'warm.m1m2m3', fold] <- rmse(all.resp.train[is.na(train.data$resp)], 
                                                    m1m2m3.mean.tens[is.na(train.data$resp)])/stddev
  
  results$mean['RMSE', 'm1', fold] <- rmse(test.data.m1$resp,
                                           m1.mean.tens[1:dim(test.data.m1$resp)[1],,])/stddev
  results$mean['RMSE', 'm2', fold] <- rmse(test.data.m2$resp,
                                           m2.mean.tens[,1:dim(test.data.m2$resp)[2],])/stddev
  results$mean['RMSE', 'm3', fold] <- rmse(test.data.m3$resp,
                                           m3.mean.tens[,,1:dim(test.data.m3$resp)[3]])/stddev
  results$mean['RMSE', 'm1m2', fold] <- rmse(test.data.m1m2$resp,
         m1m2.mean.tens[1:dim(test.data.m1m2$resp)[1], 1:dim(test.data.m1m2$resp)[2],])/stddev
  results$mean['RMSE', 'm1m3', fold] <- rmse(test.data.m1m3$resp,
         m1m3.mean.tens[1:dim(test.data.m1m3$resp)[1],, 1:dim(test.data.m1m3$resp)[3]])/stddev
  results$mean['RMSE', 'm2m3', fold] <- rmse(test.data.m2m3$resp,
         m2m3.mean.tens[,1:dim(test.data.m2m3$resp)[2], 1:dim(test.data.m2m3$resp)[3]])/stddev
  results$mean['RMSE', 'm1m2m3', fold] <- rmse(test.data.m1m2m3$resp,
         m1m2m3.mean.tens[1:dim(test.data.m1m2m3$resp)[1],
                          1:dim(test.data.m1m2m3$resp)[2],
                          1:dim(test.data.m1m2m3$resp)[3]])/stddev

  results$mean['exp.var', 'train.m1', fold] <- exp_var(train.data$resp, m1.mean.tens)
  results$mean['exp.var', 'train.m2', fold] <- exp_var(train.data$resp, m2.mean.tens)
  results$mean['exp.var', 'train.m3', fold] <- exp_var(train.data$resp, m3.mean.tens)
  results$mean['exp.var', 'train.m1m2', fold] <- exp_var(train.data$resp, m1m2.mean.tens)
  results$mean['exp.var', 'train.m1m3', fold] <- exp_var(train.data$resp, m1m3.mean.tens)
  results$mean['exp.var', 'train.m2m3', fold] <- exp_var(train.data$resp, m2m3.mean.tens)
  results$mean['exp.var', 'train.m1m2m3', fold] <- exp_var(train.data$resp, m1m2m3.mean.tens)
  
  results$mean['exp.var', 'warm.m1', fold] <- exp_var(all.resp.train[is.na(train.data$resp)],
                                                      m1.mean.tens[is.na(train.data$resp)])
  results$mean['exp.var', 'warm.m2', fold] <- exp_var(all.resp.train[is.na(train.data$resp)],
                                                      m2.mean.tens[is.na(train.data$resp)])
  results$mean['exp.var', 'warm.m3', fold] <- exp_var(all.resp.train[is.na(train.data$resp)],
                                                      m3.mean.tens[is.na(train.data$resp)])
  results$mean['exp.var', 'warm.m1m2', fold] <- exp_var(all.resp.train[is.na(train.data$resp)], 
                                                        m1m2.mean.tens[is.na(train.data$resp)])
  results$mean['exp.var', 'warm.m1m3', fold] <- exp_var(all.resp.train[is.na(train.data$resp)], 
                                                        m1m3.mean.tens[is.na(train.data$resp)])
  results$mean['exp.var', 'warm.m2m3', fold] <- exp_var(all.resp.train[is.na(train.data$resp)], 
                                                        m2m3.mean.tens[is.na(train.data$resp)])
  results$mean['exp.var', 'warm.m1m2m3', fold] <- exp_var(all.resp.train[is.na(train.data$resp)], 
                                                          m1m2m3.mean.tens[is.na(train.data$resp)])

  results$mean['exp.var', 'm1', fold] <- exp_var(test.data.m1$resp,
                                                 m1.mean.tens[1:dim(test.data.m1$resp)[1],,])
  results$mean['exp.var', 'm2', fold] <- exp_var(test.data.m2$resp,
                                                 m2.mean.tens[,1:dim(test.data.m2$resp)[2],])
  results$mean['exp.var', 'm3', fold] <- exp_var(test.data.m3$resp,
                                                 m3.mean.tens[,,1:dim(test.data.m3$resp)[3]])
  results$mean['exp.var', 'm1m2', fold] <- exp_var(test.data.m1m2$resp,
           m1m2.mean.tens[1:dim(test.data.m1m2$resp)[1], 1:dim(test.data.m1m2$resp)[2],])
  results$mean['exp.var', 'm1m3', fold] <- exp_var(test.data.m1m3$resp,
           m1m3.mean.tens[1:dim(test.data.m1m2$resp)[1],, 1:dim(test.data.m1m2$resp)[3]])
  results$mean['exp.var', 'm2m3', fold] <- exp_var(test.data.m2m3$resp,
           m2m3.mean.tens[,1:dim(test.data.m2m3$resp)[2], 1:dim(test.data.m2m3$resp)[3]])
  results$mean['exp.var', 'm1m2m3', fold] <- exp_var(test.data.m1m2m3$resp,
           m1m2m3.mean.tens[1:dim(test.data.m1m2m3$resp)[1],
                            1:dim(test.data.m1m2m3$resp)[2],
                            1:dim(test.data.m1m2m3$resp)[3]])
  
  results$mean['p.cor', 'train.m1', fold] <- cor(train.data$resp, m1.mean.tens, use='complete.obs')
  results$mean['p.cor', 'train.m2', fold] <- cor(train.data$resp, m2.mean.tens, use='complete.obs')
  results$mean['p.cor', 'train.m3', fold] <- cor(train.data$resp, m3.mean.tens, use='complete.obs')
  results$mean['p.cor', 'train.m1m2', fold] <- cor(train.data$resp, m1m2.mean.tens, use='complete.obs')
  results$mean['p.cor', 'train.m1m3', fold] <- cor(train.data$resp, m1m3.mean.tens, use='complete.obs')
  results$mean['p.cor', 'train.m2m3', fold] <- cor(train.data$resp, m2m3.mean.tens, use='complete.obs')
  results$mean['p.cor', 'train.m1m2m3', fold] <- cor(train.data$resp, m1m2m3.mean.tens, use='complete.obs')
  
  results$mean['p.cor', 'warm.m1', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)],
             m1.mean.tens[is.na(train.data$resp)], use='complete.obs'), error=function(e) NA)
  results$mean['p.cor', 'warm.m2', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)],
             m2.mean.tens[is.na(train.data$resp)], use='complete.obs'), error=function(e) NA)
  results$mean['p.cor', 'warm.m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)],
             m3.mean.tens[is.na(train.data$resp)], use='complete.obs'), error=function(e) NA)
  results$mean['p.cor', 'warm.m1m2', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
             m1m2.mean.tens[is.na(train.data$resp)], use='complete.obs'), error=function(e) NA)
  results$mean['p.cor', 'warm.m1m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
             m1m3.mean.tens[is.na(train.data$resp)], use='complete.obs'), error=function(e) NA)
  results$mean['p.cor', 'warm.m2m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
             m2m3.mean.tens[is.na(train.data$resp)], use='complete.obs'), error=function(e) NA)
  results$mean['p.cor', 'warm.m1m2m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
             m1m2m3.mean.tens[is.na(train.data$resp)], use='complete.obs'), error=function(e) NA)
  
  results$mean['p.cor', 'm1', fold] <- cor(test.data.m1$resp,
             m1.mean.tens[1:dim(test.data.m1$resp)[1],,], use='complete.obs')
  results$mean['p.cor', 'm2', fold] <- cor(test.data.m2$resp,
             m2.mean.tens[,1:dim(test.data.m2$resp)[2],], use='complete.obs')
  results$mean['p.cor', 'm3', fold] <- cor(test.data.m3$resp,
             m3.mean.tens[,,1:dim(test.data.m3$resp)[3]], use='complete.obs')
  results$mean['p.cor', 'm1m2', fold] <- cor(test.data.m1m2$resp,
             m1m2.mean.tens[1:dim(test.data.m1m2$resp)[1], 1:dim(test.data.m1m2$resp)[2],], use='complete.obs')
  results$mean['p.cor', 'm1m3', fold] <- cor(test.data.m1m3$resp,
             m1m3.mean.tens[1:dim(test.data.m1m3$resp)[1],, 1:dim(test.data.m1m3$resp)[3]], use='complete.obs')
  results$mean['p.cor', 'm2m3', fold] <- cor(test.data.m2m3$resp,
             m2m3.mean.tens[,1:dim(test.data.m2m3$resp)[2], 1:dim(test.data.m2m3$resp)[3]], use='complete.obs')
  results$mean['p.cor', 'm1m2m3', fold] <- cor(test.data.m1m2m3$resp,
             m1m2m3.mean.tens[1:dim(test.data.m1m2m3$resp)[1], 
                              1:dim(test.data.m1m2m3$resp)[2],
                              1:dim(test.data.m1m2m3$resp)[3]], use='complete.obs')
  
  results$mean['s.cor', 'train.m1', fold] <- cor(train.data$resp, m1.mean.tens, use='complete.obs', method='spearman')
  results$mean['s.cor', 'train.m2', fold] <- cor(train.data$resp, m2.mean.tens, use='complete.obs', method='spearman')
  results$mean['s.cor', 'train.m3', fold] <- cor(train.data$resp, m3.mean.tens, use='complete.obs', method='spearman')
  results$mean['s.cor', 'train.m1m2', fold] <- cor(train.data$resp, m1m2.mean.tens, use='complete.obs', method='spearman')
  results$mean['s.cor', 'train.m1m3', fold] <- cor(train.data$resp, m1m3.mean.tens, use='complete.obs', method='spearman')
  results$mean['s.cor', 'train.m2m3', fold] <- cor(train.data$resp, m2m3.mean.tens, use='complete.obs', method='spearman')
  results$mean['s.cor', 'train.m1m2m3', fold] <- cor(train.data$resp, m1m2m3.mean.tens, use='complete.obs', method='spearman')
  
  results$mean['s.cor', 'warm.m1', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)],
            m1.mean.tens[is.na(train.data$resp)], use='complete.obs', method='spearman'), error=function(e) NA)
  results$mean['s.cor', 'warm.m2', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)],
            m2.mean.tens[is.na(train.data$resp)], use='complete.obs', method='spearman'), error=function(e) NA)
  results$mean['s.cor', 'warm.m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)],
            m3.mean.tens[is.na(train.data$resp)], use='complete.obs', method='spearman'), error=function(e) NA)
  results$mean['s.cor', 'warm.m1m2', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
            m1m2.mean.tens[is.na(train.data$resp)], use='complete.obs', method='spearman'), error=function(e) NA)
  results$mean['s.cor', 'warm.m1m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
            m1m3.mean.tens[is.na(train.data$resp)], use='complete.obs', method='spearman'), error=function(e) NA)
  results$mean['s.cor', 'warm.m2m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
            m2m3.mean.tens[is.na(train.data$resp)], use='complete.obs', method='spearman'), error=function(e) NA)
  results$mean['s.cor', 'warm.m1m2m3', fold] <- tryCatch(cor(all.resp.train[is.na(train.data$resp)], 
            m1m2m3.mean.tens[is.na(train.data$resp)], use='complete.obs', method='spearman'), error=function(e) NA)
  
  results$mean['s.cor', 'm1', fold] <- cor(test.data.m1$resp,
            m1.mean.tens[1:dim(test.data.m1$resp)[1],,], use='complete.obs', method='spearman')
  results$mean['s.cor', 'm2', fold] <- cor(test.data.m2$resp,
            m2.mean.tens[,1:dim(test.data.m2$resp)[2],], use='complete.obs', method='spearman')
  results$mean['s.cor', 'm3', fold] <- cor(test.data.m3$resp,
            m3.mean.tens[,,1:dim(test.data.m3$resp)[3]], use='complete.obs', method='spearman')
  results$mean['s.cor', 'm1m2', fold] <- cor(test.data.m1m2$resp,
            m1m2.mean.tens[1:dim(test.data.m1m2$resp)[1], 1:dim(test.data.m1m2$resp)[2],], use='complete.obs', method='spearman')
  results$mean['s.cor', 'm1m3', fold] <- cor(test.data.m1m3$resp,
            m1m3.mean.tens[1:dim(test.data.m1m3$resp)[1],, 1:dim(test.data.m1m3$resp)[3]], use='complete.obs', method='spearman')
  results$mean['s.cor', 'm2m3', fold] <- cor(test.data.m2m3$resp,
            m2m3.mean.tens[,1:dim(test.data.m2m3$resp)[2], 1:dim(test.data.m2m3$resp)[3]], use='complete.obs', method='spearman')
  results$mean['p.cor', 'm1m2m3', fold] <- cor(test.data.m1m2m3$resp,
            m1m2m3.mean.tens[1:dim(test.data.m1m2m3$resp)[1], 
                             1:dim(test.data.m1m2m3$resp)[2],
                             1:dim(test.data.m1m2m3$resp)[3]], use='complete.obs', method='spearman')
  
  results$training['lower.bnd', fold, 1:trained$iter] <- trained$lower.bnd
  
  results$lower.bnd['final', fold] <- trained$lower.bnd[trained$iter]
  results$lower.bnd['max', fold] <- max(trained$lower.bnd)
  results$lower.bnd['which.max', fold] <- which.max(trained$lower.bnd)
  results$lower.bnd['not.mono', fold] <- sum((trained$lower.bnd[-1] -
                                                trained$lower.bnd[-trained$iter]) < 0)
  
  results$training['A.RMSE', fold, 1:trained$iter] <- trained$RMSE
  results$training['H.RMSE', fold, 1:trained$iter] <- trained$H.RMSE
  results$training['exp.var', fold, 1:trained$iter] <- trained$exp.var
  if(length(trained$p.cor))
    results$training['p.cor', fold, 1:trained$iter] <- trained$p.cor
  if(length(trained$s.cor))
    results$training['s.cor', fold, 1:trained$iter] <- trained$s.cor
  if('warm.RMSE' %in% names(test.results)) {
    results$training['warm.RMSE', fold, 1:trained$iter] <- test.results$warm.RMSE
    results$training['warm.exp.var', fold, 1:trained$iter] <- test.results$warm.exp.var
    if(length(test.results$warm.p.cor))
      results$training['warm.p.cor', fold, 1:trained$iter] <- test.results$warm.p.cor
    if(length(test.results$warm.s.cor))
      results$training['warm.s.cor', fold, 1:trained$iter] <- test.results$warm.s.cor
  }
  if('m1.RMSE' %in% names(test.results)) {
    results$training['m1.RMSE', fold, 1:trained$iter] <- test.results$m1.RMSE
    results$training['m1.exp.var', fold, 1:trained$iter] <- test.results$m1.exp.var
    if(length(test.results$m1.p.cor))
      results$training['m1.p.cor', fold, 1:trained$iter] <- test.results$m1.p.cor
    if(length(test.results$m1.s.cor))
      results$training['m1.s.cor', fold, 1:trained$iter] <- test.results$m1.s.cor
  }
  if('m2.RMSE' %in% names(test.results)) {
    results$training['m2.RMSE', fold, 1:trained$iter] <- test.results$m2.RMSE
    results$training['m2.exp.var', fold, 1:trained$iter] <- test.results$m2.exp.var
    if(length(test.results$m2.p.cor))
      results$training['m2.p.cor', fold, 1:trained$iter] <- test.results$m2.p.cor
    if(length(test.results$m2.s.cor))
      results$training['m2.s.cor', fold, 1:trained$iter] <- test.results$m2.s.cor
  }
  if('m3.RMSE' %in% names(test.results)) {
    results$training['m3.RMSE', fold, 1:trained$iter] <- test.results$m3.RMSE
    results$training['m3.exp.var', fold, 1:trained$iter] <- test.results$m3.exp.var
    if(length(test.results$m3.p.cor))
      results$training['m3.p.cor', fold, 1:trained$iter] <- test.results$m3.p.cor
    if(length(test.results$m3.s.cor))
      results$training['m3.s.cor', fold, 1:trained$iter] <- test.results$m3.s.cor
  }
  if('m1m2.RMSE' %in% names(test.results)) {
    results$training['m1m2.RMSE', fold, 1:trained$iter] <- test.results$m1m2.RMSE
    results$training['m1m2.exp.var', fold, 1:trained$iter] <- test.results$m1m2.exp.var
    if(length(test.results$m1m2.p.cor))
      results$training['m1m2.p.cor', fold, 1:trained$iter] <- test.results$m1m2.p.cor
    if(length(test.results$m1m2.s.cor))
      results$training['m1m2.s.cor', fold, 1:trained$iter] <- test.results$m1m2.s.cor
  }
  if('m1m3.RMSE' %in% names(test.results)) {
    results$training['m1m3.RMSE', fold, 1:trained$iter] <- test.results$m1m3.RMSE
    results$training['m1m3.exp.var', fold, 1:trained$iter] <- test.results$m1m3.exp.var
    if(length(test.results$m1m3.p.cor))
      results$training['m1m3.p.cor', fold, 1:trained$iter] <- test.results$m1m3.p.cor
    if(length(test.results$m1m3.s.cor))
      results$training['m1m3.s.cor', fold, 1:trained$iter] <- test.results$m1m3.s.cor
  }
  if('m2m3.RMSE' %in% names(test.results)) {
    results$training['m2m3.RMSE', fold, 1:trained$iter] <- test.results$m2m3.RMSE
    results$training['m2m3.exp.var', fold, 1:trained$iter] <- test.results$m2m3.exp.var
    if(length(test.results$m2m3.p.cor))
      results$training['m2m3.p.cor', fold, 1:trained$iter] <- test.results$m2m3.p.cor
    if(length(test.results$m2m3.s.cor))
      results$training['m2m3.s.cor', fold, 1:trained$iter] <- test.results$m2m3.s.cor
  }
  
  if('m1m2m3.RMSE' %in% names(test.results)) {
    results$training['m1m2m3.RMSE', fold, 1:trained$iter] <- test.results$m1m2m3.RMSE
    results$training['m1m2m3.exp.var', fold, 1:trained$iter] <- test.results$m1m2m3.exp.var
    if(length(test.results$m1m2m3.p.cor))
      results$training['m1m2m3.p.cor', fold, 1:trained$iter] <- test.results$m1m2m3.p.cor
    if(length(test.results$m1m2m3.s.cor))
      results$training['m1m2m3.s.cor', fold, 1:trained$iter] <- test.results$m1m2m3.s.cor
  }
  
  results$summaries['H', 'RMSE', fold] <- trained$H.RMSE[trained$iter]
  results$summaries['H', 'min.RMSE.iter', fold] <- which.min(trained$H.RMSE)
  results$summaries['H', 'min.RMSE', fold] <- min(trained$H.RMSE)
  
  results$summaries['train', 'RMSE', fold] <- trained$RMSE[trained$iter]
  results$summaries['train', 'min.RMSE.iter', fold] <- which.min(trained$RMSE)
  results$summaries['train', 'min.RMSE', fold] <- min(trained$RMSE)
  results$summaries['train', 'exp.var', fold] <- trained$exp.var[trained$iter]
  results$summaries['train', 'max.exp.var.iter', fold] <- which.max(trained$exp.var)
  results$summaries['train', 'max.exp.var', fold] <- max(trained$exp.var)
  if(length(trained$p.cor)) {
    results$summaries['train', 'p.cor', fold] <- trained$p.cor[trained$iter]
    results$summaries['train', 'max.p.cor.iter', fold] <- which.max(trained$p.cor)
    results$summaries['train', 'max.p.cor', fold] <- max(trained$p.cor)
    results$summaries['train', 's.cor', fold] <- trained$s.cor[trained$iter]
    results$summaries['train', 'max.s.cor.iter', fold] <- which.max(trained$s.cor)
    results$summaries['train', 'max.s.cor', fold] <- max(trained$s.cor)
  }
  
  if('warm.RMSE' %in% names(test.results)) {
    results$summaries['warm', 'RMSE', fold] <- test.results$warm.RMSE[trained$iter]
    results$summaries['warm', 'min.RMSE.iter', fold] <- which.min(test.results$warm.RMSE)
    results$summaries['warm', 'min.RMSE', fold] <- min(test.results$warm.RMSE)
    results$summaries['warm', 'clip.RMSE', fold] <- test.results$warm.RMSE.clip[trained$iter]
    results$summaries['warm', 'min.clip.RMSE.iter', fold] <- which.min(test.results$warm.RMSE.clip)
    results$summaries['warm', 'min.clip.RMSE', fold] <- min(test.results$warm.RMSE.clip)
    results$summaries['warm', 'exp.var', fold] <- test.results$warm.exp.var[trained$iter]
    results$summaries['warm', 'max.exp.var.iter', fold] <- which.max(test.results$warm.exp.var)
    results$summaries['warm', 'max.exp.var', fold] <- max(test.results$warm.exp.var)
    if(length(test.results$warm.p.cor)) {
      results$summaries['warm', 'p.cor', fold] <- test.results$warm.p.cor[trained$iter]
      results$summaries['warm', 'max.p.cor.iter', fold] <- which.max(test.results$warm.p.cor)
      results$summaries['warm', 'max.p.cor', fold] <- max(test.results$warm.p.cor)
      results$summaries['warm', 's.cor', fold] <- test.results$warm.s.cor[trained$iter]
      results$summaries['warm', 'max.s.cor.iter', fold] <- which.max(test.results$warm.s.cor)
      results$summaries['warm', 'max.s.cor', fold] <- max(test.results$warm.s.cor)
    }
  }
  
  if('m1.RMSE' %in% names(test.results)) {
    results$summaries['m1', 'RMSE', fold] <- test.results$m1.RMSE[trained$iter]
    results$summaries['m1', 'min.RMSE.iter', fold] <- which.min(test.results$m1.RMSE)
    results$summaries['m1', 'min.RMSE', fold] <- min(test.results$m1.RMSE)
    results$summaries['m1', 'clip.RMSE', fold] <- test.results$m1.RMSE.clip[trained$iter]
    results$summaries['m1', 'min.clip.RMSE.iter', fold] <- which.min(test.results$m1.RMSE.clip)
    results$summaries['m1', 'min.clip.RMSE', fold] <- min(test.results$m1.RMSE.clip)
    results$summaries['m1', 'exp.var', fold] <- test.results$m1.exp.var[trained$iter]
    results$summaries['m1', 'max.exp.var.iter', fold] <- which.max(test.results$m1.exp.var)
    results$summaries['m1', 'max.exp.var', fold] <- max(test.results$m1.exp.var)
    if(length(test.results$m1.p.cor)) {
      results$summaries['m1', 'p.cor', fold] <- test.results$m1.p.cor[trained$iter]
      results$summaries['m1', 'max.p.cor.iter', fold] <- which.max(test.results$m1.p.cor)
      results$summaries['m1', 'max.p.cor', fold] <- max(test.results$m1.p.cor)
      results$summaries['m1', 's.cor', fold] <- test.results$m1.s.cor[trained$iter]
      results$summaries['m1', 'max.s.cor.iter', fold] <- which.max(test.results$m1.s.cor)
      results$summaries['m1', 'max.s.cor', fold] <- max(test.results$m1.s.cor)
    }
  }
  
  if('m2.RMSE' %in% names(test.results)) {
    results$summaries['m2', 'RMSE', fold] <- test.results$m2.RMSE[trained$iter]
    results$summaries['m2', 'min.RMSE.iter', fold] <- which.min(test.results$m2.RMSE)
    results$summaries['m2', 'min.RMSE', fold] <- min(test.results$m2.RMSE)
    results$summaries['m2', 'clip.RMSE', fold] <- test.results$m2.RMSE.clip[trained$iter]
    results$summaries['m2', 'min.clip.RMSE.iter', fold] <- which.min(test.results$m2.RMSE.clip)
    results$summaries['m2', 'min.clip.RMSE', fold] <- min(test.results$m2.RMSE.clip)
    results$summaries['m2', 'exp.var', fold] <- test.results$m2.exp.var[trained$iter]
    results$summaries['m2', 'max.exp.var.iter', fold] <- which.max(test.results$m2.exp.var)
    results$summaries['m2', 'max.exp.var', fold] <- max(test.results$m2.exp.var)
    if(length(test.results$m2.p.cor)) {
      results$summaries['m2', 'p.cor', fold] <- test.results$m2.p.cor[trained$iter]
      results$summaries['m2', 'max.p.cor.iter', fold] <- which.max(test.results$m2.p.cor)
      results$summaries['m2', 'max.p.cor', fold] <- max(test.results$m2.p.cor)
      results$summaries['m2', 's.cor', fold] <- test.results$m2.s.cor[trained$iter]
      results$summaries['m2', 'max.s.cor.iter', fold] <- which.max(test.results$m2.s.cor)
      results$summaries['m2', 'max.s.cor', fold] <- max(test.results$m2.s.cor)
    }
  }
  
  if('m3.RMSE' %in% names(test.results)) {
    results$summaries['m3', 'RMSE', fold] <- test.results$m3.RMSE[trained$iter]
    results$summaries['m3', 'min.RMSE.iter', fold] <- which.min(test.results$m3.RMSE)
    results$summaries['m3', 'min.RMSE', fold] <- min(test.results$m3.RMSE)
    results$summaries['m3', 'clip.RMSE', fold] <- test.results$m3.RMSE.clip[trained$iter]
    results$summaries['m3', 'min.clip.RMSE.iter', fold] <- which.min(test.results$m3.RMSE.clip)
    results$summaries['m3', 'min.clip.RMSE', fold] <- min(test.results$m3.RMSE.clip)
    results$summaries['m3', 'exp.var', fold] <- test.results$m3.exp.var[trained$iter]
    results$summaries['m3', 'max.exp.var.iter', fold] <- which.max(test.results$m3.exp.var)
    results$summaries['m3', 'max.exp.var', fold] <- max(test.results$m3.exp.var)
    if(length(test.results$m3.p.cor)) {
      results$summaries['m3', 'p.cor', fold] <- test.results$m3.p.cor[trained$iter]
      results$summaries['m3', 'max.p.cor.iter', fold] <- which.max(test.results$m3.p.cor)
      results$summaries['m3', 'max.p.cor', fold] <- max(test.results$m3.p.cor)
      results$summaries['m3', 's.cor', fold] <- test.results$m3.s.cor[trained$iter]
      results$summaries['m3', 'max.s.cor.iter', fold] <- which.max(test.results$m3.s.cor)
      results$summaries['m3', 'max.s.cor', fold] <- max(test.results$m3.s.cor)
    }
  }
  
  if('m1m2.RMSE' %in% names(test.results)) {
    results$summaries['m1m2', 'RMSE', fold] <- test.results$m1m2.RMSE[trained$iter]
    results$summaries['m1m2', 'min.RMSE.iter', fold] <- which.min(test.results$m1m2.RMSE)
    results$summaries['m1m2', 'min.RMSE', fold] <- min(test.results$m1m2.RMSE)
    results$summaries['m1m2', 'clip.RMSE', fold] <- test.results$m1m2.RMSE.clip[trained$iter]
    results$summaries['m1m2', 'min.clip.RMSE.iter', fold] <- which.min(test.results$m1m2.RMSE.clip)
    results$summaries['m1m2', 'min.clip.RMSE', fold] <- min(test.results$m1m2.RMSE.clip)
    results$summaries['m1m2', 'exp.var', fold] <- test.results$m1m2.exp.var[trained$iter]
    results$summaries['m1m2', 'max.exp.var.iter', fold] <- which.max(test.results$m1m2.exp.var)
    results$summaries['m1m2', 'max.exp.var', fold] <- max(test.results$m1m2.exp.var)
    if(length(test.results$m1m2.p.cor)) {
      results$summaries['m1m2', 'p.cor', fold] <- test.results$m1m2.p.cor[trained$iter]
      results$summaries['m1m2', 'max.p.cor.iter', fold] <- which.max(test.results$m1m2.p.cor)
      results$summaries['m1m2', 'max.p.cor', fold] <- max(test.results$m1m2.p.cor)
      results$summaries['m1m2', 's.cor', fold] <- test.results$m1m2.s.cor[trained$iter]
      results$summaries['m1m2', 'max.s.cor.iter', fold] <- which.max(test.results$m1m2.s.cor)
      results$summaries['m1m2', 'max.s.cor', fold] <- max(test.results$m1m2.s.cor)
    }
  }
  
  if('m1m3.RMSE' %in% names(test.results)) {
    results$summaries['m1m3', 'RMSE', fold] <- test.results$m1m3.RMSE[trained$iter]
    results$summaries['m1m3', 'min.RMSE.iter', fold] <- which.min(test.results$m1m3.RMSE)
    results$summaries['m1m3', 'min.RMSE', fold] <- min(test.results$m1m3.RMSE)
    results$summaries['m1m3', 'clip.RMSE', fold] <- test.results$m1m3.RMSE.clip[trained$iter]
    results$summaries['m1m3', 'min.clip.RMSE.iter', fold] <- which.min(test.results$m1m3.RMSE.clip)
    results$summaries['m1m3', 'min.clip.RMSE', fold] <- min(test.results$m1m3.RMSE.clip)
    results$summaries['m1m3', 'exp.var', fold] <- test.results$m1m3.exp.var[trained$iter]
    results$summaries['m1m3', 'max.exp.var.iter', fold] <- which.max(test.results$m1m3.exp.var)
    results$summaries['m1m3', 'max.exp.var', fold] <- max(test.results$m1m3.exp.var)
    if(length(test.results$m1m3.p.cor)) {
      results$summaries['m1m3', 'p.cor', fold] <- test.results$m1m3.p.cor[trained$iter]
      results$summaries['m1m3', 'max.p.cor.iter', fold] <- which.max(test.results$m1m3.p.cor)
      results$summaries['m1m3', 'max.p.cor', fold] <- max(test.results$m1m3.p.cor)
      results$summaries['m1m3', 's.cor', fold] <- test.results$m1m3.s.cor[trained$iter]
      results$summaries['m1m3', 'max.s.cor.iter', fold] <- which.max(test.results$m1m3.s.cor)
      results$summaries['m1m3', 'max.s.cor', fold] <- max(test.results$m1m3.s.cor)
    }
  }
  
  if('m2m3.RMSE' %in% names(test.results)) {
    results$summaries['m2m3', 'RMSE', fold] <- test.results$m2m3.RMSE[trained$iter]
    results$summaries['m2m3', 'min.RMSE.iter', fold] <- which.min(test.results$m2m3.RMSE)
    results$summaries['m2m3', 'min.RMSE', fold] <- min(test.results$m2m3.RMSE)
    results$summaries['m2m3', 'clip.RMSE', fold] <- test.results$m2m3.RMSE.clip[trained$iter]
    results$summaries['m2m3', 'min.clip.RMSE.iter', fold] <- which.min(test.results$m2m3.RMSE.clip)
    results$summaries['m2m3', 'min.clip.RMSE', fold] <- min(test.results$m2m3.RMSE.clip)
    results$summaries['m2m3', 'exp.var', fold] <- test.results$m2m3.exp.var[trained$iter]
    results$summaries['m2m3', 'max.exp.var.iter', fold] <- which.max(test.results$m2m3.exp.var)
    results$summaries['m2m3', 'max.exp.var', fold] <- max(test.results$m2m3.exp.var)
    if(length(test.results$m2m3.p.cor)) {
      results$summaries['m2m3', 'p.cor', fold] <- test.results$m2m3.p.cor[trained$iter]
      results$summaries['m2m3', 'max.p.cor.iter', fold] <- which.max(test.results$m2m3.p.cor)
      results$summaries['m2m3', 'max.p.cor', fold] <- max(test.results$m2m3.p.cor)
      results$summaries['m2m3', 's.cor', fold] <- test.results$m2m3.s.cor[trained$iter]
      results$summaries['m2m3', 'max.s.cor.iter', fold] <- which.max(test.results$m2m3.s.cor)
      results$summaries['m2m3', 'max.s.cor', fold] <- max(test.results$m2m3.s.cor)
    }
  }
  
  if('m1m2m3.RMSE' %in% names(test.results)) {
    results$summaries['m1m2m3', 'RMSE', fold] <- test.results$m1m2m3.RMSE[trained$iter]
    results$summaries['m1m2m3', 'min.RMSE.iter', fold] <- which.min(test.results$m1m2m3.RMSE)
    results$summaries['m1m2m3', 'min.RMSE', fold] <- min(test.results$m1m2m3.RMSE)
    results$summaries['m1m2m3', 'clip.RMSE', fold] <- test.results$m1m2m3.RMSE.clip[trained$iter]
    results$summaries['m1m2m3', 'min.clip.RMSE.iter', fold] <- which.min(test.results$m1m2m3.RMSE.clip)
    results$summaries['m1m2m3', 'min.clip.RMSE', fold] <- min(test.results$m1m2m3.RMSE.clip)
    results$summaries['m1m2m3', 'exp.var', fold] <- test.results$m1m2m3.exp.var[trained$iter]
    results$summaries['m1m2m3', 'max.exp.var.iter', fold] <- which.max(test.results$m1m2m3.exp.var)
    results$summaries['m1m2m3', 'max.exp.var', fold] <- max(test.results$m1m2m3.exp.var)
    if(length(test.results$m1m2m3.p.cor)) {
      results$summaries['m1m2m3', 'p.cor', fold] <- test.results$m1m2m3.p.cor[trained$iter]
      results$summaries['m1m2m3', 'max.p.cor.iter', fold] <- which.max(test.results$m1m2m3.p.cor)
      results$summaries['m1m2m3', 'max.p.cor', fold] <- max(test.results$m1m2m3.p.cor)
      results$summaries['m1m2m3', 's.cor', fold] <- test.results$m1m2m3.s.cor[trained$iter]
      results$summaries['m1m2m3', 'max.s.cor.iter', fold] <- which.max(test.results$m1m2m3.s.cor)
      results$summaries['m1m2m3', 'max.s.cor', fold] <- max(test.results$m1m2m3.s.cor)
    }
  }
  
  results$sparsity['m1', fold] <- sum(1/(trained$mode1.lambda.shape * 
                                           trained$mode1.lambda.scale) > m1.rem.cut)/dim(train.data$mode1.X)[2]
  results$sparsity['m2', fold] <- sum(1/(trained$mode2.lambda.shape * 
                                           trained$mode2.lambda.scale) > m2.rem.cut)/dim(train.data$mode2.X)[2]
  results$sparsity['m3', fold] <- sum(1/(trained$mode3.lambda.shape * 
                                           trained$mode3.lambda.scale) > m3.rem.cut)/dim(train.data$mode2.X)[3]
  results$sparsity['core', fold] <- sum(1/(trained$core.lambda.shape * 
                                             trained$core.lambda.scale) > core.rem.cut)/prod(dim(trained$core.mean))
  
  return(results)
}

########### MAIN ###############

# Determine the number of runs with this prefix
n.files <- length(list.files(path = dirname(run_prefix),
                  pattern = paste0(basename(run_prefix), '.[0-9]+.out')))

f1 <- paste0(run_prefix, '.0/image.Rdata')

iters <- get_iters(f1)

training <- array(NA, dim=c(38, n.files, iters),
                  dimnames=list(c('lower.bnd', 'H.RMSE', 'A.RMSE', 'warm.RMSE', 
                                  'm1.RMSE', 'm2.RMSE', 'm3.RMSE',
                                  'm1m2.RMSE', 'm1m3.RMSE', 'm2m3.RMSE',
                                  'm1m2m3.RMSE',
                                  'exp.var', 'warm.exp.var', 
                                  'm1.exp.var', 'm2.exp.var', 'm3.exp.var',
                                  'm1m2.exp.var', 'm1m3.exp.var', 'm2m3.exp.var',
                                  'm1m2m3.exp.var',
                                  'p.cor', 'warm.p.cor', 
                                  'm1.p.cor', 'm2.p.cor', 'm3.p.cor',
                                  'm1m2.p.cor', 'm1m3.p.cor', 'm2m3.p.cor',
                                  'm1m2m3.p.cor',
                                  's.cor', 'warm.s.cor', 
                                  'm1.s.cor', 'm2.s.cor', 'm3.s.cor',
                                  'm1m2.s.cor', 'm1m3.s.cor', 'm2m3.s.cor', 
                                  'm1m2m3.s.cor'),
                                paste0('fold.', 1:n.files), 1:iters))

summaries <- array(NA, dim=c(11, 15, n.files),
                   dimnames=list(c('A' ,'H', 'train', 'warm',
                                   'm1', 'm2', 'm3', 
                                   'm1m2', 'm1m3', 'm2m3', 'm1m2m3'),
                                 c('RMSE', 'min.RMSE', 'min.RMSE.iter',
                                   'clip.RMSE', 'min.clip.RMSE', 'min.clip.RMSE.iter', 
                                   'exp.var', 'max.exp.var.iter', 'max.exp.var',
                                   'p.cor', 'max.p.cor.iter', 'max.p.cor',
                                   's.cor', 'max.s.cor.iter', 'max.s.cor'),
                                 paste0('fold.', 1:n.files)))

mean <- array(NA, dim=c(4, 21, n.files), 
              dimnames=list(c('RMSE', 'exp.var', 'p.cor', 's.cor'),
                            c('train.m1', 'train.m2', 'train.m3',
                              'train.m1m2', 'train.m1m3', 'train.m2m3', 
                              'train.m1m2m3',
                              'warm.m1', 'warm.m2', 'warm.m3', 
                              'warm.m1m2', 'warm.m1m3', 'warm.m2m3', 
                              'warm.m1m2m3',
                              'm1', 'm2', 'm3', 'm1m2', 'm1m3', 'm2m3', 'm1m2m3'),
                            paste0('fold.', 1:n.files)))

lower.bnd <- matrix(NA, 4, n.files, dimnames=
                      list(c('final', 'max', 'which.max', 'not.mono'),
                           paste0('fold.', 1:n.files)))

sparsity <- matrix(NA, 4, n.files,
                   dimnames=list(c('m1', 'm2', 'm3', 'core'),
                                 paste0('fold.', 1:n.files)))

results <- list(training=training, summaries=summaries,
                mean=mean, lower.bnd=lower.bnd, sparsity=sparsity)

rm(training, summaries, mean, lower.bnd, sparsity)

for(fold in 1:n.files) {
  # Load in the run data
  f <- paste0(run_prefix, '.', (fold-1), '/image.Rdata')
  results <- loadData(f, results)
}

# Make data frame counting how many folds peform better than the mean
better <- matrix(NA, 11, 10, dimnames=list(dimnames(results$summaries)[[1]],
                dimnames(results$summaries)[[2]][!grepl('iter', dimnames(results$summaries)[[2]])]))

for(type in c('A', 'H', 'train')) for(resp in c('RMSE', 'min.RMSE', 'clip.RMSE', 'min.clip.RMSE')) 
  better[type, resp] <- sum(results$summaries[type, resp,] < results$mean['RMSE', 'train.m1',])

for(resp in c('RMSE', 'min.RMSE', 'clip.RMSE', 'min.clip.RMSE')) 
  better['warm', resp] <- sum(results$summaries['warm', resp,] < results$mean['RMSE', 'warm.m1',])

for(type in c('m1', 'm2', 'm3', 'm1m2', 'm1m3', 'm2m3', 'm1m2m3')) 
  for (resp in c('RMSE', 'min.RMSE', 'clip.RMSE', 'min.clip.RMSE')) 
    better[type, resp] <- sum(results$summaries[type, resp,] < results$mean['RMSE', type,])

for(type in c('A', 'H', 'train')) for(resp in c('exp.var', 'max.exp.var')) 
  better[type, resp] <- sum(results$summaries[type, resp,] > results$mean['exp.var', 'train.m1',])
for(resp in c('exp.var', 'max.exp.var')) 
  better['warm', resp] <- sum(results$summaries['warm', resp,] > results$mean['exp.var', 'warm.m1',])
for(type in c('m1', 'm2', 'm3', 'm1m2', 'm1m3', 'm2m3', 'm1m2m3')) for (resp in c('exp.var', 'max.exp.var')) 
  better[type, resp] <- sum(results$summaries[type, resp,] > results$mean['exp.var', type,])
for(type in c('A', 'H', 'train')) for(resp in c('p.cor', 'max.p.cor')) 
  better[type, resp] <- sum(results$summaries[type, resp,] > results$mean['p.cor', 'train.m1',])
for(resp in c('p.cor', 'max.p.cor')) 
  better['warm', resp] <- sum(results$summaries['warm', resp,] > results$mean['p.cor', 'warm.m1',])
for(type in c('m1', 'm2', 'm3', 'm1m2', 'm1m3', 'm2m3', 'm1m2m3')) for (resp in c('p.cor', 'max.p.cor')) 
  better[type, resp] <- sum(results$summaries[type, resp,] > results$mean['p.cor', type,])
for(type in c('A', 'H', 'train')) for(resp in c('s.cor', 'max.s.cor')) 
  better[type, resp] <- sum(results$summaries[type, resp,] > results$mean['s.cor', 'train.m1',])
for(resp in c('s.cor', 'max.s.cor')) 
  better['warm', resp] <- sum(results$summaries['warm', resp,] > results$mean['s.cor', 'warm.m1',])
for(type in c('m1', 'm2', 'm3', 'm1m2', 'm1m3', 'm2m3', 'm1m2m3')) for (resp in c('s.cor', 'max.s.cor')) 
  better[type, resp] <- sum(results$summaries[type, resp,] > results$mean['s.cor', type,])

# Read log files to get run time if the runs finished
if(length(system2('grep', c('"Job terminated"', paste0(run_prefix, '.*.log')), stdout=T)>0)) {
  months <- c(31,28,31,30,31,30,31,31,30,31,30,31)
  month.days <- cumsum(months)
  log.starts <- system2('grep', c('"Job submitted"', paste0(run_prefix, '.*.log')), stdout=T)
  log.ends <- system2('grep', c('"Job terminated"', paste0(run_prefix, '.*.log')), stdout=T)
  month.start <- as.numeric(sapply(strsplit(sapply(strsplit(log.starts, split=' '), 
                                                   '[', 3), split='/'), '[', 1))
  month.end <- as.numeric(sapply(strsplit(sapply(strsplit(log.ends, split=' '), 
                                                 '[', 3), split='/'), '[', 1))
  day.start <- month.days[month.start] + 
    as.numeric(sapply(strsplit(sapply(strsplit(log.starts, split=' '), '[', 3), split='/'), '[', 2))
  day.end <- month.days[month.end] + 
    as.numeric(sapply(strsplit(sapply(strsplit(log.ends, split=' '), '[', 3), split='/'), '[', 2))
  days <- day.end - day.start
  hour.start <- as.numeric(sapply(strsplit(sapply(strsplit(log.starts, split=' '), 
                                                  '[', 4), split=':'), '[', 1))
  hour.end <- as.numeric(sapply(strsplit(sapply(strsplit(log.ends, split=' '), 
                                                '[', 4), split=':'), '[', 1))
  min.start <- as.numeric(sapply(strsplit(sapply(strsplit(log.starts, split=' '), 
                                                 '[', 4), split=':'), '[', 2))
  min.end <- as.numeric(sapply(strsplit(sapply(strsplit(log.ends, split=' '), 
                                               '[', 4), split=':'), '[', 2))
  
  hours <- days * 24 + (hour.end - hour.start) + (min.end - min.start)/60
  rm(log.starts, log.ends)
  print("Run time statistics (hours):")
  summary(hours)
}

# Save all the data
save.image(paste0(run_prefix, '_summary.Rdata'))

print("######################################################")
print('## Sparsity ##')
print(apply(results$sparsity, 1, summary))
print("######################################################")
print("Lower bounds")
print(apply(results$lower.bnd, 1, summary, na.rm=T))
print("Standard deviations")
print(apply(results$lower.bnd, 1, sd, na.rm=T))
print("######################################################")
print('## Means ##')
print(apply(results$summaries, c(1,2), mean, na.rm=T))
print('## Standard deviations ##')
print(apply(results$summaries, c(1,2), sd, na.rm=T))
print("######################################################")
print('## Better than mean ##')
print(better)

# Plotting
#####################################################################

cols <- rainbow(n.files)

pdf(file=paste0(run_prefix, '_lower_bounds.pdf'), height=14)
par(mfrow=c(2,1))
plot(results$training['lower.bnd',1,], type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['lower.bnd',,]))
for(i in 2:n.files) points(results$training['lower.bnd',i,], type='l', lwd=2, col=cols[i])
plot(results$training['lower.bnd',1,], type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['lower.bnd',,5:iters]))
for(i in 2:n.files) points(results$training['lower.bnd',i,], type='l', lwd=2, col=cols[i])
dev.off()

######### Plot RMSEs ##############

pdf(file=paste0(run_prefix, '_training_RMSE.pdf'))
plot(results$training['A.RMSE',1,] - results$mean['RMSE', 'train.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['A.RMSE',,] - results$mean['RMSE', 'train.m1',]),
     main="Training RMSEs relative \n to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['A.RMSE',i,] - results$mean['RMSE', 'train.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_warm_RMSEs.pdf'))
plot(results$training['warm.RMSE',1,] - results$mean['RMSE', 'warm.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['warm.RMSE',,] - results$mean['RMSE', 'warm.m1',]),
     main="Warm RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['warm.RMSE',i,] - results$mean['RMSE', 'warm.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1_RMSEs.pdf'))
plot(results$training['m1.RMSE',1,] - results$mean['RMSE', 'm1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1.RMSE',,] - results$mean['RMSE', 'm1',]),
     main="Cold mode 1 RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['m1.RMSE',i,] - results$mean['RMSE', 'm1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2_RMSEs.pdf'))
plot(results$training['m2.RMSE',1,] - results$mean['RMSE', 'm2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2.RMSE',,] - results$mean['RMSE', 'm2',]),
     main="Cold mode 2 RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['m2.RMSE',i,] - results$mean['RMSE', 'm2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m3_RMSEs.pdf'))
plot(results$training['m3.RMSE',1,] - results$mean['RMSE', 'm3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m3.RMSE',,] - results$mean['RMSE', 'm3',]),
     main="Cold mode 3 RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['m3.RMSE',i,] - results$mean['RMSE', 'm3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2_RMSEs.pdf'))
plot(results$training['m1m2.RMSE',1,] - results$mean['RMSE', 'm1m2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2.RMSE',,] - results$mean['RMSE', 'm1m2',], na.rm=T),
     main="Cold mode 1/mode 2 RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['m1m2.RMSE',i,] - results$mean['RMSE', 'm1m2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m3_RMSEs.pdf'))
plot(results$training['m1m3.RMSE',1,] - results$mean['RMSE', 'm1m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m3.RMSE',,] - results$mean['RMSE', 'm1m3',], na.rm=T),
     main="Cold mode 1/mode 3 RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['m1m3.RMSE',i,] - results$mean['RMSE', 'm1m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2m3_RMSEs.pdf'))
plot(results$training['m2m3.RMSE',1,] - results$mean['RMSE', 'm2m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2m3.RMSE',,] - results$mean['RMSE', 'm2m3',], na.rm=T),
     main="Cold mode 2/mode 3 RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['m2m3.RMSE',i,] - results$mean['RMSE', 'm2m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2m3_RMSEs.pdf'))
plot(results$training['m1m2m3.RMSE',1,] - results$mean['RMSE', 'm1m2m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2m3.RMSE',,] - results$mean['RMSE', 'm1m2m3',], na.rm=T),
     main="Cold mode 1,2 & 3 RMSEs relative to predicting mean response",
     xlab="Iteration", ylab="Relative RMSE")
for(i in 2:n.files)
  points(results$training['m1m2m3.RMSE',i,] - results$mean['RMSE', 'm1m2m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

########## Plot explained variances ################

pdf(file=paste0(run_prefix, '_training_exp_var.pdf'))
plot(results$training['exp.var',1,] - results$mean['exp.var', 'train.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['exp.var',,] - results$mean['exp.var', 'train.m1',]),
     main="Training explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['exp.var',i,] - results$mean['exp.var', 'train.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_warm_exp_var.pdf'))
plot(results$training['warm.exp.var',1,] - results$mean['exp.var', 'warm.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['warm.exp.var',,] - results$mean['exp.var', 'warm.m1',]),
     main="Warm explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['warm.exp.var',i,] - results$mean['exp.var', 'warm.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1_exp_var.pdf'))
plot(results$training['m1.exp.var',1,] - results$mean['exp.var', 'm1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1.exp.var',,] - results$mean['exp.var', 'm1',]),
     main="Cold mode 1 explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['m1.exp.var',i,] - results$mean['exp.var', 'm1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2_exp_var.pdf'))
plot(results$training['m2.exp.var',1,] - results$mean['exp.var', 'm2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2.exp.var',,] - results$mean['exp.var', 'm2',]),
     main="Cold mode 2 explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['m2.exp.var',i,] - results$mean['exp.var', 'm2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m3_exp_var.pdf'))
plot(results$training['m3.exp.var',1,] - results$mean['exp.var', 'm3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m3.exp.var',,] - results$mean['exp.var', 'm3',]),
     main="Cold mode 3 explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['m3.exp.var',i,] - results$mean['exp.var', 'm3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2_exp_var.pdf'))
plot(results$training['m1m2.exp.var',1,] - results$mean['exp.var', 'm1m2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2.exp.var',,] - results$mean['exp.var', 'm1m2',], na.rm=T),
     main="Cold mode 1/mode 2 explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['m1m2.exp.var',i,] - results$mean['exp.var', 'm1m2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m3_exp_var.pdf'))
plot(results$training['m1m3.exp.var',1,] - results$mean['exp.var', 'm1m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m3.exp.var',,] - results$mean['exp.var', 'm1m3',], na.rm=T),
     main="Cold mode 1/mode 3 explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['m1m3.exp.var',i,] - results$mean['exp.var', 'm1m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2m3_exp_var.pdf'))
plot(results$training['m2m3.exp.var',1,] - results$mean['exp.var', 'm2m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2m3.exp.var',,] - results$mean['exp.var', 'm2m3',], na.rm=T),
     main="Cold mode 2/mode 3 explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['m2m3.exp.var',i,] - results$mean['exp.var', 'm2m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2m3_exp_var.pdf'))
plot(results$training['m1m2m3.exp.var',1,] - results$mean['exp.var', 'm1m2m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2m3.exp.var',,] - results$mean['exp.var', 'm1m2m3',], na.rm=T),
     main="Cold mode 1,2 & 3 explained variance \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative exp. var.")
for(i in 2:n.files)
  points(results$training['m1m2m3.exp.var',i,] - results$mean['exp.var', 'm1m2m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

########## Plot Pearson correlations ################

pdf(file=paste0(run_prefix, '_training_p_cor.pdf'))
plot(results$training['p.cor',1,] - results$mean['p.cor', 'train.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['p.cor',,] - results$mean['p.cor', 'train.m1',]),
     main="Training Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['p.cor',i,] - results$mean['p.cor', 'train.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_warm_p_cor.pdf'))
plot(results$training['warm.p.cor',1,] - results$mean['p.cor', 'warm.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['warm.p.cor',,] - results$mean['p.cor', 'warm.m1',]),
     main="Warm Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['warm.p.cor',i,] - results$mean['p.cor', 'warm.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1_p_cor.pdf'))
plot(results$training['m1.p.cor',1,] - results$mean['p.cor', 'm1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1.p.cor',,] - results$mean['p.cor', 'm1',]),
     main="Cold mode 1 Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['m1.p.cor',i,] - results$mean['p.cor', 'm1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2_p_cor.pdf'))
plot(results$training['m2.p.cor',1,] - results$mean['p.cor', 'm2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2.p.cor',,] - results$mean['p.cor', 'm2',]),
     main="Cold mode 2 Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['m2.p.cor',i,] - results$mean['p.cor', 'm2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m3_p_cor.pdf'))
plot(results$training['m3.p.cor',1,] - results$mean['p.cor', 'm3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m3.p.cor',,] - results$mean['p.cor', 'm3',]),
     main="Cold mode 3 Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['m3.p.cor',i,] - results$mean['p.cor', 'm3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2_p_cor.pdf'))
plot(results$training['m1m2.p.cor',1,] - results$mean['p.cor', 'm1m2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2.p.cor',,] - results$mean['p.cor', 'm1m2',], na.rm=T),
     main="Cold mode 1/mode 2 Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['m1m2.p.cor',i,] - results$mean['p.cor', 'm1m2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m3_p_cor.pdf'))
plot(results$training['m1m3.p.cor',1,] - results$mean['p.cor', 'm1m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m3.p.cor',,] - results$mean['p.cor', 'm1m3',], na.rm=T),
     main="Cold mode 1/mode 3 Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['m1m3.p.cor',i,] - results$mean['p.cor', 'm1m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2m3_p_cor.pdf'))
plot(results$training['m2m3.p.cor',1,] - results$mean['p.cor', 'm2m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2m3.p.cor',,] - results$mean['p.cor', 'm2m3',], na.rm=T),
     main="Cold mode 1/mode 3 Pearson correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Pearson correlation")
for(i in 2:n.files)
  points(results$training['m2m3.p.cor',i,] - results$mean['p.cor', 'm2m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2m3_p_cor.pdf'))
plot(results$training['m1m2m3.p.cor',1,] , 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2m3.p.cor',,], na.rm=T),
     main="Cold mode 1,2 & 3 Pearson correlation",
     xlab="Iteration", ylab="Pearson correlation")
for(i in 2:n.files)
  points(results$training['m1m2m3.p.cor',i,],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

########## Plot Spearman correlations ################

pdf(file=paste0(run_prefix, '_training_s_cor.pdf'))
plot(results$training['s.cor',1,] - results$mean['s.cor', 'train.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['s.cor',,] - results$mean['s.cor', 'train.m1',]),
     main="Training Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['s.cor',i,] - results$mean['s.cor', 'train.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_warm_s_cor.pdf'))
plot(results$training['warm.s.cor',1,] - results$mean['s.cor', 'warm.m1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['warm.s.cor',,] - results$mean['s.cor', 'warm.m1',]),
     main="Warm Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['warm.s.cor',i,] - results$mean['s.cor', 'warm.m1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1_s_cor.pdf'))
plot(results$training['m1.s.cor',1,] - results$mean['s.cor', 'm1',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1.s.cor',,] - results$mean['s.cor', 'm1',]),
     main="Cold mode 1 Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['m1.s.cor',i,] - results$mean['s.cor', 'm1',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2_s_cor.pdf'))
plot(results$training['m2.s.cor',1,] - results$mean['s.cor', 'm2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2.s.cor',,] - results$mean['s.cor', 'm2',]),
     main="Cold mode 2 Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['m2.s.cor',i,] - results$mean['s.cor', 'm2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m3_s_cor.pdf'))
plot(results$training['m3.s.cor',1,] - results$mean['s.cor', 'm3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m3.s.cor',,] - results$mean['s.cor', 'm3',]),
     main="Cold mode 3 Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['m3.s.cor',i,] - results$mean['s.cor', 'm3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2_s_cor.pdf'))
plot(results$training['m1m2.s.cor',1,] - results$mean['s.cor', 'm1m2',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2.s.cor',,] - results$mean['s.cor', 'm1m2',], na.rm=T),
     main="Cold mode 1/mode 2 Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['m1m2.s.cor',i,] - results$mean['s.cor', 'm1m2',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m3_s_cor.pdf'))
plot(results$training['m1m3.s.cor',1,] - results$mean['s.cor', 'm1m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m3.s.cor',,] - results$mean['s.cor', 'm1m3',], na.rm=T),
     main="Cold mode 1/mode 3 Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['m1m3.s.cor',i,] - results$mean['s.cor', 'm1m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m2m3_s_cor.pdf'))
plot(results$training['m2m3.s.cor',1,] - results$mean['s.cor', 'm2m3',1], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m2m3.s.cor',,] - results$mean['s.cor', 'm2m3',], na.rm=T),
     main="Cold mode 2/mode 3 Spearman correlation \n relative to predicting mean response",
     xlab="Iteration", ylab="Relative Spearman correlation")
for(i in 2:n.files)
  points(results$training['m2m3.s.cor',i,] - results$mean['s.cor', 'm2m3',i],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()

pdf(file=paste0(run_prefix, '_cold_m1m2m3_s_cor.pdf'))
plot(results$training['m1m2m3.s.cor',1,], 
     type='l', lwd=2, col=cols[1], 
     ylim=range(results$training['m1m2m3.s.cor',,], na.rm=T),
     main="Cold mode 1,2 & 3 Spearman correlation",
     xlab="Iteration", ylab="Spearman correlation")
for(i in 2:n.files)
  points(results$training['m1m2m3.s.cor',i,],
         type='l', lwd=2, col=cols[i])
abline(h=0, lty=2, lwd=2)
dev.off()
nathanlazar/BaTFLED3D documentation built on May 23, 2019, 12:19 p.m.