pyimcom.psfutil

Utilities for PSFs and system matrices.

Classes

OutPSF: Simple target output PSF models.
PSFGrp: Group of PSFs attached to either an InStamp or a Block.
PSFOvl: Overlap between two PSFGrp instances or a PSFGrp instance and itself.
SysMatA: System matrix A attached to a coadd.Block instance.
SysMatB: System matrix B attached to a coadd.Block instance.

Attributes

_hasG4460

Classes

`PSFInterpolator`	Selector for PSF interpolations.
`OutPSF`	Simple target output PSF models (for testing or outputs).
`PSFGrp`	Group of PSFs.
`PSFOvl`	Overlap between two PSFGrp instances or a PSFGrp instance and itself.
`SysMatA`	System matrix A attached to a coadd.Block instance.
`SysMatB`	System matrix B attached to a coadd.Block instance.

Module Contents

_hasG4460 = True[source]

class PSFInterpolator[source]

Selector for PSF interpolations.

set_G4460()[source]: Set up G4460 as default interpolator (classmethod).

unset_G4460()[source]: Unset G4460 as default interpolator (classmethod).

gridC[source]

Grid interpolator, default = gridD5512C.

Type:: function

iC[source]

General interpolator, default = iD5512C.

Type:: function

iC_sym[source]

Symmetric interpolator, default = iD5512C_sym.

Type:: function

gridC[source]

iC[source]

iC_sym[source]

classmethod set_G4460()[source]: Set up G4460 as default interpolator (classmethod).

classmethod unset_G4460()[source]: Return to default interpolator (D5512)

class OutPSF[source]

Simple target output PSF models (for testing or outputs).

psf_gaussian()[source]: Gaussian spot (staticmethod).

psf_simple_airy()[source]: Airy spot (staticmethod).

psf_cplx_airy()[source]: Airy function with extra features printed on top to make it messier (staticmethod).

iD5512C_getw()[source]: No-Numba version of routine.iD5512C_getw (staticmethod).

psf_get_fwhm(): Compute full width at half maximum of a PSF (staticmethod).

psf_get_inv_width(): Compute shear-invariant width of a PSF (staticmethod).

static psf_gaussian(n: int, sigmax: float, sigmay: float) → numpy.array[source]

Gaussian spot, n x n, given sigma, centered (useful for testing).

Parameters:

n (int) – Size of output to generate.
sigmax (float) – 1 sigma width in the x-direction.
sigmay (float) – 1 sigma width in the y-direction.

Returns:

Output image, shape (n, n).

Return type:

np.array of flat

static psf_simple_airy(n: int, ldp: float, obsc: float = 0.0, tophat_conv: float = 0.0, sigma: float = 0.0) → numpy.array[source]

Airy spot, optionally obscured and convolved with a Gaussian and tophat.

Airy spot, n x n, with lambda/D = ldp pixels, and convolved with a tophat (square, full width tophat_conv) and Gaussian (sigma) and linear obscuration factor (obsc).

Parameters:

n (int) – Size of output to generate.
ldp (float) – Diffraction spot width lambda/D in pixels.
obsc (float, optional) – Fractional linear obscuration.
tophat_conv (float, optional) – Convolution with a square top-hat of the given full width.
sigma (float, optional) – Convolution with a Gaussian of the given 1 sigma width.

Returns:

The output spot. Shape (n, n).

Return type:

np.array

Notes

The result is centered on (n-1)/2, (n-1)/2 (so on a pixel if n is odd and a corner if n is even).

The image is normalized to sum to unity if analytically extended to infinity (so with a stamp enclosing 90% of the energy, will sum to 0.9).

static psf_cplx_airy(n: int, ldp: float, tophat_conv: float = 0.0, sigma: float = 0.0, features: int = 0) → numpy.array[source]

somewhat messier Airy function with a few diffraction features printed on ‘features’ is an integer that can be added. everything is band limited

Parameters:

