#' Evaluates a trained model on fasta, fastq or rds files

#'

#' Returns evaluation metric like confusion matrix, loss, AUC, AUPRC, MAE, MSE (depending on output layer).

#'

#' @inheritParams generator_fasta_lm

#' @inheritParams generator_fasta_label_folder

#' @inheritParams generator_fasta_label_header_csv

#' @param path_input Input directory where fasta, fastq or rds files are located.

#' @param model A keras model.

#' @param batch_size Number of samples per batch.

#' @param step How often to take a sample.

#' @param vocabulary Vector of allowed characters. Character outside vocabulary get encoded as specified in ambiguous_nuc.

#' @param vocabulary_label List of labels for targets of each output layer.

#' @param number_batches How many batches to evaluate.

#' @param format File format, `"fasta"`, `"fastq"` or `"rds"`.

#' @param mode Either `"lm"` for language model or `"label_header"`, `"label_csv"` or `"label_folder"` for label classification.

#' @param verbose Boolean.

#' @param target_middle Whether model is language model with separate input layers. 

#' @param evaluate_all_files Boolean, if `TRUE` will iterate over all files in \code{path_input} once. \code{number_batches} will be overwritten.

#' @param auc Whether to include AUC metric. If output layer activation is `"softmax"`, only possible for 2 targets. Computes the average if output layer has sigmoid

#' activation and multiple targets.

#' @param auprc Whether to include AUPRC metric. If output layer activation is `"softmax"`, only possible for 2 targets. Computes the average if output layer has sigmoid

#' activation and multiple targets.

#' @param path_pred_list Path to store list of predictions (output of output layers) and corresponding true labels as rds file. 

#' @param exact_num_samples Exact number of samples to evaluate. If you want to evaluate a number of samples not divisible by batch_size. Useful if you want

#' to evaluate a data set exactly ones and know the number of samples already. Should be a vector if `mode = "label_folder"` (with same length as `vocabulary_label`)

#' and else an integer.

#' @param activations List containing output formats for output layers (`softmax, sigmoid` or `linear`). If `NULL`, will be estimated from model.   

#' @param include_seq Whether to store input. Only applies if `path_pred_list` is not `NULL`.

#' @param ... Further generator options. See \code{\link{get_generator}}.

#' @examplesIf reticulate::py_module_available("tensorflow")

#' # create dummy data

#' path_input <- tempfile()

#' dir.create(path_input)

#' create_dummy_data(file_path = path_input,

#'                   num_files = 3,

#'                   seq_length = 11, 

#'                   num_seq = 5,

#'                   vocabulary = c("a", "c", "g", "t"))

#' # create model

#' model <- create_model_lstm_cnn(layer_lstm = 8, layer_dense = 4, maxlen = 10, verbose = FALSE)

#' # evaluate

#' evaluate_model(path_input = path_input,

#'   model = model,

#'   step = 11,

#'   vocabulary = c("a", "c", "g", "t"),

#'   vocabulary_label = list(c("a", "c", "g", "t")),

#'   mode = "lm",

#'   output_format = "target_right",

#'   evaluate_all_files = TRUE,

#'   verbose = FALSE)

#'   

#' @returns A list of evaluation results. Each list element corresponds to an output layer of the model.   

#' @export

