PURIFY
Next-generation radio interferometric imaging
utils.cc
Go to the documentation of this file.
1 #include <random>
2 #include "catch2/catch_all.hpp"
3 #include "purify/directories.h"
4 #include "purify/utilities.h"
5 #include "purify/uvw_utilities.h"
6 
7 using namespace purify;
8 using Catch::Approx;
9 
10 TEST_CASE("utilities [mod]", "[mod]") {
11  Array<t_real> range;
12  range.setLinSpaced(201, -100, 100);
13  Array<t_real> output(range.size());
14  for (t_int i = 0; i < range.size(); ++i) {
15  output(i) = utilities::mod(range(i), 20);
16  }
17  Array<t_real> expected(201);
18  expected << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4,
19  5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
20  12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
21  18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4,
22  5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
23  12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
24  18, 19, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0, 1, 2, 3, 4,
25  5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 0;
26  CHECK(expected.isApprox(output, 1e-13));
27 }
28 
29 TEST_CASE("utilities [reshape]", "[reshape]") {
30  // Testing if resize works the same as matlab's reshape
31  Matrix<t_real> magic(4, 4);
32  magic << 16, 2, 3, 13, 5, 11, 10, 8, 9, 7, 6, 12, 4, 14, 15, 1;
33  CAPTURE(magic);
34  Vector<t_real> magic_vector(16);
35  magic_vector << 16, 5, 9, 4, 2, 11, 7, 14, 3, 10, 6, 15, 13, 8, 12, 1;
36  Matrix<t_real> magic_matrix(16, 1);
37  magic_matrix << 16, 5, 9, 4, 2, 11, 7, 14, 3, 10, 6, 15, 13, 8, 12, 1;
38  magic_matrix.resize(4, 4);
39  CHECK(magic.isApprox(magic_matrix, 1e-13));
40  magic.resize(16, 1);
41  CAPTURE(magic);
42  CAPTURE(magic_vector);
43  CHECK(magic.isApprox(magic_vector, 1e-13));
44 }
45 TEST_CASE("utilities [variance]", "[variance]") {
46  // tests if mean and variance calculations are the same as done in matlab
47  Vector<t_complex> real_data = Vector<t_complex>::Random(1000);
48  Vector<t_complex> imag_data = Vector<t_complex>::Random(1000);
49  t_complex I(0, 1);
50  const Vector<t_complex> data = (real_data + I * imag_data);
51  const t_complex mu = utilities::mean(data);
52  CAPTURE(mu);
53  CHECK(std::abs(mu - data.mean()) < 1e-13);
54  const t_real var = utilities::variance(data);
55  const t_real expected_var =
56  std::abs(((data.array() - mu) * (data.array() - mu).conjugate()).sum()) / (data.size() - 1.);
57  CAPTURE(var);
58  CAPTURE(expected_var);
59  CHECK(std::abs(var - expected_var) < 1e-13);
60 }
61 TEST_CASE("utilities [median]", "[median]") {
62  Vector<t_real> x(100);
63  x << 0.824149356327936, 0.218232263733975, 0.0996423029616137, 0.619505816445823,
64  0.103814782651106, 0.799061803589637, 0.902925119282939, 0.312512845986667, 0.281589166159224,
65  0.00678194126805909, 0.495871665666786, 0.988482243012287, 0.737940751495519,
66  0.310719817391570, 0.600407500351151, 0.781679826661116, 0.111533124684943, 0.579329289388906,
67  0.870371220311720, 0.689776249611124, 0.242970250921629, 0.342722251869759, 0.545439787975426,
68  0.0675726926893163, 0.410448399698545, 0.237511397229394, 0.488973556424727,
69  0.806066588301788, 0.377847951671179, 0.517978144507388, 0.0945978893449516,
70  0.909095192849836, 0.207630593997342, 0.382064523540530, 0.660280265178685, 0.758373302277572,
71  0.173069412392096, 0.517379761506326, 0.995338112303466, 0.707609406590709,
72  0.0805673009239308, 0.0433080963751025, 0.491155562984260, 0.446596440169203,
73  0.486798565845346, 0.165891122085342, 0.360656504813112, 0.880721941763713, 0.744352179996200,
74  0.416771050919451, 0.907354808475253, 0.0943068737715302, 0.181326261018145,
75  0.946587232183202, 0.100849056860381, 0.388038069974207, 0.289224862024986,
76  0.0730868607975930, 0.194608005437083, 0.417485231103283, 0.292926559143940,
77  0.702138511791249, 0.239713177935322, 0.959482620708005, 0.305462202307119, 0.154915458213084,
78  0.555508419323977, 0.790543689112232, 0.443872293587844, 0.995818721727493, 0.436586698162556,
79  0.304441172973089, 0.246510340916772, 0.960825251654563, 0.222880958438666, 0.395609181466230,
80  0.224524553336017, 0.270024669015495, 0.418441479002589, 0.997735979773324, 0.911030403930781,
81  0.550426888716037, 0.596335815776380, 0.0791452393029073, 0.576636037390549,
82  0.898172311971638, 0.463310567259345, 0.398396426873099, 0.104461114798391, 0.652236388678710,
83  0.991700005353602, 0.678072815973294, 0.428477593633851, 0.654801321575878, 0.588745528605038,
84  0.745054910247788, 0.640870288169766, 0.503676186207024, 0.938041475805194, 0.605343580621484;
85  t_real const median = 0.475054566552345;
86  CHECK(std::abs(median - utilities::median(x)) < 1e-13);
87 }
88 TEST_CASE("utilities [read_write_vis]", "[read_write_vis]") {
89  // tests the read and write function for a visibility data set
90  std::string vis_file = vla_filename("at166B.3C129.c0.vis");
91  std::string out_file = output_filename("test_output.vis");
92  std::string out_w_file = output_filename("test_w_output.vis");
93  auto uv_data = utilities::read_visibility(vis_file);
94  utilities::write_visibility(uv_data, out_file);
95  auto new_uv_data = utilities::read_visibility(out_file);
96  CHECK(new_uv_data.u.isApprox(uv_data.u, 1e-13));
97  CHECK(new_uv_data.v.isApprox(uv_data.v, 1e-13));
98  CHECK(new_uv_data.vis.isApprox(uv_data.vis, 1e-13));
99  CHECK(new_uv_data.weights.isApprox(uv_data.weights, 1e-13));
100  t_int number_of_random_vis = 1e5;
101  const bool w_term = true;
102  auto random_uv_data = utilities::random_sample_density(number_of_random_vis, 0, constant::pi / 3);
103  utilities::write_visibility(random_uv_data, out_w_file, w_term);
104  auto new_random_uv_data = utilities::read_visibility(out_w_file, w_term);
105  CHECK(new_random_uv_data.u.isApprox(random_uv_data.u, 1e-8));
106  CHECK(new_random_uv_data.v.isApprox(random_uv_data.v, 1e-8));
107  CHECK(new_random_uv_data.w.isApprox(random_uv_data.w, 1e-8));
108  CHECK(new_random_uv_data.vis.isApprox(random_uv_data.vis, 1e-8));
109  CHECK(new_random_uv_data.weights.isApprox(random_uv_data.weights, 1e-8));
110 }
111 TEST_CASE("read_mutiple_vis") {
112  std::string csv_file = vla_filename("at166B.3C129.c0.vis");
113  std::string h5_file = atca_filename("0332-391.h5");
114  SECTION("one csv file") {
115  const std::vector<std::string> names = {csv_file};
116  const auto uv_data = utilities::read_visibility(csv_file);
117  const auto uv_multi = utilities::read_visibility(names);
118  CAPTURE(names.size());
119  CAPTURE(uv_data.size());
120  CHECK(uv_data.size() * names.size() == uv_multi.size());
121  for (int i = 0; i < names.size(); i++) {
122  CHECK(uv_data.u.isApprox(uv_multi.u.segment(i * uv_data.size(), uv_data.size())));
123  CHECK(uv_data.v.isApprox(uv_multi.v.segment(i * uv_data.size(), uv_data.size())));
124  CHECK(uv_data.w.isApprox(uv_multi.w.segment(i * uv_data.size(), uv_data.size())));
125  CHECK(uv_data.vis.isApprox(uv_multi.vis.segment(i * uv_data.size(), uv_data.size())));
126  CHECK(uv_data.weights.isApprox(uv_multi.weights.segment(i * uv_data.size(), uv_data.size())));
127  }
128  }
129  SECTION("many csv files") {
130  const std::vector<std::string> names = {csv_file, csv_file, csv_file};
131  const auto uv_data = utilities::read_visibility(csv_file);
132  const auto uv_multi = utilities::read_visibility(names);
133  CAPTURE(names.size());
134  CAPTURE(uv_data.size());
135  CHECK(uv_data.size() * names.size() == uv_multi.size());
136  for (int i = 0; i < names.size(); i++) {
137  CHECK(uv_data.u.isApprox(uv_multi.u.segment(i * uv_data.size(), uv_data.size())));
138  CHECK(uv_data.v.isApprox(uv_multi.v.segment(i * uv_data.size(), uv_data.size())));
139  CHECK(uv_data.w.isApprox(uv_multi.w.segment(i * uv_data.size(), uv_data.size())));
140  CHECK(uv_data.vis.isApprox(uv_multi.vis.segment(i * uv_data.size(), uv_data.size())));
141  CHECK(uv_data.weights.isApprox(uv_multi.weights.segment(i * uv_data.size(), uv_data.size())));
142  }
143  }
144 #ifdef PURIFY_H5
145  SECTION("one HDF5 file") {
146  const std::vector<std::string> names = {h5_file};
147  const auto uv_data = utilities::read_visibility(h5_file);
148  const auto uv_multi = utilities::read_visibility(names);
149  CAPTURE(names.size());
150  CAPTURE(uv_data.size());
151  CHECK(uv_data.size() * names.size() == uv_multi.size());
152  for (int i = 0; i < names.size(); i++) {
153  CHECK(uv_data.u.isApprox(uv_multi.u.segment(i * uv_data.size(), uv_data.size())));
154  CHECK(uv_data.v.isApprox(uv_multi.v.segment(i * uv_data.size(), uv_data.size())));
155  CHECK(uv_data.w.isApprox(uv_multi.w.segment(i * uv_data.size(), uv_data.size())));
156  CHECK(uv_data.vis.isApprox(uv_multi.vis.segment(i * uv_data.size(), uv_data.size())));
157  CHECK(uv_data.weights.isApprox(uv_multi.weights.segment(i * uv_data.size(), uv_data.size())));
158  }
159  }
160 #endif
161 }
162 TEST_CASE("utilities [file exists]", "[file exists]") {
163  std::string vis_file = vla_filename("at166B.3C129.c0.vis");
164  // File should exist
165  CHECK(utilities::file_exists(vis_file));
166  // File should not exist
167  CHECK(not utilities::file_exists("adfadsf"));
168 }
169 TEST_CASE("utilities [fit_fwhm]", "[fit_fwhm]") {
170  // testing that the gaussian fitting works.
171  t_int imsizex = 512;
172  t_int imsizey = 512;
173  Image<t_real> psf = Image<t_real>::Zero(imsizey, imsizex);
174  // choice of parameters
175  t_real const fwhm_x = 3;
176  t_real const fwhm_y = 6;
177  t_real const theta = 0;
178  // setting up Gaussian calculation
179  t_real const sigma_x = fwhm_x / (2 * std::sqrt(2 * std::log(2)));
180  t_real const sigma_y = fwhm_y / (2 * std::sqrt(2 * std::log(2)));
181  // calculating Gaussian
182  t_real const a = std::pow(std::cos(theta), 2) / (2 * sigma_x * sigma_x) +
183  std::pow(std::sin(theta), 2) / (2 * sigma_y * sigma_y);
184  t_real const b = -std::sin(2 * theta) / (4 * sigma_x * sigma_x) +
185  std::sin(2 * theta) / (4 * sigma_y * sigma_y);
186  t_real const c = std::pow(std::sin(theta), 2) / (2 * sigma_x * sigma_x) +
187  std::pow(std::cos(theta), 2) / (2 * sigma_y * sigma_y);
188  auto x0 = imsizex * 0.5;
189  auto y0 = imsizey * 0.5;
190  for (t_int i = 0; i < imsizex; ++i) {
191  for (t_int j = 0; j < imsizey; ++j) {
192  t_real x = i - x0;
193  t_real y = j - y0;
194  psf(j, i) = std::exp(-a * x * x + 2 * b * x * y - c * y * y);
195  }
196  }
197  auto const fit = utilities::fit_fwhm(psf, 3);
198  auto new_fwhm_x = std::get<0>(fit) * 2 * std::sqrt(2 * std::log(2));
199  auto new_fwhm_y = std::get<1>(fit) * 2 * std::sqrt(2 * std::log(2));
200  auto new_theta = std::get<2>(fit);
201  CHECK(std::abs(new_fwhm_x - fwhm_x) < 1e-13);
202  CHECK(std::abs(new_fwhm_y - fwhm_y) < 1e-13);
203  CHECK(std::abs(new_theta - theta) < 1e-13);
204 }
205 
206 TEST_CASE("utilities [sparse multiply]", "[sparse multiply]") {
207  // Checking that the parallel sparse matrix multiplication works against Eigen.
208  t_int cols = 1024 * 1024;
209  t_int rows = 1e3;
210  t_int nz_values = 16 * rows;
211  Vector<t_complex> const x = Vector<t_complex>::Random(cols);
212 
213  std::vector<t_tripletList> tripletList;
214  tripletList.reserve(nz_values);
215  Vector<t_complex> M_values = Vector<t_complex>::Random(nz_values);
216  std::random_device rd;
217  std::mt19937 gen(rd());
218  std::uniform_real_distribution<> dis_row(0, rows);
219  std::uniform_real_distribution<> dis_col(0, cols);
220  for (t_int i = 0; i < nz_values; ++i) {
221  tripletList.emplace_back(std::floor(dis_row(rd)), std::floor(dis_col(rd)), M_values(i));
222  }
223  Sparse<t_complex> M(rows, cols);
224  M.setFromTriplets(tripletList.begin(), tripletList.end());
225 
226  Vector<t_complex> const parallel_output = utilities::sparse_multiply_matrix(M, x);
227  Vector<t_complex> const correct_output = M * x;
228 
229  for (t_int i = 0; i < rows; ++i) {
230  CHECK(std::abs((correct_output(i) - parallel_output(i)) / correct_output(i)) < 1e-13);
231  }
232 }
233 
234 TEST_CASE("generate_baseline") {
235  // testing if randomly generating a uvcoverage from baseline configuration works
236  const Matrix<t_real> B = utilities::generate_antennas(4, 1);
237  CHECK(B.allFinite());
238  CHECK(B.rows() == 4);
239  CHECK(B.cols() == 3);
240  CAPTURE(B);
241  SECTION("one wavelength") {
242  const t_real frequency = constant::c;
243  const t_real times = 0.;
244  const t_real phi_dec = 0.;
245  const t_real theta_ra = 0.;
246  const t_real latitude = 0.;
247  const utilities::vis_params test_coverage =
248  utilities::antenna_to_coverage(B, frequency, times, theta_ra, phi_dec, latitude);
249  CHECK(test_coverage.units == utilities::vis_units::lambda);
250  const utilities::vis_params test_coverage_xyz = utilities::antenna_to_coverage(
251  B.col(0), B.col(1), B.col(2), frequency, times, theta_ra, phi_dec, latitude);
252  CHECK(test_coverage_xyz.u.isApprox(test_coverage.u));
253  CHECK(test_coverage_xyz.v.isApprox(test_coverage.v));
254  CHECK(test_coverage_xyz.w.isApprox(test_coverage.w));
255  CHECK(test_coverage.u.allFinite());
256  CHECK(test_coverage.v.allFinite());
257  CHECK(test_coverage.w.allFinite());
258  CHECK(test_coverage.size() == 2 * 3);
259  const Vector<t_real> R0 = B.row(0) - B.row(1);
260  const Vector<t_real> R1 = B.row(0) - B.row(2);
261  const Vector<t_real> R2 = B.row(0) - B.row(3);
262  const Vector<t_real> R3 = B.row(1) - B.row(2);
263  const Vector<t_real> R4 = B.row(1) - B.row(3);
264  const Vector<t_real> R5 = B.row(2) - B.row(3);
265  CHECK(R0(0) == Approx(test_coverage.u(0)));
266  CHECK(R0(1) == Approx(test_coverage.v(0)));
267  CHECK(R0(2) == Approx(test_coverage.w(0)));
268  CHECK(R1(0) == Approx(test_coverage.u(1)));
269  CHECK(R1(1) == Approx(test_coverage.v(1)));
270  CHECK(R1(2) == Approx(test_coverage.w(1)));
271  CHECK(R2(0) == Approx(test_coverage.u(2)));
272  CHECK(R2(1) == Approx(test_coverage.v(2)));
273  CHECK(R2(2) == Approx(test_coverage.w(2)));
274  CHECK(R3(0) == Approx(test_coverage.u(3)));
275  CHECK(R3(1) == Approx(test_coverage.v(3)));
276  CHECK(R3(2) == Approx(test_coverage.w(3)));
277  CHECK(R4(0) == Approx(test_coverage.u(4)));
278  CHECK(R4(1) == Approx(test_coverage.v(4)));
279  CHECK(R4(2) == Approx(test_coverage.w(4)));
280  CHECK(R5(0) == Approx(test_coverage.u(5)));
281  CHECK(R5(1) == Approx(test_coverage.v(5)));
282  CHECK(R5(2) == Approx(test_coverage.w(5)));
283  }
284  SECTION("more wavelengths") {
285  const auto frequency = std::vector<t_real>{constant::c, 2 * constant::c, 5 * constant::c};
286  const auto times = std::vector<t_real>{0, 2, 100};
287  const t_real phi_dec = 0.;
288  const t_real theta_ra = 0.;
289  const t_real latitude = 0.;
290  const t_int M = B.rows() * (B.rows() - 1) / 2;
291  const utilities::vis_params test_coverage =
292  utilities::antenna_to_coverage(B, frequency, times, theta_ra, phi_dec, latitude);
293  CHECK(test_coverage.u.allFinite());
294  CHECK(test_coverage.v.allFinite());
295  CHECK(test_coverage.w.allFinite());
296  CHECK(test_coverage.size() == M * frequency.size() * times.size());
297  t_int f_index = 0;
298  t_int t_index = 0;
299  for (t_int k = 0; k < frequency.size(); k++) {
300  for (t_int j = 0; j < times.size(); j++) {
301  const t_real f = frequency.at(k);
302  const t_real t = times.at(j);
303  CAPTURE(f);
304  CAPTURE(t);
305  CAPTURE(f_index);
306  CAPTURE(frequency.size());
307  const utilities::vis_params test_coverage_f =
308  utilities::antenna_to_coverage(B, f, t, theta_ra, phi_dec, latitude);
309  const t_int index = (k + j * frequency.size()) * M;
310  CAPTURE(test_coverage.u.segment(index, M).head(5));
311  CAPTURE(test_coverage_f.u.head(5));
312  CHECK(test_coverage_f.u.isApprox(test_coverage.u.segment(index, M)));
313  CHECK(test_coverage_f.v.isApprox(test_coverage.v.segment(index, M)));
314  CHECK(test_coverage_f.w.isApprox(test_coverage.w.segment(index, M)));
315  }
316  }
317  }
318 }
319 
320 TEST_CASE("generate coverage from antenna positions") {
321  const std::string pos_filename = mwa_filename("Phase2_config.txt");
322 
323  auto const B = utilities::read_ant_positions(pos_filename);
324  CHECK(B.rows() == 128);
325  CHECK(B.cols() == 3);
326  CHECK(B.allFinite());
327  SECTION("generate coverage") {
328  SECTION("one frequency") {
329  const t_real f = constant::c;
330  const t_real t = 0.;
331  const t_real latitude = 0.;
332  auto const coverage_from_file =
333  utilities::read_ant_positions_to_coverage(pos_filename, f, t, 0., 0., latitude);
334  auto const coverage_from_B = utilities::antenna_to_coverage(B, f, t, 0., 0., latitude);
335  CHECK(coverage_from_file.u.isApprox(coverage_from_B.u));
336  CHECK(coverage_from_file.v.isApprox(coverage_from_B.v));
337  CHECK(coverage_from_file.w.isApprox(coverage_from_B.w));
338  }
339  SECTION("multi frequency") {
340  const std::vector<t_real> f = {constant::c, constant::c * 4, constant::c * 2};
341  const std::vector<t_real> t = {0., 4, 24};
342  const t_real latitude = 0.;
343  auto const coverage_from_file =
344  utilities::read_ant_positions_to_coverage(pos_filename, f, t, 0., 0., latitude);
345  auto const coverage_from_B = utilities::antenna_to_coverage(B, f, t, 0., 0., latitude);
346  CHECK(coverage_from_file.u.isApprox(coverage_from_B.u));
347  CHECK(coverage_from_file.v.isApprox(coverage_from_B.v));
348  CHECK(coverage_from_file.w.isApprox(coverage_from_B.w));
349  }
350  }
351 }
352 
353 TEST_CASE("rotations") {
354  SECTION("u = 0, v = 0, w = 1") {
355  const t_real theta_0 = 0.;
356  const t_real phi_0 = 0.;
357  CHECK(0. ==
358  Approx(utilities::calculate_rotated_u(0., 0., 1., theta_0, phi_0, 0.)).margin(1e-6));
359  CHECK(0. ==
360  Approx(utilities::calculate_rotated_v(0., 0., 1., theta_0, phi_0, 0.)).margin(1e-6));
361  CHECK(1. == Approx(utilities::calculate_rotated_w(0., 0., 1., theta_0, phi_0, 0.)));
362  }
363  SECTION("u = 1, v = 0, w = 0") {
364  const t_real theta_0 = 0;
365  const t_real phi_0 = constant::pi / 2.;
366  CHECK(1. == Approx(utilities::calculate_rotated_u(0., 0., 1., theta_0, phi_0, 0.)));
367  CHECK(0. ==
368  Approx(utilities::calculate_rotated_v(0., 0., 1., theta_0, phi_0, 0.)).margin(1e-6));
369  CHECK(0. ==
370  Approx(utilities::calculate_rotated_w(0., 0., 1., theta_0, phi_0, 0.)).margin(1e-6));
371  }
372  SECTION("u = 0, v = 1, w = 0") {
373  const t_real theta_0 = constant::pi / 2.;
374  const t_real phi_0 = constant::pi / 2.;
375  CHECK(0. ==
376  Approx(utilities::calculate_rotated_u(0., 0., 1., theta_0, phi_0, 0.)).margin(1e-6));
377  CHECK(1. == Approx(utilities::calculate_rotated_v(0., 0., 1., theta_0, phi_0, 0.)));
378  CHECK(0. ==
379  Approx(utilities::calculate_rotated_w(0., 0., 1., theta_0, phi_0, 0.)).margin(1e-6));
380  }
381  SECTION("inverse") {
382  const t_real offset_alpha = constant::pi / 2.;
383  const t_real offset_beta = constant::pi / 3.;
384  const t_real offset_gamma = constant::pi / 3.;
385  const Vector<t_real> u = Vector<t_real>::Random(10);
386  const Vector<t_real> v = Vector<t_real>::Random(10);
387  const Vector<t_real> w = Vector<t_real>::Random(10);
388 
389  const Vector<t_real> rotated_u =
390  utilities::calculate_rotated_u(u, v, w, offset_alpha, offset_beta, offset_gamma);
391  const Vector<t_real> rotated_v =
392  utilities::calculate_rotated_v(u, v, w, offset_alpha, offset_beta, offset_gamma);
393  const Vector<t_real> rotated_w =
394  utilities::calculate_rotated_w(u, v, w, offset_alpha, offset_beta, offset_gamma);
395 
396  const Vector<t_real> inv_u = utilities::calculate_rotated_u(
397  rotated_u, rotated_v, rotated_w, (-offset_gamma), (-offset_beta), -offset_alpha);
398  const Vector<t_real> inv_v = utilities::calculate_rotated_v(
399  rotated_u, rotated_v, rotated_w, (-offset_gamma), (-offset_beta), -offset_alpha);
400  const Vector<t_real> inv_w = utilities::calculate_rotated_w(
401  rotated_u, rotated_v, rotated_w, (-offset_gamma), (-offset_beta), -offset_alpha);
402  CAPTURE(u);
403  CAPTURE(inv_u);
404  CHECK(u.isApprox(inv_u));
405  CHECK(v.isApprox(inv_v));
406  CHECK(w.isApprox(inv_w));
407  }
408 }
409 
410 TEST_CASE("conjugate symmetry") {
411  t_uint const number_of_vis = 100;
412  t_uint const max_w = 100;
413  t_real const sigma_m = 1000;
414  const auto uv_data = utilities::random_sample_density(number_of_vis, 0, sigma_m, max_w);
415  const auto reflected_data = utilities::conjugate_w(uv_data);
416  REQUIRE(uv_data.size() == reflected_data.size());
417  for (t_uint i = 0; i < uv_data.size(); i++) {
418  if (uv_data.w(i) < 0) {
419  REQUIRE(uv_data.u(i) == Approx(-reflected_data.u(i)).epsilon(1e-12));
420  REQUIRE(uv_data.v(i) == Approx(-reflected_data.v(i)).epsilon(1e-12));
421  REQUIRE(uv_data.w(i) == Approx(-reflected_data.w(i)).epsilon(1e-12));
422  REQUIRE(uv_data.vis(i).real() ==
423  Approx(std::conj(reflected_data.vis(i)).real()).epsilon(1e-12));
424  REQUIRE(uv_data.vis(i).imag() ==
425  Approx(std::conj(reflected_data.vis(i)).imag()).epsilon(1e-12));
426  } else {
427  REQUIRE(uv_data.u(i) == Approx(reflected_data.u(i)).epsilon(1e-12));
428  REQUIRE(uv_data.v(i) == Approx(reflected_data.v(i)).epsilon(1e-12));
429  REQUIRE(uv_data.w(i) == Approx(reflected_data.w(i)).epsilon(1e-12));
430  REQUIRE(uv_data.vis(i).real() == Approx(reflected_data.vis(i).real()).epsilon(1e-12));
431  REQUIRE(uv_data.vis(i).imag() == Approx(reflected_data.vis(i).imag()).epsilon(1e-12));
432  }
433  }
434 }
CHECK
#define CHECK(CONDITION, ERROR)
Definition: casa.cc:6
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
purify::constant::c
const t_real c
speed of light in vacuum
Definition: types.h:72
purify::constant::pi
const t_real pi
mathematical constant
Definition: types.h:70
purify::utilities::write_visibility
void write_visibility(const utilities::vis_params &uv_vis, const std::string &file_name, const bool w_term)
Writes visibilities to txt.
Definition: uvw_utilities.cc:243
purify::utilities::antenna_to_coverage
utilities::vis_params antenna_to_coverage(const Matrix< t_real > &B, const t_real frequency, const t_real times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
Provided antenna positions generate a coverage for a fixed frequency a fixed time.
Definition: uvw_utilities.cc:43
purify::utilities::sparse_multiply_matrix
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
purify::utilities::generate_antennas
Matrix< t_real > generate_antennas(const t_uint N, const t_real scale)
Generate guassianly distributed (sigma = scale) antenna positions.
Definition: uvw_utilities.cc:19
purify::utilities::median
t_real median(const Vector< t_real > &input)
Return median of real vector.
Definition: utilities.cc:221
purify::utilities::calculate_rotated_v
T calculate_rotated_v(const T &u, const T &v, const T &w, const t_real alpha, const t_real beta, const t_real gamma)
calculate the rotated v from euler angles in zyz and starting coordinates (u, v, w)
Definition: uvw_utilities.h:69
purify::utilities::fit_fwhm
std::tuple< t_real, t_real, t_real > fit_fwhm(const Image< t_real > &psf, const t_int &size)
Method to fit Gaussian to PSF.
Definition: utilities.cc:148
purify::utilities::variance
t_real variance(const K x)
Calculate variance of vector.
Definition: utilities.h:27
purify::utilities::mean
K::Scalar mean(const K x)
Calculate mean of vector.
Definition: utilities.h:21
purify::utilities::read_ant_positions
Matrix< t_real > read_ant_positions(const std::string &pos_name)
Read in a text file of antenna positions into a matrix [x, y ,z].
Definition: uvw_utilities.cc:77
purify::utilities::read_ant_positions_to_coverage
utilities::vis_params read_ant_positions_to_coverage(const std::string &pos_name, const T &frequencies, const K &times, const t_real theta_ra, const t_real phi_dec, const t_real latitude)
Definition: uvw_utilities.h:183
purify::utilities::calculate_rotated_u
T calculate_rotated_u(const T &u, const T &v, const T &w, const t_real alpha, const t_real beta, const t_real gamma)
calculate the rotated u from euler angles in zyz and starting coordinates (u, v, w)
Definition: uvw_utilities.h:59
purify::utilities::conjugate_w
utilities::vis_params conjugate_w(const utilities::vis_params &uv_vis)
reflects visibilities into the w >= 0 domain
Definition: uvw_utilities.cc:352
purify::utilities::file_exists
bool file_exists(const std::string &name)
Test to see if file exists.
Definition: utilities.cc:137
purify::utilities::read_visibility
utilities::vis_params read_visibility(const std::vector< std::string > &names, const bool w_term)
Read visibility files from name of vector.
Definition: uvw_utilities.cc:144
purify::utilities::calculate_rotated_w
T calculate_rotated_w(const T &u, const T &v, const T &w, const t_real alpha, const t_real beta, const t_real gamma)
calculate the rotated w from euler angles in zyz and starting coordinates (u, v, w)
Definition: uvw_utilities.h:79
purify::utilities::mod
t_real mod(const t_real &x, const t_real &y)
Mod function modified to wrap circularly for negative numbers.
Definition: utilities.cc:41
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
purify::output_filename
std::string output_filename(std::string const &filename)
Test output file.
Definition: directories.in.h:49
purify::atca_filename
std::string atca_filename(std::string const &filename)
Specific atca data.
Definition: directories.in.h:32
purify::mwa_filename
std::string mwa_filename(std::string const &filename)
Specific mwa data.
Definition: directories.in.h:36
purify::Sparse
Eigen::SparseMatrix< T, Eigen::RowMajor, I > Sparse
A matrix of a given type.
Definition: types.h:40
purify::vla_filename
std::string vla_filename(std::string const &filename)
Specific vla data.
Definition: directories.in.h:28
purify::utilities::vis_params
Definition: uvw_utilities.h:15
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::size
t_uint size() const
return number of measurements
Definition: uvw_utilities.h:54
purify::utilities::vis_params::units
vis_units units
Definition: uvw_utilities.h:24
TEST_CASE
TEST_CASE("utilities [mod]", "[mod]")
Definition: utils.cc:10
uvw_utilities.h
Generated by doxygen 1.9.1