From a96fdc6b6c719b0e3cfa02bc3e43f1cfe95e5265 Mon Sep 17 00:00:00 2001
From: Tammo Jan Dijkema <dijkema@astron.nl>
Date: Wed, 19 Apr 2017 11:20:59 +0000
Subject: [PATCH] Task #10121: reintegrate task branch: experimental direction
dependent solver in DPPP
---
.gitattributes | 19 +
CEP/DP3/CMakeLists.txt | 6 +
CEP/DP3/DPPP_DDECal/CMakeLists.txt | 11 +
.../include/DPPP_DDECal/CMakeLists.txt | 12 +
.../include/DPPP_DDECal/Constraint.h | 137 +++++
.../DPPP_DDECal/include/DPPP_DDECal/DDECal.h | 147 +++++
.../include/DPPP_DDECal/KLFitter.h | 36 ++
.../include/DPPP_DDECal/PiercePoint.h | 41 ++
.../include/DPPP_DDECal/Register.h | 37 ++
.../include/DPPP_DDECal/ScreenConstraint.h | 40 ++
.../include/DPPP_DDECal/multidirsolver.h | 79 +++
.../include/DPPP_DDECal/screenfitter.h | 45 ++
CEP/DP3/DPPP_DDECal/src/CMakeLists.txt | 11 +
CEP/DP3/DPPP_DDECal/src/Constraint.cc | 131 ++++
CEP/DP3/DPPP_DDECal/src/DDECal.cc | 557 ++++++++++++++++++
CEP/DP3/DPPP_DDECal/src/KLFitter.cc | 43 ++
CEP/DP3/DPPP_DDECal/src/PiercePoint.cc | 50 ++
CEP/DP3/DPPP_DDECal/src/Register.cc | 36 ++
CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc | 158 +++++
CEP/DP3/DPPP_DDECal/src/multidirsolver.cc | 209 +++++++
CEP/DP3/DPPP_DDECal/src/screenfitter.cc | 8 +
CMake/LofarPackageList.cmake | 1 +
22 files changed, 1814 insertions(+)
create mode 100644 CEP/DP3/DPPP_DDECal/CMakeLists.txt
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/CMakeLists.txt
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/KLFitter.h
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/PiercePoint.h
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Register.h
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/multidirsolver.h
create mode 100644 CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/screenfitter.h
create mode 100644 CEP/DP3/DPPP_DDECal/src/CMakeLists.txt
create mode 100644 CEP/DP3/DPPP_DDECal/src/Constraint.cc
create mode 100644 CEP/DP3/DPPP_DDECal/src/DDECal.cc
create mode 100644 CEP/DP3/DPPP_DDECal/src/KLFitter.cc
create mode 100644 CEP/DP3/DPPP_DDECal/src/PiercePoint.cc
create mode 100644 CEP/DP3/DPPP_DDECal/src/Register.cc
create mode 100644 CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc
create mode 100644 CEP/DP3/DPPP_DDECal/src/multidirsolver.cc
create mode 100644 CEP/DP3/DPPP_DDECal/src/screenfitter.cc
diff --git a/.gitattributes b/.gitattributes
index 7c3febba9f2..99b198590cd 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -932,6 +932,25 @@ CEP/DP3/DPPP_AOFlag/test/CMakeLists.txt -text
CEP/DP3/DPPP_AOFlag/test/tAOFlaggerStep.cc -text
CEP/DP3/DPPP_AOFlag/test/tAOFlaggerStep.run -text
CEP/DP3/DPPP_AOFlag/test/tAOFlaggerStep.sh -text
+CEP/DP3/DPPP_DDECal/CMakeLists.txt -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/CMakeLists.txt -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/KLFitter.h -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/PiercePoint.h -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Register.h -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/multidirsolver.h -text
+CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/screenfitter.h -text
+CEP/DP3/DPPP_DDECal/src/CMakeLists.txt -text
+CEP/DP3/DPPP_DDECal/src/Constraint.cc -text
+CEP/DP3/DPPP_DDECal/src/DDECal.cc -text
+CEP/DP3/DPPP_DDECal/src/KLFitter.cc -text
+CEP/DP3/DPPP_DDECal/src/PiercePoint.cc -text
+CEP/DP3/DPPP_DDECal/src/Register.cc -text
+CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc -text
+CEP/DP3/DPPP_DDECal/src/multidirsolver.cc -text
+CEP/DP3/DPPP_DDECal/src/screenfitter.cc -text
CEP/DP3/PythonDPPP/CMakeLists.txt -text
CEP/DP3/PythonDPPP/include/PythonDPPP/CMakeLists.txt -text
CEP/DP3/PythonDPPP/include/PythonDPPP/DPStepBase.h -text
diff --git a/CEP/DP3/CMakeLists.txt b/CEP/DP3/CMakeLists.txt
index 7c693c95b83..2ef1a8b00ff 100644
--- a/CEP/DP3/CMakeLists.txt
+++ b/CEP/DP3/CMakeLists.txt
@@ -7,3 +7,9 @@ lofar_add_package(DPPP_AOFlag)
lofar_add_package(SPW_Combine SPWCombine)
lofar_add_package(AOFlagger)
+lofar_find_package(Armadillo)
+if(${ARMADILLO_FOUND})
+ lofar_add_package(DPPP_DDECal)
+else()
+ message(WARNING "Armadillo was not found, NOT building DPPP_DDECal")
+endif()
diff --git a/CEP/DP3/DPPP_DDECal/CMakeLists.txt b/CEP/DP3/DPPP_DDECal/CMakeLists.txt
new file mode 100644
index 00000000000..4ebe58fb959
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/CMakeLists.txt
@@ -0,0 +1,11 @@
+# $Id: CMakeLists.txt 27640 2013-12-04 08:02:49Z diepen $
+
+lofar_package(DPPP_DDECal 1.0 DEPENDS DPPP)
+
+include(LofarFindPackage)
+lofar_find_package(Casacore COMPONENTS casa ms tables REQUIRED)
+lofar_find_package(Armadillo REQUIRED)
+
+add_subdirectory(include/DPPP_DDECal)
+add_subdirectory(src)
+# add_subdirectory(test)
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/CMakeLists.txt b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/CMakeLists.txt
new file mode 100644
index 00000000000..4abee5ccfdd
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/CMakeLists.txt
@@ -0,0 +1,12 @@
+# $Id: CMakeLists.txt 30990 2015-02-12 12:27:47Z diepen $
+
+# List of header files that will be installed.
+set(inst_HEADERS Register.h DDECal.h multidirsolver.h H5Parm.h Constraint.h ScreenConstraint.h PiercePoint.h)
+
+# Create symbolic link to include directory.
+execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_BINARY_DIR}/include/${PACKAGE_NAME})
+
+# Install header files.
+#install(FILES ${inst_HEADERS} DESTINATION include/${PACKAGE_NAME})
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h
new file mode 100644
index 00000000000..e3a42b03ebe
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h
@@ -0,0 +1,137 @@
+#ifndef CONSTRAINT_H
+#define CONSTRAINT_H
+
+#ifdef AOPROJECT
+#include "phasefitter.h"
+#define UPTR std::unique_ptr
+#else
+#include <DPPP/phasefitter.h>
+#define UPTR std::auto_ptr
+#endif
+
+#include <complex>
+#include <memory>
+#include <set>
+#include <vector>
+
+/**
+ * This class is the base class for classes that implement a constraint on
+ * calibration solutions. Constraints are used to increase
+ * the converge of calibration by applying these inside the solving step.
+ *
+ * The MultiDirSolver class uses this class for constrained calibration.
+ */
+class Constraint
+{
+public:
+ typedef std::complex<double> dcomplex;
+ struct Result
+ {
+ public:
+ std::vector<double> vals;
+ std::string axes; // Comma-separated string with axes names, fastest varying last
+ std::vector<size_t> dims;
+ std::string name;
+ };
+
+ virtual ~Constraint() { }
+
+ /**
+ * TODO To be removed: constraints should have their own init specific initialization methods.
+ */
+ virtual void init(size_t nAntennas, size_t nDirections,
+ size_t nChannelBlocks, const double* frequencies) = 0;
+
+ /**
+ * This method applies the constraints to the solutions.
+ * @param solutions is an array of array, such that:
+ * - solutions[ch] is a pointer for channelblock ch to antenna x directions solutions.
+ * - directions is the dimension with the fastest changing index.
+ * @param time Central time of interval.
+ */
+ virtual std::vector<Result> Apply(
+ std::vector<std::vector<dcomplex> >& solutions,
+ double time) = 0;
+};
+
+/**
+ * Awkwardly named, this class constraints the amplitudes of the solution to be unity, but
+ * keeps the phase.
+ */
+class PhaseConstraint : public Constraint
+{
+public:
+ PhaseConstraint() {};
+
+ virtual void init(size_t, size_t, size_t, const double*) {};
+
+ virtual std::vector<Result> Apply(
+ std::vector<std::vector<dcomplex> >& solutions,
+ double time);
+};
+
+class TECConstraint : public Constraint
+{
+public:
+ enum Mode {
+ /** Solve for both a (differential) TEC and an XX/YY-common scalar */
+ TECAndCommonScalarMode,
+ /** Solve only for a (differential) TEC value */
+ TECOnlyMode
+ };
+
+ TECConstraint(Mode mode, size_t nAntennas, size_t nDirections,
+ size_t nChannelBlocks, const double* frequencies);
+ TECConstraint(Mode mode);
+
+ virtual void init(size_t nAntennas, size_t nDirections,
+ size_t nChannelBlocks, const double* frequencies);
+
+ virtual std::vector<Result> Apply(
+ std::vector<std::vector<dcomplex> >& solutions,
+ double time);
+
+private:
+ Mode _mode;
+ size_t _nAntennas, _nDirections, _nChannelBlocks;
+ std::vector<PhaseFitter> _phaseFitters;
+};
+
+/**
+ * This constraint limits averages the solutions of a given list of antennas,
+ * so that they have equal solutions.
+ *
+ * The DDE solver uses this constraint to average the solutions of the core
+ * antennas. Core antennas are there defined by a given maximum distance from
+ * a reference antenna. In that case, the reference antenna is by default the first
+ * antenna. This constraint is meant to force all core stations to
+ * have the same solution, thereby decreasing the noise in their solutions.
+ */
+class CoreConstraint : public Constraint
+{
+public:
+ CoreConstraint() { }
+
+ virtual void init(size_t nAntennas, size_t nDirections,
+ size_t nChannelBlocks, const double*)
+ {
+ _nAntennas = nAntennas;
+ _nDirections = nDirections;
+ _nChannelBlocks = nChannelBlocks;
+ }
+
+ void setCoreAntennas(const std::set<size_t>& coreAntennas)
+ {
+ _coreAntennas = coreAntennas;
+ }
+
+ virtual std::vector<Result> Apply(
+ std::vector<std::vector<dcomplex> >& solutions,
+ double time);
+
+private:
+ size_t _nAntennas, _nDirections, _nChannelBlocks;
+ std::set<size_t> _coreAntennas;
+};
+
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
new file mode 100644
index 00000000000..319660eeca6
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h
@@ -0,0 +1,147 @@
+//# DDE.h: DPPP step class to calibrate direction dependent gains
+//# 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: DDECal.h 21598 2012-07-16 08:07:34Z diepen $
+//#
+//# @author Tammo Jan Dijkema
+
+#ifndef DPPP_DDECAL_H
+#define DPPP_DDECAL_H
+
+// @file
+// @brief DPPP step class to apply a calibration correction to the data
+
+#include <DPPP/DPInput.h>
+#include <DPPP/DPBuffer.h>
+#include <DPPP/H5Parm.h>
+#include <DPPP/BaselineSelection.h>
+#include <DPPP/Patch.h>
+#include <DPPP/Predict.h>
+#include <DPPP/SourceDBUtil.h>
+#include <DPPP/ApplyBeam.h>
+#include <DPPP_DDECal/multidirsolver.h>
+#include <DPPP_DDECal/Constraint.h>
+#include <StationResponse/Station.h>
+#include <StationResponse/Types.h>
+#include <ParmDB/Parm.h>
+#include <casa/Arrays/Cube.h>
+#include <casa/Quanta/MVEpoch.h>
+#include <measures/Measures/MEpoch.h>
+#include <casa/Arrays/ArrayMath.h>
+
+namespace LOFAR {
+
+ class ParameterSet;
+
+ namespace DPPP {
+ // @ingroup NDPPP
+
+ // This class is a DPStep class to calibrate (direction independent) gains.
+
+ typedef vector<Patch::ConstPtr> PatchList;
+ typedef std::pair<size_t, size_t >Baseline;
+
+ class DDECal: public DPStep
+ {
+ public:
+ // Construct the object.
+ // Parameters are obtained from the parset using the given prefix.
+ DDECal (DPInput*, const ParameterSet&, const string& prefix);
+
+ virtual ~DDECal();
+
+ // Create an DDECal object using the given parset.
+ static DPStep::ShPtr makeStep (DPInput*, const ParameterSet&,
+ const std::string&);
+
+ // Process the data.
+ // It keeps the data.
+ // When processed, it invokes the process function of the next step.
+ virtual bool process (const DPBuffer&);
+
+ // Initialize H5parm-file
+ void initH5parm();
+
+ // Write out the solutions
+ void writeSolutions();
+
+ // 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:
+ // Initialize the parmdb
+ void initH5Parm();
+
+ //# Data members.
+ DPInput* itsInput;
+ string itsName;
+ vector<DPBuffer> itsBufs;
+
+ // The time of the current buffer (in case of solint, average time)
+ double itsAvgTime;
+ std::vector<casa::Complex*> itsDataPtrs;
+
+ // For each timeslot, a vector of nDir buffers
+ std::vector<std::vector<casa::Complex*> > itsModelDataPtrs;
+
+ // For each time, for each channel block, a vector of size nAntennas * nDirections
+ std::vector<std::vector<std::vector<casa::DComplex> > > itsSols;
+
+ // For each time, for each constraint, a vector of results (e.g. tec and phase)
+ std::vector<std::vector<std::vector<Constraint::Result> > > itsConstraintSols;
+
+ casa::Cube<casa::Complex> itsModelData;
+ string itsH5ParmName;
+ H5Parm itsH5Parm;
+
+ string itsMode;
+ uint itsTimeStep;
+ uint itsSolInt;
+ uint itsStepInSolInt;
+ uint itsNChan;
+ uint itsNFreqCells;
+ vector<vector<string> > itsDirections; // For each direction, a vector of patches
+ vector<casa::CountedPtr<Constraint> > itsConstraints;
+
+ vector<Predict> itsPredictSteps;
+ vector<MultiResultStep::ShPtr> itsResultSteps; // For each directions, a multiresultstep with all times
+
+ NSTimer itsTimer;
+ NSTimer itsTimerPredict;
+ NSTimer itsTimerSolve;
+ NSTimer itsTimerWrite;
+ double itsCoreConstraint;
+
+ MultiDirSolver itsMultiDirSolver;
+ };
+
+ } //# end namespace
+}
+
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/KLFitter.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/KLFitter.h
new file mode 100644
index 00000000000..301d5686afa
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/KLFitter.h
@@ -0,0 +1,36 @@
+#ifndef KLFITTER_H
+#define KLFITTER_H
+#include<vector>
+#include <armadillo>
+#include <DPPP_DDECal/PiercePoint.h>
+
+using namespace arma;
+
+class KLFitter
+{//creates KH base and fits screens from collection of PiercePoints
+public:
+ KLFitter(double r0=1000.,double beta=5./3.,int order=9);
+ void calculateCorrMatrix(const vector<PiercePoint> pp);
+ void doFit();
+ size_t getOrder() const {return itsOrder;}
+ double* PhaseData() { return _phases.memptr(); }
+ double* WData() { return _weights.memptr(); }
+ double* ParData() { return itsPar.memptr(); }
+ double* TECFitWhiteData() { return itsTECFitWhite.memptr(); }
+ double* PPData() { return itsPiercePoints.memptr(); }
+
+private:
+ size_t itsOrder;
+ double itsR0,itsBeta;
+ Mat<double> itsPiercePoints;
+ Col<double> _phases;
+ Mat<double> _weights; //weights of the data points
+ Mat<double> itsCorrMatrix; //Correlation Matrix for KL fit
+ Mat<double> itsinvC; //for quick interpolation
+ Mat<double> itsU;
+ Mat<double> itsinvU;
+ Mat<double> itsTECFitWhite;
+ Col<double> itsPar;
+};
+
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/PiercePoint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/PiercePoint.h
new file mode 100644
index 00000000000..a39ef99a99b
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/PiercePoint.h
@@ -0,0 +1,41 @@
+#ifndef PIERCEPOINT_H
+#define PIERCEPOINT_H
+
+#include <vector>
+#include <measures/Measures/MDirection.h>
+#include <measures/Measures/MCDirection.h>
+#include <measures/Measures/MCPosition.h>
+#include <measures/Measures/MPosition.h>
+#include <measures/Measures/MeasConvert.h>
+#include <measures/Measures/MEpoch.h>
+
+#include <armadillo>
+
+using namespace arma;
+
+class PiercePoint
+{
+public:
+ PiercePoint();
+ PiercePoint(const casa::MPosition &ant,const casa::MDirection &source,const double height);
+ PiercePoint(const casa::MPosition &ant,const casa::MDirection &source);
+ void init(const casa::MPosition &ant,const casa::MDirection &source,const double height);
+ void evaluate(casa::MEpoch time);
+ Col<double> getValue() const {return itsValue;}
+ casa::MPosition getPos() const {return itsPosition;}
+ casa::MDirection getDir() const {return itsDirection;}
+private:
+ static const double IONOheight; //default height in meter
+ static const double EarthRadius; //default Earth radius in meter
+ //station position
+ casa::MPosition itsPosition;
+ //source position
+ casa::MDirection itsDirection;
+ // Ionospheric layer height.
+ double itsIonoHeight;
+ // square of length antenna vector (int ITRF) minus square of vector to piercepoint. This is constant for a assumed spherical Earth
+ double itsC;
+ Col<double> itsValue; //PiercePoint in ITRF coordinates
+};
+
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Register.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Register.h
new file mode 100644
index 00000000000..8144c212546
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Register.h
@@ -0,0 +1,37 @@
+//# Register.h: Register AOFlag steps in DPPP
+//# Copyright (C) 2015
+//# 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: AORFlagger.h 26900 2013-10-08 20:12:58Z loose $
+//#
+//# @author Ger van Diepen
+
+#ifndef DPPP_DDECAL_REGISTER_H
+#define DPPP_DDECAL_REGISTER_H
+
+// @file
+// @brief Register AOFlag steps in DPPP
+
+
+// Define the function (without name mangling) to register the 'constructor'.
+extern "C"
+{
+ void register_ddecal();
+}
+
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h
new file mode 100644
index 00000000000..6683a99582f
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h
@@ -0,0 +1,40 @@
+#ifndef SCREEN_CONSTRAINT_H
+#define SCREEN_CONSTRAINT_H
+
+#include <DPPP_DDECal/multidirsolver.h>
+#include <DPPP_DDECal/PiercePoint.h>
+#include <DPPP_DDECal/KLFitter.h>
+
+#include <cmath>
+#include <vector>
+#include <memory>
+
+class ScreenConstraint : public Constraint
+{
+public:
+ ScreenConstraint();
+ virtual void init(size_t nAntennas, size_t nDirections,
+ size_t nChannelBlocks, const double* frequencies);
+
+ virtual std::vector<Constraint::Result> Apply(std::vector<std::vector<MultiDirSolver::DComplex> >& solutions,double time);
+ virtual void CalculatePiercepoints();
+
+ void setAntennaPositions(const std::vector<std::vector<double> > antenna_pos);
+ void setDirections(const std::vector<std::pair<double, double> > source_pos);
+ void setTime(double time);
+ void initPiercePoints();
+
+private:
+ size_t _nAntennas, _nDirections, _nChannelBlocks;
+ std::vector<std::vector<double> > itsAntennaPos;
+ std::vector<std::vector<double> > itsSourcePos;
+ std::vector<double> itsFrequencies;
+ // antenna positions
+ // source positions
+ // measures instance ofzo
+ std::vector<std::vector<PiercePoint> > itsPiercePoints; //temporary hold calculated piercepoints per antenna
+ std::vector<KLFitter> _screenFitters;
+ double itsCurrentTime;
+};
+
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/multidirsolver.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/multidirsolver.h
new file mode 100644
index 00000000000..32af51b3c36
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/multidirsolver.h
@@ -0,0 +1,79 @@
+#ifndef MULTI_DIR_SOLVER_H
+#define MULTI_DIR_SOLVER_H
+
+#ifdef AOPROJECT
+#include "phasefitter.h"
+#include "Constraint.h"
+#define UPTR std::unique_ptr
+#else
+#include <DPPP/phasefitter.h>
+#include <DPPP_DDECal/Constraint.h>
+#define UPTR std::auto_ptr
+#endif
+
+#include <armadillo>
+
+#include <complex>
+#include <vector>
+#include <memory>
+
+class MultiDirSolver
+{
+public:
+ typedef std::complex<double> DComplex;
+ typedef std::complex<float> Complex;
+
+ struct SolveResult {
+ size_t iterations;
+ std::vector<std::vector<Constraint::Result> > _results;
+ };
+
+ enum CalibrationMode { CalibrateComplexGain,
+ CalibrateTEC1,
+ CalibrateTEC2,
+ CalibratePhase };
+
+ MultiDirSolver(size_t maxIterations, double accuracy, double stepSize);
+
+ void init(size_t nAntennas, size_t nDirections, size_t nChannels,
+ const std::vector<int>& ant1, const std::vector<int>& ant2);
+
+ // data[i] is een pointer naar de data voor tijdstap i, vanaf die pointer staat het in volgorde als in MS (bl, chan, pol)
+ // mdata[i] is een pointer voor tijdstap i naar arrays van ndir model data pointers (elk van die data pointers staat in zelfde volgorde als data)
+ // solutions[ch] is een pointer voor channelblock ch naar antenna x directions oplossingen.
+ SolveResult process(std::vector<Complex*>& data, std::vector<std::vector<Complex* > >& modelData,
+ std::vector<std::vector<DComplex> >& solutions, double time) const;
+
+ void set_mode(CalibrationMode mode) { _mode = mode; }
+
+ void set_channel_blocks(size_t nChannelBlocks) { _nChannelBlocks = nChannelBlocks; }
+
+ void set_max_iterations(size_t maxIterations) { _maxIterations = maxIterations; }
+
+ void set_accuracy(double accuracy) { _accuracy = accuracy; }
+
+ void set_step_size(double stepSize) { _stepSize = stepSize; }
+
+ void add_constraint(Constraint* constraint) { _constraints.push_back(constraint); }
+
+private:
+ void performSolveIteration(size_t channelBlockIndex,
+ std::vector<arma::cx_mat>& gTimesCs,
+ std::vector<arma::cx_vec>& vs,
+ const std::vector<DComplex>& solutions,
+ std::vector<DComplex>& nextSolutions,
+ const std::vector<Complex *>& data,
+ const std::vector<std::vector<Complex *> >& modelData) const;
+
+ size_t _nAntennas, _nDirections, _nChannels, _nChannelBlocks;
+ std::vector<int> _ant1, _ant2;
+
+ // Calibration setup
+ enum CalibrationMode _mode;
+ size_t _maxIterations;
+ double _accuracy;
+ double _stepSize;
+ std::vector<Constraint*> _constraints;
+};
+
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/screenfitter.h b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/screenfitter.h
new file mode 100644
index 00000000000..b148b1ec114
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/screenfitter.h
@@ -0,0 +1,45 @@
+//# screenfitter.h: Class to perform screen fitting
+//# Copyright (C) 2016
+//# 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/>.
+//#
+//#
+//# @author Maaijke Mevius
+
+/**
+ * @file screenfitter.h Implements TEC model screen filter @ref ScreenFitter.
+ * @author Maaijke Mevius
+ * @date 2017-02-01
+ */
+
+#ifndef SCREEN_FITTER_H
+#define SCREEN_FITTER_H
+#include <armadillo>
+#include <vector>
+using namespace arma;
+
+class ScreenFitter{
+ public:
+ ScreenFitter();
+ double* PhaseData() { return _phases.data(); }
+
+ private:
+ std::vector<double> _phases,_frequencies, _weights;
+ mat _corrmatrix; //correlation matrix
+
+};
+#endif
diff --git a/CEP/DP3/DPPP_DDECal/src/CMakeLists.txt b/CEP/DP3/DPPP_DDECal/src/CMakeLists.txt
new file mode 100644
index 00000000000..0e43f83595f
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/CMakeLists.txt
@@ -0,0 +1,11 @@
+# $Id: CMakeLists.txt 30439 2014-11-19 15:04:34Z dijkema $
+
+include(LofarPackageVersion)
+
+lofar_add_library(dppp_ddecal
+ Package__Version.cc
+ DDECal.cc Register.cc
+ KLFitter.cc DDECal.cc multidirsolver.cc Constraint.cc ScreenConstraint.cc PiercePoint.cc
+)
+
+lofar_add_bin_program(versiondppp_ddecal versiondppp_ddecal.cc)
diff --git a/CEP/DP3/DPPP_DDECal/src/Constraint.cc b/CEP/DP3/DPPP_DDECal/src/Constraint.cc
new file mode 100644
index 00000000000..53f1abd4177
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/Constraint.cc
@@ -0,0 +1,131 @@
+#ifdef AOPROJECT
+#include "Constraint.h"
+#include <omp.h> // for tec constraints
+#else
+#include <DPPP_DDECal/Constraint.h>
+#include <Common/OpenMP.h>
+#endif
+
+
+std::vector<Constraint::Result> PhaseConstraint::Apply(
+ std::vector<std::vector<dcomplex> >& solutions, double)
+{
+ for (uint ch=0; ch<solutions.size(); ++ch) {
+ for (uint solIndex=0; solIndex<solutions[ch].size(); ++solIndex) {
+ solutions[ch][solIndex] /= std::abs(solutions[ch][solIndex]);
+ }
+ }
+
+ return std::vector<Constraint::Result>();
+}
+
+TECConstraint::TECConstraint(Mode mode) :
+ _mode(mode),
+ _nAntennas(0)
+ ,_nDirections(0),
+ _nChannelBlocks(0),
+ _phaseFitters()
+{
+}
+
+TECConstraint::TECConstraint(Mode mode, size_t nAntennas, size_t nDirections,
+ size_t nChannelBlocks, const double* frequencies) :
+ _mode(mode)
+{
+ init(nAntennas, nDirections, nChannelBlocks, frequencies);
+}
+
+void TECConstraint::init(size_t nAntennas, size_t nDirections,
+ size_t nChannelBlocks, const double* frequencies) {
+ _nAntennas = nAntennas;
+ _nDirections = nDirections;
+ _nChannelBlocks = nChannelBlocks;
+ _phaseFitters.resize(
+#ifdef AOPROJECT
+ omp_get_max_threads()
+#else
+ LOFAR::OpenMP::maxThreads()
+#endif
+ );
+
+ for(size_t i=0; i!=_phaseFitters.size(); ++i)
+ {
+ _phaseFitters[i].SetChannelCount(_nChannelBlocks);
+ std::memcpy(_phaseFitters[i].FrequencyData(), frequencies, sizeof(double) * _nChannelBlocks);
+ }
+}
+
+std::vector<Constraint::Result> TECConstraint::Apply(
+ std::vector<std::vector<dcomplex> >& solutions, double)
+{
+ std::vector<Constraint::Result> res(2);
+
+ res[0].vals.resize(_nAntennas*_nDirections);
+ res[0].axes="ant,dir,freq";
+ res[0].name="tec";
+ res[0].dims.resize(3);
+ res[0].dims[0]=_nAntennas;
+ res[0].dims[1]=_nDirections;
+ res[0].dims[2]=1;
+ res[1]=res[0];
+ res[1].name="scalarphase";
+
+#pragma omp parallel for
+ for(size_t solutionIndex = 0; solutionIndex<_nAntennas*_nDirections; ++solutionIndex)
+ {
+ size_t thread =
+#ifdef AOPROJECT
+ omp_get_thread_num();
+#else
+ LOFAR::OpenMP::threadNum();
+#endif
+
+ for(size_t ch=0; ch!=_nChannelBlocks; ++ch) {
+ _phaseFitters[thread].PhaseData()[ch] = std::arg(solutions[ch][solutionIndex]);
+ }
+
+ double alpha, beta=0.0;
+ if(_mode == TECOnlyMode) {
+ _phaseFitters[thread].FitDataToTEC1Model(alpha);
+ } else {
+ _phaseFitters[thread].FitDataToTEC2Model(alpha, beta);
+ }
+
+ res[0].vals[solutionIndex] = alpha / 8.44797245e9;
+ res[1].vals[solutionIndex] = beta;
+
+ for(size_t ch=0; ch!=_nChannelBlocks; ++ch) {
+ solutions[ch][solutionIndex] = std::polar<double>(1.0, _phaseFitters[thread].PhaseData()[ch]);
+ }
+ }
+
+ return res;
+}
+
+std::vector<Constraint::Result> CoreConstraint::Apply(
+ std::vector<std::vector<dcomplex> >& solutions, double)
+{
+ for (uint ch=0; ch<solutions.size(); ++ch) {
+ std::vector<dcomplex> coreSolutions(_nDirections, 0.0);
+ // Calculate the sum of solutions over the core stations
+ for(std::set<size_t>::const_iterator antennaIter = _coreAntennas.begin(); antennaIter!=_coreAntennas.end(); ++antennaIter)
+ {
+ size_t startIndex = (*antennaIter)*_nDirections;
+ for(size_t direction = 0; direction != _nDirections; ++direction)
+ coreSolutions[direction] += solutions[ch][startIndex + direction];
+ }
+
+ // Divide by nr of core stations to get the mean solution
+ for(std::vector<dcomplex>::iterator solutionIter = coreSolutions.begin(); solutionIter!=coreSolutions.end(); ++solutionIter)
+ (*solutionIter) /= _coreAntennas.size();
+
+ // Assign all core stations to the mean solution
+ for(std::set<size_t>::const_iterator antennaIter = _coreAntennas.begin(); antennaIter!=_coreAntennas.end(); ++antennaIter)
+ {
+ size_t startIndex = (*antennaIter)*_nDirections;
+ for(size_t direction = 0; direction != _nDirections; ++direction)
+ solutions[ch][startIndex + direction] = coreSolutions[direction];
+ }
+ }
+ return std::vector<Constraint::Result>();
+}
diff --git a/CEP/DP3/DPPP_DDECal/src/DDECal.cc b/CEP/DP3/DPPP_DDECal/src/DDECal.cc
new file mode 100644
index 00000000000..bb094e25a01
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/DDECal.cc
@@ -0,0 +1,557 @@
+//# DDECal.cc: DPPP step class to do a direction dependent gain calibration
+//# 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: DDECal.cc 21598 2012-07-16 08:07:34Z diepen $
+//#
+//# @author Tammo Jan Dijkema
+
+#include <lofar_config.h>
+#include <DPPP_DDECal/DDECal.h>
+#include <DPPP/Simulate.h>
+#include <DPPP/ApplyCal.h>
+#include <DPPP/DPBuffer.h>
+#include <DPPP/DPInfo.h>
+#include <DPPP/SourceDBUtil.h>
+#include <DPPP/MSReader.h>
+#include <DPPP/DPLogger.h>
+#include <DPPP_DDECal/ScreenConstraint.h>
+#include <ParmDB/ParmDB.h>
+#include <ParmDB/ParmValue.h>
+#include <ParmDB/SourceDB.h>
+#include <Common/ParameterSet.h>
+#include <Common/StringUtil.h>
+#include <Common/LofarLogger.h>
+#include <Common/OpenMP.h>
+
+#include <fstream>
+#include <ctime>
+#include <utility>
+
+#include <casa/Arrays/ArrayMath.h>
+#include <casa/Arrays/MatrixMath.h>
+#include <measures/Measures/MEpoch.h>
+#include <measures/Measures/MeasConvert.h>
+#include <measures/Measures/MCDirection.h>
+#include <casa/Quanta/Quantum.h>
+#include <casa/OS/File.h>
+
+#include <vector>
+#include <algorithm>
+
+#include <limits>
+#include <iostream>
+#include <iomanip>
+
+using namespace casa;
+using namespace LOFAR::BBS;
+
+namespace LOFAR {
+ namespace DPPP {
+
+ DDECal::DDECal (DPInput* input,
+ const ParameterSet& parset,
+ const string& prefix)
+ : itsInput (input),
+ itsName (prefix),
+ itsAvgTime (0),
+ itsSols (),
+ itsH5ParmName (parset.getString (prefix + "h5parm",
+ parset.getString("msin")+
+ "/instrument.h5")),
+ itsH5Parm (itsH5ParmName),
+ itsTimeStep (0),
+ itsSolInt (parset.getInt (prefix + "solint", 1)),
+ itsStepInSolInt (0),
+ itsNChan (parset.getInt (prefix + "nchan", 0)),
+ itsNFreqCells (0),
+ itsCoreConstraint(parset.getDouble (prefix + "coreconstraint", 0.0)),
+ itsMultiDirSolver(parset.getInt (prefix + "maxiter", 50),
+ parset.getDouble (prefix + "tolerance", 1.e-5),
+ parset.getDouble (prefix + "stepsize", 0.2))
+ {
+ vector<string> strDirections = parset.getStringVector (prefix + "directions");
+ itsDirections.resize(strDirections.size());
+ for (uint i=0; i<strDirections.size(); ++i) {
+ ParameterValue dirStr(strDirections[i]);
+ itsDirections[i] = dirStr.getStringVector();
+ }
+
+ itsMode = toLower(parset.getString(prefix + "mode", "complexgain"));
+ if(itsCoreConstraint != 0.0) {
+ itsConstraints.push_back(casa::CountedPtr<Constraint>(
+ new CoreConstraint()));
+ }
+ if (itsMode == "phaseonly") {
+ itsConstraints.push_back(casa::CountedPtr<Constraint>(
+ new PhaseConstraint()));
+ } else if (itsMode == "tec") {
+ itsConstraints.push_back(casa::CountedPtr<Constraint>(
+ new TECConstraint(TECConstraint::TECOnlyMode)));
+ } else if (itsMode == "tecandphase"){
+ itsConstraints.push_back(casa::CountedPtr<Constraint>(
+ new TECConstraint(TECConstraint::TECAndCommonScalarMode)));
+ } else if (itsMode == "complexgain") {
+ // no constraints
+ } else if (itsMode == "tecscreen") {
+ itsConstraints.push_back(casa::CountedPtr<Constraint>(
+ new ScreenConstraint()));
+ } else {
+ THROW (Exception, "Unexpected mode: " << itsMode);
+ }
+
+ itsDataPtrs.resize(itsSolInt);
+ itsModelDataPtrs.resize(itsSolInt);
+ for (uint t=0; t<itsSolInt; ++t) {
+ itsModelDataPtrs[t].resize(itsDirections.size());
+ }
+ for (uint i=0;i<itsConstraints.size();i++) {
+ itsMultiDirSolver.add_constraint(itsConstraints[i].get());
+ }
+
+ const size_t nDir = itsDirections.size();
+
+ itsBufs.resize(itsSolInt);
+
+ itsPredictSteps.resize(nDir);
+
+ itsResultSteps.resize(nDir);
+ for (size_t dir=0; dir<nDir; ++dir) {
+ itsResultSteps[dir] = MultiResultStep::ShPtr(new MultiResultStep(itsSolInt));
+ itsPredictSteps[dir]=Predict(input, parset, prefix, itsDirections[dir]);
+ itsPredictSteps[dir].setNextStep(itsResultSteps[dir]);
+ }
+ }
+
+ DDECal::~DDECal()
+ {}
+
+ DPStep::ShPtr DDECal::makeStep (DPInput* input,
+ const ParameterSet& parset,
+ const std::string& prefix)
+ {
+ return DPStep::ShPtr(new DDECal(input, parset, prefix));
+ }
+
+ void DDECal::updateInfo (const DPInfo& infoIn)
+ {
+ info() = infoIn;
+ info().setNeedVisData();
+
+ const size_t nDir=itsDirections.size();
+
+ for (size_t dir=0; dir<nDir; ++dir) {
+ itsPredictSteps[dir].updateInfo(infoIn);
+ }
+
+ if (itsSolInt==0) {
+ itsSolInt=info().ntime();
+ }
+
+ if (itsNChan==0) {
+ itsNChan = info().nchan();
+ }
+ if (itsNChan>info().nchan()) {
+ itsNChan=info().nchan();
+ }
+ itsNFreqCells = info().nchan() / itsNChan;
+ if (itsNChan*itsNFreqCells<info().nchan()) { // If last freq cell is smaller
+ itsNFreqCells++;
+ }
+
+ // Convert from casa::Vector to std::vector
+ vector<int> ant1(info().getAnt1().size());
+ vector<int> ant2(info().getAnt2().size());
+ for (uint i=0; i<ant1.size(); ++i) {
+ ant1[i]=info().getAnt1()[i];
+ ant2[i]=info().getAnt2()[i];
+ }
+
+ // Fill antenna info in H5Parm, need to convert from casa types to std types
+ std::vector<std::string> antennaNames(info().antennaNames().size());
+ std::vector<std::vector<double> > antennaPos(info().antennaPos().size());
+ for (uint i=0; i<info().antennaNames().size(); ++i) {
+ antennaNames[i]=info().antennaNames()[i];
+ casa::Quantum<casa::Vector<double> > pos = info().antennaPos()[i].get("m");
+ antennaPos[i].resize(3);
+ antennaPos[i][0] = pos.getValue()[0];
+ antennaPos[i][1] = pos.getValue()[1];
+ antennaPos[i][2] = pos.getValue()[2];
+ }
+
+ itsH5Parm.addAntennas(antennaNames, antennaPos);
+
+ std::vector<std::string> sourceNames(itsDirections.size());
+ std::vector<std::pair<double, double> > sourcePositions(itsDirections.size());
+ for (uint i=0; i<itsDirections.size(); ++i) {
+ sourceNames[i]=itsDirections[i][0]; // This only gives the name of the first patch
+ sourcePositions[i]=itsPredictSteps[i].getFirstDirection();
+ }
+ itsH5Parm.addSources(sourceNames, sourcePositions);
+
+ uint nSolTimes = (info().ntime()+itsSolInt-1)/itsSolInt;
+ itsSols.resize(nSolTimes);
+ itsConstraintSols.resize(nSolTimes);
+
+ vector<double> chanFreqs(info().nchan()); //nChannelBlocks
+ for (uint i=0;i<info().chanFreqs().size();++i) {
+ chanFreqs[i]=info().chanFreqs()[i];
+ }
+
+ for (uint i=0; i<itsConstraints.size();++i) {
+ itsConstraints[i]->init(
+ info().antennaNames().size(),
+ itsDirections.size(),
+ info().nchan(), //nChannelBlocks
+ &(chanFreqs[0])
+ );
+
+ CoreConstraint* coreConstraint = dynamic_cast<CoreConstraint*>(itsConstraints[i].get());
+ if(coreConstraint != 0)
+ {
+ double
+ refX = antennaPos[i][0],
+ refY = antennaPos[i][1],
+ refZ = antennaPos[i][2];
+ std::set<size_t> coreAntennaIndices;
+ const double coreDistSq = itsCoreConstraint*itsCoreConstraint;
+ for(size_t ant=0; ant!=antennaPos.size(); ++ant)
+ {
+ double
+ dx = refX - antennaPos[ant][0],
+ dy = refY - antennaPos[ant][1],
+ dz = refZ - antennaPos[ant][2],
+ distSq = dx*dx + dy*dy + dz*dz;
+ if(distSq <= coreDistSq)
+ coreAntennaIndices.insert(ant);
+ }
+ coreConstraint->setCoreAntennas(coreAntennaIndices);
+ }
+
+ ScreenConstraint* screenConstraint = dynamic_cast<ScreenConstraint*>(itsConstraints[i].get());
+ if(screenConstraint != 0)
+ {
+ screenConstraint->setAntennaPositions(antennaPos);
+ screenConstraint->setDirections(sourcePositions);
+ screenConstraint->initPiercePoints();
+ //itsConstraints.setHeight(...);
+ }
+ }
+
+ uint nSt = info().antennaNames().size();
+ itsMultiDirSolver.init(nSt, nDir, info().nchan(), ant1, ant2);
+ for (uint i=0; i<nSolTimes; ++i) {
+ itsSols[i].resize(info().nchan());
+ }
+ }
+
+ void DDECal::show (std::ostream& os) const
+ {
+ os << "DDECal " << itsName << endl;
+ os << " H5Parm: " << itsH5ParmName <<endl;
+ os << " solint: " << itsSolInt <<endl;
+ os << " nchan: " << itsNChan <<endl;
+ os << " directions: " << itsDirections << endl;
+ os << " mode (constraints): " << itsMode << endl;
+ os << " coreconstraint: " << itsCoreConstraint << endl;
+ for (uint i=0; i<itsPredictSteps.size(); ++i) {
+ itsPredictSteps[i].show(os);
+ }
+ }
+
+ void DDECal::showTimings (std::ostream& os, double duration) const
+ {
+ double totaltime=itsTimer.getElapsed();
+ os << " ";
+ FlagCounter::showPerc1 (os, itsTimer.getElapsed(), duration);
+ os << " DDECal " << itsName << endl;
+
+ os << " ";
+ FlagCounter::showPerc1 (os, itsTimerPredict.getElapsed(), totaltime);
+ os << " of it spent in predict" << endl;
+
+ os << " ";
+ FlagCounter::showPerc1 (os, itsTimerSolve.getElapsed(), totaltime);
+ os << " of it spent in estimating gains and computing residuals" << endl;
+
+ os << " ";
+ FlagCounter::showPerc1 (os, itsTimerWrite.getElapsed(), totaltime);
+ os << " of it spent in writing gain solutions to disk" << endl;
+ }
+
+ bool DDECal::process (const DPBuffer& bufin)
+ {
+ itsTimer.start();
+
+ itsBufs[itsStepInSolInt].copy(bufin);
+ itsDataPtrs[itsStepInSolInt] = itsBufs[itsStepInSolInt].getData().data();
+
+ itsInput->fetchUVW(bufin, itsBufs[itsStepInSolInt], itsTimer);
+ itsInput->fetchWeights(bufin, itsBufs[itsStepInSolInt], itsTimer);
+ itsInput->fetchFullResFlags(bufin, itsBufs[itsStepInSolInt], itsTimer);
+
+ itsTimerPredict.start();
+
+#pragma omp parallel for
+ for (size_t dir=0; dir<itsPredictSteps.size(); ++dir) {
+ itsPredictSteps[dir].process(itsBufs[itsStepInSolInt]);
+ itsModelDataPtrs[itsStepInSolInt][dir] =
+ itsResultSteps[dir]->get()[itsStepInSolInt].getData().data();
+ }
+
+ // Handle weights and flags
+ const size_t nBl = info().nbaselines();
+ const size_t nCh = info().nchan();
+ const size_t nCr = 4;
+
+ for (size_t ch=0; ch<nCh; ++ch) {
+ for (size_t bl=0; bl<nBl; ++bl) {
+ for (size_t cr=0; cr<nCr; ++cr) {
+ if (itsBufs[itsStepInSolInt].getFlags().data()[bl*nCr*nCh+ch*nCr+cr]) {
+ // Flagged points: set data and model to 0
+ itsDataPtrs[itsStepInSolInt][bl*nCr*nCh+ch*nCr+cr] = 0;
+ for (size_t dir=0; dir<itsPredictSteps.size(); ++dir) {
+ itsModelDataPtrs[itsStepInSolInt][dir][bl*nCr*nCh+ch*nCr+cr] = 0;
+ }
+ } else {
+ // Premultiply non-flagged data with sqrt(weight)
+ double weight = sqrt(itsBufs[itsStepInSolInt].getWeights().data()[bl*nCr*nCh+ch*nCr+cr]);
+ itsDataPtrs[itsStepInSolInt][bl*nCr*nCh+ch*nCr+cr] *= weight;
+ for (size_t dir=0; dir<itsPredictSteps.size(); ++dir) {
+ itsModelDataPtrs[itsStepInSolInt][dir][bl*nCr*nCh+ch*nCr+cr] *= weight;
+ }
+ }
+ }
+ }
+ }
+
+ itsTimerPredict.stop();
+
+ itsAvgTime += itsAvgTime + bufin.getTime();
+
+ if (itsStepInSolInt==itsSolInt-1) {
+ itsTimerSolve.start();
+ MultiDirSolver::SolveResult solveResult =
+ itsMultiDirSolver.process(itsDataPtrs, itsModelDataPtrs,
+ itsSols[itsTimeStep/itsSolInt],
+ itsAvgTime / itsSolInt);
+ itsTimerSolve.stop();
+
+ // Store constraint solutions
+ if (itsMode!="complexgain" && itsMode!="phaseonly") {
+ itsConstraintSols[itsTimeStep/itsSolInt]=solveResult._results;
+ }
+
+ itsStepInSolInt=0;
+ itsAvgTime=0;
+ for (size_t dir=0; dir<itsResultSteps.size(); ++dir) {
+ itsResultSteps[dir]->clear();
+ }
+ } else {
+ itsStepInSolInt++;
+ }
+
+ itsTimeStep++;
+ itsTimer.stop();
+
+ getNextStep()->process(bufin);
+
+ return false;
+ }
+
+ void DDECal::writeSolutions() {
+ itsTimer.start();
+ itsTimerWrite.start();
+
+ uint nDir = itsDirections.size();
+
+ vector<double> weights(info().nchan()*info().nantenna()*info().ntime()*itsDirections.size(),1.);
+
+ if (itsConstraintSols[0].empty()) {
+ // Record the actual iterands of the solver
+ vector<DComplex> sols(info().nchan()*info().nantenna()*info().ntime()*nDir);
+ size_t i=0;
+
+ // Put solutions in a contiguous piece of memory
+ for (uint dir=0; dir<nDir; ++dir) {
+ for (uint ant=0; ant<info().nantenna(); ++ant) {
+ for (uint chan=0; chan<info().nchan(); ++chan) {
+ for (uint time=0; time<itsSols.size(); ++time) {
+ ASSERT(!itsSols[time].empty());
+ sols[i] = itsSols[time][chan][ant*nDir+dir];
+ ++i;
+ }
+ }
+ }
+ }
+ string axesnames="dir,ant,freq,time";
+ vector<hsize_t> dims(4);
+ dims[0]=nDir;
+ dims[1]=info().nantenna();
+ dims[2]=info().nchan();
+ dims[3]=itsSols.size();
+
+ uint numsols=(itsMode=="complexgain"?2:1);
+ for (uint solnum=0; solnum<numsols; ++solnum) {
+ string solTabName;
+ if (itsMode=="phaseonly") {
+ solTabName = "phaseonly000";
+ itsH5Parm.addSolution(solTabName, "scalarphase", axesnames,
+ dims, sols, weights, false);
+ } else if (itsMode=="complexgain") {
+ if (solnum==0) {
+ solTabName = "scalarphase000";
+ itsH5Parm.addSolution(solTabName, "scalarphase", axesnames,
+ dims, sols, weights, false);
+ } else {
+ solTabName = "scalaramplitude000";
+ itsH5Parm.addSolution(solTabName, "scalaramplitude", axesnames,
+ dims, sols, weights, true);
+ }
+ } else {
+ THROW(Exception, "Constraint should have produced output");
+ }
+ // Tell H5Parm that all antennas and directions were used
+ // TODO: do this more cleanly
+ std::vector<std::string> antennaNames(info().antennaNames().size());
+ for (uint i=0; i<info().antennaNames().size(); ++i) {
+ antennaNames[i]=info().antennaNames()[i];
+ }
+ itsH5Parm.setSolAntennas(solTabName, antennaNames);
+
+ std::vector<std::string> sourceNames(itsDirections.size());
+ for (uint i=0; i<itsDirections.size(); ++i) {
+ sourceNames[i]=itsDirections[i][0]; // This only gives the name of the first patch
+ }
+ itsH5Parm.setSolSources(solTabName, sourceNames);
+
+ // Set channel to frequency of middle channel
+ vector<double> chanFreqs(info().nchan());
+ for (uint chan=0; chan<info().nchan(); ++chan) {
+ chanFreqs[chan] = info().chanFreqs()[chan];
+ }
+ itsH5Parm.setFreqs(solTabName, chanFreqs);
+
+ uint nSolTimes = (info().ntime()+itsSolInt-1)/itsSolInt;
+ vector<double> solTimes(nSolTimes);
+ double starttime=info().startTime();
+ for (uint t=0; t<nSolTimes; ++t) {
+ solTimes[t] = starttime+t*info().timeInterval()*itsSolInt+0.5*info().timeInterval();
+ }
+ // End TODO
+ itsH5Parm.setTimes(solTabName, solTimes);
+ } // solnums loop
+ } else {
+ // Record the Constraint::Result in the H5Parm
+
+ uint nConstraints = itsConstraintSols[0].size();
+
+ for (uint constraintNum=0; constraintNum<nConstraints; ++constraintNum) {
+ // Number of solution names, e.g. 2 for "TEC" and "ScalarPhase"
+ uint nSolNames = itsConstraintSols[0][constraintNum].size();
+ for (uint solNameNum=0; solNameNum<nSolNames; ++solNameNum) {
+ // Get the result of the constraint solution at first time to get metadata
+ Constraint::Result firstResult = itsConstraintSols[0][constraintNum][solNameNum];
+
+ vector<hsize_t> dims(firstResult.dims.size()+1); // Add time dimension
+ dims[dims.size()-1]=itsConstraintSols.size();
+ size_t numSols=dims[dims.size()-1];
+ for (uint i=0; i<dims.size()-1; ++i) {
+ dims[i] = firstResult.dims[i];
+ numSols *= dims[i];
+ }
+
+ string axesnames = firstResult.axes + ",time";
+
+ // Put solutions in a contiguous piece of memory
+ vector<double> sols(numSols);
+ size_t posInFlatSol=0;
+ for (uint i=0; i<firstResult.vals.size(); ++i) {
+ for (uint time=0; time<itsSols.size(); ++time) {
+ sols[posInFlatSol++] =
+ itsConstraintSols[time][constraintNum][solNameNum].vals[i];
+ }
+ }
+
+ string solTabName = firstResult.name+"000";
+ itsH5Parm.addSolution(solTabName, firstResult.name,
+ axesnames, dims, sols, weights);
+
+ // Tell H5Parm that all antennas and directions were used
+ // TODO: do this more cleanly
+ std::vector<std::string> antennaNames(info().antennaNames().size());
+ for (uint i=0; i<info().antennaNames().size(); ++i) {
+ antennaNames[i]=info().antennaNames()[i];
+ }
+ itsH5Parm.setSolAntennas(solTabName, antennaNames);
+
+ std::vector<std::string> sourceNames(itsDirections.size());
+ for (uint i=0; i<itsDirections.size(); ++i) {
+ sourceNames[i]=itsDirections[i][0]; // This only gives the name of the first patch
+ }
+ itsH5Parm.setSolSources(solTabName, sourceNames);
+
+ // Set channel to frequency of middle channel
+ vector<double> oneFreq(1);
+ oneFreq[0] = info().chanFreqs()[info().nchan()/2];
+ itsH5Parm.setFreqs(solTabName, oneFreq);
+
+ uint nSolTimes = (info().ntime()+itsSolInt-1)/itsSolInt;
+ vector<double> solTimes(nSolTimes);
+ double starttime=info().startTime();
+ for (uint t=0; t<nSolTimes; ++t) {
+ solTimes[t] = starttime+t*info().timeInterval()*itsSolInt+0.5*info().timeInterval();
+ }
+ itsH5Parm.setTimes(solTabName, solTimes);
+ // End TODO
+ }
+ }
+ }
+
+
+ itsTimerWrite.stop();
+ itsTimer.stop();
+ }
+
+ void DDECal::finish()
+ {
+ itsTimer.start();
+
+ if (itsStepInSolInt!=0) {
+ //shrink itsDataPtrs, itsModelDataPtrs
+ std::vector<casa::Complex*>(itsDataPtrs.begin(),
+ itsDataPtrs.begin()+itsStepInSolInt).swap(itsDataPtrs);
+ std::vector<std::vector<casa::Complex*> >(itsModelDataPtrs.begin(),
+ itsModelDataPtrs.begin()+itsStepInSolInt).swap(itsModelDataPtrs);
+ itsTimerSolve.start();
+ itsMultiDirSolver.process(itsDataPtrs, itsModelDataPtrs,
+ itsSols[itsTimeStep/itsSolInt],
+ itsAvgTime/itsStepInSolInt);
+ itsTimerSolve.stop();
+ }
+
+ writeSolutions();
+
+ itsTimer.stop();
+
+ // Let the next steps finish.
+ getNextStep()->finish();
+ }
+
+ } //# end namespace
+}
diff --git a/CEP/DP3/DPPP_DDECal/src/KLFitter.cc b/CEP/DP3/DPPP_DDECal/src/KLFitter.cc
new file mode 100644
index 00000000000..2a87522e16c
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/KLFitter.cc
@@ -0,0 +1,43 @@
+#include "DPPP_DDECal/KLFitter.h"
+using namespace arma;
+
+KLFitter::KLFitter(double r0,double beta,int order):
+ itsOrder(order),
+ itsR0(r0),
+ itsBeta(beta)
+{
+
+}
+
+
+
+void KLFitter::calculateCorrMatrix(const vector<PiercePoint> pp){
+ itsPiercePoints.set_size(pp.size(),3);
+ _phases.set_size(pp.size());
+ _weights=eye<mat>(pp.size(),pp.size());
+ Mat<double> Distance=zeros<mat>(pp.size(),pp.size());
+ for(size_t i=0; i<pp.size();i++){
+ Mat<double> A(pp[i].getValue().memptr(),1,3);
+ itsPiercePoints.row(i)=A;
+ }
+ for(size_t n=0;n<pp.size();n++)
+ for(size_t m=0;m<pp.size();m++)
+ for(size_t i=0;i<3;i++)
+ Distance(n,m)+=pow(itsPiercePoints.col(i)[n]-itsPiercePoints.col(i)[m],2);
+ itsCorrMatrix=-pow((Distance / ( itsR0*itsR0 ) ),( itsBeta / 2.0 ))/2.0;
+ itsinvC=pinv(itsCorrMatrix);
+ mat V,U;
+ Col<double> S;
+ svd(U,S,V,itsCorrMatrix);
+ itsU=U(span::all,span(0,itsOrder));
+ itsinvU=inv(itsU.t()*(_weights*itsU));
+}
+
+
+void KLFitter::doFit(){
+ Mat<double> A=itsU.t()*(_weights*_phases);
+ itsPar=itsinvU* A; //TODO add weight info
+ itsTECFitWhite=(itsinvC*(itsU*itsPar));
+
+ _phases=itsCorrMatrix*itsTECFitWhite;
+}
diff --git a/CEP/DP3/DPPP_DDECal/src/PiercePoint.cc b/CEP/DP3/DPPP_DDECal/src/PiercePoint.cc
new file mode 100644
index 00000000000..5f9c399c516
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/PiercePoint.cc
@@ -0,0 +1,50 @@
+#include <DPPP_DDECal/PiercePoint.h>
+using namespace arma;
+
+const double PiercePoint::IONOheight = 300000.; //default height in meter
+const double PiercePoint::EarthRadius = 6371000.; //default Earth radius in meter
+
+PiercePoint::PiercePoint():
+ itsValue(3)
+{
+ casa::MPosition ant;//ITRF
+ casa::MDirection source;//J2000 pole
+ init(ant,source,PiercePoint::IONOheight);
+}
+
+PiercePoint::PiercePoint(const casa::MPosition &ant,const casa::MDirection &source,const double height):
+ itsValue(3)
+{
+ init(ant,source,height);
+};
+
+PiercePoint::PiercePoint(const casa::MPosition &ant,const casa::MDirection &source):
+ itsValue(3)
+ {
+ init(ant,source,PiercePoint::IONOheight);
+};
+
+void PiercePoint::init(const casa::MPosition &ant,const casa::MDirection &source,const double height){
+
+ itsPosition=casa::MPosition::Convert(ant,casa::MPosition::ITRF)();
+ itsDirection=source;
+ itsIonoHeight=height;
+ const casa::MVPosition &mPosition = itsPosition.getValue();
+ itsC = mPosition(0)*mPosition(0)+mPosition(1)*mPosition(1)+mPosition(2)*mPosition(2)-
+ (itsIonoHeight+PiercePoint::EarthRadius)*(itsIonoHeight+PiercePoint::EarthRadius);
+
+}
+
+void PiercePoint::evaluate(casa::MEpoch time){
+ //Convert direction to ITRF vector
+ casa::MeasFrame myframe(itsPosition,time);
+ casa::MDirection::Ref myref(casa::MDirection::ITRF,myframe);
+ const casa::MDirection dir = casa::MDirection::Convert(itsDirection,myref)();
+ const casa::MVDirection &mDir = dir.getValue();
+ const casa::MVPosition &mPos = itsPosition.getValue();
+ double A = mDir(0)*mDir(0)+mDir(1)*mDir(1)+mDir(2)*mDir(2);
+ double B = mDir(0)*mPos(0) + mDir(1)*mPos(1) +mDir(2)*mPos(2);
+ double alpha = (-B + sqrt(B*B - A*itsC))/A;
+ for(uword i=0;i<3;i++)
+ itsValue(i) = mPos(i) + alpha*mDir(i);
+};
diff --git a/CEP/DP3/DPPP_DDECal/src/Register.cc b/CEP/DP3/DPPP_DDECal/src/Register.cc
new file mode 100644
index 00000000000..856e6fdf318
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/Register.cc
@@ -0,0 +1,36 @@
+//# Register.cc: Register steps in DPPP
+//# Copyright (C) 2015
+//# 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: AOFlaggerStep.cc 31423 2015-04-03 14:06:21Z dijkema $
+//#
+//# @author Ger van Diepen
+
+#include <lofar_config.h>
+#include <DPPP_DDECal/Register.h>
+#include <DPPP_DDECal/DDECal.h>
+#include <DPPP/DPRun.h>
+
+// Define the function to make the DDECal 'constructor' known.
+// Its suffix must be the (lowercase) name of the package (library).
+// Also make the SlidingFlagger known.
+void register_ddecal()
+{
+ LOFAR::DPPP::DPRun::registerStepCtor ("ddecal",
+ LOFAR::DPPP::DDECal::makeStep);
+}
diff --git a/CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc b/CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc
new file mode 100644
index 00000000000..0f00d3e1de6
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc
@@ -0,0 +1,158 @@
+#include <DPPP_DDECal/ScreenConstraint.h>
+#include <Common/OpenMP.h>
+
+
+ScreenConstraint::ScreenConstraint() :
+ _nAntennas(0),
+ _nDirections(0),
+ _nChannelBlocks(0),
+ itsCurrentTime(0)
+{
+}
+
+void ScreenConstraint::init(size_t nAntennas, size_t nDirections, size_t nChannelBlocks, const double* frequencies) {
+ _nAntennas = nAntennas;
+ _nDirections = nDirections;
+ _nChannelBlocks = nChannelBlocks;
+ itsFrequencies.resize(_nChannelBlocks);
+ std::memcpy( itsFrequencies.data(),frequencies, sizeof(double) * _nChannelBlocks);
+ itsAntennaPos.resize(_nAntennas);
+ itsSourcePos.resize(_nDirections);
+ itsPiercePoints.resize(_nAntennas);
+ for(uint i=0;i<itsPiercePoints.size();i++)
+ itsPiercePoints[i].resize(_nDirections);
+ _screenFitters.resize(_nAntennas);
+
+ /* for(size_t i=0; i!=_screenFitters.size(); ++i)
+ {
+ _screenFitters[i].SetChannelCount(_nChannelBlocks);
+ std::memcpy(_screenFitters[i].FrequencyData(), frequencies, sizeof(double) * _nChannelBlocks);
+ _screenFitters[i].SetPPCount(_nDirections);
+
+ }
+ */
+}
+
+void ScreenConstraint::setAntennaPositions(const std::vector<std::vector<double> > antenna_pos) {
+ for(uint i=0;i<antenna_pos.size();i++){
+ itsAntennaPos[i].resize(3);
+ for(int j=0;j<3;j++)
+ itsAntennaPos[i][j]=antenna_pos[i][j];
+ }
+}
+
+void ScreenConstraint::setDirections(const std::vector<std::pair<double, double> > source_pos) {
+ for(uint i=0;i<source_pos.size();i++){
+ itsSourcePos[i].resize(2);
+ itsSourcePos[i][0]=source_pos[i].first;
+ itsSourcePos[i][1]=source_pos[i].second;
+ }
+
+}
+
+void ScreenConstraint::initPiercePoints(){
+ for(uint ipos=0;ipos<itsAntennaPos.size();ipos++){
+ casa::MPosition ant(casa::MVPosition(itsAntennaPos[ipos][0],itsAntennaPos[ipos][1],itsAntennaPos[ipos][2]),
+ casa::MPosition::ITRF);
+ for(uint isrc=0;isrc<itsSourcePos.size();isrc++){
+ casa::MDirection src(casa::MVDirection(itsSourcePos[isrc][0],itsSourcePos[isrc][1]),
+ casa::MDirection::J2000);
+
+ itsPiercePoints[ipos][isrc]=PiercePoint(ant,src);
+ }
+ }
+}
+
+void ScreenConstraint::setTime(double time){
+ if (itsCurrentTime!=time){
+
+ itsCurrentTime=time;
+ CalculatePiercepoints();
+#pragma omp parallel for
+ for(uint ipos=0;ipos<itsAntennaPos.size();ipos++){
+ _screenFitters[ipos].calculateCorrMatrix(itsPiercePoints[ipos]);
+ }
+
+ }
+}
+
+void ScreenConstraint::CalculatePiercepoints(){
+ casa::MEpoch time(casa::MVEpoch(itsCurrentTime/(24.*3600.))); //convert to MJD
+ for (uint i=0;i<itsPiercePoints.size();i++)
+ for (uint j=0;j<itsPiercePoints[i].size();j++)
+ itsPiercePoints[i][j].evaluate(time);
+}
+
+
+std::vector<Constraint::Result> ScreenConstraint::Apply(std::vector<std::vector<MultiDirSolver::DComplex> >& solutions,double time) {
+ //check if we need to reinitialize piercepoints
+ setTime(time);
+
+ std::vector<Result> res(4);
+ size_t numberofPar=std::min(_screenFitters[0].getOrder(),_nDirections);
+ res[0].vals.resize(_nAntennas*numberofPar);
+ res[0].axes="antenna,par";
+ res[0].dims.resize(2);
+ res[0].dims[0]=_nAntennas;
+ res[0].dims[1]=numberofPar;
+ res[0].name="screenpar";
+ res[1].vals.resize(_nAntennas*_nDirections*3);
+ res[1].axes="antenna,dir,xyz";
+ res[1].dims.resize(3);
+ res[1].dims[0]=_nAntennas;
+ res[1].dims[1]=_nDirections;
+ res[1].dims[2]=3;
+ res[1].name="piercepoints";
+ res[2].vals.resize(_nAntennas*_nDirections);
+ res[2].axes="antenna,dir";
+ res[2].dims.resize(2);
+ res[2].dims[0]=_nAntennas;
+ res[2].dims[1]=_nDirections;
+ res[2].name="TECfitwhite";
+ res[3].vals.resize(_nAntennas*_nDirections);
+ res[3].axes="antenna,dir";
+ res[3].dims.resize(2);
+ res[3].dims[0]=_nAntennas;
+ res[3].dims[1]=_nDirections;
+ res[3].name="phases";
+
+ //TODOEstimate Weights
+
+
+#pragma omp parallel for
+ for(size_t antIndex = 0; antIndex<_nAntennas; ++antIndex)
+ {
+ for(size_t dirIndex = 0; dirIndex<_nDirections; ++dirIndex){
+ double avgTEC=0;
+ size_t solutionIndex=antIndex*_nDirections+dirIndex;
+ for(size_t ch=0;ch<_nChannelBlocks; ++ch){
+ double refphase=std::arg(solutions[ch][dirIndex]);
+ //TODO: more advance frequency averaging...
+ avgTEC += std::arg(solutions[ch][solutionIndex]*std::polar<double>(1.0,-1*refphase))*itsFrequencies[ch];
+ }
+ res[3].vals[antIndex*_nDirections+dirIndex]= avgTEC/_nChannelBlocks;
+ _screenFitters[antIndex].PhaseData()[dirIndex] = avgTEC/_nChannelBlocks;
+ }
+
+ _screenFitters[antIndex].doFit();
+
+ for(size_t dirIndex = 0; dirIndex<_nDirections; ++dirIndex){
+ size_t solutionIndex=antIndex*_nDirections+dirIndex;
+ for(size_t ch=0;ch<_nChannelBlocks; ++ch)
+ solutions[ch][solutionIndex] = std::polar<double>(1.0, _screenFitters[antIndex].PhaseData()[dirIndex]/itsFrequencies[ch]);
+ //res[3].vals[antIndex*_nDirections+dirIndex]= _screenFitters[antIndex].PhaseData()[dirIndex];
+ for (size_t i=0;i<3;i++)
+ res[1].vals[antIndex*_nDirections*3+dirIndex*3+i]= _screenFitters[antIndex].PPData()[i*_nDirections+dirIndex];
+ }
+ for(size_t i=0;i<numberofPar;i++)
+ res[0].vals[antIndex*numberofPar+i]= _screenFitters[antIndex].ParData()[i];
+ for(size_t dirIndex = 0; dirIndex<_nDirections; ++dirIndex)
+ res[2].vals[antIndex*_nDirections+dirIndex]=
+ _screenFitters[antIndex].TECFitWhiteData()[dirIndex];
+
+
+ }
+
+
+ return res;
+}
diff --git a/CEP/DP3/DPPP_DDECal/src/multidirsolver.cc b/CEP/DP3/DPPP_DDECal/src/multidirsolver.cc
new file mode 100644
index 00000000000..643cb4aadc0
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/multidirsolver.cc
@@ -0,0 +1,209 @@
+#ifdef AOPROJECT
+#include "multidirsolver.h"
+#else
+#include <DPPP_DDECal/multidirsolver.h>
+#endif
+
+using namespace arma;
+
+MultiDirSolver::MultiDirSolver(size_t maxIterations, double accuracy, double stepSize) :
+ _nAntennas(0),
+ _nDirections(0),
+ _nChannels(0),
+ _nChannelBlocks(0),
+ _mode(CalibrateComplexGain),
+ _maxIterations(maxIterations),
+ _accuracy(accuracy),
+ _stepSize(stepSize)
+{
+}
+
+void MultiDirSolver::init(size_t nAntennas,
+ size_t nDirections,
+ size_t nChannels,
+ const std::vector<int>& ant1,
+ const std::vector<int>& ant2)
+{
+ _nAntennas = nAntennas;
+ _nDirections = nDirections;
+ _nChannels = nChannels;
+ _nChannelBlocks = nChannels;
+ _ant1 = ant1;
+ _ant2 = ant2;
+}
+
+MultiDirSolver::SolveResult MultiDirSolver::process(std::vector<Complex *>& data,
+ std::vector<std::vector<Complex *> >& modelData,
+ std::vector<std::vector<DComplex> >& solutions, double time) const
+{
+ const size_t nTimes = data.size();
+ SolveResult result;
+
+ std::vector<std::vector<DComplex> > nextSolutions(_nChannelBlocks);
+ if (solutions.size() != _nChannelBlocks) {
+ cout << "Error: 'solutions' parameter does not have the right shape" << endl;
+ result.iterations = 0;
+ return result;
+ }
+
+ result._results.resize(_constraints.size());
+
+ // Model matrix ant x [N x D] and visibility matrix ant x [N x 1],
+ // for each channelblock
+ // The following loop allocates all structures
+ std::vector<std::vector<cx_mat> > gTimesCs(_nChannelBlocks);
+ std::vector<std::vector<cx_vec> > vs(_nChannelBlocks);
+ for(size_t chBlock=0; chBlock!=_nChannelBlocks; ++chBlock)
+ {
+ solutions[chBlock].assign(_nDirections * _nAntennas, 1.0);
+ nextSolutions[chBlock].resize(_nDirections * _nAntennas);
+ const size_t
+ channelIndexStart = chBlock * _nChannels / _nChannelBlocks,
+ channelIndexEnd = (chBlock+1) * _nChannels / _nChannelBlocks,
+ curChannelBlockSize = channelIndexEnd - channelIndexStart;
+ gTimesCs[chBlock].resize(_nAntennas);
+ vs[chBlock].resize(_nAntennas);
+
+ for(size_t ant=0; ant!=_nAntennas; ++ant)
+ {
+ // Model matrix [N x D] and visibility matrix x [N x 1]
+ // Also space for the auto correlation is reserved, but it will be set to 0.
+ gTimesCs[chBlock][ant] = cx_mat(_nAntennas * nTimes * curChannelBlockSize,
+ _nDirections, fill::zeros);
+ vs[chBlock][ant] = cx_vec(_nAntennas * nTimes * curChannelBlockSize, fill::zeros);
+ }
+ }
+
+ // TODO the data and model data needs to be preweighted.
+ // Maybe we can get a non-const pointer from DPPP, that saves copying/allocating
+
+ ///
+ /// Start iterating
+ ///
+ size_t iteration = 0;
+ double normSum = 0.0, sum = 0.0;
+ do {
+#pragma omp parallel for
+ for(size_t chBlock=0; chBlock<_nChannelBlocks; ++chBlock)
+ {
+ performSolveIteration(chBlock, gTimesCs[chBlock], vs[chBlock],
+ solutions[chBlock], nextSolutions[chBlock],
+ data, modelData);
+ }
+
+ // Move the solutions towards nextSolutions
+ // (the moved solutions are stored in 'nextSolutions')
+ for(size_t chBlock=0; chBlock!=_nChannelBlocks; ++chBlock)
+ {
+ for(size_t i=0; i!=_nAntennas*_nDirections; ++i)
+ {
+ nextSolutions[chBlock][i] = solutions[chBlock][i]*(1.0-_stepSize) +
+ nextSolutions[chBlock][i] * _stepSize;
+ }
+ }
+ for(size_t i=0; i!=_constraints.size(); ++i)
+ {
+ result._results[i] = _constraints[i]->Apply(nextSolutions, time);
+ }
+
+ // Calculate the norm of the difference between the old and new solutions
+ for(size_t chBlock=0; chBlock!=_nChannelBlocks; ++chBlock)
+ {
+ for(size_t i=0; i!=_nAntennas*_nDirections; ++i)
+ {
+ double e = std::norm(nextSolutions[chBlock][i] - solutions[chBlock][i]);
+ normSum += e;
+ sum += std::abs(solutions[chBlock][i]);
+
+ solutions[chBlock][i] = nextSolutions[chBlock][i];
+
+ // For debug: output the solutions of the first antenna
+ if(i<_nDirections && false)
+ {
+ std::cout << " |s_" << i << "|=|" << solutions[chBlock][i] << "|="
+ << std::abs(solutions[chBlock][i]);
+ }
+ }
+ }
+ normSum /= _nChannelBlocks*_nAntennas*_nDirections;
+ sum /= _nChannelBlocks*_nAntennas*_nDirections;
+ iteration++;
+
+ } while(iteration < _maxIterations && normSum/sum > _accuracy);
+
+ if(normSum/sum <= _accuracy)
+ result.iterations = iteration;
+ else
+ result.iterations = _maxIterations+1;
+ return result;
+}
+
+void MultiDirSolver::performSolveIteration(size_t channelBlockIndex,
+ std::vector<arma::cx_mat>& gTimesCs,
+ std::vector<arma::cx_vec>& vs,
+ const std::vector<DComplex>& solutions,
+ std::vector<DComplex>& nextSolutions,
+ const std::vector<Complex *>& data,
+ const std::vector<std::vector<Complex *> >& modelData) const
+{
+ const size_t
+ channelIndexStart = channelBlockIndex * _nChannels / _nChannelBlocks,
+ channelIndexEnd = (channelBlockIndex+1) * _nChannels / _nChannelBlocks,
+ curChannelBlockSize = channelIndexEnd - channelIndexStart,
+ nTimes = data.size();
+
+ // The following loop fills the matrices for all antennas
+ for(size_t timeIndex=0; timeIndex!=nTimes; ++timeIndex)
+ {
+ std::vector<const Complex*> modelPtrs(_nDirections);
+ for(size_t baseline=0; baseline!=_ant1.size(); ++baseline)
+ {
+ size_t antenna1 = _ant1[baseline];
+ size_t antenna2 = _ant2[baseline];
+ if(antenna1 != antenna2)
+ {
+ cx_mat& gTimesC1 = gTimesCs[antenna1];
+ cx_vec& v1 = vs[antenna1];
+ cx_mat& gTimesC2 = gTimesCs[antenna2];
+ cx_vec& v2 = vs[antenna2];
+ for(size_t d=0; d!=_nDirections; ++d)
+ modelPtrs[d] = modelData[timeIndex][d] + (channelIndexStart + baseline * _nChannels) * 4;
+ const Complex* dataPtr = data[timeIndex] + (channelIndexStart + baseline * _nChannels) * 4;
+ for(size_t ch=channelIndexStart; ch!=channelIndexEnd; ++ch)
+ {
+ size_t dataIndex1 = ch-channelIndexStart + (timeIndex + antenna1 * nTimes) * curChannelBlockSize;
+ size_t dataIndex2 = ch-channelIndexStart + (timeIndex + antenna2 * nTimes) * curChannelBlockSize;
+ //std::cout << "timeindex" << timeIndex << ' ';
+ for(size_t d=0; d!=_nDirections; ++d)
+ {
+ std::complex<double> predicted = modelPtrs[d][0] + modelPtrs[d][3];
+ //std::cout << predicted << ' ';
+
+ size_t solIndex1 = antenna1*_nDirections + d;
+ size_t solIndex2 = antenna2*_nDirections + d;
+ gTimesC1(dataIndex2, d) = solutions[solIndex2] * predicted;
+ gTimesC2(dataIndex1, d) = solutions[solIndex1] * predicted;
+
+ modelPtrs[d] += 4; // Goto the next 2x2 matrix.
+ }
+ v1(dataIndex2) = dataPtr[0] + dataPtr[3]; // Solve using Stokes I
+ v2(dataIndex1) = v1(dataIndex2);
+ dataPtr += 4; // Goto the next 2x2 matrix.
+ }
+ }
+ }
+ }
+
+ // The matrices have been filled; compute the linear solution
+ // for each antenna.
+ for(size_t ant=0; ant!=_nAntennas; ++ant)
+ {
+ //std::cout << ant << '\n';
+ cx_mat& gTimesC = gTimesCs[ant];
+ cx_vec& v = vs[ant];
+ // solve [g C] x = v
+ cx_vec x = solve(gTimesC, v);
+ for(size_t d=0; d!=_nDirections; ++d)
+ nextSolutions[ant*_nDirections + d] = x(d);
+ }
+}
diff --git a/CEP/DP3/DPPP_DDECal/src/screenfitter.cc b/CEP/DP3/DPPP_DDECal/src/screenfitter.cc
new file mode 100644
index 00000000000..18563530647
--- /dev/null
+++ b/CEP/DP3/DPPP_DDECal/src/screenfitter.cc
@@ -0,0 +1,8 @@
+
+#include <DPPP/screenfitter.h>
+
+ScreenFitter::ScreenFitter():
+ _phases(),
+ _frequencies(),
+ _weights()
+{ }
diff --git a/CMake/LofarPackageList.cmake b/CMake/LofarPackageList.cmake
index a1fca4d0c42..1b171d44aad 100644
--- a/CMake/LofarPackageList.cmake
+++ b/CMake/LofarPackageList.cmake
@@ -36,6 +36,7 @@ if(NOT DEFINED LOFAR_PACKAGE_LIST_INCLUDED)
set(TestDynDPPP_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/DP3/TestDynDPPP)
set(PythonDPPP_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/DP3/PythonDPPP)
set(DPPP_AOFlag_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/DP3/DPPP_AOFlag)
+ set(DPPP_DDECal_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/DP3/DPPP_DDECal)
set(SPW_Combine_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/DP3/SPWCombine)
set(AOFlagger_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/DP3/AOFlagger)
set(LofarFT_SOURCE_DIR ${CMAKE_SOURCE_DIR}/CEP/Imager/LofarFT)
--
GitLab