evaluate_model <- function(path_input,

                           model = NULL,

                           batch_size = 100,

                           step = 1,

                           padding = FALSE,

                           vocabulary = c("a", "c", "g", "t"),

                           vocabulary_label = list(c("a", "c", "g", "t")),

                           number_batches = 10,

                           format = "fasta",

                           target_middle = FALSE,

                           mode = "lm",

                           output_format = "target_right",

                           ambiguous_nuc = "zero",

                           evaluate_all_files = FALSE,

                           verbose = TRUE,

                           max_iter = 20000,

                           target_from_csv = NULL,

                           max_samples = NULL,

                           proportion_per_seq = NULL,

                           concat_seq = NULL,

                           seed = 1234,

                           auc = FALSE,

                           auprc = FALSE,

                           path_pred_list = NULL,

                           exact_num_samples = NULL,

                           activations = NULL,

                           shuffle_file_order = FALSE,

                           include_seq = FALSE,

                           ...) {

  set.seed(seed)

  path_model <- NULL

  stopifnot(mode %in% c("lm", "label_header", "label_folder", "label_csv", "lm_rds", "label_rds"))

  stopifnot(format %in% c("fasta", "fastq", "rds"))

  stopifnot(is.null(proportion_per_seq) || proportion_per_seq <= 1)

  if (!is.null(exact_num_samples) & evaluate_all_files) {

    warning(paste("Will evaluate number of samples as specified in exact_num_samples argument. Setting evaluate_all_files to FALSE."))

    evaluate_all_files <- FALSE

  }

  eval_exact_num_samples <- !is.null(exact_num_samples) | evaluate_all_files

  if (is.null(activations)) activations <- get_output_activations(model)

  if (is.null(path_pred_list) & include_seq) {

    stop("Can only store input, if path_pred_list is specified.")

  }

  if (is.null(vocabulary_label)) vocabulary_label <- list(vocabulary)

  if (!is.list(vocabulary_label)) vocabulary_label <- list(vocabulary_label)

  if (mode == "label_folder") {

    number_batches <- rep(ceiling(number_batches/length(path_input)), length(path_input))

  }

  num_classes <- ifelse(mode == "label_folder", length(path_input), 1)

  num_out_layers <- length(activations)

  # extract maxlen from model

  num_in_layers <- length(model$inputs)

  if (num_in_layers == 1) {

    maxlen <- model$input$shape[[2]]

  } else {

    if (!target_middle) {

      maxlen <- model$input[[num_in_layers]]$shape[[2]]

    } else {

      maxlen <- model$input[[num_in_layers - 1]]$shape[[2]] + model$input[[num_in_layers]]$shape[[2]]

    }

  }

  if (evaluate_all_files & (format %in% c("fasta", "fastq"))) {

    number_batches <- NULL

    num_samples <- rep(0, length(path_input))

    for (i in 1:num_classes) {

      if (mode == "label_folder") {

        files <- list_fasta_files(path_input[[i]], format = format, file_filter = NULL)

      } else {

        files <- list_fasta_files(path_input, format = format, file_filter = NULL)

      }

      # remove files not in csv table

      if (mode == "label_csv") {

        csv_file <- utils::read.csv2(target_from_csv, header = TRUE, stringsAsFactors = FALSE)

        if (dim(csv_file)[2] == 1) {

          csv_file <- utils::read.csv(target_from_csv, header = TRUE, stringsAsFactors = FALSE)

        }

        index <- basename(files) %in% csv_file$file

        files <- files[index]

        if (length(files) == 0) {

          stop("No files from path_input have label in target_from_csv file.")

        }

      }

      for (file in files) {

        if (format == "fasta") {

          fasta_file <- microseq::readFasta(file)

        } else {

          fasta_file <- microseq::readFastq(file)

        }

        # remove entries with wrong header

        if (mode == "label_header") {

          index <- fasta_file$Header %in% vocabulary_label

          fasta_file <- fasta_file[index, ]

        }

        seq_vector <- fasta_file$Sequence

        if (!is.null(concat_seq)) {

          seq_vector <- paste(seq_vector, collapse = concat_seq)

        }

        if (!is.null(proportion_per_seq)) {

          fasta_width <- nchar(seq_vector)

          sample_range <- floor(fasta_width - (proportion_per_seq * fasta_width))

          start <- mapply(sample_range, FUN = sample, size = 1)

          perc_length <- floor(fasta_width * proportion_per_seq)

          stop <- start + perc_length

          seq_vector <- mapply(seq_vector, FUN = substr, start = start, stop = stop)

        }

        if (mode == "lm") {

          if (!padding) {

            seq_vector <- seq_vector[nchar(seq_vector) >= (maxlen + 1)]

          } else {

            length_vector <- nchar(seq_vector)

            short_seq_index <- which(length_vector < (maxlen + 1))

            for (ssi in short_seq_index) {

              seq_vector[ssi] <- paste0(paste(rep("0", (maxlen + 1) - length_vector[ssi]), collapse = ""), seq_vector[ssi])

            }

          }

        } else {

          if (!padding) {

            seq_vector <- seq_vector[nchar(seq_vector) >= (maxlen)]

          } else {

            length_vector <- nchar(seq_vector)

            short_seq_index <- which(length_vector < (maxlen))

            for (ssi in short_seq_index) {

              seq_vector[ssi] <- paste0(paste(rep("0", (maxlen) - length_vector[ssi]), collapse = ""), seq_vector[ssi])

            }

          }

        }

        if (length(seq_vector) == 0) next

        new_samples <- get_start_ind(seq_vector = seq_vector,

                                     length_vector = nchar(seq_vector),

                                     maxlen = maxlen,

                                     step = step,

                                     train_mode = ifelse(mode == "lm", "lm", "label"),

                                     discard_amb_nuc = ifelse(ambiguous_nuc == "discard", TRUE, FALSE),

                                     vocabulary = vocabulary

        ) %>% length()

        if (is.null(max_samples)) {

          num_samples[i] <- num_samples[i] + new_samples

        } else {

          num_samples[i] <- num_samples[i] + min(new_samples, max_samples)

        }

      }

      number_batches[i] <- ceiling(num_samples[i]/batch_size)

    }

    if (mode == "label_folder") {

      message_string <- paste0("Evaluate ", num_samples, " samples for class ", vocabulary_label[[1]], ".\n")

    } else {

      message_string <- paste0("Evaluate ", sum(num_samples), " samples.")

    }

    message(message_string)

  }

  if (evaluate_all_files & format == "rds") {

    rds_files <- list_fasta_files(path_corpus = path_input,

                                  format = "rds",

                                  file_filter = NULL)

    num_samples <- 0

    for (file in rds_files) {

      rds_file <- readRDS(file)

      x <- rds_file[[1]]

      while (is.list(x)) {

        x <- x[[1]]

      }

      num_samples <- dim(x)[1] + num_samples

    }

    number_batches <- ceiling(num_samples/batch_size)

    message_string <- paste0("Evaluate ", num_samples, " samples.")

    message(message_string)

  }

  if (!is.null(exact_num_samples)) {

    num_samples <- exact_num_samples

    number_batches <- ceiling(num_samples/batch_size)

  }

  overall_num_batches <- sum(number_batches)

  if (mode == "lm") {

    gen <- generator_fasta_lm(path_corpus = path_input,

                              format = format,

                              batch_size = batch_size,

                              maxlen = maxlen,

                              max_iter = max_iter,

                              vocabulary = vocabulary,

                              verbose = FALSE,

                              shuffle_file_order = shuffle_file_order,

                              step = step,

                              concat_seq = concat_seq,

                              padding = padding,

                              shuffle_input = FALSE,

                              reverse_complement = FALSE,

                              output_format = output_format,

                              ambiguous_nuc = ambiguous_nuc,

                              proportion_per_seq = proportion_per_seq,

                              max_samples = max_samples,

                              seed = seed,

                              ...)

  }

  if (mode == "label_header" | mode == "label_csv") {

    gen <- generator_fasta_label_header_csv(path_corpus = path_input,

                                            format = format,

                                            batch_size = batch_size,

                                            maxlen = maxlen,

                                            max_iter = max_iter,

                                            vocabulary = vocabulary,

                                            verbose = FALSE,

                                            shuffle_file_order = shuffle_file_order,

                                            step = step,

                                            padding = padding,

                                            shuffle_input = FALSE,

                                            concat_seq = concat_seq,

                                            vocabulary_label = vocabulary_label[[1]],

                                            reverse_complement = FALSE,

                                            ambiguous_nuc = ambiguous_nuc,

                                            target_from_csv = target_from_csv,

                                            proportion_per_seq = proportion_per_seq,

                                            max_samples = max_samples,

                                            seed = seed, ...)

  }

  if (mode == "label_rds" | mode == "lm_rds") {

    gen <- generator_rds(rds_folder = path_input, batch_size = batch_size, path_file_log = NULL, ...)

  }

  batch_index <- 1

  start_time <- Sys.time()

  ten_percent_steps <- seq(overall_num_batches/10, overall_num_batches, length.out = 10)

  percentage_index <- 1

  count <- 1

  y_conf_list <- vector("list", overall_num_batches)

  y_list <- vector("list", overall_num_batches)

  if (include_seq) {

    x_list <- vector("list", overall_num_batches)

  }

  for (k in 1:num_classes) {

    index <- NULL

    if (mode == "label_folder") {

      gen <- generator_fasta_label_folder(path_corpus = path_input[[k]],

                                          format = format,

                                          batch_size = batch_size,

                                          maxlen = maxlen,

                                          max_iter = max_iter,

                                          vocabulary = vocabulary,

                                          step = step,

                                          padding = padding,

                                          concat_seq = concat_seq,

                                          reverse_complement = FALSE,

                                          num_targets = length(path_input),

                                          ones_column = k,

                                          ambiguous_nuc = ambiguous_nuc,

                                          proportion_per_seq = proportion_per_seq,

                                          max_samples = max_samples,

                                          seed = seed, ...)

    }

    for (i in 1:number_batches[k]) {

      z <- gen()

      x <- z[[1]]

      y <- z[[2]]

      y_conf <- model(x)

      batch_index <- batch_index + 1

      # remove double predictions

      if (eval_exact_num_samples & (i == number_batches[k])) {

        double_index <- (i * batch_size) - num_samples[k]

        if (double_index > 0) {

          index <- 1:(nrow(y_conf) - double_index)

          if (is.list(y_conf)) {

            for (m in 1:length(y_conf)) {

              y_conf[[m]] <- y_conf[[m]][index, ]

              y[[m]] <- y[[m]][index, ]

            }

          } else {

            y_conf <- y_conf[index, ]

            y <- y[index, ]

          }

          # vector to matrix

          if (length(index) == 1) {

            if (is.list(y_conf)) {

              for (m in 1:length(y_conf)) {

                y_conf[[m]] <- array(as.array(y_conf[[m]]), dim = c(1, length(y_conf[[m]])))

                y[[m]] <- matrix(y[[m]], ncol = length(y[[m]]))

              }

            } else {

              y_conf <- array(as.array(y_conf), dim = c(1, length(y_conf)))

              y <- matrix(y, ncol = length(y))

            }

          }

        }

      }

      if (include_seq) {

        x_list[[count]] <- x

      }

      y_conf_list[[count]] <- y_conf

      if (batch_size == 1 | (!is.null(index) && length(index == 1))) {

        col_num <- ncol(y_conf)

        if (is.na(col_num)) col_num <- length(y_conf)

        y_list[[count]] <- matrix(y, ncol = col_num)

      } else {

        y_list[[count]] <- y

      }

      count <- count + 1

      if (verbose & (batch_index == 10)) {

        time_passed <- as.double(difftime(Sys.time(), start_time, units = "hours"))

        time_estimation <- (overall_num_batches/10) * time_passed

        cat("Evaluation will take approximately", round(time_estimation, 3), "hours. Starting time:", format(Sys.time(), "%F %R."), " \n")

      }

      if (verbose & (batch_index > ten_percent_steps[percentage_index]) & percentage_index < 10) {

        cat("Progress: ", percentage_index * 10 ,"% \n")

        time_passed <- as.double(difftime(Sys.time(), start_time, units = "hours"))

        cat("Time passed: ", round(time_passed, 3), "hours \n")

        percentage_index <- percentage_index + 1

      }

    }

  }

  if (verbose) {

    cat("Progress: 100 % \n")

    time_passed <- as.double(difftime(Sys.time(), start_time, units = "hours"))

    cat("Time passed: ", round(time_passed, 3), "hours \n")

  }

  y_conf_list <- reshape_y_list(y_conf_list, num_out_layers = num_out_layers, tf_format = TRUE)

  y_list <- reshape_y_list(y_list, num_out_layers = num_out_layers, tf_format = FALSE)

  if (!is.null(path_pred_list)) {

    if (include_seq) {

      if (is.list(x_list[[1]])) {

        num_layers <- length(x_list[[1]])

      } else {

        num_layers <- 1 

      }

      x_list <- reshape_y_list(x_list, num_out_layers = num_layers, tf_format = FALSE)

      saveRDS(list(pred = y_conf_list, true = y_list, x = x_list), path_pred_list)

    } else {

      saveRDS(list(pred = y_conf_list, true = y_list), path_pred_list)

    }

  } 

  eval_list <- list()

  for (i in 1:num_out_layers) {

    if (activations[i] == "softmax") {

      eval_list[[i]] <- evaluate_softmax(y = y_list[[i]], y_conf = y_conf_list[[i]],

                                         auc = auc, auprc = auprc,

                                         label_names = vocabulary_label[[i]])

    }

    if (activations[i] == "sigmoid") {

      eval_list[[i]] <- evaluate_sigmoid(y = y_list[[i]], y_conf = y_conf_list[[i]],

                                         auc = auc, auprc = auprc,

                                         label_names = vocabulary_label[[i]])

    }

    if (activations[i] == "linear") {

      eval_list[[i]] <- evaluate_linear(y_true = y_list[[i]], y_pred = y_conf_list[[i]], label_names = vocabulary_label[[i]])

    }

  }

  return(eval_list)

}

