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 {
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() = design.asEigen().template cast<double>();
        _getDataVector() = data.asEigen().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() = design.asEigen().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() = fisher.asEigen().template cast<double>();
        } else {
            _getFisherMatrix() = fisher.asEigen().transpose().template cast<double>();
        }
        _getRhsVector() = rhs.asEigen().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);

    // 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);

    PTR(Impl) _impl;
};

}}} // namespace lsst::afw::math

#endif // !LSST_AFW_MATH_LeastSquares_h_INCLUDED