TruncatedGaussian.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_TruncatedGaussian_h_INCLUDED
#define LSST_MEAS_MODELFIT_TruncatedGaussian_h_INCLUDED
 
#include "Eigen/Core"
#include "ndarray.h"
 
#include "lsst/base.h"
#include "lsst/afw/math/Random.h"
#include "lsst/meas/modelfit/common.h"
 
// TODO: we should really integrate this with Mixture somehow
 
namespace lsst { namespace meas { namespace modelfit {
 
class TruncatedGaussianSampler;
class TruncatedGaussianEvaluator;
class TruncatedGaussianLogEvaluator;
 
class TruncatedGaussian {
public:
 
    enum SampleStrategy {
 
        DIRECT_WITH_REJECTION, 
        ALIGN_AND_WEIGHT       
    };
 
    typedef TruncatedGaussianSampler Sampler;
    typedef TruncatedGaussianLogEvaluator LogEvaluator;
    typedef TruncatedGaussianEvaluator Evaluator;
 
    static TruncatedGaussian fromSeriesParameters(Scalar q0, Vector const & gradient, Matrix const & hessian);
 
    static TruncatedGaussian fromStandardParameters(Vector const & mean, Matrix const & covariance);
 
    Sampler sample(SampleStrategy strategy) const;
 
    Sampler sample(Scalar minRejectionEfficiency=0.1) const;
 
    LogEvaluator evaluateLog() const;
 
    Evaluator evaluate() const;
 
    int getDim() const;
 
    Vector maximize() const;
 
    Scalar getUntruncatedFraction() const;
 
    Scalar getLogPeakAmplitude() const;
 
    Scalar getLogIntegral() const;
 
    ~TruncatedGaussian();
 
private:
 
    friend class TruncatedGaussianSampler;
    friend class TruncatedGaussianLogEvaluator;
 
    class Impl;
 
    explicit TruncatedGaussian(std::shared_ptr<Impl> impl) : _impl(impl) {}
 
    std::shared_ptr<Impl> _impl;
};
 
class TruncatedGaussianLogEvaluator {
public:
 
    explicit TruncatedGaussianLogEvaluator(TruncatedGaussian const & parent);
 
    template <typename Derived>
    Scalar operator()(Eigen::MatrixBase<Derived> const & alpha) const {
        if ((alpha.array() < 0.0).any()) return std::numeric_limits<Scalar>::infinity();
        _workspace = alpha - _mu;
        return 0.5*(_rootH*_workspace).squaredNorm() + _norm;
    }
 
    Scalar operator()(ndarray::Array<Scalar const,1,1> const & alpha) const;
 
    void operator()(
        ndarray::Array<Scalar const,2,1> const & alpha,
        ndarray::Array<Scalar,1,1> const & output
    ) const;
 
protected:
    Scalar _norm;
    Vector _mu;
    mutable Vector _workspace;
    Matrix _rootH;
};
 
class TruncatedGaussianEvaluator {
public:
 
    explicit TruncatedGaussianEvaluator(TruncatedGaussian const & parent) :
        _internal(parent)
    {}
 
    template <typename Derived>
    Scalar operator()(Eigen::MatrixBase<Derived> const & alpha) const {
        return std::exp(-_internal(alpha));
    }
 
    Scalar operator()(ndarray::Array<Scalar const,1,1> const & alpha) const;
 
    void operator()(
        ndarray::Array<Scalar const,2,1> const & alpha,
        ndarray::Array<Scalar,1,1> const & output
    ) const;
 
private:
    TruncatedGaussianLogEvaluator _internal;
};
 
class TruncatedGaussianSampler {
public:
 
    explicit TruncatedGaussianSampler(
        TruncatedGaussian const & parent,
        TruncatedGaussian::SampleStrategy strategy
    );
 
    Scalar operator()(afw::math::Random & rng, ndarray::Array<Scalar,1,1> const & alpha) const;
 
    void operator()(
        afw::math::Random & rng,
        ndarray::Array<Scalar,2,1> const & alpha,
        ndarray::Array<Scalar,1,1> const & weights,
        bool multiplyWeights=false
    ) const;
 
    ~TruncatedGaussianSampler(); // defined in .cc so it can see Impl's dtor
 
    class Impl; // public so we can inherit from it in the .cc file
 
private:
    std::shared_ptr<Impl> _impl;
};
 
inline TruncatedGaussian::Sampler TruncatedGaussian::sample(SampleStrategy strategy) const {
    return Sampler(*this, strategy);
}
 
inline TruncatedGaussian::Sampler TruncatedGaussian::sample(Scalar minRejectionEfficiency) const {
    return Sampler(
        *this,
        (getUntruncatedFraction() < minRejectionEfficiency) ? ALIGN_AND_WEIGHT : DIRECT_WITH_REJECTION
    );
}
 
inline TruncatedGaussian::LogEvaluator TruncatedGaussian::evaluateLog() const { return LogEvaluator(*this); }
 
inline TruncatedGaussian::Evaluator TruncatedGaussian::evaluate() const { return Evaluator(*this); }
 
}}} // namespace lsst::meas::modelfit
 
#endif // !LSST_MEAS_MODELFIT_TruncatedGaussian_h_INCLUDED