Unified crossover step

R/GetAlgoParams.R

@@ -4,13 +4,15 @@
 #' @param n_particles The number of particles (population size), 3*n_params is the default value.
 #' @param n_iter The number of iterations to run the algorithm, 1000 is default.
 #' @param n_diff The number of mutually exclusive vector pairs to stochastically approximate the gradient.
-#' @param crossover_rate A numeric scalar on the interval (0,1]. Determines the probability a parameter on a chain is updated on a given crossover step, sampled from a Bernoulli distribution. The default value is 1.
+#' @param crossover_rate A numeric scalar on the interval [0,1]. Determines the probability a parameter on a chain is updated on a given crossover step, sampled from a Bernoulli distribution. When 0, exactly one randomly chosen parameter is updated per iteration. The default value is 1.
 #' @param init_sd A positive scalar or n_params-dimensional numeric vector, determines the standard deviation of the Gaussian initialization distribution. The default value is 0.01.
 #' @param init_center A scalar or n_params-dimensional numeric vector, determines the mean of the Gaussian initialization distribution. The default value is 0.
 #' @param n_cores_use An integer specifying the number of cores used when using parallelization. The default value is 1.
 #' @param 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).
 #' @param jitter_size A positive scalar that determines the jitter (noise) size. Noise is added during adaption step from Uniform(-jitter_size,jitter_size) distribution. 1e-6 is the default value. Set to 0 to turn off jitter.
 #' @param parallel_type A string specifying parallelization type. 'none','FORK', or 'PSOCK' are valid values. 'none' is default value. 'FORK' does not work with Windows OS.
+#' @param recovery_path A character scalar giving a file path where partial results are written via \code{saveRDS} periodically. Allows recovery of progress if the run is interrupted or crashes. Load with \code{readRDS(recovery_path)}. The saved object matches the structure of the return value of \code{optim_SQGDE} but excludes trace arrays. Set to \code{NULL} (default) to disable.
+#' @param recovery_freq A positive integer controlling how often the recovery file is written. The file is saved every \code{recovery_freq} iterations. Default is 1 (every iteration). Ignored when \code{recovery_path} is \code{NULL}.
 #' @param return_trace A boolean, if true, the function returns particle trajectories. This is helpful for assessing convergence or debugging model code. The trace will be an iteration/thin $x$ n_particles $x$ n_params array containing parameter values and an iteration/thin $x$ n_particles array containing particle weights.
 #' @param thin A positive integer. Only every 'thin'-th iteration will be stored in memory. The default value is 1. Increasing thin will reduce the memory required when running the algorithim for longer.
 #' @param purify A positive integer. On every 'purify'-th iteration the particle weights are recomputed. This is useful if the objective function is stochastic/noisy. If the objective function is deterministic, this computation is redundant. Purify is set to Inf by default, disabling it.
@@ -21,6 +23,7 @@
 #' @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 trace_print_freq A positive integer controlling how often iteration progress is printed. A message is printed every \code{trace_print_freq} iterations. Set to \code{Inf} to suppress progress messages entirely. Default is 100.
 #' @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
