## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(ParmOff) library(magicaxis) ## ----basics-drop-------------------------------------------------------------- model <- function(x, y, z) x * y + z # 't' is not a formal of model — ParmOff silently drops it params <- c(x = 1, y = 2, z = 3, t = 4) ParmOff(model, params) ## ----basics-list-vs-vec------------------------------------------------------- # Named numeric vector ParmOff(model, c(x = 2, y = 3, z = 1)) # Named list — preferred when types differ ParmOff(model, list(x = 2L, y = 3.0, z = TRUE)) ## ----basics-dots-------------------------------------------------------------- # z is in .args; the conflicting z in ... is ignored ParmOff(model, list(x = 1, y = 2, z = 3), z = 99) # z is missing from .args but supplied via ... ParmOff(model, list(x = 1, y = 2), z = 3) ## ----use-args----------------------------------------------------------------- f <- function(x, y) x + y # Only x and y are passed; the extra z never reaches f ParmOff(f, list(x = 2, y = 3, z = 99), .use_args = c("x", "y")) ## ----rem-args----------------------------------------------------------------- # Remove z before passing; f receives x and y only f_xyz <- function(x = 1, y = 2, z = 3) x * y + z ParmOff(f_xyz, list(x = 2, y = 3, z = 99), .rem_args = "z") ## ----use-rem-combined--------------------------------------------------------- # Keep x, y, z — then remove z ParmOff(f, list(x = 2, y = 3, z = 99, w = 0), .use_args = c("x", "y", "z"), .rem_args = "z") ## ----strip-------------------------------------------------------------------- # A parameter vector from, say, a Bayesian sampler with a "fit." prefix fit_params <- c(fit.x = 1, fit.y = 2, fit.z = 3, fit.t = 4) ParmOff(model, fit_params, .strip = "fit\\.") ## ----strip-with-use----------------------------------------------------------- ParmOff(f, list(p.x = 2, p.y = 3, p.z = 99), .strip = "p\\.", .use_args = c("x", "y")) ## ----logged------------------------------------------------------------------- # y is stored as log10(2) ≈ 0.301; 10^0.301 ≈ 2 ParmOff(model, list(x = 1, y = log10(2), z = 3), .logged = "y") ## ----logged-multi------------------------------------------------------------- # Both x and y are in log10 space ParmOff(model, list(x = 1, y = 1, z = 3), .logged = c("x", "y")) # x = 10^1 = 10, y = 10^1 = 10 → 10*10 + 3 = 103 ## ----bounds-basic------------------------------------------------------------- # Clamp y upward, z downward ParmOff(model, list(x = 1, y = 0.1, z = 15), .lower = c(y = 1), .upper = c(z = 10)) # y was 0.1, clamped to 1 → 1*1 + 10 = 11 ## ----bounds-raw-true---------------------------------------------------------- # y is in log10 space; clamp to [-1, 1] before back-transforming # log10(y) clamped to 1 → 10^1 = 10 ParmOff(model, list(x = 1, y = 2, z = 3), .logged = "y", .lower = c(y = -1), .upper = c(y = 1)) # Result: 1 * 10 + 3 = 13 ## ----bounds-raw-false--------------------------------------------------------- # y = 2 (log10) → 10^2 = 100; upper 5 (real) → clamp to 5 # Result: 1 * 5 + 3 = 8 ParmOff(model, list(x = 1, y = 2, z = 3), .logged = "y", .upper = c(y = 5), .bound_raw = FALSE) ## ----parmlim-pure-vectors----------------------------------------------------- # Named vector: partial named bound — only 'a' is clamped ParmLimLo(c(a = -5, b = 10, c = 3), lower = c(a = 0)) # Named vectors: both lower and upper, leaving c untouched ParmLimBoth(c(a = -5, b = 10, c = 3), lower = c(a = 0), upper = c(b = 7)) ## ----parmlim-list-partial----------------------------------------------------- # Partial named-list bound: only supply bounds for the children you need x <- list(a = -1, b = 5, c = -3) ParmLimLo(x, lower = list(a = 0, c = 0)) # b is left unchanged ## ----parmlim-scalar-broadcast------------------------------------------------- # Scalar broadcast: one value applied to every leaf, regardless of nesting deep <- list(p = list(q = list(r = -5, s = 20), t = 3), u = -1) ParmLimBoth(deep, lower = 0, upper = 10) ## ----parmlim-deep-partial----------------------------------------------------- # Deep partial bound: supply bounds only at the levels you need ParmLimBoth(deep, lower = list(p = list(q = list(r = 0), t = 0), u = 0), upper = list(p = list(q = list(s = 10), t = 5)) ) # p$q$r: 0 (was -5) p$q$s: 10 (was 20) # p$t: 3 (unchanged) u: 0 (was -1) ## ----parmlog-character-------------------------------------------------------- params <- list(amplitude = 100, scale = 10, offset = 3) # Forward-transform two parameters to log10 space logged_params <- ParmLog(params, logged = c("amplitude", "scale")) logged_params # Back-transform them ParmUnLog(logged_params, logged = c("amplitude", "scale")) ## ----parmlog-logical---------------------------------------------------------- # Logical-vector selector: transform the first two elements ParmUnLog(list(a = 2, b = 1, c = 5), logged = c(TRUE, TRUE, FALSE)) # a = 10^2 = 100, b = 10^1 = 10, c unchanged ## ----parmlog-ln--------------------------------------------------------------- # Natural-log flavour ParmLog(list(sigma = exp(3), mu = 0), logged = "sigma", log_type = 'ln') # sigma = log(exp(3)) = 3 ## ----parmlog-matrix----------------------------------------------------------- # Matrix elements retain their shape mat_params <- list(cov = matrix(c(100, 0, 0, 100), 2, 2), mu = 5) out <- ParmLog(mat_params, logged = "cov") out$cov # still a 2×2 matrix, values are log10 of originals ## ----------------------------------------------------------------------------- model_ex = function(x,y,z){x * y + z} input = c(x=1, y=2, z=3) ParmOff(model_ex, input, .logged='y', .lower=list(y=0), .upper=list(y=1), .return='args') #Make y constrained to be 2*x: constrain_func = function(x, y, z){list(x=x, y = 2 * x, z=z)} ParmOff(model_ex, input, .logged='y', .lower=list(y=0), .upper=list(y=1), .constrain=constrain_func, .return='args') ## ----dnorm-optim-------------------------------------------------------------- set.seed(42) data_obs <- rnorm(200, mean = 5, sd = 2) # Negative log-likelihood; 'sd' is optimised in log10 space neg_ll <- function(par, data) { -sum(dnorm(data, mean = par["mean"], sd = 10^par["log_sd"], log = TRUE)) } # Starting values: mean ~ 0, log10(sd) ~ 0 (i.e. sd ~ 1) start <- c(mean = 0, log_sd = 0) fit <- optim(start, neg_ll, data = data_obs, method = "BFGS") cat("Estimated mean :", round(fit$par["mean"], 3), "\n") cat("Estimated sd :", round(10^fit$par["log_sd"], 3), "\n") ## ----dnorm-parmoff------------------------------------------------------------ fitted_par <- fit$par # named c(mean=..., log_sd=...) # strip "log_" prefix → names become mean, sd # de-log "sd" → 10^log_sd ll_val <- ParmOff( function(mean, sd) sum(dnorm(data_obs, mean, sd, log = TRUE)), fitted_par, .strip = "log_", .logged = "sd" ) cat("Log-likelihood at fitted parameters:", round(ll_val, 2), "\n") ## ----galaxy-model------------------------------------------------------------- # Sérsic profile I(r) = I0 * exp(-b_n * ((r/Re)^(1/n) - 1)) # For simplicity we integrate along a 1-D radial profile sersic <- function(r, I0, Re, n) { bn <- 2 * n - 1/3 # approximation valid for n > 0.5 I0 * exp(-bn * ((r / Re)^(1 / n) - 1)) } # Exponential disc: I(r) = I0d * exp(-r / Rd) disc <- function(r, I0d, Rd) I0d * exp(-r / Rd) # Combined profile galaxy_profile <- function(r, I0, Re, n, I0d, Rd) { sersic(r, I0, Re, n) + disc(r, I0d, Rd) } # Simulate "observed" data set.seed(7) r_grid <- seq(0.1, 10, length.out = 60) true_par <- c(I0 = 100, Re = 2, n = 4, I0d = 50, Rd = 3) obs <- galaxy_profile(r_grid, I0 = 100, Re = 2, n = 4, I0d = 50, Rd = 3) * exp(rnorm(60, 0, 0.01)) # 1 % Gaussian scatter in log space # Chi-squared objective — parameters in "fit space": # log10(I0), log10(Re), n (linear), log10(I0d), log10(Rd) # Bounds: I0 in [1,1e4], Re in [0.1,20], n in [0.5,8], # I0d in [1,1e4], Rd in [0.1,20] fit_bounds_lower <- c(log10_I0 = 0, log10_Re = -1, n = 0.5, log10_I0d = 0, log10_Rd = -1) fit_bounds_upper <- c(log10_I0 = 4, log10_Re = log10(20), n = 8, log10_I0d = 4, log10_Rd = log10(20)) # Objective function — par is named in "log/linear fit space" chi_sq <- function(par, r, obs, lower, upper) { # Clamp to bounds (ParmOff does this too, but optim L-BFGS-B handles it here) # par <- pmax(pmin(par, upper[names(par)]), lower[names(par)]) # Use ParmOff to back-transform and forward to galaxy_profile model_val <- ParmOff( galaxy_profile, .args = as.list(par), .strip = "log10_", # log10_I0 -> I0, log10_Re -> Re, etc. .logged = c("I0", "Re", "I0d", "Rd"), # back-transform these four r = r # extra arg via ... ) sum((log(obs) - log(model_val))^2) } start_fit <- c(log10_I0 = 1.5, log10_Re = 0, n = 2, log10_I0d = 1.5, log10_Rd = 0.3) fit_gal <- optim( start_fit, chi_sq, r = r_grid, obs = obs, lower = fit_bounds_lower, upper = fit_bounds_upper, method = "L-BFGS-B" ) # Recover physical parameters recovered <- ParmOff( function(log10_I0, log10_Re, n, log10_I0d, log10_Rd) c(I0 = 10^log10_I0, Re = 10^log10_Re, n = n, I0d = 10^log10_I0d, Rd = 10^log10_Rd), as.list(fit_gal$par) ) cat("True :", paste(names(true_par), round(true_par, 2), sep = "=", collapse = ", "), "\n") cat("Fitted :", paste(names(recovered), round(recovered, 2), sep = "=", collapse = ", "), "\n") ## ----fig.width=8, fig.height=6------------------------------------------------ #We can then re-produce the best fit profile easily and compare: model_val <- ParmOff( galaxy_profile, .args = fit_gal$par, .strip = "log10_", # log10_I0 -> I0, log10_Re -> Re, etc. .logged = c("I0", "Re", "I0d", "Rd"), # back-transform these four r = r_grid # extra arg via ... ) magplot(r_grid, obs, type='l', log='xy', xlab='Rad', ylab='Intensity', lwd=5) lines(r_grid, model_val, col='lightgreen', lwd=3) ## ----return-args-------------------------------------------------------------- result <- ParmOff( model, list(x = 1, y = 2, z = 3, t = 99), .return = "args" ) cat("Arguments that WILL be passed:\n") print(result$current_args) cat("\nArguments that were IGNORED:\n") print(result$ignore_args) ## ----verbose-basic------------------------------------------------------------ model <- function(x, y, z) x * y + z # x is below its lower bound, y is above its upper bound # z is in log10 space and will be de-logged ParmOff(model, list(x = -1, y = 20, z = 1), .lower = list(x = 0), .upper = list(y = 10), .logged = "z", .verbose = TRUE) # Messages emitted (multi-line format): # Lower limit imposed on 'x' # before: -1 # after: 0 # Upper limit imposed on 'y' # before: 20 # after: 10 # ParmUnLog (10^x) applied to 'z' # before: 1 # after: 10 # Used arguments: # x y z # # Ignored arguments: # (blank — no ignored arguments) # Result: 0 * 10 + 10 = 10 ## ----verbose-helpers---------------------------------------------------------- # Lower limit ParmLimLo(list(alpha = -0.5, beta = 2), lower = list(alpha = 0), verbose = TRUE) # Upper limit ParmLimHi(list(sigma = 100, mu = 5), upper = list(sigma = 50), verbose = TRUE) # Forward log (e.g. before passing to an optimiser) ParmLog(list(amplitude = 500, index = 1.5), logged = "amplitude", verbose = TRUE) # De-log ParmUnLog(list(amplitude = 2.7, index = 1.5), logged = "amplitude", verbose = TRUE) ## ----perf, eval=FALSE--------------------------------------------------------- # arg_list <- list(x = 1, y = 2, z = 3) # # # Baseline: direct call # system.time(for (i in 1:1e4) model(1, 2, 3)) # # # do.call # system.time(for (i in 1:1e4) do.call(model, arg_list)) # # # ParmOff with full checking # system.time(for (i in 1:1e4) ParmOff(model, arg_list)) # # # ParmOff without checking (closer to do.call overhead) # system.time(for (i in 1:1e4) ParmOff(model, arg_list, .check = FALSE))