0.26.23

simple plot method

Author: bartzbeielstein

notebooks/00_spotPython_tests.ipynb

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

pyproject.toml

@@ -66,7 +66,7 @@ def branin(x: np.ndarray) -> np.ndarray:
 
     Examples:
         >>> import numpy as np
-        >>> from spotpython.surrogate.functions.branin import branin
+        >>> from spotpython.surrogate.functions.forr08a import branin
         >>> # Single input
         >>> print(branin(np.array([[0.5, 0.5]])))
         [26.63]

src/spotpython/surrogate/functions/forr08a.py

@@ -4,6 +4,8 @@
 from scipy.optimize import differential_evolution
 from sklearn.base import BaseEstimator, RegressorMixin
 from scipy.special import erf
+import matplotlib.pyplot as plt
+from numpy import linspace, meshgrid, array
 
 
 class Kriging(BaseEstimator, RegressorMixin):
@@ -53,6 +55,8 @@ def __init__(self, eps: float = None, penalty: float = 1e4, method="regression")
         if method not in ["interpolation", "regression", "reinterpolation"]:
             raise ValueError("method must be one of 'interpolation', 'regression', or 'reinterpolation']")
         self.method = method
+        self.return_ei = False
+        self.return_std = False
 
     def _get_eps(self) -> float:
         """
@@ -159,6 +163,8 @@ def predict(self, X: np.ndarray, return_std=False, return_ei=False) -> np.ndarra
             >>> print("Standard deviations:", sd)
             >>> print("Expected improvement:", ei)
         """
+        self.return_std = return_std
+        self.return_ei = return_ei
         X = np.atleast_2d(X)
         if return_std and return_ei:
             # Return predictions, standard deviations, and expected improvements
@@ -325,12 +331,14 @@ def _pred(self, x: np.ndarray) -> float:
         f = mu + psi @ resid_tilde
 
         # Compute ExpImp
