Skip to content

Use interpax.interp1d to support non-monotonic temperature evolution (reheating)#16

Merged
smsharma merged 1 commit intomainfrom
fix/interpax-reheating
Jan 23, 2026
Merged

Use interpax.interp1d to support non-monotonic temperature evolution (reheating)#16
smsharma merged 1 commit intomainfrom
fix/interpax-reheating

Conversation

@smsharma
Copy link
Collaborator

Summary

  • Replaces jnp.interp with interpax.interp1d in abundances.py, nuclear.py, and weak_rates.py to support reheating scenarios where photon temperature may not be monotonically decreasing
  • Arrays are now sorted by the x-coordinate before interpolation to ensure monotonicity
  • thermo.py and reactions.py are unchanged as they use pre-tabulated monotonic data

Problem

jnp.interp requires monotonically increasing x-coordinates, which breaks when T_g_vec is non-monotonic (e.g., during reheating where temperature can increase before decreasing again). The previous code used jnp.flip() assuming temperature always decreases, which doesn't work for reheating.

Solution

Use interpax.interp1d with sorted arrays:

  1. Sort arrays by the x-coordinate (temperature or energy density)
  2. Use interpax.interp1d with method='linear' and extrap=True

Test plan

  • Existing tests pass (pytest/test_abundances.py)
  • Verified interpolation works correctly with non-monotonic temperature arrays
  • Manual testing with reheating scenarios

cc @cgiovanetti (as requested)

Closes #15

🤖 Generated with Claude Code

Replace jnp.interp with interpax.interp1d in abundances.py, nuclear.py,
and weak_rates.py to support reheating scenarios where photon temperature
may not be monotonically decreasing.

jnp.interp requires monotonically increasing x-coordinates, which breaks
when T_g_vec is non-monotonic (e.g., during reheating where temperature
can increase before decreasing again). The fix sorts arrays by the
x-coordinate before interpolation.

Changes:
- Add interpax dependency to requirements.txt
- abundances.py: Sort by log(T_g) for a_start, t_start, t_end interpolations
  and by log(rho_tot) for P_tot interpolations in get_t/get_a
- weak_rates.py: Sort by Tg for neutrino-to-photon temperature ratio lookup
- nuclear.py: Sort by T_interval for weak rate interpolations

Note: thermo.py and reactions.py are unchanged as they use pre-tabulated
data that is guaranteed to be monotonic.

Closes #15

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
@cgiovanetti
Copy link
Owner

This looks great to me, thank you @smsharma ! I'm even seeing a ~25% speedup on M1 Mac (haven't tested other hardware yet)--with that in mind we should probably replace the jnp.interp's in thermo too, I see no reason not to. Passes the abundance test so I'm happy to merge pending @hongwanliu 's review.

Copy link
Owner

@cgiovanetti cgiovanetti left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Significant performance enhancement from interp1d over jnp.interp!

I'd hoped interp1d could take non-monotonic functions given their docs, but it seems like they just left that detail out; the sorting is a nice fix. Since we can still pass in T_g in any order we like this does not affect code output.

@smsharma
Copy link
Collaborator Author

Gonna merge in!

@smsharma smsharma merged commit 8d81f2b into main Jan 23, 2026
1 check passed
smsharma added a commit that referenced this pull request Jan 24, 2026
Follow-up to PR #16 which replaced jnp.interp with interpax in abundances.py.
This change applies the same improvement to thermo.py for consistency and
better performance.

Changes:
- Import interpax module
- Flip QED correction tables at load time (instead of at each call) for
  monotonically increasing x coordinates required by interpax
- Replace 6 jnp.interp calls with interpax.interp1d:
  - rho_EM_std: 2 calls for QED corrections
  - p_EM_std: 1 call for QED correction
  - rho_plus_p_EM_std: 1 call for QED correction
  - G_nue_with_me: 1 call for collision factor interpolation
  - G_numt_with_me: 1 call for collision factor interpolation

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
@cgiovanetti cgiovanetti deleted the fix/interpax-reheating branch January 26, 2026 19:52
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

jnp.interp precludes reheating scenarios

2 participants