Namespace lsst::shapelet

namespace shapelet

Typedefs

typedef afw::geom::ellipses::Quadrupole EllipseCore
typedef ndarray::Array<double, 1> Array1d

Typedef for a commonly-used array type.

Note

Serves as a workaround for ambiguities in the C++ standard itself: http://www.open-std.org/jtc1/sc22/wg21/docs/cwg_active.html#325

Enums

enum BasisTypeEnum

An enum that sets whether to use real-valued polar shapelets or Cartesian shapelets.

The conversion between the two bases is theoretically exact, but of course subject to round-off error here.

Values:

HERMITE

Cartesian shapelets or Gauss-Hermite functions, as defined in Refregier, 2003. That is, \( \psi(x, y)_{n_x, n_y} = \frac{H_{n_x}(x) H_{n_y}(y) e^{-\frac{x^2 + y^2}{2}}} {\sigma 2^{n_x + n_y} \sqrt{\pi n_x! n_y!}} \) where \(H_n(x)\) is a Hermite polynomial.

The ordering of coefficients [n_x, n_y] is (row-major packed): [0,0], [0,1], [1,0], [0,2], [1,1], [2,0], [0,3], [1,2], [2,1], [3,0], [0,4], [1,3], [2,2], [3,1], [4,0] etc.

LAGUERRE

Polar shapelets or Gauss-Laguerre functions, as defined in Bernstein and Jarvis, 2002. That is, \( \psi(x, y, \sigma)_{p, q} = (-1)^q \sqrt{\frac{q!}{p!}} (x + i y)^{p-q} e^{-\frac{x^2 + y^2}{2}} L^{(p-q)}_q(x^2 + y^2) \) where \(L^{(m)}_n(r)\) is an associated Laguerre polynomial.

