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(PTR(Impl) impl) : _impl(impl) {}

    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:
    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