Merge branch 'refactor-to-simulatorplan' into 'main'

Refactor code base into plan-based implementation

Closes #19, #22, and #25

See merge request !25

.gitlab-ci.yml

@@ -13,7 +13,7 @@ stages:
   image: ubuntu:22.04
   before_script:
     - apt-get update -qq
-    - apt-get install -y casacore-dev cmake gcc gcovr git g++ libboost-test-dev make python3 python3-pip
+    - apt-get install -y casacore-dev libblas-dev liblapack-dev cmake gcc gcovr git g++ libboost-test-dev make python3 python3-pip
     - python3 -m pip install --upgrade pip
     - python3 -m pip install clang-format==14.0.0 cmake-format pre-commit clang-tidy seaborn
 
@@ -24,6 +24,7 @@ stages:
     - module load python
     - module load py-pip
     - module load boost
+    - module load openblas
 
 # --------- Precheck: Code formatting ----------
 format:
@@ -150,6 +151,7 @@ collect-performance:
   - 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
+  - python3 ci/summarize-results.py --filter PhasesBenchmark/ComputePhases results*.json result-summary-phases
   - python3 ci/summarize-results.py --filter SpectrumBenchmark/Spectrum results*.json result-summary-spectrum
   
   

CMakeLists.txt

@@ -22,32 +22,31 @@ set(CMAKE_MODULE_PATH ${aocommon_SOURCE_DIR}/CMake)
 set(XTENSOR_LIBRARIES xtl xsimd xtensor xtensor-blas)
 include(${aocommon_SOURCE_DIR}/CMake/FetchXTensor.cmake)
 
+set(PREDICT_BUILD_ARGUMENTS
+    $<$<CONFIG:Debug>:-g> $<$<CONFIG:Release>:-O3;-march=native >
+    $<$<CONFIG:RelWithDebInfo>:-O3;-g;-march=native >)
+
 # Find or fetch dependencies
 find_package(Casacore REQUIRED casa tables)
 
