LSST Applications g063fba187b+66a50001ff,g0f08755f38+1a22dc2551,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g217e2c1bcf+12e87a5bd8,g246886dfd9+466c7b6c06,g28da252d5a+858b171e04,g2bbee38e9b+404b60ec9b,g2bc492864f+404b60ec9b,g3156d2b45e+6e55a43351,g347aa1857d+404b60ec9b,g35bb328faa+a8ce1bb630,g3a166c0a6a+404b60ec9b,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+70bea58702,g599934f4f4+b8c5400ca5,g781aacb6e4+a8ce1bb630,g7af13505b9+b5b9cefdb8,g80478fca09+c2997882f3,g82479be7b0+8974e6af0f,g858d7b2824+1a22dc2551,g89c8672015+f4add4ffd5,g8f1c07a47a+de51c9b0a5,g9125e01d80+a8ce1bb630,ga5288a1d22+b66f8cf76b,gb58c049af0+d64f4d3760,gc28159a63d+404b60ec9b,gcab2d0539d+66cf1de5d4,gcf0d15dbbd+12cb7e2563,gda6a2b7d83+12cb7e2563,gdaeeff99f8+1711a396fd,ge79ae78c31+404b60ec9b,gef2f8181fd+414189b318,gf0baf85859+c1f95f4921,gf0c06eb49c+1a22dc2551,gfa517265be+1a22dc2551,gfa999e8aa5+17cd334064,v28.0.0.rc2
LSST Data Management Base Package
|
The meas_modelfit package contains code for fitting PSF-convolved models of astronomical objects to data. This includes both traditional fitting via "greedy" derivative-based optimizers and Monte Carlo exploration using adaptive importance sampling.
The galaxy modeling code can be roughly divided into three levels:
The high- and medium-level code in meas_modelfit is discussed in the next few sections, with links there to the API documentation for individual classes. For more information on the low-level model evaluation code, see the lsst::shapelet package, and in particular the lsst::shapelet::MultiShapeletBasis and lsst::shapelet::MultiShapeletMatrixBuilder classes.
The Model and Interpreter classes define a number of different types of parameters, and understanding the relationships between these parameters is crucial in understanding how the package works.
The Model class represents a particular parametrized form and defines how it is mapped to shapelet expansions for convolution and evaluation. An instance of Model does not itself contain parameters: we have one Model instance for a particular fitting configuration, not one Model instance per fitted object. Model splits its parameters up into three separate arrays:
'nonlinear'
are the nonlinear parameters of the model that will be fit for. These usually define one or more ellipses. In the mathematical notation used throughout the meas_modelfit documentation, the symbol \(\theta\) is used for the nonlinear parameters.'amplitudes'
are linear parameters of the model, in the sense that an evaluation of the model on a pixel grid can be represented by a matrix that is multiplied by the amplitude vector (see Likelihood). The symbol \(\alpha\) is used to represent the amplitudes in mathematical notation.'fixed'
are nonlinear parameters that are not fit for, but are rather held fixed at their initial values. It's particularly common to fix the position parameters of the model. Note that we can't just subsume the fixed parameters into the definition of the Model instance itself, as the fixed parameters may change from source to source but the Model instance is defined purely by the configuration.Like Model, an Interpreter instance is defined only by configuration, but it contains information not only about the model (via a Model object attribute), but about the fitter and Prior probabilities used. A single Interpreter instance is attached to the ModelFitCatalog used to store fitting results, and it contains methods that can be used to interpret the records of that catalog. We have different Interpreter classes for Monte Carlo sampling (SamplingInterpreter and its subclasses) and greedy optimization (OptimizerInterpreter), as well as different classes to account for different ways the fitters can handle the amplitude parameters: because the model is linear in the amplitudes, and we assume Gaussian noise on the data, we know the likelihood is Gaussian in the amplitude parameters. This can be used in different ways by different fitters to reduce the dimensionality of the fit. In order to do this without specializing the fitters themselves, we introduce one more set of parameters, called simply 'parameters'
. These are some combination of the amplitude and nonlinear parameters defined by the model, and represent the actual space explored by the fitter. There are generally two options here:
'parameters'
vector is simply the concatenation of the 'nonlinear'
vector with the 'amplitudes'
vector; we explore the full space at once, and do not take advantage of the special nature of the amplitudes. This is the approach used by the DirectSamplingInterpreter and the OptimizerInterpreter.'parameters'
vector is just the 'nonlinear'
vector. At each point in the nonlinear parameter space we analytically marginalize over or otherwise project out the amplitude dimensions and store information necessary to reconstruct them with each sample. This approach is currently only used by the MarginalSamplingInterpreter, but a variant for greedy optimization may be added in the future.By using a polymorphic Interpreter class whose instances are shared by all output records, we can avoid having polymorphic output objects that must be subclassed for every kind of fitter and objective. We use a single ModelFitRecord class to store the results of any fitter. ModelFitRecord contains a number of normal, Schema-defined fields, as well three "blob" objects:
The Likelihood class represents the comparison between data and the evaluation of a Model in the coordinate and photometric system of that data. The data can correspond to a single image or multiple images to be fit simultaneously. Likelihood also handles the inverse-variance weighting of both model and data, so the two can be compared directly. It does not actually compute the likelihood or log likelihood values, however - it simply stores the weighted data and computes a pixelized model that can be compared with it.
It is the responsibility of the Objective classes to take the Likelihood and actually compare the data and model in a way that is useful for a particular fitter. We have different Objective class hierarchies with different APIs for sampling (SamplingObjective) and optimization (OptimizerObjective), reflecting the fact that we need to compute slightly different quantities in these cases. Each Objective class has a corresponding Interpreter class that can be used to interpret the outputs of that particular Objective.
In addition to the Likelihood, the Objective classes also include a Bayesian prior, represented by the Prior class. There is currenly only one concrete Prior class, MixturePrior, which uses the Mixture class to represent the nonlinear parameters along with a flat prior for the amplitudes (though the amplitudes are required to be nonnegative). Using a flat prior for the amplitudes means that the posterior probability as well as the likelihood will be Gaussian, but the nonnegative constraint means it will be truncated. As a result, MixturePrior offloads many of its more complex numerical calculations to the TruncatedGaussian class.
All parameters in meas_modelfit are defined in coordinate and photometric systems that are unique to each object (or perhaps a small group of neighboring objects).
The coordinate system is a local TAN projection centered on the object, with arcsecond-size pixels; this allows the model to be defined in a coordinate system that can be trivially converted to celestial coordinates with no changes in ellipticity, while keeping radius parameters around unity. It also allows us to directly evaluate the prior probability of a set of parameters without having to transform the prior to different coordinate systems.
We also define a custom photometric system for each source, using a crudely-estimated nominal magnitude (provided as an input to the fitting code; usually this will be a simple aperture magnitude). We set the photometric system such that this nominal magnitude corresponds to unit flux, ensuring that amplitude parameters also remain of order unity. This helps avoid numerical problems, and (along with the choice of arcseconds for radius units and a sensible definition of ellipticity) allows us to avoid rescaling parameters within our fitters when using algorithms that must define sensible norms over the parameter space.
See the UnitSystem and ProjectedLikelihood classes for more information on the unit systems and transforms.
There are three concrete command-line Tasks in meas_modelfit:
The first two run on CCD images and coadd patch images, respectively, and inherit from MeasureImageTask, which provides much of the implementation. (this is analogous to lsst.pipe.tasks.ProcessImageTask and its subclasses).
MeasureMultiTask fits to multiple CCD images simultaneously, using the output from MeasureCoaddTask as an input (the CCD images included in the fit correspond to those that overlap a particular coadd patch).
All of these tasks ultimately inherit from BaseMeasureTask, which constains most of the implementation: subclasses of BaseMeasureTask only have to define how to create load inputs, write outputs, and create a Likelihood object.
BaseMeasureTask contains two registries that act as factories for Model and Prior objects (see models.py and priors.py) and a "fitter" subclass that must be one of AdaptiveImportanceSamplerTask or OptimizerTask. Determination of the pixel region to fit is handled by the code in fitRegion.py, and the shapelet PSF approximation needed to convolve the galaxy models is currently delegated to the meas_extensions_multiShapelet package.
All of these tasks have four distinct stages:
We currently provide two very different kinds of fitters:
In both cases, we provide a Task that may be used as the "fitter" subtask of BaseMeasureTask. These tasks initialize the sample schema and create an Interpreter upon construction, which will be used by BaseMeasureTask when creating the output catalogs. They also provide two methods (initialize and adaptPrevious) that update fitter-specific quantities during the catalog-prep stage. Their run() methods do the actual fitting, using a Likelihood object and an output ModelFitRecord provided by BaseMeasureTask.
src
catalog as an input, so the Tasks can be run on arbitrary data, not just the S13 sims.