Line data Source code
1 : //# SDAlgorithmMSClean.cc: Implementation of SDAlgorithmMSClean classes
2 : //# Copyright (C) 1996,1997,1998,1999,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 addressed 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 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/OS/HostInfo.h>
30 : #include <synthesis/ImagerObjects/SDAlgorithmMSClean.h>
31 :
32 : #include <components/ComponentModels/SkyComponent.h>
33 : #include <components/ComponentModels/ComponentList.h>
34 : #include <casacore/images/Images/TempImage.h>
35 : #include <casacore/images/Images/SubImage.h>
36 : #include <casacore/images/Regions/ImageRegion.h>
37 : #include <casacore/casa/OS/File.h>
38 : #include <casacore/lattices/LEL/LatticeExpr.h>
39 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
40 : #include <casacore/lattices/Lattices/LatticeStepper.h>
41 : #include <casacore/lattices/Lattices/LatticeIterator.h>
42 : #include <casacore/lattices/Lattices/LatticeLocker.h>
43 : #include <synthesis/TransformMachines/StokesImageUtil.h>
44 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
45 : #include <casacore/casa/Exceptions/Error.h>
46 : #include <casacore/casa/BasicSL/String.h>
47 : #include <casacore/casa/Utilities/Assert.h>
48 : #include <casacore/casa/OS/Directory.h>
49 : #include <casacore/tables/Tables/TableLock.h>
50 :
51 : #include<synthesis/ImagerObjects/SIMinorCycleController.h>
52 :
53 : #include <sstream>
54 :
55 : #include <casacore/casa/Logging/LogMessage.h>
56 : #include <casacore/casa/Logging/LogIO.h>
57 : #include <casacore/casa/Logging/LogSink.h>
58 :
59 : #include <casacore/casa/System/Choice.h>
60 : #include <msvis/MSVis/StokesVector.h>
61 :
62 :
63 : using namespace casacore;
64 : namespace casa { //# NAMESPACE CASA - BEGIN
65 :
66 :
67 34 : SDAlgorithmMSClean::SDAlgorithmMSClean( Vector<Float> scalesizes,
68 : Float smallscalebias,
69 : // Int stoplargenegatives,
70 34 : Int stoppointmode ):
71 : SDAlgorithmBase(),
72 34 : itsMatPsf(), itsMatResidual(), itsMatModel(),
73 34 : itsCleaner(),
74 34 : itsScaleSizes(scalesizes),
75 34 : itsSmallScaleBias(smallscalebias),
76 : // itsStopLargeNegatives(stoplargenegatives),
77 68 : itsStopPointMode(stoppointmode)
78 : {
79 34 : itsAlgorithmName=String("multiscale");
80 34 : if( itsScaleSizes.nelements()==0 ){ itsScaleSizes.resize(1); itsScaleSizes[0]=0.0; }
81 34 : }
82 :
83 :
84 68 : SDAlgorithmMSClean::~SDAlgorithmMSClean()
85 : {
86 :
87 68 : }
88 :
89 13 : Long SDAlgorithmMSClean::estimateRAM(const vector<int>& imsize){
90 13 : Long mem=0;
91 13 : IPosition shp;
92 13 : if(itsImages)
93 0 : shp=itsImages->getShape();
94 13 : else if(imsize.size() >1)
95 13 : shp=IPosition(imsize);
96 : else
97 0 : return 0;
98 : //throw(AipsError("Deconvolver cannot estimate the memory usage at this point"));
99 :
100 : //Number of planes in memory
101 : //npsf=nscales+1
102 : //nresidual=4+nscales
103 : //nmodel=1
104 : //nmasks=2+nscales
105 : //transfer functions=nscales*2
106 13 : Long nplanes=5*itsScaleSizes.nelements()+7;
107 : //in kB
108 13 : mem=sizeof(Float)*(shp(0))*(shp(1))*nplanes/1024;
109 :
110 13 : return mem;
111 13 : }
112 :
113 : // void SDAlgorithmMSClean::initializeDeconvolver( Float &peakresidual, Float &modelflux )
114 44 : void SDAlgorithmMSClean::initializeDeconvolver()
115 : {
116 88 : LogIO os( LogOrigin("SDAlgorithmMSClean","initializeDeconvolver",WHERE) );
117 :
118 44 : AlwaysAssert( (bool) itsImages, AipsError );
119 : {
120 44 : LatticeLocker lock1 (*(itsImages->residual()), FileLocker::Read);
121 44 : LatticeLocker lock2 (*(itsImages->model()), FileLocker::Read);
122 44 : LatticeLocker lock3 (*(itsImages->psf()), FileLocker::Read);
123 44 : LatticeLocker lock4 (*(itsImages->mask()), FileLocker::Read);
124 44 : (itsImages->residual())->get( itsMatResidual , true );
125 44 : (itsImages->model())->get( itsMatModel , true );
126 44 : (itsImages->psf())->get( itsMatPsf , true );
127 44 : itsImages->mask()->get( itsMatMask, true );
128 :
129 44 : }
130 : //// Initialize the MatrixCleaner.
131 : /// ----------- do once ----------
132 : {
133 44 : itsCleaner.defineScales( itsScaleSizes );
134 :
135 44 : if(itsSmallScaleBias > 1)
136 : {
137 0 : os << LogIO::WARN << "Acceptable smallscalebias values are [-1,1].Changing smallscalebias from " << itsSmallScaleBias <<" to 1." << LogIO::POST;
138 0 : itsSmallScaleBias = 1;
139 : }
140 :
141 44 : if(itsSmallScaleBias < -1)
142 : {
143 0 : os << LogIO::WARN << "Acceptable smallscalebias values are [-1,1].Changing smallscalebias from " << itsSmallScaleBias <<" to -1." << LogIO::POST;
144 0 : itsSmallScaleBias = -1;
145 : }
146 :
147 :
148 44 : itsCleaner.setSmallScaleBias( itsSmallScaleBias );
149 : //itsCleaner.stopAtLargeScaleNegative( itsStopLargeNegatives );// In MFMSCleanImageSkyModel.cc, this is only for the first two major cycles...
150 44 : itsCleaner.stopPointMode( itsStopPointMode );
151 44 : itsCleaner.ignoreCenterBox( true ); // Clean full image
152 :
153 44 : Matrix<Float> tempMat;
154 44 : tempMat.reference( itsMatPsf );
155 44 : itsCleaner.setPsf( tempMat );
156 44 : itsCleaner.makePsfScales();
157 44 : }
158 : /// -----------------------------------------
159 :
160 : /*
161 : /// Find initial max vals..
162 : findMaxAbsMask( itsMatResidual, itsMatMask, itsPeakResidual, itsMaxPos );
163 : itsModelFlux = sum( itsMatModel );
164 :
165 : peakresidual = itsPeakResidual;
166 : modelflux = itsModelFlux;
167 : */
168 :
169 : // Parts to be repeated at each minor cycle start....
170 :
171 44 : itsCleaner.setcontrol(CleanEnums::MULTISCALE,0,0,0);/// Needs to come before makeDirtyScales
172 :
173 44 : Matrix<Float> tempmask(itsMatMask);
174 44 : itsCleaner.setMask( tempmask );
175 :
176 44 : Matrix<Float> tempMat1;
177 44 : tempMat1.reference( itsMatResidual );
178 44 : itsCleaner.setDirty( tempMat1 );
179 44 : itsCleaner.makeDirtyScales();
180 :
181 44 : }
182 :
183 44 : void SDAlgorithmMSClean::takeOneStep( Float loopgain, Int cycleNiter, Float cycleThreshold, Float &peakresidual, Float &modelflux, Int &iterdone)
184 : {
185 88 : LogIO os( LogOrigin("SDAlgorithmMSClean","takeOneStep",WHERE) );
186 :
187 44 : Quantity thresh( cycleThreshold, "Jy" );
188 : // Quantity ftthresh( 100.0, "Jy" ); /// Look at MFMSCleanImageSkyModel.cc for more.
189 44 : itsCleaner.setcontrol(CleanEnums::MULTISCALE, cycleNiter, loopgain, thresh); //, ftthresh);
190 :
191 44 : Matrix<Float> tempModel;
192 44 : tempModel.reference( itsMatModel );
193 : //save the previous model
194 44 : Matrix<Float> prevModel;
195 44 : prevModel=itsMatModel;
196 :
197 : //cout << "SDALMS, matrix shape : " << tempModel.shape() << " array shape : " << itsMatModel.shape() << endl;
198 :
199 : // retval
200 : // 1 = converged
201 : // 0 = not converged but behaving normally
202 : // -1 = not converged and stopped on cleaning consecutive smallest scale
203 : // -2 = not converged and either large scale hit negative or diverging
204 : // -3 = clean is diverging rather than converging
205 44 : itsCleaner.startingIteration( 0 );
206 44 : Int retval = itsCleaner.clean( tempModel );
207 44 : iterdone = itsCleaner.numberIterations();
208 :
209 44 : if( retval==-1 ) {os << LogIO::WARN << "MSClean minor cycle stopped on cleaning consecutive smallest scale" << LogIO::POST; }
210 44 : if( retval==-2 ) {os << LogIO::WARN << "MSClean minor cycle stopped at large scale negative or diverging" << LogIO::POST;}
211 44 : if( retval==-3 ) {os << LogIO::WARN << "MSClean minor cycle stopped because it is diverging" << LogIO::POST; }
212 :
213 : // Retrieve residual to be saved to the .residual file in finalizeDeconvolver.
214 : ////This is going to be wrong if there is no 0 scale;
215 : ///Matrix<Float> residual(itsCleaner.residual());
216 : //Matrix<Float> residual(itsCleaner.residual(tempModel-prevModel));
217 : // cout << "Max tempModel : " << max(abs(tempModel)) << " Max prevModel : " << max(abs(prevModel)) << endl;
218 44 : itsMatResidual = itsCleaner.residual(tempModel-prevModel);
219 :
220 : // account for mask as well
221 : //peakresidual = max(abs(residual*itsMatMask));
222 44 : peakresidual = max(abs(itsMatResidual*itsMatMask));
223 44 : modelflux = sum( itsMatModel ); // Performance hog ?
224 44 : }
225 :
226 44 : void SDAlgorithmMSClean::finalizeDeconvolver()
227 : {
228 : ///MatrixCleaner does not modify the original residual image matrix
229 : ///so the first line is a dummy.
230 44 : LatticeLocker lock1 (*(itsImages->residual()), FileLocker::Write);
231 44 : LatticeLocker lock2 (*(itsImages->model()), FileLocker::Write);
232 44 : (itsImages->residual())->put( itsMatResidual );
233 44 : (itsImages->model())->put( itsMatModel );
234 44 : }
235 :
236 :
237 : } //# NAMESPACE CASA - END
238 :
|
