LeastSquares.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_AFW_MATH_LeastSquares_h_INCLUDED
#define LSST_AFW_MATH_LeastSquares_h_INCLUDED

#include "lsst/base.h"
#include "ndarray/eigen.h"

namespace lsst {
namespace afw {
namespace math {

class LeastSquares final {
public:
    class Impl;  

    enum Factorization {
        NORMAL_EIGENSYSTEM, 
        NORMAL_CHOLESKY,    
        DIRECT_SVD          
    };

    template <typename T1, typename T2, int C1, int C2>
    static LeastSquares fromDesignMatrix(ndarray::Array<T1, 2, C1> const& design,
                                         ndarray::Array<T2, 1, C2> const& data,
                                         Factorization factorization = NORMAL_EIGENSYSTEM) {
        LeastSquares r(factorization, design.template getSize<1>());
        r.setDesignMatrix(design, data);
        return r;
    }

    template <typename D1, typename D2>
    static LeastSquares fromDesignMatrix(Eigen::MatrixBase<D1> const& design,
                                         Eigen::MatrixBase<D2> const& data,
                                         Factorization factorization = NORMAL_EIGENSYSTEM) {
        LeastSquares r(factorization, design.cols());
        r.setDesignMatrix(design, data);
        return r;
    }

    template <typename T1, typename T2, int C1, int C2>
    void setDesignMatrix(ndarray::Array<T1, 2, C1> const& design, ndarray::Array<T2, 1, C2> const& data) {
        // n.b. "template cast<T>" below is not a special kind of cast; it's just
        // the weird C++ syntax required for calling a templated member function
        // called "cast" in this context; see
        // http://eigen.tuxfamily.org/dox-devel/TopicTemplateKeyword.html
        _getDesignMatrix() = ndarray::asEigenMatrix(design).template cast<double>();
        _getDataVector() = ndarray::asEigenMatrix(data).template cast<double>();
        _factor(false);
    }

    template <typename D1, typename D2>
    void setDesignMatrix(Eigen::MatrixBase<D1> const& design, Eigen::MatrixBase<D2> const& data) {
        _getDesignMatrix() = design.template cast<double>();
        _getDataVector() = data.template cast<double>();
        _factor(false);
    }

    template <typename T1, int C1>
    void setDesignMatrix(ndarray::Array<T1, 2, C1> const& design) {
        _getDesignMatrix() = ndarray::asEigenMatrix(design).template cast<double>();
        _factor(false);
    }

    template <typename D1, typename D2>
    void setDesignMatrix(Eigen::MatrixBase<D1> const& design) {
        _getDesignMatrix() = design.template cast<double>();
        _factor(false);
    }

    template <typename T1, typename T2, int C1, int C2>
    static LeastSquares fromNormalEquations(ndarray::Array<T1, 2, C1> const& fisher,
                                            ndarray::Array<T2, 1, C2> const& rhs,
                                            Factorization factorization = NORMAL_EIGENSYSTEM) {
        LeastSquares r(factorization, fisher.template getSize<0>());
        r.setNormalEquations(fisher, rhs);
        return r;
    }

    template <typename D1, typename D2>
    static LeastSquares fromNormalEquations(Eigen::MatrixBase<D1> const& fisher,
                                            Eigen::MatrixBase<D2> const& rhs,
                                            Factorization factorization = NORMAL_EIGENSYSTEM) {
        LeastSquares r(factorization, fisher.rows());
        r.setNormalEquations(fisher, rhs);
        return r;
    }

    template <typename T1, typename T2, int C1, int C2>
    void setNormalEquations(ndarray::Array<T1, 2, C1> const& fisher, ndarray::Array<T2, 1, C2> const& rhs) {
        if ((C1 > 0) == bool(Eigen::MatrixXd::IsRowMajor)) {
            _getFisherMatrix() = ndarray::asEigenMatrix(fisher).template cast<double>();
        } else {
            _getFisherMatrix() = ndarray::asEigenMatrix(fisher).transpose().template cast<double>();
        }
        _getRhsVector() = ndarray::asEigenMatrix(rhs).template cast<double>();
        _factor(true);
    }

    template <typename D1, typename D2>
    void setNormalEquations(Eigen::MatrixBase<D1> const& fisher, Eigen::MatrixBase<D2> const& rhs) {
        if (bool(Eigen::MatrixBase<D1>::IsRowMajor) == bool(Eigen::MatrixXd::IsRowMajor)) {
            _getFisherMatrix() = fisher.template cast<double>();
        } else {
            _getFisherMatrix() = fisher.transpose().template cast<double>();
        }
        _getRhsVector() = rhs.template cast<double>();
        _factor(true);
    }

    void setThreshold(double threshold);

    double getThreshold() const;

    ndarray::Array<double const, 1, 1> getSolution();

    ndarray::Array<double const, 2, 2> getCovariance();

    ndarray::Array<double const, 2, 2> getFisherMatrix();

    ndarray::Array<double const, 1, 1> getDiagnostic(Factorization factorization);

    int getDimension() const;

    int getRank() const;

    Factorization getFactorization() const;

    LeastSquares(Factorization factorization, int dimension);

    LeastSquares(LeastSquares const&);
    LeastSquares(LeastSquares&&);
    LeastSquares& operator=(LeastSquares const&);
    LeastSquares& operator=(LeastSquares&&);

    // Need to define dtor in source file so it can see Impl declaration.
    ~LeastSquares();

private:
    // We want a column-major design matrix so the self-adjoint product is cache-friendly, so we always
    // copy a (possibly row-major) design matrix into a col-major one.  This is an unnecessarily and
    // cache-unfriendly operation when solver is DIRECT_SVD, but right now it doesn't seem to be worth
    // special-casing the design for that case.  In other cases it's a cache-unfriendly op that avoids
    // an even worse one, and it can always be avoided by using a column-major design matrix.
    Eigen::MatrixXd& _getDesignMatrix();
    Eigen::VectorXd& _getDataVector();

    // Storage order matters less here, so we just go with what Eigen is most accustomed to.
    // Because the Fisher matrix is symmetric, we can just transpose it before copying to avoid doing
    // any expensive copies between different storage orders.
    Eigen::MatrixXd& _getFisherMatrix();
    Eigen::VectorXd& _getRhsVector();

    // Do the actual high-level linear algrebra factorization; the appropriate input matrices and
    // vectors must already be set before this is called.  The solution and other quantities
    // will not be computed until they are first requested.
    void _factor(bool haveNormalEquations);

    std::shared_ptr<Impl> _impl;
};
}  // namespace math
}  // namespace afw
}  // namespace lsst

#endif  // !LSST_AFW_MATH_LeastSquares_h_INCLUDED