feat: add sombrero projection system implementation

- Introduce `sombrero_projection_system` with new configs and equations.
- Added tests, rollout functions, and integrated into Math-Gymnasium pipeline.
- Adjusted related configuration files and simulator entries for testing.

Author: RedLeader962

src/math_gymnasium/dynamical_systems/__init__.py

@@ -23,6 +23,9 @@
 from .non_linear_system.wavy_projection_system import (
     rollout_wavy_projection_partial_derivative,
 )
+from .non_linear_system.sombrero_projection_system import (
+    rollout_sombrero_projection_partial_derivative,
+)
 from .linear_debug_system import rollout_linear_debug_system
 
 __all__ = [
@@ -36,5 +39,6 @@
     "rollout_damped_oscillator_partial_derivative",
     "rollout_simple_limit_cycle_partial_derivative",
     "rollout_wavy_projection_partial_derivative",
+    "rollout_sombrero_projection_partial_derivative",
     "rollout_linear_debug_system",
 ]

src/math_gymnasium/dynamical_systems/non_linear_system/sombrero_projection_system.py

@@ -0,0 +1,119 @@
+# coding=utf-8
+
+import numpy as np
+
+from tools.math_tools.ndarray_tools.custom_msg import nan_infinity_console_warning
+from tools.math_tools.space_conversion_tools.time_to_delta_time import (
+    convert_state_time_to_state_delta_time,
+)
+
+
+def sombrero_projection_partial_derivative(
+    xyz: np.ndarray,
+    vx: float = 1.0,
+    vy: float = 1.0,
+    A: float = 1.0,
+    sigma: float = 5.0,
+    k: float = 2.0,
+    dtype: np.dtype = np.float64,
+    debug: bool = False,
+):
+    """
+    Computes the partial derivatives for a sombrero projection system.
+    This projects 2D linear motion onto a decaying ripple surface.
+
+    System equations:
+    dx/dt = vx
+    dy/dt = vy
+    dz/dt = dz/dr * dr/dt
+    where z = A * exp(-r^2 / sigma^2) * cos(k * r) and r = sqrt(x^2 + y^2)
+
+    :param xyz: An array containing the x, y, and z coordinates.
+    :param vx: Velocity in x direction.
+    :param vy: Velocity in y direction.
+    :param A: Amplitude of the ripples.
+    :param sigma: Decay constant (larger means slower decay).
+    :param k: Radial frequency of the ripples.
+    :param dtype: Data type for computations.
+    :param debug: Warn if nan or infinity values are encountered.
+    :return: An array containing the partial derivatives [x_dot, y_dot, z_dot].
+    """
+    xyz = np.nan_to_num(xyz)
+    x, y, z = xyz.astype(dtype)
+
+    r = np.sqrt(x**2 + y**2)
+
+    x_dot = vx
+    y_dot = vy
+
+    # dz/dt = A * exp(-r^2/sigma^2) * [ (-2r/sigma^2)*cos(kr) - k*sin(kr) ] * dr/dt
+    # dr/dt = (x*vx + y*vy) / r
+    # dz/dt = A * exp(-r^2/sigma^2) * [ (-2/sigma^2)*cos(kr) - k*sin(kr)/r ] * (x*vx + y*vy)
+
+    # Use sinc to handle r=0 safely: sinc(x) = sin(pi*x)/(pi*x)
+    # sin(kr)/r = k * sin(pi * (kr/pi)) / (pi * (kr/pi)) = k * sinc(kr/pi)
+    sinc_val = np.sinc(k * r / np.pi)
+
+    decay = np.exp(-(r**2) / (sigma**2))
+    bracket = (-2.0 / sigma**2) * np.cos(k * r) - k * k * sinc_val
+    z_dot = A * decay * bracket * (x * vx + y * vy)
+
+    xyz_dot = np.array([x_dot, y_dot, z_dot], dtype=dtype)
+
+    if debug and not np.all(np.isfinite(xyz_dot)):
+        nan_infinity_console_warning("xyz_dot")
+
+    xyz_dot = np.nan_to_num(xyz_dot)
+    return xyz_dot.squeeze()
+
+
+def rollout_sombrero_projection_partial_derivative(
+    time_space: np.ndarray,
+    vx: float = 1.0,
+    vy: float = 1.0,
+    A: float = 1.0,
+    sigma: float = 5.0,
+    k: float = 2.0,
+    initiale_coordinates=(0.0, 0.0, 1.0), # z(0,0) = A * exp(0) * cos(0) = A
+    time_space_is_delta_time: bool = False,
+    dtype: np.dtype = np.float64,
+) -> np.ndarray:
+    """
+    Calculates the trajectory of a sombrero projection system over a given time space.
+
+    Assume `time_space` values are increasing if `time_space_is_delta_time=True`
+
+    :param time_space: Array representing time steps in wallclock time or delta time.
+    :param vx: Velocity in x direction.
+    :param vy: Velocity in y direction.
+    :param A: Amplitude of the ripples.
+    :param sigma: Decay constant.
+    :param k: Radial frequency.
+    :param initiale_coordinates: The state at timestep 0
+    :param time_space_is_delta_time: Set to True if `time_space` is an array of delta time.
+    :param dtype: Data type for computations.
+    :return: Array of computed x, y, z coordinates over the given time space.
+    """
+    assert isinstance(initiale_coordinates, tuple) and len(initiale_coordinates) == 3
+    assert isinstance(time_space, np.ndarray) and time_space.ndim == 1
+
+    xyzs = np.empty((time_space.size, 3), dtype=dtype)
+    xyzs_partial_derivative = np.empty_like(xyzs)
+
+    xyzs[0] = initiale_coordinates
+    xyzs_partial_derivative[0] = (0.0, 0.0, 0.0)
+
+    if not time_space_is_delta_time:
+        delta_time = convert_state_time_to_state_delta_time(time_space.astype(dtype))
+    else:
+        delta_time = time_space.copy().astype(dtype)
+
+    for i in np.arange(time_space.size - 1):
+        xyzs_partial_derivative[i + 1] = sombrero_projection_partial_derivative(
+            xyzs[i], vx, vy, A, sigma, k, dtype
+        )
+        xyzs[i + 1] = xyzs[i] + xyzs_partial_derivative[i + 1] * delta_time[i + 1]
+
+    assert np.all(np.isfinite(xyzs_partial_derivative))
+    assert np.all(np.isfinite(xyzs))
+    return xyzs

