PURIFY
Next-generation radio interferometric imaging
index_mapping.cc
Go to the documentation of this file.
1 #include "purify/config.h"
2 #include "purify/types.h"
3 #include <numeric>
4 #include <random>
5 #include "catch2/catch_all.hpp"
6 #include "purify/IndexMapping.h"
7 #include "purify/utilities.h"
8 
9 using namespace purify;
10 
11 TEST_CASE("Index mapping operator") {
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);
16 
17  auto mapper = IndexMapping<t_int>(indices, input.size());
18  Vector<t_int> output;
19 
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);
25  }
26 
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);
33  adjoint(index) = 0;
34  }
35  CHECK(adjoint == adjoint.Zero(adjoint.size()));
36  }
37  }
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);
42 
43  auto mapper = IndexMapping<t_int>(indices, input.size());
44  Vector<t_int> output;
45 
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);
51  }
52 
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));
59  adjoint(index) = 0;
60  }
61  CHECK(adjoint == adjoint.Zero(adjoint.size()));
62  }
63  }
64 }
65 
66 TEST_CASE("Non empty outer vectors") {
67  Sparse<t_int> matrix(4, 4);
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());
70 
71  auto const indices = non_empty_outers(matrix);
72  CHECK(indices.size() == 2);
73  CHECK(indices.count(0) == 1);
74  CHECK(indices.count(3) == 1);
75 }
76 
77 TEST_CASE("Shrink sparse matrix") {
78  Sparse<t_int> matrix(4, 4);
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());
81  Sparse<t_int> shrinked_matrix = compress_outer(matrix);
82 
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));
87 }
88 
89 TEST_CASE("Vector-shrinked matrix multiplication") {
90  auto const N = 10;
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))
94  .matrix()
95  .select(Matrix<uint8_t>::Random(N, N), Matrix<uint8_t>::Zero(N, N))
96  .cast<t_int>();
97 
98  matrix.row(uniform_dist(*mersenne)).fill(0);
99  matrix.row(uniform_dist(*mersenne)).fill(0);
100  Sparse<t_int> const sparse = matrix.sparseView();
101  Vector<t_int> const input = Vector<uint8_t>::Random(N).cast<t_int>();
102 
103  Sparse<t_int> const compressed = compress_outer(sparse);
104  auto const mapper = IndexMapping<t_int>(non_empty_outers(sparse), input.size());
105  Vector<t_int> comp_input;
106  mapper(input, comp_input);
107 
108  CHECK(sparse * input == compressed * comp_input);
109  CHECK(utilities::sparse_multiply_matrix(sparse, input) == compressed * comp_input);
110  CHECK(utilities::sparse_multiply_matrix(compressed, comp_input) == compressed * comp_input);
111 }
#define CHECK(CONDITION, ERROR)
Definition: casa.cc:6
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.
Definition: utilities.h:67
std::set< STORAGE_INDEX_TYPE > non_empty_outers(Eigen::SparseMatrixBase< T0 > const &matrix)
Indices of non empty outer indices.
Definition: IndexMapping.h:61
Eigen::SparseMatrix< T, Eigen::RowMajor, I > Sparse
A matrix of a given type.
Definition: types.h:40
Sparse< typename T0::Scalar > compress_outer(T0 const &matrix)
Indices of non empty outer indices.
Definition: IndexMapping.h:71