+# Find BLAS
+find_package(BLAS REQUIRED)
+
 file(GLOB_RECURSE SOURCES src/*.cpp)
 file(GLOB_RECURSE HEADERS include/*.h)
 
 # Library target
 add_library(predict STATIC ${SOURCES})
 
-target_compile_options(
-  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>)
+target_compile_options(predict PRIVATE ${PREDICT_BUILD_ARGUMENTS})
+target_link_options(predict PRIVATE ${PREDICT_BUILD_ARGUMENTS})
 
 target_include_directories(
   predict PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
                  $<INSTALL_INTERFACE:include>)
 
-target_link_libraries(predict PUBLIC xtensor xsimd ${CASACORE_LIBRARIES})
+target_link_libraries(predict PUBLIC xtensor xtensor-blas xsimd
+                                     ${CASACORE_LIBRARIES} ${BLAS_LIBRARIES})
 
 target_include_directories(predict PRIVATE ${CASACORE_INCLUDE_DIRS})
 
@@ -93,12 +92,13 @@ if(PREDICT_BUILD_TESTING)
                                PRIVATE ${CASACORE_INCLUDE_DIRS})
 
     # Link Boost test framework to the test executable
-    target_link_libraries(${test_name}_unittests
-                          PRIVATE predict Boost::unit_test_framework)
+    target_link_libraries(
+      ${test_name}_unittests PRIVATE predict Boost::unit_test_framework xtensor
+                                     xtensor-blas)
     # Add a test that runs the test executable
 
-    target_compile_options(${test_name}_unittests PRIVATE -coverage
-                                                          $<$<CONFIG:Debug>:-g>)
+    target_compile_options(${test_name}_unittests
+                           PRIVATE -coverage ${PREDICT_BUILD_ARGUMENTS})
     target_link_options(${test_name}_unittests PRIVATE -coverage
                         $<$<CONFIG:Debug>:-g>)
 

benchmark/CMakeLists.txt

@@ -15,15 +15,13 @@ FetchContent_MakeAvailable(googlebenchmark)
 set(BENCHMARKS simulator)
 
 add_executable(microbenchmarks main.cpp predict.cpp directions.cpp spectrum.cpp
-                               spectrum_cross.cpp)
+                               phases.cpp spectrum_cross.cpp)
 
-target_compile_options(
-  microbenchmarks PRIVATE $<$<CONFIG:Debug>:-g> $<$<CONFIG:Release>:-O3>
-                          $<$<CONFIG:RelWithDebInfo>:-O3 -g>)
+target_compile_options(microbenchmarks PRIVATE ${PREDICT_BUILD_ARGUMENTS})
 
 target_link_libraries(
   microbenchmarks PRIVATE benchmark::benchmark predict # your main library
-                          xtensor xsimd ${CASACORE_LIBRARIES})
+                          xtensor xtensor-blas xsimd ${CASACORE_LIBRARIES})
 target_include_directories(
   microbenchmarks PRIVATE include ${CMAKE_SOURCE_DIR}/include
                           ${CASACORE_INCLUDE_DIRS})

benchmark/directions.cpp

-#include <Simulator.h>
+#include "Predict.h"
 #include <benchmark/benchmark.h>
 #include <random>
 #include <xtensor/xtensor.hpp>
@@ -24,7 +24,7 @@ public:
     n_directions = state.range(0);
     lmn = xt::xtensor<double, 2>({n_directions, 3});
     directions = std::make_unique<Directions>();
-    directions->reserve(n_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,
@@ -32,7 +32,7 @@ public:
       } else {
         SetUpFixedSource(test_reference_point, test_offset_point);
       }
-      directions->add(test_offset_point);
+      directions->Add(test_offset_point);
     }
   }
   void TearDown(benchmark::State &) override { directions.reset(); }
@@ -46,7 +46,7 @@ protected:
 BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsOneByOne)
 (benchmark::State &state) {
   for (auto _ : state) {
-    directions->radec2lmn<Directions::computation_strategy::SINGLE>(
+    directions->RaDec2Lmn<Directions::computation_strategy::SINGLE>(
         test_reference_point, lmn);
   }
 }
@@ -54,7 +54,7 @@ BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsOneByOne)
 BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsMulti)
 (benchmark::State &state) {
   for (auto _ : state) {
-    directions->radec2lmn<Directions::computation_strategy::MULTI>(
+    directions->RaDec2Lmn<Directions::computation_strategy::MULTI>(
         test_reference_point, lmn);
   }
 }
@@ -62,7 +62,7 @@ BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsMulti)
 BENCHMARK_DEFINE_F(DirectionsBenchmark, DirectionsMultiSIMD)
 (benchmark::State &state) {
   for (auto _ : state) {
-    directions->radec2lmn<Directions::computation_strategy::MULTI_SIMD>(
+    directions->RaDec2Lmn<Directions::computation_strategy::MULTI_SIMD>(
         test_reference_point, lmn);
   }
 }

benchmark/phases.cpp

+#include "PointSourceCollection.h"
+#include "Predict.h"
+#include "PredictPlan.h"
+#include "PredictPlanExec.h"
+#include <benchmark/benchmark.h>
+#include <random>
+#include <xtensor/xtensor.hpp>
+
+#include <Directions.h>
+#include <test/Common.h>
+
+const std::vector<int64_t> nstations = {24, 48};
+
+class PhasesBenchmark : public benchmark::Fixture {
+public:
+  void SetUp(benchmark::State &state) override {
+    n_stations = state.range(0);
+    n_channels = state.range(1);
+
+    PointSourceCollection sources;
+
+    PredictPlan plan{};
+    plan.nstations = n_stations;
+    plan.nchannels = n_channels;
+
+    plan_exec = std::make_unique<PredictPlanExec>(plan);
+
+    plan_exec->uvw = xt::xtensor<double, 2>({n_stations, 3});
+    for (size_t i = 0; i < n_stations; ++i) {
+      plan_exec->uvw(i, 0) = 100.0 + i;
+      plan_exec->uvw(i, 1) = 200.0 + i;
+      plan_exec->uvw(i, 2) = 300.0 + i;
+    }
+
+    plan_exec->frequencies.resize({n_channels});
+    for (size_t i = 0; i < n_channels; ++i) {
+      plan_exec->frequencies[i] = 1e8 + i * 1e6;
+    }
+
+    plan_exec->lmn = {{0.1, 0.2, 0.3}};
+  }
+
+  void TearDown(benchmark::State &) override { plan_exec.reset(); }
+
+protected:
+  std::unique_ptr<PredictPlanExec> plan_exec;
+  size_t n_stations;
+  size_t n_channels;
+};
+
+BENCHMARK_DEFINE_F(PhasesBenchmark, ComputePhases)(benchmark::State &state) {
+
+  for (auto _ : state) {
+    plan_exec->ComputeStationPhases();
+  }
+}
+
+BENCHMARK_REGISTER_F(PhasesBenchmark, ComputePhases)
+    ->ArgsProduct({nstations, {64, 128, 256}})
+    ->ArgNames({"#stations", "#channels"})
+    ->Unit(benchmark::kMicrosecond);

benchmark/predict.cpp

+#include "GaussianSourceCollection.h"
 #include <random>
 
 #include <GaussianSource.h>
@@ -6,7 +7,7 @@
 #include <test/Common.h>
 
 using dp3::base::GaussianSource;
-using dp3::base::InitializePredict;
+using dp3::base::MakePredictRun;
 using dp3::base::PointSource;
 using dp3::base::PredictRun;
 
@@ -35,7 +36,7 @@ public:
       SetUpFixedSource(test_reference_point, test_offset_point);
     }
 
-    predict_run = dp3::base::InitializePredict(
+    predict_run = dp3::base::MakePredictRun(
         test_reference_point, test_offset_point, stations, channels,
         !is_fullstokes, is_freqsmear);
   }
@@ -43,15 +44,16 @@ public:
   void TearDown(benchmark::State &) override { predict_run.reset(); }
 
 protected:
-  std::unique_ptr<const PredictRun> predict_run;
+  std::unique_ptr<PredictRun> predict_run;
 };
 
 BENCHMARK_DEFINE_F(PredictBenchmark, PointSource)
 (benchmark::State &state) {
-  auto source = predict_run->makeSource<PointSource>();
+  PointSourceCollection sources;
+  sources.Add(predict_run->makeSource<PointSource>());
 
   for (auto _ : state) {
-    predict_run->simulate(source);
+    predict_run->Run(sources);
   }
 }
 
@@ -67,10 +69,11 @@ BENCHMARK_REGISTER_F(PredictBenchmark, PointSource)
 
 BENCHMARK_DEFINE_F(PredictBenchmark, GaussianSource)
 (benchmark::State &state) {
-  auto source = predict_run->makeSource<GaussianSource>();
+  GaussianSourceCollection sources;
+  sources.Add(predict_run->makeSource<GaussianSource>());
 
   for (auto _ : state) {
-    predict_run->simulate(source);
+    predict_run->Run(sources);
   }
 }
 

benchmark/spectrum.cpp

-#include <Simulator.h>
+#include "Predict.h"
 #include <benchmark/benchmark.h>
 #include <random>
 #include <xtensor/xtensor.hpp>
@@ -26,26 +26,24 @@ public:
 
     spectrum = std::make_unique<Spectrum>();
     auto &spec = *spectrum;
-    spec.has_logarithmic_spectral_index = is_log;
-    spec.reference_flux = {1.0, 0.0, 0.0, 0.0};
-    spec.spectral_terms = {1.0e-6, 2.0e-6, 3.0e-6};
-    spec.reference_frequency = MHz;
-    spec.polarization_angle = 0.3;
-    spec.polarization_factor = 0.2;
-    spec.rotation_measure = 0.1;
-    spec.has_rotation_measure = compute_rotation;
+
+    spec.SetSpectralTerms(MHz, is_log, {1.0e-6, 2.0e-6, 3.0e-6});
+    spec.SetReferenceFlux({1.0, 0.0, 0.0, 0.0});
+    spec.SetRotationMeasure(0.1, compute_rotation);
+    spec.SetPolarization(0.3, 0.2);
+
     computed_spectrum = std::make_unique<StokesVector>();
 
     point_source = std::make_unique<dp3::base::PointSource>(
-        dp3::base::Direction(0.0, 0.0), spec.reference_flux);
+        dp3::base::Direction(0.0, 0.0), spec.GetReferenceFlux());
 
-    point_source->setSpectralTerms(
-        spec.reference_frequency, spec.has_logarithmic_spectral_index,
-        spec.spectral_terms.begin(), spec.spectral_terms.end());
+    point_source->SetSpectralTerms(
+        spec.GetReferenceFrequency(), spec.HasLogarithmicSpectralIndex(),
+        spec.GetSpectralTerms().begin(), spec.GetSpectralTerms().end());
     if (compute_rotation) {
-      point_source->setRotationMeasure(spec.polarization_factor,
-                                       spec.polarization_angle,
-                                       spec.rotation_measure);
+      point_source->SetRotationMeasure(spec.GetPolarizationFactor(),
+                                       spec.GetPolarizationAngle(),
+                                       spec.GetRotationMeasure());
     }
   }
   void TearDown(benchmark::State &) override { spectrum.reset(); }
@@ -69,7 +67,7 @@ BENCHMARK_DEFINE_F(SpectrumBenchmark, SpectrumReference)
 (benchmark::State &state) {
   for (auto _ : state) {
     for (size_t i = 0; i < n_frequencies; i++) {
-      benchmark::DoNotOptimize(point_source->stokes(frequencies(i)).I);
+      benchmark::DoNotOptimize(point_source->GetStokes(frequencies(i)).I);
     }
   }
 }

benchmark/spectrum_cross.cpp

-#include <Simulator.h>
 #include <benchmark/benchmark.h>
-#include <random>
 #include <xtensor/xtensor.hpp>
 
 #include <Direction.h>
@@ -25,14 +23,12 @@ public:
 
     spectrum = std::make_unique<Spectrum>();
     auto &spec = *spectrum;
-    spec.has_logarithmic_spectral_index = is_log;
-    spec.reference_flux = {1.0, 0.0, 0.0, 0.0};
-    spec.spectral_terms = {1.0e-6, 2.0e-6, 3.0e-6};
-    spec.reference_frequency = MHz;
-    spec.polarization_angle = 0.3;
-    spec.polarization_factor = 0.2;
-    spec.rotation_measure = 0.1;
-    spec.has_rotation_measure = compute_rotation;
+
+    spec.SetSpectralTerms(MHz, is_log, {1.0e-6, 2.0e-6, 3.0e-6});
+    spec.SetRotationMeasure(0.1);
+    spec.SetReferenceFlux({1.0, 0.0, 0.0, 0.0});
+    spec.SetPolarization(0.3, 0.2);
+
     computed_spectrum =
         std::make_unique<xt::xtensor<std::complex<double>, 2>>();
     computed_spectrum->resize({4, n_frequencies});

ci/das6/compile_and_run.sh

@@ -11,7 +11,7 @@ echo "Running on compiler version: ${COMPILER_VERSION}, architecture: ${ARCHITEC
 BUILD_DIR=build-${COMPILER_VERSION}-${ARCHITECTURE}
 module purge
 module load spack/20250403
-module load cmake casacore boost
+module load cmake casacore boost openblas
 
 cmake -B ${BUILD_DIR} . -DCMAKE_BUILD_TYPE=Release -DPREDICT_BUILD_BENCHMARK=ON
 make -C ${BUILD_DIR} -j

include/Directions.h

@@ -14,40 +14,44 @@
 #include <common/SmartVector.h>
 
 #include <vector>
+#include <xsimd/xsimd.hpp>
 #include <xtensor/xmath.hpp>
+#include <xtensor/xtensor.hpp>
 #include <xtensor/xview.hpp>
 
 using Direction = dp3::base::Direction;
 
+enum computation_strategy { SINGLE, MULTI, MULTI_SIMD };
+
 class Directions : ObjectCollection<Direction> {
 public:
   enum computation_strategy { SINGLE, MULTI, MULTI_SIMD };
   Directions() : ra(), dec() {}
-  void add(const Direction &direction) {
+  void Add(const Direction &direction) {
     ra.push_back(direction.ra);
     dec.push_back(direction.dec);
   }
 
-  void add(const std::span<Direction> &directions) {
-    ra.reserve(size() + directions.size());
-    dec.reserve(size() + directions.size());
+  void Add(const std::span<Direction> &directions) {
+    ra.reserve(Size() + directions.size());
+    dec.reserve(Size() + directions.size());
     for (size_t i = 0; i < directions.size(); ++i) {
       ra.push_back(directions[i].ra);
       dec.push_back(directions[i].dec);
     }
   }
 
-  void reserve(size_t new_size) {
+  void Reserve(size_t new_size) {
     ra.reserve(new_size);
     dec.reserve(new_size);
   }
 
-  void clear() {
+  void Clear() {
     ra.clear();
     dec.clear();
   }
 
-  size_t size() const { return ra.size(); }
+  size_t Size() const { return ra.size(); }
 
   SmartVector<double> ra;
   SmartVector<double> dec;
@@ -63,13 +67,13 @@ public:
    * coordinates will be written.
    */
   template <computation_strategy T = computation_strategy::MULTI_SIMD>
-  inline void radec2lmn(const Direction &reference,
-                        xt::xtensor<double, 2> &lmn);
+  inline void RaDec2Lmn(const Direction &reference,
+                        xt::xtensor<double, 2> &lmn) const;
 };
 
 template <>
