--- a
+++ b/R/predict.R
@@ -0,0 +1,1014 @@
+#' 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)
+}