Merge branch 'implement_spectrum_computation' into 'main'

Add implementation spectrum and tests

See merge request mancini/predict!26

.gitlab-ci.yml

@@ -150,6 +150,8 @@ 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 SpectrumBenchmark/Spectrum results*.json result-summary-spectrum
+  
   
   
 

benchmark/CMakeLists.txt

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

benchmark/spectrum.cpp

+#include <Simulator.h>
+#include <benchmark/benchmark.h>
+#include <random>
+#include <xtensor/xtensor.hpp>
+
+#include <Direction.h>
+#include <PointSource.h>
+#include <Spectrum.h>
+
+const std::vector<int64_t> n_frequencies = {128, 256, 512};
+
+const double kHz = 1.e3;
+const double MHz = 1.e6;
+namespace {} // namespace
+
+class SpectrumBenchmark : public benchmark::Fixture {
+public:
+  void SetUp(benchmark::State &state) override {
+    n_frequencies = state.range(0);
+    const bool is_log = state.range(1);
+    const bool compute_rotation = state.range(2);
+    frequencies.resize({n_frequencies});
+    for (int i = 0; i < n_frequencies; i++) {
+      frequencies(i) = i * 100. * kHz + MHz;
+    }
+
+    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;
+    computed_spectrum = std::make_unique<StokesVector>();
+
+    point_source = std::make_unique<dp3::base::PointSource>(
+        dp3::base::Direction(0.0, 0.0), spec.reference_flux);
+
+    point_source->setSpectralTerms(
+        spec.reference_frequency, spec.has_logarithmic_spectral_index,
+        spec.spectral_terms.begin(), spec.spectral_terms.end());
+    if (compute_rotation) {
+      point_source->setRotationMeasure(spec.polarization_factor,
+                                       spec.polarization_angle,
+                                       spec.rotation_measure);
+    }
+  }
+  void TearDown(benchmark::State &) override { spectrum.reset(); }
+
+protected:
+  std::unique_ptr<Spectrum> spectrum;
+  std::unique_ptr<StokesVector> computed_spectrum;
+  size_t n_frequencies;
+  xt::xtensor<double, 1> frequencies;
+  std::unique_ptr<dp3::base::PointSource> point_source;
+};
+
+BENCHMARK_DEFINE_F(SpectrumBenchmark, SpectrumBenchmark)
+(benchmark::State &state) {
+  for (auto _ : state) {
+    spectrum->evaluate(frequencies, *computed_spectrum);
+  }
+}
+
+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_REGISTER_F(SpectrumBenchmark, SpectrumBenchmark)
+    ->ArgsProduct({
+        n_frequencies,
+        {false, true},
+        {false, true},
+
+    })
+    ->ArgNames({"#frequencies", "is_log", "compute_rotation"})
+    ->Unit(benchmark::kNanosecond);
+
+BENCHMARK_REGISTER_F(SpectrumBenchmark, SpectrumReference)
+    ->ArgsProduct({
+        n_frequencies,
+        {false, true},
+        {false, true},
+
+    })
+    ->ArgNames({"#frequencies", "is_log", "compute_rotation"})
+    ->Unit(benchmark::kNanosecond);

include/PointSource.h

@@ -13,8 +13,8 @@
 #include "Stokes.h"
 #include <Direction.h>
 
+#include <iostream>
 #include <memory>