n (int) – Size of output to generate.
ldp (float) – Diffraction spot width lambda/D in pixels.
tophat_conv (float, optional) – Convolution with a square top-hat of the given full width.
sigma (float, optional) – Convolution with a Gaussian of the given 1 sigma width.
features (int, optional) – Flags controlling which messy features to add.

Returns:

The output spot. Shape (n, n).

Return type:

np.array

See also

psf_simple_airy: Version without messy added features.

static iD5512C_getw(w: numpy.array, fh: float) → None[source]

No-Numba version of routine.iD5512C_getw.

The Numba version may not work in a Jupyter environment.

Parameters:

w (np.array, shape : (10,)) – Interpolation weights in one direction.
fh (float) – ‘xfh’ and ‘yfh’ with 1/2 subtracted.

Return type:

None.

static get_psf_fwhm(psf: numpy.array, visualize: bool = False) → float[source]

Compute FWHM of a PSF.

This function assumes that the given PSF is azimuthally symmetric.

Parameters:

psf (np.array, shape : (ny, nx)) – PSF of which the FWHM is to be computed. Typically something returned by PSFGrp._get_outpsf.
visualize (bool, optional) – Whether to visualize PSF and FWHM. The default is False.

Returns:

FWHM of given PSF in units of pixels.

Return type:

float

static get_psf_inv_width(psf: numpy.array) → float[source]

Compute shear-invariant width of a PSF.

Parameters:: psf (np.array, shape : (ny, nx)) – PSF of which the shear-invariant width is to be computed. Typically something returned by PSFGrp._get_outpsf.
Returns:: Shear-invariant width of given PSF in units of pixels.
Return type:: float

class PSFGrp(in_or_out: bool = True, inst=None, blk=None, verbose: bool = False, visualize: bool = False)[source]

Group of PSFs.

Either a group of input PSFs attached to an coadd.InStamp instance or a group of output PSFs attached to a coadd.Block instance.

Parameters:

in_or_out (bool, optional) – True if input PSF group, False if output PSF group.
inst (pyimcom.coadd.InStamp) – Input stamp, must be provided when in_or_out is True.
blk (pyimcom.coadd.Block) – Block, must be provided when in_or_out is False.
verbose (bool, optional) – Whether to print additional information.
visualize (bool, optional) – Whether to visualize PSFs and sampling positions.

setup()[source]: Set up class attributes (classmethod).

__init__()[source]: Constructor.

visualize_psf()[source]: Visualize a PSF (staticmethod).

_sample_psf()[source]: Sample a single PSF or a set of PSFs.

_build_inpsfgrp()[source]: Build a group of input PSFs.

_get_outpsf()[source]: Get an output PSF specified by configuration (staticmethod).

_build_outpsfgrp()[source]: Build a group of output PSFs.

accel_pad_and_rfft2()[source]: Zero-padding and 2d rfft (staticmethod).

clear()[source]: Free up memory space.

nsamp = 383[source]

nc = 191[source]

nfft = 768[source]

classmethod setup(npixpsf: int = 48, oversamp: int = 8, dtheta: float = 0.025 / 3600, psfsplit: bool = False) → None[source]

Set up class attributes.

Parameters:

npixpsf (int, optional) – Size of PSF postage stamp in native pixels.
oversamp (int, optional) – PSF oversampling factor relative to native pixel scale.
dtheta (float, optional) – Output pixel scale in degrees. The default is 0.025/3600 (corresponds to 0.025 arcsec)).
psfsplit (bool, optional) – Is PSF splitting implemented in this run?

Return type:

None

in_or_out = True[source]

psf_rft[source]

static visualize_psf(psf_: numpy.array, yxco: numpy.array, xctr: float, yctr: float) → None[source]

Visualize a single PSF, together with sampling positions.

Parameters:

psf (np.array) – A single PSF. Must be 2-dimensional.
yxco (np.array) – Sampling positions, with (0, 0) at the PSF center. Must be 3-dimensional. `yxco`[0,:,:] represents x-positions and `yxco`[1,:,:] represents y-positions.
xctr (float) – Position of the PSF center.
yctr (float) – Position of the PSF center.

