From 7e8dbd2236067ba88f9da7048ec2698336224393 Mon Sep 17 00:00:00 2001 From: Aaron Robotham Date: Wed, 2 Sep 2026 17:06:49 +0800 Subject: [PATCH 1/2] Better Z helper functions --- .github/workflows/main.yml | 3 +- DESCRIPTION | 4 +- NAMESPACE | 2 +- R/Highlander.R | 181 +++++++++++++++++++++- R/Lowlander.R | 48 +++++- R/Lowlander_Z.R | 131 ++++++++++++++++ R/Specs_help.R | 299 ++++++++++++++++++++++++++++++++---- man/Highlander.Rd | 54 ++++++- man/Lowlander.Rd | 13 +- man/Lowlander_helpers.Rd | 102 ++++++++++++ man/Specs_help.Rd | 117 +++++++++++++- vignettes/LD-Algorithms.Rmd | 103 ++++++++++++- 12 files changed, 999 insertions(+), 58 deletions(-) create mode 100644 R/Lowlander_Z.R create mode 100644 man/Lowlander_helpers.Rd diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index ac5e9b6..7ca64b9 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -53,7 +53,8 @@ jobs: any::knitr, any::rmarkdown, any::lhs, - any::ParmOff + any::ParmOff, + any::testthat - uses: r-lib/actions/check-r-package@v2 with: diff --git a/DESCRIPTION b/DESCRIPTION index 4120379..870a035 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -8,7 +8,8 @@ Maintainer: Aaron Robotham Description: There Can be Only One! Finds the highest points in log-likelihood space by switching between CMA and MCMC phases. Defaults should work pretty well. License: LGPL-3 Imports: LaplacesDemon, cmaeshpc, parallel, numDeriv, ParmOff -Suggests: knitr, magicaxis, ProSpect, celestial, foreach, rmarkdown, lhs +Suggests: knitr, magicaxis, ProSpect, celestial, foreach, rmarkdown, lhs, + testthat (>= 3.0.0) Remotes: asgr/cmaeshpc, asgr/LaplacesDemon, @@ -16,3 +17,4 @@ Remotes: Encoding: UTF-8 LazyData: true VignetteBuilder: knitr +Config/testthat/edition: 3 diff --git a/NAMESPACE b/NAMESPACE index 1dd7dda..fc127d4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,5 +4,5 @@ import("LaplacesDemon") import("cmaeshpc") import("ParmOff") -importFrom("stats", "quantile", "cor", "pnorm") +importFrom("stats", "quantile", "cor", "pnorm", "approx", "rnorm") importFrom("numDeriv", "hessian") diff --git a/R/Highlander.R b/R/Highlander.R index 8c8a59f..ebde9d9 100644 --- a/R/Highlander.R +++ b/R/Highlander.R @@ -1,5 +1,6 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, prior=NULL, Algorithm='CHARM', seed=666, lower=NULL, upper=NULL, + Lowlander=NULL, applyintervals=TRUE, updateintervals=FALSE, applyconstraints=TRUE, dynlim=2, ablim=0, optim_iters=2, Niters=c(100,100), NfinalMCMC=Niters[2], walltime = Inf, @@ -13,6 +14,33 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, date = date() call = match.call(expand.dots=TRUE) + # Capture what the user actually supplied *before* anything is forced. The + # Specs and LDargs defaults reference each other (LDargs embeds Specs), so + # evaluating one lazily evaluates the other and destroys the ability to tell + # "absent" from "supplied as-is". + parm_given = !missing(parm) + specs_given = !missing(Specs) + ldargs_given = !missing(LDargs) + + # Validate the Lowlander object up front, before the parallel dispatch below. + # Each worker wraps its own call in try(), so a malformed object would + # otherwise surface only as "every chain failed" rather than one clear error. + if(!is.null(Lowlander)){ + if(!is.list(Lowlander) || is.null(Lowlander$lower) || is.null(Lowlander$upper) || + is.null(Lowlander$latin_mod) || is.null(Lowlander$keep) || is.null(Lowlander$output)){ + stop("'Lowlander' must be the return value of a Lowlander() call (needs $lower, $upper, $latin_mod, $keep, $output).") + } + Npar_ll = length(Lowlander$lower) + arglens = list(parm=parm, lower=lower, upper=upper) + for (nm in names(arglens)) { + v = arglens[[nm]] + if(!is.null(v) && length(v) != Npar_ll){ + stop(paste0("'", nm, "' has length ", length(v), " but the Lowlander object has ", + Npar_ll, " parameters.")) + } + } + } + # Validate cores cores = as.integer(cores) if(is.na(cores) || cores < 1L){ @@ -32,12 +60,23 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, run_args = list( parm=parm, Data=Data, likefunc=likefunc, likefunctype=likefunctype, - liketype=liketype, prior=prior, lower=lower, upper=upper, applyintervals=applyintervals, + liketype=liketype, prior=prior, lower=lower, upper=upper, Lowlander=Lowlander, + applyintervals=applyintervals, updateintervals=updateintervals, applyconstraints=applyconstraints, dynlim=dynlim, ablim=ablim, optim_iters=optim_iters, Niters=Niters, NfinalMCMC=NfinalMCMC, - walltime=walltime, CMAargs=CMAargs, LDargs=LDargs, parm.names=parm.names, + walltime=walltime, CMAargs=CMAargs, parm.names=parm.names, keepall=keepall, cores=1L ) + # Algorithm must be forwarded explicitly now. Previously the whole default + # LDargs was passed down and carried Algorithm (and Specs) inside it, which + # is how the worker ended up running the requested algorithm; but forwarding + # an already-evaluated LDargs would mark Specs as user-supplied inside the + # worker and disable the Lowlander seeding, including the per-call DEMC + # resize. So forward the algorithm itself, and only those of Specs/LDargs + # the user actually wrote out. + run_args$Algorithm = Algorithm + if(specs_given){run_args$Specs = Specs} + if(ldargs_given){run_args$LDargs = LDargs} if(.Platform$OS.type == "windows"){ cl = parallel::makeCluster(cores) @@ -74,7 +113,7 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, #Make sure the mon.names will match up with what we do internally Data[['mon.names']] = c("LP", Data[['mon.names']][! Data[['mon.names']] == 'LP']) best_result$LD_last_comb = try( - LaplacesDemon::Combine(lapply(results,\(x) x$LD_last), Data=Data) + LaplacesDemon::Combine(lapply(results, function(x) x$LD_last), Data=Data) ) best_result$best_job = best_idx best_result$High_jobs = results @@ -100,6 +139,37 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, # optim_iters: number of CMA / LD loops # Niters: iters per CMA and LD (so can be vector length 2, otherwise value is repeated) + # Lowlander: optional result of a Lowlander() run, used to seed the parts of + # Specs that LaplacesDemon would otherwise guess blind (its GIV() generator + # probes from very wide Normals and returns NA on restricted likelihoods). + # It also supplies parm/lower/upper, but *only* where the user left them NULL. + if(!is.null(Lowlander)){ + if(!parm_given || is.null(parm)){ + # Prefer explicit parm.names, then the names Lowlander itself recorded, + # then Data$parm.names (Lowlander ignores the latter, so an object built + # without parm.names would otherwise hand back an unnamed start). + ll_names = if(!is.null(parm.names)) parm.names + else if(!is.null(colnames(Lowlander$latin_mod))) colnames(Lowlander$latin_mod) + else Data$parm.names + parm = Lowlander_best(Lowlander, parm.names = ll_names) + if(is.null(parm)){ + stop("The Lowlander object has no $best to start from.") + } + parm_given = TRUE + }else{ + # LaplacesDemon evaluates the model at Initial.Values before it considers + # Z at all, and stops outright with "The posterior is infinite!" if that + # is not finite. For AIES the reported thinned chain is walker 1's + # trajectory and Z[1,] is ignored, so an out-of-region parm can abort the + # run no matter how good the seeded Z is. + if(any(parm < Lowlander$lower | parm > Lowlander$upper)){ + message('Highlander: parm sits outside the Lowlander plausible region; LaplacesDemon evaluates the model there first and stops if the posterior is not finite, so this can abort the run regardless of the seeded Specs. Consider parm = Lowlander(...)$best.') + } + } + if(is.null(lower)){lower = Lowlander$lower} + if(is.null(upper)){upper = Lowlander$upper} + } + if(is.null(parm) & !is.null(lower) & !is.null(upper)){ parm = (lower + upper)/2 } @@ -154,8 +224,68 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, stop('lower and upper cannot have the same values!') } + # likefunctype is detected by calling likefunc once at parm and looking at the + # length of what comes back. Two problems with doing that raw: this is the + # first evaluation in the whole run, so any error surfaces with no hint that + # type detection is what tripped (and in multicore mode not one but `cores` + # identical failures appear); and Data has not had mon.names / parm.names / N + # filled in yet, which is exactly what an LD-style likefunc will be handed + # later, so a function can fail here and then work fine when fitted. Retry + # once against the augmented Data before giving up. if(is.null(likefunctype)){ - if(length(likefunc(parm, Data)) == 1){ + .probe = function(D) try(likefunc(parm, D), silent=TRUE) + + probe_out = .probe(Data) + probe_dat = 'Data exactly as supplied' + + if(inherits(probe_out, 'try-error')){ + DataProbe = Data + if(is.null(DataProbe[['mon.names']])){ + DataProbe[['mon.names']] = 'LP' + }else{ + DataProbe[['mon.names']] = c('LP', DataProbe[['mon.names']][!DataProbe[['mon.names']] == 'LP']) + } + if(is.null(DataProbe[['parm.names']])){ + if(!is.null(parm.names) && length(parm.names) == length(parm)){ + DataProbe[['parm.names']] = parm.names + }else{ + DataProbe[['parm.names']] = letters[seq_along(parm)] + } + } + if(is.null(DataProbe[['N']])){DataProbe[['N']] = 1} + + probe_out2 = .probe(DataProbe) + if(inherits(probe_out2, 'try-error')){ + stop(paste0( + "likefunc() failed when Highlander called it once at the start parm to work out whether it ", + "returns a CMA scalar or an LD list, so the fit has not started. Both the Data as supplied and ", + "Data with mon.names/parm.names/N filled in were tried. likefunc must be able to run at parm = ", + paste(if(is.null(names(parm))) sprintf('%g', parm) + else sprintf('%s=%g', names(parm), parm), collapse = ', '), + " (an out-of-region start is the usual cause: LaplacesDemon likewise evaluates the model at ", + "Initial.Values before it does anything else, so it would stop there too). If likefunc legitimately ", + "cannot be evaluated at this parm, set likefunctype explicitly ('CMA' for a single numeric return, ", + "'LD' for the list with LP/Dev/Monitor) to skip the probe. Underlying error: ", + conditionMessage(attr(probe_out2, 'condition')))) + } + probe_out = probe_out2 + probe_dat = 'Data with mon.names/parm.names/N filled in' + } + + probe_len = length(probe_out) + if(probe_len == 0){ + stop(paste0("likefunc(", probe_dat, ") returned nothing at the start parm, so Highlander cannot ", + "tell whether it is a CMA scalar or an LD list. It must return either a single numeric ", + "value or the LaplacesDemon list (LP, Dev, Monitor, ...); check that it returns a value ", + "on every path, rather than falling off the end or returning NULL.")) + } + if(probe_dat != 'Data exactly as supplied'){ + message("Highlander: likefunc failed with the Data as supplied but worked once mon.names/parm.names/N ", + "were present, so likefunctype = '", if(probe_len == 1) 'CMA' else 'LD', "' was detected from ", + "the second call. Pass likefunctype explicitly to avoid the extra probe.") + } + + if(probe_len == 1){ likefunctype = 'CMA' }else{ likefunctype = 'LD' @@ -228,6 +358,38 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, } } + # Rebuild the default LDargs by hand rather than forcing the promise: the + # signature's default embeds Specs, and evaluating it would print a second + # round of Specs_help guidance when we re-seed below. Keep the raw Data for + # that default -- DataLD has synthesized parm.names, which would change the + # auto-derived Nc (and epsilon length) versus the Specs promise default. + # When the user supplied LDargs we deliberately leave Specs absent, as before, + # so LaplacesDemon still falls back to its own per-algorithm defaults. + if(!ldargs_given){ + LDargs = list(control=list(abstol=0.1), Iterations=Niters[2], + Algorithm=Algorithm, Thinning=1) + LDargs[['Specs']] = if(specs_given) Specs else Specs_help(Algorithm, Data) + } + + # Normalize the pieces of LDargs that Highlander itself needs to know about. + # A user-supplied LDargs may omit any of these. + if(is.null(LDargs[['Iterations']])){LDargs[['Iterations']] = Niters[2]} + if(is.null(LDargs[['Thinning']])){LDargs[['Thinning']] = 1} + if(is.null(LDargs[['Algorithm']])){LDargs[['Algorithm']] = Algorithm} + + # Specs Highlander generated itself get seeded from the Lowlander samples. + # Anything the user passed by hand is left exactly as-is, including a Specs + # element embedded in their own LDargs. + user_spec = specs_given || (ldargs_given && !is.null(LDargs[['Specs']])) + seed_specs = (!user_spec && !is.null(Lowlander)) + preloop_iters = NULL + if(seed_specs){ + preloop_iters = LDargs[['Iterations']] + LDargs[['Specs']] = Specs_help(Algorithm, DataLD, Lowlander = Lowlander, + Iterations = LDargs[['Iterations']], + Thinning = LDargs[['Thinning']]) + } + parm_out = parm CMA_out = NULL @@ -340,6 +502,17 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, if(i == optim_iters){LDargs[['Iterations']] = NfinalMCMC} + # The run length of this particular LaplacesDemon call is only known here, + # and DEMC sizes Z as floor(Iterations/Thinning)+1, so Specs seeded above + # (which used the pre-loop Iterations) are rebuilt at the true length. + # Iterations/Thinning do not otherwise feed into Specs, so this is a no-op + # in shape for every other algorithm. + if(seed_specs && !identical(LDargs[['Iterations']], preloop_iters)){ + LDargs[['Specs']] = suppressMessages(Specs_help(Algorithm, DataLD, Lowlander = Lowlander, + Iterations = LDargs[['Iterations']], + Thinning = LDargs[['Thinning']])) + } + if(LDargs[['Iterations']] > 0){ LD_out = do.call('LaplacesDemon', c(list(Model=LDfunc, Data=quote(DataLD), Initial.Values=parm_out), LDargs)) diff --git a/R/Lowlander.R b/R/Lowlander.R index 9ad03c4..35320a1 100644 --- a/R/Lowlander.R +++ b/R/Lowlander.R @@ -78,21 +78,54 @@ Lowlander = function(lower, upper, Data, likefunc, Nsamp = 'auto', pcut = 0.1, latin_mod[, j] = lower[j] + latin[, j] * (upper[j] - lower[j]) } - # Evaluate likefunc in parallel (Windows-safe) + # Evaluate likefunc in parallel (Windows-safe). The results are kept as a + # list rather than unlisted straight away, so the shape of what each + # evaluation returned can be checked first (below). if (.Platform$OS.type == "windows" && ncores > 1L) { # PSOCK cluster for Windows multi-core cl = parallel::makeCluster(ncores) on.exit(parallel::stopCluster(cl), add = TRUE) parallel::clusterExport(cl, varlist = c("latin_mod", "Data", "likefunc"), envir = environment()) - output = unlist(parallel::parLapply(cl, seq_len(Nsamp), - function(i) likefunc(latin_mod[i, ], Data))) + raw = parallel::parLapply(cl, seq_len(Nsamp), + function(i) likefunc(latin_mod[i, ], Data)) } else { - output = unlist(parallel::mclapply(seq_len(Nsamp), - function(i) likefunc(latin_mod[i, ], Data), - mc.cores = ncores)) + raw = parallel::mclapply(seq_len(Nsamp), + function(i) likefunc(latin_mod[i, ], Data), + mc.cores = ncores) } + # `output` has to be one value per row of latin_mod, because `keep` is an + # index into latin_mod's rows. unlist() quietly breaks that alignment: a + # length-0 result (NULL) shifts every later sample up a row, and a length-4 + # result (the LaplacesDemon-style list) leaves `output` Nsamp*4 long, which + # surfaces much later as an opaque subscript error. Reject both here, where + # the offending sample numbers are still known. + # Name at most a few offending rows; a likefunc that is wrong for one sample + # is almost always wrong for all of them, and a full index is just noise. + .first = function(idx, n = 5L) idx[seq_len(min(n, length(idx)))] + .samples = function(idx) paste0(if (length(idx) > 1L) "samples " else "sample ", + paste(.first(idx), collapse = ", ")) + + len = vapply(raw, length, integer(1)) + badlen = which(len != 1L) + if (length(badlen) > 0L) { + stop(paste0( + "'likefunc' must return a single numeric value for each sample, but it ", + "returned ", paste(unique(len[badlen]), collapse = "/"), " value(s) for ", + length(badlen), " of ", Nsamp, " samples (e.g. ", .samples(badlen), + "). This is usually the LaplacesDemon-style list output (LP, Dev, ", + "Monitor, ...); Lowlander needs one scalar per sample, so wrap it and ", + "pick the element you want, e.g. likefunc = function(parm, Data) ", + "mylikefunc(parm, Data)$LP.")) + } + badtype = which(!vapply(raw, is.numeric, logical(1))) + if (length(badtype) > 0L) { + stop(paste0("'likefunc' must return numeric values, but ", .samples(badtype), + " returned ", class(raw[[badtype[1]]])[1], ".")) + } + output = unlist(raw) + # Determine which samples to keep based on quantile cut good = is.finite(output) output_good = output[good] @@ -166,6 +199,7 @@ Lowlander = function(lower, upper, Data, likefunc, Nsamp = 'auto', pcut = 0.1, cutval = cutval, cor_data = cor_data, latin = latin, - latin_mod = latin_mod + latin_mod = latin_mod, + liketype = liketype ))) } diff --git a/R/Lowlander_Z.R b/R/Lowlander_Z.R new file mode 100644 index 0000000..65de9ee --- /dev/null +++ b/R/Lowlander_Z.R @@ -0,0 +1,131 @@ +Lowlander_kept = function(Lowlander = NULL, liketype = NULL, parm.names = NULL) { + # The kept (plausible) parameter samples from a Lowlander() return, ranked + # best-first. Returns NULL with a warning when the object is unusable. + # + # Lowlander restricts `keep` to finite-likelihood samples past its quantile + # cut, so these rows are in-bounds by construction. + + if (is.null(Lowlander)) { + return(NULL) + } + + if (!is.list(Lowlander) || is.null(Lowlander$latin_mod) || + is.null(Lowlander$keep) || is.null(Lowlander$output)) { + warning('Lowlander object is missing latin_mod/keep/output.') + return(NULL) + } + + latin_mod = as.matrix(Lowlander$latin_mod) + keep = as.integer(Lowlander$keep) + output = as.numeric(Lowlander$output) + + if (length(keep) == 0 || any(keep < 1L | keep > nrow(latin_mod))) { + warning('Lowlander object has no valid kept samples.') + return(NULL) + } + + # Defensive: a hand-edited or older object should not be able to smuggle a + # non-finite row through. + outk = output[keep] + ok = is.finite(outk) + if (!any(ok)) { + warning('All kept Lowlander samples have non-finite likelihood.') + return(NULL) + } + kept = latin_mod[keep[ok], , drop = FALSE] + outk = outk[ok] + + # Prefer the direction Lowlander actually used over anything passed in here. + lt = Lowlander$liketype + if (is.null(lt)) lt = liketype + if (is.null(lt)) lt = 'min' + + ord = if (identical(lt, 'max')) order(outk, decreasing = TRUE) else order(outk) + kept = kept[ord, , drop = FALSE] + + if (!is.null(parm.names) && length(parm.names) == ncol(kept)) { + colnames(kept) = parm.names + } + + return(kept) +} + +Lowlander_Z = function(Lowlander = NULL, N = 1, liketype = NULL, parm.names = NULL) { + # A design matrix of N plausible parameter samples drawn from a Lowlander() + # return, to seed ensemble samplers (AIES Z) or anything else that needs + # several starting points inside the sensible region. Returns NULL when + # unavailable so callers can fall back to their own defaults. + # + # Rows sit at evenly-spaced quantile ranks of the kept samples, ranked + # best-first. This guarantees distinct walkers and covers the plausible + # region rather than piling up on the mode; AIES proposes + # theta' = theta_s + z (theta_i - theta_s), so walkers that all sit on top of + # each other can only ever propose where they already are. Choosing a + # different rank spacing (e.g. denser near the best point) performed + # similarly in testing, so nothing here should be read as an optimised + # ranking - the benefit comes from using Lowlander's kept samples at all + # rather than leaving Z = NULL. + + kept = Lowlander_kept(Lowlander, liketype = liketype, parm.names = parm.names) + if (is.null(kept)) { + return(NULL) + } + + # One point cannot support an ensemble: the stretch move divides by the + # difference between two distinct walkers. + nk = nrow(kept) + if (nk < 2L) { + warning('Fewer than 2 usable Lowlander samples; cannot build a spread Z.') + return(NULL) + } + + N = max(1L, as.integer(N)) + + if (N <= nk) { + idx = floor(approx(seq(0, 1, length.out = nk), seq_len(nk), + xout = (seq_len(N) - 0.5) / N)$y) + idx = pmin(pmax(idx, 1L), nk) + Z = kept[idx, , drop = FALSE] + } else { + # More walkers than samples: recycle the kept set, then jitter the repeats + # (exactly-duplicated rows degenerate the stretch move for any pair sharing + # a position). Jitter is clamped back inside the plausible box. + reps = rep(seq_len(nk), length.out = N) + Z = kept[reps, , drop = FALSE] + dup = which(duplicated(reps)) + if (length(dup) > 0) { + span = apply(kept, 2, max) - apply(kept, 2, min) + flat = !is.finite(span) | span <= 0 + span[flat] = pmax(abs(kept[1, ][flat]), 1) + jit = matrix(rnorm(length(dup) * ncol(Z)), nrow = length(dup)) * + matrix(rep(span * 1e-3, each = length(dup)), nrow = length(dup)) + Z[dup, ] = Z[dup, , drop = FALSE] + jit + lo = apply(kept, 2, min) + hi = apply(kept, 2, max) + Z = matrix(pmin(pmax(as.vector(Z), rep(lo, each = nrow(Z))), + rep(hi, each = nrow(Z))), nrow = nrow(Z)) + } + } + + Z = matrix(as.numeric(Z), nrow = N) + if (!is.null(colnames(kept))) colnames(Z) = colnames(kept) + return(Z) +} + +Lowlander_best = function(Lowlander = NULL, parm.names = NULL) { + # The single most plausible point from a Lowlander() return as a plain + # numeric vector, for seeding Initial.Values / twalk SIV / IM mu. + if (is.null(Lowlander)) { + return(NULL) + } + best = Lowlander$best + if (is.null(best)) { + warning('Lowlander object has no $best element.') + return(NULL) + } + best = as.numeric(unname(best)) + if (!is.null(parm.names) && length(parm.names) == length(best)) { + names(best) = parm.names + } + return(best) +} diff --git a/R/Specs_help.R b/R/Specs_help.R index fdb7a14..e711649 100644 --- a/R/Specs_help.R +++ b/R/Specs_help.R @@ -1,11 +1,101 @@ -Specs_help = function(Algorithm = 'CHARM', Data = NULL, ...) { +Specs_help = function(Algorithm = 'CHARM', Data = NULL, Lowlander = NULL, + Iterations = NULL, Thinning = 1, ...) { # Returns the default Specs list for the given LaplacesDemon Algorithm name. # For algorithms that require user-supplied values (no built-in defaults), # a template list is returned with placeholder values and a message is # printed to guide the user. + # + # If a Lowlander() return object is supplied it is used to fill in the + # parts that LaplacesDemon would otherwise guess blind (via GIV), which is + # the usual cause of out-of-bounds failures on restricted likelihoods. dots = list(...) + # Number of parameters, from whichever source is available. + Npar = NULL + if (!is.null(Data$parm.names)) { + Npar = length(Data$parm.names) + } else if (!is.null(Lowlander$latin_mod)) { + Npar = ncol(as.matrix(Lowlander$latin_mod)) + } + if (!is.null(Npar) && Npar < 1L) { + stop('Cannot determine the number of parameters from Data or Lowlander.') + } + + .even_Nc = function(Nc, CPUs = 1) { + Nc = max(3L, as.integer(Nc)) + if (CPUs > 1 && Nc %% 2 != 0) { + Nc = Nc + 1L + message(Algorithm, ': Nc rounded up to ', Nc, ' (must be even when CPUs > 1).') + } + Nc + } + .Nc_ens = function(CPUs = 1) { + .even_Nc(if (is.null(Npar)) 3L else max(4L, 2L * Npar), CPUs) + } + + # LaplacesDemon rewrites Iterations/Thinning before it ever looks at Specs: + # Iterations becomes round(abs(Iterations)) clamped up to a minimum of 11, and + # Thinning becomes round(abs(Thinning)) but falls back to 1 unless it sits in + # [1, Iterations]. DEMC then sizes Z as floor(Iterations/Thinning)+1 rows, so + # an array built from the *raw* values can mismatch and be hard-stopped with + # "The first dimension of Z is incorrect". Mirror the same rules here. + .ld_iters = function(Iters) { + Iters = round(abs(Iters)) + max(11L, as.integer(Iters)) + } + .ld_thin = function(Thin, Iters_adj) { + Thin = round(abs(Thin)) + if (is.na(Thin) || Thin < 1L || Thin > Iters_adj) 1L else as.integer(Thin) + } + .z_rows = function(Iters, Thin) { + Iters = .ld_iters(Iters) + floor(Iters / .ld_thin(Thin, Iters)) + 1L + } + # Non-finite Iterations cannot size anything (as.integer(Inf) is NA), so the + # callers below must fall back to Z = NULL rather than build a broken array. + .sizeable = !is.null(Iterations) && is.finite(Iterations) && !is.null(Npar) + + # Plausible starting points mined from the Lowlander run, if we have one. + # With no Lowlander object these stay NULL and LaplacesDemon falls back to + # GIV(), its own blind guesser; that is the pre-existing behaviour and is not + # worth warning about. + .Zmatrix = function(Nc) { + if (is.null(Lowlander)) return(NULL) + Z = Lowlander_Z(Lowlander, N = Nc, parm.names = Data$parm.names) + if (is.null(Z)) { + message(paste0(Algorithm, ': no usable Z could be built from the Lowlander ', + 'object; LaplacesDemon will fall back to GIV(), which often ', + 'fails on out-of-bounds proposals.')) + } + Z + } + .seed = function(what) { + if (is.null(Lowlander)) return(NULL) + b = Lowlander_best(Lowlander, parm.names = Data$parm.names) + if (is.null(b)) { + message(paste0(Algorithm, ': ', what, ' not auto-filled from the Lowlander object.')) + } + b + } + # Per-parameter evaluation grids for the Griddy-Gibbs family. Grid holds + # OFFSETS from the current value (LaplacesDemon evaluates at + # prop[j] + Grid[[j]]), so the built-in default of seq(-0.1, 0.1) is a fixed + # +-0.1 probe that is far too coarse for parameters spanning several decades + # and too fine for tightly-bounded ones. When a Lowlander object gives us the + # plausible width of each parameter, scale each grid to that width. + .grid = function(n = 7, frac = 0.5) { + if (is.null(Npar)) return(NULL) + lo = Lowlander$lower; hi = Lowlander$upper + if (is.null(lo) || is.null(hi) || length(lo) != Npar || length(hi) != Npar || + any(!is.finite(lo)) || any(!is.finite(hi)) || any(lo >= hi)) { + return(NULL) + } + n = max(3L, as.integer(n)) + half = frac * (hi - lo) / 2 + lapply(seq_len(Npar), function(j) seq(-half[j], half[j], length.out = n)) + } + if (Algorithm == 'ADMG') { # Adaptive Directional Metropolis-within-Gibbs Specs_list = list(n = 0, Periodicity = 1) @@ -16,25 +106,38 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, ...) { } else if (Algorithm == 'AGG') { # Adaptive Griddy-Gibbs (Specs required – no default) - message('AGG requires user-supplied Specs. Returning template.') - Specs_list = list(Grid = NULL, dparm = NULL, smax = 0.1, CPUs = 1, + Grid = .grid() + if (is.null(Grid)) { + message('AGG requires a per-parameter Grid; pass a Lowlander object to auto-build one. Returning template.') + Grid = NULL + } + # dparm must stay NULL: it names the *discrete* parameters, and passing 0 + # makes LaplacesDemon clamp it to a real index, silently treating one + # continuous parameter as discrete. NULL is the "none" sentinel. + Specs_list = list(Grid = Grid, dparm = NULL, smax = 0.1, CPUs = 1, Packages = NULL, Dyn.libs = NULL) } else if (Algorithm == 'AHMC') { # Adaptive Hamiltonian Monte Carlo - # epsilon and m depend on the number of parameters; placeholders shown here - message('AHMC: epsilon and m length must equal length(Initial.Values).') - Specs_list = list(epsilon = 0.1, L = 2, m = NULL, Periodicity = 100) + # epsilon must be length(Initial.Values); LaplacesDemon recycles a scalar, + # so we build the per-parameter vector it is asking for. + npar = if (is.null(Npar)) 1L else Npar + message('AHMC: epsilon length must equal length(Initial.Values).') + # m = NULL makes LaplacesDemon use an identity mass matrix; a plain vector + # is rejected with "x must be a square matrix". + Specs_list = list(epsilon = rep(1 / npar, npar), L = 2, m = NULL, + Periodicity = 1) } else if (Algorithm == 'AIES') { # Affine-Invariant Ensemble Sampler (Specs required – no default) - if(is.null(Data$parm.names)){ + Nc = if (is.null(dots$Nc)) .Nc_ens(if (is.null(dots$CPUs)) 1 else dots$CPUs) + else .even_Nc(dots$Nc, if (is.null(dots$CPUs)) 1 else dots$CPUs) + if (is.null(dots$Nc) && is.null(Npar)) { message('AIES requires user-supplied Nc or Data. Returning template default Nc = 3.') - Nc = 3 - }else{ - Nc = 2*length(Data$parm.names) + } else if (is.null(dots$Nc)) { + message('AIES: using Nc = 2 * length(parm.names) = ', Nc, '.') } - Specs_list = list(Nc = Nc, Z = NULL, beta = 2, CPUs = 1, + Specs_list = list(Nc = Nc, Z = .Zmatrix(Nc), beta = 2, CPUs = 1, Packages = NULL, Dyn.libs = NULL) } else if (Algorithm == 'AM') { @@ -58,7 +161,41 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, ...) { } } else if (Algorithm == 'DEMC') { # Differential Evolution Markov Chain (Specs required – no default) - Specs_list = list(Nc = 3, Z = NULL, gamma = NULL, w = 0.1) + # LaplacesDemon wants Z as an array: floor(Iterations/Thinning)+1 x + # length(parm) x Nc. Without Iterations we cannot size it, so leave Z + # NULL and let LaplacesDemon generate it (via GIV). + Nc = .Nc_ens(1) + if (!is.null(dots$Nc)) Nc = max(3L, as.integer(dots$Nc)) + Z = NULL + if (is.null(Iterations) || !is.finite(Iterations)) { + message("DEMC: pass a finite Iterations (and Thinning) to auto-build Z; ", + "leaving Z = NULL so LaplacesDemon falls back to GIV().") + } else if (is.null(Npar)) { + message('DEMC: cannot determine the parameter count; leaving Z = NULL.') + } else { + pool = Lowlander_kept(Lowlander, parm.names = Data$parm.names) + if (is.null(Lowlander)) { + message('DEMC: pass a Lowlander object to seed Z with in-bounds samples; ', + 'leaving Z = NULL so LaplacesDemon falls back to GIV().') + } else if (is.null(pool)) { + message('DEMC: no usable Lowlander samples; Z = NULL, LaplacesDemon ', + 'will fall back to GIV().') + } else { + # Each chain gets its own history of plausible samples. LaplacesDemon + # proposes gamma * (Z[r1,,s1] - Z[r2,,s2]) from these rows, so they all + # need to be usable. Note it first overwrites Z[1,,1] with the value + # returned at Initial.Values (so that one slot is ours to fill but not + # to keep, exactly as AIES ignores Z[1,]), then evaluates Model() at + # Z[1,,i] for i = 2..Nc with no try() wrapper -- an error there aborts + # the run outright, so those rows must be in-bounds. + nrow_z = .z_rows(Iterations, Thinning) + Z = array(0, dim = c(nrow_z, Npar, Nc)) + for (i in seq_len(dim(Z)[3])) { + Z[, , i] = pool[sample.int(nrow(pool), nrow_z, replace = TRUE), , drop = FALSE] + } + } + } + Specs_list = list(Nc = Nc, Z = Z, gamma = NULL, w = 0.1) } else if (Algorithm == 'DRAM') { # Delayed Rejection Adaptive Metropolis @@ -74,8 +211,12 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, ...) { } else if (Algorithm == 'GG') { # Griddy-Gibbs (Specs required – no default) - message('GG requires user-supplied Specs. Returning template.') - Specs_list = list(Grid = NULL, dparm = NULL, CPUs = 1, + Grid = .grid() + if (is.null(Grid)) { + message('GG requires a per-parameter Grid; pass a Lowlander object to auto-build one. Returning template.') + Grid = NULL + } + Specs_list = list(Grid = Grid, dparm = NULL, CPUs = 1, Packages = NULL, Dyn.libs = NULL) } else if (Algorithm == 'Gibbs') { @@ -92,20 +233,29 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, ...) { } } else if (Algorithm == 'HMC') { # Hamiltonian Monte Carlo - # epsilon and m depend on the number of parameters; placeholders shown here - message('HMC: epsilon and m length must equal length(Initial.Values).') - Specs_list = list(epsilon = 0.1, L = 2, m = NULL) + npar = if (is.null(Npar)) 1L else Npar + message('HMC: epsilon length must equal length(Initial.Values).') + Specs_list = list(epsilon = rep(1 / npar, npar), L = 2, m = NULL) } else if (Algorithm == 'HMCDA') { # Hamiltonian Monte Carlo with Dual-Averaging (Specs required – no default) - message('HMCDA requires user-supplied Specs. Returning template.') - Specs_list = list(A = 1000, delta = 0.65, epsilon = NULL, Lmax = 1000, - lambda = 0.1) + # A NULL epsilon lets dual-averaging tune the step size adaptively. + message('HMCDA: A must be well below Iterations; epsilon = NULL lets dual-averaging tune it.') + Specs_list = list(A = max(1L, if (is.null(Iterations)) floor(1000) else floor(Iterations / 2)), + delta = 0.65, epsilon = NULL, Lmax = 1000, lambda = 0.1) } else if (Algorithm == 'IM') { # Independence Metropolis (Specs required – no default) - message('IM requires user-supplied Specs. mu must equal Initial.Values in length.') - Specs_list = list(mu = NULL) + # mu is the mean of the proposal distribution and must be exactly + # length(Initial.Values) or LaplacesDemon stops. Left NULL rather than + # guessed when there is no Lowlander object: an arbitrary mu (e.g. all + # zeros) still "runs" but silently biases where proposals land, which is + # worse than LaplacesDemon's own clear length error. + mu = .seed('mu') + if (is.null(mu)) { + message('IM requires user-supplied mu of length equal to Initial.Values (or pass a Lowlander object).') + } + Specs_list = list(mu = mu) } else if (Algorithm == 'INCA') { # Interchain Adaptation @@ -198,8 +348,19 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, ...) { Specs_list = list(epsilon = 0.1, L = 2, m = NULL, Temperature = 1) } else if (Algorithm == 'twalk') { - # t-walk - Specs_list = list(SIV = NULL, n1 = 4, at = 6, aw = 1.5) + # t-walk. LaplacesDemon stops if SIV gives a non-finite posterior, or if + # SIV equals Initial.Values after the model update. + SIV = NULL + if (!is.null(Lowlander)) { + # Take a plausible point that is *not* the best one, to reduce the chance + # of colliding with an Initial.Values that is already near the mode. + Zt = Lowlander_Z(Lowlander, N = 2, parm.names = Data$parm.names) + if (!is.null(Zt)) SIV = as.numeric(Zt[nrow(Zt), ]) + } + if (is.null(SIV)) { + message("twalk: pass a Lowlander object to auto-fill SIV; it must give a finite posterior and must differ from Initial.Values.") + } + Specs_list = list(SIV = SIV, n1 = 4, at = 6, aw = 1.5) } else if (Algorithm == 'UESS') { # Univariate Eigenvector Slice Sampler @@ -219,8 +380,94 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, ...) { stop(paste0('Unknown Algorithm: "', Algorithm, '". See LaplacesDemon documentation for supported algorithm names.')) } - if(any(names(dots) %in% names(Specs_list))){ - Specs_list[names(dots)] = dots + if (!is.null(Specs_list) && length(dots) > 0) { + hit = names(dots) %in% names(Specs_list) + Specs_list[names(dots)[hit]] = dots[hit] + # Names that match no known element of this algorithm's Specs are silently + # dropped by LaplacesDemon, so surface them here rather than let a typo + # (e.g. NC = 8) pass unnoticed. + if (any(!hit)) { + warning(paste0("Specs_help: ignoring unknown argument(s) for ", Algorithm, ": ", + paste(names(dots)[!hit], collapse = ", "), + ". Valid names are: ", paste(names(Specs_list), collapse = ", "), ".")) + } + } else if (is.null(Specs_list) && length(dots) > 0) { + # DRM/RDMH take no Specs; a bare NULL keeps LaplacesDemon's own default path. + warning(paste0('Specs_help: ', Algorithm, ' takes no Specs; the following ', + 'argument(s) are ignored: ', paste(names(dots), collapse = ', '), '.')) + } + + # User overrides above can desynchronise Z from Nc (LaplacesDemon stops on a + # mismatch), so repair the pair here. + if (Algorithm %in% c('AIES', 'DEMC') && !is.null(Specs_list$Z)) { + Nc = max(3L, abs(round(Specs_list$Nc))) + CPUs = if (is.null(Specs_list$CPUs)) 1 else Specs_list$CPUs + if (Algorithm == 'AIES' && CPUs > 1 && Nc %% 2 != 0) { + Nc = Nc + 1L + message('AIES: Nc rounded up to ', Nc, ' (must be even when CPUs > 1).') + } + Specs_list$Nc = Nc + if (Algorithm == 'AIES') { + Z = Specs_list$Z + if (is.matrix(Z) && nrow(Z) != Nc) { + # Only rebuild from Lowlander when Z is ours; a user-supplied Z is + # tiled/trimmed to match Nc rather than discarded. + if ('Z' %in% names(dots)) { + Specs_list$Z = matrix(as.numeric(Z[rep(seq_len(nrow(Z)), length.out = Nc), , drop = FALSE]), + nrow = Nc, ncol = ncol(Z)) + } else { + Zn = Lowlander_Z(Lowlander, N = Nc, parm.names = Data$parm.names) + if (!is.null(Zn)) Specs_list$Z = Zn + } + } + } else if (!is.array(Specs_list$Z) || length(dim(Specs_list$Z)) != 3 || + dim(Specs_list$Z)[3] != Nc || + # identical(), not !=: Npar can be NULL here (no Data, no + # Lowlander), and `dim(Z)[2] != NULL` is logical(0), which turns + # the whole condition into NA and errors inside if(). + !identical(dim(Specs_list$Z)[2], Npar) || + (.sizeable && dim(Specs_list$Z)[1] != .z_rows(Iterations, Thinning))) { + # DEMC wants a 3-D array: floor(Iterations/Thinning)+1 x Npar x Nc, with + # Iterations/Thinning as LaplacesDemon *adjusts* them. A bare matrix + # silently hits LaplacesDemon's rbind-recycle and dies with "subscript out + # of bounds", so rebuild from whatever we have. + Z = Specs_list$Z + Zn = NULL + if (.sizeable) { + nrow_z = .z_rows(Iterations, Thinning) + src = NULL + if ('Z' %in% names(dots)) { + if (is.matrix(Z) && ncol(Z) == Npar) src = Z + if (is.null(src) && is.array(Z) && length(dim(Z)) == 3 && + dim(Z)[2] == Npar) src = Z[, , 1] + } else { + src = Lowlander_kept(Lowlander, parm.names = Data$parm.names) + } + if (!is.null(src)) { + Zn = array(0, dim = c(nrow_z, Npar, Nc)) + for (i in seq_len(Nc)) { + Zn[, , i] = if ('Z' %in% names(dots)) { + src[rep(seq_len(nrow(src)), length.out = nrow_z), , drop = FALSE] + } else { + # Draw independently per chain: DEMC proposes gamma * + # (Z[r1,,s1] - Z[r2,,s2]), so chains sharing an identical history + # produce degenerate (zero) difference vectors. + src[sample.int(nrow(src), nrow_z, replace = TRUE), , drop = FALSE] + } + } + } + } + if (is.null(Zn)) { + warning(paste0('DEMC: Z must be a 3-D array of dim ', + if (.sizeable) .z_rows(Iterations, Thinning) else 'floor(Iterations/Thinning)+1', + ' x ', Npar, ' x ', Nc, '. Set to NULL so LaplacesDemon ', + 'regenerates it. Pass Iterations and Thinning to Specs_help ', + 'to have the array built for you.')) + Specs_list$Z = NULL + } else { + Specs_list$Z = Zn + } + } } return(Specs_list) diff --git a/man/Highlander.Rd b/man/Highlander.Rd index 858159e..b5a08ef 100644 --- a/man/Highlander.Rd +++ b/man/Highlander.Rd @@ -10,6 +10,7 @@ There Can be Only One! Finds the highest points in log-likelihood space by switc \usage{ Highlander(parm = NULL, Data, likefunc, likefunctype = NULL, liketype = NULL, prior = NULL, Algorithm = "CHARM", seed = 666, lower = NULL, upper = NULL, + Lowlander = NULL, applyintervals = TRUE, updateintervals = FALSE, applyconstraints = TRUE, dynlim = 2, ablim = 0, optim_iters = 2, Niters = c(100,100), NfinalMCMC = Niters[2], walltime = Inf, Specs = Specs_help(Algorithm, Data), @@ -22,7 +23,7 @@ Highlander(parm = NULL, Data, likefunc, likefunctype = NULL, liketype = NULL, %- maybe also 'usage' for other objects documented here. \arguments{ \item{parm}{ -Numeric vector; usual parameter vector- input to \option{likefunc} +Numeric vector; usual parameter vector- input to \option{likefunc}. May be left as \option{NULL} if either \option{lower} and \option{upper} are given (the midpoint is used) or a \option{Lowlander} object is supplied (\code{Lowlander$best} is used). } \item{Data}{ List; usual data input to \option{likefunc} This is where you should put all additional arguments needs since \code{\link[LaplacesDemon]{LaplacesDemon}} cannot be passed additional likelihood function arguments. @@ -52,10 +53,13 @@ THMC, twalk, UESS, USAMWG, USMWG. Integer scalar/vector; random seed to start the \code{Highlander} function with. If using multiple cores and the length of \option{seed} equals the number of cores, then those values will be used, else the seeds will become \option{seed} + 0:(\option{cores} - 1), i.e. the first job will use the \option{seed} provided, and the second job becomes \option{seed} + 1 etc. } \item{lower}{ -Numeric vector; lower limit allowed for \option{parm}; if \option{lower} = NULL and Data$intervals$lo = NULL it becomes the minimum of \option{parm}*(1/\option{dynlim}) - \option{ablim} or \option{parm}*\option{dynlim} - \option{ablim}; else if Data$intervals$lo is not NULL then this will be used instead. +Numeric vector; lower limit allowed for \option{parm}; if \option{lower} = NULL and a \option{Lowlander} object is supplied then \code{Lowlander$lower} is used; else if Data$intervals$lo = NULL it becomes the minimum of \option{parm}*(1/\option{dynlim}) - \option{ablim} or \option{parm}*\option{dynlim} - \option{ablim}; else if Data$intervals$lo is not NULL then this will be used instead. } \item{upper}{ -Numeric vector; upper limit allowed for \option{parm}; if \option{upper} = NULL and Data$intervals$hi = NULL it becomes the maximum of \option{parm}*(1/\option{dynlim}) + \option{ablim} or \option{parm}*\option{dynlim} + \option{ablim}; else if Data$intervals$hi is not NULL then this will be used instead. +Numeric vector; upper limit allowed for \option{parm}; if \option{upper} = NULL and a \option{Lowlander} object is supplied then \code{Lowlander$upper} is used; else if \option{upper} = NULL and Data$intervals$hi = NULL it becomes the maximum of \option{parm}*(1/\option{dynlim}) + \option{ablim} or \option{parm}*\option{dynlim} + \option{ablim}; else if Data$intervals$hi is not NULL then this will be used instead. +} + \item{Lowlander}{ +List; optional return value of a \code{\link{Lowlander}} call, used to seed the parts of \option{Specs} that \code{\link[LaplacesDemon]{LaplacesDemon}} would otherwise guess blind. When supplied, and only when \option{Specs} (and \option{LDargs$Specs}) have been left at their defaults, they are built with \code{Specs_help(Algorithm, Data, Lowlander = Lowlander, Iterations = ..., Thinning = ...)} so that ensemble and grid-based algorithms (\option{AIES}, \option{DEMC}, \option{twalk}, \option{IM}, \option{GG}, \option{AGG}) start from Lowlander's kept in-bounds samples rather than from \code{LaplacesDemon}'s \code{GIV()} generator, which probes from very wide Normal distributions and typically returns \code{NA} for likelihoods that are undefined outside a narrow region. \option{Lowlander} additionally supplies \option{parm}, \option{lower} and \option{upper}, but only where each has been left as \option{NULL}; anything passed explicitly is used as-is, so the object never overrides a user decision. For \option{parm} this means \code{Lowlander$best} is used as the starting point, which is generally what you want: \code{LaplacesDemon} evaluates the model at \option{parm} before it looks at \option{Specs$Z} at all, and for \option{AIES} \code{Z[1,]} is ignored in favour of \option{parm} (the reported thinned chain is walker 1's trajectory). If \option{parm} is outside the finite region the run can abort no matter how good the seeded \option{Z} is. Supplying \option{parm} explicitly from outside the Lowlander region produces a message rather than an override. Anything the user passes as \option{Specs} or embedded in their own \option{LDargs} is never rewritten, so a hand-built \option{Z} remains exactly that (and for \option{DEMC} must be sized for the run length itself; see \code{\link{Specs_help}}). Note the Lowlander object is validated before any parallel dispatch, so a malformed one fails with a single clear error rather than once per worker. Also note that a Lowlander run selects its samples from \option{likefunc} alone; if a \option{prior} is in play it is not reflected in the seeded region unless \option{likefunc} already incorporates it, so treat the seeded values as in-bounds starting points rather than as draws from the posterior. } \item{applyintervals}{ Logical; should the lower and upper intervals be applied by \code{Highlander} directly? This means it will create the list items Data$intervals$lo = lower and Data$intervals$hi = upper. You might want to set this to FALSE when your likelihood function deals with intervals itself using a different structure to the Data$intervals and creating such structures will break functionality inside \option{likefunc} (e.g. this is the case when running with \code{ProFit}). @@ -85,13 +89,13 @@ Integer scalar; number of iterations to run for the final MCMC run. You might wa Numeric scalar; the maximum allowed CPU time for \code{Highlander} in minutes. The function will stop with the best solution to date once the walltime has been exceeded (although it does not stop mid CMA or MCMC, so the \option{walltime} will usually be exceeded a bit). } \item{Specs}{ -List; \option{Specs} to be passed into \option{LDargs} as a convenience. \option{Specs} is by default populated with the settings obtained from \code{\link{Specs_help}}, if non-defaults are desired then extra arguments can be passed into \code{\link{Specs_help}} via \dots or a fully user-defined list can be constructed obeying the requirements of \code{\link[LaplacesDemon]{LaplacesDemon}} for the \option{Algorithm} being used. +List; \option{Specs} to be passed into \option{LDargs} as a convenience. \option{Specs} is by default populated with the settings obtained from \code{\link{Specs_help}}, if non-defaults are desired then extra arguments can be passed into \code{\link{Specs_help}} via \dots or a fully user-defined list can be constructed obeying the requirements of \code{\link[LaplacesDemon]{LaplacesDemon}} for the \option{Algorithm} being used. For algorithms that need several plausible starting points (\option{AIES}, \option{DEMC}, \option{twalk}, \option{IM}, \option{GG}, \option{AGG}), pass a \code{\link{Lowlander}} result to \code{\link{Specs_help}} so these are seeded from its kept in-bounds samples rather than left for \code{LaplacesDemon} to guess; see \code{\link{Specs_help}}. If the \option{Lowlander} argument is supplied and \option{Specs} is left at its default, this seeding happens automatically and \option{Specs} need not be built by hand at all. Supplying \option{Specs} explicitly always wins: \code{Highlander} then leaves it untouched, even when a \option{Lowlander} object is present. } \item{CMAargs}{ List; arguments to be passed into \code{\link[cmaeshpc]{cmaeshpc}}. } \item{LDargs}{ -List; arguments to be passed into \code{\link[LaplacesDemon]{LaplacesDemon}}. Note that \option{Algorithm} and \option{Specs} are passed into the \option{LDargs} by default. +List; arguments to be passed into \code{\link[LaplacesDemon]{LaplacesDemon}}. Note that \option{Algorithm} and \option{Specs} are passed into the \option{LDargs} by default. Elements the user omits fall back to \code{Niters[2]} for \option{Iterations} and 1 for \option{Thinning}. If a \option{Lowlander} object is supplied and \option{LDargs} itself is provided, any \option{Specs} already embedded in it counts as user-supplied and is likewise never overwritten. Because \option{Iterations} can differ between the CMA/MCMC passes (the final pass uses \option{NfinalMCMC}), auto-seeded \option{Specs} are rebuilt at the true length of each pass; a hand-built \option{DEMC} \code{Z} does not get that treatment and must be sized by the user for the pass it runs in. } \item{parm.names}{ Character vector; optional vector of parameter names. This will over-write the current names in \option{parm} if these also exist. If provided it ensures the \option{parm} and \option{CMA_last$par} outputs are named. @@ -105,6 +109,16 @@ Positive integer scalar; number of CPU cores to use (default 1). When \option{co } \details{ \code{Highlander} is designed to find good global solutions to highly multi-modal likelihood problems. It does this by oscillating between genetic algorithm stages, and classical MCMC stages. + +\subsection{Running with a Lowlander object}{ +\code{\link{Lowlander}} explores the parameter space with a Latin hypercube design and keeps the samples with the most plausible likelihood, returning shrunk bounds, the best point, and the kept samples themselves. Handing that result to \code{Highlander} via \option{Lowlander} does two things: it seeds \option{Specs} (when you have not supplied your own) so algorithms that need several starting points begin inside the plausible region, and it fills \option{parm}, \option{lower} and \option{upper} where they were left \option{NULL}. + +The reason this matters is that when \code{Specs$Z} (or \code{SIV}, \code{mu}, \code{Grid}) is left \code{NULL}, \code{LaplacesDemon} falls back to \code{GIV()}, which searches for usable initial values by evaluating the model at random Normal draws with means in \code{[-100, 100]} and standard deviations in \code{[0.1, 1000]}. On likelihoods that are undefined outside a much narrower region -- the usual case for physical models -- \code{GIV} gives up and returns \code{NA}, leaving walkers stacked at \code{NA} so the chain aborts or never approaches the mode. + +The interaction with \option{parm} is worth noting explicitly. \code{LaplacesDemon} evaluates the model at \option{parm} before it looks at \option{Specs$Z}, and stops outright with \code{"The posterior is infinite!"} if that is not finite. Worse, for \option{AIES} the reported thinned chain is walker 1's trajectory and \code{Z[1,]} is ignored, so a good \option{Z} cannot rescue a \option{parm} that sits outside the finite region. Defaulting \option{parm} to \code{Lowlander$best} avoids that trap; passing a \option{parm} from outside the region warns rather than overriding. + +All other behaviour is unchanged: with \option{Lowlander} = \code{NULL} (the default) nothing here is different from previous versions. +} } \value{ A list output containing: @@ -127,13 +141,13 @@ A list output containing: \item{best_job}{Integer scalar; present when \option{cores} > 1 and \option{keepall} = TRUE. Index (1-based) of the parallel chain that achieved the highest LP.} \item{High_jobs}{List; present when \option{cores} > 1 and \option{keepall} = TRUE. Contains the full \code{Highlander} output list for every parallel chain, in chain order. Use \option{best_job} to select the winning chain: \code{High_out$High_jobs[[High_out$best_job]]}.} -When run in multicore mode \option{parm} through to \option{LD_all} inclusive are the results of the 'best' run. \option{LD_last_comb}, \option{best_job and \option{High_jobs} are only present when \option{cores} > 1 and \option{keepall} = TRUE.} +When run in multicore mode \option{parm} through to \option{LD_all} inclusive are the results of the 'best' run. \option{LD_last_comb}, \option{best_job} and \option{High_jobs} are only present when \option{cores} > 1 and \option{keepall} = TRUE. } \author{ Aaron Robotham } \seealso{ -\code{\link[cmaeshpc]{cmaeshpc}}, \code{\link[LaplacesDemon]{LaplacesDemon}} +\code{\link[cmaeshpc]{cmaeshpc}}, \code{\link[LaplacesDemon]{LaplacesDemon}}, \code{\link{Lowlander}}, \code{\link{Specs_help}}, \code{\link{Lowlander_helpers}} } \examples{ # Create some fake data (note it needs to be a list, but data need not be called 'x'): @@ -176,6 +190,32 @@ highfit_prior = Highlander(parm=c(-1,2), Data=mock, likefunc=likefunc, liketype= prior=priorfunc, lower=c(-1e3,1), upper=c(1e3,1e3), parm.names=c('Mean', 'SD')) highfit_prior$parm + +\dontrun{ +# A Lowlander pre-pass finds in-bounds, plausible samples, and handing it to +# Highlander seeds Specs from those samples and fills parm/lower/upper where +# left NULL. This matters for algorithms that need several plausible starting +# points (AIES, DEMC, twalk, IM, GG, AGG), where LaplacesDemon otherwise guesses +# via GIV() and usually fails on likelihoods that are undefined out of a narrow +# region: + +intervals = Lowlander(c(-1e3, 1), c(1e3, 1e3), Data = mock, likefunc = likefunc, + liketype = 'max', parm.names = c('Mean', 'SD')) + +highfit_ll = Highlander(Data = mock, likefunc = likefunc, liketype = 'max', + Lowlander = intervals, Algorithm = 'AIES', + parm.names = c('Mean', 'SD')) + +highfit_ll$parm + +# Passing parm, lower or upper explicitly is never overridden by the Lowlander +# object, and neither is a Specs list you build yourself: + +highfit_ll2 = Highlander(parm = c(20, 100), Data = mock, likefunc = likefunc, + liketype = 'max', Lowlander = intervals, + Algorithm = 'AIES', parm.names = c('Mean', 'SD'), + lower = c(-1e3, 1), upper = c(1e3, 1e3)) +} } % Add one or more standard keywords, see file 'KEYWORDS' in the % R documentation directory. diff --git a/man/Lowlander.Rd b/man/Lowlander.Rd index 9fc4bc1..fd038b7 100644 --- a/man/Lowlander.Rd +++ b/man/Lowlander.Rd @@ -29,7 +29,15 @@ List; passed as the second argument to \option{likefunc} at every evaluation. } \item{likefunc}{ Function; takes arguments \option{parm} (numeric vector) and \option{Data} -(list) and returns a single numeric scalar. +(list) and returns a single numeric scalar. This is stricter than +\code{\link{Highlander}}, which accepts the \code{\link[LaplacesDemon]{LaplacesDemon}} +list with \code{LP}, \code{Dev} and \code{Monitor} elements: \option{likefunc} +is evaluated once per sample and its result must be one numeric value per row, +because \code{keep} indexes rows of \option{latin_mod}. A function that returns +a list or a longer vector is rejected up front with a message naming the +offending samples, so an \code{LD}-style function needs wrapping, e.g. +\code{likefunc = function(parm, Data) mylikefunc(parm, Data)$LP}. Note that +\option{liketype} below refers to that extracted scalar. } \item{Nsamp}{ Positive integer scalar; number of Latin-hypercube samples to draw. Default 'auto' will @@ -120,12 +128,13 @@ An invisible list containing: \item{cor_data}{Numeric matric; correlation matrix between all parameters and \option{output} (latter is the last row and column).} \item{latin}{Numeric matrix; the \option{Nsamp} x \option{Npar} unit-cube LHS design used.} \item{latin_mod}{Numeric matrix; the \option{Nsamp} x \option{Npar} parameter matrix mapped to physical bounds.} +\item{liketype}{Character scalar; the \option{liketype} actually used, recorded so that \code{\link{Lowlander_kept}} can rank the kept samples in the correct direction without being told again.} } \author{ Aaron Robotham } \seealso{ -\code{\link{Highlander}}, \code{\link[lhs]{randomLHS}} +\code{\link{Highlander}}, \code{\link{Lowlander_kept}}, \code{\link{Lowlander_Z}}, \code{\link{Lowlander_best}}, \code{\link[lhs]{randomLHS}} } \examples{ \dontrun{ diff --git a/man/Lowlander_helpers.Rd b/man/Lowlander_helpers.Rd new file mode 100644 index 0000000..fe05e9f --- /dev/null +++ b/man/Lowlander_helpers.Rd @@ -0,0 +1,102 @@ +\name{Lowlander_helpers} +\alias{Lowlander_helpers} +\alias{Lowlander_kept} +\alias{Lowlander_Z} +\alias{Lowlander_best} +\title{ +Extract usable starting values from a Lowlander result +} +\description{ +Helpers that pull in-bounds, plausible parameter values out of the return value +of \code{\link{Lowlander}}, for seeding algorithms that need them. These are +used by \code{\link{Specs_help}} when a \option{Lowlander} object is supplied, +and are exported so they can be used directly. + +\code{Lowlander_kept} returns the kept samples ranked best-first, \code{Lowlander_Z} +returns a design matrix of \option{N} such samples spread across the plausible +region, and \code{Lowlander_best} returns the single most plausible point. +} +\usage{ +Lowlander_kept(Lowlander = NULL, liketype = NULL, parm.names = NULL) + +Lowlander_Z(Lowlander = NULL, N = 1, liketype = NULL, parm.names = NULL) + +Lowlander_best(Lowlander = NULL, parm.names = NULL) +} +\arguments{ + \item{Lowlander}{ +List; the (invisible) return value of a \code{\link{Lowlander}} call. +\code{Lowlander_kept} and \code{Lowlander_Z} require the \code{latin_mod}, +\code{keep} and \code{output} elements; \code{Lowlander_best} requires +\code{best}. If \code{NULL}, \code{NULL} is returned so callers can fall back to +their own defaults. +} + \item{N}{ +Integer scalar; number of rows (walker positions) requested from +\code{Lowlander_Z}. If \option{N} exceeds the number of kept samples the samples +are recycled and the repeats jittered, since exactly-duplicated walker positions +degenerate the affine-invariant stretch move. +} + \item{liketype}{ +Character scalar; \option{'min'} or \option{'max'}, the direction in which +\code{output} should be read. Only used as a fallback: when the +\option{Lowlander} object records the \code{liketype} it was run with (the case +for \code{Lowlander} versions that return it), that value takes precedence so the +ranking cannot disagree with how the samples were selected. +} + \item{parm.names}{ +Character vector; optional names for the columns of the result (and for the +returned vector from \code{Lowlander_best}). Only applied when its length matches +the number of parameters. +} +} +\details{ +Lowlander restricts \code{keep} to samples with a finite likelihood that also +beat its quantile cut, so every row returned here lies inside the region +Lowlander judged plausible. This is what makes them safe as sampler seeds, in +contrast to \code{LaplacesDemon}'s own \code{GIV()} initial-value generator, +which proposes from broad Normal distributions and frequently returns +\code{NA} for likelihoods that are undefined outside a narrow region. + +The rows of \code{Lowlander_Z} are taken at evenly-spaced quantile ranks of the +kept samples rather than as the best \option{N} points, so the ensemble covers +the width of the plausible region instead of collapsing onto the mode. In +testing the choice of rank spacing made little difference; what matters is +supplying in-bounds seeds at all rather than leaving \code{Z = NULL}. + +\code{Lowlander_Z} requires at least two usable kept samples and warns (returning +\code{NULL}) when fewer are available. +} +\value{ +\code{Lowlander_kept} and \code{Lowlander_Z} return a numeric matrix with one +column per parameter, or \code{NULL} if nothing usable is available. +\code{Lowlander_best} returns a numeric vector, or \code{NULL}. +} +\author{ +Aaron Robotham +} +\seealso{ +\code{\link{Lowlander}}, \code{\link{Specs_help}}, \code{\link{Highlander}} +} +\examples{ +\dontrun{ +likefunc = function(parm, Data) { + sum((parm - c(2, 5))^2) +} +Data = list(parm.names = c('a', 'b'), N = 2) + +LL = Lowlander(c(0, 0), c(10, 10), Data = Data, likefunc = likefunc, + liketype = 'min', Nsamp = 500, pcut = 0.1) + +# The kept samples, ranked best-first: +Lowlander_kept(LL, parm.names = Data$parm.names) + +# 8 spread walker positions, suitable as Specs$Z for AIES: +Lowlander_Z(LL, N = 8, parm.names = Data$parm.names) + +# The single best point, suitable as Initial.Values: +Lowlander_best(LL, parm.names = Data$parm.names) +} +} +\concept{ optimisation } +\concept{ MCMC } diff --git a/man/Specs_help.Rd b/man/Specs_help.Rd index ad7ec18..c43cb83 100644 --- a/man/Specs_help.Rd +++ b/man/Specs_help.Rd @@ -12,7 +12,8 @@ is returned with placeholder values and an informative message is printed to guide the user. } \usage{ -Specs_help(Algorithm = 'CHARM', Data = NULL, ...) +Specs_help(Algorithm = 'CHARM', Data = NULL, Lowlander = NULL, + Iterations = NULL, Thinning = 1, ...) } \arguments{ \item{Algorithm}{ @@ -25,11 +26,37 @@ OHSS, pCN, RAM, RDMH, Refractive, RJ, RSS, RWM, SAMWG, SGLD, Slice, SMWG, THMC, twalk, UESS, USAMWG, USMWG. } \item{Data}{ -List; AIES can generate sensible default \option{Nc} if provided with -LD style \option{Data}. This defaults to 2 x length(\option{parm.names}). +List; used to determine the number of parameters from \code{Data$parm.names}, +which drives \option{Nc} for the ensemble samplers (\option{AIES}, \option{DEMC}) +and the length of \option{epsilon}/\option{m} for the Hamiltonian family. If +\code{Data$parm.names} is absent, the parameter count is taken from the +\option{Lowlander} object instead. +} + \item{Lowlander}{ +List; the return value of a \code{\link{Lowlander}} call. When supplied, the +parts of \code{Specs} that \code{LaplacesDemon} would otherwise guess blind +(via its \code{GIV} initial-value generator) are filled from Lowlander's kept, +in-bounds samples. This is the recommended way to run \option{AIES}, +\option{DEMC}, \option{twalk}, \option{IM}, \option{GG} and \option{AGG} on +likelihoods that are only well-behaved inside a restricted region, since +\code{GIV} routinely proposes wildly out-of-bounds points and the samplers then +either abort or stall. See Details. +} + \item{Iterations}{ +Integer scalar; the number of MCMC iterations that will be passed to +\code{\link[LaplacesDemon]{LaplacesDemon}}. Required to size the \option{Z} +array for \option{DEMC}, which must have +\code{floor(Iterations/Thinning) + 1} rows. If \code{NULL}, DEMC's \option{Z} is +left as \code{NULL} (and \code{LaplacesDemon} generates it) with a message. +} + \item{Thinning}{ +Integer scalar; the thinning interval corresponding to \option{Iterations}, used +only for sizing DEMC's \option{Z} array. Defaults to 1. } \item{\dots}{ -List; Named list where we wish to override the defaults for a particular algorithm. See Examples. +Named values used to override the defaults for a particular algorithm. See +Examples. Names that do not correspond to an element of that algorithm's +\code{Specs} produce a warning rather than being silently ignored. } } \details{ @@ -40,6 +67,52 @@ the user to supply all \code{Specs} values (e.g. \option{AIES}, values is returned and a message is printed advising the user of what must be filled in. The algorithms \option{DRM} and \option{RDMH} do not use a \code{Specs} argument, so \code{NULL} is returned. + +\subsection{Seeding from Lowlander}{ +When \code{Specs$Z} (or \code{SIV}, \code{mu}, \code{Grid}) is left \code{NULL}, +\code{LaplacesDemon} falls back to \code{GIV()}, which searches for usable +initial values by evaluating the model at random Normal draws with means in +\code{[-100, 100]} and standard deviations in \code{[0.1, 1000]}. On +astrophysical likelihoods that are undefined (or return non-finite values) +outside a much narrower region, \code{GIV} frequently gives up and returns +\code{NA}, which leaves ensemble walkers stacked at \code{NA} and the chain +either aborts or never approaches the mode. + +Passing the output of a \code{\link{Lowlander}} call avoids this. Lowlander has +already identified a set of in-bounds, plausible parameter samples, so: +\itemize{ +\item \option{AIES} \code{Z} becomes an \code{Nc} x \code{Npar} matrix of kept +samples taken at evenly-spaced quantile ranks (\code{\link{Lowlander_Z}}). +\item \option{DEMC} \code{Z} becomes the required 3-D array, with each chain +filled from an independent draw of the kept samples. +\item \option{twalk} \code{SIV} becomes a plausible point that is deliberately +not the best one, since \code{LaplacesDemon} stops if \code{SIV} coincides with +\code{Initial.Values}. +\item \option{IM} \code{mu} becomes the best kept point. +\item \option{GG}/\option{AGG} \code{Grid} becomes a per-parameter list of +\emph{offset} grids scaled to each parameter's own plausible width, rather than +\code{LaplacesDemon}'s fixed \code{seq(-0.1, 0.1)}. +} +} + +\subsection{A good \code{Z} cannot rescue a bad \code{Initial.Values}}{ +\code{LaplacesDemon} evaluates the model at \code{Initial.Values} before it ever +looks at \code{Z}, and stops outright with \code{"The posterior is infinite!"} if +that value is not finite. For \option{AIES} an explicitly supplied \code{Z} is +taken \emph{whole} -- rows 2 to \code{Nc} become walkers 2 to \code{Nc}, but +walker 1 starts from \code{Initial.Values} and \code{Z[1,]} is ignored. The +reported thinned chain is walker 1's trajectory, so it is \code{Initial.Values}, +not \code{Z[1,]}, that anchors what you get back. (Walker 1 is not frozen at +\code{Initial.Values}: it proposes and moves like any other walker. Only the +\emph{starting} value is taken from there.) + +\code{Highlander} passes a \code{Lowlander}-seeded \code{Z} through as-is, so if +\code{parm} is outside the region where the likelihood is finite the run can +abort regardless of how good \code{Z} is. When using Lowlander to bound the +problem, take \code{Initial.Values} from \code{Lowlander(...)$best} as well -- +which is what \code{Highlander} does automatically if \code{parm} is left +\code{NULL}. +} } \value{ A named list suitable for passing as the \code{Specs} argument of @@ -67,10 +140,40 @@ Specs_help('AM') # Metropolis-within-Gibbs: Specs_help('MWG') -# Affine-Invariant Ensemble Sampler (requires user values; prints a message): -Specs_help('AIES') - # Algorithms with no Specs return NULL: Specs_help('DRM') + +\dontrun{ +# The interesting case: seed the ensemble samplers from a Lowlander run so that +# LaplacesDemon never has to guess initial values blind. +lower = c(-4, -2) +upper = c(3, 1) +Data = list(parm.names = c('mu', 'sigma'), mon.names = 'LP', N = 100) + +# Note the two likefunc contracts: Highlander wants the LaplacesDemon-style +# list (LP, Dev, Monitor), while Lowlander collapses its likefunc to one +# numeric scalar per sample. Wrap it and pull out the LP, as here (assuming +# mylikefunc is the LD-style list-returning function; 'max' matches how +# Highlander's own liketype defaults to 'max' for LD-style outputs: +intervals = Lowlander(lower, upper, Data = Data, liketype = 'max', + likefunc = function(parm, Data) mylikefunc(parm, Data)$LP) + +# Z now holds Nc in-bounds walker positions instead of NULL: +Specs_help('AIES', Data, Lowlander = intervals) + +# DEMC needs Iterations/Thinning to size its 3-D Z array: +Specs_help('DEMC', Data, Lowlander = intervals, Iterations = 1000, Thinning = 1) + +# And take Initial.Values from Lowlander's best point, because LaplacesDemon +# evaluates the model there before it ever looks at Z: +Highlander(intervals$best, Data, mylikefunc, lower = intervals$lower, + upper = intervals$upper, Algorithm = 'AIES', + Specs = Specs_help('AIES', Data, Lowlander = intervals)) + +# Simpler: Highlander does all of the above when given the Lowlander result, +# seeding Specs and filling parm/lower/upper wherever they are left NULL. +Highlander(Data = Data, likefunc = mylikefunc, Algorithm = 'AIES', + Lowlander = intervals) +} } \concept{ MCMC } diff --git a/vignettes/LD-Algorithms.Rmd b/vignettes/LD-Algorithms.Rmd index e078abd..b8e3868 100644 --- a/vignettes/LD-Algorithms.Rmd +++ b/vignettes/LD-Algorithms.Rmd @@ -120,6 +120,105 @@ Fit <- LaplacesDemon( Each algorithm may require a different `Specs` list. The required and optional elements of `Specs` for every algorithm are documented below. +## Auto-filling `Specs` with `Specs_help()` + +`Specs_help(Algorithm, Data)` returns the correct `Specs` list for any +algorithm, so it rarely needs to be written by hand: + +```{r specs-help, eval=FALSE} +Highlander(..., Algorithm = "AIES", + Specs = Specs_help("AIES", Data)) +``` + +Any named argument is applied over the default, and a name that does not +belong to that algorithm's `Specs` raises a warning rather than being +silently dropped: + +```{r specs-help-override, eval=FALSE} +Specs_help("AIES", Data, Nc = 16, beta = 3) +``` + +### Seeding from `Lowlander` + +Several algorithms need starting values that already sit in a region where +the likelihood is finite: `Z` for `AIES` and `DEMC`, `SIV` for `twalk`, `mu` +for `IM`, and the evaluation `Grid` for `GG`/`AGG`. If these are left `NULL`, +**LaplacesDemon** falls back to its own `GIV()` generator, which probes the +model at random Normal draws with means in $[-100, 100]$ and standard +deviations in $[0.1, 1000]$. On likelihoods that are undefined outside a much +narrower region (the usual case for physical models), `GIV` gives up and +returns `NA`, leaving walkers stacked at `NA` so the chain either aborts or +never approaches the mode. + +Since `Lowlander` has already found in-bounds, plausible parameter samples, the +usual thing is simply to hand its output to `Highlander`: + +```{r specs-help-lowlander, eval=FALSE} +# Lowlander wants a single number per sample, so pull the scalar out of the +# LaplacesDemon-style list that mylikefunc returns. +intervals = Lowlander(lower, upper, Data = Data, liketype = 'min', + likefunc = function(parm, Data) mylikefunc(parm, Data)$LP) + +Highlander(Data = Data, likefunc = mylikefunc, liketype = 'min', + Lowlander = intervals, Algorithm = "AIES") +``` + +Note the two `likefunc` contracts differ, and this is the one place where the +same function cannot be passed to both calls. `Highlander` wants what +**LaplacesDemon** wants: a list with `LP`, `Dev` and `Monitor` elements. +`Lowlander` evaluates its `likefunc` once per Latin-hypercube sample and needs a +numeric scalar back (see `?Lowlander`); `liketype = 'min'` / `'max'` then says +which end of those scalars to keep. Passing `Lowlander` an LD-style list is +caught up front with a message naming the samples that misbehaved and showing +the `$LP` wrapper, so wrap the function and select the element you want, as +above. +The two `liketype` arguments are separate decisions even when they happen to +agree: `Lowlander`'s refers to the scalar you extracted, while `Highlander`'s +refers to how it should interpret `mylikefunc`'s `LP` field. + +That single `Lowlander` argument seeds `Specs` and fills in `parm`, `lower` and +`upper` wherever they were left `NULL`, so `intervals$best` and `intervals$lower` +no longer have to be repeated by hand. Anything you *do* pass explicitly wins, so +a bespoke `parm`, set of bounds, or `Specs` list is never overwritten. One +caveat: `Lowlander` selects its samples from its own `likefunc` alone, so if you +also pass a `prior` to `Highlander`, the seeded region does not reflect it unless +the likelihood already folds it in. Treat the seeds as in-bounds starting points +rather than posterior draws. Calling `Specs_help()` yourself is still available +and equivalent for the seeding step: + +```{r specs-help-manual, eval=FALSE} +Specs = Specs_help("AIES", Data, Lowlander = intervals) +``` + +Note that starting `Initial.Values` from `intervals$best` matters independently +of `Z`. **LaplacesDemon** evaluates the model at `Initial.Values` *before* it +looks at `Z` and stops outright if that posterior is not finite; and when `Z` is +supplied, **AIES** takes rows 2..$N_c$ as walkers 2..$N_c$ while walker 1 starts +from `Initial.Values` (`Z[1,]` is ignored). The reported thinned chain is walker +1's trajectory, so a start outside the finite region can abort the fit no matter +how good `Z` is. + +`DEMC` needs one extra thing. Its `Z` is a 3-D array rather than a matrix, with +`floor(Iterations / Thinning) + 1` slices, so `Specs_help` can only build it when +told how long the run will be: + +```{r specs-help-demc, eval=FALSE} +Specs = Specs_help("DEMC", Data, Lowlander = intervals, + Iterations = 1000, Thinning = 1) +``` + +**LaplacesDemon** rewrites both numbers before it sizes `Z` — `Iterations` is +rounded and clamped up to a minimum of 11, and `Thinning` falls back to 1 unless +it lies in $[1, \text{Iterations}]$ — so `Specs_help` applies the same +adjustment rather than trusting the raw values. `Highlander` handles this for you +because it knows each MCMC pass's length; it also rebuilds `Z` on the final pass, +which runs for `NfinalMCMC` rather than `Niters[2]`. A hand-built array is left +alone, so it must be sized against the *adjusted* values for whichever pass uses +it. + +Without `Iterations`, `Z` is left `NULL` (with a message) and +**LaplacesDemon** generates it as before. + --- # Metropolis–Hastings Algorithms @@ -599,7 +698,7 @@ with $\beta = 2$. | Spec | Description | |------|-------------| | `Nc` | Number of walkers (≥ 3; must be even when `CPUs > 1`). | -| `Z` | Optional $N_c \times d$ matrix of initial walker positions (`NULL` to auto-generate). | +| `Z` | Optional $N_c \times d$ matrix of initial walker positions (`NULL` to auto-generate via `GIV()`, which often fails on restricted likelihoods — see *Seeding from `Lowlander`* above). | | `beta` | Scale parameter controlling the stretch move (default 2; must be > 1). | | `CPUs` | Number of CPUs for parallel walker updates (≥ 1). | | `Packages` | Packages to load on each CPU (`NULL` if none). | @@ -638,7 +737,7 @@ to a randomly-chosen history point. | `gamma` | Scale for difference vector (defaults to $2.38/\sqrt{2d}$). | | `Nc` | Number of parallel chains (≥ 3). | | `w` | Probability of snooker update (recommended 0.1). | -| `Z` | Optional matrix/array of previous thinned samples. | +| `Z` | Optional 3-D array of previous thinned samples, of dim `floor(Iterations/Thinning) + 1` x $d$ x `Nc`, with `Iterations` and `Thinning` taken *after* `LaplacesDemon` adjusts them. A plain matrix survives only if it has at least half that many rows (it gets row-recycled, which duplicates history and weakens the proposals); below that it dies with `subscript out of bounds`, and an array whose first dimension is off is a hard stop. Use `Specs_help("DEMC", Data, Lowlander = ..., Iterations = ..., Thinning = ...)` (or pass `Lowlander` to `Highlander`) to build it. | ```{r demc} Fit <- LaplacesDemon(Model, MyData, Initial.Values, From fd8ab5bac28a0f1e7dfb6798e9f6c9a108b4d174 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 2 Sep 2026 10:39:08 +0000 Subject: [PATCH 2/2] Fix review thread logic issues Co-authored-by: asgr <5617132+asgr@users.noreply.github.com> --- R/Highlander.R | 2 +- R/Specs_help.R | 12 +++++------- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/R/Highlander.R b/R/Highlander.R index ebde9d9..41df7c4 100644 --- a/R/Highlander.R +++ b/R/Highlander.R @@ -162,7 +162,7 @@ Highlander=function(parm=NULL, Data, likefunc, likefunctype=NULL, liketype=NULL, # is not finite. For AIES the reported thinned chain is walker 1's # trajectory and Z[1,] is ignored, so an out-of-region parm can abort the # run no matter how good the seeded Z is. - if(any(parm < Lowlander$lower | parm > Lowlander$upper)){ + if(any(parm < Lowlander$lower | parm > Lowlander$upper, na.rm = TRUE)){ message('Highlander: parm sits outside the Lowlander plausible region; LaplacesDemon evaluates the model there first and stops if the posterior is not finite, so this can abort the run regardless of the seeded Specs. Consider parm = Lowlander(...)$best.') } } diff --git a/R/Specs_help.R b/R/Specs_help.R index e711649..24d0699 100644 --- a/R/Specs_help.R +++ b/R/Specs_help.R @@ -135,7 +135,7 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, Lowlander = NULL, if (is.null(dots$Nc) && is.null(Npar)) { message('AIES requires user-supplied Nc or Data. Returning template default Nc = 3.') } else if (is.null(dots$Nc)) { - message('AIES: using Nc = 2 * length(parm.names) = ', Nc, '.') + message('AIES: using default Nc = ', Nc, '.') } Specs_list = list(Nc = Nc, Z = .Zmatrix(Nc), beta = 2, CPUs = 1, Packages = NULL, Dyn.libs = NULL) @@ -422,10 +422,7 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, Lowlander = NULL, } } else if (!is.array(Specs_list$Z) || length(dim(Specs_list$Z)) != 3 || dim(Specs_list$Z)[3] != Nc || - # identical(), not !=: Npar can be NULL here (no Data, no - # Lowlander), and `dim(Z)[2] != NULL` is logical(0), which turns - # the whole condition into NA and errors inside if(). - !identical(dim(Specs_list$Z)[2], Npar) || + (!is.null(Npar) && dim(Specs_list$Z)[2] != Npar) || (.sizeable && dim(Specs_list$Z)[1] != .z_rows(Iterations, Thinning))) { # DEMC wants a 3-D array: floor(Iterations/Thinning)+1 x Npar x Nc, with # Iterations/Thinning as LaplacesDemon *adjusts* them. A bare matrix @@ -433,7 +430,7 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, Lowlander = NULL, # of bounds", so rebuild from whatever we have. Z = Specs_list$Z Zn = NULL - if (.sizeable) { + if (.sizeable && !is.null(Npar)) { nrow_z = .z_rows(Iterations, Thinning) src = NULL if ('Z' %in% names(dots)) { @@ -458,9 +455,10 @@ Specs_help = function(Algorithm = 'CHARM', Data = NULL, Lowlander = NULL, } } if (is.null(Zn)) { + npar_msg = if (is.null(Npar)) 'Npar' else Npar warning(paste0('DEMC: Z must be a 3-D array of dim ', if (.sizeable) .z_rows(Iterations, Thinning) else 'floor(Iterations/Thinning)+1', - ' x ', Npar, ' x ', Nc, '. Set to NULL so LaplacesDemon ', + ' x ', npar_msg, ' x ', Nc, '. Set to NULL so LaplacesDemon ', 'regenerates it. Pass Iterations and Thinning to Specs_help ', 'to have the array built for you.')) Specs_list$Z = NULL