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