Line data    Source code

       1             : //# VisCalSolver.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: casa-feedback@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_VISCALSOL_H
      29             : #define SYNTHESIS_VISCALSOL_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/MeasurementEquations/VisEquation.h>
      40             : #include <msvis/MSVis/VisBuffGroupAcc.h>
      41             : #include <msvis/MSVis/CalVisBuffer.h>
      42             : 
      43             : namespace casa { //# NAMESPACE CASA - BEGIN
      44             : 
      45             : // <summary> 
      46             : // VisCalSolver: Default solver for calibration using visibility data
      47             : // </summary>
      48             : 
      49             : // <use visibility=export>
      50             : 
      51             : // <reviewed reviewer="" date="" tests="" demos="">
      52             : 
      53             : // <prerequisite>
      54             : //   <li> <linkto class="MeasurementComponents">MeasurementComponents</linkto> module
      55             : //   <li> <linkto class="VisEquation">VisEquation</linkto> module
      56             : // </prerequisite>
      57             : //
      58             : // <etymology>
      59             : // VisCal for visibility calibration (meaning solved from visibilities), Solver for
      60             : // solving. 
      61             : // </etymology>
      62             : //
      63             : // <synopsis> 
      64             : //
      65             : // VisCalSolver describes an interface for solving for calibration from visibility
      66             : // data in a VisBuffer.  It hosts the communication of visibility data with otherwise
      67             : // generic solving mechanims.  A Levenberg-Marquardt solver is supplied here by
      68             : // default, but it is intended that this class be the template for interfaces to
      69             : // third-party solving mechanisms.
      70             : //
      71             : // </synopsis> 
      72             : //
      73             : // <example>
      74             : // <srcblock>
      75             : 
      76             : // </srcblock>
      77             : // </example>
      78             : //
      79             : // <motivation>
      80             : // It is desirable to establish the distinct communicative boundary between generic
      81             : // solving mechanims and the particulars of visibility data and calibration
      82             : // component descriptions.  This class is intended to serve this purpose, by providing
      83             : // access to visibility data and calibration in terms of quantities necessary for
      84             : // least-squares (and other) style solvers.  
      85             : // </motivation>
      86             : //
      87             : // <todo asof="97/10/01">
      88             : // </todo>
      89             : 
      90             : // **********************************************************
      91             : //  VisCalSolver
      92             : //
      93             : 
      94             : class VisCalSolver {
      95             : public:
      96             : 
      97             :   // Constructor currently generic
      98             :   VisCalSolver();
      99             :   
     100             :   // Destructor
     101             :   ~VisCalSolver();
     102             : 
     103             :   // Do the solve
     104             :   casacore::Bool solve(VisEquation& viseq, SolvableVisCal& svc, VisBuffer& svb);
     105             :   casacore::Bool solve(VisEquation& viseq, SolvableVisCal& svc, VisBuffGroupAcc& vbga);
     106             : 
     107             : protected:
     108             : 
     109             :   // Access to fundamental external objects:
     110           0 :   inline VisBuffer&       svb()  { return *svb_; };
     111           0 :   inline VisBuffGroupAcc& vbga() { return *vbga_; };
     112           0 :   inline VisEquation&     ve()   { return *ve_; };
     113           0 :   inline SolvableVisCal&  svc()  { return *svc_; };
     114             : 
     115             :   // Accessors to current svb's (differentiated) Residuals
     116           0 :   inline casacore::Cube<casacore::Complex>&    R()    { return R_; };
     117           0 :   inline casacore::Array<casacore::Complex>&   dR()   { return dR_; };
     118           0 :   inline casacore::Matrix<casacore::Bool>&     Rflg() { return Rflg_; };
     119             : 
     120           0 :   inline casacore::Array<casacore::Complex>&   dSrc()   { return dSrc_; };
     121             : 
     122             :   // Access to maxIter_
     123           0 :   inline casacore::Int&    maxIter() { return maxIter_; };
     124             : 
     125             :   // Access to chi2
     126           0 :   inline casacore::Double& chiSq()     { return chiSq_; };
     127           0 :   inline casacore::Vector<casacore::Double>& chiSqV()     { return chiSqV_; };
     128           0 :   inline casacore::Double& lastChiSq() { return lastChiSq_; };
     129           0 :   inline casacore::Double& dChiSq()    { return dChiSq_; };
     130           0 :   inline casacore::Double& sumWt()     { return sumWt_; };
     131           0 :   inline casacore::Int&    nWt()       { return nWt_; };
     132             : 
     133             :   // Access to parameters, & grad,hess,dp
     134           0 :   inline casacore::Int&              nTotalPar() { return nTotalPar_; };
     135           0 :   inline casacore::Int&              nCalPar()   { return nCalPar_; };
     136           0 :   inline casacore::Int&              nSrcPar()   { return nSrcPar_; };
     137           0 :   inline casacore::Vector<casacore::Complex>&  par()       { return par_; };
     138           0 :   inline casacore::Vector<casacore::Bool>&     parOK()     { return parOK_; };
     139           0 :   inline casacore::Vector<casacore::Float>&    parErr()    { return parErr_; };
     140           0 :   inline casacore::Vector<casacore::Complex>&  srcPar()    { return srcPar_; };
     141           0 :   inline casacore::Vector<casacore::DComplex>& grad()      { return grad_; };
     142           0 :   inline casacore::Vector<casacore::Double>&   hess()      { return hess_; };
     143           0 :   inline casacore::Vector<casacore::Complex>&  dpar()      { return dpar_; };
     144           0 :   inline casacore::Vector<casacore::Complex>&  dCalPar()   { return dsrcpar_; };
     145           0 :   inline casacore::Vector<casacore::Complex>&  dSrcPar()   { return dcalpar_; };
     146           0 :   inline casacore::Vector<casacore::Complex>&  lastCalPar()   { return lastCalPar_; };
     147           0 :   inline casacore::Vector<casacore::Complex>&  lastSrcPar()   { return lastSrcPar_; };
     148             : 
     149           0 :   inline casacore::Double&           lambda()    { return lambda_; };
     150             : 
     151             :   // Initialize solving data
     152             :   void initSolve();
     153             : 
     154             :   // Obtain trial residuals w.r.t svc's current pars
     155             :   void residualate();
     156             :   void residualate2();
     157             : 
     158             :   // Differentiate the svb w.r.t svc's pars
     159             :   void differentiate();
     160             :   void differentiate2();
     161             : 
     162             :   // Calculate residuals (incl. diff'd) and chi2 
     163             :   void chiSquare();
     164             :   void chiSquare2();
     165             : 
     166             :   // Check for convergence
     167             :   casacore::Bool converged();
     168             : 
     169             :   // Internal solving methods
     170             :   void accGradHess();
     171             :   void accGradHess2();
     172             :   void revert();
     173             :   void solveGradHess();
     174             :   void updatePar();
     175             : 
     176             :   // Optimize the step parabolically
     177             :   void optStepSize();
     178             :   void optStepSize2();
     179             : 
     180             :   // Get and print par errors
     181             :   void getErrors();
     182             : 
     183             :   void printPar(const casacore::Int& iter);
     184             : 
     185             : private:
     186             :   
     187             :   // Diagnostic print level
     188           0 :   inline casacore::Int& prtlev() { return prtlev_; };
     189             : 
     190             :   // VisBuffer (from outside)
     191             :   VisBuffer* svb_;
     192             :   VisBuffGroupAcc* vbga_;
     193             : 
     194             :   // VisEquation (from outside)
     195             :   VisEquation* ve_;
     196             : 
     197             :   // SVC (from outside)
     198             :   SolvableVisCal* svc_;
     199             : 
     200             :   // Total Number of parameters
     201             :   casacore::Int nTotalPar_;
     202             :   casacore::Int nCalPar_;
     203             :   casacore::Int nSrcPar_;
     204             : 
     205             :   // Residual/Differentiation caches
     206             :   casacore::Cube<casacore::Complex>  R_;
     207             :   casacore::Array<casacore::Complex> dR_;
     208             :   casacore::Matrix<casacore::Bool>    Rflg_;
     209             : 
     210             :   // Derivative wrt Q and U
     211             :   casacore::Array<casacore::Complex> dSrc_;
     212             : 
     213             :   // Maximum number of solve iterations to attempt
     214             :   casacore::Int maxIter_;
     215             : 
     216             :   // Chi2, sum wts
     217             :   casacore::Double chiSq_;
     218             :   casacore::Vector<casacore::Double> chiSqV_;
     219             :   casacore::Double lastChiSq_;
     220             :   casacore::Double dChiSq_;
     221             :   casacore::Double sumWt_;
     222             :   casacore::Int    nWt_;
     223             :   casacore::Int    cvrgcount_;
     224             : 
     225             :   // Parameter storage
     226             :   // (these are casacore::Complex to match the VisCal solvePar)
     227             :   casacore::Vector<casacore::Complex> par_;
     228             :   casacore::Vector<casacore::Bool>    parOK_;
     229             :   casacore::Vector<casacore::Float>   parErr_;
     230             :   casacore::Vector<casacore::Complex> srcPar_;  
     231             :   casacore::Vector<casacore::Complex> lastCalPar_,lastSrcPar_;
     232             : 
     233             :   // Parameter update
     234             :   casacore::Vector<casacore::Complex> dpar_;
     235             :   casacore::Vector<casacore::Complex> dcalpar_,dsrcpar_;
     236             : 
     237             :   // Gradient, Hessian
     238             :   //  (these are casacore::Double for precision in accumulation
     239             :   casacore::Vector<casacore::DComplex> grad_;
     240             :   casacore::Vector<casacore::Double> hess_;
     241             : 
     242             :   // LM factor
     243             :   casacore::Double lambda_;
     244             : 
     245             :   // Step optimization toggle
     246             :   casacore::Bool optstep_;
     247             : 
     248             :   // Diagnostic print level
     249             :   casacore::Int prtlev_;
     250             : 
     251             : };
     252             : 
     253             : }
     254             : #endif