diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 4be7730e0284abfcd06788553dccdea8f62cbc2b..330f0ba9de41f56fdd74aa8899163ff660ea07ba 100644
--- a/.gitlab-ci.yml
+++ b/.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
+
diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
index d739d25f3db8d7ffb0592e2ff4f7f5721893e466..27847c8f4643b5990fdbef202df6118d1afb67da 100644
--- a/benchmark/CMakeLists.txt
+++ b/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>
diff --git a/benchmark/spectrum.cpp b/benchmark/spectrum.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..295b2601c166bae5144563a2c9e7646e39d1d586
--- /dev/null
+++ b/benchmark/spectrum.cpp
@@ -0,0 +1,95 @@
+#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);
diff --git a/include/PointSource.h b/include/PointSource.h
index 3bd6cfcc20b3a2af7cb569147da9d312105975d2..8703442a71941afb3db1664d6eccf469513ce0c8 100644
--- a/include/PointSource.h
+++ b/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:
diff --git a/include/Spectrum.h b/include/Spectrum.h
index c464b90855712299eee74a38d197795fde2a4a56..94f2154f019d8f836dedaee610570549afe2052b 100644
--- a/include/Spectrum.h
+++ b/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> ¶meters) 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> ¶meters) 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> ¶meters) 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> ¶meters) const {
+ double x = frequency / reference_frequency - 1.0;
+
+ return I0 + evaluate_polynomial(x, parameters) * x;
}
};
diff --git a/tests/test_spectrum.cpp b/tests/test_spectrum.cpp
index ba1fa273c3e2af7dba1bca80bbedb94fb92aef0c..74bdaeb95946af02140ec90317d7a021ae3a4362 100644
--- a/tests/test_spectrum.cpp
+++ b/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