#' Make prediction for nucleotide sequence or entries in fasta/fastq file

#'

#' @description Removes layers (optional) from pretrained model and calculates states of fasta/fastq file or nucleotide sequence.

#' Writes states to h5 or csv file (access content of h5 output with \code{\link{load_prediction}} function).

#' There are several options on how to process an input file:

#' \itemize{

#' \item If `"one_seq"`, computes prediction for sequence argument or fasta/fastq file.

#' Combines fasta entries in file to one sequence. This means predictor sequences can contain elements from more than one fasta entry.

#' \item If `"by_entry"`, will output a separate file for each fasta/fastq entry.

#' Names of output files are: `output_dir` + "Nr" + i + `filename` + `output_type`, where i is the number of the fasta entry.

#' \item If `"by_entry_one_file"`, will store prediction for all fasta entries in one h5 file.

#' \item If `"one_pred_per_entry"`, will make one prediction for each entry by either picking random sample for long sequences

#' or pad sequence for short sequences.

#' }

#' 

#' @inheritParams get_generator

#' @inheritParams train_model

#' @inheritParams get_generator 

#' @inheritParams train_model

#' @param layer_name Name of layer to get output from. If `NULL`, will use the last layer.

#' @param path_input Path to fasta file.

#' @param sequence Character string, ignores path_input if argument given.

#' @param round_digits Number of decimal places.

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

#' @param include_seq Whether to include input sequence in h5 file.

#' @param output_format Either `"one_seq"`, `"by_entry"`, `"by_entry_one_file"`, `"one_pred_per_entry"`.

#' @param output_type `"h5"` or `"csv"`. If `output_format`` is `"by_entries_one_file", "one_pred_per_entry"` can only be `"h5"`.

#' @param return_states Return predictions as data frame. Only supported for output_format `"one_seq"`.

#' @param padding Either `"none"`, `"maxlen"`, `"standard"` or `"self"`.

#' \itemize{

#' \item If `"none"`, apply no padding and skip sequences that are too short.

#' \item If `"maxlen"`, pad with maxlen number of zeros vectors.

#' \item If `"standard"`, pad with zero vectors only if sequence is shorter than maxlen. Pads to minimum size required for one prediction.

#' \item If `"self"`, concatenate sequence with itself until sequence is long enough for one prediction.

#' Example: if sequence is "ACGT" and maxlen is 10, make prediction for "ACGTACGTAC". 

#' Only applied if sequence is shorter than maxlen.

#' }

#' @param verbose Boolean.

#' @param filename Filename to store states in. No file output if argument is `NULL`.

#' If `output_format = "by_entry"`, adds "_nr_" + "i" after name, where i is entry number.

#' @param output_dir Directory for file output.

#' @param use_quality Whether to use quality scores.

#' @param quality_string String for encoding with quality scores (as used in fastq format).

#' @param lm_format Either `"target_right"`, `"target_middle_lstm"`, `"target_middle_cnn"` or `"wavenet"`.

#' @param ... Further arguments for sequence encoding with \code{\link{seq_encoding_label}}.

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

#' # make prediction for single sequence and write to h5 file

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

#' vocabulary <- c("a", "c", "g", "t")

#' sequence <- paste(sample(vocabulary, 200, replace = TRUE), collapse = "")

#' output_file <- tempfile(fileext = ".h5")

#' predict_model(output_format = "one_seq", model = model, step = 10,

#'              sequence = sequence, filename = output_file, mode = "label")

#' 

#' # make prediction for fasta file with multiple entries, write output to separate h5 files

#' fasta_path <- tempfile(fileext = ".fasta")

#' create_dummy_data(file_path = fasta_path, num_files = 1,

#'                  num_seq = 5, seq_length = 100,

#'                  write_to_file_path = TRUE)

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

#' output_dir <- tempfile()

#' dir.create(output_dir)

#' predict_model(output_format = "by_entry", model = model, step = 10, verbose = FALSE,

#'                output_dir = output_dir, mode = "label", path_input = fasta_path)

#' list.files(output_dir)

#' 

#' @returns If `return_states = TRUE` returns a list of model predictions and position of corresponding sequences.

#' If additionally `include_seq = TRUE`, list contains sequence strings.

#' If `return_states = FALSE` returns nothing, just writes output to file(s).  

#' @export

