Ira development: add/generate random noise to data

CMakeLists.txt

@@ -155,6 +155,7 @@ add_library(IDGPredict OBJECT
 set(IDGPREDICT_OBJECT $<TARGET_OBJECTS:IDGPredict>)
 
 add_library(DPPP_OBJ OBJECT
+  DPPP/AddNoise.cc
   DPPP/Apply.cc
   DPPP/ApplyCal.cc
   DPPP/Averager.cc

DPPP/AddNoise.cc

@@ -0,0 +1,222 @@
+// AddNoise.cc: DPPP step class to add HBA/LBA random noise to data
+// Copyright (C) 2013
+// 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: GainCal.cc 21598 2012-07-16 08:07:34Z diepen $
+//
+// @author Claudio Gheller, Henrik Edler
+
+#include "AddNoise.h"
+
+#include <iostream>
+
+#include "../Common/ParameterSet.h"
+#include "../Common/Timer.h"
+
+#include <stddef.h>
+#include <random>
+#include <sstream>
+#include <string>
+#include <utility>
+#include <vector>
+#include "Exceptions.h"
+
+using namespace casacore;
+
+namespace DP3 {
+namespace DPPP {
+
+AddNoise::AddNoise(DPInput* input, const ParameterSet& parset,
+                   const string& prefix)
+    : itsInput(input) {
+  mode = parset.getInt(prefix + "mode", 0);
+  factor = parset.getDouble(prefix + "factor", 1.0);
+}
+
+AddNoise::~AddNoise() {}
+
+void AddNoise::updateInfo(const DPInfo& infoIn) {
+  info() = infoIn;
+  info().setNeedVisData();
+  info().setWriteData();
+}
+
+void AddNoise::show(std::ostream& os) const {
+  os << "AddNoise " << itsName << endl;
+}
+
+void AddNoise::showTimings(std::ostream& os, double duration) const {
+  os << "  ";
+  FlagCounter::showPerc1(os, itsTimer.getElapsed(), duration);
+  os << " AddNoise " << itsName << endl;
+}
+
+bool AddNoise::process(const DPBuffer& buf) {
+  itsTimer.start();
+
+  // Name of the column to add the noise (at the moment not used, just a
+  // placeholder)
+  string column = "DATA";
+  DPBuffer itsBuf;
+  itsBuf.copy(buf);
+
+  Array<Complex>::contiter indIter = itsBuf.getData().cbegin();
+
+  // Set the exposure
+  double exposure = buf.getExposure();
+
+  // Load the Antenna columns
+  Vector<Int> antenna1 = info().getAnt1();
+  Vector<Int> antenna2 = info().getAnt2();
+
+  // Load the Antenna type
+  antennaSet = info().antennaSet();
+
+  // Set Number of baselines
+  int n_baselines = antenna1.size();
+
+  // Set the number of correlations
+  int n_corr = info().ncorr();
+
+  // Set the Antenna names
+  const Vector<String>& allNames = getInfo().antennaNames();
+
+  // Set the Polynomial coefficients
+  double* coeff;
+  if (antennaSet.compare("LBA") == 0 || antennaSet.compare("LBA_ALL") == 0)
+    coeff = coeffs_all_LBA;
+  if (antennaSet.compare("LBA_INNER") == 0) coeff = coeffs_inner_LBA;
+  if (antennaSet.compare("LBA_OUTER") == 0) coeff = coeffs_outer_LBA;
+
+  // Set the number of frequency channels
+  int n_freq = getInfo().chanFreqs().size();
+  Vector<Double> chan_freq = getInfo().chanFreqs();
+  Vector<Double> chan_width = getInfo().chanWidths();
+
+  unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+  std::default_random_engine generator(seed);
+
+  DPBuffer runBuf1;
+  runBuf1.copy(buf);
+  Array<Complex>::contiter run1Iter = runBuf1.getData().cbegin();
+  DPBuffer runBuf2;
+  runBuf2.copy(buf);
+  Array<Complex>::contiter run2Iter = runBuf2.getData().cbegin();
+
+  // Add noise
+
+  int ibegin = 0;
+  int iend = n_baselines;
+  long icount = 0;
+  double nu;
+  double stddev;
+  double sefd;
+  double sefd1;
+  double sefd2;
+
+  //if (antennaSet.compare("LBA") == 0 || antennaSet.compare("LBA_INNER") == 0 || antennaSet.compare("LBA_OUTER") == 0 ||
+  //    antennaSet.compare("LBA_ALL") == 0) {
+  if (antennaSet.substr(0, 3) == "LBA") {
+    for (int ibase = ibegin; ibase < iend; ibase++) {
+      for (int ifreq = 0; ifreq < n_freq; ifreq++) {
+        nu = chan_freq(ifreq);
+        sefd = coeff[0] + coeff[1] * nu + coeff[2] * pow(nu, 2.0) +
+               coeff[3] * pow(nu, 3.0) + coeff[4] * pow(nu, 4.0) +
+               coeff[5] * pow(nu, 5.0);
+        stddev = factor * sefd / eta;
+        stddev = stddev / sqrt(2.0 * exposure * chan_width[ifreq]);
+        std::normal_distribution<double> distribution(0.0, stddev);
+
+        for (int icorr = 0; icorr < n_corr; icorr++) {
+          double noise_real = distribution(generator);
+          double noise_img = distribution(generator);
+          std::complex<float> c_noise((float)noise_real, (float)noise_img);
+          *run1Iter = c_noise;
+          *run2Iter = *indIter + c_noise;
+          run1Iter++;
+          run2Iter++;
+          indIter++;
+          icount++;
+        }
+      }
+    }
+  } else if (antennaSet.substr(0, 3) == "HBA") {
+    double* coeff1;
+    double* coeff2;
+    for (int ibase = ibegin; ibase < iend; ibase++) {
+      if (allNames[antenna1[ibase]].substr(0, 2) == "RS") {
+        coeff1 = coeffs_rs_HBA;
+      } else {
+        coeff1 = coeffs_cs_HBA;
+      }
+      if (allNames[antenna2[ibase]].substr(0, 2) == "RS") {
+        coeff2 = coeffs_rs_HBA;
+      } else {
+        coeff2 = coeffs_cs_HBA;
+      }
+      for (int ifreq = 0; ifreq < n_freq; ifreq++) {
+        nu = chan_freq(ifreq);
+        sefd1 = coeff1[0] + coeff1[1] * nu + coeff1[2] * pow(nu, 2.0) +
+                coeff1[3] * pow(nu, 3.0) + coeff1[4] * pow(nu, 4.0) +
+                coeff1[5] * pow(nu, 5.0);
+        sefd2 = coeff2[0] + coeff2[1] * nu + coeff2[2] * pow(nu, 2.0) +
+                coeff2[3] * pow(nu, 3.0) + coeff2[4] * pow(nu, 4.0) +
+                coeff2[5] * pow(nu, 5.0);
+        sefd = sqrt(sefd1 * sefd2);
+        stddev = factor * sefd / eta;
+        stddev = stddev / sqrt(2.0 * exposure * chan_width[ifreq]);
+        std::normal_distribution<double> distribution(0.0, stddev);
+
+        for (int icorr = 0; icorr < n_corr; icorr++) {
+          double noise_real = distribution(generator);
+          double noise_img = distribution(generator);
+          std::complex<float> c_noise((float)noise_real, (float)noise_img);
+          *run1Iter = c_noise;
+          *run2Iter = *indIter + c_noise;
+          run1Iter++;
+          run2Iter++;
+          indIter++;
+          icount++;
+        }
+      }
+    }
+  } else {
+    throw Exception(
+        "Antenna should be HBA, LBA, LBA_INNER, LBA_OUTER, LBA_ALL");
+  }
+
+  if (mode == 0) {
+    itsBuf.setData(runBuf1.getData());
+  } else if (mode == 1) {
+    itsBuf.setData(runBuf2.getData());
+  } else {
+    cout << "Mode not supported" << endl;
+    exit(100);
+  }
+
+  itsTimer.stop();
+  getNextStep()->process(itsBuf);
+  return false;
+}
+
+void AddNoise::finish() {
+  // Let the next steps finish.
+  getNextStep()->finish();
+}
+}  // namespace DPPP
+}  // namespace DP3

DPPP/AddNoise.h

@@ -0,0 +1,105 @@
+// AddNoise.h: DPPP step class to add HBA/LBA random noise to data
+// Copyright (C) 2013
+// 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/>.
+
+/// @file
+/// @brief DPPP step class to add HBA/LBA random noise to data
+/// @author Tammo Jan Dijkema
+//
+
+#ifndef DPPP_AddNoise_H
+#define DPPP_AddNoise_H
+
+#include "DPBuffer.h"
+#include "DPInput.h"
+
+#include <utility>
+
+#define POL_DEGREE 5
+
+namespace DP3 {
+
+class ParameterSet;
+
+namespace DPPP {
+/// @brief DPPP step class to AddNoise visibilities from a source model
+
+class AddNoise : public DPStep {
+ public:
+  /// Construct the object.
+  /// Parameters are obtained from the parset using the given prefix.
+  AddNoise(DPInput*, const ParameterSet&, const string& prefix);
+
+  virtual ~AddNoise();
+
+  /// Process the data calculating or adding to DATA the estimated random
+  /// noise (van Haarlem et al. 2013). When processed, it invokes the process
+  /// function of the next step.
+  virtual bool process(const DPBuffer&);
+
+  /// Finish the processing of this step and subsequent steps.
+  virtual void finish();
+
+  /// Update the general info.
+  virtual void updateInfo(const DPInfo&);
+
+  /// Show the step parameters.
+  virtual void show(std::ostream&) const;
+
+  /// Show the timings.
+  virtual void showTimings(std::ostream&, double duration) const;
+
+ private:
+  DPInput* itsInput;
+  string itsName;
+  DPBuffer itsBuffer;
+
+  NSTimer itsTimer;
+  /// The mode parameter select the kind of output from the function
+  /// SET mode = 0 data are modified with the noise > data = data + noise
+  /// ADD mode = 1 a new array is created > newdata = data + noise
+  int mode;
+  /// The antennaSet parameter selects between: LBA=LBA_ALL, LBA_INNER
+  ///                                        LBA_OUTER and HBA (read from MS)
+  string antennaSet;
+  /// Coefficients for polynomial interpolation (from constant -first- to highet
+  /// order -last-)
+  double coeffs_all_LBA[POL_DEGREE + 1] = {
+      4.194759691372669e+05,  -0.040624470842582,    1.657038099833047e-09,
+      -3.318591685264548e-17, 3.220530883859981e-25, -1.199723767939448e-33};
+  double coeffs_inner_LBA[POL_DEGREE + 1] = {
+      6.468156199342838e+05,  -0.065296541139271,    2.752773309538937e-09,
+      -5.659747549065881e-17, 5.612743945799180e-25, -2.133417264334753e-33};
+  double coeffs_outer_LBA[POL_DEGREE + 1] = {
+      4.716452746313004e+05,  -0.042301679603444,    1.632924288392071e-09,
+      -3.129891572910005e-17, 2.925471926740372e-25, -1.048279929083482e-33};
+  double coeffs_cs_HBA[POL_DEGREE + 1] = {
+      1.403173283860732e+06,  -0.044309811184301,    5.564568767538230e-10,
+      -3.461658816150442e-18, 1.064770242682354e-26, -1.292413137320037e-35};
+  double coeffs_rs_HBA[POL_DEGREE + 1] = {
+      2.643667639489899e+05,  -0.007554769559170,    8.554267734878638e-11,
+      -4.722267097362212e-19, 1.251625317661630e-27, -1.236074872829482e-36};
+  /// system efficiency: roughly 1.0
+  double eta = 0.95;
+  double factor;
+};
+
+}  // namespace DPPP
+}  // namespace DP3
+
+#endif

DPPP/DPRun.cc

@@ -25,6 +25,7 @@
 
 #include <boost/algorithm/string.hpp>
 
+#include "AddNoise.h"
 #include "ApplyBeam.h"
 #include "ApplyCal.h"
 #include "Averager.h"
@@ -326,6 +327,8 @@ DPStep::ShPtr DPRun::makeSteps(const ParameterSet& parset, const string& prefix,
       step = std::make_shared<DDECal>(reader, parset, prefix);
     } else if (type == "interpolate") {
       step = std::make_shared<Interpolate>(reader, parset, prefix);
+    } else if (type == "addnoise") {
+      step = std::make_shared<AddNoise>(reader, parset, prefix);
     } else if (type == "out" || type == "output" || type == "msout") {
       step = makeOutputStep(reader, parset, prefix, currentMSName,
                             lastStep->outputs() == DPStep::MSType::BDA);

cmakecommand.txt

@@ -0,0 +1 @@
+cmake .. -DCFITSIO_INCLUDE_DIR="/homes/gheller/Work/cfitsio-3.48/build/include"