From b50b79000274e2ee740477e1381bdf9ac8807705 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Tue, 17 Dec 2019 12:41:58 -0500
Subject: [PATCH 01/25] added extra checks for number of outputs matching
 layers of validation data

---
 R/validate.R | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/R/validate.R b/R/validate.R
index daebb36e..cc319f76 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -147,6 +147,18 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   
   number_of_years <- ceiling(time_length(duration, "year"))
   
+  if (output_frequency == "week") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "week"))
+  } else if (output_frequency == "month") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "month"))
+  } else if (output_frequency == "day") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "day"))
+  } else if (output_frequency == "year") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "year"))
+  } else if (output_frequency == "time_step") {
+    number_of_outputs <- number_of_time_steps
+  }
+  
   infected <- raster::raster(infected_file)
   infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
@@ -375,7 +387,7 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   infection_years <- raster::reclassify(infection_years, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   
   num_layers_infected_years <- raster::nlayers(infection_years)
-  if (num_layers_infected_years < number_of_time_steps) {
+  if (num_layers_infected_years < number_of_outputs) {
     return(paste("The infection years file must have enough layers to match the number of outputs from the model. The number of layers of your infected year file is", num_layers_infected_years, "and the number of outputs is", number_of_time_steps))
   }
   

From cd05a483f8afe132a4aca6d2a55d5d7fbe64579f Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Wed, 18 Dec 2019 15:46:32 -0500
Subject: [PATCH 02/25] fixed small bug in is_last_week_of_month function

---
 inst/include/date.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/inst/include/date.hpp b/inst/include/date.hpp
index 7cc2dd88..432848ce 100644
--- a/inst/include/date.hpp
+++ b/inst/include/date.hpp
@@ -162,11 +162,11 @@ bool Date::is_last_day_of_year()
 bool Date::is_last_week_of_month()
 {
   if (this->is_leap_year()){
-    if ((day_ + 7) == day_in_month[1][month_])
+    if ((day_ + 7) >= day_in_month[1][month_])
       return true;
     return false;
   } else {
-    if ((day_ + 7)== day_in_month[0][month_])
+    if ((day_ + 7) >= day_in_month[0][month_])
       return true;
     return false;
   }

From 6e4001f85fad2de59b48607ad0d61965aa6f2bef Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Wed, 18 Dec 2019 15:47:10 -0500
Subject: [PATCH 03/25] updated to bayesion updating of parameters and added
 anthropogenic distance and percent natural dispersal as calibrated parameters

---
 R/calibrate.R    | 389 ++++++++++++++++++++++++++++++++++++++++++++---
 man/calibrate.Rd |  56 ++++---
 2 files changed, 391 insertions(+), 54 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 9abb1b00..4f03d1f6 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -9,11 +9,13 @@
 #' @inheritParams pops
 #' @param infected_years_file Raster file with years of initial infection/infestation as individual locations of a pest or pathogen
 #' @param num_iterations how many iterations do you want to run to allow the calibration to converge (recommend a minimum of at least 100,000 but preferably 1 million).
-#' @param start_reproductive_rate starting reproductive rate for MCMC calibration (affects how quickly a series converges) numeric value > 0
-#' @param start_natural_distance_scale starting short distance scale parameter for MCMC calibration (affects how quickly a series converges) numeric value > 0
+#' @param number_of_observations the number of observations used for this calibartion 
+#' @param prior_number_of_observations the number of total observations from previous calibrations used to weight the posterior distributions (if this is a new calibration this value takes the form of a prior weight (0 - 1))
+#' @param prior_reproductive_rate the prior reproductive rate for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1
+#' @param prior_natural_distance_scale the prior natural distance scale for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1
+#' @param prior_percent_natural_dispersal the prior percent natural distance for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1
+#' @param prior_anthropogenic_distance_scale the prior anthropogenic distance scale for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1
 #' @param number_of_cores number of cores to use for calibration (defaults to the number of cores available on the machine) integer value >= 1
-#' @param sd_reproductive_rate starting standard deviation for reproductive rate for MCMC calibration (can affect how quickly and if a series converges) numeric value > 0
-#' @param sd_natural_distance_scale starting standard deviation for short distance scale for MCMC calibration (can affect how quickly and if a series converges) numeric value > 0
 #' @param success_metric Choose which success metric to use for calibration. Choices are "quantity", "quantity and configuration", "residual error" and "odds ratio". Default is "quantity"
 #' @param 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). 
 #'
@@ -33,12 +35,16 @@
 #' \dontrun{
 #' }
 
-calibrate <- function(infected_years_file, num_iterations, start_reproductive_rate, number_of_cores = NA,
-                      start_natural_distance_scale, sd_reproductive_rate, sd_natural_distance_scale,
+calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA,
+                      number_of_observations, prior_number_of_observations,
+                      prior_reproductive_rate,
+                      prior_natural_distance_scale,
+                      prior_percent_natural_dispersal = c(1.0,0), 
+                      prior_anthropogenic_distance_scale = c(1000,0),
                       infected_file, host_file, total_plants_file, 
                       temp = FALSE, temperature_coefficient_file = "", 
                       precip = FALSE, precipitation_coefficient_file = "", 
-                      time_step = "month", reproductive_rate = 3.0,
+                      time_step = "month", 
                       season_month_start = 1, season_month_end = 12, 
                       start_date = '2008-01-01', end_date = '2008-12-31', 
                       use_lethal_temperature = FALSE, temperature_file = "",
@@ -46,13 +52,11 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
                       mortality_on = FALSE, mortality_rate = 0, mortality_time_lag = 0, 
                       management = FALSE, treatment_dates = c(0), treatments_file = "",
                       treatment_method = "ratio",
-                      percent_natural_dispersal = 1.0,
                       natural_kernel_type = "cauchy", anthropogenic_kernel_type = "cauchy",
-                      natural_distance_scale = 21, anthropogenic_distance_scale = 0.0,
                       natural_dir = "NONE", natural_kappa = 0, 
                       anthropogenic_dir = "NONE", anthropogenic_kappa = 0,
                       pesticide_duration = c(0), pesticide_efficacy = 1.0,
-                      mask = NULL, success_metric = "quantity", output_frequency = "year"){ 
+                      mask = NULL, success_metric = "quantity", output_frequency = "year") { 
   
   if (success_metric == "quantity") {
     configuration = FALSE
@@ -74,7 +78,7 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
     return("Infected file does not exist") 
   }
   
-  if (!(extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
+  if (!(raster::extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
     return("Infected file is not one of '.grd', '.tif', '.img'")
   }
   
@@ -82,7 +86,7 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
     return("Host file does not exist") 
   }
   
-  if (!(extension(host_file) %in% c(".grd", ".tif", ".img"))) {
+  if (!(raster::extension(host_file) %in% c(".grd", ".tif", ".img"))) {
     return("Host file is not one of '.grd', '.tif', '.img'")
   }
   
@@ -90,7 +94,7 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
     return("Total plants file does not exist") 
   }
   
-  if (!(extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
+  if (!(raster::extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
@@ -130,6 +134,114 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
   
   number_of_years <- ceiling(lubridate::time_length(duration, "year"))
   
+  if (output_frequency == "week") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "week"))
+  } else if (output_frequency == "month") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "month"))
+  } else if (output_frequency == "day") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "day"))
+  } else if (output_frequency == "year") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "year"))
+  } else if (output_frequency == "time_step") {
+    number_of_outputs <- number_of_time_steps
+  }
+  
+  ## Setup for reproductive rate to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
+  if (class(prior_reproductive_rate) == "numeric" && length(prior_reproductive_rate) == 2) {
+    prior_reproductive_rate <- matrix(prior_reproductive_rate, ncol = 2)
+  } 
+  
+  if (class(prior_reproductive_rate) %in% c("matrix", "data.frame") && ncol(prior_reproductive_rate) == 2) {
+    if (class(prior_reproductive_rate) == "matrix" && nrow(prior_reproductive_rate) == 1) {
+      start_reproductive_rate <- prior_reproductive_rate[1]
+      sd_reproductive_rate <- prior_reproductive_rate[2]
+    } else if (class(prior_reproductive_rate) == "data.frame" && nrow(prior_reproductive_rate) == 1) {
+      start_reproductive_rate <- prior_reproductive_rate[[1]]
+      sd_reproductive_rate <- 0
+    } else if (class(prior_reproductive_rate) %in% c("matrix", "data.frame") && nrow(prior_reproductive_rate) > 1) {
+      names(prior_reproductive_rate) <- c('var', 'prob')
+      start_reproductive_rate <- prior_reproductive_rate$var[prior_reproductive_rate$prob == max(prior_reproductive_rate$prob)]
+      if(length(start_reproductive_rate) > 1) {
+        start_reproductive_rate <- mean(start_reproductive_rate)
+      }
+      sd_reproductive_rate <- sd(prior_reproductive_rate$var)
+    }
+  } else {
+    return("Incorrect format for prior reproductive rate")
+  }
+  
+  ## Setup for natural dispersal scale to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
+  if (class(prior_natural_distance_scale) == "numeric" && length(prior_natural_distance_scale) == 2) {
+    prior_natural_distance_scale <- matrix(prior_natural_distance_scale, ncol = 2)
+  } 
+  
+  if (class(prior_natural_distance_scale) %in% c("matrix", "data.frame") && ncol(prior_natural_distance_scale) == 2) {
+    if (class(prior_natural_distance_scale) == "matrix" && nrow(prior_natural_distance_scale) == 1) {
+      start_natural_distance_scale <- prior_natural_distance_scale[1]
+      sd_natural_distance_scale <- prior_natural_distance_scale[2]
+    } else if (class(prior_natural_distance_scale) == "data.frame" && nrow(prior_natural_distance_scale) == 1) {
+      start_natural_distance_scale <- prior_natural_distance_scale[[1]]
+      sd_natural_distance_scale <- 0
+    } else if (class(prior_natural_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_natural_distance_scale) > 1) {
+      names(prior_natural_distance_scale) <- c('var', 'prob')
+      start_natural_distance_scale <- prior_natural_distance_scale$var[prior_natural_distance_scale$prob == max(prior_natural_distance_scale$prob)]
+      if(length(start_natural_distance_scale) > 1) {
+        start_natural_distance_scale <- mean(start_natural_distance_scale)
+      }
+      sd_natural_distance_scale <- sd(prior_natural_distance_scale$var)
+    }
+  } else {
+    return("Incorrect format for prior natural distance scale")
+  }
+  
+  ## Setup for anthropogenic dispersal scale to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
+  if (class(prior_anthropogenic_distance_scale) == "numeric" && length(prior_anthropogenic_distance_scale) == 2) {
+    prior_anthropogenic_distance_scale <- matrix(prior_anthropogenic_distance_scale, ncol = 2)
+  } 
+  
+  if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && ncol(prior_anthropogenic_distance_scale) == 2) {
+    if (class(prior_anthropogenic_distance_scale) == "matrix" && nrow(prior_anthropogenic_distance_scale) == 1) {
+      start_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale[1]
+      sd_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale[2]
+    } else if (class(prior_anthropogenic_distance_scale) == "data.frame" && nrow(prior_anthropogenic_distance_scale) == 1) {
+      start_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale[[1]]
+      sd_anthropogenic_distance_scale <- 0
+    } else if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_anthropogenic_distance_scale) > 1) {
+      names(prior_anthropogenic_distance_scale) <- c('var', 'prob')
+      start_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale$var[prior_anthropogenic_distance_scale$prob == max(prior_anthropogenic_distance_scale$prob)]
+      if(length(start_anthropogenic_distance_scale) > 1) {
+        start_anthropogenic_distance_scale <- mean(start_anthropogenic_distance_scale)
+      }
+      sd_anthropogenic_distance_scale <- sd(prior_anthropogenic_distance_scale$var)
+    }
+  } else {
+    return("Incorrect format for prior athropogenic distance scale")
+  }
+    
+  ## Setup for percent natural distance to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
+  if (class(prior_percent_natural_dispersal) == "numeric" && length(prior_percent_natural_dispersal) == 2) {
+    prior_percent_natural_dispersal <- matrix(prior_percent_natural_dispersal, ncol = 2)
+  } 
+    
+  if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && ncol(prior_percent_natural_dispersal) == 2) {
+    if (class(prior_percent_natural_dispersal) == "matrix" && nrow(prior_percent_natural_dispersal) == 1) {
+      start_percent_natural_dispersal <- prior_percent_natural_dispersal[1]
+      sd_percent_natural_dispersal <- prior_percent_natural_dispersal[2]
+    } else if (class(prior_percent_natural_dispersal) == "data.frame" && nrow(prior_percent_natural_dispersal) == 1) {
+      start_percent_natural_dispersal <- prior_percent_natural_dispersal[[1]]
+      sd_percent_natural_dispersal <- 0
+    } else if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && nrow(prior_percent_natural_dispersal) > 1) {
+      names(prior_percent_natural_dispersal) <- c('var', 'prob')
+      start_percent_natural_dispersal <- prior_percent_natural_dispersal$var[prior_percent_natural_dispersal$prob == max(prior_percent_natural_dispersal$prob)]
+      if(length(start_percent_natural_dispersal) > 1) {
+        start_percent_natural_dispersal <- mean(start_percent_natural_dispersal)
+      }
+      sd_percent_natural_dispersal <- sd(prior_percent_natural_dispersal$var)
+    }
+  } else {
+    return("Incorrect format for prior percent natural distance")
+  }
+  
   infected <- raster::raster(infected_file)
   infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
@@ -322,9 +434,9 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
     treatment_maps <- list(raster::as.matrix(treatment_map))
   }
   
