PURIFY
Next-generation radio interferometric imaging
operators.cc
Go to the documentation of this file.
1 #include "purify/operators.h"
2 #include "purify/config.h"
3 #include "purify/types.h"
4 #include "catch2/catch_all.hpp"
5 #include "purify/directories.h"
6 #include "purify/fly_operators.h"
7 #include "purify/kernels.h"
8 #include "purify/logging.h"
9 #include "purify/test_data.h"
10 #include "purify/wproj_operators.h"
11 #include <sopt/power_method.h>
12 
13 using namespace purify;
14 
15 namespace operators_test {
16 const std::string test_dir = "expected/operators/";
18 const std::vector<t_real> u =
19  read_data<t_real>(data_filename(test_dir + "u_data"));
20 const std::vector<t_real> v = read_data<t_real>(data_filename(test_dir + "v_data"));
22 const std::vector<t_complex> direct_input =
23  read_data<t_complex>(data_filename(test_dir + "direct_input_data"));
25 const std::vector<t_complex> expected_direct =
26  read_data<t_complex>(data_filename(test_dir + "expected_direct_data"));
28 const std::vector<t_complex> indirect_input =
29  read_data<t_complex>(data_filename(test_dir + "indirect_input_data"));
31 const std::vector<t_complex> expected_indirect =
32  read_data<t_complex>(data_filename(test_dir + "expected_indirect_data"));
33 
35 const std::vector<t_complex> expected_S =
36  read_data(read_data<t_real>(data_filename(test_dir + "expected_S_data")));
37 } // namespace operators_test
38 
39 TEST_CASE("Operators") {
40  const t_uint M = 10;
41  const t_real oversample_ratio = 2;
42  const t_uint imsizex = 16;
43  const t_uint imsizey = 16;
44  const t_uint ftsizev = std::floor(imsizey * oversample_ratio);
45  const t_uint ftsizeu = std::floor(imsizex * oversample_ratio);
46  const t_uint Ju = 4;
47  const t_uint Jv = 4;
48  const t_uint power_iters = 1e4;
49  const t_real power_tol = 1e-9;
50  const kernels::kernel kernel = kernels::kernel::kb;
51  const std::string &weighting_type = "natural";
52  utilities::vis_params uv_vis;
53  uv_vis.u = Vector<t_real>::Map(operators_test::u.data(), operators_test::u.size());
54  uv_vis.v = Vector<t_real>::Map(operators_test::v.data(), operators_test::v.size());
55  uv_vis.w = uv_vis.u * 0.;
56  uv_vis.weights = Vector<t_complex>::Random(M);
57  uv_vis.vis = Vector<t_complex>::Random(M);
58  uv_vis.units = utilities::vis_units::pixels;
59  std::function<t_real(t_real)> kbu, kbv, ftkbu, ftkbv;
60  std::tie(kbu, kbv, ftkbu, ftkbv) =
61  create_kernels(kernel, Ju, Jv, imsizey, imsizex, oversample_ratio);
62  SECTION("Gridding on the fly") {
63  sopt::OperatorFunction<Vector<t_complex>> directG, indirectG;
64  std::tie(directG, indirectG) = operators::init_on_the_fly_gridding_matrix_2d<Vector<t_complex>>(
65  uv_vis.u, uv_vis.v, Vector<t_complex>::Constant(M, 1.), imsizey, imsizex, oversample_ratio,
66  kbu, Ju, 4e5);
67  Vector<t_complex> direct_output;
68  directG(direct_output,
69  Vector<t_complex>::Map(operators_test::direct_input.data(), ftsizeu * ftsizev));
70  CHECK(direct_output.size() == M);
71  CHECK(direct_output.size() == operators_test::expected_direct.size());
72  CAPTURE(direct_output);
73  CAPTURE(operators_test::expected_direct);
74  REQUIRE(direct_output.isApprox(Vector<t_complex>::Map(operators_test::expected_direct.data(),
75  operators_test::expected_direct.size()),
76  1e-5));
77  Vector<t_complex> indirect_output;
78  indirectG(indirect_output, Vector<t_complex>::Map(operators_test::indirect_input.data(),
79  operators_test::indirect_input.size()));
80  CHECK(indirect_output.size() == ftsizev * ftsizeu);
81  CAPTURE(Vector<t_complex>::Map(operators_test::expected_indirect.data(),
82  operators_test::expected_indirect.size())
83  .head(5));
84  CAPTURE((indirect_output - Vector<t_complex>::Map(operators_test::expected_indirect.data(),
85  operators_test::expected_indirect.size()))
86  .head(5));
87  REQUIRE(
88  indirect_output.isApprox(Vector<t_complex>::Map(operators_test::expected_indirect.data(),
89  operators_test::expected_indirect.size()),
90  1e-5));
91  }
92  SECTION("Gridding") {
93  sopt::OperatorFunction<Vector<t_complex>> directG, indirectG;
94  std::tie(directG, indirectG) = operators::init_gridding_matrix_2d<Vector<t_complex>>(
95  uv_vis.u, uv_vis.v, Vector<t_complex>::Constant(M, 1.), imsizey, imsizex, oversample_ratio,
96  kbv, kbu, Ju, Jv);
97  Vector<t_complex> direct_output;
98  directG(direct_output,
99  Vector<t_complex>::Map(operators_test::direct_input.data(), ftsizeu * ftsizev));
100  CHECK(direct_output.size() == M);
101  CHECK(direct_output.size() == operators_test::expected_direct.size());
102  REQUIRE(direct_output.isApprox(Vector<t_complex>::Map(operators_test::expected_direct.data(),
103  operators_test::expected_direct.size()),
104  1e-5));
105  Vector<t_complex> indirect_output;
106  indirectG(indirect_output, Vector<t_complex>::Map(operators_test::indirect_input.data(),
107  operators_test::indirect_input.size()));
108  CHECK(indirect_output.size() == ftsizev * ftsizeu);
109  CAPTURE(Vector<t_complex>::Map(operators_test::expected_indirect.data(),
110  operators_test::expected_indirect.size())
111  .head(5));
112  CAPTURE((indirect_output - Vector<t_complex>::Map(operators_test::expected_indirect.data(),
113  operators_test::expected_indirect.size()))
114  .head(5));
115  REQUIRE(
116  indirect_output.isApprox(Vector<t_complex>::Map(operators_test::expected_indirect.data(),
117  operators_test::expected_indirect.size()),
118  1e-5));
119  }
120  SECTION("Zero Padding") {
121  const Image<t_complex> S =
122  details::init_correction2d(oversample_ratio, imsizey, imsizex, ftkbu, ftkbv, 0, 0, 0);
123  CHECK(imsizex == S.cols());
124  CHECK(imsizey == S.rows());
125  INFO(S(0) / operators_test::expected_S.at(0));
126  INFO(S(5) / operators_test::expected_S.at(5));
127  REQUIRE(S.isApprox(Image<t_complex>::Map(operators_test::expected_S.data(), imsizey, imsizex),
128  1e-6));
129  sopt::OperatorFunction<Vector<t_complex>> directZ, indirectZ;
130  std::tie(directZ, indirectZ) =
131  operators::init_zero_padding_2d<Vector<t_complex>>(S, oversample_ratio);
132  const Vector<t_complex> direct_input = Vector<t_complex>::Random(imsizex * imsizey);
133  Vector<t_complex> direct_output;
134  directZ(direct_output, direct_input);
135  CHECK(direct_output.size() == ftsizeu * ftsizev);
136  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(ftsizeu * ftsizev);
137  Vector<t_complex> indirect_output;
138  indirectZ(indirect_output, indirect_input);
139  CHECK(indirect_output.size() == imsizex * imsizey);
140  }
141  SECTION("FFT") {
142  sopt::OperatorFunction<Vector<t_complex>> directFFT, indirectFFT;
143  std::tie(directFFT, indirectFFT) =
144  operators::init_FFT_2d<Vector<t_complex>>(imsizey, imsizex, oversample_ratio);
145  const Vector<t_complex> direct_input = Vector<t_complex>::Random(ftsizev * ftsizeu);
146  Vector<t_complex> direct_output;
147  directFFT(direct_output, direct_input);
148  CHECK(direct_output.size() == ftsizeu * ftsizev);
149  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(ftsizev * ftsizeu);
150  Vector<t_complex> indirect_output;
151  indirectFFT(indirect_output, indirect_input);
152  CHECK(indirect_output.size() == ftsizev * ftsizeu);
153  Vector<t_complex> inverse_check;
154  Vector<t_complex> buff;
155  directFFT(buff, direct_input);
156  indirectFFT(inverse_check, buff);
157  CHECK(inverse_check.isApprox(direct_input, 1e-4));
158  }
159  SECTION("Create Weighted Measurement Operator") {
160  const t_uint M_measures = 1e4;
161  const Vector<t_real> u = Vector<t_real>::Random(M_measures);
162  const Vector<t_real> v = Vector<t_real>::Random(M_measures);
163  const Vector<t_real> w = Vector<t_real>::Random(M_measures);
164  const Vector<t_complex> weights = Vector<t_complex>::Random(M_measures);
165  const auto measure_op = measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
166  u, v, w, weights, imsizey, imsizex, oversample_ratio, kernel, Ju, Jv);
167  const Vector<t_complex> direct_input = Vector<t_complex>::Random(imsizex * imsizey);
168  const Vector<t_complex> direct_output = *measure_op * direct_input;
169  CHECK(direct_output.size() == M_measures);
170  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(M_measures).eval();
171  const Vector<t_complex> indirect_output = measure_op->adjoint() * indirect_input;
172  CHECK(indirect_output.size() == imsizex * imsizey);
173  }
174 }
175 TEST_CASE("Degridding from multiple Images") {
176  const t_int imsizex = 32;
177  const t_int imsizey = 32;
178  const t_int M = 5;
179  const t_real oversample_ratio = 2;
180  const t_int ftsizeu = std::floor(imsizex * oversample_ratio);
181  const t_int ftsizev = std::floor(imsizey * oversample_ratio);
182  auto const kernelu = [](t_real) -> t_real { return 1.; };
183  auto const kernelv = [](t_real) -> t_real { return 1.; };
184  const t_int Ju = 1;
185  const t_int Jv = 1;
186  for (t_int images : {1, 2, 3, 5, 10}) {
187  Vector<t_real> u = Vector<t_real>::Zero(M * images);
188  Vector<t_real> v = Vector<t_real>::Zero(M * images);
189  Vector<t_complex> weights = Vector<t_complex>::Zero(M * images);
190  Vector<t_complex> vis = Vector<t_complex>::Zero(M * images);
191  Vector<t_complex> image = Vector<t_complex>::Zero(ftsizeu * ftsizev * images);
192  std::vector<t_int> image_index(M * images, 0);
193  // double up vector
194  u.segment(0, M) = Vector<t_real>::Random(M);
195  v.segment(0, M) = Vector<t_real>::Random(M);
196  weights.segment(0, M) = Vector<t_complex>::Random(M);
197  vis.segment(0, M) = Vector<t_complex>::Random(M);
198  image.segment(0, ftsizeu * ftsizev) = Vector<t_complex>::Random(ftsizeu * ftsizev);
199  for (t_int i = 1; i < images; i++) {
200  image.segment(ftsizeu * ftsizev * i, ftsizeu * ftsizev) =
201  image.segment(0, ftsizeu * ftsizev) * i;
202  u.segment(M * i, M) = u.segment(0, M);
203  v.segment(M * i, M) = v.segment(0, M);
204  weights.segment(M * i, M) = weights.segment(0, M);
205  vis.segment(M * i, M) = vis.segment(0, M) * i;
206  for (t_int c = 0; c < M; c++) image_index[M * i + c] = i;
207  }
208  const auto G_tuple = operators::init_gridding_matrix_2d<Vector<t_complex>>(
209  images, image_index, u, v, weights, imsizey, imsizex, oversample_ratio, kernelu, kernelv,
210  Ju, Jv);
211  sopt::LinearTransform<Vector<t_complex>> const G = {
212  std::get<0>(G_tuple), {0, 1, M}, std::get<1>(G_tuple), {0, 1, ftsizeu * ftsizev * images}};
213  const Vector<t_complex> forward = G * image;
214  for (t_int i = 1; i < images; i++)
215  REQUIRE((i * forward.segment(0, M)).isApprox(forward.segment(M * i, M), 1e-7));
216  const Vector<t_complex> backward = G.adjoint() * vis;
217  for (t_int i = 1; i < images; i++)
218  REQUIRE((i * backward.segment(0, ftsizeu * ftsizev))
219  .isApprox(backward.segment(ftsizeu * ftsizev * i, ftsizeu * ftsizev), 1e-7));
220  }
221 }
222 TEST_CASE("Degridding from multiple Images wproj") {
223  const t_int imsizex = 32;
224  const t_int imsizey = 32;
225  const t_int M = 5;
226  const t_real oversample_ratio = 2;
227  const t_int ftsizeu = std::floor(imsizex * oversample_ratio);
228  const t_int ftsizev = std::floor(imsizey * oversample_ratio);
229  auto const kernelu = [](t_real) -> t_real { return 1.; };
230  auto const kernelv = [](t_real) -> t_real { return 1.; };
231  const t_int Ju = 1;
232  const t_int Jv = 1;
233  for (t_int images : {1, 2, 3, 5, 10}) {
234  Vector<t_real> u = Vector<t_real>::Zero(M * images);
235  Vector<t_real> v = Vector<t_real>::Zero(M * images);
236  Vector<t_complex> weights = Vector<t_complex>::Zero(M * images);
237  Vector<t_complex> vis = Vector<t_complex>::Zero(M * images);
238  Vector<t_complex> image = Vector<t_complex>::Zero(ftsizeu * ftsizev * images);
239  std::vector<t_int> image_index(M * images, 0);
240  // double up vector
241  u.segment(0, M) = Vector<t_real>::Random(M);
242  v.segment(0, M) = Vector<t_real>::Random(M);
243  weights.segment(0, M) = Vector<t_complex>::Random(M);
244  vis.segment(0, M) = Vector<t_complex>::Random(M);
245  image.segment(0, ftsizeu * ftsizev) = Vector<t_complex>::Random(ftsizeu * ftsizev);
246  for (t_int i = 1; i < images; i++) {
247  image.segment(ftsizeu * ftsizev * i, ftsizeu * ftsizev) =
248  image.segment(0, ftsizeu * ftsizev) * i;
249  u.segment(M * i, M) = u.segment(0, M);
250  v.segment(M * i, M) = v.segment(0, M);
251  weights.segment(M * i, M) = weights.segment(0, M);
252  vis.segment(M * i, M) = vis.segment(0, M) * i;
253  for (t_int c = 0; c < M; c++) image_index[M * i + c] = i;
254  }
255  const auto G_tuple = operators::init_gridding_matrix_2d<Vector<t_complex>>(
256  images, image_index, std::vector<t_real>(images, 0), u, v, u * 0, weights, imsizey, imsizex,
257  oversample_ratio, kernelu, kernelv, Ju, 10, 1, 1, 1e-6, 1e-6, dde_type::wkernel_radial);
258  sopt::LinearTransform<Vector<t_complex>> const G = {
259  std::get<0>(G_tuple), {0, 1, M}, std::get<1>(G_tuple), {0, 1, ftsizeu * ftsizev * images}};
260  const Vector<t_complex> forward = G * image;
261  for (t_int i = 1; i < images; i++)
262  REQUIRE((i * forward.segment(0, M)).isApprox(forward.segment(M * i, M), 1e-7));
263  const Vector<t_complex> backward = G.adjoint() * vis;
264  for (t_int i = 1; i < images; i++)
265  REQUIRE((i * backward.segment(0, ftsizeu * ftsizev))
266  .isApprox(backward.segment(ftsizeu * ftsizev * i, ftsizeu * ftsizev), 1e-7));
267  }
268 }
269 TEST_CASE("on the fly with more presamples") {
270  const t_uint M = 1e4;
271  const t_real oversample_ratio = 2;
272  const t_uint imsizex = 16;
273  const t_uint imsizey = 16;
274  const t_uint ftsizev = std::floor(imsizey * oversample_ratio);
275  const t_uint ftsizeu = std::floor(imsizex * oversample_ratio);
276  const t_uint Ju = 4;
277  const t_uint Jv = 4;
278  const t_uint power_iters = 1e4;
279  const t_real power_tol = 1e-9;
280  const kernels::kernel kernel = kernels::kernel::kb;
281  const std::string &weighting_type = "natural";
282  utilities::vis_params uv_vis;
283  uv_vis.u = Vector<t_real>::Random(M) * ftsizeu * 0.5;
284  uv_vis.v = Vector<t_real>::Random(M) * ftsizev * 0.5;
285  uv_vis.w = uv_vis.u * 0.;
286  uv_vis.weights = Vector<t_complex>::Random(M);
287  uv_vis.vis = Vector<t_complex>::Random(M);
288  uv_vis.units = utilities::vis_units::pixels;
289  std::function<t_real(t_real)> kbu, kbv, ftkbu, ftkbv;
290  std::tie(kbu, kbv, ftkbu, ftkbv) =
291  create_kernels(kernel, Ju, Jv, imsizey, imsizex, oversample_ratio);
292  SECTION("Gridding on the fly") {
293  sopt::OperatorFunction<Vector<t_complex>> directG, indirectG;
294  std::tie(directG, indirectG) = operators::init_gridding_matrix_2d<Vector<t_complex>>(
295  uv_vis.u, uv_vis.v, Vector<t_complex>::Constant(M, 1.), imsizey, imsizex, oversample_ratio,
296  kbv, kbu, Ju, Jv);
297  sopt::OperatorFunction<Vector<t_complex>> flydirectG, flyindirectG;
298  std::tie(flydirectG, flyindirectG) =
299  operators::init_on_the_fly_gridding_matrix_2d<Vector<t_complex>>(
300  uv_vis.u, uv_vis.v, Vector<t_complex>::Constant(M, 1.), imsizey, imsizex,
301  oversample_ratio, kbu, Ju, 4e5);
302  SECTION("direct") {
303  Vector<t_complex> direct_output;
304  Vector<t_complex> flydirect_output;
305  const Vector<t_complex> direct_input = Vector<t_complex>::Random(ftsizev * ftsizeu);
306  directG(direct_output, direct_input);
307  flydirectG(flydirect_output, direct_input);
308  CHECK(direct_output.size() == M);
309  CHECK(direct_output.size() == flydirect_output.size());
310  CAPTURE(direct_output);
311  CAPTURE(flydirect_output);
312  REQUIRE(direct_output.isApprox(flydirect_output, 1e-5));
313  }
314  SECTION("indirect") {
315  Vector<t_complex> indirect_output;
316  Vector<t_complex> flyindirect_output;
317  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(M);
318  indirectG(indirect_output, indirect_input);
319  flyindirectG(flyindirect_output, indirect_input);
320  CHECK(indirect_output.size() == ftsizev * ftsizeu);
321  CAPTURE(flyindirect_output.head(5));
322  CAPTURE((indirect_output - flyindirect_output).head(5));
323  REQUIRE(indirect_output.isApprox(flyindirect_output, 1e-5));
324  }
325  }
326 }
CHECK
#define CHECK(CONDITION, ERROR)
Definition: casa.cc:6
fly_operators.h
operators_test
Definition: operators.cc:15
operators_test::test_dir
const std::string test_dir
Definition: operators.cc:16
operators_test::expected_indirect
const std::vector< t_complex > expected_indirect
data for gridding output
Definition: operators.cc:31
operators_test::indirect_input
const std::vector< t_complex > indirect_input
data for gridding input
Definition: operators.cc:28
operators_test::u
const std::vector< t_real > u
data for u coordinate
Definition: operators.cc:18
operators_test::v
const std::vector< t_real > v
data for v coordinate
Definition: operators.cc:20
operators_test::expected_direct
const std::vector< t_complex > expected_direct
data for degridding output
Definition: operators.cc:25
operators_test::expected_S
const std::vector< t_complex > expected_S
data for gridding correction
Definition: operators.cc:35
operators_test::direct_input
const std::vector< t_complex > direct_input
data for degridding input
Definition: operators.cc:22
purify::constant::c
const t_real c
speed of light in vacuum
Definition: types.h:72
purify::details::init_correction2d
Image< t_complex > init_correction2d(const t_real &oversample_ratio, const t_uint &imsizey_, const t_uint &imsizex_, const std::function< t_real(t_real)> ftkernelu, const std::function< t_real(t_real)> ftkernelv, const t_real &w_mean, const t_real &cellx, const t_real &celly)
Definition: operators.cc:47
purify::data_filename
std::string data_filename(std::string const &filename)
Holds data and such.
Definition: directories.in.h:16
purify::create_kernels
std::tuple< std::function< t_real(t_real)>, std::function< t_real(t_real)>, std::function< t_real(t_real)>, std::function< t_real(t_real)> > create_kernels(const kernels::kernel kernel_name_, const t_uint Ju_, const t_uint Jv_, const t_real imsizey_, const t_real imsizex_, const t_real oversample_ratio)
Definition: kernels.cc:249
purify::read_data
std::enable_if< std::is_scalar< T >::value, std::vector< T > >::type read_data(const std::string &filename)
read real values from data file
Definition: data.in.h:14
purify::utilities::vis_params
Definition: uvw_utilities.h:15
purify::utilities::vis_params::vis
Vector< t_complex > vis
Definition: uvw_utilities.h:22
purify::utilities::vis_params::v
Vector< t_real > v
Definition: uvw_utilities.h:17
purify::utilities::vis_params::u
Vector< t_real > u
Definition: uvw_utilities.h:16
purify::utilities::vis_params::w
Vector< t_real > w
Definition: uvw_utilities.h:18
purify::utilities::vis_params::weights
Vector< t_complex > weights
Definition: uvw_utilities.h:23
purify::utilities::vis_params::units
vis_units units
Definition: uvw_utilities.h:24
TEST_CASE
TEST_CASE("Operators")
Definition: operators.cc:39
wproj_operators.h
Generated by doxygen 1.9.1