From ab5c1dee3d39e5c1e64784cbf8b16fce7790f2d0 Mon Sep 17 00:00:00 2001
From: Brendan Matthew Galdo <brendan.m.galdo@gmail.com>
Date: Sun, 10 May 2026 21:05:14 -0400
Subject: [PATCH 1/5] create unified function for adapation sequences

---
 R/SQG_DE_bin_1.R      | 108 ++++++++++++++++++++++++++++++++++++++++++
 R/SQG_DE_bin_1_best.R |  94 ------------------------------------
 R/SQG_DE_bin_1_curr.R |  93 ------------------------------------
 R/SQG_DE_bin_1_rand.R |  94 ------------------------------------
 R/optim_SQGDE.R       |  29 ++++--------
 5 files changed, 118 insertions(+), 300 deletions(-)
 create mode 100644 R/SQG_DE_bin_1.R
 delete mode 100644 R/SQG_DE_bin_1_best.R
 delete mode 100644 R/SQG_DE_bin_1_curr.R
 delete mode 100644 R/SQG_DE_bin_1_rand.R

diff --git a/R/SQG_DE_bin_1.R b/R/SQG_DE_bin_1.R
new file mode 100644
index 0000000..b57ba97
--- /dev/null
+++ b/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(len_param_use, -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, ]))
+}
diff --git a/R/SQG_DE_bin_1_best.R b/R/SQG_DE_bin_1_best.R
deleted file mode 100644
index 3bfdb46..0000000
--- a/R/SQG_DE_bin_1_best.R
+++ /dev/null
@@ -1,94 +0,0 @@
-#' SQG_DE_bin_1_best
-#'
-#' @param pmem_index  Index of population (particle) member you are you are updating
-#' @param current_params Current parameter values for partcle (numeric vector)
-#' @param current_weight  weights for current population
-#' @param objFun  function we want to minimize
-#' @param n_particles number of particles
-#' @param crossover_rate i.i.d. bernouli probability for updating a parameter
-#' @param n_diff number of particles used to estimate the gradient
-#' @param ... additional arguments for objective function
-#' @noRd
-#' 
-
-SQG_DE_bin_1_best=function(pmem_index,
-                           current_params,
-                           current_weight,
-                           objFun,
-                           step_size = .8,
-                           jitter_size = 1e-6,
-                           n_particles,
-                           crossover_rate = 1,
-                           lower = -Inf,
-                           upper = Inf,
-                           bounds_type = 'reflect',
-                           n_diff, ... ){
-
-  # get statistics about particle
-  weight_use = current_weight[pmem_index]
-  best_pmem_index = which.min(current_weight)
-  params_use = current_params[pmem_index,]
-  len_param_use = length(params_use)
-  
-  # sample parent chains
-  parent_indices = sample(c(1:n_particles)[-best_pmem_index],size=2*n_diff,replace=FALSE)
-  
-  # use binomial to sample which params to update matching crossover rate frequency
-  param_idices_bool = stats::rbinom(len_param_use, prob = crossover_rate, size=1)
-
-  # if no params selected, randomly sample 1 parameter
-  if(all(param_idices_bool==0)){
-    param_idices_bool[sample(x=1:len_param_use,size=1)] = 1
-  }
-
-  # indices of parameters to be updated
-  param_indices = seq(1,len_param_use,by=1)[as.logical(param_idices_bool)]
-
-
-  # get rid of this for loop later
-  vec_diff_sum = grad_approx = numeric(length(param_indices))
-  for(d in 1:n_diff){
-    # difference in vector pairs
-    vec_diff_temp = (current_params[parent_indices[d],param_indices] -
-                     current_params[parent_indices[d+n_diff],param_indices])
-
-    # sum up all vector differences for normalization step
-    vec_diff_sum = vec_diff_sum + vec_diff_temp
-
-    # difference in function values
-    weight_diff_temp = (current_weight[parent_indices[d]] -
-                        current_weight[parent_indices[d+n_diff]])
-
-    # sum the approximate gradient vectors
-    grad_approx = grad_approx + vec_diff_temp*(weight_diff_temp/sqrt(sum(vec_diff_temp^2)))
-  }
-
-  # calculate normalization factor for algorithm self scaling
-  psi_num = sqrt(sum(vec_diff_sum^2))*(1/n_diff)
-  psi_den = sqrt(sum(grad_approx^2))
-  psi = psi_num/psi_den
-
-
-  if(all(is.finite(grad_approx)) & is.finite(psi) & (psi>0)){
-    # mate parents for proposal
-    params_use[param_indices] = current_params[best_pmem_index,param_indices] -
-      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 = apply_bounds(params_use, lower, upper, bounds_type)
-  params_use = matrix(params_use,1,len_param_use)
-
-  weight_proposal = NA
-  # get weight
-  if(all(is.finite(params_use)))weight_proposal = objFun(params_use,...)
-  if(is.na(weight_proposal))weight_proposal = Inf
-
-  #negative greedy acceptance rule
-  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,]))
-
-}
diff --git a/R/SQG_DE_bin_1_curr.R b/R/SQG_DE_bin_1_curr.R
deleted file mode 100644
index 5ac6bf9..0000000
--- a/R/SQG_DE_bin_1_curr.R
+++ /dev/null
@@ -1,93 +0,0 @@
-#' SQG_DE_bin_1_curr
-#'
-#' @param pmem_index  Index of population (particle) member you are you are updating
-#' @param current_params Current parameter values for partcle (numeric vector)
-#' @param current_weight  weights for current population
-#' @param objFun  function we want to minimize
-#' @param n_particles number of particles
-#' @param crossover_rate i.i.d. bernouli probability for updating a parameter
-#' @param n_diff number of particles used to estimate the gradient
-#' @param ... additional arguments for objective function
-#' @noRd
-#' 
-SQG_DE_bin_1_curr=function(pmem_index,
-                           current_params,
-                           current_weight,
-                           objFun,
-                           step_size = .8,
-                           jitter_size = 1e-6,
-                           n_particles,
-                           crossover_rate = 1,
-                           lower = -Inf,
-                           upper = Inf,
-                           bounds_type = 'reflect',
-                           n_diff, ... ){
-
-  # sample parents
-  parent_indices = sample(c(1:n_particles)[-pmem_index],size=2*n_diff,replace=F)
-
-  # get statistics about sparticle
-  weight_use = current_weight[pmem_index]
-  # best_pmem_index = which.max(current_weight)
-  params_use = current_params[pmem_index,]
-  len_param_use = length(params_use)
-
-  # use binomial to sample which params to update matching crossover rate frequency
-  param_idices_bool = stats::rbinom(len_param_use, prob = crossover_rate, size=1)
-
-  # if no pars selected, randomly sample 1 parameter
-  if(all(param_idices_bool==0)){
-    param_idices_bool[sample(x=1:len_param_use,size=1)] = 1
-  }
-
-  # indices of parameters to be updated
-  param_indices = seq(1,len_param_use,by=1)[as.logical(param_idices_bool)]
-
-
-  # get rid of this for loop later
-  vec_diff_sum = grad_approx = numeric(length(param_indices))
-  for(d in 1:n_diff){
-    # difference in vector pairs
-    vec_diff_temp = (current_params[parent_indices[d],param_indices] -
-                     current_params[parent_indices[d+n_diff],param_indices])
-
-    # sum up all vector differences for normalization step
-    vec_diff_sum = vec_diff_sum+vec_diff_temp
-
-    # difference in function values
-    weight_diff_temp = (current_weight[parent_indices[d]] -
-                        current_weight[parent_indices[d+n_diff]])
-
-    # sum the approximate gradient vectors
-    grad_approx = grad_approx + vec_diff_temp*(weight_diff_temp/sqrt(sum(vec_diff_temp^2)))
-  }
-
-  # calculate normalization factor for algorithm self scaling
-  psi_num = sqrt(sum(vec_diff_sum^2))*(1/n_diff)
-  psi_den = sqrt(sum(grad_approx^2))
-  psi=psi_num/psi_den
-
-
-  if(all(is.finite(grad_approx)) & is.finite(psi) & (psi>0)){
-    # mate parents for proposal
-    params_use[param_indices] = current_params[pmem_index, param_indices] -
-      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 = apply_bounds(params_use, lower, upper, bounds_type)
-  params_use = matrix(params_use,1,len_param_use)
-
-  weight_proposal = NA
-  # get weight
-  if(all(is.finite(params_use)))weight_proposal = objFun(params_use,...)
-  if(is.na(weight_proposal))weight_proposal = Inf
-
-  #negative greedy acceptance rule
-  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,]))
-
-}
diff --git a/R/SQG_DE_bin_1_rand.R b/R/SQG_DE_bin_1_rand.R
deleted file mode 100644
index 99da953..0000000
--- a/R/SQG_DE_bin_1_rand.R
+++ /dev/null
@@ -1,94 +0,0 @@
-#' SQG_DE_bin_1_rand
-#'
-#' @param pmem_index  Index of population (particle) member you are you are updating
-#' @param current_params Current parameter values for partcle (numeric vector)
-#' @param current_weight  weights for current population
-#' @param objFun  function we want to minimize
-#' @param n_particles number of particles
-#' @param crossover_rate i.i.d. bernouli probability for updating a parameter
-#' @param n_diff number of particles used to estimate the gradient
-#' @param ... additional arguments for objective function
-#' 
-#' @noRd
-#' 
-SQG_DE_bin_1_rand=function(pmem_index,
-                           current_params,
-                           current_weight,
-                           objFun,
-                           step_size=.8,
-                           jitter_size=1e-6,
-                           n_particles,
-                           crossover_rate=1,
-                           lower=-Inf,
-                           upper=Inf,
-                           bounds_type='reflect',
-                           n_diff, ... ){
-
-
-  # get current particle info
-  weight_use = current_weight[pmem_index]
-  params_use = current_params[pmem_index,]
-  len_param_use = length(params_use)
-  
-  # sample parent 
-  parent_indices = sample(c(1:n_particles),size=2*n_diff+1,replace=FALSE)
-  
-  # use binomial to sample which params to update matching crossover rate frequency
-  param_idices_bool = stats::rbinom(len_param_use, prob = crossover_rate, size=1)
-
-  # if no params selected, randomly sample 1 parameter
-  if(all(param_idices_bool==0)){
-    param_idices_bool[sample(x=1:len_param_use,size=1)] = 1
-  }
-
-  # indices of parameters to be updated
-  param_indices = seq(1,len_param_use,by=1)[as.logical(param_idices_bool)]
-
-
-  # get rid of this for loop later
-  vec_diff_sum = grad_approx = numeric(length(param_indices))
-  for(d in 1:n_diff){
-    # difference in vector pairs
-    vec_diff_temp=(current_params[parent_indices[d],param_indices] -
-                     current_params[parent_indices[d+n_diff],param_indices])
-
-    # sum up all vector differences for normalization step
-    vec_diff_sum = vec_diff_sum+vec_diff_temp
-
-    # difference in function values
-    weight_diff_temp = (current_weight[parent_indices[d]] -
-                        current_weight[parent_indices[d+n_diff]])
-
-    # sum the approximate gradient vectors
-    grad_approx = grad_approx+vec_diff_temp*(weight_diff_temp/sqrt(sum(vec_diff_temp^2)))
-  }
-
-  # calculate normalization factor for algorithm self scaling
-  psi_num = sqrt(sum(vec_diff_sum^2))*(1/n_diff)
-  psi_den = sqrt(sum(grad_approx^2))
-  psi = psi_num/psi_den
-
-
-  if(all(is.finite(grad_approx)) & is.finite(psi) & (psi>0)){
-    # mate parents for proposal
-    params_use[param_indices] = current_params[parent_indices[2*n_diff+1],param_indices] -
-      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 = apply_bounds(params_use, lower, upper, bounds_type)
-  params_use = matrix(params_use,1,len_param_use)
-
-  weight_proposal = NA
-  # get weight
-  if(all(is.finite(params_use)))weight_proposal = objFun(params_use,...)
-  if(is.na(weight_proposal))weight_proposal = Inf
-
-  #negative greedy acceptance rule
-  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,]))
-
-}
diff --git a/R/optim_SQGDE.R b/R/optim_SQGDE.R
index 7fd8243..9fdb9f5 100644
--- a/R/optim_SQGDE.R
+++ b/R/optim_SQGDE.R
@@ -122,17 +122,6 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), warm_start = NU
     message('population initialization complete  :)')
   }
 