reshape_y_list <- function(y, num_out_layers, tf_format = TRUE) {

  if (num_out_layers > 1) {

    y <- do.call(c, y)

  }

  reshaped_list <- vector("list", num_out_layers)

  for (i in 1:num_out_layers) {

    index <- seq(i, length(y), by = num_out_layers)

    if (tf_format) {

      reshaped_list[[i]] <- y[index] %>%

        tensorflow::tf$concat(axis = 0L) %>%

        keras::k_eval()

    } else {

      reshaped_list[[i]] <- do.call(rbind, y[index])

    }

  }

  return(reshaped_list)

}

#' Evaluate matrices of true targets and predictions from layer with softmax activation. 

#' 

#' Compute confusion matrix, accuracy, categorical crossentropy and (optionally) AUC or AUPRC, given predictions and

#' true targets. AUC and AUPRC only possible for 2 targets. 

#' 

#' @param y Matrix of true target.

#' @param y_conf Matrix of predictions.

#' @param auc Whether to include AUC metric. Only possible for 2 targets. 

#' @param auprc Whether to include AUPRC metric. Only possible for 2 targets. 

#' @param label_names Names of corresponding labels. Length must be equal to number of columns of \code{y}.

#' @examplesIf reticulate::py_module_available("tensorflow")

#' y <- matrix(c(1, 0, 0, 0, 1, 1), ncol = 2)