tests/tests_tools/tests_non_linear_system/test_sombrero_projection_system.py

@@ -0,0 +1,64 @@
+# coding=utf-8
+from typing import Tuple
+
+import numpy as np
+import pytest
+
+from math_gymnasium.dynamical_systems.non_linear_system.sombrero_projection_system import (
+    rollout_sombrero_projection_partial_derivative,
+)
+
+
+@pytest.fixture
+def setup_sombrero_projection_sys_spaces() -> Tuple[
+    np.ndarray,
+    Tuple[float, float, float],
+]:
+    t_time_space = np.arange(10) * 0.1
+    t_init_coord = (0.0, 0.0, 1.0)
+    return t_time_space, t_init_coord
+
+
+class TestRolloutSombreroProjectionPartialDerivative:
+    def test_default(self, setup_sombrero_projection_sys_spaces):
+        (t_time_space, t_init_coord) = setup_sombrero_projection_sys_spaces
+        t_time_space_freeze = t_time_space.copy()
+        coords = rollout_sombrero_projection_partial_derivative(
+            t_time_space, vx=1.0, vy=1.0, initiale_coordinates=t_init_coord,
+        )
+        assert np.array_equal(
+            t_time_space_freeze, t_time_space
+        ), "original array was overridden"
+        assert coords.shape == (t_time_space.shape[0], 3)
+        
+        # Check if x and y are linear
+        assert np.allclose(coords[:, 0], t_time_space)
+        assert np.allclose(coords[:, 1], t_time_space)
+        
+        # Initial z should be A = 1.0
+        assert np.allclose(coords[0, 2], 1.0)
+        
+        # Check that it's not linear
+        assert not np.allclose(coords[:, 2], coords[0, 2] + (coords[1, 2] - coords[0, 2]) * np.arange(10))
+
+    def test_output_coord_using_delta_time(self, setup_sombrero_projection_sys_spaces):
+        (t_time_space, t_init_coord) = setup_sombrero_projection_sys_spaces
+        t_dt_space = np.ones_like(t_time_space) * 0.1
+        t_dt_space[0] = 0.0
+        
+        coord = rollout_sombrero_projection_partial_derivative(
+            t_dt_space, vx=1.0, vy=1.0, initiale_coordinates=t_init_coord, time_space_is_delta_time=True
+        )
+        assert coord.shape == (t_time_space.shape[0], 3)
+        # Initial z is 1.0.
+        # At r=0, dz/dt = 0 (as calculated in thoughts, since x*vx+y*vy = 0).
+        # Actually, let's check the first step from origin.
+        # If x=0, y=0, then x*vx + y*vy = 0, so dz/dt = 0.
+        # So the second z should still be 1.0 if starting EXACTLY at origin.
+        assert np.allclose(coord[1, 2], 1.0)
+        
+        # If we start offset, it should change.
+        coord_offset = rollout_sombrero_projection_partial_derivative(
+            t_dt_space, vx=1.0, vy=1.0, initiale_coordinates=(0.1, 0.1, 1.0), time_space_is_delta_time=True
+        )
+        assert not np.allclose(coord_offset[1, 2], 1.0)