CModel.h

Go to the documentation of this file.

// -*- lsst-c++ -*-
/*
 * LSST Data Management System
 * Copyright 2008-2013 LSST Corporation.
 *
 * This product includes software developed by the
 * LSST Project (http://www.lsst.org/).
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the LSST License Statement and
 * the GNU General Public License along with this program.  If not,
 * see <http://www.lsstcorp.org/LegalNotices/>.
 */

#ifndef LSST_MEAS_MODELFIT_CModelFit_h_INCLUDED
#define LSST_MEAS_MODELFIT_CModelFit_h_INCLUDED

#include <bitset>
#include <string>

#include "ndarray.h"

#include "lsst/geom.h"
#include "lsst/pex/config.h"
#include "lsst/meas/base/exceptions.h"
#include "lsst/afw/table/Source.h"
#include "lsst/shapelet/RadialProfile.h"
#include "lsst/meas/modelfit/Model.h"
#include "lsst/meas/modelfit/Prior.h"
#include "lsst/meas/modelfit/MixturePrior.h"
#include "lsst/meas/modelfit/SoftenedLinearPrior.h"
#include "lsst/meas/modelfit/SemiEmpiricalPrior.h"
#include "lsst/meas/modelfit/UnitTransformedLikelihood.h"
#include "lsst/meas/modelfit/optimizer.h"
#include "lsst/meas/modelfit/PixelFitRegion.h"

namespace lsst { namespace meas { namespace modelfit {

struct CModelStageControl {

    CModelStageControl() :
        profileName("lux"),
        priorSource("EMPIRICAL"),
        priorName(),
        nComponents(8),
        maxRadius(0),
        usePixelWeights(false),
        weightsMultiplier(1.0),
        doRecordHistory(true),
        doRecordTime(true)
    {}

    shapelet::RadialProfile const & getProfile() const {
        return shapelet::RadialProfile::get(profileName);
    }

    PTR(Model) getModel() const;

    PTR(Prior) getPrior() const;

    LSST_CONTROL_FIELD(
        profileName, std::string,
        "Name of the shapelet.RadialProfile that defines the model to fit"
    );

    LSST_CONTROL_FIELD(
        priorSource, std::string,
        "One of 'FILE', 'LINEAR', 'EMPIRICAL', or 'NONE', indicating whether the prior should be loaded "
        "from disk, created from one of the nested prior config/control objects, or None"
    );

    LSST_CONTROL_FIELD(
        priorName, std::string,
        "Name of the Prior that defines the model to fit (a filename in $MEAS_MODELFIT_DIR/data, "
        "with no extension), if priorSource='FILE'.  Ignored for forced fitting."
    );

    LSST_NESTED_CONTROL_FIELD(
        linearPriorConfig, lsst.meas.modelfit.priors, SoftenedLinearPriorControl,
        "Configuration for a linear prior, used if priorSource='LINEAR'."
    );

    LSST_NESTED_CONTROL_FIELD(
        empiricalPriorConfig, lsst.meas.modelfit.priors, SemiEmpiricalPriorControl,
        "Configuration for an empirical prior, used if priorSource='EMPIRICAL'."
    );

    LSST_CONTROL_FIELD(nComponents, int, "Number of Gaussian used to approximate the profile");

    LSST_CONTROL_FIELD(
        maxRadius,
        int,
        "Maximum radius used in approximating profile with Gaussians (0=default for this profile)"
    );

    LSST_CONTROL_FIELD(
        usePixelWeights,
        bool,
        "Use per-pixel variances as weights in the nonlinear fit (the final linear fit for"
        " flux never uses per-pixel variances)"
    );

    LSST_CONTROL_FIELD(
        weightsMultiplier,
        double,
        "Scale the likelihood by this factor to artificially reweight it w.r.t. the prior."
    );

