Rate vs readout#

coronagraphoto exposes two parallel families of simulation functions for every source type (star, planet, disk, zodi, and the whole-scene aggregator system):

Family	Returns	Stochastic?	Differentiable?
`*_rate`	Per-pixel count rate (e/s/pixel)	No	Yes (through the whole forward model)
`*_readout`	Per-pixel electron count (e) for one exposure	Yes (Poisson + binomial QE)	No

This page explains when to reach for which.

The mathematics#

For a single source \(S\) through the optical path \(\mathcal{O}\), the per-pixel expected electron count during an exposure of length \(t\) on a detector with quantum efficiency \(\eta\) is:

\[ \mathbb{E}[N_{\text{electrons}}] \;=\; r_S(\mathcal{O}) \cdot t \cdot \eta \]

where \(r_S(\mathcal{O})\) is the count rate (electrons per second per pixel). The *_rate functions return \(r_S(\mathcal{O})\).

The actual observed electron count is a Poisson realisation of that expectation. The *_readout functions return:

\[ N_{\text{electrons}} \;\sim\; \mathrm{Poisson}(r_S(\mathcal{O}) \cdot t \cdot \eta) \]

For the high-count regime coronagraphoto also offers a “thinned” variant that uses the Poisson thinning theorem (\(\mathrm{Binomial}(\mathrm{Poisson}(\lambda), p) \sim \mathrm{Poisson}(\lambda p)\)) to fuse the QE binomial draw into the Poisson draw, which is faster without changing the distribution. That’s the readout_source_electrons_thinned path inside the detector model.

Use `*_rate` for fitting#

Anything that needs gradients should differentiate through the rate pipeline. Examples:

Likelihood evaluation – the rate map IS the model expectation; compare against observed counts via a Poisson or Gaussian likelihood.
MAP / HMC retrieval – gradient-based parameter inference (planet position, atmospheric composition, orbital elements).
Sensitivity studies – \(\partial(\text{observable}) / \partial(\text{parameter})\) for, e.g., the dependence of detection significance on telescope diameter or exposure time.
Yield calculator forward model – the deterministic rate-summed-over-targets that downstream optimisers query.

import jax
import equinox as eqx
from coronagraphoto import system_rate

# Differentiable rate pipeline -- gradient flows end-to-end
@eqx.filter_jit
def total_electrons(wavelength_nm, scene, optical_path):
    rate_map = system_rate(
        scene, optical_path,
        start_time_jd=2_460_000.0,
        wavelength_nm=wavelength_nm,
        bin_width_nm=50.0,
        telescope_pa_deg=0.0,
        ecliptic_lat_deg=0.0, solar_lon_deg=135.0,
    )
    return rate_map.sum() * EXPOSURE_S

grad_fn = eqx.filter_grad(total_electrons)

Differentiable parameters can be any JAX array inside scene or optical_path – use eqx.tree_at to swap a field, or differentiate directly w.r.t. a kwarg like wavelength_nm. The Richardson-converged gradient tests in the test suite cover the wavelength case end-to-end.

The system_rate orchestrator sums all per-source rates – star, every planet, disk if present, zodi if present – into one differentiable map for the entire scene.

Use `*_readout` for data generation#

Anything that needs a noisy detector image should use the readout family. Examples:

Generating an HWO survey dataset – thousands of frames, each with independent Poisson noise.
End-to-end ADI / RDI demos – generate noisy science frames at varying telescope rotations, then run KLIP on the cube.
Comparing post-processing algorithms – the algorithm should see realistic noise, not the noiseless rate.

from coronagraphoto import system_readout
import jax

key = jax.random.PRNGKey(0)
image = system_readout(
    scene, optical_path, key,
    start_time_jd=..., exposure_time_s=..., wavelength_nm=...,
    bin_width_nm=..., telescope_pa_deg=...,
    ecliptic_lat_deg=..., solar_lon_deg=...,
)

system_readout independently draws Poisson noise per source (each gets its own PRNG subkey) and sums the realisations. The sum of Poissons is Poisson with summed mean, so this is statistically equivalent to one combined draw, but the per-source split gives you deterministic noise for the same key + scene combination.

What you cannot do#

You cannot differentiate through *_readout. The Poisson draw is a discrete random integer and its gradient is zero almost everywhere. JAX will not raise; it will silently return zero gradients. If you find yourself doing jax.grad(simulate_frame) where simulate_frame calls *_readout, the result is meaningless. Two acceptable workarounds:

Differentiate the rate instead. The expected count \(r_S \cdot t \cdot \eta\) is what almost every inference problem actually wants, and it’s differentiable.
Use a Gaussian approximation at high counts: \(\mathrm{Poisson}(\lambda) \approx \mathcal{N}(\lambda, \sqrt{\lambda})\). The Gaussian distribution is reparameterisable and differentiable through the noise sample. Build it yourself if you need it.

Performance#

Both families share most of the work – the readout family is just the rate family followed by a single Poisson draw via the detector’s readout method. The Poisson draw is a few percent of the per-frame cost on GPU. See performance for the closure-vs-args trap that dominates compile time and how to avoid it.

Rate vs readout

Contents