AST-1527 Fix AVX-enabled EveryBeam

.gitlab-ci.common.yml

@@ -114,16 +114,20 @@ format-2204:
     - git submodule update --init external/aocommon
     - ./scripts/run-format.sh
 
+.download-wsrt-measures:
+  before_script:
+    - wget -q ftp://ftp.astron.nl/outgoing/Measures/WSRT_Measures.ztar
+    - tar -xf WSRT_Measures.ztar -C /var/lib/casacore/data/
+    - rm -f WSRT_Measures.ztar
+
 # Build a debug version of EveryBeam from the base image
 test-and-coverage-2204:
-  extends: .needs-base-2204
+  extends: [ .needs-base-2204, .download-wsrt-measures ]
   stage: build
   variables:
     GIT_SUBMODULE_STRATEGY: normal
   script:
     - WORKDIR=$PWD
-    # Download casacore wsrt measures
-    - wget -q ftp://ftp.astron.nl/outgoing/Measures/WSRT_Measures.ztar && tar -xf WSRT_Measures.ztar -C /var/lib/casacore/data/ && rm -f WSRT_Measures.ztar
     # Build in Debug mode
     - mkdir -p build/coeffs/lobes && cd build
     - mv /coeffs/lobes/* ./coeffs/lobes/
@@ -161,14 +165,12 @@ test-and-coverage-2204:
 # Python tests do not work with the sanitizer, since the main executable
 # (python) was not built with the sanitizer.
 sanitize-2204:
-  extends: [ ".needs-base-2204", ".dind-requester" ]
+  extends: [ .needs-base-2204, .dind-requester, .download-wsrt-measures ]
   stage: build
   image: $BASE_IMAGE_2204
   variables:
     GIT_SUBMODULE_STRATEGY: normal
   script:
-    # Download casacore wsrt measures
-    - wget -q ftp://ftp.astron.nl/outgoing/Measures/WSRT_Measures.ztar && tar -xf WSRT_Measures.ztar -C /var/lib/casacore/data/ && rm -f WSRT_Measures.ztar
     # Build and run unit tests.
     - mkdir -p build/coeffs/lobes && cd build
     - mv /coeffs/lobes/* ./coeffs/lobes/
@@ -183,14 +185,15 @@ sanitize-2204:
   script:
     - mkdir -p build/coeffs/lobes && cd build
     - mv /coeffs/lobes/* ./coeffs/lobes/
-    - echo Configuring with CXX flags \"$CXX_FLAGS\"
     - cmake -DBUILD_TESTING=ON -DBUILD_WITH_PYTHON=ON -DCMAKE_CXX_FLAGS="$CXX_FLAGS" ..
-    - make -j`nproc`
+    - make -j`nproc` install
+    - if [ -n "$RUN_TESTS" ]; then ctest -j`nproc` --output-on-failure; fi
 
-build-everybeam-2004:
-  extends: [.needs-base-2004,.build-everybeam]
+build-test-2004:
+  extends: [ .needs-base-2004, .build-everybeam, .download-wsrt-measures ]
   variables:
     CXX_FLAGS: -mavx2 -mfma
+    RUN_TESTS: yes please run all tests with AVX support
 
 build-everybeam-2204:
   extends: [ .needs-base-2204, .build-everybeam ]
@@ -212,7 +215,7 @@ build-doc-2204:
     - build/doc/html
 
 build-package-2204:
-  extends: [.needs-base-2204]
+  extends: .needs-base-2204
   stage: package
   variables:
     GIT_SUBMODULE_STRATEGY: normal
@@ -236,7 +239,7 @@ build-package-2204:
 
 build-compare-oskar-2204:
   stage: oskar-comparison
-  extends: [.needs-base-2204]
+  extends: .needs-base-2204
   image: $BASE_IMAGE_2204
   variables:
     GIT_SUBMODULE_STRATEGY: normal
@@ -302,7 +305,7 @@ deploy-image-2204:
       - dist/everybeam*.whl
 
 cibuild-python-wheels:
-  extends: [ ".needs-base-2204", ".cibuildwheels-custom" ]
+  extends: [ .needs-base-2204, .cibuildwheels-custom ]
   stage: deploy
   image: python:3.8
   variables:

cpp/beamformer.cc

@@ -84,7 +84,7 @@ std::vector<aocommon::MC2x2Diag> BeamFormer::ComputeWeightedResponses(
       ComputeGeometricResponse(delta_phase_reference_positions_, pointing);
 
   // Initialize and fill result
-  double weight_sum[2] = {0.0, 0.0};
+  std::array<double, 2> weight_sum = {0.0, 0.0};
   std::vector<aocommon::MC2x2Diag> result(antennas_.size());
   for (size_t idx = 0; idx < antennas_.size(); ++idx) {
     // Get geometric response at index
@@ -98,8 +98,7 @@ std::vector<aocommon::MC2x2Diag> BeamFormer::ComputeWeightedResponses(
 
   // Normalize the weight by the number of antennas
   for (auto& entry : result) {
-    entry[0] /= weight_sum[0];
-    entry[1] /= weight_sum[1];
+    entry = {entry[0] / weight_sum[0], entry[1] / weight_sum[1]};
   }
   return result;
 }

cpp/beamformerlofar.cc

@@ -27,21 +27,21 @@ aocommon::MC2x2Diag BeamFormerLofar::FieldArrayFactor(
   aocommon::UVector<std::complex<double>> geometric_response =
       BeamFormer::ComputeGeometricResponse(antenna_positions, delta_direction);
 
-  double weight_sum[2] = {0.0, 0.0};
-  aocommon::MC2x2Diag result(0., 0.);
+  std::array<std::complex<double>, 2> response_sum = {0.0, 0.0};
+  std::array<double, 2> weight_sum = {0.0, 0.0};
 
   for (size_t idx = 0; idx < antenna_positions.size(); ++idx) {
-    result[0] += geometric_response[idx] * (1.0 * antenna_enabled[idx][0]);
-    result[1] += geometric_response[idx] * (1.0 * antenna_enabled[idx][1]);
+    response_sum[0] +=
+        geometric_response[idx] * (1.0 * antenna_enabled[idx][0]);
+    response_sum[1] +=
+        geometric_response[idx] * (1.0 * antenna_enabled[idx][1]);
     weight_sum[0] += (1.0 * antenna_enabled[idx][0]);
     weight_sum[1] += (1.0 * antenna_enabled[idx][1]);
   }
 
   // Normalize the weight by the number of enabled tiles
-  result[0] /= weight_sum[0];
-  result[1] /= weight_sum[1];
-
-  return result;
+  return aocommon::MC2x2Diag(response_sum[0] / weight_sum[0],
+                             response_sum[1] / weight_sum[1]);
 }
 
 aocommon::MC2x2 BeamFormerLofar::LocalResponse(

cpp/hamaker/hamakerelementresponse.cc

@@ -99,10 +99,8 @@ aocommon::MC2x2 HamakerElementResponse::Response(double freq, double theta,
     const double caz = std::cos(angle);
     const double saz = std::sin(angle);
 
-    response[0] += caz * P.first;
-    response[1] += -saz * P.second;
-    response[2] += saz * P.first;
-    response[3] += caz * P.second;
+    response += aocommon::MC2x2(caz * P.first, -saz * P.second, saz * P.first,
+                                caz * P.second);
 
     // Update sign and kappa.
     sign = -sign;

cpp/pointresponse/dishpoint.cc

@@ -86,15 +86,7 @@ aocommon::MC2x2 DishPoint::Response(BeamMode beam_mode, size_t station_idx,
 
   vp.Render(buffer, pointing_direction_angle, 0.0, 0.0, 0.0, freq);
 
-  // Because of the conversion from std::complex<float> to std::complex<double>
-  // It is not possible to directly return the buffer by a
-  //    return aocommon::MC2x2(buffer);
-  // Instead, cast the buffer element by element before returning
-  aocommon::MC2x2 result;
-  for (int i = 0; i < 4; i++) {
-    result[i] = buffer[i];
-  }
-  return result;
+  return aocommon::MC2x2(aocommon::MC2x2F(buffer));
 }
 
 void DishPoint::ResponseAllStations(BeamMode beam_mode,

cpp/sphericalharmonicsresponse.cc

@@ -220,10 +220,10 @@ aocommon::MC2x2 SphericalHarmonicsResponse::ComputeResponse(
     std::complex<double> q3;
     std::tie(q2, q3) =
         common::F4far_new(nms_(i, 2), nms_(i, 1), nms_(i, 0), theta, phi);
-    response[0] += q2 * c0;  // xx
-    response[1] += q3 * c0;  // xy
-    response[2] += q2 * c1;  // yx
-    response[3] += q3 * c1;  // yy
+
+    //                                         xx, xy, yx, yy
+    response += ElementProduct(aocommon::MC2x2(q2, q3, q2, q3),
+                               aocommon::MC2x2(c0, c0, c1, c1));
   }
 
   return response;

cpp/sphericalharmonicsresponsefixeddirection.cc

@@ -66,10 +66,10 @@ aocommon::MC2x2 SphericalHarmonicsResponseFixedDirection::ComputeResponse(
         coefficients(1, frequency_index, element_index, i);
     const std::complex<double> q2 = base_functions_(i, 0);
     const std::complex<double> q3 = base_functions_(i, 1);
-    response[0] += q2 * c0;  // xx
-    response[1] += q3 * c0;  // xy
-    response[2] += q2 * c1;  // yx
-    response[3] += q3 * c1;  // yy
+
+    //                                         xx, xy, yx, yy
+    response += ElementProduct(aocommon::MC2x2(q2, q3, q2, q3),
+                               aocommon::MC2x2(c0, c0, c1, c1));
   }
 
   return response;

cpp/test/tatermconfig.cc

@@ -228,17 +228,14 @@ BOOST_AUTO_TEST_CASE(combine_aterms) {
                  n_channels, n_times, coord_system);
 
   // Combine the aterms generated in the beam-only and diagonal-only modes
-  std::vector<std::complex<float>> aterms_product(aterm_buffer_size);
   for (size_t j = 0; j != aterm_buffer_size; j += kAtermMatrixSize) {
-    aocommon::MC2x2F scratch =
+    const aocommon::MC2x2F product =
         aocommon::MC2x2F(aterm_buffer_diagonal.data() + j) *
         aocommon::MC2x2F(aterm_buffer_beam.data() + j);
-    scratch.AssignTo(aterms_product.data() + j);
+    for (size_t i = 0; i < 4; ++i) {
+      BOOST_CHECK_CLOSE(product[i], aterm_buffer_combined[j + i], 1.0e-4);
+    }
   }
-
-  BOOST_CHECK_EQUAL_COLLECTIONS(aterms_product.begin(), aterms_product.end(),
-                                aterm_buffer_combined.begin(),
-                                aterm_buffer_combined.end());
 }
 
 BOOST_AUTO_TEST_CASE(fourierfit) {

cpp/test/tdish.cc

@@ -45,7 +45,7 @@ BOOST_AUTO_TEST_CASE(load_dish) {
       dynamic_cast<everybeam::pointresponse::DishPoint*>(point_response.get());
   BOOST_REQUIRE(dish_point_response);
 
-  std::vector<std::vector<std::complex<float>>> reference_response = {
+  const std::vector<std::vector<std::complex<float>>> kReferenceResponse = {
       {{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}},
       {{0.382599, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.382599, 0.0}}};
 
@@ -64,9 +64,8 @@ BOOST_AUTO_TEST_CASE(load_dish) {
                                   kFrequency, station_id, field_id);
 
     for (std::size_t i = 0; i < 4; ++i) {
-      std::cout << point_response_buffer[i] << std::endl;
-      BOOST_CHECK(
-          std::abs(point_response_buffer[i] - reference_response[j][i]) < 1e-6);
+      BOOST_CHECK_CLOSE(point_response_buffer[i], kReferenceResponse[j][i],
+                        2.0e-4);
     }
 
     casacore::MDirection pointing(casacore::Quantity(ra, "rad"),
@@ -76,13 +75,12 @@ BOOST_AUTO_TEST_CASE(load_dish) {
     const coords::ItrfConverter itrf_converter(kTime);
     vector3r_t direction = itrf_converter.ToItrf(pointing);
 
-    aocommon::MC2x2 response = dish_point_response->Response(
+    const aocommon::MC2x2 response = dish_point_response->Response(
         everybeam::BeamMode::kFull, station_id, kFrequency, direction);
 
     for (std::size_t i = 0; i < 4; ++i) {
-      std::cout << response[i] << "  " << reference_response[j][i] << std::endl;
-      BOOST_CHECK(std::abs(response[i] - std::complex<double>(
-                                             reference_response[j][i])) < 1e-6);
+      BOOST_CHECK_CLOSE(response[i],
+                        std::complex<double>(kReferenceResponse[j][i]), 2.0e-4);
     }
   }
 }

cpp/test/tlofar_lba.cc

@@ -30,19 +30,26 @@ using everybeam::griddedresponse::PhasedArrayGrid;
 using everybeam::telescope::LOFAR;
 using everybeam::telescope::Telescope;
 
-BOOST_AUTO_TEST_SUITE(tlofar_lba)
+BOOST_AUTO_TEST_SUITE(lofar_lba)
 
+namespace {
 // Properties extracted from MS
-double time = 4.92183348e+09;
-double frequency = 57812500.;
-double ra(-1.44194878), dec(0.85078091);
+const double kTime = 4.92183348e+09;
+const double kFrequency = 57812500.0;
+const double kRa = -1.44194878;
+const double kDec = 0.85078091;
 
 // Properties of grid
-std::size_t width(4), height(4);
-double dl(0.5 * M_PI / 180.), dm(0.5 * M_PI / 180.), shift_l(0.), shift_m(0.);
+const std::size_t kWidth = 4;
+const std::size_t kHeight = 4;
+const double kDl = 0.5 * M_PI / 180.0;
+const double kDm = 0.5 * M_PI / 180.0;
+const double kShiftL = 0.0;
+const double kShiftM = 0.0;
 
-CoordinateSystem coord_system = {width, height, ra,      dec,
-                                 dl,    dm,     shift_l, shift_m};
+const CoordinateSystem kCoordinateSystem = {kWidth, kHeight, kRa,     kDec,
+                                            kDl,    kDm,     kShiftL, kShiftM};
+}  // namespace
 
 BOOST_AUTO_TEST_CASE(test_hamaker) {
   Options options;
@@ -71,7 +78,7 @@ BOOST_AUTO_TEST_CASE(test_hamaker) {
 
   const Station& station = lofartelescope.GetStation(19);
   aocommon::MC2x2 element_response =
-      station.ComputeElementResponse(time, frequency, direction, false, true);
+      station.ComputeElementResponse(kTime, kFrequency, direction, false, true);
 
   for (size_t i = 0; i != 4; ++i) {
     BOOST_CHECK(std::abs(element_response[i] - target_element_response[i]) <
@@ -86,7 +93,7 @@ BOOST_AUTO_TEST_CASE(test_hamaker) {
   vector3r_t direction_s31 = {0.667806, -0.0770635, 0.740335};
   vector3r_t station0_dir = {0.655743, -0.0670973, 0.751996};
   aocommon::MC2x2 station31_response = station31.Response(
-      time, freq4, direction_s31, freq4, station0_dir, station0_dir);
+      kTime, freq4, direction_s31, freq4, station0_dir, station0_dir);
 
   aocommon::MC2x2 target_station_response(
       {-0.71383788, 0.00612506}, {-0.4903527, 0.00171652},
@@ -99,7 +106,7 @@ BOOST_AUTO_TEST_CASE(test_hamaker) {
 
   // Gridded response
   std::unique_ptr<GriddedResponse> grid_response =
-      telescope->GetGriddedResponse(coord_system);
+      telescope->GetGriddedResponse(kCoordinateSystem);
   BOOST_CHECK(nullptr != dynamic_cast<PhasedArrayGrid*>(grid_response.get()));
 
   // Define buffer and get gridded responses
@@ -107,7 +114,8 @@ BOOST_AUTO_TEST_CASE(test_hamaker) {
       grid_response->GetStationBufferSize(1));
 
   grid_response->Response(everybeam::BeamMode::kFull,
-                          antenna_buffer_single.data(), time, frequency, 31, 0);
+                          antenna_buffer_single.data(), kTime, kFrequency, 31,
+                          0);
 
   // Compare with everybeam at pixel (1, 3), reference solution obtained with
   // everybeam at commit sha 70a286e7dace4616417b0e973a624477f15c9ce3
@@ -117,10 +125,10 @@ BOOST_AUTO_TEST_CASE(test_hamaker) {
       {-0.541472, 0.0053149},
       {0.775696, -0.00369905}};
 
-  std::size_t offset_13 = (3 + 1 * width) * 4;
+  std::size_t offset_13 = (3 + 1 * kWidth) * 4;
   for (std::size_t i = 0; i < 4; ++i) {
     BOOST_CHECK(std::abs(antenna_buffer_single[offset_13 + i] -
-                         everybeam_ref_p13[i]) < 1e-6);
+                         everybeam_ref_p13[i]) < 1.0e-6);
   }
 }
 
@@ -142,7 +150,7 @@ BOOST_AUTO_TEST_CASE(test_lobes) {
 
   // Gridded response
   std::unique_ptr<GriddedResponse> grid_response =
-      telescope->GetGriddedResponse(coord_system);
+      telescope->GetGriddedResponse(kCoordinateSystem);
 
   // Define buffer and get gridded responses
   std::vector<std::complex<float>> antenna_buffer_single(
@@ -150,7 +158,8 @@ BOOST_AUTO_TEST_CASE(test_lobes) {
 
   // Get the gridded response for station 20 (of course!)
   grid_response->Response(everybeam::BeamMode::kFull,
-                          antenna_buffer_single.data(), time, frequency, 20, 0);
+                          antenna_buffer_single.data(), kTime, kFrequency, 20,
+                          0);
 
   // Compare with everybeam at pixel (1, 3). This solution only is a "reference"
   // certainly not a "ground-truth"
@@ -160,11 +169,11 @@ BOOST_AUTO_TEST_CASE(test_lobes) {
       {0.0086051673, 0.034685627},
       {1.5801408, 2.0045171}};
 
-  std::size_t offset_13 = (3 + 1 * width) * 4;
+  std::size_t offset_13 = (3 + 1 * kWidth) * 4;
   for (std::size_t i = 0; i < 4; ++i) {
-    // relative tolerance is 2e-1%, so 2e-3
+    // relative tolerance is 5.0e-1%, so 5.0e-3
     BOOST_CHECK_CLOSE(antenna_buffer_single[offset_13 + i],
-                      everybeam_ref_p13[i], 2e-1);
+                      everybeam_ref_p13[i], 5.0e-1);
   }
 }
 

docker/ubuntu_20_04_base

@@ -30,6 +30,8 @@ RUN export DEBIAN_FRONTEND="noninteractive" && \
 		wget \
 	&& \
 	rm -rf /var/lib/apt/lists/*
+# pytest-lazy-fixture does not work with pytest-8.x -> Use pytest-7.x for now.
+# TODO(AST-1501): Investigate removing pytest-lazy-fixture.
 RUN pip3 install \
 		aptly-api-client \
 		astropy \
@@ -42,6 +44,8 @@ RUN pip3 install \
 		myst-parser \
 		numpy==1.19.0 \
 		pandas \
+		pytest~=7.0 \
+		pytest-lazy-fixture \
 		scipy \
 		sphinx \
 		sphinx_rtd_theme \