reflect when hitting bounds

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,

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

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

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

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))
+}

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)