PURIFY
Next-generation radio interferometric imaging
Functions
operators_gpu.cc File Reference
#include "purify/operators_gpu.h"
#include "purify/config.h"
#include "purify/types.h"
#include <arrayfire.h>
#include <iostream>
#include "catch2/catch_all.hpp"
#include "purify/directories.h"
#include "purify/kernels.h"
#include "purify/logging.h"
#include "purify/operators.h"
#include "purify/wproj_operators_gpu.h"
#include <sopt/power_method.h>
+ Include dependency graph for operators_gpu.cc:

Go to the source code of this file.

Functions

 TEST_CASE ("GPU Operators")
 
 TEST_CASE ("wprojection")
 

Function Documentation

◆ TEST_CASE() [1/2]

TEST_CASE ( "GPU Operators"  )

Definition at line 14 of file operators_gpu.cc.

14  {
15  af::setDevice(0);
16  af::info();
17  // sopt::logging::set_level("debug");
18  // purify::logging::set_level("debug");
19  const t_uint M = 1000;
20  const t_real oversample_ratio = 2;
21  const t_uint imsizex = 128;
22  const t_uint imsizey = 128;
23  const t_uint N = imsizex * imsizey;
24  const t_uint ftsizev = std::floor(imsizey * oversample_ratio);
25  const t_uint ftsizeu = std::floor(imsizex * oversample_ratio);
26  const t_uint Ju = 4;
27  const t_uint Jv = 4;
28  const t_uint power_iters = 100;
29  const t_real power_tol = 1e-9;
30  const kernels::kernel kernel = kernels::kernel::kb;
31  const std::string& weighting_type = "natural";
32  const t_real R = 0;
33  auto u = Vector<t_real>::Random(M);
34  auto v = Vector<t_real>::Random(M);
35  utilities::vis_params uv_vis;
36  uv_vis.u = u;
37  uv_vis.v = v;
38  uv_vis.w = u * 0.;
39  uv_vis.weights = Vector<t_complex>::Random(M);
40  uv_vis.vis = Vector<t_complex>::Random(M);
41  uv_vis.units = utilities::vis_units::pixels;
42  std::function<t_real(t_real)> kbu, kbv, ftkbu, ftkbv;
43  std::tie(kbu, kbv, ftkbu, ftkbv) =
44  create_kernels(kernel, Ju, Jv, imsizey, imsizex, oversample_ratio);
45  SECTION("FFT") {
46  const Vector<t_complex> input = Vector<t_complex>::Random(ftsizev * ftsizeu);
47  Vector<t_complex> output;
48  Vector<t_complex> new_input;
49  sopt::OperatorFunction<af::array> direct_gpu_fft, indirect_gpu_fft;
50  std::tie(direct_gpu_fft, indirect_gpu_fft) =
51  gpu::operators::init_af_FFT_2d(imsizey, imsizex, oversample_ratio);
52  sopt::OperatorFunction<Vector<t_complex>> fft =
53  gpu::host_wrapper(direct_gpu_fft, input.size(), input.size());
54  sopt::OperatorFunction<Vector<t_complex>> ifft =
55  gpu::host_wrapper(indirect_gpu_fft, input.size(), input.size());
56  fft(output, input);
57  ifft(new_input, output);
58  CAPTURE(input.transpose().array() / new_input.transpose().array());
59  CHECK(input.isApprox(new_input, 1e-6));
60  sopt::OperatorFunction<Vector<t_complex>> direct_fft, indirect_fft;
61  std::tie(direct_fft, indirect_fft) =
62  operators::init_FFT_2d<Vector<t_complex>>(imsizey, imsizex, oversample_ratio);
63  const Vector<t_complex> direct_input = Vector<t_complex>::Random(ftsizev * ftsizeu);
64  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(ftsizev * ftsizeu);
65  Vector<t_complex> output_new;
66  Vector<t_complex> output_old;
67  fft(output_new, direct_input);
68  direct_fft(output_old, direct_input);
69  CHECK(output_old.isApprox(output_new, 1e-6));
70  ifft(output_new, indirect_input);
71  indirect_fft(output_old, indirect_input);
72  CHECK(output_old.isApprox(output_new, 1e-6));
73  }
74  SECTION("Gridding") {
75  sopt::OperatorFunction<af::array> direct_gpu_G, indirect_gpu_G;
76  std::tie(direct_gpu_G, indirect_gpu_G) = gpu::operators::init_af_gridding_matrix_2d(
77  uv_vis.u, uv_vis.v, Vector<t_complex>::Constant(M, 1.), imsizey, imsizex, oversample_ratio,
78  kbu, kbv, Ju, Jv);
79  const sopt::OperatorFunction<Vector<t_complex>> gpu_direct_G =
80  gpu::host_wrapper(direct_gpu_G, ftsizeu * ftsizev, M);
81  const sopt::OperatorFunction<Vector<t_complex>> gpu_indirect_G =
82  gpu::host_wrapper(indirect_gpu_G, M, ftsizeu * ftsizev);
83  sopt::OperatorFunction<Vector<t_complex>> direct_G, indirect_G;
84  std::tie(direct_G, indirect_G) = operators::init_gridding_matrix_2d<Vector<t_complex>>(
85  uv_vis.u, uv_vis.v, Vector<t_complex>::Constant(M, 1.), imsizey, imsizex, oversample_ratio,
86  kbv, kbu, Ju, Jv);
87  const Vector<t_complex> direct_input = Vector<t_complex>::Random(ftsizev * ftsizeu);
88  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(M);
89  Vector<t_complex> output_new;
90  Vector<t_complex> output_old;
91  gpu_direct_G(output_new, direct_input);
92  direct_G(output_old, direct_input);
93  CHECK(output_new.isApprox(output_old, 1e-6));
94  gpu_indirect_G(output_new, indirect_input);
95  indirect_G(output_old, indirect_input);
96  CHECK(output_new.isApprox(output_old, 1e-6));
97  }
98  SECTION("Zero Padding") {
99  const Image<t_complex> S =
100  details::init_correction2d(oversample_ratio, imsizey, imsizex, ftkbu, ftkbv, 0, 1, 1);
101  CHECK(imsizex == S.cols());
102  CHECK(imsizey == S.rows());
103  sopt::OperatorFunction<Vector<t_complex>> directZ, indirectZ;
104  std::tie(directZ, indirectZ) =
105  operators::init_zero_padding_2d<Vector<t_complex>>(S, oversample_ratio);
106  sopt::OperatorFunction<af::array> directZgpu, indirectZgpu;
107  std::tie(directZgpu, indirectZgpu) =
108  gpu::operators::init_af_zero_padding_2d(S.cast<t_complexf>(), oversample_ratio);
109  sopt::OperatorFunction<Vector<t_complex>> directZ_gpu =
110  gpu::host_wrapper(directZgpu, imsizex * imsizey, ftsizeu * ftsizev);
111  sopt::OperatorFunction<Vector<t_complex>> indirectZ_gpu =
112  gpu::host_wrapper(indirectZgpu, ftsizeu * ftsizev, imsizex * imsizey);
113  const Vector<t_complex> direct_input = Vector<t_complex>::Random(imsizex * imsizey);
114  Vector<t_complex> direct_output_old;
115  Vector<t_complex> direct_output_new;
116  directZ(direct_output_old, direct_input);
117  directZ_gpu(direct_output_new, direct_input);
118  CHECK(direct_output_new.size() == ftsizeu * ftsizev);
119  CHECK(direct_output_old.isApprox(direct_output_new, 1e-6));
120  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(ftsizeu * ftsizev);
121  Vector<t_complex> indirect_output_old;
122  Vector<t_complex> indirect_output_new;
123  indirectZ(indirect_output_old, indirect_input);
124  indirectZ_gpu(indirect_output_new, indirect_input);
125  CHECK(indirect_output_new.size() == imsizex * imsizey);
126  CHECK(indirect_output_old.isApprox(indirect_output_new, 1e-6));
127  }
128  SECTION("Serial Operator") {
129  const auto measure_op = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
130  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
131  uv_vis.u, uv_vis.v, uv_vis.w, uv_vis.weights, imsizey, imsizex, oversample_ratio,
132  kernel, Ju, Jv),
133  power_iters, power_tol, Vector<t_complex>::Random(imsizex * imsizey)));
134 
135  const auto measure_op_gpu = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
136  gpu::measurementoperator::init_degrid_operator_2d(uv_vis.u, uv_vis.v, uv_vis.w,
137  uv_vis.weights, imsizey, imsizex,
138  oversample_ratio, kernel, Ju, Jv),
139  power_iters, power_tol, Vector<t_complex>::Random(imsizex * imsizey)));
140  const Vector<t_complex> direct_input = Vector<t_complex>::Random(imsizex * imsizey);
141  const Vector<t_complex> direct_output = *measure_op_gpu * direct_input;
142  CHECK(direct_output.size() == M);
143  const Vector<t_complex> indirect_input = Vector<t_complex>::Random(M);
144  const Vector<t_complex> indirect_output = measure_op_gpu->adjoint() * indirect_input;
145  CHECK(indirect_output.size() == imsizex * imsizey);
146  SECTION("Power Method") {
147  auto op_norm = std::get<0>(sopt::algorithm::power_method<Vector<t_complex>>(
148  *measure_op, power_iters, power_tol, Vector<t_complex>::Random(imsizex * imsizey)));
149  CHECK(std::abs(op_norm - 1.) < power_tol);
150  }
151  SECTION("Degrid") {
152  const Vector<t_complex> input = Vector<t_complex>::Random(imsizex * imsizey);
153  const Vector<t_complex> expected_output = *measure_op * input;
154  const Vector<t_complex> actual_output = *measure_op_gpu * input;
155  CHECK(expected_output.size() == actual_output.size());
156  CHECK(actual_output.isApprox(expected_output, 1e-4));
157  }
158  SECTION("Grid") {
159  const Vector<t_complex> input = Vector<t_complex>::Random(M);
160  const Vector<t_complex> expected_output = measure_op->adjoint() * input;
161  const Vector<t_complex> actual_output = measure_op_gpu->adjoint() * input;
162  CHECK(expected_output.size() == actual_output.size());
163  CHECK(actual_output.isApprox(expected_output, 1e-4));
164  }
165  }
166 }
CHECK
#define CHECK(CONDITION, ERROR)
Definition: casa.cc:6
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::direct_input
const std::vector< t_complex > direct_input
data for degridding input
Definition: operators.cc:22
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::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::t_complexf
std::complex< float > t_complexf
Definition: types.h:21
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::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

