diff --git a/.Rprofile b/.Rprofile deleted file mode 100644 index 0f860de..0000000 --- a/.Rprofile +++ /dev/null @@ -1,2 +0,0 @@ -library(roxygen2) -library(devtools) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index f36e31f..e2f46c8 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -8,27 +8,16 @@ name: R-CMD-check jobs: R-CMD-check: - runs-on: ${{ matrix.config.os }} - name: ${{ matrix.config.os }} (R ${{ matrix.config.r }}) - - strategy: - fail-fast: false - matrix: - config: - - {os: macos-latest, r: 'release'} - - {os: windows-latest, r: 'release'} - - {os: ubuntu-latest, r: 'release'} - - {os: ubuntu-latest, r: 'devel'} + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - uses: actions/checkout@v4 - - uses: r-lib/actions/setup-pandoc@v2 - - uses: r-lib/actions/setup-r@v2 with: - r-version: ${{ matrix.config.r }} - use-public-rspm: true + r-version: 'release' - uses: r-lib/actions/setup-r-dependencies@v2 with: @@ -37,4 +26,5 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: + args: 'c("--no-manual", "--as-cran")' upload-snapshots: true diff --git a/.gitignore b/.gitignore index ee875fc..b4aaf6d 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,7 @@ tests/testthat/_snaps/ # devtools / roxygen /doc/ /Meta/ + +# local R session config (machine-specific) +.Rprofile +.Renviron diff --git a/DESCRIPTION b/DESCRIPTION index 45826ab..0267e70 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -19,5 +19,4 @@ Imports: parallel Suggests: testthat (>= 3.0.0) -Roxygen: list(markdown = TRUE) -RoxygenNote: 7.1.2 +RoxygenNote: 7.3.2 diff --git a/R/GetAlgoParams.R b/R/GetAlgoParams.R index f239ad6..b999e63 100644 --- a/R/GetAlgoParams.R +++ b/R/GetAlgoParams.R @@ -18,6 +18,8 @@ #' @param give_up_init An integer for how many failed initialization attempts before stopping the optimization routine. 100 is the default value. #' @param stop_check An integer for how often to check the convergence criterion. The default is 10 iterations. #' @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 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 @@ -27,6 +29,8 @@ GetAlgoParams = function(n_params, n_iter = 1000, init_sd = 0.01, init_center = 0, + lower = -Inf, + upper = Inf, n_cores_use = 1, step_size = NULL, jitter_size = 1e-6, @@ -95,6 +99,28 @@ GetAlgoParams = function(n_params, stop('ERROR: init_center vector length must be 1 or n_params') } + # lower + lower = as.numeric(lower) + if(any(is.nan(lower))){ + stop('ERROR: lower contains NaN') + } else if(!(length(lower) == 1 | length(lower) == n_params)){ + stop('ERROR: lower vector length must be 1 or n_params') + } + if(length(lower) == 1) lower = rep(lower, n_params) + + # upper + upper = as.numeric(upper) + if(any(is.nan(upper))){ + stop('ERROR: upper contains NaN') + } else if(!(length(upper) == 1 | length(upper) == n_params)){ + stop('ERROR: upper vector length must be 1 or n_params') + } + if(length(upper) == 1) upper = rep(upper, n_params) + + if(any(lower >= upper)){ + stop('ERROR: lower must be strictly less than upper for all parameters') + } + # n_cores_use ### assign NULL value default if(is.null(n_cores_use)){ @@ -264,6 +290,8 @@ GetAlgoParams = function(n_params, 'n_iter' = n_iter, 'init_sd' = init_sd, 'init_center' = init_center, + 'lower' = lower, + 'upper' = upper, '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 bcbfbee..5b595f3 100644 --- a/R/SQG_DE_bin_1_best.R +++ b/R/SQG_DE_bin_1_best.R @@ -17,8 +17,10 @@ SQG_DE_bin_1_best=function(pmem_index, objFun, step_size = .8, jitter_size = 1e-6, - n_particles, + n_particles, crossover_rate = 1, + lower = -Inf, + upper = Inf, n_diff, ... ){ # get statistics about particle @@ -72,6 +74,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 = 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 f88545a..10b7210 100644 --- a/R/SQG_DE_bin_1_curr.R +++ b/R/SQG_DE_bin_1_curr.R @@ -18,6 +18,8 @@ SQG_DE_bin_1_curr=function(pmem_index, jitter_size = 1e-6, n_particles, crossover_rate = 1, + lower = -Inf, + upper = Inf, n_diff, ... ){ # sample parents @@ -71,6 +73,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 = 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 a448cc2..2a700e4 100644 --- a/R/SQG_DE_bin_1_rand.R +++ b/R/SQG_DE_bin_1_rand.R @@ -19,6 +19,8 @@ SQG_DE_bin_1_rand=function(pmem_index, jitter_size=1e-6, n_particles, crossover_rate=1, + lower=-Inf, + upper=Inf, n_diff, ... ){ @@ -72,6 +74,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 = matrix(params_use,1,len_param_use) weight_proposal = NA diff --git a/R/optim_SQGDE.R b/R/optim_SQGDE.R index f53a98f..195abf2 100644 --- a/R/optim_SQGDE.R +++ b/R/optim_SQGDE.R @@ -77,9 +77,11 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), ...){ for(pmem_index in 1:control_params$n_particles){ count = 0 # establish a count variable to avoid infinite run time while(weights[1,pmem_index]==Inf) { - particles[1, pmem_index, ] = stats::rnorm(control_params$n_params, - control_params$init_center, - control_params$init_sd) + particles[1, pmem_index, ] = pmax(control_params$lower, + pmin(control_params$upper, + stats::rnorm(control_params$n_params, + control_params$init_center, + control_params$init_sd))) weights[1, pmem_index] = ObjFun(particles[1, pmem_index, ], ...) @@ -139,6 +141,8 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), ...){ jitter_size = control_params$jitter_size, n_particles = control_params$n_particles, crossover_rate = control_params$crossover_rate, + lower = control_params$lower, + upper = control_params$upper, n_diff = control_params$n_diff, ...)), control_params$n_particles, control_params$n_params+1, byrow=TRUE) @@ -152,6 +156,8 @@ optim_SQGDE = function(ObjFun, control_params = GetAlgoParams(), ...){ jitter_size = control_params$jitter_size, n_particles = control_params$n_particles, crossover_rate = control_params$crossover_rate, + lower = control_params$lower, + upper = control_params$upper, 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 1bdee6e..6fa70b9 100644 --- a/man/GetAlgoParams.Rd +++ b/man/GetAlgoParams.Rd @@ -11,6 +11,8 @@ GetAlgoParams( n_iter = 1000, init_sd = 0.01, init_center = 0, + lower = -Inf, + upper = Inf, n_cores_use = 1, step_size = NULL, jitter_size = 1e-06, @@ -39,6 +41,10 @@ GetAlgoParams( \item{init_center}{A scalar or n_params-dimensional numeric vector, determines the mean of the Gaussian initialization distribution. The default value is 0.} +\item{lower}{A numeric scalar or n_params-dimensional vector specifying lower bounds for each parameter. Default is -Inf (no lower bound).} + +\item{upper}{A numeric scalar or n_params-dimensional vector specifying upper bounds for each parameter. Default is Inf (no upper bound).} + \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).} diff --git a/tests/testthat/test-bounds.R b/tests/testthat/test-bounds.R new file mode 100644 index 0000000..a9804da --- /dev/null +++ b/tests/testthat/test-bounds.R @@ -0,0 +1,148 @@ +library(graDiEnt) + +# ── GetAlgoParams: bounds validation ───────────────────────────────────────── + +test_that("GetAlgoParams rejects lower >= upper", { + expect_error(GetAlgoParams(n_params = 2, lower = 1, upper = 0), + "lower must be strictly less than upper") + expect_error(GetAlgoParams(n_params = 2, lower = c(0, 5), upper = c(1, 5)), + "lower must be strictly less than upper") +}) + +test_that("GetAlgoParams rejects wrong-length bounds", { + expect_error(GetAlgoParams(n_params = 2, lower = c(0, 0, 0)), + "lower vector length must be 1 or n_params") + expect_error(GetAlgoParams(n_params = 2, upper = c(1, 1, 1)), + "upper vector length must be 1 or n_params") +}) + +test_that("GetAlgoParams stores bounds as length-n_params vectors", { + cp <- GetAlgoParams(n_params = 3, lower = -1, upper = 2) + expect_length(cp$lower, 3) + expect_length(cp$upper, 3) + expect_equal(cp$lower, c(-1, -1, -1)) + expect_equal(cp$upper, c(2, 2, 2)) +}) + +test_that("GetAlgoParams accepts per-parameter bounds", { + cp <- GetAlgoParams(n_params = 3, lower = c(-1, 0, -5), upper = c(1, 3, 5)) + expect_equal(cp$lower, c(-1, 0, -5)) + expect_equal(cp$upper, c( 1, 3, 5)) +}) + +# ── optim_SQGDE: solution stays within symmetric scalar bounds ──────────────── + +test_that("solution stays within scalar bounds when optimum is inside", { + set.seed(42) + out <- suppressMessages( + optim_SQGDE( + ObjFun = function(x) sum((x - c(1, 2))^2), + control_params = GetAlgoParams( + n_params = 2, + n_iter = 400, + n_particles = 10, + n_diff = 2, + init_sd = 1, + init_center = 0, + lower = -5, + upper = 5 + ) + ) + ) + + expect_true(all(out$solution >= -5)) + expect_true(all(out$solution <= 5)) + expect_lt(max(abs(out$solution - c(1, 2))), 0.15) +}) + +# ── optim_SQGDE: per-parameter asymmetric bounds ───────────────────────────── + +test_that("solution stays within per-parameter asymmetric bounds", { + set.seed(7) + true_opt <- c(0.5, 2.5, -0.5) + out <- suppressMessages( + optim_SQGDE( + ObjFun = function(x) sum((x - true_opt)^2), + control_params = GetAlgoParams( + n_params = 3, + n_iter = 500, + n_particles = 12, + n_diff = 2, + init_sd = 0.3, + init_center = c(0.5, 2, -0.5), + lower = c( 0, 1, -1), + upper = c( 1, 3, 0) + ) + ) + ) + + lower <- c(0, 1, -1) + upper <- c(1, 3, 0) + expect_true(all(out$solution >= lower)) + expect_true(all(out$solution <= upper)) + expect_lt(max(abs(out$solution - true_opt)), 0.15) +}) + +# ── optim_SQGDE: optimum outside bounds → clips to boundary ────────────────── + +test_that("solution clips to boundary when true optimum is outside bounds", { + set.seed(99) + out <- suppressMessages( + optim_SQGDE( + ObjFun = function(x) sum((x - c(-10, -10))^2), + control_params = GetAlgoParams( + n_params = 2, + n_iter = 400, + n_particles = 10, + n_diff = 2, + init_sd = 0.5, + init_center = 1, + lower = 0, + upper = 2 + ) + ) + ) + + expect_true(all(out$solution >= 0)) + expect_true(all(out$solution <= 2)) + # constrained optimum is at lower boundary (0, 0) + expect_lt(max(abs(out$solution - c(0, 0))), 0.1) +}) + +# ── optim_SQGDE: unbounded by default (no behaviour change) ────────────────── + +test_that("default bounds (-Inf / Inf) do not alter unbounded behaviour", { + obj <- function(x) sum((x - c(-2, 3))^2) + + set.seed(5) + out_nobound <- suppressMessages( + optim_SQGDE( + ObjFun = obj, + control_params = GetAlgoParams( + n_params = 2, + n_iter = 400, + n_particles = 10, + n_diff = 2, + init_sd = 1 + ) + ) + ) + + set.seed(5) + out_explicitinf <- suppressMessages( + optim_SQGDE( + ObjFun = obj, + control_params = GetAlgoParams( + n_params = 2, + n_iter = 400, + n_particles = 10, + n_diff = 2, + init_sd = 1, + lower = -Inf, + upper = Inf + ) + ) + ) + + expect_equal(out_nobound$solution, out_explicitinf$solution) +}) diff --git a/tests/testthat/test-optim_SQGDE.R b/tests/testthat/test-optim_SQGDE.R index a2dce99..a7228b8 100644 --- a/tests/testthat/test-optim_SQGDE.R +++ b/tests/testthat/test-optim_SQGDE.R @@ -5,7 +5,8 @@ library(graDiEnt) test_that("MLE factorizable MVN recovers sample means", { set.seed(42) true_mu <- c(-1, 1, 0, 2) - data <- matrix(rnorm(500 * 4, mean = true_mu, sd = 1), nrow = 500, ncol = 4) + data <- matrix(rnorm(500 * 4, mean = true_mu, sd = 1), + nrow = 500, ncol = 4 ,byrow = TRUE) analytic_mu <- colMeans(data) neg_log_lik <- function(x, data) { @@ -44,7 +45,7 @@ test_that("MLE factorizable MVN recovers sample means", { test_that("univariate objective recovers sample mean", { set.seed(42) - obs <- rnorm(300, mean = 2.5, sd = 1) + obs <- rnorm(300, mean = 4, sd = 1) analytic_mu <- mean(obs) neg_log_lik_1d <- function(x, obs) { @@ -59,9 +60,9 @@ test_that("univariate objective recovers sample mean", { n_params = 1, n_iter = 500, n_particles = 6, - n_diff = 1, + n_diff = 2, init_center = 0, - init_sd = 2 + init_sd = 1 ), obs = obs ) @@ -155,8 +156,9 @@ test_that("PSOCK parallel execution returns valid solution", { skip_on_cran() set.seed(42) - true_mu <- c(-1, 1, 0, 2) - data_par <- matrix(rnorm(500 * 4, mean = true_mu, sd = 1), nrow = 500, ncol = 4) + true_mu <- c(-1, 8, 0, 2) + data_par <- matrix(rnorm(500 * 4, mean = true_mu, sd = 1), + nrow = 500, ncol = 4 ,byrow = TRUE) analytic_mu <- colMeans(data_par) neg_log_lik <- function(x, data_par) { @@ -201,7 +203,8 @@ test_that("FORK parallel execution returns valid solution", { set.seed(42) true_mu <- c(-1, 1, 0, 2) - data_par <- matrix(rnorm(500 * 4, mean = true_mu, sd = 1), nrow = 500, ncol = 4) + data_par <- matrix(rnorm(500 * 4, mean = true_mu, sd = 1), + nrow = 500, ncol = 4 ,byrow = TRUE) analytic_mu <- colMeans(data_par) neg_log_lik <- function(x, data_par) {