predict_model <- function(model, output_format = "one_seq", layer_name = NULL, sequence = NULL, path_input = NULL,

                          round_digits = NULL, filename = "states.h5", step = 1, vocabulary = c("a", "c", "g", "t"),

                          batch_size = 256, verbose = TRUE, return_states = FALSE, 

                          output_type = "h5", padding = "none", use_quality = FALSE, quality_string = NULL,

                          mode = "label", lm_format = "target_right", output_dir = NULL,

                          format = "fasta", include_seq = FALSE, reverse_complement_encoding = FALSE,

                          ambiguous_nuc = "zero", ...) {

  stopifnot(padding %in% c("standard", "self", "none", "maxlen"))

  stopifnot(output_format %in% c("one_seq", "by_entry", "by_entry_one_file", "one_pred_per_entry"))

  if (output_format %in% c("by_entry_one_file", "one_pred_per_entry") & output_type == "csv") {

    message("by_entry_one_file or one_pred_per_entry only implemented for h5 output.

            Setting output_type to h5")

    output_type <- "h5"

  }

  if (output_format == "one_seq") {

    output_list <- predict_model_one_seq(layer_name = layer_name, sequence = sequence, path_input = path_input,

                                         round_digits = round_digits, filename = filename, step = step, vocabulary = vocabulary,

                                         batch_size = batch_size, verbose = verbose, return_states = return_states, 

                                         padding = padding, quality_string = quality_string, use_quality = use_quality,

                                         output_type = output_type, model = model, mode = mode, lm_format = lm_format,

                                         format = format, include_seq = include_seq, ambiguous_nuc = ambiguous_nuc,

                                         reverse_complement_encoding = reverse_complement_encoding, ...)

    return(output_list)

  }

  if (output_format == "by_entry") {

    predict_model_by_entry(layer_name = layer_name, path_input = path_input,

                           round_digits = round_digits, filename = filename, step = step, 

                           # vocabulary = vocabulary, quality_string = quality_string, ambiguous_nuc = ambiguous_nuc,

                           batch_size = batch_size, verbose = verbose, use_quality = use_quality,

                           output_type = output_type, model = model, mode = mode, lm_format = lm_format,

                           output_dir = output_dir, format = format, include_seq = include_seq,

                           padding = padding, reverse_complement_encoding = reverse_complement_encoding, ...)

  }

  if (output_format == "by_entry_one_file") {

    predict_model_by_entry_one_file(layer_name = layer_name, path_input = path_input,

                                    round_digits = round_digits, filename = filename, step = step, 

                                    # vocabulary = vocabulary, quality_string = quality_string, ambiguous_nuc = ambiguous_nuc, 

                                    batch_size = batch_size, verbose = verbose, use_quality = use_quality,

                                    model = model, mode = mode, lm_format = lm_format, format = format,

                                    padding = padding, include_seq = include_seq, 

                                    reverse_complement_encoding = reverse_complement_encoding, ...)

  }

  if (output_format == "one_pred_per_entry") {

    if (mode == "lm") {

      stop("one_pred_per_entry only implemented for label classification")

    }

    predict_model_one_pred_per_entry(layer_name = layer_name, path_input = path_input,

                                     round_digits = round_digits, filename = filename, 

                                     batch_size = batch_size, verbose = verbose, model = model, format = format,

                                     # ambiguous_nuc = ambiguous_nuc, use_quality = use_quality, vocabulary = vocabulary,

                                     reverse_complement_encoding = reverse_complement_encoding, ...)

  }

}

#' Write output of specific model layer to h5 or csv file.

#'

#' Removes layers (optional) from pretrained model and calculates states of fasta file, writes states to h5/csv file.

#' Function combines fasta entries in file to one sequence. This means predictor sequences can contain elements from more than one fasta entry.

#' h5 file also contains sequence and positions of targets corresponding to states.

#'

#' @inheritParams generator_fasta_lm

#' @param layer_name Name of layer to get output from. If `NULL`, will use the last layer.

#' @param path_input Path to fasta file.

#' @param sequence Character string, ignores path_input if argument given.

#' @param round_digits Number of decimal places.

#' @param batch_size Number of samples to evaluate at once. Does not change output, only relevant for speed and memory.

#' @param step Frequency of sampling steps.

#' @param filename Filename to store states in. No file output if argument is `NULL`.

#' @param vocabulary Vector of allowed characters, character outside vocabulary get encoded as 0-vector.

#' @param return_states Logical scalar, return states matrix.

#' @param ambiguous_nuc `"zero"` or `"equal"`. 

#' @param verbose Whether to print model before and after removing layers.

#' @param output_type Either `"h5"` or `"csv"`.

#' @param model A keras model. 

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

#' @param format Either `"fasta"` or `"fastq"`.

#' @param include_seq Whether to include input sequence in h5 file.

#' @param ... Further arguments for sequence encoding with \code{\link{seq_encoding_label}}.

#' @noRd

predict_model_one_seq <- function(model, layer_name = NULL, sequence = NULL, path_input = NULL, round_digits = 2,

                                  filename = "states.h5", step = 1, vocabulary = c("a", "c", "g", "t"), batch_size = 256, verbose = TRUE,

                                  return_states = FALSE, target_len = 1, use_quality = FALSE, quality_string = NULL,

                                  output_type = "h5", mode = "lm", lm_format = "target_right",

                                  ambiguous_nuc = "zero", padding = "none", format = "fasta", output_dir = NULL,

                                  include_seq = TRUE, reverse_complement_encoding = FALSE, ...) {

  vocabulary <- stringr::str_to_lower(vocabulary)

  stopifnot(mode %in% c("lm", "label"))

  stopifnot(output_type %in% c("h5", "csv"))

  if (!is.null(quality_string)) use_quality <- TRUE

  # if (!is.null(quality_string) & !is.null(sequence)) {

  #   stopifnot(length(sequence) == length(quality_to_probability(quality_string)))

  # }

  file_output <- !is.null(filename)

  if (!file_output) {

    if (!return_states) stop("If filename is NULL, return_states must be TRUE; otherwise function produces no output.")

    filename <- tempfile(fileext = paste0(".", output_type))

  }

  stopifnot(batch_size > 0)

  stopifnot(!file.exists(filename))

  if (reverse_complement_encoding) {

    test_len <- length(vocabulary) != 4

    if (test_len || all(sort(stringr::str_to_lower(vocabulary)) != c("a", "c", "g", "t"))) {

      stop("reverse_complement_encoding only implemented for A,C,G,T vocabulary")

    }

  }

  # token for ambiguous nucleotides

  for (i in letters) {

    if (!(i %in% stringr::str_to_lower(vocabulary))) {

      amb_nuc_token <- i

      break

    }

  }

  tokenizer <- keras::fit_text_tokenizer(keras::text_tokenizer(char_level = TRUE, lower = TRUE, oov_token = "0"), c(vocabulary, amb_nuc_token))

  if (is.null(layer_name)) {

    layer_name <- model$output_names

    if (verbose) message(paste("layer_name not specified. Using layer", layer_name))

  }

  if (!is.null(sequence) && (!missing(sequence) & sequence != "")) {

    nt_seq <- sequence %>% stringr::str_to_lower()

  } else {

    if (format == "fasta") {

      fasta.file <- microseq::readFasta(path_input)

    }

    if (format == "fastq") {

      fasta.file <- microseq::readFastq(path_input)

    }

    if (nrow(fasta.file) > 1 & verbose) {

      text_1 <- paste("Your file has", nrow(fasta.file), "entries. 'one_seq'  output_format will concatenate them to a single sequence.\n")

      text_2 <- "Use 'by_entry' or 'by_entry_one_file' output_format to evaluate them separately."

      message(paste0(text_1, text_2))

    }

    nt_seq <- paste(fasta.file$Sequence, collapse = "") %>% stringr::str_to_lower()

  }

  # tokenize ambiguous nt

  pattern <- paste0("[^", paste0(vocabulary, collapse = ""), "]")

  nt_seq <- stringr::str_replace_all(string = nt_seq, pattern = pattern, amb_nuc_token)

  if (use_quality & is.null(quality_string)) {

    quality_string <- paste(fasta.file$Quality, collapse = "")

  }

  # extract maxlen

  target_middle <- ifelse(mode == "lm" && (lm_format %in% c("target_middle_lstm", "target_middle_cnn")), TRUE, FALSE)

  if (!target_middle) {

    if (reverse_complement_encoding) {

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

    } else {

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

    }

  } else {

    maxlen_1 <- model$input[[1]]$shape[[2]]

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

    maxlen <- maxlen_1 + maxlen_2

  }

  total_seq_len <- ifelse(mode == "lm", maxlen + target_len, maxlen)

  # pad sequence

  unpadded_seq_len <- nchar(nt_seq)

  pad_len <- 0

  if (padding == "maxlen") {

    pad_len <- maxlen

  }

  if (padding == "standard" & (unpadded_seq_len < total_seq_len)) {

    pad_len <- total_seq_len - unpadded_seq_len

  }

  if (padding == "self"  & (unpadded_seq_len < total_seq_len)) {

    nt_seq <- strrep(nt_seq, ceiling(total_seq_len / unpadded_seq_len))

    nt_seq <- substr(nt_seq, 1, total_seq_len)

    if (use_quality) {

      quality_string <- strrep(quality_string, ceiling(total_seq_len / unpadded_seq_len))

      quality_string <- substr(quality_string, 1, total_seq_len) 

    }

  } else {

    nt_seq <- paste0(strrep("0", pad_len), nt_seq)

    if (use_quality) quality_string <- paste0(strrep("0", pad_len), quality_string)

  }

  if (nchar(nt_seq) < total_seq_len) {

    stop(paste0("Input sequence is shorter than required length (", total_seq_len, "). Use padding argument to pad sequence to bigger size."))

  }

  if (use_quality) {

    quality_vector <- quality_string %>% quality_to_probability()

  } else {

    quality_vector <- NULL

  }

  # start of samples

  start_indices <- seq(1, nchar(nt_seq) - total_seq_len + 1, by = step)

  num_samples <- length(start_indices)

  check_layer_name(model, layer_name)

  model <- tensorflow::tf$keras$Model(model$input, model$get_layer(layer_name)$output)

  if (verbose) {

    cat("Computing output for model at layer", layer_name,  "\n")

    print(model)

  }

  # extract number of neurons in last layer

  if (length(model$output$shape$dims) == 3) {

    if (!("lstm" %in% stringr::str_to_lower(model$output_names))) {

      stop("Output dimension of layer is > 1, format not supported yet")

    }

    layer.size <- model$output$shape[[3]]

  } else {

    layer.size <- model$output$shape[[2]]

  }

  # tokenize sequence

  nt_seq <- stringr::str_to_lower(nt_seq)

  tokSeq <- keras::texts_to_sequences(tokenizer, nt_seq)[[1]] - 1

  # seq end position

  pos_arg <- start_indices + total_seq_len - pad_len - 1

  if (include_seq) {

    output_seq <- substr(nt_seq, pad_len + 1, nchar(nt_seq))

  }

  # create h5 file to store states

  if (output_type == "h5") {

    h5_file <- hdf5r::H5File$new(filename, mode = "w")

    h5_file[["multi_entries"]] <- FALSE

    h5_file[["sample_end_position"]] <- pos_arg

    if (include_seq) h5_file[["sequence"]] <- output_seq

  }

  number_batches <- ceiling(length(start_indices)/batch_size)

  pred_list <- vector("list", number_batches)

  col_names <- c(as.character(1:layer.size), "sample_end_position")

  # subset input for target middle

  if (mode == "lm" && lm_format %in% c("target_middle_lstm", "target_middle_cnn")) {

    index_x_1 <- 1:ceiling((total_seq_len - target_len)/2)

    index_x_2 <- (max(index_x_1) + target_len + 1) : total_seq_len

  } 

  for (i in 1:number_batches) {

    index_start <- ((i - 1) * batch_size) + 1

    index_end <- min(c(num_samples + 1, index_start + batch_size)) - 1

    index <- index_start : index_end 

    x <- seq_encoding_label(sequence = tokSeq, 

                            maxlen = total_seq_len,

                            vocabulary = vocabulary,

                            start_ind = start_indices[index],

                            ambiguous_nuc = ambiguous_nuc,

                            tokenizer = NULL,

                            adjust_start_ind = FALSE,

                            quality_vector = quality_vector,

                            ...

    )

    if (mode == "lm" && lm_format == "target_middle_lstm") {

      x1 <- x[ , index_x_1, ]

      x2 <- x[ , index_x_2, ]

      if (length(index_x_1) == 1 | dim(x)[1] == 1) {

        x1 <- array(x1, dim = c(1, dim(x1)))

      }

      if (length(index_x_2) == 1 | dim(x)[1] == 1) {

        x2 <- array(x2, dim = c(1, dim(x2)))

      }

      x2 <- x2[ , dim(x2)[2]:1, ] # reverse order

      if (length(dim(x2)) == 2) {

        x2 <- array(x2, dim = c(1, dim(x2)))

      }

      x <- list(x1, x2)

    } 

    if (mode == "lm" && lm_format == "target_middle_cnn") {

      x <- x[ , c(index_x_1, index_x_2), ]

    } 

    if (reverse_complement_encoding) x <- list(x, reverse_complement_tensor(x))

    y <- stats::predict(model, x, verbose = 0)

    if (!is.null(round_digits)) y <- round(y, round_digits)

    pred_list[[i]] <- y

  }

  states <- do.call(rbind, pred_list)

  if (file_output) {

    if (output_type == "h5") {

      h5_file[["states"]] <- states

      h5_file$close_all()

    } else {

      col_names <- paste0("N", 1:ncol(states))

      colnames(states) <- col_names 

      utils::write.csv(x = states, file = filename, row.names = FALSE)

    }

  }

  if (return_states) {

    output_list <- list()

    output_list$states <- states

    output_list$sample_end_position <- pos_arg

    if (include_seq) output_list$sequence <- output_seq

    return(output_list)

  } else {

    return(NULL)

  }

}

#' Write states to h5 file

#'

#' @description Removes layers (optional) from pretrained model and calculates states of fasta file, writes a separate

#' h5 file for every fasta entry in fasta file. h5 files also contain the nucleotide sequence and positions of targets corresponding to states.

#' Names of output files are: file_path + "Nr" + i + filename + output_type, where i is the number of the fasta entry.

#'

#' @param filename Filename to store states, function adds "_nr_" + "i" after name, where i is entry number.

#' @param output_dir Path to folder, where to write output.

#' @noRd

predict_model_by_entry <- function(model, layer_name = NULL, path_input, round_digits = 2,

                                   filename = "states.h5", output_dir = NULL, step = 1, vocabulary = c("a", "c", "g", "t"),

                                   batch_size = 256, output_type = "h5", mode = "lm",

                                   lm_format = "target_right", format = "fasta", use_quality = FALSE,

                                   reverse_complement_encoding = FALSE, padding = "none", 

                                   verbose = FALSE, include_seq = FALSE, ambiguous_nuc = "zero", ...) {

  stopifnot(mode %in% c("lm", "label"))

  stopifnot(!is.null(filename))

  stopifnot(!is.null(output_dir))

  if (endsWith(filename, paste0(".", output_type))) {

    filename <- stringr::str_remove(filename, paste0(".", output_type, "$"))

    filename <- basename(filename)

  }

  # extract maxlen

  target_middle <- ifelse(mode == "lm" && (lm_format %in% c("target_middle_lstm", "target_middle_cnn")), TRUE, FALSE)

  if (!target_middle) {

    if (reverse_complement_encoding) {

      model$input[[1]]$shape[[2]]

    } else {

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

    }

  } else {

    maxlen_1 <- model$input[[1]]$shape[[2]]

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

    maxlen <- maxlen_1 + maxlen_2

  }

  # load fasta file

  if (format == "fasta") {

    fasta.file <- microseq::readFasta(path_input)

  }

  if (format == "fastq") {

    fasta.file <- microseq::readFastq(path_input)

  }

  df <- fasta.file[ , c("Sequence", "Header")]

  names(df) <- c("seq", "header")

  rownames(df) <- NULL

  num_skipped_seq <- 0

  for (i in 1:nrow(df)) {

    # skip entry if too short

    if ((nchar(df[i, "seq"]) < maxlen) & padding == "none") {

      num_skipped_seq <- num_skipped_seq + 1

      next

    } 

    if (use_quality) {

      quality_string <- fasta.file$Quality[i]

    } else {

      quality_string <- NULL

    }

    current_file <- paste0(output_dir, "/", filename, "_nr_", as.character(i), ".", output_type)

    predict_model_one_seq(layer_name = layer_name, sequence = df[i, "seq"],

                          round_digits = round_digits, path_input = path_input,

                          filename = current_file, quality_string = quality_string,

                          step = step, vocabulary = vocabulary, batch_size = batch_size,

                          verbose = ifelse(i > 1, FALSE, verbose), 

                          output_type = output_type, mode = mode,

                          lm_format = lm_format, model = model, include_seq = include_seq,

                          padding = padding, ambiguous_nuc = "zero", 

                          reverse_complement_encoding = reverse_complement_encoding, ...)

  }

  if (verbose & num_skipped_seq > 0) {

    message(paste0("Skipped ", num_skipped_seq,

                   ifelse(num_skipped_seq == 1, " entry", " entries"),

                   ". Use different padding option to evaluate all."))

  }  

}

#' Write states to h5 file

#'

#' @description Removes layers (optional) from pretrained model and calculates states of fasta file,

#' writes separate states matrix in one .h5 file for every fasta entry.

#' h5 file also contains the nucleotide sequences and positions of targets corresponding to states.

#' @noRd

predict_model_by_entry_one_file <- function(model, path_input, round_digits = 2, filename = "states.h5",

                                            step = 1,  vocabulary = c("a", "c", "g", "t"), batch_size = 256, layer_name = NULL,

                                            verbose = TRUE, mode = "lm", use_quality = FALSE,

                                            lm_format = "target_right", padding = "none",

                                            format = "fasta", include_seq = TRUE, reverse_complement_encoding = FALSE, ...) {

  vocabulary <- stringr::str_to_lower(vocabulary)

  stopifnot(mode %in% c("lm", "label"))

  target_middle <- ifelse(mode == "lm" && (lm_format %in% c("target_middle_lstm", "target_middle_cnn")), TRUE, FALSE)

  # extract maxlen

  if (!target_middle) {

    if (reverse_complement_encoding) {

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

    } else {

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

    }

  } else {

    maxlen_1 <- model$input[[1]]$shape[[2]]

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

    maxlen <- maxlen_1 + maxlen_2

  }

  # extract number of neurons in last layer

  if (length(model$output$shape$dims) == 3) {

    if (!("lstm" %in% stringr::str_to_lower(model$output_names))) {

      stop("Output dimension of layer is > 1, format not supported yet")

    }

    layer.size <- model$output$shape[[3]]

  } else {

    layer.size <- model$output$shape[[2]]

  }

  # load fasta file

  if (format == "fasta") {

    fasta.file <- microseq::readFasta(path_input)

  }

  if (format == "fastq") {

    fasta.file <- microseq::readFastq(path_input)

  }

  df <- fasta.file[ , c("Sequence", "Header")]

  names(df) <- c("seq", "header")

  rownames(df) <- NULL

  if (verbose) {

    # check if names are unique

    if (length(df$header) != length(unique(df$header))) {

      message("Header names are not unique, adding '_header_x' to names (x being the header number)")

      df$header <- paste0(df$header, paste0("_header_", 1:length(df$header)))

    }

  }

  # create h5 file to store states

  h5_file <- hdf5r::H5File$new(filename, mode = "w")

  h5_file[["multi_entries"]] <- TRUE

  states.grp <- h5_file$create_group("states")

  sample_end_position.grp <- h5_file$create_group("sample_end_position")

  if (include_seq) seq.grp <- h5_file$create_group("sequence")

  num_skipped_seq <- 0

  for (i in 1:nrow(df)) {

    #seq_name <- df$header[i]

    seq_name <- paste0("entry_", i)

    temp_file <- tempfile(fileext = ".h5")

    # skip entry if too short

    if ((nchar(df[i, "seq"]) < maxlen) & padding == "none") {

      num_skipped_seq <- num_skipped_seq + 1

      next

    } 

    if (use_quality) {

      quality_string <- fasta.file$Quality[i]

    } else {

      quality_string <- NULL

    }

    output_list <- predict_model_one_seq(layer_name = layer_name, sequence = df$seq[i], path_input = path_input,

                                         round_digits = round_digits, filename = temp_file, step = step, vocabulary = vocabulary,

                                         batch_size = batch_size, return_states = TRUE, quality_string = quality_string,

                                         output_type = "h5", model = model, mode = mode, lm_format = lm_format,

                                         ambiguous_nuc = "zero", verbose = ifelse(i > 1, FALSE, verbose), 

                                         padding = padding, format = format, include_seq = include_seq,

                                         reverse_complement_encoding = reverse_complement_encoding, ...)

    states.grp[[seq_name]] <- output_list$states

    sample_end_position.grp[[seq_name]] <- output_list$sample_end_position

    if (include_seq) seq.grp[[seq_name]] <- output_list$sequence

  }

  if (verbose & num_skipped_seq > 0) {

    message(paste0("Skipped ", num_skipped_seq,

                   ifelse(num_skipped_seq == 1, " entry", " entries"),

                   ". Use different padding option to evaluate all."))

  }  

  h5_file$close_all()

}

#' Get states for label classification model.

#'

#' Computes output at specified model layer. Forces every fasta entry to have length maxlen by either padding sequences shorter than maxlen or taking random subsequence for

#' longer sequences.

#'

#' @inheritParams predict_model_one_seq

#' @noRd

predict_model_one_pred_per_entry <- function(model, layer_name = NULL, path_input, round_digits = 2, format = "fasta",

                                             ambiguous_nuc = "zero", filename = "states.h5", padding = padding,

                                             vocabulary = c("a", "c", "g", "t"), batch_size = 256, verbose = TRUE,

                                             return_states = FALSE, reverse_complement_encoding = FALSE, 

                                             include_seq = FALSE, use_quality = FALSE, ...) {

  vocabulary <- stringr::str_to_lower(vocabulary)

  file_type <- "h5"

  stopifnot(batch_size > 0)

  stopifnot(!file.exists(filename))

  # token for ambiguous nucleotides

  for (i in letters) {

    if (!(i %in% stringr::str_to_lower(vocabulary))) {

      amb_nuc_token <- i

      break

    }

  }

  tokenizer <- keras::fit_text_tokenizer(keras::text_tokenizer(char_level = TRUE, lower = TRUE, oov_token = "N"), vocabulary)

  if (is.null(layer_name)) {

    layer_name <- model$output_names

    if (verbose) message(paste("layer_name not specified. Using layer", layer_name))

  }

  # extract maxlen

  if (reverse_complement_encoding) {

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

  } else {

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

  }

  if (format == "fasta") {

    fasta.file <- microseq::readFasta(path_input)

  }

  if (format == "fastq") {

    fasta.file <- microseq::readFastq(path_input)

  }

  num_samples <- nrow(fasta.file)

  nucSeq <- as.character(fasta.file$Sequence)

  seq_length <- nchar(fasta.file$Sequence)

  if (use_quality) {

    quality_string <- vector("character", length(nucSeq)) 

  } else {

    quality_string <- NULL

  }

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

    # take random subsequence

    if (seq_length[i] > maxlen) {

      start <- sample(1 : (seq_length[i] - maxlen + 1) , size = 1)

      nucSeq[i] <- substr(nucSeq[i], start = start, stop = start + maxlen - 1)

      if (use_quality) {

        quality_string[i] <- substr(fasta.file$Quality[i], start = start, stop = start + maxlen - 1)

      }  

    }

    # pad sequence

    if (seq_length[i] < maxlen) {

      nucSeq[i] <- paste0(paste(rep("N", maxlen - seq_length[i]), collapse = ""), nucSeq[i])

      if (use_quality) {

        quality_string[i] <- paste0(paste(rep("0", maxlen - seq_length[i]), collapse = ""), fasta.file$Quality[i])

      }  

    }

  }

  model <- tensorflow::tf$keras$Model(model$input, model$get_layer(layer_name)$output)

  if (verbose) {

    cat("Computing output for model at layer", layer_name,  "\n")

    print(model)

  } 

  # extract number of neurons in last layer

  if (length(model$output$shape$dims) == 3) {

    if (!("lstm" %in% stringr::str_to_lower(model$output_names))) {

      stop("Output dimension of layer is > 1, format not supported yet")

    }

    layer.size <- model$output$shape[[3]]

  } else {

    layer.size <- model$output$shape[[2]]

  }

  if (file_type == "h5") {

    # create h5 file to store states

    h5_file <- hdf5r::H5File$new(filename, mode = "w")

    if (!missing(path_input)) h5_file[["fasta_file"]] <- path_input

    h5_file[["header_names"]] <- fasta.file$Header

    if (include_seq) h5_file[["sequences"]] <- nucSeq

    h5_file[["states"]] <- array(0, dim = c(0, layer.size))

    h5_file[["multi_entries"]] <- FALSE

    writer <- h5_file[["states"]]

  }

  rm(fasta.file)

  #nucSeq <- paste(nucSeq, collapse = "") %>% stringr::str_to_lower()

  number_batches <- ceiling(num_samples/batch_size)

  if (verbose) cat("Evaluating", number_batches, ifelse(number_batches > 1, "batches", "batch"), "\n")

  row <- 1

  string_start_index <- 1 

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

  percentage_index <- 1

  if (number_batches > 1) {

    for (i in 1:(number_batches - 1)) {

      string_end_index <-string_start_index + batch_size - 1 

      char_seq <- nucSeq[string_start_index : string_end_index] %>% paste(collapse = "") 

      if (use_quality) {

        quality_string_subset <- quality_string[string_start_index : string_end_index] %>% paste(collapse = "") 

      } else {

        quality_string_subset <- NULL

      }

      if (i == 1) start_ind <- seq(1, nchar(char_seq), maxlen)

      one_hot_batch <- seq_encoding_label(sequence = NULL, maxlen = maxlen, vocabulary = vocabulary,

                                          start_ind = start_ind, ambiguous_nuc = ambiguous_nuc, 

                                          char_sequence = char_seq, quality_vector = quality_string_subset,

                                          tokenizer = tokenizer, adjust_start_ind = TRUE, ...) 

      if (reverse_complement_encoding) one_hot_batch <- list(one_hot_batch, reverse_complement_tensor(one_hot_batch))

      activations <- keras::predict_on_batch(model, one_hot_batch)

      writer[row : (row + batch_size - 1), ] <- activations

      row <- row + batch_size

      string_start_index <- string_end_index + 1

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

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

        percentage_index <- percentage_index + 1

      }

    }

  }

  # last batch might be shorter

  char_seq <- nucSeq[string_start_index : length(nucSeq)] %>% paste(collapse = "") 

  if (use_quality) {

    quality_string_subset <- quality_string[string_start_index : length(nucSeq)] %>% paste(collapse = "") 

  } else {

    quality_string_subset <- NULL

  }

  one_hot_batch <- seq_encoding_label(sequence = NULL, maxlen = maxlen, vocabulary = vocabulary,

                                      start_ind = seq(1, nchar(char_seq), maxlen), ambiguous_nuc = "zero", nuc_dist = NULL,

                                      quality_vector = quality_string_subset, use_coverage = FALSE, max_cov = NULL,

                                      cov_vector = NULL, n_gram = NULL, n_gram_stride = 1, char_sequence = char_seq,

                                      tokenizer = tokenizer, adjust_start_ind = TRUE, ...) 

  if (reverse_complement_encoding) one_hot_batch <- list(one_hot_batch, reverse_complement_tensor(one_hot_batch))

  activations <- keras::predict_on_batch(model, one_hot_batch)

  writer[row : num_samples, ] <- activations[1 : length(row:num_samples), ]

  if (verbose) cat("Progress: 100 % \n")

  if (return_states & (file_type == "h5")) states <- writer[ , ]

  if (file_type == "h5") h5_file$close_all()

  if (return_states) return(states)

}

