Skip to content

Commit e4afb2d

Browse files
eddelbuetteleddelbuettel
committed
tinytest conversion step two: conversion
1 parent 297667c commit e4afb2d

File tree

8 files changed

+523
-548
lines changed

8 files changed

+523
-548
lines changed

‎ChangeLog

+14
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,17 @@
1+
2019-10-30 Dirk Eddelbuettel <[email protected]>
2+
3+
* test_sparse.R: Converted to tinytest
4+
* cpp/sparse.cpp: Added
5+
6+
2019-10-29 Dirk Eddelbuettel <[email protected]>
7+
8+
* test_fastLm.R: Converted to tinytest
9+
* test_RcppEigen.R: Idem
10+
* test_solution.R: Idem
11+
12+
* cpp/rcppeigen.cpp: Added
13+
* cpp/solution.cpp: Idem
14+
115
2019-10-28 Dirk Eddelbuettel <[email protected]>
216

317
* tests/tinytest.R: Renamed from tests/doRUnit.R

‎inst/tinytest/cpp/rcppeigen.cpp

+157
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,157 @@
1+
2+
#include <RcppEigen.h>
3+
4+
// [[Rcpp::depends(RcppEigen)]]
5+
6+
// [[Rcpp::export]]
7+
Rcpp::List fx() {
8+
Rcpp::List vecs = Rcpp::List::create(
9+
Rcpp::_["Vec<complex>"] = Eigen::VectorXcd::Zero(5),
10+
Rcpp::_["Vec<double>"] = Eigen::VectorXd::Zero(5),
11+
Rcpp::_["Vec<float>"] = Eigen::VectorXf::Zero(5),
12+
Rcpp::_["Vec<int>"] = Eigen::VectorXi::Zero(5)
13+
);
14+
15+
// A VectorX<T> behaves as a matrix with one column but is converted to
16+
// a vector object in R, not a matrix of one column. The distinction is
17+
// that VectorX<T> objects are defined at compile time to have one column,
18+
// whereas a MatrixX<T> has a dynamic number of columns that is set to 1
19+
// during execution of the code. A MatrixX<T> object can be resized to have
20+
// a different number of columns. A VectorX<T> object cannot.
21+
Rcpp::List cols = Rcpp::List::create(
22+
Rcpp::_["Col<complex>"] = Eigen::MatrixXcd::Zero(5, 1),
23+
Rcpp::_["Col<double>"] = Eigen::MatrixXd::Zero(5, 1),
24+
Rcpp::_["Col<float>"] = Eigen::MatrixXf::Zero(5, 1),
25+
Rcpp::_["Col<int>"] = Eigen::MatrixXi::Zero(5, 1)
26+
);
27+
28+
Rcpp::List rows = Rcpp::List::create(
29+
Rcpp::_["Row<complex>"] = Eigen::RowVectorXcd::Zero(5),
30+
Rcpp::_["Row<double>"] = Eigen::RowVectorXd::Zero(5),
31+
Rcpp::_["Row<float>"] = Eigen::RowVectorXf::Zero(5),
32+
Rcpp::_["Row<int>"] = Eigen::RowVectorXi::Zero(5)
33+
);
34+
35+
Rcpp::List matrices = Rcpp::List::create(
36+
Rcpp::_["Mat<complex>"] = Eigen::MatrixXcd::Identity(3, 3),
37+
Rcpp::_["Mat<double>"] = Eigen::MatrixXd::Identity(3, 3),
38+
Rcpp::_["Mat<float>"] = Eigen::MatrixXf::Identity(3, 3),
39+
Rcpp::_["Mat<int>"] = Eigen::MatrixXi::Identity(3, 3)
40+
);
41+
42+
// ArrayXX<t> objects have the same structure as matrices but allow
43+
// componentwise arithmetic. A * B is matrix multiplication for
44+
// matrices and componentwise multiplication for arrays.
45+
Rcpp::List arrays2 = Rcpp::List::create(
46+
Rcpp::_["Arr2<complex>"] = Eigen::ArrayXXcd::Zero(3, 3),
47+
Rcpp::_["Arr2<double>"] = Eigen::ArrayXXd::Zero(3, 3),
48+
Rcpp::_["Arr2<float>"] = Eigen::ArrayXXf::Zero(3, 3),
49+
Rcpp::_["Arr2<int>"] = Eigen::ArrayXXi::Zero(3, 3)
50+
);
51+
52+
// ArrayX<t> objects have the same structure as VectorX<T> objects
53+
// but allow componentwise arithmetic, including functions like exp, log,
54+
// sqrt, ...
55+
Rcpp::List arrays1 = Rcpp::List::create(
56+
Rcpp::_["Arr1<complex>"] = Eigen::ArrayXcd::Zero(5),
57+
Rcpp::_["Arr1<double>"] = Eigen::ArrayXd::Zero(5),
58+
Rcpp::_["Arr1<float>"] = Eigen::ArrayXf::Zero(5),
59+
Rcpp::_["Arr1<int>"] = Eigen::ArrayXi::Zero(5)
60+
);
61+
62+
Rcpp::List operations = Rcpp::List::create(
63+
Rcpp::_["Op_seq"] = Eigen::ArrayXd::LinSpaced(6, 1, 10), // arguments are length.out, start, end
64+
Rcpp::_["Op_log"] = Eigen::ArrayXd::LinSpaced(6, 1, 10).log(),
65+
Rcpp::_["Op_exp"] = Eigen::ArrayXd::LinSpaced(6, 1, 10).exp(),
66+
Rcpp::_["Op_sqrt"] = Eigen::ArrayXd::LinSpaced(6, 1, 10).sqrt(),
67+
Rcpp::_["Op_cos"] = Eigen::ArrayXd::LinSpaced(6, 1, 10).cos()
68+
);
69+
70+
Rcpp::List output = Rcpp::List::create(
71+
Rcpp::_["vectors : VectorX<T>"] = vecs,
72+
Rcpp::_["matrices : MatrixX<T>"] = matrices,
73+
Rcpp::_["rows : RowVectorX<T>"] = rows,
74+
Rcpp::_["columns : MatrixX<T>"] = cols,
75+
Rcpp::_["arrays2d : ArrayXX<T>"] = arrays2,
76+
Rcpp::_["arrays1d : ArrayX<T>"] = arrays1,
77+
Rcpp::_["operations : ArrayXd"] = operations
78+
);
79+
80+
return output ;
81+
}
82+
83+
// [[Rcpp::export]]
84+
Rcpp::List fx2(Rcpp::List input) {
85+
Eigen::VectorXi m1 = input[0] ; /* implicit as */
86+
Eigen::VectorXd m2 = input[1] ; /* implicit as */
87+
Eigen::Matrix<unsigned int, Eigen::Dynamic, 1> m3 = input[0] ; /* implicit as */
88+
Eigen::VectorXf m4 = input[1] ; /* implicit as */
89+
90+
Rcpp::List res = Rcpp::List::create(m1.sum(), m2.sum(), m3.sum(), m4.sum());
91+
92+
return res ;
93+
}
94+
95+
96+
// [[Rcpp::export]]
97+
Rcpp::List fx3(Rcpp::List input) {
98+
99+
const Eigen::Map<Eigen::VectorXi> m1 = input[0] ; // maps share storage and do not allow conversion
100+
const Eigen::Map<Eigen::VectorXd> m2 = input[1] ;
101+
102+
Rcpp::List res = Rcpp::List::create(m1.sum(), m2.sum());
103+
104+
return res ;
105+
}
106+
107+
// [[Rcpp::export]]
108+
Rcpp::List fx4(Rcpp::List input) {
109+
110+
const Eigen::Map<Eigen::RowVectorXi> m1 = input[0] ; // maps share storage, do not allow conversion
111+
const Eigen::Map<Eigen::RowVectorXd> m2 = input[1] ;
112+
113+
Rcpp::List res = Rcpp::List::create(m1.sum(), m2.sum());
114+
115+
return res ;
116+
}
117+
118+
119+
// [[Rcpp::export]]
120+
Rcpp::List fx5(Rcpp::List input) {
121+
const Eigen::Map<Eigen::MatrixXi> m1 = input[0]; // maps share storage, do not allow conversion
122+
const Eigen::Map<Eigen::MatrixXd> m2 = input[1] ;
123+
// FIXME: Write a version of as specifically for complex matrices.
124+
// const Eigen::Map<Eigen::MatrixXcd> m3 = input[2] ;
125+
126+
Rcpp::List res = Rcpp::List::create(m1.sum(), m2.sum());//, m3.sum());
127+
128+
return res ;
129+
}
130+
131+
132+
// [[Rcpp::export]]
133+
Rcpp::List fx6(Rcpp::List input) {
134+
const Eigen::MappedSparseMatrix<double> m1 = input[0]; // maps share storage and do not allow conversion
135+
136+
Rcpp::List res = Rcpp::List::create(Rcpp::_["nnz"] = double(m1.nonZeros()),
137+
Rcpp::_["nr"] = double(m1.rows()),
138+
Rcpp::_["nc"] = double(m1.cols()),
139+
Rcpp::_["inSz"] = double(m1.innerSize()),
140+
Rcpp::_["outSz"] = double(m1.outerSize()),
141+
Rcpp::_["sum"] = m1.sum());
142+
143+
return res ;
144+
}
145+
146+
147+
// [[Rcpp::export]]
148+
Rcpp::List fx7(Rcpp::List input) {
149+
const Eigen::SparseMatrix<double> m1 = input[0];
150+
Rcpp::List res = Rcpp::List::create(Rcpp::_["nnz"] = double(m1.nonZeros()),
151+
Rcpp::_["nr"] = double(m1.rows()),
152+
Rcpp::_["nc"] = double(m1.cols()),
153+
Rcpp::_["inSz"] = double(m1.innerSize()),
154+
Rcpp::_["outSz"] = double(m1.outerSize()),
155+
Rcpp::_["sum"] = m1.sum());
156+
return res ;
157+
}