@@ -38,14 +41,17 @@ GetAlgoParams = function(n_params,
                          jitter_size = 1e-6,
                          crossover_rate = 1,
                          parallel_type = 'none',
+                         recovery_path = NULL,
+                         recovery_freq = 1,
                          return_trace = FALSE,
                          thin = 1,
                          purify = Inf,
                          adapt_scheme = NULL,
                          give_up_init = 100,
                          stop_check = 10,
                          stop_tol = 1e-4,
-                         converge_crit = 'stdev'){
+                         converge_crit = 'stdev',
+                         trace_print_freq = 100){
   # n_params
   ### catch errors
   n_params = as.integer(n_params)
@@ -181,8 +187,8 @@ GetAlgoParams = function(n_params,
   crossover_rate = as.numeric(crossover_rate)
   if(any(!is.finite(crossover_rate))){
     stop('ERROR: crossover_rate is not finite')
-  } else if(any(crossover_rate>1) | any(crossover_rate<= 0) | length(crossover_rate)>1){
-    stop('ERROR: crossover_rate must be a numeric scalar on the interval (0,1]')
+  } else if(any(crossover_rate>1) | any(crossover_rate< 0) | length(crossover_rate)>1){
+    stop('ERROR: crossover_rate must be a numeric scalar on the interval [0,1]')
   } else if(is.complex(crossover_rate)){
     stop('ERROR: crossover_rate cannot be complex')
   }
@@ -198,6 +204,22 @@ GetAlgoParams = function(n_params,
     stop(paste('ERROR: invalid parallel_type.'))
   }
 
+  # recovery_path
+  if(!is.null(recovery_path)){
+    if(!is.character(recovery_path) | length(recovery_path) != 1){
+      stop('ERROR: recovery_path must be NULL or a character scalar file path')
+    }
+    if(!grepl('\\.rds$', recovery_path, ignore.case = TRUE)){
+      stop('ERROR: recovery_path must end in .rds')
+    }
+  }
+
+  # recovery_freq
+  recovery_freq = as.integer(recovery_freq)
+  if(length(recovery_freq) != 1 || !is.finite(recovery_freq) || recovery_freq < 1){
+    stop('ERROR: recovery_freq must be a positive integer scalar')
+  }
+
   #converge_crit
   validConType = c('stdev','percent')
   ### assign NULL value default
@@ -294,6 +316,18 @@ GetAlgoParams = function(n_params,
     stop('ERROR: stop_tol must be a scalar positive')
   }
 
+  # trace_print_freq
+  if(is.null(trace_print_freq)) trace_print_freq = 100
+  if(!(length(trace_print_freq) == 1)){
+    stop('ERROR: trace_print_freq must be a scalar')
+  }
+  if(is.finite(trace_print_freq)){
+    trace_print_freq = as.integer(trace_print_freq)
+    if(trace_print_freq < 1){
+      stop('ERROR: trace_print_freq must be a positive integer or Inf')
+    }
+  }
+
   out = list('n_params' = n_params,
              'n_particles' = n_particles,
              'n_iter' = n_iter,
@@ -307,6 +341,8 @@ GetAlgoParams = function(n_params,
              'crossover_rate' = crossover_rate,
              'jitter_size' = jitter_size,
              'parallel_type' = parallel_type,
+             'recovery_path' = recovery_path,
+             'recovery_freq' = recovery_freq,
              'thin' = thin,
              'purify' = purify,
              'n_iters_per_particle' = n_iters_per_particle,
@@ -316,7 +352,8 @@ GetAlgoParams = function(n_params,
              'give_up_init'= give_up_init,
              'stop_tol' = stop_tol,
              'stop_check' = stop_check,
-             'converge_crit' = converge_crit)
+             'converge_crit' = converge_crit,
+             'trace_print_freq' = trace_print_freq)
 
   return(out)
 }

R/SQG_DE_bin_1.R

@@ -0,0 +1,108 @@
+# Internal helper: stochastic quasi-gradient approximation.
+# Returns psi (self-scaling factor) and grad_approx (approximate gradient vector).
+# Skips any pair where the parameter difference or weight difference is zero
+# to avoid division by zero.
+grad_approx_fn = function(param_indices, n_diff, current_params,
+                           parent_indices, current_weight) {
+  vec_diff_sum = grad_approx = numeric(length(param_indices))
+  for (d in 1:n_diff) {
+    vec_diff_temp = (current_params[parent_indices[d],     param_indices] -
+                     current_params[parent_indices[d+n_diff], param_indices])
+    vec_diff_sum  = vec_diff_sum + vec_diff_temp
+    weight_diff_temp = (current_weight[parent_indices[d]] -
+                        current_weight[parent_indices[d+n_diff]])
+    if (!(all(vec_diff_temp == 0) | weight_diff_temp == 0)) {
+      grad_approx = grad_approx +
+        vec_diff_temp * (weight_diff_temp / sqrt(sum(vec_diff_temp^2)))
+    }
+  }
+  psi_num = sqrt(sum(vec_diff_sum^2)) * (1/n_diff)
+  psi_den = sqrt(sum(grad_approx^2))
+  list(psi = psi_num / psi_den, grad_approx = grad_approx)
+}
+
+
+#' SQG_DE_bin_1
+#'
+#' @param pmem_index Index of population member being updated.
+#' @param current_params Current parameter matrix (n_particles x n_params).
+#' @param current_weight Current weight vector (length n_particles).
+#' @param objFun Objective function to minimize.
+#' @param scheme Adaptation scheme: 'rand', 'best', or 'current'.
+#' @param step_size Scaling factor for gradient step.
+#' @param jitter_size Half-width of uniform jitter added to proposals.
+#' @param n_particles Number of particles.
+#' @param crossover_rate Bernoulli probability of updating each parameter.
+#' @param lower Lower bound vector (length n_params).
+#' @param upper Upper bound vector (length n_params).
+#' @param bounds_type Bounds enforcement: 'reflect' or 'clip'.
+#' @param n_diff Number of particle pairs used for gradient approximation.
+#' @param ... Additional arguments passed to objFun.
+#' @noRd
+
+SQG_DE_bin_1 = function(pmem_index,
+                         current_params,
+                         current_weight,
+                         objFun,
+                         scheme       = 'rand',
+                         step_size    = .8,
+                         jitter_size  = 1e-6,
+                         n_particles,
+                         crossover_rate = 1,
+                         lower        = -Inf,
+                         upper        =  Inf,
+                         bounds_type  = 'reflect',
+                         n_diff, ...) {
+
+  weight_use       = current_weight[pmem_index]
+  params_use       = current_params[pmem_index, ]
+  len_param_use    = length(params_use)
+  best_pmem_index  = which.min(current_weight)
+
+  # scheme determines: parent pool, sample size, and base position for proposal
+  if (scheme == 'best') {
+    parent_indices = sample(c(1:n_particles)[-best_pmem_index],
+                            size = 2*n_diff, replace = FALSE)
+    base_index = best_pmem_index
+  } else if (scheme == 'current') {
+    parent_indices = sample(c(1:n_particles)[-pmem_index],
+                            size = 2*n_diff, replace = FALSE)
+    base_index = pmem_index
+  } else {
+    parent_indices = sample(1:n_particles, size = 2*n_diff+1, replace = FALSE)
+    base_index = parent_indices[2*n_diff+1]
+  }
+
+  # crossover: select which parameters to update
+  param_idices_bool = stats::rbinom(len_param_use, prob = crossover_rate, size = 1)
+  if (all(param_idices_bool == 0)) {
+    param_idices_bool[sample(x = 1:len_param_use, size = 1)] = 1
+  }
+  param_indices = seq(1, len_param_use, by = 1)[as.logical(param_idices_bool)]
+
+  # approximate gradient and self-scaling factor
+  ga           = grad_approx_fn(param_indices, n_diff, current_params,
+                                parent_indices, current_weight)
+  grad_approx  = ga$grad_approx
+  psi          = ga$psi
+
+  if (all(is.finite(grad_approx)) & is.finite(psi) & (psi > 0)) {
+    params_use[param_indices] = current_params[base_index, param_indices] -
+      step_size * psi * grad_approx +
+      stats::runif(length(param_indices), -jitter_size, jitter_size)
+  }
+
+  params_use = apply_bounds(params_use, lower, upper, bounds_type)
+  params_use = matrix(params_use, 1, len_param_use)
+
+  weight_proposal = NA
+  if (all(is.finite(params_use))) weight_proposal = objFun(params_use, ...)
+  if (is.na(weight_proposal)) weight_proposal = Inf
+
+  if (weight_proposal < weight_use) {
+    current_params[pmem_index, ] = params_use
+    current_weight[pmem_index]   = weight_proposal
+  }
+
+  return(c(current_weight[pmem_index], current_params[pmem_index, ]))
+}