tranformed GPs, updated AI skill, new light forced docs, new tests

MarcusMNoack · MarcusMNoack · commit 90b598571aee · 2026-06-04T16:40:36.000-07:00
diff --git a/CLAUDE.md b/CLAUDE.md
@@ -20,6 +20,7 @@ You are helping scientists design autonomous experiments using gpCAM, a Gaussian
 - **Cost functions (travel time)**: `skills/cost-functions/SKILL.md`
 - **Large-scale (>10k points)**: `skills/gp2scale-advanced/SKILL.md`
 - **Multi-task/multi-output**: `skills/multi-task-advanced/SKILL.md`
+- **Constrained observations (y>0 or y∈[0,1])**: `skills/transformed-optimizers-advanced/SKILL.md`
 
 These skills also ship as a Claude Code plugin marketplace (see README.md), so they are available outside this repo once installed.
 
diff --git a/README.md b/README.md
@@ -69,6 +69,7 @@ To update later, run `/plugin marketplace update gpcam`; to remove, `/plugin uni
 | **cost-functions** | Account for motor travel time, settling, directional costs, and zone-based penalties. |
 | **gp2scale-advanced** | Large-scale experiments (>10k points) using sparse kernels and Dask distributed computing. |
 | **multi-task-advanced** | Multi-output / function-valued experiments with `fvGPOptimizer`. |
+| **transformed-optimizers-advanced** | Constrained observations: `LogGPOptimizer` for `y > 0` and `LogitGPOptimizer` for `y ∈ [0, 1]` — predictions and credible intervals stay inside the constrained range. |
 
 Once installed, the skills activate automatically when you describe an experiment design problem to Claude, or you can invoke one explicitly (e.g. _"use the experiment-designer skill to set up an adaptive XRD scan"_).
 
diff --git a/docs/source/_static/custom.css b/docs/source/_static/custom.css
@@ -1,3 +1,9 @@
+/* Keep native browser controls (scrollbars, form widgets) in their light
+   variants -- pydata's default_mode="light" handles the rest. */
+:root {
+    color-scheme: light;
+}
+
 /* ------------------------------------------------------------------ */
 /* Headshot grid used on contributor/team pages                        */
 /* ------------------------------------------------------------------ */
diff --git a/docs/source/api/overview.md b/docs/source/api/overview.md
@@ -6,6 +6,7 @@
 
 gpOptimizer.md
 fvgpOptimizer.md
+transformedOptimizers.md
 gpMCMC.md
 kernels.md
 logging.md
diff --git a/docs/source/api/transformedOptimizers.md b/docs/source/api/transformedOptimizers.md
@@ -0,0 +1,34 @@
+# Transformed Optimizers
+
+Single-task optimizers that fit a Gaussian process to observations after a fixed
+output transform (link function), and use
+:py:meth:`gpcam.GPOptimizer.evaluate_posterior` to push the latent Gaussian
+posterior back through the inverse link to the original observation space.
+
+Use :py:class:`gpcam.LogGPOptimizer` for strictly positive observations in
+``(0, inf)`` — predictions and credible intervals are guaranteed positive, and
+the original-scale mean/std are available in closed form (lognormal).
+
+Use :py:class:`gpcam.LogitGPOptimizer` for observations bounded in ``[0, 1]`` —
+predictions and credible intervals are guaranteed inside ``(0, 1)``; the
+original-scale mean/std are estimated by Monte Carlo (logistic-normal has no
+closed form). Observations are clipped to ``[eps, 1 - eps]`` because
+``logit(0)``/``logit(1)`` are infinite.
+
+The inherited :py:meth:`posterior_mean` / :py:meth:`posterior_covariance`
+operate in the GP's modeling space (log / logit); use
+:py:meth:`evaluate_posterior` to obtain a posterior on the original scale.
+
+## LogGPOptimizer
+
+```{eval-rst}
+.. autoclass:: gpcam.gp_optimizer.LogGPOptimizer
+    :members:
+```
+
+## LogitGPOptimizer
+
+```{eval-rst}
+.. autoclass:: gpcam.gp_optimizer.LogitGPOptimizer
+    :members:
+```
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -71,6 +71,11 @@ def setup(app):
 # -- Options for HTML output -------------------------------------------------
 html_theme = 'pydata_sphinx_theme'
 
+# Force light mode regardless of the visitor's OS preference. pydata-sphinx-theme
+# sets data-theme="light" on <html> when default_mode == "light", and the navbar
+# theme-switcher is excluded via navbar_end so visitors can't toggle to dark.
+html_context = {'default_mode': 'light'}
+
 html_theme_options = {
     'logo': {
         'image_light': '_static/gpCAM_bright_bg.png',
diff --git a/docs/source/examples/index.md b/docs/source/examples/index.md
@@ -11,6 +11,7 @@ GPOptimizer_NonEuclideanInputSpaces
 GPOptimizer_gp2ScaleTest
 GPOptimizer_MultiTaskTest
 GPOptimizer_Optimization
+GPOptimizer_LogAndLogit
 ```
 
 The notebooks below walk through common gpCAM use cases.