‎inst/tinytest/cpp/solution.cpp

+37
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,37 @@
1+
2+
#include <RcppEigen.h>
3+
4+
// [[Rcpp::depends(RcppEigen)]]
5+
6+
typedef Eigen::ArrayXd Ar1;
7+
typedef Eigen::Map<Ar1> MAr1;
8+
typedef Eigen::ArrayXXd Ar2;
9+
typedef Eigen::Map<Ar2> MAr2;
10+
typedef Eigen::MatrixXd Mat;
11+
typedef Eigen::Map<Mat> MMat;
12+
typedef Eigen::VectorXd Vec;
13+
typedef Eigen::Map<Vec> MVec;
14+
typedef Eigen::PartialPivLU<Mat> PPLU;
15+
typedef Eigen::ColPivHouseholderQR<Mat> CPQR;
16+
17+
18+
// [[Rcpp::export]]
19+
Rcpp::List dense_PPLU(MMat A, MVec b) {
20+
PPLU lu(A);
21+
Mat Ainv(lu.inverse());
22+
Vec x(lu.solve(b));
23+
24+
return Rcpp::List::create(Rcpp::Named("A", A),
25+
Rcpp::Named("Ainv", Ainv),
26+
Rcpp::Named("b", b),
27+
Rcpp::Named("x", x));
28+
}
29+
30+
// [[Rcpp::export]]
31+
Rcpp::List dense_CPQR(MMat A, MVec b) {
32+
CPQR qr(A);
33+
Mat Ainv(qr.inverse());
34+
Vec x(qr.solve(b));
35+
return Rcpp::List::create(Rcpp::Named("Ainv", Ainv),
36+
Rcpp::Named("x", x));
37+
}

