Merge branch 'port_direction_related_routines' into 'main'

Initial porting for the radec2lmn

See merge request mancini/predict!20

.gitlab-ci.yml

@@ -149,5 +149,7 @@ collect-performance:
   - ls -la
   - python3 ci/summarize-results.py --filter PredictBenchmark/PointSource results*.json result-summary-pointsource
   - python3 ci/summarize-results.py --filter PredictBenchmark/GaussianSource results*.json result-summary-gaussian
+  - python3 ci/summarize-results.py --filter DirectionsBenchmark/Directions results*.json result-summary-directions
+  
   
 

CMakeLists.txt

@@ -32,8 +32,14 @@ file(GLOB_RECURSE HEADERS include/*.h)
 add_library(predict STATIC ${SOURCES})
 
 target_compile_options(
-  predict PRIVATE $<$<CONFIG:Debug>:-g> $<$<CONFIG:Release>:-O3>
-                  $<$<CONFIG:RelWithDebInfo>:-O3 -g>)
+  predict
+  PRIVATE $<$<CONFIG:Debug>:-g
+          -O3>
+          $<$<CONFIG:Release>:-O3
+          -march=native>
+          $<$<CONFIG:RelWithDebInfo>:-O3
+          -g
+          -march=native>)
 target_link_options(predict PRIVATE $<$<CONFIG:Debug>:-g>
                     $<$<CONFIG:Release>:-O3> $<$<CONFIG:RelWithDebInfo>:-O3 -g>)
 
@@ -80,8 +86,7 @@ if(PREDICT_BUILD_TESTING)
     get_filename_component(test_name ${test} NAME_WE)
     string(REPLACE "test_" "" test_name ${test_name})
     message(STATUS "Adding test: ${test_name}")
-    add_test(NAME ${test_name}_unittests COMMAND ${test_name}_unittests
-                                                 ${test_name})
+    add_test(NAME ${test_name}_unittests COMMAND ${test_name}_unittests)
     add_executable(${test_name}_unittests ${test})
 
     target_include_directories(${test_name}_unittests

benchmark/CMakeLists.txt

@@ -14,7 +14,7 @@ FetchContent_MakeAvailable(googlebenchmark)
 
 set(BENCHMARKS simulator)
 
-add_executable(microbenchmarks main.cpp predict.cpp)
+add_executable(microbenchmarks main.cpp predict.cpp directions.cpp)
 
 target_compile_options(
   microbenchmarks PRIVATE $<$<CONFIG:Debug>:-g> $<$<CONFIG:Release>:-O3>

benchmark/directions.cpp

+#include <Simulator.h>
+#include <benchmark/benchmark.h>
+#include <random>
+#include <xtensor/xtensor.hpp>
+
+#include <Directions.h>
+#include <test/Common.h>
+
+extern bool do_randomized_run;
+extern int randomized_run_seed;
+
+const std::vector<int64_t> nstations = {24, 48};
+const std::vector<int64_t> nsources = {64, 128};
+
+namespace {
+Direction test_reference_point{0.5, 0.1};
+Direction test_offset_point{test_reference_point.ra + 0.02,
+                            test_reference_point.dec + 0.02};
+} // namespace
+
+class DirectionsBenchmark : public benchmark::Fixture {
+public:
+  void SetUp(benchmark::State &state) override {
+    n_directions = state.range(0);
+    lmn = xt::xtensor<double, 2>({n_directions, 3});
+    directions = std::make_unique<Directions>();
+    directions->reserve(n_directions);
+    for (size_t i = 0; i < n_directions; i++) {
+      if (do_randomized_run) {
+        SetUpRandomizedSource(test_reference_point, test_offset_point,
+                              randomized_run_seed);
+      } else {
+        SetUpFixedSource(test_reference_point, test_offset_point);
+      }
+      directions->add(test_offset_point);
+    }
+  }
+  void TearDown(benchmark::State &) override { directions.reset(); }
+
+protected:
+  std::unique_ptr<Directions> directions;
+  size_t n_directions;
+  xt::xtensor<double, 2> lmn;
+};
+
+BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsOneByOne)
+(benchmark::State &state) {
+  for (auto _ : state) {
+    directions->radec2lmn<Directions::computation_strategy::SINGLE>(
+        test_reference_point, lmn);
+  }
+}
+
+BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsMulti)
+(benchmark::State &state) {
+  for (auto _ : state) {
+    directions->radec2lmn<Directions::computation_strategy::MULTI>(
+        test_reference_point, lmn);
+  }
+}
+
+BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsMultiSIMD)
+(benchmark::State &state) {
+  for (auto _ : state) {
+    directions->radec2lmn<Directions::computation_strategy::MULTI_SIMD>(
+        test_reference_point, lmn);
+  }
+}
+
+BENCHMARK_REGISTER_F(DirectionsBenchmark, DirectionsOneByOne)
+    ->ArgsProduct({{1024, 2048, 4096, 8192}})
+    ->ArgNames({
+        "#dir",
+    })
+    ->Unit(benchmark::kMillisecond);
+
+BENCHMARK_REGISTER_F(DirectionsBenchmark, DirectionsMulti)
+    ->ArgsProduct({{1024, 2048, 4096, 8192}})
+    ->ArgNames({
+        "#dir",
+    })
+    ->Unit(benchmark::kMillisecond);
+
+BENCHMARK_REGISTER_F(DirectionsBenchmark, DirectionsMultiSIMD)
+    ->ArgsProduct({{1024, 2048, 4096, 8192}})
+    ->ArgNames({
+        "#dir",
+    })
+    ->Unit(benchmark::kMillisecond);
\ No newline at end of file

benchmark/main.cpp

@@ -10,7 +10,9 @@ void ParseCustomFlags(int argc, char **argv) {
   extern int randomized_run_seed;
 
   for (int i = 0; i < argc; ++i) {
-    do_randomized_run = (std::string(argv[i]) == "--randomize");
+    if (std::string(argv[i]) == std::string("--randomize")) {
+      do_randomized_run = true;
+    }
 
     if (std::string(argv[i]) == "--seed") {
       if (i + 1 < argc) {

ci/das6/compile_and_run.sh

@@ -10,7 +10,7 @@ COMPILER_VERSION=$2
 echo "Running on compiler version: ${COMPILER_VERSION}, architecture: ${ARCHITECTURE}, hostname: $(hostname)"
 BUILD_DIR=build-${COMPILER_VERSION}-${ARCHITECTURE}
 module purge
-module load spack/20250403/${COMPILER_VERSION}
+module load spack/20250403
 module load cmake casacore boost
 
 cmake -B ${BUILD_DIR} . -DCMAKE_BUILD_TYPE=Release -DPREDICT_BUILD_BENCHMARK=ON

ci/summarize-results.py

@@ -35,12 +35,17 @@ def store_combined(outfile, obj):
 
 def create_summary_plot(metrics, outplot_name):
     time_unit = metrics.time_unit[0]
+    metrics["test_name"] = metrics["name"].str.split("/", n=2).str[1] 
+    metrics["parameter"] = metrics["name"].str.split("/", n=2).str[2]
     
-    grid = seaborn.FacetGrid(metrics, col="architecture")
-    grid.map(seaborn.barplot, "compiler_version", "cpu_time", "name")
+    grid = seaborn.FacetGrid(metrics, col="architecture", row="compiler_version")
+    grid.map(seaborn.scatterplot, "parameter", "cpu_time", "test_name")
     grid.set_ylabels(f"CPU time({time_unit})")
+    grid.set_xlabels("Parameter")
     grid.add_legend()
     plt.setp(grid._legend.get_texts(), fontsize=9)
+    plt.xticks(rotation=90)
+    
     grid.savefig(outplot_name + ".png")
 
 def main():

include/Directions.h

@@ -12,12 +12,16 @@
 #include <Direction.h>
 #include <ObjectCollection.h>
 #include <common/SmartVector.h>
+
 #include <vector>
+#include <xtensor/xmath.hpp>
+#include <xtensor/xview.hpp>
 
 using Direction = dp3::base::Direction;
 
 class Directions : ObjectCollection<Direction> {
 public:
+  enum computation_strategy { SINGLE, MULTI, MULTI_SIMD };
   Directions() : ra(), dec() {}
   void add(const Direction &direction) {
     ra.push_back(direction.ra);
@@ -47,5 +51,153 @@ public:
 
   SmartVector<double> ra;
   SmartVector<double> dec;
+
+  /**
+   * Compute LMN coordinates of \p direction relative to \p reference.
+   * \param[in]   reference
+   * Reference direction on the celestial sphere.
+   * \param[in]   direction
+   * Direction of interest on the celestial sphere.
+   * \param[out]   lmn
+   * Pointer to a buffer of (at least) length three into which the computed LMN
+   * coordinates will be written.
+   */
+  template <computation_strategy T = computation_strategy::MULTI_SIMD>
+  inline void radec2lmn(const Direction &reference,
+                        xt::xtensor<double, 2> &lmn);
+};
+
+template <>
+inline void Directions::radec2lmn<Directions::computation_strategy::MULTI>(
+    const Direction &reference, xt::xtensor<double, 2> &lmn) {
+  /**
+   * \f{eqnarray*}{
+   *   \ell &= \cos(\delta) \sin(\alpha - \alpha_0) \\
+   *      m &= \sin(\delta) \cos(\delta_0) - \cos(\delta) \sin(\delta_0)
+   *                                         \cos(\alpha - \alpha_0)
+   * \f}
+   */
+
+  const auto ra_view = ra.view();
+  const auto dec_view = dec.view();
+
+  const xt::xtensor<double, 1> delta_ra = ra_view - reference.ra;
+  const xt::xtensor<double, 1> sin_delta_ra = xt::sin(delta_ra);
+  const xt::xtensor<double, 1> cos_delta_ra = xt::cos(delta_ra);
+  const xt::xtensor<double, 1> sin_dec = xt::sin(dec_view);
+  const xt::xtensor<double, 1> cos_dec = xt::cos(dec_view);
+  const double sin_dec0 = std::sin(reference.dec);
+  const double cos_dec0 = std::cos(reference.dec);
+
+  xt::view(lmn, xt::all(), 0) = cos_dec * sin_delta_ra;
+  xt::view(lmn, xt::all(), 1) =
+      sin_dec * cos_dec0 - cos_dec * sin_dec0 * cos_delta_ra;
+  // Normally, n is calculated using n = std::sqrt(1.0 - l * l - m * m).
+  // However the sign of n is lost that way, so a calculation is used that
+  // avoids losing the sign. This formula can be found in Perley (1989).
+  // Be aware that we asserted that the sign is wrong in Perley (1989),
+  // so this is corrected.
+  xt::view(lmn, xt::all(), 2) =
+      sin_dec * sin_dec0 + cos_dec * cos_dec0 * cos_delta_ra;
+}
+
+template <>
+inline void Directions::radec2lmn<Directions::computation_strategy::SINGLE>(
+    const Direction &reference, xt::xtensor<double, 2> &lmn) {
+  /**
+   * \f{eqnarray*}{
+   *   \ell &= \cos(\delta) \sin(\alpha - \alpha_0) \\
+   *      m &= \sin(\delta) \cos(\delta_0) - \cos(\delta) \sin(\delta_0)
+   *                                         \cos(\alpha - \alpha_0)
+   * \f}
+   */
+  for (size_t i = 0; i < ra.size(); i++) {
+    const double delta_ra = ra[i] - reference.ra;
+    const double sin_delta_ra = std::sin(delta_ra);
+    const double cos_delta_ra = std::cos(delta_ra);
+    const double sin_dec = std::sin(dec[i]);
+    const double cos_dec = std::cos(dec[i]);
+    const double sin_dec0 = std::sin(reference.dec);
+    const double cos_dec0 = std::cos(reference.dec);
+
+    lmn(i, 0) = cos_dec * sin_delta_ra;
+    lmn(i, 1) = sin_dec * cos_dec0 - cos_dec * sin_dec0 * cos_delta_ra;
+    // Normally, n is calculated using n = std::sqrt(1.0 - l * l - m * m).
+    // However the sign of n is lost that way, so a calculation is used that
+    // avoids losing the sign. This formula can be found in Perley (1989).
+    // Be aware that we asserted that the sign is wrong in Perley (1989),
+    // so this is corrected.
+    lmn(i, 2) = sin_dec * sin_dec0 + cos_dec * cos_dec0 * cos_delta_ra;
+  }
+}
+
+#include <cstddef>
+#include <xsimd/xsimd.hpp>
+
+namespace xsimd {
+
+struct radec2lmn {
+  template <class C, class Tag, class Arch>
+  void operator()(Arch, const Direction &reference, const C &ra, const C &dec,
+                  xt::xtensor<double, 2> &lmn, Tag);
 };
+
+template <class C, class Tag, class Arch>
+void radec2lmn::operator()(Arch, const Direction &reference, const C &ra,
+                           const C &dec, xt::xtensor<double, 2> &lmn, Tag) {
+  using b_type = xsimd::batch<double, Arch>;
+  std::size_t inc = b_type::size;
+  std::size_t size = ra.size();
+  double sin_dec0 = std::sin(reference.dec);
+  double cos_dec0 = std::cos(reference.dec);
+  // size for which the vectorization is possible
+  std::size_t vec_size = size - size % inc;
+  for (std::size_t i = 0; i < vec_size; i += inc) {
+    const b_type ra_vec = b_type::load(&ra[i], Tag());
+    const b_type dec_vec = b_type::load(&dec[i], Tag());
+
+    const b_type delta_ra = ra_vec - reference.ra;
+
+    const std::pair<b_type, b_type> sin_cos_delta_ra = xsimd::sincos(delta_ra);
+    const std::pair<b_type, b_type> sin_cos_dec = xsimd::sincos(dec_vec);
+
+    const b_type &cos_dec = sin_cos_dec.second;
+    const b_type &sin_dec = sin_cos_dec.first;
+    const b_type &sin_delta_ra = sin_cos_delta_ra.first;
+    const b_type &cos_delta_ra = sin_cos_delta_ra.second;
+
+    const b_type result_vec_l = cos_dec * sin_delta_ra;
+    const b_type result_vec_m =
+        sin_dec * cos_dec0 - cos_dec * sin_dec0 * cos_delta_ra;
+    const b_type result_vec_n =
+        sin_dec * sin_dec0 + cos_dec * cos_dec0 * cos_delta_ra;
+
+    xsimd::store(&lmn(0, i), result_vec_l, Tag());
+    xsimd::store(&lmn(1, i), result_vec_m, Tag());
+    xsimd::store(&lmn(2, i), result_vec_n, Tag());
+  }
+  // Remaining part that cannot be vectorize
+  for (std::size_t i = vec_size; i < size; ++i) {
+    double delta_ra = ra[i] - reference.ra;
+    double sin_delta_ra = std::sin(delta_ra);
+    double cos_delta_ra = std::cos(delta_ra);
+    double sin_dec = std::sin(dec[i]);
+    double cos_dec = std::cos(dec[i]);
+
+    lmn(i, 0) = cos_dec * sin_delta_ra;
+    lmn(i, 1) = sin_dec * cos_dec0 - cos_dec * sin_dec0 * cos_delta_ra;
+    lmn(i, 2) = sin_dec * sin_dec0 + cos_dec * cos_dec0 * cos_delta_ra;
+  }
+}
+
+} // namespace xsimd
+
+template <>
+inline void Directions::radec2lmn<Directions::computation_strategy::MULTI_SIMD>(
+    const Direction &reference, xt::xtensor<double, 2> &lmn) {
+  xt::xtensor<double, 2> lmn_tmp({3, ra.size()});
+  xsimd::dispatch(xsimd::radec2lmn{})(reference, ra, dec, lmn_tmp,
+                                      xsimd::aligned_mode());
+  lmn = xt::transpose(lmn_tmp);
+}
 #endif