The ordering of coefficients [p, q] is (row-major packed): [0,0], Re([1,0]), Im([1,0]), Re([2,0]), Im([2,0]), [1,1], Re([3,0]), Im([3,0], Re([2,1]), Im([2,1]), Re([4,0]), Im([4,0], Re([3,1]), Im([3,1]), [2,2] etc.

Elements with p < q are redundant in representing real-valued functions, while those with p == q are inherently real.

Functions

int computeOffset(int order)

Return the offset of the given order in a coefficient vector.

int computeSize(int order)

Return the size of the coefficient vector for the given order.

int computeOrder(int size)

Infer the order of a shapelet expansion from the number of coefficients.

Exceptions

  • InvalidParameterError: if the number of coefficients does not correspond to any shapelet order.

double intSqrt(int n)

Compute the square root of an integer number.

double rationalSqrt(int n, int d)

Compute the square root of a rational number i.e. sqrt(n/d)

Variables

double const BASIS_NORMALIZATION

Normalization factor for 1-d orthonormal shapelets: pi^(-1/4)

class BasisEvaluator
#include <BasisEvaluator.h>

Evaluates a standard shapelet Basis.

class ConversionMatrix
#include <ConversionMatrix.h>

Conversions between shapelet basis types.

The basis conversion matrix is block-diagonal and only needs to be computed once, so we cache the blocks in a hidden singleton and provide operations that act on shapelet matrices while taking advantage of the sparseness of the conversion.

class GaussHermiteConvolution
#include <GaussHermiteConvolution.h>

A parametrized matrix that performs a convolution in shapelet space.

GaussHermiteConvolution is defined only for the HERMITE basis type.

class GaussHermiteEvaluator
#include <GaussHermiteEvaluator.h>

A class to evaluate HERMITE shapelet-related quantities.

class HermiteTransformMatrix
#include <HermiteTransformMatrix.h>

A class that computes a matrix that applies a linear transform to a 2-d Hermite polynomial expansion.

Let

\[ Z_{\boldsymbol{n}}\!(\boldsymbol{x}) \equiv \mathcal{H}_{n_0}\!(x_0)\;\mathcal{H}_{n_1}\!(x_1) \]
where
\[ \mathcal{H}_n(x)=(2^n \pi^{1/2} n!)^{-1/2}H_n(x) \]
is the \(i\)th “alternate” Hermite polynomial. This function computes the matrix \(\boldsymbol{Q}(\boldsymbol{U})\) given a linear transform \(\boldsymbol{U}\) such that
\[ Z_{\boldsymbol{m}}\!(\boldsymbol{U}\boldsymbol{x}) = \sum_{\boldsymbol{n}} Q_{\boldsymbol{m},\boldsymbol{n}}\!(\boldsymbol{U})\,Z_{\boldsymbol{n}}\!(\boldsymbol{x}) \]

template<typename T>
class MatrixBuilder
#include <MatrixBuilder.h>

Class that evaluates a (multi-)shapelet basis at predefined points.

The output “matrix” has pixels along the rows, and basis elements along columns; this is the design matrix involved in a linear least squares fit for the basis coefficients.

A MatrixBuilder can be constructed in two ways: via one of its own constructors, or via a MatrixBuilderFactory. Using the latter allows the workspace arrays used by the MatrixBuilder to be shared between instances. See MatrixBuilderFactory and MatrixBuilderWorkspace for more information.

template<typename T>
class MatrixBuilderFactory
#include <MatrixBuilder.h>

A factory class for MatrixBuilder, providing more control over workspace memory.

To allocate workspace manually for a MatrixBuilder, we use the following pattern:; 

MatrixBuilderFactory<T> factory(...);
MatrixBuilderWorkspace<T> workspace(factory.computeWorkspace());
MatrixBuilder<T> builder = factory(workspace);
Because we aren’t doing anything special with the workspace, however, this is actually exactly equivalent to just constructing the MatrixBuilder directly, i.e.: 
MatrixBuilder<T> builder(...);

A more interesting case is if we want to use the same workspace for a pair of MatrixBuilders: 

MatrixBuilderFactory<T> factory1(...);
MatrixBuilderFactory<T> factory2(...);
MatrixBuilderWorkspace<T> workspace1(
   std::max(factory1.computeWorkspace(), factory2.computeWorkspace())
);
MatrixBuilderWorkspace<T> workspace2(workspace1);
MatrixBuilder<T> builder1 = factory1(workspace1);
MatrixBuilder<T> builder2 = factory2(workspace2);

template<typename T>
class MatrixBuilderWorkspace
#include <MatrixBuilder.h>

Reusable, shareable workspace for MatrixBuilder.

Multiple MatrixBuilders are often used together to evaluate a multi-shapelet model, and in this case it’s more efficient for the MatrixBuilders to use the same memory for temporary arrays rather than have them each allocate their own workspace. At other times, it may be useful to use one workspace for a sequence of throwaway MatrixBuilders, to avoid unnecessary memory allocations. This class manages the memory used for a MatrixBuilder’s workspace arrays, and provides methods for tracking it and sharing it between multple MatrixBuilders. See MatrixBuilderFactory for examples.

MatrixBuilderWorkspace holds a ndarray::Manager::Ptr that “owns” the block of memory. This is passed on to any MatrixBuilders constructed with it, ensuring that the block of memory remains alive with the MatrixBuilder even if the workspace object goes out of scope - so it is not necessary to keep the workspace object alive manually for the duration of its users.

In addition, MatrixBuilderWorkspace tracks both the current point in memory, and increments this as workspace matrices and vectors are created from it. It also checks that any created arrays do not exceed the bounds of the allocated memory. In order to share workspace memory between MatrixBuilders, the workspace object is simply copied - copied objects share the same memory, but maintain different “current” pointers, and hence create arrays from the same block of memory.

class MultiShapeletBasis
#include <MultiShapeletBasis.h>

A basis formed from a linear combination of shapelet bases that differ only in radius.

A MultiShapeletBasis can have many “components” (shapelet bases with different orders and radii), which are mapped via matrices into one or more “elements”. It’s common for a basis to have only one or two elements, representing a galaxy model that is a linear combination of one or two radial profiles. It’s also common for most components to be zeroth order (Gaussians), as higher- order shapelet terms don’t provide much of an advantage when approximating axisymmetric functions.

MultiShapeletBasis itself provides the interface to define these multi-Gaussian approximations and combine and refine them, and delegates the work of defining them to MultiShapeletFunction (via the makeFunction() method) and the MultiShapeletMatrixBuilder class. MultiShapeletFunction is a more user-friendly and versatile approach, intended for debugging and testing, while the MultiShapletMatrixBuilder approach is intended for performance-critical evaluation of PSF-convolved MultiShapeletBasis objects.

class MultiShapeletBasisComponent
#include <MultiShapeletBasis.h>

Simple struct that represents one shapelet expansion in a MultiShapeletBasis.

A MultiShapeletBasis is formed by the linear combination of several shapelet bases with different radii and common ellipticity; this represents a single shapelet basis within the MultiShapeletBasis.

Note

This really ought to be an inner class, and should generally be referred to via the MultiShapeletBasis::Component typedef, but Swig doesn’t handle inner classes.

class MultiShapeletFunction
#include <MultiShapeletFunction.h>

A multi-scale shapelet function.

class MultiShapeletFunctionEvaluator
#include <MultiShapeletFunction.h>

Evaluates a MultiShapeletFunction.

This is distinct from MultiShapeletFunction to allow the evaluator to construct temporaries and allocate workspace that will be reused when evaluating at multiple points.

A MultiShapeletFunctionEvaluator is invalidated whenever the MultiShapeletFunction it was constructed from is modified.

class MultiShapeletFunctionKey : public lsst::afw::table::FunctorKey<MultiShapeletFunction>
#include <FunctorKeys.h>

Class that maps MultiShapeletFunction objects to fields in afw::table objects.

A MultiShapeletFunctionKey holds a sequnece of ShapeletFunctionKey, with an numerical prefix in front of each component. A two-component MultiShapeletFunctionKey would thus be associated with the following keys:

  • ”<prefix>_0_xx”

  • ”<prefix>_0_yy”

  • ”<prefix>_0_xy”

  • ”<prefix>_0_x”

  • ”<prefix>_0_y”

  • ”<prefix>_1_xx”

  • ”<prefix>_1_yy”

  • ”<prefix>_1_xy”

  • ”<prefix>_1_x”

  • ”<prefix>_1_y”

As with all FunctorKeys, a MultiShapeletFunctorKey can be used to directly get or set objects on an afw::table::BaseRecord, just as with a true Key.

class PackedIndex
#include <GaussHermiteEvaluator.h>

An iterator-like object to help in traversing “packed” shapelet or Hermite polynomial matrix or vector dimensions.

A pair of indices (x,y) is mapped to the packed position i = (x+y)(x+y+1)/2 + x.

Typical usage is in a nested loop of the form: 

for (PackedIndex i; i.getOrder() <= order; ++i) {
    // utilize i
}

class RadialProfile
#include <RadialProfile.h>

Registry and utility class for multi-Gaussian approximations to radial profiles.

RadialProfile provides C++ access to the saved multi-Gaussian approximations stored in the pickle files in the data subdirectory of the shapelet package, by maintaining a singleton registry of these that is populated at import time by tractor.py. It also provides access to the true (exact) profiles, (mostly useful for educational/diagnostic purposes) and information about the relationship between different radii for this profile (see getMomentsRadiusFactor()).

Each RadialProfile object represents a particular “true” profile, and can hold multiple multi-Gaussian (or multi-Shapelet) approximations with different degrees of fidelity, as controlled by the number of components in the approximation and the maximum radius at which the approximation was fit. For more information about these approximations, see tractor.py and references therein. All RadialProfiles are defined in units of the half-light radius of the true profile, even for profiles for which this is not the radius typically used.

Several predefined subclasses of RadialProfile are defined privately, and can be accessed by name through the static get() method:

  • ”gaussian”: a single Gaussian profile

  • ”exp”: Exponential profile (Sersic n=1)

  • ”ser2”: Sersic profile with n=2

  • ”ser3”: Sersic profile with n=3

  • ”dev”: de Vaucouleur profile (Sersic n=4)

  • ”ser5”: Sersic profile with n=5

  • ”lux”: SDSS truncated Exponential profile

  • ”luv”: SDSS truncated and softened de Vaucouleur profile

class ShapeletFunction
#include <ShapeletFunction.h>

A 2-d function defined by an expansion onto a Gauss-Laguerre or Gauss-Hermite basis.

The coefficients are in units of flux, not surface brightness; increasing the area of the basis ellipse while leaving the coefficients unchanged will decrease the surface brightness by the ratio of the areas of the new and old ellipses. This convention is necessary to ensure that convolution with a delta function is sane. Because the BasisEvaluator class does not deal with ellipses, it is necessary to divide its output by R^2 to get the same result as a ShapeletFunctionEvaluator whose ellipse has determinant radius of R.

The coefficient normalization is not identical to that of a Gaussian, however; a zeroth-order ShapeletFunction with its only coefficient value set to 1 has a flux of 2.0 * pi^(1/2). This value is defined as ShapeletFunction::FLUX_FACTOR. Of course, to get the flux of a more complex shapelet expansion you have to use ShapeletFunctionEvaluator::integrate().

Note that the units of the coefficients would have to be radius for basis functions with the same ellipse to be orthonormal, but this orthonormality isn’t very useful, because the basis functions aren’t even orthogonal in the more common case that the ellipses differ. However, basis functions defined on the unit circle are still orthonormal.

class ShapeletFunctionEvaluator
#include <ShapeletFunction.h>

Evaluates a ShapeletFunction.

This is distinct from ShapeletFunction to allow the evaluator to construct temporaries and allocate workspace that will be reused when evaluating at multiple points.

A ShapeletFunctionEvaluator is invalidated whenever the ShapeletFunction it was constructed from is modified.

class ShapeletFunctionKey : public lsst::afw::table::FunctorKey<ShapeletFunction>
#include <FunctorKeys.h>

Class that maps ShapeletFunction objects to fields in afw::table objects.

A ShapeletFunctionKey manages a set of fields with a common prefix and the following suffixes:

  • ”x”, “y”, “xx”, “yy”, “xy”: ellipse

  • ”0”, “1”, “2”, …: coefficients.

As with all FunctorKeys, a ShapeletFunctorKey can be used to directly get or set objects on an afw::table::BaseRecord, just as with a true Key.