@@ -23,3 +24,4 @@ Each can be downloaded and run locally after `pip install gpcam`.
 - **[gp2Scale](GPOptimizer_gp2ScaleTest.ipynb)** — large-scale sparse GP with Dask distributed computation
 - **[Multi-task GP](GPOptimizer_MultiTaskTest.ipynb)** — multi-output regression with the `fvGPOptimizer` class
 - **[Function optimization](GPOptimizer_Optimization.ipynb)** — minimizing a known function with the `optimize` loop
+- **[Log and Logit GP optimizers](GPOptimizer_LogAndLogit.ipynb)** — fitting strictly-positive data with `LogGPOptimizer` and bounded `[0, 1]` data with `LogitGPOptimizer`
diff --git a/examples/GPOptimizer_LogAndLogit.ipynb b/examples/GPOptimizer_LogAndLogit.ipynb
diff --git a/gpcam/__init__.py b/gpcam/__init__.py
@@ -7,10 +7,13 @@
 
 from .gp_optimizer import GPOptimizer
 from .gp_optimizer import fvGPOptimizer
+from .gp_optimizer import LogGPOptimizer
+from .gp_optimizer import LogitGPOptimizer
 from .gp_mcmc import gpMCMC, ProposalDistribution
 from .autonomous_experimenter import AutonomousExperimenterGP
 from .autonomous_experimenter import AutonomousExperimenterFvGP
 
-__all__ = ['GPOptimizer', 'fvGPOptimizer','AutonomousExperimenterGP','AutonomousExperimenterFvGP', 'gpMCMC', 'ProposalDistribution']
+__all__ = ['GPOptimizer', 'fvGPOptimizer', 'LogGPOptimizer', 'LogitGPOptimizer',
+           'AutonomousExperimenterGP', 'AutonomousExperimenterFvGP', 'gpMCMC', 'ProposalDistribution']
 
 logger.disable('gpcam')
diff --git a/gpcam/gp_optimizer.py b/gpcam/gp_optimizer.py
@@ -1,5 +1,7 @@
 #!/usr/bin/env python
+import warnings
 import numpy as np
+from scipy.special import expit, logit
 from fvgp import fvGP
 from .gp_optimizer_base import GPOptimizerBase
 
@@ -687,3 +689,111 @@ def __init__(
                          logging=logging,
                          multi_task=True,
                          args=args)
