diff --git a/R/GetAlgoParams.R b/R/GetAlgoParams.R index b999e63..bed2b4a 100644 --- a/R/GetAlgoParams.R +++ b/R/GetAlgoParams.R @@ -20,6 +20,7 @@ #' @param stop_tol A convergence metric must be less than value to be labeled as converged. The default is 1e-4. #' @param lower A numeric scalar or n_params-dimensional vector specifying lower bounds for each parameter. Default is -Inf (no lower bound). #' @param upper A numeric scalar or n_params-dimensional vector specifying upper bounds for each parameter. Default is Inf (no upper bound). +#' @param bounds_type A string specifying how parameter bounds are enforced on proposals. 'clip' truncates proposals to [lower, upper]. 'reflect' mirrors proposals back across the boundary (with clip fallback for very large steps). Default is 'reflect'. #' @param converge_crit A string denoting the convergence metric used, valid metrics are 'stdev' (standard deviation of population weight in the last stop_check iterations) and 'percent' (percent improvement in median particle weight in the last stop_check iterations). 'stdev' is the default. #' @return A list of control parameters for the optim_SQGDE function. #' @export @@ -31,6 +32,7 @@ GetAlgoParams = function(n_params, init_center = 0, lower = -Inf, upper = Inf, + bounds_type = 'reflect', n_cores_use = 1, step_size = NULL, jitter_size = 1e-6, @@ -121,6 +123,13 @@ GetAlgoParams = function(n_params, stop('ERROR: lower must be strictly less than upper for all parameters') } + # bounds_type + validBoundsType = c('clip', 'reflect') + if(is.null(bounds_type)) bounds_type = 'reflect' + if(!bounds_type %in% validBoundsType){ + stop('ERROR: bounds_type must be "clip" or "reflect"') + } + # n_cores_use ### assign NULL value default if(is.null(n_cores_use)){ @@ -292,6 +301,7 @@ GetAlgoParams = function(n_params, 'init_center' = init_center, 'lower' = lower, 'upper' = upper, + 'bounds_type' = bounds_type, 'n_cores_use' = n_cores_use, 'step_size' = step_size, 'crossover_rate' = crossover_rate, diff --git a/R/SQG_DE_bin_1_best.R b/R/SQG_DE_bin_1_best.R index 5b595f3..3bfdb46 100644 --- a/R/SQG_DE_bin_1_best.R +++ b/R/SQG_DE_bin_1_best.R @@ -21,6 +21,7 @@ SQG_DE_bin_1_best=function(pmem_index, crossover_rate = 1, lower = -Inf, upper = Inf, + bounds_type = 'reflect', n_diff, ... ){ # get statistics about particle @@ -74,7 +75,7 @@ SQG_DE_bin_1_best=function(pmem_index, step_size*psi*(grad_approx) + # move in the direction against the gradient stats::runif(len_param_use,-jitter_size,jitter_size) # a little noise } - params_use = pmax(lower, pmin(upper, params_use)) + params_use = apply_bounds(params_use, lower, upper, bounds_type) params_use = matrix(params_use,1,len_param_use) weight_proposal = NA diff --git a/R/SQG_DE_bin_1_curr.R b/R/SQG_DE_bin_1_curr.R index 10b7210..5ac6bf9 100644 --- a/R/SQG_DE_bin_1_curr.R +++ b/R/SQG_DE_bin_1_curr.R @@ -20,6 +20,7 @@ SQG_DE_bin_1_curr=function(pmem_index, crossover_rate = 1, lower = -Inf, upper = Inf, + bounds_type = 'reflect', n_diff, ... ){ # sample parents @@ -73,7 +74,7 @@ SQG_DE_bin_1_curr=function(pmem_index, step_size*psi*(grad_approx) + # move in the direction against the gradient stats::runif(len_param_use, -jitter_size, jitter_size) # a little noise } - params_use = pmax(lower, pmin(upper, params_use)) + params_use = apply_bounds(params_use, lower, upper, bounds_type) params_use = matrix(params_use,1,len_param_use) weight_proposal = NA diff --git a/R/SQG_DE_bin_1_rand.R b/R/SQG_DE_bin_1_rand.R index 2a700e4..99da953 100644 --- a/R/SQG_DE_bin_1_rand.R +++ b/R/SQG_DE_bin_1_rand.R @@ -21,6 +21,7 @@ SQG_DE_bin_1_rand=function(pmem_index, crossover_rate=1, lower=-Inf, upper=Inf, + bounds_type='reflect', n_diff, ... ){ @@ -74,7 +75,7 @@ SQG_DE_bin_1_rand=function(pmem_index, step_size*psi*(grad_approx) + # move in the direction against the gradient stats::runif(len_param_use,-jitter_size,jitter_size) # a little noise } - params_use = pmax(lower, pmin(upper, params_use)) + params_use = apply_bounds(params_use, lower, upper, bounds_type) params_use = matrix(params_use,1,len_param_use) weight_proposal = NA diff --git a/R/apply_bounds.R b/R/apply_bounds.R new file mode 100644 index 0000000..cdfaa17 --- /dev/null +++ b/R/apply_bounds.R @@ -0,0 +1,12 @@ +# Internal helper: enforce box constraints on a parameter vector. +# "clip" — truncate to [lower, upper] +# "reflect" — reflect off each boundary once, then clip (handles large steps) +apply_bounds <- function(x, lower, upper, type) { + if (type == "reflect") { + below <- is.finite(lower) & x < lower + x[below] <- 2 * lower[below] - x[below] + above <- is.finite(upper) & x > upper + x[above] <- 2 * upper[above] - x[above] + } + pmax(lower, pmin(upper, x)) +} diff --git a/R/optim_SQGDE.R b/R/optim_SQGDE.R index 195abf2..2ff385a 100644 --- a/R/optim_SQGDE.R +++ b/R/optim_SQGDE.R @@ -143,6 +143,7 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), ...){ crossover_rate = control_params$crossover_rate, lower = control_params$lower, upper = control_params$upper, + bounds_type = control_params$bounds_type, n_diff = control_params$n_diff, ...)), control_params$n_particles, control_params$n_params+1, byrow=TRUE) @@ -158,6 +159,7 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), ...){ crossover_rate = control_params$crossover_rate, lower = control_params$lower, upper = control_params$upper, + bounds_type = control_params$bounds_type, n_diff = control_params$n_diff, ...)), control_params$n_particles, control_params$n_params+1, byrow=TRUE) diff --git a/man/GetAlgoParams.Rd b/man/GetAlgoParams.Rd index 6fa70b9..7e08409 100644 --- a/man/GetAlgoParams.Rd +++ b/man/GetAlgoParams.Rd @@ -13,6 +13,7 @@ GetAlgoParams( init_center = 0, lower = -Inf, upper = Inf, + bounds_type = "reflect", n_cores_use = 1, step_size = NULL, jitter_size = 1e-06, @@ -45,6 +46,8 @@ GetAlgoParams( \item{upper}{A numeric scalar or n_params-dimensional vector specifying upper bounds for each parameter. Default is Inf (no upper bound).} +\item{bounds_type}{A string specifying how parameter bounds are enforced on proposals. 'clip' truncates proposals to [lower, upper]. 'reflect' mirrors proposals back across the boundary (with clip fallback for very large steps). Default is 'reflect'.} + \item{n_cores_use}{An integer specifying the number of cores used when using parallelization. The default value is 1.} \item{step_size}{A positive scalar, jump size or "F" in the DE crossover step notation. The default value is 2.38/sqrt(2*n_params).}