    LSST_NESTED_CONTROL_FIELD(
        optimizer, lsst.meas.modelfit.optimizer, OptimizerControl,
        "Configuration for how the objective surface is explored.  Ignored for forced fitting"
    );

    LSST_CONTROL_FIELD(
        doRecordHistory, bool,
        "Whether to record the steps the optimizer takes (or just the number, if running as a plugin)"
    );

    LSST_CONTROL_FIELD(
        doRecordTime, bool,
        "Whether to record the time spent in this stage"
    );

};

struct CModelControl {

    CModelControl() :
        psfName("modelfit_DoubleShapeletPsfApprox"),
        minInitialRadius(0.1),
        fallbackInitialMomentsPsfFactor(1.5)
    {
        initial.nComponents = 3; // use very rough model in initial fit
        initial.optimizer.gradientThreshold = 1E-3; // with slightly coarser convergence criteria
        initial.optimizer.minTrustRadiusThreshold = 1E-2;
        initial.usePixelWeights = true;
        dev.profileName = "luv";
        exp.nComponents = 6;
        exp.optimizer.maxOuterIterations = 250;
    }

    LSST_CONTROL_FIELD(
        psfName,
        std::string,
        "Field name prefix of the Shapelet PSF approximation used to convolve the galaxy model; "
        "must contain a set of fields matching the schema defined by shapelet.MultiShapeletFunctionKey."
    );

    LSST_NESTED_CONTROL_FIELD(
        region, lsst.meas.modelfit.pixelFitRegion, PixelFitRegionControl,
        "Configuration parameters related to the determination of the pixels to include in the fit."
    );

    LSST_NESTED_CONTROL_FIELD(
        initial, lsst.meas.modelfit.cmodel, CModelStageControl,
        "An initial fit (usually with a fast, approximate model) used to warm-start the exp and dev fits, "
        "convolved with only the zeroth-order terms in the multi-shapelet PSF approximation."
    );

    LSST_NESTED_CONTROL_FIELD(
        exp, lsst.meas.modelfit.cmodel, CModelStageControl,
        "Independent fit of the exponential component"
    );

    LSST_NESTED_CONTROL_FIELD(
        dev, lsst.meas.modelfit.cmodel, CModelStageControl,
        "Independent fit of the de Vaucouleur component"
    );

    LSST_CONTROL_FIELD(
        minInitialRadius, double,
        "Minimum initial radius in pixels (used to regularize initial moments-based PSF deconvolution)"
    );

    LSST_CONTROL_FIELD(
        fallbackInitialMomentsPsfFactor, double,
        "If the 2nd-moments shape used to initialize the fit failed, use the PSF moments multiplied by this."
        "  If <= 0.0, abort the fit early instead."
    );

};

struct CModelStageResult {

    enum FlagBit {
        FAILED=0,        
        TR_SMALL,        
        MAX_ITERATIONS,  
        NUMERIC_ERROR,   
        BAD_REFERENCE,   
        N_FLAGS          
    };

    CModelStageResult();

    PTR(Model) model;    
    PTR(Prior) prior;    
    PTR(OptimizerObjective) objfunc;  
    PTR(UnitTransformedLikelihood) likelihood; 
    Scalar instFlux;         
    Scalar instFluxErr;    
    Scalar instFluxInner;    
    Scalar objective;    
    Scalar time;         
    afw::geom::ellipses::Quadrupole ellipse;  

    ndarray::Array<Scalar const,1,1> nonlinear;  
    ndarray::Array<Scalar const,1,1> amplitudes; 
    ndarray::Array<Scalar const,1,1> fixed;      

    afw::table::BaseCatalog history;  
    std::bitset<N_FLAGS> flags; 
};

struct CModelResult {

