From 61a057e9158f73c792ff8b22d460b4b6ecdbbbbe Mon Sep 17 00:00:00 2001
From: mancini <mancini@astron.nl>
Date: Tue, 20 May 2025 14:20:04 +0200
Subject: [PATCH] Refactor
---
include/Directions.h | 2 ++
include/GaussianSource.h | 7 ++--
include/GaussianSourceCollection.h | 17 +++++++++
include/PointSource.h | 32 ++++-------------
include/PointSourceCollection.h | 21 +++++++----
include/Predict.h | 7 +++-
include/Spectrum.h | 4 +--
src/GaussianSource.cpp | 10 ++++++
src/PointSource.cpp | 47 ++++++++++++++-----------
src/Predict.cpp | 35 ++++++++++++++++++
tests/test_gaussiansourcecollection.cpp | 8 ++---
tests/test_pointsourcecollection.cpp | 16 ++++-----
tests/test_spectrum.cpp | 4 +--
13 files changed, 136 insertions(+), 74 deletions(-)
diff --git a/include/Directions.h b/include/Directions.h
index f172324..58a982b 100644
--- a/include/Directions.h
+++ b/include/Directions.h
@@ -51,6 +51,8 @@ public:
dec.clear();
}
+ Direction operator[](size_t i) const { return Direction(ra[i], dec[i]); }
+
size_t Size() const { return ra.size(); }
SmartVector<double> ra;
diff --git a/include/GaussianSource.h b/include/GaussianSource.h
index e12c799..17065a5 100644
--- a/include/GaussianSource.h
+++ b/include/GaussianSource.h
@@ -6,7 +6,7 @@
#ifndef DPPP_GAUSSIANSOURCE_H
#define DPPP_GAUSSIANSOURCE_H
-#include "PointSource.h"
+#include <PointSource.h>
namespace dp3 {
namespace base {
@@ -15,7 +15,7 @@ namespace base {
/// @{
-class GaussianSource : public PointSource {
+class GaussianSource : public dp3::base::PointSource {
public:
typedef std::shared_ptr<GaussianSource> Ptr;
typedef std::shared_ptr<const GaussianSource> ConstPtr;
@@ -23,6 +23,9 @@ public:
GaussianSource(const Direction &direction);
GaussianSource(const Direction &direction, const Stokes &stokes,
size_t beam_id = 0);
+ GaussianSource(const Direction &direction, const Spectrum &spectrum,
+ double position_angle, bool is_position_angle_absolute,
+ double minor_axis, double major_axis, size_t beam_id = 0);
/// Set position angle in radians. The position angle is the smallest angle
/// between the major axis and North, measured positively North over East.
diff --git a/include/GaussianSourceCollection.h b/include/GaussianSourceCollection.h
index 9c3e2b5..3c0c26f 100644
--- a/include/GaussianSourceCollection.h
+++ b/include/GaussianSourceCollection.h
@@ -38,6 +38,23 @@ public:
PointSourceCollection::Reserve(new_size);
}
+ GaussianSource operator[](size_t i) const {
+ return GaussianSource(direction_vector[i], spectrums[i], beam_id[i],
+ position_angle[i], major_axis[i], minor_axis[i],
+ position_angle_is_absolute[i]);
+ }
+
+ std::unique_ptr<GaussianSourceCollection> SelectBeamID(size_t beam_id) {
+ auto selected = std::make_unique<GaussianSourceCollection>();
+
+ for (size_t i = 0; i < Size(); ++i) {
+ if (this->beam_id[i] == beam_id) {
+ selected->Add(operator[](i));
+ }
+ }
+ return std::move(selected);
+ }
+
SmartVector<double> position_angle;
SmartVector<double> major_axis;
SmartVector<double> minor_axis;
diff --git a/include/PointSource.h b/include/PointSource.h
index 154503f..0b41548 100644
--- a/include/PointSource.h
+++ b/include/PointSource.h
@@ -10,7 +10,6 @@
#include "Direction.h"
#include "Spectrum.h"
#include "Stokes.h"
-#include "StokesVector.h"
#include <memory>
#include <set>
namespace dp3 {
@@ -28,6 +27,8 @@ public:
PointSource(const Direction &direction, const Stokes &stokes,
const size_t beam_id = 0);
+ PointSource(const Direction &direction, const Spectrum &stokes,
+ const size_t beam_id = 0);
const Direction &GetDirection() const { return direction_; }
void SetDirection(const Direction &position);
@@ -37,8 +38,7 @@ public:
const Spectrum &GetSpectrum() const { return spectrum_; }
Stokes GetStokes(double freq) const;
- const Stokes &GetStokes() const { return stokes_; }
- const StokesVector &GetStokesVector() const { return stokes_vector_; }
+ const Stokes &GetStokes() const { return spectrum_.GetReferenceFlux(); }
const size_t &GetBeamId() const { return beam_id_; }
void SetBeamId(size_t beam_id) { beam_id_ = beam_id; }
@@ -54,19 +54,7 @@ public:
protected:
Direction direction_;
Spectrum spectrum_;
- StokesVector stokes_vector_;
- Stokes stokes_;
size_t beam_id_;
- // 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_;
};
/// @}
@@ -77,18 +65,10 @@ void PointSource::SetSpectralTerms(double refFreq, bool isLogarithmic, T first,
// 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));
-}
+ spectrum_.SetSpectralTerms(refFreq, isLogarithmic,
+ xt::adapt(std::vector<double>(first, last)));
+};
} // namespace base
} // namespace dp3
diff --git a/include/PointSourceCollection.h b/include/PointSourceCollection.h
index e4cfebf..8171a83 100644
--- a/include/PointSourceCollection.h
+++ b/include/PointSourceCollection.h
@@ -22,15 +22,12 @@ using PointSource = dp3::base::PointSource;
class PointSourceCollection : public ObjectCollection<PointSource> {
public:
Directions direction_vector;
- std::vector<StokesVector> spectrum_vector;
std::vector<Spectrum> spectrums;
std::vector<size_t> beam_id;
std::set<size_t> beam_ids;
- StokesVector stokes_vector;
void Add(const PointSource &point_source) {
direction_vector.Add(point_source.GetDirection());
- stokes_vector.Add(point_source.GetStokes());
spectrums.push_back(point_source.GetSpectrum());
beam_id.push_back(point_source.GetBeamId());
beam_ids.insert(point_source.GetBeamId());
@@ -38,20 +35,32 @@ public:
void Reserve(size_t new_size) {
direction_vector.Reserve(new_size);
- stokes_vector.Reserve(new_size);
spectrums.reserve(new_size);
beam_id.reserve(new_size);
}
void Clear() {
- PointSourceCollection::Clear();
direction_vector.Clear();
- spectrum_vector.clear();
spectrums.clear();
beam_id.clear();
beam_ids.clear();
}
+ PointSource operator[](size_t i) const {
+ return PointSource(direction_vector[i], spectrums[i], beam_id[i]);
+ }
+
+ std::unique_ptr<PointSourceCollection> SelectBeamID(size_t beam_id) {
+ auto selected = std::make_unique<PointSourceCollection>();
+
+ for (size_t i = 0; i < Size(); ++i) {
+ if (this->beam_id[i] == beam_id) {
+ selected->Add(operator[](i));
+ }
+ }
+ return std::move(selected);
+ }
+
size_t Size() const { return direction_vector.Size(); }
};
diff --git a/include/Predict.h b/include/Predict.h
index 0851580..08f2579 100644
--- a/include/Predict.h
+++ b/include/Predict.h
@@ -86,7 +86,7 @@ inline void spectrum(const PointSource &component, size_t nChannel,
/// @{
/**
- * @brief Simulator to compute visibilities given a sky model
+ * @brief Predict class to compute visibilities given a sky model
*
* This class computes visibilities given model components in the sky.
* Effectively, it evaluates the equation:
@@ -126,6 +126,11 @@ public:
xt::xtensor<dcomplex, 3> &buffer) const;
void run(PredictPlanExec &plan, const GaussianSourceCollection &sources,
xt::xtensor<dcomplex, 3> &buffer) const;
+
+ void run(PredictPlanExec &plan, const GaussianSourceCollection &sources,
+ xt::xtensor<dcomplex, 4> &buffer) const;
+ void run(PredictPlanExec &plan, const PointSourceCollection &sources,
+ xt::xtensor<dcomplex, 4> &buffer) const;
};
/// @}
diff --git a/include/Spectrum.h b/include/Spectrum.h
index 3fabc95..369a4dd 100644
--- a/include/Spectrum.h
+++ b/include/Spectrum.h
@@ -11,7 +11,6 @@
#include <Stokes.h>
#include <StokesVector.h>
-
#include <casacore/casa/BasicSL/Constants.h>
#include <complex>
#include <xtensor/xcomplex.hpp>
@@ -91,7 +90,7 @@ public:
return has_logarithmic_spectral_index_;
}
bool HasSpectralTerms() const { return spectral_terms_.size() > 0; }
-
+ bool HasRotationMeasure() const { return has_rotation_measure_; }
void SetPolarization(double angle, double factor) {
polarization_angle_ = angle;
polarization_factor_ = factor;
@@ -207,7 +206,6 @@ private:
return EvaluatePolynomial(x, parameters) * x;
}
-private:
Stokes reference_flux_;
xt::xtensor<double, 1> spectral_terms_;
diff --git a/src/GaussianSource.cpp b/src/GaussianSource.cpp
index 2e558aa..3d67e40 100644
--- a/src/GaussianSource.cpp
+++ b/src/GaussianSource.cpp
@@ -20,6 +20,16 @@ GaussianSource::GaussianSource(const Direction &direction, const Stokes &stokes,
: PointSource(direction, stokes, beam_id), position_angle_(0.0),
is_position_angle_absolute_(true), major_axis_(0.0), minor_axis_(0.0) {}
+GaussianSource::GaussianSource(const Direction &direction,
+ const Spectrum &spectrum, double position_angle,
+ bool is_position_angle_absolute,
+ double minor_axis, double major_axis,
+ size_t beam_id)
+ : PointSource(direction, spectrum, beam_id),
+ position_angle_(position_angle),
+ is_position_angle_absolute_(is_position_angle_absolute),
+ major_axis_(major_axis), minor_axis_(minor_axis) {}
+
void GaussianSource::SetPositionAngle(double angle) { position_angle_ = angle; }
void GaussianSource::SetMajorAxis(double fwhm) { major_axis_ = fwhm; }
diff --git a/src/PointSource.cpp b/src/PointSource.cpp
index 569edda..88c1ef9 100644
--- a/src/PointSource.cpp
+++ b/src/PointSource.cpp
@@ -18,13 +18,14 @@ namespace base {
PointSource::PointSource(const Direction &position, const Stokes &stokes,
const size_t beam_id)
- : direction_(position), stokes_(stokes), reference_frequency_(0.0),
- polarizated_fraction_(0.0), polarization_angle_(0.0),
- rotation_measure_(0.0), has_rotation_measure_(false),
- has_logarithmic_si_(true), beam_id_(beam_id) {
+ : direction_(position), beam_id_(beam_id) {
spectrum_.SetReferenceFlux(stokes);
}
+PointSource::PointSource(const Direction &position, const Spectrum &spectrum,
+ const size_t beam_id)
+ : direction_(position), beam_id_(beam_id), spectrum_(spectrum) {}
+
void PointSource::SetDirection(const Direction &direction) {
direction_ = direction;
}
@@ -36,31 +37,30 @@ void PointSource::ComputeSpectrum(
}
void PointSource::SetRotationMeasure(double fraction, double angle, double rm) {
- polarizated_fraction_ = fraction;
- polarization_angle_ = angle;
- rotation_measure_ = rm;
- has_rotation_measure_ = true;
+ spectrum_.SetPolarization(angle, fraction);
+ spectrum_.SetRotationMeasure(rm);
}
// FIXME: legacy code, remove it later
Stokes PointSource::GetStokes(double freq) const {
- Stokes stokes(stokes_);
+ Stokes stokes(spectrum_.GetReferenceFlux());
if (HasSpectralTerms()) {
- if (has_logarithmic_si_) {
+ if (spectrum_.HasLogarithmicSpectralIndex()) {
// Compute spectral index as:
// (v / v0) ^ (c0 + c1 * log10(v / v0) + c2 * log10(v / v0)^2 + ...)
// Where v is the frequency and v0 is the reference frequency.
// Compute log10(v / v0).
- double base = log10(freq) - log10(reference_frequency_);
+ double base = log10(freq) - log10(spectrum_.GetReferenceFrequency());
// Compute c0 + log10(v / v0) * c1 + log10(v / v0)^2 * c2 + ...
// using Horner's rule.
double exponent = 0.0;
- typedef std::vector<double>::const_reverse_iterator iterator_type;
- for (iterator_type it = spectral_terms_.rbegin(),
- end = spectral_terms_.rend();
+ typedef xt::xtensor<double, 1>::const_reverse_iterator iterator_type;
+ auto spectral_terms = spectrum_.GetSpectralTerms();
+ for (iterator_type it = spectral_terms.rbegin(),
+ end = spectral_terms.rend();
it != end; ++it) {
exponent = exponent * base + *it;
}
@@ -70,9 +70,10 @@ Stokes PointSource::GetStokes(double freq) const {
stokes.I *= pow(10., base * exponent);
stokes.V *= pow(10., base * exponent);
} else {
- double x = freq / reference_frequency_ - 1.0;
- typedef std::vector<double>::const_reverse_iterator iterator_type;
+ double x = freq / spectrum_.GetReferenceFrequency() - 1.0;
+ typedef xt::xtensor<double, 1>::const_reverse_iterator iterator_type;
double val = 0.0;
+ auto spectral_terms_ = spectrum_.GetSpectralTerms();
for (iterator_type it = spectral_terms_.rbegin(),
end = spectral_terms_.rend();
it != end; ++it) {
@@ -85,18 +86,22 @@ Stokes PointSource::GetStokes(double freq) const {
if (HasRotationMeasure()) {
double lambda = casacore::C::c / freq;
- double chi =
- 2.0 * (polarization_angle_ + rotation_measure_ * lambda * lambda);
- double stokesQU = stokes.I * polarizated_fraction_;
+ double chi = 2.0 * (spectrum_.GetPolarizationAngle() +
+ spectrum_.GetRotationMeasure() * lambda * lambda);
+ double stokesQU = stokes.I * spectrum_.GetPolarizationFactor();
stokes.Q = stokesQU * cos(chi);
stokes.U = stokesQU * sin(chi);
}
return stokes;
}
-bool PointSource::HasSpectralTerms() const { return !spectral_terms_.empty(); }
+bool PointSource::HasSpectralTerms() const {
+ return spectrum_.HasSpectralTerms();
+}
-bool PointSource::HasRotationMeasure() const { return has_rotation_measure_; }
+bool PointSource::HasRotationMeasure() const {
+ return spectrum_.HasRotationMeasure();
+}
} // namespace base
} // namespace dp3
diff --git a/src/Predict.cpp b/src/Predict.cpp
index a05045a..1ed984d 100644
--- a/src/Predict.cpp
+++ b/src/Predict.cpp
@@ -34,5 +34,40 @@ void Predict::run(PredictPlanExec &plan,
plan.Compute(sources, buffer);
}
+void Predict::run(PredictPlanExec &plan,
+ const GaussianSourceCollection &sources,
+ xt::xtensor<dcomplex, 4> &buffer) const {
+ plan.Precompute(sources);
+
+ xt::xtensor<dcomplex, 3>::shape_type tmp_buffer_shape = {
+ buffer.shape(1), buffer.shape(2), buffer.shape(3)};
+ for (auto &beam_id : sources.beam_ids) {
+ xt::xtensor<dcomplex, 3> direction_buffer(tmp_buffer_shape);
+
+ // TODO compute should be extended to accumulate per beam_id
+ plan.Compute(sources, direction_buffer);
+ xt::view(buffer, beam_id, xt::all(), xt::all(), xt::all()) =
+ direction_buffer;
+ }
+}
+
+void Predict::run(PredictPlanExec &plan, const PointSourceCollection &sources,
+ xt::xtensor<dcomplex, 4> &buffer) const {
+ // buffer dimensions are (nBean, nFreq, nBaselines, nSources)
+ plan.Precompute(sources);
+
+ xt::xtensor<dcomplex, 3>::shape_type tmp_buffer_shape = {
+ buffer.shape(1), buffer.shape(2), buffer.shape(3)};
+ for (auto &beam_id : sources.beam_ids) {
+ xt::xtensor<dcomplex, 3> direction_buffer(tmp_buffer_shape);
+
+ // TODO compute should be extended to accumulate per beam_id
+ plan.Compute(sources, direction_buffer);
+ xt::view(buffer, beam_id, xt::all(), xt::all(), xt::all()) =
+ direction_buffer;
+ }
+}
+// plan.Compute(sources, buffer);}
+
} // namespace base
} // namespace dp3
diff --git a/tests/test_gaussiansourcecollection.cpp b/tests/test_gaussiansourcecollection.cpp
index e797036..3749a1c 100644
--- a/tests/test_gaussiansourcecollection.cpp
+++ b/tests/test_gaussiansourcecollection.cpp
@@ -14,10 +14,10 @@ BOOST_AUTO_TEST_CASE(add_single) {
collection.Add(source);
BOOST_CHECK_EQUAL(collection.Size(), 1);
BOOST_CHECK_EQUAL(collection.direction_vector.Size(), 1);
- BOOST_CHECK_EQUAL(collection.stokes_vector.Size(), 1);
+
BOOST_CHECK_EQUAL(collection.direction_vector.ra[0], 0.1);
BOOST_CHECK_EQUAL(collection.direction_vector.dec[0], 0.2);
- BOOST_CHECK_EQUAL(collection.stokes_vector.I[0], 1.0);
+ BOOST_CHECK_EQUAL(collection.spectrums[0].GetReferenceFlux().I, 1.0);
}
BOOST_AUTO_TEST_CASE(reserve) {
@@ -27,13 +27,13 @@ BOOST_AUTO_TEST_CASE(reserve) {
PointSourceCollection collection;
collection.Reserve(2);
- auto ptr = collection.stokes_vector.I.data();
+ auto ptr = collection.direction_vector.ra.data();
collection.Add(sources[0]);
collection.Add(sources[1]);
BOOST_CHECK_EQUAL(collection.Size(), 2);
- BOOST_CHECK_EQUAL(collection.stokes_vector.I.data(), ptr);
+ BOOST_CHECK_EQUAL(collection.direction_vector.ra.data(), ptr);
}
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
diff --git a/tests/test_pointsourcecollection.cpp b/tests/test_pointsourcecollection.cpp
index dcddb1e..2d48acc 100644
--- a/tests/test_pointsourcecollection.cpp
+++ b/tests/test_pointsourcecollection.cpp
@@ -14,10 +14,10 @@ BOOST_AUTO_TEST_CASE(add_single) {
collection.Add(source);
BOOST_CHECK_EQUAL(collection.Size(), 1);
BOOST_CHECK_EQUAL(collection.direction_vector.Size(), 1);
- BOOST_CHECK_EQUAL(collection.stokes_vector.Size(), 1);
+ BOOST_CHECK_EQUAL(collection.spectrums.size(), 1);
BOOST_CHECK_EQUAL(collection.direction_vector.ra[0], 0.1);
BOOST_CHECK_EQUAL(collection.direction_vector.dec[0], 0.2);
- BOOST_CHECK_EQUAL(collection.stokes_vector.I[0], 1.0);
+ BOOST_CHECK_EQUAL(collection.spectrums[0].GetReferenceFlux().I, 1.0);
}
BOOST_AUTO_TEST_CASE(add_single_unspecified_beam_id) {
@@ -26,10 +26,10 @@ BOOST_AUTO_TEST_CASE(add_single_unspecified_beam_id) {
collection.Add(source);
BOOST_CHECK_EQUAL(collection.Size(), 1);
BOOST_CHECK_EQUAL(collection.direction_vector.Size(), 1);
- BOOST_CHECK_EQUAL(collection.stokes_vector.Size(), 1);
+ BOOST_CHECK_EQUAL(collection.spectrums.size(), 1);
BOOST_CHECK_EQUAL(collection.direction_vector.ra[0], 0.1);
BOOST_CHECK_EQUAL(collection.direction_vector.dec[0], 0.2);
- BOOST_CHECK_EQUAL(collection.stokes_vector.I[0], 1.0);
+
BOOST_CHECK_EQUAL(collection.beam_id[0], 0);
}
@@ -39,10 +39,10 @@ BOOST_AUTO_TEST_CASE(add_single_with_beam_id) {
collection.Add(source);
BOOST_CHECK_EQUAL(collection.Size(), 1);
BOOST_CHECK_EQUAL(collection.direction_vector.Size(), 1);
- BOOST_CHECK_EQUAL(collection.stokes_vector.Size(), 1);
+ BOOST_CHECK_EQUAL(collection.spectrums.size(), 1);
BOOST_CHECK_EQUAL(collection.direction_vector.ra[0], 0.1);
BOOST_CHECK_EQUAL(collection.direction_vector.dec[0], 0.2);
- BOOST_CHECK_EQUAL(collection.stokes_vector.I[0], 1.0);
+ BOOST_CHECK_EQUAL(collection.spectrums[0].GetReferenceFlux().I, 1.0);
BOOST_CHECK_EQUAL(collection.beam_id[0], 15);
BOOST_CHECK_EQUAL(collection.beam_ids.size(), 1);
}
@@ -69,13 +69,13 @@ BOOST_AUTO_TEST_CASE(reserve) {
PointSourceCollection collection;
collection.Reserve(2);
- auto ptr = collection.stokes_vector.I.data();
+ auto ptr = collection.direction_vector.ra.data();
collection.Add(sources[0]);
collection.Add(sources[1]);
BOOST_CHECK_EQUAL(collection.Size(), 2);
- BOOST_CHECK_EQUAL(collection.stokes_vector.I.data(), ptr);
+ BOOST_CHECK_EQUAL(collection.direction_vector.ra.data(), ptr);
}
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
diff --git a/tests/test_spectrum.cpp b/tests/test_spectrum.cpp
index 10badfb..2131121 100644
--- a/tests/test_spectrum.cpp
+++ b/tests/test_spectrum.cpp
@@ -24,13 +24,11 @@ struct SpectrumFixture {
};
BOOST_FIXTURE_TEST_CASE(normal_no_rotation, SpectrumFixture) {
- dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0),
- spectrum.GetReferenceFlux());
+ dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0), spectrum);
point_source.SetSpectralTerms(
spectrum.GetReferenceFrequency(), spectrum.HasLogarithmicSpectralIndex(),
spectrum.GetSpectralTerms().begin(), spectrum.GetSpectralTerms().end());
- point_source.SetRotationMeasure(0.0, 0.0, 0.0);
Stokes stokes_nu0(0, 0, 0, 0);
Stokes stokes_nu1(0, 0, 0, 0);
--
GitLab