#' y_conf <- matrix(c(0.3, 0.5, 0.1, 0.7, 0.5, 0.9), ncol = 2)

#' evaluate_softmax(y, y_conf, auc = TRUE, auprc = TRUE, label_names = c("A", "B")) 

#' 

#' @returns A list of evaluation results. 

#' @export    

evaluate_softmax <- function(y, y_conf, auc = FALSE, auprc = FALSE, label_names = NULL) {

  if (ncol(y) != 2 & (auc | auprc)) {

    message("Can only compute AUC or AUPRC if output layer with softmax acticvation has two neurons.")

    auc <- FALSE

    auprc <- FALSE

  }

  y_pred <- apply(y_conf, 1, which.max)

  y_true <- apply(y, 1, FUN = which.max) - 1

  df_true_pred <- data.frame(

    true = factor(y_true + 1, levels = 1:(length(label_names)), labels = label_names),

    pred = factor(y_pred, levels = 1:(length(label_names)), labels = label_names)

  )

  loss_per_class <- list()

  for (i in 1:ncol(y)) {

    index <- (y_true + 1) == i

    if (any(index)) {

      cce_loss_class <- tensorflow::tf$keras$losses$categorical_crossentropy(y[index, ], y_conf[index, ])

      loss_per_class[[i]] <- cce_loss_class$numpy()

    } else {

      loss_per_class[[i]] <- NA

    }

  }

  cm <- yardstick::conf_mat(df_true_pred, true, pred)

  confMat <- cm[[1]]

  acc <- sum(diag(confMat))/sum(confMat)

  loss <- mean(unlist(loss_per_class))

  for (i in 1:length(loss_per_class)) {

    loss_per_class[[i]] <- mean(unlist(loss_per_class[[i]]), na.rm = TRUE)

  }

  loss_per_class <- unlist(loss_per_class)

  m <- as.matrix(confMat)

  class_acc <- vector("numeric")

  for (i in 1:ncol(m)) {

    if (sum(m[ , i]) == 0) {

      class_acc[i] <- NA

    } else {

      class_acc[i] <- m[i, i]/sum(m[ , i])

    }

  }

  names(class_acc) <- label_names

  names(loss_per_class) <- label_names

  balanced_acc <- mean(class_acc)

  if (auc) {

    auc_list <- PRROC::roc.curve(

      scores.class0 = y_conf[ , 2],

      weights.class0 = y_true)

  } else {

    auc_list <- NULL

  }

  if (auprc) {

    auprc_list <- PRROC::pr.curve(

      scores.class0 = y_conf[ , 2],

      weights.class0 = y_true)

  } else {

    auprc_list <- NULL

  }

  return(list(confusion_matrix = confMat,

              accuracy = acc,

              categorical_crossentropy_loss = loss,

              #balanced_accuracy = balanced_acc,

              #loss_per_class = loss_per_class,

              #accuracy_per_class = class_acc,

              AUC = auc_list$auc,

              AUPRC = auprc_list$auc.integral))

}

