run_offset_simulation_routines <- function(simulation_inputs, current_simulation_params, global_params, parcels,
decline_rates_initial, dev_probability_list, offset_probability_list,
scenario_ind, realisation_ind){
# run simulation with identical realisation instantiation
# list used to govern feature_layers rate changes
current_data_dir = write_folder(paste0(global_params$output_folder,
'scenario_', formatC(scenario_ind, width = 3, format = "d", flag = "0"),
'/realisation_', formatC(realisation_ind, width = 3, format = "d", flag = "0"), '/'))
flog.info('current data dir is %s', current_data_dir)
# run the model and return outputs
simulation_outputs <- run_simulation(simulation_inputs,
global_params,
current_simulation_params,
parcels,
decline_rates_initial,
dev_probability_list,
offset_probability_list,
current_data_dir)
# save raw simulation data
if (global_params$save_simulation_outputs == TRUE){
saveRDS(simulation_outputs, paste0(current_data_dir, 'realisation_',
formatC(realisation_ind, width = 3, format = "d", flag = "0"),
'_outputs.rds'))
}
for (feature_ind in seq(current_simulation_params$feature_num)){
# select current layer index
current_feature = current_simulation_params$features_to_use_in_simulation[feature_ind]
# read current feature layer files over time series and arrange into data stack
current_data_stack = form_data_stack(current_data_dir,
feature_string = formatC(current_feature, width = 3, format = "d", flag = "0"),
land_parcels = parcels$land_parcels,
current_simulation_params$time_steps)
# run series of routines used to calculate gains and losses at multiple scales over current feature layer
current_collated_realisation = run_collate_routines(simulation_outputs,
current_data_stack,
decline_rates_initial,
simulation_inputs$current_feature_layers,
current_data_dir,
current_simulation_params,
realisation_ind,
feature_ind)
file_prefix = paste0(global_params$collated_folder,
'collated_scenario_', formatC(scenario_ind, width = 3, format = "d", flag = "0"),
'_realisation_', formatC(realisation_ind, width = 3, format = "d", flag = "0"))
saveRDS(current_collated_realisation, paste0(file_prefix, '_feature_',
formatC(current_feature, width = 3, format = "d", flag = "0"), '.rds'))
if ((realisation_ind == 1) && (feature_ind == 1)){
mov_folder = paste0(global_params$collated_folder, '/mov_',
formatC(scenario_ind, width = 3, format = "d", flag = "0"), '/')
if( (global_params$write_movie == TRUE) || (global_params$write_offset_layer == TRUE) ){
if(!(file.exists(mov_folder))){
dir.create(mov_folder)
}
}
if ( global_params$write_movie == TRUE){
write_frames(current_data_stack, filetype = 'png', mov_folder, parcels, global_params, current_simulation_params,
offset_site_indexes = unlist(simulation_outputs$offsets$site_indexes),
offset_yrs = unlist(simulation_outputs$offsets$offset_yrs))
}
}
}
# delete current temporary files and folder
if (global_params$save_simulation_outputs == FALSE){
unlink(current_data_dir, recursive = TRUE)
}
}
# main engine for code - returns all simulation outputs including developments, offsets etc.
run_simulation <- function(simulation_outputs, global_params, current_simulation_params, parcels,
decline_rates_initial, dev_probability_list, offset_probability_list, current_data_dir){
#run through main time loop
for (yr in seq_len(current_simulation_params$time_steps)){
flog.info('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~')
flog.info('t = %s', yr)
# if running multiple regions with distinct policies cycle through each region
#when running in offset banking select out current set of sites to add to bank
if (current_simulation_params$use_offset_bank == TRUE){
simulation_outputs <- perform_banking_routine(simulation_outputs, current_simulation_params, yr)
}
# determine current set of available offset sites and calculate gains structure as detailed in current policy params
simulation_outputs$offset_pool_object <- prepare_offset_pool(simulation_outputs,
current_simulation_params,
parcels,
decline_rates_initial,
yr)
# cycle through number of developments and associated offsets
for (current_dev_index in seq_len(current_simulation_params$intervention_vec[yr])){
if (current_simulation_params$allow_developments_from_credit == TRUE){
###### TODO REINSTATE OFFSETBANK ROUTINES ####
# if (current_simulation_params$use_offset_bank == TRUE){
# simulation_outputs$current_credit = assess_banking_credit(simulation_outputs, current_simulation_params)
# }
# attempt to develop from current available credit
flog.info(cat('current credit ', paste(simulation_outputs$current_credit), '\n'))
credit_match_object = develop_from_credit(simulation_outputs$current_feature_layers,
simulation_outputs$current_credit,
dev_probability_list,
current_simulation_params,
intervention_vec = current_simulation_params$intervention_vec,
dev_indexes_to_use = simulation_outputs$index_object$indexes_to_use$devs,
decline_rates_initial,
parcels$land_parcels,
yr,
current_simulation_params$offset_time_horizon)
# if development was permitted add current site to current group of developments, set all feature_layers of development site to zero
if (credit_match_object$match_flag == TRUE){
simulation_outputs$current_credit = credit_match_object$current_credit
flog.info(cat('developed site with value', paste(round(unlist(credit_match_object$development_object$parcel_vals_used), 1)),
'from credit, remaining =', paste(round(unlist(credit_match_object$current_credit), 1)), '\n'))
simulation_outputs <- perform_clearing_routine(simulation_outputs,
simulation_outputs$index_object,
simulation_outputs$decline_rates,
current_dev_object = credit_match_object$development_object,
clearing_type = 'develop_from_credit',
current_simulation_params)
}
} else {
credit_match_object = false_match()
}
#if insufficient credits accumulated to allow development attempt development with offset match.
if ( (credit_match_object$match_flag == FALSE && current_simulation_params$use_parcel_sets == TRUE)){
#perform the matching routine - i.e. find a matching development/offset set
# this is the trickiest set of routines and may need series of robustness checks.
match_object <- match_sites(offset_pool_object = simulation_outputs$offset_pool_object,
simulation_outputs$current_credit,
dev_probability_list,
offset_probability_list,
current_simulation_params,
intervention_vec = current_simulation_params$intervention_vec,
indexes_to_use = simulation_outputs$index_object$indexes_to_use$devs,
simulation_outputs$current_feature_layers,
decline_rates = simulation_outputs$decline_rates,
parcels$land_parcels,
yr,
current_simulation_params$offset_time_horizon)
# if a match was found record current development and associated offsets and update site parameters.
if (match_object$match_flag == TRUE){
flog.info(cat('developed site', paste(match_object$development_object$site_indexes),
'with value', paste(round(unlist(match_object$development_object$parcel_vals_used), 1)),
'and offset with sites', paste(match_object$offset_object$site_indexes),
'of value', paste(round(unlist(match_object$offset_object$parcel_vals_used), 1)), '\n'))
#update available credit
simulation_outputs$current_credit = match_object$current_credit
# remove selected offset sites from available pool, save offset site characteristics, update decline rates to implement offset.
simulation_outputs <- perform_offset_routine(simulation_outputs,
index_object = simulation_outputs$index_object,
decline_rates = simulation_outputs$decline_rates,
current_offset_object = match_object$offset_object,
current_offset_indexes = match_object$offset_object$site_indexes,
match_object,
current_simulation_params)
# remove selected development sites from available pool, save development site characteristics,
# update decline rates to implement development loss.
simulation_outputs <- perform_clearing_routine(simulation_outputs,
simulation_outputs$index_object,
simulation_outputs$decline_rates,
current_dev_object = match_object$development_object,
clearing_type = 'development',
current_simulation_params)
}
}
}
#select sites for clearing
unregulated_loss_object <- perform_unregulated_loss(simulation_outputs$current_feature_layers,
simulation_outputs$index_object,
yr,
current_simulation_params,
decline_rates_initial,
current_simulation_params$offset_time_horizon)
# remove selected sites from available pool, save site characteristics, update decline rates to implement loss.
if (!is.null(unregulated_loss_object)){
simulation_outputs <- perform_clearing_routine(simulation_outputs,
index_object = simulation_outputs$index_object,
decline_rates = simulation_outputs$decline_rates,
current_dev_object = unregulated_loss_object,
clearing_type = 'unregulated',
current_simulation_params)
}
if ( (length(unlist(simulation_outputs$offsets_object$site_indexes)) > 0) & (current_simulation_params$offset_action_type == 'restore')){
# assess whether offsets have achieved their gains as determined in initial offset gains calculation
assessed_parcel_sets_object <- assess_parcel_sets(simulation_outputs$current_feature_layers,
simulation_outputs$offsets_object,
simulation_outputs$index_object$site_indexes_used$offsets,
current_simulation_params,
decline_rates_initial,
current_simulation_params$offset_time_horizon,
yr)
# update feature_layers change parameters
simulation_outputs$decline_rates <- update_decline_rates(simulation_outputs$decline_rates,
restoration_rate = current_simulation_params$restoration_rate,
restoration_rate_std = current_simulation_params$restoration_rate_std,
features_to_use_in_offset_intervention = current_simulation_params$features_to_use_in_offset_intervention,
feature_num = current_simulation_params$feature_num,
decline_rate_type = 'offset',
action_type = 'maintain',
site_indexes = assessed_parcel_sets_object$site_success_inds)
}
# implement ecological loss for developed sites (those with decline_rates = 0)
simulation_outputs$current_feature_layers <- destroy_development_sites(simulation_outputs$current_feature_layers,
simulation_outputs$decline_rates,
current_simulation_params$feature_num)
for (feature_ind in seq(current_simulation_params$feature_num)){
feature_layers_to_save = lapply(seq_along(simulation_outputs$current_feature_layers), function(i) simulation_outputs$current_feature_layers[[i]][[feature_ind]])
saveRDS(feature_layers_to_save, paste0(current_data_dir,
'feature_', formatC(current_simulation_params$features_to_use_in_simulation[feature_ind], width = 3, format = "d", flag = "0"),
'_yr_', formatC(yr, width = 3, format = "d", flag = "0"), '.rds'))
}
# update feature_layers for all sites (including those just developed/offset)
simulation_outputs$current_feature_layers = project_sites(simulation_outputs$current_feature_layers,
simulation_outputs$decline_rates,
current_time_horizons = rep(list(1), length(simulation_outputs$current_feature_layers)),
current_simulation_params,
features_to_project = seq_along(current_simulation_params$features_to_use_in_simulation),
time_fill = FALSE)
}
return(simulation_outputs)
}
# build nested list of feature_layers change parameters of used to determine how the feature_layers evolves in the absence of the
# development and offset intervention. List length is determined by how many sites and list depth is determined by
# how many features are currently in use
#' @export
simulate_decline_rates <- function(parcel_num, sample_decline_rate, mean_decline_rates, decline_rate_std){
feature_num = length(mean_decline_rates)
if (sample_decline_rate == TRUE){
# sample change rate from normal distribution
decline_rates = lapply(seq(parcel_num), function(i) lapply(seq(feature_num),
function(j) rnorm(1, mean_decline_rates[[j]], decline_rate_std[[j]])))
} else {
# copy same rate to all sites
decline_rates = lapply(seq(parcel_num), function(i) lapply(seq(feature_num),
function(j) mean_decline_rates[[j]]))
}
return(decline_rates)
}
# for multuiple regions with multiple policies select current policy
select_current_policy <- function(current_simulation_params_set, region_num){
if (region_num > 1){
current_simulation_params = current_simulation_params_set
} else {
current_simulation_params = current_simulation_params_set
}
return(current_simulation_params)
}
# used after simulation has completed. Takes time series layer files for each feature and stacks
form_data_stack <- function(current_data_dir, feature_string, land_parcels, time_steps){
current_filenames <- list.files(path = current_data_dir,
pattern = paste0('feature_', feature_string), all.files = FALSE,
full.names = FALSE, recursive = FALSE, ignore.case = FALSE,
include.dirs = FALSE, no.. = FALSE)
data_stack = lapply(seq_along(land_parcels), function(i) array(0, c(time_steps, length(land_parcels[[i]]))))
for (yr in seq(time_steps)){
current_feature_layers = readRDS(paste0(current_data_dir, current_filenames[yr]))
data_stack <- lapply(seq_along(data_stack), function(i) stack_current_yr(data_stack[[i]], current_feature_layers[[i]], yr))
}
return(data_stack)
}
stack_current_yr <- function(current_parcel, current_parcel_feature_layers, yr){
current_parcel[yr, ] = current_parcel_feature_layers
return(current_parcel)
}
write_frames <- function(data_stack, filetype, mov_folder, parcels, global_params, current_simulation_params, offset_site_indexes, offset_yrs){
# gray.colors(n = 1024, start = 0, end = 1, gamma = 2.2, alpha = NULL)
graphics.off()
rgb.palette <- colorRampPalette(c("black", "green"), space = "rgb")
if (!dir.exists(mov_folder)){
dir.create(mov_folder)
}
filename = paste0(mov_folder, "tmp%03d.", filetype, sep = '')
if (filetype == 'png'){
png(filename, height = parcels$landscape_dims[1], width = parcels$landscape_dims[2])
} else if (filetype == 'jpg'){
jpeg(filename, height = parcels$landscape_dims[1], width = parcels$landscape_dims[2])
}
for (yr in seq(current_simulation_params$time_steps)){
current_feature_layers = lapply(seq_along(data_stack), function(i) data_stack[[i]][yr, ])
landscape = array(0, parcels$landscape_dims)
landscape[unlist(parcels$land_parcels)] = unlist(current_feature_layers)
if (global_params$write_offset_layer == TRUE){
offset_sites_to_use = which(offset_yrs <= yr)
landscape[unlist(parcels$land_parcels[offset_site_indexes[offset_sites_to_use] ])] = current_simulation_params$max_eco_val
}
image(landscape, zlim = c(0, current_simulation_params$max_eco_val), col = rgb.palette(512)) #, col = grey(seq(0, 1, length = 256))
cat('\n image ', yr, ' of ', current_simulation_params$time_steps, 'done\n')
}
dev.off()
}
#write all offset parcels to single layer to output as image
write_site_mask <- function(output_filename, landscape_dims, land_parcels, current_site_indexes){
site_mask = array(0, landscape_dims)
site_mask[ unlist(land_parcels[unlist(current_site_indexes)])] = 1
# rgb.palette <- colorRampPalette(color_vec, space = "rgb")
png(filename = output_filename, height = dim(site_mask)[1], width = dim(site_mask)[2])
image(site_mask, zlim = c(0,1), col = rgb.palette(512))
dev.off()
return(site_mask)
}
# remove site from current offset pool
remove_site_from_pool <- function(offset_pool_object, current_site_indexes){
subset_pool_to_remove <- list_intersect(offset_pool_object$site_indexes, current_site_indexes)
if (length(subset_pool_to_remove$match_ind) > 0){
subset_pool_to_use <- seq_along(offset_pool_object$site_indexes)[-subset_pool_to_remove$match_ind]
offset_pool_object <- select_pool_subset(offset_pool_object, subset_pool_to_use)
}
return(offset_pool_object)
}
#sample over uniform random vector, indicies less than the threshold level are selected for clearing
select_sites_to_clear <- function(site_indexes, current_simulation_params){
clearing_thresh <- rep(current_simulation_params$unregulated_loss_prob, length(site_indexes))
discrim <- (runif(length(clearing_thresh)) < clearing_thresh)
inds_to_clear <- site_indexes[discrim]
return(inds_to_clear)
}
perform_unregulated_loss <- function(current_feature_layers, index_object, yr, current_simulation_params, decline_rates_initial, time_horizon){
if (current_simulation_params$unregulated_loss_prob == 0){
# return null object when clearing is inactive
return()
}
available_site_indexes = setdiff(index_object$global_indexes, unique(unlist(index_object$site_indexes_used)))
inds_to_clear <- select_sites_to_clear(available_site_indexes, current_simulation_params)
if (length(inds_to_clear) == 0){ #return null for no sites selected for clearing
return()
} else {
flog.info('unregulated losses %s' , as.vector(inds_to_clear))
}
# current remaining sites - used for calculation of cfac
parcel_num_remaining = length(available_site_indexes)
# store group of site characteristics in site characteristics object
unregulated_loss_object <- record_current_parcel_set(current_feature_layers[inds_to_clear],
inds_to_clear,
parcel_num_remaining,
yr)
# record characteristics of cleared site
unregulated_loss_object <- assess_current_pool(pool_object = unregulated_loss_object,
pool_type = 'devs',
calc_type = current_simulation_params$dev_calc_type,
cfacs_flag = current_simulation_params$dev_cfacs_flag,
adjust_cfacs_flag = current_simulation_params$adjust_dev_cfacs_flag,
action_type = current_simulation_params$offset_action_type,
include_potential_developments = current_simulation_params$include_potential_developments_in_dev_calc,
include_potential_offsets = current_simulation_params$include_potential_offsets_in_dev_calc,
include_unregulated_loss = current_simulation_params$include_unregulated_loss_in_dev_calc,
time_horizon_type = 'future',
current_simulation_params,
decline_rates_initial,
time_horizon,
yr) #determine future development parcel attributes
return(unregulated_loss_object)
}
# series of routines to implement offset
perform_offset_routine <- function(simulation_outputs, index_object, decline_rates, current_offset_object, current_offset_indexes, match_object,
current_simulation_params){
# if running in banking mode remove offset site from available bank
if (current_simulation_params$use_offset_bank == TRUE){
banked_offset_pool = simulation_outputs$offset_bank_object$site_indexes
banked_offset_inds_used = list_intersect(banked_offset_pool, current_offset_indexes)
# determine parcels used in matching routine and remove from available pool
simulation_outputs$offset_bank_object$site_indexes = remove_index(banked_offset_pool, banked_offset_inds_used$match_ind)
} else {
# determine parcels used in matching routine and remove from available pool
simulation_outputs$index_object <- update_index_object(index_object, update_type = 'offset', current_offset_indexes)
simulation_outputs$decline_rates <- update_decline_rates(decline_rates,
restoration_rate = current_simulation_params$restoration_rate,
restoration_rate_std = current_simulation_params$restoration_rate_std,
sample_restoration_rate = current_simulation_params$sample_restoration_rate,
features_to_use_in_offset_intervention = current_simulation_params$features_to_use_in_offset_intervention,
feature_num = current_simulation_params$feature_num,
decline_rate_type = 'offset',
action_type = current_simulation_params$offset_action_type,
current_offset_indexes) # set elements in decline rates array corresponding to offsets to restoration rates
}
#record current offset site characteristics
simulation_outputs$offsets_object <- append_current_group(simulation_outputs$offsets_object, current_offset_object, append_routine = 'standard')
#remove offset sites from available pool
simulation_outputs$offset_pool_object <- remove_parcel_from_current_pool(simulation_outputs$offset_pool_object,
current_site_indexes = current_offset_indexes)
return(simulation_outputs)
}
# remove site characteristics from current pool of available sites
remove_parcel_from_current_pool <- function(offset_pool_object, current_site_indexes){
# work out indexes of current sites in pool
sites_to_remove <- list_intersect(offset_pool_object$site_indexes, current_site_indexes)
if (length(sites_to_remove$match_ind) > 0){
# remove site index from available pool
subset_pool_to_use <- seq_along(offset_pool_object$site_indexes)[-sites_to_remove$match_ind]
# redefine pool with current used sites removed
offset_pool_object <- select_pool_subset(offset_pool_object, subset_pool_to_use)
}
return(offset_pool_object)
}
# routines to mark and destroy feature_layers in cleared sites e.g. Development or unregulated
perform_clearing_routine <- function(simulation_outputs, index_object, decline_rates, current_dev_object, clearing_type, current_simulation_params){
#remove development parcels from available pool
simulation_outputs$index_object <- update_index_object(index_object, clearing_type, current_dev_object$site_indexes)
if (clearing_type == 'development'){
#record current development site characteristics
simulation_outputs$dev_object <- append_current_group(simulation_outputs$dev_object, current_dev_object, append_routine = 'standard')
} else if (clearing_type == 'develop_from_credit'){
#record current development site characteristics
simulation_outputs$credit_object <- append_current_group(simulation_outputs$credit_object, current_dev_object, append_routine = 'standard')
} else if (clearing_type == 'unregulated'){
#record current cleared site characteristics
simulation_outputs$unregulated_loss_object <- append_current_group(simulation_outputs$unregulated_loss_object, current_dev_object, append_routine = 'standard')
}
# set elements corresponding to developed parcels in decline rates array to zero
simulation_outputs$decline_rates <- update_decline_rates(decline_rates,
restoration_rate = vector(),
restoration_rate_std = current_simulation_params$restoration_rate_std,
sample_restoration_rate = vector(),
features_to_use_in_offset_intervention = current_simulation_params$features_to_use_in_offset_intervention,
feature_num = current_simulation_params$feature_num,
decline_rate_type = 'development',
action_type = vector(),
current_dev_object$site_indexes)
if (length(current_dev_object$site_indexes) > 0){
# if any development was allowed remove developed site from available offset pool
simulation_outputs$offset_pool_object <- remove_parcel_from_current_pool(simulation_outputs$offset_pool_object,
current_site_indexes = current_dev_object$site_indexes)
}
return(simulation_outputs)
}
# determine current available credit accumulated through offsets for development
assess_banking_credit <- function(simulation_outputs, current_simulation_params){
features_to_use_in_offset_calc = current_simulation_params$features_to_use_in_offset_calc
# determine total offset gains
offset_credit = nested_list_sum(simulation_outputs$offset_pool_object$parcel_vals_used)
dev_list = append(simulation_outputs$credit_object$parcel_vals_used, simulation_outputs$dev_object$parcel_vals_used)
if (length(dev_list) > 0){
# determine total development losses
dev_sum = nested_list_sum(dev_list)
current_credit = subtract_nested_lists(offset_credit, dev_sum)
} else{
current_credit = offset_credit
}
return(current_credit)
}
# determine characteristics of potential offset sites
prepare_offset_pool <- function(simulation_outputs, current_simulation_params,
parcels, decline_rates_initial, yr){
# if no developments or banked offsets for the current year return null object
if (current_simulation_params$intervention_vec[yr] == 0){
offset_pool_object = list()
return(offset_pool_object)
}
# select current set of available offset sites
current_offset_pool <- simulation_outputs$index_object$indexes_to_use$offsets
# if pool is empty return null object and print error
if (length(current_offset_pool) == 0){
flog.error('empty offset pool flag')
offset_pool_object = list()
return(offset_pool_object)
}
if (current_simulation_params$use_offset_bank == TRUE){
# if running in offset bank mode select sites from current region
flog.error('offset bank in development')
stop()
subset_pool = simulation_outputs$offset_bank_object$site_indexes
# find set of offset characteristics that apply to current set of available sites
offset_pool_object <- select_pool_subset(simulation_outputs$offset_bank_object, subset_pool)
# find set of current cumulative site vals, record as projected val as calculation is from time
# of original offset to current time for banking
offset_pool_object$projected_vals <- find_current_parcel_sums(simulation_outputs$current_feature_layers[unlist(simulation_outputs$offset_bank_object$site_indexes[subset_pool])])
offset_pool_type = 'offset_bank'
} else {
# when running in standard mode record current offsite site characteristics
offset_pool_object <- record_current_parcel_set(simulation_outputs$current_feature_layers[current_offset_pool],
current_offset_pool,
parcel_num_remaining = length(current_offset_pool),
yr) #arrange available parcel pool into form to use in parcel set determination
offset_pool_type = 'offsets'
}
# determine current gains characteristics
offset_pool_object <- assess_current_pool(pool_object = offset_pool_object,
pool_type = offset_pool_type,
calc_type = current_simulation_params$offset_calc_type,
cfacs_flag = current_simulation_params$offset_cfacs_flag,
adjust_cfacs_flag = current_simulation_params$adjust_offset_cfacs_flag,
action_type = current_simulation_params$offset_action_type,
include_potential_developments = current_simulation_params$include_potential_developments_in_offset_calc,
include_potential_offsets = current_simulation_params$include_potential_offsets_in_offset_calc,
include_unregulated_loss = current_simulation_params$include_unregulated_loss_in_offset_calc,
time_horizon_type = current_simulation_params$offset_time_horizon_type,
current_simulation_params,
decline_rates_initial,
time_horizon = current_simulation_params$offset_time_horizon,
yr) #determine available parcel values, depending on what particular offset policy is in use using counterfactuals etc.
return(offset_pool_object)
}
# routines to perform offset banking
perform_banking_routine <- function(simulation_outputs, current_simulation_params, yr){
# how many offsets to be added in current year
offset_bank_num = unlist(current_simulation_params$banked_offset_vec[yr])
if (offset_bank_num == 0){
return(simulation_outputs)
}
# select current number of offset sites from current available pool to add to banked offset pool
current_banked_offset_pool <- sample(simulation_outputs$index_object$indexes_to_use$offsets, offset_bank_num)
# number of sites to potentially offset
parcel_num_remaining = length(c(unlist(simulation_outputs$index_object$indexes_to_use$offsets),
unlist(simulation_outputs$index_object$indexes_to_use$devs)))
# record current offset pool characteristics
current_banked_object <- record_current_parcel_set(simulation_outputs$current_feature_layers[current_banked_offset_pool],
current_pool = current_banked_offset_pool,
parcel_num_remaining,
yr) # arrange current parcel data
simulation_outputs$offset_bank_object = append_current_group(object_to_append = simulation_outputs$offset_bank_object,
current_object = current_banked_object,
append_routine = 'banked_offset')
# remove current group of sites from available pool
simulation_outputs$index_object <- update_index_object(simulation_outputs$index_object,
update_type = 'banking',
current_banked_offset_pool)
simulation_outputs$decline_rates <- update_decline_rates(simulation_outputs$decline_rates, # set decline rate parameters to offset
restoration_rate = current_simulation_params$restoration_rate,
restoration_rate_std = current_simulation_params$restoration_rate_std,
sample_restoration_rate = current_simulation_params$sample_restoration_rate,
features_to_use_in_offset_intervention = current_simulation_params$features_to_use_in_offset_intervention,
feature_num = current_simulation_params$feature_num,
decline_rate_type = 'offset',
action_type = current_simulation_params$offset_action_type,
current_banked_offset_pool)
return(simulation_outputs)
}
#assess whether parcel sets, ie. developments and offsets have achieved NNL according to gains structures
# current_feature_layers = simulation_outputs$current_feature_layers
# offsets_object = simulation_outputs$offsets_object
# offset_parcel_sets = simulation_outputs$index_object$site_indexes_used$offsets
# time_horizon = current_simulation_params$offset_time_horizon
assess_parcel_sets <- function(current_feature_layers, offsets_object, offset_parcel_sets, current_simulation_params, decline_rates_initial, time_horizon, yr){
assessed_parcel_sets_object <- list()
parcel_set_count = length((offset_parcel_sets))
discriminator_metric <- vector('list', parcel_set_count)
for (parcel_set_ind in seq_len(parcel_set_count)){
# select current set of sites from offsets object
current_parcel_set <- which(unlist(offsets_object$site_indexes) %in% unlist(offset_parcel_sets[[parcel_set_ind]]))
# select characteristics of current set of sites
current_offset_object <- lapply(seq_along(offsets_object), function(i) offsets_object[[i]][current_parcel_set])
names(current_offset_object) <- names(offsets_object)
parcel_vals_achieved <- assess_current_gain_pool(current_feature_layers,
pool_object = current_offset_object,
pool_type = "offset_bank",
calc_type = current_simulation_params$offset_calc_type,
cfacs_flag = current_simulation_params$offset_cfacs_flag,
adjust_cfacs_flag = current_simulation_params$adjust_offset_cfacs_flag,
include_potential_developments = current_simulation_params$include_potential_developments_in_offset_calc,
include_potential_offsets = current_simulation_params$include_potential_offsets_in_offset_calc,
include_unregulated_loss = current_simulation_params$include_unregulated_loss_in_offset_calc,
time_horizon_type = current_simulation_params$offset_time_horizon_type,
current_simulation_params,
decline_rates_initial,
time_horizon,
yr)
for (feature_ind in seq_len(current_simulation_params$feature_num)){
discriminator_metric[[parcel_set_ind]] = nested_list_sum(subtract_nested_lists(parcel_vals_achieved, current_offset_object$parcel_vals_used))
}
}
# determine whether site has reached gains for subset of features used in offset calculations
parcel_set_success_inds <- unlist(lapply(seq_along(discriminator_metric),
function(i) (all(unlist(discriminator_metric[[i]][current_simulation_params$features_to_use_in_offset_calc]) > 0))))
site_success_inds = which(unlist(offsets_object$site_indexes) %in% unlist(offset_parcel_sets[parcel_set_success_inds]))
assessed_parcel_sets_object$site_success_inds <- offsets_object$site_indexes[site_success_inds]
assessed_parcel_sets_object$discriminator_metric <- discriminator_metric
return(assessed_parcel_sets_object)
}
assess_current_gain_pool <- function(current_feature_layers, pool_object, pool_type, calc_type, cfacs_flag, adjust_cfacs_flag,
include_potential_developments,include_potential_offsets,include_unregulated_loss,
time_horizon_type, current_simulation_params, decline_rates_initial, time_horizon, yr){
current_pool = unlist(pool_object$site_indexes)
parcel_count = length(current_pool)
offset_yrs = unlist(pool_object$offset_yrs)
time_horizons <- generate_time_horizons(project_type = 'current', yr, offset_yrs, time_horizon, parcel_count)
if (cfacs_flag == TRUE){
cfacs_object = calc_cfacs(parcel_feature_layers = pool_object$parcel_feature_layers,
parcel_num_remaining = pool_object$parcel_num_remaining,
current_simulation_params,
decline_rates = decline_rates_initial[current_pool],
time_horizons,
offset_yrs,
include_potential_developments,
include_potential_offsets,
include_unregulated_loss,
adjust_cfacs_flag = current_simulation_params$adjust_offset_cfacs_flag,
features_to_project = current_simulation_params$features_to_use_in_offset_calc)
cfac_vals = nested_list_tail(cfacs_object$cfacs_to_use)
}
projected_vals = current_feature_layers[current_pool]
projected_vals = lapply(seq_along(projected_vals), function(i) lapply(seq_along(current_simulation_params$features_to_use_in_offset_calc), function(j) sum(projected_vals[[i]][[j]] )))
parcel_vals_achieved = evaluate_parcel_vals(calc_type, pool_object$parcel_sums_at_offset, projected_vals, cfac_vals)
return(parcel_vals_achieved)
}
#construct a list of zero arrays with identical dimensions defined by array_dims
list_of_zeros <- function(list_dims, array_dims){
list_object = vector('list', list_dims)
for (list_ind in seq_len(list_dims)){
list_object[[list_ind]] = array(0, array_dims)
}
return(list_object)
}
#used to break up array into smaller set of sub arrays defined by vx and vy that fit together to give input array
mcell <- function(x, vx, vy){
rowsizes = vy;
colsizes = vx;
rows = length(rowsizes);
cols = length(colsizes);
a = 1
# make an array composed of lists with dimenisons that define the land parcels/regions.
# The list format allows arrays of different sizes to be stored
B = vector('list', rows*cols)
colStart = 0
# run down through the columns of input array
for (i in seq_len(cols)){
rowStart = 0
for (j in seq_len(rows)){
#group elements of input array into sub arrays and assign to B
B[[a]] = x[rowStart+(1:rowsizes[j]), colStart+(1:colsizes[i])]
rowStart = rowStart + rowsizes[j]
a = a + 1
}
colStart = colStart + colsizes[i]
}
parcel = list()
parcel$dims = c(length(vy), length(vx))
parcel$elements = B
return(parcel)
}
user_projection <- function(parcel_vals, min_eco_val, max_eco_val, decline_rate, time_horizon, time_fill){
if (time_fill == TRUE){
time_vec = 0:time_horizon
} else {time_vec = time_horizon}
# If decline_rate is set to zero return same value for all time steps
if (decline_rate == 0){
eco_projected = rep(parcel_vals, length(time_vec))
} else if ( (decline_rate != 0) & (parcel_vals > min_eco_val)){
# define logistic curve given logistic parameter decline_rate.
t_sh = -1/decline_rate * log( ((parcel_vals - min_eco_val)/(max_eco_val - parcel_vals)))
eco_projected = min_eco_val + (max_eco_val - min_eco_val)/(1 + exp(-decline_rate*(time_vec - t_sh)))
}
return(eco_projected)
}
# define feature_layers dynamics by logistic curve
logistic_projection <- function(parcel_vals, min_eco_val, max_eco_val, current_dec_rate, time_vec){
t_sh = -1/current_dec_rate * log( ((parcel_vals - min_eco_val)/(max_eco_val - parcel_vals)))
# define logistic curve given logistic parameter set.
eco_projected = min_eco_val + (max_eco_val - min_eco_val)/(1 + exp(-current_dec_rate*(time_vec - t_sh)))
return(eco_projected)
}
user_projection <- function(parcel_vals, min_eco_val, max_eco_val, current_dec_rate, time_vec){
lo <- loess(y ~ time_vec)
return(eco_projected)
}
# calculate the expected feature_layers under maintenaince, restoration
project_feature_layers <- function(projection_type, current_parcel_feature_layers, min_eco_val, max_eco_val, current_dec_rate, time_horizon, time_fill){
# for maintain feature_layers or project development vals, copy current feature_layers to matrix of depth (time_horizon + 1)
if ( (current_dec_rate == 1) || (current_dec_rate == 0) || ( sum(current_parcel_feature_layers) == 0) ){
# for maintain feature_layers copy current feature_layers to matrix of depth (time_horizon + 1)
if (time_fill == TRUE){
# return all feature_layers states over all time steps
current_parcel_feature_layers = matrix(rep(current_parcel_feature_layers, (time_horizon + 1)), ncol = length(current_parcel_feature_layers), byrow = TRUE)
} else {
# project for single time step
return(current_parcel_feature_layers)
}
} else {
if (time_fill == TRUE){
time_vec = 0:time_horizon
} else {
time_vec = time_horizon
}
# update feature_layers according to function defined in project_feature_layers function
# function parameters are contained in decline_rates array
if (projection_type == 'logistic_function'){
current_parcel_feature_layers = sapply(current_parcel_feature_layers, logistic_projection, min_eco_val, max_eco_val,
current_dec_rate, time_vec)
}
}
if (time_horizon == 0){
dim(current_parcel_feature_layers) = c(1, length(current_parcel_feature_layers))
}
return(current_parcel_feature_layers)
}
# project sites through all features - returns nested list of arrays where each nested array has length
# defined by current_time_horizons and depth defined by feature number for all sites
project_sites <- function(current_parcel_feature_layers, current_decline_rates, current_time_horizons, current_simulation_params, features_to_project, time_fill){
projected_features = lapply(seq_along(current_parcel_feature_layers),
function(i) lapply(features_to_project,
function(j) project_feature_layers(current_simulation_params$projection_type,
current_parcel_feature_layers[[i]][[j]],
current_simulation_params$min_eco_val,
current_simulation_params$max_eco_val,
current_dec_rate = current_decline_rates[[i]][[j]],
time_horizon = unlist(current_time_horizons[i]),
time_fill)))
return(projected_features)
}
# find sum nested list of depth j with identical structure
nested_list_sum <- function(nested_list){
if (length(nested_list) > 0){
summed_list <- lapply(seq_along(nested_list[[1]]),
function(j) Reduce('+', lapply(seq_along(nested_list), function(i) nested_list[[i]][[j]])))
return(summed_list)
} else {
return(NULL)
}
}
#build nested list with length outer_dim and depth defined by inner_dim
generate_nested_list <- function(outer_dim, inner_dim){
if (outer_dim > 0){
nested_list <- vector('list', outer_dim)
} else {
nested_list = list()
}
for (outer_ind in seq_len(outer_dim)){
nested_list[[outer_ind]] <- vector('list', inner_dim)
}
return(nested_list)
}
remove_index <- function(object_list, ind_to_remove){
if (length(ind_to_remove) > 0){
object_list <- object_list[-ind_to_remove]
}
return(object_list)
}
recalculate_probabilities <- function(current_probability_list){
current_probability_list = lapply(seq_along(current_probability_list), function(i) current_probability_list[[i]]/sum(unlist(current_probability_list)))
return(current_probability_list)
}
match_sites <- function(offset_pool_object, current_credit, dev_probability_list, offset_probability_list, current_simulation_params,
intervention_vec, indexes_to_use, current_feature_layers, decline_rates_initial,
land_parcels, yr, time_horizon){
match_object = false_match()
current_pool_vals = offset_pool_object$parcel_vals_used
current_pool_indexes = offset_pool_object$site_indexes
parcel_num_remaining = length(indexes_to_use)
current_match_pool = indexes_to_use
if ((length(offset_pool_object$site_indexes) == 0) | (parcel_num_remaining == 0)){
return(match_object)
}
current_pool_vals_array <- matrix(unlist(current_pool_vals), nrow = length(current_pool_vals), byrow=TRUE)
zero_inds <- which(apply(current_pool_vals_array, MARGIN = 1, sum) == 0)
###### TODO check zero development routines the code below may not be the most efficient method for this ###
current_pool_vals = remove_index(current_pool_vals, zero_inds)
current_pool_indexes = remove_index(current_pool_indexes, zero_inds)
# store group of site characteristics in site characteristics object
dev_pool_object <- record_current_parcel_set(current_feature_layers[current_match_pool],
current_match_pool,
parcel_num_remaining,
yr)
dev_pool_object <- assess_current_pool(pool_object = dev_pool_object,
pool_type = 'devs',
calc_type = current_simulation_params$dev_calc_type,
cfacs_flag = current_simulation_params$dev_cfacs_flag,
adjust_cfacs_flag = current_simulation_params$adjust_dev_cfacs_flag,
action_type = current_simulation_params$offset_action_type,
include_potential_developments = current_simulation_params$include_potential_developments_in_dev_calc,
include_potential_offsets = current_simulation_params$include_potential_offsets_in_dev_calc,
include_unregulated_loss = current_simulation_params$include_unregulated_loss_in_dev_calc,
time_horizon_type = 'future',
current_simulation_params,
decline_rates_initial,
time_horizon,
yr) #determine future development parcel attributes
current_match_vals_pool = dev_pool_object$parcel_vals_used
if (sum(unlist(current_match_vals_pool))== 0){
flog.info('all projected developments are zero - blocking all developments')
match_object = false_match()
match_object$offset_object = list()
match_object$current_credit = current_credit
return(match_object)
}
#any potential parcel set match requires a minimum of two sites
while( (match_object$match_flag == FALSE) & (length(current_match_pool) > 1) ){
if (current_simulation_params$development_selection_type == 'random'){
sample_ind = sample(seq_along(current_match_pool), size = 1)
} else if (current_simulation_params$development_selection_type == 'weighted'){
current_dev_probability_list = recalculate_probabilities(dev_probability_list[unlist(current_match_pool)])
sample_ind = sample(seq_along(current_match_pool), size = 1, prob = current_dev_probability_list)
}
current_test_index = current_match_pool[sample_ind]
vals_to_match = current_match_vals_pool[[sample_ind]]
if (current_simulation_params$use_offset_bank == FALSE){
dev_ind = list_intersect(current_pool_indexes, current_test_index) #find and remove index that corresponds to potiential development index
match_pool_to_use = remove_index(current_pool_indexes, dev_ind$match_ind)
pool_vals_to_use = remove_index(current_pool_vals, dev_ind$match_ind)
} else {
match_pool_to_use = current_pool_indexes #if performing offset banking use any of the available banked offset pool
pool_vals_to_use = current_pool_vals
}
match_object <- match_from_pool(match_type = 'offset',
current_pool = match_pool_to_use,
pool_vals_to_use,
current_credit,
offset_probability_list,
vals_to_match_initial = vals_to_match,
current_simulation_params,
yr) #perform matching routine
if (match_object$match_flag == FALSE){
# only try to match sites with smaller ecological value than current site
inds_to_keep = which(lapply(seq_along(dev_pool_object$parcel_vals_used),
function(i) all(unlist(subtract_nested_lists(dev_pool_object$parcel_vals_used[[i]], vals_to_match)) < 0) ) == TRUE)
current_match_pool = dev_pool_object$site_indexes[inds_to_keep]
current_match_vals_pool = dev_pool_object$parcel_vals_used[inds_to_keep]
}
}
if (match_object$match_flag == TRUE){
dev_match_index = which(unlist(dev_pool_object$site_indexes) == current_test_index)
match_object$development_object = select_pool_subset(dev_pool_object, unlist(dev_match_index))
subset_pool = list_intersect(offset_pool_object$site_indexes, match_object$match_indexes)
offset_object <- select_pool_subset(pool_object = offset_pool_object, subset_pool = subset_pool$match_ind)
match_object$offset_object = offset_object
} else if (match_object$match_flag == FALSE){
match_object$offset_object = list()
match_object$current_credit = current_credit
}
return(match_object)
}
develop_from_credit <- function(current_feature_layers, current_credit, dev_probability_list, current_simulation_params,
intervention_vec, dev_indexes_to_use, decline_rates_initial, land_parcels, yr, time_horizon){
parcel_num_remaining = length(dev_indexes_to_use)
# store group of site characteristics in site characteristics object
dev_pool_object <- record_current_parcel_set(current_feature_layers[dev_indexes_to_use], dev_indexes_to_use, parcel_num_remaining, yr)
dev_pool_object <- assess_current_pool(pool_object = dev_pool_object,
pool_type = 'devs',
calc_type = current_simulation_params$dev_calc_type,
cfacs_flag = current_simulation_params$dev_cfacs_flag,
adjust_cfacs_flag = current_simulation_params$adjust_dev_cfacs_flag,
action_type = current_simulation_params$offset_action_type,
include_potential_developments = current_simulation_params$include_potential_developments_in_dev_calc,
include_potential_offsets = current_simulation_params$include_potential_offsets_in_dev_calc,
include_unregulated_loss = current_simulation_params$include_unregulated_loss_in_dev_calc,
time_horizon_type = 'future',
current_simulation_params,
decline_rates_initial,
time_horizon,
yr)
# subset_pool = list_intersect(dev_pool_object$site_indexes, match_object$match_indexes)
# match_object$development_object = select_pool_subset(dev_pool_object, subset_pool = subset_pool$match_ind)
if (length(unlist(dev_pool_object$site_indexes)) > 0){
pool_vals_to_use = dev_pool_object$parcel_vals_use
match_object <- match_from_pool(match_type = 'development',
current_pool = dev_pool_object$site_indexes,
pool_vals_to_use,
current_credit,
dev_probability_list,
vals_to_match_initial = current_credit,
current_simulation_params,
yr)
} else{
match_object = false_match()
}
if (match_object$match_flag == TRUE){
subset_pool = list_intersect(dev_pool_object$site_indexes, match_object$match_indexes)
match_object$development_object = select_pool_subset(dev_pool_object, subset_pool = subset_pool$match_ind)
} else{
match_object$development_object = list()
match_object$current_credit = current_credit
}
return(match_object)
}
evaluate_parcel_vals <- function(calc_type, current_condition_vals, projected_vals, cfac_vals){
if (calc_type == 'current_condition'){
projected_vals = current_condition_vals
cfac_vals = list_of_zeros(length(current_condition_vals), 1)
} else if (calc_type == 'restoration_gains'){
cfac_vals = current_condition_vals
} else if (calc_type == 'restoration_condition_value'){
cfac_vals = list_of_zeros(length(projected_vals), 1)
} else if (calc_type == 'avoided_condition_decline'){
projected_vals = current_condition_vals
} else if ((calc_type == 'future_condition')){
projected_vals = cfac_vals
cfac_vals = list_of_zeros(length(current_condition_vals), 1)
}
parcel_vals_pool = subtract_nested_lists(projected_vals, cfac_vals)
return(parcel_vals_pool)
}
subtract_nested_lists <- function(list_a, list_b){
length_a = unlist(lapply(seq_along(list_a), function(i) length(list_a[[i]])))
length_b = unlist(lapply(seq_along(list_b), function(i) length(list_b[[i]])))
stopifnot(all(length_a == length_b))
list_c = lapply( seq_along(list_a), function(i) mapply('-', list_a[[i]], list_b[[i]], SIMPLIFY = FALSE))
return(list_c)
}
sum_nested_lists <- function(list_to_sum){
summed_list = list_to_sum[[1]]
if (length(list_to_sum) > 1){
for (list_ind in 2:length(list_to_sum)){
current_list = list_to_sum[[list_ind]]
summed_list = lapply(seq_along(summed_list), function(i) mapply('+', summed_list[[i]], current_list[[i]], SIMPLIFY = FALSE))
}
}
return(summed_list)
}
subtract_lists <- function(list_a, list_b){
list_c = mapply('-', list_a, list_b, SIMPLIFY = FALSE)
return(list_c)
}
sum_lists <- function(list_a, list_b){
if ( (length(list_a) == length(list_b)) & (length(list_a) >0)){
list_c = mapply('+', list_a, list_b, SIMPLIFY = FALSE)
} else if (length(list_a) == 0){
if (length(list_b) == 0){
list_c = list()
} else{
list_c = list_b
}
} else if (length(list_b) == 0){
if (length(list_a) == 0){
list_c = list()
} else{
list_c = list_a
}
}
return(list_c)
}
# store group of site characteristics in parcel_set_object
record_current_parcel_set <- function(current_feature_layers, current_pool, parcel_num_remaining, yr){
parcel_set_object = list()
parcel_set_object$offset_yrs = rep(list(yr), length(current_pool))
parcel_set_object$parcel_feature_layers = current_feature_layers
parcel_set_object$parcel_sums_at_offset = sum_sites(current_feature_layers)
parcel_set_object$site_indexes = as.list(current_pool)
parcel_set_object$parcel_num_remaining = rep(list(parcel_num_remaining), length(current_pool))
return(parcel_set_object)
}
# determine last element of all nested list vectors
nested_list_tail <- function(list_a){
last_elements <- lapply(seq_along(list_a),
function(i) unlist(lapply(seq_along(list_a[[i]]),
function(j) list_a[[i]][[j]][ length(list_a[[i]][[j]]) ] )))
return(last_elements)
}
##### TODO probability list should be passed to this - otherwise underlying assumption is equal weights for all sites which is wrong #####
find_intervention_probability <- function(intervention_vec, offset_yrs, calc_type, offset_intervention_scale,
time_horizons, parcel_num, parcel_num_remaining, time_steps){
intervention_probs = vector('list', parcel_num)
parcel_num_remaining = unlist(parcel_num_remaining)
offset_yrs = unlist(offset_yrs)
for (parcel_ind in seq_len(parcel_num)){
time_horizon = time_horizons[parcel_ind]
offset_yr = offset_yrs[parcel_ind]
current_prob = array(0, (time_horizon + 1))
current_vec = intervention_vec[offset_yr:time_steps]
if (length(current_vec) < (time_horizon + 1)){
current_vec = c(current_vec, array(0, ((time_horizon + 1) - length(current_vec))))
}
current_vec = current_vec[1:(time_horizon + 1)]
current_parcel_num_remaining = parcel_num_remaining[parcel_ind]
if (current_parcel_num_remaining > 0) {
current_prob = current_prob + current_vec/current_parcel_num_remaining
}
if (calc_type == 'offset'){
current_prob = offset_intervention_scale*current_prob
}
intervention_probs[[parcel_ind]] = current_prob
}
return(intervention_probs)
}
# test to determine nearest neighbour by Euclidean norm given a pool of potential candidates (parcel_vals_pool)
# and the criterion to match to (vals_to_match)
euclidean_norm_match <- function(parcel_vals_pool, vals_to_match){
euclidean_norm = lapply(seq_along(parcel_vals_pool), function(i) parcel_vals_pool[[i]] - vals_to_match)
#build euclidean metric to perform euclidean match
euclidean_norm = unlist(lapply(seq_along(euclidean_norm), function(i) sum(euclidean_norm[[i]]^2)))
match_ind = which(euclidean_norm == min(euclidean_norm))
if (length(match_ind) > 1){
#list where the duplicate occurred
flog.error(cat('duplicate match flag on ', parcel_vals_pool[match_ind]), '\n')
# chose first site of the multiple match
match_ind = sample(match_ind, 1)
}
match_vals = parcel_vals_pool[match_ind]
match_object = list()
match_object$match_vals = match_vals
match_object$match_ind = match_ind
return(match_object)
}
#return false match object
false_match <- function(){
match_object = list()
match_object$match_flag = FALSE
return(match_object)
}
select_cols <- function(arr_to_use, col_inds){
arr_to_use <- arr_to_use[, col_inds]
arr_to_use <- t(t(arr_to_use))
return(arr_to_use)
}
select_pool_to_match <- function(vals_to_match, current_features_to_use_in_offset_calc, use_specified_offset_metric, metric_type, thresh, pool_vals_to_use,
max_parcel_num, current_pool, match_type, screen_site_zeros){
pool_object = list()
if (length(unlist(pool_vals_to_use)) == 0){
pool_object$break_flag = TRUE
return(pool_object)
}
pool_vals_to_use <- lapply(seq_along(pool_vals_to_use), function(i) pool_vals_to_use[[i]][current_features_to_use_in_offset_calc])
if (use_specified_offset_metric == TRUE){
pool_vals_to_use = transform_features_to_offset_metric(pool_vals_to_use, metric_type)
}
zero_inds <- which(unlist(lapply(seq_along(pool_vals_to_use), function(i) sum(pool_vals_to_use[[i]]) == 0)))
if (screen_site_zeros == TRUE){
if (length(zero_inds) > 0){
current_pool <- remove_index(current_pool, zero_inds)
pool_vals_to_use <- remove_index(pool_vals_to_use, zero_inds)
}
if (length(current_pool) == 0){
cat('\nall sites yield zero assessment')
pool_object$break_flag = TRUE
return(pool_object)
}
}
##### TODO check robustness #######
if (max_parcel_num > 1){
pool_vals_to_test = Reduce('+', pool_vals_to_use)
} else {
pool_vals_to_test = pool_vals_to_use
}
if (match_type == 'offset'){
pool_condition = unlist(lapply(seq_along(pool_vals_to_test), function(i) all( (vals_to_match - pool_vals_to_test[[i]]) <= thresh)))
} else if (match_type == 'development'){
pool_condition = unlist(lapply(seq_along(pool_vals_to_test), function(i) all( (vals_to_match - pool_vals_to_test[[i]]) >= -thresh)))
}
if (all(pool_condition == FALSE)){
pool_object$break_flag = TRUE
return(pool_object)
}
if (max_parcel_num == 1){
pool_vals_to_use <- pool_vals_to_use[pool_condition]
current_pool <- current_pool[pool_condition]
}
pool_object$break_flag = FALSE
pool_object$pool_vals_to_use = pool_vals_to_use
pool_object$current_pool = current_pool
pool_object$zero_inds = zero_inds
return(pool_object)
}
transform_features_to_offset_metric <- function(pool_vals, metric_type){
if (metric_type == 'euclidian_norm'){
if (class(pool_vals) == 'list'){
pool_vals = lapply(seq_along(pool_vals), function(i) sqrt(sum(pool_vals[[i]]^2)))
} else {
pool_vals = sqrt(sum(pool_vals^2))
}
}
return(pool_vals)
}
match_from_pool <- function(match_type, current_pool, pool_vals_to_use, current_credit, current_probability_list,
vals_to_match_initial, current_simulation_params, yr){
if (length(unlist(current_pool)) == 0){
match_object = false_match()
return(match_object)
}
metric_type = current_simulation_params$offset_metric_type
current_features_to_use_in_offset_calc = current_simulation_params$features_to_use_in_offset_calc
if (match_type == 'offset'){
max_parcel_num = current_simulation_params$max_offset_parcel_num
match_procedure = current_simulation_params$offset_selection_type
screen_site_zeros = current_simulation_params$screen_offset_zeros
if (current_simulation_params$use_specified_offset_metric == FALSE){
current_features_to_use_in_offset_calc = which(vals_to_match_initial > 0)
if ( length(current_features_to_use_in_offset_calc) == 0 ) {
match_object = false_match()
return(match_object)
}
} else {
current_features_to_use_in_offset_calc = current_simulation_params$features_to_use_in_offset_calc
}
offset_multiplier = current_simulation_params$offset_multiplier
vals_to_match = offset_multiplier * vals_to_match_initial[current_features_to_use_in_offset_calc]
if (current_simulation_params$use_specified_offset_metric == TRUE){
vals_to_match = transform_features_to_offset_metric(vals_to_match, metric_type = current_simulation_params$offset_metric_type)
}
if ( (current_simulation_params$allow_developments_from_credit == TRUE)){
if (current_simulation_params$use_specified_offset_metric == FALSE){
vals_to_match = vals_to_match - current_credit[current_features_to_use_in_offset_calc]
} else {
vals_to_match = vals_to_match - current_credit
}
}
} else if (match_type == 'development'){
screen_site_zeros = current_simulation_params$screen_dev_zeros
# force only single development for development routine
max_parcel_num = 1
if (screen_site_zeros == FALSE){
# if zeros are allowed use random selection
match_procedure = 'random'
} else {
match_procedure = current_simulation_params$development_selection_type
}
# when developing from credit, inverse offset multiplier
vals_to_match = 1/current_simulation_params$offset_multiplier*vals_to_match_initial
}
#create an array of threshold values defined by user based proportion
thresh = array(current_simulation_params$match_threshold_ratio * vals_to_match)
pool_object <- select_pool_to_match(vals_to_match, current_features_to_use_in_offset_calc, current_simulation_params$use_specified_offset_metric, metric_type, thresh,
pool_vals_to_use, max_parcel_num, current_pool, match_type, screen_site_zeros)
if (pool_object$break_flag == TRUE){
match_object = false_match()
return(match_object)
}
# initialise parcel_vals for matching procedure
parcel_vals_pool = pool_object$pool_vals_to_use
current_pool = pool_object$current_pool
match_flag = FALSE
match_vals = list()
match_indexes = list()
while(match_flag == FALSE){
if (length(current_pool) == 0){
break
}
if (match_procedure == 'greedy'){
match_params = euclidean_norm_match(parcel_vals_pool, vals_to_match)
} else if (match_procedure == 'random'){
match_params = list()
match_params$match_ind = sample(length(current_pool), 1)
match_params$match_vals = parcel_vals_pool[match_params$match_ind]
} else if (match_procedure == 'weighted'){
match_params = list()
probability_list_to_use = recalculate_probabilities(current_probability_list[unlist(current_pool)])
match_params$match_ind = sample(length(current_pool), 1, probability_list_to_use)
match_params$match_vals = parcel_vals_pool[match_params$match_ind]
}
current_match_val = unlist(match_params$match_vals)
current_match_index = current_pool[match_params$match_ind]
vals_to_match = vals_to_match - current_match_val
if (max_parcel_num > 1){
ind_to_remove = list_intersect(current_pool, current_match_index)
current_pool = remove_index(current_pool, ind_to_remove$match_ind)
parcel_vals_pool = remove_index(parcel_vals_pool, ind_to_remove$match_ind)
match_vals = append(match_vals, current_match_val)
match_indexes = append(match_indexes, current_match_index)
if (length(unlist(match_indexes)) > max_parcel_num){
match_flag = FALSE
break
}
} else {
match_indexes = list(current_match_index)
match_vals = list(current_match_val)
}
if (match_type == 'offset'){
match_flag = all(vals_to_match <= thresh)
} else if (match_type == 'development'){
match_flag = all(vals_to_match >= -thresh)
}
}
#### TODO NOTE if running with offset_metric specified current credit is wrong as it copies all entries to the same value
if (match_type == 'offset'){
# switch sign for any additional credit from offset
current_credit_to_update = -vals_to_match
} else if (match_type == 'development'){
current_credit_to_update = vals_to_match
}
if (current_simulation_params$use_specified_offset_metric == FALSE){
current_credit[current_features_to_use_in_offset_calc] = current_credit_to_update
} else {
current_credit = current_credit_to_update
}
match_object = list()
match_object$match_indexes = match_indexes
match_object$match_vals = match_vals
match_object$match_flag = match_flag
match_object$current_credit = current_credit
return(match_object)
}
# determine cumulative value of all sites within parcel feature_layers for multiple features
sum_sites <- function(current_feature_layers){
parcel_sums = lapply(seq_along(current_feature_layers),
function(i) unlist(lapply(seq_along(current_feature_layers[[i]]),
function(j) sum(current_feature_layers[[i]][[j]]) )))
return(parcel_sums)
}
# determine cumulative sum of all sites within parcel feature_layers for 1D features
sum_feature_layers <- function(parcel_feature_layers){
parcel_sums <- lapply(seq_along(parcel_feature_layers), function(i) sum(parcel_feature_layers[[i]] ))
}
#remove site from available pool for offsets and developments. This is a two stage process to cover when offsets and developments may not overlap
update_index_object <- function(index_object, update_type, site_indexes){
#remove parcel from available list
index_object$indexes_to_use$offsets = setdiff(index_object$indexes_to_use$offsets, site_indexes)
index_object$indexes_to_use$devs = setdiff(index_object$indexes_to_use$devs, site_indexes)
if (update_type == 'offset'){
index_object$site_indexes_used$offsets = append(index_object$site_indexes_used$offsets, list(site_indexes))
} else if (update_type == 'development'){
index_object$site_indexes_used$devs = append(index_object$site_indexes_used$devs, list(site_indexes))
} else if (update_type == 'unregulated'){
index_object$site_indexes_used$unregulated = append(index_object$site_indexes_used$unregulated, list(site_indexes))
} else if (update_type == 'develop_from_credit'){
index_object$site_indexes_used$dev_credit = append(index_object$site_indexes_used$dev_credits, list(site_indexes))
} else if (update_type == 'banking'){
index_object$site_indexes_used$banking = append(index_object$site_indexes_used$banking, list(site_indexes))
}
return(index_object)
}
# update feature_layers change parameters
# routine to label offset and development sites - updates decline_rates with 0 entry for developments,
# and 1/restoration_parameter for maintain/restoration offset. Protection offset leaves decline_rates unaffected
update_decline_rates <- function(decline_rates, restoration_rate, restoration_rate_std, sample_restoration_rate,
features_to_use_in_offset_intervention, feature_num, decline_rate_type, action_type, site_indexes){
for (current_parcel_ind in unlist(site_indexes)){
if (decline_rate_type == 'development'){
decline_rates[[current_parcel_ind]] <- rep(list(0), feature_num)
} else if (decline_rate_type == 'offset'){
if (action_type == 'maintain'){
decline_rates[[current_parcel_ind]][features_to_use_in_offset_intervention] = rep(list(1), length(features_to_use_in_offset_intervention))
} else if (action_type == 'restore'){
if (sample_restoration_rate == TRUE){
# spread restoration rates about mean
decline_rates[[current_parcel_ind]][features_to_use_in_offset_intervention] = as.list(rnorm(length(features_to_use_in_offset_intervention), mean = restoration_rate, sd = restoration_rate_std))
} else{
# use identical restoration rates for each site
decline_rates[[current_parcel_ind]][features_to_use_in_offset_intervention] = rep(list(restoration_rate), length(features_to_use_in_offset_intervention))
}
}
}
}
return(decline_rates)
}
#intersection routine for lists with catches on null lists
list_intersect <- function(list_a, list_b){
list_match = list()
vec_a <- unlist(list_a)
vec_b <- unlist(list_b)
if ( (length(vec_a) == 0) || (length(vec_b) == 0)){
return(list_match)
}
match_ind <- which(vec_a %in% vec_b)
match_val <- vec_a[match_ind]
list_match$match_ind = match_ind
list_match$match_val = match_val
return(list_match)
}
# Determine subset of object (eg development or offset group) containing a nested set of components with the same structure
select_pool_subset <- function(pool_object, subset_pool){
object_subset = lapply(seq_along(pool_object), function(i) pool_object[[i]][subset_pool])
names(object_subset) <- names(pool_object)
return(object_subset)
}
#function to work out vector of time intervals used in gains calculations
generate_time_horizons <- function(project_type, yr, offset_yrs, time_horizon, parcel_count){
if (project_type == 'current'){
# work out time intervals from time of offset to current year
time_horizons = rep(yr, parcel_count) - offset_yrs
} else if (project_type == 'future'){
# work out time intervals from current year to projected year defined by time_horizon
time_horizons = rep(time_horizon, parcel_count)
} else if (project_type == 'zeros'){
#set up dummy times
time_horizons = rep(0, parcel_count)
}
return(time_horizons)
}
assess_current_pool <- function(pool_object, pool_type, calc_type, cfacs_flag, adjust_cfacs_flag, action_type,
include_potential_developments, include_potential_offsets, include_unregulated_loss,
time_horizon_type, current_simulation_params, decline_rates_initial, time_horizon, yr){
current_pool = unlist(pool_object$site_indexes)
parcel_count = length(current_pool)
offset_yrs = unlist(pool_object$offset_yrs)
feature_num = length(current_simulation_params$features_to_use_in_simulation)
current_condition_vals = lapply(seq_along(pool_object$parcel_sums_at_offset),
function(i) pool_object$parcel_sums_at_offset[[i]][current_simulation_params$features_to_use_in_offset_calc])
if (calc_type == 'current_condition') {
projected_vals = current_condition_vals
cfac_vals = list_of_zeros(length(pool_object$parcel_sums_at_offset), length(current_simulation_params$features_to_use_in_offset_calc))
} else {
time_horizons <- generate_time_horizons(project_type = 'future', yr, offset_yrs, time_horizon, parcel_count)
if (cfacs_flag == TRUE){
cfacs_object = calc_cfacs(parcel_feature_layers = pool_object$parcel_feature_layers,
parcel_num_remaining = pool_object$parcel_num_remaining,
current_simulation_params,
decline_rates_initial[current_pool],
time_horizons,
offset_yrs,
include_potential_developments,
include_potential_offsets,
include_unregulated_loss,
adjust_cfacs_flag,
features_to_project = current_simulation_params$features_to_use_in_offset_calc)
cfac_vals = nested_list_tail(cfacs_object$cfacs_to_use)
} else if (calc_type == 'restoration_gains'){
cfac_vals = current_condition_vals
} else {
cfac_vals = list_of_zeros(length(pool_object$parcel_sums_at_offset), length(current_simulation_params$features_to_use_in_offset_calc))
}
if (pool_type == 'offsets') {
if (action_type == 'maintain'){
current_decline_rates = simulate_decline_rates(length(pool_object$parcel_feature_layers),
sample_decline_rate = FALSE,
mean_decline_rates = rep(list(1), feature_num),
decline_rate_std = vector())
} else if (action_type == 'restore'){
current_decline_rates = simulate_decline_rates(length(pool_object$parcel_feature_layers),
sample_decline_rate = FALSE,
mean_decline_rates = rep(list(current_simulation_params$restoration_rate), feature_num),
decline_rate_std = rep(list(current_simulation_params$restoration_rate_std), feature_num))
} else if (action_type == 'protect'){
current_decline_rates = decline_rates_initial[current_pool]
}
projected_feature_layers = project_sites(current_parcel_feature_layers = pool_object$parcel_feature_layers,
current_decline_rates,
time_horizons,
current_simulation_params,
features_to_project = current_simulation_params$features_to_use_in_offset_calc,
time_fill = FALSE)
#sum_sites
projected_vals = sum_sites(projected_feature_layers)
} else if (pool_type == 'devs') {
projected_vals = cfac_vals
cfac_vals = list_of_zeros(length(pool_object$parcel_sums_at_offset), length(current_simulation_params$features_to_use_in_offset_calc))
}
}
pool_object$parcel_vals_used = mapply('-', projected_vals, cfac_vals, SIMPLIFY = FALSE)
return(pool_object)
}
# function to append current object characteristics to group
append_current_group <- function(object_to_append, current_object, append_routine){
if (length(object_to_append) == 0){
#append null object with same fields to maintain order
object_to_append = vector('list', length(current_object))
names(object_to_append) = names(current_object)
}
appended_object <- append_current_object(object_to_append,
current_object,
append_type = 'as_list',
inds_to_append = seq_along(object_to_append))
# if (append_routine == 'banked_offset'){
# appended_object <- append_current_object(object_to_append,
# current_object,
# append_type = 'as_list',
# inds_to_append = seq_along(object_to_append)) #record current offset parcels in offsets object containing all offsets info
# } else {
#record current offset parcels in offsets object containing all offset characteristics
# needs a two step process as site_indexes are nested and the others are not
# appended_object <- append_current_object(object_to_append,
# current_object,
# append_type = 'as_list',
# inds_to_append = which(!(names(current_object) == 'site_indexes')))
# appended_object <- append_current_object(appended_object,
# current_object,
# append_type = 'as_group',
# inds_to_append = which(names(current_object) == 'site_indexes'))
# }
return(appended_object)
}
append_current_object <- function(parcel_set_object, current_parcel_set_object, append_type, inds_to_append){
if (append_type == 'as_group'){
#routine to append by characteristic for nested object
parcel_set_object[inds_to_append] <- lapply(inds_to_append, function(i) append(parcel_set_object[[i]], list(current_parcel_set_object[[i]])))
} else if (append_type == 'as_list'){
#routine to append by characteristic for non nested object
parcel_set_object[inds_to_append] <- lapply(inds_to_append, function(i) append(parcel_set_object[[i]], current_parcel_set_object[[i]]))
}
#set names of group to current object name
names(parcel_set_object) = names(current_parcel_set_object)
return(parcel_set_object)
}
destroy_development_sites <- function(current_feature_layers, decline_rates, feature_num){
parcel_num = length(current_feature_layers)
for (parcel_ind in seq_len(parcel_num)){
for (feature_ind in seq_len(feature_num)){
#select all development sites identified by decline rates = 0
current_decline_rate = decline_rates[[parcel_ind]][[feature_ind]]
if (current_decline_rate == 0){
#set all features of development sites to zero
current_feature_layers[[parcel_ind]][[feature_ind]] = array(0, length(current_feature_layers[[parcel_ind]][[feature_ind]]))
}
}
}
return(current_feature_layers)
}
calc_cfacs <- function(parcel_feature_layers, parcel_num_remaining, current_simulation_params,
decline_rates, time_horizons, offset_yrs, include_potential_developments, include_potential_offsets, include_unregulated_loss,
adjust_cfacs_flag, features_to_project){
cfacs_object = list()
if (length(decline_rates) != length(parcel_feature_layers)){
flog.error('calc cfacs length error')
}
cfacs_object$cfacs = project_sites(parcel_feature_layers,
decline_rates,
time_horizons,
current_simulation_params,
features_to_project,
time_fill = TRUE)
if (adjust_cfacs_flag == TRUE){
cfacs_object$adjusted_cfacs = adjust_cfacs(cfacs_object$cfacs,
include_potential_developments,
include_potential_offsets,
include_unregulated_loss,
current_simulation_params,
parcel_num_remaining,
decline_rates,
time_horizons,
offset_yrs)
cfacs_object$cfacs_to_use = sum_trajectories(cfacs_object$adjusted_cfacs)
} else{
cfacs_object$cfacs_to_use = sum_trajectories(cfacs_object$cfacs)
}
return(cfacs_object)
}
adjust_cfacs <- function(current_cfacs, include_potential_developments,include_potential_offsets, include_unregulated_loss,
current_simulation_params, parcel_num_remaining, decline_rates, time_horizons, offset_yrs){
time_horizons = unlist(time_horizons)
parcel_num = length(current_cfacs)
counter_weights = rep(list(1), parcel_num)
if (include_unregulated_loss == TRUE){
unregulated_loss_weights <- generate_weights(include_unregulated_loss,
calc_type = 'unregulated_loss',
current_simulation_params$max_offset_parcel_num,
current_simulation_params$intervention_vec,
offset_yrs,
time_horizons,
parcel_num,
parcel_num_remaining,
current_simulation_params$time_steps,
current_simulation_params$unregulated_loss_prob)
counter_weights <- lapply(seq_len(parcel_num), function(i) counter_weights[[i]] - unregulated_loss_weights$weights[[i]])
}
if (include_potential_developments == TRUE){
dev_probability_list <- generate_weights(include_potential_developments,
calc_type = 'development',
current_simulation_params$max_offset_parcel_num,
current_simulation_params$intervention_vec,
offset_yrs,
time_horizons,
parcel_num,
parcel_num_remaining,
current_simulation_params$time_steps,
current_simulation_params$unregulated_loss_prob)
counter_weights <- lapply(seq_len(parcel_num), function(i) counter_weights[[i]] - dev_probability_list$weights[[i]])
}
if (include_potential_offsets == TRUE){
current_offset_probability_list <- generate_weights(include_potential_offsets,
calc_type = 'offset',
current_simulation_params$max_offset_parcel_num,
current_simulation_params$intervention_vec,
offset_yrs,
time_horizons,
parcel_num,
parcel_num_remaining,
current_simulation_params$time_steps,
current_simulation_params$unregulated_loss_prob)
counter_weights <- lapply(seq_len(parcel_num), function(i) counter_weights[[i]] - current_offset_probability_list$weights[[i]])
}
inds_to_accept = lapply(seq_along(counter_weights), function(i) counter_weights[[i]] >= 0)
counter_weights <- remove_neg_probs(counter_weights, inds_to_accept)
adjusted_cfacs = lapply(seq_along(current_cfacs), function(i) lapply(seq_along(current_cfacs[[i]]),
function(j) current_cfacs[[i]][[j]]*matrix(rep(counter_weights[[i]], dim(current_cfacs[[i]][[j]])[2]),
nrow = dim(current_cfacs[[i]][[j]])[1], byrow = FALSE)))
if (include_potential_offsets == TRUE){
offset_intervention_probs <- remove_neg_probs(current_offset_probability_list$weighted_probs, inds_to_accept)
offset_projections <- calc_offset_projections(current_simulation_params$projection_type,
current_cfacs,
offset_intervention_probs,
current_simulation_params$restoration_rate,
time_horizons,
current_simulation_params$feature_num,
current_simulation_params$min_eco_val,
current_simulation_params$max_eco_val)
summed_offset_projections <- sum_offset_projs(offset_projections,
offset_intervention_probs,
current_simulation_params$feature_num,
time_horizons)
adjusted_cfacs = sum_clearing_offsets(adjusted_cfacs, summed_offset_projections, current_simulation_params$feature_num)
}
return(adjusted_cfacs)
}
remove_neg_probs <- function(weight_list, inds_to_accept){
weight_list <- lapply(seq_along(weight_list), function(i) weight_list[[i]]*inds_to_accept[[i]])
return(weight_list)
}
generate_weights <- function(perform_weight, calc_type, offset_intervention_scale, intervention_vec, offset_yrs, time_horizons,
parcel_num, parcel_num_remaining, time_steps, unregulated_loss_prob){
if (perform_weight == TRUE){
if (calc_type == 'unregulated_loss'){
weighted_probs <- lapply(seq_len(parcel_num), function(i) rep(unregulated_loss_prob, (time_horizons[i] + 1))) #runif(n = (time_horizon + 1), min = 0, max = current_simulation_params$unregulated_loss_prob)
} else {
weighted_probs <- find_intervention_probability(intervention_vec,
offset_yrs,
calc_type,
offset_intervention_scale,
time_horizons,
parcel_num,
parcel_num_remaining,
time_steps)
}
} else {
weighted_probs <- lapply(seq_len(parcel_num), function(i) rep(0, (time_horizons[i] + 1)))
}
weights <- lapply(weighted_probs, cumsum)
weighted_object = list()
weighted_object$weighted_probs = weighted_probs
weighted_object$weights = weights
return(weighted_object)
}
calc_offset_projections <- function(projection_type,current_cfacs, offset_probs, restoration_rate, time_horizons, feature_num, min_eco_val, max_eco_val){
parcel_num = length(current_cfacs)
offset_projections = vector('list', parcel_num)
for (parcel_ind in seq_len(parcel_num)){
time_horizon = time_horizons[parcel_ind] + 1
current_offset_probs = offset_probs[[parcel_ind]]
current_offset_projections = generate_nested_list(outer_dim = feature_num, inner_dim = time_horizon)
for (feature_ind in seq_len(feature_num)){
current_cfac = current_cfacs[[parcel_ind]][[feature_ind]]
for (proj_yr in seq_len(time_horizon)){
current_offset_projections[[feature_ind]][[proj_yr]] = array(0, dim(current_cfac))
if (current_offset_probs[proj_yr] > 0){
current_offset_proj = project_feature_layers(projection_type,
current_cfac[proj_yr, ],
min_eco_val,
max_eco_val,
restoration_rate,
(time_horizon - proj_yr),
time_fill = TRUE)
current_offset_projections[[feature_ind]][[proj_yr]][proj_yr:time_horizon, ] = current_offset_proj
}
}
}
offset_projections[[parcel_ind]] = current_offset_projections
}
return(offset_projections)
}
sum_offset_projs <- function(offset_projections, offset_probs, feature_num, time_horizons){
parcel_num = length(offset_projections)
summed_offset_projections = vector('list', parcel_num)
for (parcel_ind in seq_len(parcel_num)){
summed_offset_projections[[parcel_ind]] = vector('list', feature_num)
current_offset_prob = offset_probs[[parcel_ind]]
current_offset_prob <- current_offset_prob*(current_offset_prob > 0)
current_offset_proj = offset_projections[[parcel_ind]]
for (feature_ind in seq_len(feature_num)){
current_offset_projections <- current_offset_proj[[feature_ind]]
current_offset_projections <- lapply(seq_along(current_offset_projections), function(i) current_offset_projections[[i]]*current_offset_prob[i])
summed_offset_projections[[parcel_ind]][[feature_ind]] = Reduce('+', current_offset_projections)
}
}
return(summed_offset_projections)
}
sum_clearing_offsets <- function(cfacs_include_clearing, summed_offset_projections, feature_num){
parcel_num = length(cfacs_include_clearing)
cfacs_include_clearing_offsets = vector('list', parcel_num)
for (parcel_ind in 1:parcel_num){
cfacs_include_clearing_offsets[[parcel_ind]] = vector('list', feature_num)
}
for (parcel_ind in seq_len(parcel_num)){
if (length(summed_offset_projections[[parcel_ind]]) > 0 ){
for (feature_ind in seq_len(feature_num)){
cfacs_include_clearing_offsets[[parcel_ind]][[feature_ind]] = summed_offset_projections[[parcel_ind]][[feature_ind]] + cfacs_include_clearing[[parcel_ind]][[feature_ind]]
}
}
}
return(cfacs_include_clearing_offsets)
}
sum_cols_multi <- function(arr_to_use){
if (length(dim(arr_to_use)) <= 1){
arr_out = t(t(arr_to_use))
} else if (length(dim(arr_to_use)) == 2){
arr_out = apply(arr_to_use, MARGIN = 1, sum)
arr_out = t(t(arr_out))
} else if (length(dim(arr_to_use)) == 3){
arr_out = apply(arr_to_use, MARGIN = c(1, 3), sum)
dim(arr_out) = c(dim(arr_out), 1)
arr_out = aperm(arr_out, c(1, 3, 2))
}
return(arr_out)
}
# sum through features for each site to yield single dimensional site state vector for each feature
sum_trajectories <- function(traj_list, features_to_use_in_offset_calc){
parcel_traj_list = lapply(seq_along(traj_list), function(i) lapply(seq_along(traj_list[[i]]),
function(j) sum_cols_multi(traj_list[[i]][[j]])))
return(parcel_traj_list)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.