Task #10121: reintegrate task branch: experimental direction dependent solver in DPPP

.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

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()

CEP/DP3/DPPP_DDECal/CMakeLists.txt

+# $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)

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/CMakeLists.txt

+# $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})

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Constraint.h

+#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

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/DDECal.h

+//# 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

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/KLFitter.h

+#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

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/PiercePoint.h

+#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

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/Register.h

+//# 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

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/ScreenConstraint.h

+#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

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/multidirsolver.h

+#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

CEP/DP3/DPPP_DDECal/include/DPPP_DDECal/screenfitter.h

+//# 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

CEP/DP3/DPPP_DDECal/src/CMakeLists.txt

+# $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)

CEP/DP3/DPPP_DDECal/src/Constraint.cc

+#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>();
+}

CEP/DP3/DPPP_DDECal/src/KLFitter.cc

+#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;
+}

CEP/DP3/DPPP_DDECal/src/PiercePoint.cc

+#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);
+};

CEP/DP3/DPPP_DDECal/src/Register.cc

+//# 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);
+}

CEP/DP3/DPPP_DDECal/src/ScreenConstraint.cc

+#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;
+}

CEP/DP3/DPPP_DDECal/src/multidirsolver.cc

+#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);
+  }
+}