From 6493f3978729f8f428397f6ecdc7000cbdc20f35 Mon Sep 17 00:00:00 2001 From: Claudio Gheller Date: Fri, 9 Oct 2020 14:47:52 +0200 Subject: [PATCH 01/10] commit basic functions --- CMakeLists.txt | 20 +----- DPPP/AddNoiseLBA.cc | 170 ++++++++++++++++++++++++++++++++++++++++++++ DPPP/AddNoiseLBA.h | 94 ++++++++++++++++++++++++ DPPP/DPRun.cc | 3 + cmakecommand.txt | 1 + 5 files changed, 271 insertions(+), 17 deletions(-) create mode 100644 DPPP/AddNoiseLBA.cc create mode 100644 DPPP/AddNoiseLBA.h create mode 100644 cmakecommand.txt diff --git a/CMakeLists.txt b/CMakeLists.txt index 3adb7e6c..4b9cfb52 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -155,6 +155,7 @@ add_library(IDGPredict OBJECT set(IDGPREDICT_OBJECT $) add_library(DPPP_OBJ OBJECT + DPPP/AddNoiseLBA.cc DPPP/Apply.cc DPPP/ApplyCal.cc DPPP/Averager.cc @@ -185,9 +186,6 @@ add_library(DPPP_OBJ OBJECT DPPP/MedFlagger.cc DPPP/ModelComponent.cc DPPP/ModelComponentVisitor.cc - DPPP/MS.cc - DPPP/MSBDAReader.cc - DPPP/MSBDAWriter.cc DPPP/MSReader.cc DPPP/MSUpdater.cc DPPP/MSWriter.cc @@ -222,7 +220,8 @@ add_library(DPPP_OBJ OBJECT ) set(DPPP_OBJECT $) -find_package(AOFlagger 3 REQUIRED) +#find_package(AOFlagger 3 REQUIRED) +#find_package(aoflagger 2 REQUIRED) include_directories(${AOFLAGGER_INCLUDE_DIR}) add_library(AOFlaggerStep_OBJ OBJECT @@ -308,7 +307,6 @@ endif() set(TEST_FILENAMES DPPP/test/runtests.cc DPPP/test/unit/mock/MockStep.cc - DPPP/test/unit/fixtures/fDirectory.cc DPPP/test/unit/tApplyCal.cc DPPP/test/unit/tApplyCalH5.cc DPPP/test/unit/tAverager.cc @@ -338,9 +336,6 @@ if(Boost_VERSION_STRING VERSION_GREATER_EQUAL "1.59") list(APPEND TEST_FILENAMES DPPP/test/unit/tBDAAverager.cc DPPP/test/unit/tDPInfo.cc - DPPP/test/unit/tDPInput.cc - DPPP/test/unit/tMSBDAReader.cc - DPPP/test/unit/tMSBDAWriter.cc DPPP/test/unit/tScaleDataBDA.cc ) else() @@ -373,15 +368,6 @@ set_tests_properties( PROPERTIES WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} FIXTURES_SETUP dppp ) -add_test( - NAME extract_NDPPP_BDA - COMMAND ${CMAKE_COMMAND} -E tar xzf ${CMAKE_SOURCE_DIR}/DPPP/test/resources/tNDPPP_bda.in_MS.tgz -) -set_tests_properties( - extract_NDPPP_BDA - PROPERTIES WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} FIXTURES_SETUP dppp -) - add_executable( unittests EXCLUDE_FROM_ALL ${TEST_FILENAMES} ${OS_SPECIFIC_TESTS} diff --git a/DPPP/AddNoiseLBA.cc b/DPPP/AddNoiseLBA.cc new file mode 100644 index 00000000..78993d2a --- /dev/null +++ b/DPPP/AddNoiseLBA.cc @@ -0,0 +1,170 @@ +// AddNoiseLBA.h: DPPP step class to add 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 "AddNoiseLBA.h" + +#include + +#include "../Common/ParameterSet.h" +#include "../Common/Timer.h" + +#include +#include +#include +#include +#include +#include + +using namespace casacore; + +namespace DP3 { +namespace DPPP { + +AddNoiseLBA::AddNoiseLBA(DPInput* input, const ParameterSet& parset, + const string& prefix) + : itsInput(input) { + mode = parset.getInt(prefix + "mode", 0); +} + +AddNoiseLBA::~AddNoiseLBA() {} + +void AddNoiseLBA::updateInfo(const DPInfo& infoIn) { + info() = infoIn; + info().setNeedVisData(); + info().setWriteData(); +} + +void AddNoiseLBA::show(std::ostream& os) const { + os << "AddNoiseLBA " << itsName << endl; +} + +void AddNoiseLBA::showTimings(std::ostream& os, double duration) const { + os << " "; + FlagCounter::showPerc1(os, itsTimer.getElapsed(), duration); + os << " AddNoiseLBA " << itsName << endl; +} + +bool AddNoiseLBA::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(); + + // Set Number of baselines + int n_baselines = antenna1.size(); + // cout << "Number of baselines = " << n_baselines << endl; + + // Set the number of correlations + int n_corr = info().ncorr(); + + // Set the LOFAR_ANTENNA_SET + string antennaSet1 = "LBA_OUTER"; + string antennaSet2 = "LBA_INNER"; + string antennaSet = getInfo().antennaSet(); + ; + + // Set the Polynomial coefficients + double* coeff; + if (antennaSet.compare(antennaSet1)) coeff = coeffs_outer; + if (antennaSet.compare(antennaSet2)) coeff = coeffs_inner; + + // Set the number of frequency channels + int n_freq = getInfo().chanFreqs().size(); + // cout << "Number of frequencies = " << n_freq << endl; + 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; + 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 = eta * sefd; + 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++; + } + } + } + + 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 AddNoiseLBA::finish() { + // Let the next steps finish. + getNextStep()->finish(); +} +} // namespace DPPP +} // namespace DP3 diff --git a/DPPP/AddNoiseLBA.h b/DPPP/AddNoiseLBA.h new file mode 100644 index 00000000..1a618f83 --- /dev/null +++ b/DPPP/AddNoiseLBA.h @@ -0,0 +1,94 @@ +// AddNoiseLBA.h: DPPP step class to add 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 LBA random noise to data +/// @author Tammo Jan Dijkema +// + +#ifndef DPPP_AddNoiseLBA_H +#define DPPP_AddNoiseLBA_H + +#include "DPBuffer.h" +#include "DPInput.h" + +#include + +#define POL_DEGREE 5 + +namespace DP3 { + +class ParameterSet; + +namespace DPPP { +/// @brief DPPP step class to AddNoiseLBA visibilities from a source model + +class AddNoiseLBA : public DPStep { + public: + /// Construct the object. + /// Parameters are obtained from the parset using the given prefix. + AddNoiseLBA(DPInput*, const ParameterSet&, const string& prefix); + + virtual ~AddNoiseLBA(); + + /// Process the data calculating or adding to DATA the LBA 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; + /// Coefficients for polynomial interpolation (from constant -first- to highet + /// order -last-) + double coeffs_outer[POL_DEGREE + 1] = {4.46492043e+05, -4.04156579e-02, + 1.58636639e-09, -3.09364148e-17, + 2.93955326e-25, -1.06998148e-33}; + double coeffs_inner[POL_DEGREE + 1] = {8.32889327e+05, -8.93829326e-02, + 3.90153820e-09, -8.23245656e-17, + 8.35181243e-25, -3.25202160e-33}; + /// lba_mode can be: lba_inner or lba_outer. Read from the parset file + string lba_mode; + /// system efficiency: for the moment set to 1.0 + double eta = 1.0; +}; + +} // namespace DPPP +} // namespace DP3 + +#endif diff --git a/DPPP/DPRun.cc b/DPPP/DPRun.cc index 3e2e41bb..a1bba2dc 100644 --- a/DPPP/DPRun.cc +++ b/DPPP/DPRun.cc @@ -25,6 +25,7 @@ #include +#include "AddNoiseLBA.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 == "addnoiselba") { + step = DPStep::ShPtr(new AddNoiseLBA (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" From 999e671eebd015b5620acca729ccd62d6a1c1b9b Mon Sep 17 00:00:00 2001 From: Claudio Gheller Date: Fri, 9 Oct 2020 14:53:30 +0200 Subject: [PATCH 02/10] correction to CMakeList.txt --- CMakeLists.txt | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4b9cfb52..1bf2376c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -155,7 +155,7 @@ add_library(IDGPredict OBJECT set(IDGPREDICT_OBJECT $) add_library(DPPP_OBJ OBJECT - DPPP/AddNoiseLBA.cc + DPPP/AddNoiseLBA.cc DPPP/Apply.cc DPPP/ApplyCal.cc DPPP/Averager.cc @@ -186,6 +186,9 @@ add_library(DPPP_OBJ OBJECT DPPP/MedFlagger.cc DPPP/ModelComponent.cc DPPP/ModelComponentVisitor.cc + DPPP/MS.cc + DPPP/MSBDAReader.cc + DPPP/MSBDAWriter.cc DPPP/MSReader.cc DPPP/MSUpdater.cc DPPP/MSWriter.cc @@ -220,8 +223,7 @@ add_library(DPPP_OBJ OBJECT ) set(DPPP_OBJECT $) -#find_package(AOFlagger 3 REQUIRED) -#find_package(aoflagger 2 REQUIRED) +find_package(AOFlagger 3 REQUIRED) include_directories(${AOFLAGGER_INCLUDE_DIR}) add_library(AOFlaggerStep_OBJ OBJECT @@ -307,6 +309,7 @@ endif() set(TEST_FILENAMES DPPP/test/runtests.cc DPPP/test/unit/mock/MockStep.cc + DPPP/test/unit/fixtures/fDirectory.cc DPPP/test/unit/tApplyCal.cc DPPP/test/unit/tApplyCalH5.cc DPPP/test/unit/tAverager.cc @@ -336,6 +339,9 @@ if(Boost_VERSION_STRING VERSION_GREATER_EQUAL "1.59") list(APPEND TEST_FILENAMES DPPP/test/unit/tBDAAverager.cc DPPP/test/unit/tDPInfo.cc + DPPP/test/unit/tDPInput.cc + DPPP/test/unit/tMSBDAReader.cc + DPPP/test/unit/tMSBDAWriter.cc DPPP/test/unit/tScaleDataBDA.cc ) else() @@ -368,6 +374,15 @@ set_tests_properties( PROPERTIES WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} FIXTURES_SETUP dppp ) +add_test( + NAME extract_NDPPP_BDA + COMMAND ${CMAKE_COMMAND} -E tar xzf ${CMAKE_SOURCE_DIR}/DPPP/test/resources/tNDPPP_bda.in_MS.tgz +) +set_tests_properties( + extract_NDPPP_BDA + PROPERTIES WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} FIXTURES_SETUP dppp +) + add_executable( unittests EXCLUDE_FROM_ALL ${TEST_FILENAMES} ${OS_SPECIFIC_TESTS} From dd54432e02617a69d549e1d186f4605fcf9f48d5 Mon Sep 17 00:00:00 2001 From: "claudio.gheller@gmail.com" Date: Wed, 14 Oct 2020 08:22:44 +0000 Subject: [PATCH 03/10] few style issues fixed --- DPPP/AddNoiseLBA.cc | 3 ++- DPPP/AddNoiseLBA.h | 6 +++--- DPPP/DPRun.cc | 2 +- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/DPPP/AddNoiseLBA.cc b/DPPP/AddNoiseLBA.cc index 78993d2a..a0ca09b4 100644 --- a/DPPP/AddNoiseLBA.cc +++ b/DPPP/AddNoiseLBA.cc @@ -68,7 +68,8 @@ bool AddNoiseLBA::process(const DPBuffer& buf) { itsTimer.start(); - // Name of the column to add the noise (at the moment not used, just a placeholder) + // Name of the column to add the noise (at the moment not used, just a + // placeholder) string column = "DATA"; DPBuffer itsBuf; itsBuf.copy(buf); diff --git a/DPPP/AddNoiseLBA.h b/DPPP/AddNoiseLBA.h index 1a618f83..09cfcf15 100644 --- a/DPPP/AddNoiseLBA.h +++ b/DPPP/AddNoiseLBA.h @@ -47,9 +47,9 @@ class AddNoiseLBA : public DPStep { virtual ~AddNoiseLBA(); - /// Process the data calculating or adding to DATA the LBA estimated random noise - /// (van Haarlem et al. 2013). - /// When processed, it invokes the process function of the next step. + /// Process the data calculating or adding to DATA the LBA 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. diff --git a/DPPP/DPRun.cc b/DPPP/DPRun.cc index a1bba2dc..67ab42a4 100644 --- a/DPPP/DPRun.cc +++ b/DPPP/DPRun.cc @@ -328,7 +328,7 @@ DPStep::ShPtr DPRun::makeSteps(const ParameterSet& parset, const string& prefix, } else if (type == "interpolate") { step = std::make_shared(reader, parset, prefix); } else if (type == "addnoiselba") { - step = DPStep::ShPtr(new AddNoiseLBA (reader, parset, prefix)); + 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); From a0976823cb2584a49303420c14a690297980e94e Mon Sep 17 00:00:00 2001 From: root Date: Mon, 9 Nov 2020 11:15:38 +0000 Subject: [PATCH 04/10] HBA implementation completed --- DPPP/AddNoise.cc | 227 +++++++++++++++++++++++++++++++++++++++++++++++ DPPP/AddNoise.h | 106 ++++++++++++++++++++++ DPPP/DPRun.cc | 6 +- 3 files changed, 336 insertions(+), 3 deletions(-) create mode 100644 DPPP/AddNoise.cc create mode 100644 DPPP/AddNoise.h diff --git a/DPPP/AddNoise.cc b/DPPP/AddNoise.cc new file mode 100644 index 00000000..54099f5e --- /dev/null +++ b/DPPP/AddNoise.cc @@ -0,0 +1,227 @@ +// 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); + antenna = parset.getString(prefix + "antenna", "LBA"); +} + +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(); + + // Set Number of baselines + int n_baselines = antenna1.size(); + // cout << "Number of baselines = " << n_baselines << endl; + + // Set the number of correlations + int n_corr = info().ncorr(); + + // Set the LOFAR_ANTENNA_SET + string antennaSet = getInfo().antennaSet(); + const Vector& allNames = getInfo().antennaNames(); + + // Set the Polynomial coefficients + double* coeff; + if (antennaSet.compare("LBA") || antennaSet.compare("LBA_ALL")) coeff = coeffs_all_LBA; + if (antennaSet.compare("LBA_INNER")) coeff = coeffs_inner_LBA; + if (antennaSet.compare("LBA_OUTER")) coeff = coeffs_outer_LBA; + + // Set the number of frequency channels + int n_freq = getInfo().chanFreqs().size(); + // cout << "Number of frequencies = " << n_freq << endl; + 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 (antenna == "LBA" || + antenna == "LBA_INNER" || + antenna == "LBA_OUTER" || + antenna == "LBA_ALL") + { + 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 = eta * sefd; + 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 (antenna == "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 = eta * sefd; + 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..f9b8bd00 --- /dev/null +++ b/DPPP/AddNoise.h @@ -0,0 +1,106 @@ +// 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 antenna parameter selects between: LBA=LBA_ALL, LBA_INNER + /// LBA_OUTER and HBA + string antenna; + /// 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}; + /// lba_mode can be: lba_inner or lba_outer. Read from the parset file + string lba_mode; + /// system efficiency: roughly 1.0 + double eta = 0.95; +}; + +} // namespace DPPP +} // namespace DP3 + +#endif diff --git a/DPPP/DPRun.cc b/DPPP/DPRun.cc index 67ab42a4..6f9a9a66 100644 --- a/DPPP/DPRun.cc +++ b/DPPP/DPRun.cc @@ -25,7 +25,7 @@ #include -#include "AddNoiseLBA.h" +#include "AddNoise.h" #include "ApplyBeam.h" #include "ApplyCal.h" #include "Averager.h" @@ -327,8 +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 == "addnoiselba") { - 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); From decf5bd177c32c8bee9344c16f1a1d51725022f5 Mon Sep 17 00:00:00 2001 From: "claudio.gheller@gmail.com" Date: Mon, 9 Nov 2020 12:08:16 +0000 Subject: [PATCH 05/10] class AddNoise added to the list --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1bf2376c..5d7866db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -155,7 +155,7 @@ add_library(IDGPredict OBJECT set(IDGPREDICT_OBJECT $) add_library(DPPP_OBJ OBJECT - DPPP/AddNoiseLBA.cc + DPPP/AddNoise.cc DPPP/Apply.cc DPPP/ApplyCal.cc DPPP/Averager.cc From 8d133982a10687f133e356facd8e1ae0f28bc6c4 Mon Sep 17 00:00:00 2001 From: Claudio Gheller Date: Tue, 10 Nov 2020 10:37:56 +0100 Subject: [PATCH 06/10] Deleted the NoiseLBA function from the git repository. Substituted by AddNoise, comprising both LBA and HBA --- DPPP/AddNoiseLBA.cc | 171 -------------------------------------------- DPPP/AddNoiseLBA.h | 94 ------------------------ 2 files changed, 265 deletions(-) delete mode 100644 DPPP/AddNoiseLBA.cc delete mode 100644 DPPP/AddNoiseLBA.h diff --git a/DPPP/AddNoiseLBA.cc b/DPPP/AddNoiseLBA.cc deleted file mode 100644 index a0ca09b4..00000000 --- a/DPPP/AddNoiseLBA.cc +++ /dev/null @@ -1,171 +0,0 @@ -// AddNoiseLBA.h: DPPP step class to add 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 "AddNoiseLBA.h" - -#include - -#include "../Common/ParameterSet.h" -#include "../Common/Timer.h" - -#include -#include -#include -#include -#include -#include - -using namespace casacore; - -namespace DP3 { -namespace DPPP { - -AddNoiseLBA::AddNoiseLBA(DPInput* input, const ParameterSet& parset, - const string& prefix) - : itsInput(input) { - mode = parset.getInt(prefix + "mode", 0); -} - -AddNoiseLBA::~AddNoiseLBA() {} - -void AddNoiseLBA::updateInfo(const DPInfo& infoIn) { - info() = infoIn; - info().setNeedVisData(); - info().setWriteData(); -} - -void AddNoiseLBA::show(std::ostream& os) const { - os << "AddNoiseLBA " << itsName << endl; -} - -void AddNoiseLBA::showTimings(std::ostream& os, double duration) const { - os << " "; - FlagCounter::showPerc1(os, itsTimer.getElapsed(), duration); - os << " AddNoiseLBA " << itsName << endl; -} - -bool AddNoiseLBA::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(); - - // Set Number of baselines - int n_baselines = antenna1.size(); - // cout << "Number of baselines = " << n_baselines << endl; - - // Set the number of correlations - int n_corr = info().ncorr(); - - // Set the LOFAR_ANTENNA_SET - string antennaSet1 = "LBA_OUTER"; - string antennaSet2 = "LBA_INNER"; - string antennaSet = getInfo().antennaSet(); - ; - - // Set the Polynomial coefficients - double* coeff; - if (antennaSet.compare(antennaSet1)) coeff = coeffs_outer; - if (antennaSet.compare(antennaSet2)) coeff = coeffs_inner; - - // Set the number of frequency channels - int n_freq = getInfo().chanFreqs().size(); - // cout << "Number of frequencies = " << n_freq << endl; - 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; - 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 = eta * sefd; - 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++; - } - } - } - - 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 AddNoiseLBA::finish() { - // Let the next steps finish. - getNextStep()->finish(); -} -} // namespace DPPP -} // namespace DP3 diff --git a/DPPP/AddNoiseLBA.h b/DPPP/AddNoiseLBA.h deleted file mode 100644 index 09cfcf15..00000000 --- a/DPPP/AddNoiseLBA.h +++ /dev/null @@ -1,94 +0,0 @@ -// AddNoiseLBA.h: DPPP step class to add 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 LBA random noise to data -/// @author Tammo Jan Dijkema -// - -#ifndef DPPP_AddNoiseLBA_H -#define DPPP_AddNoiseLBA_H - -#include "DPBuffer.h" -#include "DPInput.h" - -#include - -#define POL_DEGREE 5 - -namespace DP3 { - -class ParameterSet; - -namespace DPPP { -/// @brief DPPP step class to AddNoiseLBA visibilities from a source model - -class AddNoiseLBA : public DPStep { - public: - /// Construct the object. - /// Parameters are obtained from the parset using the given prefix. - AddNoiseLBA(DPInput*, const ParameterSet&, const string& prefix); - - virtual ~AddNoiseLBA(); - - /// Process the data calculating or adding to DATA the LBA 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; - /// Coefficients for polynomial interpolation (from constant -first- to highet - /// order -last-) - double coeffs_outer[POL_DEGREE + 1] = {4.46492043e+05, -4.04156579e-02, - 1.58636639e-09, -3.09364148e-17, - 2.93955326e-25, -1.06998148e-33}; - double coeffs_inner[POL_DEGREE + 1] = {8.32889327e+05, -8.93829326e-02, - 3.90153820e-09, -8.23245656e-17, - 8.35181243e-25, -3.25202160e-33}; - /// lba_mode can be: lba_inner or lba_outer. Read from the parset file - string lba_mode; - /// system efficiency: for the moment set to 1.0 - double eta = 1.0; -}; - -} // namespace DPPP -} // namespace DP3 - -#endif From 1c6f190e5fbc8c64ceb5ddc79aabe8c29a4f59ae Mon Sep 17 00:00:00 2001 From: Claudio Gheller Date: Tue, 10 Nov 2020 11:23:20 +0100 Subject: [PATCH 07/10] AntennaSet added --- DPPP/AddNoise.cc | 147 +++++++++++++++++++++++------------------------ DPPP/AddNoise.h | 40 ++++++------- 2 files changed, 90 insertions(+), 97 deletions(-) diff --git a/DPPP/AddNoise.cc b/DPPP/AddNoise.cc index 54099f5e..4894c5d4 100644 --- a/DPPP/AddNoise.cc +++ b/DPPP/AddNoise.cc @@ -42,10 +42,9 @@ namespace DP3 { namespace DPPP { AddNoise::AddNoise(DPInput* input, const ParameterSet& parset, - const string& prefix) + const string& prefix) : itsInput(input) { mode = parset.getInt(prefix + "mode", 0); - antenna = parset.getString(prefix + "antenna", "LBA"); } AddNoise::~AddNoise() {} @@ -67,10 +66,9 @@ void AddNoise::showTimings(std::ostream& os, double duration) const { } bool AddNoise::process(const DPBuffer& buf) { - itsTimer.start(); - // Name of the column to add the noise (at the moment not used, just a + // Name of the column to add the noise (at the moment not used, just a // placeholder) string column = "DATA"; DPBuffer itsBuf; @@ -85,6 +83,9 @@ bool AddNoise::process(const DPBuffer& buf) { Vector antenna1 = info().getAnt1(); Vector antenna2 = info().getAnt2(); + // Load the Antenna type + antennaSet = info().antennaSet(); + // Set Number of baselines int n_baselines = antenna1.size(); // cout << "Number of baselines = " << n_baselines << endl; @@ -92,13 +93,13 @@ bool AddNoise::process(const DPBuffer& buf) { // Set the number of correlations int n_corr = info().ncorr(); - // Set the LOFAR_ANTENNA_SET - string antennaSet = getInfo().antennaSet(); + // Set the Antenna names const Vector& allNames = getInfo().antennaNames(); // Set the Polynomial coefficients double* coeff; - if (antennaSet.compare("LBA") || antennaSet.compare("LBA_ALL")) coeff = coeffs_all_LBA; + if (antennaSet.compare("LBA") || antennaSet.compare("LBA_ALL")) + coeff = coeffs_all_LBA; if (antennaSet.compare("LBA_INNER")) coeff = coeffs_inner_LBA; if (antennaSet.compare("LBA_OUTER")) coeff = coeffs_outer_LBA; @@ -129,80 +130,74 @@ bool AddNoise::process(const DPBuffer& buf) { double sefd1; double sefd2; - if (antenna == "LBA" || - antenna == "LBA_INNER" || - antenna == "LBA_OUTER" || - antenna == "LBA_ALL") - { - 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 = eta * sefd; - 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++; + if (antennaSet.compare("LBA") || antennaSet.compare("LBA_INNER") || antennaSet.compare("LBA_OUTER") || + antennaSet.compare("LBA_ALL") { + 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 = eta * sefd; + 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 (antenna == "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; + } 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 = eta * sefd; + 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++; } - 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 = eta * sefd; - 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"); + } else { + throw Exception( + "Antenna should be HBA, LBA, LBA_INNER, LBA_OUTER, LBA_ALL"); } if (mode == 0) { diff --git a/DPPP/AddNoise.h b/DPPP/AddNoise.h index f9b8bd00..31b5cfa3 100644 --- a/DPPP/AddNoise.h +++ b/DPPP/AddNoise.h @@ -47,7 +47,7 @@ class AddNoise : public DPStep { virtual ~AddNoise(); - /// Process the data calculating or adding to DATA the estimated random + /// 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&); @@ -74,28 +74,26 @@ class AddNoise : public DPStep { /// 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 antenna parameter selects between: LBA=LBA_ALL, LBA_INNER - /// LBA_OUTER and HBA - string antenna; + /// The antennaSet parameter selects between: LBA=LBA_ALL, LBA_INNER + /// LBA_OUTER and HBA (read from MS) + const 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}; - /// lba_mode can be: lba_inner or lba_outer. Read from the parset file - string lba_mode; + 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; }; From 5bde439cc771921295e949f93cf8f623f0763d7b Mon Sep 17 00:00:00 2001 From: "claudio.gheller@gmail.com" Date: Tue, 10 Nov 2020 10:48:02 +0000 Subject: [PATCH 08/10] Some bugs fixed --- DPPP/AddNoise.cc | 4 +--- DPPP/AddNoise.h | 2 +- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/DPPP/AddNoise.cc b/DPPP/AddNoise.cc index 4894c5d4..2db28abb 100644 --- a/DPPP/AddNoise.cc +++ b/DPPP/AddNoise.cc @@ -88,7 +88,6 @@ bool AddNoise::process(const DPBuffer& buf) { // Set Number of baselines int n_baselines = antenna1.size(); - // cout << "Number of baselines = " << n_baselines << endl; // Set the number of correlations int n_corr = info().ncorr(); @@ -105,7 +104,6 @@ bool AddNoise::process(const DPBuffer& buf) { // Set the number of frequency channels int n_freq = getInfo().chanFreqs().size(); - // cout << "Number of frequencies = " << n_freq << endl; Vector chan_freq = getInfo().chanFreqs(); Vector chan_width = getInfo().chanWidths(); @@ -131,7 +129,7 @@ bool AddNoise::process(const DPBuffer& buf) { double sefd2; if (antennaSet.compare("LBA") || antennaSet.compare("LBA_INNER") || antennaSet.compare("LBA_OUTER") || - antennaSet.compare("LBA_ALL") { + antennaSet.compare("LBA_ALL")) { for (int ibase = ibegin; ibase < iend; ibase++) { for (int ifreq = 0; ifreq < n_freq; ifreq++) { nu = chan_freq(ifreq); diff --git a/DPPP/AddNoise.h b/DPPP/AddNoise.h index 31b5cfa3..0d5346d0 100644 --- a/DPPP/AddNoise.h +++ b/DPPP/AddNoise.h @@ -76,7 +76,7 @@ class AddNoise : public DPStep { int mode; /// The antennaSet parameter selects between: LBA=LBA_ALL, LBA_INNER /// LBA_OUTER and HBA (read from MS) - const string antennaSet; + string antennaSet; /// Coefficients for polynomial interpolation (from constant -first- to highet /// order -last-) double coeffs_all_LBA[POL_DEGREE + 1] = { From 05240de196eae21667da0ba278553b8667ea2d53 Mon Sep 17 00:00:00 2001 From: Claudio Gheller Date: Fri, 20 Nov 2020 10:43:41 +0100 Subject: [PATCH 09/10] optional multiplicative factor for variance added --- DPPP/AddNoise.cc | 10 ++++++---- DPPP/AddNoise.h | 1 + 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/DPPP/AddNoise.cc b/DPPP/AddNoise.cc index 2db28abb..3dcd9f2f 100644 --- a/DPPP/AddNoise.cc +++ b/DPPP/AddNoise.cc @@ -45,6 +45,7 @@ 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() {} @@ -128,15 +129,16 @@ bool AddNoise::process(const DPBuffer& buf) { double sefd1; double sefd2; - if (antennaSet.compare("LBA") || antennaSet.compare("LBA_INNER") || antennaSet.compare("LBA_OUTER") || - antennaSet.compare("LBA_ALL")) { + //if (antennaSet.compare("LBA") || antennaSet.compare("LBA_INNER") || antennaSet.compare("LBA_OUTER") || + // antennaSet.compare("LBA_ALL")) { + 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 = eta * sefd; + stddev = factor * sefd / eta; stddev = stddev / sqrt(2.0 * exposure * chan_width[ifreq]); std::normal_distribution distribution(0.0, stddev); @@ -176,7 +178,7 @@ bool AddNoise::process(const DPBuffer& buf) { coeff2[3] * pow(nu, 3.0) + coeff2[4] * pow(nu, 4.0) + coeff2[5] * pow(nu, 5.0); sefd = sqrt(sefd1 * sefd2); - stddev = eta * sefd; + stddev = factor * sefd / eta; stddev = stddev / sqrt(2.0 * exposure * chan_width[ifreq]); std::normal_distribution distribution(0.0, stddev); diff --git a/DPPP/AddNoise.h b/DPPP/AddNoise.h index 0d5346d0..251d04b5 100644 --- a/DPPP/AddNoise.h +++ b/DPPP/AddNoise.h @@ -96,6 +96,7 @@ class AddNoise : public DPStep { -4.722267097362212e-19, 1.251625317661630e-27, -1.236074872829482e-36}; /// system efficiency: roughly 1.0 double eta = 0.95; + double factor; }; } // namespace DPPP From 57f340977c45d02efc9ac01ecff5a10d3cd349f3 Mon Sep 17 00:00:00 2001 From: Claudio Gheller Date: Fri, 20 Nov 2020 11:10:28 +0100 Subject: [PATCH 10/10] bug fixed --- DPPP/AddNoise.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/DPPP/AddNoise.cc b/DPPP/AddNoise.cc index 3dcd9f2f..756adb3f 100644 --- a/DPPP/AddNoise.cc +++ b/DPPP/AddNoise.cc @@ -98,10 +98,10 @@ bool AddNoise::process(const DPBuffer& buf) { // Set the Polynomial coefficients double* coeff; - if (antennaSet.compare("LBA") || antennaSet.compare("LBA_ALL")) + if (antennaSet.compare("LBA") == 0 || antennaSet.compare("LBA_ALL") == 0) coeff = coeffs_all_LBA; - if (antennaSet.compare("LBA_INNER")) coeff = coeffs_inner_LBA; - if (antennaSet.compare("LBA_OUTER")) coeff = coeffs_outer_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(); @@ -129,8 +129,8 @@ bool AddNoise::process(const DPBuffer& buf) { double sefd1; double sefd2; - //if (antennaSet.compare("LBA") || antennaSet.compare("LBA_INNER") || antennaSet.compare("LBA_OUTER") || - // antennaSet.compare("LBA_ALL")) { + //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++) {