Resolve "Integrate primary beam functionality I - III"

CMake/FindFFTW3.cmake

+# - 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

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)

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

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")

cpp/common/fftresampler.cc

+#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];
+  }
+}

cpp/common/fftresampler.h

+#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_

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,
@@ -36,3 +42,55 @@ void DishGrid::CalculateAllStations(std::complex<float>* buffer, double,
                 buffer + i * width_ * height_ * 4);
   }
 }
+
+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

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

cpp/griddedresponse/griddedresponse.cc

+#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

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

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;

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)
    */

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()

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);

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);
 

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)
 

aocommon @ 100ae51f

-Subproject commit 05917d37e32a1bcd3fc985bc775b0b7b13d12532
+Subproject commit 100ae51f7a13ea4dfc453513970eaa91ab44b987