0.33.13

Author: bartzbeielstein

notebooks/00_spotPython_tests.ipynb

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "spotpython"
-version = "0.33.12"
+version = "0.33.13"
 authors = [
   { name="T. Bartz-Beielstein", email="tbb@bartzundbartz.de" }
 ]

pyproject.toml

@@ -451,7 +451,20 @@ def predict(self, X: np.ndarray, return_std=False, return_val: str = "y") -> np.
             >>> print("Predictions:", y_pred)
         """
         self.return_std = return_std
-        X = np.atleast_2d(X)
+        X = np.asarray(X)
+
+        # Ensure X has shape (n_samples, n_features=self.k)
+        if X.ndim == 1:
+            X = X.reshape(-1, self.k)
+        else:
+            if X.shape[1] != self.k:
+                if X.shape[0] == self.k:  # common case: row/col swap for 1D
+                    X = X.T
+                elif self.k == 1:
+                    X = X.reshape(-1, 1)
+                else:
+                    raise ValueError(f"X has shape {X.shape}, expected (*, {self.k}).")
+
         if return_std:
             # Return predictions and standard deviations
             # Compatibility with scikit-learn
@@ -549,8 +562,7 @@ def likelihood(self, x: np.ndarray) -> Tuple[float, np.ndarray, np.ndarray]:
 
         Args:
             x (np.ndarray):
-                1D array of log(theta) parameters of length k. If self.method is "regression" or
-                "reinterpolation", length is k+1 and the last element of x is the log(noise) parameter.
+                1D array of log(theta), log(Lambda) (if method is "regression" or "reinterpolation"), and p values (if optim_p is True).
 
         Returns:
             (float, np.ndarray, np.ndarray):
@@ -762,11 +774,24 @@ def max_likelihood(self, bounds: List[Tuple[float, float]]) -> Tuple[np.ndarray,
 
         Returns:
             (np.ndarray, float): (best_x, best_fun) where best_x is the
-            optimal log(theta) array and best_fun is the minimized negative log-likelihood.
+            optimal [log(theta) log(lambda) p] array and best_fun is the minimized negative log-likelihood.
+
+        Examples:
+            >>> import numpy as np
+                from spotpython.surrogate.kriging import Kriging
+                # Training data
+                X_train = np.array([[0.0, 0.0], [0.5, 0.5], [1.0, 1.0]])
+                y_train = np.array([0.1, 0.2, 0.3])
+                # Fit the Kriging model
+                model = Kriging().fit(X_train, y_train)
+                bounds = [(-3.0, 2.0), (-3.0, 2.0), (-9.0, 2.0)] # Example bounds for log(theta) and log(lambda)
+                best_x, best_fun = model.max_likelihood(bounds)
+                print("Optimal parameters (log(theta),log(theta), log(lambda)):", best_x)
+                print("Minimized negative log-likelihood:", best_fun)
         """
 
-        def objective(logtheta_lambda):
-            neg_ln_like, _, _ = self.likelihood(logtheta_lambda)
+        def objective(logtheta_loglambda_p: np.ndarray) -> float:
+            neg_ln_like, _, _ = self.likelihood(logtheta_loglambda_p)
             return neg_ln_like
 
         result = differential_evolution(objective, bounds)
