diff --git a/CMake/FindFFTW3.cmake b/CMake/FindFFTW3.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..8b5b4c27e1671695253a0e90c480ee3effb53a53
--- /dev/null
+++ b/CMake/FindFFTW3.cmake
@@ -0,0 +1,95 @@
+# - Try to find FFTW3.
+# Usage: find_package(FFTW3 [COMPONENTS [single double long-double threads]])
+#
+# Variables used by this module:
+# FFTW3_ROOT_DIR - FFTW3 root directory
+# Variables defined by this module:
+# FFTW3_FOUND - system has FFTW3
+# FFTW3_INCLUDE_DIR - the FFTW3 include directory (cached)
+# FFTW3_INCLUDE_DIRS - the FFTW3 include directories
+# (identical to FFTW3_INCLUDE_DIR)
+# FFTW3[FL]?_LIBRARY - the FFTW3 library - double, single(F),
+# long-double(L) precision (cached)
+# FFTW3[FL]?_THREADS_LIBRARY - the threaded FFTW3 library - double, single(F),
+# long-double(L) precision (cached)
+# FFTW3_LIBRARIES - list of all FFTW3 libraries found
+
+# Copyright (C) 2009-2010
+# ASTRON (Netherlands Institute for Radio Astronomy)
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This file is part of the LOFAR software suite.
+# The LOFAR software suite is free software: you can redistribute it and/or
+# modify it under the terms of the GNU General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The LOFAR software suite is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+#
+# $Id: FindFFTW3.cmake 15918 2010-06-25 11:12:42Z loose $
+
+# Use double precision by default.
+if(FFTW3_FIND_COMPONENTS MATCHES "^$")
+ set(_components double)
+else()
+ set(_components ${FFTW3_FIND_COMPONENTS})
+endif()
+
+# Loop over each component.
+set(_libraries)
+foreach(_comp ${_components})
+ if(_comp STREQUAL "single")
+ list(APPEND _libraries fftw3f)
+ elseif(_comp STREQUAL "double")
+ list(APPEND _libraries fftw3)
+ elseif(_comp STREQUAL "long-double")
+ list(APPEND _libraries fftw3l)
+ elseif(_comp STREQUAL "threads")
+ set(_use_threads ON)
+ else(_comp STREQUAL "single")
+ message(FATAL_ERROR "FindFFTW3: unknown component `${_comp}' specified. "
+ "Valid components are `single', `double', `long-double', and `threads'.")
+ endif(_comp STREQUAL "single")
+endforeach(_comp ${_components})
+
+# If using threads, we need to link against threaded libraries as well.
+if(_use_threads)
+ set(_thread_libs)
+ foreach(_lib ${_libraries})
+ list(APPEND _thread_libs ${_lib}_threads)
+ endforeach(_lib ${_libraries})
+ set(_libraries ${_thread_libs} ${_libraries})
+endif(_use_threads)
+
+# Keep a list of variable names that we need to pass on to
+# find_package_handle_standard_args().
+set(_check_list)
+
+# Search for all requested libraries.
+foreach(_lib ${_libraries})
+ string(TOUPPER ${_lib} _LIB)
+ find_library(${_LIB}_LIBRARY ${_lib}
+ HINTS ${FFTW3_ROOT_DIR} PATH_SUFFIXES lib)
+ mark_as_advanced(${_LIB}_LIBRARY)
+ list(APPEND FFTW3_LIBRARIES ${${_LIB}_LIBRARY})
+ list(APPEND _check_list ${_LIB}_LIBRARY)
+endforeach(_lib ${_libraries})
+
+# Search for the header file.
+find_path(FFTW3_INCLUDE_DIR fftw3.h
+ HINTS ${FFTW3_ROOT_DIR} PATH_SUFFIXES include)
+mark_as_advanced(FFTW3_INCLUDE_DIR)
+list(APPEND _check_list FFTW3_INCLUDE_DIR)
+
+set(FFTW3_INCLUDE_DIRS ${FFTW3_INCLUDE_DIR})
+
+# Handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if
+# all listed variables are TRUE
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(FFTW3 DEFAULT_MSG ${_check_list})
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 903af813301dac1b83ffc2fc9b39bfa49799b422..996273fdc7302c7498c9e762ab9b3bf958d0e203 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -53,6 +53,9 @@ find_package(OpenMP REQUIRED)
find_package(Boost REQUIRED)
include_directories(${Boost_INCLUDE_DIR})
+# Find and include FFTW3 float libraries
+find_package(FFTW3 COMPONENTS single REQUIRED)
+
#------------------------------------------------------------------------------
# Set CMake and compiler options
if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
@@ -89,10 +92,6 @@ if("${isSystemDir}" STREQUAL "-1")
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
endif("${isSystemDir}" STREQUAL "-1")
-#------------------------------------------------------------------------------
-# Find/load top level dependencies
-# add_subdirectory(external)
-
#------------------------------------------------------------------------------
# Add source
add_subdirectory(cpp)
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index 182a1ef51fd877de6619cb890c8fff7649e6f04a..9505ade91e2d98ce999e5e7e240e34c97d804ad8 100644
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
@@ -16,6 +16,7 @@ add_library(everybeam SHARED
beamformeridenticalantennas.cc
element.cc
load.cc
+ common/fftresampler.cc
coords/itrfconverter.cc
coords/itrfdirection.cc
lofarreadutils.cc
@@ -25,6 +26,7 @@ add_library(everybeam SHARED
telescope/dish.cc
telescope/mwa.cc
telescope/oskar.cc
+ griddedresponse/griddedresponse.cc
griddedresponse/lofargrid.cc
griddedresponse/dishgrid.cc
griddedresponse/mwagrid.cc
@@ -43,7 +45,7 @@ target_include_directories(everybeam PUBLIC
target_include_directories(everybeam PUBLIC ${CASACORE_INCLUDE_DIR})
target_link_libraries(everybeam PUBLIC hamaker lobes oskar)
-target_link_libraries(everybeam PUBLIC ${CASACORE_LIBRARIES} ${HDF5_LIBRARIES})
+target_link_libraries(everybeam PUBLIC ${CASACORE_LIBRARIES} ${HDF5_LIBRARIES} ${FFTW3F_LIBRARY})
install (
TARGETS everybeam
diff --git a/cpp/common/CMakeLists.txt b/cpp/common/CMakeLists.txt
index 6ba321b5db43d6568326f2d8b4ce5509502ee4a0..d61c842af59936672cf332818564afffd5893bbd 100644
--- a/cpp/common/CMakeLists.txt
+++ b/cpp/common/CMakeLists.txt
@@ -4,4 +4,5 @@ install (FILES
mathutils.h
mutable_ptr.h
types.h
+ fftresampler.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/common")
diff --git a/cpp/common/fftresampler.cc b/cpp/common/fftresampler.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c6b8a82f91cd270d1d876b5201f8a8bab536dad4
--- /dev/null
+++ b/cpp/common/fftresampler.cc
@@ -0,0 +1,171 @@
+#include "fftresampler.h"
+
+#include <complex>
+#include <iostream>
+
+using aocommon::WindowFunction;
+using everybeam::common::FFTResampler;
+
+FFTResampler::FFTResampler(size_t width_in, size_t height_in, size_t width_out,
+ size_t height_out)
+ : width_in_(width_in),
+ height_in_(height_in),
+ width_out_(width_out),
+ height_out_(height_out),
+ fft_width_(std::max(width_in, width_out)),
+ fft_height_(std::max(height_in, height_out)),
+ window_function_(WindowFunction::Rectangular),
+ tukey_inset_size_(0.0),
+ correct_window_(false) {
+ float* input_data = reinterpret_cast<float*>(
+ fftwf_malloc(fft_width_ * fft_height_ * sizeof(float)));
+ fftwf_complex* fft_data = reinterpret_cast<fftwf_complex*>(
+ fftwf_malloc(fft_width_ * fft_height_ * sizeof(fftwf_complex)));
+ in_to_f_plan_ = fftwf_plan_dft_r2c_2d(height_in_, width_in_, input_data,
+ fft_data, FFTW_ESTIMATE);
+ f_to_out_plan_ = fftwf_plan_dft_c2r_2d(height_out_, width_out_, fft_data,
+ input_data, FFTW_ESTIMATE);
+ fftwf_free(fft_data);
+ fftwf_free(input_data);
+}
+
+FFTResampler::~FFTResampler() {
+ fftwf_destroy_plan(in_to_f_plan_);
+ fftwf_destroy_plan(f_to_out_plan_);
+}
+
+void FFTResampler::RunSingle(const Task& task, bool skip_window) const {
+ float* ptr_end = task.input + width_in_ * height_in_;
+ for (float* i = task.input; i != ptr_end; ++i) {
+ if (!std::isfinite(*i)) *i = 0.0;
+ }
+
+ if (window_function_ != WindowFunction::Rectangular && !skip_window)
+ ApplyWindow(task.input);
+
+ size_t fft_in_width = width_in_ / 2 + 1;
+ std::complex<float>* fft_data = reinterpret_cast<std::complex<float>*>(
+ fftwf_malloc(fft_in_width * height_in_ * sizeof(std::complex<float>)));
+ fftwf_execute_dft_r2c(in_to_f_plan_, task.input,
+ reinterpret_cast<fftwf_complex*>(fft_data));
+
+ size_t fft_out_width = width_out_ / 2 + 1;
+ // TODO this can be done without allocating more mem!
+ std::complex<float>* fft_data_new = reinterpret_cast<std::complex<float>*>(
+ fftwf_malloc(fft_out_width * height_out_ * sizeof(std::complex<float>)));
+ std::uninitialized_fill_n(fft_data_new, fft_out_width * height_out_,
+ std::complex<float>(0));
+
+ size_t old_mid_x = width_in_ / 2;
+ size_t new_mid_x = width_out_ / 2;
+
+ size_t min_width = std::min(width_in_, width_out_);
+ size_t min_height = std::min(height_in_, height_out_);
+
+ size_t min_mid_x = min_width / 2;
+ size_t min_mid_y = min_height / 2;
+
+ float factor = 1.0 / (min_width * min_height);
+
+ for (size_t y = 0; y != min_height; ++y) {
+ size_t y_old = y - min_mid_y + height_in_;
+ size_t y_new = y - min_mid_y + height_out_;
+ if (y_old >= height_in_) y_old -= height_in_;
+ if (y_new >= height_out_) y_new -= height_out_;
+
+ // The last dimension is stored half
+ for (size_t x = 0; x != min_mid_x; ++x) {
+ size_t index_old = x + y_old * (old_mid_x + 1);
+ size_t index_new = x + y_new * (new_mid_x + 1);
+
+ fft_data_new[index_new] = fft_data[index_old] * factor;
+ }
+ if (width_in_ >= width_out_) {
+ size_t index_old = width_in_ / 2 + y_old * (old_mid_x + 1);
+ size_t index_new = width_out_ / 2 + y_new * (new_mid_x + 1);
+ fft_data_new[index_new] = fft_data[index_old] * factor;
+ }
+ }
+
+ fftwf_free(fft_data);
+
+ fftwf_execute_dft_c2r(f_to_out_plan_,
+ reinterpret_cast<fftwf_complex*>(fft_data_new),
+ task.output);
+
+ fftwf_free(fft_data_new);
+
+ if (correct_window_ && window_function_ != WindowFunction::Rectangular &&
+ !skip_window)
+ UnapplyWindow(task.output);
+}
+
+void FFTResampler::MakeWindow(aocommon::UVector<float>& data,
+ size_t width) const {
+ if (window_function_ == WindowFunction::Tukey)
+ MakeTukeyWindow(data, width);
+ else {
+ data.resize(width);
+ for (size_t x = 0; x != width; ++x)
+ data[x] = WindowFunction::Evaluate(window_function_, width, x) + 1e-5;
+ }
+}
+
+void FFTResampler::MakeTukeyWindow(aocommon::UVector<float>& data,
+ size_t width) const {
+ // Make a Tukey window, which consists of
+ // left: a cosine going from 0 to 1
+ // mid: all 1
+ // right: a cosine going from 1 to 0
+ data.resize(width);
+ for (size_t x = 0; x != width; ++x) {
+ // left part of Tukey window
+ double x_sh = (0.5 + x) * 2;
+ if (x_sh < width - tukey_inset_size_) {
+ double pos = x_sh / (width - tukey_inset_size_);
+ data[x] = (std::cos((pos + 1.0) * M_PI) + 1.0) * 0.5;
+ } else if (x_sh < width + tukey_inset_size_) {
+ data[x] = 1.0;
+ } else {
+ double pos =
+ (x_sh - (width + tukey_inset_size_)) / (width - tukey_inset_size_);
+ data[x] = (std::cos(pos * M_PI) + 1.0) * 0.5;
+ }
+ }
+}
+
+void FFTResampler::ApplyWindow(float* data) const {
+ if (window_row_in_.empty()) {
+ MakeWindow(window_row_in_, width_in_);
+ MakeWindow(window_col_in_, height_in_);
+ if (correct_window_) {
+ aocommon::UVector<float> window_img_in(width_in_ * height_in_);
+ float* ptr_in = window_img_in.data();
+ for (size_t y = 0; y != height_in_; ++y) {
+ for (size_t x = 0; x != width_in_; ++x) {
+ *ptr_in = window_row_in_[x] * window_col_in_[y];
+ ++ptr_in;
+ }
+ }
+
+ window_out_.resize(width_out_ * height_out_);
+ Task task;
+ task.input = window_img_in.data();
+ task.output = window_out_.data();
+ RunSingle(task, true);
+ }
+ }
+ for (size_t y = 0; y != height_in_; ++y) {
+ for (size_t x = 0; x != width_in_; ++x) {
+ *data *= window_row_in_[x] * window_col_in_[y];
+ ++data;
+ }
+ }
+}
+
+void FFTResampler::UnapplyWindow(float* data) const {
+ size_t n = width_out_ * height_out_;
+ for (size_t i = 0; i != n; ++i) {
+ data[i] /= window_out_[i];
+ }
+}
diff --git a/cpp/common/fftresampler.h b/cpp/common/fftresampler.h
new file mode 100644
index 0000000000000000000000000000000000000000..f6627e7000b8949956d2e5ea771d58619b805f3e
--- /dev/null
+++ b/cpp/common/fftresampler.h
@@ -0,0 +1,95 @@
+#ifndef EVERYBEAM_COMMON_FFT_RESAMPLE_H_
+#define EVERYBEAM_COMMON_FFT_RESAMPLE_H_
+
+#include <aocommon/windowfunction.h>
+#include <aocommon/uvector.h>
+
+#include <vector>
+#include <thread>
+
+#include <fftw3.h>
+
+namespace everybeam {
+namespace common {
+
+/**
+ * @brief FFT resampling from (coarse) input grid to (high resolution) output
+ * grid
+ *
+ */
+class FFTResampler {
+ private:
+ struct Task {
+ float *input, *output;
+ };
+
+ public:
+ /**
+ * @brief Construct a new FFTResampler object
+ *
+ * @param width_in input image width (int)
+ * @param height_in input image height (int)
+ * @param width_out output image width (int)
+ * @param height_out output image height (int)
+ */
+ FFTResampler(size_t width_in, size_t height_in, size_t width_out,
+ size_t height_out);
+
+ ~FFTResampler();
+
+ /**
+ * @brief Do the FFT resampling
+ *
+ * @param input Input image buffer
+ * @param output Output image buffer
+ */
+ void Resample(float* input, float* output) {
+ Task task;
+ task.input = input;
+ task.output = output;
+ RunSingle(task, false);
+ }
+
+ /**
+ * Only to be used with SingleFT (it makes resampling thread unsafe!)
+ */
+ void SetTukeyWindow(double inset_size, bool correct_window) {
+ window_function_ = aocommon::WindowFunction::Tukey;
+ tukey_inset_size_ = inset_size;
+ correct_window_ = correct_window;
+ window_row_in_.clear();
+ window_col_in_.clear();
+ window_out_.clear();
+ }
+
+ void SetWindowFunction(aocommon::WindowFunction::Type window,
+ bool correct_window) {
+ window_function_ = window;
+ correct_window_ = correct_window;
+ window_row_in_.clear();
+ window_col_in_.clear();
+ window_out_.clear();
+ }
+
+ private:
+ void RunSingle(const Task& task, bool skip_window) const;
+ void ApplyWindow(float* data) const;
+ void UnapplyWindow(float* data) const;
+ void MakeWindow(aocommon::UVector<float>& data, size_t width) const;
+ void MakeTukeyWindow(aocommon::UVector<float>& data, size_t width) const;
+
+ size_t width_in_, height_in_;
+ size_t width_out_, height_out_;
+ size_t fft_width_, fft_height_;
+ aocommon::WindowFunction::Type window_function_;
+ double tukey_inset_size_;
+ mutable aocommon::UVector<float> window_row_in_;
+ mutable aocommon::UVector<float> window_col_in_;
+ mutable aocommon::UVector<float> window_out_;
+ bool correct_window_;
+
+ fftwf_plan in_to_f_plan_, f_to_out_plan_;
+};
+} // namespace common
+} // namespace everybeam
+#endif // EVERYBEAM_COMMON_FFT_RESAMPLE_H_
diff --git a/cpp/griddedresponse/dishgrid.cc b/cpp/griddedresponse/dishgrid.cc
index c5bd5a74ff2deb964c0ab67e6e57ebf9e54e114f..ce0f81291bbd04fe79ec3a33d9b9367a137997c5 100644
--- a/cpp/griddedresponse/dishgrid.cc
+++ b/cpp/griddedresponse/dishgrid.cc
@@ -4,9 +4,15 @@
#include "../circularsymmetric/vlabeam.h"
#include <aocommon/uvector.h>
-
+#include <aocommon/matrix2x2.h>
+#include <aocommon/matrix4x4.h>
+#include <aocommon/hmatrix4x4.h>
#include <algorithm>
+using aocommon::HMC4x4;
+using aocommon::MC2x2;
+using aocommon::MC4x4;
+using aocommon::UVector;
using everybeam::griddedresponse::DishGrid;
void DishGrid::CalculateStation(std::complex<float>* buffer, double,
@@ -35,4 +41,56 @@ void DishGrid::CalculateAllStations(std::complex<float>* buffer, double,
std::copy_n(buffer, width_ * height_ * 4,
buffer + i * width_ * height_ * 4);
}
-}
\ No newline at end of file
+}
+
+void DishGrid::CalculateIntegratedResponse(double* buffer, double,
+ double frequency, size_t field_id,
+ size_t undersampling_factor,
+ const std::vector<double>&) {
+ // Copy coordinate members
+ size_t width_original = width_, height_original = height_;
+ double dl_original = dl_, dm_original = dm_;
+
+ width_ /= undersampling_factor;
+ height_ /= undersampling_factor;
+ dl_ *= (double(width_original) / double(width_));
+ dm_ *= (double(width_original) / double(width_));
+
+ // Init (Hermitian) Mueller matrix for every pixel in the coarse grid
+ size_t npixels = width_ * height_;
+ std::vector<HMC4x4> matrices(npixels, HMC4x4::Zero());
+ MakeIntegratedSnapshot(matrices, frequency, field_id);
+
+ DoFFTResampling(buffer, width_, height_, width_original, height_original,
+ matrices);
+
+ // Reset coordinate members to original values
+ width_ = width_original;
+ height_ = height_original;
+ dl_ = dl_original;
+ dm_ = dm_original;
+};
+
+void DishGrid::MakeIntegratedSnapshot(std::vector<aocommon::HMC4x4>& matrices,
+ double frequency, size_t field_id) {
+ size_t nstations = telescope_->GetNrStations();
+ UVector<std::complex<float>> buffer_undersampled(
+ GetStationBufferSize(nstations));
+ // Assert that buffer size can accomodate Jones matrix on pixels
+ assert(buffer_undersampled.size() >= width_ * height_ * 4);
+
+ CalculateAllStations(buffer_undersampled.data(), 0., frequency, field_id);
+ // Loop over the pixels just once, and compute auto correlation
+ for (size_t y = 0; y != height_; ++y) {
+ for (size_t x = 0; x != width_; ++x) {
+ size_t offset = (y * width_ + x) * 4;
+ MC2x2 A(buffer_undersampled[offset], buffer_undersampled[offset + 1],
+ buffer_undersampled[offset + 2], buffer_undersampled[offset + 3]);
+ // Compute Mueller matrix and apply vec trick, see
+ // https://en.wikipedia.org/wiki/Kronecker_product#Matrix_equations
+ // No need to add (+) and average (*0.5) for auto-correlation
+ matrices[y * width_ + x] =
+ HMC4x4::KroneckerProduct(A.HermTranspose().Transpose(), A);
+ }
+ }
+};
\ No newline at end of file
diff --git a/cpp/griddedresponse/dishgrid.h b/cpp/griddedresponse/dishgrid.h
index e0a9bf5a681fbf774572c0824218b86c7140f259..7f2234e59e310326320622217eaa0a8ecc1878a1 100644
--- a/cpp/griddedresponse/dishgrid.h
+++ b/cpp/griddedresponse/dishgrid.h
@@ -41,6 +41,33 @@ class DishGrid final : public GriddedResponse {
void CalculateAllStations(std::complex<float>* buffer, double time,
double frequency, size_t field_id) override;
+
+ virtual void CalculateIntegratedResponse(
+ double* buffer, double time, double frequency, size_t field_id,
+ size_t undersampling_factor,
+ const std::vector<double>& baseline_weights) override;
+
+ virtual void CalculateIntegratedResponse(
+ double* buffer, const std::vector<double>& time_array, double frequency,
+ size_t field_id, size_t undersampling_factor,
+ const std::vector<double>& baseline_weights) override {
+ // Time does not play a role in the integrated response of a dish telescope,
+ // so call CalculateIntegratedResponse as if it were one time step
+ CalculateIntegratedResponse(buffer, 0., frequency, field_id,
+ undersampling_factor, std::vector<double>{0});
+ };
+
+ private:
+ /**
+ * @brief Make integrated snapshot, specialized/simplified for dish
+ * telescopes.
+ *
+ * @param matrices Vector of Mueller matrices
+ * @param frequency Frequency (Hz)
+ * @param field_id Field id
+ */
+ void MakeIntegratedSnapshot(std::vector<aocommon::HMC4x4>& matrices,
+ double frequency, size_t field_id);
};
} // namespace griddedresponse
} // namespace everybeam
diff --git a/cpp/griddedresponse/griddedresponse.cc b/cpp/griddedresponse/griddedresponse.cc
new file mode 100644
index 0000000000000000000000000000000000000000..cbe70aba13fde4b7d740eada10cb7a0a40625ade
--- /dev/null
+++ b/cpp/griddedresponse/griddedresponse.cc
@@ -0,0 +1,170 @@
+#include "griddedresponse.h"
+#include "../common/fftresampler.h"
+
+#include "../telescope/telescope.h"
+
+#include <aocommon/matrix2x2.h>
+#include <aocommon/matrix4x4.h>
+#include <vector>
+#include <complex>
+#include <stdexcept>
+
+using aocommon::HMC4x4;
+using aocommon::MC2x2;
+using aocommon::MC4x4;
+using aocommon::UVector;
+using everybeam::common::FFTResampler;
+using everybeam::griddedresponse::GriddedResponse;
+
+void GriddedResponse::CalculateIntegratedResponse(
+ double* buffer, double time, double frequency, size_t field_id,
+ size_t undersampling_factor, const std::vector<double>& baseline_weights) {
+ size_t nstations = telescope_->GetNrStations();
+ size_t nbaselines = nstations * (nstations + 1) / 2;
+ if (baseline_weights.size() != nbaselines) {
+ throw std::runtime_error("baseline_weights vector has incorrect size.");
+ }
+ // Scaling factor
+ double baseline_total_weight =
+ std::accumulate(baseline_weights.begin(), baseline_weights.end(), 0);
+
+ // Copy coordinate members
+ size_t width_original = width_, height_original = height_;
+ double dl_original = dl_, dm_original = dm_;
+
+ width_ /= undersampling_factor;
+ height_ /= undersampling_factor;
+ dl_ *= (double(width_original) / double(width_));
+ dm_ *= (double(width_original) / double(width_));
+
+ // Init (Hermitian) Mueller matrix for every pixel in the coarse grid
+ size_t npixels = width_ * height_;
+ std::vector<HMC4x4> matrices(npixels, HMC4x4::Zero());
+ MakeIntegratedSnapshot(matrices, time, frequency, field_id,
+ baseline_weights.data());
+
+ for (HMC4x4& matrix : matrices) {
+ matrix /= baseline_total_weight;
+ }
+
+ DoFFTResampling(buffer, width_, height_, width_original, height_original,
+ matrices);
+
+ // Reset coordinate members to original values
+ width_ = width_original;
+ height_ = height_original;
+ dl_ = dl_original;
+ dm_ = dm_original;
+}
+
+void GriddedResponse::CalculateIntegratedResponse(
+ double* buffer, const std::vector<double>& time_array, double frequency,
+ size_t field_id, size_t undersampling_factor,
+ const std::vector<double>& baseline_weights) {
+ size_t nstations = telescope_->GetNrStations();
+ size_t nbaselines = nstations * (nstations + 1) / 2;
+ if (baseline_weights.size() != time_array.size() * nbaselines) {
+ throw std::runtime_error("baseline_weights vector has incorrect size.");
+ }
+
+ // Scaling factor
+ double baseline_total_weight =
+ std::accumulate(baseline_weights.begin(), baseline_weights.end(), 0);
+
+ // Copy coordinate members
+ size_t width_original = width_, height_original = height_;
+ double dl_original = dl_, dm_original = dm_;
+
+ width_ /= undersampling_factor;
+ height_ /= undersampling_factor;
+ dl_ *= (double(width_original) / double(width_));
+ dm_ *= (double(width_original) / double(width_));
+
+ // Init (Hermitian) Mueller matrix for every pixel in the coarse grid
+ size_t npixels = width_ * height_;
+ std::vector<HMC4x4> matrices(npixels, HMC4x4::Zero());
+
+ for (std::size_t tstep = 0; tstep != time_array.size(); ++tstep) {
+ MakeIntegratedSnapshot(matrices, time_array[tstep], frequency, field_id,
+ baseline_weights.data() + tstep * nbaselines);
+ }
+
+ for (HMC4x4& matrix : matrices) {
+ matrix /= baseline_total_weight;
+ }
+
+ DoFFTResampling(buffer, width_, height_, width_original, height_original,
+ matrices);
+
+ // Reset coordinate members to original values
+ width_ = width_original;
+ height_ = height_original;
+ dl_ = dl_original;
+ dm_ = dm_original;
+}
+
+void GriddedResponse::MakeIntegratedSnapshot(
+ std::vector<HMC4x4>& matrices, double time, double frequency,
+ size_t field_id, const double* baseline_weights_interval) {
+ size_t nstations = telescope_->GetNrStations();
+ UVector<std::complex<float>> buffer_undersampled(
+ GetStationBufferSize(nstations));
+ CalculateAllStations(buffer_undersampled.data(), time, frequency, field_id);
+
+ size_t npixels = width_ * height_;
+ for (size_t y = 0; y != height_; ++y) {
+ for (size_t x = 0; x != width_; ++x) {
+ size_t index = 0;
+ HMC4x4 gain = HMC4x4::Zero();
+ for (size_t s1 = 0; s1 != nstations; ++s1) {
+ size_t offset_s1 = (s1 * npixels + y * width_ + x) * 4;
+ MC2x2 A(buffer_undersampled[offset_s1],
+ buffer_undersampled[offset_s1 + 1],
+ buffer_undersampled[offset_s1 + 2],
+ buffer_undersampled[offset_s1 + 3]);
+ for (size_t s2 = s1; s2 != nstations; ++s2) {
+ size_t offset_s2 = offset_s1 + (s2 - s1) * npixels * 4;
+ MC2x2 B(buffer_undersampled[offset_s2],
+ buffer_undersampled[offset_s2 + 1],
+ buffer_undersampled[offset_s2 + 2],
+ buffer_undersampled[offset_s2 + 3]);
+
+ // Compute Mueller matrix and apply vec trick, see
+ // https://en.wikipedia.org/wiki/Kronecker_product#Matrix_equations
+ MC4x4 baseline_mueller =
+ MC4x4::KroneckerProduct(B.HermTranspose().Transpose(), A) +
+ MC4x4::KroneckerProduct(A.HermTranspose().Transpose(), B);
+
+ // Note: keep scalar factor in brackets to avoid a redundant
+ // element-wise multiplication of the HMC4x4 matrix. Factor 0.5 is to
+ // account for Kronecker delta additions.
+ gain += HMC4x4(baseline_mueller) *
+ (0.5 * baseline_weights_interval[index]);
+ ++index;
+ }
+ }
+ matrices[y * width_ + x] += gain;
+ }
+ }
+}
+
+void GriddedResponse::DoFFTResampling(
+ double* buffer, int width_in, int height_in, int width_out, int height_out,
+ const std::vector<aocommon::HMC4x4>& matrices) {
+ // (FFT) resampling, run multi-threaded?
+ FFTResampler resampler(width_in, height_in, width_out, height_out);
+ UVector<float> lowres_input(width_in * height_in);
+ // Loop over the "double" representation of the HMC4x4 Hermitian matrix
+ for (size_t p = 0; p != 16; ++p) {
+ for (int i = 0; i != width_in * height_in; ++i) {
+ lowres_input[i] = matrices[i].Data(p);
+ }
+ // Resample and write to the output buffer
+ UVector<float> hires_output(width_out * height_out);
+ resampler.Resample(lowres_input.data(), hires_output.data());
+ // Copy hires_output (float) to buffer (double), conversion is implicit in
+ // std::copy
+ std::copy(hires_output.begin(), hires_output.end(),
+ buffer + p * width_out * height_out);
+ }
+}
\ No newline at end of file
diff --git a/cpp/griddedresponse/griddedresponse.h b/cpp/griddedresponse/griddedresponse.h
index ce070ec3bb57ffdb748817ca265a957fade34ed9..f1dc5068638fbe95dfbcd975f213dfd69b5c6861 100644
--- a/cpp/griddedresponse/griddedresponse.h
+++ b/cpp/griddedresponse/griddedresponse.h
@@ -31,7 +31,11 @@
#include <memory>
#include <vector>
#include <thread>
+
#include <aocommon/lane.h>
+#include <aocommon/uvector.h>
+#include <aocommon/hmatrix4x4.h>
+
#include <casacore/measures/Measures/MDirection.h>
namespace everybeam {
@@ -52,7 +56,7 @@ class GriddedResponse {
* @brief Compute the Beam response for a single station
*
* @param buffer Output buffer, compute and set size with
- * GriddedResponse::GetBufferSize(1)
+ * GriddedResponse::GetStationBufferSize(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).
@@ -66,15 +70,54 @@ class GriddedResponse {
* @brief Compute the Beam response for all stations in a Telescope
*
* @param buffer Output buffer, compute and set size with
- * GriddedResponse::GetBufferSize()
+ * GriddedResponse::GetStationBufferSize()
* @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
* @param frequency Frequency (Hz)
*/
virtual void CalculateAllStations(std::complex<float>* buffer, double time,
double frequency, size_t field_id) = 0;
- std::size_t GetBufferSize(std::size_t nstations) const {
- return std::size_t(nstations * width_ * height_ * 2 * 2);
+ /**
+ * @brief Calculate integrated beam for a single time step.
+ *
+ * @param buffer Buffer for storing the result, should have size width *
+ * height * 16
+ * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * @param frequency Frequency (Hz)
+ * @param field_id Field id
+ * @param undersampling_factor Undersampling factor
+ * @param baseline_weights Baseline weights, size should equal
+ * Telescope::GetNrStations() * (Telescope::GetNrStations() + 1)/2
+ */
+ virtual void CalculateIntegratedResponse(
+ double* buffer, double time, double frequency, size_t field_id,
+ size_t undersampling_factor, const std::vector<double>& baseline_weights);
+
+ /**
+ * @brief Calculate integrated beam over multiple time steps.
+ *
+ * @param buffer Buffer for storing the result, should have size width *
+ * height * 16
+ * @param time_array Vector with probing times, modified Julian date, UTC, in
+ * seconds (MJD(UTC), s).
+ * @param frequency Frequency (Hz)
+ * @param field_id Field id
+ * @param undersampling_factor Undersampling factor
+ * @param baseline_weights Baseline weights, size should equal
+ * (Telescope::GetNrStations() * (Telescope::GetNrStations() + 1)/2) *
+ * time_array.size()
+ */
+ virtual void CalculateIntegratedResponse(
+ double* buffer, const std::vector<double>& time_array, double frequency,
+ size_t field_id, size_t undersampling_factor,
+ const std::vector<double>& baseline_weights);
+
+ std::size_t GetStationBufferSize(std::size_t nstations) const {
+ return std::size_t(nstations * width_ * height_ * 4);
+ }
+
+ std::size_t GetIntegratedBufferSize() const {
+ return std::size_t(width_ * height_ * 16);
}
protected:
@@ -96,9 +139,18 @@ class GriddedResponse {
phase_centre_dl_(coordinate_system.phase_centre_dl),
phase_centre_dm_(coordinate_system.phase_centre_dm){};
+ static void DoFFTResampling(double* buffer, int width_in, int height_in,
+ int width_out, int height_out,
+ const std::vector<aocommon::HMC4x4>& matrices);
+
const telescope::Telescope* const telescope_;
size_t width_, height_;
double ra_, dec_, dl_, dm_, phase_centre_dl_, phase_centre_dm_;
+
+ private:
+ void MakeIntegratedSnapshot(std::vector<aocommon::HMC4x4>& matrices,
+ double time, double frequency, size_t field_id,
+ const double* baseline_weights_interval);
};
} // namespace griddedresponse
} // namespace everybeam
diff --git a/cpp/griddedresponse/lofargrid.cc b/cpp/griddedresponse/lofargrid.cc
index 0afee4becd9568d49b7a8c8ca3f9100c44f389e7..27dbd6078bad04dba9a84ca91c00d99bac372b66 100644
--- a/cpp/griddedresponse/lofargrid.cc
+++ b/cpp/griddedresponse/lofargrid.cc
@@ -160,7 +160,7 @@ void LOFARGrid::CalcThread(std::complex<float>* buffer, double time,
itrf_direction[2] =
l * l_vector_itrf_[2] + m * m_vector_itrf_[2] + n * n_vector_itrf_[2];
- std::complex<float>* base_buffer = buffer + (x + job.y * height_) * 4;
+ std::complex<float>* base_buffer = buffer + (x + job.y * width_) * 4;
std::complex<float>* ant_buffer_ptr =
base_buffer + job.buffer_offset * values_per_ant;
diff --git a/cpp/griddedresponse/lofargrid.h b/cpp/griddedresponse/lofargrid.h
index 28ada8044524046b0312a0a885aa4aabecabcd65..40d41de119d8100cc9a11fea7936039e7705d3af 100644
--- a/cpp/griddedresponse/lofargrid.h
+++ b/cpp/griddedresponse/lofargrid.h
@@ -54,7 +54,7 @@ class LOFARGrid final : public GriddedResponse {
* @brief Compute the Beam response for a single station
*
* @param buffer Output buffer, compute and set size with
- * GriddedResponse::GetBufferSize(1)
+ * GriddedResponse::GetStationBufferSize(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).
@@ -68,7 +68,7 @@ class LOFARGrid final : public GriddedResponse {
* @brief Compute the Beam response for all stations in a Telescope
*
* @param buffer Output buffer, compute and set size with
- * GriddedResponse::GetBufferSize()
+ * GriddedResponse::GetStationBufferSize()
* @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
* @param freq Frequency (Hz)
*/
diff --git a/cpp/test/tlofar.cc b/cpp/test/tlofar.cc
index 1f0433eb39659ff463bd05ca25b518743e0f3a88..87f769416f6c07eb8da7208bfc086a7a7a896607 100644
--- a/cpp/test/tlofar.cc
+++ b/cpp/test/tlofar.cc
@@ -9,6 +9,7 @@
#include "config.h"
#include <complex>
#include <cmath>
+#include <iostream>
using everybeam::ElementResponseModel;
using everybeam::Load;
@@ -62,7 +63,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
// Define buffer and get gridded responses
std::vector<std::complex<float>> antenna_buffer_single(
- grid_response->GetBufferSize(1));
+ grid_response->GetStationBufferSize(1));
grid_response->CalculateStation(antenna_buffer_single.data(), time, frequency,
23, 0);
BOOST_CHECK_EQUAL(antenna_buffer_single.size(),
@@ -94,8 +95,9 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
1e-4);
}
+ // All stations
std::vector<std::complex<float>> antenna_buffer_all(
- grid_response->GetBufferSize(telescope->GetNrStations()));
+ grid_response->GetStationBufferSize(telescope->GetNrStations()));
grid_response->CalculateAllStations(antenna_buffer_all.data(), time,
frequency, 0);
BOOST_CHECK_EQUAL(
@@ -114,7 +116,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
telescope_diff_beam->GetGriddedResponse(coord_system);
std::vector<std::complex<float>> antenna_buffer_diff_beam(
- grid_response_diff_beam->GetBufferSize(1));
+ grid_response_diff_beam->GetStationBufferSize(1));
grid_response_diff_beam->CalculateStation(antenna_buffer_diff_beam.data(),
time, frequency, 15, 0);
@@ -124,10 +126,87 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
}
BOOST_CHECK(std::abs(norm_jones_mat - 2.) < 1e-6);
- // Print to np array
- const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
- npy::SaveArrayAsNumpy("lofar_station_responses.npy", false, 4, leshape,
- antenna_buffer_single);
+ // Primary beam tests
+ //
+ // Just check whether CalculateIntegratedResponse does run and reproduces
+ // results One time interval
+ std::vector<double> antenna_buffer_integrated(
+ grid_response->GetIntegratedBufferSize());
+ std::vector<double> baseline_weights(
+ telescope->GetNrStations() * (telescope->GetNrStations() + 1) / 2, 1.);
+ grid_response->CalculateIntegratedResponse(antenna_buffer_integrated.data(),
+ time, frequency, 0, 2,
+ baseline_weights);
+
+ // Just check whether some (rather arbitrary) numbers are reproduced
+ BOOST_CHECK(std::abs(antenna_buffer_integrated[10] - 0.0262708) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer_integrated[20] - 0.127972) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer_integrated.back() - 0.00847742) < 1e-6);
+
+ // Two time intervals, should give same output as single time interval
+ std::fill(antenna_buffer_integrated.begin(), antenna_buffer_integrated.end(),
+ 0);
+ std::vector<double> tarray = {time, time};
+ baseline_weights.resize(baseline_weights.size() * tarray.size());
+ std::fill(baseline_weights.begin(), baseline_weights.end(), 1.);
+ grid_response->CalculateIntegratedResponse(antenna_buffer_integrated.data(),
+ tarray, frequency, 0, 2,
+ baseline_weights);
+
+ BOOST_CHECK(std::abs(antenna_buffer_integrated[10] - 0.0262708) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer_integrated[20] - 0.127972) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer_integrated.back() - 0.00847742) < 1e-6);
+
+ // Primary beam response on 40 x 40 grid.
+ // Validated results were obtained with the following wsclean command
+ //
+ // wsclean -size 40 40 -scale 900asec -apply-primary-beam LOFAR_MOCK.ms
+ //
+ // where LOFAR_MOCK.ms the MS available from
+ // https://www.astron.nl/citt/EveryBeam/L258627-one-timestep.tar.bz2
+ //
+ // PLEASE NOTE: for the sake of testing, the baseline weights were set to 1
+ // in wsclean::lbeamimagemaker, i.e.
+ //
+ // --- a/lofar/lbeamimagemaker.cpp
+ // +++ b/lofar/lbeamimagemaker.cpp
+ // @@ -380,7 +380,7 @@ void LBeamImageMaker::makeBeamSnapshot(
+ // MC4x4::KroneckerProduct(
+ // stationGains[a1].HermTranspose().Transpose(),
+ // stationGains[a2]);
+ // - double w = weights.Value(a1, a2);
+ // + double w = 1.;
+ //
+ std::size_t width_pb = 40, height_pb = 40;
+ CoordinateSystem coord_system_pb = {.width = width_pb,
+ .height = height_pb,
+ .ra = ra,
+ .dec = dec,
+ // 900asec
+ .dl = (0.25 * M_PI / 180.),
+ .dm = (0.25 * M_PI / 180.),
+ .phase_centre_dl = shift_l,
+ .phase_centre_dm = shift_m};
+ std::unique_ptr<GriddedResponse> grid_response_pb =
+ telescope->GetGriddedResponse(coord_system_pb);
+
+ std::vector<double> antenna_buffer_pb(
+ grid_response_pb->GetIntegratedBufferSize());
+ std::vector<double> baseline_weights_pb(
+ telescope->GetNrStations() * (telescope->GetNrStations() + 1) / 2, 1.);
+ grid_response_pb->CalculateIntegratedResponse(
+ antenna_buffer_pb.data(), time, frequency, 0, 8, baseline_weights_pb);
+ // Check diagonal and off-diagonal term in component 0 and 5 of HMC4x4
+ // representation of Mueller matrix
+ std::size_t offset_01616 = 16 * width_pb + 16,
+ offset_02310 = 23 * width_pb + 10,
+ offset_52020 = 5 * width_pb * height_pb + 20 * width_pb + 20,
+ offset_51825 = 5 * width_pb * height_pb + 18 * width_pb + 25;
+
+ BOOST_CHECK(std::abs(antenna_buffer_pb[offset_01616] - 0.020324793) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer_pb[offset_02310] - 0.0059926948) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer_pb[offset_52020] - 0.00018088287) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer_pb[offset_51825] - 0.00013052078) < 1e-6);
}
BOOST_AUTO_TEST_SUITE_END()
diff --git a/cpp/test/tmwa.cc b/cpp/test/tmwa.cc
index 82840e1798d2a66ef921df83c43e87209700e70d..5a8fa5ce2b1edc25858a72dd2f06aa0ff9d73256 100644
--- a/cpp/test/tmwa.cc
+++ b/cpp/test/tmwa.cc
@@ -54,7 +54,7 @@ BOOST_AUTO_TEST_CASE(load_mwa) {
BOOST_CHECK(nullptr != dynamic_cast<MWAGrid*>(grid_response.get()));
std::vector<std::complex<float>> antenna_buffer(
- grid_response->GetBufferSize(telescope->GetNrStations()));
+ grid_response->GetStationBufferSize(telescope->GetNrStations()));
grid_response->CalculateAllStations(antenna_buffer.data(), time, frequency,
0);
diff --git a/cpp/test/tvla.cc b/cpp/test/tvla.cc
index 76364a394d16161a7f795f9c0f935f3d796a9ae1..903c6374f14e9a80c0e146d9616b752bda13cc37 100644
--- a/cpp/test/tvla.cc
+++ b/cpp/test/tvla.cc
@@ -53,7 +53,7 @@ BOOST_AUTO_TEST_CASE(load_vla) {
BOOST_CHECK(nullptr != dynamic_cast<DishGrid*>(grid_response.get()));
std::vector<std::complex<float>> antenna_buffer(
- grid_response->GetBufferSize(telescope->GetNrStations()));
+ grid_response->GetStationBufferSize(telescope->GetNrStations()));
grid_response->CalculateAllStations(antenna_buffer.data(), time, frequency,
0);
diff --git a/external/README.md b/external/README.md
index 18b169d710a50d55c412eccf60083686ccf4f981..bedab49b789c7059603473b5704cf4f1905f3287 100644
--- a/external/README.md
+++ b/external/README.md
@@ -3,7 +3,7 @@ External dependencies are included via git submodules, see .gitmodules in the ro
Please note that the dependencies are fixed to a specific commit, to make sure tests do not break as a
result of silently updating the submodules. Dependencies are checked out on the following commits:
-- `aocommon`: 05917d37e32a1bcd3fc985bc775b0b7b13d12532 (https://gitlab.com/aroffringa/aocommon)
+- `aocommon`: dd3722b9c3c3b8f0eb14f75db4b0bfb80cc91a68 (https://gitlab.com/aroffringa/aocommon)
- `eigen`: `3.3.7` (https://gitlab.com/libeigen/eigen/-/tags/3.3.7)
- `pybind11`: `v2.5.0` (https://github.com/pybind/pybind11/releases/tag/v2.5.0)
diff --git a/external/aocommon b/external/aocommon
index 05917d37e32a1bcd3fc985bc775b0b7b13d12532..100ae51f7a13ea4dfc453513970eaa91ab44b987 160000
--- a/external/aocommon
+++ b/external/aocommon
@@ -1 +1 @@
-Subproject commit 05917d37e32a1bcd3fc985bc775b0b7b13d12532
+Subproject commit 100ae51f7a13ea4dfc453513970eaa91ab44b987