Implementing and Testing Jacobian Functions in Stingray in ``models.py``

[kashish2210]

Contribution: Implementing and Testing Jacobian Functions in Stingray

Overview

This contribution implements and tests the Jacobian functions for GeneralizedLorentzianJacobian and SmoothBrokenPowerLawJacobian in stingray/simulator/models.py. These were listed as TODO and have now been completed and verified with test plots.

Implemented Functions

The following Jacobian functions were added: [file_change]

Testing and Verification

The functions were tested using the following code to ensure correctness and visualize the Jacobian components.

x_values = np.linspace(0.1, 10, 100)

# Test values for Generalized Lorentzian Jacobian
params_lorentz = (5.0, 1.0, 1.0, 2.0)  # (x_0, fwhm, value, power_coeff)
lorentzian_results = generalizedLorentzianJacobian(x_values, *params_lorentz)

# Test values for Smooth Broken Power Law Jacobian
params_powerlaw = (1.0, 1.5, 2.5, 3.0)  # (norm, gamma_low, gamma_high, break_freq)
powerlaw_results = SmoothBrokenPowerLawJacobian(x_values, *params_powerlaw)

# Plot results
fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# Plot Lorentzian Jacobian components
axes[0].plot(x_values, lorentzian_results[:, 0], label='d_x0')
axes[0].plot(x_values, lorentzian_results[:, 1], label='d_fwhm')
axes[0].plot(x_values, lorentzian_results[:, 2], label='d_value')
axes[0].plot(x_values, lorentzian_results[:, 3], label='d_power_coeff')
axes[0].set_title("Generalized Lorentzian Jacobian Components")
axes[0].set_xlabel("x")
axes[0].set_ylabel("Jacobian Component Values")
axes[0].legend()
axes[0].grid()

# Plot Smooth Broken Power Law Jacobian components
axes[1].plot(x_values, powerlaw_results[:, 0], label='d_norm')
axes[1].plot(x_values, powerlaw_results[:, 1], label='d_gamma_low')
axes[1].plot(x_values, powerlaw_results[:, 2], label='d_gamma_high')
axes[1].plot(x_values, powerlaw_results[:, 3], label='d_break_freq')
axes[1].set_title("Smooth Broken Power Law Jacobian Components")
axes[1].set_xlabel("x")
axes[1].set_ylabel("Jacobian Component Values")
axes[1].legend()
axes[1].grid()

plt.tight_layout()
plt.show()

Results

Helper_Notes:

• I have tried implementing @custom_model from astropy.modeling.Model and found an error

self._initialize_parameters(args, kwargs)
raise TypeError(
    ...<2 lines>...
    )

• then I found out Why Not Use @custom_model?

@custom_model expects a function returning a single array, while the Jacobian returns multiple derivatives.
@custom_model does not directly handle gradients.

This contribution addresses @matteobachetti's # TODO: Add Jacobian in stingray/simulator/models.py by implementing and testing the required Jacobian functions. The code has been validated with test plots demonstrating correct behavior.

guide me for my workflow:

[kashish2210]

for further testing with some complications:

testing code:

# Define test values
x = np.logspace(-2, 2, 100)

# Parameters for GeneralizedLorentz1D
x_0 = 1.0
fwhm = 1.0
value = 1.0
power_coeff = 2.0

# Parameters for SmoothBrokenPowerLaw
norm = 1.0
gamma_low = 1.0
gamma_high = 2.0
break_freq = 1.0

# Compute Jacobians using provided functions
jac_lorentz = GeneralizedLorentz1DJacobian(x, x_0, fwhm, value, power_coeff)
jac_smooth_bkn_po = SmoothBrokenPowerLawJacobian(x, norm, gamma_low, gamma_high, break_freq)

# Plot Jacobians
fig, ax = plt.subplots(figsize=(10, 5))
ax.loglog(x, np.abs(jac_lorentz[:, 0]), label='d/dx_0 (Lorentz)', linestyle='dashed')
ax.loglog(x, np.abs(jac_lorentz[:, 1]), label='d/dFWHM (Lorentz)', linestyle='dotted')
ax.loglog(x, np.abs(jac_lorentz[:, 2]), label='d/dValue (Lorentz)', linestyle='dashdot')
ax.loglog(x, np.abs(jac_lorentz[:, 3]), label='d/dPowerCoeff (Lorentz)', linestyle='solid')
ax.loglog(x, np.abs(jac_smooth_bkn_po[:, 0]), label='d/dNorm (Smooth BPL)', linestyle='dashed', color='red')
ax.loglog(x, np.abs(jac_smooth_bkn_po[:, 1]), label='d/dGammaLow (Smooth BPL)', linestyle='dotted', color='green')
ax.loglog(x, np.abs(jac_smooth_bkn_po[:, 2]), label='d/dGammaHigh (Smooth BPL)', linestyle='dashdot', color='purple')
ax.loglog(x, np.abs(jac_smooth_bkn_po[:, 3]), label='d/dBreakFreq (Smooth BPL)', linestyle='solid', color='brown')
ax.set_xlabel('Frequency (x)')
ax.set_ylabel('Jacobian Values (Absolute)')
ax.legend()
ax.set_title('Jacobian Values for Generalized Lorentzian and Smooth Broken Power Law')
ax.grid(True, which='both', linestyle='--', linewidth=0.5)

plt.show()

# Print Jacobian values for verification
print("Jacobian for Generalized Lorentzian (first few values):\n", jac_lorentz[:5])
print("Jacobian for Smooth Broken Power Law (first few values):\n", jac_smooth_bkn_po[:5])

output:

Jacobian for Generalized Lorentzian (first few values):
 [[ 0.327133    0.57168971  0.20323551 -0.22191759]
 [ 0.32783851  0.57284784  0.2035548  -0.22213299]
 [ 0.32861489  0.57412271  0.20390602 -0.22236968]
 [ 0.32946947  0.57552649  0.20429244 -0.22262979]
 [ 0.33041041  0.57707272  0.20471769 -0.22291567]]
Jacobian for Smooth Broken Power Law (first few values):
 [[ 1.00005000e+02  4.60540044e+02  5.00000000e-03 -9.99950004e-03]
 [ 9.11217629e+01  4.11153827e+02  5.48749382e-03 -1.09743267e-02]
 [ 8.30277791e+01  3.66908283e+02  6.02251770e-03 -1.20441617e-02]
 [ 7.56529422e+01  3.27279845e+02  6.60970573e-03 -1.32182566e-02]
 [ 6.89333748e+01  2.91797403e+02  7.25414388e-03 -1.45067611e-02]]

let me know about bugs

[codecov[bot]]

Codecov Report

Attention: Patch coverage is 20.00000% with 16 lines in your changes missing coverage. Please review.

Project coverage is 94.09%. Comparing base (a91c031) to head (d513fe6).

Files with missing lines	Patch %	Lines
stingray/simulator/models.py	20.00%	16 Missing :warning:

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #898      +/-   ##
==========================================
- Coverage   96.03%   94.09%   -1.95%     
==========================================
  Files          48       48              
  Lines        9770     9788      +18     
==========================================
- Hits         9383     9210     -173     
- Misses        387      578     +191     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.