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