CorrelatedObservation#

class scarlet2.CorrelatedObservation(data, psf=None, wcs=None, channels=None, renderer=None, name='', correlation_function=None, mask=None, native_scale=None, n_eff=None)[source]#

Bases: Observation

Content and definition of an observation with pixel correlations

The noise model is still assumed to be Gaussian, but with correlations between pixels. The implementation computes the goodness of fit in Fourier space from the noise power spectrum to avoid the expensive computation of/with an inverse banded matrix in configuration space.

Create an observation with a correlated-noise Gaussian likelihood.

The noise model is given by the measured correlation_function; the power spectrum used by the likelihood is derived from it. In most cases it is easier to use from_observation() (for a coadd already on the desired grid) or from_resampling() (to resample the data here), or stack_observations() to combine several correlated observations that share a grid.

Parameters:
  • data (array) – Observed data, 2D or 3D (C, H, W).

  • psf – As for Observation.

  • wcs – As for Observation.

  • channels – As for Observation.

  • renderer – As for Observation.

  • name – As for Observation.

  • correlation_function (dict) – 2D pixel correlation function keyed by integer (dy, dx) lag, as returned by scarlet2.measure.correlation_function(). Required.

  • mask (array, optional) – Boolean array marking invalid pixels. Defaults to all-valid.

  • native_scale (astropy.units.Quantity, optional) – Angular size of a native, pre-resampling pixel, e.g. 0.17 * u.arcsec. Set this if the data sit on a grid finer than their true resolution (resampled by an upstream pipeline or by from_resampling()). Only 1 / (native_scale / pixel_scale) ** 2 of the pixels are then independent, and chi^2 is restricted to the spatial frequencies the native sampling supports. Requires wcs. Mutually exclusive with n_eff.

  • n_eff (sequence of int, optional) – Per-channel effective degrees of freedom, one entry per channel. Overrides the count derived from mask and native_scale. Used by stack_observations() to carry through the (possibly heterogeneous) counts of the stacked observations.

property N#

Effective number of degrees of freedom of the noise model

For an observation without resampling this is the number of unmasked pixels, as in Observation.N. When the correlations come from resampling onto a finer grid, the noise covariance is rank-deficient: there are only as many independent noise values as unmasked pixels in the original observation. Using the pixel count there would dilute goodness_of_fit() and the likelihood normalization by the surplus modes, which carry no information.

This is the sum of the per-channel n_eff.

__call__()#

Evaluate the model

check_set_renderer(frame)#

Check existence of renderer, or set it by calling match()

Parameters:

frame (Frame) – The frame to match

Return type:

None

correlation_function: dict#

Measured data correlation function, from py:meth:scarlet2.measure.correlation_function

data: Array#

Observed data

eval_chi_square_in_box_and_border(scene, border_width=3)#

Evaluate the weighted mean (weighted by the inverse variance weights) of the squared residuals for each source. Chi square is also computed for the perimeter outside the box of with border_width.

Parameters:
  • scene (Scene) – Scene containing the sources

  • border_width (int) – width of the border around the source box

Return type:

Dict of sources indices and their corresponding Dict of residuals inside and outside source box.

frame: Frame#

Metadata to describe what view of the sky data amounts to

classmethod from_observation(obs, patch_size=50, maxlength=12, native_scale=None)[source]#

Create a CorrelatedObservation from data that already sit on the desired grid.

Use this for a coadd (drizzled, warped, stacked, …) whose pixel grid matches the model frame. The method finds a patch of size patch_size with as few sources as possible, measures the 2D pixel correlation function in it (out to maxlength), and turns that into the noise power spectrum used by the likelihood. Truncating the correlation function at maxlength biases the power spectrum, severely so once the correlation length approaches maxlength.

To resample the data onto the model grid here instead, use from_resampling().

Parameters:
  • obs (Observation) – Observation containing the data and weight map. The weights are used only to locate bad and bright pixels for the correlation measurement.

  • patch_size (int) – Linear size (in pixels) of the source-free patch used to measure the correlation function. Clamped to the image size; if it reaches the image size the whole image is used.

  • maxlength (int) – Maximum lag (in pixels) of the 2D correlation function. Large enough to cover the extent of the correlations, and small compared to patch_size so every lag is averaged over many pixel pairs (the assertion enforces 4 * maxlength <= patch_size).

  • native_scale (astropy.units.Quantity, optional) – Angular size of a native, pre-resampling pixel, e.g. 0.17 * u.arcsec. Set this only if an upstream pipeline already resampled the data onto a grid finer than their true resolution: the pixels are then correlated and not all independent, so chi^2 is restricted to the spatial frequencies the native sampling supports. See CorrelatedObservation.__init__().

Return type:

CorrelatedObservation

