PURIFY
Next-generation radio interferometric imaging
mpi_measurement_operator.cc
Go to the documentation of this file.
1 #include <iomanip>
2 #include <random>
3 #include "catch2/catch_all.hpp"
4 #include "purify/distribute.h"
5 #include "purify/logging.h"
6 #include "purify/mpi_utilities.h"
7 #include "purify/operators.h"
8 #include "purify/utilities.h"
9 #include "purify/wproj_operators.h"
10 #include <sopt/mpi/communicator.h>
11 #include <sopt/power_method.h>
12 #ifdef PURIFY_ARRAYFIRE
13 #include "purify/operators_gpu.h"
14 #include "purify/wproj_operators_gpu.h"
15 #endif
16 using namespace purify;
17 
18 TEST_CASE("Serial vs Distributed Operator") {
19  auto const world = sopt::mpi::Communicator::World();
20 
21  auto const N = 1000;
22  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
23  uv_serial.u = world.broadcast(uv_serial.u);
24  uv_serial.v = world.broadcast(uv_serial.v);
25  uv_serial.w = world.broadcast(uv_serial.w);
26  uv_serial.units = utilities::vis_units::radians;
27  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
28  uv_serial.weights =
29  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
30 
31  utilities::vis_params uv_mpi;
32  if (world.is_root()) {
33  auto const order =
34  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
35  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
36  } else
37  uv_mpi = utilities::scatter_visibilities(world);
38 
39  auto const over_sample = 2;
40  auto const J = 4;
41  auto const kernel = kernels::kernel::kb;
42  auto const width = 128;
43  auto const height = 128;
44  const Vector<t_complex> power_init =
45  world.broadcast(Vector<t_complex>::Random(height * width).eval());
46  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
47  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
48  uv_serial.u, uv_serial.v, uv_serial.w, uv_serial.weights, height, width, over_sample),
49  100, 1e-4, power_init));
50 
51  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
52  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
53  world, uv_mpi.u, uv_mpi.v, uv_mpi.w, uv_mpi.weights, height, width, over_sample),
54  100, 1e-4, power_init));
55 
56  if (uv_serial.u.size() == uv_mpi.u.size()) {
57  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
58  CHECK(uv_serial.v.isApprox(uv_mpi.v));
59  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
60  }
61  SECTION("Degridding") {
62  Vector<t_complex> const image =
63  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
64 
65  auto uv_degrid = uv_serial;
66  if (world.is_root()) {
67  uv_degrid.vis = *op_serial * image;
68  auto const order =
69  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
70  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
71  } else
72  uv_degrid = utilities::scatter_visibilities(world);
73  Vector<t_complex> const degridded = *op * image;
74  REQUIRE(degridded.size() == uv_degrid.vis.size());
75  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
76  }
77  SECTION("Gridding") {
78  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
79  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
80  REQUIRE(gridded.size() == gridded_serial.size());
81  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
82  }
83 }
84 TEST_CASE("Serial vs Distributed Fourier Grid Operator") {
85  auto const world = sopt::mpi::Communicator::World();
86 
87  auto const N = 1000;
88  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
89  uv_serial.u = world.broadcast(uv_serial.u);
90  uv_serial.v = world.broadcast(uv_serial.v);
91  uv_serial.w = world.broadcast(uv_serial.w);
92  uv_serial.units = utilities::vis_units::radians;
93  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
94  uv_serial.weights =
95  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
96 
97  utilities::vis_params uv_mpi;
98  if (world.is_root()) {
99  auto const order =
100  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
101  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
102  } else
103  uv_mpi = utilities::scatter_visibilities(world);
104 
105  auto const over_sample = 2;
106  auto const J = 4;
107  auto const kernel = kernels::kernel::kb;
108  auto const width = 128;
109  auto const height = 128;
110  const Vector<t_complex> power_init =
111  world.broadcast(Vector<t_complex>::Random(height * width).eval());
112  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
113  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
114  uv_serial.u, uv_serial.v, uv_serial.w, uv_serial.weights, height, width, over_sample),
115  100, 1e-4, power_init));
116 
117  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
118  purify::measurementoperator::init_degrid_operator_2d_mpi<Vector<t_complex>>(
119  world, uv_mpi.u, uv_mpi.v, uv_mpi.w, uv_mpi.weights, height, width, over_sample),
120  100, 1e-4, power_init));
121 
122  if (uv_serial.u.size() == uv_mpi.u.size()) {
123  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
124  CHECK(uv_serial.v.isApprox(uv_mpi.v));
125  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
126  }
127  SECTION("Degridding") {
128  Vector<t_complex> const image =
129  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
130 
131  auto uv_degrid = uv_serial;
132  if (world.is_root()) {
133  uv_degrid.vis = *op_serial * image;
134  auto const order =
135  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
136  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
137  } else
138  uv_degrid = utilities::scatter_visibilities(world);
139  Vector<t_complex> const degridded = *op * image;
140  REQUIRE(degridded.size() == uv_degrid.vis.size());
141  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
142  }
143  SECTION("Gridding") {
144  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
145  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
146  REQUIRE(gridded.size() == gridded_serial.size());
147  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
148  }
149 }
150 
151 TEST_CASE("Serial vs Distributed Fourier Grid Operator weighted") {
152  // sopt::logging::set_level("debug");
153  // purify::logging::set_level("debug");
154  auto const world = sopt::mpi::Communicator::World();
155 
156  auto const N = 1000;
157  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
158  uv_serial.u = world.broadcast(uv_serial.u);
159  uv_serial.v = world.broadcast(uv_serial.v);
160  uv_serial.w = world.broadcast(uv_serial.w);
161  uv_serial.units = utilities::vis_units::radians;
162  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
163  uv_serial.weights =
164  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
165 
166  utilities::vis_params uv_mpi;
167  if (world.is_root()) {
168  auto const order =
169  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
170  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
171  } else
172  uv_mpi = utilities::scatter_visibilities(world);
173 
174  auto const over_sample = 2;
175  auto const J = 4;
176  auto const kernel = kernels::kernel::kb;
177  auto const width = 128;
178  auto const height = 128;
179  const Vector<t_complex> power_init =
180  world.broadcast(Vector<t_complex>::Random(height * width).eval());
181  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
182  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
183  uv_serial.u, uv_serial.v, uv_serial.w, uv_serial.weights, height, width, over_sample),
184  100, 1e-4, power_init));
185  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
186  purify::measurementoperator::init_degrid_operator_2d_mpi<Vector<t_complex>>(
187  world, uv_mpi.u, uv_mpi.v, uv_mpi.w, uv_mpi.weights, height, width, over_sample),
188  100, 1e-4, power_init));
189 
190  if (world.size() == 1) {
191  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
192  CHECK(uv_serial.v.isApprox(uv_mpi.v));
193  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
194  }
195  SECTION("Degridding") {
196  Vector<t_complex> const image =
197  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
198 
199  auto uv_degrid = uv_serial;
200  if (world.is_root()) {
201  uv_degrid.vis = *op_serial * image;
202  auto const order =
203  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
204  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
205  } else
206  uv_degrid = utilities::scatter_visibilities(world);
207  Vector<t_complex> const degridded = *op * image;
208  REQUIRE(degridded.size() == uv_degrid.vis.size());
209  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
210  }
211  SECTION("Gridding") {
212  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
213  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
214  REQUIRE(gridded.size() == gridded_serial.size());
215  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
216  }
217 }
218 
219 TEST_CASE("Serial vs All to All Fourier Grid Operator weighted") {
220  // sopt::logging::set_level("debug");
221  // purify::logging::set_level("debug");
222  auto const world = sopt::mpi::Communicator::World();
223 
224  auto const N = 1000;
225  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
226  uv_serial.u = world.broadcast(uv_serial.u);
227  uv_serial.v = world.broadcast(uv_serial.v);
228  uv_serial.w = world.broadcast(uv_serial.w);
229  uv_serial.units = utilities::vis_units::radians;
230  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
231  uv_serial.weights =
232  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
233 
234  utilities::vis_params uv_mpi;
235  if (world.is_root()) {
236  auto const order =
237  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
238  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
239  } else
240  uv_mpi = utilities::scatter_visibilities(world);
241 
242  auto const over_sample = 2;
243  auto const J = 4;
244  auto const kernel = kernels::kernel::kb;
245  auto const width = 128;
246  auto const height = 128;
247  const Vector<t_complex> power_init =
248  world.broadcast(Vector<t_complex>::Random(height * width).eval());
249  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
250  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
251  uv_serial.u, uv_serial.v, uv_serial.w, uv_serial.weights, height, width, over_sample),
252  100, 1e-4, power_init));
253  // First create an instance of an engine.
254  std::random_device rnd_device;
255  // Specify the engine and distribution.
256  std::mt19937 mersenne_engine(rnd_device()); // Generates random integers
257  std::uniform_int_distribution<t_int> dist(0, world.size() - 1);
258 
259  auto gen = [&dist, &mersenne_engine]() { return dist(mersenne_engine); };
260  std::vector<t_int> image_index(uv_mpi.size());
261  std::generate(image_index.begin(), image_index.end(), gen);
262  std::vector<t_real> w_stacks(world.size(), 0.);
263 
264  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
265  purify::measurementoperator::init_degrid_operator_2d_all_to_all<Vector<t_complex>>(
266  world, image_index, w_stacks, uv_mpi.u, uv_mpi.v, uv_mpi.w, uv_mpi.weights, height, width,
267  over_sample),
268  100, 1e-4, power_init));
269 
270  if (world.size() == 1) {
271  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
272  CHECK(uv_serial.v.isApprox(uv_mpi.v));
273  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
274  }
275  SECTION("Degridding") {
276  Vector<t_complex> const image =
277  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
278 
279  auto uv_degrid = uv_serial;
280  if (world.is_root()) {
281  uv_degrid.vis = *op_serial * image;
282  auto const order =
283  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
284  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
285  } else
286  uv_degrid = utilities::scatter_visibilities(world);
287  Vector<t_complex> const degridded = *op * image;
288  REQUIRE(degridded.size() == uv_degrid.vis.size());
289  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
290  }
291  SECTION("Gridding") {
292  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
293  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
294  REQUIRE(gridded.size() == gridded_serial.size());
295  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
296  }
297 }
298 
299 TEST_CASE("Standard vs All to All stacking") {
300  // sopt::logging::set_level("debug");
301  // purify::logging::set_level("debug");
302  auto const world = sopt::mpi::Communicator::World();
303 
304  auto const N = 1000;
305  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3, 100);
306  uv_serial.u = world.broadcast(uv_serial.u);
307  uv_serial.v = world.broadcast(uv_serial.v);
308  uv_serial.w = world.broadcast(uv_serial.w);
309  uv_serial.units = utilities::vis_units::radians;
310  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
311  uv_serial.weights =
312  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
313 
314  utilities::vis_params uv_mpi;
315  if (world.is_root()) {
316  auto const order =
317  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
318  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
319  } else
320  uv_mpi = utilities::scatter_visibilities(world);
321 
322  auto const over_sample = 2;
323  auto const J = 4;
324  auto const kernel = kernels::kernel::kb;
325  auto const width = 128;
326  auto const height = 128;
327  auto const cell_size = 1;
328  // First create an instance of an engine.
329  const auto kmeans = distribute::kmeans_algo(uv_mpi.w, world.size(), 100, world);
330  const std::vector<t_int> image_index = std::get<0>(kmeans);
331  const std::vector<t_real> w_stacks = std::get<1>(kmeans);
332  const Vector<t_complex> power_init =
333  world.broadcast(Vector<t_complex>::Random(height * width).eval());
334 
335  const auto uv_stacks = utilities::regroup_and_all_to_all(uv_mpi, image_index, world);
336  // standard operator
337  const auto op_stack = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
338  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
339  world, uv_stacks, height, width, cell_size, cell_size, over_sample, kernel, J, J, true),
340  100, 1e-4, power_init));
341  // all to all operator
342  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
343  purify::measurementoperator::init_degrid_operator_2d_all_to_all<Vector<t_complex>>(
344  world, image_index, w_stacks, uv_mpi, height, width, cell_size, cell_size, over_sample,
345  kernel, J, J, true),
346  100, 1e-4, power_init));
347  if (world.size() == 1) {
348  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
349  CHECK(uv_serial.v.isApprox(uv_mpi.v));
350  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
351  }
352  SECTION("Degridding") {
353  Vector<t_complex> const image =
354  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
355 
356  const Vector<t_complex> degridded = *op_stack * image;
357  auto uv_degrid = uv_mpi;
358  uv_degrid.vis = *op * image;
359  uv_degrid = utilities::regroup_and_all_to_all(uv_degrid, image_index, world);
360  REQUIRE(degridded.size() == uv_degrid.vis.size());
361  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
362  }
363  SECTION("Gridding") {
364  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
365  Vector<t_complex> const gridded_serial = op_stack->adjoint() * uv_stacks.vis;
366  REQUIRE(gridded.size() == gridded_serial.size());
367  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
368  }
369 }
370 TEST_CASE("Standard vs All to All wproj") {
371  // sopt::logging::set_level("debug");
372  // purify::logging::set_level("debug");
373  auto const world = sopt::mpi::Communicator::World();
374 
375  auto const N = 1000;
376  const t_real w_rms = 100;
377  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3, w_rms);
378  uv_serial.u = world.broadcast(uv_serial.u);
379  uv_serial.v = world.broadcast(uv_serial.v);
380  uv_serial.w = world.broadcast(uv_serial.w);
381  uv_serial.units = utilities::vis_units::radians;
382  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
383  uv_serial.weights =
384  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
385 
386  utilities::vis_params uv_mpi;
387  if (world.is_root()) {
388  auto const order =
389  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
390  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
391  } else
392  uv_mpi = utilities::scatter_visibilities(world);
393 
394  auto const over_sample = 2;
395  auto const J = 4;
396  auto const kernel = kernels::kernel::kb;
397  auto const width = 128;
398  auto const height = 128;
399  auto const cell_size = 60;
400  // First create an instance of an engine.
401  const auto kmeans = distribute::kmeans_algo(uv_mpi.w, world.size(), 100, world);
402  const std::vector<t_int> image_index = std::get<0>(kmeans);
403  const std::vector<t_real> w_stacks = std::get<1>(kmeans);
404  const Vector<t_complex> power_init =
405  world.broadcast(Vector<t_complex>::Random(height * width).eval());
406 
407  const auto uv_stacks = utilities::regroup_and_all_to_all(uv_mpi, image_index, world);
408  // standard operator
409  const auto op_wproj = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
410  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
411  world, uv_stacks, height, width, cell_size, cell_size, over_sample, kernel, J, 100, true,
412  1e-8, 1e-8, dde_type::wkernel_radial),
413  100, 1e-4, power_init));
414  // all to all operator
415  const auto op_wproj_all = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
416  purify::measurementoperator::init_degrid_operator_2d_all_to_all<Vector<t_complex>>(
417  world, image_index, w_stacks, uv_mpi, height, width, cell_size, cell_size, over_sample,
418  kernel, J, 100, true, 1e-8, 1e-8, dde_type::wkernel_radial),
419  100, 1e-4, power_init));
420  if (world.size() == 1) {
421  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
422  CHECK(uv_serial.v.isApprox(uv_mpi.v));
423  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
424  }
425  SECTION("Degridding") {
426  Vector<t_complex> const image =
427  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
428 
429  const Vector<t_complex> degridded = *op_wproj * image;
430  auto uv_degrid = uv_mpi;
431  uv_degrid.vis = *op_wproj_all * image;
432  uv_degrid = utilities::regroup_and_all_to_all(uv_degrid, image_index, world);
433  REQUIRE(degridded.size() == uv_degrid.vis.size());
434  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
435  }
436  SECTION("Gridding") {
437  Vector<t_complex> const gridded = op_wproj_all->adjoint() * uv_mpi.vis;
438  Vector<t_complex> const gridded_serial = op_wproj->adjoint() * uv_stacks.vis;
439  REQUIRE(gridded.size() == gridded_serial.size());
440  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
441  }
442 }
443 #ifdef PURIFY_ARRAYFIRE
444 TEST_CASE("Serial vs Distributed GPU Fourier Grid Operator weighted") {
445  // sopt::logging::set_level("debug");
446  // purify::logging::set_level("debug");
447  af::setDevice(0);
448  auto const world = sopt::mpi::Communicator::World();
449  auto const N = 1000;
450  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
451  uv_serial.u = world.broadcast(uv_serial.u);
452  uv_serial.v = world.broadcast(uv_serial.v);
453  uv_serial.w = world.broadcast(uv_serial.w);
454  uv_serial.units = utilities::vis_units::radians;
455  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
456  uv_serial.weights =
457  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
458 
459  utilities::vis_params uv_mpi;
460  if (world.is_root()) {
461  auto const order =
462  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
463  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
464  } else
465  uv_mpi = utilities::scatter_visibilities(world);
466 
467  auto const over_sample = 2;
468  auto const J = 4;
469  auto const kernel = kernels::kernel::kb;
470  auto const width = 128;
471  auto const height = 128;
472  const Vector<t_complex> power_init =
473  world.broadcast(Vector<t_complex>::Random(height * width).eval());
474  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
475  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
476  uv_serial.u, uv_serial.v, uv_serial.w, uv_serial.weights, height, width, over_sample),
477  100, 1e-4, power_init));
478  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
479  purify::gpu::measurementoperator::init_degrid_operator_2d_mpi(
480  world, uv_mpi.u, uv_mpi.v, uv_mpi.w, uv_mpi.weights, height, width, over_sample),
481  100, 1e-4, power_init));
482 
483  if (world.size() == 1) {
484  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
485  CHECK(uv_serial.v.isApprox(uv_mpi.v));
486  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
487  }
488  SECTION("Power Method") {
489  auto op_norm = std::get<0>(sopt::algorithm::power_method<Vector<t_complex>>(
490  *op, 100, 1e-4, Vector<t_complex>::Random(width * height)));
491  CHECK(std::abs(op_norm - 1.) < 1e-4);
492  auto op_norm_old = std::get<0>(sopt::algorithm::power_method<Vector<t_complex>>(
493  *op_serial, 100, 1e-4, Vector<t_complex>::Random(width * height)));
494  CHECK(std::abs(op_norm_old - 1.) < 1e-4);
495  }
496  SECTION("Degridding") {
497  Vector<t_complex> const image =
498  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
499 
500  auto uv_degrid = uv_serial;
501  if (world.is_root()) {
502  uv_degrid.vis = *op_serial * image;
503  auto const order =
504  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
505  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
506  } else
507  uv_degrid = utilities::scatter_visibilities(world);
508  Vector<t_complex> const degridded = *op * image;
509  REQUIRE(degridded.size() == uv_degrid.vis.size());
510  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
511  }
512  SECTION("Gridding") {
513  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
514  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
515  REQUIRE(gridded.size() == gridded_serial.size());
516  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
517  }
518 }
519 TEST_CASE("Serial vs Distributed GPU Operator weighted") {
520  // sopt::logging::set_level("debug");
521  // purify::logging::set_level("debug");
522  af::setDevice(0);
523  auto const world = sopt::mpi::Communicator::World();
524  auto const N = 1000;
525  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
526  uv_serial.u = world.broadcast(uv_serial.u);
527  uv_serial.v = world.broadcast(uv_serial.v);
528  uv_serial.w = world.broadcast(uv_serial.w);
529  uv_serial.units = utilities::vis_units::radians;
530  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
531  uv_serial.weights =
532  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
533 
534  utilities::vis_params uv_mpi;
535  if (world.is_root()) {
536  auto const order =
537  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
538  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
539  } else
540  uv_mpi = utilities::scatter_visibilities(world);
541 
542  auto const over_sample = 2;
543  auto const J = 4;
544  auto const kernel = kernels::kernel::kb;
545  auto const width = 128;
546  auto const height = 128;
547  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
548  purify::measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
549  uv_serial.u, uv_serial.v, uv_serial.w, uv_serial.weights, height, width, over_sample),
550  100, 1e-4, Vector<t_complex>::Random(height * width)));
551  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
552  purify::gpu::measurementoperator::init_degrid_operator_2d(
553  world, uv_mpi.u, uv_mpi.v, uv_mpi.w, uv_mpi.weights, height, width, over_sample),
554  100, 1e-4, world.broadcast(Vector<t_complex>::Random(height * width).eval())));
555 
556  if (world.size() == 1) {
557  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
558  CHECK(uv_serial.v.isApprox(uv_mpi.v));
559  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
560  }
561  SECTION("Power Method") {
562  auto op_norm = std::get<0>(sopt::algorithm::power_method<Vector<t_complex>>(
563  *op, 100, 1e-4, Vector<t_complex>::Random(width * height)));
564  CHECK(std::abs(op_norm - 1.) < 1e-4);
565  auto op_norm_old = std::get<0>(sopt::algorithm::power_method<Vector<t_complex>>(
566  *op_serial, 100, 1e-4, Vector<t_complex>::Random(width * height)));
567  CHECK(std::abs(op_norm_old - 1.) < 1e-4);
568  }
569  SECTION("Degridding") {
570  Vector<t_complex> const image =
571  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
572 
573  auto uv_degrid = uv_serial;
574  if (world.is_root()) {
575  uv_degrid.vis = *op_serial * image;
576  auto const order =
577  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
578  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
579  } else
580  uv_degrid = utilities::scatter_visibilities(world);
581  Vector<t_complex> const degridded = *op * image;
582  REQUIRE(degridded.size() == uv_degrid.vis.size());
583  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
584  }
585  SECTION("Gridding") {
586  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
587  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
588  REQUIRE(gridded.size() == gridded_serial.size());
589  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
590  }
591 }
592 TEST_CASE("Serial vs Distributed GPU Operator Radial WProjection") {
593  auto const world = sopt::mpi::Communicator::World();
594 
595  auto const N = 1000;
596  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
597  uv_serial.u = world.broadcast(uv_serial.u);
598  uv_serial.v = world.broadcast(uv_serial.v);
599  uv_serial.w = world.broadcast(Vector<t_real>::Zero(uv_serial.size()).eval());
600  uv_serial.units = utilities::vis_units::radians;
601  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
602  uv_serial.weights =
603  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
604 
605  utilities::vis_params uv_mpi;
606  if (world.is_root()) {
607  auto const order =
608  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
609  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
610  } else
611  uv_mpi = utilities::scatter_visibilities(world);
612 
613  auto const over_sample = 2;
614  auto const J = 4;
615  auto const kernel = kernels::kernel::kb;
616  auto const width = 128;
617  auto const height = 128;
618  const t_real cellx = 1;
619  const t_real celly = 1;
620  const t_real abs_error = 1e-6;
621  const t_real rel_error = 1e-8;
622  const Vector<t_complex> power_init =
623  world.broadcast(Vector<t_complex>::Random(height * width).eval());
624  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
625  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
626  uv_serial, height, width, cellx, celly, over_sample, kernel, J, 4, true, abs_error,
627  rel_error, dde_type::wkernel_radial),
628  10000, 1e-5, power_init);
629 
630  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
631  gpu::measurementoperator::init_degrid_operator_2d(world, uv_mpi, height, width, cellx, celly,
632  over_sample, kernel, J, 4, true, abs_error,
633  rel_error, dde_type::wkernel_radial),
634  10000, 1e-5, power_init);
635 
636  if (uv_serial.u.size() == uv_mpi.u.size()) {
637  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
638  CHECK(uv_serial.v.isApprox(uv_mpi.v));
639  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
640  }
641  SECTION("Degridding") {
642  Vector<t_complex> const image =
643  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
644 
645  auto uv_degrid = uv_serial;
646  if (world.is_root()) {
647  uv_degrid.vis = *op_serial * image;
648  auto const order =
649  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
650  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
651  } else
652  uv_degrid = utilities::scatter_visibilities(world);
653  Vector<t_complex> const degridded = *op * image;
654  REQUIRE(degridded.size() == uv_degrid.vis.size());
655  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
656  }
657  SECTION("Gridding") {
658  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
659  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
660  REQUIRE(gridded.size() == gridded_serial.size());
661  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
662  }
663 }
664 #endif
665 
666 TEST_CASE("Serial vs Distributed Operator Radial WProjection") {
667  auto const world = sopt::mpi::Communicator::World();
668 
669  auto const N = 1000;
670  auto uv_serial = utilities::random_sample_density(N, 0, constant::pi / 3);
671  uv_serial.u = world.broadcast(uv_serial.u);
672  uv_serial.v = world.broadcast(uv_serial.v);
673  uv_serial.w = world.broadcast(Vector<t_real>::Zero(uv_serial.size()).eval());
674  uv_serial.units = utilities::vis_units::radians;
675  uv_serial.vis = world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
676  uv_serial.weights =
677  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(uv_serial.u.size()));
678 
679  utilities::vis_params uv_mpi;
680  if (world.is_root()) {
681  auto const order =
682  distribute::distribute_measurements(uv_serial, world, distribute::plan::radial);
683  uv_mpi = utilities::regroup_and_scatter(uv_serial, order, world);
684  } else
685  uv_mpi = utilities::scatter_visibilities(world);
686 
687  auto const over_sample = 2;
688  auto const J = 4;
689  auto const kernel = kernels::kernel::kb;
690  auto const width = 128;
691  auto const height = 128;
692  const t_real cellx = 1;
693  const t_real celly = 1;
694  const t_real abs_error = 1e-8;
695  const t_real rel_error = 1e-8;
696  const Vector<t_complex> power_init =
697  world.broadcast(Vector<t_complex>::Random(height * width).eval());
698  const auto op_serial = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
699  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
700  uv_serial, height, width, cellx, celly, over_sample, kernel, J, 4, true, abs_error,
701  rel_error, dde_type::wkernel_radial),
702  10000, 1e-5, power_init));
703 
704  const auto op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
705  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
706  world, uv_mpi, height, width, cellx, celly, over_sample, kernel, J, 4, true, abs_error,
707  rel_error, dde_type::wkernel_radial),
708  10000, 1e-5, power_init));
709 
710  if (uv_serial.u.size() == uv_mpi.u.size()) {
711  REQUIRE(uv_serial.u.isApprox(uv_mpi.u));
712  CHECK(uv_serial.v.isApprox(uv_mpi.v));
713  CHECK(uv_serial.weights.isApprox(uv_mpi.weights));
714  }
715  SECTION("Degridding") {
716  Vector<t_complex> const image =
717  world.broadcast<Vector<t_complex>>(Vector<t_complex>::Random(width * height));
718 
719  auto uv_degrid = uv_serial;
720  if (world.is_root()) {
721  uv_degrid.vis = *op_serial * image;
722  auto const order =
723  distribute::distribute_measurements(uv_degrid, world, distribute::plan::radial);
724  uv_degrid = utilities::regroup_and_scatter(uv_degrid, order, world);
725  } else
726  uv_degrid = utilities::scatter_visibilities(world);
727  Vector<t_complex> const degridded = *op * image;
728  REQUIRE(degridded.size() == uv_degrid.vis.size());
729  REQUIRE(degridded.isApprox(uv_degrid.vis, 1e-4));
730  }
731  SECTION("Gridding") {
732  Vector<t_complex> const gridded = op->adjoint() * uv_mpi.vis;
733  Vector<t_complex> const gridded_serial = op_serial->adjoint() * uv_serial.vis;
734  REQUIRE(gridded.size() == gridded_serial.size());
735  REQUIRE(gridded.isApprox(gridded_serial, 1e-4));
736  }
737 }
CHECK
#define CHECK(CONDITION, ERROR)
Definition: casa.cc:6
TEST_CASE
TEST_CASE("Serial vs Distributed Operator")
Definition: mpi_measurement_operator.cc:18
mpi_utilities.h
purify::constant::pi
const t_real pi
mathematical constant
Definition: types.h:70
purify::distribute::kmeans_algo
std::tuple< std::vector< t_int >, std::vector< t_real > > kmeans_algo(const Vector< t_real > &w, const t_int number_of_nodes, const t_int iters, const std::function< t_real(t_real)> &cost, const t_real rel_diff)
patition w terms using k-means
Definition: distribute.cc:103
purify::distribute::distribute_measurements
std::vector< t_int > distribute_measurements(Vector< t_real > const &u, Vector< t_real > const &v, Vector< t_real > const &w, t_int const number_of_nodes, distribute::plan const distribution_plan, t_int const &grid_size)
Distribute visiblities into groups.
Definition: distribute.cc:6
purify::measurementoperator::init_degrid_operator_2d
std::shared_ptr< sopt::LinearTransform< T > > init_degrid_operator_2d(const Vector< t_real > &u, const Vector< t_real > &v, const Vector< t_real > &w, const Vector< t_complex > &weights, const t_uint &imsizey, const t_uint &imsizex, const t_real &oversample_ratio=2, const kernels::kernel kernel=kernels::kernel::kb, const t_uint Ju=4, const t_uint Jv=4, const bool w_stacking=false, const t_real &cellx=1, const t_real &celly=1)
Returns linear transform that is the standard degridding operator.
Definition: operators.h:608
purify::utilities::regroup_and_all_to_all
std::tuple< vis_params, std::vector< t_int > > regroup_and_all_to_all(vis_params const &params, const std::vector< t_int > &image_index, std::vector< t_int > const &groups, sopt::mpi::Communicator const &comm)
Definition: mpi_utilities.cc:94
purify::utilities::scatter_visibilities
vis_params scatter_visibilities(vis_params const &params, std::vector< t_int > const &sizes, sopt::mpi::Communicator const &comm)
Definition: mpi_utilities.cc:137
purify::utilities::regroup_and_scatter
vis_params regroup_and_scatter(vis_params const &params, std::vector< t_int > const &groups, sopt::mpi::Communicator const &comm)
Definition: mpi_utilities.cc:77
purify::utilities::random_sample_density
utilities::vis_params random_sample_density(const t_int vis_num, const t_real mean, const t_real standard_deviation, const t_real rms_w)
Generates a random visibility coverage.
Definition: uvw_utilities.cc:104
operators_gpu.h
purify::utilities::vis_params
Definition: uvw_utilities.h:15
wproj_operators.h
wproj_operators_gpu.h
Generated by doxygen 1.9.1