diff --git a/skills/acquisition-functions/SKILL.md b/skills/acquisition-functions/SKILL.md
index 9e0c917..c6fd654 100644
--- a/skills/acquisition-functions/SKILL.md
+++ b/skills/acquisition-functions/SKILL.md
@@ -1,334 +1,113 @@
 ---
 name: acquisition-functions
-description: Use when choosing or designing the acquisition function for a gpCAM experiment — deciding which built-in actually fits the scientist's goal (and confirming it with them rather than defaulting to expected improvement), or writing a custom one for exploration/exploitation balance, multi-objective targets, constrained search regions, cost-aware moves, change-mapping via gradient or radical gradient, knowledge gradient and noisy EI under observation noise, UCB/LCB, target-value, or threshold-finding criteria.
-gpcam_version: "8.4.x"
-fvgp_version: "4.8.x"
-last_verified: "2026-08-04 (gpCAM 23c886d + joint candidate-batch selection); multi-task pointer updated for gpCAM 8.7.0 on 2026-10-28"
+description: Use when choosing or writing an acquisition function for gpCAM — selecting a built-in criterion for exploration, optimization, batch, or multi-task runs, or implementing a custom callable for UCB/LCB, minimization EI, constraints, radical gradient, or threshold finding.
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: gpCAM Acquisition Functions
 
-*Verified against gpCAM 8.4.x / fvgp 4.8.x — last checked 2026-08-04 (gpCAM `23c886d` + joint candidate-batch selection).*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Choose — and where needed, write — the acquisition function that controls where gpCAM measures next.
 
 ## When to Use
 
-Any time an experiment script needs an acquisition function. That includes simply
-picking a built-in string, not only writing a custom callable:
-- Deciding what to pass to `ask(acquisition_function=...)`
-- Balancing exploration and exploitation
-- Multi-objective optimization
-- Constrained search (avoid forbidden regions)
-- Cost-aware acquisition (expensive moves)
-- Upper/lower confidence bounds
-- Probability of improvement
+- Deciding which built-in string or custom callable to pass to `ask(acquisition_function=...)`
+- Balancing exploration and exploitation, or switching strategies mid-run
+- Minimization, constrained search, multi-objective targets, change-mapping (`"gradient"` / `radical_gradient`), or level-set (`threshold_finder` / `"target probability"`) criteria
 
 ## Choose the Acquisition Function *With* the User
 