-        yBest = np.min(y)
-        EITermOne = (yBest - f) * (0.5 + 0.5 * erf((1 / np.sqrt(2)) * ((yBest - f) / s)))
-        EITermTwo = s * (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((yBest - f) ** 2 / SSqr))
-        ExpImp = np.log10(EITermOne + EITermTwo + self.eps)
-
-        return float(f), float(s), float(-ExpImp)
+        if self.return_ei:
+            yBest = np.min(y)
+            EITermOne = (yBest - f) * (0.5 + 0.5 * erf((1 / np.sqrt(2)) * ((yBest - f) / s)))
+            EITermTwo = s * (1 / np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((yBest - f) ** 2 / SSqr))
+            ExpImp = np.log10(EITermOne + EITermTwo + self.eps)
+            return float(f), float(s), float(-ExpImp)
+        else:
+            return float(f), float(s)
 
     def max_likelihood(self, bounds: List[Tuple[float, float]]) -> Tuple[np.ndarray, float]:
         """
@@ -351,3 +359,103 @@ def objective(logtheta_lambda):
 
         result = differential_evolution(objective, bounds)
         return result.x, result.fun
+
+    def plot(self, show: Optional[bool] = True, alpha=0.8) -> None:
+        """
+        This function plots 1D and 2D surrogates.
+
+        Args:
+            show (bool):
+                If `True`, the plots are displayed.
+                If `False`, `plt.show()` should be called outside this function.
+
+        Returns:
+            None
+
+        Examples:
+            >>> model = Kriging()
+            >>> model.fit(X_train, y_train)
+            >>> model.plot()
+        """
+        if self.X_ is None or self.y_ is None:
+            raise ValueError("The model must be fitted before calling the plot method.")
+
+        k = self.X_.shape[1]  # Number of dimensions
+
+        if k == 1:
+            # 1D Plot
+            fig = plt.figure(figsize=(9, 6))
+            n_grid = 100
+            x = linspace(self.X_[:, 0].min(), self.X_[:, 0].max(), num=n_grid).reshape(-1, 1)
+            y_pred, y_std = self.predict(x, return_std=True)
+
+            plt.plot(x, y_pred, "k", label="Prediction")
+            plt.fill_between(
+                x.ravel(),
+                y_pred - 1.96 * y_std,
+                y_pred + 1.96 * y_std,
+                alpha=0.2,
+                label="95% Confidence Interval",
+            )
+            plt.scatter(self.X_, self.y_, color="red", label="Training Data")
+            plt.xlabel("X")
+            plt.ylabel("Prediction")
+            plt.title("1D Kriging Surrogate")
+            plt.legend()
+            if show:
+                plt.show()
+
+        elif k == 2:
+            # 2D Plot
+            fig = plt.figure(figsize=(12, 10))
+            n_grid = 100
+            x1 = linspace(self.X_[:, 0].min(), self.X_[:, 0].max(), num=n_grid)
+            x2 = linspace(self.X_[:, 1].min(), self.X_[:, 1].max(), num=n_grid)
+            X1, X2 = meshgrid(x1, x2)
+            grid_points = array([X1.ravel(), X2.ravel()]).T
+
+            y_pred, y_std = self.predict(grid_points, return_std=True)
+            Z_pred = y_pred.reshape(X1.shape)
+            Z_std = y_std.reshape(X1.shape)
+
+            # Plot predicted values
+            ax1 = fig.add_subplot(221, projection="3d")
+            ax1.plot_surface(X1, X2, Z_pred, cmap="viridis", alpha=alpha)
+            ax1.set_title("Prediction Surface")
+            ax1.set_xlabel("X1")
+            ax1.set_ylabel("X2")
+            ax1.set_zlabel("Prediction")
+
+            # Plot prediction error
+            ax2 = fig.add_subplot(222, projection="3d")
+            ax2.plot_surface(X1, X2, Z_std, cmap="viridis", alpha=alpha)
+            ax2.set_title("Prediction Error Surface")
+            ax2.set_xlabel("X1")
+            ax2.set_ylabel("X2")
+            ax2.set_zlabel("Error")
+
+            # Contour plot of predicted values
+            ax3 = fig.add_subplot(223)
+            contour = ax3.contourf(X1, X2, Z_pred, cmap="viridis", levels=30)
+            plt.colorbar(contour, ax=ax3)
+            ax3.scatter(self.X_[:, 0], self.X_[:, 1], color="red", label="Training Data")
+            ax3.set_title("Prediction Contour")
+            ax3.set_xlabel("X1")
+            ax3.set_ylabel("X2")
+            ax3.legend()
+
+            # Contour plot of prediction error
+            ax4 = fig.add_subplot(224)
+            contour = ax4.contourf(X1, X2, Z_std, cmap="viridis", levels=30)
+            plt.colorbar(contour, ax=ax4)
+            ax4.scatter(self.X_[:, 0], self.X_[:, 1], color="red", label="Training Data")
+            ax4.set_title("Error Contour")
+            ax4.set_xlabel("X1")
+            ax4.set_ylabel("X2")
+            ax4.legend()
+
+            if show:
+                plt.show()
+
+        else:
+            raise ValueError("Plotting is only supported for 1D or 2D input data.")

src/spotpython/surrogate/kriging.py

@@ -620,7 +620,7 @@ def mmlhs(X_start: np.ndarray, population: int, iterations: int, q: Optional[flo
     return X_best
 
 
-def bestlh(n: int, k: int, population: int, iterations: int, p=1, plot=False) -> np.ndarray:
+def bestlh(n: int, k: int, population: int, iterations: int, p=1, plot=False, verbosity=0) -> np.ndarray:
     """
     Generates an optimized Latin hypercube by evolving the Morris-Mitchell
     criterion across multiple exponents (q values) and selecting the best plan.
@@ -640,6 +640,8 @@ def bestlh(n: int, k: int, population: int, iterations: int, p=1, plot=False) ->
         plot (bool, optional):
             If True, a scatter plot of the optimized plan in the first two dimensions
             will be displayed. Only if k>=2.  Defaults to False.
+        verbosity (int, optional):
+            Verbosity level. 0 is silent, 1 prints the best q value found. Defaults to 0.
 
     Returns:
         np.ndarray:
@@ -681,15 +683,17 @@ def bestlh(n: int, k: int, population: int, iterations: int, p=1, plot=False) ->
 
     # Evolve the plan for each q in q_list
     for i, q_val in enumerate(q_list):
-        print(f"Now optimizing for q={q_val}...")
+        if verbosity > 0:
+            print(f"Now optimizing for q={q_val}...")
         X3D[:, :, i] = mmlhs(X_start, population, iterations, q_val)
 
     # Sort the set of evolved plans according to the Morris-Mitchell criterion
     index_order = mmsort(X3D, p=p)
 
     # index_order is a 1-based array of plan indices; the first element is the best
     best_idx = index_order[0] - 1
-    print(f"Best lh found using q={q_list[best_idx]}...")
+    if verbosity > 0:
+        print(f"Best lh found using q={q_list[best_idx]}...")
 
     # The best plan in 3D array order
     X = X3D[:, :, best_idx]