From 167ae7730c85bf8d88a56cca143567f2e04ad608 Mon Sep 17 00:00:00 2001
From: Jakob Maljaars <jakob.maljaars@stcorp.nl>
Date: Tue, 14 Jul 2020 17:24:24 +0200
Subject: [PATCH] Resolve "Add GriddedResponse test"
---
cpp/LofarMetaDataUtil.cc | 2 +-
cpp/antenna.h | 4 +-
cpp/beamformer.cc | 41 +++---
cpp/common/types.h | 8 +-
cpp/gridded_response/griddedresponse.h | 10 +-
cpp/gridded_response/lofargrid.cc | 10 +-
cpp/gridded_response/lofargrid.h | 8 +-
cpp/telescope/lofar.cc | 155 +-------------------
cpp/test/tload_lofar.cc | 78 ++++++++--
demo/comparison-oskar/main.cpp | 3 +-
{demo/comparison-oskar => external}/npy.hpp | 0
11 files changed, 121 insertions(+), 198 deletions(-)
rename {demo/comparison-oskar => external}/npy.hpp (100%)
diff --git a/cpp/LofarMetaDataUtil.cc b/cpp/LofarMetaDataUtil.cc
index 55756c04..826467e5 100644
--- a/cpp/LofarMetaDataUtil.cc
+++ b/cpp/LofarMetaDataUtil.cc
@@ -219,7 +219,7 @@ BeamFormer::Ptr readAntennaField(const Table &table, unsigned int id,
assert(aips_flag.shape().isEqual(IPosition(2, 2, aips_offset.ncolumn())));
TileConfig tile_config;
- if (name != "LBA") readTileConfig(table, id);
+ if (name != "LBA") tile_config = readTileConfig(table, id);
transformToFieldCoordinates(tile_config, coordinate_system.axes);
for (size_t i = 0; i < aips_offset.ncolumn(); ++i) {
diff --git a/cpp/antenna.h b/cpp/antenna.h
index b4b8ed54..029b840e 100644
--- a/cpp/antenna.h
+++ b/cpp/antenna.h
@@ -135,9 +135,9 @@ class Antenna {
.rotate = options.rotate,
.east = TransformToLocalDirection(options.east),
.north = TransformToLocalDirection(options.north)};
- matrix22c_t Response =
+ matrix22c_t response =
LocalResponse(time, freq, local_direction, local_options);
- return Response;
+ return response;
}
/**
diff --git a/cpp/beamformer.cc b/cpp/beamformer.cc
index f9adb94d..66e260d9 100644
--- a/cpp/beamformer.cc
+++ b/cpp/beamformer.cc
@@ -25,16 +25,17 @@ std::vector<std::complex<double>> BeamFormer::ComputeGeometricResponse(
const double freq, const vector3r_t &direction) const {
// Initialize and fill result vector by looping over antennas
std::vector<std::complex<double>> result(antennas_.size());
- for (auto &antenna : antennas_) {
- const double dl = direction[0] * (antenna->m_phase_reference_position[0] -
- local_phase_reference_position_[0]) +
- direction[1] * (antenna->m_phase_reference_position[1] -
- local_phase_reference_position_[1]) +
- direction[2] * (antenna->m_phase_reference_position[2] -
- local_phase_reference_position_[2]);
+ for (std::size_t idx = 0; idx < antennas_.size(); ++idx) {
+ const double dl =
+ direction[0] * (antennas_[idx]->m_phase_reference_position[0] -
+ local_phase_reference_position_[0]) +
+ direction[1] * (antennas_[idx]->m_phase_reference_position[1] -
+ local_phase_reference_position_[1]) +
+ direction[2] * (antennas_[idx]->m_phase_reference_position[2] -
+ local_phase_reference_position_[2]);
double phase = -2 * M_PI * dl / (common::c / freq);
- result.push_back({std::sin(phase), std::cos(phase)});
+ result[idx] = {std::sin(phase), std::cos(phase)};
}
return result;
}
@@ -48,23 +49,20 @@ BeamFormer::ComputeWeights(const vector3r_t &pointing, double freq) const {
double weight_sum[2] = {0.0, 0.0};
std::vector<std::pair<std::complex<double>, std::complex<double>>> result(
geometric_response.size());
- for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size();
- ++antenna_idx) {
+ for (std::size_t idx = 0; idx < antennas_.size(); ++idx) {
// Compute conjugate of geometric response
- auto phasor_conj = std::conj(geometric_response[antenna_idx]);
+ auto phasor_conj = std::conj(geometric_response[idx]);
// Compute the delays in x/y direction
- result.push_back(
- {phasor_conj * (1.0 * antennas_[antenna_idx]->m_enabled[0]),
- phasor_conj * (1.0 * antennas_[antenna_idx]->m_enabled[1])});
- weight_sum[0] += (1.0 * antennas_[antenna_idx]->m_enabled[0]);
- weight_sum[1] += (1.0 * antennas_[antenna_idx]->m_enabled[1]);
+ result[idx] = {phasor_conj * (1.0 * antennas_[idx]->m_enabled[0]),
+ phasor_conj * (1.0 * antennas_[idx]->m_enabled[1])};
+ weight_sum[0] += (1.0 * antennas_[idx]->m_enabled[0]);
+ weight_sum[1] += (1.0 * antennas_[idx]->m_enabled[1]);
}
// Normalize the weight by the number of antennas
- for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size();
- ++antenna_idx) {
- result[antenna_idx].first /= weight_sum[0];
- result[antenna_idx].second /= weight_sum[1];
+ for (std::size_t idx = 0; idx < antennas_.size(); ++idx) {
+ result[idx].first /= weight_sum[0];
+ result[idx].second /= weight_sum[1];
}
return result;
@@ -79,7 +77,7 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
auto geometric_response = ComputeGeometricResponse(freq, direction);
matrix22c_t result = {0};
- for (unsigned int antenna_idx = 0; antenna_idx < antennas_.size();
+ for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size();
++antenna_idx) {
auto antenna = antennas_[antenna_idx];
auto antenna_weight = weights[antenna_idx];
@@ -87,7 +85,6 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
matrix22c_t antenna_response =
antenna->Response(time, freq, direction, options);
-
result[0][0] += antenna_weight.first * antenna_geometric_reponse *
antenna_response[0][0];
result[0][1] += antenna_weight.first * antenna_geometric_reponse *
diff --git a/cpp/common/types.h b/cpp/common/types.h
index 0384b2ca..35a440fc 100644
--- a/cpp/common/types.h
+++ b/cpp/common/types.h
@@ -96,7 +96,13 @@ typedef std::array<vector3r_t, 16> TileConfig;
template <typename T, size_t N>
std::ostream &operator<<(std::ostream &out, const std::array<T, N> &obj) {
- print(out, obj.begin(), obj.end());
+ // print(out, obj.begin(), obj.end());
+ out << "[";
+ for (auto it : obj) {
+ out << it;
+ if (it != *obj.rbegin()) out << ", ";
+ }
+ out << "]\n";
return out;
}
diff --git a/cpp/gridded_response/griddedresponse.h b/cpp/gridded_response/griddedresponse.h
index 35abf9d5..b152850c 100644
--- a/cpp/gridded_response/griddedresponse.h
+++ b/cpp/gridded_response/griddedresponse.h
@@ -51,7 +51,8 @@ class GriddedResponse {
/**
* @brief Compute the Beam response for a single station
*
- * @param buffer Output buffer
+ * @param buffer Output buffer, compute and set size with
+ * GriddedResponse::GetBufferSize(1)
* @param station_idx Station index, must be smaller than number of stations
* in the Telescope
* @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
@@ -63,13 +64,18 @@ class GriddedResponse {
/**
* @brief Compute the Beam response for all stations in a Telescope
*
- * @param buffer Output buffer
+ * @param buffer Output buffer, compute and set size with
+ * GriddedResponse::GetBufferSize()
* @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
* @param frequency Frequency (Hz)
*/
virtual bool CalculateAllStations(std::complex<float>* buffer, double time,
double frequency) = 0;
+ std::size_t GetBufferSize(std::size_t nstations) {
+ return std::size_t(nstations * width_ * height_ * 2 * 2);
+ }
+
protected:
/**
* @brief Construct a new Gridded Response object
diff --git a/cpp/gridded_response/lofargrid.cc b/cpp/gridded_response/lofargrid.cc
index e42bfc40..b6884670 100644
--- a/cpp/gridded_response/lofargrid.cc
+++ b/cpp/gridded_response/lofargrid.cc
@@ -4,7 +4,7 @@
#include <aocommon/imagecoordinates.h>
#include <cmath>
-// #include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
+#include <iostream>
using namespace everybeam::gridded_response;
@@ -41,7 +41,7 @@ bool LOFARGrid::CalculateStation(std::complex<float>* buffer, double time,
}
for (size_t y = 0; y != height_; ++y) {
- lane.write(Job{.y = y, .antenna_idx = station_idx});
+ lane.write(Job{.y = y, .antenna_idx = station_idx, .buffer_offset = 0});
}
lane.write_end();
@@ -73,7 +73,8 @@ bool LOFARGrid::CalculateAllStations(std::complex<float>* buffer, double time,
for (size_t y = 0; y != height_; ++y) {
for (size_t antenna_idx = 0; antenna_idx != telescope_->GetNrStations();
++antenna_idx) {
- lane.write(Job{.y = y, .antenna_idx = antenna_idx});
+ lane.write(Job{
+ .y = y, .antenna_idx = antenna_idx, .buffer_offset = antenna_idx});
}
}
@@ -135,7 +136,8 @@ void LOFARGrid::CalcThread(std::complex<float>* buffer, double time,
std::complex<float>* base_buffer = buffer + (x + job.y * height_) * 4;
std::complex<float>* ant_buffer_ptr =
- base_buffer + job.antenna_idx * values_per_ant;
+ base_buffer + job.buffer_offset * values_per_ant;
+
matrix22c_t gain_matrix = telescope_->GetStation(job.antenna_idx)
->Response(time, frequency, itrf_direction,
sb_freq, station0_, tile0_);
diff --git a/cpp/gridded_response/lofargrid.h b/cpp/gridded_response/lofargrid.h
index ea524aa6..6ab3eb99 100644
--- a/cpp/gridded_response/lofargrid.h
+++ b/cpp/gridded_response/lofargrid.h
@@ -53,7 +53,8 @@ class LOFARGrid final : public GriddedResponse {
/**
* @brief Compute the Beam response for a single station
*
- * @param buffer Output buffer
+ * @param buffer Output buffer, compute and set size with
+ * GriddedResponse::GetBufferSize(1)
* @param stationIdx Station index, must be smaller than number of stations
* in the Telescope
* @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
@@ -65,7 +66,8 @@ class LOFARGrid final : public GriddedResponse {
/**
* @brief Compute the Beam response for all stations in a Telescope
*
- * @param buffer Output buffer
+ * @param buffer Output buffer, compute and set size with
+ * GriddedResponse::GetBufferSize()
* @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
* @param freq Frequency (Hz)
*/
@@ -82,7 +84,7 @@ class LOFARGrid final : public GriddedResponse {
std::vector<std::thread> threads_;
struct Job {
- size_t y, antenna_idx;
+ size_t y, antenna_idx, buffer_offset;
};
aocommon::Lane<Job>* lane_;
diff --git a/cpp/telescope/lofar.cc b/cpp/telescope/lofar.cc
index e1b02a6b..3cb55809 100644
--- a/cpp/telescope/lofar.cc
+++ b/cpp/telescope/lofar.cc
@@ -2,6 +2,7 @@
#include "./../gridded_response/lofargrid.h"
#include "./../common/math_utils.h"
#include "./../common/casa_utils.h"
+#include "../LofarMetaDataUtil.h"
#include <aocommon/banddata.h>
#include <cassert>
@@ -11,112 +12,6 @@ using namespace everybeam;
using namespace everybeam::telescope;
using namespace casacore;
-// LOFAR Telescope specific utils
-namespace {
-
-constexpr Antenna::CoordinateSystem::Axes lofar_antenna_orientation = {
- {
- -std::sqrt(.5),
- -std::sqrt(.5),
- 0.0,
- },
- {
- std::sqrt(.5),
- -std::sqrt(.5),
- 0.0,
- },
- {0.0, 0.0, 1.0},
-};
-
-// TODO: move to common utils?
-vector3r_t transformToFieldCoordinates(
- const vector3r_t &position, const Antenna::CoordinateSystem::Axes &axes) {
- const vector3r_t result{dot(position, axes.p), dot(position, axes.q),
- dot(position, axes.r)};
- return result;
-}
-
-// // Seems LOFAR specific?
-// void transformToFieldCoordinates(TileConfig &config,
-// const Antenna::CoordinateSystem::Axes &axes)
-// {
-// for (unsigned int i = 0; i < config.size(); ++i) {
-// const vector3r_t position = config[i];
-// config[i][0] = dot(position, axes.p);
-// config[i][1] = dot(position, axes.q);
-// config[i][2] = dot(position, axes.r);
-// }
-// }
-
-// LOFAR specific? Or is this generic for each telescope?
-BeamFormer::Ptr readAntennaField(const Table &table, std::size_t id,
- ElementResponse::Ptr element_response) {
- Antenna::CoordinateSystem coordinate_system =
- common::readCoordinateSystem(table, id);
-
- BeamFormer::Ptr beam_former(new BeamFormer(coordinate_system));
-
- ROScalarColumn<String> c_name(table, "NAME");
- ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");
- ROArrayColumn<Bool> c_flag(table, "ELEMENT_FLAG");
-
- const string &name = c_name(id);
-
- // Read element offsets and flags.
- Matrix<Quantity> aips_offset = c_offset(id);
- assert(aips_offset.shape().isEqual(IPosition(2, 3, aips_offset.ncolumn())));
-
- Matrix<Bool> aips_flag = c_flag(id);
- assert(aips_flag.shape().isEqual(IPosition(2, 2, aips_offset.ncolumn())));
-
- // TileConfig tile_config;
- // if(name != "LBA") readTileConfig(table, id);
- // transformToFieldCoordinates(tile_config, coordinate_system.axes);
-
- for (size_t i = 0; i < aips_offset.ncolumn(); ++i) {
- vector3r_t antenna_position = {aips_offset(0, i).getValue(),
- aips_offset(1, i).getValue(),
- aips_offset(2, i).getValue()};
- antenna_position =
- transformToFieldCoordinates(antenna_position, coordinate_system.axes);
-
- Antenna::Ptr antenna;
- Antenna::CoordinateSystem antenna_coordinate_system{
- antenna_position, lofar_antenna_orientation};
- if (name == "LBA") {
- antenna = Element::Ptr(
- new Element(antenna_coordinate_system, element_response, id));
- } else {
- antenna = Element::Ptr(
- new Element(antenna_coordinate_system, element_response, id));
-
- // TODO
- // HBA, HBA0, HBA1
- // antenna = make_tile(name, coordinate_system, tile_config,
- // element_response);
- }
-
- antenna->m_enabled[0] = !aips_flag(0, i);
- antenna->m_enabled[1] = !aips_flag(1, i);
- beam_former->AddAntenna(antenna);
- }
- return beam_former;
-}
-
-vector3r_t readStationPhaseReference(const Table &table, unsigned int id) {
- vector3r_t phase_reference = {0.0, 0.0, 0.0};
- const string columnName("LOFAR_PHASE_REFERENCE");
- if (common::hasColumn(table, columnName)) {
- ROScalarMeasColumn<MPosition> c_reference(table, columnName);
- MPosition mReference =
- MPosition::Convert(c_reference(id), MPosition::ITRF)();
- MVPosition mvReference = mReference.getValue();
- phase_reference = {mvReference(0), mvReference(1), mvReference(2)};
- }
- return phase_reference;
-}
-} // namespace
-
LOFAR::LOFAR(MeasurementSet &ms, const ElementResponseModel model,
const Options &options)
: Telescope(ms, model, options) {
@@ -163,52 +58,6 @@ std::unique_ptr<gridded_response::GriddedResponse> LOFAR::GetGriddedResponse(
Station::Ptr LOFAR::ReadStation(const MeasurementSet &ms, std::size_t id,
const ElementResponseModel model) const {
- ROMSAntennaColumns antenna(ms.antenna());
- assert(nstations_ > id && !antenna.flagRow()(id));
-
- // Get station name.
- const string name(antenna.name()(id));
-
- // Get station position (ITRF).
- MPosition mPosition =
- MPosition::Convert(antenna.positionMeas()(id), MPosition::ITRF)();
- MVPosition mvPosition = mPosition.getValue();
- const vector3r_t position = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
-
- // Create station.
- Station::Ptr station(new Station(name, position, model));
-
- // Read phase reference position (if available).
- station->setPhaseReference(readStationPhaseReference(ms.antenna(), id));
-
- Table tab_field = common::getSubTable(ms, "LOFAR_ANTENNA_FIELD");
- tab_field = tab_field(tab_field.col("ANTENNA_ID") == static_cast<Int>(id));
-
- // The Station will consist of a BeamFormer that combines the fields
- // coordinate system is ITRF
- // phase reference is station position
- auto beam_former = BeamFormer::Ptr(new BeamFormer(
- Antenna::IdentityCoordinateSystem, station->phaseReference()));
-
- for (std::size_t i = 0; i < tab_field.nrow(); ++i) {
- beam_former->AddAntenna(
- readAntennaField(tab_field, i, station->get_element_response()));
- }
-
- // TODO
- // If There is only one field, the top level beamformer is not needed
- // and the station antenna can be set the the beamformer of the field
- station->SetAntenna(beam_former);
-
- size_t field_id = 0;
- size_t element_id = 0;
- Antenna::CoordinateSystem coordinate_system =
- common::readCoordinateSystem(tab_field, field_id);
-
- // TODO: rotate coordinate system for antenna
- auto element = Element::Ptr(new Element(
- coordinate_system, station->get_element_response(), element_id));
- station->SetElement(element);
-
+ Station::Ptr station = readStation(ms, id);
return station;
}
\ No newline at end of file
diff --git a/cpp/test/tload_lofar.cc b/cpp/test/tload_lofar.cc
index 6babadd5..d79b589d 100644
--- a/cpp/test/tload_lofar.cc
+++ b/cpp/test/tload_lofar.cc
@@ -4,8 +4,11 @@
#include "./../options.h"
#include "./../gridded_response/lofargrid.h"
#include "./../element_response.h"
+#include "../../external/npy.hpp"
#include "config.h"
+#include <complex>
+#include <cmath>
using namespace everybeam;
@@ -15,7 +18,8 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
casacore::MeasurementSet ms(LOFAR_MOCK_MS);
// Load LOFAR Telescope
- auto telescope = Load(ms, response_model, options);
+ std::unique_ptr<telescope::Telescope> telescope =
+ Load(ms, response_model, options);
// Assert if we indeed have a LOFAR pointer
BOOST_CHECK(nullptr != dynamic_cast<telescope::LOFAR*>(telescope.get()));
@@ -27,14 +31,70 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
// Assert if GetStation(stationd_id) behaves properly
BOOST_CHECK_EQUAL(telescope->GetStation(0)->name(), "CS001HBA0");
- // Get gridded response
- coords::CoordinateSystem coord_system;
- auto grrp = telescope->GetGriddedResponse(coord_system);
+ // Properties extracted from MS
+ double time = 4929192878.008341;
+ double frequency = 138476562.5;
+ std::size_t width(4), height(4);
+ double ra(2.15374123), dec(0.8415521), dl(0.5 * M_PI / 180.),
+ dm(0.5 * M_PI / 180.), shift_l(0.), shift_m(0.);
+
+ coords::CoordinateSystem coord_system = {.width = width,
+ .height = height,
+ .ra = ra,
+ .dec = dec,
+ .dl = dl,
+ .dm = dm,
+ .phase_centre_dl = shift_l,
+ .phase_centre_dm = shift_m};
+ std::unique_ptr<gridded_response::GriddedResponse> grid_response =
+ telescope->GetGriddedResponse(coord_system);
BOOST_CHECK(nullptr !=
- dynamic_cast<gridded_response::LOFARGrid*>(grrp.get()));
+ dynamic_cast<gridded_response::LOFARGrid*>(grid_response.get()));
+
+ // Define buffer and get gridded responses
+ std::vector<std::complex<float>> antenna_buffer_single(
+ grid_response->GetBufferSize(1));
+ grid_response->CalculateStation(antenna_buffer_single.data(), time, frequency,
+ 23);
+ BOOST_CHECK_EQUAL(antenna_buffer_single.size(),
+ std::size_t(width * height * 2 * 2));
+
+ // LOFARBeam output at pixel (2,2):
+ std::vector<std::complex<float>> lofar_p22 = {{-0.175908, -0.000478397},
+ {-0.845988, -0.00121503},
+ {-0.89047, -0.00125383},
+ {0.108123, -5.36076e-05}};
+
+ // Compare with everybeam
+ std::size_t offset_22 = (2 + 2 * height) * 4;
+ for (std::size_t i = 0; i < 4; ++i) {
+ BOOST_CHECK(std::abs(antenna_buffer_single[offset_22 + i] - lofar_p22[i]) <
+ 1e-4);
+ }
+
+ // LOFARBeam output at pixel (1,3):
+ std::vector<std::complex<float>> lofar_p13 = {{-0.158755, -0.000749433},
+ {-0.816165, -0.00272568},
+ {-0.863389, -0.00283979},
+ {0.0936919, 0.000110673}};
+
+ // Compare with everybeam
+ std::size_t offset_13 = (1 + 3 * height) * 4;
+ for (std::size_t i = 0; i < 4; ++i) {
+ BOOST_CHECK(std::abs(antenna_buffer_single[offset_13 + i] - lofar_p13[i]) <
+ 1e-4);
+ }
+
+ std::vector<std::complex<float>> antenna_buffer_all(
+ grid_response->GetBufferSize(telescope->GetNrStations()));
+ grid_response->CalculateAllStations(antenna_buffer_all.data(), time,
+ frequency);
+ BOOST_CHECK_EQUAL(
+ antenna_buffer_all.size(),
+ std::size_t(telescope->GetNrStations() * width * height * 2 * 2));
- // TODO: add test
- // grrp->CalculateStation(0);
- // grrp->CalculateStation();
- // grrp->CalculateAllStations();
+ // Print to np array
+ // const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
+ // npy::SaveArrayAsNumpy("station_responses.npy", false, 4, leshape,
+ // antenna_buffer_single);
}
\ No newline at end of file
diff --git a/demo/comparison-oskar/main.cpp b/demo/comparison-oskar/main.cpp
index e765ccf0..9f9410c1 100644
--- a/demo/comparison-oskar/main.cpp
+++ b/demo/comparison-oskar/main.cpp
@@ -4,7 +4,8 @@
#include <OSKARElementResponse.h>
-#include "npy.hpp" // to save arrays in numpy format
+#include "../../external/npy.hpp"
+// #include "npy.hpp" // to save arrays in numpy format
int main(int argc, char** argv){
// int main() {
diff --git a/demo/comparison-oskar/npy.hpp b/external/npy.hpp
similarity index 100%
rename from demo/comparison-oskar/npy.hpp
rename to external/npy.hpp
--
GitLab