-  # assign adaption scheme
-  if(control_params$adapt_scheme=='rand'){
-    AdaptSQGDE = SQG_DE_bin_1_rand
-  }
-  if(control_params$adapt_scheme=='best'){
-    AdaptSQGDE = SQG_DE_bin_1_best
-  }
-  if(control_params$adapt_scheme=='current'){
-    AdaptSQGDE = SQG_DE_bin_1_curr
-  }
-
   # cluster initialization
   if(!control_params$parallel_type=='none'){
 
@@ -155,10 +144,11 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), warm_start = NU
 
     if(control_params$parallel_type=='none'){
       # adapt particles using SQG DE sequentially
-      temp=matrix(unlist(lapply(1:control_params$n_particles, AdaptSQGDE,
-                                current_params = matrix(particles[iter_idx, , ], nrow = control_params$n_particles, ncol = control_params$n_params),   # current parameter values (numeric matrix)
-                                current_weight = weights[iter_idx, ],  # corresponding weights (numeric vector)
-                                objFun = ObjFun,  # objective function (returns scalar)
+      temp=matrix(unlist(lapply(1:control_params$n_particles, SQG_DE_bin_1,
+                                current_params = matrix(particles[iter_idx, , ], nrow = control_params$n_particles, ncol = control_params$n_params),
+                                current_weight = weights[iter_idx, ],
+                                objFun = ObjFun,
+                                scheme = control_params$adapt_scheme,
                                 step_size = control_params$step_size,
                                 jitter_size = control_params$jitter_size,
                                 n_particles = control_params$n_particles,
@@ -171,10 +161,11 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), warm_start = NU
                   control_params$n_params+1, byrow=TRUE)
     } else {
       # adapt particles using SQG DE in parallel
-      temp=matrix(unlist(parallel::parLapplyLB(cl_use, 1:control_params$n_particles, AdaptSQGDE,
-                                             current_params = matrix(particles[iter_idx, , ], nrow = control_params$n_particles, ncol = control_params$n_params),   # current parameter values (numeric matrix)
-                                             current_weight = weights[iter_idx, ],  # corresponding weight (numeric vector)
-                                             objFun = ObjFun,  # function we want to minimize (returns scalar)
+      temp=matrix(unlist(parallel::parLapplyLB(cl_use, 1:control_params$n_particles, SQG_DE_bin_1,
+                                             current_params = matrix(particles[iter_idx, , ], nrow = control_params$n_particles, ncol = control_params$n_params),
+                                             current_weight = weights[iter_idx, ],
+                                             objFun = ObjFun,
+                                             scheme = control_params$adapt_scheme,
                                              step_size = control_params$step_size,
                                              jitter_size = control_params$jitter_size,
                                              n_particles = control_params$n_particles,

From ba055efe389abf76dccd4b191a0da4ffe8550f63 Mon Sep 17 00:00:00 2001
From: Brendan Matthew Galdo <brendan.m.galdo@gmail.com>
Date: Sun, 10 May 2026 21:10:53 -0400
Subject: [PATCH 2/5] Create test-adapt_scheme.R

---
 tests/testthat/test-adapt_scheme.R | 76 ++++++++++++++++++++++++++++++
 1 file changed, 76 insertions(+)
 create mode 100644 tests/testthat/test-adapt_scheme.R

diff --git a/tests/testthat/test-adapt_scheme.R b/tests/testthat/test-adapt_scheme.R
new file mode 100644
index 0000000..bfbd38c
--- /dev/null
+++ b/tests/testthat/test-adapt_scheme.R
@@ -0,0 +1,76 @@
+library(graDiEnt)
+
+obj <- function(x) sum((x - c(1, -1, 2))^2)
+true_opt <- c(1, -1, 2)
+
+cp_base <- function(scheme) {
+  GetAlgoParams(
+    n_params    = 3,
+    n_iter      = 600,
+    n_particles = 12,
+    n_diff      = 2,
+    init_sd     = 1,
+    init_center = 0,
+    adapt_scheme = scheme
+  )
+}
+
+# ── each scheme converges ─────────────────────────────────────────────────────
+
+test_that("adapt_scheme 'rand' converges to known optimum", {
+  set.seed(42)
+  out <- suppressMessages(optim_SQGDE(obj, cp_base("rand")))
+
+  expect_true(all(is.finite(out$solution)))
+  expect_lt(max(abs(out$solution - true_opt)), 0.1)
+})
+
+test_that("adapt_scheme 'best' converges to known optimum", {
+  set.seed(42)
+  out <- suppressMessages(optim_SQGDE(obj, cp_base("best")))
+
+  expect_true(all(is.finite(out$solution)))
+  expect_lt(max(abs(out$solution - true_opt)), 0.1)
+})
+
+test_that("adapt_scheme 'current' converges to known optimum", {
+  set.seed(42)
+  out <- suppressMessages(optim_SQGDE(obj, cp_base("current")))
+
+  expect_true(all(is.finite(out$solution)))
+  expect_lt(max(abs(out$solution - true_opt)), 0.1)
+})
+
+# ── each scheme returns correct output structure ──────────────────────────────
+
+test_that("all schemes return valid list structure", {
+  for (scheme in c("rand", "best", "current")) {
+    set.seed(1)
+    out <- suppressMessages(
+      optim_SQGDE(obj, GetAlgoParams(n_params = 3, n_iter = 5,
+                                     n_particles = 12, n_diff = 2,
+                                     adapt_scheme = scheme))
+    )
+    expect_type(out, "list")
+    expect_true(all(is.finite(out$solution)))
+    expect_true(is.finite(out$weight))
+    expect_length(out$solution, 3)
+  }
+})
+
+# ── invalid scheme errors ─────────────────────────────────────────────────────
+
+test_that("GetAlgoParams rejects invalid adapt_scheme", {
+  expect_error(
+    GetAlgoParams(n_params = 2, adapt_scheme = "gradient"),
+    "invalid adaption scheme"
+  )
+  expect_error(
+    GetAlgoParams(n_params = 2, adapt_scheme = "random"),
+    "invalid adaption scheme"
+  )
+  expect_error(
+    GetAlgoParams(n_params = 2, adapt_scheme = ""),
+    "invalid adaption scheme"
+  )
+})

From b1844af575025ced3fd829b6610bdf5505eef268 Mon Sep 17 00:00:00 2001
From: Brendan Matthew Galdo <brendan.m.galdo@gmail.com>
Date: Sun, 10 May 2026 21:32:13 -0400
Subject: [PATCH 3/5] give user the ability to control test frequency

---
 R/GetAlgoParams.R                      |  25 ++++-
 R/SQG_DE_bin_1.R                       |   2 +-
 R/optim_SQGDE.R                        |   2 +-
 tests/testthat/test-bug_fixes.R        | 138 +++++++++++++++++++++++++
 tests/testthat/test-trace_print_freq.R |  75 ++++++++++++++
 5 files changed, 235 insertions(+), 7 deletions(-)
 create mode 100644 tests/testthat/test-bug_fixes.R
 create mode 100644 tests/testthat/test-trace_print_freq.R

diff --git a/R/GetAlgoParams.R b/R/GetAlgoParams.R
index bed2b4a..15ec509 100644
--- a/R/GetAlgoParams.R
+++ b/R/GetAlgoParams.R
@@ -4,7 +4,7 @@
 #' @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.
@@ -21,6 +21,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
@@ -45,7 +46,8 @@ GetAlgoParams = function(n_params,
                          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 +183,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')
   }
@@ -294,6 +296,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,
@@ -316,7 +330,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)
 }
diff --git a/R/SQG_DE_bin_1.R b/R/SQG_DE_bin_1.R
index b57ba97..346de08 100644
--- a/R/SQG_DE_bin_1.R
+++ b/R/SQG_DE_bin_1.R
@@ -89,7 +89,7 @@ SQG_DE_bin_1 = function(pmem_index,
   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(len_param_use, -jitter_size, jitter_size)
+      stats::runif(length(param_indices), -jitter_size, jitter_size)
   }
 
   params_use = apply_bounds(params_use, lower, upper, bounds_type)
diff --git a/R/optim_SQGDE.R b/R/optim_SQGDE.R
index 9fdb9f5..067537a 100644
--- a/R/optim_SQGDE.R
+++ b/R/optim_SQGDE.R
@@ -243,7 +243,7 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), warm_start = NU
 
 
 
-    if(iter%%100==0){
+    if(is.finite(control_params$trace_print_freq) && iter%%control_params$trace_print_freq==0){
       message(paste0('iter ', iter, '/', control_params$n_iter))
     }
     if(iter%%control_params$thin==0 & !(iter==control_params$n_iter)){
diff --git a/tests/testthat/test-bug_fixes.R b/tests/testthat/test-bug_fixes.R
new file mode 100644
index 0000000..0e0d4bf
--- /dev/null
+++ b/tests/testthat/test-bug_fixes.R
@@ -0,0 +1,138 @@
+library(graDiEnt)
+
+# ── grad_approx_fn: zero-division guards ──────────────────────────────────────
+
+test_that("grad_approx_fn returns finite values when two particles are identical", {
+  # identical particles → vec_diff_temp == 0 → original code divides by zero
+  params <- matrix(c(1, 2, 1, 2, 3, 4, 5, 6), nrow = 4, ncol = 2)
+  params[2, ] <- params[1, ]   # particles 1 and 2 are identical
+  weights <- c(1.0, 0.5, 2.0, 1.5)
+
+  result <- graDiEnt:::grad_approx_fn(
+    param_indices  = 1:2,
+    n_diff         = 1,
+    current_params = params,
+    parent_indices = c(1, 2),   # identical pair
+    current_weight = weights
+  )
+
+  expect_true(all(is.finite(result$grad_approx)))
+  expect_true(is.finite(result$psi) | is.nan(result$psi))  # psi may be NaN when grad is zero vector
+})
+
+test_that("grad_approx_fn returns finite grad when two particles have identical weights", {
+  # identical weights → weight_diff_temp == 0 → original skips this implicitly,
+  # but the guard makes it explicit and prevents any numerical issues
+  params <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8), nrow = 4, ncol = 2)
+  weights <- c(1.0, 1.0, 2.0, 1.5)   # particles 1 and 2 have same weight
+
+  result <- graDiEnt:::grad_approx_fn(
+    param_indices  = 1:2,
+    n_diff         = 1,
+    current_params = params,
+    parent_indices = c(1, 2),
+    current_weight = weights
+  )
+
+  expect_true(all(is.finite(result$grad_approx)))
+})
+
+test_that("grad_approx_fn handles mixed valid and zero-diff pairs without NaN", {
+  # n_diff = 2: one pair identical, one pair distinct — guard skips bad pair
+  params <- matrix(c(1,2, 1,2, 3,4, 5,6, 7,8), nrow = 5, ncol = 2, byrow = TRUE)
+  weights <- c(1.0, 1.0, 0.5, 2.0, 1.5)
+
+  result <- graDiEnt:::grad_approx_fn(
+    param_indices  = 1:2,
+    n_diff         = 2,
+    current_params = params,
+    parent_indices = c(1, 3, 2, 4),   # pair (1,2) identical; pair (3,4) distinct
+    current_weight = weights
+  )
+
+  expect_true(all(is.finite(result$grad_approx)))
+})
+
+# ── degenerate population: optimizer doesn't crash ───────────────────────────
+
+test_that("optimizer runs without error when particles initialise very close together", {
+  # extremely tight init forces near-identical particles early on,
+  # triggering the zero-division guard repeatedly
+  set.seed(99)
+  expect_no_error(
+    suppressMessages(
+      optim_SQGDE(
+        ObjFun = function(x) sum((x - c(1, 2))^2),
+        control_params = GetAlgoParams(
+          n_params    = 2,
+          n_iter      = 50,
+          n_particles = 8,
+          n_diff      = 2,
+          init_sd     = 1e-8,
+          init_center = 0
+        )
+      )
+    )
+  )
+})
+
+# ── rand scheme sampling: correct number of parents ───────────────────────────
+
+test_that("rand scheme works with minimum viable n_particles and n_diff", {
+  # n_particles=4, n_diff=1 for rand requires sampling 2*1+1=3 parents from 4
+  # the original SQG_DE_bin_1_pos bug would sample only 2 and then access index 3
+  set.seed(7)
+  out <- suppressMessages(
+    optim_SQGDE(
+      ObjFun = function(x) sum((x - 1)^2),
+      control_params = GetAlgoParams(
+        n_params     = 1,
+        n_iter       = 100,
+        n_particles  = 6,
+        n_diff       = 1,
+        adapt_scheme = "rand",
+        init_sd      = 1
+      )
+    )
+  )
+  expect_true(is.finite(out$weight))
+  expect_true(all(is.finite(out$solution)))
+})
+
+test_that("best scheme works with minimum viable n_particles and n_diff", {
+  set.seed(7)
+  out <- suppressMessages(
+    optim_SQGDE(
+      ObjFun = function(x) sum((x - 1)^2),
+      control_params = GetAlgoParams(
+        n_params     = 1,
+        n_iter       = 100,
+        n_particles  = 6,
+        n_diff       = 1,
+        adapt_scheme = "best",
+        init_sd      = 1
+      )
+    )
+  )
+  expect_true(is.finite(out$weight))
+  expect_true(all(is.finite(out$solution)))
+})
+
+test_that("current scheme works with minimum viable n_particles and n_diff", {
+  set.seed(7)
+  out <- suppressMessages(
+    optim_SQGDE(
+      ObjFun = function(x) sum((x - 1)^2),
+      control_params = GetAlgoParams(
+        n_params     = 1,
+        n_iter       = 100,
+        n_particles  = 6,
+        n_diff       = 1,
+        adapt_scheme = "current",
+        init_sd      = 1
+      )
+    )
+  )
+  expect_true(is.finite(out$weight))
+  expect_true(all(is.finite(out$solution)))
+})
diff --git a/tests/testthat/test-trace_print_freq.R b/tests/testthat/test-trace_print_freq.R
new file mode 100644
index 0000000..7304c99
--- /dev/null
+++ b/tests/testthat/test-trace_print_freq.R
@@ -0,0 +1,75 @@
+library(graDiEnt)
+
+obj <- function(x) sum((x - 1)^2)
+
+# helper: run optimizer and return only "iter N/M" progress messages
+capture_iter_messages <- function(control_params) {
+  msgs <- character(0)
+  withCallingHandlers(
+    optim_SQGDE(obj, control_params),
+    message = function(m) {
+      txt <- conditionMessage(m)
+      if (grepl("^iter ", txt)) msgs <<- c(msgs, txt)
+      invokeRestart("muffleMessage")
+    }
+  )
+  msgs
+}
+
+# ── Inf suppresses all iter messages ─────────────────────────────────────────
+
+test_that("trace_print_freq = Inf produces no iter progress messages", {
+  set.seed(1)
+  msgs <- capture_iter_messages(
+    GetAlgoParams(n_params = 1, n_iter = 200, n_particles = 6,
+                  n_diff = 1, trace_print_freq = Inf)
+  )
+  expect_length(msgs, 0)
+})
+
+# ── custom freq prints at correct multiples ───────────────────────────────────
+
+test_that("trace_print_freq = 10 prints at iterations 10 and 20", {
+  set.seed(1)
+  msgs <- capture_iter_messages(
+    GetAlgoParams(n_params = 1, n_iter = 25, n_particles = 6,
+                  n_diff = 1, trace_print_freq = 10, stop_check = 30)
+  )
+  expect_length(msgs, 2)
+  expect_match(msgs[1], "^iter 10/")
+  expect_match(msgs[2], "^iter 20/")
+})
+
+test_that("trace_print_freq = 50 prints at iterations 50 and 100", {
+  set.seed(1)
+  msgs <- capture_iter_messages(
+    GetAlgoParams(n_params = 1, n_iter = 110, n_particles = 6,
+                  n_diff = 1, trace_print_freq = 50, stop_check = 120)
+  )
+  expect_length(msgs, 2)
+  expect_match(msgs[1], "^iter 50/")
+  expect_match(msgs[2], "^iter 100/")
+})
+
+# ── validation errors ─────────────────────────────────────────────────────────
+
+test_that("GetAlgoParams rejects trace_print_freq = 0", {
+  expect_error(
+    GetAlgoParams(n_params = 2, trace_print_freq = 0),
+    "positive integer or Inf"
+  )
+})
+
+test_that("GetAlgoParams rejects negative trace_print_freq", {
+  expect_error(
+    GetAlgoParams(n_params = 2, trace_print_freq = -10),
+    "positive integer or Inf"
+  )
+})
+
+test_that("GetAlgoParams rejects non-scalar trace_print_freq", {
+  expect_error(
+    GetAlgoParams(n_params = 2, trace_print_freq = c(10, 20)),
+    "scalar"
+  )
+})

From 35c745f8a043007cff7f7ce8ffd9f7eae253e62d Mon Sep 17 00:00:00 2001
From: Brendan Matthew Galdo <brendan.m.galdo@gmail.com>
Date: Mon, 11 May 2026 21:28:01 -0400
Subject: [PATCH 4/5] add recovery feature

---
 R/GetAlgoParams.R              |  22 +++
 R/optim_SQGDE.R                |  17 +++
 tests/testthat/test-recovery.R | 235 +++++++++++++++++++++++++++++++++
 3 files changed, 274 insertions(+)
 create mode 100644 tests/testthat/test-recovery.R

diff --git a/R/GetAlgoParams.R b/R/GetAlgoParams.R
index 15ec509..9438630 100644
--- a/R/GetAlgoParams.R
+++ b/R/GetAlgoParams.R
@@ -11,6 +11,8 @@
 #' @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.
@@ -39,6 +41,8 @@ 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,
@@ -200,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
@@ -321,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,
diff --git a/R/optim_SQGDE.R b/R/optim_SQGDE.R
index 067537a..53f836a 100644
--- a/R/optim_SQGDE.R
+++ b/R/optim_SQGDE.R
@@ -134,6 +134,9 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), warm_start = NU
     parallel::clusterSetRNGStream(cl_use)
   }
 
+  if(!is.null(control_params$recovery_path)){
+    message(paste0('recovery enabled: writing partial results to "', control_params$recovery_path, '" each iteration'))
+  }
   message("running SQG-DE...")
 
 
@@ -220,6 +223,20 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), warm_start = NU
       }
     }
 
+    if(!is.null(control_params$recovery_path) && iter %% control_params$recovery_freq == 0){
+      .minIdx_r <- which.min(weights[iter_idx, ])
+      saveRDS(list(
+        solution       = particles[iter_idx, .minIdx_r, ],
+        weight         = weights[iter_idx, .minIdx_r],
+        last_particles = matrix(particles[iter_idx, , ],
+                                nrow = control_params$n_particles,
+                                ncol = control_params$n_params),
+        last_weights   = weights[iter_idx, ],
+        converged      = FALSE,
+        iter_completed = iter
+      ), file = control_params$recovery_path)
+    }
+
     if(iter %% control_params$stop_check==0){
       # assign convergence method
 
diff --git a/tests/testthat/test-recovery.R b/tests/testthat/test-recovery.R
new file mode 100644
index 0000000..8b13b51
--- /dev/null
+++ b/tests/testthat/test-recovery.R
@@ -0,0 +1,235 @@
+library(graDiEnt)
+
+obj <- function(x) sum((x - 1)^2)
+
+# ── GetAlgoParams validation ──────────────────────────────────────────────────
+
+test_that("GetAlgoParams accepts recovery_path = NULL", {
+  cp <- GetAlgoParams(n_params = 2, recovery_path = NULL)
+  expect_null(cp$recovery_path)
+})
+
+test_that("GetAlgoParams accepts a character scalar recovery_path", {
+  cp <- GetAlgoParams(n_params = 2, recovery_path = tempfile(fileext = ".rds"))
+  expect_true(is.character(cp$recovery_path))
+})
+
+test_that("GetAlgoParams rejects non-character recovery_path", {
+  expect_error(
+    GetAlgoParams(n_params = 2, recovery_path = 123),
+    "character scalar"
+  )
+})
+
+test_that("GetAlgoParams rejects non-scalar recovery_path", {
+  expect_error(
+    GetAlgoParams(n_params = 2, recovery_path = c("a.rds", "b.rds")),
+    "character scalar"
+  )
+})
+
+test_that("GetAlgoParams rejects recovery_path without .rds extension", {
+  expect_error(
+    GetAlgoParams(n_params = 2, recovery_path = "output.txt"),
+    "\\.rds"
+  )
+})
+
+test_that("GetAlgoParams rejects recovery_path with no extension", {
+  expect_error(
+    GetAlgoParams(n_params = 2, recovery_path = "myfile"),
+    "\\.rds"
+  )
+})
+
+test_that("GetAlgoParams accepts .RDS extension (case-insensitive)", {
+  expect_no_error(GetAlgoParams(n_params = 2, recovery_path = "output.RDS"))
+})
+
+# ── recovery_freq validation ──────────────────────────────────────────────────
+
+test_that("GetAlgoParams accepts recovery_freq = 1", {
+  cp <- GetAlgoParams(n_params = 2, recovery_freq = 1)
+  expect_equal(cp$recovery_freq, 1L)
+})
+
+test_that("GetAlgoParams accepts recovery_freq = 10", {
+  cp <- GetAlgoParams(n_params = 2, recovery_freq = 10)
+  expect_equal(cp$recovery_freq, 10L)
+})
+
+test_that("GetAlgoParams rejects recovery_freq = 0", {
+  expect_error(GetAlgoParams(n_params = 2, recovery_freq = 0), "positive integer")
+})
+
+test_that("GetAlgoParams rejects negative recovery_freq", {
+  expect_error(GetAlgoParams(n_params = 2, recovery_freq = -5), "positive integer")
+})
+
+test_that("GetAlgoParams rejects non-scalar recovery_freq", {
+  expect_error(GetAlgoParams(n_params = 2, recovery_freq = c(1, 2)), "positive integer")
+})
+
+# ── recovery file written during run ─────────────────────────────────────────
+
+test_that("recovery file is created and readable after run", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 1, n_iter = 50, n_particles = 6,
+                              n_diff = 1, recovery_path = f))
+  )
+
+  expect_true(file.exists(f))
+  rec <- readRDS(f)
+  expect_true(is.list(rec))
+})
+
+test_that("recovery file has correct structure", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 2, n_iter = 50, n_particles = 8,
+                              n_diff = 1, recovery_path = f))
+  )
+
+  rec <- readRDS(f)
+  expect_named(rec, c("solution", "weight", "last_particles", "last_weights",
+                      "converged", "iter_completed"),
+               ignore.order = TRUE)
+})
+
+test_that("recovery file solution and weight are finite", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 2, n_iter = 50, n_particles = 8,
+                              n_diff = 1, recovery_path = f))
+  )
+
+  rec <- readRDS(f)
+  expect_true(all(is.finite(rec$solution)))
+  expect_true(is.finite(rec$weight))
+})
+
+test_that("recovery last_particles has correct dimensions", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 3, n_iter = 50, n_particles = 12,
+                              n_diff = 1, recovery_path = f))
+  )
+
+  rec <- readRDS(f)
+  expect_equal(dim(rec$last_particles), c(12L, 3L))
+  expect_equal(length(rec$last_weights), 12L)
+})
+
+test_that("recovery iter_completed equals n_iter when run completes", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 1, n_iter = 30, n_particles = 6,
+                              n_diff = 1, recovery_path = f,
+                              stop_check = 40))  # prevent early stop
+  )
+
+  rec <- readRDS(f)
+  expect_equal(rec$iter_completed, 30L)
+})
+
+test_that("recovery converged field is FALSE", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 1, n_iter = 30, n_particles = 6,
+                              n_diff = 1, recovery_path = f,
+                              stop_check = 40))
+  )
+
+  rec <- readRDS(f)
+  expect_false(rec$converged)
+})
+
+# ── recovery file does not contain trace arrays ───────────────────────────────
+
+test_that("recovery file excludes trace arrays even when return_trace = TRUE", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 1, n_iter = 30, n_particles = 6,
+                              n_diff = 1, recovery_path = f,
+                              return_trace = TRUE, stop_check = 40))
+  )
+
+  rec <- readRDS(f)
+  expect_false("particles_trace" %in% names(rec))
+  expect_false("weights_trace"   %in% names(rec))
+})
+
+# ── recovery_freq controls save timing ───────────────────────────────────────
+
+test_that("recovery file not written when n_iter < recovery_freq", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 1, n_iter = 5, n_particles = 6,
+                              n_diff = 1, recovery_path = f, recovery_freq = 10,
+                              stop_check = 20))
+  )
+  expect_false(file.exists(f))
+})
+
+test_that("recovery file written when iter reaches recovery_freq multiple", {
+  f <- tempfile(fileext = ".rds")
+  on.exit(unlink(f))
+
+  set.seed(1)
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 1, n_iter = 20, n_particles = 6,
+                              n_diff = 1, recovery_path = f, recovery_freq = 10,
+                              stop_check = 30))
+  )
+  expect_true(file.exists(f))
+  rec <- readRDS(f)
+  expect_true(rec$iter_completed %% 10 == 0)
+})
+
+# ── no file written when recovery_path = NULL ─────────────────────────────────
+
+test_that("no rds file written when recovery_path is NULL", {
+  set.seed(1)
+  before <- list.files(tempdir(), pattern = "\\.rds$")
+  suppressMessages(
+    optim_SQGDE(obj,
+                GetAlgoParams(n_params = 1, n_iter = 20, n_particles = 6,
+                              n_diff = 1, recovery_path = NULL))
+  )
+  after <- list.files(tempdir(), pattern = "\\.rds$")
+  expect_equal(before, after)
+})