References CHECK, purify::create_kernels(), operators_test::direct_input, operators_test::indirect_input, purify::details::init_correction2d(), purify::measurementoperator::init_degrid_operator_2d(), purify::kernels::kb, purify::utilities::pixels, purify::utilities::vis_params::u, operators_test::u, purify::utilities::vis_params::units, purify::utilities::vis_params::v, operators_test::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, and purify::utilities::vis_params::weights.

◆ TEST_CASE() [2/2]

TEST_CASE ( "wprojection"  )

Definition at line 167 of file operators_gpu.cc.

167  {
168  Image<t_complex> const M31 = Image<t_complex>::Random(256, 256);
169  Vector<t_complex> const input = Vector<t_complex>::Map(M31.data(), M31.size());
170  const Vector<t_complex> vis = Vector<t_complex>::Random(10);
171  const t_uint imsizex = M31.cols();
172  const t_uint imsizey = M31.rows();
173  const t_uint Jw = 30;
174  const t_real cell_x = 1;
175  const t_real cell_y = 1;
176  const t_uint power_iters = 1000;
177  const t_real power_tol = 1e-4;
178  const t_real abs_error = 1e-9;
179  const t_real rel_error = 1e-9;
180  const bool w_term = false;
181  const t_uint M = 10;
182  utilities::vis_params uv_data;
183  uv_data.u = Vector<t_real>::Random(M);
184  uv_data.v = Vector<t_real>::Random(M);
185  uv_data.w = Vector<t_real>::Zero(M);
186  uv_data.weights = Vector<t_complex>::Ones(M);
187  uv_data.vis = Vector<t_complex>::Ones(M);
188  SECTION("oversample 2") {
189  const kernels::kernel kernel = kernels::kernel::kb;
190  const t_real oversample_ratio = 2;
191  const t_uint Ju = 4;
192  const t_uint Jv = 4;
193  auto mop = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
194  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
195  uv_data, imsizey, imsizex, cell_x, cell_y, oversample_ratio, kernel, Ju, Jv, w_term),
196  power_iters, power_tol, Vector<t_complex>::Random(imsizex * imsizey)));
197  auto mop_wproj = std::get<2>(sopt::algorithm::normalise_operator<Vector<t_complex>>(
198  measurementoperator::init_degrid_operator_2d<Vector<t_complex>>(
199  uv_data, imsizey, imsizex, cell_x, cell_y, oversample_ratio, kernel, Ju, Jw, w_term,
200  abs_error, rel_error, dde_type::wkernel_radial),
201  power_iters, power_tol, Vector<t_complex>::Random(imsizex * imsizey)));
202  REQUIRE((mop_wproj->adjoint() * vis).size() == imsizex * imsizey);
203  REQUIRE((mop->adjoint() * vis).size() == imsizex * imsizey);
204  REQUIRE((*mop * input).size() == M);
205  REQUIRE((*mop_wproj * input).size() == M);
206  const t_real op_norm = std::get<0>(sopt::algorithm::power_method<Vector<t_complex>>(
207  {[=](Vector<t_complex>& output, const Vector<t_complex>& inp) {
208  output = (*mop * inp).eval() - (*mop_wproj * inp).eval();
209  },
210  [=](Vector<t_complex>& output, const Vector<t_complex>& inp) {
211  output = (mop->adjoint() * inp).eval() - (mop_wproj->adjoint() * inp).eval();
212  }},
213  power_iters, power_tol, input));
214  REQUIRE(op_norm == Approx(0).margin(1e-3));
215  }
216 }
purify::read_measurements::format::vis
@ vis

References purify::measurementoperator::init_degrid_operator_2d(), purify::kernels::kb, purify::utilities::vis_params::u, purify::utilities::vis_params::v, purify::utilities::vis_params::vis, purify::utilities::vis_params::w, purify::utilities::vis_params::weights, and purify::wkernel_radial.

Generated by doxygen 1.9.1