-inline void Directions::radec2lmn<Directions::computation_strategy::MULTI>(
-    const Direction &reference, xt::xtensor<double, 2> &lmn) {
+inline void Directions::RaDec2Lmn<Directions::computation_strategy::MULTI>(
+    const Direction &reference, xt::xtensor<double, 2> &lmn) const {
   /**
    * \f{eqnarray*}{
    *   \ell &= \cos(\delta) \sin(\alpha - \alpha_0) \\
@@ -102,8 +106,8 @@ inline void Directions::radec2lmn<Directions::computation_strategy::MULTI>(
 }
 
 template <>
-inline void Directions::radec2lmn<Directions::computation_strategy::SINGLE>(
-    const Direction &reference, xt::xtensor<double, 2> &lmn) {
+inline void Directions::RaDec2Lmn<Directions::computation_strategy::SINGLE>(
+    const Direction &reference, xt::xtensor<double, 2> &lmn) const {
   /**
    * \f{eqnarray*}{
    *   \ell &= \cos(\delta) \sin(\alpha - \alpha_0) \\
@@ -131,9 +135,6 @@ inline void Directions::radec2lmn<Directions::computation_strategy::SINGLE>(
   }
 }
 
-#include <cstddef>
-#include <xsimd/xsimd.hpp>
-
 namespace xsimd {
 
 struct radec2lmn {
@@ -152,6 +153,7 @@ void radec2lmn::operator()(Arch, const Direction &reference, const C &ra,
   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());
@@ -193,9 +195,10 @@ void radec2lmn::operator()(Arch, const Direction &reference, const C &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()});
+inline void Directions::RaDec2Lmn<Directions::computation_strategy::MULTI_SIMD>(
+    const Direction &reference, xt::xtensor<double, 2> &lmn) const {
+  xt::xtensor<double, 2> lmn_tmp;
+  lmn_tmp.resize({3, ra.size()});
   xsimd::dispatch(xsimd::radec2lmn{})(reference, ra, dec, lmn_tmp,
                                       xsimd::unaligned_mode());
   lmn = xt::transpose(lmn_tmp);

include/GaussianSource.h

@@ -25,37 +25,37 @@ public:
 
   /// Set position angle in radians. The position angle is the smallest angle
   /// between the major axis and North, measured positively North over East.
-  void setPositionAngle(double angle);
-  double getPositionAngle() const { return itsPositionAngle; }
+  void SetPositionAngle(double angle);
+  double GetPositionAngle() const { return position_angle_; }
 
   /// Set whether the position angle (orientation) is absolute, see
   /// documentation of class member)
-  void setPositionAngleIsAbsolute(bool positionAngleIsAbsolute) {
-    itsPositionAngleIsAbsolute = positionAngleIsAbsolute;
+  void SetPositionAngleIsAbsolute(bool positionAngleIsAbsolute) {
+    is_position_angle_absolute_ = positionAngleIsAbsolute;
   }
 
   /// Return whether the position angle (orientation) is absolute, see
   /// documentation of class member.
-  bool getPositionAngleIsAbsolute() const { return itsPositionAngleIsAbsolute; }
+  bool GetPositionAngleIsAbsolute() const {
+    return is_position_angle_absolute_;
+  }
 
   /// Set the major axis length (FWHM in radians).
-  void setMajorAxis(double fwhm);
-  double getMajorAxis() const { return itsMajorAxis; }
+  void SetMajorAxis(double fwhm);
+  double GetMajorAxis() const { return major_axis_; }
 
   /// Set the minor axis length (FWHM in radians).
-  void setMinorAxis(double fwhm);
-  double getMinorAxis() const { return itsMinorAxis; }
-
-  void accept(ModelComponentVisitor &visitor) const override;
+  void SetMinorAxis(double fwhm);
+  double GetMinorAxis() const { return minor_axis_; }
 
 private:
-  double itsPositionAngle;
+  double position_angle_;
   /// Whether the position angle (also refered to as orientation) is absolute
   /// (w.r.t. to the local declination axis) or with respect to the declination
   /// axis at the phase center (the default until 2022, it was fixed in 5.3.0)
-  bool itsPositionAngleIsAbsolute;
-  double itsMajorAxis;
-  double itsMinorAxis;
+  bool is_position_angle_absolute_;
+  double major_axis_;
+  double minor_axis_;
 };
 
 /// @}

include/GaussianSourceCollection.h

@@ -9,6 +9,8 @@
 #ifndef GAUSSIAN_SOURCE_COLLECTION_H_
 #define GAUSSIAN_SOURCE_COLLECTION_H_
 
+#include <vector>
+
 #include <GaussianSource.h>
 #include <ObjectCollection.h>
 #include <PointSourceCollection.h>
@@ -18,25 +20,21 @@ using GaussianSource = dp3::base::GaussianSource;
 
 class GaussianSourceCollection : public PointSourceCollection {
 public:
-  void add(const GaussianSource &gaussian_source) {
-    PointSourceCollection::add(gaussian_source.direction(),
-                               gaussian_source.stokes());
+  void Add(const GaussianSource &gaussian_source) {
+    PointSourceCollection::Add(gaussian_source);
 
-    position_angle.push_back(gaussian_source.getPositionAngle());
-    major_axis.push_back(gaussian_source.getMajorAxis());
-    minor_axis.push_back(gaussian_source.getMinorAxis());
+    position_angle.push_back(gaussian_source.GetPositionAngle());
+    major_axis.push_back(gaussian_source.GetMajorAxis());
+    minor_axis.push_back(gaussian_source.GetMinorAxis());
     position_angle_is_absolute.push_back(
-        gaussian_source.getPositionAngleIsAbsolute());
+        gaussian_source.GetPositionAngleIsAbsolute());
   }
 
-  size_t size() const { return position_angle.size(); }
-
-  void reserve(size_t new_size) {
+  void Reserve(size_t new_size) {
     position_angle.reserve(new_size);
     major_axis.reserve(new_size);
     minor_axis.reserve(new_size);
     position_angle_is_absolute.reserve(new_size);
-    PointSourceCollection::reserve(new_size);
   }
 
   SmartVector<double> position_angle;

include/ModelComponent.h

-// ModelComponent.h: Base class for model components.
-//
-// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
-// SPDX-License-Identifier: GPL-3.0-or-later
-
-#ifndef DPPP_MODELCOMPONENT_H
-#define DPPP_MODELCOMPONENT_H
-
-#include <memory>
-
-namespace dp3 {
-namespace base {
-
-class ModelComponentVisitor;
-struct Direction;
-
-/// \brief Base class for model components.
-
-/// @{
-
-class ModelComponent {
-public:
-  virtual ~ModelComponent() {}
-  virtual const Direction &direction() const = 0;
-  virtual void accept(ModelComponentVisitor &) const = 0;
-};
-
-/// @}
-
-} // namespace base
-} // namespace dp3
-
-#endif

include/ModelComponentVisitor.h

-// ModelComponentVisitor.h: Base class for visitors that visit model component
-// hierarchies.
-//
-// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
-// SPDX-License-Identifier: GPL-3.0-or-later
-
-#ifndef DPPP_MODELCOMPONENTVISITOR_H
-#define DPPP_MODELCOMPONENTVISITOR_H
-
-namespace dp3 {
-namespace base {
-
-class PointSource;
-class GaussianSource;
-
-/// \brief Base class for visitors that visit model component hierarchies.
-
-/// @{
-
-class ModelComponentVisitor {
-public:
-  virtual ~ModelComponentVisitor() {}
-
-  virtual void visit(const PointSource &) = 0;
-  virtual void visit(const GaussianSource &) = 0;
-};
-
-/// @}
-
-} // namespace base
-} // namespace dp3
-
-#endif

include/ObjectCollection.h

@@ -13,11 +13,11 @@
 
 template <class T> class ObjectCollection {
 public:
-  virtual void add(const T &object) = 0;
-  virtual void add(const std::span<T> &object) = 0;
-  virtual void clear() = 0;
-  virtual size_t size() const = 0;
-  virtual void reserve(size_t size) = 0;
+  virtual void Add(const T &object) = 0;
+  virtual void Add(const std::span<T> &object) = 0;
+  virtual void Clear() = 0;
+  virtual size_t Size() const = 0;
+  virtual void Reserve(size_t size) = 0;
 };
 
 #endif // OBJECT_COLLECTION_H_

include/PointSource.h

@@ -7,13 +7,10 @@
 #ifndef DPPP_POINTSOURCE_H
 #define DPPP_POINTSOURCE_H
 
-#include <vector>
-
-#include "ModelComponent.h"
+#include "Direction.h"
+#include "Spectrum.h"
 #include "Stokes.h"
-#include <Direction.h>
-
-#include <iostream>
+#include "StokesVector.h"
 #include <memory>
 namespace dp3 {
 namespace base {
@@ -23,59 +20,68 @@ namespace base {
 
 /// @{
 
-class PointSource : public ModelComponent {
+class PointSource {
 public:
   typedef std::shared_ptr<PointSource> Ptr;
   typedef std::shared_ptr<const PointSource> ConstPtr;
 
-  PointSource(const Direction &position);
-  PointSource(const Direction &position, const Stokes &stokes);
+  PointSource(const Direction &direction, const Stokes &stokes);
 
-  const Direction &direction() const override { return itsDirection; }
-  void setDirection(const Direction &position);
+  const Direction &GetDirection() const { return direction_; }
+  void SetDirection(const Direction &position);
 
-  void setStokes(const Stokes &stokes);
+  void ComputeSpectrum(const xt::xtensor<double, 1> &frequencies,
+                       xt::xtensor<std::complex<double>, 2> &result) const;
+  const Spectrum &GetSpectrum() const { return spectrum_; }
+
+  Stokes GetStokes(double freq) const;
+  const Stokes &GetStokes() const { return stokes_; }
+  const StokesVector &GetStokesVector() const { return stokes_vector_; }
 
   template <typename T>
-  void setSpectralTerms(double refFreq, bool isLogarithmic, T first, T last);
+  void SetSpectralTerms(double refFreq, bool isLogarithmic, T first, T last);
 
-  void setRotationMeasure(double fraction, double angle, double rm);
+  void SetRotationMeasure(double fraction, double angle, double rm);
 
-  Stokes stokes(double freq) const;
-  Stokes stokes() const;
-  const std::vector<double> &spectrum() const { return itsSpectralTerms; }
-  double referenceFreq() const { return itsRefFreq; }
-  double spectralTerm(size_t i) const { return itsSpectralTerms[i]; }
-  void accept(ModelComponentVisitor &visitor) const override;
+  bool HasRotationMeasure() const;
+  bool HasSpectralTerms() const;
 
 private:
-  bool hasSpectralTerms() const;
-  bool hasRotationMeasure() const;
-
-  Direction itsDirection;
-  Stokes itsStokes;
-  double itsRefFreq;
-  std::vector<double> itsSpectralTerms;
-  double itsPolarizedFraction;
-  double itsPolarizationAngle;
-  double itsRotationMeasure;
-  bool itsHasRotationMeasure;
-  bool itsHasLogarithmicSI;
+  Direction direction_;
+  Spectrum spectrum_;
+  StokesVector stokes_vector_;
+  Stokes stokes_;
+
+  // FIXME: remove the following declerations.
+  // They are not needed anymore, but are kept for
+  // backward compatibility for the tests.
+  double reference_frequency_;
+  std::vector<double> spectral_terms_;
+  double polarizated_fraction_;
+  double polarization_angle_;
+  double rotation_measure_;
+  bool has_rotation_measure_;
+  bool has_logarithmic_si_;
 };
 
 /// @}
 
-// -------------------------------------------------------------------------- //
-// - Implementation: PointSource                                            - //
-// -------------------------------------------------------------------------- //
-
 template <typename T>
-void PointSource::setSpectralTerms(double refFreq, bool isLogarithmic, T first,
+void PointSource::SetSpectralTerms(double refFreq, bool isLogarithmic, T first,
                                    T last) {
-  itsRefFreq = refFreq;
-  itsHasLogarithmicSI = isLogarithmic;
-  itsSpectralTerms.clear();
-  itsSpectralTerms.insert(itsSpectralTerms.begin(), first, last);
+  // FIXME: the following four statements can be removed later.
+  // They are not needed anymore, but are kept for
+  // backward compatibility for the tests.
+  reference_frequency_ = refFreq;
+  has_logarithmic_si_ = isLogarithmic;
+
+  spectral_terms_.clear();
+  spectral_terms_.insert(spectral_terms_.begin(), first, last);
+  auto new_terms = xt::adapt(std::vector<double>(first, last));
+
+  spectrum_.SetSpectralTerms(
+      refFreq, isLogarithmic,
+      xt::concatenate(xt::xtuple(spectrum_.GetSpectralTerms(), new_terms), 0));
 }
 
 } // namespace base

include/PointSourceCollection.h

@@ -9,6 +9,8 @@
 #ifndef POINT_SOURCE_COLLECTION_H_
 #define POINT_SOURCE_COLLECTION_H_
 
+#include <vector>
+
 #include <Directions.h>
 #include <ObjectCollection.h>
 #include <PointSource.h>
@@ -19,45 +21,55 @@ using PointSource = dp3::base::PointSource;
 class PointSourceCollection : public ObjectCollection<PointSource> {
 public:
   Directions direction_vector;
+  std::vector<StokesVector> spectrum_vector;
+  std::vector<Spectrum> spectrums;
+
   StokesVector stokes_vector;
 
-  void add(const PointSource &point_source) {
-    direction_vector.add(point_source.direction());
-    stokes_vector.add(point_source.stokes());
+  void Add(const PointSource &point_source) {
+    direction_vector.Add(point_source.GetDirection());
+    stokes_vector.Add(point_source.GetStokes());
+    spectrums.push_back(point_source.GetSpectrum());
   }
 
-  void add(const Direction &direction, const Stokes &stokes) {
-    direction_vector.add(direction);
-    stokes_vector.add(stokes);
+  void Add(const Direction &direction, const Stokes &stokes) {
+    direction_vector.Add(direction);
+    stokes_vector.Add(stokes);
+    spectrums.push_back(Spectrum());
   }
 
-  void add(const std::span<Direction> &directions,
+  void Add(const std::span<Direction> &directions,
            const std::span<Stokes> &stokes) {
-    direction_vector.add(directions);
-    stokes_vector.add(stokes);
+    direction_vector.Add(directions);
+    stokes_vector.Add(stokes);
+    spectrums.push_back(Spectrum());
   }
 
-  void reserve(size_t new_size) {
-    direction_vector.reserve(new_size);
-    stokes_vector.reserve(new_size);
+  void Reserve(size_t new_size) {
+    direction_vector.Reserve(new_size);
+    stokes_vector.Reserve(new_size);
+    spectrums.reserve(new_size);
   }
 
-  void add(const std::span<PointSource> &point_sources) {
-    const size_t expectedSize = size() + point_sources.size();
-    direction_vector.reserve(expectedSize);
-    stokes_vector.reserve(expectedSize);
+  void Add(const std::span<PointSource> &point_sources) {
+    const size_t expectedSize = Size() + point_sources.size();
+    direction_vector.Reserve(expectedSize);
     std::for_each(point_sources.begin(), point_sources.end(),
                   [this](const PointSource &point_source) {
-                    direction_vector.add(point_source.direction());
-                    stokes_vector.add(point_source.stokes());
+                    direction_vector.Add(point_source.GetDirection());
+                    stokes_vector.Add(point_source.GetStokes());
+                    spectrums.push_back(point_source.GetSpectrum());
                   });
   }
 
-  void clear() {
-    direction_vector.clear();
-    stokes_vector.clear();
+  void Clear() {
+    PointSourceCollection::Clear();
+    direction_vector.Clear();
+    spectrum_vector.clear();
+    spectrums.clear();
   }
-  size_t size() const { return direction_vector.size(); }
+
+  size_t Size() const { return direction_vector.Size(); }
 };
 
 #endif // POINT_SOURCE_COLLECTION_H_

include/Simulator.h → include/Predict.h

-// Simulator.h: Compute visibilities for different model components types
+// Predict.h: Compute visibilities for different model components types
 // (implementation of ModelComponentVisitor).
 //
 // Copyright (C) 2023 ASTRON (Netherlands Institute for Radio Astronomy)
 // SPDX-License-Identifier: GPL-3.0-or-later
 
-#ifndef DP3_BASE_SIMULATOR_H_
-#define DP3_BASE_SIMULATOR_H_
-
-#include <vector>
+#ifndef PREDICT_H_
+#define PREDICT_H_
 
 #include <casacore/casa/Arrays/Cube.h>
+#include <casacore/measures/Measures.h>
 
 #include <xtensor/xtensor.hpp>
 
-#include "Baseline.h"
-#include "ModelComponent.h"
-#include "ModelComponentVisitor.h"
+#include "GaussianSourceCollection.h"
+#include "PredictPlanExec.h"
 #include <Direction.h>
+#include <PointSource.h>
+#include <Stokes.h>
 
 namespace dp3 {
 namespace base {
@@ -58,9 +58,32 @@ inline void radec2lmn(const Direction &reference, const Direction &direction,
   lmn[2] = sin_dec * sin_dec0 + cos_dec * cos_dec0 * cos_delta_ra;
 }
 
-/// @{
+// Compute component spectrum.
+inline void spectrum(const PointSource &component, size_t nChannel,
+                     const xt::xtensor<double, 1> &freq,
+                     xt::xtensor<std::complex<double>, 2> &spectrum,
+                     bool stokesIOnly = false) {
+#pragma GCC ivdep
+  for (size_t ch = 0; ch < nChannel; ++ch) {
+    Stokes stokes = component.GetStokes(freq(ch));
+
+    if (stokesIOnly) {
+      spectrum(0, ch).real(stokes.I);
+      spectrum(0, ch).imag(0.0);
+    } else {
+      spectrum(0, ch).real(stokes.I + stokes.Q);
+      spectrum(0, ch).imag(0.0);
+      spectrum(1, ch).real(stokes.U);
+      spectrum(1, ch).imag(stokes.V);
+      spectrum(2, ch).real(stokes.U);
+      spectrum(2, ch).imag(-stokes.V);
+      spectrum(3, ch).real(stokes.I - stokes.Q);
+      spectrum(3, ch).imag(0.0);
+    }
+  }
+}
 
-typedef std::complex<double> dcomplex;
+/// @{
 
 /**
  * @brief Simulator to compute visibilities given a sky model
@@ -97,88 +120,12 @@ typedef std::complex<double> dcomplex;
  * to the position of the source.
  */
 
-class Simulator : public ModelComponentVisitor {
+class Predict {
 public:
-  /**
-   * @brief Construct a new Simulator object
-   *
-   * @param reference Reference direction (phase center)
-   * @param baselines Vector of Baselines
-   * @param freq Channel frequencies (Hz)
-   * @param chanWidths Channel widths (Hz)
-   * @param stationUVW Station UVW coordinates. This structure has to remain
-   * valid during the lifetime of the Simulator.
-   * @param buffer Output buffer, should be of shape (nCor, nFreq, nBaselines),
-   * where nCor should be 1 if stokesIOnly is true, else 4
-   * @param correctFreqSmearing Correct for frequency smearing
-   * @param stokesIOnly Stokes I only, to avoid a loop over correlations
-   */
-  Simulator(const Direction &reference, size_t nStation,
-            const std::vector<Baseline> &baselines,
-            const std::vector<double> &freq,
-            const std::vector<double> &chanWidths,
-            const xt::xtensor<double, 2> &stationUVW,
-            xt::xtensor<dcomplex, 3> &buffer, bool correctFreqSmearing,
-            bool stokesIOnly);
-
-  // Note DuoMatrix is actually two T matrices
-  // T: floating point type, ideally float, double, or long double.
-  template <typename T> class DuoMatrix {
-  public:
-    DuoMatrix() = default;
-
-    DuoMatrix(size_t nrows, size_t ncols) { resize(nrows, ncols); }
-
-    void resize(size_t nrows, size_t ncols) {
-      itsNRows = nrows;
-      itsData_real.resize(nrows * ncols);
-      itsData_imag.resize(nrows * ncols);
-    }
-
-    size_t nRows() { return itsNRows; }
-    size_t nCols() { return itsData_real.size() / itsNRows; }
-
-    T &real(size_t row, size_t col) {
-      return itsData_real[col * itsNRows + row];
-    }
-    T &imag(size_t row, size_t col) {
-      return itsData_imag[col * itsNRows + row];
-    }
-    T *realdata() { return &itsData_real[0]; }
-    T *imagdata() { return &itsData_imag[0]; }
-
-  private:
-    // Use separate storage for real/imag parts
-    std::vector<T> itsData_real;
-    std::vector<T> itsData_imag;
-    size_t itsNRows{0};
-  };
-
-  enum StokesComponent { I = 0, Q = 1, U = 2, V = 3 };
-
-  void simulate(const std::shared_ptr<const ModelComponent> &component);
-
-private:
-  void visit(const PointSource &component) override;
-  void visit(const GaussianSource &component) override;
-
-private:
-  Direction itsReference;
-  size_t itsNStation, itsNBaseline, itsNChannel;
-  bool itsCorrectFreqSmearing;
-  bool itsStokesIOnly;
-  std::vector<Baseline> itsBaselines;
-  std::vector<double> itsFreq;
-  std::vector<double> itsChanWidths;
-  /// Non-owning pointer to UVW values for each station. The user of Simulator
-  /// supplies them in the constructor, and ensures they remain valid.
-  /// Using a pointer avoids copying the values.
-  const xt::xtensor<double, 2> *itsStationUVW;
-  xt::xtensor<dcomplex, 3> *itsBuffer;
-
-  std::vector<double> itsStationPhases;
-  DuoMatrix<double> itsShiftBuffer;
-  DuoMatrix<double> itsSpectrumBuffer;
+  void run(PredictPlanExec &plan, const PointSourceCollection &sources,
+           xt::xtensor<dcomplex, 3> &buffer) const;
+  void run(PredictPlanExec &plan, const GaussianSourceCollection &sources,
+           xt::xtensor<dcomplex, 3> &buffer) const;
 };
 
 /// @}

include/PredictPlan.h

+// PredictPlan.h: PredictPlan
+//
+// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: Apache2.0
+
+/// \file
+/// \brief PredictPlan
+
+#ifndef PREDICT_PLAN_H_
+#define PREDICT_PLAN_H_
+
+#include <Baseline.h>
+#include <Direction.h>
+#include <xtensor/xlayout.hpp>
+#include <xtensor/xtensor.hpp>
+
+#include "GaussianSourceCollection.h"
+#include "PointSourceCollection.h"
+
+using Direction = dp3::base::Direction;
+using Baseline = dp3::base::Baseline;
+
+typedef std::complex<double> dcomplex;
+
+struct PredictPlan {
+  enum StokesComponent : uint_fast8_t { I = 0, Q = 1, U = 2, V = 3 };
+
+  Direction reference;
+  bool correct_frequency_smearing;
+  bool compute_stokes_I_only;
+  std::vector<Baseline> baselines;
+  std::vector<double> channel_widths;
+  size_t nchannels;
+  xt::xtensor<double, 2, xt::layout_type::column_major> uvw;
+  xt::xtensor<double, 1> frequencies;
+  size_t nstations;
+  size_t nbaselines;
+  size_t nstokes;
+
+  // FIXME: We probably want to implement this as a pure virtual function, i.e.,
+  // precompute() = 0, although this would not allow us to construct a
+  // standalone PredictPlan object
+  virtual void Precompute(const PointSourceCollection &) {}
+
+  virtual void Compute(const PointSourceCollection &sources,
+                       xt::xtensor<std::complex<double>, 3> &buffer) {}
+
+  virtual void Compute(const GaussianSourceCollection &sources,
+                       xt::xtensor<std::complex<double>, 3> &buffer) {}
+};
+
+#endif // PREDICT_PLAN_H_

include/PredictPlanExec.h

+// PredictPlanExec.h: PredictPlanExec
+//
+// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: Apache2.0
+
+/// \file
+/// \brief PredictPlanExec
+
+#ifndef PREDICT_PLAN_EXEC_H_
+#define PREDICT_PLAN_EXEC_H_
+
+#include "GaussianSourceCollection.h"
+#include "PointSourceCollection.h"
+#include "PredictPlan.h"
+#include <casacore/measures/Measures.h>
+#include <complex>
+#include <xtensor/xtensor.hpp>
+class PredictPlanExec : public PredictPlan {
+public:
+  explicit PredictPlanExec(const PredictPlan &plan) : PredictPlan(plan) {
+    // Initialize the station phases and shifts
+    Initialize();
+  }
+
+  virtual void Precompute(const PointSourceCollection &sources) final override;
+
+  virtual void
+  Compute(const PointSourceCollection &sources,
+          xt::xtensor<std::complex<double>, 3> &buffer) final override;
+
+  virtual void
+  Compute(const GaussianSourceCollection &sources,
+          xt::xtensor<std::complex<double>, 3> &buffer) final override;
+
+  void ComputeStationPhases();
+
+  const xt::xtensor<double, 1> &GetStationPhases() const {
+    return station_phases;
+  }
+  const xt::xtensor<std::complex<double>, 2> &GetShiftData() const {
+    return shift_data;
+  }
+
+  void FillEulerMatrix(xt::xtensor<double, 2> &mat, const double ra,
+                       const double dec);
+
+  // FIXME: move this to a more appropriate location
+  inline float ComputeSmearterm(double uvw, double halfwidth) {
+    float smearterm = static_cast<float>(uvw) * static_cast<float>(halfwidth);
+    return (smearterm == 0.0f) ? 1.0f
+                               : std::fabs(std::sin(smearterm) / smearterm);
+  }
+
+  const xt::xtensor<double, 2> &get_lmn() const { return lmn; }
+  const xt::xtensor<double, 2, xt::layout_type::column_major> &get_uvw() const {
+    return uvw;
+  }
+  const xt::xtensor<double, 1> &get_frequencies() const { return frequencies; }
+
+  xt::xtensor<double, 2> lmn;
+  xt::xtensor<double, 1> station_phases;
+  xt::xtensor<std::complex<double>, 2> shift_data;
+
+private:
+  void Initialize();
+
+private:
+  const double kCInv = 2.0 * M_PI / casacore::C::c;
+};
+
+#endif // PREDICT_PLAN_EXEC_H_