Skip to content

Add optika.systems.LinearSystem, a linear forward-model approximation to SequentialSystem - #152

Open
roytsmart wants to merge 25 commits into
mainfrom
feature/interpolated-system
Open

Add optika.systems.LinearSystem, a linear forward-model approximation to SequentialSystem#152
roytsmart wants to merge 25 commits into
mainfrom
feature/interpolated-system

Conversation

@roytsmart

@roytsmart roytsmart commented Apr 3, 2026

Copy link
Copy Markdown
Collaborator

Summary

Adds optika.systems.LinearSystem (and its abstract base AbstractLinearSystem), a fast linear forward model that approximates a SequentialSystem.

Instead of raytracing every scene, a LinearSystem is characterized once by a set of precomputed models:

  • a distortion model mapping object-plane coordinates onto the detector,
  • an area_effective (effective area) model,
  • optional vignetting and field_stop models,
  • and a sensor.

image() maps a scene's energy spectral radiance to detector electrons by conservatively regridding it through the distortion model, weighting each cell by the effective area, vignetting, and field stop, then applying the sensor's exposure model. The expensive conservative-regridding operator is built once in weights() and reused by image_from_weights(), so many scenes can be imaged through the same optics cheaply.

Supporting changes

  • Shaped on the calibration models. AbstractDistortionModel, AbstractVignettingModel, and AbstractEffectiveAreaModel now inherit the Shaped mixin, with shape implemented on the concrete models (SimpleDistortionModel, PolynomialDistortionModel, PolynomialVignettingModel, InterpolatedEffectiveAreaModel). Each returns the model's batch shape (excluding the wavelength/field domain axes), which is what lets LinearSystem.shape report the system's parametrization.
  • Moved Plottable/DxfWritable from AbstractSystem to AbstractSequentialSystem. A LinearSystem has no geometry to plot or export to CAD, so the common AbstractSystem base no longer forces those mixins; the geometric SequentialSystem keeps them.

Tests

_linear_test.py adds TestLinearSystem with two fixtures (a minimal system, and one with vignetting and a field stop), driving image() end-to-end with noise on and off, plus the model interface. The distortion/vignetting/effective-area test bases gained shape coverage.

Docs

LinearSystem has a runnable .. jupyter-execute:: example that assembles the component models and images a USAF-1951 target, mirroring the SequentialSystem example. AbstractLinearSystem documents the full model, and a Features bullet was added to the landing page.

Follow-ups (known limitations)

  • A scene with only a single wavelength cell currently fails in image(): the size-1 wavelength axis is dropped and the regridding cannot align the weights.
  • LinearSystem.image expects energy spectral radiance, whereas SequentialSystem.image takes photon radiance; the two image() methods disagree on units.

🤖 Generated with Claude Code

https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ

@codecov

codecov Bot commented Apr 3, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.57%. Comparing base (71e7cd6) to head (3b75e6e).

Additional details and impacted files
@@            Coverage Diff             @@
##             main     #152      +/-   ##
==========================================
+ Coverage   99.55%   99.57%   +0.02%     
==========================================
  Files         116      118       +2     
  Lines        6457     6805     +348     
==========================================
+ Hits         6428     6776     +348     
  Misses         29       29              