classmethod from_resampling(obs, to_frame, lanczos_order=9, resample_psf=True, patch_size=50, maxlength=12)[source]#

Create a CorrelatedObservation by resampling obs onto the grid of to_frame.

The Lanczos resampling correlates the pixels, so the result needs a correlated-noise likelihood. This method measures the correlation function of the resampled data.

If the spatial grid of to_frame is finer than obs, the resampled noise is rank-deficient: chi^2 will automatically be restricted to the frequencies the native (obs) sampling supports.

Parameters:
  • obs (Observation) – Observation to resample. Its pixel scale is taken as the native resolution.

  • to_frame (Frame) – Frame describing the desired spatial sampling, typically the model frame.

  • lanczos_order (int) – Order of the Lanczos resampling kernel.

  • resample_psf (bool, optional) – Whether to resample obs.psf to to_frame. Set to False only if the PSF is already sampled at the resolution of to_frame.

  • patch_size (int) – Passed to from_observation().

  • maxlength (int) – Passed to from_observation().

Return type:

CorrelatedObservation

get(name=None)#

Get parameter(s) from this module

Parameters:

name (str, optional) – Name of parameter. If not set, returns all parameters.

Returns:

requested data arrays for parameters

Return type:

dict

goodness_of_fit(model)#

Evaluate the goodness of the model fit to the data

For a Gaussian noise model, the gof is defined as the averaged squared deviation of the model from the data, scaled by the variance of the data, aka mean chi squared \(\frac{1}{N}\sum_i=1^N w_i (m_i - d_i)^2\) with inverse variance weights \(w_i\).

Up to a normalization, the gof is identical to log-likelihood.

Parameters:

model (array) – The (pre-rendered) predicted data cube, typically from evaluating Scene

Return type:

float

log_likelihood(model)#

The logarithm the likelihood of data given model

Parameters:

model (array) – The (pre-rendered) predicted data cube, typically from evaluating Scene

Return type:

float

mask: Array#

Mask for invalid pixels

match(frame, renderer=None)[source]#

Construct the mapping between frame (from the model) and this observation frame

See Observation.match(). In addition, this method determines whether chi^2 can be evaluated on the padded grid of the renderer.

n_eff: tuple#

Per-channel effective number of degrees of freedom of the noise model, see N

One entry per channel. For a channel without resampling this is its number of unmasked pixels; for a channel resampled onto a finer grid it is the smaller count of independent noise values on the native grid. N is the sum over channels.

name: str#

Name to describe the observation

property parameters#

Parameters defined for this module

Returns:

name: (node, param) mapping for all parameters

Return type:

dict

render(model, **kwargs)#

Render model in the frame of this observation

Parameters:
  • model (array) – The (pre-rendered) predicted data cube, typically from evaluating Scene

  • kwargs (dict) – Additional keyword arguments to pass to the renderer

Returns:

Prediction of the observation given the model. Has the same shape as data.

Return type:

array

renderer: Renderer#

Renderer to translate from the model frame the observation frame

replace(name, value)#

Replace member attribuge name with value

Parameters:
  • name (str) – Name of member to replace

  • value (any) – Value to replace member with

Returns:

The modified module.

Return type:

Module

replace_subtree(where, value, parameters=None)#

Replace the subtree where(self) with value and carry over all other parameters

Unlike replace(), value does not need to have the same structure as the subtree it replaces, e.g. an array morphology can be replaced by a SersicMorphology. Parameters defined for this module remain attached to the same nodes, except those inside the replaced subtree, which are dropped.

Parameters:
  • where (callable) – Function module -> node that selects the subtree to replace, as in equinox.tree_at()

  • value (any) – Replacement for the subtree. If this module has a frame, attributes of value given in sky coordinates are converted to pixel coordinates of that frame.

  • parameters (callable, optional) – Function module -> None that is called with the new module inside of a Parameters context to define parameters for the replacement, e.g. by calling Parameter on nodes of value.

Returns:

The modified module.

Return type:

Module

Examples

>>> sersic = SersicMorphology(n=jnp.array(2.), size=jnp.array(2.), shape=morph.shape)
>>> def sersic_parameters(morph):
>>>     Parameter(morph.n, name="n:0", prior=dist.Uniform(0.5, 6))
>>>     Parameter(morph.size, name="size:0", constraint=constraints.positive, stepsize=0.1)
>>>
>>> scene_ = scene.replace_subtree(lambda sc: sc.sources[0].morphology, sersic, sersic_parameters)

See also

remap()

set(values)#

Set parameter(s) from this module with values

Parameters:

values (dict[str,jnp.array]) – values to replace parameters with, identified by their name

Returns:

new module with parameter(s) replaced by values

Return type:

Module

weights: Array#

Statistical weights (usually inverse variance) for log_likelihood()