diff --git a/src/prediction/abundance_prediction_bootstrap_per_boot.R b/src/prediction/abundance_prediction_bootstrap_per_boot.R index f66f9667cc80d62b763bd0f2d82211e19b0651f3..718e801a71211f2a809350cb66baa8bf40c2f1d2 100644 --- a/src/prediction/abundance_prediction_bootstrap_per_boot.R +++ b/src/prediction/abundance_prediction_bootstrap_per_boot.R @@ -134,31 +134,35 @@ all_boots <- readRDS(all_boots_path) %>% # HERE it already starts to be year specific years_to_predict <- 1999:2015 +#preds in model +predictor_names <- c("year", "temp_mean", "rain_var", "rain_dry", "dom_T_OC", + "peatswamp", "lowland_forest", + "lower_montane_forest", "deforestation_hansen", + "human_pop_dens", "ou_killing_prediction", + "perc_muslim" ) + +quadratic_terms_names <- c("rain_dry") + -pred_per_boot <- foreach(year_to_predict = years_to_predict, .combine = rbind) %do% { +pred_per_boot <- foreach(year_to_predict = years_to_predict, .combine = cbind) %dopar% { suppressPackageStartupMessages(library(stringr)) +suppressMessages(library(doParallel)) + predictors_grid_path <- path.to.current(indir, paste0("predictors_grid_scaled_", year_to_predict), "rds") -if(is_verbose){print(paste("this is predictors_grid_path", predictors_grid_path))} -predictors_grid <- readRDS(predictors_grid_path) - # predictors used in model -predictor_names <- c("year", "temp_mean", "rain_var", "rain_dry", "dom_T_OC", - "peatswamp", "lowland_forest", - "lower_montane_forest", "deforestation_hansen", - "human_pop_dens", "ou_killing_prediction", - "perc_muslim" ) +predictors_grid <- readRDS(predictors_grid_path) -quadratic_terms_names <- c("rain_dry") intercept <- rep(1, nrow(predictors_grid)) + predictor_estimates <- cbind( intercept, predictors_grid[ , predictor_names], - predictors_grid[ ,quadratic_terms_names] * predictors_grid[ ,quadratic_terms_names]) - + predictors_grid[ ,quadratic_terms_names] * predictors_grid[ ,quadratic_terms_names]) + names(predictor_estimates) <- c("intercept", predictor_names, paste0("I(", quadratic_terms_names, "^2)")) @@ -174,49 +178,30 @@ names(predictor_estimates) <- c("intercept", predictor_names, if(is_verbose){print(paste("1. start pred_per_cell", Sys.time()))} - - #pred_per_cell <- foreach(cell = 1:nrow(predictor_estimates), .combine = c) %dopar% { - pred_per_cell <- foreach(cell = 1:10, .combine = c) %dopar% { + pred_per_cell <- foreach(cell = 1:nrow(predictor_estimates), .combine = c) %do% { t_predictor_estimates <- t( predictor_estimates[cell, ]) pred_estimates_wcoeffs <- data.frame(mapply(`*`, all_boots[nr_boot, ], t_predictor_estimates, SIMPLIFY = F)) pred_estimate_cell <- apply(pred_estimates_wcoeffs, 1, sum) pred_estimate_cell <- exp(pred_estimate_cell) return(pred_estimate_cell) } - return(pred_per_boot) + return(pred_per_cell) } - if(is_verbose){print(paste(Sys.time(), "2. finished dopar loop"))} pred_per_boot <- as.data.frame(pred_per_boot) -sum_years <- as.data.frame(cbind(years_to_predict, pred_per_boot$V1)) -names(sum_years) <- c("year", "sum") -min_years <- as.data.frame(cbind(years_to_predict, pred_per_boot$V2)) -names(min_years) <- c("year", "min") -max_years <- as.data.frame(cbind(years_to_predict, pred_per_boot$V3)) -names(max_years) <- c("year", "max") +names(pred_per_boot) <- paste0("year_", c(1999:2015)) -saveRDS(sum_years, +saveRDS(pred_per_boot, file = file.path(outdir, - paste0("bootstrapped_abundance_sum_boot_", + paste0("bootstrapped_pred_per_pixel_", nr_boot, "_", Sys.Date(), ".rds"))) -saveRDS(min_years, - file = file.path(outdir, - paste0("bootstrapped_abundance_min_boot_", - nr_boot, "_", - Sys.Date(), ".rds"))) - -saveRDS(max_years, - file = file.path(outdir, - paste0("bootstrapped_abundance_max_boot_", - nr_boot, "_", - Sys.Date(), ".rds"))) @@ -224,3 +209,7 @@ print(paste("Finish model_prediction script for boot", nr_boot, Sys.time())) +if (options$worker > 1) { + stopCluster(cl) + } +