diff --git a/base/PredictBuffer.h b/base/PredictBuffer.h index 5ec50cf1bd2b5e1b3481966642817ce107ec99a1..a29dd7efae19b544474934f6d247211f75694193 100644 --- a/base/PredictBuffer.h +++ b/base/PredictBuffer.h @@ -19,7 +19,7 @@ namespace base { class PredictBuffer { public: void resize(size_t n_threads, size_t n_correlations, size_t n_channels, - size_t n_baselines, size_t n_stations, bool include_beam, + size_t n_baselines, size_t n_stations_beam, bool include_beam, bool full_beam) { if (include_beam) { // The full buffer is not used when Stokes I is used -- conditionally @@ -32,9 +32,9 @@ class PredictBuffer { for (size_t i = 0; i != n_threads; ++i) { if (full_beam) { - full_beam_values_[i].resize(n_stations * n_channels); + full_beam_values_[i].resize(n_stations_beam * n_channels); } else { - scalar_beam_values_[i].resize(n_stations * n_channels); + scalar_beam_values_[i].resize(n_stations_beam * n_channels); } } } diff --git a/base/Telescope.cc b/base/Telescope.cc index b3e23f823ee9cc5a8763775b211220f67b2be531..9d0c3eb1fa7491657bb0555d1aa388a924e3b9f1 100644 --- a/base/Telescope.cc +++ b/base/Telescope.cc @@ -3,8 +3,6 @@ #include "Telescope.h" -#include <EveryBeam/telescope/phasedarray.h> - // Support both old and new return value types of PhasedArray::GetStation(). // TODO(AST-1001): Remove support for the old type in 2023. namespace { @@ -25,14 +23,14 @@ namespace base { std::vector<size_t> SelectStationIndices( const everybeam::telescope::Telescope& telescope, const std::vector<std::string>& station_names) { + if (IsDish(telescope)) { + // For a dish telescope, the beam of all stations is assumed to be identical + return {0}; + } + + // If the telescope is not a dish, it is a phased array auto phased_array = dynamic_cast<const everybeam::telescope::PhasedArray*>(&telescope); - if (phased_array == nullptr) { - throw std::runtime_error( - "Currently, only PhasedArray telescopes accepted as input, i.e. " - "OSKAR or LOFAR. Support for other telescopes may become available " - "soon."); - } // Copy only those stations for which the name matches. // Note: the order of the station names in both vectors match, diff --git a/base/Telescope.h b/base/Telescope.h index 2d59ae933f5631abb1069c4fe856b4691abcb209..a1cd208a34e5bbe939fe370f45b3bc1944a41eb2 100644 --- a/base/Telescope.h +++ b/base/Telescope.h @@ -5,6 +5,8 @@ #define DP3_BASE_TELESCOPE_H_ #include <EveryBeam/load.h> +#include <EveryBeam/telescope/phasedarray.h> +#include <EveryBeam/telescope/dish.h> namespace dp3 { namespace base { @@ -24,10 +26,20 @@ inline std::unique_ptr<everybeam::telescope::Telescope> GetTelescope( return telescope; } +inline bool IsPhasedArray(const everybeam::telescope::Telescope& telescope) { + auto phased_array = + dynamic_cast<const everybeam::telescope::PhasedArray*>(&telescope); + return phased_array != nullptr; +} + +inline bool IsDish(const everybeam::telescope::Telescope& telescope) { + auto dish = dynamic_cast<const everybeam::telescope::Dish*>(&telescope); + return dish != nullptr; +} + /** * Find stations in a telescope by name and return their indices. * @param telescope The telescope, which contains antennae / stations. - * Only PhasedArray telescopes are currently supported. * @param station_names A list of station names. The order of the names in this * list should match the order in which they occur in the telescope. * @return The indices corresponding to the station names. Because of the diff --git a/docker/py_wheel.docker b/docker/py_wheel.docker index d3cf06547d3b23a5f02caa4debc91ac45ef4c55f..9771425779da90bcb90b53265c2755ae0b3faa9e 100644 --- a/docker/py_wheel.docker +++ b/docker/py_wheel.docker @@ -8,7 +8,7 @@ ARG PYMINOR FROM quay.io/casacore/casacore:py${PYMAJOR}${PYMINOR}_master ARG AOFLAGGER_VERSION=3.2.0 -ARG EVERYBEAM_VERSION=0.5.6 +ARG EVERYBEAM_VERSION=0.5.7 ARG HDF5_VERSION=1.12.2 ARG FFTW_VERSION=3.3.8 ARG LUA_VERSION=5.3.6 diff --git a/docker/ubuntu_20_04_base b/docker/ubuntu_20_04_base index 3bfc8f9bb8b33062e9947762a183f7170ec00720..5f20ff4d922f0f5284f8f6203173546b0c132f49 100644 --- a/docker/ubuntu_20_04_base +++ b/docker/ubuntu_20_04_base @@ -2,7 +2,7 @@ FROM ubuntu:20.04 # TODO: needs to be bumped before next DP3 release # ENV IDG_VERSION=0.8 -ENV EVERYBEAM_VERSION=v0.5.6 +ENV EVERYBEAM_VERSION=v0.5.7 ENV IDG_VERSION=6b61c038883ad3f807d20047c4f9e1a1f0b8d98a ENV AOFLAGGER_VERSION=65d5fba4f4c12797386d3fd9cd76734956a8b233 diff --git a/docker/ubuntu_22_04_base b/docker/ubuntu_22_04_base index 5954af52d4daf0fe24e1a9f844e9ef685db8d061..f3134dd23bb67cd552b0c02aec27f61872b9129a 100644 --- a/docker/ubuntu_22_04_base +++ b/docker/ubuntu_22_04_base @@ -2,7 +2,7 @@ FROM ubuntu:22.04 # TODO: needs to be bumped before next DP3 release # ENV IDG_VERSION=0.8 -ENV EVERYBEAM_VERSION=v0.5.6 +ENV EVERYBEAM_VERSION=v0.5.7 ENV IDG_VERSION=6b61c038883ad3f807d20047c4f9e1a1f0b8d98a ENV AOFLAGGER_VERSION=65d5fba4f4c12797386d3fd9cd76734956a8b233 diff --git a/resources/tDish.MS.tgz b/resources/tDish.MS.tgz new file mode 100644 index 0000000000000000000000000000000000000000..66f9f2172d49155a6aea928755743d129ef6928d Binary files /dev/null and b/resources/tDish.MS.tgz differ diff --git a/steps/ApplyBeam.cc b/steps/ApplyBeam.cc index c4b8d785ca8bd28422f40fa618dea07484ed963b..9ef1f24d9b52b22598d4ce712265fe6df0461a1f 100644 --- a/steps/ApplyBeam.cc +++ b/steps/ApplyBeam.cc @@ -148,8 +148,8 @@ void ApplyBeam::updateInfo(const DPInfo& infoIn) { static_cast<int>(everybeam::CorrectionMode::kNone)); } - const size_t nSt = info().nantenna(); - const size_t nCh = info().nchan(); + const size_t n_stations = info().nantenna(); + const size_t n_channels = info().nchan(); const size_t nThreads = aocommon::ThreadPool::GetInstance().NThreads(); itsBeamValues.resize(nThreads); @@ -161,7 +161,7 @@ void ApplyBeam::updateInfo(const DPInfo& infoIn) { telescopes_.resize(nThreads); for (size_t thread = 0; thread < nThreads; ++thread) { - itsBeamValues[thread].resize(nSt * nCh); + itsBeamValues[thread].resize(n_stations * n_channels); itsMeasFrames[thread].set(info().arrayPosCopy()); itsMeasFrames[thread].set( MEpoch(MVEpoch(info().startTime() / 86400), MEpoch::UTC)); @@ -267,13 +267,12 @@ template <typename T> void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0, float* weight0, const everybeam::vector3r_t& srcdir, const everybeam::telescope::Telescope* telescope, - std::vector<aocommon::MC2x2>& beamValues, bool invert, - everybeam::CorrectionMode mode, bool doUpdateWeights, - std::mutex* mutex) { + std::vector<aocommon::MC2x2>& beam_values, + bool invert, everybeam::CorrectionMode mode, + bool doUpdateWeights, std::mutex* mutex) { // Get the beam values for each station. - const size_t nCh = info.chanFreqs().size(); - const size_t nSt = beamValues.size() / nCh; - const size_t nBl = info.nbaselines(); + const size_t n_channels = info.chanFreqs().size(); + const size_t n_baselines = info.nbaselines(); std::unique_ptr<everybeam::pointresponse::PointResponse> point_response = telescope->GetPointResponse(time); @@ -281,26 +280,29 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0, const std::vector<size_t> station_indices = base::SelectStationIndices(*telescope, info.antennaNames()); + const size_t n_stations = station_indices.size(); + // Apply the beam values of both stations to the ApplyBeamed data. - for (size_t ch = 0; ch < nCh; ++ch) { + for (size_t ch = 0; ch < n_channels; ++ch) { switch (mode) { case everybeam::CorrectionMode::kFull: case everybeam::CorrectionMode::kElement: - // Fill beamValues for channel ch - for (size_t st = 0; st < nSt; ++st) { - beamValues[nCh * st + ch] = point_response->Response( + // Fill beam_values for channel ch + for (size_t st = 0; st < n_stations; ++st) { + beam_values[n_channels * st + ch] = point_response->Response( mode, station_indices[st], info.chanFreqs()[ch], srcdir, mutex); if (invert) { // Terminate if the matrix is not invertible. - [[maybe_unused]] bool status = beamValues[nCh * st + ch].Invert(); + [[maybe_unused]] bool status = + beam_values[n_channels * st + ch].Invert(); assert(status); } } break; case everybeam::CorrectionMode::kArrayFactor: { aocommon::MC2x2 af_tmp; - // Fill beamValues for channel ch - for (size_t st = 0; st < nSt; ++st) { + // Fill beam_values for channel ch + for (size_t st = 0; st < n_stations; ++st) { af_tmp = point_response->Response( mode, station_indices[st], info.chanFreqs()[ch], srcdir, mutex); @@ -308,13 +310,13 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0, af_tmp[0] = 1. / af_tmp[0]; af_tmp[3] = 1. / af_tmp[3]; } - beamValues[nCh * st + ch] = af_tmp; + beam_values[n_channels * st + ch] = af_tmp; } break; } case everybeam::CorrectionMode::kNone: // this should not happen - for (size_t st = 0; st < nSt; ++st) { - beamValues[nCh * st + ch] = aocommon::MC2x2::Unity(); + for (size_t st = 0; st < n_stations; ++st) { + beam_values[n_channels * st + ch] = aocommon::MC2x2::Unity(); } break; } @@ -322,15 +324,23 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0, // Apply beam for channel ch on all baselines // For mode=ARRAY_FACTOR, too much work is done here because we know // that r and l are diagonal - for (size_t bl = 0; bl < nBl; ++bl) { - T* data = data0 + bl * 4 * nCh + ch * 4; + for (size_t bl = 0; bl < n_baselines; ++bl) { + T* data = data0 + bl * 4 * n_channels + ch * 4; const aocommon::MC2x2F mat(data); - const aocommon::MC2x2F left(beamValues[nCh * info.getAnt1()[bl] + ch]); - const aocommon::MC2x2F right(beamValues[nCh * info.getAnt2()[bl] + ch]); + + // If the beam is the same for all stations (i.e. when n_stations = 1), + // all baselines will have the same beam values + size_t index_left = + (n_stations == 1 ? ch : n_channels * info.getAnt1()[bl] + ch); + size_t index_right = + (n_stations == 1 ? ch : n_channels * info.getAnt2()[bl] + ch); + const aocommon::MC2x2F left(beam_values[index_left]); + const aocommon::MC2x2F right(beam_values[index_right]); const aocommon::MC2x2F result = left * mat.MultiplyHerm(right); result.AssignTo(data); if (doUpdateWeights) { - ApplyCal::ApplyWeights(left, right, weight0 + bl * 4 * nCh + ch * 4); + ApplyCal::ApplyWeights(left, right, + weight0 + bl * 4 * n_channels + ch * 4); } } } @@ -340,7 +350,7 @@ template void ApplyBeam::applyBeam( const DPInfo& info, double time, std::complex<double>* data0, float* weight0, const everybeam::vector3r_t& srcdir, const everybeam::telescope::Telescope* telescope, - std::vector<aocommon::MC2x2>& beamValues, bool invert, + std::vector<aocommon::MC2x2>& beam_values, bool invert, everybeam::CorrectionMode mode, bool doUpdateWeights, std::mutex* mutex); void ApplyBeam::applyBeam(const DPInfo& info, double time, @@ -355,13 +365,13 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, bool do_update_weights, std::mutex* mutex) { // Get the beam values for each station. const size_t n_channels = info.chanFreqs().size(); - const size_t n_stations = beam_values.size() / n_channels; std::unique_ptr<everybeam::pointresponse::PointResponse> point_response = telescope->GetPointResponse(time); const std::vector<size_t> station_indices = base::SelectStationIndices(*telescope, info.antennaNames()); + const size_t n_stations = station_indices.size(); // Apply the beam values of both stations to the ApplyBeamed data. @@ -419,10 +429,13 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, for (size_t bl = baseline_range.first; bl < baseline_range.second; ++bl) { std::complex<double>* data = data0 + bl * 4 * n_channels + ch * 4; const aocommon::MC2x2F mat(data); - const aocommon::MC2x2F left( - beam_values[n_channels * info.getAnt1()[bl] + ch]); - const aocommon::MC2x2F right( - beam_values[n_channels * info.getAnt2()[bl] + ch]); + size_t index_left = + (n_stations == 1 ? ch : n_channels * info.getAnt1()[bl] + ch); + size_t index_right = + (n_stations == 1 ? ch : n_channels * info.getAnt2()[bl] + ch); + const aocommon::MC2x2F left(beam_values[index_left]); + const aocommon::MC2x2F right(beam_values[index_right]); + const aocommon::MC2x2F result = left * mat.MultiplyHerm(right); result.AssignTo(data); if (do_update_weights) { @@ -440,37 +453,42 @@ void ApplyBeam::applyBeamStokesIArrayFactor( const DPInfo& info, double time, T* data0, const everybeam::vector3r_t& srcdir, const everybeam::telescope::Telescope* telescope, - std::vector<everybeam::complex_t>& beamValues, bool invert, + std::vector<everybeam::complex_t>& beam_values, bool invert, everybeam::CorrectionMode mode, std::mutex* mutex) { assert(mode == everybeam::CorrectionMode::kArrayFactor); // Get the beam values for each station. - const size_t nCh = info.chanFreqs().size(); - const size_t nSt = beamValues.size() / nCh; - const size_t nBl = info.nbaselines(); + const size_t n_channels = info.chanFreqs().size(); + const size_t n_baselines = info.nbaselines(); std::unique_ptr<everybeam::pointresponse::PointResponse> point_response = telescope->GetPointResponse(time); const std::vector<size_t> station_indices = base::SelectStationIndices(*telescope, info.antennaNames()); + const size_t n_stations = station_indices.size(); // Apply the beam values of both stations to the ApplyBeamed data. - for (size_t ch = 0; ch < nCh; ++ch) { - // Fill beamValues for channel ch - for (size_t st = 0; st < nSt; ++st) { - beamValues[nCh * st + ch] = point_response->Response( + for (size_t ch = 0; ch < n_channels; ++ch) { + // Fill beam_values for channel ch + for (size_t st = 0; st < n_stations; ++st) { + beam_values[n_channels * st + ch] = point_response->Response( everybeam::BeamMode::kArrayFactor, station_indices[st], info.chanFreqs()[ch], srcdir, mutex)[0]; if (invert) { - beamValues[nCh * st + ch] = 1. / beamValues[nCh * st + ch]; + beam_values[n_channels * st + ch] = + 1. / beam_values[n_channels * st + ch]; } } // Apply beam for channel ch on all baselines - for (size_t bl = 0; bl < nBl; ++bl) { - T* data = data0 + bl * nCh + ch; - everybeam::complex_t* left = &(beamValues[nCh * info.getAnt1()[bl]]); - everybeam::complex_t* right = &(beamValues[nCh * info.getAnt2()[bl]]); + for (size_t bl = 0; bl < n_baselines; ++bl) { + T* data = data0 + bl * n_channels + ch; + size_t index_left = + (n_stations == 1 ? 0 : n_channels * info.getAnt1()[bl]); + size_t index_right = + (n_stations == 1 ? 0 : n_channels * info.getAnt2()[bl]); + everybeam::complex_t* left = &(beam_values[index_left]); + everybeam::complex_t* right = &(beam_values[index_right]); data[0] = left[ch] * std::complex<double>(data[0]) * conj(right[ch]); // TODO: update weights? @@ -482,7 +500,7 @@ template void ApplyBeam::applyBeamStokesIArrayFactor( const DPInfo& info, double time, std::complex<double>* data0, const everybeam::vector3r_t& srcdir, const everybeam::telescope::Telescope* telescope, - std::vector<everybeam::complex_t>& beamValues, bool invert, + std::vector<everybeam::complex_t>& beam_values, bool invert, everybeam::CorrectionMode mode, std::mutex* mutex); void ApplyBeam::applyBeamStokesIArrayFactor( @@ -496,13 +514,13 @@ void ApplyBeam::applyBeamStokesIArrayFactor( assert(mode == everybeam::CorrectionMode::kArrayFactor); // Get the beam values for each station. const size_t n_channels = info.chanFreqs().size(); - const size_t n_stations = beam_values.size() / n_channels; std::unique_ptr<everybeam::pointresponse::PointResponse> point_response = telescope->GetPointResponse(time); const std::vector<size_t> station_indices = base::SelectStationIndices(*telescope, info.antennaNames()); + const size_t n_stations = station_indices.size(); // Apply the beam values of both stations to the ApplyBeamed data. if (station_range.first < n_stations) { @@ -525,10 +543,12 @@ void ApplyBeam::applyBeamStokesIArrayFactor( // Apply beam for channel ch on the baselines handeled by this thread for (size_t bl = baseline_range.first; bl < baseline_range.second; ++bl) { std::complex<double>* data = data0 + bl * n_channels + ch; - everybeam::complex_t* left = - &(beam_values[n_channels * info.getAnt1()[bl]]); - everybeam::complex_t* right = - &(beam_values[n_channels * info.getAnt2()[bl]]); + size_t index_left = + (n_stations == 1 ? 0 : n_channels * info.getAnt1()[bl]); + size_t index_right = + (n_stations == 1 ? 0 : n_channels * info.getAnt2()[bl]); + everybeam::complex_t* left = &(beam_values[index_left]); + everybeam::complex_t* right = &(beam_values[index_right]); data[0] = left[ch] * std::complex<double>(data[0]) * conj(right[ch]); // TODO: update weights? diff --git a/steps/OnePredict.cc b/steps/OnePredict.cc index fe4afc05e5e59dd0876dd65021eccc1f64352442..8817738e214e12609277c2c1e3e48e5fc4408e10 100644 --- a/steps/OnePredict.cc +++ b/steps/OnePredict.cc @@ -222,11 +222,14 @@ void OnePredict::initializeThreadData() { if (!predict_buffer_) { predict_buffer_ = std::make_shared<base::PredictBuffer>(); } + bool is_dish_telescope = false; if (apply_beam_) { telescope_ = base::GetTelescope(info().msName(), element_response_model_, use_channel_freq_); + is_dish_telescope = base::IsDish(*telescope_); } - predict_buffer_->resize(nThreads, nCr, nCh, nBl, nSt, apply_beam_, + predict_buffer_->resize(nThreads, nCr, nCh, nBl, + (is_dish_telescope ? 1 : nSt), apply_beam_, !stokes_i_only_); // Create the Measure ITRF conversion info given the array position. // The time and direction are filled in later. diff --git a/steps/test/integration/tApplyBeam.py b/steps/test/integration/tApplyBeam.py index c06c34d62c1344bf03b5ada733be230d0bae120c..278f339da1f4e59ccd49a809dd3b68514303f7ee 100644 --- a/steps/test/integration/tApplyBeam.py +++ b/steps/test/integration/tApplyBeam.py @@ -13,7 +13,7 @@ import sys sys.path.append(".") import testconfig as tcf -from utils import assert_taql, untar_ms +from utils import assert_taql, untar_ms, get_taql_result """ Tests for applying the beam model. @@ -25,9 +25,24 @@ Script can be invoked in two ways: """ MSIN = "tNDPPP-generic.MS" +DISH_MSIN = "tDish.MS" MSAPPLYBEAM = "tApplyBeam.tab" CWD = os.getcwd() +""" +The tDish.MS is a reduced version of a MEERKAT dataset, generated with the following command: + +DP3 \ +msin=MKT_CDFS_2_5_856_963MHz.ms \ +msout=tDish.MS \ +msout.overwrite=true \ +steps=[filter] \ +msin.starttime='29-Jun-2019/05:08:00.581' \ +msin.ntimes=5 \ +msin.nchan=5 \ +filter.blrange="[0,100,0,100]" +""" + @pytest.fixture(autouse=True) def source_env(): @@ -37,6 +52,7 @@ def source_env(): os.chdir(tmpdir) untar_ms(f"{tcf.RESOURCEDIR}/{MSIN}.tgz") + untar_ms(f"{tcf.RESOURCEDIR}/{DISH_MSIN}.tgz") untar_ms(f"{tcf.SRCDIR}/{MSAPPLYBEAM}.tgz") # Tests are executed here @@ -111,3 +127,46 @@ def test_with_updateweights(): ) taql_command = f"select from {MSIN} where all(near(WEIGHT_SPECTRUM, NEW_WEIGHT_SPECTRUM))" assert_taql(taql_command) + + +from testconfig import TAQLEXE + + +def test_dish_beam(): + # Apply the beam with an offset of [0.01,-0.02] to the phase center + check_call( + [ + tcf.DP3EXE, + f"msin={DISH_MSIN}", + f"msout=beam_applied.ms", + "msout.overwrite=true", + "steps=[applybeam]", + "applybeam.usechannelfreq=true", + "applybeam.direction=[0.91848969rad,-0.50149271rad]", + ] + ) + + # Assert that there are 25 elements which are zeroes before and after applying the beam + zero_elements = f"select t1.DATA[0,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[0,0])==0.0 AND abs(t2.DATA[0,0])==0.0" + assert_taql(zero_elements, 25) + + # Assert that all other elements have changed after applying the beam + different_elements = f"select t1.DATA[0,0] from {DISH_MSIN} t1, beam_applied.ms t2 where (abs(t1.DATA[0,0]-t2.DATA[0,0])>1e-1) AND abs(t2.DATA[0,0])!=0.0" + assert_taql(different_elements, 425) + + # Assert that the beam value is correct per each channel + assert_taql( + f"select t1.DATA[0,0]/t2.DATA[0,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[0,0] / t2.DATA[0,0] - 0.146382) > 1e-6 AND abs(t2.DATA[0,0])!=0.0" + ) + assert_taql( + f"select t1.DATA[1,0]/t2.DATA[1,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[1,0] / t2.DATA[1,0] - 0.146003) > 1e-6 AND abs(t2.DATA[1,0])!=0.0" + ) + assert_taql( + f"select t1.DATA[2,0]/t2.DATA[2,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[2,0] / t2.DATA[2,0] - 0.145628) > 1e-6 AND abs(t2.DATA[2,0])!=0.0" + ) + assert_taql( + f"select t1.DATA[3,0]/t2.DATA[3,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[3,0] / t2.DATA[3,0] - 0.145258) > 1e-6 AND abs(t2.DATA[3,0])!=0.0" + ) + assert_taql( + f"select t1.DATA[4,0]/t2.DATA[4,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[4,0] / t2.DATA[4,0] - 0.144933) > 1e-6 AND abs(t2.DATA[4,0])!=0.0" + )