Skip to content

Fix the effective-area unit stripped by the linear-system weights - #194

Merged
roytsmart merged 1 commit into
feature/interpolated-systemfrom
fix/linear-weights-area-unit
Jul 31, 2026
Merged

Fix the effective-area unit stripped by the linear-system weights#194
roytsmart merged 1 commit into
feature/interpolated-systemfrom
fix/linear-weights-area-unit

Conversation

@jacobdparker

Copy link
Copy Markdown
Contributor

Summary

AbstractLinearSystem.weights and weights_transposed strip the unit of the effective area with .value before handing the weights to na.regridding, while image_from_weights/backproject_from_weights restore that unit as the hardcoded weights_unit (u.cm**2). If the effective-area model happens to evaluate in any other unit, every forward image is silently wrong by the unit ratio — the ESIS-II instrument model evaluates in mm², making LinearSystem.image exactly 100× too bright relative to the raytraced SequentialSystem.image.

This converts with area.to_value(self.weights_unit) instead, so the stripped number is always in the unit that gets restored.

Test

Adds test_image_area_unit_invariance to AbstractTestAbstractLinearSystem: imaging the same scene must give identical results when the effective-area model is re-expressed in mm². Fails at 100× on the base branch, passes with the fix.

After the fix, LinearSystem.image and SequentialSystem.image agree to within the degree-2 vignetting-fit residual (~20%) on an ESIS-II test scene, instead of two orders of magnitude.

🤖 Generated with Claude Code

`weights` and `weights_transposed` stripped the unit of the effective
area with `.value`, while `image_from_weights` and
`backproject_from_weights` restore it as the hardcoded `weights_unit`
(cm^2). If the effective area model evaluates in any other unit (mm^2
for the ESIS-II model), every image is wrong by the unit ratio (x100).

Convert with `.to_value(self.weights_unit)` instead, and add a
regression test checking that the image is invariant under a unit
change of the effective area model.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@roytsmart
roytsmart merged commit 532e288 into feature/interpolated-system Jul 31, 2026
4 of 10 checks passed
@roytsmart
roytsmart deleted the fix/linear-weights-area-unit branch July 31, 2026 18:50
roytsmart added a commit that referenced this pull request Aug 2, 2026
…on to `SequentialSystem` (#152)

* Added `optika.systems.InterpolatedSystem` as an approximation to `SequentialSystem`.

* lots of improvements

* lots of fixes

* add tests

* Document `LinearSystem` with a runnable end-to-end example

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

* Add `SequentialSystem.linearize()` to build a `LinearSystem`

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

* fixes

* coverage

* Add `LinearSystem.backproject` and accept photon or energy radiance

`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

* Invert `expose` in `LinearSystem.backproject`

`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

* Fold volume and sensor response into `image_from_weights`/`backproject_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

* Add an `uncertainty` flag to `LinearSystem.image`

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

* revert to old ruff rules

* Apply the noise model per readout: add `integrate` to the sensor

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

* Reconcile `image`/`backproject` signatures across systems

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

* Express the distortion and sensor grid in pixel coordinates

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

* Fit the polynomial distortion and vignetting models per channel

`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

* Fix docs build and input handling for the pixel-coordinate LinearSystem

- 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

* Cover axis inference in image_from_weights/backproject_from_weights

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

* Let `backproject` express the radiance in photon or energy units

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

* Use `SpectralPositionalVectorArray.volume_cell` for the voxel volume

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

* Accept a Doppler object-plane grid in image and backproject

`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

* Use `.spectral_positional` to accept Doppler grids

Replace the hand-written `na.SpectralPositionalVectorArray(wavelength=...,
position=...)` coercions in `LinearSystem` (`weights`, `weights_transposed`,
`image_from_weights`, `backproject_from_weights`) with the
`.spectral_positional` property added in named-arrays 2.4.0, and bump the
pin accordingly. Behavior is unchanged; the property centralizes the
spectral-/Doppler-positional normalization.

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

* Convert the effective area to `weights_unit` before stripping its unit (#194)

`weights` and `weights_transposed` stripped the unit of the effective
area with `.value`, while `image_from_weights` and
`backproject_from_weights` restore it as the hardcoded `weights_unit`
(cm^2). If the effective area model evaluates in any other unit (mm^2
for the ESIS-II model), every image is wrong by the unit ratio (x100).

Convert with `.to_value(self.weights_unit)` instead, and add a
regression test checking that the image is invariant under a unit
change of the effective area model.

Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

* Raise `NotImplementedError` for `RayVectorArray.volume_cell`

named-arrays 2.4.0 added `volume_cell` to
`AbstractSpectralPositionalVectorArray`, which `RayVectorArray`
subclasses. A ray bundle is a scattered collection of rays rather than a
logically-rectangular grid, so the per-voxel volume is undefined and the
inherited implementation crashes on scalar-wavelength ray fixtures.
Override it to raise `NotImplementedError`, matching the existing
`type_matrix`/`test_matrix` precedent and the test's expectation.

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

* Address code-review findings on the linear-system PR

- backproject_from_weights: re-attach the full target wavelength grid
  before `photons_absorbed`, so the integrated readout is spread across
  the correct number of wavelength bins. The detector readout only
  carries the two band edges, so `num_wavelength` was always 1 and the
  per-bin spreading (and per-bin quantum efficiency) never happened.
- AbstractSystem.backproject: make it a concrete default that raises
  `NotImplementedError` instead of an `@abc.abstractmethod`, since
  backprojection is an optional capability (a `SequentialSystem` has no
  transpose). This stops it from breaking unrelated `AbstractSystem`
  subclasses that never backproject.
- image/backproject and the *_from_weights variants: make the optional
  arguments keyword-only, so a positional `image(scene, pupil)` can no
  longer silently bind `pupil` to `axis_wavelength` after the parameter
  reorder.
- _collapse_wavelength: use min/max rather than the first/last edge, so
  a non-ascending wavelength grid still reports sorted band edges.

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

---------

Co-authored-by: Claude Opus 4.8 <noreply@anthropic.com>
Co-authored-by: jacobdparker <jacobdparker@gmail.com>
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