#' Evaluate matrices of true targets and predictions from layer with sigmoid activation. 

#' 

#' Compute accuracy, binary crossentropy and (optionally) AUC or AUPRC, given predictions and

#' true targets. Outputs columnwise average.  

#' 

#' @inheritParams evaluate_model

#' @inheritParams evaluate_softmax

#' @param auc Whether to include AUC metric.

#' @param auprc Whether to include AUPRC metric. 

#' @examplesIf reticulate::py_module_available("tensorflow")

#' y <- matrix(sample(c(0, 1), 30, replace = TRUE), ncol = 3)

#' y_conf <- matrix(runif(n = 30), ncol = 3)

#' evaluate_sigmoid(y, y_conf, auc = TRUE, auprc = TRUE)

#' 

#' @returns A list of evaluation results. 

#' @export    

evaluate_sigmoid <- function(y, y_conf, auc = FALSE, auprc = FALSE, label_names = NULL) {

  y_pred <- ifelse(y_conf > 0.5, 1, 0)

  loss_per_class <- list()

  for (i in 1:ncol(y)) {

    bce_loss_class <- tensorflow::tf$keras$losses$binary_crossentropy(y[ , i], y_conf[ , i])

    loss_per_class[[i]] <- bce_loss_class$numpy()

  }

  loss_per_class <- unlist(loss_per_class)

  names(loss_per_class) <- label_names

  loss <- mean(unlist(loss_per_class))

  class_acc <- vector("numeric", ncol(y))

  for (i in 1:ncol(y)) {

    num_true_pred <-  sum(y[ , i] == y_pred[ , i])

    class_acc[i] <- num_true_pred /nrow(y)

  }

  names(class_acc) <- label_names

  acc <- mean(class_acc)

  if (auc) {

    auc_list <- purrr::map(1:ncol(y_conf), ~PRROC::roc.curve(

      scores.class0 = y_conf[ , .x],

      weights.class0 = y[ , .x]))

    auc_vector <- vector("numeric", ncol(y))

    for (i in 1:length(auc_vector)) {

      auc_vector[i] <- auc_list[[i]]$auc

    }

    na_count <- sum(is.na(auc_vector))

    if (na_count > 0) {

      message(paste(sum(na_count), ifelse(na_count > 1, "columns", "column"),

                    "removed from AUC evaluation since they contain only one label"))

    }

    AUC <- mean(auc_vector, na.rm = TRUE)

  } else {

    AUC <- NULL

  }

  if (auprc) {

    auprc_list <- purrr::map(1:ncol(y_conf), ~PRROC::pr.curve(

      scores.class0 = y_conf[ , .x],

      weights.class0 = y[ , .x]))

    auprc_vector <- vector("numeric", ncol(y))

    for (i in 1:length(auprc_vector)) {

      auprc_vector[i] <- auprc_list[[i]]$auc.integral

    }

    AUPRC <- mean(auprc_vector, na.rm = TRUE) 

  } else {

    AUPRC <- NULL

  }

  return(list(accuracy = acc,

              binary_crossentropy_loss = loss,

              #loss_per_class = loss_per_class,

              #accuracy_per_class = class_acc,

              AUC = AUC,

              AUPRC = AUPRC))

}