Return type:

None

_sample_psf(idx, psf, outpix2world2inpix=None, visualize=False)[source]

Perform interpolations to sample a single PSF or a set of PSFs.

Parameters:

idx (int or None) – Index of the single PSF (int). If None, sample a set of PSFs at the same time.
psf (np.array) – PSF(s) to be sampled. Shape (n_psf, ny, nx).
outpix2world2inpix (method, optional) – InImage outpix2world2inpix of the appropriate InImage instance. The default is None, meaning not to rotate the sampling positions.
visualize (bool, optional) – Whether to visualize PSFs and sampling positions.

Return type:

None

_build_inpsfgrp(visualize: bool = False) → None[source]

Build a group of input PSFs.

Parameters:: visualize (bool, optional) – Whether to visualize PSFs and sampling positions.
Return type:: None

static _get_outpsf(outpsf: str = 'AIRYOBSC', extrasmooth: float = 0.0, use_filter: int = 4) → numpy.array[source]

Get an output PSF specified by configuration.

Parameters:

outpsf (str, optional) – Target output PSF type. Options are ‘GAUSSIAN’, ‘AIRYOBSC’, and ‘AIRYUNOBSC’.
extrasmooth (float, optional) – Target output PSF extra smearing. The default is 0.0.
use_filter (int, optional) – Which filter to use (0..10).

Returns:

np.array, shape – Output PSF specified by configuration.

Return type:

(PSFGrp.nsamp+1, PSFGrp.nsamp+1)

_build_outpsfgrp(visualize: bool = False) → None[source]

Build a group of output PSFs.

What output PSFs to include is set by the configuration file.

Parameters:: visualize (bool, optional) – Whether to visualize PSFs and sampling positions.
Return type:: None

static accel_pad_and_rfft2(psf_arr: numpy.array) → numpy.array[source]

Accelerated version of zero-padding and rfft2.

For nsamp = 537 and nfft = 1280, this saves about 20% of the time compared to simply using rfft2.