    enum FlagBit {
        FAILED=0,                
        REGION_MAX_AREA,                  
        REGION_MAX_BAD_PIXEL_FRACTION,    
        REGION_USED_FOOTPRINT_AREA,       
        REGION_USED_PSF_AREA,             
        REGION_USED_INITIAL_ELLIPSE_MIN,  
        REGION_USED_INITIAL_ELLIPSE_MAX,  
        NO_SHAPE,                
        SMALL_SHAPE,             
        NO_SHAPELET_PSF,         
        BAD_CENTROID,            
        BAD_REFERENCE,           
        N_FLAGS                  
    };

    CModelResult();

    Scalar instFlux;       
    Scalar instFluxErr;  
    Scalar instFluxInner;  
    Scalar fracDev;    
    Scalar objective;  

    CModelStageResult initial; 
    CModelStageResult exp;     
    CModelStageResult dev;     

    afw::geom::ellipses::Quadrupole initialFitRegion;  
    afw::geom::ellipses::Quadrupole finalFitRegion;    

    LocalUnitTransform fitSysToMeasSys; 
    std::bitset<N_FLAGS> flags; 
};

class CModelAlgorithm {
public:

    typedef CModelControl Control; 
    typedef CModelResult Result;   

    CModelAlgorithm(
        std::string const & name,
        Control const & ctrl,
        afw::table::Schema & schema
    );

    CModelAlgorithm(
        std::string const & name,
        Control const & ctrl,
        afw::table::SchemaMapper & schemaMapper
    );

    explicit CModelAlgorithm(Control const & ctrl);

    Control const & getControl() const { return _ctrl; }

    Result apply(
        afw::image::Exposure<Pixel> const & exposure,
        shapelet::MultiShapeletFunction const & psf,
        geom::Point2D const & center,
        afw::geom::ellipses::Quadrupole const & moments,
        Scalar approxFlux=-1,
        Scalar kronRadius=-1,
        int footprintArea=-1
    ) const;

    Result applyForced(
        afw::image::Exposure<Pixel> const & exposure,
        shapelet::MultiShapeletFunction const & psf,
        geom::Point2D const & center,
        Result const & reference,
        Scalar approxFlux=-1
    ) const;

     void measure(
        afw::table::SourceRecord & measRecord,
        afw::image::Exposure<Pixel> const & exposure
    ) const;

    void measure(
        afw::table::SourceRecord & measRecord,
        afw::image::Exposure<Pixel> const & exposure,
        afw::table::SourceRecord const & refRecord
    ) const;

    void fail(
        afw::table::SourceRecord & measRecord,
        meas::base::MeasurementError * error
    ) const;

    void writeResultToRecord(Result const & result, afw::table::BaseRecord & record) const;

private:

    friend class CModelAlgorithmControl;

    // Actual implementations go here; we use an output argument for the result so we can get partial
    // results to the plugin version when we throw.
    void _applyImpl(
        Result & result,
        afw::image::Exposure<Pixel> const & exposure,
        shapelet::MultiShapeletFunction const & psf,
        geom::Point2D const & center,
        afw::geom::ellipses::Quadrupole const & moments,
        Scalar approxFlux,
        Scalar kronRadius=-1,
        int footprintArea=-1
    ) const;

    // Actual implementations go here; we use an output argument for the result so we can get partial
    // results to the SourceRecord version when we throw.
    void _applyForcedImpl(
        Result & result,
        afw::image::Exposure<Pixel> const & exposure,
        shapelet::MultiShapeletFunction const & psf,
        geom::Point2D const & center,
        Result const & reference,
        Scalar approxFlux
    ) const;

    // gets/checks inputs from SourceRecord that are needed by both apply and applyForced
    template <typename PixelT>
    shapelet::MultiShapeletFunction _processInputs(
        afw::table::SourceRecord & source,
        afw::image::Exposure<PixelT> const & exposure
    ) const;

    class Impl;

    Control _ctrl;
    PTR(Impl) _impl;
};

}}} // namespace lsst::meas::modelfit

#endif // !LSST_MEAS_MODELFIT_CModelFit_h_INCLUDED