#' Evaluate matrices of true targets and predictions from layer with linear activation. 

#' 

#' Compute MAE and MSE, given predictions and

#' true targets. Outputs columnwise average.  

#' 

#' @inheritParams evaluate_model

#' @inheritParams evaluate_softmax

#' @param y_true Matrix of true labels.

#' @param y_pred Matrix of predictions.

#' @examplesIf reticulate::py_module_available("tensorflow")

#' y_true <- matrix(rnorm(n = 12), ncol = 3)

#' y_pred <- matrix(rnorm(n = 12), ncol = 3)

#' evaluate_linear(y_true, y_pred)

#' 

#' @returns A list of evaluation results. 

#' @export    

evaluate_linear <- function(y_true, y_pred, label_names = NULL) {

  loss_per_class_mse <- list()

  loss_per_class_mae <- list()

  for (i in 1:ncol(y_true)) {

    mse_loss_class <- tensorflow::tf$keras$losses$mean_squared_error(y_true[ ,i], y_pred[ , i])

    mae_loss_class <- tensorflow::tf$keras$losses$mean_absolute_error(y_true[ ,i], y_pred[ , i])

    loss_per_class_mse[[i]] <- mse_loss_class$numpy()

    loss_per_class_mae[[i]] <- mae_loss_class$numpy()

  }

  return(list(mse = mean(unlist(loss_per_class_mse)),

              mae = mean(unlist(loss_per_class_mae))))

}