-
 namespace dp3 {
 namespace base {
 
@@ -45,7 +45,7 @@ public:
   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;
 
 private:

include/Spectrum.h

@@ -12,11 +12,13 @@
 #include <Stokes.h>
 #include <StokesVector.h>
 
+#include <casacore/casa/BasicSL/Constants.h>
 #include <xtensor/xtensor.hpp>
 
 using Stokes = dp3::base::Stokes;
 
-struct Spectrum {
+class Spectrum {
+public:
   Stokes reference_flux;
   xt::xtensor<double, 1> spectral_terms;
   double reference_frequency;
@@ -27,8 +29,76 @@ struct Spectrum {
   bool has_rotation_measure;
   bool has_logarithmic_spectral_index;
 
-  StokesVector evaluate(const xt::xtensor<double, 1> &frequencies) const {
-    return StokesVector();
+  void evaluate(const xt::xtensor<double, 1> &frequencies,
+                StokesVector &result) const {
+    result.I.resize(frequencies.size());
+    result.Q.resize(frequencies.size());
+    result.U.resize(frequencies.size());
+    result.V.resize(frequencies.size());
+
+    double I0 = reference_flux.I;
+
+    if (spectral_terms.size() > 0 && has_logarithmic_spectral_index) {
+      for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
+        result.I[f_idx] = evaluate_log_spectrum(
+            I0, frequencies[f_idx], reference_frequency, spectral_terms);
+      }
+    } else if (spectral_terms.size()) {
+      for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
+        result.I[f_idx] = evaluate_spectrum(
+            I0, frequencies[f_idx], reference_frequency, spectral_terms);
+      }
+    } else {
+      for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
+        result.I[f_idx] = I0;
+      }
+    }
+
+    if (has_rotation_measure) {
+      // This way of computing is faster then using the xtensor extension
+      for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
+        const double lambda = casacore::C::c / frequencies[f_idx];
+        const double chi =
+            2.0 * (polarization_angle + (rotation_measure * lambda * lambda));
+        const double stokesQU = result.I[f_idx] * polarization_factor;
+        result.Q[f_idx] = stokesQU * cos(chi);
+        result.U[f_idx] = stokesQU * sin(chi);
+      }
+    }
+  }
+
+private:
+  inline double
+  evaluate_polynomial(double x,
+                      const xt::xtensor<double, 1> &parameters) const {
+    double partial = 0.0;
+    for (auto it = spectral_terms.rbegin(), end = spectral_terms.rend();
+         it != end; ++it) {
+      partial = std::fma(partial, x, *it);
+    }
+    return partial;
+  }
+  inline double
+  evaluate_log_polynomial(double x,
+                          const xt::xtensor<double, 1> &parameters) const {
+    double base = log10(x);
+    return evaluate_polynomial(base, parameters);
+  }
+
+  inline double
+  evaluate_log_spectrum(double I0, const double frequency,
+                        const double reference_frequency,
+                        const xt::xtensor<double, 1> &parameters) const {
+    double x = frequency / reference_frequency;
+    return I0 * pow(10, x * evaluate_log_polynomial(x, parameters));
+  }
+  inline double
+  evaluate_spectrum(double I0, const double frequency,
+                    const double reference_frequency,
+                    const xt::xtensor<double, 1> &parameters) const {
+    double x = frequency / reference_frequency - 1.0;
+
+    return I0 + evaluate_polynomial(x, parameters) * x;
   }
 };
 

tests/test_spectrum.cpp

@@ -3,11 +3,175 @@
 // SPDX-License-Identifier: Apache2.0
 
 #define BOOST_TEST_MODULE SPECTRUM_TEST
+#include <Direction.h>
+#include <PointSource.h>
 #include <Spectrum.h>
 #include <boost/test/unit_test.hpp>
 
 BOOST_AUTO_TEST_SUITE(SpectrumTestSuite)
 
-BOOST_AUTO_TEST_CASE(test_Spectrum) { BOOST_CHECK(true); }
+struct SpectrumFixture {
+  SpectrumFixture() {
+    spectrum.reference_flux = {1.0, 0.0, 0.0, 0.0};
+    spectrum.spectral_terms = {1.0, 2.0, 3.0};
+    spectrum.reference_frequency = 1.0e6;
+    spectrum.polarization_angle = 0.3;
+    spectrum.polarization_factor = 0.2;
+    spectrum.rotation_measure = 0.1;
+    spectrum.has_rotation_measure = false;
+    spectrum.has_logarithmic_spectral_index = false;
+  }
+  ~SpectrumFixture() { BOOST_TEST_MESSAGE("teardown fixture"); }
+
+  Spectrum spectrum;
+  xt::xtensor<double, 1> frequencies = {15.0e6, 20.0e6, 30.0e6};
+};
+
+BOOST_FIXTURE_TEST_CASE(normal_no_rotation, SpectrumFixture) {
+  dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0),
+                                      spectrum.reference_flux);
+
+  point_source.setSpectralTerms(
+      spectrum.reference_frequency, spectrum.has_logarithmic_spectral_index,
+      spectrum.spectral_terms.begin(), spectrum.spectral_terms.end());
+  point_source.setRotationMeasure(0.0, 0.0, 0.0);
+
+  Stokes stokes_nu0(0, 0, 0, 0);
+  Stokes stokes_nu1(0, 0, 0, 0);
+  Stokes stokes_nu2(0, 0, 0, 0);
+
+  stokes_nu0 = point_source.stokes(frequencies[0]);
+  stokes_nu1 = point_source.stokes(frequencies[1]);
+  stokes_nu2 = point_source.stokes(frequencies[2]);
+
+  StokesVector result;
+  spectrum.evaluate(frequencies, result);
+
+  BOOST_CHECK_CLOSE(result.I[0], stokes_nu0.I, 1.e-6);
+  BOOST_CHECK_CLOSE(result.I[1], stokes_nu1.I, 1.e-6);
+  BOOST_CHECK_CLOSE(result.I[2], stokes_nu2.I, 1.e-6);
+  BOOST_CHECK_CLOSE(result.Q[0], stokes_nu0.Q, 1.e-6);
+  BOOST_CHECK_CLOSE(result.Q[1], stokes_nu1.Q, 1.e-6);
+  BOOST_CHECK_CLOSE(result.Q[2], stokes_nu2.Q, 1.e-6);
+  BOOST_CHECK_CLOSE(result.U[0], stokes_nu0.U, 1.e-6);
+  BOOST_CHECK_CLOSE(result.U[1], stokes_nu1.U, 1.e-6);
+  BOOST_CHECK_CLOSE(result.U[2], stokes_nu2.U, 1.e-6);
+  BOOST_CHECK_CLOSE(result.V[0], stokes_nu0.V, 1.e-6);
+  BOOST_CHECK_CLOSE(result.V[1], stokes_nu1.V, 1.e-6);
+  BOOST_CHECK_CLOSE(result.V[2], stokes_nu2.V, 1.e-6);
+}
+
+BOOST_FIXTURE_TEST_CASE(normal_with_rotation, SpectrumFixture) {
+  spectrum.has_rotation_measure = true;
+
+  dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0),
+                                      spectrum.reference_flux);
+
+  point_source.setSpectralTerms(
+      spectrum.reference_frequency, spectrum.has_logarithmic_spectral_index,
+      spectrum.spectral_terms.begin(), spectrum.spectral_terms.end());
+  point_source.setRotationMeasure(spectrum.polarization_factor,
+                                  spectrum.polarization_angle,
+                                  spectrum.rotation_measure);
+  Stokes stokes_nu0(0, 0, 0, 0);
+  Stokes stokes_nu1(0, 0, 0, 0);
+  Stokes stokes_nu2(0, 0, 0, 0);
+
+  stokes_nu0 = point_source.stokes(frequencies[0]);
+  stokes_nu1 = point_source.stokes(frequencies[1]);
+  stokes_nu2 = point_source.stokes(frequencies[2]);
+
+  StokesVector result;
+  spectrum.evaluate(frequencies, result);
+
+  BOOST_CHECK_CLOSE(result.I[0], stokes_nu0.I, 1.e-6);
+  BOOST_CHECK_CLOSE(result.I[1], stokes_nu1.I, 1.e-6);
+  BOOST_CHECK_CLOSE(result.I[2], stokes_nu2.I, 1.e-6);
+  BOOST_CHECK_CLOSE(result.Q[0], stokes_nu0.Q, 1.e-6);
+  BOOST_CHECK_CLOSE(result.Q[1], stokes_nu1.Q, 1.e-6);
+  BOOST_CHECK_CLOSE(result.Q[2], stokes_nu2.Q, 1.e-6);
+  BOOST_CHECK_CLOSE(result.U[0], stokes_nu0.U, 1.e-6);
+  BOOST_CHECK_CLOSE(result.U[1], stokes_nu1.U, 1.e-6);
+  BOOST_CHECK_CLOSE(result.U[2], stokes_nu2.U, 1.e-6);
+  BOOST_CHECK_CLOSE(result.V[0], stokes_nu0.V, 1.e-6);
+  BOOST_CHECK_CLOSE(result.V[1], stokes_nu1.V, 1.e-6);
+  BOOST_CHECK_CLOSE(result.V[2], stokes_nu2.V, 1.e-6);
+}
+
+BOOST_FIXTURE_TEST_CASE(normal_log_no_rotation, SpectrumFixture) {
+  spectrum.spectral_terms = {1.0e-6, 2.0e-6, 3.0e-6};
+  spectrum.has_logarithmic_spectral_index = true;
+
+  dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0),
+                                      spectrum.reference_flux);
+
+  point_source.setSpectralTerms(
+      spectrum.reference_frequency, spectrum.has_logarithmic_spectral_index,
+      spectrum.spectral_terms.begin(), spectrum.spectral_terms.end());
+  Stokes stokes_nu0(0, 0, 0, 0);
+  Stokes stokes_nu1(0, 0, 0, 0);
+  Stokes stokes_nu2(0, 0, 0, 0);
+
+  stokes_nu0 = point_source.stokes(frequencies[0]);
+  stokes_nu1 = point_source.stokes(frequencies[1]);
+  stokes_nu2 = point_source.stokes(frequencies[2]);
+
+  StokesVector result;
+  spectrum.evaluate(frequencies, result);
+
+  BOOST_CHECK_CLOSE_FRACTION(result.I[0], stokes_nu0.I, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.I[1], stokes_nu1.I, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.I[2], stokes_nu2.I, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.Q[0], stokes_nu0.Q, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.Q[1], stokes_nu1.Q, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.Q[2], stokes_nu2.Q, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.U[0], stokes_nu0.U, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.U[1], stokes_nu1.U, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.U[2], stokes_nu2.U, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.V[0], stokes_nu0.V, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.V[1], stokes_nu1.V, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.V[2], stokes_nu2.V, 1.e-3);
+}
+
+BOOST_FIXTURE_TEST_CASE(normal_log_with_rotation, SpectrumFixture) {
+  spectrum.spectral_terms = {1.0e-6, 2.0e-6, 3.0e-6};
+  spectrum.has_logarithmic_spectral_index = true;
+  spectrum.has_rotation_measure = true;
+  spectrum.has_logarithmic_spectral_index = true;
+
+  dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0),
+                                      spectrum.reference_flux);
+
+  point_source.setSpectralTerms(
+      spectrum.reference_frequency, spectrum.has_logarithmic_spectral_index,
+      spectrum.spectral_terms.begin(), spectrum.spectral_terms.end());
+
+  point_source.setRotationMeasure(spectrum.polarization_factor,
+                                  spectrum.polarization_angle,
+                                  spectrum.rotation_measure);
+  Stokes stokes_nu0(0, 0, 0, 0);
+  Stokes stokes_nu1(0, 0, 0, 0);
+  Stokes stokes_nu2(0, 0, 0, 0);
+
+  stokes_nu0 = point_source.stokes(frequencies[0]);
+  stokes_nu1 = point_source.stokes(frequencies[1]);
+  stokes_nu2 = point_source.stokes(frequencies[2]);
+
+  StokesVector result;
+  spectrum.evaluate(frequencies, result);
+
+  BOOST_CHECK_CLOSE_FRACTION(result.I[0], stokes_nu0.I, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.I[1], stokes_nu1.I, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.I[2], stokes_nu2.I, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.Q[0], stokes_nu0.Q, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.Q[1], stokes_nu1.Q, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.Q[2], stokes_nu2.Q, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.U[0], stokes_nu0.U, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.U[1], stokes_nu1.U, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.U[2], stokes_nu2.U, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.V[0], stokes_nu0.V, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.V[1], stokes_nu1.V, 1.e-3);
+  BOOST_CHECK_CLOSE_FRACTION(result.V[2], stokes_nu2.V, 1.e-3);
+}
 
 BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file