Namespace lsst::meas::modelfit¶

namespace modelfit

Unnamed Group

typedef float Pixel: Typedefs to be used for pixel values

Unnamed Group

typedef double Scalar: Typedefs to be used for probability and parameter values

typedef Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> Matrix

typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> Vector

typedef afw::table::Key<Scalar> ScalarKey

typedef afw::table::Key<afw::table::Array<Scalar>> ArrayKey

Typedefs

typedef std::vector<PTR(Model)> ModelVector

Functions

void solveTrustRegion(ndarray::Array<Scalar, 1, 1> const &x, ndarray::Array<Scalar const, 2, 1> const &F, ndarray::Array<Scalar const, 1, 1> const &g, double r, double tolerance)

Solve a symmetric quadratic matrix equation with a ball constraint.

This computes a near-exact solution to the “trust region subproblem” necessary in trust-region-based nonlinear optimizers:

\[ \min_x{\quad g^T x + \frac{1}{2}x^T F x}\quad\quad\quad \text{s.t.} ||x|| \le r \]

The tolerance parameter sets how close to \(r\) we require the norm of the solution to be when it lies on the constraint, as a fraction of \(r\) itself.

This implementation is based on the algorithm described in Section 4.3 of “Nonlinear Optimization” by Nocedal and Wright.

class AdaptiveImportanceSampler : public lsst::meas::modelfit::Sampler

#include <AdaptiveImportanceSampler.h>

Sampler class that performs Monte Carlo sampling, while iteratively updating the analytic distribution from which points are drawn.

Between the iterations defined in the control object, the prior is applied to the samples, and the mixture distribution is updated using expectation-maximization to match the samples.

class CModelAlgorithm

#include <CModel.h>

Main public interface class for CModel algorithm.

See CModel Magnitudes for a full description of the algorithm.

This class provides the methods that actually execute the algorithm, and (depending on how it is constructed) holds the Key objects necessary to use SourceRecords for input and output.

struct CModelControl: #include <CModel.h>
The main control object for CModel, containing parameters for the final linear fit and aggregating the other control objects.

struct CModelResult: #include <CModel.h>
Master result object for CModel, containing results for the final linear fit and three nested CModelStageResult objects for the results of the previous stages.

struct CModelStageControl: #include <CModel.h>
Nested control object for CModel that configures one of the three (“initial”, “exp”, “dev”) nonlinear fitting stages.

struct CModelStageResult: #include <CModel.h>
Result object for a single nonlinear fitting stage of the CModel algorithm

class DoubleShapeletPsfApproxAlgorithm : public lsst::meas::base::SimpleAlgorithm

#include <DoubleShapeletPsfApprox.h>

An algorithm that fits a 2-component shapelet approximation to the PSF model.

This algorithm fits a similar model to the one that has been used on the HSC side for several data releases, but uses the improved optimizer in meas_modelfit instead of the one in (the now defunct) meas_extensions_multiShapelet and using moments to reduce the number of parameters in the fit. It is faster and more robust than GeneralShapeletPsfApprox, but much less flexible.

The model is not a fully general one, even for two shapelet components. We tie the ellipticities of the two components together, hold the center fixed, and keep their radii fixed at a value set in the configuration; this means we fit only one set of ellipse parameters, instead of two.

class DoubleShapeletPsfApproxControl: #include <DoubleShapeletPsfApprox.h>
Control object used to configure a 2-shapelet fit to a PSF model; see DoubleShapeletPsfApproxAlgorithm.

class EpochFootprint

#include <UnitTransformedLikelihood.h>

An image at one epoch of a galaxy, plus associated info

Includes one image of a galaxy and and associated footprint and multi-shapelet PSF model

class GeneralPsfFitter

#include <GeneralPsfFitter.h>

Class for fitting multishapelet models to PSF images.

This class fits up to four shapelet expansions simultaneously to a PSF image, with the relative radii and number of shapelet coefficients for each expansion separately configurable. These expansions are also named; this allows us to map different fits with some expansions disabled to each other, in order to first fit an approximate model and follow this up with a more complete model, using the approximate model as a starting point.