#' Read states from h5 file

#'

#' Reads h5 file created by  \code{\link{predict_model}} function.

#'

#' @param h5_path Path to h5 file.

#' @param rows Range of rows to read. If `NULL` read everything.

#' @param get_sample_position Return position of sample corresponding to state if `TRUE`.

#' @param get_seq Return nucleotide sequence if `TRUE`.

#' @param verbose Boolean.

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

#' # make prediction for single sequence and write to h5 file

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

#' vocabulary <- c("a", "c", "g", "t")

#' sequence <- paste(sample(vocabulary, 200, replace = TRUE), collapse = "")

#' output_file <- tempfile(fileext = ".h5")

#' predict_model(output_format = "one_seq", model = model, step = 10,

#'               sequence = sequence, filename = output_file, mode = "label")

#' load_prediction(h5_path = output_file)

#' 

#' @returns A list of data frames, containing model predictions and sequence positions.

#' @export

load_prediction <- function(h5_path, rows = NULL, verbose = FALSE,

                            get_sample_position = FALSE, get_seq = FALSE) {

  if (is.null(rows)) complete <- TRUE

  h5_file <- hdf5r::H5File$new(h5_path, mode = "r")

  multi_entries <- ifelse(h5_file[["multi_entries"]][], TRUE, FALSE)

  if (!multi_entries) {

    number_entries <- 1

  } else {

    entry_names <- names(h5_file[["states"]])

    number_entries <- length(entry_names)

    output_list <- list()

  }

  train_mode <- "label"

  if (get_sample_position & !any(c("sample_end_position", "target_position") %in% names(h5_file))) {

    get_sample_position <- FALSE

    message("File does not contain target positions.")

  }

  if (!multi_entries) {

    read_states <- h5_file[["states"]]

    if (get_sample_position) {

      read_targetPos <- h5_file[["sample_end_position"]]

    }

    if (verbose) {

      cat("states matrix has", dim(read_states[ , ])[1], "rows and",  dim(read_states[ , ])[2], "columns \n")

    }

    if (complete) {

      states <- read_states[ , ]

      if (get_sample_position) {

        targetPos <- read_targetPos[ ]

      }

    } else {

      states <- read_states[rows, ]

      if (get_sample_position) {

        targetPos <- read_targetPos[rows]

      }

    }

    if (is.null(dim(states))) {

      states <- matrix(states, nrow = 1)

    }

    contains_seq <- FALSE

    if (get_seq) {

      if ("sequence" %in% names(h5_file)) {

        contains_seq <- TRUE

        sequence <- h5_file[["sequence"]][]

      } else {

        contains_seq <- FALSE

        message("File does not contain sequence.")

      }

    }

    h5_file$close_all()

    output_list <- list(states = states)

    if (get_sample_position) {

      if (train_mode == "label") {

        output_list$sample_end_position <- targetPos

      } else {

        output_list$target_position <- targetPos

      }

    }

    if (get_seq && contains_seq) {

      output_list$sequence <- sequence

    }

    return(output_list)

    # multi entries

  } else {

    if (verbose) {

      cat("file contains", number_entries, "entries \n")

    }

    if (get_sample_position) {

      target_name <- "sample_end_position"

    }

    if (get_seq & !("sequence" %in% names(h5_file))) {

      message("File does not contain sequence.")

      get_seq <- FALSE

    }

    for (i in 1:number_entries) {

      entry_name <- entry_names[i]

      states <- h5_file[["states"]][[entry_name]][ , ]

      if (is.null(dim(states))) {

        states <- matrix(states, nrow = 1)

      }

      if (get_seq) {

        sequence <- h5_file[["sequence"]][[entry_name]][ ]

      }

      if (get_sample_position) {

        targetPos <- h5_file[[target_name]][[entry_name]][ ]

      }

      if (!complete) {

        states <- states[rows, ]

        if (is.null(dim(states))) {

          states <- matrix(states, nrow = 1)

        }

        if (get_sample_position) {

          targetPos <- hdf5r::h5attr(read_states, target_name)[rows]

        }

      }

      if (get_sample_position) {

        l <- list(states = states, sample_end_position = targetPos)

      } else {

        l <- list(states = states)

      }

      if (get_seq) {

        l[["sequence"]] <- sequence

      }

      output_list[[entry_name]] <- l

    }

    h5_file$close_all()

    names(output_list) <- entry_names

    return(output_list)

  }

}

