diff --git a/DESCRIPTION b/DESCRIPTION
index 46041320..fe5939f9 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -31,6 +31,7 @@ Imports:
     MASS,
     aws.s3,
     Metrics,
+    scales,
     methods
 License: GPL-3 | file LICENSE
 BugReports: https://github.com/ncsu-landscape-dynamics/rpops/issues
diff --git a/NAMESPACE b/NAMESPACE
index 77da2b05..a84d575e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -5,6 +5,7 @@ export(configuration)
 export(pops)
 export(pops_model)
 export(pops_multirun)
+export(pops_optimize)
 export(quantity_allocation_disagreement)
 export(validate)
 importFrom(MASS,mvrnorm)
@@ -33,9 +34,12 @@ importFrom(parallel,makeCluster)
 importFrom(parallel,stopCluster)
 importFrom(stats,cov)
 importFrom(stats,median)
+importFrom(stats,quantile)
 importFrom(stats,rnorm)
 importFrom(stats,runif)
 importFrom(stats,sd)
+importFrom(stats,setNames)
+importFrom(stats,weighted.mean)
 importFrom(terra,app)
 importFrom(terra,as.matrix)
 importFrom(terra,as.points)
diff --git a/R/pops_optimize.R b/R/pops_optimize.R
new file mode 100644
index 00000000..bb7ef5fe
--- /dev/null
+++ b/R/pops_optimize.R
@@ -0,0 +1,760 @@
+quarantine_distance <-
+  function(quarantine_areas, quarantine_directions) {
+    bboxrast <-
+      terra::subst(terra::trim(quarantine_areas, value = 0), 1, 0)
+    ncols <- terra::ncol(bboxrast)
+    nrows <- terra::nrow(bboxrast)
+    if (is.null(quarantine_directions) ||
+        grepl("N", quarantine_directions)) {
+      bboxrast[terra::cellFromRowCol(bboxrast, 1, 1:ncols)] <- 1
+    }
+    if (is.null(quarantine_directions) ||
+        grepl("S", quarantine_directions)) {
+      bboxrast[terra::cellFromRowCol(bboxrast, nrows, 1:ncols)] <- 1
+    }
+    if (is.null(quarantine_directions) ||
+        grepl("E", quarantine_directions)) {
+      bboxrast[terra::cellFromRowCol(bboxrast, 1:nrows, ncols)] <- 1
+    }
+    if (is.null(quarantine_directions) ||
+        grepl("W", quarantine_directions)) {
+      bboxrast[terra::cellFromRowCol(bboxrast, 1:nrows, 1)] <- 1
+    }
+    return(terra::gridDist(bboxrast, target = 1))
+  }
+
+prior_weight <- function(cost_column,
+                         potential_column,
+                         quarantine_column,
+                         score_weights) {
+  minmax_scale <- function(column) {
+    if (max(column) - min(column) > 0) {
+      return(scales::rescale(column, to = c(0, 1)))
+    } else {
+      return(0)
+    }
+  }
+  cost_norm <- minmax_scale(cost_column)
+  potential_norm <- minmax_scale(potential_column)
+  if (!is.null(quarantine_column)) {
+    quarantine_column <- minmax_scale(quarantine_column)
+    weights <- score_weights / sum(score_weights)
+    return((
+      weights[1] * potential_norm +
+        1 - cost_norm + weights[2] * (1 - quarantine_column)
+    ) / 3)
+  }
+  return((potential_norm + 1 - cost_norm) / 2)
+}
+
+pixel_treatments <- function(points, treatments_raster) {
+  raster <-
+    terra::rasterize(points, treatments_raster, background = 0)
+  return(list(raster = raster, cost = sum(points$cost)))
+}
+
+buffer_treatments <-
+  function(points, treatments_raster, cost_raster) {
+    rasterization_factor <- 10
+    buffers <- terra::buffer(points, width = points$buffer_size)
+    high_res_raster <-
+      terra::disagg(treatments_raster, fact = rasterization_factor)
+    buffers_raster <- terra::rasterize(buffers, high_res_raster)
+    count_raster <- terra::aggregate(
+      buffers_raster,
+      fact = rasterization_factor,
+      fun = function(x) {
+        length(x[!is.na(x)])
+      }
+    )
+    treatments_raster <-
+      count_raster / terra::global(count_raster, "max")[[1]]
+    treatment_cost_raster <- treatments_raster * cost_raster
+    actual_cost <-
+      terra::global(treatment_cost_raster, "sum", na.rm = TRUE)[[1]]
+    return(list(raster = treatments_raster, cost = actual_cost))
+  }
+
+treatments <- function(points, treatments_raster, cost_raster) {
+  if ("buffer_size" %in% names(points)) {
+    return(buffer_treatments(points,
+                             treatments_raster,
+                             cost_raster))
+  }
+  return(pixel_treatments(points, treatments_raster))
+}
+
+
+run_pops <- function(config, treatment_raster = NULL) {
+  if (!is.null(treatment_raster)) {
+    config$treatment_maps <- list(terra::as.matrix(treatment_raster * config$pesticide_efficacy,
+                                                   wide = TRUE))
+  } else {
+    config$treatment_maps <- list(matrix(0, nrow = 1, ncol = 1))
+  }
+
+  cl <- parallel::makeCluster(config$core_count)
+  doParallel::registerDoParallel(cl)
+  infected_stack <-
+    foreach::foreach(
+      i = seq_len(config$number_of_iterations),
+      .combine = c,
+      .packages = c("PoPS", "terra")
+    ) %dopar% {
+      config$random_seed <- round(stats::runif(1, 1, 1000000))
+      data <- PoPS::pops_model(
+        random_seed = config$random_seed,
+        multiple_random_seeds = config$multiple_random_seeds,
+        random_seeds = as.matrix(config$random_seeds[1, ])[1, ],
+        use_lethal_temperature = config$use_lethal_temperature,
+        lethal_temperature = config$lethal_temperature,
+        lethal_temperature_month = config$lethal_temperature_month,
+        use_survival_rates = config$use_survival_rates,
+        survival_rate_month = config$survival_rate_month,
+        survival_rate_day = config$survival_rate_day,
+        infected = config$infected,
+        total_exposed = config$total_exposed,
+        exposed = config$exposed,
+        susceptible = config$susceptible,
+        total_populations = config$total_populations,
+        total_hosts = config$total_hosts,
+        mortality_on = config$mortality_on,
+        mortality_tracker = config$mortality_tracker,
+        mortality = config$mortality,
+        quarantine_areas = config$quarantine_areas,
+        treatment_maps = config$treatment_maps,
+        treatment_dates = config$treatment_dates,
+        pesticide_duration = config$pesticide_duration,
+        resistant = config$resistant,
+        use_movements = config$use_movements,
+        movements = config$movements,
+        movements_dates = config$movements_dates,
+        weather = config$weather,
+        temperature = config$temperature,
+        survival_rates = config$survival_rates,
+        weather_coefficient = config$weather_coefficient,
+        weather_coefficient_sd = config$weather_coefficient_sd,
+        res = config$res,
+        rows_cols = config$rows_cols,
+        time_step = config$time_step,
+        reproductive_rate = config$reproductive_rate,
+        spatial_indices = config$spatial_indices,
+        season_month_start_end = config$season_month_start_end,
+        soil_reservoirs = config$soil_reservoirs,
+        mortality_rate = config$mortality_rate,
+        mortality_time_lag = config$mortality_time_lag,
+        start_date = config$start_date,
+        end_date = config$end_date,
+        treatment_method = config$treatment_method,
+        natural_kernel_type = config$natural_kernel_type,
+        anthropogenic_kernel_type = config$anthropogenic_kernel_type,
+        use_anthropogenic_kernel = config$use_anthropogenic_kernel,
+        percent_natural_dispersal = config$percent_natural_dispersal,
+        natural_distance_scale = config$natural_distance_scale,
+        anthropogenic_distance_scale = config$anthropogenic_distance_scale,
+        natural_dir = config$natural_dir,
+        natural_kappa = config$natural_kappa,
+        anthropogenic_dir = config$anthropogenic_dir,
+        anthropogenic_kappa = config$anthropogenic_kappa,
+        output_frequency = config$output_frequency,
+        output_frequency_n = config$output_frequency_n,
+        quarantine_frequency = config$quarantine_frequency,
+        quarantine_frequency_n = config$quarantine_frequency_n,
+        use_quarantine = config$use_quarantine,
+        quarantine_directions = config$quarantine_directions,
+        spreadrate_frequency = config$spreadrate_frequency,
+        spreadrate_frequency_n = config$spreadrate_frequency_n,
+        mortality_frequency = config$mortality_frequency,
+        mortality_frequency_n = config$mortality_frequency_n,
+        use_spreadrates = config$use_spreadrates,
+        model_type_ = config$model_type,
+        latency_period = config$latency_period,
+        generate_stochasticity = config$generate_stochasticity,
+        establishment_stochasticity = config$establishment_stochasticity,
+        movement_stochasticity = config$movement_stochasticity,
+        dispersal_stochasticity = config$dispersal_stochasticity,
+        establishment_probability = config$establishment_probability,
+        dispersal_percentage = config$dispersal_percentage,
+        use_overpopulation_movements = config$use_overpopulation_movements,
+        overpopulation_percentage = config$overpopulation_percentage,
+        leaving_percentage = config$leaving_percentage,
+        leaving_scale_coefficient = config$leaving_scale_coefficient,
+        bbox = config$bounding_box,
+        network_min_distance = config$network_min_distance,
+        network_max_distance = config$network_max_distance,
+        network_filename = config$network_filename,
+        network_movement = config$network_movement,
+        weather_size = config$weather_size,
+        weather_type = config$weather_type,
+        dispersers_to_soils_percentage = config$dispersers_to_soils_percentage,
+        use_soils = config$use_soils
+      )
+      run <- c()
+      run$infected_area <- data$area_infected
+      run$quarantine_escape_distance <-
+        data$quarantine_escape_distance
+      run
+    }
+  stopCluster(cl)
+
+  area_infected_runs <-
+    infected_stack[seq(1, length(infected_stack), 2)]
+  infected_area <-
+    data.frame(t(rep(0, length(
+      area_infected_runs[[1]]
+    ))))
+  if (config$use_quarantine) {
+    quarantine_escape_distance_runs <-
+      infected_stack[seq(2, length(infected_stack), 2)]
+    quarantine_escape_distances <-
+      data.frame(t(rep(0, length(
+        area_infected_runs[[1]]
+      ))))
+  }
+  for (p in seq_len(length(area_infected_runs))) {
+    infected_area[p, ] <- area_infected_runs[[p]]
+    if (config$use_quarantine) {
+      quarantine_escape_distances[p, ] <-
+        quarantine_escape_distance_runs[[p]]
+    }
+  }
+  results <- list()
+  infected_areas <-
+    round(sapply(infected_area, mean, na.rm = TRUE), digits = 0)
+  results$infected_area <- infected_areas[[length(infected_areas)]]
+  if (config$use_quarantine) {
+    quarantine_escape_distances[is.na(quarantine_escape_distances)] <- 0
+    quarantine_distance <-
+      round(sapply(quarantine_escape_distances, mean), digits = 0)
+    results$quarantine_distance <-
+      quarantine_distance[[length(quarantine_distance)]]
+  }
+  return(results)
+}
+
+pops_init <- function(config) {
+  config$function_name <- "multirun"
+  config$management <- FALSE
+  config <- PoPS::configuration(config)
+
+  if (!is.null(config$failure)) {
+    stop(config$failure)
+  }
+  config$crs <- terra::crs(config$host)
+  config$infected <- config$infected_mean
+  exposed2 <- config$exposed_mean
+  exposed <- config$exposed
+  exposed[[config$latency_period + 1]] <- exposed2
+  config$exposed <- exposed
+  config$host <- config$host_mean
+  susceptible <- config$host - config$infected - exposed2
+  susceptible[susceptible < 0] <- 0
+  config$susceptible <- susceptible
+  config$total_hosts <- config$host
+  config$total_exposed <- exposed2
+  if (config$mortality_on) {
+    mortality_tracker2 <- config$mortality_tracker
+    mortality_tracker2[[length(mortality_tracker2)]] <-
+      config$infected
+    config$mortality_tracker <- mortality_tracker2
+  }
+  return(config)
+}
+
+best_guess <- function(points,
+                       weight_column,
+                       treatments_raster,
+                       cost_raster,
+                       budget,
+                       config) {
+  sorted <-
+    points[order(points[[weight_column]][[1]], decreasing = TRUE), ]
+  candidate <- sorted[cumsum(sorted$cost) <= budget, ]
+
+  treatment <- treatments(candidate,
+                          treatments_raster,
+                          cost_raster)
+  result <- run_pops(config, treatment$raster)
+  return(list(candidate = candidate, result = result))
+}
+
+estimate_baseline <- function(config) {
+  # run with 0 spread to get initial quarantine distance
+  config_no_spread <- config
+  config_no_spread$parameter_means[1] <- 0
+  config_no_spread$reproductive_rate <- 0
+  baseline <- run_pops(config_no_spread)
+  quarantine_distance <- baseline$quarantine_distance
+  baseline <- run_pops(config)
+  baseline$quarantine_distance <- quarantine_distance
+  return(baseline)
+}
+
+estimate_initial_threshold <- function(points,
+                                       weight_column,
+                                       treatments_raster,
+                                       cost_raster,
+                                       budget,
+                                       baseline,
+                                       score_weights,
+                                       config) {
+  runs <- 10
+  results_list <- list()
+  for (run in seq(runs)) {
+    candidate <- sample_candidate(points, weight_column, budget)
+    treatment <- treatments(candidate,
+                            treatments_raster,
+                            cost_raster)
+    results_list[[run]] <- run_pops(config, treatment$raster)
+  }
+  scores <- sapply(results_list, scoring, baseline, score_weights)
+  threshold <- quantile(scores, probs = 0.1)
+  threshold_step <- abs(threshold - quantile(scores, probs = 0.2))
+  return(list(threshold = threshold, threshold_step = threshold_step))
+}
+
+scoring <- function(simulated, baseline, weights = c(1, 1)) {
+  scores <- c(NA, NA)
+
+  if (!is.null(simulated$infected_area)) {
+    scores[1] <- simulated$infected_area / baseline$infected_area
+  }
+  if (!is.null(simulated$quarantine_distance)) {
+    scores[2] <-
+      (baseline$quarantine_distance - simulated$quarantine_distance) /
+      baseline$quarantine_distance
+  }
+  return(weighted.mean(scores, w = weights, na.rm = TRUE))
+}
+
+sample_candidate <- function(points, weight_column, budget) {
+  # to allow sample vector size 1
+  sample.vec <- function(x, ...)
+    x[sample(length(x), ...)]
+  tmp_budget <- 0
+  tries <- 0
+  count <- 1
+  candidate_cats <- c()
+  cats <- points$cat
+  weights <- points[[weight_column]][[1]]
+  costs <- points$cost
+  while (length(cats)) {
+    sample_cat <- sample.vec(cats, 1, prob = weights)
+    selected_index <- which(cats == sample_cat)
+    cost <- costs[[selected_index]]
+    if (tmp_budget + cost <= budget) {
+      tmp_budget <- tmp_budget + cost
+      candidate_cats[count] <- sample_cat
+      count <- count + 1
+      cats <- cats[-selected_index]
+      weights <- weights[-selected_index]
+      costs <- costs[-selected_index]
+    } else {
+      tries <- tries + 1
+      if (tries > 50) {
+        break
+      }
+    }
+  }
+  return(points[points$cat %in% candidate_cats, ])
+}
+
+
+generation <- function(points,
+                       weight_column,
+                       treatments_raster,
+                       cost_raster,
+                       budget,
+                       min_particles,
+                       threshold_percentile,
+                       threshold,
+                       threshold_step,
+                       baseline,
+                       score_weights,
+                       config) {
+  score_list <- numeric(min_particles)
+  particle_count <- 0
+  tested <- 0
+  best <- list()
+  best$simulation <- baseline
+  best$score <- 1
+  new_weights <-
+    setNames(as.list(rep(0, length(points$cat))), points$cat)
+  while (particle_count < min_particles) {
+    tested <- tested + 1
+    candidate <- sample_candidate(points, weight_column, budget)
+    treatment <- treatments(candidate,
+                            treatments_raster,
+                            cost_raster)
+    simulation <- run_pops(config, treatment$raster)
+    score <- scoring(simulation, baseline, score_weights)
+    # success, add particle
+    if (score < threshold) {
+      particle_count <- particle_count + 1
+      for (cat in as.character(candidate$cat)) {
+        new_weights[[cat]] <- new_weights[[cat]] + 1
+      }
+      score_list[particle_count] <- score
+      if (score < best$score) {
+        best$simulation <- simulation
+        best$candidate <- candidate
+        best$score <- score
+        best$cost <- treatment$cost
+        best$treatment <- treatment$raster
+      }
+    }
+    acceptance_rate <- particle_count / tested
+    particles_left <- min_particles - particle_count
+    # print intermediate
+    if (tested %% 10 == 0) {
+      message(
+        "Rate: ",
+        acceptance_rate,
+        ", particles left: ",
+        particles_left,
+        ", best score: ",
+        best$score
+      )
+    }
+    if ((tested == 50) &&
+        (acceptance_rate < 0.1 || acceptance_rate > 0.2)) {
+      if (acceptance_rate < 0.1) {
+        threshold <- threshold + threshold_step
+      } else {
+        threshold <- threshold - threshold_step
+      }
+      message("Adjust threshold to: ", threshold)
+      particle_count <- 0
+      tested <- 0
+      best$score <- 1
+      new_weights <-
+        setNames(as.list(rep(0, length(points$cat))), points$cat)
+    }
+  }
+  new_threshold <-
+    quantile(score_list, probs = threshold_percentile / 100)
+  threshold_step <-
+    abs(new_threshold - quantile(score_list, probs = (threshold_percentile + 10) / 100))
+  output <- list(
+    weights = new_weights,
+    acceptance_rate = acceptance_rate,
+    threshold = new_threshold,
+    threshold_step = threshold_step,
+    best = best
+  )
+  return(output)
+}
+
+filter_particles <-
+  function(points,
+           weight_column,
+           iteration,
+           percentile) {
+    # filter unsuccessful ones
+    filtered <- points[points[[weight_column]] > 0, ]
+    weights_percentile <- quantile(filtered[[weight_column]][[1]],
+                                   percentile / 100,
+                                   names = FALSE)
+    # filter unlikely ones
+    points[points[[weight_column]] <= weights_percentile, weight_column] <-
+      0
+    return(points)
+  }
+
+#' @title Optimize treatments using ABC
+#'
+#' @description Optimizes treatment locations using ABC based on infestation potential and cost.
+#' Treatments are pixel based with the option to specify constant or variable buffer size.
+#' The criteria can be infected area, distance to quarantine boundary or both.
+#'
+#'
+#' @inheritParams pops_multirun
+#' @param infestation_potential_file Raster file with infestation potential.
+#' @param cost_file Raster file with cost of treatment
+#'  (constant value or spatially variable). When buffers are used,
+#'  the cost must represent the cost of the buffer around a pixel.
+#' @param buffer_size_file Raster file with size of buffers
+#' (constant values or spatially variable).
+#' @param min_particles Number of successful combinations per generation.
+#' @param budget Budget to spend. Limits the number of pixels that can be treated.
+#' @param score_weights Weights to compute weighted average of infected_area and
+#' distance to quarantine boundary success metrics. For example, c(1, 0) means only
+#' infected area is used, c(2, 1) means both metrics are averaged with provided weights.
+#' @param filter_percentile Lower value removes fewer pixels from the pool,
+#' resulting in more iterations and lower acceptance rate, but potentially better results.
+#' @param threshold_percentile Determines threshold for next iteration.
+#' @param output_frequency Default is final_step, no need for yearly output
+
+#' @importFrom stats quantile setNames weighted.mean
+#' @useDynLib PoPS, .registration = TRUE
+#' @return results
+#' @export
+#'
+pops_optimize <- function(infestation_potential_file,
+                          cost_file,
+                          buffer_size_file = NULL,
+                          min_particles,
+                          budget,
+                          score_weights = c(1, 1),
+                          filter_percentile = 15,
+                          threshold_percentile = 10,
+                          infected_file,
+                          host_file,
+                          total_populations_file,
+                          parameter_means,
+                          parameter_cov_matrix,
+                          temp = FALSE,
+                          temperature_coefficient_file = "",
+                          precip = FALSE,
+                          precipitation_coefficient_file = "",
+                          model_type = "SI",
+                          latency_period = 0,
+                          time_step = "month",
+                          season_month_start = 1,
+                          season_month_end = 12,
+                          start_date = "2008-01-01",
+                          end_date = "2008-12-31",
+                          use_survival_rates = FALSE,
+                          survival_rate_month = 3,
+                          survival_rate_day = 15,
+                          survival_rates_file = "",
+                          use_lethal_temperature = FALSE,
+                          temperature_file = "",
+                          lethal_temperature = -12.87,
+                          lethal_temperature_month = 1,
+                          mortality_on = FALSE,
+                          mortality_rate = 0,
+                          mortality_time_lag = 0,
+                          mortality_frequency = "year",
+                          mortality_frequency_n = 1,
+                          treatment_dates = c(""),
+                          treatment_method = "ratio",
+                          natural_kernel_type = "cauchy",
+                          anthropogenic_kernel_type = "cauchy",
+                          natural_dir = "NONE",
+                          anthropogenic_dir = "NONE",
+                          number_of_iterations = 100,
+                          number_of_cores = NA,
+                          pesticide_duration = 0,
+                          pesticide_efficacy = 1.0,
+                          random_seed = NULL,
+                          output_frequency = "final_step",
+                          output_frequency_n = 1,
+                          movements_file = "",
+                          use_movements = FALSE,
+                          start_exposed = FALSE,
+                          generate_stochasticity = TRUE,
+                          establishment_stochasticity = TRUE,
+                          movement_stochasticity = TRUE,
+                          dispersal_stochasticity = TRUE,
+                          establishment_probability = 0.5,
+                          dispersal_percentage = 0.99,
+                          quarantine_areas_file = "",
+                          use_quarantine = FALSE,
+                          use_spreadrates = FALSE,
+                          use_overpopulation_movements = FALSE,
+                          overpopulation_percentage = 0,
+                          leaving_percentage = 0,
+                          leaving_scale_coefficient = 1,
+                          exposed_file = "",
+                          mask = NULL,
+                          write_outputs = "None",
+                          output_folder_path = "",
+                          network_filename = "",
+                          network_movement = "walk",
+                          use_initial_condition_uncertainty = FALSE,
+                          use_host_uncertainty = FALSE,
+                          weather_type = "deterministic",
+                          temperature_coefficient_sd_file = "",
+                          precipitation_coefficient_sd_file = "",
+                          dispersers_to_soils_percentage = 0,
+                          quarantine_directions = "",
+                          multiple_random_seeds = FALSE,
+                          file_random_seeds = NULL,
+                          use_soils = FALSE,
+                          soil_starting_pest_file = "",
+                          start_with_soil_populations = FALSE) {
+  # parameters for pops
+  config <- as.list(environment())
+
+  # extract infected points and associated cost and potential
+  infected_raster <- terra::rast(infected_file)
+  infected_points <- terra::as.points(infected_raster)
+  names(infected_points) <- "infected"
+  infected_points <- infected_points[infected_points$infected > 0]
+  # cat
+  infected_points$cat <- seq_len(nrow(infected_points))
+  # cost
+  cost_raster <- terra::rast(cost_file)
+  names(cost_raster) <- "cost"
+  infected_points <-
+    terra::extract(cost_raster, infected_points, bind = TRUE)
+  # infestation potential
+  potential_raster <- terra::rast(infestation_potential_file)
+  names(potential_raster) <- "potential"
+  infected_points <-
+    terra::extract(potential_raster, infected_points, bind = TRUE)
+
+  if (!is.null(buffer_size_file)) {
+    buffer_size_raster <- terra::rast(buffer_size_file)
+    names(buffer_size_raster) <- "buffer_size"
+    infected_points <-
+      terra::extract(buffer_size_raster, infected_points, bind = TRUE)
+  }
+  buffer_size_raster <- NULL
+  if (!is.null(quarantine_areas_file)) {
+    bboxrast <- quarantine_distance(terra::rast(quarantine_areas_file),
+                                    config$quarantine_directions)
+    names(bboxrast) <- "quarantine_distance"
+    infected_points <-
+      terra::extract(bboxrast, infected_points, bind = TRUE)
+  }
+
+  # prior weights
+  iteration <- 1
+  weight_column <- paste0("weight_", iteration)
+  infected_points[[weight_column]] <- prior_weight(
+    infected_points$cost,
+    infected_points$potential,
+    infected_points$quarantine_distance,
+    score_weights
+  )
+
+  # prepare treatment raster
+  treatments_raster <- terra::rast(terra::ext(infected_raster),
+                                   resolution = terra::res(infected_raster))
+  terra::crs(treatments_raster) <- terra::crs(infected_raster)
+  config <- pops_init(config)
+
+  # baseline
+  baseline <- estimate_baseline(config)
+  message("Baseline area:", baseline$infected_area)
+  if (score_weights[2] > 0) {
+    message("Initial distance to quarantine boundary:",
+            baseline$quarantine_distance)
+  }
+
+  # best guess
+  best_guess <- best_guess(infected_points,
+                           weight_column,
+                           treatments_raster,
+                           cost_raster,
+                           budget,
+                           config)
+  message("Best guess infected area:", best_guess$result$infected_area)
+  if (score_weights[2] > 0) {
+    message(
+      "Best guess distance to quarantine boundary:",
+      best_guess$result$quarantine_distance
+    )
+  }
+
+  # initial threshold
+  thresholds <- c()
+
+  initial_threshold <- estimate_initial_threshold(
+    infected_points,
+    weight_column,
+    treatments_raster,
+    cost_raster,
+    budget,
+    baseline,
+    score_weights,
+    config
+  )
+  thresholds[1] <- initial_threshold$threshold
+  threshold_step <- initial_threshold$threshold_step
+
+  filtered_points <- infected_points
+  acceptance_rates <- c()
+  while (TRUE) {
+    message("Iteration ", iteration,
+            ", threshold ", thresholds[length(thresholds)])
+    results <- generation(
+      filtered_points,
+      weight_column,
+      treatments_raster,
+      cost_raster,
+      budget,
+      min_particles,
+      threshold_percentile,
+      thresholds[length(thresholds)],
+      threshold_step,
+      baseline,
+      score_weights,
+      config
+    )
+    acceptance_rates <-
+      append(acceptance_rates, results$acceptance_rate)
+    thresholds <- append(thresholds, results$threshold)
+    threshold_step <- results$threshold_step
+    new_weight_column <- paste0("weight_", iteration + 1)
+    tmp_match <- match(filtered_points$cat, as.integer(names(results$weights)))
+    filtered_points[[new_weight_column]] <- results$weights[tmp_match]
+    # filtering
+    tmp <- filter_particles(filtered_points,
+                            new_weight_column,
+                            iteration,
+                            filter_percentile)
+    before <-
+      filtered_points[filtered_points[[new_weight_column]] > 0, ]
+    after <- tmp[tmp[[new_weight_column]] > 0, ]
+    if (sum(after$cost) >= budget) {
+      message("Filtered ",
+              length(before) - length(after),
+              ", remaining: ",
+              length(after))
+      filtered_points <- tmp
+      weight_column <- new_weight_column
+    } else {
+      break
+    }
+    iteration <- iteration + 1
+  }
+
+  output <- list(
+    best_guess = best_guess,
+    best = results$best,
+    weights = filtered_points,
+    acceptance_rates = acceptance_rates,
+    thresholds = thresholds
+  )
+
+  cat("Baseline area:", baseline$infected_area, "\n")
+  cat("Best guess infected area:",
+      best_guess$result$infected_area,
+      "\n")
+  cat("Optimized infected area:",
+      results$best$simulation$infected_area,
+      "\n")
+
+  if (score_weights[2] > 0) {
+    cat("Initial distance to quarantine boundary:",
+        baseline$quarantine_distance,
+        "\n")
+    cat(
+      "Best guess distance to quarantine boundary:",
+      best_guess$result$quarantine_distance,
+      "\n"
+    )
+    cat(
+      "Optimized distance to quarantine boundary:",
+      results$best$simulation$quarantine_distance,
+      "\n"
+    )
+  }
+  cat("Acceptance rates: ", acceptance_rates, "\n")
+  cat("Thresholds: ", thresholds, "\n")
+  if (file.exists(output_folder_path)) {
+    file_name <- file.path(output_folder_path, "best_candidate.gpkg")
+    terra::writeVector(results$best$candidate, file_name, overwrite = TRUE)
+    file_name <- file.path(output_folder_path, "best_treatment.tif")
+    terra::writeRaster(results$best$treatment, file_name, overwrite = TRUE)
+    file_name <-
+      file.path(output_folder_path, "best_guess_candidate.gpkg")
+    terra::writeVector(best_guess$candidate, file_name, overwrite = TRUE)
+    file_name <- file.path(output_folder_path, "output.rdata")
+    save(output, file = file_name)
+  }
+  return(output)
+}
diff --git a/man/pops_optimize.Rd b/man/pops_optimize.Rd
new file mode 100644
index 00000000..e6c6fc34
--- /dev/null
+++ b/man/pops_optimize.Rd
@@ -0,0 +1,381 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pops_optimize.R
+\name{pops_optimize}
+\alias{pops_optimize}
+\title{Optimize treatments using ABC}
+\usage{
+pops_optimize(
+  infestation_potential_file,
+  cost_file,
+  buffer_size_file = NULL,
+  min_particles,
+  budget,
+  score_weights = c(1, 1),
+  filter_percentile = 15,
+  threshold_percentile = 10,
+  infected_file,
+  host_file,
+  total_populations_file,
+  parameter_means,
+  parameter_cov_matrix,
+  temp = FALSE,
+  temperature_coefficient_file = "",
+  precip = FALSE,
+  precipitation_coefficient_file = "",
+  model_type = "SI",
+  latency_period = 0,
+  time_step = "month",
+  season_month_start = 1,
+  season_month_end = 12,
+  start_date = "2008-01-01",
+  end_date = "2008-12-31",
+  use_survival_rates = FALSE,
+  survival_rate_month = 3,
+  survival_rate_day = 15,
+  survival_rates_file = "",
+  use_lethal_temperature = FALSE,
+  temperature_file = "",
+  lethal_temperature = -12.87,
+  lethal_temperature_month = 1,
+  mortality_on = FALSE,
+  mortality_rate = 0,
+  mortality_time_lag = 0,
+  mortality_frequency = "year",
+  mortality_frequency_n = 1,
+  treatment_dates = c(""),
+  treatment_method = "ratio",
+  natural_kernel_type = "cauchy",
+  anthropogenic_kernel_type = "cauchy",
+  natural_dir = "NONE",
+  anthropogenic_dir = "NONE",
+  number_of_iterations = 100,
+  number_of_cores = NA,
+  pesticide_duration = 0,
+  pesticide_efficacy = 1,
+  random_seed = NULL,
+  output_frequency = "final_step",
+  output_frequency_n = 1,
+  movements_file = "",
+  use_movements = FALSE,
+  start_exposed = FALSE,
+  generate_stochasticity = TRUE,
+  establishment_stochasticity = TRUE,
+  movement_stochasticity = TRUE,
+  dispersal_stochasticity = TRUE,
+  establishment_probability = 0.5,
+  dispersal_percentage = 0.99,
+  quarantine_areas_file = "",
+  use_quarantine = FALSE,
+  use_spreadrates = FALSE,
+  use_overpopulation_movements = FALSE,
+  overpopulation_percentage = 0,
+  leaving_percentage = 0,
+  leaving_scale_coefficient = 1,
+  exposed_file = "",
+  mask = NULL,
+  write_outputs = "None",
+  output_folder_path = "",
+  network_filename = "",
+  network_movement = "walk",
+  use_initial_condition_uncertainty = FALSE,
+  use_host_uncertainty = FALSE,
+  weather_type = "deterministic",
+  temperature_coefficient_sd_file = "",
+  precipitation_coefficient_sd_file = "",
+  dispersers_to_soils_percentage = 0,
+  quarantine_directions = "",
+  multiple_random_seeds = FALSE,
+  file_random_seeds = NULL,
+  use_soils = FALSE,
+  soil_starting_pest_file = "",
+  start_with_soil_populations = FALSE
+)
+}
+\arguments{
+\item{infestation_potential_file}{Raster file with infestation potential.}
+
+\item{cost_file}{Raster file with cost of treatment
+(constant value or spatially variable). When buffers are used,
+the cost must represent the cost of the buffer around a pixel.}
+
+\item{buffer_size_file}{Raster file with size of buffers
+(constant values or spatially variable).}
+
+\item{min_particles}{Number of successful combinations per generation.}
+
+\item{budget}{Budget to spend. Limits the number of pixels that can be treated.}
+
+\item{score_weights}{Weights to compute weighted average of infected_area and
+distance to quarantine boundary success metrics. For example, c(1, 0) means only
+infected area is used, c(2, 1) means both metrics are averaged with provided weights.}
+
+\item{filter_percentile}{Lower value removes fewer pixels from the pool,
+resulting in more iterations and lower acceptance rate, but potentially better results.}
+
+\item{threshold_percentile}{Determines threshold for next iteration.}
+
+\item{infected_file}{Raster file with initial infections. Units for infections are based on data
+availability and the way the units used for your host file is created (e.g. percent area, # of
+hosts per cell, etc.).}
+
+\item{host_file}{path to raster files with number of hosts and standard deviation on those
+estimates can be based in 3 formats (a single file with number of hosts, a single file with 2
+layers number of hosts and standard deviation, or two files 1 with number of hosts and the other
+with standard deviation of those estimates). The units for this can be of many formats the two
+most common that we use are either percent area (0 to 100) or # of hosts in the cell. Usually
+depends on data available and estimation methods.}
+
+\item{total_populations_file}{path to raster file with number of total populations of all hosts
+and non-hosts. This depends on how your host data is set up. If host is percent area then this
+should be a raster with values that are 100 anywhere with host. If host file is # of hosts in a
+cell then this should be a raster with values that are the max of the host raster any where the
+# of hosts is greater than 0.}
+
+\item{parameter_means}{A vector of the means of the model parameters (reproductive_rate,
+natural_dispersal_distance, percent_natural_dispersal, anthropogenic_dispersal_distance, natural
+kappa, anthropogenic kappa, network_min_distance, and network_max_distance). 1x8 vector.}
+
+\item{parameter_cov_matrix}{A covariance matrix from the previous years posterior parameter
+estimation ordered from (reproductive_rate, natural_dispersal_distance,
+percent_natural_dispersal, anthropogenic_dispersal_distance, natural kappa, anthropogenic kappa,
+network_min_distance, and network_max_distance) Should be 8x8 matrix.}
+
+\item{temp}{boolean that allows the use of temperature coefficients to modify spread
+(TRUE or FALSE)}
+
+\item{temperature_coefficient_file}{path to raster file with temperature coefficient data for the
+timestep and and time period specified (e.g. if timestep = week and start_date = 2017_01_01 and
+end_date = 2019_12_31 this file would have 52 * 3 bands = 156 bands with data being weekly
+precipitation coefficients). We convert raw precipitation values to coefficients that affect the
+reproduction and survival of the pest all values in the raster are between 0 and 1.}
+
+\item{precip}{boolean that allows the use of precipitation coefficients to modify spread
+(TRUE or FALSE)}
+
+\item{precipitation_coefficient_file}{Raster file with precipitation coefficient data for the
+timestep and time period specified (e.g. if timestep = week and start_date = 2017_01_01 and
+end_date = 2019_12_31 this file would have 52 * 3 bands = 156 bands with data being weekly
+precipitation coefficients). We convert raw precipitation values to coefficients that affect the
+reproduction and survival of the pest all values in the raster are between 0 and 1.}
+
+\item{model_type}{What type of model most represents your system. Options are "SEI"
+(Susceptible - Exposed - Infected/Infested) or "SI" (Susceptible - Infected/Infested). Default
+value is "SI".}
+
+\item{latency_period}{How many times steps does it take to for exposed populations become
+infected/infested. This is an integer value and must be greater than 0 if model type is SEI.}
+
+\item{time_step}{How often should spread occur options: ('day', 'week', 'month').}
+
+\item{season_month_start}{When does spread first start occurring in the year for your pest or
+pathogen (integer value between 1 and 12)}
+
+\item{season_month_end}{When does spread end during the year for your pest
+or pathogen (integer value between 1 and 12)}
+
+\item{start_date}{Date to start the simulation with format ('YYYY_MM_DD')}
+
+\item{end_date}{Date to end the simulation with format ('YYYY_MM_DD')}
+
+\item{use_survival_rates}{Boolean to indicate if the model will use survival rates to limit the
+survival or emergence of overwintering generations.}
+
+\item{survival_rate_month}{What month do over wintering generations emerge. We suggest using the
+month before for this parameter as it is when the survival rates raster will be applied.}
+
+\item{survival_rate_day}{What day should the survival rates be applied}
+
+\item{survival_rates_file}{Raster file with survival rates from 0 to 1 representing the
+percentage of emergence for a cell.}
+
+\item{use_lethal_temperature}{A boolean to answer the question: does your pest or pathogen have
+a temperature at which it cannot survive? (TRUE or FALSE)}
+
+\item{temperature_file}{Path to raster file with temperature data for minimum temperature}
+
+\item{lethal_temperature}{The temperature in degrees C at which lethal temperature related
+mortality occurs for your pest or pathogen (-50 to 60)}
+
+\item{lethal_temperature_month}{The month in which lethal temperature related mortality occurs
+for your pest or pathogen integer value between 1 and 12}
+
+\item{mortality_on}{Boolean to turn host mortality on and off (TRUE or FALSE)}
+
+\item{mortality_rate}{Rate at which mortality occurs value between 0 and 1}
+
+\item{mortality_time_lag}{Time lag from infection until mortality can occur in time steps
+integer >= 1}
+
+\item{mortality_frequency}{Sets the frequency of mortality calculations occur either ('year',
+'month', week', 'day', 'time step', or 'every_n_steps')}
+
+\item{mortality_frequency_n}{Sets number of units from mortality_frequency in which to run the
+mortality calculation if mortality_frequency is 'every_n_steps'. Must be an integer >= 1.}
+
+\item{treatment_dates}{Dates in which to apply treatment list with format ('YYYY_MM_DD')
+(needs to be the same length as treatment_file and pesticide_duration)}
+
+\item{treatment_method}{What method to use when applying treatment one of ("ratio" or "all
+infected"). ratio removes a portion of all infected and susceptibles, all infected removes all
+infected a portion of susceptibles.}
+
+\item{natural_kernel_type}{What type of dispersal kernel should be used for natural dispersal.
+Current dispersal kernel options are ('Cauchy', 'exponential', 'uniform',
+'deterministic neighbor','power law', 'hyperbolic secant', 'gamma', 'weibull', 'logistic')}
+
+\item{anthropogenic_kernel_type}{What type of dispersal kernel should be used for anthropogenic
+dispersal. Current dispersal kernel options are ('cauchy', 'exponential', 'uniform',
+'deterministic neighbor','power law', 'hyperbolic secant', 'gamma', 'weibull', 'logistic',
+'network')}
+
+\item{natural_dir}{Sets the predominate direction of natural dispersal usually due to wind values
+('N', 'NW', 'W', 'SW', 'S', 'SE', 'E', 'NE', 'NONE')}
+
+\item{anthropogenic_dir}{Sets the predominate direction of anthropogenic dispersal usually due
+to human movement typically over long distances (e.g. nursery trade, movement of firewood, etc..)
+('N', 'NW', 'W', 'SW', 'S', 'SE', 'E', 'NE', 'NONE')}
+
+\item{number_of_iterations}{how many iterations do you want to run to allow the calibration to
+converge at least 10}
+
+\item{number_of_cores}{enter how many cores you want to use (default = NA). If not set uses the
+# of CPU cores - 1. must be an integer >= 1}
+
+\item{pesticide_duration}{How long does the pesticide (herbicide, vaccine, etc..) last before the
+host is susceptible again. If value is 0 treatment is a culling (i.e. host removal) not a
+pesticide treatment. (needs to be the same length as treatment_dates and treatment_file)}
+
+\item{pesticide_efficacy}{How effective is the pesticide at preventing the disease or killing the
+pest (if this is 0.70 then when applied it successfully treats 70 percent of the plants or
+animals).}
+
+\item{random_seed}{Sets the random seed for the simulation used for reproducibility}
+
+\item{output_frequency}{Default is final_step, no need for yearly output}
+
+\item{output_frequency_n}{Sets number of units from output_frequency in which to export model
+results if mortality_frequency is 'every_n_steps'. Must be an integer >= 1.}
+
+\item{movements_file}{This is a csv file with columns lon_from, lat_from, lon_to, lat_to, number
+of animals, and date.}
+
+\item{use_movements}{This is a boolean to turn on use of the movement module.}
+
+\item{start_exposed}{Do your initial conditions start as exposed or infected (only used if
+model_type is "SEI"). Default False. If this is TRUE need to have both an infected_file (this
+can be a raster of all 0's) and exposed_file}
+
+\item{generate_stochasticity}{Boolean to indicate whether to use stochasticity in reproductive
+functions default is TRUE}
+
+\item{establishment_stochasticity}{Boolean to indicate whether to use stochasticity in
+establishment functions default is TRUE}
+
+\item{movement_stochasticity}{Boolean to indicate whether to use stochasticity in movement
+functions default is TRUE}
+
+\item{dispersal_stochasticity}{Boolean to indicate whether to use a stochasticity in the
+dispersal kernel default is TRUE}
+
+\item{establishment_probability}{Threshold to determine establishment if
+establishment_stochasticity is FALSE (range 0 to 1, default = 0.5)}
+
+\item{dispersal_percentage}{Percentage of dispersal used to calculate the bounding box for
+deterministic dispersal}
+
+\item{quarantine_areas_file}{Path to raster file with quarantine boundaries used in calculating
+likelihood of quarantine escape if use_quarantine is TRUE}
+
+\item{use_quarantine}{Boolean to indicate whether or not there is a quarantine area if TRUE must
+pass in a raster file indicating the quarantine areas (default = FALSE)}
+
+\item{use_spreadrates}{Boolean to indicate whether or not to calculate spread rates}
+
+\item{use_overpopulation_movements}{Boolean to indicate whether to use the overpopulation pest
+movement module (driven by the natural kernel with its scale parameter modified by a coefficient)}
+
+\item{overpopulation_percentage}{Percentage of occupied hosts when the cell is considered to be
+overpopulated}
+
+\item{leaving_percentage}{Percentage of pests leaving an overpopulated cell}
+
+\item{leaving_scale_coefficient}{Coefficient to multiply scale parameter of the natural kernel
+(if applicable)}
+
+\item{exposed_file}{A file with the exposed for the current}
+
+\item{mask}{Raster file used to provide a mask to remove 0's that are not true negatives from
+comparisons (e.g. mask out lakes and oceans from statics if modeling terrestrial species). This
+can also be used to mask out areas that can't be managed in the auto_manage function.}
+
+\item{write_outputs}{Either c("summary_outputs", "all_simulations", or "None"). If not
+"None" output folder path must be provided.}
+
+\item{output_folder_path}{this is the full path with either / or \\ (e.g.,
+"C:/user_name/desktop/pops_sod_2020_2023/outputs/")}
+
+\item{network_filename}{The entire file path for the network file. Used if
+anthropogenic_kernel_type = 'network'.}
+
+\item{network_movement}{What movement type do you want to use in the network kernel either
+"walk", "jump", or "teleport". "walk" allows dispersing units to leave the network at any cell
+along the edge. "jump" automatically moves to the nearest node when moving through the network.
+"teleport" moves from node to node most likely used for airport and seaport networks.}
+
+\item{use_initial_condition_uncertainty}{Boolean to indicate whether or not to propagate and
+partition uncertainty from initial conditions. If TRUE the infected_file needs to have 2 layers
+one with the mean value and one with the standard deviation. If an SEI model is used the
+exposed_file needs to have 2 layers one with the mean value and one with the standard
+deviation}
+
+\item{use_host_uncertainty}{Boolean to indicate whether or not to propagate and partition
+uncertainty from host data. If TRUE the host_file needs to have 2 layers one with the mean value
+and one with the standard deviation.}
+
+\item{weather_type}{string indicating how the weather data is passed in  either as a mean and
+standard deviation to represent uncertainty ("probabilistic") or as a time series
+("deterministic")}
+
+\item{temperature_coefficient_sd_file}{Raster file with temperature coefficient standard
+deviation data for the timestep and time period specified (e.g. if timestep = week this file
+would have 52 bands with data being weekly temperature coefficient standard deviations). We
+convert raw temperature values to coefficients that affect the reproduction and survival of
+the pest all values in the raster are between 0 and 1.}
+
+\item{precipitation_coefficient_sd_file}{Raster file with precipitation coefficient standard
+deviation data for the timestep and time period specified (e.g. if timestep = week this file
+would have 52 bands with data being weekly precipitation coefficient standard deviations). We
+convert raw precipitation values to coefficients that affect the reproduction and survival of
+the pest all values in the raster are between 0 and 1.}
+
+\item{dispersers_to_soils_percentage}{Range from 0 to 1 representing the percentage of dispersers
+that fall to the soil and survive.}
+
+\item{quarantine_directions}{String with comma separated directions to include in the quarantine
+direction analysis, e.g., 'N,E'. By default all directions (N, S, E, W) are considered}
+
+\item{multiple_random_seeds}{Boolean to indicate if the model should use multiple random seeds
+(allows for performing uncertainty partitioning) or a single random seed (backwards
+compatibility option). Default is FALSE.}
+
+\item{file_random_seeds}{A file path to the file with the .csv file containing random_seeds
+table. Use if you are trying to recreate an exact analysis otherwise we suggest leaving the
+default. Default is Null which draws the seed numbers for each.}
+
+\item{use_soils}{Boolean to indicate if pests establish in the soil and spread out from there.
+Typically used for soil borne pathogens.}
+
+\item{soil_starting_pest_file}{path to the raster file with the starting amount of pest or
+pathogen.}
+
+\item{start_with_soil_populations}{Boolean to indicate whether to use a starting soil pest or
+pathogen population if TRUE then soil_starting_pest_file is required.}
+}
+\value{
+results
+}
+\description{
+Optimizes treatment locations using ABC based on infestation potential and cost.
+Treatments are pixel based with the option to specify constant or variable buffer size.
+The criteria can be infected area, distance to quarantine boundary or both.
+}