1 #include "purify/config.h"
5 #include "catch2/catch_all.hpp"
12 SECTION(
"No overlap") {
13 std::vector<t_int> indices{1, 3, 5, 2};
14 Vector<t_int> input(10);
15 std::iota(input.data(), input.data() + input.size(), 0);
20 mapper(input, output);
21 SECTION(
"Direct application") {
22 CHECK(output.size() == indices.size());
23 decltype(output)::Index i(0);
24 for (
auto const &index : indices)
CHECK(output(i++) == index);
27 Vector<t_int> adjoint;
28 mapper.adjoint(output, adjoint);
29 SECTION(
"Adjoint application") {
30 CHECK(adjoint.size() == input.size());
31 for (
auto const &index : indices) {
32 CHECK(adjoint(index) == index);
35 CHECK(adjoint == adjoint.Zero(adjoint.size()));
38 SECTION(
"With overlap") {
39 std::vector<t_int> indices{1, 2, 5, 2};
40 Vector<t_int> input(10);
41 std::iota(input.data(), input.data() + input.size(), 0);
46 mapper(input, output);
47 SECTION(
"Direct application") {
48 CHECK(output.size() == indices.size());
49 decltype(output)::Index i(0);
50 for (
auto const &index : indices)
CHECK(output(i++) == index);
53 SECTION(
"Adjoint application") {
54 Vector<t_int> adjoint;
55 mapper.adjoint(output, adjoint);
56 CHECK(adjoint.size() == input.size());
57 for (
auto const &index : std::set<t_int>(indices.begin(), indices.end())) {
58 CHECK(adjoint(index) == (index == 2 ? 2 * index : index));
61 CHECK(adjoint == adjoint.Zero(adjoint.size()));
68 std::vector<Eigen::Triplet<t_int>>
const triplets = {{0, 0, 1}, {0, 3, 1}, {2, 3, 1}};
69 matrix.setFromTriplets(triplets.begin(), triplets.end());
72 CHECK(indices.size() == 2);
73 CHECK(indices.count(0) == 1);
74 CHECK(indices.count(3) == 1);
79 std::vector<Eigen::Triplet<t_int>>
const triplets = {{0, 0, 1}, {0, 3, 2}, {2, 3, 3}};
80 matrix.setFromTriplets(triplets.begin(), triplets.end());
83 CHECK(shrinked_matrix.nonZeros() == matrix.nonZeros());
84 CHECK(shrinked_matrix.coeffRef(0, 0) == matrix.coeffRef(0, 0));
85 CHECK(shrinked_matrix.coeffRef(0, 1) == matrix.coeffRef(0, 3));
86 CHECK(shrinked_matrix.coeffRef(2, 1) == matrix.coeffRef(2, 3));
89 TEST_CASE(
"Vector-shrinked matrix multiplication") {
91 extern std::unique_ptr<std::mt19937_64>
mersenne;
92 std::uniform_int_distribution<sopt::t_int> uniform_dist(0, N - 1);
93 Matrix<t_int> matrix = (Image<uint8_t>::Random(N, N) < uint8_t(64))
95 .select(Matrix<uint8_t>::Random(N, N), Matrix<uint8_t>::Zero(N, N))
98 matrix.row(uniform_dist(*
mersenne)).fill(0);
99 matrix.row(uniform_dist(*
mersenne)).fill(0);
101 Vector<t_int>
const input = Vector<uint8_t>::Random(N).cast<t_int>();
105 Vector<t_int> comp_input;
106 mapper(input, comp_input);
108 CHECK(sparse * input == compressed * comp_input);
#define CHECK(CONDITION, ERROR)
std::unique_ptr< std::mt19937_64 > mersenne(new std::mt19937_64(0))
TEST_CASE("Index mapping operator")
std::enable_if< std::is_same< typename T0::Scalar, typename T1::Scalar >::value and T0::IsRowMajor, Vector< typename T0::Scalar > >::type sparse_multiply_matrix(const Eigen::SparseMatrixBase< T0 > &M, const Eigen::MatrixBase< T1 > &x)
Parallel multiplication with a sparse matrix and vector.
std::set< STORAGE_INDEX_TYPE > non_empty_outers(Eigen::SparseMatrixBase< T0 > const &matrix)
Indices of non empty outer indices.
Eigen::SparseMatrix< T, Eigen::RowMajor, I > Sparse
A matrix of a given type.
Sparse< typename T0::Scalar > compress_outer(T0 const &matrix)
Indices of non empty outer indices.