diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 436667f21..3f9b71930 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,3 +1,38 @@ +v2.2.7.dev43+gd513fe6e.d20250312 (2025-03-12) +--------------------------------------------- + +New Features +^^^^^^^^^^^^ + +- Better support for Chandra event files (`#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 `__) + + +Bug Fixes +^^^^^^^^^ + +- Fixed issue with counting ``0`` bin occupancy in ``fold_events`` with ``mode='pdm'`` (`#872 `__) + + +Documentation +^^^^^^^^^^^^^ + +- Added `stingray.varenergyspectrum.CountSpectrum` to docs. (`#884 `__) + + +Internal Changes +^^^^^^^^^^^^^^^^ + +- Bump jinja2 version to 3.1.5 (`#878 `__) +- Added a test to check ``profile``, when ``fold_events`` method called with ``mode= "pdm"`` after the bug fix #872. (`#880 `__) +- Suppressed deprecation warnings related to `numpy.core.einsumfunc` as it is causing test failure described in (`#886 `__). + + v2.3 (under development) ------------------------ diff --git a/stingray/simulator/models.py b/stingray/simulator/models.py index c209d5a29..82b06d839 100644 --- a/stingray/simulator/models.py +++ b/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): diff --git a/stingray/simulator/tests/test_simulator.py b/stingray/simulator/tests/test_simulator.py index 053f37b06..d2755dbf6 100644 --- a/stingray/simulator/tests/test_simulator.py +++ b/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))