inner_procedure.r

#' Perform Habitat Sampling and Probability Mapping
#'
#'This is the function that performs: initiate sampling and model building
#'
#' @param raster satellite time series stack (rasterBrickObject) or just any type of image (*rasterObject)
#' @param nb_samples number of spatial locations
#' @param sample_type distribution of spatial locations c("random","regular")
#' @param nb_models number of models (independent classifiers) to collect
#' @param nb_mean number of iterations for model accuracy
#' @param buffer distance (in m) for new sample collection around initial samples (depends on pixel size)
#' @param reference reference spectra either SpatialPointsDataFrame (shape file) or data.frame with lines = classes, column = predictors]
#' @param model which machine learning classifier to use c("rf", "svm") for random forest or suppurt vector machine implementation
#' @param mtry number of predictor used at random forest splitting nodes (mtry << n predictors)
#' @param mod.error threshold for model error until which iteration is being executed
#' @param in.memory boolean for raster processing (memory = "TRUE", from disk = "FALSE")
#' @param last only true for one class classifier c("FALSE", TRUE")
#' @param seed set seed for reproducable results
#' @param init.seed "sample" for new or use run1@seeds to reproduce previous steps
#' @param save_runs if the user wants to save the runs, if TRUE the complete Habitat Class object is returned
#' @param parallel_mode run loops in parallel
#' @param max_num_cores maximum number of cores for parallelism
#' @param progress_bar if true use a normal progress bar, otherwise a shiny progress bar
#' @param session shiny session, needed for progress bar update
#'
#' @return a list with 5 elements:
#' 1) returns 0 succeeded, 1 increase init.samples, or 2 increase init.samples and nb_models
#' 2) An index
#' 3) num_models number of models selected
#' 4) Accuracy vector
#' 5) A vector with a Habitat objects, each consisting of 7 slots: \cr
#' run1@models - list of selected classifiers (only if save_runs is TRUE) \cr
#' run1@ref_samples - list of SpatialPointsDataFrames with same length as run1@models holding reference labels [1,2] for each selected model (only if save_runs is TRUE) \cr
#' run1@switch - vector of length run1@models indicating if target class equals 2, if not NA the labels need to be switched (only if save_runs is TRUE) \cr
#' run1@layer - raster map of habitat type probability \cr
#' run1@mod_all - list of all classifiers (equals nb_models) (only if save_runs is TRUE) \cr
#' run1@class_ind - vector of predictive distance measure for all habitats (only if save_runs is TRUE) \cr
#' run1@seeds - vector of seeds for random sampling (only if save_runs is TRUE) \cr
#' all files are saved with step number, the *.tif file is additionally saved with class names
#' @keywords internal
#' @export