#' Plot ROC

#' 

#' Compute ROC and AUC from target and prediction matrix and plot ROC. Target/prediction matrix should 

#' have one column if output of layer with sigmoid activation and two columns for softmax activation. 

#' 

#' @inheritParams evaluate_softmax

#' @inheritParams evaluate_linear

#' @param path_roc_plot Where to store ROC plot.

#' @param return_plot Whether to return plot.

#' @examples

#' y_true <- matrix(c(1, 0, 0, 0, 1, 1), ncol = 1)

#' y_conf <- matrix(runif(n = nrow(y_true)), ncol = 1)

#' p <- plot_roc(y_true, y_conf, return_plot = TRUE)

#' p

#' 

#' @returns A ggplot of ROC curve.

#' @export    

plot_roc <- function(y_true, y_conf, path_roc_plot = NULL,

                     return_plot = TRUE) {

  if (!all(y_true == 0 | y_true == 1)) {

    stop("y_true should only contain 0 and 1 entries")

  }

  if (is.matrix(y_true) && ncol(y_true) > 2) {

    stop("y_true can contain 1 or 2 columns")

  }

  if (is.matrix(y_true) && ncol(y_true) == 2) {

    y_true <- y_true[ , 1] 

    y_conf <- y_conf[ , 2]

  }

  if (stats::var(y_true) == 0) {

    stop("y_true contains just one label")

  }

  y_true <- as.vector(y_true)

  y_conf <- as.vector(y_conf)

  rocobj <-  pROC::roc(y_true, y_conf, quiet = TRUE)

  auc <- round(pROC::auc(y_true, y_conf, quiet = TRUE), 4)

  p <- pROC::ggroc(rocobj,  size = 1, color = "black")

  p <- p + ggplot2::theme_classic() + ggplot2::theme(aspect.ratio = 1) 

  p <- p + ggplot2::ggtitle(paste0('ROC Curve ', '(AUC = ', auc, ')'))

  p <- p + ggplot2::geom_abline(intercept = 1, linetype = 2, color = "grey50")

  p <- p + ggplot2::geom_vline(xintercept = 1, linetype = 2, color = "grey50")

  p <- p + ggplot2::geom_hline(yintercept = 1,  linetype = 2, color = "grey50")

  if (!is.null(path_roc_plot)) {

    ggplot2::ggsave(path_roc_plot, p)

  }

  if (return_plot) {

    return(p)

  } else {

    return(NULL)

  }

}

