Implementing and Testing Jacobian Functions in Stingray in models.py

CHANGELOG.rst

@@ -1,3 +1,38 @@
+v2.2.7.dev43+gd513fe6e.d20250312 (2025-03-12)
+---------------------------------------------
+
+New Features
+^^^^^^^^^^^^
+
+- Better support for Chandra event files (`#877 <https://github.com/StingraySoftware/stingray/pull/877>`__)
+- Added Jacobian functions for Generalized Lorentzian and Smooth Broken Power Law models.
+
+  - Introduced `GeneralizedLorentz1DJacobian` to compute analytical derivatives for the Generalized Lorentzian function.
+  - Implemented `SmoothBrokenPowerLawJacobian` to enhance model fitting for Smooth Broken Power Law.
+  - Improved numerical stability in `stingray.simulator.models` by providing Jacobian support.
+  - Added new test cases ensuring the correctness of Jacobian computations. (`#898 <https://github.com/StingraySoftware/stingray/pull/898>`__)
+
+
+Bug Fixes
+^^^^^^^^^
+
+- Fixed issue with counting ``0`` bin occupancy in ``fold_events`` with ``mode='pdm'`` (`#872 <https://github.com/StingraySoftware/stingray/pull/872>`__)
+
+
+Documentation
+^^^^^^^^^^^^^
+
+- Added `stingray.varenergyspectrum.CountSpectrum` to docs. (`#884 <https://github.com/StingraySoftware/stingray/pull/884>`__)
+
+
+Internal Changes
+^^^^^^^^^^^^^^^^
+
+- Bump jinja2 version to 3.1.5 (`#878 <https://github.com/StingraySoftware/stingray/pull/878>`__)
+- Added a test to check ``profile``, when ``fold_events`` method called with ``mode= "pdm"`` after the bug fix #872. (`#880 <https://github.com/StingraySoftware/stingray/pull/880>`__)
+- Suppressed deprecation warnings related to `numpy.core.einsumfunc` as it is causing test failure described in (`#886 <https://github.com/StingraySoftware/stingray/pull/886>`__).
+
+
 v2.3 (under development)
 ------------------------
 

stingray/simulator/models.py

@@ -1,4 +1,54 @@
 from astropy.modeling.models import custom_model
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# TODO: Added Jacobian functions
+def GeneralizedLorentz1DJacobian(
+    x: np.ndarray, x_0: float, fwhm: float, value: float, power_coeff: float
+) -> np.ndarray:
+    """
+    Compute the Jacobian matrix for the Generalized Lorentzian function.
+
+    Parameters
+    ----------
+    x : numpy.ndarray
+        Non-zero frequencies.
+    x_0 : float
+        Peak central frequency.
+    fwhm : float
+        Full width at half maximum (FWHM) of the peak.
+    value : float
+        Peak value at x = x_0.
+    power_coeff : float
+        Power coefficient.
+
+    Returns
+    -------
+    numpy.ndarray
+        The computed Jacobian matrix of shape (len(x), 4).
+    """
+    dx = x - x_0
+    gamma_pow = (fwhm / 2) ** power_coeff
+    denom = dx**power_coeff + gamma_pow
+    denom_sq = denom**2
+
+    d_x0 = power_coeff * value * gamma_pow * dx ** (power_coeff - 1) / denom_sq
+    d_fwhm = (
+        power_coeff
+        * value
+        * (fwhm / 2) ** (power_coeff - 1)
+        / 2
+        * (1 + power_coeff * gamma_pow / denom)
+        / denom
+    )
+    d_value = gamma_pow / denom
+    d_power_coeff = (
+        value * gamma_pow * np.log(fwhm / 2) / denom
+        - value * gamma_pow * np.log(abs(dx) + (fwhm / 2)) / denom
+    )
+
+    return np.vstack([-d_x0, d_fwhm, d_value, d_power_coeff]).T
 
 
 # TODO: Add Jacobian
@@ -39,6 +89,48 @@ def GeneralizedLorentz1D(x, x_0=1.0, fwhm=1.0, value=1.0, power_coeff=1.0):
     )
 
 
+# TODO: Added Jacobian functions
+def SmoothBrokenPowerLawJacobian(
+    x: np.ndarray, norm: float, gamma_low: float, gamma_high: float, break_freq: float
+) -> np.ndarray:
+    """
+    Compute the Jacobian matrix for the Smooth Broken Power Law function.
+
+    Parameters
+    ----------
+    x : numpy.ndarray
+        Non-zero frequencies.
+    norm : float
+        Normalization frequency.
+    gamma_low : float
+        Power law index for f → zero.
+    gamma_high : float
+        Power law index for f → infinity.
+    break_freq : float
+        Break frequency.
+
+    Returns
+    -------
+    numpy.ndarray
+        The computed Jacobian matrix of shape (len(x), 4).
+    """
+    x_bf2 = (x / break_freq) ** 2
+    denom = (1.0 + x_bf2) ** (-(gamma_low - gamma_high) / 2)
+    d_norm = x ** (-gamma_low) * denom
+    d_gamma_low = -norm * np.log(x) * x ** (-gamma_low) * denom
+    d_gamma_high = norm * x ** (-gamma_low) * denom * np.log(1 + x_bf2) / 2
+    d_break_freq = (
+        norm
+        * x ** (-gamma_low)
+        * denom
+        * (gamma_low - gamma_high)
+        * x_bf2
+        / (break_freq * (1 + x_bf2))
+    )
+
+    return np.vstack([d_norm, d_gamma_low, d_gamma_high, d_break_freq]).T
+
+
 # TODO: Add Jacobian
 @custom_model
 def SmoothBrokenPowerLaw(x, norm=1.0, gamma_low=1.0, gamma_high=1.0, break_freq=1.0):

