0.29.13

Author: bartzbeielstein

pyproject.toml

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "spotpython"
-version = "0.29.12"
+version = "0.29.13"
 authors = [
   { name="T. Bartz-Beielstein", email="tbb@bartzundbartz.de" }
 ]
@@ -42,6 +42,8 @@ dependencies = [
   "nbformat",
   "pandas",
   "plotly",
+  "pytest",
+  "pytest-mock",
   "PyQt6",
   "python-markdown-math",
   "pytorch-lightning>=1.4",

src/spotpython/pinns/nn/fcn.py

@@ -0,0 +1,106 @@
+import torch
+import torch.nn as nn
+
+
+class FCN(nn.Module):
+    """A Fully Connected Network (FCN).
+
+    This network consists of an input layer, a specified number of hidden layers,
+    and an output layer. All hidden layers use the Tanh activation function.
+
+    Attributes:
+        fcs (nn.Sequential): Sequential container for the first linear layer
+                             (input to hidden) and its activation.
+        fch (nn.Sequential): Sequential container for the hidden layers. Each hidden
+                             layer consists of a linear transformation and an activation.
+        fce (nn.Linear): The final linear layer (hidden to output).
+
+    References:
+        - Solving differential equations using physics informed deep learning: a hand-on tutorial with benchmark tests. Baty, Hubert and Baty, Leo. April 2023.
+    """
+
+    def __init__(self, N_INPUT: int, N_OUTPUT: int, N_HIDDEN: int, N_LAYERS: int):
+        """Initializes the FCN.
+
+        Args:
+            N_INPUT (int): The number of input features.
+            N_OUTPUT (int): The number of output features.
+            N_HIDDEN (int): The number of neurons in each hidden layer.
+            N_LAYERS (int): The total number of layers, including the input layer
+                            (which is N_INPUT -> N_HIDDEN), hidden layers, but
+                            not counting the final output layer transformation.
+                            A N_LAYERS=1 means only input to hidden, then hidden to output.
+                            A N_LAYERS=2 means input to hidden, one hidden to hidden, then hidden to output.
+                            The number of hidden-to-hidden layers is N_LAYERS - 1.
+                            If N_LAYERS is 1, there are no fch layers.
+
+        Examples:
+            >>> # Example of creating a FCN
+            >>> from spotpython.pinns.nn.fcn import FCN
+            >>> model = FCN(N_INPUT=1, N_OUTPUT=1, N_HIDDEN=10, N_LAYERS=3)
+            >>> print(model)
+            FCN(
+              (fcs): Sequential(
+                (0): Linear(in_features=1, out_features=10, bias=True)
+                (1): Tanh()
+              )
+              (fch): Sequential(
+                (0): Sequential(
+                  (0): Linear(in_features=10, out_features=10, bias=True)
+                  (1): Tanh()
+                )
+                (1): Sequential(
+                  (0): Linear(in_features=10, out_features=10, bias=True)
+                  (1): Tanh()
+                )
+              )
+              (fce): Linear(in_features=10, out_features=1, bias=True)
+            )
+            >>> # Example of a forward pass
+            >>> input_tensor = torch.randn(5, 1) # Batch of 5, 1 input feature
+            >>> output_tensor = model(input_tensor)
+            >>> print(output_tensor.shape)
+            torch.Size([5, 1])
+
+            >>> # Example with N_LAYERS = 1 (no hidden-to-hidden layers)
+            >>> model_simple = FCN(N_INPUT=2, N_OUTPUT=1, N_HIDDEN=5, N_LAYERS=1)
+            >>> print(model_simple)
+            FCN(
+              (fcs): Sequential(
+                (0): Linear(in_features=2, out_features=5, bias=True)
+                (1): Tanh()
+              )
+              (fch): Sequential()
+              (fce): Linear(in_features=5, out_features=1, bias=True)
+            )
+        """
+        super().__init__()
+        activation = nn.Tanh
+
+        # Input layer: N_INPUT -> N_HIDDEN
+        self.fcs = nn.Sequential(nn.Linear(N_INPUT, N_HIDDEN), activation())
+
+        # Hidden layers: N_HIDDEN -> N_HIDDEN, (N_LAYERS - 1) times
+        # If N_LAYERS is 1, range(N_LAYERS-1) is range(0), so fch will be empty.
+        hidden_layers = []
+        if N_LAYERS > 1:
+            for _ in range(N_LAYERS - 1):
+                hidden_layers.append(nn.Sequential(nn.Linear(N_HIDDEN, N_HIDDEN), activation()))
+        self.fch = nn.Sequential(*hidden_layers)
+
+        # Output layer: N_HIDDEN -> N_OUTPUT
+        self.fce = nn.Linear(N_HIDDEN, N_OUTPUT)
+
+    def forward(self, x: torch.Tensor) -> torch.Tensor:
+        """Performs the forward pass of the network.
+
+        Args:
+            x (torch.Tensor): The input tensor. Shape (batch_size, N_INPUT).
+
+        Returns:
+            torch.Tensor: The output tensor. Shape (batch_size, N_OUTPUT).
+        """
+        x = self.fcs(x)
+        x = self.fch(x)
+        x = self.fce(x)
+        return x

src/spotpython/pinns/plot/result.py

@@ -0,0 +1,139 @@
+import torch
+import matplotlib.pyplot as plt
+from typing import Optional, Union, List
+import numpy as np
+
+
+def plot_result(
+    x: Union[torch.Tensor, List[float], "np.ndarray"],
+    y: Union[torch.Tensor, List[float], "np.ndarray"],
+    x_data: Union[torch.Tensor, List[float], "np.ndarray"],
+    y_data: Union[torch.Tensor, List[float], "np.ndarray"],
+    yh: Union[torch.Tensor, List[float], "np.ndarray"],
+    current_step: int,
+    xp: Optional[Union[torch.Tensor, List[float], "np.ndarray"]] = None,
+    figure_size: tuple = (8, 4),
+    xlims: Optional[tuple] = (-1.25, 31.05),
+    ylims: Optional[tuple] = (-0.65, 2.25),
+    show_plot: bool = True,
+    save_path: Optional[str] = None,
+) -> None:
+    """Plots the results of a PINN training, comparing predictions with exact solutions.
+
+    Displays the neural network's prediction, the exact solution, training data points,
+    and optionally, collocation points.
+
+    Args:
+        x (Union[torch.Tensor, List[float], "np.ndarray"]):
+            The x-coordinates for the continuous plots (e.g., time points).
+        y (Union[torch.Tensor, List[float], "np.ndarray"]):
+            The y-coordinates of the exact solution corresponding to `x`.
+        x_data (Union[torch.Tensor, List[float], "np.ndarray"]):
+            The x-coordinates of the training data points.
+        y_data (Union[torch.Tensor, List[float], "np.ndarray"]):
+            The y-coordinates of the training data points.
+        yh (Union[torch.Tensor, List[float], "np.ndarray"]):
+            The y-coordinates of the neural network's prediction corresponding to `x`.
+        current_step (int):
+            The current training step or epoch number to display on the plot.
+        xp (Optional[Union[torch.Tensor, List[float], "np.ndarray"]], optional):
+            The x-coordinates of the collocation points. If None, these are not plotted.
+            Defaults to None.
+        figure_size (tuple, optional):
+            Size of the matplotlib figure. Defaults to (8, 4).
+        xlims (Optional[tuple], optional):
+            Tuple defining the x-axis limits. If None, matplotlib's default is used.
+            Defaults to (-1.25, 31.05).
+        ylims (Optional[tuple], optional):
+            Tuple defining the y-axis limits. If None, matplotlib's default is used.
+            Defaults to (-0.65, 2.25).
+        show_plot (bool, optional):
+            Whether to display the plot using `plt.show()`. Defaults to True.
+        save_path (Optional[str], optional):
+            If provided, the path to save the figure to. If None, the figure is not saved.
+            Defaults to None.
+
+    Examples:
+        >>> from spotpython.pinns.plot.result import plot_result
+        >>> import torch
+        >>> import numpy as np
+        >>> # Generate some dummy data
+        >>> x_plot = torch.linspace(0, 30, 100)
+        >>> y_exact = torch.sin(x_plot / 5)
+        >>> y_pred = torch.sin(x_plot / 5 + 0.1) # Slightly off prediction
+        >>> x_train = torch.rand(10) * 30
+        >>> y_train = torch.sin(x_train / 5)
+        >>> collocation_points = torch.rand(50) * 30
+        >>> current_training_step = 1000
+        >>> # plot_result( # This would show a plot if run in an interactive environment
+        ... #     x_plot, y_exact, x_train, y_train, y_pred,
+        ... #     current_training_step, xp=collocation_points,
+        ... #     show_plot=False, save_path="temp_plot.png"
+        ... # )
+        >>> # To avoid actual plotting in doctest, we'll just confirm it runs
+        >>> try:
+        ...     plot_result(
+        ...         x_plot.numpy(), y_exact.numpy(), x_train.numpy(), y_train.numpy(), y_pred.numpy(),
+        ...         current_training_step, xp=collocation_points.numpy(),
+        ...         show_plot=False
+        ...     )
+        ... except Exception as e:
+        ...     print(f"Plotting failed: {e}")
+
+    Note:
+        If using PyTorch tensors as input, they will be detached and moved to CPU
+        before plotting. Consider converting to NumPy arrays beforehand if preferred.
+
+    References:
+        - Solving differential equations using physics informed deep learning: a hand-on tutorial with benchmark tests. Baty, Hubert and Baty, Leo. April 2023.
+    """
+
+    # Convert tensors to numpy arrays for plotting if they are tensors
+    def to_numpy(data):
+        if isinstance(data, torch.Tensor):
+            return data.detach().cpu().numpy()
+        return data
+
+    x_np = to_numpy(x)
+    y_np = to_numpy(y)
+    x_data_np = to_numpy(x_data)
+    y_data_np = to_numpy(y_data)
+    yh_np = to_numpy(yh)
+    if xp is not None:
+        xp_np = to_numpy(xp)
+    else:
+        xp_np = None
+
+    plt.figure(figsize=figure_size)
+    plt.plot(x_np, yh_np, color="tab:red", linewidth=2, alpha=0.8, label="NN prediction")
+    plt.plot(x_np, y_np, color="blue", linewidth=2, alpha=0.8, linestyle="--", label="Exact solution")
+    plt.scatter(x_data_np, y_data_np, s=60, color="tab:red", alpha=0.4, label="Training data")
+
+    if xp_np is not None:
+        # Create y-values for collocation points at y=0 or a specified level
+        # Original code used -0*torch.ones_like(xp), which is just zeros.
+        xp_y_values = np.zeros_like(xp_np)
+        plt.scatter(xp_np, xp_y_values, s=30, color="tab:green", alpha=0.4, label="Collocation points")
+
+    legend_handle = plt.legend(loc=(0.67, 0.62), frameon=False, fontsize="large")
+    plt.setp(legend_handle.get_texts(), color="k")
+
+    if xlims:
+        plt.xlim(xlims)
+    if ylims:
+        plt.ylim(ylims)
+
+    plt.text(0.05, 0.95, f"Training step: {current_step}", fontsize="xx-large", color="k", transform=plt.gca().transAxes, ha="left", va="top")
+
+    plt.ylabel("y", fontsize="xx-large")
+    plt.xlabel("Time", fontsize="xx-large")
+    plt.axis("on")
+    plt.grid(True, linestyle="--", alpha=0.7)
+
+    if save_path:
+        plt.savefig(save_path, dpi=300, bbox_inches="tight")
+
+    if show_plot:
+        plt.show()
+    else:
+        plt.close()  # Close the figure if not shown to free memory

src/spotpython/pinns/solvers.py

@@ -0,0 +1,82 @@
+import numpy as np
+import torch
+from typing import Tuple
+
+
+def oscillatorb(n_steps: int = 3000, t_min: float = 0.0, t_max: float = 30.0, y0: float = 1.0, alpha: float = 0.1, omega: float = np.pi / 2) -> Tuple[torch.Tensor, torch.Tensor]:
+    """Solves the first-order ODE y' = -alpha*y + sin(omega*t) using a
+    two-stage explicit Runge-Kutta method as described in the reference.
+    The ODE represents a damped harmonic oscillator with a sine forcing term.
+
+    The specific numerical scheme used is:
+    1. y_intermediate = y_current + (t_step/2) * f(t_current + t_step/2, y_current)
+    2. y_next = y_current + t_step * f(t_current + t_step, y_intermediate)
+    where f(t,y) = -alpha*y + sin(omega*t). This is a second-order method.
+
+    Args:
+        n_steps (int): Number of time points in the discretized time domain,
+            including the initial point. Defaults to 3000.
+        t_min (float): Initial time. Defaults to 0.0.
+        t_max (float): Defines the nominal end of the time interval. The time step
+            is calculated as (t_max - t_min) / n_steps. The actual last
+            time point will be t_min + (n_steps - 1) * t_step. Defaults to 30.0.
+        y0 (float): Initial condition for y at t_min. Defaults to 1.0.
+        alpha (float): Damping coefficient in the ODE. Defaults to 0.1.
+        omega (float): Angular frequency for the sine forcing term in the ODE.
+            Defaults to np.pi / 2.
+
+    Returns:
+        Tuple[torch.Tensor, torch.Tensor]: A tuple containing two PyTorch tensors:
+            - t_points_tensor: Tensor of time points, shape (n_steps, 1).
+            - y_values_tensor: Tensor of corresponding y values, shape (n_steps, 1).
+
+    Examples:
+        >>> from spotpython.pinns.solvers import oscillatorb
+        >>> import numpy as np
+        >>> import torch
+        >>> t_vals, y_vals = oscillatorb(n_steps=100, t_min=0.0, t_max=10.0, y0=1.0, alpha=0.1, omega=np.pi/2)
+        >>> print(t_vals.shape, y_vals.shape)
+        torch.Size([100, 1]) torch.Size([100, 1])
+        >>> print(f"Initial t: {t_vals[0].item():.2f}, Initial y: {y_vals[0].item():.2f}")
+        Initial t: 0.00, Initial y: 1.00
+        >>> # Last t will be t_max - t_step for this configuration
+        >>> # t_step = (10.0 - 0.0) / 100 = 0.1
+        >>> # Last t = 0.0 + (100-1)*0.1 = 9.9
+        >>> print(f"Final t: {t_vals[-1].item():.2f}, Final y: {y_vals[-1].item():.2f}")
+        Final t: 9.90, Final y: ...
+
+    References:
+        - Solving differential equations using physics informed deep learning: a hand-on tutorial with benchmark tests. Baty, Hubert and Baty, Leo. April 2023.
+    """
+    t_step = (t_max - t_min) / n_steps  # Time step
+    # Time points: t_min, t_min + t_step, ..., t_min + (n_steps-1)*t_step
+    t_points = np.arange(t_min, t_min + n_steps * t_step, t_step)[:n_steps]
+
+    y = [y0]  # List to store y values, starting with initial condition
+
+    # Solve for the time evolution
+    # t_points[0] corresponds to y0. Loop starts from t_points[1].
+    for t_current_step_end in t_points[1:]:
+        # t_midpoint is the midpoint of the current integration interval
+        # Interval: [t_current_step_end - t_step, t_current_step_end]
+        # Midpoint: (t_current_step_end - t_step) + t_step/2 = t_current_step_end - t_step/2
+        t_midpoint = t_current_step_end - t_step / 2.0
+        # y_prev is the last computed value of y
+        y_prev = y[-1]
+
+        # Stage 1: Calculate intermediate y value (y_intermediate)
+        # Uses slope at t_midpoint, with y_prev
+        # f(t,y) = -alpha*y + sin(omega*t)
+        slope_at_t_mid_using_y_prev = -alpha * y_prev + np.sin(omega * t_midpoint)
+        y_intermediate = y_prev + (t_step / 2.0) * slope_at_t_mid_using_y_prev
+
+        # Stage 2: Calculate y at t_current_step_end
+        # Uses slope at t_current_step_end, with y_intermediate
+        slope_at_t_end_using_y_intermediate = -alpha * y_intermediate + np.sin(omega * t_current_step_end)
+        y_next = y_prev + t_step * slope_at_t_end_using_y_intermediate
+        y.append(y_next)
+
+    t_points_tensor = torch.tensor(t_points, dtype=torch.float32).view(-1, 1)
+    y_values_tensor = torch.tensor(y, dtype=torch.float32).view(-1, 1)
+
+    return t_points_tensor, y_values_tensor