@@ -793,17 +818,17 @@ def plot(self, i: int = 0, j: int = 1, show: Optional[bool] = True, add_points:
         Returns:
             None
 
-        Note:
-            * This method provides only a basic plot. For more advanced plots,
-                use the `plot_contour()` method of the `Spot` class.
+        Notes:
+            * This method is a wrapper around the `plotkd` function for 2D plots.
+            * For 1D plots, it generates a line plot of the surrogate model.
 
         Examples:
             >>> import numpy as np
                 from spotpython.fun.objectivefunctions import Analytical
                 from spotpython.spot import spot
                 from spotpython.utils.init import fun_control_init, design_control_init
                 # 1-dimensional example
-                fun = analytical().fun_sphere
+                fun = Analytical().fun_sphere
                 fun_control=fun_control_init(lower = np.array([-1]),
                                             upper = np.array([1]),
                                             noise=False)
@@ -816,7 +841,7 @@ def plot(self, i: int = 0, j: int = 1, show: Optional[bool] = True, add_points:
                 S.fit_surrogate()
                 S.surrogate.plot()
                 # 2-dimensional example
-                fun = analytical().fun_sphere
+                fun = Analytical().fun_sphere
                 fun_control=fun_control_init(lower = np.array([-1, -1]),
                                             upper = np.array([1, 1]),
                                             noise=False)
@@ -831,10 +856,15 @@ def plot(self, i: int = 0, j: int = 1, show: Optional[bool] = True, add_points:
         """
         if self.k == 1:
             n_grid = 100
-            x = linspace(self.min_X[0], self.max_X[0], num=n_grid)
+            x = linspace(self.min_X[0], self.max_X[0], num=n_grid).reshape(-1, 1)
             y = self.predict(x)
             plt.figure()
-            plt.plot(x, y, "k")
+            plt.plot(x.ravel(), y, "k")
+            # add the original points
+            if add_points:
+                plt.plot(self.X_[:, 0], self.y_, "ro")
+            # add grid
+            plt.grid()
             if show:
                 plt.show()
         else:

src/spotpython/surrogate/kriging.py

@@ -0,0 +1,236 @@
+import numpy as np
+import pytest
+
+from spotpython.surrogate.kriging import Kriging
+
+
+def _make_xy(k=2):
+    X = np.array([[0.0] * k, [0.5] * k, [1.0] * k], dtype=float)
+    y = np.array([0.1, 0.2, 0.3], dtype=float)
+    return X, y
+
+
+def test_fit_regression_default_bounds_and_state_updated(monkeypatch):
+    X, y = _make_xy(k=2)
+    model = Kriging(method="regression")
+
+    seen = {"bounds": None, "x_likelihood": None}
+
+    def fake_max_likelihood(bounds):
+        seen["bounds"] = bounds
+        # k=2 thetas + 1 lambda
+        return np.array([0.1, -0.2, -6.0], dtype=float), -999.0
+
+    def fake_likelihood(x):
+        seen["x_likelihood"] = x
+        n = X.shape[0]
+        return 3.14, np.eye(n), np.eye(n)
+
+    monkeypatch.setattr(model, "max_likelihood", fake_max_likelihood, raising=True)
+    monkeypatch.setattr(model, "likelihood", fake_likelihood, raising=True)
+
+    ret = model.fit(X, y)
+
+    assert ret is model
+
+    expected_bounds = [(model.min_theta, model.max_theta)] * X.shape[1] + [
+        (model.min_Lambda, model.max_Lambda)
+    ]
+    assert seen["bounds"] == expected_bounds
+
+    assert np.allclose(model.theta, np.array([0.1, -0.2]))
+    assert np.allclose(model.Lambda, np.array([-6.0]))
+    assert np.allclose(model.logtheta_lambda_, np.array([0.1, -0.2, -6.0]))
+    assert model.negLnLike == pytest.approx(3.14)
+    assert np.allclose(model.Psi_, np.eye(X.shape[0]))
+    assert np.allclose(model.U_, np.eye(X.shape[0]))
+    assert np.allclose(seen["x_likelihood"], np.array([0.1, -0.2, -6.0]))
+
+
+def test_fit_interpolation_bounds_no_lambda_and_lambda_none(monkeypatch):
+    X, y = _make_xy(k=2)
+    model = Kriging(method="interpolation")
+
+    seen = {"bounds": None, "x_likelihood": None}
+
+    def fake_max_likelihood(bounds):
+        seen["bounds"] = bounds
+        # Only k=2 thetas
+        return np.array([0.25, -0.75], dtype=float), -1.0
+
+    def fake_likelihood(x):
+        seen["x_likelihood"] = x
+        n = X.shape[0]
+        return 1.23, np.eye(n), np.eye(n)
+
+    monkeypatch.setattr(model, "max_likelihood", fake_max_likelihood, raising=True)
+    monkeypatch.setattr(model, "likelihood", fake_likelihood, raising=True)
+
+    model.fit(X, y)
+
+    expected_bounds = [(model.min_theta, model.max_theta)] * X.shape[1]
+    assert seen["bounds"] == expected_bounds
+
+    assert np.allclose(model.theta, np.array([0.25, -0.75]))
+    assert model.Lambda is None
+    assert np.allclose(model.logtheta_lambda_, np.array([0.25, -0.75]))
+    assert model.negLnLike == pytest.approx(1.23)
+    assert np.allclose(model.Psi_, np.eye(X.shape[0]))
+    assert np.allclose(model.U_, np.eye(X.shape[0]))
+    assert np.allclose(seen["x_likelihood"], np.array([0.25, -0.75]))
+
+
+def test_fit_reinterpolation_bounds_and_state_updated(monkeypatch):
+    X, y = _make_xy(k=3)
+    model = Kriging(method="reinterpolation")
+
+    seen = {"bounds": None}
+
+    def fake_max_likelihood(bounds):
+        seen["bounds"] = bounds
+        # k=3 thetas + 1 lambda
+        return np.array([0.0, 0.1, 0.2, -3.0], dtype=float), -5.0
+
+    def fake_likelihood(x):
+        n = X.shape[0]
+        return 0.77, np.eye(n), np.eye(n)
+
+    monkeypatch.setattr(model, "max_likelihood", fake_max_likelihood, raising=True)
+    monkeypatch.setattr(model, "likelihood", fake_likelihood, raising=True)
+
+    model.fit(X, y)
+
+    expected_bounds = [(model.min_theta, model.max_theta)] * X.shape[1] + [
+        (model.min_Lambda, model.max_Lambda)
+    ]
+    assert seen["bounds"] == expected_bounds
+
+    assert np.allclose(model.theta, np.array([0.0, 0.1, 0.2]))
+    assert np.allclose(model.Lambda, np.array([-3.0]))
+    assert np.allclose(model.logtheta_lambda_, np.array([0.0, 0.1, 0.2, -3.0]))
+    assert model.negLnLike == pytest.approx(0.77)
+    assert np.allclose(model.Psi_, np.eye(X.shape[0]))
+    assert np.allclose(model.U_, np.eye(X.shape[0]))
+
+
+def test_fit_isotropic_uses_k_theta_bounds_but_sets_n_theta_1(monkeypatch):
+    X, y = _make_xy(k=3)
+    model = Kriging(method="regression", isotropic=True)
+
+    seen = {"bounds": None}
+
+    def fake_max_likelihood(bounds):
+        seen["bounds"] = bounds
+        # despite isotropic, code builds k thetas + 1 lambda bounds
+        return np.array([0.5, 0.4, 0.3, -4.0], dtype=float), -2.0
+
+    def fake_likelihood(x):
+        n = X.shape[0]
+        return 2.22, np.eye(n), np.eye(n)
+
+    monkeypatch.setattr(model, "max_likelihood", fake_max_likelihood, raising=True)
+    monkeypatch.setattr(model, "likelihood", fake_likelihood, raising=True)
+
+    model.fit(X, y)
+
+    expected_bounds = [(model.min_theta, model.max_theta)] * X.shape[1] + [
+        (model.min_Lambda, model.max_Lambda)
+    ]
+    assert seen["bounds"] == expected_bounds
+    assert model.n_theta == 1
+    # Only first theta retained in model.theta for isotropic
+    assert np.allclose(model.theta, np.array([0.5]))
+    # Lambda is taken from the first param after n_theta
+    assert np.allclose(model.Lambda, np.array([0.4]))
+    # logtheta_lambda_ still holds all optimized parameters
+    assert np.allclose(model.logtheta_lambda_, np.array([0.5, 0.4, 0.3, -4.0]))
+
+
+def test_fit_respects_explicit_bounds_override(monkeypatch):
+    X, y = _make_xy(k=2)
+    model = Kriging(method="regression")
+
+    seen = {"bounds": None, "x_likelihood": None}
+
+    def fake_max_likelihood(bounds):
+        seen["bounds"] = bounds
+        # Must match provided bounds length (2 theta + 1 lambda)
+        return np.array([1.0, -1.0, -2.0], dtype=float), -7.0
+
+    def fake_likelihood(x):
+        seen["x_likelihood"] = x
+        n = X.shape[0]
+        return 4.56, np.eye(n), np.eye(n)
+
+    monkeypatch.setattr(model, "max_likelihood", fake_max_likelihood, raising=True)
+    monkeypatch.setattr(model, "likelihood", fake_likelihood, raising=True)
+
+    custom_bounds = [(-1.0, 1.0), (-2.0, 2.0), (-8.0, -1.0)]
+    model.fit(X, y, bounds=custom_bounds)
+
+    assert seen["bounds"] == custom_bounds
+    assert np.allclose(model.theta, np.array([1.0, -1.0]))
+    assert np.allclose(model.Lambda, np.array([-2.0]))
+    assert model.negLnLike == pytest.approx(4.56)
+    assert np.allclose(seen["x_likelihood"], np.array([1.0, -1.0, -2.0]))
+
+
+def test_fit_regression_with_optim_p_adds_p_bounds_and_sets_p(monkeypatch):
+    X, y = _make_xy(k=2)
+    model = Kriging(method="regression", optim_p=True, n_p=2, min_p=1.1, max_p=1.9)
+
+    seen = {"bounds": None}
+
+    def fake_max_likelihood(bounds):
+        seen["bounds"] = bounds
+        # k=2 thetas + 1 lambda + n_p=2 p-values
+        return np.array([0.0, 0.2, -5.0, 1.3, 1.7], dtype=float), -10.0
+
+    def fake_likelihood(x):
+        n = X.shape[0]
+        return 0.5, np.eye(n), np.eye(n)
+
+    monkeypatch.setattr(model, "max_likelihood", fake_max_likelihood, raising=True)
+    monkeypatch.setattr(model, "likelihood", fake_likelihood, raising=True)
+
+    model.fit(X, y)
+
+    expected_bounds = (
+        [(model.min_theta, model.max_theta)] * X.shape[1]
+        + [(model.min_Lambda, model.max_Lambda)]
+        + [(model.min_p, model.max_p)] * model.n_p
+    )
+    assert seen["bounds"] == expected_bounds
+    assert np.allclose(model.theta, np.array([0.0, 0.2]))
+    assert np.allclose(model.Lambda, np.array([-5.0]))
+    assert np.allclose(model.p_val, np.array([1.3, 1.7]))
+
+
+def test_fit_interpolation_with_optim_p_adds_p_bounds_and_sets_p(monkeypatch):
+    X, y = _make_xy(k=2)
+    model = Kriging(method="interpolation", optim_p=True, n_p=2, min_p=1.2, max_p=1.8)
+
+    seen = {"bounds": None}
+
+    def fake_max_likelihood(bounds):
+        seen["bounds"] = bounds
+        # k=2 thetas + n_p=2 p-values
+        return np.array([0.3, -0.4, 1.25, 1.55], dtype=float), -3.0
+
+    def fake_likelihood(x):
+        n = X.shape[0]
+        return 1.0, np.eye(n), np.eye(n)
+
+    monkeypatch.setattr(model, "max_likelihood", fake_max_likelihood, raising=True)
+    monkeypatch.setattr(model, "likelihood", fake_likelihood, raising=True)
+
+    model.fit(X, y)
+
+    expected_bounds = (
+        [(model.min_theta, model.max_theta)] * X.shape[1]
+        + [(model.min_p, model.max_p)] * model.n_p
+    )
+    assert seen["bounds"] == expected_bounds
+    assert model.Lambda is None
+    assert np.allclose(model.theta, np.array([0.3, -0.4]))
+    assert np.allclose(model.p_val, np.array([1.25, 1.55]))