[Rcpp-commits] r3217 - pkg/RcppEigen/src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Oct 22 00:49:04 CEST 2011
Author: dmbates
Date: 2011-10-22 00:49:03 +0200 (Sat, 22 Oct 2011)
New Revision: 3217
Modified:
pkg/RcppEigen/src/fastLm.cpp
pkg/RcppEigen/src/fastLm.h
Log:
Clean up code.
Modified: pkg/RcppEigen/src/fastLm.cpp
===================================================================
--- pkg/RcppEigen/src/fastLm.cpp 2011-10-21 22:48:32 UTC (rev 3216)
+++ pkg/RcppEigen/src/fastLm.cpp 2011-10-21 22:49:03 UTC (rev 3217)
@@ -22,263 +22,201 @@
#include "fastLm.h"
-using namespace Rcpp;
+namespace lmsol {
+ using Rcpp::_;
+ using Rcpp::as;
+ using Rcpp::CharacterVector;
+ using Rcpp::clone;
+ using Rcpp::List;
+ using Rcpp::NumericMatrix;
+ using Rcpp::NumericVector;
+ using Rcpp::RObject;
+ using Rcpp::wrap;
-lm::lm(const MMatrixXd &X, const MVectorXd &y)
- : m_n(X.rows()),
- m_p(X.cols()),
- m_coef(m_p),
- m_r(NA_INTEGER),
- m_df(m_n - m_p),
- m_perm(m_p),
- m_fitted(m_n),
- m_unsc(m_p, m_p),
- m_usePrescribedThreshold(false) {
-}
+ using std::invalid_argument;
+ using std::numeric_limits;
-lm& lm::setThreshold(const RealScalar& threshold)
-{
- m_usePrescribedThreshold = true;
- m_prescribedThreshold = threshold;
- return *this;
-}
+ lm::lm(const Map<MatrixXd> &X, const Map<VectorXd> &y)
+ : m_n(X.rows()),
+ m_p(X.cols()),
+ m_coef(m_p),
+ m_r(NA_INTEGER),
+ m_df(m_n - m_p),
+ m_perm(m_p),
+ m_fitted(m_n),
+ m_unsc(MatrixXd::Identity(m_p, m_p)),
+ m_usePrescribedThreshold(false) {
+ for (Index i = 0; i < m_p; ++i) m_perm[i] = i;
+ }
-/** Returns the threshold that will be used by certain methods such as rank().
- *
- */
-RealScalar lm::threshold() const
-{
- return m_usePrescribedThreshold ? m_prescribedThreshold
-// this formula comes from experimenting (see "LU precision tuning" thread on the list)
-// and turns out to be identical to Higham's formula used already in LDLt.
- : std::numeric_limits<double>::epsilon() * m_p;
-}
+ lm& lm::setThreshold(const RealScalar& threshold)
+ {
+ m_usePrescribedThreshold = true;
+ m_prescribedThreshold = threshold;
+ return *this;
+ }
-ColPivQR::ColPivQR(const MMatrixXd &X, const MVectorXd &y) : lm(X, y) {
- PivQRType PQR(X); // decompose the model matrix
- PermutationType Pmat = PQR.colsPermutation();
- m_perm = Pmat.indices();
- m_r = PQR.rank();
- MatrixXd R = PQR.matrixQR().topLeftCorner(m_p, m_p);
+ /** Returns the threshold that will be used by certain methods such as rank().
+ *
+ */
+ RealScalar lm::threshold() const
+ {
+ return m_usePrescribedThreshold ? m_prescribedThreshold
+ // this formula comes from experimenting (see "LU precision tuning" thread on the list)
+ // and turns out to be identical to Higham's formula used already in LDLt.
+ : numeric_limits<double>::epsilon() * m_p;
+ }
- if (m_r < (int)m_p) { // The rank-deficient case
- m_df = m_n - m_r;
- int nsing = m_p - m_r;
- MatrixXd Atrunc = (X * Pmat).leftCols(m_r);
- QRType QR(Atrunc);
- VectorXd coefTrunc = QR.solve(y);
- m_fitted = Atrunc * coefTrunc;
- m_coef.topRows(m_r) = coefTrunc;
- std::fill(m_coef.data() + m_r, m_coef.data() + m_p, NA_REAL);
- m_coef = Pmat * m_coef;
- MatrixXd Rtrunc = R.topLeftCorner(m_r, m_r);
- MatrixXd Rinvtrunc = Rtrunc.triangularView<Eigen::Upper>().solve(MatrixXd::Identity(m_r, m_r));
- m_unsc.topLeftCorner(m_r, m_r) = Rtrunc.setZero().selfadjointView<Eigen::Upper>().rankUpdate(Rinvtrunc);
- m_unsc.rightCols(nsing).fill(NA_REAL);
- m_unsc.bottomRows(nsing).fill(NA_REAL);
- } else { // full rank
- MatrixXd Rinv = R.triangularView<Eigen::Upper>().solve(MatrixXd::Identity(m_p, m_p));
- m_unsc = R.setZero().selfadjointView<Eigen::Upper>().rankUpdate(Rinv);
- m_coef = PQR.solve(y);
- m_fitted = X * m_coef;
+ inline MatrixXd AAt(const MatrixXd& A) {
+ return MatrixXd(MatrixXd(A.cols(), A.cols()).
+ setZero().selfadjointView<Lower>().rankUpdate(A));
}
-}
-QR::QR(const MMatrixXd &X, const MVectorXd &y) : lm(X, y) {
- QRType QR(X);
- m_coef = QR.solve(y);
- m_fitted = X * m_coef;
- MatrixXd R = QR.matrixQR().topLeftCorner(m_p, m_p);
- MatrixXd Rinv = R.triangularView<Eigen::Upper>().solve(MatrixXd::Identity(m_p, m_p));
- m_unsc = m_unsc.setZero().selfadjointView<Eigen::Upper>().rankUpdate(Rinv);
- for (Index i = 0; i < m_p; ++i) m_perm[i] = i;
-}
+ ColPivQR::ColPivQR(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
+ ColPivHouseholderQR<MatrixXd> PQR(X); // decompose the model matrix
+ Permutation Pmat(PQR.colsPermutation());
+ m_perm = Pmat.indices();
+ m_r = PQR.rank();
+ MatrixXd R(PQR.matrixQR().topRows(m_p).triangularView<Upper>());
+
+ if (m_r < (int)m_p) { // X is rank-deficient
+ m_df = m_n - m_r;
+ int nsing(m_p - m_r);
+ MatrixXd Atrunc((X * Pmat).leftCols(m_r));
+ HouseholderQR<MatrixXd> QR(Atrunc);
+ VectorXd coefTrunc(QR.solve(y));
+ m_fitted = Atrunc * coefTrunc;
+ m_coef.fill(::NA_REAL);
+ m_coef.topRows(m_r) = coefTrunc;
+ m_coef = Pmat * m_coef;
+ MatrixXd Rtrunc(MatrixXd(R).topLeftCorner(m_r, m_r));
+ MatrixXd Rinvtrunc(Rtrunc.triangularView<Upper>().
+ solve(MatrixXd::Identity(m_r, m_r)));
+ m_unsc.topLeftCorner(m_r,m_r) = AAt(Rinvtrunc);
+ m_unsc.rightCols(nsing).fill(NA_REAL);
+ m_unsc.bottomRows(nsing).fill(NA_REAL);
+ } else { // full rank X
+ MatrixXd Rinv(R.triangularView<Upper>().solve(m_unsc));
+ m_unsc = AAt(Rinv);
+ m_coef = PQR.solve(y);
+ m_fitted = X * m_coef;
+ }
+ }
+
+ QR::QR(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
+ HouseholderQR<MatrixXd> QR(X);
+ MatrixXd Rinv(QR.matrixQR().topRows(m_p).triangularView<Upper>().solve(m_unsc));
+ m_unsc = AAt(Rinv);
+ m_coef = QR.solve(y);
+ m_fitted = X * m_coef;
+ }
+
+ Llt::Llt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
+ LLT<MatrixXd> Ch(m_unsc.setZero().selfadjointView<Lower>().rankUpdate(X.adjoint()));
+ m_coef = Ch.solve(X.adjoint() * y);
+ m_fitted = X * m_coef;
+ m_unsc = Ch.solve(MatrixXd::Identity(m_p, m_p));
+ }
+
+ Ldlt::Ldlt(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
+ LDLT<MatrixXd> Ch(m_unsc.setZero().selfadjointView<Lower>().rankUpdate(X.adjoint()));
+ m_coef = Ch.solve(X.adjoint() * y);
+ m_fitted = X * m_coef;
+ m_unsc = Ch.solve(MatrixXd::Identity(m_p, m_p));
+ }
+
+ typedef Eigen::DiagonalMatrix<double, Eigen::Dynamic> DiagXd;
-LLT::LLT(const MMatrixXd &X, const MVectorXd &y) : lm(X, y) {
- LLTType Ch(m_unsc.setZero().selfadjointView<Eigen::Lower>().rankUpdate(X.adjoint()));
- m_coef = Ch.solve(X.adjoint() * y);
- m_fitted = X * m_coef;
- m_unsc = Ch.solve(MatrixXd::Identity(m_p, m_p));
- for (Index i = 0; i < m_p; ++i) m_perm[i] = i;
-}
+ inline DiagXd Dplus(const ArrayXd& D, Index r) {
+ VectorXd Di(VectorXd::Constant(D.size(), 0.));
+ for (Index i = 0; i < r; ++i) Di[i] = 1. / D[i];
+ return DiagXd(Di);
+ }
-LDLT::LDLT(const MMatrixXd &X, const MVectorXd &y) : lm(X, y) {
- LDLTType Ch(m_unsc.setZero().selfadjointView<Eigen::Lower>().rankUpdate(X.adjoint()));
- m_coef = Ch.solve(X.adjoint() * y);
- for (Index i = 0; i < m_p; ++i) m_perm[i] = i; // for now. Maybe add a boolean flag to indicate if unsc is already reordered
-// m_perm = PermutationType(Ch.transpositionsP()).indices();
- m_fitted = X * m_coef;
- m_unsc = Ch.solve(MatrixXd::Identity(m_p, m_p));
-}
+ SVD::SVD(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
+ JacobiSVD<MatrixXd> UDV(X.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV));
+ ArrayXd D(UDV.singularValues());
+ m_r = (D > D[0] * threshold()).count();
+ MatrixXd VDi(UDV.matrixV() * Dplus(D, m_r));
+ m_coef = VDi * UDV.matrixU().adjoint() * y;
+ m_fitted = X * m_coef;
+ m_unsc = AAt(VDi);
+ }
-// SVD method
-MatrixXd pseudoInverse(const MatrixXd& X, double tolerance) {
- SVDType UDV = X.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV);
- VectorXd D = UDV.singularValues();
- double tol = D[0] * tolerance;
- // There are better ways of doing this
- double *Dpt = D.data();
- for (int i = 0; i < D.size(); ++i)
- Dpt[i] = Dpt[i] < tol ? 0. : 1/Dpt[i];
-// Eigen2 code
-//UDV.matrixV() * (D.cwise() > tol).select(D.cwise().inverse(), 0).
-// asDiagonal() * UDV.matrixU().adjoint();
- return UDV.matrixV() * D.asDiagonal() * UDV.matrixU().adjoint();
-}
-
-SVD::SVD(const MMatrixXd &X, const MVectorXd &y) : lm(X, y) {
- SVDType UDV = X.jacobiSvd(Eigen::ComputeThinU|Eigen::ComputeThinV);
- VectorXd D = UDV.singularValues();
- m_r = (D.array() > threshold() * D[0]).count();
- m_coef = UDV.solve(y);
- m_fitted = X * m_coef;
- MatrixXd VDi = UDV.matrixV() * DiagType(UDV.singularValues().array().inverse().matrix());
- m_unsc = m_unsc.selfadjointView<Eigen::Lower>().rankUpdate(VDi);
- for (Index i = 0; i < m_p; ++i) m_perm[i] = i;
-}
-
-
-SymmEigen::SymmEigen(const MMatrixXd &X, const MVectorXd &y) : lm(X, y) {
- EigvType eig(m_unsc.setZero().selfadjointView<Eigen::Lower>().rankUpdate(X.adjoint()));
- ArrayXd D = eig.eigenvalues().array().sqrt(); // singular values of X
- MatrixXd VDi = eig.eigenvectors() * DiagType(D.inverse().matrix());
- m_coef = VDi * VDi.adjoint() * X.adjoint() * y;
- m_r = std::count_if(D.data(), D.data() + m_p,
- std::bind2nd(std::greater<double>(), D[0] * threshold()));
- m_fitted = X * m_coef;
- m_unsc = m_unsc.setZero().selfadjointView<Eigen::Lower>().rankUpdate(VDi);
- for (Index i = 0; i < m_p; ++i) m_perm[i] = i;
-}
-
-enum {ColPivQR_t = 0, QR_t, LDLT_t, LLT_t, SVD_t, SymmEigen_t};
-
-static inline lm do_lm(const MMatrixXd &X, const MVectorXd &y, int type)
- throw (std::invalid_argument) {
- switch(type) {
- case ColPivQR_t:
- return ColPivQR(X, y);
- case QR_t:
- return QR(X, y);
- case LLT_t:
- return LLT(X, y);
- case LDLT_t:
- return LDLT(X, y);
- case SVD_t:
- return SVD(X, y);
- case SymmEigen_t:
- return SymmEigen(X, y);
+ SymmEigen::SymmEigen(const Map<MatrixXd> &X, const Map<VectorXd> &y) : lm(X, y) {
+ SelfAdjointEigenSolver<MatrixXd>
+ eig(m_unsc.setZero().selfadjointView<Lower>().rankUpdate(X.adjoint()));
+ ArrayXd D(eig.eigenvalues().array().sqrt());
+ m_r = (D > D[0] * threshold()).count();
+ MatrixXd VDi(eig.eigenvectors() * Dplus(D, m_r));
+ m_unsc = AAt(VDi);
+ m_coef = m_unsc * X.adjoint() * y;
+ m_fitted = X * m_coef;
}
- throw std::invalid_argument("invalid type");
- return ColPivQR(X, y); // -Wall
-}
-extern "C" SEXP fastLm(SEXP Xs, SEXP ys, SEXP type) {
- try {
- const MMatrixXd X(as<MMatrixXd>(Xs));
- const MVectorXd y(as<MVectorXd>(ys));
- Index n = X.rows(), p = X.cols();
- if ((Index)y.size() != n)
- throw std::invalid_argument("size mismatch");
+ enum {ColPivQR_t = 0, QR_t, LLT_t, LDLT_t, SVD_t, SymmEigen_t};
- lm ans = do_lm(X, y, ::Rf_asInteger(type));
- // Copy coefficients and install names, if available
- NumericVector coef = wrap(ans.coef());
- List dimnames = NumericMatrix(Xs).attr("dimnames");
- if (dimnames.size() > 1) {
- RObject colnames = dimnames[1];
- if (!(colnames).isNULL())
- coef.attr("names") = clone(CharacterVector(colnames));
+ static inline lm do_lm(const Map<MatrixXd> &X, const Map<VectorXd> &y, int type) {
+ switch(type) {
+ case ColPivQR_t:
+ return ColPivQR(X, y);
+ case QR_t:
+ return QR(X, y);
+ case LLT_t:
+ return Llt(X, y);
+ case LDLT_t:
+ return Ldlt(X, y);
+ case SVD_t:
+ return SVD(X, y);
+ case SymmEigen_t:
+ return SymmEigen(X, y);
}
+ throw invalid_argument("invalid type");
+ return ColPivQR(X, y); // -Wall
+ }
+
+ extern "C" SEXP fastLm(SEXP Xs, SEXP ys, SEXP type) {
+ try {
+ const Map<MatrixXd> X(as<Map<MatrixXd> >(Xs));
+ const Map<VectorXd> y(as<Map<VectorXd> >(ys));
+ Index n = X.rows(), p = X.cols();
+ if ((Index)y.size() != n)
+ throw invalid_argument("size mismatch");
+ // Select and apply the least squares method
+ lm ans(do_lm(X, y, ::Rf_asInteger(type)));
+ // Copy coefficients
+ NumericVector coef(wrap(ans.coef()));
+ // Install names, if available
+ List dimnames(NumericMatrix(Xs).attr("dimnames"));
+ if (dimnames.size() > 1) {
+ RObject colnames = dimnames[1];
+ if (!(colnames).isNULL())
+ coef.attr("names") = clone(CharacterVector(colnames));
+ }
- VectorXd resid = y - ans.fitted();
- double s2 = resid.squaredNorm()/ans.df();
+ VectorXd resid = y - ans.fitted();
+ double s2 = resid.squaredNorm()/ans.df();
// Create the standard errors
- PermutationType Pmat = PermutationType(p);
- Pmat.indices() = ans.perm();
- VectorXd dd = Pmat * ans.unsc().diagonal();
- ArrayXd se = (dd.array() * s2).sqrt();
+ Permutation Pmat = Permutation(p);
+ Pmat.indices() = ans.perm();
+ VectorXd dd = Pmat * ans.unsc().diagonal();
+ ArrayXd se = (dd.array() * s2).sqrt();
- return List::create(_["coefficients"] = coef,
- _["se"] = se,
- _["rank"] = ans.rank(),
- _["df.residual"] = ans.df(),
- _["perm"] = ans.perm(),
- _["residuals"] = resid,
- _["s2"] = s2,
- _["fitted.values"] = ans.fitted(),
- _["unsc"] = ans.unsc());
+ return List::create(_["coefficients"] = coef,
+ _["se"] = se,
+ _["rank"] = ans.rank(),
+ _["df.residual"] = ans.df(),
+ _["perm"] = ans.perm(),
+ _["residuals"] = resid,
+ _["s2"] = s2,
+ _["fitted.values"] = ans.fitted(),
+ _["unsc"] = ans.unsc());
- } catch( std::exception &ex ) {
- forward_exception_to_r( ex );
- } catch(...) {
- ::Rf_error( "c++ exception (unknown reason)" );
+ } catch( std::exception &ex ) {
+ forward_exception_to_r( ex );
+ } catch(...) {
+ ::Rf_error( "c++ exception (unknown reason)" );
+ }
+ return R_NilValue; // -Wall
}
- return R_NilValue; // -Wall
}
-
-extern "C" SEXP crossprod(SEXP Xs) {
- try {
- const NumericMatrix X(Xs);
- Index n = X.nrow(), p = X.ncol();
- const MMatrixXd Xe(X.begin(), n, p);
- MatrixXd XtX(p, p);
- XtX = XtX.setZero().selfadjointView<Eigen::Lower>().rankUpdate(Xe.adjoint());
- return wrap(XtX);
- } catch( std::exception &ex ) {
- forward_exception_to_r( ex );
- } catch(...) {
- ::Rf_error( "c++ exception (unknown reason)" );
- }
- return R_NilValue; // -Wall
-}
-
-extern "C" SEXP crossprod1(SEXP Xs) {
- try {
- const NumericMatrix X(Xs);
- Index n = X.nrow(), p = X.ncol();
- const MMatrixXd Xe(X.begin(), n, p);
- MatrixXd XtX(p, p);
- XtX = XtX.setZero().selfadjointView<Eigen::Lower>().rankUpdate(Xe.adjoint());
- NumericMatrix ans(p, p);
- std::copy(XtX.data(), XtX.data() + XtX.size(), ans.begin());
- return ans;
- } catch( std::exception &ex ) {
- forward_exception_to_r( ex );
- } catch(...) {
- ::Rf_error( "c++ exception (unknown reason)" );
- }
- return R_NilValue; // -Wall
-}
-
-extern "C" SEXP tcrossprod(SEXP Xs) {
- try {
- const NumericMatrix X(Xs);
- Index n = X.nrow(), p = X.ncol();
- MatrixXd XXt(n, n);
- XXt = XXt.setZero().selfadjointView<Eigen::Lower>().rankUpdate(MMatrixXd(X.begin(), n, p));
- return wrap(XXt);
- } catch( std::exception &ex ) {
- forward_exception_to_r( ex );
- } catch(...) {
- ::Rf_error( "c++ exception (unknown reason)" );
- }
- return R_NilValue; // -Wall
-}
-
-extern "C" SEXP tcrossprod1(SEXP Xs) {
- try {
- const NumericMatrix X(Xs);
- Index n = X.nrow(), p = X.ncol();
- MatrixXd XXt(n, n);
- XXt = XXt.setZero().selfadjointView<Eigen::Lower>().rankUpdate(MMatrixXd(X.begin(), n, p));
- NumericMatrix ans(n, n);
- std::copy(XXt.data(), XXt.data() + XXt.size(), ans.begin());
- return ans;
- } catch( std::exception &ex ) {
- forward_exception_to_r( ex );
- } catch(...) {
- ::Rf_error( "c++ exception (unknown reason)" );
- }
- return R_NilValue; // -Wall
-}
-
Modified: pkg/RcppEigen/src/fastLm.h
===================================================================
--- pkg/RcppEigen/src/fastLm.h 2011-10-21 22:48:32 UTC (rev 3216)
+++ pkg/RcppEigen/src/fastLm.h 2011-10-21 22:49:03 UTC (rev 3217)
@@ -23,104 +23,84 @@
#define RCPPEIGEN_FASTLM_H
#include <RcppEigen.h>
- // Basic matrix, vector and array types for double
-typedef Eigen::MatrixXd MatrixXd;
-typedef Eigen::VectorXd VectorXd;
-typedef Eigen::ArrayXd ArrayXd;
-typedef Eigen::ArrayXXd ArrayXXd;
- // integer
-typedef Eigen::MatrixXi MatrixXi;
-typedef Eigen::VectorXi VectorXi;
-typedef Eigen::ArrayXi ArrayXi;
-typedef Eigen::ArrayXXi ArrayXXi;
- // complex
-typedef Eigen::MatrixXcd MatrixXcd;
-typedef Eigen::VectorXcd VectorXcd;
-typedef Eigen::ArrayXcd ArrayXcd;
-typedef Eigen::ArrayXXcd ArrayXXcd;
- // these should be defined for each base type - we use double
-typedef MatrixXd::Index Index;
-typedef MatrixXd::Scalar Scalar;
-typedef MatrixXd::RealScalar RealScalar;
- // Mapped matrix and vector types (share storage)
-typedef Eigen::Map<MatrixXd> MMatrixXd;
-typedef Eigen::Map<VectorXd> MVectorXd;
- // Views
-typedef Eigen::TriangularView<MatrixXd, Eigen::Upper> UpperTri;
-typedef Eigen::TriangularView<MatrixXd, Eigen::Lower> LowerTri;
-typedef Eigen::SelfAdjointView<MatrixXd, Eigen::Upper> UpperSym;
-typedef Eigen::SelfAdjointView<MatrixXd, Eigen::Lower> LowerSym;
- // Decomposition types
-typedef Eigen::LLT<MatrixXd> LLTType;
-typedef Eigen::LDLT<MatrixXd> LDLTType;
-typedef Eigen::ColPivHouseholderQR<MatrixXd> PivQRType;
-typedef Eigen::HouseholderQR<MatrixXd> QRType;
-typedef Eigen::JacobiSVD<MatrixXd> SVDType;
-typedef Eigen::DiagonalMatrix<double, Eigen::Dynamic> DiagType;
-typedef Eigen::SelfAdjointEigenSolver<MatrixXd> EigvType;
- // Types derived from decompositions
-typedef PivQRType::PermutationType PermutationType;
-typedef PermutationType::IndicesType IndicesType;
-class lm {
-protected:
- Index m_n; /**< number of rows */
- Index m_p; /**< number of columns */
- VectorXd m_coef; /**< coefficient vector */
- int m_r; /**< computed rank or NA_INTEGER */
- int m_df; /**< residual degrees of freedom */
- IndicesType m_perm; /**< column permutation */
- VectorXd m_fitted; /**< vector of fitted values */
- MatrixXd m_unsc; /**< unscaled variance-covariance matrix */
- RealScalar m_prescribedThreshold; /**< user specified tolerance */
- bool m_usePrescribedThreshold;
-public:
- lm(const MMatrixXd&, const MVectorXd&);
- lm& setThreshold(const RealScalar&); // patterned after ColPivHouseholderQR
- RealScalar threshold() const;
- const VectorXd& coef() const {return m_coef;}
- int rank() const {return m_r;}
- int df() const {return m_df;}
- const IndicesType& perm() const {return m_perm;}
- const VectorXd& fitted() const {return m_fitted;}
- const MatrixXd& unsc() const {return m_unsc;}
-};
+namespace lmsol {
+ using Eigen::ArrayXd;
+ using Eigen::ColPivHouseholderQR;
+ using Eigen::HouseholderQR;
+ using Eigen::JacobiSVD;
+ using Eigen::LDLT;
+ using Eigen::LLT;
+ using Eigen::Lower;
+ using Eigen::Map;
+ using Eigen::MatrixXd;
+ using Eigen::SelfAdjointEigenSolver;
+ using Eigen::TriangularView;
+ using Eigen::VectorXd;
+ using Eigen::Upper;
-class ColPivQR : public lm {
-public:
- ColPivQR(const MMatrixXd&, const MVectorXd&);
-};
+ typedef MatrixXd::Index Index;
+ typedef MatrixXd::Scalar Scalar;
+ typedef MatrixXd::RealScalar RealScalar;
+ typedef ColPivHouseholderQR<MatrixXd>::PermutationType Permutation;
+ typedef Permutation::IndicesType Indices;
-class LLT : public lm {
-public:
- LLT(const MMatrixXd&, const MVectorXd&);
-};
+ class lm {
+ protected:
+ Index m_n; /**< number of rows */
+ Index m_p; /**< number of columns */
+ VectorXd m_coef; /**< coefficient vector */
+ int m_r; /**< computed rank or NA_INTEGER */
+ int m_df; /**< residual degrees of freedom */
+ Indices m_perm; /**< column permutation */
+ VectorXd m_fitted; /**< vector of fitted values */
+ MatrixXd m_unsc; /**< unscaled variance-covariance */
+ RealScalar m_prescribedThreshold; /**< user specified tolerance */
+ bool m_usePrescribedThreshold;
+ public:
+ lm(const Map<MatrixXd>&, const Map<VectorXd>&);
+ lm& setThreshold(const RealScalar&); // patterned after ColPivHouseholderQR code
+ RealScalar threshold() const;
+ const VectorXd& coef() const {return m_coef;}
+ int rank() const {return m_r;}
+ int df() const {return m_df;}
+ const Indices& perm() const {return m_perm;}
+ const VectorXd& fitted() const {return m_fitted;}
+ const MatrixXd& unsc() const {return m_unsc;}
+ };
-class LDLT : public lm {
-public:
- LDLT(const MMatrixXd&, const MVectorXd&);
-};
+ class ColPivQR : public lm {
+ public:
+ ColPivQR(const Map<MatrixXd>&, const Map<VectorXd>&);
+ };
-class QR : public lm {
-public:
- QR(const MMatrixXd&, const MVectorXd&);
-};
+ class Llt : public lm {
+ public:
+ Llt(const Map<MatrixXd>&, const Map<VectorXd>&);
+ };
-class SVD : public lm {
-public:
- SVD(const MMatrixXd&, const MVectorXd&);
-};
+ class Ldlt : public lm {
+ public:
+ Ldlt(const Map<MatrixXd>&, const Map<VectorXd>&);
+ };
-class SymmEigen : public lm {
-public:
- SymmEigen(const MMatrixXd&, const MVectorXd&);
-};
+ class QR : public lm {
+ public:
+ QR(const Map<MatrixXd>&, const Map<VectorXd>&);
+ };
-extern "C" SEXP fastLm(SEXP Xs, SEXP ys, SEXP types);
+ class SVD : public lm {
+ public:
+ SVD(const Map<MatrixXd>&, const Map<VectorXd>&);
+ };
-extern "C" SEXP crossprod(SEXP Xs);
+ class SymmEigen : public lm {
+ public:
+ SymmEigen(const Map<MatrixXd>&, const Map<VectorXd>&);
+ };
+}
-extern "C" SEXP tcrossprod(SEXP Xs);
+extern "C" SEXP fastLm(SEXP Xs, SEXP ys, SEXP types);
#endif
More information about the Rcpp-commits
mailing list