Line data    Source code

       1             : //# VisCalSolver2.h: Default solver for calibration using visibilities
       2             : //# Copyright (C) 1996,1997,2000,2001,2002,2003
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be adressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //#
      27             : 
      28             : #ifndef SYNTHESIS_VISCALSOL2_H
      29             : #define SYNTHESIS_VISCALSOL2_H
      30             : 
      31             : #include <casacore/casa/aips.h>
      32             : #include <casacore/casa/BasicSL/Complex.h>
      33             : #include <casacore/casa/BasicSL/Constants.h>
      34             : #include <casacore/casa/Arrays/Array.h>
      35             : #include <casacore/casa/Arrays/Vector.h>
      36             : #include <casacore/casa/Arrays/Matrix.h>
      37             : #include <synthesis/MeasurementComponents/VisCal.h>
      38             : #include <synthesis/MeasurementComponents/SolvableVisCal.h>
      39             : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
      40             : #include <synthesis/MeasurementEquations/VisEquation.h>
      41             : 
      42             : namespace casa { //# NAMESPACE CASA - BEGIN
      43             : 
      44             : // <summary> 
      45             : // VisCalSolver2: Default solver for calibration using visibility data
      46             : // </summary>
      47             : 
      48             : // <use visibility=export>
      49             : 
      50             : // <reviewed reviewer="" date="" tests="" demos="">
      51             : 
      52             : // <prerequisite>
      53             : //   <li> <linkto class="MeasurementComponents">MeasurementComponents</linkto> module
      54             : //   <li> <linkto class="VisEquation">VisEquation</linkto> module
      55             : // </prerequisite>
      56             : //
      57             : // <etymology>
      58             : // VisCal for visibility calibration (meaning solved from visibilities), Solver for
      59             : // solving. 
      60             : // </etymology>
      61             : //
      62             : // <synopsis> 
      63             : //
      64             : // VisCalSolver2 describes an interface for solving for calibration from visibility
      65             : // data in a VisBuffer.  It hosts the communication of visibility data with otherwise
      66             : // generic solving mechanims.  A Levenberg-Marquardt solver is supplied here by
      67             : // default, but it is intended that this class be the template for interfaces to
      68             : // third-party solving mechanisms.
      69             : //
      70             : // </synopsis> 
      71             : //
      72             : // <example>
      73             : // <srcblock>
      74             : 
      75             : // </srcblock>
      76             : // </example>
      77             : //
      78             : // <motivation>
      79             : // It is desirable to establish the distinct communicative boundary between generic
      80             : // solving mechanims and the particulars of visibility data and calibration
      81             : // component descriptions.  This class is intended to serve this purpose, by providing
      82             : // access to visibility data and calibration in terms of quantities necessary for
      83             : // least-squares (and other) style solvers.  
      84             : // </motivation>
      85             : //
      86             : // <todo asof="97/10/01">
      87             : // </todo>
      88             : 
      89             : // **********************************************************
      90             : //  VisCalSolver2
      91             : //
      92             : 
      93             : class VisCalSolver2 {
      94             : public:
      95             : 
      96             :   // Constructor
      97             :   VisCalSolver2();
      98             :   VisCalSolver2(casacore::String solmode,
      99             :                 casacore::Vector<casacore::Float>& rmsthresh);
     100             :   
     101             :   // Destructor
     102             :   ~VisCalSolver2();
     103             : 
     104             :   // Do the solve
     105             :   bool solve(VisEquation& viseq, SolvableVisCal& svc, SDBList& sdbs);
     106             : 
     107             :   // L1R-capable version
     108             :   bool solveL1R(VisEquation& viseq, SolvableVisCal& svc, SDBList& sdbs);
     109             : 
     110             : protected:
     111             : 
     112             :   // Access to fundamental external objects:
     113     2076889 :   inline SDBList&         sdbs() { return *SDBs_; };
     114      328675 :   inline VisEquation&     ve()   { return *ve_; };
     115      153230 :   inline SolvableVisCal&  svc()  { return *svc_; };
     116             : 
     117             :   // Access to maxIter_
     118       67430 :   inline int&    maxIter() { return maxIter_; };
     119             : 
     120             :   // Access to chi2
     121     1227737 :   inline double& chiSq()     { return chiSq_; };
     122    60020214 :   inline casacore::Vector<double>& chiSqV()     { return chiSqV_; };
     123      237557 :   inline double& lastChiSq() { return lastChiSq_; };
     124      391870 :   inline double& dChiSq()    { return dChiSq_; };
     125      523640 :   inline double& sumWt()     { return sumWt_; };
     126    60020214 :   inline casacore::Vector<double>& sumWtV()     { return sumWtV_; };
     127    59793524 :   inline int&    nWt()       { return nWt_; };
     128             : 
     129             :   // Access to parameters, & grad,hess,dp
     130     1422724 :   inline int&              nPar() { return nPar_; };
     131      614523 :   inline casacore::Vector<casacore::Complex>&  par()       { return par_; };
     132     1102826 :   inline casacore::Vector<bool>&     parOK()     { return parOK_; };
     133      168242 :   inline casacore::Vector<float>&    parErr()    { return parErr_; };
     134    48845597 :   inline casacore::Vector<casacore::DComplex>& grad()      { return grad_; };
     135    50158903 :   inline casacore::Vector<double>&   hess()      { return hess_; };
     136     2689482 :   inline casacore::Vector<casacore::Complex>&  dpar()      { return dpar_; };
     137      261683 :   inline casacore::Vector<casacore::Complex>&  lastPar()   { return lastPar_; };
     138             : 
     139      152020 :   inline double&           lambda()    { return lambda_; };
     140             : 
     141             :   // Initialize solving data
     142             :   void initSolve();
     143             : 
     144             :   // Obtain trial residuals w.r.t svc's current pars
     145             :   void residualate2();
     146             : 
     147             :   // Differentiate w.r.t svc's pars
     148             :   void differentiate2();
     149             : 
     150             :   // Calculate residuals (incl. diff'd) and chi2 
     151             :   void chiSquare2();
     152             : 
     153             :   // Apply RMS threshold to current residuals
     154             :   void RMSThresh(casacore::Int RejIter);
     155             : 
     156             :   // Check for convergence
     157             :   bool converged();
     158             : 
     159             :   // Internal solving methods
     160             :   void accGradHess2();
     161             :   void revert();
     162             :   void solveGradHess();
     163             :   void updatePar();
     164             : 
     165             :   // Optimize the step parabolically
     166             :   void optStepSize2();
     167             : 
     168             :   // Get and print par errors
     169             :   void getErrors();
     170             : 
     171             :   void printPar(const int& iter);
     172             : 
     173             : private:
     174             :   
     175             :   // Diagnostic print level
     176     1271101 :   inline int& prtlev() { return prtlev_; };
     177             : 
     178             :   // VisBuffer (from outside)
     179             :   SDBList* SDBs_;
     180             : 
     181             :   // VisEquation (from outside)
     182             :   VisEquation* ve_;
     183             : 
     184             :   // SVC (from outside)
     185             :   SolvableVisCal* svc_;
     186             : 
     187             :   // Total Number of parameters
     188             :   int nPar_;
     189             : 
     190             :   // Maximum number of solve iterations to attempt
     191             :   int maxIter_;
     192             : 
     193             :   // Chi2, sum wts
     194             :   double chiSq_;
     195             :   casacore::Vector<double> chiSqV_;
     196             :   double lastChiSq_;
     197             :   double dChiSq_;
     198             :   double sumWt_;
     199             :   casacore::Vector<double> sumWtV_;
     200             :   int    nWt_;
     201             :   int    cvrgcount_;
     202             : 
     203             :   // Parameter storage
     204             :   // (these are casacore::Complex to match the VisCal solvePar)
     205             :   casacore::Vector<casacore::Complex> par_;
     206             :   casacore::Vector<bool>    parOK_;
     207             :   casacore::Vector<float>   parErr_;
     208             :   casacore::Vector<casacore::Complex> lastPar_;
     209             : 
     210             :   // Parameter update
     211             :   casacore::Vector<casacore::Complex> dpar_;
     212             : 
     213             :   // Gradient, Hessian
     214             :   //  (these are double for precision in accumulation
     215             :   casacore::Vector<casacore::DComplex> grad_;
     216             :   casacore::Vector<double> hess_;
     217             : 
     218             :   // LM factor
     219             :   double lambda_;
     220             : 
     221             :   // Step optimization toggle
     222             :   bool optstep_;
     223             : 
     224             :   // Control L1 solution
     225             :   bool doL1_;
     226             :   casacore::Vector<casacore::Float> L1clamp_;
     227             : 
     228             :   // Control iterative rejection
     229             :   bool doRMSThresh_;
     230             :   casacore::Vector<casacore::Float> RMSThresh_;
     231             :   int nRMSThresh_;
     232             : 
     233             :   // Diagnostic print level
     234             :   int prtlev_;
     235             : 
     236             : };
     237             : 
     238             : }
     239             : #endif