Flag Coverage Δ
unittests 99.57% <100.00%> (+0.02%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 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.

@roytsmart
roytsmart force-pushed the feature/interpolated-system branch from 4b5ad3a to c625f4b Compare July 6, 2026 16:36
@roytsmart
roytsmart force-pushed the feature/interpolated-system branch from 23ec8b0 to 33b7cd0 Compare July 21, 2026 22:16
Rewrite the `AbstractLinearSystem` docstring to describe the full linear
forward model (distortion, effective area, and optional vignetting/field
stop) instead of the stale "define a distortion/vignetting method" text, and
add a `LinearSystem` class docstring with a `jupyter-execute` example that
assembles the component models and images a USAF-1951 target, mirroring the
`SequentialSystem` example.

Also document the `shape` property, clarify that `image()` expects energy
spectral radiance, and add a Features bullet for the linear forward model.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
@roytsmart roytsmart changed the title Added optika.systems.InterpolatedSystem as an approximation to SequentialSystem. Add optika.systems.LinearSystem, a linear forward-model approximation to SequentialSystem Jul 21, 2026
roytsmart and others added 4 commits July 21, 2026 20:22
Fit the system's distortion, vignetting, and effective area models and
assemble them (with the sensor) into a `LinearSystem`. The `pupil` argument is
interpreted as cell vertices: the effective area fit uses the vertices to
compute pupil cell areas, while the distortion and vignetting fits trace at the
corresponding cell centers. The field stop is not modeled, so `field_stop` is
left `None`.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
`LinearSystem.image` no longer hardcodes an energy-to-photon conversion; it
passes the flux through in whatever units the scene radiance is given (energy
or photon) and lets the sensor convert, matching `SequentialSystem`.

Add `backproject` / `backproject_from_weights`, the transpose of the optical
forward model, which projects a detector-plane image back onto the object
plane and recovers the scene radiance (relying on the corrected
`transpose_weights_conservative`).

Fix the `direction` computed by `SequentialSystem.linearize`: use the refracted
cosine from `rayfunction_default` (as `sensor.collect` does) rather than the
uncalled `rayfunction` method and the raw direction vector, averaging only over
the grid axes actually present.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
roytsmart and others added 6 commits July 23, 2026 10:00
`backproject` now accepts a detector image of electrons and inverts the
sensor response with `ImagingSensor.photons_absorbed` before applying the
optical transpose, making `image`/`backproject` a matched pair. Add a
`test_backproject` roundtrip.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
…t_from_weights`

`image_from_weights` now takes a scene and returns electrons (integrating
the radiance over each voxel, regridding, and exposing the sensor);
`backproject_from_weights` takes the detector image plus object-plane
coordinates and returns the backprojected radiance (inverting the sensor
response, regridding, and dividing out the voxel volume). `image` and
`backproject` become thin wrappers that build the weights and delegate,
so callers holding precomputed weights (e.g. ctis) can reuse the full
forward/transpose model without recomputing the regridding operator.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
When `uncertainty=True`, `image`/`image_from_weights` attach the standard
deviation of the measurement noise to the result as a
`NormalUncertainScalarArray`, using the sensor's `uncertainty` method. This
lets a caller obtain the expected image and its per-pixel noise width in a
single pass, which is needed to compute the uncertainty of a
wavelength-integrated image before the integration discards the
per-wavelength electrons.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
roytsmart and others added 9 commits July 23, 2026 18:50
Read noise is a per-readout effect, but the CTIS `LinearSystem` forward
model calls `expose` per wavelength, so adding read noise there
over-counted it by a factor of sqrt(N_wavelength).

Give `expose`, `uncertainty`, `photons_absorbed`, and `measure` an
`integrate` keyword (default `True`): the shot/Fano/QY conversion stays
per wavelength (exact, and correct for large EUV quantum yields), then the
electrons are summed over wavelength into a single readout and the read
noise is applied once. `expose` also gains an `uncertainty` flag so it can
return a `NormalUncertainScalarArray` in one pass (evaluating `signal`
once). `LinearSystem.image`/`backproject` forward `integrate`;
`SequentialSystem.image` now delegates its wavelength integration to
`expose` instead of collapsing the wavelength grid itself.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
Declare the shared forward/transpose interface on `AbstractSystem`: both
`image` and `backproject` take the common `scene`/`image`,
`axis_wavelength`, `axis_field`, `integrate` (and `noise` for `image`)
arguments in the same order. `SequentialSystem.image`, `LinearSystem.image`,
and `LinearSystem.image_from_weights`/`backproject_from_weights` are all
reordered to that common prefix (in both signatures and docstrings), with
their extra arguments (`pupil`/`axis_pupil`, `uncertainty`, `weights`) moved
to the end. `SequentialSystem` gains a `backproject` that raises
`NotImplementedError`, since a ray-traced system is not a linear operator.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
Add `AbstractImagingSensor.pixels()`, mapping an in-plane sensor position to
fractional pixel coordinates (pixel 0 at the lower edge of the
light-sensitive area). `SequentialSystem.distortion` now fits the sensor-side
coordinates in pixels, and `LinearSystem.coordinates_sensor` returns the
pixel edges in `u.pix`, so the distortion (plate scale in arcsec/pix,
dispersion in nm/pix) and the regrid grid share pixel units. Ray tracing and
`collect` keep the physical units they naturally produce.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
`PolynomialDistortionModel.fit`/`fit_inverse` and
`PolynomialVignettingModel.fit` did not pass `axis_polynomial` to
`named_arrays.PolynomialFitFunctionArray.from_degree`, so the
least-squares sums ran over every axis of the calibration points,
including axes orthogonal to the scene such as the channel axis of a
multi-channel instrument. All channels were collapsed into a single
averaged polynomial that fit none of them. Passing the scene axes
restricts the fit sums to wavelength/field, giving an independent
polynomial per orthogonal element.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HuYXL1BkWsWGAxYnpD8Kqz
- Update the LinearSystem and SequentialSystem docstring examples for the
  pixel-based sensor grid: express the SimpleDistortionModel in pix instead of
  mm (so it matches the pixel `coordinates_sensor`) and plot the image
  directly, since the `integrate=True` default already collapses the
  wavelength axis (dropping the now-invalid `.sum("wavelength")`).
- Coerce `coordinates` to a plain SpectralPositionalVectorArray in
  `weights()`/`weights_transposed()` so richer scene vectors (e.g. Doppler
  scenes) are normalized to wavelength+position before distortion.
- Normalize `AbstractImagingSensor.pixels()` output to pix units.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
Codecov flagged 6 uncovered lines in optika/systems/_linear.py: the
`axis_field`/`axis_wavelength` inference branches of `image_from_weights`
and `backproject_from_weights`, which never ran because `image()`/
`backproject()` always pass those axes explicitly. Add
`test_from_weights_infer_axes`, which calls the `*_from_weights` methods
directly with the axes left to infer from the scene, bringing the module
to 100%.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
The forward model, `image`, accepts a scene in either photon or energy
units (the sensor detects and converts), so make the transpose symmetric: add
a `unit` parameter to `AbstractSystem.backproject`, `LinearSystem.backproject`,
and `backproject_from_weights` so the caller can request the backprojected
spectral radiance in whichever units match the scene.

The conversion (`_radiance_to_unit`) scales the natural photon radiance by the
energy per photon (hc/lambda) when an energy unit is requested, using the
wavelength already available in `backproject_from_weights`. The default
(`unit=None`) leaves the radiance in photon units, unchanged from before.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
Replace the hand-written `wavelength.volume_cell(...) * position.volume_cell(...)`
(with the manual cell-center alignment) in `image_from_weights` and
`backproject_from_weights` with the new `coordinates.volume_cell(...)` from
named-arrays 2.3. Bump the named-arrays lower bound to ~=2.3 accordingly.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
`image_from_weights` and `backproject_from_weights` computed the voxel
volume with `coordinates.volume_cell(...)`, which raised
`NotImplementedError` for a `DopplerPositionalVectorArray` scene (as
produced by the ctis `IdealInstrument`) since `volume_cell` is only
defined for the spectral-positional vector.

Coerce to a plain `SpectralPositionalVectorArray` just for the volume
computation, mirroring the ctis `_volume_scene` helper. The coercion is
local to the volume, so `image` is unchanged for a normal scene and
`backproject` returns the radiance on the original grid (a Doppler grid
stays a Doppler grid, preserving `wavelength_rest`).

Add `test_image_doppler` and `test_backproject_doppler` (with a
`_scene_doppler` helper) covering the Doppler path across both linear
system configurations.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CtdKmedevWkDab6BWupXqQ
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.

2 participants