CEMemModel.h

Go to the documentation of this file.
00001 //# CEMemModel.h: this defines CEMemModel
00002 //# Copyright (C) 1996,1997,1998,1999,2000
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //# $Id$
00027 
00028 #ifndef SYNTHESIS_CEMEMMODEL_H
00029 #define SYNTHESIS_CEMEMMODEL_H
00030 
00031 
00032 #include <casa/aips.h>
00033 #include <casa/Arrays/Matrix.h>
00034 #include <casa/Arrays/Vector.h>
00035 #include <casa/Arrays/Array.h>
00036 #include <lattices/Lattices/Lattice.h>
00037 #include <images/Images/PagedImage.h>
00038 #include <synthesis/MeasurementEquations/Entropy.h>
00039 #include <synthesis/MeasurementEquations/LinearEquation.h>
00040 #include <synthesis/MeasurementEquations/LinearModel.h>
00041 #include <casa/Logging/LogIO.h>
00042 
00043 
00044 namespace casa { //# NAMESPACE CASA - BEGIN
00045 
00046 // Forward declaration
00047 class LatConvEquation;
00048 class CEMemProgress;
00049 
00050 
00051 // <summary> Implements the Cornwell & Evans MEM Algorithm on Lattices </summary>
00052 
00053 // <use visibility=export>
00054 
00055 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
00056 // </reviewed>
00057 
00058 // <prerequisite> 
00059 // <li> ResidualEquation/ConvolutionEquation 
00060 // <li> LinearModel/LinearEquation Paradigm 
00061 // </prerequisite>
00062 //
00063 // <etymology>
00064 // This class is called CEMemModel because it uses the Cornwell and
00065 // Evans MEM algorithm to deconvolve the model. 
00066 // </etymology>
00067 //
00068 // <synopsis>
00069 // This class is used to perform the  Cornwell and Evans MEM Algorithm on an
00070 // Array. Only the deconvolved model of the sky are directly stored by this
00071 // class. The point spread function (psf) and convolved (dirty) image are
00072 // stored in a companion class which is must be derived from
00073 // ResidualEquation. 
00074 // 
00075 // The deconvolution works like this. The user constructs a CEMemModel by
00076 // specifying an initial model of the sky. This can by be
00077 // one,two,three... dimensional depending on the dimension of the psf (see
00078 // below). The user then constructs a class which implements the forward
00079 // equation between the model and the dirty image. Typically this will be
00080 // the ConvolutionEquation class, although any class which has a
00081 // ResidualEquation interface will work (but perhaps very slowly, as the
00082 // ConvolutionEquation class has member functions optimised for CLEAN and MEM)
00083 //
00084 // The user then calls the solve() function (with the appropriate equation
00085 // class as an arguement), and this class will perform the MEM deconvolution.
00086 // The various MEM parameters are set (prior to calling solve) using the
00087 // functions derived from the Iterate class, in particular
00088 // setNumberIterations().
00089 // 
00090 // The solve() function does not return either the deconvolved model or the
00091 // residuals. The solved model can be obtained using the getModel() function
00092 // (derived from ArrayModel()) and the residual can be obtained using the
00093 // residual() member function of the Convolution/Residual Equation Class.
00094 // 
00095 // The size and shape of the model used in this class MUST be the same as
00096 // the convolved data (Dirty Image), stored in the companion
00097 // ResidualEquation Class. However the model (and convolved data) can have
00098 // more dimensions than the psf, as well as a different size (either larger
00099 // or smaller). When the dimensionality is different the deconvolution is done
00100 // independendtly in each "plane" of the model. (Note this has not
00101 // been implemented yet but is relatively simple to do if necessary). 
00102 //
00103 // StokesVectors are not yet implemented.
00104 //
00105 // A companion class to this one is MaskedCEMemModel. This provides
00106 // the same functionality but is used with MaskedArrays which indicate which
00107 // regions of the NewtonRaphson (residual) image to apply when forming the
00108 // step image (MaskedCEMemModel is not yet implemented).
00109 //
00110 // </synopsis>
00111 //
00112 // <example>
00113 // <srcblock>
00114 // Matrix<Float> psf(12,12), dirty(10,10), initialModel(10,10);
00115 // ...put appropriate values into psf, dirty, & initialModel....
00116 // CEMemModel<Float> deconvolvedModel(initialModel); 
00117 // ConvolutionEquation convEqn(psf, dirty);
00118 // deconvolvedModel.setSigma(0.001); 
00119 // deconvolvedModel.setTargetFlux(-2.500); 
00120 // deconvolvedModel.setNumberIterations(30);
00121 // Bool convWorked = deconvolvedModel.solve(convEqn);
00122 // Array<Float> finalModel, residuals;
00123 // if (convWorked){
00124 //   finalModel = deconvolvedModel.getModel();
00125 //   ConvEqn.residual(deconvolvedModel, finalResidual);
00126 // }
00127 // </srcblock> 
00128 // </example>
00129 //
00130 // <motivation>
00131 // This class is needed to deconvolve extended images.  
00132 // In time, the MEM algorithm will be a principle player in the 
00133 // mosaicing stuff.
00134 // </motivation>
00135 //
00136 // <templating arg=T>
00137 //    For testing:
00138 //    <li> Float: lets try it first
00139 //    <li> StokesVector: will require lots more work
00140 // </templating>
00141 //
00142 // <todo asof="1998/12/02">
00143 //   <li> We need to implement soft Masking
00144 // </todo>
00145 
00146 class CEMemModel: public LinearModel< Lattice<Float> > 
00147 {
00148 
00149   // Any new entropies derived from Entropy must sign the friend list:
00150 friend class Entropy;
00151 friend class EntropyI;
00152 friend class EntropyEmptiness;
00153 
00154 
00155 public:
00156 
00157   // Construct the CEMemModel object and initialise the model.
00158   CEMemModel(Entropy &ent, 
00159              Lattice<Float> & model,
00160              uInt nIntegrations = 10,
00161              Float sigma = 0.001,
00162              Float targetFlux = 1.0,
00163              Bool useFluxConstraint = False,
00164              Bool initializeModel = True,
00165              Bool imagePlane = False);
00166 
00167   // Construct the CEMemModel object, initialise the model and Prior images.
00168   CEMemModel(Entropy &ent, 
00169              Lattice<Float> & model,
00170              Lattice<Float> & prior,
00171              uInt nIntegrations = 10,
00172              Float sigma = 0.001,
00173              Float targetFlux = 1.0,
00174              Bool useFluxConstraint = False,
00175              Bool initializeModel = True,
00176              Bool imagePlane = False);
00177 
00178   // Construct the CEMemModel object, initialise the model, Prior,
00179   // and mask images.
00180   CEMemModel(Entropy &ent, 
00181              Lattice<Float> & model, 
00182              Lattice<Float> & prior,
00183              Lattice<Float> & mask,
00184              uInt nIntegrations = 10,
00185              Float sigma = 0.001,
00186              Float targetFlux = 1.0,
00187              Bool useFluxConstraint = False,
00188              Bool initializeModel = True,
00189              Bool imagePlane = False);
00190 
00191 
00192   // destructor
00193   virtual ~CEMemModel();
00194 
00195   
00196   // solve the convolution equation
00197   // returns True if converged
00198 
00199   // Gives information about the state of the CEMem
00200   // (ie, using mask image, using prior image; more work here!)
00201   void state();
00202 
00203 
00204   //  This needs to be "ResidualEquation", using LatConvEquation as
00205   //  polymorphism is broken
00206   Bool solve(ResidualEquation<Lattice<Float> >  & eqn);
00207 
00208   // Set and get various state images and classes
00209   // <group>
00210   // set or get the Entropy class
00211   void setEntropy(Entropy &myEntropy ) {itsEntropy_ptr = &myEntropy;}
00212   void getEntropy(Entropy &myEntropy ) {myEntropy = *itsEntropy_ptr;}
00213 
00214   // set or get the Model image
00215   Lattice<Float>& getModel() const 
00216     { return (*(itsModel_ptr->clone())); }
00217   void setModel(const Lattice<Float> & model)
00218     { itsModel_ptr = model.clone(); }
00219 
00220   // set or get the Prior image
00221   Lattice<Float>& getPrior() const 
00222     { return (*(itsPrior_ptr->clone())); }
00223   void setPrior(Lattice<Float> & prior);
00224 
00225   // set or get the Mask image
00226   Lattice<Float>& getMask() const 
00227     { return (*(itsMask_ptr->clone())); }
00228   void setMask(Lattice<Float> & mask);
00229 
00230   // get the Residual image
00231   Lattice<Float>& getResidual() const 
00232     { return (*(itsResidual_ptr->clone())); }
00233 
00234   // </group>
00235 
00236 
00237   // set and get alpha and beta
00238   // <group>
00239   Float getAlpha() const { return itsAlpha; }
00240   Float getBeta() const { return itsBeta; }
00241   void setAlpha(Float alpha) {itsAlpha = alpha; }
00242   void setBeta(Float beta) {itsBeta = beta; }
00243   // </group>
00244 
00245   // Set various controlling parameters
00246   // (The most popular controlling parameters are 
00247   // set in the constructor)
00248   // <group>
00249   Float getTolerance() {return itsTolerance; }
00250   void  setTolerance(Float x) { itsTolerance = x; }
00251   Float getQ() {return itsQ; }
00252   void  setQ(Float x) { itsQ= x; }
00253   Float getGain() {return itsGain; }
00254   void  setGain(Float x) { itsGain = x; }
00255   Float getMaxNormGrad() {return itsMaxNormGrad; }
00256   void  setMaxNormGrad(Float x) { itsMaxNormGrad = x; }
00257   Int getInitialNumberIterations() {return itsFirstIteration; }
00258   void  setInitialNumberIterations(Int x) { itsFirstIteration = x; }
00259   // </group>
00260 
00261   // The convergence can also be in terms of the maximum residual
00262   // (ie, for use in stopping in a major cycle).
00263   void setThreshold (const Float x ) { itsThreshold0 = x; }
00264   // Thresh doubles in iter iterations
00265   void setThresholdSpeedup (const Float iter) {itsThresholdSpeedup = iter; }
00266   Float getThreshold();
00267 
00268   // Set/get the progress display 
00269   // <group>
00270   virtual void setProgress(CEMemProgress& ccp) { itsProgressPtr = &ccp; }
00271   virtual CEMemProgress& getProgress() { return *itsProgressPtr; }
00272   // </group>
00273 
00274   // return the number of iterations completed
00275   Int numberIterations () { return itsIteration; }
00276 
00277   // if this MEMModel is constructed in a MF loop, we may need to
00278   // increment the flux by the last iterations flux
00279   void setCycleFlux(Float x) { itsCycleFlux = x; }
00280   Float getCycleFlux() { return itsCycleFlux; }
00281 
00282 protected:
00283 
00284 
00285   // Perform One Iteration
00286   void oneIteration();
00287 
00288   // apply mask to a lattice; returns True if mask is available, 
00289   // False if not
00290   Bool applyMask( Lattice<Float> & lat );
00291 
00292   // Helper functions that interface with Entropy routines
00293   // Access to the entropy should be through this interface; the
00294   // points at which the Entropy is mentioned is then limited to
00295   // these lines right here, and to the constructors which set the
00296   // Entropy.  The entropy should not ever change the private
00297   // 
00298   // <group>
00299   void formEntropy() { itsEntropy = itsEntropy_ptr->formEntropy(); }
00300 
00301   void  formGDG() { itsEntropy_ptr->formGDG( itsGDG ); }
00302 
00303   void  formGDGStep() { itsEntropy_ptr->formGDGStep( itsGDG ); }
00304 
00305   void formGDS() { itsGradDotStep1=itsEntropy_ptr->formGDS(); }
00306 
00307   void entropyType(String &str) { itsEntropy_ptr->entropyType(str); }
00308 
00309   void relaxMin() { itsRequiredModelMin = itsEntropy_ptr->relaxMin(); }
00310 
00311   Bool testConvergence() { return itsEntropy_ptr->testConvergence(); }
00312 
00313   // </group>
00314 
00315 
00316   // protected generic constrcutor: DON'T USE IT!
00317   CEMemModel();
00318 
00319   // Set entropy pointer to zero: called by Entropy's destructor
00320   void letEntropyDie() { itsEntropy_ptr = 0; }
00321 
00322   // initialize itsStep and itsResidual and other stuff
00323   Bool initStuff();
00324 
00325 
00326   // controls how to change Alpha and Beta
00327   // <group>
00328   void changeAlphaBeta ();
00329   // initialize Alpha-Beta (called by changeAlphaBeta)
00330   void initializeAlphaBeta();
00331   // update Alpha-Beta (called by changeAlphaBeta)
00332   void updateAlphaBeta();
00333   // </group>
00334 
00335 
00336 
00337 
00338   // Generic utility functions
00339   // <group>
00340   // check that a single image is onf plausible shape
00341   Bool checkImage(const Lattice<Float> *);
00342   // check that the lattices and the underlying tile sizes are consistent
00343   Bool checkImages(const Lattice<Float> *one, const Lattice<Float> *other);
00344   // check that all is well in Denmark:
00345   //     ensure all images are the same size,
00346   //     ensure we have an entropy,
00347   //     ensure state variables have reasonable values
00348   Bool ok();
00349   // </group>
00350   
00351 
00352 
00353 
00354   // Helper functions for oneIteration:
00355   // <group>
00356   // calculate size of step
00357   void calculateStep();
00358 
00359   // take one step: clipped addition of
00360   // wt1*itsModel + wt2*itsStep
00361   void takeStep(Float wt1, Float wt2);
00362 
00363   // Calculate the flux, itsModMin, and itsModMax
00364   Float formFlux();
00365   // </group>
00366 
00367   // Determine if the peak residual is less than the getThreshold()
00368   Bool testConvergenceThreshold();
00369 
00370   // ------------Private Member Data---------------------
00371   // functional form of the entropy
00372   Entropy  *itsEntropy_ptr;   
00373 
00374   // form of the Residual method
00375   ResidualEquation< Lattice<Float> > * itsResidualEquation_ptr;
00376 
00377   // Images:
00378   Lattice<Float> * itsModel_ptr;
00379   Lattice<Float> * itsPrior_ptr;
00380   Lattice<Float> * itsMask_ptr;
00381   // Our OWN internal temp images; delete these upon destruction
00382   Lattice<Float> * itsStep_ptr;
00383   Lattice<Float> * itsResidual_ptr;
00384 
00385 
00386   // Controlling parameters
00387   // <group>
00388   Bool itsInitializeModel;
00389   uInt itsNumberIterations;
00390   Bool  itsDoInit;
00391   Float itsSigma;
00392   Float itsTargetFlux;  
00393   Float itsDefaultLevel;
00394   Float itsTargetChisq;
00395   // tolerance for convergence
00396   Float itsTolerance;   
00397   // N points per beam
00398   Float itsQ;           
00399   // gain for adding step image
00400   Float itsGain;        
00401   Float itsMaxNormGrad;
00402   // constrain flux or not?
00403   Bool  itsUseFluxConstraint;  
00404   // is this an image plane problem (like Single Dish or Optical?)
00405   Bool  itsDoImagePlane;
00406   Float itsThreshold0;
00407   Float itsThresholdSpeedup;
00408   Float itsCycleFlux; // flux from previous cycles
00409   // </group>
00410 
00411   // State variables
00412   // <group>  
00413   Float itsAlpha;
00414   Float itsBeta;
00415   Float itsNormGrad;
00416   Float itsFlux;
00417   Float itsTotalFlux;      // itsCycleFlux + itsFlux
00418   Float itsChisq;
00419   // sqrt( chi^2/target_chi^2 )
00420   Float itsFit;        
00421   // sqrt( chi^2/ Npixels )
00422   Float itsAFit;       
00423   // numerical value of entropy
00424   Float itsEntropy;    
00425   // Model is constrained to be >= itsNewMin
00426   Float itsRequiredModelMin; 
00427   // maximum pixel value in model
00428   Float itsModelMax;   
00429   // minimum pixel value n model
00430   Float itsModelMin;   
00431   Float itsLength;
00432   Double itsGradDotStep0;
00433   Double itsGradDotStep1;
00434   uInt   itsIteration;
00435   uInt   itsFirstIteration;
00436   // matrix of gradient dot products
00437   Matrix<Double> itsGDG;  
00438   Float itsCurrentPeakResidual;
00439   // </group>
00440   Float itsNumberPixels;
00441 
00442 
00443   // Accesories
00444   // <group>  
00445   Bool itsChoose;
00446   LogIO itsLog;
00447   // </group>  
00448 
00449   // Enumerate the different Gradient subscript types
00450   enum GradientType {
00451     // Entropy
00452     H=0,
00453     // Chi_sq
00454     C=1,
00455     // Flux
00456     F=2,
00457     // Objective function J
00458     J=3
00459   };
00460 
00461   CEMemProgress* itsProgressPtr;
00462  
00463 
00464 };
00465 
00466 
00467 
00468 } //# NAMESPACE CASA - END
00469 
00470 #endif
00471 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

Generated on 31 Aug 2016 for casa by  doxygen 1.6.1