The configuration also defines a simple Bayesian prior for the fit, defined using simple independent Gaussians for the ellipse parameters of each component. The priors can be disabled by setting their width (xxPriorSigma in the control object) to infinity, and those parameters can be held fixed at their input values by setting the prior width to zero. The priors are always centered at the input value, meaning that it may be more appropriate to think of the priors as a form of regularization, rather than a rigorous prior. In fact, it’s impossible to use a prior here rigorously without a noise model for the PSF image, which is something the LSST Psf class doesn’t provide, and here is just provided as a constant noise sigma to be provided by the user (who generally just has to chose a small number arbitrarily). Decreasing the noise sigma will of course decrease the effect of the priors (and vice versa). In any case, having some sort of regularization is probably a good idea, as this is a very high-dimensional fit.

Subclassed by lsst::meas::modelfit::GeneralPsfFitterAlgorithm

class GeneralPsfFitterComponentControl: #include <GeneralPsfFitter.h>
Control object used to define one piece of multishapelet fit to a PSF model; see GeneralPsfFitterControl

class GeneralPsfFitterControl

#include <GeneralPsfFitter.h>

Control object used to configure a multishapelet fit to a PSF model; see GeneralPsfFitter.

The default configuration corresponds to fitting an elliptical double-Gaussian, in which each component can have different radii, positions, and ellipticities. While the fitter can support much more complex models, at present, fitting these is prohibitively slow, and is not recommended in production environments (use DoubleShapeletPsfApprox instead).

class ImportanceSamplerControl

#include <AdaptiveImportanceSampler.h>

Control object for one iteration of adaptive importance sampling.

See: AdaptiveImportanceSampler, AdaptiveImportanceSamplerTask

class Likelihood

#include <Likelihood.h>

Base class for optimizer/sampler likelihood functions that compute likelihood at a point.

Likelihood abstracts the problem of computing the likelihood over different kinds of data. It is responsible for creating a “model matrix” that maps amplitudes to data values, and maintaining a vector of scaled, weighted data values that corresponds to it. Its components can be represented best in the mathematical formula for a -log likelihood assuming Gaussian data and a model with both nonlinear parameters \(\theta\) and linear (“amplitude”) parameters \(\alpha\):

\[ L(\alpha,\theta) = \frac{1}{2}\left(y - A(\theta)\alpha\right)^T\, \Sigma^{-1}\,\left(y - A(\theta)\alpha\right) \]

where \(y\) is the data vector, \(\Sigma\) is the data covariance matrix (assumed to be diagonal), and \(A(\theta)\) is the “true” model matrix (parametrized on the nonlinear parameters).

When fitting or sampling from the likelihood, however, we don’t want to use these quantities directly, and they aren’t what the Likelihood class provides. Instead, we reparametrize with:

\[ w_i \equiv \Sigma_{i,i}^{-1/2} \]

\[ z_i = w_i y_i \]

\[ B_{i,j} = w_i A_{i,j} \]

resulting in the equivalent formula:

\[ L(\alpha,\theta) = \frac{1}{2}\left(z-B(\theta)\alpha\right)^T\,\left(z-B(\theta)\alpha\right) \]

The \(w_i\) are the weights, which are applied to both the data vector and the model matrix to account for the noise in the data. In some cases, we may choose to use a constant weight rather than per-pixel weights, but will will still use a vector to represent it.

Subclassed by lsst::meas::modelfit::MultiShapeletPsfLikelihood, lsst::meas::modelfit::UnitTransformedLikelihood

struct LocalUnitTransform

#include <UnitSystem.h>

A local mapping between two UnitSystems.

LocalUnitTransform is “local” because it linearizes the Wcs and evaluates the PhotoCalib transform at a particular predifined point, allowing it to represent the geometric transform as an AffineTransform and the photometric transform as a simple scaling.

class MixtureComponent: #include <Mixture.h>
A weighted Student’s T or Gaussian distribution used as a component in a Mixture.

class MixturePrior : public lsst::meas::modelfit::Prior: #include <MixturePrior.h>
A prior that’s flat in amplitude parameters, and uses a Mixture for nonlinear parameters.

class MixtureUpdateRestriction

#include <Mixture.h>

Helper class used to define restrictions to the form of the component parameters in Mixture::updateEM.

The base class implementation does not apply any restrictions.

class Model

#include <Model.h>

Abstract base class and concrete factories that define multi-shapelet galaxy models.

A Model provides a mapping from its parameters to ellipses and to a realization using shapelet objects. A Model does not “hold” its parameters; parameters are always stored in separate arrays.

Model parameters are split into three categories: nonlinear, amplitudes, and fixed. These are described more fully in modelfitParameters.