# plot_roc_auprc <- function(y_true, y_conf, path_roc_plot = NULL, path_auprc_plot = NULL,

#                            return_plot = TRUE, layer_activation = "softmax") {

#   

#   if (layer_activation == "softmax") {

#     

#     if (!all(y_true == 0 | y_true == 1)) {

#       stop("y_true should only contain 0 and 1 entries")

#     }

#     

#     if (ncol(y_true) != 2 & (auc | auprc)) {

#       message("Can only compute AUC or AUPRC if output layer with softmax acticvation has two neurons.")

#     }

#     

#     auc_list <- PRROC::roc.curve(

#       scores.class0 = y_conf[ , 2],

#       weights.class0 = y_true[ , 2], curve = TRUE)

#     

#     

#     auprc_list <- PRROC::pr.curve(

#       scores.class0 = y_conf[ , 2],

#       weights.class0 = y_true[ , 2], curve = TRUE)

#     

#     #auc_plot <- NULL

#     #auprc_plot <- NULL  

#     

#   }

#   

#   if (layer_activation == "sigmoid") {

#     

#     auc_list <- purrr::map(1:ncol(y_conf), ~PRROC::roc.curve(

#       scores.class0 = y_conf[ , .x],

#       weights.class0 = y[ , .x], curve = TRUE))

#     auc_vector <- vector("numeric", ncol(y))

#     

#     

#     auprc_list <- purrr::map(1:ncol(y_conf), ~PRROC::pr.curve(

#       scores.class0 = y_conf[ , .x],

#       weights.class0 = y[ , .x], curve = TRUE))

#     auprc_vector <- vector("numeric", ncol(y))

#     

#   }

#   

#   if (!is.null(path_roc_plot)) {

#     

#   }

#   

#   if (!is.null(path_auprc_plot)) {

#     

#   }

#   

# }