Line data Source code
1 : //# SDAlgorithmMEM.cc: Implementation of SDAlgorithmMEM 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 :
31 : #include <components/ComponentModels/SkyComponent.h>
32 : #include <components/ComponentModels/ComponentList.h>
33 : #include <casacore/images/Images/TempImage.h>
34 : #include <casacore/images/Images/SubImage.h>
35 : #include <casacore/images/Regions/ImageRegion.h>
36 : #include <casacore/casa/OS/File.h>
37 : #include <casacore/lattices/LEL/LatticeExpr.h>
38 : #include <casacore/lattices/Lattices/TiledLineStepper.h>
39 : #include <casacore/lattices/Lattices/LatticeStepper.h>
40 : #include <casacore/lattices/Lattices/LatticeIterator.h>
41 : #include <synthesis/TransformMachines/StokesImageUtil.h>
42 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
43 : #include <casacore/casa/Exceptions/Error.h>
44 : #include <casacore/casa/BasicSL/String.h>
45 : #include <casacore/casa/Utilities/Assert.h>
46 : #include <casacore/casa/OS/Directory.h>
47 : #include <casacore/tables/Tables/TableLock.h>
48 :
49 : #include<synthesis/ImagerObjects/SIMinorCycleController.h>
50 :
51 : #include <sstream>
52 :
53 : #include <casacore/casa/Logging/LogMessage.h>
54 : #include <casacore/casa/Logging/LogIO.h>
55 : #include <casacore/casa/Logging/LogSink.h>
56 :
57 : #include <casacore/casa/System/Choice.h>
58 : #include <msvis/MSVis/StokesVector.h>
59 :
60 : #include <synthesis/ImagerObjects/SDAlgorithmMEM.h>
61 : #include <synthesis/MeasurementEquations/CEMemModel.h>
62 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
63 : #include <synthesis/MeasurementEquations/IPLatConvEquation.h>
64 :
65 :
66 : using namespace casacore;
67 : namespace casa { //# NAMESPACE CASA - BEGIN
68 :
69 0 : SDAlgorithmMEM::SDAlgorithmMEM(String entropy):
70 0 : SDAlgorithmBase()
71 : {
72 0 : LogIO os( LogOrigin("SDAlgorithmMEM","Constructor",WHERE) );
73 0 : itsAlgorithmName=String("Mem");
74 0 : itsMatDeltaModel.resize();
75 :
76 0 : if(entropy=="entropy") {
77 : // os << "Deconvolving image using maximum entropy algorithm" << LogIO::POST;
78 0 : itsEnt = new EntropyI;
79 : }
80 0 : else if (entropy=="emptiness") {
81 0 : itsEnt = new EntropyEmptiness;
82 : }
83 : else {
84 : // Put check in Deconvolver parameter object.
85 0 : os << " Known MEM entropies: entropy | emptiness " << LogIO::POST;
86 0 : os << LogIO::SEVERE << "Unknown MEM entropy: " << entropy << LogIO::POST;
87 : // return false;
88 : }
89 :
90 0 : }
91 :
92 0 : SDAlgorithmMEM::~SDAlgorithmMEM()
93 : {
94 :
95 0 : }
96 :
97 : // void SDAlgorithmMEM::initializeDeconvolver( Float &peakresidual, Float &modelflux )
98 0 : void SDAlgorithmMEM::initializeDeconvolver()
99 : {
100 0 : LogIO os( LogOrigin("SDAlgorithmMEM","initializeDeconvolver",WHERE) );
101 :
102 0 : itsImages->residual()->get( itsMatResidual, true );
103 0 : itsImages->model()->get( itsMatModel, true );
104 0 : itsImages->psf()->get( itsMatPsf, true );
105 0 : itsImages->mask()->get( itsMatMask, true );
106 :
107 :
108 : // cout << "initDecon : " << itsImages->residual()->shape() << " : " << itsMatResidual.shape()
109 : // << itsImages->model()->shape() << " : " << itsMatModel.shape()
110 : // << itsImages->psf()->shape() << " : " << itsMatPsf.shape()
111 : // << endl;
112 :
113 : /*
114 : findMaxAbs( itsMatResidual, itsPeakResidual, itsMaxPos );
115 : itsModelFlux = sum( itsMatModel );
116 :
117 : peakresidual = itsPeakResidual;
118 : modelflux = itsModelFlux;
119 : */
120 :
121 : // Initialize the Delta Image model. Resize if needed.
122 0 : if ( itsMatDeltaModel.shape().nelements() != itsMatModel.shape().nelements() )
123 0 : { itsMatDeltaModel.resize ( itsMatModel.shape() ); }
124 :
125 :
126 0 : }
127 :
128 : // Code obtained from Deconvolver.cc
129 0 : void SDAlgorithmMEM::takeOneStep( Float /*loopgain*/,
130 : Int cycleNiter,
131 : Float cycleThreshold,
132 : Float &peakresidual,
133 : Float &modelflux,
134 : Int &iterdone)
135 : {
136 0 : LogIO os( LogOrigin("SDAlgorithmMEM","takeOneStep",WHERE) );
137 : // tmp
138 0 : itsImages->residual()->get( itsMatResidual, true );
139 0 : itsImages->mask()->get( itsMatMask, true );
140 0 : findMaxAbsMask( itsMatResidual, itsMatMask, itsPeakResidual, itsMaxPos );
141 0 : cout << "Peak Res at start of step : " << itsPeakResidual << endl;
142 : // tmp
143 :
144 :
145 : // Store current model in this matrix.
146 0 : itsImages->model()->get( itsMatDeltaModel, true );
147 0 : itsMatModel.assign( itsMatDeltaModel ); // This should make an explicit copy
148 0 : cout << "Flux at start of step : " << sum(itsMatModel) << endl;
149 :
150 : // Set model to zero
151 :
152 0 : LatticeLocker lockmod(*(itsImages->model()), FileLocker::Write);
153 0 : itsImages->model()->set( 0.0 );
154 :
155 : // Add to construction params
156 0 : Float targetFlux=1.0;
157 0 : Bool constrainTargetFlux=false;
158 0 : Bool initializeModel=true; // Use incremental model ?
159 0 : Bool imagePlane=true; // Use full image plane. Otherwise, use inner quarter.
160 :
161 0 : CEMemModel memer( *itsEnt, *(itsImages->model()), cycleNiter, cycleThreshold,
162 : targetFlux, constrainTargetFlux,
163 0 : initializeModel, imagePlane);
164 :
165 0 : if (!initializeModel) {
166 0 : Record info=itsImages->model()->miscInfo();
167 : try {
168 0 : Float alpha = 0.0;
169 0 : Float beta = 0.0;
170 0 : info.get("ALPHA", alpha);
171 0 : memer.setAlpha(alpha);
172 0 : info.get("BETA", beta);
173 0 : memer.setBeta(beta);
174 0 : } catch (AipsError x) {
175 : // could not get Alpha and Beta for initialization
176 : // continue
177 : os << "Could not retrieve Alpha and Beta from previously initialized model"
178 0 : << LogIO::POST;
179 0 : }
180 0 : }
181 :
182 : /// Set the Prior
183 : /// if(prior != 0){
184 : /// memer.setPrior(priorSub);
185 : /// }
186 :
187 0 : memer.setMask(*(itsImages->mask()));
188 :
189 0 : CountedPtr<ResidualEquation<Lattice<Float> > > residEqn;
190 :
191 0 : if (imagePlane) {
192 0 : residEqn = new IPLatConvEquation (*(itsImages->psf()), *(itsImages->residual()));
193 : } else {
194 0 : residEqn = new LatConvEquation (*(itsImages->psf()), *(itsImages->residual()));
195 : }
196 :
197 : //Bool result=
198 0 : memer.solve(*residEqn);
199 :
200 0 : Record info=itsImages->model()->miscInfo();
201 0 : info.define("ALPHA", memer.getAlpha());
202 0 : info.define("BETA", memer.getBeta());
203 0 : itsImages->model()->setMiscInfo(info);
204 :
205 0 : iterdone = memer.numberIterations();
206 :
207 0 : LatticeExprNode maxres( max( memer.getResidual() ) );
208 0 : cout << "MAX RES at end : " << maxres.getFloat() << endl;
209 :
210 : // Retrieve residual before major cycle
211 0 : itsImages->residual()->copyData( memer.getResidual() );
212 :
213 : // Add delta model to old model
214 : //Bool ret2 =
215 0 : itsImages->model()->get( itsMatDeltaModel, true );
216 0 : itsMatModel += itsMatDeltaModel;
217 :
218 : //---------------------------------
219 :
220 : // Find Peak Residual
221 0 : itsImages->residual()->get( itsMatResidual, true );
222 0 : itsImages->mask()->get( itsMatMask, true );
223 0 : findMaxAbsMask( itsMatResidual, itsMatMask, itsPeakResidual, itsMaxPos );
224 0 : peakresidual = itsPeakResidual;
225 :
226 : // Find Total Model flux
227 0 : modelflux = sum( itsMatModel ); // Performance hog ?
228 0 : (itsImages->model())->put( itsMatModel );
229 :
230 0 : cout << "peakres : " << peakresidual << " model : " << modelflux << endl;
231 :
232 0 : }
233 :
234 0 : void SDAlgorithmMEM::finalizeDeconvolver()
235 : {
236 0 : LatticeLocker lock1(*(itsImages->residual()), FileLocker::Write);
237 0 : LatticeLocker lock2(*(itsImages->model()), FileLocker::Write);
238 0 : (itsImages->residual())->put( itsMatResidual );
239 0 : (itsImages->model())->put( itsMatModel );
240 0 : }
241 :
242 :
243 : } //# NAMESPACE CASA - END
244 :
|