Parameters:: psf_arr (np.array) – PSF array to be padded and Fourier transformed. The shape should be (…, PSFGrp.nsamp, PSFGrp.nsamp).
Returns:: Real Fourier transform results of psf_arr. The shape is (…, PSFGrp.nfft, PSFGrp.nfft//2+1).
Return type:: np.array

Notes

Diagram of the steps involved here:

#                                       +---+      +---+
#                                       |   | fft  |***|
#  +---+ pad  +-------+ rfft +---+ pad  |   | ===> |***|
#  |***| ===> |***    | ===> |***| ===> |***|      |***|
#  +---+      +-------+      +---+      +---+      +---+
#  psf_arr     pad_m1                   pad_m2      res

Note: m1 stands for axis=-1, and m2 stands for axis=-2.

clear(verbose: bool = False) → None[source]

Free up memory space.

Use this in addition to the default __del__ method to make sure that the PSF array is removed when it is no longer used.

Parameters:: verbose (bool, optional) – Whether to print additional information.
Return type:: None

class PSFOvl(psfgrp1: PSFGrp, psfgrp2: PSFGrp = None, verbose: bool = False, visualize: bool = False)[source]

Overlap between two PSFGrp instances or a PSFGrp instance and itself.

Parameters:

psfgrp1 (pyimcom.psfutil.PSFGrp) – The first PSFGrp instance.
psfgrp2 (pyimcom.psfutil.PSFGrp or None, optional) – The second PSFGrp instance. The default is None, indicating the self-overlap of psfgrp1.
verbose (bool, optional) – Whether to print additional information.
visualize (bool, optional) – Whether to visualize the PSF overlap array.

setup()[source]: Set up class attribute (classmethod).

__init__()[source]: Constructor.

_idx_square2triangle()[source]: Convert a 2D square index to a 1D triangle index.

accel_irfft2_and_extract()[source]: irfft2 and extraction (staticmethod).

_build_psfovl()[source]: Build the PSF overlap array.

visualize_psfovl()[source]: Visualize the PSF overlap array.

__call__()[source]: Wrapper for the C interpolators.

_call_ii_cross()[source]: Interpolations in the input-input cross-overlap case.

_call_io_cross()[source]: Interpolations in the input-output cross-overlap case.

_call_ii_self()[source]: Interpolations in the input-input self-overlap case.

clear()[source]: Free up memory space.

flat_penalty = 1e-07[source]

classmethod setup(flat_penalty: float = 1e-07) → None[source]

Set up class attribute.

This must be run after PSFGrp.setup.

Parameters:: flat_penalty (float, optional) – Amount by which to penalize having different contributions to the output from different input images.
Return type:: None

grp1[source]

grp2 = None[source]

_idx_square2triangle(idx1: int, idx2: int) → int[source]

Convert a 2D square index to a 1D triangle index.

Parameters:

idx1 (int, int) – 2D square index.
idx2 (int, int) – 2D square index.

Returns:

1D triangle index.

Return type:

int

Notes

The motivation is that, in the input self-overlap case, the overlap between psf_arr[idx2] and psf_arr[idx1] is simply the flipped version of the overlap between psf_arr[idx1] and psf_arr[idx2]. Therefore, we only need to compute the latter, and the results are stored in a np.array of shape (n_psf*(n_psf+1)//2, nsamp, nsamp).

For example, in the n_psf == 3 case:

#  \idx2  0   1   2
# idx1  +---+---+---+
#   0   | 0 | 1 | 2 |
#       +---+---+---+
#   1   |   | 3 | 4 |
#       +---+---+---+
#   2   |   |   | 5 |
#       +---+---+---+

static accel_irfft2_and_extract(ovl_rft: numpy.array) → numpy.array[source]

Accelerated version of irfft2 and extraction (used in FFT-based convolution).

Parameters:: ovl_rft (np.array) – Real Fourier transform of the PSF overlap array we want. The shape is (…, PSFGrp.nfft, PSFGrp.nfft//2+1).
Returns:: The PSF overlap array we want. The shape is (…, nsamp, nsamp)
Return type:: np.array

Notes

This acceleration is based on the fact that we only need less than (1/2) x (1/2) of the irfft2 results. For nsamp = 537 and nfft = 1280, this saves about 40% of the time compared to simply using irfft2 and ifftshift.

Because of the ifftshift part, this function is not the inverse function of PSFGrp.accel_pad_and_rfft2.

Diagram of the steps involved here (ext. stands for extraction):

# +---+      +---+
# |***| ifft |***|
# |***| ===> |***| ext. +---+ irfft+-------+ ext. +---+
# |***|      |***| ===> |***| ===> |*******| ===> |***|
# +---+      +---+      +---+      +-------+      +---+
# ovl_rft    ift_m2     ovl_m2      ift_m1        ovl_m1

Note: m2 stands for axis=-2, and m1 stands for axis=-1.

If PSFOvl.nsamp is at least half of PSFGrp.nfft, then reverts to using irfft2 and ifftshift.

_build_psfovl(visualize: bool = False) → None[source]

Build the PSF overlap array.

Parameters:: visualize (bool, optional) – Whether to visualize the PSF overlap array.
Return type:: None

visualize_psfovl() → None[source]: Visualize the PSF overlap array.

__call__(st1, st2=None, visualize=False) → numpy.array[source]

Wrapper for the C interpolators.

Parameters:

st1 (pyimcom.coadd.InStamp) – First input stamp.
st2 (pyimcom.coadd.InStamp or coadd.OutStamp or None) – Second input stamp. The default is None, indicating diagonal blocks of the A matrix (2nd stamp same as 1st).
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.

Returns:

The shape depends on the nature of this PSFOvl instance and the input.

Return type:

np.array

_call_ii_cross(st1, st2, visualize=False) → numpy.array[source]

Interpolations in the input-input cross-overlap case.

Parameters:

st1 (pyimcom.coadd.InStamp) – First input stamp.
st2 (pyimcom.coadd.InStamp) – Second input stamp.
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.

Returns:

A 2D array of the cross-correlation between the two input stamps. The shape is (st1.pix_cumsum[-1], st2.pix_cumsum[-1]).

Return type:

np.array

_call_io_cross(st1, st2, visualize=False)[source]

Interpolations in the input-output cross-overlap case.

Parameters:

st1 (pyimcom.coadd.InStamp) – Input stamp.
st2 (pyimcom.coadd.OutStamp) – Output stamp.
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.

Returns:

A 3D image of the cross-correlation of the input and output PSFs, with axis=0 indicating which output PSF is used. The shape is either (self.grp2.n_psf, n_outpix, st1.pix_cumsum[-1]) or (self.grp2.n_psf, n_outpix, selection.shape[0]).

Return type:

np.array

_call_ii_self(st1, st2, visualize=False)[source]

Interpolations in the input-input self-overlap case.

Note that the self-overlap of an input PSF group can be used to compute either the block-diagonal submatrices for a single input stamp or the submatrices for a pair of input stamps within a 2x2 group.

Parameters:

st1 (pyimcom.coadd.InStamp) – First input stamp.
st2 (pyimcom.coadd.InStamp or None, optional) – Second input stamp. The default is None, indicating diagonal blocks of the A matrix.
visualize (bool, optional) – Whether to visualize overlap arrays and sampling positions.

Returns:

A 2D image of the cross-correlation. The shape is (st1.pix_cumsum[-1], st2.pix_cumsum[-1]).

Return type:

np.array

clear(verbose: bool = False) → None[source]

Free up memory space.

Use this in addition to the default __del__ method to make sure that the huge PSF overlap array is removed when it is no longer used.

Parameters:: verbose (bool, optional) – Whether to print additional information.
Return type:: None

class SysMatA(blk)[source]

System matrix A attached to a coadd.Block instance.

The symtem matrix A is defined in Rowe+ 2011 Equation (17).

Parameters:: blk (pyimcom.coadd.Block) – The Block instance to which this SysMatA instance is attached.

__init__()[source]: Constructor.

ji_st2psf()[source]: Convert InStamp index to PSFGrp index (staticmethod).

shift_ji_st()[source]: Tool function for tuple addition (staticmethod).

iisubmat_dist()[source]: Calculate the index for iisubmats_ref (staticmethod).

_compute_iisubmats()[source]: Make input-input PSFOvl and compute A submatrices.

get_iisubmat()[source]: Return a requested A submatrix.

clear()[source]: Free up memory space.

blk[source]

iisubmats[source]

iisubmats_ref[source]

static ji_st2psf(ji_st: int, int)[source]

Convert InStamp index to PSFGrp index.

More precisely, convert index of the InStamp of interest to that of its neighbor to which the PSFGrp in that 2x2 group is attached. For example, any element of {(0, 0), (0, 1), (1, 0), (1, 1)} -> (0, 0).

Parameters:: ji_st ((int, int)) – Index of the InStamp of interest.
Returns:: Index of the InStamp which harbors the PSFGrp.
Return type:: (int, int)

static shift_ji_st(ji_st: int, int, dji_st: int, int)[source]

Tool function for tuple addition.

For our purposes, this function shifts an InStamp index in each direction by either 0 or 1.

Parameters:

ji_st ((int, int)) – InStamp index to be shifted.
dji_st ((int, int)) – Shift in InStamp index.

Returns:

Shifted InStamp index.

Return type:

(int, int)

static iisubmat_dist(ji_st1: int, int, ji_st2: int, int)[source]

Calculate the index for iisubmats_ref.

Specifically, calculate the “distance” between two input postage stamps in order to store the reference counts in the 3D array iisubmats_ref using indices (`ji_st1`[0], `ji_st1`[1], dist), which this static method returns.

Parameters:

ji_st1 ((int, int)) – Index of the first InStamp.
ji_st2 ((int, int)) – Index of the second InStamp.

Returns:

Index of the first InStamp combined with the “distance,” which can be directly used as the index of iisubmats_ref. Returns None when the “distance” is out of range.

Return type:

(int, int, int) or None

Notes

The “distance” is defined as follows:

# +---+---+---+---+---+
# | 8 | 9 | 10| 11| 12|
# +---+---+---+---+---+
# | 3 | 4 | 5 | 6 | 7 |
# +---+---+---+---+---+
# |   |   | 0*| 1 | 2 |
# +---+---+---+---+---+

where “*” denotes ji_st1. Note that ji_st1 must precede ji_st2. If the “distance” is out of range, this function returns None.

_compute_iisubmats(ji_st1: int, int, ji_st2: int, int, sim_mode: bool = False, verbose: bool = False, visualize: bool = False) → None[source]

Make input-input PSFOvl and compute A submatrices.

This method is only called in SysMatA.get_iisubmat. When an A submatrix is requested from an OutStamp but not in iisubmats, get_iisubmat calls this method, which makes the input-input PSFOvl, computes all the needed A submatrices which are supposed to be computed using that particular PSFOvl instance, and clears the PSFOvl.

Note that removing the huge PSF overlap arrays is important! Otherwise pyimcom would demand much more memory.

To determine whether a specific system submatrix is needed, a Block instance first runs the postage coaddition framework in sim_mode to count references to all possible system submatrices.

Parameters:

ji_st1 ((int, int)) – Index of the first InStamp.
ji_st2 ((int, int)) – Index of the second InStamp.
sim_mode (bool, optional) – Whether to count references without actually computing submatrices.
verbose (bool, optional) – Whether to print additional information.
visualize (bool, optional) – Whether to make more plots.

Return type:

None

get_iisubmat(ji_st1: int, int, ji_st2: int, int, sim_mode: bool = False, ji_st_out: int, int = None, visualize: bool = False) → numpy.array[source]

Return the requested A submatrix.

This is the only public interface of this class. For the sim_mode, see the docstring of SysMatA._compute_iisubmats.

Parameters:

ji_st1 ((int, int)) – Index of the first InStamp.
ji_st2 ((int, int)) – Index of the second InStamp.
sim_mode (bool, optional) – Whether to count references without actually computing submatrices.
ji_st_out ((int, int)) – Index of the OutStamp. Needed for virtual memory.
visualize (bool, optional) – Whether to make plots.

Returns:

The requested A submatrix. The shape is (st1.pix_cumsum[-1], st2.pix_cumsum[-1]).

Return type:

np.array

clear() → None[source]: Free up memory space.

class SysMatB(blk)[source]

System matrix B attached to a coadd.Block instance.

The symtem matrix B is defined in Rowe+ 2011 Equation (18). IMPORTANT: The -2 coefficient is NOT included in this program.

Parameters:: blk (pyimcom.coadd.Block) – The Block instance to which this SysMatB instance is attached.

__init__()[source]: Constructor.

get_iosubmat()[source]: Return a requested B submatrix.

clear()[source]: Free up memory space.

blk[source]

iopsfovls[source]

iopsfovls_ref[source]

get_iosubmat(ji_st_in: int, int, ji_st_out: int, int, sim_mode: bool = False, visualize: bool = False) → numpy.array[source]

Return a requested B submatrix.

Note that this is the courterpart of the combination of SysMatA.get_iisubmat and the _compute_iisubmats method behind it, since SysMatB is much simpler than SysMatA.

Parameters:

ji_st_in ((int, int)) – Index of the InStamp.
ji_st_out ((int, int)) – Index of the OutStamp.
sim_mode (bool, optional) – Whether to count references without actually computing submatrices.
visualize (bool, optional) – Whether to visualize the in-out overlap.

Returns:

np.array, shape – The requested B submatrix.

Return type:

same as PSFOvl._call_io_cross