+
+
+class LogGPOptimizer(GPOptimizer):
+    """
+    A single-task :py:class:`GPOptimizer` for strictly positive observations in (0, inf).
+
+    Observations are modeled in log-space (the GP sees ``log(y)``), and posterior
+    predictions are mapped back to the original scale with ``exp`` via
+    :py:meth:`evaluate_posterior`, which guarantees strictly positive predictions and
+    credible intervals. ``exp`` of a Gaussian is lognormal, so the original-scale mean
+    and standard deviation are available in closed form.
+
+    All constructor arguments are identical to :py:class:`GPOptimizer`. Note that the
+    inherited :py:meth:`posterior_mean` / :py:meth:`posterior_covariance` operate in
+    log-space; use :py:meth:`evaluate_posterior` for the original (positive) scale.
+
+    Acquisition functions: :py:meth:`ask` optimizes the GP in log-space. Because ``log``
+    is monotone increasing, ranking acquisitions (``variance``, ``ucb``, ``lcb``,
+    ``maximum``, ``minimum``) still identify the same locations as on the original scale.
+    For ``target probability``, pass bounds already in log-space (``np.log(a)``, ``np.log(b)``).
+    """
+
+    def _prepare(self, y):
+        y = np.asarray(y, dtype=float)
+        if np.any(y <= 0.0):
+            raise ValueError("LogGPOptimizer requires strictly positive observations (y > 0).")
+        return y
+
+    def _forward(self, y):
+        return np.log(y)
+
+    def _inverse(self, z):
+        return np.exp(z)
+
+    def _forward_deriv(self, y):
+        return 1.0 / y
+
+    def _moments(self, mu, var):
+        # exp(Normal(mu, var)) is lognormal
+        mean = np.exp(mu + var / 2.0)
+        std = np.sqrt((np.exp(var) - 1.0) * np.exp(2.0 * mu + var))
+        return mean, std
+
+
+class LogitGPOptimizer(GPOptimizer):
+    """
+    A single-task :py:class:`GPOptimizer` for observations bounded in [0, 1].
+
+    Observations are modeled in logit (log-odds) space (the GP sees ``logit(y)``), and
+    posterior predictions are mapped back with the logistic/sigmoid via
+    :py:meth:`evaluate_posterior`, which guarantees predictions and credible intervals
+    inside (0, 1). Because ``logit(0)`` / ``logit(1)`` are infinite, observations are
+    clipped to ``[eps, 1 - eps]`` (a warning is emitted when clipping occurs). The
+    logistic-normal distribution has no closed-form moments, so the original-scale mean
+    and standard deviation are estimated by Monte-Carlo.
+
+    Parameters
+    ----------
+    eps : float, optional
+        Clipping margin for the open interval; observations are clipped to
+        ``[eps, 1 - eps]`` before the logit transform. The default is 1e-6.
+    n_samples : int, optional
+        Number of Monte-Carlo samples used to estimate the original-scale mean/std in
+        :py:meth:`evaluate_posterior`. The default is 10000.
+
+    Notes
+    -----
+    All other constructor arguments are identical to :py:class:`GPOptimizer`. The
+    inherited :py:meth:`posterior_mean` / :py:meth:`posterior_covariance` operate in
+    logit-space; use :py:meth:`evaluate_posterior` for the original (0, 1) scale. The
+    acquisition-function note for :py:class:`LogGPOptimizer` applies here too (pass
+    ``target probability`` bounds in logit-space).
+    """
+
+    def __init__(self, x_data=None, y_data=None, eps=1e-6, n_samples=10000, **kwargs):
+        self.eps = eps
+        self.n_samples = n_samples
+        super().__init__(x_data=x_data, y_data=y_data, **kwargs)
+
+    def _prepare(self, y):
+        y = np.asarray(y, dtype=float)
+        if np.any(y < self.eps) or np.any(y > 1.0 - self.eps):
+            warnings.warn("LogitGPOptimizer clipped observations to "
+                          f"[{self.eps}, {1.0 - self.eps}] before the logit transform.")
+        return np.clip(y, self.eps, 1.0 - self.eps)
+
+    def _forward(self, y):
+        return logit(y)
+
+    def _inverse(self, z):
+        return expit(z)
+
+    def _forward_deriv(self, y):
+        return 1.0 / (y * (1.0 - y))
+
+    def _moments(self, mu, var):
+        # sigmoid(Normal(mu, var)) is logistic-normal -> no closed form, estimate by MC
+        mu = np.asarray(mu).reshape(-1)
+        sd = np.sqrt(np.asarray(var).reshape(-1))
+        samples = expit(np.random.normal(loc=mu[:, None], scale=sd[:, None],
+                                         size=(mu.shape[0], self.n_samples)))
+        return samples.mean(axis=1), samples.std(axis=1)
+
+    def __getstate__(self):
+        state = super().__getstate__()
+        state["eps"] = self.eps
+        state["n_samples"] = self.n_samples
+        return state
diff --git a/gpcam/gp_optimizer_base.py b/gpcam/gp_optimizer_base.py
@@ -7,6 +7,7 @@
 import warnings
 import random
 from distributed import Client
+from scipy.stats import norm
 
 
 # TODO (for gpCAM):
