Mixture.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_Mixture_h_INCLUDED
#define LSST_MEAS_MODELFIT_Mixture_h_INCLUDED
 
#include <limits>
 
#include "Eigen/Cholesky"
#include "Eigen/StdVector"
#include "Eigen/Dense"
 
#include "ndarray.h"
 
#include "lsst/base.h"
#include "lsst/afw/math/Random.h"
#include "lsst/afw/table/io/Persistable.h"
#include "lsst/meas/modelfit/common.h"
#include "lsst/afw/table/io/python.h"  // for declarePersistableFacade
 
 
namespace lsst { namespace meas { namespace modelfit {
 
class MixtureComponent {
public:
 
    int getDimension() const { return _mu.size(); }
 
    Scalar weight;
 
    Vector getMu() const { return _mu; }
    void setMu(Vector const & mu) { _mu = mu; }
 
 
    Matrix getSigma() const { return _sigmaLLT.reconstructedMatrix(); }
    void setSigma(Matrix const & sigma);
 
    MixtureComponent project(int dim) const;
 
    MixtureComponent project(int dim1, int dim2) const;
 
    explicit MixtureComponent(int dim);
 
    MixtureComponent(Scalar weight_, Vector const & mu, Matrix const & sigma);
 
    MixtureComponent & operator=(MixtureComponent const & other);
 
    friend std::ostream & operator<<(std::ostream & os, MixtureComponent const & self) {
        self._stream(os);
        return os;
    }
 
private:
 
    friend class Mixture;
 
    void _stream(std::ostream & os, int offset=0) const;
 
    Scalar _sqrtDet;
    Vector _mu;
    Eigen::LLT<Matrix> _sigmaLLT;
};
 
class MixtureUpdateRestriction {
public:
 
    int getDimension() const { return _dim; }
 
    virtual void restrictMu(Vector & mu) const {}
 
    virtual void restrictSigma(Matrix & sigma) const {}
 
    virtual ~MixtureUpdateRestriction() {}
 
    explicit MixtureUpdateRestriction(int dim) : _dim(dim) {}
 
private:
    int _dim;
};
 
class Mixture : public afw::table::io::PersistableFacade<Mixture>, public afw::table::io::Persistable {
public:
 
    typedef MixtureComponent Component;
    typedef MixtureUpdateRestriction UpdateRestriction;
    typedef std::vector<Component> ComponentList;
    typedef ComponentList::iterator iterator;
    typedef ComponentList::const_iterator const_iterator;
 
 
    iterator begin() { return _components.begin(); }
    iterator end() { return _components.end(); }
 
    const_iterator begin() const { return _components.begin(); }
    const_iterator end() const { return _components.end(); }
 
    Component & operator[](std::size_t i) { return _components[i]; }
    Component const & operator[](std::size_t i) const { return _components[i]; }
 
    std::size_t size() const { return _components.size(); }
 
    virtual int getComponentCount() const { return size(); }
 
    PTR(Mixture) project(int dim) const;
 
    PTR(Mixture) project(int dim1, int dim2) const;
 
    int getDimension() const { return _dim; }
 
    void normalize();
 
    void shift(int dim, Scalar offset);
 
    std::size_t clip(Scalar threshold=0.0);
 
    Scalar getDegreesOfFreedom() const { return _df; }
 
    void setDegreesOfFreedom(Scalar df=std::numeric_limits<Scalar>::infinity());
 
    template <typename Derived>
    Scalar evaluate(Component const & component, Eigen::MatrixBase<Derived> const & x) const {
        Scalar z = _computeZ(component, x);
        return component.weight * _evaluate(z) / component._sqrtDet;
    }
 
    template <typename Derived>
    Scalar evaluate(Eigen::MatrixBase<Derived> const & x) const {
        Scalar p = 0.0;
        for (const_iterator i = begin(); i != end(); ++i) {
            p += evaluate(*i, x);
        }
        return p;
    }
 
    void evaluate(
        ndarray::Array<Scalar const,2,1> const & x,
        ndarray::Array<Scalar,1,0> const & p
    ) const;
 
    void evaluateComponents(
        ndarray::Array<Scalar const,2,1> const & x,
        ndarray::Array<Scalar,2,1> const & p
    ) const;
 
    void evaluateDerivatives(
        ndarray::Array<Scalar const,1,1> const & x,
        ndarray::Array<Scalar,1,1> const & gradient,
        ndarray::Array<Scalar,2,1> const & hessian
    ) const;
 
    void evaluateDerivatives(
        Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & x,
        Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & gradient,
        Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> & hessian
    ) const;
 
   void evaluateDerivatives(
        Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & x,
        Eigen::Matrix<Scalar, Eigen::Dynamic, 1> & gradient
    ) const;
 
    void draw(afw::math::Random & rng, ndarray::Array<Scalar,2,1> const & x) const;
 
    void updateEM(
        ndarray::Array<Scalar const,2,1> const & x,
        ndarray::Array<Scalar const,1,0> const & w,
        Scalar tau1=0.0, Scalar tau2=0.5
    );
 
    void updateEM(
        ndarray::Array<Scalar const,2,1> const & x,
        ndarray::Array<Scalar const,1,0> const & w,
        UpdateRestriction const & restriction,
        Scalar tau1=0.0, Scalar tau2=0.5
    );
 
    void updateEM(
        ndarray::Array<Scalar const,2,1> const & x,
        UpdateRestriction const & restriction,
        Scalar tau1=0.0, Scalar tau2=0.5
    );
 
    virtual PTR(Mixture) clone() const;
 
    explicit Mixture(int dim, ComponentList & components, Scalar df=std::numeric_limits<Scalar>::infinity());
 
    friend std::ostream & operator<<(std::ostream & os, Mixture const & self) {
        self._stream(os);
        return os;
    }
 
    virtual bool isPersistable() const noexcept override { return true; }
 
protected:
 
    std::string getPythonModule() const override { return "lsst.meas.modelfit"; }
 
    std::string getPersistenceName() const override;
 
    void write(OutputArchiveHandle & handle) const override;
 
private:
    template <typename A, typename B, typename C>
    void _evaluateDerivativesImpl(A const & x,
                                  B & gradient,
                                  C * hessian,
                                  bool computeHessian = true) const;
 
    template <typename Derived>
    Scalar _computeZ(Component const & component, Eigen::MatrixBase<Derived> const & x) const {
        _workspace = x - component._mu;
        component._sigmaLLT.matrixL().solveInPlace(_workspace);
        return _workspace.squaredNorm();
    }
 
    // Helper function used in updateEM
    void updateDampedSigma(int k, Matrix const & sigma, double tau1, double tau2);
 
    Scalar _evaluate(Scalar z) const;
 
    void _stream(std::ostream & os) const;
 
    bool _isGaussian;
    int _dim;
    Scalar _df;
    Scalar _norm;
    mutable Vector _workspace;
    ComponentList _components;
};
 
}}} // namespace lsst::meas::modelfit
 
#endif // !LSST_MEAS_MODELFIT_Mixture_h_INCLUDED