From c4a5503d417d7a960b26ff4aa7e1c78c7b3253f7 Mon Sep 17 00:00:00 2001
From: Brendan Matthew Galdo <brendan.m.galdo@gmail.com>
Date: Mon, 11 May 2026 21:28:40 -0400
Subject: [PATCH 5/5] Update GetAlgoParams.Rd

update documentation
---
 man/GetAlgoParams.Rd | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/man/GetAlgoParams.Rd b/man/GetAlgoParams.Rd
index 7e08409..b069be5 100644
--- a/man/GetAlgoParams.Rd
+++ b/man/GetAlgoParams.Rd
@@ -19,6 +19,8 @@ GetAlgoParams(
   jitter_size = 1e-06,
   crossover_rate = 1,
   parallel_type = "none",
+  recovery_path = NULL,
+  recovery_freq = 1,
   return_trace = FALSE,
   thin = 1,
   purify = Inf,
@@ -26,7 +28,8 @@ GetAlgoParams(
   give_up_init = 100,
   stop_check = 10,
   stop_tol = 1e-04,
-  converge_crit = "stdev"
+  converge_crit = "stdev",
+  trace_print_freq = 100
 )
 }
 \arguments{
@@ -54,10 +57,14 @@ GetAlgoParams(
 
 \item{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.}
 
-\item{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.}
+\item{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.}
 
 \item{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.}
 
+\item{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.}
+
+\item{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}.}
+
 \item{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.}
 
 \item{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.}
@@ -73,6 +80,8 @@ GetAlgoParams(
 \item{stop_tol}{A convergence metric must be less than value to be labeled as converged. The default is 1e-4.}
 
 \item{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.}
+
+\item{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.}
 }
 \value{
 A list of control parameters for the optim_SQGDE function.