Line data Source code
1 : //# SDAlgorithmAAspClean.cc: Implementation of SDAlgorithmAAspClean 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: 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 : //# $Id$
27 :
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/OS/HostInfo.h>
30 : #include <synthesis/ImagerObjects/SDAlgorithmAAspClean.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 <synthesis/TransformMachines/StokesImageUtil.h>
43 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
44 : #include <casacore/casa/Exceptions/Error.h>
45 : #include <casacore/casa/BasicSL/String.h>
46 : #include <casacore/casa/Utilities/Assert.h>
47 : #include <casacore/casa/OS/Directory.h>
48 : #include <casacore/tables/Tables/TableLock.h>
49 :
50 : #include<synthesis/ImagerObjects/SIMinorCycleController.h>
51 :
52 : #include <sstream>
53 :
54 : #include <casacore/casa/Logging/LogMessage.h>
55 : #include <casacore/casa/Logging/LogIO.h>
56 : #include <casacore/casa/Logging/LogSink.h>
57 :
58 : #include <casacore/casa/System/Choice.h>
59 : #include <msvis/MSVis/StokesVector.h>
60 :
61 : using namespace casacore;
62 : namespace casa { //# NAMESPACE CASA - BEGIN
63 :
64 0 : SDAlgorithmAAspClean::SDAlgorithmAAspClean(Float fusedThreshold, bool isSingle, Int largestScale, Int stoppointmode):
65 : SDAlgorithmBase(),
66 : itsMatPsf(), itsMatResidual(), itsMatModel(),
67 : itsCleaner(),
68 : itsStopPointMode(stoppointmode),
69 : itsMCsetup(true),
70 : itsFusedThreshold(fusedThreshold),
71 : itsPrevPsfWidth(0),
72 : itsIsSingle(isSingle),
73 0 : itsUserLargestScale(largestScale)
74 : {
75 0 : itsAlgorithmName = String("asp");
76 0 : }
77 :
78 0 : SDAlgorithmAAspClean::~SDAlgorithmAAspClean()
79 : {
80 :
81 0 : }
82 :
83 0 : void SDAlgorithmAAspClean::initializeDeconvolver()
84 : {
85 0 : LogIO os(LogOrigin("SDAlgorithmAAspClean", "initializeDeconvolver", WHERE));
86 0 : AlwaysAssert((bool)itsImages, AipsError);
87 :
88 0 : itsImages->residual()->get( itsMatResidual, true );
89 0 : itsImages->model()->get( itsMatModel, true );
90 0 : itsImages->psf()->get( itsMatPsf, true );
91 0 : itsImages->mask()->get( itsMatMask, true );
92 :
93 : // Initialize the AspMatrixCleaner.
94 : // If it's single channel, this only needs to be computed once.
95 : // Otherwise, it needs to be called repeatedly at each minor cycle start to
96 : // get psf for each channel
97 0 : if (itsMCsetup)
98 : {
99 0 : Matrix<Float> tempMat(itsMatPsf);
100 0 : itsCleaner.setPsf(tempMat);
101 : // Initial scales are unchanged and only need to be
102 : // computed when psf width is updated
103 0 : const Float width = itsCleaner.getPsfGaussianWidth(*(itsImages->psf()));
104 : //if user does not provide the largest scale, we calculate it internally.
105 0 : itsCleaner.setUserLargestScale(itsUserLargestScale);
106 : // we do not use the shortest baseline approach below b/c of reasons in CAS-940 dated around Jan 2022
107 : // itsCleaner.getLargestScaleSize(*(itsImages->psf()));
108 :
109 0 : if (itsPrevPsfWidth != width)
110 : {
111 0 : itsPrevPsfWidth = width;
112 0 : itsCleaner.setInitScaleXfrs(width);
113 : }
114 :
115 0 : itsCleaner.stopPointMode( itsStopPointMode );
116 0 : itsCleaner.ignoreCenterBox( true ); // Clean full image
117 : // If it's single channel, we do not do the expensive set up repeatedly
118 0 : if (itsIsSingle)
119 0 : itsMCsetup = false;
120 : // Not used. Kept for unit test
121 : //Matrix<Float> tempMat1(itsMatResidual);
122 : //itsCleaner.setOrigDirty( tempMat1 );
123 :
124 0 : if (itsFusedThreshold < 0)
125 : {
126 0 : os << LogIO::WARN << "Acceptable fusedthreshld values are >= 0. Changing fusedthreshold from " << itsFusedThreshold << " to 0." << LogIO::POST;
127 0 : itsFusedThreshold = 0.0;
128 : }
129 :
130 0 : itsCleaner.setFusedThreshold(itsFusedThreshold);
131 : }
132 :
133 : // Parts to be repeated at each minor cycle start....
134 0 : itsCleaner.setInitScaleMasks(itsMatMask);
135 0 : itsCleaner.setaspcontrol(0, 0, 0, Quantity(0.0, "%"));/// Needs to come before the rest
136 :
137 0 : Matrix<Float> tempMat1;
138 0 : tempMat1.reference(itsMatResidual);
139 0 : itsCleaner.setDirty( tempMat1 );
140 : // InitScaleXfrs and InitScaleMasks should already be set
141 0 : itsScaleSizes.clear();
142 0 : itsScaleSizes = itsCleaner.getActiveSetAspen();
143 0 : itsScaleSizes.push_back(0.0); // put 0 scale
144 0 : itsCleaner.defineAspScales(itsScaleSizes);
145 0 : }
146 :
147 :
148 0 : void SDAlgorithmAAspClean::takeOneStep( Float loopgain,
149 : Int cycleNiter,
150 : Float cycleThreshold,
151 : Float &peakresidual,
152 : Float &modelflux,
153 : Int &iterdone)
154 : {
155 0 : LogIO os( LogOrigin("SDAlgorithmAAspClean","takeOneStep", WHERE) );
156 :
157 0 : Quantity thresh(cycleThreshold, "Jy");
158 0 : itsCleaner.setaspcontrol(cycleNiter, loopgain, thresh, Quantity(0.0, "%"));
159 0 : Matrix<Float> tempModel;
160 0 : tempModel.reference( itsMatModel );
161 : //save the previous model
162 0 : Matrix<Float> prevModel;
163 0 : prevModel = itsMatModel;
164 :
165 : //cout << "AAspALMS, matrix shape : " << tempModel.shape() << " array shape : " << itsMatModel.shape() << endl;
166 :
167 : // retval
168 : // 1 = converged
169 : // 0 = not converged but behaving normally
170 : // -1 = not converged and stopped on cleaning consecutive smallest scale
171 : // -2 = not converged and either large scale hit negative or diverging
172 : // -3 = clean is diverging rather than converging
173 0 : itsCleaner.startingIteration( 0 );
174 0 : Int retval = itsCleaner.aspclean( tempModel );
175 0 : iterdone = itsCleaner.numberIterations();
176 :
177 0 : if( retval==-1 ) {os << LogIO::WARN << "AspClean minor cycle stopped on cleaning consecutive smallest scale" << LogIO::POST; }
178 0 : if( retval==-2 ) {os << LogIO::WARN << "AspClean minor cycle stopped at large scale negative or diverging" << LogIO::POST;}
179 0 : if( retval==-3 ) {os << LogIO::WARN << "AspClean minor cycle stopped because it is diverging" << LogIO::POST; }
180 :
181 : // update residual - this is critical
182 0 : itsMatResidual = itsCleaner.getterResidual();
183 :
184 0 : peakresidual = itsCleaner.getterPeakResidual();
185 : //cout << "SDAlg: peakres " << peakresidual << endl;
186 0 : modelflux = sum( itsMatModel );
187 0 : }
188 :
189 0 : void SDAlgorithmAAspClean::finalizeDeconvolver()
190 : {
191 0 : (itsImages->residual())->put( itsMatResidual );
192 0 : (itsImages->model())->put( itsMatModel );
193 0 : }
194 :
195 : } //# NAMESPACE CASA - END
|
