LSSTApplications  17.0+124,17.0+14,17.0+73,18.0.0+37,18.0.0+80,18.0.0-4-g68ffd23+4,18.1.0-1-g0001055+12,18.1.0-1-g03d53ef+5,18.1.0-1-g1349e88+55,18.1.0-1-g2505f39+44,18.1.0-1-g5315e5e+4,18.1.0-1-g5e4b7ea+14,18.1.0-1-g7e8fceb+4,18.1.0-1-g85f8cd4+48,18.1.0-1-g8ff0b9f+4,18.1.0-1-ga2c679d+1,18.1.0-1-gd55f500+35,18.1.0-10-gb58edde+2,18.1.0-11-g0997b02+4,18.1.0-13-gfe4edf0b+12,18.1.0-14-g259bd21+21,18.1.0-19-gdb69f3f+2,18.1.0-2-g5f9922c+24,18.1.0-2-gd3b74e5+11,18.1.0-2-gfbf3545+32,18.1.0-26-g728bddb4+5,18.1.0-27-g6ff7ca9+2,18.1.0-3-g52aa583+25,18.1.0-3-g8ea57af+9,18.1.0-3-gb69f684+42,18.1.0-3-gfcaddf3+6,18.1.0-32-gd8786685a,18.1.0-4-gf3f9b77+6,18.1.0-5-g1dd662b+2,18.1.0-5-g6dbcb01+41,18.1.0-6-gae77429+3,18.1.0-7-g9d75d83+9,18.1.0-7-gae09a6d+30,18.1.0-9-gc381ef5+4,w.2019.45
LSSTDataManagementBasePackage
lsst::meas::modelfit; Source Measurement via Model Fitting

Introduction

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.

Models, Interpreters, and Parameters

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:

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:

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:

, Objective, and Prior

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.

Units and Coordinate Systems

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.

Command-line Tasks

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:

  1. We read inputs using the Butler.
  2. We "prep" the entire output catalog, creating empty records containing starting parameters and the Footprint that defines the pixels to fit. We can do this either by bootstrapping from a SourceCatalog produced by running one of the ProcessImageTasks on the data, or by using the outputs of a previous run of one of the meas_modelfit tasks as inputs. Unlike most tasks, multiple meas_modelfit Task runs with different configurations may coexist in the same data repository, as long as each run has a unique value for the 'tag' configuration option. These tags are included in the data ID for the output products, so we can chain multiple runs with different tags, or simply compare independent runs without having to create multiple output repositories and use multiple Butlers in our analysis code.
  3. We iterate over the output catalog, and fit each record. We first create a Likelihood object, then pass it to the 'fitter' subtask, which creates the appropriate Objective object, fits it, and sets the fields of the output record.
  4. We write the output catalog using the Butler.

Samplers and Optimizers

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.

Future Plans and Known Issues