‎inst/tinytest/cpp/sparse.cpp

+110
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
2+
#include <RcppEigen.h>
3+
4+
// [[Rcpp::depends(RcppEigen)]]
5+
6+
// [[Rcpp::export]]
7+
SEXP wrapSparseDouble() {
8+
Eigen::SparseMatrix<double> mm(9,3);
9+
mm.reserve(9);
10+
for (int j = 0; j < 3; ++j) {
11+
mm.startVec(j);
12+
for (int i = 3 * j; i < 3 * (j + 1); ++i)
13+
mm.insertBack(i, j) = 1.;
14+
}
15+
mm.finalize();
16+
return Rcpp::wrap(mm);
17+
}
18+
19+
// [[Rcpp::export]]
20+
SEXP wrapSparseDoubleColumnMajor() {
21+
Eigen::SparseMatrix<double, Eigen::ColMajor> mm(9,3);
22+
mm.reserve(9);
23+
for (int j = 0; j < 3; ++j) {
24+
mm.startVec(j);
25+
for (int i = 3 * j; i < 3 * (j + 1); ++i)
26+
mm.insertBack(i, j) = 1.;
27+
}
28+
mm.finalize();
29+
return Rcpp::wrap(mm);
30+
}
31+
32+
// [[Rcpp::export]]
33+
SEXP wrapSparseDoubleRowMajor() {
34+
Eigen::SparseMatrix<double, Eigen::RowMajor> mm(9,3);
35+
mm.reserve(9);
36+
for (int irow = 0; irow < 9; ++irow) {
37+
mm.startVec(irow);
38+
mm.insertBack(irow, irow / 3) = static_cast<double>( 9 - irow );
39+
}
40+
mm.finalize();
41+
return Rcpp::wrap(mm);
42+
}
43+
44+
// [[Rcpp::export]]
45+
SEXP asSparseDoubleColumnMajor(SEXP R_mm) {
46+
Eigen::SparseMatrix<double, Eigen::ColMajor> mm =
47+
Rcpp::as<Eigen::SparseMatrix<double, Eigen::ColMajor> >( R_mm );
48+
return Rcpp::wrap(mm);
49+
}
50+
51+
// [[Rcpp::export]]
52+
SEXP asMappedSparseDoubleColMajor(SEXP R_mm) {
53+
typedef Eigen::Map<Eigen::SparseMatrix<double, Eigen::ColMajor> > MapMat;
54+
MapMat mm = Rcpp::as<MapMat>( R_mm );
55+
return Rcpp::wrap(mm);
56+
}
57+
58+
// [[Rcpp::export]]
59+
SEXP asMappedSparseDeprecatedDoubleColMajor(SEXP R_mm) {
60+
// Deprecated
61+
typedef Eigen::MappedSparseMatrix<double, Eigen::ColMajor> MapMat;
62+
MapMat mm = Rcpp::as<MapMat>( R_mm );
63+
return Rcpp::wrap(mm);
64+
}
65+
66+
// [[Rcpp::export]]
67+
SEXP asSparseDoubleRowMajor(SEXP R_mm) {
68+
Eigen::SparseMatrix<double, Eigen::RowMajor> mm =
69+
Rcpp::as<Eigen::SparseMatrix<double, Eigen::RowMajor> >( R_mm );
70+
return Rcpp::wrap(mm);
71+
}
72+
73+
// [[Rcpp::export]]
74+
SEXP asMappedSparseDoubleRowMajor(SEXP R_mm) {
75+
typedef Eigen::Map<Eigen::SparseMatrix<double, Eigen::RowMajor> > MapMat;
76+
MapMat mm = Rcpp::as<MapMat>( R_mm );
77+
return Rcpp::wrap(mm);
78+
}
79+
80+
// [[Rcpp::export]]
81+
SEXP asMappedSparseDeprecatedDoubleRowMajor(SEXP R_mm) {
82+
// Deprecated
83+
typedef Eigen::MappedSparseMatrix<double, Eigen::RowMajor> MapMat;
84+
MapMat mm = Rcpp::as<MapMat>( R_mm );
85+
return Rcpp::wrap(mm);
86+
}
87+
88+
// [[Rcpp::export]]
89+
Rcpp::List sparseCholesky(Rcpp::List input) {
90+
using Eigen::VectorXd;
91+
using Eigen::MatrixXd;
92+
using Eigen::Lower;
93+
using Eigen::Map;
94+
using Eigen::SparseMatrix;
95+
using Eigen::SimplicialLDLT;
96+
using Eigen::Success;
97+
98+
const Map<SparseMatrix<double> > m1 = input[0];
99+
const Map<VectorXd> v1 = input[1];
100+
SparseMatrix<double> m2(m1.cols(), m1.cols());
101+
m2.selfadjointView<Lower>().rankUpdate(m1.adjoint());
102+
103+
SimplicialLDLT<SparseMatrix<double> > ff(m2);
104+
VectorXd res = ff.solve(m1.adjoint() * v1);
105+
106+
return Rcpp::List::create(Rcpp::Named("res") = res,
107+
Rcpp::Named("rows") = double(ff.rows()),
108+
Rcpp::Named("cols") = double(ff.cols()));
109+
110+
}

0 commit comments

Comments
 (0)