-  if(percent_natural_dispersal == 1.0) {
+  if(start_percent_natural_dispersal == 1.0) {
     use_anthropogenic_kernel = FALSE
-  } else if (percent_natural_dispersal < 1.0  && percent_natural_dispersal >= 0.0) {
+  } else if (start_percent_natural_dispersal < 1.0  && start_percent_natural_dispersal >= 0.0) {
     use_anthropogenic_kernel = TRUE
   } else {
     return("Percent natural dispersal must be between 0.0 and 1.0")
@@ -361,12 +473,12 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
   infection_years <- raster::reclassify(infection_years, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   num_layers_infected_years <- raster::nlayers(infection_years)
   
-  if (num_layers_infected_years < number_of_time_steps) {
+  if (num_layers_infected_years < number_of_outputs) {
     return(paste("The infection years file must have enough layers to match the number of outputs from the model. The number of layers of your infected year file is", num_layers_infected_years, "and the number of outputs is", number_of_time_steps))
   }
   
   ## set the parameter function to only need the parameters that chanage
-  param_func <- function(reproductive_rate, natural_distance_scale) {
+  param_func <- function(reproductive_rate, natural_distance_scale, anthropogenic_distance_scale, percent_natural_dispersal) {
     random_seed <- round(runif(1, 1, 1000000))
     data <- pops_model(random_seed = random_seed, 
                        use_lethal_temperature = use_lethal_temperature, 
@@ -400,7 +512,7 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
     return(data)
   }
   
-  data <- param_func(start_reproductive_rate, start_natural_distance_scale)
+  data <- param_func(start_reproductive_rate, start_natural_distance_scale, start_anthropogenic_distance_scale, start_percent_natural_dispersal)
   
   ## set up comparison
   
@@ -412,7 +524,7 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
   }
   
   ## save current state of the system
-  current <- best <- data.frame(t(all_disagreement), reproductive_rate = start_reproductive_rate, natural_distance_scale = start_natural_distance_scale)
+  current <- best <- data.frame(t(all_disagreement), reproductive_rate = start_reproductive_rate, natural_distance_scale = start_natural_distance_scale, anthropogenic_distance_scale = start_anthropogenic_distance_scale, percent_natural_dispersal = start_percent_natural_dispersal)
   
   ## create parallel environment
   if (is.na(number_of_cores) || number_of_cores > parallel::detectCores()) {
@@ -433,11 +545,21 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
     proposed_natural_distance_scale <- round(rnorm(1, mean = best$natural_distance_scale, sd = sd_natural_distance_scale), digits = 0)
   }
   
+  proposed_anthropogenic_distance_scale <- 0
+  while (proposed_anthropogenic_distance_scale <= 0) {
+    proposed_anthropogenic_distance_scale <- round(rnorm(1, mean = best$anthropogenic_distance_scale, sd = sd_anthropogenic_distance_scale)/10, digits = 0) * 10
+  }
+  
+  proposed_percent_natural_dispersal <- 0
+  while (proposed_percent_natural_dispersal <= 0 || proposed_percent_natural_dispersal > 1.000) {
+    proposed_percent_natural_dispersal <- round(rnorm(1, mean = best$percent_natural_dispersal, sd = sd_percent_natural_dispersal), digits = 3)
+  }
+  
   params <- foreach(icount(num_iterations), .combine = rbind, .packages = c("raster", "PoPS", "foreach", "iterators"), .inorder = TRUE) %do% {
     average_disagreements_odds_ratio <- foreach(p = 1:10, .combine = rbind, .packages = c("raster", "PoPS", "foreach"), .final = colMeans) %dopar% {
-      disagreements_odds_ratio <- data.frame(reproductive_rate = 0, natural_distance_scale = 0, total_disagreement = 0, quantity_disagreement = 0, allocation_disagreement = 0, odds_ratio = 0)
+      # disagreements_odds_ratio <- data.frame(reproductive_rate = 0, natural_distance_scale = 0, total_disagreement = 0, quantity_disagreement = 0, allocation_disagreement = 0, odds_ratio = 0)
       
-      data <- param_func(proposed_reproductive_rate, proposed_natural_distance_scale)
+      data <- param_func(proposed_reproductive_rate, proposed_natural_distance_scale, proposed_anthropogenic_distance_scale, proposed_percent_natural_dispersal)
       
       # set up comparison
       all_disagreement <- foreach(q = 1:length(data$infected), .combine = rbind, .packages =c("raster", "PoPS"), .final = colSums) %dopar% {
@@ -449,7 +571,7 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
       to.average_disagreements_odds_ratio <- all_disagreement
     }
     
-    proposed <- data.frame(t(average_disagreements_odds_ratio), reproductive_rate = proposed_reproductive_rate, natural_distance_scale = proposed_natural_distance_scale)
+    proposed <- data.frame(t(average_disagreements_odds_ratio), reproductive_rate = proposed_reproductive_rate, natural_distance_scale = proposed_natural_distance_scale, anthropogenic_distance_scale = proposed_anthropogenic_distance_scale, percent_natural_dispersal = proposed_percent_natural_dispersal)
     
     # make sure no proposed statistics are 0 or the calculation fails instead set them all to the lowest possible non-zero value
     if (proposed$allocation_disagreement == 0) {proposed$allocation_disagreement <- 1}
@@ -483,10 +605,22 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
         while (proposed_reproductive_rate <= 0) {
           proposed_reproductive_rate <- round(rnorm(1, mean = best$reproductive_rate, sd = sd_reproductive_rate), digits = 1)
         }
+        
         proposed_natural_distance_scale <- 0
         while (proposed_natural_distance_scale <= 0) {
           proposed_natural_distance_scale <- round(rnorm(1, mean = best$natural_distance_scale, sd = sd_natural_distance_scale), digits = 0)
         }
+        
+        proposed_anthropogenic_distance_scale <- 0
+        while (proposed_anthropogenic_distance_scale <= 0) {
+          proposed_anthropogenic_distance_scale <- round(rnorm(1, mean = best$anthropogenic_distance_scale, sd = sd_anthropogenic_distance_scale)/10, digits = 0) * 10
+        }
+        
+        proposed_percent_natural_dispersal <- 0
+        while (proposed_percent_natural_dispersal <= 0 || proposed_percent_natural_dispersal > 1.000) {
+          proposed_percent_natural_dispersal <- round(rnorm(1, mean = best$percent_natural_dispersal, sd = sd_percent_natural_dispersal), digits = 3)
+        }
+        
         to.params <- param
       }
     } else if (success_metric == "quantity and configuration") {
@@ -500,10 +634,21 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
         while (proposed_reproductive_rate <= 0) {
           proposed_reproductive_rate <- round(rnorm(1, mean = best$reproductive_rate, sd_reproductive_rate), digits = 1)
         }
+        
         proposed_natural_distance_scale <- 0.0
         while (proposed_natural_distance_scale <= 0) {
           proposed_natural_distance_scale <- round(rnorm(1, mean = best$natural_distance_scale, sd_natural_distance_scale), digits = 0)
         }
+        proposed_anthropogenic_distance_scale <- 0
+        while (proposed_anthropogenic_distance_scale <= 0) {
+          proposed_anthropogenic_distance_scale <- round(rnorm(1, mean = best$anthropogenic_distance_scale, sd = sd_anthropogenic_distance_scale)/10, digits = 0) * 10
+        }
+        
+        proposed_percent_natural_dispersal <- 0
+        while (proposed_percent_natural_dispersal <= 0 || proposed_percent_natural_dispersal > 1.000) {
+          proposed_percent_natural_dispersal <- round(rnorm(1, mean = best$percent_natural_dispersal, sd = sd_percent_natural_dispersal), digits = 3)
+        }
+        
         to.params <- param
       } else if (configuration_pass && quantity_pass == FALSE) {
         current <- proposed
@@ -515,6 +660,16 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
         while (proposed_natural_distance_scale <= 0) {
           proposed_natural_distance_scale <- round(rnorm(1, mean = best$natural_distance_scale, sd_natural_distance_scale), digits = 0)
         }
+        
+        proposed_anthropogenic_distance_scale <- 0
+        while (proposed_anthropogenic_distance_scale <= 0) {
+          proposed_anthropogenic_distance_scale <- round(rnorm(1, mean = best$anthropogenic_distance_scale, sd = sd_anthropogenic_distance_scale)/10, digits = 0) * 10
+        }
+        
+        proposed_percent_natural_dispersal <- 0
+        while (proposed_percent_natural_dispersal <= 0 || proposed_percent_natural_dispersal > 1.000) {
+          proposed_percent_natural_dispersal <- round(rnorm(1, mean = best$percent_natural_dispersal, sd = sd_percent_natural_dispersal), digits = 3)
+        }
         to.params <- param
       } else if (quantity_pass && configuration_pass == FALSE) {
         current <- proposed
@@ -544,6 +699,16 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
         while (proposed_natural_distance_scale <= 0) {
           proposed_natural_distance_scale <- round(rnorm(1, mean = best$natural_distance_scale, sd = sd_natural_distance_scale), digits = 0)
         }
+        proposed_anthropogenic_distance_scale <- 0
+        while (proposed_anthropogenic_distance_scale <= 0) {
+          proposed_anthropogenic_distance_scale <- round(rnorm(1, mean = best$anthropogenic_distance_scale, sd = sd_anthropogenic_distance_scale), digits = 3)
+        }
+        
+        proposed_percent_natural_dispersal <- 0
+        while (proposed_percent_natural_dispersal <= 0 || proposed_percent_natural_dispersal > 1.000) {
+          proposed_percent_natural_dispersal <- round(rnorm(1, mean = best$percent_natural_dispersal, sd = sd_percent_natural_dispersal), digits = 3)
+        }
+        
         to.params <- param
       } 
     } else if (success_metric == "residual error") {
@@ -562,6 +727,16 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
         while (proposed_natural_distance_scale <= 0) {
           proposed_natural_distance_scale <- round(rnorm(1, mean = best$natural_distance_scale, sd = sd_natural_distance_scale), digits = 0)
         }
+        proposed_anthropogenic_distance_scale <- 0
+        while (proposed_anthropogenic_distance_scale <= 0) {
+          proposed_anthropogenic_distance_scale <- round(rnorm(1, mean = best$anthropogenic_distance_scale, sd = sd_anthropogenic_distance_scale)/10, digits = 0) * 10
+        }
+        
+        proposed_percent_natural_dispersal <- 0
+        while (proposed_percent_natural_dispersal <= 0 || proposed_percent_natural_dispersal > 1.000) {
+          proposed_percent_natural_dispersal <- round(rnorm(1, mean = best$percent_natural_dispersal, sd = sd_percent_natural_dispersal), digits = 3)
+        }
+        
         to.params <- param
       } 
     } else {
@@ -571,5 +746,171 @@ calibrate <- function(infected_years_file, num_iterations, start_reproductive_ra
   }
   stopCluster(cl)
 
-  return(params)
-}
\ No newline at end of file
+  if (prior_number_of_observations < 1) {
+    prior_weight <- prior_number_of_observations
+    total_number_of_observations <- number_of_observations + round(number_of_observations * prior_number_of_observations)
+    weight <- 1 - prior_weight
+  } else if (prior_number_of_observations >= 1) {
+    total_number_of_observations <- prior_number_of_observations + number_of_observations
+    prior_weight <- prior_number_of_observations/total_number_of_observations
+    weight <- 1 - prior_weight
+  }
+  
+  ## Use prior and calibrated parameters to update to posteriors
+  count <- 10000000
+  # Reproductive Rate
+  if (class(prior_reproductive_rate) == "matrix" && nrow(prior_reproductive_rate) == 1) {
+    prior_reproductive_rates <- round(rnorm(count, start_reproductive_rate, sd_reproductive_rate), digits = 1)
+    prior_reproductive_rates <- as.data.frame(table(prior_reproductive_rates))
+    prior_reproductive_rates$prior_reproductive_rates <- as.numeric(as.character(prior_reproductive_rates$prior_reproductive_rates))
+    prior_reproductive_rates$prob <- round(prior_reproductive_rates$Freq/count, digits = 3)
+    prior_reproductive_rates <- prior_reproductive_rates[prior_reproductive_rates$prob > 0.00,]
+  } else if (class(prior_reproductive_rate) %in% c("matrix", "data.frame") && nrow(prior_reproductive_rate) > 1) {
+    prior_reproductive_rates <- prior_reproductive_rate[ , 1:2]
+    names(prior_reproductive_rates) <- c('prior_reproductive_rates', 'prob')
+    
+  }
+
+  calibration_count <- nrow(params)
+  calibrated_reproductive_rates <- as.data.frame(table(params$reproductive_rate))
+  calibrated_reproductive_rates$Var1 <- as.numeric(as.character(calibrated_reproductive_rates$Var1))
+  calibrated_reproductive_rates$prob <- round(calibrated_reproductive_rates$Freq/calibration_count, digits = 3)
+  
+  min_reproductive_rate <- min(min(prior_reproductive_rates$prior_reproductive_rates), min(calibrated_reproductive_rates$Var1))
+  max_reproductive_rate <- max(max(prior_reproductive_rates$prior_reproductive_rates), max(calibrated_reproductive_rates$Var1))
+  
+  reproductive_rates <- data.frame(reproductive_rate = round(seq(min_reproductive_rate, max_reproductive_rate, 0.1), digits = 1), 
+                                   prior_probability = rep(0, length(seq(min_reproductive_rate, max_reproductive_rate, 0.1))),
+                                   calibrated_probability = rep(0, length(seq(min_reproductive_rate, max_reproductive_rate, 0.1))),
+                                   posterior_probability = rep(0, length(seq(min_reproductive_rate, max_reproductive_rate, 0.1))))
+  
+  for (i in 1:nrow(reproductive_rates)) {
+    if (length(prior_reproductive_rates$prob[prior_reproductive_rates$prior_reproductive_rates == reproductive_rates$reproductive_rate[i]]) > 0) {
+      reproductive_rates$prior_probability[i] <- prior_reproductive_rates$prob[prior_reproductive_rates$prior_reproductive_rates == reproductive_rates$reproductive_rate[i]]
+    }
+    if (length(calibrated_reproductive_rates$prob[calibrated_reproductive_rates$Var1 == reproductive_rates$reproductive_rate[i]]) > 0) {
+      reproductive_rates$calibrated_probability[i] <- calibrated_reproductive_rates$prob[calibrated_reproductive_rates$Var1 == reproductive_rates$reproductive_rate[i]]
+    }
+  }
+  reproductive_rates$posterior_probability <- round(reproductive_rates$prior_probability*prior_weight + reproductive_rates$calibrated_probability * weight, digits = 3)
+  colSums(reproductive_rates)
+  posterior_reproductive_rates <- reproductive_rates[,c(1,4)]
+  
+  # Natural Distance Scale
+  if (class(prior_natural_distance_scale) == "matrix" && nrow(prior_natural_distance_scale) == 1) {
+    prior_natural_distance_scales <- round(rnorm(count, start_natural_distance_scale, sd_natural_distance_scale), digits = 0)
+    prior_natural_distance_scales <- as.data.frame(table(prior_natural_distance_scales))
+    prior_natural_distance_scales$prior_natural_distance_scales <- as.numeric(as.character(prior_natural_distance_scales$prior_natural_distance_scales))
+    prior_natural_distance_scales$prob <- round(prior_natural_distance_scales$Freq/count, digits = 3)
+    prior_natural_distance_scales <- prior_natural_distance_scales[prior_natural_distance_scales$prob > 0.000,]
+  } else if (class(prior_natural_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_natural_distance_scale) > 1) {
+    prior_natural_distance_scales <- prior_natural_distance_scale[ , 1:2]
+    names(prior_natural_distance_scales) <- c('prior_natural_distance_scales', 'prob')
+  }
+  
+  calibration_count <- nrow(params)
+  calibrated_natural_distance_scales <- as.data.frame(table(params$natural_distance_scale))
+  calibrated_natural_distance_scales$Var1 <- as.numeric(as.character(calibrated_natural_distance_scales$Var1))
+  calibrated_natural_distance_scales$prob <- round(calibrated_natural_distance_scales$Freq/calibration_count, digits = 3)
+  
+  min_natural_distance_scale <- min(min(prior_natural_distance_scales$prior_natural_distance_scales), min(calibrated_natural_distance_scales$Var1))
+  max_natural_distance_scale <- max(max(prior_natural_distance_scales$prior_natural_distance_scales), max(calibrated_natural_distance_scales$Var1))
+  
+  natural_distance_scales <- data.frame(natural_distance_scale = round(seq(min_natural_distance_scale, max_natural_distance_scale, 1), digits = 1), 
+                                   prior_probability = rep(0, length(seq(min_natural_distance_scale, max_natural_distance_scale, 1))),
+                                   calibrated_probability = rep(0, length(seq(min_natural_distance_scale, max_natural_distance_scale, 1))),
+                                   posterior_probability = rep(0, length(seq(min_natural_distance_scale, max_natural_distance_scale, 1))))
+  
+  for (i in 1:nrow(natural_distance_scales)) {
+    if (length(prior_natural_distance_scales$prob[prior_natural_distance_scales$prior_natural_distance_scales == natural_distance_scales$natural_distance_scale[i]]) > 0) {
+      natural_distance_scales$prior_probability[i] <- prior_natural_distance_scales$prob[prior_natural_distance_scales$prior_natural_distance_scales == natural_distance_scales$natural_distance_scale[i]]
+    }
+    if (length(calibrated_natural_distance_scales$prob[calibrated_natural_distance_scales$Var1 == natural_distance_scales$natural_distance_scale[i]]) > 0) {
+      natural_distance_scales$calibrated_probability[i] <- calibrated_natural_distance_scales$prob[calibrated_natural_distance_scales$Var1 == natural_distance_scales$natural_distance_scale[i]]
+    }
+  }
+  natural_distance_scales$posterior_probability <- round(natural_distance_scales$prior_probability*prior_weight + natural_distance_scales$calibrated_probability * weight, digits = 3)
+  colSums(natural_distance_scales)
+  posterior_natural_distance_scales <- natural_distance_scales[,c(1,4)]
+  
+  # Anthropogenic Distance Scale
+  if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_anthropogenic_distance_scale) == 1) {
+    prior_anthropogenic_distance_scales <- round(rnorm(count, start_anthropogenic_distance_scale, sd_anthropogenic_distance_scale)/10, digits = 0)*10
+    prior_anthropogenic_distance_scales <- as.data.frame(table(prior_anthropogenic_distance_scales))
+    prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales <- as.numeric(as.character(prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales))
+    prior_anthropogenic_distance_scales$prob <- round(prior_anthropogenic_distance_scales$Freq/count, digits = 3)
+    prior_anthropogenic_distance_scales <- prior_anthropogenic_distance_scales[prior_anthropogenic_distance_scales$prob > 0.000,]
+  } else if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_anthropogenic_distance_scale) > 1) {
+    prior_anthropogenic_distance_scales <- prior_anthropogenic_distance_scale[ , 1:2]
+    names(prior_anthropogenic_distance_scales) <- c('prior_anthropogenic_distance_scales', 'prob')
+  }
+  
+  calibration_count <- nrow(params)
+  calibrated_anthropogenic_distance_scales <- as.data.frame(table(params$anthropogenic_distance_scale))
+  calibrated_anthropogenic_distance_scales$Var1 <- as.numeric(as.character(calibrated_anthropogenic_distance_scales$Var1))
+  calibrated_anthropogenic_distance_scales$prob <- round(calibrated_anthropogenic_distance_scales$Freq/calibration_count, digits = 3)
+  
+  min_anthropogenic_distance_scale <- min(min(prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales), min(calibrated_anthropogenic_distance_scales$Var1))
+  max_anthropogenic_distance_scale <- max(max(prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales), max(calibrated_anthropogenic_distance_scales$Var1))
+  
+  anthropogenic_distance_scales <- data.frame(anthropogenic_distance_scale = round(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10), digits = 1), 
+                                        prior_probability = rep(0, length(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10))),
+                                        calibrated_probability = rep(0, length(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10))),
+                                        posterior_probability = rep(0, length(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10))))
+  
+  for (i in 1:nrow(anthropogenic_distance_scales)) {
+    if (length(prior_anthropogenic_distance_scales$prob[prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]) > 0) {
+      anthropogenic_distance_scales$prior_probability[i] <- prior_anthropogenic_distance_scales$prob[prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]
+    }
+    if (length(calibrated_anthropogenic_distance_scales$prob[calibrated_anthropogenic_distance_scales$Var1 == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]) > 0) {
+      anthropogenic_distance_scales$calibrated_probability[i] <- calibrated_anthropogenic_distance_scales$prob[calibrated_anthropogenic_distance_scales$Var1 == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]
+    }
+  }
+  anthropogenic_distance_scales$posterior_probability <- round(anthropogenic_distance_scales$prior_probability*prior_weight + anthropogenic_distance_scales$calibrated_probability * weight, digits = 3)
+  colSums(anthropogenic_distance_scales)
+  posterior_anthropogenic_distance_scales <- anthropogenic_distance_scales[,c(1,4)]
+  
+  # Percent Natural Dispersal
+  if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && nrow(prior_percent_natural_dispersal) == 1) {
+    prior_percent_natural_dispersals <- round(rnorm(count, start_percent_natural_dispersal, sd_percent_natural_dispersal), digits = 3)
+    prior_percent_natural_dispersals <- as.data.frame(table(prior_percent_natural_dispersals))
+    prior_percent_natural_dispersals$prior_percent_natural_dispersals <- as.numeric(as.character(prior_percent_natural_dispersals$prior_percent_natural_dispersals))
+    prior_percent_natural_dispersals$prob <- round(prior_percent_natural_dispersals$Freq/count, digits = 3)
+    prior_percent_natural_dispersals <- prior_percent_natural_dispersals[prior_percent_natural_dispersals$prob > 0.000,]
+    prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals == 1.000] <- sum(prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals >= 1.0])
+    prior_percent_natural_dispersals <- prior_percent_natural_dispersals[prior_percent_natural_dispersals$prior_percent_natural_dispersals <= 1.0,]
+  } else if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && nrow(prior_percent_natural_dispersal) > 1) {
+    prior_percent_natural_dispersals <- prior_percent_natural_dispersal[ , 1:2]
+    names(prior_percent_natural_dispersals) <- c('prior_percent_natural_dispersals', 'prob')
+  }
+  
+  calibration_count <- nrow(params)
+  calibrated_percent_natural_dispersals <- as.data.frame(table(params$percent_natural_dispersal))
+  calibrated_percent_natural_dispersals$Var1 <- as.numeric(as.character(calibrated_percent_natural_dispersals$Var1))
+  calibrated_percent_natural_dispersals$prob <- round(calibrated_percent_natural_dispersals$Freq/calibration_count, digits = 3)
+  
+  min_percent_natural_dispersal <- min(min(prior_percent_natural_dispersals$prior_percent_natural_dispersals), min(calibrated_percent_natural_dispersals$Var1))
+  max_percent_natural_dispersal <- max(max(prior_percent_natural_dispersals$prior_percent_natural_dispersals), max(calibrated_percent_natural_dispersals$Var1))
+  
+  percent_natural_dispersals <- data.frame(percent_natural_dispersal = round(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001), digits = 3), 
+                                        prior_probability = rep(0, length(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001))),
+                                        calibrated_probability = rep(0, length(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001))),
+                                        posterior_probability = rep(0, length(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001))))
+  
+  for (i in 1:nrow(percent_natural_dispersals)) {
+    if (length(prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals == percent_natural_dispersals$percent_natural_dispersal[i]]) > 0) {
+      percent_natural_dispersals$prior_probability[i] <- prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals == percent_natural_dispersals$percent_natural_dispersal[i]]
+    }
+    if (length(calibrated_percent_natural_dispersals$prob[calibrated_percent_natural_dispersals$Var1 == percent_natural_dispersals$percent_natural_dispersal[i]]) > 0) {
+      percent_natural_dispersals$calibrated_probability[i] <- calibrated_percent_natural_dispersals$prob[calibrated_percent_natural_dispersals$Var1 == percent_natural_dispersals$percent_natural_dispersal[i]]
+    }
+  }
+  percent_natural_dispersals$posterior_probability <- round(percent_natural_dispersals$prior_probability*prior_weight + percent_natural_dispersals$calibrated_probability * weight, digits = 3)
+  colSums(percent_natural_dispersals)
+  posterior_percent_natural_dispersals <- percent_natural_dispersals[,c(1,4)]
+  
+  outputs <- list(posterior_reproductive_rates, posterior_natural_distance_scales, posterior_anthropogenic_distance_scales, posterior_percent_natural_dispersals, total_number_of_observations, reproductive_rates, natural_distance_scales, anthropogenic_distance_scales, percent_natural_dispersals, params)
+  names(outputs) <- c('posterior_reproductive_rates', 'posterior_natural_distance_scales', 'posterior_anthropogenic_distance_scales', 'posterior_percent_natural_dispersals', 'total_number_of_observations', 'reproductive_rates', 'natural_distance_scales', 'anthropogenic_distance_scales', 'percent_natural_dispersals', 'raw_calibration_data')
+  
+  return(outputs)
+}
diff --git a/man/calibrate.Rd b/man/calibrate.Rd
index 5cc1b6b0..69a753cc 100644
--- a/man/calibrate.Rd
+++ b/man/calibrate.Rd
@@ -4,41 +4,45 @@
 \alias{calibrate}
 \title{Calibrates the reproductive rate and dispersal scales of the pops model.}
 \usage{
-calibrate(infected_years_file, num_iterations, start_reproductive_rate,
-  number_of_cores = NA, start_natural_distance_scale,
-  sd_reproductive_rate, sd_natural_distance_scale, infected_file,
+calibrate(infected_years_file, num_iterations, number_of_cores = NA,
+  number_of_observations, prior_number_of_observations,
+  prior_reproductive_rate, prior_natural_distance_scale,
+  prior_percent_natural_dispersal = c(1, 0),
+  prior_anthropogenic_distance_scale = c(1000, 0), infected_file,
   host_file, total_plants_file, temp = FALSE,
   temperature_coefficient_file = "", precip = FALSE,
   precipitation_coefficient_file = "", time_step = "month",
-  reproductive_rate = 3, season_month_start = 1,
-  season_month_end = 12, start_date = "2008-01-01",
-  end_date = "2008-12-31", use_lethal_temperature = FALSE,
-  temperature_file = "", lethal_temperature = -12.87,
-  lethal_temperature_month = 1, mortality_on = FALSE,
-  mortality_rate = 0, mortality_time_lag = 0, management = FALSE,
-  treatment_dates = c(0), treatments_file = "",
-  treatment_method = "ratio", percent_natural_dispersal = 1,
-  natural_kernel_type = "cauchy", anthropogenic_kernel_type = "cauchy",
-  natural_distance_scale = 21, anthropogenic_distance_scale = 0,
-  natural_dir = "NONE", natural_kappa = 0,
-  anthropogenic_dir = "NONE", anthropogenic_kappa = 0,
-  pesticide_duration = c(0), pesticide_efficacy = 1, mask = NULL,
-  success_metric = "quantity", output_frequency = "year")
+  season_month_start = 1, season_month_end = 12,
+  start_date = "2008-01-01", end_date = "2008-12-31",
+  use_lethal_temperature = FALSE, temperature_file = "",
+  lethal_temperature = -12.87, lethal_temperature_month = 1,
+  mortality_on = FALSE, mortality_rate = 0, mortality_time_lag = 0,
+  management = FALSE, treatment_dates = c(0), treatments_file = "",
+  treatment_method = "ratio", natural_kernel_type = "cauchy",
+  anthropogenic_kernel_type = "cauchy", natural_dir = "NONE",
+  natural_kappa = 0, anthropogenic_dir = "NONE",
+  anthropogenic_kappa = 0, pesticide_duration = c(0),
+  pesticide_efficacy = 1, mask = NULL, success_metric = "quantity",
+  output_frequency = "year")
 }
 \arguments{
 \item{infected_years_file}{Raster file with years of initial infection/infestation as individual locations of a pest or pathogen}
 
 \item{num_iterations}{how many iterations do you want to run to allow the calibration to converge (recommend a minimum of at least 100,000 but preferably 1 million).}
 
-\item{start_reproductive_rate}{starting reproductive rate for MCMC calibration (affects how quickly a series converges) numeric value > 0}
-
 \item{number_of_cores}{number of cores to use for calibration (defaults to the number of cores available on the machine) integer value >= 1}
 
-\item{start_natural_distance_scale}{starting short distance scale parameter for MCMC calibration (affects how quickly a series converges) numeric value > 0}
+\item{number_of_observations}{the number of observations used for this calibartion}
+
+\item{prior_number_of_observations}{the number of total observations from previous calibrations used to weight the posterior distributions (if this is a new calibration this value takes the form of a prior weight (0 - 1))}
+
+\item{prior_reproductive_rate}{the prior reproductive rate for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1}
+
+\item{prior_natural_distance_scale}{the prior natural distance scale for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1}
 
-\item{sd_reproductive_rate}{starting standard deviation for reproductive rate for MCMC calibration (can affect how quickly and if a series converges) numeric value > 0}
+\item{prior_percent_natural_dispersal}{the prior percent natural distance for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1}
 
-\item{sd_natural_distance_scale}{starting standard deviation for short distance scale for MCMC calibration (can affect how quickly and if a series converges) numeric value > 0}
+\item{prior_anthropogenic_distance_scale}{the prior anthropogenic distance scale for MCMC calibration as a list with mean and standard deviation ( e.g. c(mean, sd)) with mean and sd being numeric values or as a 2-column matrix with value and probability as columns probabilites must sum to 1}
 
 \item{infected_file}{path to raster file with initial infections}
 
@@ -56,8 +60,6 @@ calibrate(infected_years_file, num_iterations, start_reproductive_rate,
 
 \item{time_step}{how often should spread occur options: ('day', 'week', 'month')}
 
-\item{reproductive_rate}{number of spores or pest units produced by a single host under optimal weather conditions}
-
 \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)}
@@ -88,16 +90,10 @@ calibrate(infected_years_file, num_iterations, start_reproductive_rate,
 
 \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{percent_natural_dispersal}{what percentage of dispersal is natural range versus anthropogenic range value between 0 and 1}
-
 \item{natural_kernel_type}{what type of dispersal kernel should be used for natural dispersal can be ('cauchy', 'exponential')}
 
 \item{anthropogenic_kernel_type}{what type of dispersal kernel should be used for anthropogenic dispersal can be ('cauchy', 'exponential')}
 
-\item{natural_distance_scale}{distance scale parameter for natural range dispersal kernel numeric value > 0}
-
-\item{anthropogenic_distance_scale}{distance scale parameter for anthropogenic range dispersal kernel numeric value > 0}
-
 \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{natural_kappa}{sets the strength of the natural direction in the von-mises distribution numeric value between 0.01 and 12}

From 0cdf02df682606377e16245f4a95a9ec5d2f7c71 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 10:22:18 -0500
Subject: [PATCH 04/25] improved checks

---
 R/pops.r | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/R/pops.r b/R/pops.r
index d366a67f..7efc8252 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -1,11 +1,11 @@
 #' PoPS (Pest or Pathogen Spread) model
 #' 
 #' A dynamic species distribution model for pest or pathogen spread in forest or agricultural ecosystems. The model is process based
-#' meaning that it uses understanding of the effect of weather on reproduction and survival of the pest/pathogen in order to forecast
+#' meaning that it uses understanding of the effect of weather and other environmental factors on reproduction and survival of the pest/pathogen in order to forecast
 #' spread of the pest/pathogen into the future. 
 #'
 #' @param infected_file path to raster file with initial infections 
-#' @param host_file path to raster file with number of hosts
+#' @param 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)
 #' @param total_plants_file path to raster file with number of total plants
 #' @param temp boolean that allows the use of temperature coefficients to modify spread (TRUE or FALSE)
 #' @param temperature_coefficient_file path to raster file with temperature coefficient data for the timestep and number of years specified
@@ -125,7 +125,7 @@ pops <- function(infected_file, host_file, total_plants_file,
   }
   
   if (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d"))){
-    return("End time and/or start time not of type numeric and/or in format YYYY")
+    return("End time and/or start time not of type numeric and/or in format YYYY-MM-DD")
   }
   
   if (!(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
@@ -162,14 +162,11 @@ pops <- function(infected_file, host_file, total_plants_file,
   
   infected <- raster::raster(infected_file)
   infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  # infected[is.na(infected)] <- 0
   host <- raster::raster(host_file)
   host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  # host[is.na(host)] <- 0
   total_plants <- raster::raster(total_plants_file)
   total_plants <- raster::reclassify(total_plants, matrix(c(NA, 0), ncol = 2, byrow = TRUE), right = NA)
-  # total_plants[is.na(total_plants)] <- 0
-  
+
   if (!(raster::extent(infected) == raster::extent(host) && raster::extent(infected) == raster::extent(total_plants))) {
     return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
   }
@@ -184,7 +181,6 @@ pops <- function(infected_file, host_file, total_plants_file,
   
   susceptible <- host - infected
   susceptible <- raster::reclassify(susceptible, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  # susceptible[is.na(susceptible)] <- 0
   susceptible[susceptible < 0] <- 0
   
   if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
@@ -198,8 +194,7 @@ pops <- function(infected_file, host_file, total_plants_file,
   if (use_lethal_temperature == TRUE) {
     temperature_stack <- raster::stack(temperature_file)
     temperature_stack <- raster::reclassify(temperature_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-    # temperature_stack[is.na(temperature_stack)] <- 0
-    
+
     if (!(raster::extent(infected) == raster::extent(temperature_stack))) {
       return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
     }

From 2a9052eb6866ef112e1e1de09125c9550178d644 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 11:49:30 -0500
Subject: [PATCH 05/25] updated documentation

---
 man/calibrate.Rd     | 2 +-
 man/pops.Rd          | 4 ++--
 man/pops_multirun.Rd | 2 +-
 man/validate.Rd      | 2 +-
 4 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/man/calibrate.Rd b/man/calibrate.Rd
index 69a753cc..8b4be88c 100644
--- a/man/calibrate.Rd
+++ b/man/calibrate.Rd
@@ -46,7 +46,7 @@ calibrate(infected_years_file, num_iterations, number_of_cores = NA,
 
 \item{infected_file}{path to raster file with initial infections}
 
-\item{host_file}{path to raster file with number of hosts}
+\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)}
 
 \item{total_plants_file}{path to raster file with number of total plants}
 
diff --git a/man/pops.Rd b/man/pops.Rd
index cab99b48..7a1ea377 100644
--- a/man/pops.Rd
+++ b/man/pops.Rd
@@ -25,7 +25,7 @@ pops(infected_file, host_file, total_plants_file, temp = FALSE,
 \arguments{
 \item{infected_file}{path to raster file with initial infections}
 
-\item{host_file}{path to raster file with number of hosts}
+\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)}
 
 \item{total_plants_file}{path to raster file with number of total plants}
 
@@ -102,6 +102,6 @@ list of infected and susceptible per year
 }
 \description{
 A dynamic species distribution model for pest or pathogen spread in forest or agricultural ecosystems. The model is process based
-meaning that it uses understanding of the effect of weather on reproduction and survival of the pest/pathogen in order to forecast
+meaning that it uses understanding of the effect of weather and other environmental factors on reproduction and survival of the pest/pathogen in order to forecast
 spread of the pest/pathogen into the future.
 }
diff --git a/man/pops_multirun.Rd b/man/pops_multirun.Rd
index 7fe6c1c5..8664a5c7 100644
--- a/man/pops_multirun.Rd
+++ b/man/pops_multirun.Rd
@@ -26,7 +26,7 @@ pops_multirun(infected_file, host_file, total_plants_file, temp = FALSE,
 \arguments{
 \item{infected_file}{path to raster file with initial infections}
 
-\item{host_file}{path to raster file with number of hosts}
+\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)}
 
 \item{total_plants_file}{path to raster file with number of total plants}
 
diff --git a/man/validate.Rd b/man/validate.Rd
index d4b32a7b..301d2250 100644
--- a/man/validate.Rd
+++ b/man/validate.Rd
@@ -32,7 +32,7 @@ validate(infected_years_file, num_iterations, number_of_cores = NA,
 
 \item{infected_file}{path to raster file with initial infections}
 
-\item{host_file}{path to raster file with number of hosts}
+\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)}
 
 \item{total_plants_file}{path to raster file with number of total plants}
 

From fdf6f6bac7bed938165df144099ca8ef34623637 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 11:52:31 -0500
Subject: [PATCH 06/25] add checks to keep main functions more compact and have
 less copy and paste code checks

---
 R/checks.R | 102 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 102 insertions(+)
 create mode 100644 R/checks.R

diff --git a/R/checks.R b/R/checks.R
new file mode 100644
index 00000000..880eb2df
--- /dev/null
+++ b/R/checks.R
@@ -0,0 +1,102 @@
+## These functions are designed to improve data format checks and reduce copy and pasting of code across functions
+
+initial_raster_check <- function(x) {
+  checks_passed <- TRUE
+  
+  if (!file.exists(x)) {
+    checks_passed <- FALSE
+    failed_check <- "file does not exist" 
+  }
+  
+  if (!(raster::extension(x) %in% c(".grd", ".tif", ".img"))) {
+    checks_passed <- FALSE
+    failed_check <- "file is not one of '.grd', '.tif', '.img'"
+  }
+  
+  if (checks_passed) {
+    r<- raster::raster(x)
+    r <- raster::reclassify(r, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, r)
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    return(outs)
+  }
+}
+
+
+
+
+
+time_checks <- function(end_date, start_date, time_step, output_frequency) {
+  checks_passed <- TRUE
+  
+  if (!(time_step %in% list("week", "month", "day"))) {
+    checks_passed <- FALSE
+    failed_check <- "Time step must be one of 'week', 'month' or 'day'"
+  }
+  
+  if (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d"))){
+    checks_passed <- FALSE
+    failed_check <- "End time and/or start time not of type numeric and/or in format YYYY-MM-DD"
+  }
+  
+  if (!(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
+    checks_passed <- FALSE
+    failed_check <- "Time step must be one of 'week', 'month' or 'day'"
+  }
+  
+  if (output_frequency == "day") {
+    if (time_step == "week" || time_step == "month") {
+      checks_passed <- FALSE
+      failed_check <- "Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible"
+    }
+  }
+  
+  if (output_frequency == "week") {
+    if (time_step == "month") {
+      checks_passed <- FALSE
+      failed_check <- "Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible"
+    }
+  }
+  
+  
+  duration <- lubridate::interval(start_date, end_date)
+  
+  if (time_step == "week") {
+    number_of_time_steps <- ceiling(lubridate::time_length(duration, "week"))
+  } else if (time_step == "month") {
+    number_of_time_steps <- ceiling(lubridate::time_length(duration, "month"))
+  } else if (time_step == "day") {
+    number_of_time_steps <- ceiling(lubridate::time_length(duration, "day"))
+  }
+  
+  number_of_years <- ceiling(lubridate::time_length(duration, "year"))
+  
+  if (output_frequency == "week") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "week"))
+  } else if (output_frequency == "month") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "month"))
+  } else if (output_frequency == "day") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "day"))
+  } else if (output_frequency == "year") {
+    number_of_outputs <- ceiling(lubridate::time_length(duration, "year"))
+  } else if (output_frequency == "time_step") {
+    number_of_outputs <- number_of_time_steps
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, number_of_time_steps, number_of_years, number_of_outputs)
+    names(outs) <- c('checks_passed', 'number_of_time_steps', 'number_of_years', 'number_of_outputs')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+  
+}
+

From 1e5026f520f52de3005243da1e9ea0d68316fc2e Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 11:53:09 -0500
Subject: [PATCH 07/25] added time checks to all r wrapper functions