#' Create summary of predictions 

#'

#' Create summary data frame for confidence predictions over 1 or several state files or a data frame.

#' Columns in file or data frame should be confidence predictions for one class,

#' i.e. each rows should sum to 1 and have nonnegative entries. 

#' Output data frame contains average confidence scores, max score and percentage of votes for each class.

#'

#' @param states_path Folder containing state files or a single file with same ending as `file_type`.

#' @param label_names Names of predicted classes.

#' @param file_type `"h5"` or `"csv"`.

#' @param df A states data frame. Ignore `states_dir` argument if not `NULL`.

#' @examples 

#' m <- c(0.9,  0.1, 0.2, 0.01,

#'        0.05, 0.7, 0.2, 0,

#'        0.05, 0.2, 0.6, 0.99) %>% matrix(ncol = 3)

#' 

#' label_names <- paste0("class_", 1:3)

#' df <- as.data.frame(m)

#' pred_summary <- summarize_states(label_names = label_names, df = df)

#' pred_summary

#' 

#' @returns A data frame of predictions summaries.

#' @export

summarize_states <- function(states_path = NULL, label_names = NULL, file_type = "h5", df = NULL) {

  if (!is.null(df)) {

    states_path <- NULL

  }

  if (is.null(states_path)) {

    state_files <- 1

  } else {

    if (endsWith(states_path, file_type)) {

      state_files <- states_path

    } else {

      state_files <- list.files(states_path, full.names = TRUE)

    }

  }

  if (!is.null(label_names)) {

    num_labels <- length(label_names)

  }

  summary_list <- list()

  for (state_file in state_files) {

    if (is.null(df)) {

      if (file_type == "h5") {

        df <- load_prediction(h5_path = state_file, get_sample_position = FALSE, verbose = FALSE)

        df <- as.data.frame(df$states)

      }

      if (file_type == "csv") {

        df <- utils::read.csv(state_file)

        if (ncol(df) != num_labels) {

          df <- utils::read.csv2(state_file)

        }

      }

    }

    if (state_file == state_files[1] & is.null(label_names)) {

      label_names <- paste0("X_", 1:ncol(df))

      num_labels <- length(label_names)

    }

    stopifnot(ncol(df) == num_labels)

    names(df) <- c(label_names)

    mean_df <- data.frame(matrix(0, nrow = 1, ncol = num_labels))

    names(mean_df) <- paste0("mean_conf_", label_names)

    max_df <- data.frame(matrix(0, nrow = 1, ncol = num_labels))

    names(max_df) <- paste0("max_conf_", label_names)

    for (label in label_names) {

      mean_df[[paste0("mean_conf_", label)]] <-  mean(df[[label]])

      max_df[[paste0("max_conf_", label)]] <-  max(df[[label]])

    }

    vote_distribution <- apply(df[label_names], 1, which.max)

    vote_perc <- table(factor(vote_distribution, levels = 1:length(label_names)))/length(vote_distribution)

    votes_df <- data.frame(matrix(vote_perc, nrow = 1, ncol = num_labels))

    names(votes_df) <- paste0("vote_perc_", label_names)

    mean_prediction <- label_names[which.max(unlist(mean_df))]

    max_prediction <- label_names[which.max(unlist(max_df))]

    vote_prediction <- label_names[which.max(vote_perc)]

    if (is.null(states_path)) {

      file_name <- NA

    } else {

      file_name <- basename(state_file)

    }

    summary_list[[state_file]] <- data.frame(file_name, mean_df, max_df, votes_df,

                                             mean_prediction, max_prediction, vote_prediction,

                                             num_prediction = nrow(df))

  }

  summary_df <- data.table::rbindlist(summary_list)

  return(summary_df)

}