###################################################################################
sample_nb <- function(raster,
                      nb_samples,
                      sample_type,
                      nb_mean,
                      nb_it,
                      buffer,
                      reference,
                      model,
                      mtry,
                      mod.error,
                      last,
                      seed,
                      init.seed,
                      in.memory,
                      save_runs,
                      parallel_mode,
                      max_num_cores,
                      progress_bar = TRUE,
                      session) {
  print(paste(paste("init.samples = ", nb_samples[1]),
              paste("models = ", nb_mean)))
  ###
  n_channel <- length(names(raster))
  ###raster in MEMORY or NOT
  if (in.memory == TRUE &&
      raster::fromDisk(raster) == TRUE) {
    rast <- raster::readAll(raster)
  } else {
    rast <- raster
  }
  rast <- raster::mask(rast, raster::calc(rast, fun = sum))
  ###
  l <- 1        ###6. opt=260
  pbtn1 <- matrix(1, nrow = 1, ncol = 1)
  m <- vector("numeric", length = length(nb_samples))
  layer <- list()
  for (r in nb_samples) {
    ############################################################################
    if (last == T) {
      reference <- rbind(reference, rep(3000, ncol(reference)))
    }
    n <- nb_it
    sample_size <- r
    max_samples_per_class <- sample_size * 5
    if (class(init.seed) == "character" && init.seed == "sample") {
      seed2 <- sample(c(1:1000000), size = nb_mean, replace = F)
    } else {
      seed2 <- init.seed
    }
    oobe <- matrix(NA, nrow = n, ncol = nb_mean)
    models_list <- list()
    points_list <- list()
    dif <- matrix(NA, nrow = nb_mean, ncol = nrow(reference))
    channel <- matrix(NA, nrow = nb_mean, ncol = nrow(reference))
    switch <- matrix(NA, nrow = nb_mean, ncol = nrow(reference))
    if (progress_bar == TRUE) {
        pb <- utils::txtProgressBar(min = 1,
                                    max = nb_mean,
                                    style = 3)
    }

    if (parallel_mode == TRUE) {
      cores = parallel::detectCores( logical = TRUE)
      if (cores > max_num_cores) {
        cores <- max_num_cores
      } else {
        cores <- cores - 1
      }

      res <- parallel::mclapply(
        1:nb_mean,
        model_opt_r,
        raster = rast,
        col_names = names(raster),
        sample_type = sample_type,
        buffer = buffer,
        model = model,
        seed = seed,
        n = n,
        sample_size = sample_size,
        n_channel = n_channel,
        seed2 = seed2,
        mtry = mtry,
        mod.error = mod.error,
        pbtn1 = pbtn1,
        max_samples_per_class = max_samples_per_class,
        mc.cores = cores,
        mc.preschedule = TRUE,
        mc.cleanup = TRUE
      )
      for (k in 1:nb_mean) {
        points_list[[k]] <- res[["k" = k]][["points"]]
        models_list[[k]] <- res[["k" = k]][["models"]]
        oobe[, k] <- res[["k" = k]][["oobe"]][, 1]
      }
    } else {
      for (k in 1:nb_mean) {
        res <- model_opt_r(
          raster = rast,
          col_names = names(raster),
          sample_type = sample_type,
          buffer = buffer,
          model = model,
          seed = seed,
          k = k,
          n = n,
          sample_size = sample_size,
          n_channel = n_channel,
          seed2 = seed2,
          mtry = mtry,
          mod.error = mod.error,
          pbtn1 = pbtn1,
          max_samples_per_class = max_samples_per_class
        )
        points_list[[k]] <- res$points
        models_list[[k]] <- res$models
        oobe[, k] <- res$oobe[, 1]
        if (progress_bar) {
            setTxtProgressBar(pb, k)
        } else {
            shinyWidgets::updateProgressBar(
              session = session,
              id = "pbar",
              value = k,
              total = nb_mean,
              title = paste("Doing sampling for model", k)
            )
        }
      }
    }

    #which_models_null <- which(models == "NULL")
    which_models_null <- vapply(models_list, is.not.null <- function(x){!is.null(x)}, FALSE)
    if (length(models_list) == 0 |
        length(which_models_null[which_models_null == FALSE]) == length(models_list)) {
      remove(points_list)
      remove(models_list)
      out <- list(
        returns = 1,
        index = NULL,
        num_models = 0,
        acc = NULL,
        obj = NULL
      )
      return(out)
    }
    if (length(which_models_null[which_models_null == FALSE]) > 0) {
      models <- models_list[which_models_null]
    } else {
      models <- models_list
    }
    remove(models_list)

    for (jj in 1:nrow(reference)) {
      ref <- jj
      rr = 3
      for (i in 1:length(models)) {
        if (i == 1) {
          dummy <- as.numeric(
            as.character(stats::predict(models[[i]], newdata = reference)))
          if (dummy[ref] != 2) {
            dummy[dummy == 1] <- 3
            dummy[dummy == 2] <- 1
            dummy[dummy == 3] <- 2
            switch[i, jj] <- i
          }
        } else {
          dummy2 <- as.numeric(
            as.character(stats::predict(models[[i]], newdata = reference)))

          if (dummy2[ref] != 2) {
            dummy2[dummy2 == 1] <- 3
            dummy2[dummy2 == 2] <- 1
            dummy2[dummy2 == 3] <- 2
            switch[i, jj] <- i
          }

          dummy_set <- dummy
          dummy <- dummy + dummy2

          dummy3 <- dummy / dummy[ref]
          dif[i, jj] <- dummy3[ref] - max(dummy3[-ref])

          if (i > 2) {
            if (dummy3[ref] - max(dummy3[-ref]) > dif[(rr - 1), jj]) {
              dummy <- dummy
              channel[i, jj] <- i
              dif[(rr - 1), jj] <- dif[i, jj]
            } else {
              dummy <- dummy_set
            }
          }
        }
      }
      m[l] <- max(dif[2,], na.rm = T)
      if (!progress_bar) {
        shinyWidgets::updateProgressBar(
          session = session,
          id = "pbar",
          value = jj,
          total = nrow(reference),
          title = paste("Doing prediction for class", jj)
        )
      }
    }
    index <- which.max(dif[2,])
    ch <- as.numeric(na.omit(channel[, index]))
    if (length(ch) == 0) {
      out <- list(
        returns = 2,
        index = NULL,
        num_models = 0,
        acc = NULL,
        obj = NULL
      )
      return(out)
    }
    acc <- (round(m[l] ^ 2, 2) / 0.25)

    cat("\n")
    print(paste("class=", index, "  difference=", (round(m[l] ^ 2, 2) / 0.25),
                sep = ""))
    l <- l + 1
  }
  if (progress_bar == TRUE) {
    close(pb)
  }

  if (save_runs == TRUE) {
    mod_all <- models
  } else {
    mod_all = list()
  }
  models <- models[ch]
  num_models <- length(models)
  print(paste("n_models =", num_models))
  flush(stdout())
  switch <- switch[ch, index]
  points <- points_list[ch]
  dif <- dif[2,]
  ##############################################################################
  ###Vohersage
  ch <- c(1:length(models))
  switch <- switch

  j <- 1
  for (i in ch) {
    if (j == 1) {
      result1 <- raster::predict(object = raster, model = models[[i]])

      if (is.na(switch[i]) == F) {
        result1 <-
          raster::reclassify(result1, rbind(c(0.5, 1.5, 2), c(1.6, 2.5, 1)))
      }
    } else {
      result1 <-
        raster::stack( result1, raster::predict(object = raster, model = models[[i]]))

      if (is.na(switch[i]) == F) {
        result1[[j]] <-
          raster::reclassify(result1[[j]], rbind(c(0.5, 1.5, 2), c(1.6, 2.5, 1)))
      }
    }
    #print(j)
    j <- j + 1
  }
  ###
  dummy <- raster::brick(result1)
  dummy <- raster::calc(dummy, fun = sum)
  layer[[1]] <- dummy

  if (save_runs == TRUE) {
    obj <- new(
      "Habitat",
      models = models,
      ref_samples = points,
      switch = switch,
      layer = layer,
      mod_all = mod_all,
      class_ind = dif,
      seeds = seed2
    )
  } else {
    obj <- new(
      "Habitat",
      models = list(),
      ref_samples = points,
      switch = switch,
      layer = layer,
      mod_all = list(),
      class_ind = 0,
      seeds = 0
    )
  }
  remove(models)
  remove(points)
  remove(switch)
  remove(layer)
  remove(mod_all)
  remove(dif)
  remove(seed2)
  gc(full = TRUE)

  out <- list(
    returns = 0,
    index = index,
    num_models = num_models,
    acc = acc,
    obj = obj
  )
  return(out)
}