@@ -106,6 +107,9 @@ def _initializeGP(self, x_data, y_data, noise_variances=None):
         If data is prided at initialization this function is NOT needed.
         It has the same parameters as the initialization of the class.
         """
+        y_prepared = self._prepare(np.asarray(y_data))
+        noise_variances = self._transform_noise_variances(y_prepared, noise_variances)
+        y_data = self._forward(y_prepared)
         if self.multi_task: self.x_out = np.arange(y_data.shape[1])
         else: self.x_out = None
 
@@ -145,6 +149,7 @@ def get_data(self):
                 "input dim": self.input_space_dimension,
                 "x data": self.x_data,
                 "y data": self.y_data,
+                "original y data": self._inverse(self.y_data),
                 "measurement variances": self.likelihood.V,
                 "hyperparameters": self.hyperparameters,
                 "cost function": self.cost_function}
@@ -161,6 +166,99 @@ def get_data(self):
 
         ############################################################################
 
+    # ---- Output-transform hooks (identity by default; overridden by transformed optimizers) ----
+    def _prepare(self, y):
+        """Validate/condition original-space observations before the forward transform."""
+        return y
+
+    def _forward(self, y):
+        """Map original-space observations into the GP's modeling space."""
+        return y
+
+    def _inverse(self, z):
+        """Map GP-space values back to the original observation space."""
+        return z
+
+    def _forward_deriv(self, y):
+        """Derivative of the forward transform w.r.t. y (for delta-method noise propagation)."""
+        return 1.0
+
+    def _moments(self, mu, var):
+        """Mean and std in the original space given GP-space mean and variance."""
+        return mu, np.sqrt(var)
+
+    def _samples(self, mu, var, n_samples):
+        """
+        Draw `n_samples` posterior samples in the original space at each evaluation point.
+
+        Default: sample the GP-space Gaussian and push through :py:meth:`_inverse`.
+        Identity transform -> Gaussian samples; Log -> lognormal; Logit -> logistic-normal.
+        Subclasses with an analytic sampler may override.
+        Returns an array of shape ``(n_points, n_samples)``.
+        """
+        mu = np.asarray(mu).reshape(-1)
+        sd = np.sqrt(np.asarray(var).reshape(-1))
+        z = np.random.normal(loc=mu[:, None], scale=sd[:, None],
+                             size=(mu.shape[0], int(n_samples)))
+        return self._inverse(z)
+
+    def _transform_noise_variances(self, y_prepared, noise_variances):
+        """Delta-method transform of observation variances: var_z ~= (g'(y))**2 * var_y."""
+        if noise_variances is None: return None
+        d = self._forward_deriv(y_prepared)
+        return np.asarray(noise_variances) * np.asarray(d) ** 2
+
+    def evaluate_posterior(self, x, x_out=None, level=0.95, return_samples=False, n_samples=10000):
+        """
+        Posterior of the original-space observations at points `x`.
+
+        Unlike :py:meth:`posterior_mean` / :py:meth:`posterior_covariance`, which act in
+        the GP's modeling space, this pushes the Gaussian posterior through the inverse
+        output transform. For the default (identity) transform it simply bundles the
+        Gaussian posterior. For monotone transforms the median and credible bounds map
+        exactly; the mean/std are exact when a closed form exists (e.g. lognormal) and
+        otherwise estimated (see the specific optimizer).
+
+        Parameters
+        ----------
+        x : np.ndarray | list
+            Points at which to evaluate the posterior. np.ndarray of shape (V x D) or list.
+        x_out : np.ndarray, optional
+            Output-space positions; only used in the multi-task setting.
+        level : float, optional
+            Central credible-interval mass in (0, 1). The default is 0.95.
+        return_samples : bool, optional
+            If True, also draw `n_samples` posterior samples in the original space at each
+            point in `x` and include them under the key ``"samples"``. The default is False.
+        n_samples : int, optional
+            Number of samples to draw when ``return_samples=True``. The default is 10000.
+
+        Return
+        ------
+        Posterior summary : dict
+            Keys `"median"`, `"mean"`, `"std"`, `"lower"`, `"upper"`, `"level"`. When
+            ``return_samples=True``, an additional `"samples"` array of shape
+            ``(n_points, n_samples)`` is included, where ``samples[i]`` are the draws at
+            the i-th input point.
+        """
+        assert self.gp, "GP not yet initialized; tell() data!"
+        if x_out is None: x_out = self.x_out
+        mu = np.atleast_1d(np.asarray(self.posterior_mean(x, x_out=x_out)["m(x)"]))
+        var = np.atleast_1d(np.asarray(self.posterior_covariance(x, x_out=x_out, variance_only=True)["v(x)"]))
+        sd = np.sqrt(var)
+        z = norm.ppf(0.5 + level / 2.0)
+        mean, std = self._moments(mu, var)
+        result = {
+            "median": self._inverse(mu),
+            "mean": mean,
+            "std": std,
+            "lower": self._inverse(mu - z * sd),
+            "upper": self._inverse(mu + z * sd),
+            "level": level}
+        if return_samples:
+            result["samples"] = self._samples(mu, var, n_samples)
+        return result
+
     def evaluate_acquisition_function(self, x, x_out=None, acquisition_function="variance", origin=None, args=None):
         """
         Function to evaluate the acquisition function.
@@ -239,7 +337,9 @@ def tell(self, x, y, noise_variances=None, append=True, rank_n_update=None):
 
         if rank_n_update is None: rank_n_update = append
         if self.gp:
-            self.update_gp_data(x, y, noise_variances_new=noise_variances,
+            y_prepared = self._prepare(np.asarray(y))
+            noise_variances = self._transform_noise_variances(y_prepared, noise_variances)
+            self.update_gp_data(x, self._forward(y_prepared), noise_variances_new=noise_variances,
                                 append=append, rank_n_update=rank_n_update)
         else:
             self._initializeGP(x, y, noise_variances=noise_variances)
diff --git a/skills/experiment-designer/SKILL.md b/skills/experiment-designer/SKILL.md
diff --git a/skills/gp2scale-advanced/SKILL.md b/skills/gp2scale-advanced/SKILL.md
diff --git a/skills/transformed-optimizers-advanced/SKILL.md b/skills/transformed-optimizers-advanced/SKILL.md
diff --git a/tests/test_gpCAM.py b/tests/test_gpCAM.py