---
 R/calibrate.R     | 53 +++++++----------------------------------------
 R/pops.r          | 41 +++++++-----------------------------
 R/pops_multirun.R | 41 +++++++-----------------------------
 R/validate.R      | 53 +++++++----------------------------------------
 4 files changed, 28 insertions(+), 160 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 4f03d1f6..0d2711c4 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -98,52 +98,13 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  if (!(time_step %in% list("week", "month", "day"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d"))){
-    return("End time and/or start time not of type numeric and/or in format YYYY")
-  }
-  
-  if (!(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (output_frequency == "day") {
-    if (time_step == "week" || time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  if (output_frequency == "week") {
-    if (time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  duration <- lubridate::interval(start_date, end_date)
-  
-  if (time_step == "week") {
-    number_of_time_steps <- ceiling(lubridate::time_length(duration, "week"))
-  } else if (time_step == "month") {
-    number_of_time_steps <- ceiling(lubridate::time_length(duration, "month"))
-  } else if (time_step == "day") {
-    number_of_time_steps <- ceiling(lubridate::time_length(duration, "day"))
-  }
-  
-  number_of_years <- ceiling(lubridate::time_length(duration, "year"))
-  
-  if (output_frequency == "week") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "week"))
-  } else if (output_frequency == "month") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "month"))
-  } else if (output_frequency == "day") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "day"))
-  } else if (output_frequency == "year") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "year"))
-  } else if (output_frequency == "time_step") {
-    number_of_outputs <- number_of_time_steps
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  if(time_check$checks_passed) {
+    number_of_time_steps <- time_check$number_of_time_steps
+    number_of_years <- time_check$number_of_years
+    number_of_outputs <- time_check$number_of_outputs
+  } else {
+    return(time_check$failed_check)
   }
   
   ## Setup for reproductive rate to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
diff --git a/R/pops.r b/R/pops.r
index 7efc8252..a561946a 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -120,42 +120,15 @@ pops <- function(infected_file, host_file, total_plants_file,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  if (!(time_step %in% list("week", "month", "day"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d"))){
-    return("End time and/or start time not of type numeric and/or in format YYYY-MM-DD")
-  }
-  
-  if (!(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (output_frequency == "day") {
-    if (time_step == "week" || time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  if (output_frequency == "week") {
-    if (time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  duration <- lubridate::interval(start_date, end_date)
-  
-  if (time_step == "week") {
-    number_of_time_steps <- ceiling(time_length(duration, "week"))
-  } else if (time_step == "month") {
-    number_of_time_steps <- ceiling(time_length(duration, "month"))
-  } else if (time_step == "day") {
-    number_of_time_steps <- ceiling(time_length(duration, "day"))
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  if(time_check$checks_passed) {
+    number_of_time_steps <- time_check$number_of_time_steps
+    number_of_years <- time_check$number_of_years
+    number_of_outputs <- time_check$number_of_outputs
+  } else {
+    return(time_check$failed_check)
   }
   
-  number_of_years <- ceiling(time_length(duration, "year"))
-  
   if (is.null(random_seed)) {
     random_seed = round(stats::runif(1, 1, 1000000))
   }
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index 03db037f..6a01dff5 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -89,42 +89,15 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  if (!(time_step %in% list("week", "month", "day"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d"))){
-    return("End time and/or start time not of type numeric and/or in format YYYY")
-  }
-  
-  if (!(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (output_frequency == "day") {
-    if (time_step == "week" || time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  if (output_frequency == "week") {
-    if (time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  duration <- lubridate::interval(start_date, end_date)
-  
-  if (time_step == "week") {
-    number_of_time_steps <- ceiling(time_length(duration, "week"))
-  } else if (time_step == "month") {
-    number_of_time_steps <- ceiling(time_length(duration, "month"))
-  } else if (time_step == "day") {
-    number_of_time_steps <- ceiling(time_length(duration, "day"))
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  if(time_check$checks_passed) {
+    number_of_time_steps <- time_check$number_of_time_steps
+    number_of_years <- time_check$number_of_years
+    number_of_outputs <- time_check$number_of_outputs
+  } else {
+    return(time_check$failed_check)
   }
   
-  number_of_years <- ceiling(time_length(duration, "year"))
-  
   infected <- raster::raster(infected_file)
   infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
diff --git a/R/validate.R b/R/validate.R
index cc319f76..ed375aad 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -111,52 +111,13 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  if (!(time_step %in% list("week", "month", "day"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d"))){
-    return("End time and/or start time not of type numeric and/or in format YYYY")
-  }
-  
-  if (!(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
-    return("Time step must be one of 'week', 'month' or 'day'")
-  }
-  
-  if (output_frequency == "day") {
-    if (time_step == "week" || time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  if (output_frequency == "week") {
-    if (time_step == "month") {
-      return("Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
-    }
-  }
-  
-  duration <- lubridate::interval(start_date, end_date)
-  
-  if (time_step == "week") {
-    number_of_time_steps <- ceiling(time_length(duration, "week"))
-  } else if (time_step == "month") {
-    number_of_time_steps <- ceiling(time_length(duration, "month"))
-  } else if (time_step == "day") {
-    number_of_time_steps <- ceiling(time_length(duration, "day"))
-  }
-  
-  number_of_years <- ceiling(time_length(duration, "year"))
-  
-  if (output_frequency == "week") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "week"))
-  } else if (output_frequency == "month") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "month"))
-  } else if (output_frequency == "day") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "day"))
-  } else if (output_frequency == "year") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "year"))
-  } else if (output_frequency == "time_step") {
-    number_of_outputs <- number_of_time_steps
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  if(time_check$checks_passed) {
+    number_of_time_steps <- time_check$number_of_time_steps
+    number_of_years <- time_check$number_of_years
+    number_of_outputs <- time_check$number_of_outputs
+  } else {
+    return(time_check$failed_check)
   }
   
   infected <- raster::raster(infected_file)

From 11b983ac9af26ae7322486faede728f9fa275c23 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 12:18:47 -0500
Subject: [PATCH 08/25] added percent checks and implemented them in r wrappers

---
 R/calibrate.R     | 57 +++++++++++++++++++++++------------------------
 R/checks.R        | 23 ++++++++++++++++++-
 R/pops.r          | 33 +++++++++++++--------------
 R/pops_multirun.R | 33 +++++++++++++--------------
 R/validate.R      | 33 +++++++++++++--------------
 5 files changed, 98 insertions(+), 81 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 0d2711c4..0c273ede 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -74,30 +74,6 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     return("treatment method is not one of the valid treatment options")
   }
   
-  if (!file.exists(infected_file)) {
-    return("Infected file does not exist") 
-  }
-  
-  if (!(raster::extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Infected file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (!file.exists(host_file)) {
-    return("Host file does not exist") 
-  }
-  
-  if (!(raster::extension(host_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Host file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (!file.exists(total_plants_file)) {
-    return("Total plants file does not exist") 
-  }
-  
-  if (!(raster::extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Total plants file is not one of '.grd', '.tif', '.img'")
-  }
-  
   time_check <- time_checks(end_date, start_date, time_step, output_frequency)
   if(time_check$checks_passed) {
     number_of_time_steps <- time_check$number_of_time_steps
@@ -203,6 +179,30 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     return("Incorrect format for prior percent natural distance")
   }
   
+  if (!file.exists(infected_file)) {
+    return("Infected file does not exist") 
+  }
+  
+  if (!(raster::extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
+    return("Infected file is not one of '.grd', '.tif', '.img'")
+  }
+  
+  if (!file.exists(host_file)) {
+    return("Host file does not exist") 
+  }
+  
+  if (!(raster::extension(host_file) %in% c(".grd", ".tif", ".img"))) {
+    return("Host file is not one of '.grd', '.tif', '.img'")
+  }
+  
+  if (!file.exists(total_plants_file)) {
+    return("Total plants file does not exist") 
+  }
+  
+  if (!(raster::extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
+    return("Total plants file is not one of '.grd', '.tif', '.img'")
+  }
+  
   infected <- raster::raster(infected_file)
   infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
@@ -395,12 +395,11 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     treatment_maps <- list(raster::as.matrix(treatment_map))
   }
   
-  if(start_percent_natural_dispersal == 1.0) {
-    use_anthropogenic_kernel = FALSE
-  } else if (start_percent_natural_dispersal < 1.0  && start_percent_natural_dispersal >= 0.0) {
-    use_anthropogenic_kernel = TRUE
+  percent_check <- percent_checks(start_percent_natural_dispersal)
+  if (percent_check$checks_passed){
+    use_anthropogenic_kernel <- percent_check$use_anthropogenic_kernel
   } else {
-    return("Percent natural dispersal must be between 0.0 and 1.0")
+    return(percent_check$failed_check)
   }
   
   ew_res <- xres(susceptible)
diff --git a/R/checks.R b/R/checks.R
index 880eb2df..b5092e1e 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -28,7 +28,28 @@ initial_raster_check <- function(x) {
 }
 
 
-
+percent_checks <- function(percent_natural_dispersal) {
+  checks_passed <- TRUE
+  
+  if(percent_natural_dispersal == 1.0) {
+    use_anthropogenic_kernel = FALSE
+  } else if (percent_natural_dispersal < 1.0  && percent_natural_dispersal >= 0.0) {
+    use_anthropogenic_kernel = TRUE
+  } else {
+    checks_passed <- FALSE
+    failed_check <- "Percent natural dispersal must be between 0.0 and 1.0"
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, use_anthropogenic_kernel)
+    names(outs) <- c('checks_passed', 'use_anthropogenic_kernel')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+}
 
 
 time_checks <- function(end_date, start_date, time_step, output_frequency) {
diff --git a/R/pops.r b/R/pops.r
index a561946a..269c8c7f 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -96,6 +96,22 @@ pops <- function(infected_file, host_file, total_plants_file,
     return("treatment method is not one of the valid treatment options")
   }
   
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  if(time_check$checks_passed) {
+    number_of_time_steps <- time_check$number_of_time_steps
+    number_of_years <- time_check$number_of_years
+    number_of_outputs <- time_check$number_of_outputs
+  } else {
+    return(time_check$failed_check)
+  }
+  
+  percent_check <- percent_checks(percent_natural_dispersal)
+  if (percent_check$checks_passed){
+    use_anthropogenic_kernel <- percent_check$use_anthropogenic_kernel
+  } else {
+    return(percent_check$failed_check)
+  }
+  
   if (!file.exists(infected_file)) {
     return("Infected file does not exist") 
   }
@@ -120,15 +136,6 @@ pops <- function(infected_file, host_file, total_plants_file,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
-  if(time_check$checks_passed) {
-    number_of_time_steps <- time_check$number_of_time_steps
-    number_of_years <- time_check$number_of_years
-    number_of_outputs <- time_check$number_of_outputs
-  } else {
-    return(time_check$failed_check)
-  }
-  
   if (is.null(random_seed)) {
     random_seed = round(stats::runif(1, 1, 1000000))
   }
@@ -330,14 +337,6 @@ pops <- function(infected_file, host_file, total_plants_file,
     treatment_maps <- list(raster::as.matrix(treatment_map))
   }
   
-  if(percent_natural_dispersal == 1.0) {
-    use_anthropogenic_kernel = FALSE
-  } else if (percent_natural_dispersal < 1.0  && percent_natural_dispersal >= 0.0) {
-    use_anthropogenic_kernel = TRUE
-  } else {
-    return("Percent natural dispersal must be between 0.0 and 1.0")
-  }
-  
   ew_res <- raster::xres(susceptible)
   ns_res <- raster::yres(susceptible)
   num_cols <- raster::ncol(susceptible)
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index 6a01dff5..e60db8fb 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -65,6 +65,22 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     return("treatment method is not one of the valid treatment options")
   }
   
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  if(time_check$checks_passed) {
+    number_of_time_steps <- time_check$number_of_time_steps
+    number_of_years <- time_check$number_of_years
+    number_of_outputs <- time_check$number_of_outputs
+  } else {
+    return(time_check$failed_check)
+  }
+  
+  percent_check <- percent_checks(percent_natural_dispersal)
+  if (percent_check$checks_passed){
+    use_anthropogenic_kernel <- percent_check$use_anthropogenic_kernel
+  } else {
+    return(percent_check$failed_check)
+  }
+  
   if (!file.exists(infected_file)) {
     return("Infected file does not exist") 
   }
@@ -89,15 +105,6 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
-  if(time_check$checks_passed) {
-    number_of_time_steps <- time_check$number_of_time_steps
-    number_of_years <- time_check$number_of_years
-    number_of_outputs <- time_check$number_of_outputs
-  } else {
-    return(time_check$failed_check)
-  }
-  
   infected <- raster::raster(infected_file)
   infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
@@ -291,14 +298,6 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     treatment_maps <- list(raster::as.matrix(treatment_map))
   }
   
-  if(percent_natural_dispersal == 1.0) {
-    use_anthropogenic_kernel = FALSE
-  } else if (percent_natural_dispersal < 1.0  && percent_natural_dispersal >= 0.0) {
-    use_anthropogenic_kernel = TRUE
-  } else {
-    return("Percent natural dispersal must be between 0.0 and 1.0")
-  }
-  
   ew_res <- raster::xres(susceptible)
   ns_res <- raster::yres(susceptible)
   num_cols <- raster::ncol(susceptible)
diff --git a/R/validate.R b/R/validate.R
index ed375aad..4065c8f4 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -87,6 +87,22 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     return("treatment method is not one of the valid treatment options")
   }
   
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  if(time_check$checks_passed) {
+    number_of_time_steps <- time_check$number_of_time_steps
+    number_of_years <- time_check$number_of_years
+    number_of_outputs <- time_check$number_of_outputs
+  } else {
+    return(time_check$failed_check)
+  }
+  
+  percent_check <- percent_checks(percent_natural_dispersal)
+  if (percent_check$checks_passed){
+    use_anthropogenic_kernel <- percent_check$use_anthropogenic_kernel
+  } else {
+    return(percent_check$failed_check)
+  }
+  
   if (!file.exists(infected_file)) {
     return("Infected file does not exist") 
   }
@@ -111,15 +127,6 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
-  if(time_check$checks_passed) {
-    number_of_time_steps <- time_check$number_of_time_steps
-    number_of_years <- time_check$number_of_years
-    number_of_outputs <- time_check$number_of_outputs
-  } else {
-    return(time_check$failed_check)
-  }
-  
   infected <- raster::raster(infected_file)
   infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
@@ -313,14 +320,6 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     treatment_maps <- list(raster::as.matrix(treatment_map))
   }
   
-  if(percent_natural_dispersal == 1.0) {
-    use_anthropogenic_kernel = FALSE
-  } else if (percent_natural_dispersal < 1.0  && percent_natural_dispersal >= 0.0) {
-    use_anthropogenic_kernel = TRUE
-  } else {
-    return("Percent natural dispersal must be between 0.0 and 1.0")
-  }
-  
   ew_res <- raster::xres(susceptible)
   ns_res <- raster::yres(susceptible)
   num_cols <- raster::ncol(susceptible)

From 3a66ee106503d80485052d2e16bc1ef8ac95e9cc Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 12:43:49 -0500
Subject: [PATCH 09/25] added treatment and metric checks and implemented in
 wrappers

---
 R/calibrate.R     | 18 +++++++-----------
 R/checks.R        | 47 ++++++++++++++++++++++++++++++++++++++++++++++-
 R/pops.r          |  5 +++--
 R/pops_multirun.R |  5 +++--
 R/validate.R      | 16 +++++++---------
 5 files changed, 66 insertions(+), 25 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 0c273ede..ed39860f 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -58,20 +58,16 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
                       pesticide_duration = c(0), pesticide_efficacy = 1.0,
                       mask = NULL, success_metric = "quantity", output_frequency = "year") { 
   
-  if (success_metric == "quantity") {
-    configuration = FALSE
-  } else if (success_metric == "quantity and configuration") {
-    configuration = TRUE
-  } else if (success_metric == "odds ratio") {
-    configuration = FALSE
-  } else if (success_metric == "residual error"){
-    configuration = FALSE
+  metric_check <- metric_checks(success_metric)
+  if (metric_check$checks_passed){
+    configuration <- metric_check$configuration
   } else {
-    return("Success metric must be one of 'quantity', 'quantity and configuration', 'residual error', or 'odds ratio'")
+    return(metric_check$failed_check)
   }
   
-  if (!treatment_method %in% c("ratio", "all infected")) {
-    return("treatment method is not one of the valid treatment options")
+  treatment_metric_check <- treatment_metric_checks(treatment_method)
+  if (!treatment_metric_check$checks_passed) {
+    return(treatment_metric_check$failed_check)
   }
   
   time_check <- time_checks(end_date, start_date, time_step, output_frequency)
diff --git a/R/checks.R b/R/checks.R
index b5092e1e..7dc72f67 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -28,6 +28,52 @@ initial_raster_check <- function(x) {
 }
 
 
+
+treatment_metric_checks <- function(treatment_method) {
+  checks_passed <- TRUE
+  
+  if (!treatment_method %in% c("ratio", "all infected")) {
+    checks_passed <- FALSE
+    failed_check <- "treatment method is not one of the valid treatment options"
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed)
+    names(outs) <- c('checks_passed')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+}
+
+metric_checks <- function(success_metric) {
+  checks_passed <- TRUE
+  
+  if (success_metric == "quantity") {
+    configuration <- FALSE
+  } else if (success_metric == "quantity and configuration") {
+    configuration <- TRUE
+  } else if (success_metric == "odds_ratio") {
+    configuration <- FALSE
+  } else {
+    checks_passed <- FALSE
+    failed_check <- "Success metric must be one of 'quantity', 'quantity and configuration', or 'odds_ratio'"
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, configuration)
+    names(outs) <- c('checks_passed', 'configuration')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+  
+}
+
 percent_checks <- function(percent_natural_dispersal) {
   checks_passed <- TRUE
   
@@ -51,7 +97,6 @@ percent_checks <- function(percent_natural_dispersal) {
   }
 }
 
-
 time_checks <- function(end_date, start_date, time_step, output_frequency) {
   checks_passed <- TRUE
   
diff --git a/R/pops.r b/R/pops.r
index 269c8c7f..1a4adbf8 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -92,8 +92,9 @@ pops <- function(infected_file, host_file, total_plants_file,
                  pesticide_duration = c(0), pesticide_efficacy = 1.0,
                  random_seed = NULL, output_frequency = "year"){ 
 
-  if (!treatment_method %in% c("ratio", "all infected")) {
-    return("treatment method is not one of the valid treatment options")
+  treatment_metric_check <- treatment_metric_checks(treatment_method)
+  if (!treatment_metric_check$checks_passed) {
+    return(treatment_metric_check$failed_check)
   }
   
   time_check <- time_checks(end_date, start_date, time_step, output_frequency)
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index e60db8fb..c5579fee 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -61,8 +61,9 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
                  pesticide_duration = 0, pesticide_efficacy = 1.0,
                  random_seed = NULL, output_frequency = "year"){ 
   
-  if (!treatment_method %in% c("ratio", "all infected")) {
-    return("treatment method is not one of the valid treatment options")
+  treatment_metric_check <- treatment_metric_checks(treatment_method)
+  if (!treatment_metric_check$checks_passed) {
+    return(treatment_metric_check$failed_check)
   }
   
   time_check <- time_checks(end_date, start_date, time_step, output_frequency)
diff --git a/R/validate.R b/R/validate.R
index 4065c8f4..c5efb8a3 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -73,18 +73,16 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
                      mask = NULL, success_metric = "quantity", output_frequency = "year"
                      ){ 
   
-  if (success_metric == "quantity") {
-    configuration = FALSE
-  } else if (success_metric == "quantity and configuration") {
-    configuration = TRUE
-  } else if (success_metric == "odds_ratio") {
-    configuration = FALSE
+  metric_check <- metric_checks(success_metric)
+  if (metric_check$checks_passed){
+    configuration <- metric_check$configuration
   } else {
-    return("Success metric must be one of 'quantity', 'quantity and configuration', or 'odds_ratio'")
+    return(metric_check$failed_check)
   }
   
-  if (!treatment_method %in% c("ratio", "all infected")) {
-    return("treatment method is not one of the valid treatment options")
+  treatment_metric_check <- treatment_metric_checks(treatment_method)
+  if (!treatment_metric_check$checks_passed) {
+    return(treatment_metric_check$failed_check)
   }
   
   time_check <- time_checks(end_date, start_date, time_step, output_frequency)

From f4c5c60e9517b2d92e2f47456031dea4919ae420 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 13:43:04 -0500
Subject: [PATCH 10/25] added initial raster and secondary raster checks

---
 R/checks.R | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 58 insertions(+), 4 deletions(-)

diff --git a/R/checks.R b/R/checks.R
index 7dc72f67..bf2e045a 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -1,32 +1,86 @@
 ## These functions are designed to improve data format checks and reduce copy and pasting of code across functions
 
-initial_raster_check <- function(x) {
+initial_raster_checks <- function(x) {
   checks_passed <- TRUE
   
-  if (!file.exists(x)) {
+  if (!all(file.exists(x))) {
     checks_passed <- FALSE
     failed_check <- "file does not exist" 
   }
   
-  if (!(raster::extension(x) %in% c(".grd", ".tif", ".img"))) {
+  if (!all((raster::extension(x) %in% c(".grd", ".tif", ".img")))) {
     checks_passed <- FALSE
     failed_check <- "file is not one of '.grd', '.tif', '.img'"
   }
   
   if (checks_passed) {
-    r<- raster::raster(x)
+    r<- raster::stack(x)
     r <- raster::reclassify(r, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   }
   
+  if (raster::nlayers(r) > 1) {
+    r <- output_from_raster_mean_and_sd(r)
+  }
+  
   if (checks_passed) {
     outs <- list(checks_passed, r)
+    names(outs) <- c('checks_passed', 'raster')
     return(outs)
   } else {
     outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
     return(outs)
   }
 }
 
+## adds checks to test for raster extent, resolution, and crs x2 being the raster already through initial checks for comparison
+secondary_raster_checks <- function(x, x2) {
+  checks_passed <- TRUE
+  
+  if (!all(file.exists(x))) {
+    checks_passed <- FALSE
+    failed_check <- "file does not exist" 
+  }
+  
+  if (!all((raster::extension(x) %in% c(".grd", ".tif", ".img")))) {
+    checks_passed <- FALSE
+    failed_check <- "file is not one of '.grd', '.tif', '.img'"
+  }
+  
+  if (checks_passed) {
+    r<- raster::stack(x)
+    r <- raster::reclassify(r, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+  }
+  
+  if (raster::nlayers(r) > 1) {
+    r <- output_from_raster_mean_and_sd(r)
+  }
+  
+  if (!(raster::extent(x2) == raster::extent(r))) {
+    checks_passed <- FALSE
+    failed_check <- "Extents of input rasters do not match. Ensure that all of your input rasters have the same extent"
+  }
+  
+  if (!(raster::xres(x2) == raster::xres(r) && raster::yres(x2) == raster::yres(r))) {
+    checks_passed <- FALSE
+    failed_check <- "Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution"
+  }
+  
+  if (!raster::compareCRS(r,x2)) {
+    checks_passed <- FALSE
+    failed_check <- "Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs"
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, r)
+    names(outs) <- c('checks_passed', 'raster')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+}
 
 
 treatment_metric_checks <- function(treatment_method) {

From 19100a429df174d502d6a54832e3b0a9ff5069da Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 13:43:42 -0500
Subject: [PATCH 11/25] add initial raster checks to wrappers

---
 R/calibrate.R     | 13 +++++--------
 R/pops.r          | 28 +++++++++++-----------------
 R/pops_multirun.R | 13 +++++--------
 R/validate.R      | 13 +++++--------
 4 files changed, 26 insertions(+), 41 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index ed39860f..42bdce49 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -175,12 +175,11 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     return("Incorrect format for prior percent natural distance")
   }
   
-  if (!file.exists(infected_file)) {
-    return("Infected file does not exist") 
-  }
-  
-  if (!(raster::extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Infected file is not one of '.grd', '.tif', '.img'")
+  infected_check <- initial_raster_checks(infected_file)
+  if (infected_check$checks_passed) {
+    infected <- infected_check$raster
+  } else {
+    return(infected_check$failed_check)
   }
   
   if (!file.exists(host_file)) {
@@ -199,8 +198,6 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  infected <- raster::raster(infected_file)
-  infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
   host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   total_plants <- raster::raster(total_plants_file)
diff --git a/R/pops.r b/R/pops.r
index 1a4adbf8..8ebb2438 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -113,12 +113,15 @@ pops <- function(infected_file, host_file, total_plants_file,
     return(percent_check$failed_check)
   }
   
-  if (!file.exists(infected_file)) {
-    return("Infected file does not exist") 
+  if (is.null(random_seed)) {
+    random_seed = round(stats::runif(1, 1, 1000000))
   }
   
-  if (!(raster::extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Infected file is not one of '.grd', '.tif', '.img'")
+  infected_check <- initial_raster_checks(infected_file)
+  if (infected_check$checks_passed) {
+    infected <- infected_check$raster
+  } else {
+    return(infected_check$failed_check)
   }
   
   if (!file.exists(host_file)) {
@@ -129,6 +132,9 @@ pops <- function(infected_file, host_file, total_plants_file,
     return("Host file is not one of '.grd', '.tif', '.img'")
   }
   
+  host <- raster::raster(host_file)
+  host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+  
   if (!file.exists(total_plants_file)) {
     return("Total plants file does not exist") 
   }
@@ -137,14 +143,6 @@ pops <- function(infected_file, host_file, total_plants_file,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  if (is.null(random_seed)) {
-    random_seed = round(stats::runif(1, 1, 1000000))
-  }
-  
-  infected <- raster::raster(infected_file)
-  infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  host <- raster::raster(host_file)
-  host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   total_plants <- raster::raster(total_plants_file)
   total_plants <- raster::reclassify(total_plants, matrix(c(NA, 0), ncol = 2, byrow = TRUE), right = NA)
 
@@ -269,7 +267,6 @@ pops <- function(infected_file, host_file, total_plants_file,
   }
   
   if (weather == TRUE){
-    # weather_coefficient_stack[is.na(weather_coefficient_stack)] <- 0
     weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
     weather_coefficient <- list(raster::as.matrix(weather_coefficient_stack[[1]]))
     for(i in 2:number_of_time_steps) {
@@ -292,8 +289,7 @@ pops <- function(infected_file, host_file, total_plants_file,
   if (management == TRUE) {
     treatment_stack <- raster::stack(treatments_file)
     treatment_stack <- raster::reclassify(treatment_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-    # treatment_stack[is.na(treatment_stack)] <- 0
-    
+
     if (!(raster::extent(infected) == raster::extent(treatment_stack))) {
       return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
     }
@@ -330,8 +326,6 @@ pops <- function(infected_file, host_file, total_plants_file,
         }
       }
     }
-    
-    # treatment_dates <- treatment_dates
   } else {
     treatment_map <- host
     raster::values(treatment_map) <- 0
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index c5579fee..9ed52d82 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -82,12 +82,11 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     return(percent_check$failed_check)
   }
   
-  if (!file.exists(infected_file)) {
-    return("Infected file does not exist") 
-  }
-  
-  if (!(raster::extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Infected file is not one of '.grd', '.tif', '.img'")
+  infected_check <- initial_raster_checks(infected_file)
+  if (infected_check$checks_passed) {
+    infected <- infected_check$raster
+  } else {
+    return(infected_check$failed_check)
   }
   
   if (!file.exists(host_file)) {
@@ -106,8 +105,6 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  infected <- raster::raster(infected_file)
-  infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
   host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   total_plants <- raster::raster(total_plants_file)
diff --git a/R/validate.R b/R/validate.R
index c5efb8a3..07b0681e 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -101,12 +101,11 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     return(percent_check$failed_check)
   }
   
-  if (!file.exists(infected_file)) {
-    return("Infected file does not exist") 
-  }
-  
-  if (!(raster::extension(infected_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Infected file is not one of '.grd', '.tif', '.img'")
+  infected_check <- initial_raster_checks(infected_file)
+  if (infected_check$checks_passed) {
+    infected <- infected_check$raster
+  } else {
+    return(infected_check$failed_check)
   }
   
   if (!file.exists(host_file)) {
@@ -125,8 +124,6 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     return("Total plants file is not one of '.grd', '.tif', '.img'")
   }
   
-  infected <- raster::raster(infected_file)
-  infected <- raster::reclassify(infected, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   host <- raster::raster(host_file)
   host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   total_plants <- raster::raster(total_plants_file)

From b6590eca2320c8ec66dce4e56c8ada60dd6be0ba Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 13:52:57 -0500
Subject: [PATCH 12/25] implemented secondary checks in wrapper functions

---
 R/calibrate.R     | 40 ++++++++++------------------------------
 R/pops.r          | 41 ++++++++++-------------------------------
 R/pops_multirun.R | 40 ++++++++++------------------------------
 R/validate.R      | 40 ++++++++++------------------------------
 4 files changed, 40 insertions(+), 121 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 42bdce49..8971f7c4 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -182,41 +182,21 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     return(infected_check$failed_check)
   }
   
-  if (!file.exists(host_file)) {
-    return("Host file does not exist") 
-  }
-  
-  if (!(raster::extension(host_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Host file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (!file.exists(total_plants_file)) {
-    return("Total plants file does not exist") 
-  }
-  
-  if (!(raster::extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Total plants file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  host <- raster::raster(host_file)
-  host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  total_plants <- raster::raster(total_plants_file)
-  total_plants <- raster::reclassify(total_plants, matrix(c(NA, 0), ncol = 2, byrow = TRUE), right = NA)
-  
-  if (!(extent(infected) == extent(host) && extent(infected) == extent(total_plants))) {
-    return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-  }
-  
-  if (!(xres(infected) == xres(host) && xres(infected) == xres(total_plants) && yres(infected) == yres(host) && yres(infected) == yres(total_plants))) {
-    return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
+  host_check <- secondary_raster_checks(host_file, infected)
+  if (host_check$checks_passed) {
+    host <- host_check$raster
+  } else {
+    return(host_check$failed_check)
   }
   
-  if (!(compareCRS(host,infected) && compareCRS(host, total_plants))) {
-    return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+  total_plants_check <- secondary_raster_checks(total_plants_file, infected)
+  if (total_plants_check$checks_passed) {
+    total_plants <- total_plants_check$raster
+  } else {
+    return(total_plants_check$failed_check)
   }
   
   susceptible <- host - infected
-  susceptible <- raster::reclassify(susceptible, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   susceptible[susceptible < 0] <- 0
   
   if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
diff --git a/R/pops.r b/R/pops.r
index 8ebb2438..499b69bc 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -124,42 +124,21 @@ pops <- function(infected_file, host_file, total_plants_file,
     return(infected_check$failed_check)
   }
   
-  if (!file.exists(host_file)) {
-    return("Host file does not exist") 
-  }
-  
-  if (!(raster::extension(host_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Host file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  host <- raster::raster(host_file)
-  host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  
-  if (!file.exists(total_plants_file)) {
-    return("Total plants file does not exist") 
-  }
-  
-  if (!(raster::extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Total plants file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  total_plants <- raster::raster(total_plants_file)
-  total_plants <- raster::reclassify(total_plants, matrix(c(NA, 0), ncol = 2, byrow = TRUE), right = NA)
-
-  if (!(raster::extent(infected) == raster::extent(host) && raster::extent(infected) == raster::extent(total_plants))) {
-    return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-  }
-  
-  if (!(raster::xres(infected) == raster::xres(host) && raster::xres(infected) == raster::xres(total_plants) && raster::yres(infected) == raster::yres(host) && raster::yres(infected) == raster::yres(total_plants))) {
-    return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
+  host_check <- secondary_raster_checks(host_file, infected)
+  if (host_check$checks_passed) {
+    host <- host_check$raster
+  } else {
+    return(host_check$failed_check)
   }
   
-  if (!(raster::compareCRS(host,infected) && raster::compareCRS(host, total_plants))) {
-    return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+  total_plants_check <- secondary_raster_checks(total_plants_file, infected)
+  if (total_plants_check$checks_passed) {
+    total_plants <- total_plants_check$raster
+  } else {
+    return(total_plants_check$failed_check)
   }
   
   susceptible <- host - infected
-  susceptible <- raster::reclassify(susceptible, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   susceptible[susceptible < 0] <- 0
   
   if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index 9ed52d82..a7775029 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -89,41 +89,21 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     return(infected_check$failed_check)
   }
   
-  if (!file.exists(host_file)) {
-    return("Host file does not exist") 
-  }
-  
-  if (!(raster::extension(host_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Host file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (!file.exists(total_plants_file)) {
-    return("Total plants file does not exist") 
-  }
-  
-  if (!(raster::extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Total plants file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  host <- raster::raster(host_file)
-  host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  total_plants <- raster::raster(total_plants_file)
-  total_plants <- raster::reclassify(total_plants, matrix(c(NA, 0), ncol = 2, byrow = TRUE), right = NA)
-  
-  if (!(raster::extent(infected) == raster::extent(host) && raster::extent(infected) == raster::extent(total_plants))) {
-    return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-  }
-  
-  if (!(raster::xres(infected) == raster::xres(host) && raster::xres(infected) == raster::xres(total_plants) && raster::yres(infected) == raster::yres(host) && raster::yres(infected) == raster::yres(total_plants))) {
-    return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
+  host_check <- secondary_raster_checks(host_file, infected)
+  if (host_check$checks_passed) {
+    host <- host_check$raster
+  } else {
+    return(host_check$failed_check)
   }
   
-  if (!(raster::compareCRS(host,infected) && raster::compareCRS(host, total_plants))) {
-    return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+  total_plants_check <- secondary_raster_checks(total_plants_file, infected)
+  if (total_plants_check$checks_passed) {
+    total_plants <- total_plants_check$raster
+  } else {
+    return(total_plants_check$failed_check)
   }
   
   susceptible <- host - infected
-  susceptible <- raster::reclassify(susceptible, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   susceptible[susceptible < 0] <- 0
   
   if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
diff --git a/R/validate.R b/R/validate.R
index 07b0681e..568d8246 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -108,41 +108,21 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     return(infected_check$failed_check)
   }
   
-  if (!file.exists(host_file)) {
-    return("Host file does not exist") 
-  }
-  
-  if (!(raster::extension(host_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Host file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (!file.exists(total_plants_file)) {
-    return("Total plants file does not exist") 
-  }
-  
-  if (!(raster::extension(total_plants_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Total plants file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  host <- raster::raster(host_file)
-  host <- raster::reclassify(host, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-  total_plants <- raster::raster(total_plants_file)
-  total_plants <- raster::reclassify(total_plants, matrix(c(NA, 0), ncol = 2, byrow = TRUE), right = NA)
-  
-  if (!(raster::extent(infected) == raster::extent(host) && raster::extent(infected) == raster::extent(total_plants))) {
-    return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-  }
-  
-  if (!(raster::xres(infected) == raster::xres(host) && raster::xres(infected) == raster::xres(total_plants) && raster::yres(infected) == raster::yres(host) && raster::yres(infected) == raster::yres(total_plants))) {
-    return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
+  host_check <- secondary_raster_checks(host_file, infected)
+  if (host_check$checks_passed) {
+    host <- host_check$raster
+  } else {
+    return(host_check$failed_check)
   }
   
-  if (!(raster::compareCRS(host,infected) && raster::compareCRS(host, total_plants))) {
-    return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+  total_plants_check <- secondary_raster_checks(total_plants_file, infected)
+  if (total_plants_check$checks_passed) {
+    total_plants <- total_plants_check$raster
+  } else {
+    return(total_plants_check$failed_check)
   }
   
   susceptible <- host - infected
-  susceptible <- raster::reclassify(susceptible, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   susceptible[susceptible < 0] <- 0
   
   if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {

From e7e89bb51ba495f1862922e2086f547c8aa647d3 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 14:04:46 -0500
Subject: [PATCH 13/25] implemented secondary check for lethal temperatures

---
 R/calibrate.R     | 35 ++++++++++++++---------------------
 R/checks.R        |  4 ----
 R/pops.r          | 37 +++++++++++++++----------------------
 R/pops_multirun.R | 35 ++++++++++++++---------------------
 R/validate.R      | 35 ++++++++++++++---------------------
 5 files changed, 57 insertions(+), 89 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 8971f7c4..0c944514 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -178,6 +178,9 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   infected_check <- initial_raster_checks(infected_file)
   if (infected_check$checks_passed) {
     infected <- infected_check$raster
+    if (raster::nlayers(infected) > 1) {
+      infected <- output_from_raster_mean_and_sd(infected)
+    }
   } else {
     return(infected_check$failed_check)
   }
@@ -185,6 +188,9 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   host_check <- secondary_raster_checks(host_file, infected)
   if (host_check$checks_passed) {
     host <- host_check$raster
+    if (raster::nlayers(host) > 1) {
+      host <- output_from_raster_mean_and_sd(host)
+    }
   } else {
     return(host_check$failed_check)
   }
@@ -192,6 +198,9 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   total_plants_check <- secondary_raster_checks(total_plants_file, infected)
   if (total_plants_check$checks_passed) {
     total_plants <- total_plants_check$raster
+    if (raster::nlayers(total_plants) > 1) {
+      total_plants <- output_from_raster_mean_and_sd(total_plants)
+    }
   } else {
     return(total_plants_check$failed_check)
   }
@@ -199,28 +208,12 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   susceptible <- host - infected
   susceptible[susceptible < 0] <- 0
   
-  if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
-    return("Temperature file does not exist")
-  }
-  
-  if (use_lethal_temperature == TRUE  && !(extension(temperature_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature file is not one of '.grd', '.tif', '.img'")
-  }
-  
   if (use_lethal_temperature == TRUE) {
-    temperature_stack <- raster::stack(temperature_file)
-    temperature_stack <- raster::reclassify(temperature_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-    
-    if (!(extent(infected) == extent(temperature_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(xres(infected) == xres(temperature_stack) && yres(infected) == yres(temperature_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(compareCRS(infected, temperature_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    temperature_check <- secondary_raster_checks(temperature_file, infected)
+    if (temperature_check$checks_passed) {
+      temperature_stack <- temperature_check$raster
+    } else {
+      return(temperature_check$failed_check)
     }
     
     temperature <- list(as.matrix(temperature_stack[[1]]))
diff --git a/R/checks.R b/R/checks.R
index bf2e045a..79248cb7 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -52,10 +52,6 @@ secondary_raster_checks <- function(x, x2) {
     r <- raster::reclassify(r, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   }
   
-  if (raster::nlayers(r) > 1) {
-    r <- output_from_raster_mean_and_sd(r)
-  }
-  
   if (!(raster::extent(x2) == raster::extent(r))) {
     checks_passed <- FALSE
     failed_check <- "Extents of input rasters do not match. Ensure that all of your input rasters have the same extent"
diff --git a/R/pops.r b/R/pops.r
index 499b69bc..95e20739 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -120,6 +120,9 @@ pops <- function(infected_file, host_file, total_plants_file,
   infected_check <- initial_raster_checks(infected_file)
   if (infected_check$checks_passed) {
     infected <- infected_check$raster
+    if (raster::nlayers(infected) > 1) {
+      infected <- output_from_raster_mean_and_sd(infected)
+    }
   } else {
     return(infected_check$failed_check)
   }
@@ -127,6 +130,9 @@ pops <- function(infected_file, host_file, total_plants_file,
   host_check <- secondary_raster_checks(host_file, infected)
   if (host_check$checks_passed) {
     host <- host_check$raster
+    if (raster::nlayers(host) > 1) {
+      host <- output_from_raster_mean_and_sd(host)
+    }
   } else {
     return(host_check$failed_check)
   }
@@ -134,35 +140,22 @@ pops <- function(infected_file, host_file, total_plants_file,
   total_plants_check <- secondary_raster_checks(total_plants_file, infected)
   if (total_plants_check$checks_passed) {
     total_plants <- total_plants_check$raster
+    if (raster::nlayers(total_plants) > 1) {
+      total_plants <- output_from_raster_mean_and_sd(total_plants)
+    }
   } else {
     return(total_plants_check$failed_check)
   }
   
   susceptible <- host - infected
   susceptible[susceptible < 0] <- 0
-  
-  if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
-    return("Temperature file does not exist")
-  }
-  
-  if (use_lethal_temperature == TRUE  && !(raster::extension(temperature_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (use_lethal_temperature == TRUE) {
-    temperature_stack <- raster::stack(temperature_file)
-    temperature_stack <- raster::reclassify(temperature_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
 
-    if (!(raster::extent(infected) == raster::extent(temperature_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(temperature_stack) && raster::yres(infected) == raster::yres(temperature_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, temperature_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+  if (use_lethal_temperature == TRUE) {
+    temperature_check <- secondary_raster_checks(temperature_file, infected)
+    if (temperature_check$checks_passed) {
+      temperature_stack <- temperature_check$raster
+    } else {
+      return(temperature_check$failed_check)
     }
     
     temperature <- list(raster::as.matrix(temperature_stack[[1]]))
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index a7775029..be6a86f9 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -85,6 +85,9 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
   infected_check <- initial_raster_checks(infected_file)
   if (infected_check$checks_passed) {
     infected <- infected_check$raster
+    if (raster::nlayers(infected) > 1) {
+      infected <- output_from_raster_mean_and_sd(infected)
+    }
   } else {
     return(infected_check$failed_check)
   }
@@ -92,6 +95,9 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
   host_check <- secondary_raster_checks(host_file, infected)
   if (host_check$checks_passed) {
     host <- host_check$raster
+    if (raster::nlayers(host) > 1) {
+      host <- output_from_raster_mean_and_sd(host)
+    }
   } else {
     return(host_check$failed_check)
   }
@@ -99,6 +105,9 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
   total_plants_check <- secondary_raster_checks(total_plants_file, infected)
   if (total_plants_check$checks_passed) {
     total_plants <- total_plants_check$raster
+    if (raster::nlayers(total_plants) > 1) {
+      total_plants <- output_from_raster_mean_and_sd(total_plants)
+    }
   } else {
     return(total_plants_check$failed_check)
   }
@@ -106,28 +115,12 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
   susceptible <- host - infected
   susceptible[susceptible < 0] <- 0
   
-  if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
-    return("Temperature file does not exist")
-  }
-  
-  if (use_lethal_temperature == TRUE  && !(raster::extension(temperature_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature file is not one of '.grd', '.tif', '.img'")
-  }
-  
   if (use_lethal_temperature == TRUE) {
-    temperature_stack <- raster::stack(temperature_file)
-    temperature_stack <- raster::reclassify(temperature_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-    
-    if (!(raster::extent(infected) == raster::extent(temperature_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(temperature_stack) && raster::yres(infected) == raster::yres(temperature_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, temperature_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    temperature_check <- secondary_raster_checks(temperature_file, infected)
+    if (temperature_check$checks_passed) {
+      temperature_stack <- temperature_check$raster
+    } else {
+      return(temperature_check$failed_check)
     }
     
     temperature <- list(raster::as.matrix(temperature_stack[[1]]))
diff --git a/R/validate.R b/R/validate.R
index 568d8246..6007ca05 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -104,6 +104,9 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   infected_check <- initial_raster_checks(infected_file)
   if (infected_check$checks_passed) {
     infected <- infected_check$raster
+    if (raster::nlayers(infected) > 1) {
+      infected <- output_from_raster_mean_and_sd(infected)
+    }
   } else {
     return(infected_check$failed_check)
   }
@@ -111,6 +114,9 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   host_check <- secondary_raster_checks(host_file, infected)
   if (host_check$checks_passed) {
     host <- host_check$raster
+    if (raster::nlayers(host) > 1) {
+      host <- output_from_raster_mean_and_sd(host)
+    }
   } else {
     return(host_check$failed_check)
   }
@@ -118,6 +124,9 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   total_plants_check <- secondary_raster_checks(total_plants_file, infected)
   if (total_plants_check$checks_passed) {
     total_plants <- total_plants_check$raster
+    if (raster::nlayers(total_plants) > 1) {
+      total_plants <- output_from_raster_mean_and_sd(total_plants)
+    }
   } else {
     return(total_plants_check$failed_check)
   }
@@ -125,28 +134,12 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   susceptible <- host - infected
   susceptible[susceptible < 0] <- 0
   
-  if (use_lethal_temperature == TRUE  && !file.exists(temperature_file)) {
-    return("Temperature file does not exist")
-  }
-  
-  if (use_lethal_temperature == TRUE  && !(raster::extension(temperature_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature file is not one of '.grd', '.tif', '.img'")
-  }
-  
   if (use_lethal_temperature == TRUE) {
-    temperature_stack <- raster::stack(temperature_file)
-    temperature_stack <- raster::reclassify(temperature_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-    
-    if (!(raster::extent(infected) == raster::extent(temperature_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(temperature_stack) && raster::yres(infected) == raster::yres(temperature_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, temperature_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    temperature_check <- secondary_raster_checks(temperature_file, infected)
+    if (temperature_check$checks_passed) {
+      temperature_stack <- temperature_check$raster
+    } else {
+      return(temperature_check$failed_check)
     }
     
     temperature <- list(raster::as.matrix(temperature_stack[[1]]))

From 94c649f2eb6d5b48f619d749e37a30f42b5cda55 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 14:19:52 -0500
Subject: [PATCH 14/25] updated weather to use new checks in all wrapper
 functions

---
 R/calibrate.R     | 75 ++++++++++++-----------------------------------
 R/pops.r          | 69 ++++++++++---------------------------------
 R/pops_multirun.R | 69 ++++++++++---------------------------------
 R/validate.R      | 69 ++++++++++---------------------------------
 4 files changed, 67 insertions(+), 215 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 0c944514..a7ea070f 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -226,70 +226,33 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     temperature <- list(as.matrix(temperature))
   }
   
-  if (temp == TRUE  && !file.exists(temperature_coefficient_file)) {
-    return("Temperature coefficient file does not exist")
-  }
-  
-  if (temp == TRUE  && !(extension(temperature_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (precip == TRUE  && !file.exists(precipitation_coefficient_file)) {
-    return("Precipitation coefficient file does not exist")
-  }
-  
-  if (precip == TRUE  && !(extension(precipitation_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Precipitation coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
   weather <- FALSE
   if (temp == TRUE) {
-    temperature_coefficient <- stack(temperature_coefficient_file)
-    
-    if (!(extent(infected) == extent(temperature_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(xres(infected) == xres(temperature_coefficient) && yres(infected) == yres(temperature_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(compareCRS(infected, temperature_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    temperature_coefficient_check <- secondary_raster_checks(temperature_coefficient_file, infected)
+    if (temperature_coefficient_check$checks_passed) {
+      temperature_coefficient <- temperature_coefficient_check$raster
+    } else {
+      return(temperature_coefficient_check$failed_check)
     }
     
     weather <- TRUE
     weather_coefficient_stack <- temperature_coefficient
     if (precip ==TRUE){
-      precipitation_coefficient <- stack(precipitation_coefficient_file)
-      
-      if (!(extent(infected) == extent(precipitation_coefficient))) {
-        return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-      }
-      
-      if (!(xres(infected) == xres(precipitation_coefficient) && yres(infected) == yres(precipitation_coefficient))) {
-        return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-      }
-      
-      if (!(compareCRS(infected, precipitation_coefficient))) {
-        return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+      precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+      if (precipitation_coefficient_check$checks_passed) {
+        precipitation_coefficient <- precipitation_coefficient_check$raster
+      } else {
+        return(precipitation_coefficient_check$failed_check)
       }
       
       weather_coefficient_stack <- weather_coefficient_stack * precipitation_coefficient
     }
   } else if(precip == TRUE){
-    precipitation_coefficient <- stack(precipitation_coefficient_file)
-    
-    if (!(extent(infected) == extent(precipitation_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(xres(infected) == xres(precipitation_coefficient) && yres(infected) == yres(precipitation_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(compareCRS(infected, precipitation_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+    if (precipitation_coefficient_check$checks_passed) {
+      precipitation_coefficient <- precipitation_coefficient_check$raster
+    } else {
+      return(precipitation_coefficient_check$failed_check)
     }
     
     weather <- TRUE
@@ -297,15 +260,15 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   }
   
   if (weather == TRUE){
-    weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+    # weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
     weather_coefficient <- list(raster::as.matrix(weather_coefficient_stack[[1]]))
     for(i in 2:number_of_time_steps) {
-      weather_coefficient[[i]] <- as.matrix(weather_coefficient_stack[[i]])
+      weather_coefficient[[i]] <- raster::as.matrix(weather_coefficient_stack[[i]])
     }
   } else {
     weather_coefficient <- host
-    weather_coefficient[] <- 1
-    weather_coefficient <- list(as.matrix(weather_coefficient))
+    raster::values(weather_coefficient) <- 1
+    weather_coefficient <- list(raster::as.matrix(weather_coefficient))
   }
   
   if (management == TRUE  && !file.exists(treatments_file)) {
diff --git a/R/pops.r b/R/pops.r
index 95e20739..ad3823d8 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -168,70 +168,33 @@ pops <- function(infected_file, host_file, total_plants_file,
     temperature <- list(raster::as.matrix(temperature))
   }
   
-  if (temp == TRUE  && !file.exists(temperature_coefficient_file)) {
-    return("Temperature coefficient file does not exist")
-  }
-  
-  if (temp == TRUE  && !(raster::extension(temperature_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (precip == TRUE  && !file.exists(precipitation_coefficient_file)) {
-    return("Precipitation coefficient file does not exist")
-  }
-  
-  if (precip == TRUE  && !(raster::extension(precipitation_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Precipitation coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
   weather <- FALSE
   if (temp == TRUE) {
-    temperature_coefficient <- raster::stack(temperature_coefficient_file)
-    
-    if (!(raster::extent(infected) == raster::extent(temperature_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(temperature_coefficient) && raster::yres(infected) == raster::yres(temperature_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, temperature_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    temperature_coefficient_check <- secondary_raster_checks(temperature_coefficient_file, infected)
+    if (temperature_coefficient_check$checks_passed) {
+      temperature_coefficient <- temperature_coefficient_check$raster
+    } else {
+      return(temperature_coefficient_check$failed_check)
     }
     
     weather <- TRUE
     weather_coefficient_stack <- temperature_coefficient
     if (precip ==TRUE){
-      precipitation_coefficient <- raster::stack(precipitation_coefficient_file)
-      
-      if (!(raster::extent(infected) == raster::extent(precipitation_coefficient))) {
-        return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-      }
-      
-      if (!(raster::xres(infected) == raster::xres(precipitation_coefficient) && raster::yres(infected) == raster::yres(precipitation_coefficient))) {
-        return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-      }
-      
-      if (!(raster::compareCRS(infected, precipitation_coefficient))) {
-        return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+      precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+      if (precipitation_coefficient_check$checks_passed) {
+        precipitation_coefficient <- precipitation_coefficient_check$raster
+      } else {
+        return(precipitation_coefficient_check$failed_check)
       }
       
       weather_coefficient_stack <- weather_coefficient_stack * precipitation_coefficient
     }
   } else if(precip == TRUE){
-    precipitation_coefficient <- raster::stack(precipitation_coefficient_file)
-    
-    if (!(raster::extent(infected) == raster::extent(precipitation_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(precipitation_coefficient) && raster::yres(infected) == raster::yres(precipitation_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, precipitation_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+    if (precipitation_coefficient_check$checks_passed) {
+      precipitation_coefficient <- precipitation_coefficient_check$raster
+    } else {
+      return(precipitation_coefficient_check$failed_check)
     }
     
     weather <- TRUE
@@ -239,7 +202,7 @@ pops <- function(infected_file, host_file, total_plants_file,
   }
   
   if (weather == TRUE){
-    weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+    # weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
     weather_coefficient <- list(raster::as.matrix(weather_coefficient_stack[[1]]))
     for(i in 2:number_of_time_steps) {
       weather_coefficient[[i]] <- raster::as.matrix(weather_coefficient_stack[[i]])
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index be6a86f9..0ad50122 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -133,70 +133,33 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     temperature <- list(raster::as.matrix(temperature))
   }
   
-  if (temp == TRUE  && !file.exists(temperature_coefficient_file)) {
-    return("Temperature coefficient file does not exist")
-  }
-  
-  if (temp == TRUE  && !(raster::extension(temperature_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (precip == TRUE  && !file.exists(precipitation_coefficient_file)) {
-    return("Precipitation coefficient file does not exist")
-  }
-  
-  if (precip == TRUE  && !(raster::extension(precipitation_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Precipitation coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
   weather <- FALSE
   if (temp == TRUE) {
-    temperature_coefficient <- raster::stack(temperature_coefficient_file)
-    
-    if (!(raster::extent(infected) == raster::extent(temperature_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(temperature_coefficient) && raster::yres(infected) == raster::yres(temperature_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, temperature_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    temperature_coefficient_check <- secondary_raster_checks(temperature_coefficient_file, infected)
+    if (temperature_coefficient_check$checks_passed) {
+      temperature_coefficient <- temperature_coefficient_check$raster
+    } else {
+      return(temperature_coefficient_check$failed_check)
     }
     
     weather <- TRUE
     weather_coefficient_stack <- temperature_coefficient
     if (precip ==TRUE){
-      precipitation_coefficient <- raster::stack(precipitation_coefficient_file)
-      
-      if (!(raster::extent(infected) == raster::extent(precipitation_coefficient))) {
-        return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-      }
-      
-      if (!(raster::xres(infected) == raster::xres(precipitation_coefficient) && raster::yres(infected) == raster::yres(precipitation_coefficient))) {
-        return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-      }
-      
-      if (!(raster::compareCRS(infected, precipitation_coefficient))) {
-        return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+      precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+      if (precipitation_coefficient_check$checks_passed) {
+        precipitation_coefficient <- precipitation_coefficient_check$raster
+      } else {
+        return(precipitation_coefficient_check$failed_check)
       }
       
       weather_coefficient_stack <- weather_coefficient_stack * precipitation_coefficient
     }
   } else if(precip == TRUE){
-    precipitation_coefficient <- raster::stack(precipitation_coefficient_file)
-    
-    if (!(raster::extent(infected) == raster::extent(precipitation_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(precipitation_coefficient) && raster::yres(infected) == raster::yres(precipitation_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, precipitation_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+    if (precipitation_coefficient_check$checks_passed) {
+      precipitation_coefficient <- precipitation_coefficient_check$raster
+    } else {
+      return(precipitation_coefficient_check$failed_check)
     }
     
     weather <- TRUE
@@ -204,7 +167,7 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
   }
   
   if (weather == TRUE){
-    weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+    # weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
     weather_coefficient <- list(raster::as.matrix(weather_coefficient_stack[[1]]))
     for(i in 2:number_of_time_steps) {
       weather_coefficient[[i]] <- raster::as.matrix(weather_coefficient_stack[[i]])
diff --git a/R/validate.R b/R/validate.R
index 6007ca05..c453d446 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -152,70 +152,33 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     temperature <- list(raster::as.matrix(temperature))
   }
   
-  if (temp == TRUE  && !file.exists(temperature_coefficient_file)) {
-    return("Temperature coefficient file does not exist")
-  }
-  
-  if (temp == TRUE  && !(raster::extension(temperature_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Temperature coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  if (precip == TRUE  && !file.exists(precipitation_coefficient_file)) {
-    return("Precipitation coefficient file does not exist")
-  }
-  
-  if (precip == TRUE  && !(raster::extension(precipitation_coefficient_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Precipitation coefficient file is not one of '.grd', '.tif', '.img'")
-  }
-  
   weather <- FALSE
   if (temp == TRUE) {
-    temperature_coefficient <- raster::stack(temperature_coefficient_file)
-    
-    if (!(raster::extent(infected) == raster::extent(temperature_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(temperature_coefficient) && raster::yres(infected) == raster::yres(temperature_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, temperature_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    temperature_coefficient_check <- secondary_raster_checks(temperature_coefficient_file, infected)
+    if (temperature_coefficient_check$checks_passed) {
+      temperature_coefficient <- temperature_coefficient_check$raster
+    } else {
+      return(temperature_coefficient_check$failed_check)
     }
     
     weather <- TRUE
     weather_coefficient_stack <- temperature_coefficient
     if (precip ==TRUE){
-      precipitation_coefficient <- raster::stack(precipitation_coefficient_file)
-      
-      if (!(raster::extent(infected) == raster::extent(precipitation_coefficient))) {
-        return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-      }
-      
-      if (!(raster::xres(infected) == raster::xres(precipitation_coefficient) && raster::yres(infected) == raster::yres(precipitation_coefficient))) {
-        return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-      }
-      
-      if (!(raster::compareCRS(infected, precipitation_coefficient))) {
-        return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+      precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+      if (precipitation_coefficient_check$checks_passed) {
+        precipitation_coefficient <- precipitation_coefficient_check$raster
+      } else {
+        return(precipitation_coefficient_check$failed_check)
       }
       
       weather_coefficient_stack <- weather_coefficient_stack * precipitation_coefficient
     }
   } else if(precip == TRUE){
-    precipitation_coefficient <- raster::stack(precipitation_coefficient_file)
-    
-    if (!(raster::extent(infected) == raster::extent(precipitation_coefficient))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(precipitation_coefficient) && raster::yres(infected) == raster::yres(precipitation_coefficient))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, precipitation_coefficient))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+    precipitation_coefficient_check <- secondary_raster_checks(precipitation_coefficient_file, infected)
+    if (precipitation_coefficient_check$checks_passed) {
+      precipitation_coefficient <- precipitation_coefficient_check$raster
+    } else {
+      return(precipitation_coefficient_check$failed_check)
     }
     
     weather <- TRUE
@@ -223,7 +186,7 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   }
   
   if (weather == TRUE){
-    weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+    # weather_coefficient_stack <- raster::reclassify(weather_coefficient_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
     weather_coefficient <- list(raster::as.matrix(weather_coefficient_stack[[1]]))
     for(i in 2:number_of_time_steps) {
       weather_coefficient[[i]] <- raster::as.matrix(weather_coefficient_stack[[i]])

From 56bfe45ebc65b58334f90d36b8a0f9d16d83ffb6 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Thu, 19 Dec 2019 14:44:37 -0500
Subject: [PATCH 15/25] added treatment checks and implemented them in wrapper
 functions

---
 R/calibrate.R     | 51 ++++++++--------------------------------------
 R/checks.R        | 40 ++++++++++++++++++++++++++++++++++++
 R/pops.r          | 52 ++++++++---------------------------------------
 R/pops_multirun.R | 52 ++++++++---------------------------------------
 R/validate.R      | 52 ++++++++---------------------------------------
 5 files changed, 76 insertions(+), 171 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index a7ea070f..0fc22aa3 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -271,52 +271,19 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
     weather_coefficient <- list(raster::as.matrix(weather_coefficient))
   }
   
-  if (management == TRUE  && !file.exists(treatments_file)) {
-    return("Treatments file does not exist")
-  }
-  
-  if (management == TRUE  && !(extension(treatments_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Treatments file is not one of '.grd', '.tif', '.img'")
-  }
-  
   if (management == TRUE) {
-    treatment_stack <- raster::stack(treatments_file)
-    treatment_stack <- raster::reclassify(treatment_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-    
-    if (!(raster::extent(infected) == raster::extent(treatment_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(treatment_stack) && raster::yres(infected) == raster::yres(treatment_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, treatment_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
-    }
-    
-    if (length(treatments_file) != length(treatment_dates)) {
-      return("Length of list for treatment dates and treatments_file must be equal")
-    }
-    
-    if (length(pesticide_duration) != length(treatment_dates)) {
-      return("Length of list for treatment dates and pesticide_duration must be equal")
+    treatments_check <- secondary_raster_checks(treatments_file, infected)
+    if (treatments_check$checks_passed) {
+      treatment_stack <- treatments_check$raster
+    } else {
+      return(treatments_check$failed_check)
     }
     
-    if (pesticide_duration[1] > 0) {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]] * pesticide_efficacy))
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    if (treatment_check$checks_passed) {
+      treatment_maps <- treatment_check$treatment_maps
     } else {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]]))
-    }
-    if (raster::nlayers(treatment_stack) >= 2) {
-      for(i in 2:raster::nlayers(treatment_stack)) {
-        if (pesticide_duration[i] > 0) {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]] * pesticide_efficacy)
-        } else {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]])
-          
-        }
-      }
+      return(treatment_check$failed_check)
     }
   } else {
     treatment_map <- host
diff --git a/R/checks.R b/R/checks.R
index 79248cb7..fdc4b880 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -78,6 +78,46 @@ secondary_raster_checks <- function(x, x2) {
   }
 }
 
+treatment_checks <- function(treatment_stack, treatments_file, pesticide_duration, treatment_dates) {
+  checks_passed <- TRUE
+  
+  if (length(treatments_file) != length(treatment_dates)) {
+    checks_passed <- FALSE
+    failed_check <- "Length of list for treatment dates and treatments_file must be equal"
+  }
+  
+  if (length(pesticide_duration) != length(treatment_dates)) {
+    checks_passed <- FALSE
+    failed_check <- "Length of list for treatment dates and pesticide_duration must be equal"
+  }
+  
+  if (pesticide_duration[1] > 0) {
+    treatment_maps <- list(raster::as.matrix(treatment_stack[[1]] * pesticide_efficacy))
+  } else {
+    treatment_maps <- list(raster::as.matrix(treatment_stack[[1]]))
+  }
+  
+  if (raster::nlayers(treatment_stack) >= 2) {
+    for(i in 2:raster::nlayers(treatment_stack)) {
+      if (pesticide_duration[i] > 0) {
+        treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]] * pesticide_efficacy)
+      } else {
+        treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]])
+        
+      }
+    }
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, treatment_maps)
+    names(outs) <- c('checks_passed', 'treatment_maps')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+}
 
 treatment_metric_checks <- function(treatment_method) {
   checks_passed <- TRUE
diff --git a/R/pops.r b/R/pops.r
index ad3823d8..65db2400 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -213,53 +213,19 @@ pops <- function(infected_file, host_file, total_plants_file,
     weather_coefficient <- list(raster::as.matrix(weather_coefficient))
   }
   
-  if (management == TRUE  && !file.exists(treatments_file)) {
-    return("Treatments file does not exist")
-  }
-  
-  if (management == TRUE  && !(raster::extension(treatments_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Treatments file is not one of '.grd', '.tif', '.img'")
-  }
-  
   if (management == TRUE) {
-    treatment_stack <- raster::stack(treatments_file)
-    treatment_stack <- raster::reclassify(treatment_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-
-    if (!(raster::extent(infected) == raster::extent(treatment_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(treatment_stack) && raster::yres(infected) == raster::yres(treatment_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, treatment_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
-    }
-    
-    if (length(treatments_file) != length(treatment_dates)) {
-      return("Length of list for treatment dates and treatments_file must be equal")
-    }
-    
-    if (length(pesticide_duration) != length(treatment_dates)) {
-      return("Length of list for treatment dates and pesticide_duration must be equal")
-    }
-    
-    if (pesticide_duration[1] > 0) {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]] * pesticide_efficacy))
+    treatments_check <- secondary_raster_checks(treatments_file, infected)
+    if (treatments_check$checks_passed) {
+      treatment_stack <- treatments_check$raster
     } else {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]]))
+      return(treatments_check$failed_check)
     }
     
-    if (raster::nlayers(treatment_stack) >= 2) {
-      for(i in 2:raster::nlayers(treatment_stack)) {
-        if (pesticide_duration[i] > 0) {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]] * pesticide_efficacy)
-        } else {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]])
-          
-        }
-      }
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    if (treatment_check$checks_passed) {
+      treatment_maps <- treatment_check$treatment_maps
+    } else {
+      return(treatment_check$failed_check)
     }
   } else {
     treatment_map <- host
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index 0ad50122..17766301 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -178,53 +178,19 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
     weather_coefficient <- list(raster::as.matrix(weather_coefficient))
   }
   
-  if (management == TRUE  && !file.exists(treatments_file)) {
-    return("Treatments file does not exist")
-  }
-  
-  if (management == TRUE  && !(raster::extension(treatments_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Treatments file is not one of '.grd', '.tif', '.img'")
-  }
-  
   if (management == TRUE) {
-    treatment_stack <- raster::stack(treatments_file)
-    treatment_stack <- raster::reclassify(treatment_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-
-    if (!(raster::extent(infected) == raster::extent(treatment_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(treatment_stack) && raster::yres(infected) == raster::yres(treatment_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, treatment_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
-    }
-    
-    if (length(treatments_file) != length(treatment_dates)) {
-      return("Length of list for treatment dates and treatments_file must be equal")
-    }
-    
-    if (length(pesticide_duration) != length(treatment_dates)) {
-      return("Length of list for treatment dates and pesticide_duration must be equal")
-    }
-    
-    if (pesticide_duration[1] > 0) {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]] * pesticide_efficacy))
+    treatments_check <- secondary_raster_checks(treatments_file, infected)
+    if (treatments_check$checks_passed) {
+      treatment_stack <- treatments_check$raster
     } else {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]]))
+      return(treatments_check$failed_check)
     }
     
-    if (raster::nlayers(treatment_stack) >= 2) {
-      for(i in 2:raster::nlayers(treatment_stack)) {
-        if (pesticide_duration[i] > 0) {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]] * pesticide_efficacy)
-        } else {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]])
-          
-        }
-      }
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    if (treatment_check$checks_passed) {
+      treatment_maps <- treatment_check$treatment_maps
+    } else {
+      return(treatment_check$failed_check)
     }
   } else {
     treatment_map <- host
diff --git a/R/validate.R b/R/validate.R
index c453d446..c6220648 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -197,53 +197,19 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     weather_coefficient <- list(raster::as.matrix(weather_coefficient))
   }
   
-  if (management == TRUE  && !file.exists(treatments_file)) {
-    return("Treatments file does not exist")
-  }
-  
-  if (management == TRUE  && !(raster::extension(treatments_file) %in% c(".grd", ".tif", ".img"))) {
-    return("Treatments file is not one of '.grd', '.tif', '.img'")
-  }
-  
-  
   if (management == TRUE) {
-    treatment_stack <- raster::stack(treatments_file)
-    treatment_stack <- raster::reclassify(treatment_stack, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
-    
-    if (!(raster::extent(infected) == raster::extent(treatment_stack))) {
-      return("Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
-    }
-    
-    if (!(raster::xres(infected) == raster::xres(treatment_stack) && raster::yres(infected) == raster::yres(treatment_stack))) {
-      return("Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
-    }
-    
-    if (!(raster::compareCRS(infected, treatment_stack))) {
-      return("Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
-    }
-    
-    if (length(treatments_file) != length(treatment_dates)) {
-      return("Length of list for treatment dates and treatments_file must be equal")
-    }
-    
-    if (length(pesticide_duration) != length(treatment_dates)) {
-      return("Length of list for treatment dates and pesticide_duration must be equal")
+    treatments_check <- secondary_raster_checks(treatments_file, infected)
+    if (treatments_check$checks_passed) {
+      treatment_stack <- treatments_check$raster
+    } else {
+      return(treatments_check$failed_check)
     }
     
-    if (pesticide_duration[1] > 0) {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]] * pesticide_efficacy))
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    if (treatment_check$checks_passed) {
+      treatment_maps <- treatment_check$treatment_maps
     } else {
-      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]]))
-    }
-    if (raster::nlayers(treatment_stack) >= 2) {
-      for(i in 2:raster::nlayers(treatment_stack)) {
-        if (pesticide_duration[i] > 0) {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]] * pesticide_efficacy)
-        } else {
-          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]])
-          
-        }
-      }
+      return(treatment_check$failed_check)
     }
   } else {
     treatment_map <- host

From 88a00a92478f271431a7b42e18ff0ecbd882bd21 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Fri, 20 Dec 2019 07:57:31 -0500
Subject: [PATCH 16/25] added function calculating calibrated and posterior
 probabilities for output variables.

---
 R/calibrate.R | 275 ++++++++++----------------------------------------
 R/checks.R    |  92 +++++++++++++++++
 2 files changed, 144 insertions(+), 223 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 0fc22aa3..0d62034d 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -80,100 +80,44 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   }
   
   ## Setup for reproductive rate to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
-  if (class(prior_reproductive_rate) == "numeric" && length(prior_reproductive_rate) == 2) {
-    prior_reproductive_rate <- matrix(prior_reproductive_rate, ncol = 2)
-  } 
-  
-  if (class(prior_reproductive_rate) %in% c("matrix", "data.frame") && ncol(prior_reproductive_rate) == 2) {
-    if (class(prior_reproductive_rate) == "matrix" && nrow(prior_reproductive_rate) == 1) {
-      start_reproductive_rate <- prior_reproductive_rate[1]
-      sd_reproductive_rate <- prior_reproductive_rate[2]
-    } else if (class(prior_reproductive_rate) == "data.frame" && nrow(prior_reproductive_rate) == 1) {
-      start_reproductive_rate <- prior_reproductive_rate[[1]]
-      sd_reproductive_rate <- 0
-    } else if (class(prior_reproductive_rate) %in% c("matrix", "data.frame") && nrow(prior_reproductive_rate) > 1) {
-      names(prior_reproductive_rate) <- c('var', 'prob')
-      start_reproductive_rate <- prior_reproductive_rate$var[prior_reproductive_rate$prob == max(prior_reproductive_rate$prob)]
-      if(length(start_reproductive_rate) > 1) {
-        start_reproductive_rate <- mean(start_reproductive_rate)
-      }
-      sd_reproductive_rate <- sd(prior_reproductive_rate$var)
-    }
+  prior_reproductive_rate_check <- prior_checks(prior_reproductive_rate)
+  if(prior_reproductive_rate_check$checks_passed) {
+    prior_reproductive_rate <- prior_reproductive_rate_check$priors
+    start_reproductive_rate <- prior_reproductive_rate_check$start_priors
+    sd_reproductive_rate <- prior_reproductive_rate_check$sd_priors
   } else {
-    return("Incorrect format for prior reproductive rate")
+    return(prior_reproductive_rate_check$failed_check)
   }
   
   ## Setup for natural dispersal scale to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
-  if (class(prior_natural_distance_scale) == "numeric" && length(prior_natural_distance_scale) == 2) {
-    prior_natural_distance_scale <- matrix(prior_natural_distance_scale, ncol = 2)
-  } 
-  
-  if (class(prior_natural_distance_scale) %in% c("matrix", "data.frame") && ncol(prior_natural_distance_scale) == 2) {
-    if (class(prior_natural_distance_scale) == "matrix" && nrow(prior_natural_distance_scale) == 1) {
-      start_natural_distance_scale <- prior_natural_distance_scale[1]
-      sd_natural_distance_scale <- prior_natural_distance_scale[2]
-    } else if (class(prior_natural_distance_scale) == "data.frame" && nrow(prior_natural_distance_scale) == 1) {
-      start_natural_distance_scale <- prior_natural_distance_scale[[1]]
-      sd_natural_distance_scale <- 0
-    } else if (class(prior_natural_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_natural_distance_scale) > 1) {
-      names(prior_natural_distance_scale) <- c('var', 'prob')
-      start_natural_distance_scale <- prior_natural_distance_scale$var[prior_natural_distance_scale$prob == max(prior_natural_distance_scale$prob)]
-      if(length(start_natural_distance_scale) > 1) {
-        start_natural_distance_scale <- mean(start_natural_distance_scale)
-      }
-      sd_natural_distance_scale <- sd(prior_natural_distance_scale$var)
-    }
+  prior_natural_distance_scale_check <- prior_checks(prior_natural_distance_scale)
+  if(prior_natural_distance_scale_check$checks_passed) {
+    prior_natural_distance_scale <- prior_natural_distance_scale_check$priors
+    start_natural_distance_scale <- prior_natural_distance_scale_check$start_priors
+    sd_natural_distance_scale <- prior_natural_distance_scale_check$sd_priors
   } else {
-    return("Incorrect format for prior natural distance scale")
+    return(prior_natural_distance_scale_check$failed_check)
   }
-  
+
   ## Setup for anthropogenic dispersal scale to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
-  if (class(prior_anthropogenic_distance_scale) == "numeric" && length(prior_anthropogenic_distance_scale) == 2) {
-    prior_anthropogenic_distance_scale <- matrix(prior_anthropogenic_distance_scale, ncol = 2)
-  } 
-  
-  if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && ncol(prior_anthropogenic_distance_scale) == 2) {
-    if (class(prior_anthropogenic_distance_scale) == "matrix" && nrow(prior_anthropogenic_distance_scale) == 1) {
-      start_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale[1]
-      sd_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale[2]
-    } else if (class(prior_anthropogenic_distance_scale) == "data.frame" && nrow(prior_anthropogenic_distance_scale) == 1) {
-      start_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale[[1]]
-      sd_anthropogenic_distance_scale <- 0
-    } else if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_anthropogenic_distance_scale) > 1) {
-      names(prior_anthropogenic_distance_scale) <- c('var', 'prob')
-      start_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale$var[prior_anthropogenic_distance_scale$prob == max(prior_anthropogenic_distance_scale$prob)]
-      if(length(start_anthropogenic_distance_scale) > 1) {
-        start_anthropogenic_distance_scale <- mean(start_anthropogenic_distance_scale)
-      }
-      sd_anthropogenic_distance_scale <- sd(prior_anthropogenic_distance_scale$var)
-    }
+  prior_anthropogenic_distance_scale_check <- prior_checks(prior_anthropogenic_distance_scale)
+  if(prior_anthropogenic_distance_scale_check$checks_passed) {
+    prior_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale_check$priors
+    start_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale_check$start_priors
+    sd_anthropogenic_distance_scale <- prior_anthropogenic_distance_scale_check$sd_priors
   } else {
-    return("Incorrect format for prior athropogenic distance scale")
+    return(prior_anthropogenic_distance_scale_check$failed_check)
   }
-    
+
   ## Setup for percent natural distance to be passed in as either mean and sd or a 2 column data.frame with value and probability as columns 1 and 2
-  if (class(prior_percent_natural_dispersal) == "numeric" && length(prior_percent_natural_dispersal) == 2) {
-    prior_percent_natural_dispersal <- matrix(prior_percent_natural_dispersal, ncol = 2)
-  } 
-    
-  if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && ncol(prior_percent_natural_dispersal) == 2) {
-    if (class(prior_percent_natural_dispersal) == "matrix" && nrow(prior_percent_natural_dispersal) == 1) {
-      start_percent_natural_dispersal <- prior_percent_natural_dispersal[1]
-      sd_percent_natural_dispersal <- prior_percent_natural_dispersal[2]
-    } else if (class(prior_percent_natural_dispersal) == "data.frame" && nrow(prior_percent_natural_dispersal) == 1) {
-      start_percent_natural_dispersal <- prior_percent_natural_dispersal[[1]]
-      sd_percent_natural_dispersal <- 0
-    } else if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && nrow(prior_percent_natural_dispersal) > 1) {
-      names(prior_percent_natural_dispersal) <- c('var', 'prob')
-      start_percent_natural_dispersal <- prior_percent_natural_dispersal$var[prior_percent_natural_dispersal$prob == max(prior_percent_natural_dispersal$prob)]
-      if(length(start_percent_natural_dispersal) > 1) {
-        start_percent_natural_dispersal <- mean(start_percent_natural_dispersal)
-      }
-      sd_percent_natural_dispersal <- sd(prior_percent_natural_dispersal$var)
-    }
+  prior_percent_natural_dispersal_check <- prior_checks(prior_percent_natural_dispersal)
+  if(prior_percent_natural_dispersal_check$checks_passed) {
+    prior_percent_natural_dispersal <- prior_percent_natural_dispersal_check$priors
+    start_percent_natural_dispersal <- prior_percent_natural_dispersal_check$start_priors
+    sd_percent_natural_dispersal <- prior_percent_natural_dispersal_check$sd_priors
   } else {
-    return("Incorrect format for prior percent natural distance")
-  }
+    return(prior_percent_natural_dispersal_check$failed_check)
+  } 
   
   infected_check <- initial_raster_checks(infected_file)
   if (infected_check$checks_passed) {
@@ -615,158 +559,43 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   ## Use prior and calibrated parameters to update to posteriors
   count <- 10000000
   # Reproductive Rate
-  if (class(prior_reproductive_rate) == "matrix" && nrow(prior_reproductive_rate) == 1) {
-    prior_reproductive_rates <- round(rnorm(count, start_reproductive_rate, sd_reproductive_rate), digits = 1)
-    prior_reproductive_rates <- as.data.frame(table(prior_reproductive_rates))
-    prior_reproductive_rates$prior_reproductive_rates <- as.numeric(as.character(prior_reproductive_rates$prior_reproductive_rates))
-    prior_reproductive_rates$prob <- round(prior_reproductive_rates$Freq/count, digits = 3)
-    prior_reproductive_rates <- prior_reproductive_rates[prior_reproductive_rates$prob > 0.00,]
-  } else if (class(prior_reproductive_rate) %in% c("matrix", "data.frame") && nrow(prior_reproductive_rate) > 1) {
-    prior_reproductive_rates <- prior_reproductive_rate[ , 1:2]
-    names(prior_reproductive_rates) <- c('prior_reproductive_rates', 'prob')
-    
+  reproductive_rate_checks <- bayesian_checks(prior_reproductive_rate, start_reproductive_rate, sd_reproductive_rate, params$reproductive_rate, count, prior_weight, weight, step_size = 0.1, bounds = c(0, Inf), round_to = 1, round_to_digits = 1)
+  if(reproductive_rate_checks$checks_passed) {
+    reproductive_rates <- reproductive_rate_checks$rates
+    posterior_reproductive_rates <- reproductive_rate_checks$posterior_rates
+  } else {
+    return(reproductive_rate_checks$failed_check)
   }
 
-  calibration_count <- nrow(params)
-  calibrated_reproductive_rates <- as.data.frame(table(params$reproductive_rate))
-  calibrated_reproductive_rates$Var1 <- as.numeric(as.character(calibrated_reproductive_rates$Var1))
-  calibrated_reproductive_rates$prob <- round(calibrated_reproductive_rates$Freq/calibration_count, digits = 3)
-  
-  min_reproductive_rate <- min(min(prior_reproductive_rates$prior_reproductive_rates), min(calibrated_reproductive_rates$Var1))
-  max_reproductive_rate <- max(max(prior_reproductive_rates$prior_reproductive_rates), max(calibrated_reproductive_rates$Var1))
-  
-  reproductive_rates <- data.frame(reproductive_rate = round(seq(min_reproductive_rate, max_reproductive_rate, 0.1), digits = 1), 
-                                   prior_probability = rep(0, length(seq(min_reproductive_rate, max_reproductive_rate, 0.1))),
-                                   calibrated_probability = rep(0, length(seq(min_reproductive_rate, max_reproductive_rate, 0.1))),
-                                   posterior_probability = rep(0, length(seq(min_reproductive_rate, max_reproductive_rate, 0.1))))
-  
-  for (i in 1:nrow(reproductive_rates)) {
-    if (length(prior_reproductive_rates$prob[prior_reproductive_rates$prior_reproductive_rates == reproductive_rates$reproductive_rate[i]]) > 0) {
-      reproductive_rates$prior_probability[i] <- prior_reproductive_rates$prob[prior_reproductive_rates$prior_reproductive_rates == reproductive_rates$reproductive_rate[i]]
-    }
-    if (length(calibrated_reproductive_rates$prob[calibrated_reproductive_rates$Var1 == reproductive_rates$reproductive_rate[i]]) > 0) {
-      reproductive_rates$calibrated_probability[i] <- calibrated_reproductive_rates$prob[calibrated_reproductive_rates$Var1 == reproductive_rates$reproductive_rate[i]]
-    }
-  }
-  reproductive_rates$posterior_probability <- round(reproductive_rates$prior_probability*prior_weight + reproductive_rates$calibrated_probability * weight, digits = 3)
-  colSums(reproductive_rates)
-  posterior_reproductive_rates <- reproductive_rates[,c(1,4)]
-  
   # Natural Distance Scale
-  if (class(prior_natural_distance_scale) == "matrix" && nrow(prior_natural_distance_scale) == 1) {
-    prior_natural_distance_scales <- round(rnorm(count, start_natural_distance_scale, sd_natural_distance_scale), digits = 0)
-    prior_natural_distance_scales <- as.data.frame(table(prior_natural_distance_scales))
-    prior_natural_distance_scales$prior_natural_distance_scales <- as.numeric(as.character(prior_natural_distance_scales$prior_natural_distance_scales))
-    prior_natural_distance_scales$prob <- round(prior_natural_distance_scales$Freq/count, digits = 3)
-    prior_natural_distance_scales <- prior_natural_distance_scales[prior_natural_distance_scales$prob > 0.000,]
-  } else if (class(prior_natural_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_natural_distance_scale) > 1) {
-    prior_natural_distance_scales <- prior_natural_distance_scale[ , 1:2]
-    names(prior_natural_distance_scales) <- c('prior_natural_distance_scales', 'prob')
-  }
-  
-  calibration_count <- nrow(params)
-  calibrated_natural_distance_scales <- as.data.frame(table(params$natural_distance_scale))
-  calibrated_natural_distance_scales$Var1 <- as.numeric(as.character(calibrated_natural_distance_scales$Var1))
-  calibrated_natural_distance_scales$prob <- round(calibrated_natural_distance_scales$Freq/calibration_count, digits = 3)
-  
-  min_natural_distance_scale <- min(min(prior_natural_distance_scales$prior_natural_distance_scales), min(calibrated_natural_distance_scales$Var1))
-  max_natural_distance_scale <- max(max(prior_natural_distance_scales$prior_natural_distance_scales), max(calibrated_natural_distance_scales$Var1))
-  
-  natural_distance_scales <- data.frame(natural_distance_scale = round(seq(min_natural_distance_scale, max_natural_distance_scale, 1), digits = 1), 
-                                   prior_probability = rep(0, length(seq(min_natural_distance_scale, max_natural_distance_scale, 1))),
-                                   calibrated_probability = rep(0, length(seq(min_natural_distance_scale, max_natural_distance_scale, 1))),
-                                   posterior_probability = rep(0, length(seq(min_natural_distance_scale, max_natural_distance_scale, 1))))
-  
-  for (i in 1:nrow(natural_distance_scales)) {
-    if (length(prior_natural_distance_scales$prob[prior_natural_distance_scales$prior_natural_distance_scales == natural_distance_scales$natural_distance_scale[i]]) > 0) {
-      natural_distance_scales$prior_probability[i] <- prior_natural_distance_scales$prob[prior_natural_distance_scales$prior_natural_distance_scales == natural_distance_scales$natural_distance_scale[i]]
-    }
-    if (length(calibrated_natural_distance_scales$prob[calibrated_natural_distance_scales$Var1 == natural_distance_scales$natural_distance_scale[i]]) > 0) {
-      natural_distance_scales$calibrated_probability[i] <- calibrated_natural_distance_scales$prob[calibrated_natural_distance_scales$Var1 == natural_distance_scales$natural_distance_scale[i]]
-    }
+  natural_distance_scale_checks <- bayesian_checks(prior_natural_distance_scale, start_natural_distance_scale, sd_natural_distance_scale, params$natural_distance_scale, count, prior_weight, weight, step_size = 1, bounds = c(0, Inf), round_to = 1, round_to_digits = 0)
+  if(natural_distance_scale_checks$checks_passed) {
+    natural_distance_scales <- natural_distance_scale_checks$rates
+    posterior_natural_distance_scales <- natural_distance_scale_checks$posterior_rates
+  } else {
+    return(natural_distance_scale_checks$failed_check)
   }
-  natural_distance_scales$posterior_probability <- round(natural_distance_scales$prior_probability*prior_weight + natural_distance_scales$calibrated_probability * weight, digits = 3)
-  colSums(natural_distance_scales)
-  posterior_natural_distance_scales <- natural_distance_scales[,c(1,4)]
   
   # Anthropogenic Distance Scale
-  if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_anthropogenic_distance_scale) == 1) {
-    prior_anthropogenic_distance_scales <- round(rnorm(count, start_anthropogenic_distance_scale, sd_anthropogenic_distance_scale)/10, digits = 0)*10
-    prior_anthropogenic_distance_scales <- as.data.frame(table(prior_anthropogenic_distance_scales))
-    prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales <- as.numeric(as.character(prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales))
-    prior_anthropogenic_distance_scales$prob <- round(prior_anthropogenic_distance_scales$Freq/count, digits = 3)
-    prior_anthropogenic_distance_scales <- prior_anthropogenic_distance_scales[prior_anthropogenic_distance_scales$prob > 0.000,]
-  } else if (class(prior_anthropogenic_distance_scale) %in% c("matrix", "data.frame") && nrow(prior_anthropogenic_distance_scale) > 1) {
-    prior_anthropogenic_distance_scales <- prior_anthropogenic_distance_scale[ , 1:2]
-    names(prior_anthropogenic_distance_scales) <- c('prior_anthropogenic_distance_scales', 'prob')
-  }
-  
-  calibration_count <- nrow(params)
-  calibrated_anthropogenic_distance_scales <- as.data.frame(table(params$anthropogenic_distance_scale))
-  calibrated_anthropogenic_distance_scales$Var1 <- as.numeric(as.character(calibrated_anthropogenic_distance_scales$Var1))
-  calibrated_anthropogenic_distance_scales$prob <- round(calibrated_anthropogenic_distance_scales$Freq/calibration_count, digits = 3)
-  
-  min_anthropogenic_distance_scale <- min(min(prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales), min(calibrated_anthropogenic_distance_scales$Var1))
-  max_anthropogenic_distance_scale <- max(max(prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales), max(calibrated_anthropogenic_distance_scales$Var1))
-  
-  anthropogenic_distance_scales <- data.frame(anthropogenic_distance_scale = round(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10), digits = 1), 
-                                        prior_probability = rep(0, length(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10))),
-                                        calibrated_probability = rep(0, length(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10))),
-                                        posterior_probability = rep(0, length(seq(min_anthropogenic_distance_scale, max_anthropogenic_distance_scale, 10))))
-  
-  for (i in 1:nrow(anthropogenic_distance_scales)) {
-    if (length(prior_anthropogenic_distance_scales$prob[prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]) > 0) {
-      anthropogenic_distance_scales$prior_probability[i] <- prior_anthropogenic_distance_scales$prob[prior_anthropogenic_distance_scales$prior_anthropogenic_distance_scales == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]
-    }
-    if (length(calibrated_anthropogenic_distance_scales$prob[calibrated_anthropogenic_distance_scales$Var1 == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]) > 0) {
-      anthropogenic_distance_scales$calibrated_probability[i] <- calibrated_anthropogenic_distance_scales$prob[calibrated_anthropogenic_distance_scales$Var1 == anthropogenic_distance_scales$anthropogenic_distance_scale[i]]
-    }
+  anthropogenic_distance_scale_checks <- bayesian_checks(prior_anthropogenic_distance_scale, start_anthropogenic_distance_scale, sd_anthropogenic_distance_scale, params$anthropogenic_distance_scale, count, prior_weight, weight, step_size = 10, bounds = c(0, Inf), round_to = 10, round_to_digits = 0)
+  if(anthropogenic_distance_scale_checks$checks_passed) {
+    anthropogenic_distance_scales <- anthropogenic_distance_scale_checks$rates
+    posterior_anthropogenic_distance_scales <- anthropogenic_distance_scale_checks$posterior_rates
+  } else {
+    return(anthropogenic_distance_scale_checks$failed_check)
   }
-  anthropogenic_distance_scales$posterior_probability <- round(anthropogenic_distance_scales$prior_probability*prior_weight + anthropogenic_distance_scales$calibrated_probability * weight, digits = 3)
-  colSums(anthropogenic_distance_scales)
-  posterior_anthropogenic_distance_scales <- anthropogenic_distance_scales[,c(1,4)]
   
   # Percent Natural Dispersal
-  if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && nrow(prior_percent_natural_dispersal) == 1) {
-    prior_percent_natural_dispersals <- round(rnorm(count, start_percent_natural_dispersal, sd_percent_natural_dispersal), digits = 3)
-    prior_percent_natural_dispersals <- as.data.frame(table(prior_percent_natural_dispersals))
-    prior_percent_natural_dispersals$prior_percent_natural_dispersals <- as.numeric(as.character(prior_percent_natural_dispersals$prior_percent_natural_dispersals))
-    prior_percent_natural_dispersals$prob <- round(prior_percent_natural_dispersals$Freq/count, digits = 3)
-    prior_percent_natural_dispersals <- prior_percent_natural_dispersals[prior_percent_natural_dispersals$prob > 0.000,]
-    prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals == 1.000] <- sum(prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals >= 1.0])
-    prior_percent_natural_dispersals <- prior_percent_natural_dispersals[prior_percent_natural_dispersals$prior_percent_natural_dispersals <= 1.0,]
-  } else if (class(prior_percent_natural_dispersal) %in% c("matrix", "data.frame") && nrow(prior_percent_natural_dispersal) > 1) {
-    prior_percent_natural_dispersals <- prior_percent_natural_dispersal[ , 1:2]
-    names(prior_percent_natural_dispersals) <- c('prior_percent_natural_dispersals', 'prob')
-  }
-  
-  calibration_count <- nrow(params)
-  calibrated_percent_natural_dispersals <- as.data.frame(table(params$percent_natural_dispersal))
-  calibrated_percent_natural_dispersals$Var1 <- as.numeric(as.character(calibrated_percent_natural_dispersals$Var1))
-  calibrated_percent_natural_dispersals$prob <- round(calibrated_percent_natural_dispersals$Freq/calibration_count, digits = 3)
-  
-  min_percent_natural_dispersal <- min(min(prior_percent_natural_dispersals$prior_percent_natural_dispersals), min(calibrated_percent_natural_dispersals$Var1))
-  max_percent_natural_dispersal <- max(max(prior_percent_natural_dispersals$prior_percent_natural_dispersals), max(calibrated_percent_natural_dispersals$Var1))
-  
-  percent_natural_dispersals <- data.frame(percent_natural_dispersal = round(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001), digits = 3), 
-                                        prior_probability = rep(0, length(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001))),
-                                        calibrated_probability = rep(0, length(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001))),
-                                        posterior_probability = rep(0, length(seq(min_percent_natural_dispersal, max_percent_natural_dispersal, 0.001))))
-  
-  for (i in 1:nrow(percent_natural_dispersals)) {
-    if (length(prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals == percent_natural_dispersals$percent_natural_dispersal[i]]) > 0) {
-      percent_natural_dispersals$prior_probability[i] <- prior_percent_natural_dispersals$prob[prior_percent_natural_dispersals$prior_percent_natural_dispersals == percent_natural_dispersals$percent_natural_dispersal[i]]
-    }
-    if (length(calibrated_percent_natural_dispersals$prob[calibrated_percent_natural_dispersals$Var1 == percent_natural_dispersals$percent_natural_dispersal[i]]) > 0) {
-      percent_natural_dispersals$calibrated_probability[i] <- calibrated_percent_natural_dispersals$prob[calibrated_percent_natural_dispersals$Var1 == percent_natural_dispersals$percent_natural_dispersal[i]]
-    }
+  percent_natural_dispersal_checks <- bayesian_checks(prior_percent_natural_dispersal, start_percent_natural_dispersal, sd_percent_natural_dispersal, params$percent_natural_dispersal, count, prior_weight, weight, step_size = 0.001, bounds = c(0, 1.000), round_to = 1, round_to_digits = 3)
+  if(percent_natural_dispersal_checks$checks_passed) {
+    percent_natural_dispersals <- percent_natural_dispersal_checks$rates
+    posterior_percent_natural_dispersals <- percent_natural_dispersal_checks$posterior_rates
+  } else {
+    return(percent_natural_dispersal_checks$failed_check)
   }
-  percent_natural_dispersals$posterior_probability <- round(percent_natural_dispersals$prior_probability*prior_weight + percent_natural_dispersals$calibrated_probability * weight, digits = 3)
-  colSums(percent_natural_dispersals)
-  posterior_percent_natural_dispersals <- percent_natural_dispersals[,c(1,4)]
   
   outputs <- list(posterior_reproductive_rates, posterior_natural_distance_scales, posterior_anthropogenic_distance_scales, posterior_percent_natural_dispersals, total_number_of_observations, reproductive_rates, natural_distance_scales, anthropogenic_distance_scales, percent_natural_dispersals, params)
   names(outputs) <- c('posterior_reproductive_rates', 'posterior_natural_distance_scales', 'posterior_anthropogenic_distance_scales', 'posterior_percent_natural_dispersals', 'total_number_of_observations', 'reproductive_rates', 'natural_distance_scales', 'anthropogenic_distance_scales', 'percent_natural_dispersals', 'raw_calibration_data')
-  
+
   return(outputs)
 }
diff --git a/R/checks.R b/R/checks.R
index fdc4b880..5be3fd44 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -256,3 +256,95 @@ time_checks <- function(end_date, start_date, time_step, output_frequency) {
   
 }
 
+## check for making sure priors are in the proper format and output the mean where the mean serves as the starting point for the calibration
+prior_checks <- function(priors) {
+  checks_passed <- TRUE
+  
+  if (class(priors) == "numeric" && length(priors) == 2) {
+    priors <- matrix(priors, ncol = 2)
+  } 
+  
+  if (class(priors) %in% c("matrix", "data.frame") && ncol(priors) == 2) {
+    if (class(priors) == "matrix" && nrow(priors) == 1) {
+      start_priors <- priors[1]
+      sd_priors <- priors[2]
+    } else if (class(priors) == "data.frame" && nrow(priors) == 1) {
+      start_priors <- priors[[1]]
+      sd_priors <- 0
+    } else if (class(priors) %in% c("matrix", "data.frame") && nrow(priors) > 1) {
+      names(priors) <- c('var', 'prob')
+      start_priors <- priors$var[priors$prob == max(priors$prob)]
+      if(length(start_priors) > 1) {
+        start_priors <- mean(start_priors)
+      }
+      sd_priors <- sd(priors$var)
+    }
+  } else {
+    checks_passed <- FALSE
+    failed_check <- "Incorrect format for priors"
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, priors, start_priors, sd_priors)
+    names(outs) <- c('checks_passed', 'priors', 'start_priors', 'sd_priors')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+}
+
+## helper for taking priors and calibartion to posteriors
+bayesian_checks <- function(prior, start_priors, sd_priors, params, count, prior_weight, weight, step_size, bounds = c(0, Inf), round_to = 1, round_to_digits = 1) {
+  checks_passed <- TRUE
+  
+  if (class(prior) == "matrix" && nrow(prior) == 1) {
+    priors <- round(rnorm(count, start_priors, sd_priors)/round_to, digits = round_to_digits)*round_to
+    priors <- as.data.frame(table(priors))
+    priors$priors <- as.numeric(as.character(priors$priors))
+    priors$prob <- round(priors$Freq/count, digits = 3)
+    priors$prob[priors$priors == bounds[2]] <- sum(priors$prob[priors$priors >= bounds[2]])
+    priors$prob[priors$priors == bounds[1]] <- sum(priors$prob[priors$priors <= bounds[1]])
+    priors <- priors[priors$prob > 0.000,]
+    priors <- priors[priors$priors <= bounds[2] & priors$priors >= bounds[1],]
+  } else if (class(prior) %in% c("matrix", "data.frame") && nrow(prior) > 1) {
+    priors <- prior[ , 1:2]
+    names(priors) <- c('priors', 'prob')
+    
+  }
+  
+  calibration_count <- length(params)
+  calibrated_rates <- as.data.frame(table(params))
+  calibrated_rates$params <- as.numeric(as.character(calibrated_rates$params))
+  calibrated_rates$prob <- round(calibrated_rates$Freq/calibration_count, digits = 3)
+  
+  min_rate <- min(min(priors$priors), min(calibrated_rates$params))
+  max_rate <- max(max(priors$priors), max(calibrated_rates$params))
+  
+  rates <- data.frame(rate = round(seq(min_rate, max_rate, step_size), digits = round_to_digits), 
+                                   prior_probability = rep(0, length(seq(min_rate, max_rate, step_size))),
+                                   calibrated_probability = rep(0, length(seq(min_rate, max_rate, step_size))),
+                                   posterior_probability = rep(0, length(seq(min_rate, max_rate, step_size))))
+  
+  for (i in 1:nrow(rates)) {
+    if (length(priors$prob[priors$priors == rates$rate[i]]) > 0) {
+      rates$prior_probability[i] <- priors$prob[priors$priors == rates$rate[i]]
+    }
+    if (length(calibrated_rates$prob[calibrated_rates$params == rates$rate[i]]) > 0) {
+      rates$calibrated_probability[i] <- calibrated_rates$prob[calibrated_rates$params == rates$rate[i]]
+    }
+  }
+  rates$posterior_probability <- round(rates$prior_probability*prior_weight + rates$calibrated_probability * weight, digits = 3)
+  posterior_rates <- rates[,c(1,4)]
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, rates, posterior_rates)
+    names(outs) <- c('checks_passed', 'rates', 'posterior_rates')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+}

From 5212128aa92c7938d7625926a0ccbd6a3bec0525 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Fri, 20 Dec 2019 10:20:54 -0500
Subject: [PATCH 17/25] added extra imports to clear up build notes for helpers

---
 NAMESPACE         | 3 +++
 R/calibrate.R     | 4 ++--
 R/checks.R        | 8 ++------
 R/pops.r          | 4 ++--
 R/pops_multirun.R | 4 ++--
 R/validate.R      | 4 ++--
 6 files changed, 13 insertions(+), 14 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index 2517cb25..b601b3c8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,16 +26,19 @@ importFrom(parallel,detectCores)
 importFrom(parallel,makeCluster)
 importFrom(parallel,stopCluster)
 importFrom(raster,as.matrix)
+importFrom(raster,calc)
 importFrom(raster,cellStats)
 importFrom(raster,cellsFromExtent)
 importFrom(raster,compareCRS)
 importFrom(raster,extension)
 importFrom(raster,extent)
+importFrom(raster,extract)
 importFrom(raster,getValues)
 importFrom(raster,ncol)
 importFrom(raster,nlayers)
 importFrom(raster,nrow)
 importFrom(raster,raster)
+importFrom(raster,rasterToPoints)
 importFrom(raster,reclassify)
 importFrom(raster,stack)
 importFrom(raster,values)
diff --git a/R/calibrate.R b/R/calibrate.R
index 0d62034d..13e75c10 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -19,7 +19,7 @@
 #' @param success_metric Choose which success metric to use for calibration. Choices are "quantity", "quantity and configuration", "residual error" and "odds ratio". Default is "quantity"
 #' @param 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). 
 #'
-#' @importFrom raster raster values as.matrix xres yres stack reclassify cellStats nlayers extent extension compareCRS getValues
+#' @importFrom raster raster values as.matrix xres yres stack reclassify cellStats nlayers extent extension compareCRS getValues calc extract rasterToPoints
 #' @importFrom stats runif rnorm
 #' @importFrom doParallel registerDoParallel
 #' @importFrom foreach  registerDoSEQ %dopar% %do% %:% foreach
@@ -223,7 +223,7 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
       return(treatments_check$failed_check)
     }
     
-    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
     if (treatment_check$checks_passed) {
       treatment_maps <- treatment_check$treatment_maps
     } else {
diff --git a/R/checks.R b/R/checks.R
index 5be3fd44..e1c6d9e4 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -78,7 +78,7 @@ secondary_raster_checks <- function(x, x2) {
   }
 }
 
-treatment_checks <- function(treatment_stack, treatments_file, pesticide_duration, treatment_dates) {
+treatment_checks <- function(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy) {
   checks_passed <- TRUE
   
   if (length(treatments_file) != length(treatment_dates)) {
@@ -342,9 +342,5 @@ bayesian_checks <- function(prior, start_priors, sd_priors, params, count, prior
     outs <- list(checks_passed, rates, posterior_rates)
     names(outs) <- c('checks_passed', 'rates', 'posterior_rates')
     return(outs)
-  } else {
-    outs <- list(checks_passed, failed_check)
-    names(outs) <- c('checks_passed', 'failed_check')
-    return(outs)
-  }
+  } 
 }
diff --git a/R/pops.r b/R/pops.r
index 65db2400..9675e073 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -43,7 +43,7 @@
 #' @param output_frequency sets when outputs occur either ('year', 'month' or 'time step')
 #' 
 #' @useDynLib PoPS, .registration = TRUE
-#' @importFrom raster raster values as.matrix xres yres stack extent
+#' @importFrom raster raster values as.matrix xres yres stack extent calc extract rasterToPoints
 #' @importFrom Rcpp sourceCpp evalCpp
 #' @importFrom  stats runif
 #' @importFrom lubridate interval time_length
@@ -221,7 +221,7 @@ pops <- function(infected_file, host_file, total_plants_file,
       return(treatments_check$failed_check)
     }
     
-    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
     if (treatment_check$checks_passed) {
       treatment_maps <- treatment_check$treatment_maps
     } else {
diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index 17766301..0e3e60ea 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -8,7 +8,7 @@
 #' @param num_iterations how many iterations do you want to run to allow the calibration to converge at least 10 
 #' @param 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
 #'
-#' @importFrom raster raster values as.matrix xres yres stack reclassify cellStats nlayers
+#' @importFrom raster raster values as.matrix xres yres stack reclassify cellStats nlayers calc extract rasterToPoints
 #' @importFrom stats runif rnorm median sd
 #' @importFrom doParallel registerDoParallel
 #' @importFrom foreach  registerDoSEQ %dopar%
@@ -186,7 +186,7 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
       return(treatments_check$failed_check)
     }
     
-    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
     if (treatment_check$checks_passed) {
       treatment_maps <- treatment_check$treatment_maps
     } else {
diff --git a/R/validate.R b/R/validate.R
index c6220648..dd35c105 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -12,7 +12,7 @@
 #' @param success_metric Choose which success metric to use for calibration. Choices are "quantity", "quantity and configuration", "residual error" and "odds ratio". Default is "quantity"
 #' @param 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).
 #'
-#' @importFrom raster raster values as.matrix xres yres stack reclassify cellStats nlayers
+#' @importFrom raster raster values as.matrix xres yres stack reclassify cellStats nlayers calc extract rasterToPoints
 #' @importFrom stats runif rnorm
 #' @importFrom doParallel registerDoParallel
 #' @importFrom foreach  registerDoSEQ %dopar%
@@ -205,7 +205,7 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
       return(treatments_check$failed_check)
     }
     
-    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates)
+    treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
     if (treatment_check$checks_passed) {
       treatment_maps <- treatment_check$treatment_maps
     } else {

From db1e824d4a698df72c4012b4e0b4664688b4b387 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Fri, 20 Dec 2019 12:23:07 -0500
Subject: [PATCH 18/25] added functions for propogating uncertainty into the
 model

---
 R/uncertainty_propogation.R | 68 +++++++++++++++++++++++++++++++++++++
 1 file changed, 68 insertions(+)
 create mode 100644 R/uncertainty_propogation.R

diff --git a/R/uncertainty_propogation.R b/R/uncertainty_propogation.R
new file mode 100644
index 00000000..6a1e5a23
--- /dev/null
+++ b/R/uncertainty_propogation.R
@@ -0,0 +1,68 @@
+## Uncertainty propogation checks and passing
+
+output_from_raster_mean_and_sd <- function(x) {
+  x[[1]] <- reclassify(x[[1]], matrix(c(-Inf, 0, 0), ncol = 3, byrow = TRUE))
+  x[[2]] <- reclassify(x[[2]], matrix(c(-Inf, 0, 0), ncol = 3, byrow = TRUE))
+  # x <- reclassify(x, c(-Inf, 0, 0))
+  fun <- function(x) {round(rnorm(1, mean = x[1], sd = x[2]), digits = 0)}
+  x2 <- suppressWarnings(calc(x, fun))
+  return(x2)
+}
+
+choose_variable_based_on_probability <- function(x, n = 1) {
+  names(x) <- c('var', 'prob')
+  for (i in 1:nrow(x)){
+    x$bin[i] <- sum(x$prob[1:i])
+  }
+    rn <- runif(n)
+    value <- rep(0,n)
+    for (j in 1:length(rn)){
+      for (i in 1:nrow(x)) {
+        if (i <= 1) {
+          if (rn[j] <= x$bin[i]) {
+            value[j] <- x$var[i]
+          }
+        } else {
+          if (rn[j] <= x$bin[i] && rn > x$bin[i-1]) {
+            value[j] <- x$var[i]
+          }
+        }
+      }
+    }
+  return(value)
+}
+
+uncertainty_check <- function(priors, round_to = 0, n = 1) {
+  checks_passed <- TRUE
+  
+  if (class(priors) == "numeric" && length(priors) == 2) {
+    priors <- matrix(priors, ncol = 2)
+  } 
+  
+  if (class(priors) %in% c("matrix", "data.frame") && ncol(priors) == 2) {
+    if (class(priors) == "matrix" && nrow(priors) == 1) {
+      value <- round(rnorm(n, priors[1], priors[2]), digits = round_to)
+    } else if (class(priors) %in% c("matrix", "data.frame") && nrow(priors) >= 1) {
+      if (class(priors) == "matrix") {
+        priors <- data.frame(priors)
+      }
+      value <- choose_variable_based_on_probability(priors, n)
+    }
+  } else if (class(priors) == "numeric" && length(priors) == 1) {
+    value <- rep(priors, n)
+  } else {
+    checks_passed <- FALSE
+    failed_check <- "Incorrect format for priors"
+  }
+  
+  if (checks_passed) {
+    outs <- list(checks_passed, value)
+    names(outs) <- c('checks_passed', 'value')
+    return(outs)
+  } else {
+    outs <- list(checks_passed, failed_check)
+    names(outs) <- c('checks_passed', 'failed_check')
+    return(outs)
+  }
+}
+

From c6a2f1f5230531f29cf6470b5f581ab10689885d Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Fri, 20 Dec 2019 13:53:06 -0500
Subject: [PATCH 19/25] edited uncertainty checks to ensure that the sum of
 probabilities adds up to 1

---
 R/uncertainty_propogation.R | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/R/uncertainty_propogation.R b/R/uncertainty_propogation.R
index 6a1e5a23..f3e97a38 100644
--- a/R/uncertainty_propogation.R
+++ b/R/uncertainty_propogation.R
@@ -13,6 +13,11 @@ choose_variable_based_on_probability <- function(x, n = 1) {
   names(x) <- c('var', 'prob')
   for (i in 1:nrow(x)){
     x$bin[i] <- sum(x$prob[1:i])
+    if (i == nrow(x)) {
+      if ((x$bin[i] < 1.00 && x$bin[i] > 0.985) || (x$bin[i] > 1.00 && x$bin[i] < 1.01)) {
+        x$bin[i] <- 1.00
+      }
+    }
   }
     rn <- runif(n)
     value <- rep(0,n)
@@ -46,6 +51,10 @@ uncertainty_check <- function(priors, round_to = 0, n = 1) {
       if (class(priors) == "matrix") {
         priors <- data.frame(priors)
       }
+      if (round(sum(priors[,2]), 3) < 0.985 || round(sum(priors[,2]), 3) > 1.01) {
+        checks_passed <- FALSE
+        failed_check <- "probabilities do not add up to 1"
+      }
       value <- choose_variable_based_on_probability(priors, n)
     }
   } else if (class(priors) == "numeric" && length(priors) == 1) {

From 818ab24ff8da10c55568c09d3a2639616ab350c1 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Fri, 20 Dec 2019 13:53:38 -0500
Subject: [PATCH 20/25] added uncertainty checks and propogation to single
 model wrapper

---
 R/pops.r | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/R/pops.r b/R/pops.r
index 9675e073..61f011dc 100644
--- a/R/pops.r
+++ b/R/pops.r
@@ -248,6 +248,34 @@ pops <- function(infected_file, host_file, total_plants_file,
   mortality <- mortality_tracker
   resistant <- mortality_tracker
   
+  reproductive_rate_check <- uncertainty_check(reproductive_rate, round_to = 1, n = 1)
+  if (reproductive_rate_check$checks_passed) {
+    reproductive_rate <- reproductive_rate_check$value
+  } else {
+    return(reproductive_rate_check$failed_check)
+  }
+  
+  natural_distance_scale_check <- uncertainty_check(natural_distance_scale, round_to = 0, n = 1)
+  if (natural_distance_scale_check$checks_passed) {
+    natural_distance_scale <- natural_distance_scale_check$value
+  } else {
+    return(natural_distance_scale_check$failed_check)
+  }
+    
+  anthropogenic_distance_scale_check <- uncertainty_check(anthropogenic_distance_scale, round_to = 0, n = 1)
+  if (anthropogenic_distance_scale_check$checks_passed) {
+    anthropogenic_distance_scale <- anthropogenic_distance_scale_check$value
+  } else {
+    return(anthropogenic_distance_scale_check$failed_check)
+  }
+  
+  percent_natural_dispersal_check <- uncertainty_check(percent_natural_dispersal, round_to = 3, n = 1)
+  if (percent_natural_dispersal_check$checks_passed) {
+    percent_natural_dispersal <- percent_natural_dispersal_check$value
+  } else {
+    return(percent_natural_dispersal_check$failed_check)
+  }
+  
   data <- PoPS::pops_model(random_seed = random_seed, 
                      use_lethal_temperature = use_lethal_temperature, 
                      lethal_temperature = lethal_temperature, lethal_temperature_month = lethal_temperature_month,

From 1617bf81679652eb43a225535fd3d9c79b151665 Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Fri, 20 Dec 2019 14:56:09 -0500
Subject: [PATCH 21/25] added uncertainty propogation into pops multirun and
 validate functions

---
 R/pops_multirun.R | 34 +++++++++++++++++++++++++++++++---
 R/validate.R      | 34 +++++++++++++++++++++++++++++++---
 2 files changed, 62 insertions(+), 6 deletions(-)

diff --git a/R/pops_multirun.R b/R/pops_multirun.R
index 0e3e60ea..a1d3c640 100644
--- a/R/pops_multirun.R
+++ b/R/pops_multirun.R
@@ -213,6 +213,34 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
   mortality <- mortality_tracker
   resistant <- mortality_tracker
   
+  reproductive_rate_check <- uncertainty_check(reproductive_rate, round_to = 1, n = num_iterations)
+  if (reproductive_rate_check$checks_passed) {
+    reproductive_rate <- reproductive_rate_check$value
+  } else {
+    return(reproductive_rate_check$failed_check)
+  }
+  
+  natural_distance_scale_check <- uncertainty_check(natural_distance_scale, round_to = 0, n = num_iterations)
+  if (natural_distance_scale_check$checks_passed) {
+    natural_distance_scale <- natural_distance_scale_check$value
+  } else {
+    return(natural_distance_scale_check$failed_check)
+  }
+  
+  anthropogenic_distance_scale_check <- uncertainty_check(anthropogenic_distance_scale, round_to = 0, n = num_iterations)
+  if (anthropogenic_distance_scale_check$checks_passed) {
+    anthropogenic_distance_scale <- anthropogenic_distance_scale_check$value
+  } else {
+    return(anthropogenic_distance_scale_check$failed_check)
+  }
+  
+  percent_natural_dispersal_check <- uncertainty_check(percent_natural_dispersal, round_to = 3, n = num_iterations)
+  if (percent_natural_dispersal_check$checks_passed) {
+    percent_natural_dispersal <- percent_natural_dispersal_check$value
+  } else {
+    return(percent_natural_dispersal_check$failed_check)
+  }
+  
   if (is.na(number_of_cores) || number_of_cores > parallel::detectCores()) {
     core_count <- parallel::detectCores() - 1
   } else {
@@ -243,14 +271,14 @@ pops_multirun <- function(infected_file, host_file, total_plants_file,
                        temperature = temperature,
                        weather_coefficient = weather_coefficient,
                        ew_res = ew_res, ns_res = ns_res, num_rows = num_rows, num_cols = num_cols,
-                       time_step = time_step, reproductive_rate = reproductive_rate,
+                       time_step = time_step, reproductive_rate = reproductive_rate[i],
                        mortality_rate = mortality_rate, mortality_time_lag = mortality_time_lag,
                        season_month_start = season_month_start, season_month_end = season_month_end,
                        start_date = start_date, end_date = end_date,
                        treatment_method = treatment_method,
                        natural_kernel_type = natural_kernel_type, anthropogenic_kernel_type = anthropogenic_kernel_type, 
-                       use_anthropogenic_kernel = use_anthropogenic_kernel, percent_natural_dispersal = percent_natural_dispersal,
-                       natural_distance_scale = natural_distance_scale, anthropogenic_distance_scale = anthropogenic_distance_scale, 
+                       use_anthropogenic_kernel = use_anthropogenic_kernel, percent_natural_dispersal = percent_natural_dispersal[i],
+                       natural_distance_scale = natural_distance_scale[i], anthropogenic_distance_scale = anthropogenic_distance_scale[i], 
                        natural_dir = natural_dir, natural_kappa = natural_kappa,
                        anthropogenic_dir = anthropogenic_dir, anthropogenic_kappa = anthropogenic_kappa,
                        output_frequency = output_frequency
diff --git a/R/validate.R b/R/validate.R
index dd35c105..a352e94d 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -232,6 +232,34 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   mortality <- mortality_tracker
   resistant <- mortality_tracker
   
+  reproductive_rate_check <- uncertainty_check(reproductive_rate, round_to = 1, n = num_iterations)
+  if (reproductive_rate_check$checks_passed) {
+    reproductive_rate <- reproductive_rate_check$value
+  } else {
+    return(reproductive_rate_check$failed_check)
+  }
+  
+  natural_distance_scale_check <- uncertainty_check(natural_distance_scale, round_to = 0, n = num_iterations)
+  if (natural_distance_scale_check$checks_passed) {
+    natural_distance_scale <- natural_distance_scale_check$value
+  } else {
+    return(natural_distance_scale_check$failed_check)
+  }
+  
+  anthropogenic_distance_scale_check <- uncertainty_check(anthropogenic_distance_scale, round_to = 0, n = num_iterations)
+  if (anthropogenic_distance_scale_check$checks_passed) {
+    anthropogenic_distance_scale <- anthropogenic_distance_scale_check$value
+  } else {
+    return(anthropogenic_distance_scale_check$failed_check)
+  }
+  
+  percent_natural_dispersal_check <- uncertainty_check(percent_natural_dispersal, round_to = 3, n = num_iterations)
+  if (percent_natural_dispersal_check$checks_passed) {
+    percent_natural_dispersal <- percent_natural_dispersal_check$value
+  } else {
+    return(percent_natural_dispersal_check$failed_check)
+  }
+  
   # reference <- raster(infected_years_file)
   ## Load observed data on occurence
   infection_years <- stack(infected_years_file)
@@ -276,14 +304,14 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
                              temperature = temperature,
                              weather_coefficient = weather_coefficient,
                              ew_res = ew_res, ns_res = ns_res, num_rows = num_rows, num_cols = num_cols,
-                             time_step = time_step, reproductive_rate = reproductive_rate,
+                             time_step = time_step, reproductive_rate = reproductive_rate[i],
                              mortality_rate = mortality_rate, mortality_time_lag = mortality_time_lag,
                              season_month_start = season_month_start, season_month_end = season_month_end,
                              start_date = start_date, end_date = end_date,
                              treatment_method = treatment_method,
                              natural_kernel_type = natural_kernel_type, anthropogenic_kernel_type = anthropogenic_kernel_type, 
-                             use_anthropogenic_kernel = use_anthropogenic_kernel, percent_natural_dispersal = percent_natural_dispersal,
-                             natural_distance_scale = natural_distance_scale, anthropogenic_distance_scale = anthropogenic_distance_scale, 
+                             use_anthropogenic_kernel = use_anthropogenic_kernel, percent_natural_dispersal = percent_natural_dispersal[i],
+                             natural_distance_scale = natural_distance_scale[i], anthropogenic_distance_scale = anthropogenic_distance_scale[i], 
                              natural_dir = natural_dir, natural_kappa = natural_kappa,
                              anthropogenic_dir = anthropogenic_dir, anthropogenic_kappa = anthropogenic_kappa,
                              output_frequency = output_frequency

From 5305db1a62f5cc7ef632e93144f11b9d2520d72f Mon Sep 17 00:00:00 2001
From: Chris Michael Jones <cjones1688@gmail.com>
Date: Fri, 20 Dec 2019 14:57:56 -0500
Subject: [PATCH 22/25] added helper functions for automated management

---
 R/helpers.R | 256 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 256 insertions(+)
 create mode 100644 R/helpers.R

diff --git a/R/helpers.R b/R/helpers.R
new file mode 100644
index 00000000..754606c0
--- /dev/null
+++ b/R/helpers.R
@@ -0,0 +1,256 @@
+## These functions are designed to reduce code complexity and the need to copy and past code across main functions
+
+getAllInfected <- function(rast, direction = 4) {
+  p <- rasterToPoints(rast, fun=function(x){x>0}, spatial = TRUE)
+  infections <- data.frame(extract(rast, p, cellnumbers=TRUE))
+  if (direction == 4) {
+    cellnumbersb <- data.frame(extract(rast, p, buf = xres(rast), cellnumbers=TRUE))
+  } else if (direction == 8) {
+    cellnumbersb <- data.frame(extract(rast, p, buf = 1.5 * xres(rast), cellnumbers=TRUE))
+  } else {
+    return("direction should be either of 4 or 8")
+  }
+  infections$i <- infections$cells %/% ncol(rast) + 1
+  infections$j <- infections$cells %% ncol(rast)
+  infections$group <- seq(1, nrow(infections),1)
+  cell_cols <- seq(1, ncol(cellnumbersb)-1,2)
+  for (i in cell_cols) {
+    infections$group[infections$cells %in% cellnumbersb[,i]] <- min(infections$group[infections$cells %in% cellnumbersb[,i]])
+    
+  }
+  groups <- unique(infections$group)
+  group_map <- seq(1, length(groups), 1)
+  for (m in 1:length(groups)) {
+    infections$group_size[infections$group == groups[m]] <- nrow(infections[infections$group == groups[m],])
+    infections$group[infections$group == groups[m]] <- group_map[m]
+    
+  }
+  names(infections) <- c('cell_number', 'detections', 'i', 'j', 'group', 'group_size')
+  return(infections)
+}
+
+getFoci <- function(rast) {
+  indexes <- getAllInfected(rast)
+  center <- data.frame(i = floor(mean(indexes$i)), j = floor(mean(indexes$j)))
+  return(center)
+}
+
+getInfectionBorder <- function(rast) {
+  s <- getAllInfected(rast)
+  min_max_col <- data.frame(row_number = seq(1, nrow(rast),1), min_j = rep(0,nrow(rast)), max_j = rep(0,nrow(rast)))
+  min_max_row <- data.frame(column_number = seq(1, ncol(rast), 1), min_i = rep(0,ncol(rast)), max_i = rep(0,ncol(rast)))
+  for (i in 1:nrow(rast)) {
+    if (length(s$j[s$i %in% i]) > 0) {
+      min_max_col$max_j[i] <- max(s$j[s$i %in% i])
+      min_max_col$min_j[i] <- min(s$j[s$i %in% i])
+    } else {
+      min_max_col$max_j[i] <- 0
+      min_max_col$min_j[i] <- 0
+    }
+  }
+  
+  for (j in 1:ncol(rast)) {
+    if (length(s$i[s$j %in% j]) > 0) {
+      min_max_row$max_i[j] <- max(s$i[s$j %in% j])
+      min_max_row$min_i[j] <- min(s$i[s$j %in% j])
+    } else {
+      min_max_row$max_i[j] <- 0
+      min_max_row$min_i[j] <- 0
+    }
+  }
+  
+  border1 <- data.frame(rows = rep(min_max_col$row_number,2), cols = c(min_max_col$min_j, min_max_col$max_j))
+  border2 <- data.frame(rows = c(min_max_row$min_i, min_max_row$max_i), cols = rep(min_max_row$column_number))
+  border <- rbind(border1, border2)
+  borders <- data.frame(table(border[,1:2]))
+  borders <- borders[borders$Freq > 0, ]
+  borders <- borders[borders$rows != 0, ]
+  borders <- borders[borders$cols != 0, ]
+  borders$rows <- as.numeric(as.character(borders$rows))
+  borders$cols <- as.numeric(as.character(borders$cols))
+  
+  names(borders) <- c('i','j', 'freq')
+  return(borders)
+}
+
+## options for get 
+getInfectionDistances <- function(rast, method = 'Foci', points = c()) {
+  infections <- getAllInfected(rast)
+  if (method == "Foci") {
+    points <- getFoci(rast)
+  } else if (method == "Border") {
+    points <- getInfectionBorder(rast)
+  } else if (method == "Points") {
+    points <- points
+  } else {
+    return("method must be one of 'Foci', 'Border', or'Points'")
+  }
+  # minimum_distance <- Inf
+  for (i in 1:nrow(infections)) {
+    minimum_distance <- Inf
+    for (l in 1:nrow(points)) {
+      new_distance <- sqrt((abs(infections$i[i] - points$i[l]))^2 + (abs(infections$j[i] - points$j[l]))^2)
+      if (new_distance < minimum_distance) {
+        minimum_distance <- new_distance
+      }
+    }
+    infections$distance[i] <- minimum_distance
+  }
+  
+  return(infections)
+}
+
+treatmentAuto <- function(rast, rast2, method = 'Foci', priority = 'group size', number_of_locations = 1, points = c(), treatment_efficacy = 1, buffer_cells = 1.5) {
+  ## get distances and groups and group size
+  if (method == "Points") {
+    points <- points
+  }
+  infections <- getInfectionDistances(rast = rast, method = method, points = points)
+  infections$host <- 0
+  for (p in 1:length(infections)) {
+    infections$host <- rast2[infections$i[p], infections$j[p]]
+  }
+  cells_treated <- 0
+  treatments <- data.frame(i = 0, j = 0, value = 0)
+  treatment <- rast
+  treatment[] <- 0
+  
+  if (priority == 'group size') {
+    infections <- infections[order(infections$group_size, infections$distance, decreasing = c(TRUE, FALSE)),]
+    groups_unique <- unique(infections$group)
+    for (group in groups_unique) {
+      managed_group <- infections[infections$group == group,]
+      group_size <- nrow(managed_group)
+      if (group_size + cells_treated >= number_of_locations) {
+        for (m in 1:nrow(managed_group)) {
+          i = managed_group$i[m]
+          j = managed_group$j[m]
+          if (treatment[i,j] < 1 & (rast[i, j] | rast2[i, j])) {
+            value <- min(1, treatment[i, j] + 1)
+            if (value > treatment[i, j]) {
+              cells_treated <- cells_treated + value - treatment[i, j]
+            } else {
+              cells_treated <- cells_treated + value
+            }
+            treatment[i, j] <- value
+            if (cells_treated >= number_of_locations) {break}
+          } 
+        }
+      } else if (group_size + cells_treated < number_of_locations) {
+        for (m in 1:nrow(managed_group)) {
+          i = managed_group$i[m]
+          j = managed_group$j[m]
+          if (treatment[i,j] < 1 & (rast[i, j] | rast2[i, j])) {
+            value <- min(1, treatment[i, j] + 1)
+            if (value > treatment[i, j]) {
+              cells_treated <- cells_treated + value - treatment[i, j]
+            } else {
+              cells_treated <- cells_treated + value
+            }
+            treatment[i, j] <- value
+            if (cells_treated >= number_of_locations) {break}
+          } 
+        }
+        for (m in 1:nrow(managed_group)) {
+          i = managed_group$i[m]
+          j = managed_group$j[m]
+          i_s <- seq(floor(i - buffer_cells), ceiling(i + buffer_cells), 1)
+          j_s <- seq(floor(j - buffer_cells), ceiling(j + buffer_cells), 1)
+          for (s in 1:length(i_s)) {
+            for(n in 1:length(j_s)) {
+              if (treatment[i_s[s], j_s[n]] < 1 & (rast[i_s[s], j_s[n]] | rast2[i_s[s], j_s[n]])) {
+                if (abs(i - i_s[s]) > buffer_cells | abs(j - j_s[n]) > buffer_cells) {
+                  value <- min(1, treatment[i_s[s], j_s[n]] + (buffer_cells - floor(buffer_cells)))
+                  if (value > treatment[i_s[s], j_s[n]]) {
+                    cells_treated <- cells_treated + value - treatment[i_s[s], j_s[n]]
+                  } else {
+                    cells_treated <- cells_treated + value
+                  }
+                  treatment[i_s[s], j_s[n]] <- value
+                } 
+              } 
+              if (cells_treated >= number_of_locations) {break}
+            }
+          }
+          if (cells_treated >= number_of_locations) {break}
+        }
+      }
+      if (cells_treated >= number_of_locations) {break}
+    }
+  } else if (priority == 'host') {
+    infections <- infections[order(infections$host, infections$distance, decreasing = c(TRUE, FALSE)),]
+    for (t in 1:nrow(infections)) {
+      i <- infections$i[t]
+      j <- infections$j[t]
+      if (treatment[i, j] < 1) {
+        cells_treated <- cells_treated + (1 - treatment[i, j])
+        treatment[i, j] <- 1
+      }
+      i_s <- seq(floor(i - buffer_cells), ceiling(i + buffer_cells), 1)
+      j_s <- seq(floor(j - buffer_cells), ceiling(j + buffer_cells), 1)
+      for (s in 1:length(i_s)) {
+        for(n in 1:length(j_s)) {
+          if (cells_treated >= number_of_locations) {break}
+          if (treatment[i_s[s], j_s[n]] < 1 & (rast[i_s[s], j_s[n]] | rast2[i_s[s], j_s[n]])) {
+            if (abs(i - i_s[s]) > buffer_cells | abs(j - j_s[n]) > buffer_cells) {
+              value <- min(1, treatment[i_s[s], j_s[n]] + (buffer_cells - floor(buffer_cells)))
+              if (value > treatment[i_s[s], j_s[n]]) {
+                cells_treated <- cells_treated + value - treatment[i_s[s], j_s[n]]
+              } else {
+                cells_treated <- cells_treated + value
+              }
+              treatment[i_s[s], j_s[n]] <- value
+            } else if (rast[i_s[s], j_s[n]] | rast2[i_s[s], j_s[n]]) {
+              if (treatment[i_s[s], j_s[n]] < 1) {
+                cells_treated <- cells_treated + (1 - treatment[i_s[s], j_s[n]])
+                treatment[i_s[s], j_s[n]] <- 1
+              }
+            }
+          } 
+        }
+      }
+      if (cells_treated >= number_of_locations) {break}
+    }
+  } else if (priority == 'infected') {
+    infections <- infections[order(infections$detections, infections$distance, decreasing = c(TRUE, FALSE)),]
+    for (t in 1:nrow(infections)) {
+      i <- infections$i[t]
+      j <- infections$j[t]
+      if (treatment[i, j] < 1) {
+        cells_treated <- cells_treated + (1 - treatment[i, j])
+        treatment[i, j] <- 1
+      }
+      i_s <- seq(floor(i - buffer_cells), ceiling(i + buffer_cells), 1)
+      j_s <- seq(floor(j - buffer_cells), ceiling(j + buffer_cells), 1)
+      for (s in 1:length(i_s)) {
+        for(n in 1:length(j_s)) {
+          if (cells_treated >= number_of_locations) {break}
+          if (treatment[i_s[s], j_s[n]] < 1 & (rast[i_s[s], j_s[n]] | rast2[i_s[s], j_s[n]])) {
+            if (abs(i - i_s[s]) > buffer_cells | abs(j - j_s[n]) > buffer_cells) {
+              value <- min(1, treatment[i_s[s], j_s[n]] + (buffer_cells - floor(buffer_cells)))
+              if (value > treatment[i_s[s], j_s[n]]) {
+                cells_treated <- cells_treated + value - treatment[i_s[s], j_s[n]]
+              } else {
+                cells_treated <- cells_treated + value
+              }
+              treatment[i_s[s], j_s[n]] <- value
+            } else if (rast[i_s[s], j_s[n]] | rast2[i_s[s], j_s[n]]) {
+              if (treatment[i_s[s], j_s[n]] < 1) {
+                cells_treated <- cells_treated + (1 - treatment[i_s[s], j_s[n]])
+                treatment[i_s[s], j_s[n]] <- 1
+              }
+            }
+          } 
+          
+        }
+      }
+      if (cells_treated >= number_of_locations) {break}
+    }
+  } else {
+    return('priority needs to be one of "group size", "host", or "infected"')
+  }
+  
+  treatment <- treatment * treatment_efficacy
+  
+  return(treatment)
+}
\ No newline at end of file

From 21b11c9a27e0170f76abcd37d32a81442c9f8182 Mon Sep 17 00:00:00 2001
From: ChrisJones687 <cjones1688@gmail.com>
Date: Fri, 3 Jan 2020 08:09:06 -0500
Subject: [PATCH 23/25] updated checks to not fail when early checks indicated
 a problem with inputs

---
 R/checks.R | 115 ++++++++++++++++++++++++++++-------------------------
 1 file changed, 60 insertions(+), 55 deletions(-)

diff --git a/R/checks.R b/R/checks.R
index e1c6d9e4..3825ae75 100644
--- a/R/checks.R
+++ b/R/checks.R
@@ -8,7 +8,7 @@ initial_raster_checks <- function(x) {
     failed_check <- "file does not exist" 
   }
   
-  if (!all((raster::extension(x) %in% c(".grd", ".tif", ".img")))) {
+  if (checks_passed && !all((raster::extension(x) %in% c(".grd", ".tif", ".img")))) {
     checks_passed <- FALSE
     failed_check <- "file is not one of '.grd', '.tif', '.img'"
   }
@@ -16,11 +16,11 @@ initial_raster_checks <- function(x) {
   if (checks_passed) {
     r<- raster::stack(x)
     r <- raster::reclassify(r, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
+    if (raster::nlayers(r) > 1) {
+      r <- output_from_raster_mean_and_sd(r)
+    }
   }
   
-  if (raster::nlayers(r) > 1) {
-    r <- output_from_raster_mean_and_sd(r)
-  }
   
   if (checks_passed) {
     outs <- list(checks_passed, r)
@@ -42,7 +42,7 @@ secondary_raster_checks <- function(x, x2) {
     failed_check <- "file does not exist" 
   }
   
-  if (!all((raster::extension(x) %in% c(".grd", ".tif", ".img")))) {
+  if (checks_passed && !all((raster::extension(x) %in% c(".grd", ".tif", ".img")))) {
     checks_passed <- FALSE
     failed_check <- "file is not one of '.grd', '.tif', '.img'"
   }
@@ -52,17 +52,17 @@ secondary_raster_checks <- function(x, x2) {
     r <- raster::reclassify(r, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   }
   
-  if (!(raster::extent(x2) == raster::extent(r))) {
+  if (checks_passed && !(raster::extent(x2) == raster::extent(r))) {
     checks_passed <- FALSE
     failed_check <- "Extents of input rasters do not match. Ensure that all of your input rasters have the same extent"
   }
   
-  if (!(raster::xres(x2) == raster::xres(r) && raster::yres(x2) == raster::yres(r))) {
+  if (checks_passed && !(raster::xres(x2) == raster::xres(r) && raster::yres(x2) == raster::yres(r))) {
     checks_passed <- FALSE
     failed_check <- "Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution"
   }
   
-  if (!raster::compareCRS(r,x2)) {
+  if (checks_passed && !raster::compareCRS(r,x2)) {
     checks_passed <- FALSE
     failed_check <- "Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs"
   }
@@ -81,29 +81,31 @@ secondary_raster_checks <- function(x, x2) {
 treatment_checks <- function(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy) {
   checks_passed <- TRUE
   
-  if (length(treatments_file) != length(treatment_dates)) {
+  if (checks_passed && length(treatments_file) != length(treatment_dates)) {
     checks_passed <- FALSE
     failed_check <- "Length of list for treatment dates and treatments_file must be equal"
   }
   
-  if (length(pesticide_duration) != length(treatment_dates)) {
+  if (checks_passed && length(pesticide_duration) != length(treatment_dates)) {
     checks_passed <- FALSE
     failed_check <- "Length of list for treatment dates and pesticide_duration must be equal"
   }
   
-  if (pesticide_duration[1] > 0) {
-    treatment_maps <- list(raster::as.matrix(treatment_stack[[1]] * pesticide_efficacy))
-  } else {
-    treatment_maps <- list(raster::as.matrix(treatment_stack[[1]]))
-  }
-  
-  if (raster::nlayers(treatment_stack) >= 2) {
-    for(i in 2:raster::nlayers(treatment_stack)) {
-      if (pesticide_duration[i] > 0) {
-        treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]] * pesticide_efficacy)
-      } else {
-        treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]])
-        
+  if (checks_passed) {
+    if (pesticide_duration[1] > 0) {
+      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]] * pesticide_efficacy))
+    } else {
+      treatment_maps <- list(raster::as.matrix(treatment_stack[[1]]))
+    }
+    
+    if (raster::nlayers(treatment_stack) >= 2) {
+      for(i in 2:raster::nlayers(treatment_stack)) {
+        if (pesticide_duration[i] > 0) {
+          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]] * pesticide_efficacy)
+        } else {
+          treatment_maps[[i]] <- raster::as.matrix(treatment_stack[[i]])
+          
+        }
       }
     }
   }
@@ -145,11 +147,13 @@ metric_checks <- function(success_metric) {
     configuration <- FALSE
   } else if (success_metric == "quantity and configuration") {
     configuration <- TRUE
-  } else if (success_metric == "odds_ratio") {
+  } else if (success_metric == "odds ratio") {
+    configuration <- FALSE
+  } else if (success_metric == "residual error") {
     configuration <- FALSE
   } else {
     checks_passed <- FALSE
-    failed_check <- "Success metric must be one of 'quantity', 'quantity and configuration', or 'odds_ratio'"
+    failed_check <- "Success metric must be one of 'quantity', 'quantity and configuration', 'residual error', or 'odds ratio'"
   }
   
   if (checks_passed) {
@@ -190,58 +194,59 @@ percent_checks <- function(percent_natural_dispersal) {
 time_checks <- function(end_date, start_date, time_step, output_frequency) {
   checks_passed <- TRUE
   
-  if (!(time_step %in% list("week", "month", "day"))) {
+  if (checks_passed && !(time_step %in% list("week", "month", "day"))) {
     checks_passed <- FALSE
     failed_check <- "Time step must be one of 'week', 'month' or 'day'"
   }
   
-  if (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d"))){
+  if (checks_passed && (class(end_date) != "character" || class(start_date) != "character" || class(as.Date(end_date, format="%Y-%m-%d")) != "Date" || class(as.Date(start_date, format="%Y-%m-%d")) != "Date" || is.na(as.Date(end_date, format="%Y-%m-%d")) || is.na(as.Date(start_date, format="%Y-%m-%d")))){
     checks_passed <- FALSE
     failed_check <- "End time and/or start time not of type numeric and/or in format YYYY-MM-DD"
   }
   
-  if (!(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
+  if (checks_passed && !(output_frequency %in% list("week", "month", "day", "year", "time_step"))) {
     checks_passed <- FALSE
-    failed_check <- "Time step must be one of 'week', 'month' or 'day'"
+    failed_check <- "Output frequency must be one of 'week', 'month' or 'day'"
   }
   
-  if (output_frequency == "day") {
+  if (checks_passed && output_frequency == "day") {
     if (time_step == "week" || time_step == "month") {
       checks_passed <- FALSE
       failed_check <- "Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible"
     }
   }
   
-  if (output_frequency == "week") {
+  if (checks_passed && output_frequency == "week") {
     if (time_step == "month") {
       checks_passed <- FALSE
       failed_check <- "Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible"
     }
   }
   
-  
-  duration <- lubridate::interval(start_date, end_date)
-  
-  if (time_step == "week") {
-    number_of_time_steps <- ceiling(lubridate::time_length(duration, "week"))
-  } else if (time_step == "month") {
-    number_of_time_steps <- ceiling(lubridate::time_length(duration, "month"))
-  } else if (time_step == "day") {
-    number_of_time_steps <- ceiling(lubridate::time_length(duration, "day"))
-  }
-  
-  number_of_years <- ceiling(lubridate::time_length(duration, "year"))
-  
-  if (output_frequency == "week") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "week"))
-  } else if (output_frequency == "month") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "month"))
-  } else if (output_frequency == "day") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "day"))
-  } else if (output_frequency == "year") {
-    number_of_outputs <- ceiling(lubridate::time_length(duration, "year"))
-  } else if (output_frequency == "time_step") {
-    number_of_outputs <- number_of_time_steps
+  if (checks_passed) {
+    duration <- lubridate::interval(start_date, end_date)
+    
+    if (time_step == "week") {
+      number_of_time_steps <- ceiling(lubridate::time_length(duration, "week"))
+    } else if (time_step == "month") {
+      number_of_time_steps <- ceiling(lubridate::time_length(duration, "month"))
+    } else if (time_step == "day") {
+      number_of_time_steps <- ceiling(lubridate::time_length(duration, "day"))
+    }
+    
+    number_of_years <- ceiling(lubridate::time_length(duration, "year"))
+    
+    if (output_frequency == "week") {
+      number_of_outputs <- ceiling(lubridate::time_length(duration, "week"))
+    } else if (output_frequency == "month") {
+      number_of_outputs <- ceiling(lubridate::time_length(duration, "month"))
+    } else if (output_frequency == "day") {
+      number_of_outputs <- ceiling(lubridate::time_length(duration, "day"))
+    } else if (output_frequency == "year") {
+      number_of_outputs <- ceiling(lubridate::time_length(duration, "year"))
+    } else if (output_frequency == "time_step") {
+      number_of_outputs <- number_of_time_steps
+    }
   }
   
   if (checks_passed) {
@@ -253,7 +258,7 @@ time_checks <- function(end_date, start_date, time_step, output_frequency) {
     names(outs) <- c('checks_passed', 'failed_check')
     return(outs)
   }
-  
+
 }
 
 ## check for making sure priors are in the proper format and output the mean where the mean serves as the starting point for the calibration

From 215c744b7cc81c02802c43f93e9e20fe8561bedd Mon Sep 17 00:00:00 2001
From: ChrisJones687 <cjones1688@gmail.com>
Date: Sat, 4 Jan 2020 07:23:51 -0500
Subject: [PATCH 24/25] added tests for new refactored checks

---
 tests/testthat/test-checks.R | 281 +++++++++++++++++++++++++++++++++++
 1 file changed, 281 insertions(+)
 create mode 100644 tests/testthat/test-checks.R

diff --git a/tests/testthat/test-checks.R b/tests/testthat/test-checks.R
new file mode 100644
index 00000000..662520ad
--- /dev/null
+++ b/tests/testthat/test-checks.R
@@ -0,0 +1,281 @@
+context("test-checks")
+
+test_that("Initial raster checks report correct errors and return a raster", {
+  infected_file <- ""
+  infected <- initial_raster_checks(infected_file)
+  expect_equal(infected$checks_passed, FALSE)
+  expect_equal(infected$failed_check, "file does not exist")
+  
+  infected_file <-  system.file("extdata", "simple2x2", "infected.csv", package = "PoPS")
+  infected <- initial_raster_checks(infected_file)
+  expect_equal(infected$checks_passed, FALSE)
+  expect_equal(infected$failed_check, "file is not one of '.grd', '.tif', '.img'")
+  
+  infected_file <- system.file("extdata", "simple2x2", "infected.tif", package = "PoPS")
+  infected <- initial_raster_checks(infected_file)
+  expect_equal(infected$checks_passed, TRUE)
+  expect_equal(infected$raster[], raster(infected_file)[])
+})
+
+test_that("Initial raster checks report correct errors and return a raster", {
+  infected_file <- system.file("extdata", "simple2x2", "infected.tif", package = "PoPS")
+  infected <- raster(infected_file)
+  host_file <- ""
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, FALSE)
+  expect_equal(host$failed_check, "file does not exist")
+  
+  host_file <-  system.file("extdata", "simple2x2", "infected.csv", package = "PoPS")
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, FALSE)
+  expect_equal(host$failed_check, "file is not one of '.grd', '.tif', '.img'")
+  
+  host_file <- system.file("extdata", "simple5x5", "total_plants.tif", package = "PoPS")
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, FALSE)
+  expect_equal(host$failed_check, "Extents of input rasters do not match. Ensure that all of your input rasters have the same extent")
+  
+  host_file <- system.file("extdata", "simple2x2", "total_plants_diff_res.tif", package = "PoPS")
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, FALSE)
+  expect_equal(host$failed_check, "Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
+  
+  host_file <- system.file("extdata", "simple2x2", "total_plants_diff_xres.tif", package = "PoPS")
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, FALSE)
+  expect_equal(host$failed_check, "Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
+  
+  host_file <- system.file("extdata", "simple2x2", "total_plants_diff_yres.tif", package = "PoPS")
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, FALSE)
+  expect_equal(host$failed_check, "Resolution of input rasters do not match. Ensure that all of your input rasters have the same resolution")
+  
+  host_file <- system.file("extdata", "simple2x2", "total_plants_with_crs.tif", package = "PoPS")
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, FALSE)
+  expect_equal(host$failed_check, "Coordinate reference system (crs) of input rasters do not match. Ensure that all of your input rasters have the same crs")
+  
+  host_file <-  system.file("extdata", "simple2x2", "total_plants.tif", package = "PoPS")
+  host <- secondary_raster_checks(host_file, infected)
+  expect_equal(host$checks_passed, TRUE)
+  expect_equal(host$raster[], raster(host_file)[])
+  
+})
+
+test_that("Treatment checks report correct errors and return a matrices of treatments if all checks pass", {
+  treatments_file <- system.file("extdata", "simple2x2", "treatments.tif", package = "PoPS")
+  treatment_stack <- stack(treatments_file)
+  pesticide_duration <- c(60)
+  treatment_dates <- c('2016-01-01', '2016-03-15')
+  pesticide_efficacy <- 1.0
+  treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
+  expect_equal(treatment_check$checks_passed, FALSE)
+  expect_equal(treatment_check$failed_check, "Length of list for treatment dates and treatments_file must be equal")
+  
+  treatments_file <- system.file("extdata", "simple2x2", "treatments.tif", package = "PoPS")
+  treatment_stack <- stack(treatments_file)
+  pesticide_duration <- c(60,45,60)
+  treatment_dates <- c('2016-01-01')
+  pesticide_efficacy <- 1.0
+  treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
+  expect_equal(treatment_check$checks_passed, FALSE)
+  expect_equal(treatment_check$failed_check, "Length of list for treatment dates and pesticide_duration must be equal")
+  
+  treatments_file <- system.file("extdata", "simple2x2", "treatments.tif", package = "PoPS")
+  treatment_stack <- stack(treatments_file)
+  pesticide_duration <- c(60)
+  treatment_dates <- c('2016-01-01')
+  pesticide_efficacy <- 1.0
+  treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
+  expect_equal(treatment_check$checks_passed, TRUE)
+  expect_equal(treatment_check$treatment_maps[[1]], as.matrix(treatment_stack[[1]]))
+  
+  treatments_file <- c(system.file("extdata", "simple2x2", "treatments.tif", package = "PoPS"), system.file("extdata", "simple2x2", "treatments.tif", package = "PoPS"))
+  treatment_stack <- stack(treatments_file)
+  pesticide_duration <- c(60, 45)
+  treatment_dates <- c('2016-01-01', '2016-10-02')
+  pesticide_efficacy <- 1.0
+  treatment_check <- treatment_checks(treatment_stack, treatments_file, pesticide_duration, treatment_dates, pesticide_efficacy)
+  expect_equal(treatment_check$checks_passed, TRUE)
+  expect_equal(treatment_check$treatment_maps[[1]], as.matrix(treatment_stack[[1]]))
+})
+
+
+test_that("Treatment method checks report correct errors and TRUE otherwise", {
+  treatment_method <- "Test"
+  methods_check <- treatment_metric_checks(treatment_method)
+  expect_equal(methods_check$checks_passed, FALSE)
+  expect_equal(methods_check$failed_check, "treatment method is not one of the valid treatment options")
+  treatment_method <- "ratio"
+  methods_check <- treatment_metric_checks(treatment_method)
+  expect_equal(methods_check$checks_passed, TRUE)
+  treatment_method <- "all infected"
+  methods_check <- treatment_metric_checks(treatment_method)
+  expect_equal(methods_check$checks_passed, TRUE)
+  
+})
+
+test_that("Success metric checks report correct errors and return configuration disagreement boolean otherwise", {
+  success_metric <- "Test"
+  metric_check <- metric_checks(success_metric)
+  expect_equal(metric_check$checks_passed, FALSE)
+  expect_equal(metric_check$failed_check, "Success metric must be one of 'quantity', 'quantity and configuration', 'residual error', or 'odds ratio'")
+  success_metric <- "quantity"
+  metric_check <- metric_checks(success_metric)
+  expect_equal(metric_check$checks_passed, TRUE)
+  expect_equal(metric_check$configuration, FALSE)
+  success_metric <- "quantity and configuration"
+  metric_check <- metric_checks(success_metric)
+  expect_equal(metric_check$checks_passed, TRUE)
+  expect_equal(metric_check$configuration, TRUE)
+  success_metric <- "odds ratio"
+  metric_check <- metric_checks(success_metric)
+  expect_equal(metric_check$checks_passed, TRUE)
+  expect_equal(metric_check$configuration, FALSE)
+  success_metric <- "residual error"
+  metric_check <- metric_checks(success_metric)
+  expect_equal(metric_check$checks_passed, TRUE)
+  expect_equal(metric_check$configuration, FALSE)
+})
+
+test_that("Percent checks report correct errors and return whether or not to use the anthropogenic distance kernel", {
+  percent_natural_dispersal <- 1.01
+  percent_check <- percent_checks(percent_natural_dispersal)
+  expect_equal(percent_check$checks_passed, FALSE)
+  expect_equal(percent_check$failed_check, "Percent natural dispersal must be between 0.0 and 1.0")
+  percent_natural_dispersal <- -0.2
+  percent_check <- percent_checks(percent_natural_dispersal)
+  expect_equal(percent_check$checks_passed, FALSE)
+  expect_equal(percent_check$failed_check, "Percent natural dispersal must be between 0.0 and 1.0")
+  percent_natural_dispersal <- 1.0
+  percent_check <- percent_checks(percent_natural_dispersal)
+  expect_equal(percent_check$checks_passed, TRUE)
+  expect_equal(percent_check$use_anthropogenic_kernel, FALSE)
+  percent_natural_dispersal <- 0.80
+  percent_check <- percent_checks(percent_natural_dispersal)
+  expect_equal(percent_check$checks_passed, TRUE)
+  expect_equal(percent_check$use_anthropogenic_kernel, TRUE)
+
+})
+
+test_that("Time checks report correct errors and return number of time steps, number of years, and number of outputs otherwise", {
+  end_date <- 2017-01-01
+  start_date <- 2016-01-01
+  time_step <- 'hour'
+  output_frequency <- 'week'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, FALSE)
+  expect_equal(time_check$failed_check, "Time step must be one of 'week', 'month' or 'day'")
+  
+  end_date <- 2017-01-01
+  start_date <- 2016-01-01
+  time_step <- 'week'
+  output_frequency <- 'week'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, FALSE)
+  expect_equal(time_check$failed_check, "End time and/or start time not of type numeric and/or in format YYYY-MM-DD")
+  
+  end_date <- '2017-01-01'
+  start_date <- '2016-01-01'
+  time_step <- 'month'
+  output_frequency <- 'hour'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, FALSE)
+  expect_equal(time_check$failed_check, "Output frequency must be one of 'week', 'month' or 'day'")
+  
+  end_date <- '2017-01-01'
+  start_date <- '2016-01-01'
+  time_step <- 'month'
+  output_frequency <- 'week'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, FALSE)
+  expect_equal(time_check$failed_check, "Output frequency is more frequent than time_step. The minimum output_frequency you can use is the time_step of your simulation. You can set the output_frequency to 'time_step' to default to most frequent output possible")
+
+  end_date <- '2017-01-01'
+  start_date <- '2016-01-01'
+  time_step <- 'week'
+  output_frequency <- 'week'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, TRUE)
+  expect_equal(time_check$number_of_time_steps, 53)
+  expect_equal(time_check$number_of_years, 1)
+  expect_equal(time_check$number_of_outputs, 53)
+  
+  end_date <- '2016-02-01'
+  start_date <- '2016-01-01'
+  time_step <- 'week'
+  output_frequency <- 'week'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, TRUE)
+  expect_equal(time_check$number_of_time_steps, 5)
+  expect_equal(time_check$number_of_years, 1)
+  expect_equal(time_check$number_of_outputs, 5)
+  
+  end_date <- '2017-01-01'
+  start_date <- '2016-01-01'
+  time_step <- 'week'
+  output_frequency <- 'month'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, TRUE)
+  expect_equal(time_check$number_of_time_steps, 53)
+  expect_equal(time_check$number_of_years, 1)
+  expect_equal(time_check$number_of_outputs, 12)
+  
+  end_date <- '2017-01-01'
+  start_date <- '2016-01-01'
+  time_step <- 'week'
+  output_frequency <- 'year'
+  time_check <- time_checks(end_date, start_date, time_step, output_frequency)
+  expect_equal(time_check$checks_passed, TRUE)
+  expect_equal(time_check$number_of_time_steps, 53)
+  expect_equal(time_check$number_of_years, 1)
+  expect_equal(time_check$number_of_outputs, 1)
+  
+})
+
+test_that("Prior checks report correct errors and return reformatted priors, and mean and sd of priors for starting location for MCMC", {
+  priors <- 2.1
+  prior_check <- prior_checks(priors)
+  expect_equal(prior_check$checks_passed, FALSE)
+  expect_equal(prior_check$failed_check, "Incorrect format for priors")
+  
+  priors <- data.frame(data = c(1.5,1.6,1.7,1.8,1.9,2.0,2.1), per = c(0.1,0.12,0.13,0.30,0.13,0.12,0.1))
+  prior_check <- prior_checks(priors)
+  names(priors) <- c('var', 'prob')
+  expect_equal(prior_check$checks_passed, TRUE)
+  expect_equal(prior_check$priors, priors)
+  expect_equal(prior_check$start_priors, 1.8)
+  expect_equal(prior_check$sd_priors, sd(priors$var))
+  
+  priors <- c(2.1, 0.2)
+  prior_check <- prior_checks(priors)
+  expect_equal(prior_check$checks_passed, TRUE)
+  expect_equal(prior_check$priors, matrix(priors, ncol = 2))
+  expect_equal(prior_check$start_priors, 2.1)
+  expect_equal(prior_check$sd_priors, 0.2)
+  
+})
+
+test_that("Bayesian checks report correct errors and return all rates as a data frame and posterior rates as a data frame", {
+  priors <- c(2.1, 0.2)
+  prior_check <- prior_checks(priors)
+  prior <- prior_check$priors
+  start_priors <- prior_check$start_priors
+  sd_priors <- prior_check$sd_priors
+  params <- data.frame(data = c(rep(1.5,10), rep(1.5,10), rep(1.6,15), rep(1.7,15), rep(1.8,30), rep(1.9,15), rep(2.0,15), rep(2.1,15), rep(2.2, 15)))
+  count <- 10000000
+  prior_weight <- 0.05
+  weight <- 1 - prior_weight
+  step_size <- 0.1
+  bounds <- c(0, Inf)
+  round_to <- 1
+  round_to_digits <- 1
+  bayesian_check <- bayesian_checks(prior, start_priors, sd_priors, params, count, prior_weight, weight, step_size, bounds = c(0, Inf), round_to = 1, round_to_digits = 1)
+  expect_equal(bayesian_check$checks_passed, TRUE)
+  expect_length(bayesian_check$rates, 4)
+  expect_length(bayesian_check$posterior_rates, 2)
+  expect_lte(mean(bayesian_check$posterior_rates[,1]), 2.1)
+  
+  
+  
+})

From 84c83ad6760c30bb5d87d4e5f3505650d06974a8 Mon Sep 17 00:00:00 2001
From: ChrisJones687 <cjones1688@gmail.com>
Date: Sat, 4 Jan 2020 09:45:36 -0500
Subject: [PATCH 25/25] setup quantity allocation disagreemt to handle
 reclassification so that residual error can be treated indepent of other
 success metrics

---
 R/calibrate.R                        |  7 -----
 R/quantity_allocation_disagreement.R | 41 +++++++---------------------
 R/validate.R                         | 16 ++---------
 3 files changed, 12 insertions(+), 52 deletions(-)

diff --git a/R/calibrate.R b/R/calibrate.R
index 13e75c10..4b6c0c46 100644
--- a/R/calibrate.R
+++ b/R/calibrate.R
@@ -264,11 +264,6 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   if (length(total_infections) > number_of_years){
     total_infections <- total_infections[1:number_of_years]
   }
-  ## set up reclassification matrix for binary reclassification
-  rcl <- c(1, Inf, 1, 0, 0.99, NA)
-  rclmat <- matrix(rcl, ncol=3, byrow=TRUE)
-  ## reclassify to binary values
-  infection_years <- reclassify(infection_years, rclmat)
   ## Get rid of NA values to make comparisons
   infection_years <- raster::reclassify(infection_years, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   num_layers_infected_years <- raster::nlayers(infection_years)
@@ -319,7 +314,6 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
   comp_year <- raster(infected_file)
   all_disagreement <- foreach(q = 1:length(data$infected), .combine = rbind, .packages =c("raster", "PoPS"), .final = colSums) %do% {
     comp_year[] <- data$infected[[q]]
-    comp_year <- reclassify(comp_year, rclmat)
     to.all_disagreement <- quantity_allocation_disagreement(infection_years[[q]], comp_year, configuration, mask)
   }
   
@@ -364,7 +358,6 @@ calibrate <- function(infected_years_file, num_iterations,  number_of_cores = NA
       # set up comparison
       all_disagreement <- foreach(q = 1:length(data$infected), .combine = rbind, .packages =c("raster", "PoPS"), .final = colSums) %dopar% {
         comp_year[] <- data$infected[[q]]
-        comp_year <- reclassify(comp_year, rclmat)
         to.all_disagreement <- quantity_allocation_disagreement(infection_years[[q]], comp_year, configuration, mask)
       }
       
diff --git a/R/quantity_allocation_disagreement.R b/R/quantity_allocation_disagreement.R
index 87609797..203c7564 100644
--- a/R/quantity_allocation_disagreement.R
+++ b/R/quantity_allocation_disagreement.R
@@ -34,7 +34,15 @@ quantity_allocation_disagreement <- function(reference, comparison, configuratio
   }
   ## test that the comparison raster is the same extent, resolution, and crs as the reference to ensure that they can be compared accurately
   raster::compareRaster(reference, comparison)
-
+  # save initial reference and comparison to use for residual error calculation and then reclassify the original reference and comparison to binary values for other
+  # classifications as we are only concerned with locational accuracy of infection not exact population predictions (residual error is a comparison of exact population numbers)
+  comp <- comparison
+  ref <- reference
+  rcl <- c(1, Inf, 1, 0, 0.99, 0)
+  rclmat <- matrix(rcl, ncol=3, byrow=TRUE)
+  reference <- reclassify(reference, rclmat)
+  comparison <- reclassify(comparison, rclmat)
+  
   if (configuration == TRUE) {
     # calculate number of infected patches
     NP_ref <- landscapemetrics::lsm_c_np(reference, directions = 8)$value[2]
@@ -113,35 +121,6 @@ quantity_allocation_disagreement <- function(reference, comparison, configuratio
   false_negative <- sum(comparison[reference == 1] == 0)
   true_negative <- sum(comparison[reference == 0] == 0)
   
-  # ## calculate probabilities of each class based on Death to Kappa (Pontius et al. 2011)
-  # if (negatives_in_comparison == 0 && total_in_reference == 0) {
-  #   probability_00 <- (true_negative / 1) * (negatives_in_reference / 1)
-  #   probability_01 <- (false_negative / 1) * (negatives_in_reference / 1)
-  # } else if (total_in_reference == 0) {
-  #   probability_00 <- (true_negative / negatives_in_comparison) * (negatives_in_reference / 1)
-  #   probability_01 <- (false_negative / negatives_in_comparison) * (negatives_in_reference / 1)
-  # } else if (negatives_in_comparison == 0) {
-  #   probability_00 <- (true_negative / 1) * (negatives_in_reference / total_in_reference)
-  #   probability_01 <- (false_negative / 1) * (negatives_in_reference / total_in_reference)
-  # } else {
-  #   probability_00 <- (true_negative / negatives_in_comparison) * (negatives_in_reference / total_in_reference)
-  #   probability_01 <- (false_negative / negatives_in_comparison) * (negatives_in_reference / total_in_reference)
-  # }
-  # 
-  # if (positives_in_comparison == 0 && total_in_reference == 0) {
-  #   probability_10 <- (false_positive / 1) * (positives_in_reference / 1)
-  #   probability_11 <- (true_positive / 1) * (positives_in_reference / 1)
-  # } else if (total_in_reference == 0) {
-  #   probability_10 <- (false_positive / positives_in_comparison) * (positives_in_reference / 1)
-  #   probability_11 <- (true_positive / positives_in_comparison) * (positives_in_reference / 1)
-  # } else if (positives_in_comparison == 0) {
-  #   probability_10 <- (false_positive / 1) * (positives_in_reference / total_in_reference)
-  #   probability_11 <- (true_positive / 1) * (positives_in_reference / total_in_reference)
-  # } else {
-  #   probability_10 <- (false_positive / positives_in_comparison) * (positives_in_reference / total_in_reference)
-  #   probability_11 <- (true_positive / positives_in_comparison) * (positives_in_reference / total_in_reference)
-  # }
-  
   ## calculate quantity and allocation disagreements for infected/infested from probabilities based on Death to Kappa (Pontius et al. 2011)
   # quantity_disagreement <- abs((probability_11 + probability_01) - (probability_11 + probability_10))
   # allocation_disagreement <- min((probability_11 + probability_10) - probability_11, (probability_11 + probability_01) - probability_11)
@@ -169,7 +148,7 @@ quantity_allocation_disagreement <- function(reference, comparison, configuratio
   output$total_disagreement <- total_disagreement
   output$configuration_disagreement <- configuration_disagreement
   output$odds_ratio <- odds_ratio
-  output$residual_error <- cellStats(abs(reference - comparison), 'sum')
+  output$residual_error <- cellStats(abs(ref - comp), 'sum')
   
   return(output)
 }
\ No newline at end of file
diff --git a/R/validate.R b/R/validate.R
index a352e94d..9c98445c 100644
--- a/R/validate.R
+++ b/R/validate.R
@@ -260,14 +260,9 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     return(percent_natural_dispersal_check$failed_check)
   }
   
-  # reference <- raster(infected_years_file)
   ## Load observed data on occurence
   infection_years <- stack(infected_years_file)
-  ## set up reclassification matrix for binary reclassification
-  rcl <- c(1, Inf, 1, 0, 0.99, NA)
-  rclmat <- matrix(rcl, ncol=3, byrow=TRUE)
-  ## reclassify to binary values
-  infection_years <- reclassify(infection_years, rclmat)
+
   ## Get rid of NA values to make comparisons
   infection_years <- raster::reclassify(infection_years, matrix(c(NA,0), ncol = 2, byrow = TRUE), right = NA)
   
@@ -285,7 +280,7 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
   registerDoParallel(cl)
   
 
-  qa <- foreach::foreach (icount(num_iterations), .combine = rbind, .packages = c("raster", "PoPS", "foreach")) %dopar% {
+  qa <- foreach::foreach (i = 1:num_iterations, .combine = rbind, .packages = c("raster", "PoPS", "foreach")) %dopar% {
     random_seed <- round(stats::runif(1, 1, 1000000))
     data <- PoPS::pops_model(random_seed = random_seed, 
                              use_lethal_temperature = use_lethal_temperature, 
@@ -320,16 +315,9 @@ validate <- function(infected_years_file, num_iterations, number_of_cores = NA,
     comp_year <- raster(infected_file)
     all_disagreement <- foreach(q = 1:length(data$infected), .combine = rbind, .packages =c("raster", "PoPS", "foreach"), .final = colSums) %dopar% {
       comp_year[] <- data$infected[[q]]
-      comp_year <- reclassify(comp_year, rclmat)
       to.all_disagreement <- quantity_allocation_disagreement(infection_years[[q]], comp_year, configuration, mask)
     }
     
-    # comparison <- reference
-    # raster::values(comparison) <- data$infected[[1]]
-    # 
-    # comparison <- raster::reclassify(comparison, rclmat)
-    
-    # to.qa <- PoPS::quantity_allocation_disagreement(reference, comparison)
     to.qa <- data.frame(t(all_disagreement))
   }