diff --git a/CMakeLists.txt b/CMakeLists.txt index 3adb7e6c..5d7866db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -155,6 +155,7 @@ add_library(IDGPredict OBJECT set(IDGPREDICT_OBJECT $) add_library(DPPP_OBJ OBJECT + DPPP/AddNoise.cc DPPP/Apply.cc DPPP/ApplyCal.cc DPPP/Averager.cc diff --git a/DPPP/AddNoise.cc b/DPPP/AddNoise.cc new file mode 100644 index 00000000..756adb3f --- /dev/null +++ b/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 . +// +// $Id: GainCal.cc 21598 2012-07-16 08:07:34Z diepen $ +// +// @author Claudio Gheller, Henrik Edler + +#include "AddNoise.h" + +#include + +#include "../Common/ParameterSet.h" +#include "../Common/Timer.h" + +#include +#include +#include +#include +#include +#include +#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::contiter indIter = itsBuf.getData().cbegin(); + + // Set the exposure + double exposure = buf.getExposure(); + + // Load the Antenna columns + Vector antenna1 = info().getAnt1(); + Vector 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& 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 chan_freq = getInfo().chanFreqs(); + Vector 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::contiter run1Iter = runBuf1.getData().cbegin(); + DPBuffer runBuf2; + runBuf2.copy(buf); + Array::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 distribution(0.0, stddev); + + for (int icorr = 0; icorr < n_corr; icorr++) { + double noise_real = distribution(generator); + double noise_img = distribution(generator); + std::complex 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 distribution(0.0, stddev); + + for (int icorr = 0; icorr < n_corr; icorr++) { + double noise_real = distribution(generator); + double noise_img = distribution(generator); + std::complex 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 diff --git a/DPPP/AddNoise.h b/DPPP/AddNoise.h new file mode 100644 index 00000000..251d04b5 --- /dev/null +++ b/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 . + +/// @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 + +#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 diff --git a/DPPP/DPRun.cc b/DPPP/DPRun.cc index 3e2e41bb..6f9a9a66 100644 --- a/DPPP/DPRun.cc +++ b/DPPP/DPRun.cc @@ -25,6 +25,7 @@ #include +#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(reader, parset, prefix); } else if (type == "interpolate") { step = std::make_shared(reader, parset, prefix); + } else if (type == "addnoise") { + step = std::make_shared(reader, parset, prefix); } else if (type == "out" || type == "output" || type == "msout") { step = makeOutputStep(reader, parset, prefix, currentMSName, lastStep->outputs() == DPStep::MSType::BDA); diff --git a/cmakecommand.txt b/cmakecommand.txt new file mode 100644 index 00000000..f0b14499 --- /dev/null +++ b/cmakecommand.txt @@ -0,0 +1 @@ +cmake .. -DCFITSIO_INCLUDE_DIR="/homes/gheller/Work/cfitsio-3.48/build/include"