Better Z helper functions

.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:

DESCRIPTION

@@ -8,11 +8,13 @@ Maintainer: Aaron Robotham <aaron.robotham@uwa.edu.au>
 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,
   asgr/ProSpect
 Encoding: UTF-8
 LazyData: true
 VignetteBuilder: knitr
+Config/testthat/edition: 3

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")

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, 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.')
+      }
+    }
+    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))

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
   )))
 }