Merge branch 'vdtol-performance1' into 'master'

Increase computational performance

See merge request !55

cpp/CMakeLists.txt

@@ -13,6 +13,7 @@ add_library(everybeam SHARED
   antenna.cc
   elementresponse.cc
   beamformer.cc
+  beamformeridenticalantennas.cc
   element.cc
   load.cc
   coords/itrfconverter.cc
@@ -44,6 +45,7 @@ install (
 install (FILES
   antenna.h
   beamformer.h
+  beamformeridenticalantennas.h
   element.h
   elementresponse.h
   lofarreadutils.h

cpp/beamformer.cc

@@ -43,7 +43,8 @@ std::vector<std::complex<double>> BeamFormer::ComputeGeometricResponse(
 std::vector<std::pair<std::complex<double>, std::complex<double>>>
 BeamFormer::ComputeWeights(const vector3r_t &pointing, double freq) const {
   // Get geometric response for pointing direction
-  auto geometric_response = ComputeGeometricResponse(freq, pointing);
+  std::vector<std::complex<double>> geometric_response =
+      ComputeGeometricResponse(freq, pointing);
 
   // Initialize and fill result
   double weight_sum[2] = {0.0, 0.0};
@@ -51,7 +52,7 @@ BeamFormer::ComputeWeights(const vector3r_t &pointing, double freq) const {
       geometric_response.size());
   for (std::size_t idx = 0; idx < antennas_.size(); ++idx) {
     // Compute conjugate of geometric response
-    auto phasor_conj = std::conj(geometric_response[idx]);
+    std::complex<double> phasor_conj = std::conj(geometric_response[idx]);
     // Compute the delays in x/y direction
     result[idx] = {phasor_conj * (1.0 * antennas_[idx]->enabled_[0]),
                    phasor_conj * (1.0 * antennas_[idx]->enabled_[1])};
@@ -72,16 +73,20 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
                                       const vector3r_t &direction,
                                       const Options &options) const {
   // Weights based on pointing direction of beam
-  auto weights = ComputeWeights(options.station0, options.freq0);
+  std::vector<std::pair<std::complex<double>, std::complex<double>>> weights =
+      ComputeWeights(options.station0, options.freq0);
   // Weights based on direction of interest
-  auto geometric_response = ComputeGeometricResponse(freq, direction);
+  std::vector<std::complex<double>> geometric_response =
+      ComputeGeometricResponse(freq, direction);
 
   matrix22c_t result = {0};
   for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size();
        ++antenna_idx) {
-    auto antenna = antennas_[antenna_idx];
-    auto antenna_weight = weights[antenna_idx];
-    auto antenna_geometric_reponse = geometric_response[antenna_idx];
+    Antenna::Ptr antenna = antennas_[antenna_idx];
+    std::pair<std::complex<double>, std::complex<double>> antenna_weight =
+        weights[antenna_idx];
+    std::complex<double> antenna_geometric_reponse =
+        geometric_response[antenna_idx];
 
     matrix22c_t antenna_response =
         antenna->Response(time, freq, direction, options);
@@ -96,4 +101,30 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
   }
   return result;
 }
+
+diag22c_t BeamFormer::LocalArrayFactor(real_t time, real_t freq,
+                                       const vector3r_t &direction,
+                                       const Options &options) const {
+  // Weights based on pointing direction of beam
+  std::vector<std::pair<std::complex<double>, std::complex<double>>> weights =
+      ComputeWeights(options.station0, options.freq0);
+  // Weights based on direction of interest
+  std::vector<std::complex<double>> geometric_response =
+      ComputeGeometricResponse(freq, direction);
+
+  diag22c_t result = {0};
+  for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size();
+       ++antenna_idx) {
+    Antenna::Ptr antenna = antennas_[antenna_idx];
+    std::pair<std::complex<double>, std::complex<double>> antenna_weight =
+        weights[antenna_idx];
+    std::complex<double> antenna_geometric_reponse =
+        geometric_response[antenna_idx];
+
+    result[0] += antenna_weight.first * antenna_geometric_reponse;
+    result[1] += antenna_weight.second * antenna_geometric_reponse;
+  }
+  return result;
+}
+
 }  // namespace everybeam

cpp/beamformer.h

@@ -59,7 +59,7 @@ class BeamFormer : public Antenna {
    */
   void AddAntenna(Antenna::Ptr antenna) { antennas_.push_back(antenna); }
 
- private:
+ protected:
   vector3r_t
       local_phase_reference_position_;  // in coordinate system of Antenna
 
@@ -76,9 +76,7 @@ class BeamFormer : public Antenna {
   // "representation" of Jones matrix
   virtual diag22c_t LocalArrayFactor(real_t time, real_t freq,
                                      const vector3r_t &direction,
-                                     const Options &options) const override {
-    return {1.0, 1.0};
-  }
+                                     const Options &options) const override;
 
   // Compute the geometric response for all the antennas in the BeamFormer based
   // on the probing frequency and a specified direction (either pointing dir or

cpp/beamformeridenticalantennas.cc

+#include "beamformeridenticalantennas.h"
+
+#include "common/constants.h"
+#include "common/mathutils.h"
+
+#include <cmath>
+
+namespace everybeam {
+
+matrix22c_t BeamFormerIdenticalAntennas::LocalResponse(
+    real_t time, real_t freq, const vector3r_t &direction,
+    const Options &options) const {
+  matrix22c_t result;
+
+  auto antenna = antennas_[0];
+
+  matrix22c_t antenna_response =
+      antenna->Response(time, freq, direction, options);
+
+  diag22c_t array_factor = LocalArrayFactor(time, freq, direction, options);
+
+  result[0][0] = array_factor[0] * antenna_response[0][0];
+  result[0][1] = array_factor[0] * antenna_response[0][1];
+  result[1][0] = array_factor[1] * antenna_response[1][0];
+  result[1][1] = array_factor[1] * antenna_response[1][1];
+
+  return result;
+}
+
+}  // namespace everybeam

cpp/beamformeridenticalantennas.h

+#ifndef EVERYBEAM_BEAMFORMERIDENTICALANTENNAS_H
+#define EVERYBEAM_BEAMFORMERIDENTICALANTENNAS_H
+
+#include "beamformer.h"
+
+namespace everybeam {
+class BeamFormerIdenticalAntennas : public BeamFormer {
+ public:
+  /**
+   * @brief Construct a new BeamFormerIdenticalAntennas object
+   *
+   */
+  BeamFormerIdenticalAntennas() : BeamFormer() {}
+
+  /**
+   * @brief Construct a new BeamFormerIdenticalAntennas object given a
+   * coordinate system.
+   *
+   * @param coordinate_system
+   */
+  BeamFormerIdenticalAntennas(const CoordinateSystem &coordinate_system)
+      : BeamFormer(coordinate_system) {}
+
+  /**
+   * @brief Construct a new BeamFormer object given a coordinate system and a
+   * phase reference position
+   *
+   * @param coordinate_system
+   * @param phase_reference_position
+   */
+  BeamFormerIdenticalAntennas(CoordinateSystem coordinate_system,
+                              vector3r_t phase_reference_position)
+      : BeamFormer(coordinate_system, phase_reference_position) {}
+
+  BeamFormerIdenticalAntennas(vector3r_t phase_reference_position)
+      : BeamFormer(phase_reference_position) {}
+
+ private:
+  // Compute the BeamFormer response in certain direction of arrival (ITRF, m)
+  // and return (Jones) matrix of response
+  virtual matrix22c_t LocalResponse(real_t time, real_t freq,
+                                    const vector3r_t &direction,
+                                    const Options &options) const override;
+};
+}  // namespace everybeam
+#endif

cpp/lofarreadutils.cc

@@ -22,6 +22,7 @@
 // $Id$
 
 #include "lofarreadutils.h"
+#include "beamformeridenticalantennas.h"
 #include "common/mathutils.h"
 #include "common/casautils.h"
 
@@ -146,7 +147,8 @@ vector3r_t TransformToFieldCoordinates(
 BeamFormer::Ptr MakeTile(unsigned int id, const vector3r_t &position,
                          const TileConfig &tile_config,
                          ElementResponse::Ptr element_response) {
-  BeamFormer::Ptr tile = BeamFormer::Ptr(new BeamFormer(position));
+  BeamFormer::Ptr tile =
+      BeamFormer::Ptr(new BeamFormerIdenticalAntennas(position));
 
   for (unsigned int id = 0; id < tile_config.size(); id++) {
     vector3r_t antenna_position = tile_config[id];
@@ -175,7 +177,8 @@ BeamFormer::Ptr ReadAntennaField(const Table &table, unsigned int id,
   //     coordinate_system.axes.q[2] << std::endl; std::cout << "  axes.r: " <<
   //     coordinate_system.axes.r[0] << ", " << coordinate_system.axes.r[1] <<
   //     ", " << coordinate_system.axes.r[2] << std::endl;
-  BeamFormer::Ptr beam_former(new BeamFormer(coordinate_system));
+  BeamFormer::Ptr beam_former(
+      new BeamFormerIdenticalAntennas(coordinate_system));
 
   ROScalarColumn<String> c_name(table, "NAME");
   ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");