-**Never pick an acquisition function silently.** It is the single choice that most
-determines where the instrument actually goes, and the right answer depends on what
-the scientist wants to walk away with — which is rarely what their first sentence says.
-
-Before generating a script:
-
-1. **Name your recommendation and the one tradeoff it makes.** "I'd use `variance`,
-   which will give you an even-quality map but won't spend extra shots refining the
-   peak once it finds it."
-2. **Resolve the ambiguity in "find the best conditions."** That phrase covers two
-   different experiments. Ask which one:
-   - *the single best point* — they will set the instrument there and run production;
-     the rest of the space is disposable → exploitation-leaning (`ucb`, EI late)
-   - *a trustworthy map that contains the best point* — they need to see the landscape,
-     report it, or find secondary optima → exploration-leaning (`variance`, gradient
-     forms, `relative information entropy`)
-3. **Ask what they'd do with a wrong answer.** If over-exploiting a noise spike would
-   waste a shift of beam time, weight toward exploration.
-4. **Confirm before writing the loop**, and put the reasoning in a comment in the
-   generated script so the choice is revisitable.
-
-Do **not** default to `"expected improvement"`. It is the most famous acquisition
-function in the literature and the most over-applied in practice; it is also
-maximization-only and anchors to a single noisy observation — see
-[gpCAM-Specific Behavior](#gpcam-specific-behavior-you-must-know) below. And on gpCAM
-≤ 8.4.2 it degenerates to pure exploration (caveat a), so **ask which version they run
-before recommending it**.
-
-### Decision Table — start from what they want to *learn*
-
-| They say… | Start from | Tradeoff to state out loud |
-|---|---|---|
-| "map the whole space" / "replace my raster scan" | `"variance"` | even posterior confidence everywhere; won't refine the peak |
-| "find where the signal *changes*" / "find the edges, the transition" | `radical_gradient` callable (below), or `"gradient"` | targets steep regions; ignores flat regions even if unmeasured |
-| "find the single best point, I have limited time" | `"ucb"` | exploits toward the max; `beta` is fixed at 3.0 for the built-in string |
-| "find the best point, but **my measurements are noisy**" | `"noisy expected improvement"` | the correct EI under noise; costs more per evaluation than `"ucb"` |
-| "it keeps re-measuring the same spot" / "EI has gone flat" | `"knowledge gradient"` | lookahead — values a point by how much it would improve your best *decision*, so it keeps exploring where EI has stalled; the most expensive built-in |
-| "find where the signal crosses X" | `threshold_finder` callable (below) | concentrates on one level set; the rest of the space stays coarse |
-| "find conditions where output ≈ X" | `"target probability"` + `args={'a':lo,'b':hi}` | same — a band, not a global map |
-| "**minimize** this" | `"lcb"`, or a custom EI-for-minimization | **not** EI, PI, KG, NEI, or `"ucb"` — all are hardcoded for maximization |
-| "give me N points per round" (batch) | `"total correlation"` (or `"relative information entropy"`) | the only two that score a *set*, so the only two that can pick a batch jointly — see caveat (e) |
-| multi-task / spectra — *exploring* | `"relative information entropy set"` or `"variance"` | — |
-| multi-task / spectra — *optimizing* | `"knowledge gradient"` or `"noisy expected improvement"` | all three act on the task-summed objective, but KG/NEI use the exact cross-task covariance where plain EI assumes independent tasks — see caveat (d) |
-| "I already know roughly where the peak is, refine it" | `"expected improvement"` if noise is low, `"noisy expected improvement"` if not | greedy; plain EI assumes the incumbent is known exactly |
-| "I don't know yet" | `"variance"` for the first ~20 points, then revisit | you can change the acquisition between `ask()` calls — it is not fixed for the run |
-
-That last row is worth saying to the scientist explicitly: **the acquisition function
-is an argument to `ask()`, not a property of the GP.** A common good design is to
-explore first and switch to an exploitative acquisition once the posterior is
-informative.
-
-For worked examples of each archetype — how the scientist phrases it, what to ask back,
-and the resulting loop — see
-[references/choosing-an-acquisition-function.md](references/choosing-an-acquisition-function.md).
-
-## Acquisition Function Contract
-
-```python
-def my_acquisition(x, gp_optimizer):
-    """
-    Parameters
-    ----------
-    x : np.ndarray, shape (V, D)
-        Candidate points to evaluate.
-    gp_optimizer : GPOptimizer
-        The GP model. Use its posterior_mean() and posterior_covariance() methods.
-    
-    Returns
-    -------
-    scores : np.ndarray, shape (V,)
-        Score for each candidate. HIGHER = more desirable (function is MAXIMIZED).
-    """
-```
+**Never pick an acquisition function silently.** Before generating a script:
 
-**Key rule:** The acquisition function is **maximized**. Return higher values for points you want to measure.
-
-### Reading the GP's data inside an acquisition function
-
-`gpo.y_data` and `gpo.x_data` are public properties and are the right thing to use
-for an incumbent value (`y_best`). They return the data **in the GP's own modeling
-space**, which is the same space `posterior_mean()` returns — so the two are always
-directly comparable.
-
-`get_data()` returns the same arrays in a dict along with hyperparameters and
-metadata. Use `get_data()` for reporting and inspection; use the `x_data` / `y_data`
-properties inside acquisition functions, where you want the cheap attribute access
-and guaranteed space-consistency.
-
-> **With a transformed optimizer** (`LogGPOptimizer`, `LogitGPOptimizer`), `gpo.y_data`
-> holds the *transformed* observations — `log(y)` or `logit(y)` — matching
-> `posterior_mean()`. Every recipe below therefore stays internally consistent with no
-> changes. What you must **not** do is mix spaces: `get_data()["original y data"]`
-> is on the original scale and comparing it against `posterior_mean()` is a silent
-> units bug. Any **absolute constant** you hardcode — a threshold, a target value, a
-> forbidden output level — is also in latent space and must be transformed
-> (`np.log(threshold)`, `scipy.special.logit(threshold)`). See the
-> `transformed-optimizers-advanced` skill.
-
-## Built-in Options
-
-Pass these as strings to `gpo.ask(acquisition_function=...)`:
-
-| Name | String key | Best for | Avoid when |
-|------|-----------|----------|------------|
-| Variance | `"variance"` | Pure exploration / mapping | They need the peak refined, not the map |
-| Expected Improvement | `"expected improvement"` | Late refinement of a known maximum, low noise, sequential, single-task | Minimizing; noisy data; batch (`n>1`); correlated multi-task (caveat d); **any run on gpCAM ≤ 8.4.2** (caveat a) |
-| Probability of Improvement | `"probability of improvement"` | Risk-averse maximization | Minimizing; it is even greedier than EI and stalls easily |
-| Upper Confidence Bound | `"ucb"` | Maximization; the most reliable general-purpose optimizer here | Minimizing; you need `beta ≠ 3.0` (write the callable — caveat f) |
-| Lower Confidence Bound | `"lcb"` | Minimization | You need `beta ≠ 3.0` (write the callable) |
-| Predicted Maximum | `"maximum"` | Pure exploitation — mean only, no uncertainty | Almost always: with no uncertainty term it re-measures the same point forever |
-| Predicted Minimum | `"minimum"` | Pure exploitation for minimization | Same |
-| Gradient | `"gradient"` | Seek steepest regions of the posterior mean | Signal is flat but interesting; you want uniform coverage |
-| Target Probability | `"target probability"` | Find points with output near a target value | `args={'a':…, 'b':…}` not set — it raises |
-| Relative Information Entropy | `"relative information entropy"` | Information-theoretic exploration; scores a **set** | Speed matters — forces `vectorized=False` |
-| RIE Set | `"relative information entropy set"` | Batch acquisition; multi-task exploration | — |
-| Total Correlation | `"total correlation"` | Batch acquisition; scores a **set** | Speed matters — forces `vectorized=False` |
-| Knowledge Gradient | `"knowledge gradient"` | Maximization; lookahead, one-step-optimal, robust to noise; keeps exploring when EI stalls | Minimizing; large datasets or tight acquisition budgets — it loops per candidate over a joint posterior |
-| Noisy Expected Improvement | `"noisy expected improvement"` | Maximization on noisy measurements — the right EI when the incumbent isn't known exactly | Minimizing; noiseless data (plain EI is cheaper); tight compute budgets |
-
-### The two set-valued acquisitions, and how they differ
-
-`"total correlation"` and `"relative information entropy"` are the only built-ins that
-score a **whole set** with a single number rather than scoring each point on its own.
-That is what lets them judge a batch, and it is why they are the answer when a scientist
-asks "which `M` experiments should I run together". Both are KL divergences, but between
-different pairs of distributions, so they answer different questions.
-
-**`"total correlation"`** is the divergence between the joint distribution over the
-existing data together with the candidates, and that same distribution with the
-data-to-candidate cross-covariance removed and the candidates made mutually independent.
-So it measures how much statistical dependence is present: how redundant the candidates
-are with what has already been measured, and with each other. gpCAM negates it, so
-maximizing the acquisition *minimizes* dependence and prefers points that are novel
-relative to the existing data **and** mutually distinct.
-
-**`"relative information entropy"`** is the divergence between the prior and the
-posterior at the candidate points alone. It measures how much the data already collected
-has changed the belief there. Negated the same way, maximizing it prefers points where
-prior and posterior still agree — locations the data has not reached.
-
-Which to recommend: for "which measurements should I take *together*", `"total
-correlation"` is the more direct criterion, because it explicitly penalizes
-candidate-to-candidate correlation. `"relative information entropy"` also spreads points
-out, but only as a side effect of them all lying in unexplored regions, so it will
-tolerate neighbors that are individually uninformed yet redundant with one another.
-
-To sanity-check any built-in or custom acquisition on a grid of candidates without calling `ask()`:
-```python
-scores = gpo.evaluate_acquisition_function(x_grid, acquisition_function="ucb")
-```
+1. **Name your recommendation and its main tradeoff** (e.g., `"variance"` gives uniform coverage across the space, whereas `"ucb"` concentrates shots near the peak).
+2. **Disambiguate "find the best conditions":** ask whether they need the single best operating point (`"ucb"`, `"noisy expected improvement"`, `"knowledge gradient"`) or a trustworthy map of the landscape (`"variance"`, `"gradient"`, `radical_gradient`), and confirm whether they are **maximizing or minimizing**.
+3. **Remember the acquisition is an `ask()` argument, not a fixed GP property** — exploring with `"variance"` for the first ~20 points before switching to `"ucb"` or `"noisy expected improvement"` is often the best design.
 
-Plotting these scores over a 1-D or 2-D grid next to the posterior mean is the fastest
-way to show a scientist *why* one acquisition sends the instrument somewhere different
-from another. Do this when they're unsure which to pick.
+Do **not** default to `"expected improvement"`: it is maximization-only, anchors to a single noisy observation (`np.max(y_data)`), and in legacy gpCAM ≤ 8.4.2 degenerated to pure exploration below the incumbent.
 
-## gpCAM-Specific Behavior You Must Know
+For per-archetype dialogue prompts ("Ask back...") and loop templates across all 9 experimental archetypes, see [references/choosing-an-acquisition-function.md](references/choosing-an-acquisition-function.md).
 
-These are properties of gpCAM's implementations, not of the textbook definitions. Tell
-the scientist about any that apply to their choice — several will otherwise silently do
-something other than what they were told.
+## Built-in Options & Decision Guide
 
-**(a) EI degenerated to pure exploration in gpCAM ≤ 8.4.2 — check which version they run.**
-Fixed by [#55](https://github.com/lbl-camera/gpCAM/issues/55); `surrogate_model.py` now
-standardizes first and needs no clip at all (the `np.maximum` is only a guard against
-floating-point cancellation when `σ` is at or below the epsilon), in the shared
-`_expected_improvement` helper:
+Pass these strings to `gpo.ask(acquisition_function=...)`:
 
+| Goal / Scenario | String key (or callable) | Key behavior & tradeoffs |
+|---|---|---|
+| Pure exploration / replace raster scan | `"variance"` | Uniform posterior confidence everywhere; won't refine peaks |
+| Maximize (low noise, fast) | `"ucb"` | `mean + 3.0 * std`; reliable general-purpose maximizer (`beta=3.0` fixed — use a callable to tune `beta`) |
+| **Minimize** | `"lcb"` (or `expected_improvement_minimize`) | `-(mean - 3.0 * std)`; **never** use EI, PI, UCB, KG, or NEI for minimization |
+| Maximize under **observation noise** | `"noisy expected improvement"` | Integrates over the incumbent posterior (`nei_samples`, `nei_reference_set_size`, `nei_seed` via `args`); slower than `"ucb"` |
+| Maximize when greedy EI/UCB stalls | `"knowledge gradient"` | One-step lookahead expected gain in `max posterior_mean` (`kg_reference_set_size`, `kg_seed` via `args`); heaviest built-in |
+| Late local refinement (low noise) | `"expected improvement"` | Maximization-only; anchors to `np.max(y_data)` (single-task) or `_incumbent_task_sum` (multi-task) |
+| Risk-averse maximization | `"probability of improvement"` | Maximization-only; greedier than EI and stalls easily |
+| Map steep transitions / edges | `"gradient"` or `radical_gradient` callable | `"gradient"` scores `‖∇m(x)‖ · σ(x)` (sharpest ridge); `radical_gradient` uses `√‖∇m(x)‖ · σ(x)` to spread across multiple transitions |
+| Find where output ∈ `[a, b]` | `"target probability"` | Requires `args={"a": lo, "b": hi}` on `GPOptimizer(...)` or `gpo.ask(..., args=...)` |
+| Find level set `f(x) = threshold` | `threshold_finder` callable | Concentrates on one contour; alternate with `"variance"` if a full map is also needed |
+| Batch (`n > 1` points jointly) | `"total correlation"` or `"relative information entropy"` | The only two built-ins that score a **set** jointly (forces `vectorized=False`; `"total correlation"` unavailable in `"kronecker"` mode) |
+| Multi-task exploration | `"relative information entropy set"` or `"variance"` | Vectorized exploration across `x_out` tasks |
+| Multi-task maximization (`Σ_t f(x, t)`) | `"knowledge gradient"` or `"noisy expected improvement"` | Uses exact cross-task covariance (`_scalarized_blocks`), whereas multi-task `"expected improvement"` assumes independent tasks |
+| Pure exploitation (no uncertainty) | `"maximum"` / `"minimum"` | Posterior mean only (`σ` ignored) — almost always stalls at one point |
+
+To compare any built-in or custom acquisition on a grid before running the loop:
 ```python
-z = imp / (std + 1e-9)
-return np.maximum(imp * norm.cdf(z) + std * norm.pdf(z), 0.)
+scores = gpo.evaluate_acquisition_function(x_grid, acquisition_function="ucb")
 ```
+(`evaluate_acquisition_function` flips the internal minimization sign back at `gp_optimizer_base.py:301` so higher scores mean higher priority.)
 
-The old code clipped `imp` to 0 *first* and only then divided by `σ`, so at every
-candidate whose posterior mean sat **below** the incumbent the score collapsed to the
-constant `σ · φ(0) = 0.3989 · σ` — the `"variance"` acquisition up to a positive factor.
-Until the posterior mean somewhere exceeded the best observed value, EI was pure
-exploration. This matters because **8.3.9 is still the stable line and 8.4.2 the latest
-release**, so most scientists you talk to are running the broken version. If they are,
-tell them plainly that their EI run was exploration, and recommend upgrading or using
-`"knowledge gradient"` / `"noisy expected improvement"` (never affected).
+### Set-valued batch acquisitions (`"total correlation"` vs `"relative information entropy"`)
 
-Check which behavior their install has, in four lines:
+`"total correlation"` and `"relative information entropy"` are the only built-ins that score a **whole candidate set** with a single number rather than scoring each point independently:
+- **`"total correlation"`** measures the KL divergence between the joint distribution over existing data + candidates and the factored distribution with data-to-candidate and candidate-to-candidate correlations removed. Negated for maximization, it penalizes both redundancy with existing data **and** mutual redundancy among candidates in the batch. (Not available when `multi_task_mode="kronecker"`, `gp_optimizer_base.py:443-445`.)
+- **`"relative information entropy"`** measures the KL divergence between the prior and posterior at the candidates alone, preferring unexplored regions without explicitly penalizing candidate-to-candidate proximity.
+
+## gpCAM-Specific Behavior You Must Know
 
+**(a) Legacy gpCAM ≤ 8.4.2 EI bug (fixed in #55 / gpCAM ≥ 8.5).**
+In `surrogate_model.py:409-412`, `_expected_improvement` standardizes `z = imp / (std + 1e-9)` before evaluating `imp * norm.cdf(z) + std * norm.pdf(z)`. In gpCAM ≤ 8.4.2, `imp` was clipped to 0 *before* dividing by `std`, collapsing EI to `0.3989 · σ(x)` wherever `mean < y_best`. If a user is on an older install, check with:
 ```python
 ei  = gpo.evaluate_acquisition_function(x_grid, acquisition_function="expected improvement")
 var = gpo.evaluate_acquisition_function(x_grid, acquisition_function="variance")
 below = gpo.posterior_mean(x_grid)["m(x)"] < np.max(gpo.y_data)
-print(np.unique(np.round(ei[below] / var[below], 6)))
-# one value == [0.398942] -> affected build; a spread of values below 0.398942 -> fixed
+print(np.unique(np.round(ei[below] / var[below], 6)))  # [0.398942] on <=8.4.2; spread below 0.398942 on >=8.5
 ```
 
-(The public `gpo.evaluate_acquisition_function` returns the higher-is-better score —
-it flips the internal sign back at `gp_optimizer_base.py:297`. Don't negate it again.)
-
-Note this fixes only the *degeneration*. Caveats (b) and (c) below — maximization-only,
-and anchoring to a noisy observation — are properties of plain EI and still apply.
-
-**(b) EI, PI, KG, and NEI are all maximization-only.** `last_best = np.max(gpo.y_data)`
-is hardcoded for single-task EI/PI (`surrogate_model.py:251, 258`), multi-task EI takes
-the max over the task-summed observations, and KG/NEI take `np.max(mu_ref)` over the
-reference set. If the scientist is minimizing, every one of
-them produces meaningless behavior with no error. Use `"lcb"` or the
-`expected_improvement_minimize` callable below.
-
-**(c) Plain EI/PI anchor to a noisy observation.** `np.max(y_data)` is the largest
-*measured* value, not the posterior mean at that location, so on noisy data EI
-over-exploits around a noise spike. **This is what `"noisy expected improvement"` is
-for** — it integrates over the posterior of the incumbent instead of treating it as
-known. Prefer NEI (or `"knowledge gradient"`) on noisy measurements; reach for `"ucb"`
-only if you also need the cheaper evaluation.
-
-**(d) Multi-task plain EI assumes the tasks are independent.** The `x_out` branch
-(`surrogate_model.py:312-321`) scalarizes to the task-summed objective
-`g(x) = Σ_t f(x, t)`, and since #55 it does so coherently: the incumbent is the best
-task-*summed* observation (`_observed_task_sums`) and the spread is
-`sqrt(Σ_t Var[f(x,t)])`. That last step drops the cross-task covariance, so the spread
-is right only if the tasks are uncorrelated — usually they are not, which is the whole
-reason to fit a multi-task GP. `"knowledge gradient"` and `"noisy expected improvement"`
-scalarize the same objective but take the exact cross-task covariance from
-`_scalarized_blocks`, at O((N·T)²) cost per evaluation. So for correlated tasks prefer
-KG/NEI; plain multi-task EI is the cheap approximation, not the accurate one.
-
-In gpCAM ≤ 8.4.2 this branch was worse than approximate — it *summed the standard
-deviations* (not the variances) and compared the task-sum against `np.max(gpo.y_data)`,
-a single `(point, task)` entry of the product-space vector. On that build, do not
-recommend plain multi-task EI at all.
-
-**(e) Batch behavior depends on whether you passed bounds or candidates.**
-
-With an `np.ndarray` `input_set`, `n > 1`, `method != "hgdl"`, and a string acquisition,
-`gpcam/gp_optimizer_base.py:528-536` replaces whatever was requested with
-`"total correlation"` (a warning is emitted). Batch EI is not EI — and neither is batch
-KG or batch NEI. Either use a genuinely batch-aware acquisition, or pass a callable and
-`method="hgdl"` with a Dask client.
-
-With a **candidate list** and `n > 1` the split is different, and it is the one that
-decides whether a batch is any good. `"total correlation"` and `"relative information
-entropy"` score a whole set with one number, so gpCAM selects the batch jointly by
-greedy forward selection and `'f_a(x)'` comes back as a **single value for the set**.
-Every other acquisition — including any callable, which gpCAM cannot introspect — is
-scored one point at a time, `n` truncates the sorted list, and `ask()` warns that the
-points are not mutually optimal. Take that warning seriously: nothing stops the highest
-individual scorers from sitting on top of each other, so a point-wise top-`n` from a
-dense candidate pool is very often several copies of the same experiment. If the
-scientist wants "the `M` most valuable experiments to run together", that is the
-set-valued path, not the point-wise one.
-
-**(f) `"ucb"` and `"lcb"` hardcode `beta = 3.0`** (`surrogate_model.py:232, 236`). The
-strings expose no tuning. To change the exploration/exploitation balance you must pass
-the callable form below.
-
-**(g) KG and NEI are the expensive ones.** Both build a joint posterior over the
-reference set plus candidates on every evaluation, and KG loops in Python over
-candidates. On a large dataset or a tight acquisition budget this dominates the loop.
-Keep `kg_reference_set_size` / `nei_reference_set_size` modest (default 100, which
-subsamples the data), and use `method="global"` — the `"local"` finite-difference
-optimizer is weak on their flat score surfaces, exactly as with plain EI.
-
-### Knowledge gradient and noisy expected improvement
-
-Both are optimization acquisitions (find the maximum) that reason about the posterior
-of the underlying **function** rather than the raw noisy observations. They are the
-right choice over plain `"expected improvement"` when measurements are noisy, and
-they work for both single-task and multi-task GPs.
+**(b) EI, PI, UCB, KG, and NEI are all maximization-only.**
+Single-task EI/PI hardcode `last_best = np.max(gp.y_data)` (`surrogate_model.py:306, 313`), multi-task EI maximizes the task sum (`_incumbent_task_sum`, `surrogate_model.py:415-427`), and KG/NEI take `np.max(mu_ref)` (`surrogate_model.py:462, 513`). For minimization, use `"lcb"` or `expected_improvement_minimize` below.
+
+**(c) Single-task plain EI/PI anchor to a noisy observation.**
+`np.max(gpo.y_data)` is the largest raw reading, not the posterior mean, so on noisy data plain EI over-exploits noise spikes. Prefer `"noisy expected improvement"` (or `"knowledge gradient"`), which integrates over the posterior at the observed reference points (`surrogate_model.py:495-520`).
+
+**(d) Multi-task plain EI assumes uncorrelated tasks.**
+When `x_out` is set (`surrogate_model.py:375-389, 415-427`), multi-task `"expected improvement"` scalarizes to the task-summed objective `g(x) = Σ_t f(x, t)` using `_incumbent_task_sum` (the maximum task-summed **posterior mean** at the observed inputs, so missing `NaN` task entries in `"cartesian"` mode do not poison EI) and `std = sqrt(Σ_t Var[f(x, t)])`. That spread drops cross-task covariance, so it is exact only when tasks are independent. `"knowledge gradient"` and `"noisy expected improvement"` scalarize the same `Σ_t f(x, t)` objective using the full cross-task covariance via `_scalarized_blocks` (`surrogate_model.py:430-452`).
+
+**(e) Batch (`n > 1`) behavior depends on `input_set` type.**
+- With a continuous `np.ndarray` bounds `input_set`, `n > 1`, `method != "hgdl"`, and a string acquisition other than `"total correlation"` or `"relative information entropy"`, `gp_optimizer_base.py:583-594` warns and overrides the string to `"total correlation"` (or `"relative information entropy"` when `multi_task_mode == "kronecker"`). To get `n > 1` local optima of a point-wise acquisition over bounds, pass a callable with `method="hgdl"` and a `dask_client`.
+- With a discrete **candidate list** and `n > 1`, `"total correlation"` and `"relative information entropy"` select the batch jointly via greedy forward selection (`gp_optimizer_base.py:618-625`), returning a single scalar score in `"f_a(x)"`. All other strings and custom callables are scored point-by-point and truncated to the top `n` (`gp_optimizer_base.py:626-636`), which warns and can pick clustered duplicates.
 
+**(f) `"ucb"` and `"lcb"` fix `beta = 3.0`** (`surrogate_model.py:291, 295, 366, 372`). Write the 5-line callable below to change `beta`.
+
+**(g) KG and NEI cost O((N_ref + V)²) per call.**
+Keep `kg_reference_set_size` / `nei_reference_set_size` modest (default `100`) and pass via `GPOptimizer(..., args=...)` or `gpo.ask(..., args=...)`:
 ```python
-# Single-task
-gpo.ask(bounds, acquisition_function="knowledge gradient")
-gpo.ask(bounds, acquisition_function="noisy expected improvement")
-
-# Multi-task — both act on the task-summed objective sum_t f(x, t),
-# matching the built-in multi-task "expected improvement" / "ucb"
-gpo.ask(bounds, x_out=np.array([0, 1, 2]),
-        acquisition_function="knowledge gradient")
+gpo.ask(bounds, acquisition_function="noisy expected improvement",
+        args={"nei_samples": 128, "nei_reference_set_size": 100, "nei_seed": 0})
 ```
 
-- **Knowledge gradient (KG)** — the expected increase in the maximum of the posterior
-  *mean* after a hypothetical measurement at the candidate (KGCP scheme: the reference
-  set is the data points plus the candidate). It is a one-step-optimal, lookahead
-  criterion: it values a measurement by how much it is expected to improve your best
-  decision, so it keeps exploring even when simple EI has gone to zero. KG is ~0 on
-  top of existing data and largest where a measurement would most raise the achievable
-  maximum.
-- **Noisy expected improvement (NEI)** — expected improvement averaged over the
-  posterior of the incumbent (the best value so far). Plain EI treats the incumbent as
-  known exactly, which is wrong under observation noise; NEI integrates over it
-  (Monte-Carlo over posterior samples of the objective at the observed points, with
-  common random numbers so the score is smooth for the optimizer).
-
-Both accept optional tuning through the constructor `args` dict:
+## Acquisition Function Contract
 
 ```python
-gpo = GPOptimizer(x_data, y_data, args={
-    "kg_reference_set_size": 100,   # cap on reference points (subsampled if larger)
-    "kg_seed": 0,                   # seed for the reference-set subsample
-    "nei_samples": 128,             # Monte-Carlo incumbent samples for NEI
-    "nei_reference_set_size": 100,
-    "nei_seed": 0,                  # seed -> common random numbers, smooth objective
-})
+def my_acquisition(x, gp_optimizer):
+    """
+    x : np.ndarray, shape (V, D) — candidate points to evaluate.
+    gp_optimizer : GPOptimizer — use posterior_mean() and posterior_covariance().
+    Returns: scores : np.ndarray, shape (V,) — HIGHER = more desirable (MAXIMIZED).
+    """
 ```
 
-Cost note: both are substantially heavier than `"variance"` or `"ucb"` — see caveat (g)
-above for what to do about it.
+- **Reading GP data:** Use `gpo.x_data` and `gpo.y_data` inside custom acquisitions (not `get_data()["original y data"]`). For `LogGPOptimizer` and `LogitGPOptimizer`, `gpo.y_data` is already in latent (`log` or `logit`) space matching `posterior_mean()`; any hardcoded threshold or target constant must also be transformed into latent space (`np.log(threshold)` or `scipy.special.logit((threshold - lo) / (hi - lo))` — see `transformed-optimizers-advanced`).
 
 ## Custom Acquisition Recipes
 
@@ -535,19 +314,14 @@ result = gpo.ask(
 | `method="global"\|"local"\|"hgdl"\|"hgdlAsync"` | Inner optimizer that searches for the argmax of the acquisition over `input_set`. `hgdl` requires `dask_client=`; `hgdlAsync` starts a background search and returns an `opt_obj` you can poll or `kill_client()`. |
 | `dask_client=client` | Distribute the inner optimization across Dask workers. |
 | `batch_size=B` | When candidates are a list, evaluate them in chunks of `B` on the cluster. |
+| `args={...}` | Updates `gpo.args` before evaluating the acquisition (`gp_optimizer_base.py:373, 536-537, 559`) — used by `"target probability"` (`{'a': lo, 'b': hi}`) and KG/NEI tuning keys. |
 | `max_iter`, `pop_size`, `info=True` | Inner optimizer controls. |
 
 `input_set` can be continuous bounds (`np.array([[lo,hi], ...])`), a list of candidate points (discrete finite set), or a list of arbitrary objects (non-Euclidean — strings, graphs — provided your kernel handles them).
 
-## Hyperparameter Coordination
-
-Acquisition functions don't add hyperparameters — they read the GP state via `gpo.posterior_mean()` and `gpo.posterior_covariance()`. However:
+## Hyperparameter & `args` Coordination
 
-- `gpo.y_data` reflects every `tell()` immediately, so `y_best` is current without extra bookkeeping. It is in the GP's modeling space (see the note at the top of this skill).
-- Acquisition functions receive no `args` dict — unlike kernel/mean/noise callables, the signature is exactly `f(x, gp_optimizer)`. Bind extra parameters with a closure or `functools.partial`.
-- The GP must be trained before acquisition makes sense — always call `train()` first
-- For `variance_only=True`: faster, returns just diagonal variances (usually what you want)
-- For full covariance: use `variance_only=False` but this is O(V²) memory
+Acquisition functions don't add hyperparameters — they read the trained GP state via `gpo.posterior_mean(x)` and `gpo.posterior_covariance(x, variance_only=True)` (prefer `variance_only=True` to avoid O(V²) covariance matrices). Custom acquisition callables always have the 2-argument signature `f(x, gp_optimizer)`; bind extra parameters with `functools.partial` or read `gp_optimizer.args` (see the `experiment-designer` skill for the full `args` dispatch table).
 
 ## Common Pitfalls
 
diff --git a/skills/acquisition-functions/references/choosing-an-acquisition-function.md b/skills/acquisition-functions/references/choosing-an-acquisition-function.md
index 6b41cea..e1c7538 100644
--- a/skills/acquisition-functions/references/choosing-an-acquisition-function.md
+++ b/skills/acquisition-functions/references/choosing-an-acquisition-function.md
@@ -28,12 +28,10 @@ measurement?"*
 
 **Why not plain EI by default.** Its EI anchors to `np.max(y_data)` — the largest
 *noisy observation*, not the posterior mean — so on noisy data it locks onto a spike;
-this is exactly what NEI was added to fix. And on gpCAM ≤ 8.4.2 there is a second,
-bigger reason: until the posterior mean somewhere exceeded that incumbent, EI was
-algebraically `0.3989 · σ(x)` over the whole sub-incumbent region — pure exploration,
-identical in ranking to `"variance"` there. Fixed in
-[#55](https://github.com/lbl-camera/gpCAM/issues/55), so check the version before
-ruling EI in or out on those grounds.
+this is exactly what NEI was added to fix. (In legacy gpCAM ≤ 8.4.2, plain EI also
+clipped `imp` before dividing by `σ`, collapsing to `0.3989 · σ(x)` wherever the
+posterior mean sat below the incumbent — fixed in
+[#55](https://github.com/lbl-camera/gpCAM/issues/55).)
 
 ```python
 # Explore first, then exploit — a good default for peak finding.
@@ -41,9 +39,9 @@ for i in range(N_ITERATIONS):
     acq = "variance" if i < 20 else "ucb"
     next_x = gpo.ask(input_set=parameter_bounds, acquisition_function=acq)["x"]
     y, nv = measure(next_x)
-    gpo.tell(next_x, np.array([y]), noise_variances=np.array([nv]), append=True)
+    gpo.tell(next_x, np.array([y]), noise_variances=np.array([nv]))
     if (i + 1) % RETRAIN_EVERY == 0:
-        gpo.train(hyperparameter_bounds=hp_bounds, method="local", max_iter=100)
+        gpo.train(hyperparameter_bounds=hp_bounds, method="local", max_iter=100, n_restarts=1)
 ```
 
 ---
@@ -106,11 +104,12 @@ on the same data.
 
 ```python
 def threshold_finder(x, gpo):
-    threshold = 0.5  # EDIT
+    threshold = 0.5  # EDIT — in the GP's modeling space
     mean = gpo.posterior_mean(x)["m(x)"]
     std = np.sqrt(np.maximum(
         gpo.posterior_covariance(x, variance_only=True)["v(x)"], 1e-10))
-    return std / (np.abs(mean - threshold) + 0.01)
+    eps = 0.01 * np.std(gpo.y_data) + 1e-12  # scale-aware softening
+    return std / (np.abs(mean - threshold) + eps)
 
 next_x = gpo.ask(input_set=parameter_bounds,
                  acquisition_function=threshold_finder)["x"]
@@ -126,14 +125,12 @@ want a usable map, alternate this with `"variance"`.
 Common in sample synthesis and alignment: not a maximum, a specification.
 
 ```python
-gpo = GPOptimizer(x_data, y_data, args={"a": 0.95, "b": 1.05})  # target band
 next_x = gpo.ask(input_set=parameter_bounds,
-                 acquisition_function="target probability")["x"]
+                 acquisition_function="target probability",
+                 args={"a": 0.95, "b": 1.05})["x"]  # or pass args={...} to GPOptimizer(...)
 ```
 
-`args` must carry `'a'` (lower) and `'b'` (upper) or the call raises. For
-`LogGPOptimizer` / `LogitGPOptimizer`, pass the bounds already transformed —
-`np.log(0.95)`, `np.log(1.05)`.
+`args` must carry `'a'` (lower) and `'b'` (upper) — passed to `GPOptimizer(..., args=...)` or `gpo.ask(..., args=...)` — or the call raises. For `LogGPOptimizer` / `LogitGPOptimizer`, pass the bounds in latent space (`np.log(a)` for `LogGPOptimizer`; `scipy.special.logit((a - lo) / (hi - lo))` for `LogitGPOptimizer` with `range=(lo, hi)`).
 
 ---
 
@@ -177,25 +174,28 @@ next_x = result["x"]        # shape (5, D)
 ```
 
 **Critical caveat.** With an array `input_set`, `n > 1`, `method != "hgdl"`, and a
-*string* acquisition, gpCAM replaces whatever was requested with `"total correlation"`
-and emits a warning (`gp_optimizer_base.py:517-525`). So `ask(n=5,
-acquisition_function="expected improvement")` does not run EI. Either ask for a
-batch-aware acquisition deliberately, or pass a callable with `method="hgdl"` and a
-Dask client:
+*string* acquisition other than `"total correlation"` or `"relative information entropy"`,
+gpCAM replaces whatever was requested with `"total correlation"` (or `"relative information
+entropy"` when `multi_task_mode="kronecker"`) and emits a warning
+(`gp_optimizer_base.py:583-594`). So `ask(n=5, acquisition_function="expected improvement")`
+does not run EI. Either ask for a batch-aware acquisition deliberately, or pass a callable
+with `method="hgdl"` and a Dask client:
 
 ```python
 from distributed import Client
-client = Client()
-result = gpo.ask(input_set=parameter_bounds, n=5,
-                 acquisition_function=my_callable,
-                 method="hgdl", dask_client=client)
+
+if __name__ == "__main__":
+    client = Client()
+    result = gpo.ask(input_set=parameter_bounds, n=5,
+                     acquisition_function=my_callable,
+                     method="hgdl", dask_client=client)
 ```
 
 ---
 
 ## 8. Multi-task — "each measurement returns a spectrum"
 
-Use `GPOptimizer` with `y_data` of shape (V, No) and a `multi_task_mode` (see the
+Use `GPOptimizer` with `y_data` of shape `(V, No)` and a `multi_task_mode` (see the
 `multi-task-advanced` skill), and `x_out` to pick outputs.
 
 **Ask back:** *"Are you mapping all the channels, or optimizing something derived from
@@ -204,12 +204,13 @@ them?"*
 - Mapping / exploring → `"relative information entropy set"` or `"variance"`
 - Optimizing → `"knowledge gradient"` or `"noisy expected improvement"`
 
-**Do not recommend plain EI here.** The `x_out` branch (`surrogate_model.py:319`) sums
-the posterior mean across tasks and sums the standard deviations across tasks, then
-compares that sum against `np.max(gpo.y_data)` of a scalar. Summed std is not the std
-of the sum — that needs the cross-task covariance — so the result is not expected
-improvement of anything. KG and NEI scalarize the same task-summed objective
-`g(x) = Σ_t f(x, t)` but build the joint posterior properly (`_scalarized_blocks`).
+**Why KG or NEI over plain multi-task EI.** The `x_out` branch (`surrogate_model.py:375-389`)
+scalarizes to the task-summed objective `g(x) = Σ_t f(x, t)` using the best task-summed
+posterior mean at the observed inputs (`_incumbent_task_sum`, `surrogate_model.py:415-427`)
+and `std = sqrt(Σ_t Var[f(x, t)])`. That spread drops the cross-task covariance, so plain
+multi-task EI is accurate only when tasks are uncorrelated. KG and NEI scalarize the same
+task-summed objective `g(x) = Σ_t f(x, t)` using the exact cross-task covariance
+(`_scalarized_blocks`, `surrogate_model.py:430-452`).
 
 ```python
 # Exploring all channels
@@ -237,7 +238,7 @@ argument. `evaluate_acquisition_function` divides the acquisition by the cost fu
 so an expensive move must earn its score.
 
 ```python
-def my_cost(origin, x, arguments=None):
+def my_cost(origin, x):
     """Minutes to move from `origin` to each row of `x`. Must be > 0."""
     return 1.0 + 4.0 * np.abs(x[:, 0] - origin[0])
 
@@ -248,7 +249,7 @@ next_x = gpo.ask(input_set=parameter_bounds,
 ```
 
 `position` must be passed to `ask()` — `evaluate_acquisition_function` only applies the
-cost when `origin is not None` (`surrogate_model.py:196`), so omitting it silently
+cost when `origin is not None` (`surrogate_model.py:264-265`), so omitting it silently
 disables cost-awareness. See the `cost-functions` skill for the full contract and for
 fitting cost parameters from observed timings.
 
diff --git a/skills/cost-functions/SKILL.md b/skills/cost-functions/SKILL.md
index 234dbf4..4b57a5e 100644
--- a/skills/cost-functions/SKILL.md
+++ b/skills/cost-functions/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: cost-functions
 description: Use when modeling the real expense of moving between gpCAM measurement points — motor travel time, settling, directional costs, sample damage, beam time, or zone-based penalties.
-gpcam_version: "8.4.x"
-fvgp_version: "4.8.x"
-last_verified: "2026-07-23 (gpCAM dadeb65)"
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: gpCAM Cost Functions
 
-*Verified against gpCAM 8.4.x / fvgp 4.8.x — last checked 2026-07-23 (gpCAM `dadeb65`).*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Design cost functions that account for the real expense of moving between measurement points — motor travel time, sample damage, beam time, etc.
 
@@ -53,128 +53,57 @@ def my_cost(origin, x):
 ```
 
 **Key rules:**
-- Cost must be **positive** (> 0) — zero cost causes division by zero
-- Cost is used as a **divisor** — higher cost = lower effective acquisition score
-- The function receives a batch of candidate points, not a single point
-- **The cost function does NOT receive the `args` dict.** Unlike kernel / prior mean /
-  noise callables, there is no arity-based dispatch here — the call site is literally
-  `cost_function(origin, x)`. A signature like `my_cost(origin, x, arguments=None)`
-  will run without error but `arguments` is always `None`. To parameterize a cost
-  function, use a closure or `functools.partial` (see below).
-- **`cost_function` only takes effect when `origin` is not None**, i.e. when you call
-  `ask(..., position=current_position)`. Without `position=`, the cost is ignored
-  entirely and the acquisition is used unmodified. This is the single most common
-  reason a cost function "does nothing".
+- Cost must be **strictly positive** (`> 0`) and return a 1-D array of shape `(V,)` — it is used as a divisor (`surrogate_model.py:265`), so zero cost causes division by zero. Always include a positive base offset (e.g., measurement overhead time).
+- **Pass `position=current_position` to `ask()`:** gpCAM only evaluates `cost_function(origin, x)` when `origin is not None` (`surrogate_model.py:264`). Omitting `position=` in `ask()` silently disables the cost function.
+- **No `args` dict or trainable hyperparameters:** `cost_function` is always called with 2 arguments `cost_function(origin, x)` and is not optimized during `train()`. Parameterize it with `functools.partial` (which survives `pickle.dumps(gpo)` for checkpointing; see the `experiment-designer` skill for the full `args` dispatch table).
 
 ## Recipes
 
-### L2 (Euclidean) Distance Cost
-Simple travel time proportional to straight-line distance:
+### Euclidean / Anisotropic Travel Cost (with `functools.partial`)
 ```python
-def l2_cost(origin, x):
-    """Cost proportional to Euclidean distance."""
-    offset = 1.0   # minimum cost (prevents div-by-zero, represents measurement time)
-    speed = 1.0     # cost per unit distance
-    distance = np.linalg.norm(x - origin, axis=1)
-    return offset + speed * distance
-```
-
-### L1 (Manhattan) Distance Cost
-For stage systems that move one axis at a time:
-```python
-def l1_cost(origin, x):
-    """Cost proportional to Manhattan distance (axis-by-axis motion)."""
-    offset = 1.0
-    speed = 1.0
-    distance = np.sum(np.abs(x - origin), axis=1)
-    return offset + speed * distance
-```
+from functools import partial
+import numpy as np
 
-### Anisotropic Cost
-Different axes have different speeds. Since the cost function receives no `args`,
-bind the speeds with a factory:
-```python
-def make_anisotropic_cost(speeds, offset=1.0):
+def travel_cost(origin, x, offset=1.0, weights=1.0):
     """
-    Returns a two-argument cost function with per-axis speeds captured.
-    speeds : array of shape (D,) — cost per unit distance along each axis.
+    Base overhead `offset` (> 0) plus weighted Euclidean travel time.
+    Pass `weights` as a scalar or shape-(D,) per-axis array via functools.partial.
     """
-    speeds = np.asarray(speeds, dtype=float)
+    diff = (x - origin) * np.asarray(weights, dtype=float)
+    return offset + np.linalg.norm(diff, axis=1)
 
-    def anisotropic_cost(origin, x):
-        return offset + np.sum(np.abs(x - origin) * speeds, axis=1)
-    return anisotropic_cost
+my_cost = partial(travel_cost, offset=1.0, weights=[1.0, 3.0])
 ```
+For axis-by-axis (Manhattan) stages, replace `np.linalg.norm(diff, axis=1)` with `np.sum(np.abs(diff), axis=1)`.
 
-### Directional Cost
-Moving in one direction is cheaper (e.g., gravity-assisted, or always-increasing scans):
+### Directional, Settling, or Zone-Based Cost
+Combine asymmetric step costs, step-function settling times, or expensive-region penalties into the same positive 1-D return vector:
 ```python
-def directional_cost(origin, x):
-    """Cheaper to move in +x direction than -x."""
-    offset = 1.0
-    diff = x - origin
-    # Forward motion (positive) is cheap, backward is expensive
-    forward_cost = np.maximum(diff[:, 0], 0) * 1.0   # cost going forward
-    backward_cost = np.maximum(-diff[:, 0], 0) * 5.0  # 5x cost going backward
-    lateral_cost = np.sum(np.abs(diff[:, 1:]), axis=1) * 1.0
-    return offset + forward_cost + backward_cost + lateral_cost
-```
-
-### Settling Time Cost
-Fast moves need more settling time:
-```python
-def settling_cost(origin, x):
-    """
-    Short moves are fast; long moves need extra settling time.
-    Models: cost = base + travel + settling * (distance > threshold)
-    """
-    base = 1.0
-    travel_rate = 0.5
-    settle_time = 3.0
-    settle_threshold = 2.0  # distance above which settling kicks in
-    
-    distance = np.linalg.norm(x - origin, axis=1)
-    settling = np.where(distance > settle_threshold, settle_time, 0.0)
-    return base + travel_rate * distance + settling
-```
-
-### Zone-Based Cost
-Some regions of the parameter space are more expensive:
-```python
-def zone_cost(origin, x):
-    """
-    Higher cost to measure in certain zones.
-    E.g., cryogenic sample region requires cooldown.
-    """
+def stage_cost(origin, x):
+    """Asymmetric axis-0 motion + long-move settling time + zone penalty."""
     base = 1.0
-    distance = np.linalg.norm(x - origin, axis=1)
-    
-    # Expensive zone: x[:, 0] > 8.0
+    diff = x - origin
+    axis0_cost = np.maximum(diff[:, 0], 0.0) * 1.0 + np.maximum(-diff[:, 0], 0.0) * 5.0
+    settling = np.where(np.linalg.norm(diff, axis=1) > 2.0, 3.0, 0.0)
     zone_penalty = np.where(x[:, 0] > 8.0, 10.0, 0.0)
-    
-    return base + distance + zone_penalty
+    return base + axis0_cost + settling + zone_penalty
 ```
 
-## Usage
+## Usage in the Loop
 
 ```python
 gpo = GPOptimizer(
     x_data=x_data,
     y_data=y_data,
-    cost_function=l2_cost,
-    # cost_function=make_anisotropic_cost([1.0, 3.0, 0.5]),   # parameterized variant
+    cost_function=my_cost,
 )
-```
 
-**You must pass `position=` to `ask()` for the cost to be applied:**
-
-```python
 current_position = np.array([[0.0, 0.0]])       # where the instrument is now
 
 for i in range(N_ITERATIONS):
     next_x = gpo.ask(
         input_set=parameter_bounds,
-        position=current_position,               # REQUIRED — without it cost is ignored
+        position=current_position,               # REQUIRED — sets `origin`
         acquisition_function="variance",
     )["x"]
     y = measure(next_x)
@@ -182,35 +111,6 @@ for i in range(N_ITERATIONS):
     current_position = next_x                    # advance the origin
 ```
 
-Forgetting `position=` is the usual reason a cost function appears to have no effect:
-gpCAM applies the cost only when `origin is not None`, and silently skips it otherwise.
-
-## Parameterizing a Cost Function
-
-The `args` dict passed to `GPOptimizer(..., args=...)` reaches kernel, prior mean, and
-noise functions — **not** the cost function. Use a closure (as in
-`make_anisotropic_cost` above) or `functools.partial`:
-
-```python
-from functools import partial
-
-def scaled_l2_cost(origin, x, offset, speed):
-    return offset + speed * np.linalg.norm(x - origin, axis=1)
-
-gpo = GPOptimizer(x_data, y_data,
-                  cost_function=partial(scaled_l2_cost, offset=1.0, speed=2.5))
-```
-
-`partial` of a module-level function pickles by reference, so this form survives
-`pickle.dumps(gpo)` for checkpointing; a closure defined inside another function does
-not.
-
-## Cost Functions Do NOT Add Hyperparameters
-
-Unlike kernel/mean/noise functions, cost function parameters are **fixed** — they are
-not optimized during training. To tune them, refit them yourself and rebuild the
-callable.
-
 ## Calibrating a Cost Function From Observed Moves
 
 If you don't know the cost parameters ahead of time, record the real cost of each move
@@ -238,7 +138,7 @@ keep = np.abs(rate - rate.mean()) <= 2.0 * rate.std()
                                distances[keep], costs[keep], p0=[1.0, 1.0])
 print(f"offset={offset:.3f} s, speed={speed:.3f} s per unit distance")
 
-gpo.cost_function = partial(scaled_l2_cost, offset=offset, speed=speed)
+gpo.cost_function = partial(travel_cost, offset=offset, weights=speed)
 ```
 
 Fitting on measured seconds makes the cost units meaningful, which matters because the
diff --git a/skills/experiment-designer/SKILL.md b/skills/experiment-designer/SKILL.md
index 4327ad5..e0a3ae9 100644
--- a/skills/experiment-designer/SKILL.md
+++ b/skills/experiment-designer/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: experiment-designer
 description: Use for end-to-end autonomous experiment design with gpCAM. Translates a scientist's description of their measurement into a complete, runnable gpCAM script — useful for replacing raster scans with adaptive sampling, peak-finding, or parameter optimization.
-gpcam_version: "8.6.x"
-fvgp_version: "4.9.1"
-last_verified: "2026-10-14 (fvgp 4.9.1)"
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: gpCAM Experiment Designer
 
-*Verified against gpCAM 8.6.x / fvgp 4.9.1 — last checked 2026-10-14.*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Design complete autonomous experiment scripts using gpCAM. You translate a scientist's description of their measurement into a runnable Python script.
 
@@ -20,13 +20,6 @@ When a user wants to:
 - Find optimal experimental conditions (peak finding, parameter optimization)
 - Explore a parameter space efficiently
 
-## Your Role
-
-You are helping **beamline scientists** who may not know GP math or the gpCAM API. Your job is to:
-1. Understand their experiment (what they measure, what they control, what they want to find)
-2. Generate a complete, well-commented Python script they can adapt
-3. Explain the key choices you made in plain language
-
 ## Conversation Flow
 
 ### Step 1: Understand the Experiment
@@ -55,7 +48,7 @@ Based on their answers, decide:
 
 | Choice | Guidance |
 |--------|----------|
-| **Optimizer class** | Pick by the support of the observations: `GPOptimizer` for unconstrained or negative-allowed `y`; `LogGPOptimizer` if `y > 0` (intensities, rates, concentrations); `LogitGPOptimizer` if `y ∈ [0, 1]` (yields, fractions, probabilities). A plain GP on positive-only or bounded data can predict invalid values — see `transformed-optimizers-advanced` skill. |
+| **Optimizer class** | Pick by the support and dimensionality of the observations: `GPOptimizer` handles both scalar (`y_data` shape `(N,)`) and multi-output/multi-task (`y_data` shape `(N, No)`) problems via `multi_task_mode` (`"independent"` default, `"cartesian"` for joint task correlations or missing task entries, `"kronecker"` for separable product kernels — see `multi-task-advanced` skill; legacy `fvGPOptimizer` is deprecated since 8.7.0); `LogGPOptimizer` if scalar `y > 0` (intensities, rates, concentrations); `LogitGPOptimizer` if scalar `y ∈ (a, b)` (`range=(0.0, 1.0)` default — yields, fractions, bounded intervals). A plain GP on positive-only or bounded data can predict invalid values — see `transformed-optimizers-advanced` skill. |
 | **Kernel** | Default Matérn-3/2 ARD is good for most cases. Use periodic kernel if periodicity is known. Use Matérn-1/2 for rough/discontinuous data, Matérn-5/2 or SE for very smooth. See `kernel-designer` skill for custom kernels. |
 | **Acquisition function** | **Do not pick this one on your own — read the `acquisition-functions` skill and confirm the choice with the scientist.** It's the argument that decides where the instrument actually goes, and several built-ins have gpCAM-specific behavior that surprises people (plain EI is maximization-only, and on gpCAM ≤ 8.4.2 it degenerates to pure exploration; `ask(n>1)` silently overrides string acquisitions). Starting points, all subject to confirmation: `'variance'` for exploration/mapping; `'ucb'` for maximization (`beta` fixed at 3.0 — write a callable to tune it); `'lcb'` for **minimization**; `'noisy expected improvement'` when measurements are noisy, or `'knowledge gradient'` for lookahead that keeps exploring when EI stalls — both also work multi-task; `'gradient'` or a radical-gradient callable to map where the signal changes; a custom callable for multi-objective or constraints. Do **not** default to `'expected improvement'`. |
 | **Prior mean** | Default is a **constant equal to `mean(y_data)`** — not zero. Away from data the posterior reverts to that constant. Override with `prior_mean_function=` only if they have a physical model. See `prior-mean-functions` skill. |
@@ -111,7 +104,7 @@ Note `acq_func='lcb'` is the default, which **minimizes** `func` — pass
 iteration, and `break_condition=f(x_data, y_data)` returning `True` to stop early.
 
 Use it when the scientist has a simulator or instrument wrapper they can hand gpCAM as
-a Python function. For fvGP/multi-task, pass `x_out=np.array([...])` and note the
+a Python function. For multi-task `GPOptimizer` (2-D `y_data`), pass `x_out=np.array([...])` and note the
 return shape differs — see the `multi-task-advanced` skill.
 
 #### B. Full template (adaptive loop)
@@ -121,7 +114,7 @@ Output a complete Python script with this structure:
 ```python
 """
 Autonomous Experiment: [description]
-Generated for gpCAM v8.4.x
+Generated for gpCAM v8.7.x
 
 Input space: [dimensions and ranges]
 Output: [what is measured]
@@ -324,7 +317,7 @@ if __name__ == "__main__":
 4. **Hyperparameter bounds matter** — set sensible defaults based on the parameter ranges and expected signal scale.
 5. **Default kernel is usually fine** — only suggest custom kernels when there's a clear reason (known periodicity, symmetry, etc.).
 6. **The acquisition function is not a default — it's a conversation.** State your recommendation, name the one tradeoff it makes, and confirm before writing the loop. Record the reason as a comment in the script. Remind them it's an `ask()` argument, so it can change between iterations — "explore first, then exploit" is often the right design. See the `acquisition-functions` skill.
-7. **Training schedule** — train globally once at the start, then locally every N iterations. Don't train every iteration (too slow).
+7. **Training schedule** — run multi-restart local training (`gpo.train(hyperparameter_bounds=hp_bounds)`, which defaults to `method="local", n_restarts=5`) once at the start, then warm-started single-start local training (`method="local", n_restarts=1`) every N iterations. Don't train every iteration (too slow).
 8. **Initial points** — always start with random initial sampling before the adaptive loop.
 
 ## Hyperparameter Coordination
@@ -332,9 +325,9 @@ if __name__ == "__main__":
 This is critical and often the source of bugs:
 
 - The hyperparameter vector is shared across kernel, mean, and noise functions, and the index ranges each one reads **must be disjoint**. fvGP assumes a hyperparameter belonging to the mean has zero kernel derivative and vice versa, so an overlap corrupts the training gradients as well as the model.
-- **Standard ordering: kernel, then mean, then noise.** Every gpCAM skill uses this layout.
-- For the default kernel with D input dimensions: `hps[0]` = signal variance, `hps[1:D+1]` = length scales. So custom mean/noise hyperparameters start at `K = 1 + D`.
-- **Derive `K = 1 + D`; never hardcode it.** `K = 3` is correct only for 2-D input. In 1-D it aliases a length scale; in 3-D two components silently share a hyperparameter.
+- **Standard ordering: kernel (`hps[:K]`), then prior mean (`hps[K : K + M]`), then noise (`hps[K + M : K + M + L]`).** Every gpCAM skill uses this layout.
+- For the single-task (or `"kronecker"`) default kernel with `D` input dimensions: `hps[0]` = signal variance, `hps[1:D+1]` = length scales, so `K = 1 + D`. For a custom kernel, `K` is the number of parameters that kernel consumes. For multi-task default kernels, `K = 1 + (D + 1) = D + 2` under `multi_task_mode="cartesian"` (where `x` inside kernel/mean/noise callables has shape `(N, D + 1)` with the task index in the last column) and `K = No * (1 + D)` under `multi_task_mode="independent"`.
+- **Derive `K` from `D` (or `x.shape[1]`); never hardcode it.** `K = 3` is correct only for 2-D single-task input. In 1-D it aliases a length scale; in 3-D two components silently share a hyperparameter. (Under `multi_task_mode="cartesian"`, remember `x.shape[1]` inside callables is `D + 1`.)
 - **Never read `hps[-1]` in a prior mean function** — the end of the vector belongs to the noise function.
 - Assert `len(gpo.hyperparameters) == len(hp_bounds)` after construction; it catches most layout mistakes for one line.
 - **Always document the hyperparameter layout** in a comment at the top of the script
@@ -350,7 +343,7 @@ Example with custom noise:
 # Total: D + 2 hyperparameters
 
 def my_noise(x, hps):
-    K = 1 + x.shape[1]              # = D + 1; derived, not hardcoded
+    K = 1 + x.shape[1]              # = D + 1 for single-task; derived, not hardcoded
     return np.full(len(x), hps[K]**2)
 
 hp_bounds = np.array(
@@ -362,7 +355,7 @@ hp_bounds = np.array(
 
 ## Advanced Options (mention only if needed)
 
-- **Incremental data updates**: `gpo.tell(x_new, y_new, append=True)` adds points without overwriting; `append=False` replaces. Default is to replace — use `append=True` in a streaming instrument loop.
+- **Incremental data updates**: `gpo.tell(x_new, y_new)` appends new points by default (`append=True`); pass `append=False` to replace the entire dataset instead.
 - **Async training**: `opt_obj = gpo.train(..., asynchronous=True, method="hgdl", dask_client=client)` returns immediately; later call `gpo.update_hyperparameters(opt_obj)` to pull current best hyperparameters, and `gpo.stop_training(opt_obj)` to finish. Useful when training is expensive and the loop shouldn't block.
 - **Async ask**: `method="hgdlAsync"` starts a background search for the next point; the result dict contains an `opt_obj` you can `kill_client()` once you've used the suggestion.
 - **Checkpointing**: `GPOptimizer` instances are picklable — `pickle.dumps(gpo)` before a long run lets you reload state later.
@@ -371,7 +364,7 @@ hp_bounds = np.array(
 ### User-function arguments (`args`)
 
 `GPOptimizer(..., args={"a": 1.5})` plumbs a dict to **kernel, prior mean, and noise
-functions only**. Dispatch is by **parameter count** (`inspect.signature`), not by
+functions** (and `gpo.ask(..., args=...)` updates `gpo.args` before evaluating acquisitions). Dispatch for kernel, mean, and noise is by **parameter count** (`inspect.signature`), not by
 name — declare one extra positional parameter and gpCAM passes the dict; omit it and
 it doesn't. The parameter name is yours to choose.
 
@@ -381,7 +374,7 @@ it doesn't. The parameter name is yours to choose.
 | Prior mean | `f(x, hps)` | `f(x, hps, args)` |
 | Noise | `f(x, hps)` | `f(x, hps, args)` |
 | **Cost** | `f(origin, x)` | **not supported** — see below |
-| **Acquisition** | `f(x, gp_obj)` | **not supported** — see below |
+| **Acquisition** | `f(x, gp_obj)` | **not supported as 3rd param** — read `gp_obj.args` inside `f` |
 
 ```python
 def my_kernel(x1, x2, hps, args):        # 4 params -> gets args
diff --git a/skills/gp2scale-advanced/SKILL.md b/skills/gp2scale-advanced/SKILL.md
index 7e2afdd..d4bbc2d 100644
--- a/skills/gp2scale-advanced/SKILL.md
+++ b/skills/gp2scale-advanced/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: gp2scale-advanced
 description: Use for large-scale gpCAM experiments (>10k points up to millions) using sparse compactly-supported kernels and Dask distributed computing for exact GP computation at scale.
-gpcam_version: "8.6.x"
-fvgp_version: "4.9.1"
-last_verified: "2026-10-14 (fvgp 4.9.1)"
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: gp2Scale — Large-Scale GPs
 
-*Verified against gpCAM 8.6.x / fvgp 4.9.1 — last checked 2026-10-14.*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Design experiments with tens of thousands to millions of data points using gpCAM's gp2Scale mode for exact GP computation at scale.
 
@@ -30,28 +30,34 @@ gp2Scale uses:
 ## Basic Setup
 
 ```python
+import numpy as np
 from distributed import Client
 from gpcam import GPOptimizer
 
-# Start a local Dask cluster
-client = Client()                  # uses all available cores
-client.wait_for_workers(4)         # good practice: wait for workers before constructing
-
-gpo = GPOptimizer(
-    x_data=x_data,
-    y_data=y_data,
-    gp2Scale=True,
-    dask_client=client,
-    # gp2Scale_batch_size defaults to 10000 and gp2Scale_distribution to "blockwise".
-    # Leave both alone unless profiling says otherwise — the batch size means
-    # different things under the two distributions. See "Choosing the Batch Size".
-    init_hyperparameters=np.array([0.73, 0.0014]),  # signal var, length scale
-)
-
-# Under gp2Scale, train() runs MCMC: its stochastic log-determinant rules out the
-# gradient-based default. max_iter=25 keeps this example quick; for a real fit leave it at
-# the default and let the convergence check end the chain (from 400 steps on).
-gpo.train(hyperparameter_bounds=hps_bounds, max_iter=25, info=True)
+if __name__ == "__main__":
+    # Start a local Dask cluster (must be inside __main__ guard for worker spawning)
+    client = Client()                  # uses all available cores
+    client.wait_for_workers(4)         # good practice: wait for workers before constructing
+
+    D = x_data.shape[1]
+    init_hps = np.array([0.73] + [0.01] * D)
+    hps_bounds = np.array([[0.01, 10.0]] + [[1e-4, 0.1]] * D)
+
+    gpo = GPOptimizer(
+        x_data=x_data,
+        y_data=y_data,
+        gp2Scale=True,
+        dask_client=client,
+        # gp2Scale_batch_size defaults to 10000 and gp2Scale_distribution to "blockwise".
+        # Leave both alone unless profiling says otherwise — the batch size means
+        # different things under the two distributions. See "Choosing the Batch Size".
+        init_hyperparameters=init_hps,  # signal var + D length scales
+    )
+
+    # Under gp2Scale, train() runs MCMC: its stochastic log-determinant rules out the
+    # gradient-based default. max_iter=25 keeps this example quick; for a real fit leave it at
+    # the default and let the convergence check end the chain (from 400 steps on).
+    gpo.train(hyperparameter_bounds=hps_bounds, max_iter=25, info=True)
 ```
 
 ## Choosing the Batch Size
@@ -120,15 +126,20 @@ gpo.get_gp2Scale_exec_time(time_per_worker_execution=0.05, number_of_workers=32)
 
 ## Kernel Requirement
 
-When `gp2Scale=True`, the kernel MUST produce a sparse matrix. The default switches to an anisotropic Wendland kernel automatically if no custom kernel is provided.
+When `gp2Scale=True`, the kernel MUST produce a sparse matrix. If `kernel_function=None`, gpCAM automatically uses the default compactly supported Wendland kernel (`gp_optimizer.py:126-131, 224-225`).
+
+`gpcam.kernels` re-exports all `fvgp.kernels` Wendland helpers (`gpcam/kernels.py:1`, `gp_optimizer.py:128-130`):
+- `wendland_anisotropic_gp2Scale_cpu(x1, x2, hps)` — CPU Wendland kernel for `gp2Scale`
+- `wendland_anisotropic_gp2Scale_cpu_sparse(x1, x2, hps)` — builds a `scipy.sparse` block directly without materializing the dense array
+- `wendland_anisotropic_gp2Scale_gpu(x1, x2, hps)` — GPU Wendland kernel for `gp2Scale`
+- `wendland_anisotropic(x1, x2, hps)` and the distance-matrix primitive `wendland_kernel(d, l)`
 
-If providing a custom kernel, it must have compact support:
 ```python
-from gpcam.kernels import wendland_anisotropic
+from gpcam.kernels import wendland_anisotropic_gp2Scale_cpu_sparse
 
 def my_gp2scale_kernel(x1, x2, hps):
-    """Custom kernel with compact support for gp2Scale."""
-    return wendland_anisotropic(x1, x2, hps)
+    """Custom kernel with compact support returning a sparse block for gp2Scale."""
+    return wendland_anisotropic_gp2Scale_cpu_sparse(x1, x2, hps)
 ```
 
 ## Hyperparameters for Wendland Kernel
@@ -219,13 +230,14 @@ To check convergence properly, run several chains and compare them with
 from dask_jobqueue import SLURMCluster
 from distributed import Client
 
-cluster = SLURMCluster(
-    cores=32,
-    memory="64GB",
-    walltime="01:00:00",
-)
-cluster.scale(jobs=4)  # 4 nodes × 32 cores
-client = Client(cluster)
+if __name__ == "__main__":
+    cluster = SLURMCluster(
+        cores=32,
+        memory="64GB",
+        walltime="01:00:00",
+    )
+    cluster.scale(jobs=4)  # 4 nodes × 32 cores
+    client = Client(cluster)
 ```
 
 ## Common Pitfalls
diff --git a/skills/kernel-designer/SKILL.md b/skills/kernel-designer/SKILL.md
index ae51b82..eeafc14 100644
--- a/skills/kernel-designer/SKILL.md
+++ b/skills/kernel-designer/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: kernel-designer
-description: Use when designing or composing custom kernel (covariance) functions for gpCAM that encode domain knowledge — smoothness, periodicity, symmetry, anisotropy, or non-Euclidean input spaces.
-gpcam_version: "8.4.x"
-fvgp_version: "4.8.x"
-last_verified: "2026-07-23 (gpCAM dadeb65); multi-task line updated for gpCAM 8.7.0 on 2026-10-28"
+description: Use when designing or composing custom kernel (covariance) functions for gpCAM — smoothness, periodicity, symmetry, anisotropy, non-stationary amplitude/length-scale kernels (`non_stat_kernel`, Gibbs), non-Euclidean inputs, or neural-network-warped deep kernels (`Network`).
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: gpCAM Kernel Designer
 
-*Verified against gpCAM 8.4.x / fvgp 4.8.x — last checked 2026-07-23 (gpCAM `dadeb65`).*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Design custom kernel (covariance) functions for gpCAM that encode domain knowledge about the experiment.
 
@@ -42,10 +42,9 @@ def my_kernel(x1, x2, hyperparameters):
     """
 ```
 
-- `x1` and `x2` are 2D arrays even for 1D inputs
-- The output must be an `(N1, N2)` matrix
-- The kernel must be symmetric: `k(x1, x2) = k(x2, x1).T`
-- The kernel must be positive semi-definite
+- `x1` and `x2` are 2D arrays even for 1D inputs (for `multi_task_mode="cartesian"`, `x1` and `x2` have shape `(N, D + 1)` with the task index in the last column — see `multi-task-advanced`)
+- Declare a 4th positional parameter `my_kernel(x1, x2, hyperparameters, args)` to receive `GPOptimizer(..., args=...)` (`gp_optimizer.py:113-114`; see `experiment-designer` for the full `args` dispatch table)
+- The output must be an `(N1, N2)` symmetric positive semi-definite matrix (`k(x1, x2) = k(x2, x1).T`)
 - Use vectorized numpy operations — avoid Python loops over data points
 
 ## Building Blocks
@@ -344,7 +343,7 @@ hp_bounds = np.array([[0.01, 10.0], [0.01, 10.0]] + [[-2.0, 2.0]] * N_NN)
 
 gp = GPOptimizer(x_data, y_data, init_hyperparameters=init_hps,
                  kernel_function=deep_kernel)
-gp.train(hyperparameter_bounds=hp_bounds, method="global", max_iter=100)
+gp.train(hyperparameter_bounds=hp_bounds, method="adam", max_iter=100)
 ```
 
 Sanity check that the unpacking is live — evaluate the kernel at two different NN hyperparameter draws and confirm the matrices differ:
@@ -369,31 +368,11 @@ Training notes: `method="mcmc"` or `method="adam"` scale better than `global`/`l
 
 ## Hyperparameter Coordination
 
-**Critical:** When designing a custom kernel, you must:
+Kernel hyperparameters occupy the leading slice `hps[:K]`, followed by prior-mean (`hps[K : K + M]`) and noise (`hps[K + M :]`) parameters — see `experiment-designer` ("Hyperparameter Coordination") for the canonical index convention.
 
-1. **Document the hyperparameter layout** — which index maps to what
-2. **Set matching bounds** — the `hyperparameter_bounds` array must have one row per hyperparameter
-3. **Coordinate with noise/mean functions** — all three share the same hyperparameter vector
-
-```python
-# ALWAYS include a comment block like this:
-#
-# Hyperparameter layout for my_kernel:
-# hps[0]     = signal variance          bounds: [0.001, 100]
-# hps[1]     = length scale dim 0       bounds: [0.01, range_dim0 * 10]
-# hps[2]     = length scale dim 1       bounds: [0.01, range_dim1 * 10]
-# hps[3]     = period (if periodic)     bounds: [expected_period * 0.5, expected_period * 2]
-```
-
-### Setting Initial Hyperparameters
-- **Signal variance**: start near `np.var(y_data)` or `1.0`
-- **Length scales**: start near `0.1 * range(x_dim)` — not too small (overfitting) or too large (underfitting)
-- **Period**: start near the expected period if known
-
-### Setting Bounds
-- **Signal variance**: `[0.001, 10 * np.std(y_data)]`
-- **Length scales**: `[0.01, 10 * range(x_dim)]` — lower bound prevents overfitting
-- **Period**: `[0.5 * expected, 2 * expected]` — tight if well-known
+- **Signal variance**: initialize near `np.var(y_data)` (or `1.0`), bounds `[0.001, 10 * np.var(y_data)]`
+- **Length scales**: initialize near `0.1 * range(x_dim)`, bounds `[0.01 * range(x_dim), 10 * range(x_dim)]`
+- **Period**: initialize near the expected period, bounds `[0.5 * expected, 2 * expected]`
 
 ## Common Pitfalls
 
diff --git a/skills/multi-task-advanced/SKILL.md b/skills/multi-task-advanced/SKILL.md
index f0d0ced..e4017e1 100644
--- a/skills/multi-task-advanced/SKILL.md
+++ b/skills/multi-task-advanced/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: multi-task-advanced
-description: Use for multi-output, vector-valued, or function-valued gpCAM experiments using GPOptimizer with several output columns and multi_task_mode — useful when a single measurement returns multiple correlated quantities (e.g., spectra, multi-channel detectors).
+description: Use for multi-output, vector-valued, or function-valued gpCAM experiments using GPOptimizer with several output columns and multi_task_mode (or migrating legacy fvGPOptimizer code) — useful when a single measurement returns multiple correlated quantities (e.g., spectra, multi-channel detectors).
 gpcam_version: "8.7.x"
 fvgp_version: "4.11.x"
-last_verified: "2026-10-28 (gpCAM 8.7.0)"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: Multi-Task GPs with GPOptimizer
 
-*Verified against gpCAM 8.7.x / fvgp 4.11.x — last checked 2026-10-28.*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Design experiments with vector-valued or function-valued outputs using gpCAM's multi-task GP.
 
@@ -68,23 +68,8 @@ gpo.train()
 gpo = GPOptimizer(x_data, y_data, multi_task_mode="kronecker")
 gpo.train()
 
-# Cartesian with a custom kernel — supply init_hyperparameters and hp bounds.
-# The layout below matches `multi_task_kernel` in "Multi-Task Kernel Design" below:
-# D + 2 hyperparameters for D input dimensions. (The default kernel in this mode treats
-# the task index as one more input dimension — fine for tests, rarely for science.)
-gpo = GPOptimizer(
-    x_data=x_data,
-    y_data=y_data,
-    init_hyperparameters=np.ones(4) / 10.0,
-    kernel_function=multi_task_kernel,
-    multi_task_mode="cartesian",
-)
-gpo.train(hyperparameter_bounds=np.array([
-    [0.01, 100.0],  # hps[0] signal variance
-    [0.01, 10.0],   # hps[1] length scale dim 0
-    [0.01, 10.0],   # hps[2] length scale dim 1
-    [0.0, 1.0],     # hps[3] task correlation  (index D+1, D=2)
-]))
+# For multi_task_mode="cartesian" with a custom task-correlation kernel, see
+# "Multi-Task Kernel Design" below for the complete definition and bounds.
 ```
 
 ### Predictions and ask — the `x_out` argument
@@ -110,17 +95,12 @@ gpo.ask(parameter_bounds, x_out=np.array([0, 1]), n=4,
 ```
 
 **Plain `"expected improvement"` with `x_out` assumes the tasks are uncorrelated.** The
-multi-task branch of `evaluate_gp_acquisition_function` scalarizes to the task-summed
-objective and takes the spread as `sqrt(Σ_t Var[f(x,t)])`, which drops the cross-task
-covariance. Its incumbent is the best task-summed *posterior mean* at the data points
-(gpCAM ≥ 8.7; before, a missing output made every EI value NaN). It is cheap and coherent but not exact; for correlated tasks — the usual
-reason to fit a multi-task GP at all — prefer KG/NEI below.
-
-**On gpCAM ≤ 8.4.2, do not use plain `"expected improvement"` with `x_out` at all.**
-There it *summed the standard deviations* rather than the variances and compared the
-task-sum against `np.max(gpo.y_data)`, a single `(point, task)` entry of the flattened
-product-space vector — quantities that aren't comparable, with no error raised. Fixed in
-[#55](https://github.com/lbl-camera/gpCAM/issues/55).
+multi-task branch of `evaluate_gp_acquisition_function` (`surrogate_model.py:375-389, 415-427`)
+scalarizes to the task-summed objective and takes the spread as `sqrt(Σ_t Var[f(x,t)])`, which
+drops the cross-task covariance. Its incumbent is the best task-summed *posterior mean* at the
+data points (`_incumbent_task_sum`; before 8.7, a missing output made every EI value NaN, and on
+gpCAM ≤ 8.4.2 it summed standard deviations instead of variances — fixed in #55). For correlated
+tasks — the usual reason to fit a multi-task GP at all — prefer KG/NEI below.
 
 Choose by what the scientist is doing:
 
@@ -238,7 +218,7 @@ while reading a length scale as a task correlation.
 
 ## Important Notes
 
-1. **Multi-task acquisition**: Use `"relative information entropy set"` / `"relative information entropy"` / `"variance"` / `"total correlation"` for batch acquisition across tasks. For optimization, `"expected improvement"`, `"knowledge gradient"`, and `"noisy expected improvement"` all work multi-task — they act on the **task-summed** objective `sum_t f(x, t)` (`"noisy expected improvement"` is the one to reach for when measurements are noisy). Pass `x_out=...` to the `ask()` call. Custom callables are supported and often advisable when you care about a specific task or combination rather than the plain sum.
+1. **Multi-task acquisition**: Use `"relative information entropy set"` / `"relative information entropy"` / `"variance"` / `"total correlation"` for batch acquisition across tasks (`"total correlation"` is unavailable in `"kronecker"` mode). For optimization of the task-summed objective `sum_t f(x, t)`, prefer `"knowledge gradient"` or `"noisy expected improvement"`, which use the exact cross-task covariance (whereas multi-task `"expected improvement"` assumes independent tasks). Pass `x_out=...` to `ask()`, or write a custom callable if you care about a weighted combination or a single task.
 2. **Missing task observations** (`"cartesian"` only): `y_data` can have `np.nan` entries (e.g., task 1 wasn't measured at some x); the corresponding `noise_variances` entry **must also** be `np.nan`. The GP just ignores those entries — no imputation needed. `"independent"` and `"kronecker"` reject NaN.
 3. **Default kernel**: `"independent"` uses the anisotropic Matérn on the inputs; `"kronecker"` the coregionalization kernel (a Matérn on the inputs times a learned output covariance), with data-driven defaults and bounds; `"cartesian"` treats the task index as one more input dimension, which is rarely what a scientist wants — write a kernel (below). If you supply a custom `kernel_function`, you become responsible for the full `init_hyperparameters` + `hyperparameter_bounds` layout.
 4. **Deep kernel via NN warping**: For harder multi-task structure, `from gpcam.deep_kernel_network import Network` gives you an MLP you parametrize from `hps` — see the `kernel-designer` skill for the pattern.
diff --git a/skills/noise-functions/SKILL.md b/skills/noise-functions/SKILL.md
index 4904230..738270f 100644
--- a/skills/noise-functions/SKILL.md
+++ b/skills/noise-functions/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: noise-functions
 description: Use when modeling position-dependent, heteroscedastic, or otherwise structured noise in gpCAM — e.g., detector characteristics, count-rate-dependent variance, or non-uniform measurement uncertainty.
-gpcam_version: "8.4.x"
-fvgp_version: "4.8.x"
-last_verified: "2026-07-23 (gpCAM dadeb65)"
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: gpCAM Noise Functions
 
-*Verified against gpCAM 8.4.x / fvgp 4.8.x — last checked 2026-07-23 (gpCAM `dadeb65`).*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Design custom noise models for experiments with non-uniform or structured noise.
 
@@ -34,50 +34,31 @@ def my_noise(x, hyperparameters):
     
     Returns
     -------
-    noise : np.ndarray
-        Either shape (N,) for diagonal noise (independent per point),
-        or shape (N, N) for full noise covariance matrix.
+    noise : np.ndarray, shape (N,)
+        1-D array of strictly positive noise variances (one per point).
+        (Returning a 2-D (N, N) noise covariance matrix is deprecated via FutureWarning.)
     """
 ```
 
-An optional third argument receives the constructor's `args` dict — gpCAM dispatches
-on parameter count, so declaring it is what turns it on:
-
-```python
-def my_noise(x, hps, args):      # GPOptimizer(..., args={"bounds": ...})
-    ...
-```
-
-Pass it with `noise_function=my_noise`. For analytic gradients, also pass
-`noise_function_grad=`.
+- **Optional `args` parameter:** Declare a 3rd positional parameter `my_noise(x, hps, args)` to receive `GPOptimizer(..., args=...)` (`gp_optimizer.py:179-180`; see `experiment-designer` for the full `args` dispatch table).
+- **Analytic gradient (`noise_function_grad`):** Optional callable `f(x, hyperparameters)` (or `f(x, hyperparameters, args)`) returning a 2-D array of shape `(len(hyperparameters), N)` of partial derivatives of the noise variances (`0.0` for kernel/mean hyperparameter rows; or with `gradient_per_direction=True`, `f(x, hyperparameters, direction)` returning shape `(N,)`, `gp_optimizer.py:188-202`).
+- **Multi-task shapes (`gp_optimizer.py:100-102, 181-182`):**
+  - `noise_variances` array: shape `(V,)` for single-task, `"independent"`, and `"kronecker"` modes; shape `(V, No)` for `multi_task_mode="cartesian"` (with `np.nan` wherever `y_data` is `np.nan`).
+  - `noise_function(x, hps)`: in `"independent"` and `"kronecker"` modes `x` has shape `(N, D)`; in `"cartesian"` mode `x` has shape `(N, D + 1)` with the task index in the last column (`x[:, -1]`), and the function returns a 1-D vector of length `N`.
 
 ## When to Use What
 
 | Scenario | Approach |
 |----------|----------|
 | Known, uniform noise | Use `noise_variances=np.full(N, sigma**2)` — no noise function needed |
-| Known per-point noise | Use `noise_variances=my_array` — no noise function needed |
+| Known per-point noise | Use `noise_variances=my_array` (`(N,)` or `(N, No)` in `"cartesian"` mode) |
 | Unknown uniform noise | Use a noise function with a learnable hyperparameter |
 | Position-dependent noise | Use a noise function that depends on `x` |
 | No noise info at all | Don't provide either — gpCAM defaults to `(0.01 * mean|y|)²` |
 
-## The Index Convention
+## Hyperparameter Offset (`K`)
 
-Every recipe below starts its hyperparameters at index `K`. **Derive `K`; never
-hardcode it.** With the default ARD Matérn kernel the kernel owns `hps[0]` (signal
-variance) plus one length scale per input dimension:
-
-```python
-K = 1 + D    # D = number of input dimensions = x.shape[1]
-```
-
-`K = 3` is only correct for a 2-D input space. In 1-D it points at a length scale;
-in 3-D the noise function returns a squared length scale as a variance while training
-optimizes one parameter that two components are reading. All of that is silent.
-
-The standard ordering across gpCAM skills is **kernel, then mean, then noise** — so
-if you also use a custom prior mean, the noise indices start after the mean's, and
-`K` grows accordingly. See the `prior-mean-functions` skill.
+See `experiment-designer` ("Hyperparameter Coordination") for the canonical layout (`hps[:K]` kernel, `hps[K : K + M]` prior mean, `hps[K + M :]` noise). With the single-task default ARD kernel and no custom prior mean, `K = 1 + D = 1 + x.shape[1]`; with a custom kernel or prior mean, start noise indices at `K + M` (and under `multi_task_mode="cartesian"`, remember `x.shape[1]` is `D + 1`).
 
 ## Recipes
 
@@ -203,33 +184,10 @@ def two_detector_noise(x, hps):
     return noise
 ```
 
-## Hyperparameter Coordination
-
-Noise function hyperparameters come from the **same vector** as kernel and mean
-hyperparameters, and the index ranges must be **disjoint** — fvGP assumes a
-hyperparameter belonging to the noise function has zero kernel derivative, so an
-overlap corrupts the training gradients as well as the model.
+## Choosing Noise Bounds
 
-Standard ordering: **kernel, then mean, then noise.** Noise is last, which is exactly
-why a prior mean function must never read `hps[-1]`.
-
-```python
-# Example layout (2-D input, default kernel, no custom mean):
-# hps[0]     = signal variance     (kernel)
-# hps[1:3]   = length scales       (kernel)
-# hps[3]     = noise std dev       (noise function) ← K = 1 + D = 3
-#
-# Total: 4 hyperparameters
-
-hp_bounds = np.array([
-    [0.001, 100.0],  # signal variance
-    [0.01, 50.0],    # length scale dim 0
-    [0.01, 50.0],    # length scale dim 1
-    [0.001, 10.0],   # noise std dev
-])
-```
+Noise hyperparameters occupy the trailing slice `hps[K + M :]` of the shared hyperparameter vector (see `experiment-designer` for the full layout and disjoint-slice rules).
 
-### Choosing Noise Bounds
 - **Lower bound**: Never 0 — use `0.001` minimum (prevents singular matrices)
 - **Upper bound**: `10 * std(y_data)` — noise shouldn't be larger than the signal
 - **Initial value**: `0.01 * mean(|y_data|)` (gpCAM's own default)
@@ -254,4 +212,4 @@ gpo = GPOptimizer(x_data, y_data, noise_function=learnable_noise)
 4. **Zero noise**: Causes singular matrix errors. Always ensure noise > 0.
 5. **Providing both**: Don't pass `noise_variances` AND `noise_function` — pick one.
 6. **Forgetting bounds**: The noise hyperparameter needs bounds in the `hyperparameter_bounds` array.
-7. **Wrong shape**: Return shape `(N,)` for independent noise, `(N, N)` for correlated.
+7. **Wrong shape**: Return a 1-D vector of shape `(N,)` from `noise_function` (matrix-valued `(N, N)` noise is deprecated via `FutureWarning`).
diff --git a/skills/prior-mean-functions/SKILL.md b/skills/prior-mean-functions/SKILL.md
index 218c60e..3601122 100644
--- a/skills/prior-mean-functions/SKILL.md
+++ b/skills/prior-mean-functions/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: prior-mean-functions
-description: Use when encoding known physics, theoretical models, or expected trends as prior mean functions for gpCAM — useful when there's a baseline expectation the GP should regress against rather than a flat zero prior.
-gpcam_version: "8.4.x"
-fvgp_version: "4.8.x"
-last_verified: "2026-07-23 (gpCAM dadeb65)"
+description: Use when encoding known physics, theoretical models, or expected trends as prior mean functions for gpCAM — useful when there's a baseline expectation or physical model the GP should regress against rather than the default constant data-mean prior (mean(y_data)).
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: gpCAM Prior Mean Functions
 
-*Verified against gpCAM 8.4.x / fvgp 4.8.x — last checked 2026-07-23 (gpCAM `dadeb65`).*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM a21385a).*
 
 Design prior mean functions that encode known physics or expected trends.
 
@@ -54,22 +54,12 @@ gpo = GPOptimizer(
 gpo.train(hyperparameter_bounds=hp_bounds)
 ```
 
-An optional third argument receives the constructor's `args` dict — declare it and
-gpCAM passes it positionally (dispatch is by parameter count, so the name is yours
-to choose):
-
-```python
-def my_mean(x, hps, args):          # GPOptimizer(..., args={"d_spacing": 3.14})
-    return np.full(len(x), args["d_spacing"] * hps[K])
-```
-
-If you also want analytic gradients, pass `prior_mean_function_grad=` returning shape
-`(len(hps), N)`. Without it gpCAM uses finite differences, which is fine but slower.
-
-**If no prior mean is provided**, gpCAM uses **the average of `y_data`** as a constant
-mean (`fvgp.GP._default_mean_function`) — not a flat zero. This matters: away from
-data the posterior reverts to `mean(y_data)`, not to zero. It is a sensible default
-for most cases.
+- **Optional `args` parameter:** Declare a 3rd positional parameter `my_mean(x, hps, args)` to receive `GPOptimizer(..., args=...)` (`gp_optimizer.py:153-154`; see `experiment-designer` for the full `args` dispatch table).
+- **Analytic gradient (`prior_mean_function_grad`):** Optional callable returning shape `(len(hps), N)` (`gp_optimizer.py:165-171`).
+- **Default prior mean:** When `prior_mean_function` is omitted, gpCAM uses `fvgp.GP._default_mean_function`, which is a constant equal to each output's average of `y_data` (`gp_optimizer.py:157-160`) — not zero.
+- **Multi-task shapes (`gp_optimizer.py:153-160`):**
+  - In `"independent"` and `"kronecker"` modes, `x` has shape `(N, D)` and `prior_mean_function` may return `(N,)` or `(N, No)`.
+  - In `multi_task_mode="cartesian"`, `x` has shape `(N, D + 1)` with the task index in the last column (`x[:, -1]`), and `prior_mean_function` returns a 1-D vector of shape `(N,)` (remember `x.shape[1]` is `D + 1`).
 
 ## When NOT to Use a Prior Mean
 
@@ -77,26 +67,9 @@ for most cases.
 - Your model might be wrong → a bad prior mean biases the GP and can hurt more than help
 - You're purely exploring → the GP will learn the mean from data
 
-## The Index Convention
-
-Every recipe below starts its hyperparameters at index `K`. **Derive `K`; never
-hardcode it.** With the default ARD Matérn kernel the kernel owns `hps[0]` (signal
-variance) plus one length scale per input dimension, so:
-
-```python
-K = 1 + D    # D = number of input dimensions = x.shape[1]
-```
-
-`K = 3` is only correct for a 2-D input space. In 1-D it points at a length scale;
-in 3-D it points at the wrong slot and two components silently share a
-hyperparameter. Compute `K` inside the function from `x.shape[1]`, or define it once
-as a module-level constant next to your bounds array so the two cannot drift apart.
-
-**Never use negative indices** (`hps[-1]`). The end of the vector belongs to the
-noise function under the standard layout — see the ordering rule below — so `hps[-1]`
-silently ties your prior mean to the noise amplitude.
+## Hyperparameter Offset (`K`)
 
-With a custom kernel, `K` is however many hyperparameters your kernel reads.
+See `experiment-designer` ("Hyperparameter Coordination") for the canonical layout (`hps[:K]` kernel, `hps[K : K + M]` prior mean, `hps[K + M :]` noise). With the single-task default ARD kernel, `K = 1 + D = 1 + x.shape[1]`; with a custom kernel, `K` is the number of parameters that kernel consumes. Never index with `hps[-1]` in a prior mean function — the trailing slice belongs to the noise function.
 
 ## Recipes
 
@@ -175,51 +148,9 @@ def bragg_mean(x, hps):
     return expected
 ```
 
-## Hyperparameter Coordination
-
-**This is where most bugs happen.** The prior mean function receives the same
-hyperparameter vector as the kernel and noise functions, and fvGP requires the index
-ranges each callable reads to be **disjoint**. The gradient computation depends on
-it: when a hyperparameter index belongs to the mean function, its kernel derivative
-is assumed zero, and vice versa. Overlapping indices therefore produce not just a
-wrong model but wrong gradients.
-
-**The standard ordering used across all gpCAM skills is: kernel, then mean, then
-noise.**
-
-```
-[ kernel hps ][ mean hps ][ noise hps ]
-  0 .. K-1     K .. M-1     M .. end
-```
-
-You must:
-
-1. Decide which indices the mean function uses — derive `K`, don't hardcode it
-2. Add bounds for those hyperparameters
-3. Document the full layout
-
-```python
-# Example layout (2-D input, linear mean, learnable noise):
-# hps[0]     = signal variance     (kernel)
-# hps[1:3]   = length scales       (kernel, 2D input)
-# hps[3]     = intercept           (mean function)   <- K = 1 + D = 3
-# hps[4:6]   = slopes              (mean function)
-# hps[6]     = noise amplitude     (noise function)
-#
-# Total: 7 hyperparameters
-
-hp_bounds = np.array([
-    [0.001, 100.0],   # signal variance
-    [0.01, 50.0],     # length scale dim 0
-    [0.01, 50.0],     # length scale dim 1
-    [-10.0, 10.0],    # intercept
-    [-5.0, 5.0],      # slope dim 0
-    [-5.0, 5.0],      # slope dim 1
-    [0.001, 10.0],    # noise amplitude
-])
-```
+## Complete Worked Example (Kernel + Mean + Noise)
 
-### Complete worked example
+Prior mean hyperparameters sit at `hps[K : K + M]` between the kernel (`hps[:K]`) and noise (`hps[K + M :]`) slices — see `experiment-designer` for the full coordination rules.
 
 ```python
 import numpy as np
diff --git a/skills/transformed-optimizers-advanced/SKILL.md b/skills/transformed-optimizers-advanced/SKILL.md
index f1876bb..77347a0 100644
--- a/skills/transformed-optimizers-advanced/SKILL.md
+++ b/skills/transformed-optimizers-advanced/SKILL.md
@@ -1,16 +1,16 @@
 ---
 name: transformed-optimizers-advanced
-description: Use when observations are constrained — strictly positive (intensities, rates, concentrations) or bounded in [0, 1] (fractions, probabilities). LogGPOptimizer and LogitGPOptimizer fit a GP on transformed observations and push the posterior back through the inverse link, so predictions and credible intervals stay inside the constrained range. Includes how to get raw posterior samples for histograms.
-gpcam_version: "8.4.x"
-fvgp_version: "4.8.x"
-last_verified: "2026-07-23 (gpCAM dadeb65)"
+description: Use when observations are constrained — strictly positive (intensities, rates, concentrations) or bounded in [a, b] (fractions, probabilities, yields). LogGPOptimizer and LogitGPOptimizer fit a GP on transformed observations and push the posterior back through the inverse link, so predictions and credible intervals stay inside the constrained range. Includes how to get raw posterior samples for histograms.
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: Transformed-Output GP Optimizers
 
-*Verified against gpCAM 8.4.x / fvgp 4.8.x — last checked 2026-07-23 (gpCAM `dadeb65`).*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM `a21385a`).*
 
-Use this skill when measurements are guaranteed positive (`y > 0`) or bounded (`y ∈ [0, 1]`). A plain GP doesn't know about these constraints and will happily predict negative intensities or probabilities outside `[0, 1]`. The transformed optimizers fit the GP in an unconstrained "link" space (log or logit) and push the Gaussian posterior back through the inverse link via `evaluate_posterior(x)`, so the original-scale predictions and credible intervals are guaranteed to respect the constraint.
+Use this skill when measurements are guaranteed positive (`y > 0`) or bounded (`y ∈ [a, b]`, default `[0, 1]`). A plain GP doesn't know about these constraints and will happily predict negative intensities or values outside `[a, b]`. The transformed optimizers fit the GP in an unconstrained "link" space (log or logit) and push the Gaussian posterior back through the inverse link via `evaluate_posterior(x)`, so the original-scale predictions and credible intervals are guaranteed to respect the constraint.
 
 ## When to use which
 
@@ -99,7 +99,7 @@ Distributions are: Gaussian for plain `GPOptimizer`, lognormal for `LogGPOptimiz
 
 Two cases that need care:
 
-- **`"target probability"`**: pass already-transformed bounds. For `LogGPOptimizer`, use `args={"a": np.log(a), "b": np.log(b)}`. For `LogitGPOptimizer`, use `args={"a": scipy.special.logit(a), "b": scipy.special.logit(b)}`.
+- **`"target probability"`**: pass already-transformed bounds (via `GPOptimizer(..., args=...)` or `gpo.ask(..., args=...)`). For `LogGPOptimizer`, use `args={"a": np.log(a), "b": np.log(b)}`. For `LogitGPOptimizer` with `range=(lower, upper)`, first rescale `a_norm = (a - lower) / (upper - lower)` and `b_norm = (b - lower) / (upper - lower)`, then pass `args={"a": scipy.special.logit(a_norm), "b": scipy.special.logit(b_norm)}`.
 - **`"expected improvement"` / `"probability of improvement"`**: these compare against the best observed transformed value, which is also the best original value, so semantics carry over. The score *magnitude* differs from a plain GP fit but the *ranking* is consistent.
 
 ### Writing custom acquisition functions against a transformed optimizer
@@ -158,8 +158,8 @@ When you provide `noise_variances` (observation noise) at construction or via `t
 
 ## Pickling
 
-`LogGPOptimizer` and `LogitGPOptimizer` pickle the same way as `GPOptimizer`. `LogitGPOptimizer`'s extra attributes (`eps`, `n_samples`) survive the round-trip.
+`LogGPOptimizer` and `LogitGPOptimizer` pickle the same way as `GPOptimizer`. `LogitGPOptimizer`'s extra attributes (`eps`, `n_samples`, `range`) survive the round-trip.
 
 ## Example notebook
 
-See [`examples/GPOptimizer_LogAndLogit.ipynb`](../../examples/GPOptimizer_LogAndLogit.ipynb) for a side-by-side visual tour — latent-space GP, back-transformed band, plain-GP comparison that violates the constraint, and posterior-sample histograms (lognormal vs. logistic-normal).
+See [`GPOptimizer_LogAndLogit.ipynb`](../../examples/GPOptimizer_LogAndLogit.ipynb) for a side-by-side visual tour — latent-space GP, back-transformed band, plain-GP comparison that violates the constraint, and posterior-sample histograms (lognormal vs. logistic-normal).
diff --git a/skills/troubleshooting/SKILL.md b/skills/troubleshooting/SKILL.md
index 68bf3ef..27c0281 100644
--- a/skills/troubleshooting/SKILL.md
+++ b/skills/troubleshooting/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: troubleshooting
 description: Use when a gpCAM script crashes or misbehaves — Cholesky/non-positive-definite errors, singular covariance, hyperparameter bounds mismatches, NaNs, kernel signature errors, Dask hangs, pickling failures, or an ask/tell loop that keeps picking the same point. Maps the actual error text to the cause and the fix.
-gpcam_version: "8.6.x"
-fvgp_version: "4.9.1"
-last_verified: "2026-10-14 (fvgp 4.9.1); multi-task note updated for gpCAM 8.7.0 on 2026-10-28"
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: Troubleshooting gpCAM
 
-*Verified against gpCAM 8.6.x / fvgp 4.9.1 — last checked 2026-10-14.*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM `a21385a`).*
 
 A decision tree from the error you actually saw to the thing that is actually wrong.
 The audience is a scientist who does not want to learn GP internals to get unstuck.
@@ -110,17 +110,7 @@ assert len(gpo.hyperparameters) == len(hp_bounds), (
     f"{len(hp_bounds)} bounds rows but {len(gpo.hyperparameters)} hyperparameters")
 ```
 
-Count what each callable reads and confirm the total:
-
-| Component | Hyperparameters (default ARD kernel, D input dims) |
-|---|---|
-| Kernel | `1 + D` — `hps[0]` signal variance, `hps[1:D+1]` length scales |
-| Prior mean | however many it reads, starting at `K = 1 + D` |
-| Noise | however many it reads, starting after the mean's |
-
-**The usual root cause is a hardcoded `K = 3`**, which is correct only for `D = 2`.
-Derive it: `K = 1 + x.shape[1]`. See the `prior-mean-functions` and `noise-functions`
-skills.
+**The usual root cause is a hardcoded `K = 3`**, which is correct only for `D = 2` with the default single-task ARD kernel. Derive `K` from the input dimensionality (`K = 1 + x.shape[1]` for single-task, or `K = 1 + (D + 1)` in `"cartesian"` multi-task mode) and allocate disjoint slices in order: kernel (`0..K-1`), prior mean (`K..K+M-1`), noise (`K+M..`). See the `experiment-designer` skill ("Hyperparameter Coordination") for the canonical index table.
 
 Related and nastier: if the indices *overlap* rather than run past the end, there is
 **no error at all** — two components silently share a hyperparameter. See
@@ -163,22 +153,7 @@ Exception: No valid kernel function signature
 ```
 (or the prior-mean / noise equivalent)
 
-gpCAM dispatches on **parameter count**, via `inspect.signature`:
-
-| Callable | Without `args` | With `args` |
-|---|---|---|
-| Kernel | `f(x1, x2, hps)` | `f(x1, x2, hps, args)` |
-| Prior mean | `f(x, hps)` | `f(x, hps, args)` |
-| Noise | `f(x, hps)` | `f(x, hps, args)` |
-| Cost | `f(origin, x)` | **not supported** |
-| Acquisition | `f(x, gp_obj)` | **not supported** |
-
-Any other count raises. Note that a default value still counts as a parameter, so
-`def my_kernel(x1, x2, hps, args=None)` is a **4**-parameter kernel and will be handed
-the `args` dict.
-
-Cost and acquisition functions never receive `args` — bind their parameters with
-`functools.partial` or a closure.
+gpCAM dispatches on **parameter count** via `inspect.signature`: kernels take 3 (`f(x1, x2, hps)`) or 4 (`f(x1, x2, hps, args)`) parameters, mean and noise functions take 2 (`f(x, hps)`) or 3 (`f(x, hps, args)`), and cost (`f(origin, x)`) and custom acquisition (`f(x, gp_obj)`) functions always take 2 parameters and never receive `args` (bind extra parameters with `functools.partial`). A default value still counts toward `len(inspect.signature(f).parameters)`, so `def my_kernel(x1, x2, hps, args=None)` is dispatched as a 4-parameter kernel. See the `experiment-designer` skill ("User-function arguments (`args`)") for the canonical dispatch table.
 
 ---
 
@@ -236,7 +211,7 @@ With a transformed optimizer, pass the bounds **already transformed**
 Other acquisition strings that need attention:
 - Multi-task (`x_out is not None`) supports a **smaller set** than single-task.
   `"maximum"`, `"minimum"`, `"gradient"`, `"probability of improvement"`, and
-  `"target probability"` are single-task only.
+  `"target probability"` are single-task only, and `"total correlation"` is not available when `multi_task_mode="kronecker"`.
 - A custom callable must return a **1-D array of length `len(x)`**. Returning a scalar
   or a `(V, 1)` column raises further down in the optimizer.
 
@@ -292,19 +267,7 @@ AttributeError: Can't pickle local object 'make_cost.<locals>.cost'
 
 `GPOptimizer` instances are picklable, but **the callables you attached must pickle by
 reference**, which means they must be defined at **module level**. A closure or a
-lambda defined inside another function cannot pickle.
-
-```python
-# Breaks pickling:
-def make_cost(speed):
-    def cost(origin, x): return 1.0 + speed * np.linalg.norm(x - origin, axis=1)
-    return cost
-
-# Pickles fine:
-from functools import partial
-def cost(origin, x, speed): return 1.0 + speed * np.linalg.norm(x - origin, axis=1)
-gpo.cost_function = partial(cost, speed=2.5)
-```
+lambda defined inside another function cannot pickle — replace it with a module-level function bound via `functools.partial(cost, speed=2.5)` (see the `cost-functions` skill for a complete example).
 
 When reloading saved arrays that contain Python objects, remember
 `np.load(path, allow_pickle=True)`.
diff --git a/skills/uncertainty-calibration/SKILL.md b/skills/uncertainty-calibration/SKILL.md
index 44c9c99..a9527e4 100644
--- a/skills/uncertainty-calibration/SKILL.md
+++ b/skills/uncertainty-calibration/SKILL.md
@@ -1,14 +1,14 @@
 ---
 name: uncertainty-calibration
 description: Use to validate that a gpCAM posterior's error bars are trustworthy before publishing them or letting an autonomous loop act on them — scoring rules (NLPD, CRPS, MSLL), calibration/coverage curves, PICP and interval width, and what to do when the model is over- or under-confident.
-gpcam_version: "8.4.x"
-fvgp_version: "4.8.x"
-last_verified: "2026-07-23 (gpCAM dadeb65)"
+gpcam_version: "8.7.x"
+fvgp_version: "4.11.x"
+last_verified: "2026-10-28 (gpCAM 8.7.0, a21385a)"
 ---
 
 # Skill: Validating gpCAM Uncertainty
 
-*Verified against gpCAM 8.4.x / fvgp 4.8.x — last checked 2026-07-23 (gpCAM `dadeb65`).*
+*Verified against gpCAM 8.7.0 / fvgp 4.11.1 — last checked 2026-10-28 (gpCAM `a21385a`).*
 
 A GP's headline product is not the prediction — it is the *uncertainty on* the
 prediction. An autonomous loop acts on that uncertainty at every step: `"variance"`
@@ -46,9 +46,7 @@ If `PICP ≈ 0.95` and `NLPD` is small, the posterior is roughly trustworthy. If
 `PICP` is far from `0.95`, stop and read the rest of this skill — nothing downstream
 is reliable.
 
-**Shapes must match exactly.** `rmse`, `nlpd`, and `crps` assert
-`posterior_mean(x_test).shape == y_test.shape`. For a single-task GP, `y_test` must
-be shape `(V,)`, not `(V, 1)`.
+**Shapes must match `posterior_mean(x_test)["m(x)"].shape` (`gp_optimizer.py:615`).** `rmse`, `nlpd`, and `crps` assert `posterior_mean(x_test)["m(x)"].shape == y_test.shape`: for a single-task GP, `y_test` must be 1-D `(V,)` (not `(V, 1)`); for a multi-task GP (including `multi_task_mode="cartesian"`), `x_test` is `(V, D)` and `y_test` must be `(V, No)` across all `No = len(x_out)` tasks.
 
 ## The Metrics
 
@@ -198,7 +196,7 @@ history = []
 for i in range(N_ITERATIONS):
     ...                          # ask / measure / tell
     if (i + 1) % RETRAIN_EVERY == 0:
-        gpo.train(hyperparameter_bounds=hp_bounds, method="local")
+        gpo.train(hyperparameter_bounds=hp_bounds, method="local", n_restarts=1)
         history.append({
             "n": len(gpo.y_data),
             "rmse": gpo.rmse(holdout_x, holdout_y),
@@ -210,7 +208,7 @@ for i in range(N_ITERATIONS):
 
 Expect `rmse` and `nlpd` to fall and `picp95` to settle near 0.95. If `picp95` drifts
 steadily *down* as points accumulate, the hyperparameters are overfitting the growing
-dataset — retrain globally rather than locally, or tighten the length-scale lower
+dataset — re-run multi-restart local training (`n_restarts=5`, the `train()` default) or `method="mcmc"`, or tighten the length-scale lower
 bound.
 
 Never score against points the GP has been told; `rmse` on training data measures
@@ -229,7 +227,7 @@ Points should fall on the diagonal with scatter consistent with the error bars.
 2. **Unpacking `crps` as a scalar.** It returns `(mean, std)`.
 3. **`picp` without `mpiw`.** Coverage is trivially achievable with absurdly wide intervals; always report both, or use `interval_score`, which combines them.
 4. **Testing on adaptively-acquired points.** They are not i.i.d.; hold out a random subset before the loop starts.
-5. **Shape mismatch.** `y_test` must be `(V,)` for a single-task GP — the assertion message tells you the expected shape.
+5. **Shape mismatch.** `y_test` must match `posterior_mean(x_test)["m(x)"].shape` — `(V,)` for a single-task GP (not `(V, 1)`) and `(V, No)` for a multi-task GP.
 6. **Judging a transformed optimizer by latent-space metrics.** Report original-scale coverage via `evaluate_posterior`.
 7. **Calibrating once and never again.** Calibration drifts as the loop adds data and retrains.