A few private concrete subclasses of Model have been provided that will meet most needs; instances can be constructed via the make() and makeGaussian()

Subclassed by lsst::meas::modelfit::MultiModel

class MultiModel : public lsst::meas::modelfit::Model: #include <MultiModel.h>
A concrete Model class that simply concatenates several other Models.

class MultiShapeletPsfLikelihood : public lsst::meas::modelfit::Likelihood: #include <GeneralPsfFitter.h>
Likelihood object used to fit multishapelet models to PSF model images; mostly for internal use by GeneralPsfFitter.

class Optimizer

#include <optimizer.h>

A numerical optimizer customized for least-squares problems with Bayesian priors.

The algorithm used by Optimizer combines the Gauss-Newton approach of approximating the second-derivative (Hessian) matrix as the inner product of the Jacobian of the residuals, while maintaining a matrix of corrections to this to account for large residuals, which is updated using a symmetric rank-1 (SR1) secant formula. We assume the prior has analytic first and second derivatives, but use numerical derivatives to compute the Jacobian of the residuals at every step. A trust region approach is used to ensure global convergence.

We consider the function \(f(x)\) we wish to optimize to have two terms, which correspond to negative log likelihood ( \(\chi^2/2=\|r(x)|^2\), where \(r(x)\) is the vector of residuals at \(x\)) and negative log prior \(q(x)=-\ln P(x)\):

\[ f(x) = \frac{1}{2}\|r(x)\|^2 + q(x) \]

At each iteration \(k\), we expand \(f(x)\) in a Taylor series in \(s=x_{k+1}-x_{k}\):

\[ f(x) \approx m(s) = f(x_k) + g_k^T s + \frac{1}{2}s^T H_k s \]

where

\[ g_k \equiv \left.\frac{\partial f}{\partial x}\right|_{x_k} = J_k^T r_k + \nabla q_k;\quad\quad J_k \equiv \left.\frac{\partial r}{\partial x}\right|_{x_k} \]

\[ H_k = J_k^T J_k + \nabla^2 q_k + B_k \]

Here, \(B_k\) is the SR1 approximation term to the second derivative term:

\[ B_k \approx \sum_i \frac{\partial^2 r^{(i)}_k}{\partial x^2}r^{(i)}_k \]

which we initialize to zero and then update with the following formula:

\[ B_{k+1} = B_{k} + \frac{v v^T}{v^T s};\quad\quad v\equiv J^T_{k+1} r_{k+1} - J^T_k r_k \]

Unlike the more common rank-2 BFGS update formula, SR1 updates are not guaranteed to produce a positive definite Hessian. This can result in more accurate approximations of the Hessian (and hence more accurate covariance matrices), but it rules out line-search methods and the simple dog-leg approach to the trust region problem. As a result, we should require fewer steps to converge, but spend more time computing each step; this is ideal when we expect the time spent in function evaluation to dominate the time per step anyway.

class OptimizerControl

#include <optimizer.h>

Configuration object for Optimizer.

Many of these configuration options pertain to how the trust region is updated. It’s easiest to understand these together rather than separately. At each iteration, a quadratic model of the objective function is formed. We can use this model to predict how we expect the objective function to behave over a step, and compare it to how the actual objective function behaves. To do this, we’ll use the ratio of the actual reduction in the objective function to the predicted reduction in the objective function, and call this \(\rho\). Then,

the step is accepted, and the parameters updated, when \(\rho >\) stepAcceptThreshold.
if \(\rho > \) trustRegionGrowReductionRatio and the length of the step is greater than trustRegionGrowStepFraction times the current trust region radius, the trust region radius will be multiplied by trustRegionGrowFactor.
if trustRegionShrinkMinReductionRatio \(< \rho < \) trustRegionShrinkMaxReductionRatio, the trust region radius will be multiplied by trustRegionShrinkFactor.

class OptimizerObjective: #include <optimizer.h>
Base class for objective functions for Optimizer.

class Prior

#include <Prior.h>

Base class for Bayesian priors.

Subclassed by lsst::meas::modelfit::MixturePrior, lsst::meas::modelfit::SemiEmpiricalPrior, lsst::meas::modelfit::SoftenedLinearPrior

class Sampler: Subclassed by lsst::meas::modelfit::AdaptiveImportanceSampler