stingray/simulator/tests/test_simulator.py

@@ -8,6 +8,8 @@
 from stingray import Lightcurve, Crossspectrum, sampledata, Powerspectrum
 from stingray.simulator import Simulator
 from stingray.simulator import models
+from stingray.simulator.models import GeneralizedLorentz1D, SmoothBrokenPowerLaw
+from stingray.simulator.models import GeneralizedLorentz1DJacobian, SmoothBrokenPowerLawJacobian
 
 _H5PY_INSTALLED = True
 
@@ -612,3 +614,80 @@ def test_io_with_unsupported_format(self):
         with pytest.raises(KeyError):
             sim.read("sim.pickle", fmt="hdf5")
         os.remove("sim.pickle")
+
+    def test_GeneralizedLorentz1DJacobian(self):
+        """
+        Test the Jacobian function for GeneralizedLorentz1D model.
+        """
+        x = np.array([0.5, 1.0, 2.0, 5.0])
+        x_0, fwhm, value, power_coeff = 10, 1.0, 10.0, 2
+        jacobian = GeneralizedLorentz1DJacobian(x, x_0, fwhm, value, power_coeff)
+
+        assert jacobian.shape == (len(x), 4), "Jacobian shape mismatch"
+        assert np.all(np.isfinite(jacobian)), "Jacobian contains NaN/Inf"
+
+    def test_SmoothBrokenPowerLawJacobian(self):
+        """
+        Test the Jacobian function for SmoothBrokenPowerLaw model.
+        """
+        x = np.array([0.1, 1.0, 10.0, 100.0])
+        norm, gamma_low, gamma_high, break_freq = 1.0, 1.0, 2.0, 1.0
+        jacobian = SmoothBrokenPowerLawJacobian(x, norm, gamma_low, gamma_high, break_freq)
+
+        assert jacobian.shape == (len(x), 4), "Jacobian shape mismatch"
+        assert np.all(np.isfinite(jacobian)), "Jacobian contains NaN/Inf"
+
+    def test_simulate_GeneralizedLorentz1D_str(self):
+        """
+        Simulate a light curve using the GeneralizedLorentz1D model
+        called as a string
+        """
+        assert len(
+            self.simulator.simulate(
+                "GeneralizedLorentz1D", {"x_0": 10, "fwhm": 1.0, "value": 10.0, "power_coeff": 2}
+            )
+        ), 1024
+
+    def test_simulate_SmoothBrokenPowerLaw_str(self):
+        """
+        Simulate a light curve using the SmoothBrokenPowerLaw model
+        called as a string
+        """
+        assert len(
+            self.simulator.simulate(
+                "SmoothBrokenPowerLaw",
+                {"norm": 1.0, "gamma_low": 1.0, "gamma_high": 2.0, "break_freq": 1.0},
+            )
+        ), 1024
+
+    @pytest.mark.parametrize("poisson", [True, False])
+    def test_compare_composite(self, poisson):
+        """
+        Compare the PSD of a light curve simulated using a composite model
+        (using SmoothBrokenPowerLaw plus GeneralizedLorentz1D)
+        with the actual model
+        """
+        N = 50000
+        dt = 0.01
+        m = 30000.0
+
+        self.simulator = Simulator(N=N, mean=m, dt=dt, rms=self.rms, poisson=poisson)
+        smoothbknpo = SmoothBrokenPowerLaw(
+            norm=1.0, gamma_low=1.0, gamma_high=2.0, break_freq=1.0
+        )
+        lorentzian = GeneralizedLorentz1D(x_0=10, fwhm=1.0, value=10.0, power_coeff=2.0)
+        myModel = smoothbknpo + lorentzian
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            lc = [self.simulator.simulate(myModel) for i in range(1, 50)]
+
+        simulated = self.simulator.powerspectrum(lc, lc[0].tseg)
+
+        w = np.fft.rfftfreq(N, d=dt)[1:]
+        actual = myModel(w)[:-1]
+
+        actual_prob = actual / float(sum(actual))
+        simulated_prob = simulated / float(sum(simulated))
+
+        assert np.all(np.abs(actual_prob - simulated_prob) < 3 * np.sqrt(actual_prob))