\ No newline at end of file

tests/test_directions.cpp

@@ -3,7 +3,9 @@
 
 #define BOOST_TEST_MODULE DIRECTIONS_TEST
 #include <Directions.h>
+#include <Simulator.h>
 #include <boost/test/unit_test.hpp>
+#include <xsimd/xsimd.hpp>
 #include <xtensor/xadapt.hpp>
 
 BOOST_AUTO_TEST_SUITE(DirectionsTestSuite)
@@ -26,10 +28,6 @@ BOOST_AUTO_TEST_CASE(test_add_multi_direction) {
 
   BOOST_CHECK_EQUAL(directions.ra[1], 0.3);
   BOOST_CHECK_EQUAL(directions.dec[1], 0.4);
-
-  std::vector<double> ra{0.60, 0.70};
-
-  auto bla = xt::adapt(ra.data(), ra.size(), xt::no_ownership());
 }
 
 BOOST_AUTO_TEST_CASE(test_reserve) {
@@ -46,4 +44,74 @@ BOOST_AUTO_TEST_CASE(test_reserve) {
   auto bla = xt::adapt(ra.data(), ra.size(), xt::no_ownership());
 }
 
+BOOST_AUTO_TEST_CASE(test_radec_to_lmn_single) {
+  Direction reference(0.2, 0.3);
+  Direction dir_0(0.1, 0.2);
+  Direction dir_1(0.2, 0.3);
+
+  xt::xtensor<double, 2> lmn_expected({2, 3});
+
+  dp3::base::radec2lmn(reference, dir_0, &lmn_expected(0, 0));
+  dp3::base::radec2lmn(reference, dir_1, &lmn_expected(1, 0));
+
+  Directions directions;
+  directions.add(dir_0);
+  directions.add(dir_1);
+
+  xt::xtensor<double, 2> lmn_actual({2, 3});
+  directions.radec2lmn<Directions::SINGLE>(reference, lmn_actual);
+  for (size_t i; i < 2; ++i) {
+    BOOST_CHECK_CLOSE(lmn_expected(i, 0), lmn_actual(i, 0), 1.e-6);
+    BOOST_CHECK_CLOSE(lmn_expected(i, 1), lmn_actual(i, 1), 1.e-6);
+    BOOST_CHECK_CLOSE(lmn_expected(i, 2), lmn_actual(i, 2), 1.e-6);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(test_radec_to_lmn_multi) {
+  Direction reference(0.2, 0.3);
+  Direction dir_0(0.1, 0.2);
+  Direction dir_1(0.2, 0.3);
+
+  xt::xtensor<double, 2> lmn_expected({2, 3});
+
+  dp3::base::radec2lmn(reference, dir_0, &lmn_expected(0, 0));
+  dp3::base::radec2lmn(reference, dir_1, &lmn_expected(1, 0));
+
+  Directions directions;
+  directions.add(dir_0);
+  directions.add(dir_1);
+
+  xt::xtensor<double, 2> lmn_actual({2, 3});
+  directions.radec2lmn<Directions::MULTI>(reference, lmn_actual);
+  for (size_t i; i < 2; ++i) {
+    BOOST_CHECK_CLOSE(lmn_expected(i, 0), lmn_actual(i, 0), 1.e-6);
+    BOOST_CHECK_CLOSE(lmn_expected(i, 1), lmn_actual(i, 1), 1.e-6);
+    BOOST_CHECK_CLOSE(lmn_expected(i, 2), lmn_actual(i, 2), 1.e-6);
+  }
+}
+
+BOOST_AUTO_TEST_CASE(test_radec_to_lmn_xsimd) {
+  Direction reference(0.2, 0.3);
+  Direction dir_0(0.1, 0.2);
+  Direction dir_1(0.2, 0.3);
+
+  xt::xtensor<double, 2> lmn_expected({2, 3});
+
+  dp3::base::radec2lmn(reference, dir_0, &lmn_expected(0, 0));
+  dp3::base::radec2lmn(reference, dir_1, &lmn_expected(1, 0));
+
+  Directions directions;
+  directions.add(dir_0);
+  directions.add(dir_1);
+
+  xt::xtensor<double, 2> lmn_actual({2, 3});
+  directions.radec2lmn<Directions::computation_strategy::MULTI_SIMD>(
+      reference, lmn_actual);
+  for (size_t i; i < 2; ++i) {
+    BOOST_CHECK_CLOSE(lmn_expected(i, 0), lmn_actual(i, 0), 1.e-6);
+    BOOST_CHECK_CLOSE(lmn_expected(i, 1), lmn_actual(i, 1), 1.e-6);
+    BOOST_CHECK_CLOSE(lmn_expected(i, 2), lmn_actual(i, 2), 1.e-6);
+  }
+}
+
 BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file

tests/test_gaussiansourcecollection.cpp

@@ -9,8 +9,8 @@
 BOOST_AUTO_TEST_SUITE(GaussianSourceCollectionTestSuite)
 
 BOOST_AUTO_TEST_CASE(add_single) {
-  PointSource source(Direction(0.1, 0.2), Stokes(1.0, 0.0, 0.0, 0.0));
-  PointSourceCollection collection;
+  GaussianSource source(Direction(0.1, 0.2), Stokes(1.0, 0.0, 0.0, 0.0));
+  GaussianSourceCollection collection;
   collection.add(source);
   BOOST_CHECK_EQUAL(collection.size(), 1);
   BOOST_CHECK_EQUAL(collection.direction_vector.size(), 1);