class SemiEmpiricalPrior : public lsst::meas::modelfit::Prior: #include <SemiEmpiricalPrior.h>
A piecewise prior motivated by both real distributions and practical considerations.

class SoftenedLinearPrior : public lsst::meas::modelfit::Prior: #include <SoftenedLinearPrior.h>
A prior that’s linear in radius and flat in ellipticity, with a cubic roll-off at the edges.

class TruncatedGaussian

#include <TruncatedGaussian.h>

Represents a multidimensional Gaussian function truncated at zero.

This is typically used to represent the posterior probability of amplitude parameters, given a flat prior; we require that the amplitudes each be positive but otherwise do not modify the Gaussian likelihood.

Currently only 1 and 2 dimensions are supported, and all dimensions must be truncated. Computing integrals is the only operation for which > 2-d is not implemented, but the integrals must be computed upon construction, so we can’t support any other operations for > 2-d either.

Many operations on TruncatedGaussians are defined in -log space, as underflow/overflow problems will often occur in the non-log forms.

See modelfitTruncatedGaussianMath for implementation notes

class TruncatedGaussianEvaluator: #include <TruncatedGaussian.h>
Helper class for evaluating the -log of a TruncatedGaussian.

class TruncatedGaussianLogEvaluator: #include <TruncatedGaussian.h>
Helper class for evaluating the -log of a TruncatedGaussian.

class TruncatedGaussianSampler: #include <TruncatedGaussian.h>
Helper class for drawing samples from a TruncatedGaussian.

struct UnitSystem: #include <UnitSystem.h>
A simple struct that combines a Wcs and a PhotoCalib.

class UnitTransformedLikelihood : public lsst::meas::modelfit::Likelihood

#include <UnitTransformedLikelihood.h>

A concrete Likelihood class that does not require its parameters and data to be in the same UnitSystem.

This is the main concrete Likelihood class using when fitting sources, even in the case where the measurement UnitSystem is the same as that of the data; we always prefer to fit in a special UnitSystem (see modelfitUnits), using this class to transform the model realization when comparing to the data. This makes forced photometry and modelfit measurements just as easy as single-frame measurements (aside from data access); one can simply initialize a UnitTransformedLikelihood with multiple exposures instead of a single exposure to fit simultaneously to multiple exposures.

class UnitTransformedLikelihoodControl

#include <UnitTransformedLikelihood.h>

Control object used to initialize a UnitTransformedLikelihood.

Translated to Python as UnitTransformedLikelihoodConfig; the Swig-wrapped C++ Control object can be created from the config object via the makeControl() method (see lsst.pex.config.wrap).

namespace detail

Functions

Eigen::Vector4d solveRampPoly(double v0, double v1, double x0, double x1, double s0, double s1): Solve for the coefficients of a cubic polynomial p(x) that goes from p(x0)=v0 to p(x1)=v1, with p’(x0)=s0 and p’(x1)=s1.

double phid(double z)

Compute univariate normal probabilities.

This function computes

\[ \frac{1}{2\pi}\int_z^{\infty}dx\;e^{-x^2/(2z^2)} \]

This is just a simple wrapper around boost::math::erfc, used to provide the particular form expected by bvnu.

double bvnu(double h, double k, double rho)

Compute bivariate normal probabilities.

This function computes

\[ \frac{1}{2\pi\sqrt{1-\rho^2}}\int_h^{\infty}dx\int_k^{\infty}dy \;e^{-(x^2 - 2\rho x y + y^2)/(2(1-\rho^2))} \]

It is a reimplementation of the “bvnu” MatLab routine by Alan Genz. The original implementation can be found at http://www.math.wsu.edu/faculty/genz/homepage. It is based on the algorithm described in:

Drezner & Wesolowsky (1989), “On the computation of the bivariate normal integral”, Journal of Statist. Comput. Simul. 35, pp. 101-107.

A copy of Genz’s FORTRAN routine of the same is included in the tests directory, as it has been used to generate the test reference data there. It should generally not need to be compiled by users.

Most of the time, this function should be called via the methods in the TruncatedGaussian class; it is exposed publically only for testing purposes.

template<int N> class Vandermonde: #include <polynomials.h>
Class that computes rows of the Vandermonde matrix and related matrices; the dot product of these row vectors with the polynomial coefficient vectors evaluates the polynomial (or computes a derivative).

Navigation

Namespace lsst::meas::modelfit¶