Line data Source code
1 : //# NNLSImageSkyModel.cc: Implementation of NNLSImageSkyModel class
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,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/Matrix.h>
29 : #include <casacore/casa/Arrays/Cube.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/ArrayPosIter.h>
32 : #include <synthesis/MeasurementComponents/NNLSImageSkyModel.h>
33 : #include <casacore/images/Images/PagedImage.h>
34 : #include <casacore/casa/OS/File.h>
35 : #include <casacore/casa/OS/HostInfo.h>
36 : #include <casacore/lattices/Lattices/LatticeStepper.h>
37 : #include <casacore/lattices/Lattices/LatticeIterator.h>
38 : #include <synthesis/MeasurementEquations/SkyEquation.h>
39 : #include <casacore/scimath/Mathematics/NNLSMatrixSolver.h>
40 : #include <casacore/casa/Exceptions/Error.h>
41 : #include <casacore/casa/BasicSL/String.h>
42 : #include <casacore/casa/Utilities/Assert.h>
43 :
44 : #include <sstream>
45 :
46 : #include <casacore/casa/Logging/LogMessage.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 :
49 : using namespace casacore;
50 : namespace casa { //# NAMESPACE CASA - BEGIN
51 :
52 : // NNLS solver: This could be a whole lot smarter about memory use!
53 0 : Bool NNLSImageSkyModel::solve(SkyEquation& se) {
54 :
55 0 : LogMessage message(LogOrigin("NNLSImageSkyModel","solve"));
56 :
57 0 : if(numberOfModels()>1) {
58 0 : message.message("Cannot process more than one field");
59 0 : logSink().post(message);
60 0 : return false;
61 : }
62 :
63 : // Zero Stokes I for the masked region
64 0 : maskedZeroI();
65 :
66 : // Make the dirty image
67 0 : makeNewtonRaphsonStep(se);
68 :
69 : //Make the PSF
70 0 : makeApproxPSFs(se);
71 :
72 0 : if(isSolveable(0)) {
73 :
74 0 : Int nx=image(0).shape()(0);
75 0 : Int ny=image(0).shape()(1);
76 0 : Int npol=image(0).shape()(2);
77 0 : Int nchan=image(0).shape()(3);
78 :
79 0 : NNLSMatrixSolver matrixSolver; // Matrix solver to be used
80 0 : matrixSolver.setGain(gain());
81 0 : matrixSolver.setMaxIters(numberIterations());
82 0 : matrixSolver.setTolerance(tolerance());
83 :
84 : // Loop over all channels
85 0 : LatticeStepper psfls(PSF(0).shape(), IPosition(4,nx,ny,1,1),
86 0 : IPosition(4,0,1,2,3));
87 0 : RO_LatticeIterator<Float> psfli(PSF(0),psfls);
88 0 : LatticeStepper ls(image(0).shape(), IPosition(4, nx, ny, npol, 1),
89 0 : IPosition(4, 0, 1, 2, 3));
90 0 : LatticeIterator<Float> imageli(image(0), ls);
91 :
92 : // Get I plane of image
93 0 : Matrix<Float> limage;
94 0 : if(npol>1) {
95 0 : limage=imageli.cubeCursor().xyPlane(0);
96 : }
97 : else {
98 0 : limage=imageli.matrixCursor();
99 : }
100 0 : LatticeIterator<Float> imageStepli(residual(0), ls);
101 0 : Matrix<Float> limageStep;
102 0 : if(npol>1) {
103 0 : limageStep=imageStepli.cubeCursor().xyPlane(0);
104 : }
105 : else {
106 0 : limageStep=imageStepli.matrixCursor();
107 : }
108 : // Get the flux mask?
109 0 : ArrayPositionIterator fai(limageStep.shape(),0);
110 0 : Int lFluxMask=nx*ny;
111 0 : Matrix<Float> lfluxmask;
112 0 : if(hasFluxMask(0)) {
113 0 : AlwaysAssert(fluxMask(0).shape()(0)==nx, AipsError);
114 0 : AlwaysAssert(fluxMask(0).shape()(1)==ny, AipsError);
115 0 : LatticeStepper mls(fluxMask(0).shape(),
116 0 : IPosition(4, nx, ny, npol, 1),
117 0 : IPosition(4, 0, 1, 2, 3));
118 0 : RO_LatticeIterator<Float> fluxmaskli(fluxMask(0), mls);
119 0 : fluxmaskli.reset();
120 0 : if(npol>1) {
121 0 : lfluxmask=fluxmaskli.cubeCursor().xyPlane(0);
122 : }
123 : else {
124 0 : lfluxmask=fluxmaskli.matrixCursor();
125 : }
126 0 : lFluxMask=0;
127 0 : ArrayPositionIterator fai(limage.shape(),0);
128 0 : for(fai.origin();!fai.pastEnd();fai.next()) {
129 0 : if(lfluxmask(fai.pos())>0.0) lFluxMask++;
130 : }
131 0 : AlwaysAssert(lFluxMask>0, AipsError);
132 : }
133 : else {
134 0 : lfluxmask.resize(nx,ny);
135 0 : lfluxmask=1.0;
136 : }
137 : // Get the data mask?
138 0 : Int lMask=nx*ny;
139 0 : Matrix<Float> lmask;
140 0 : ArrayPositionIterator dai(limageStep.shape(),0);
141 0 : if(hasMask(0)) {
142 0 : AlwaysAssert(mask(0).shape()(0)==nx, AipsError);
143 0 : AlwaysAssert(mask(0).shape()(1)==ny, AipsError);
144 0 : LatticeStepper mls(mask(0).shape(),
145 0 : IPosition(4, nx, ny, npol, 1),
146 0 : IPosition(4, 0, 1, 2, 3));
147 0 : RO_LatticeIterator<Float> maskli(mask(0), mls); maskli.reset();
148 0 : if(npol>1) {
149 0 : lmask=maskli.cubeCursor().xyPlane(0);
150 : }
151 : else {
152 0 : lmask=maskli.matrixCursor();
153 : }
154 0 : lMask=0;
155 0 : ArrayPositionIterator dai(limageStep.shape(),0);
156 0 : for(dai.origin();!dai.pastEnd();dai.next()) {
157 0 : if(lmask(dai.pos())>0.0) lMask++;
158 : }
159 0 : AlwaysAssert(lMask>0, AipsError);
160 : }
161 : else {
162 0 : lmask.resize(nx,ny);
163 0 : lmask=1.0;
164 : }
165 :
166 : {
167 0 : ostringstream o; o <<""<<lMask<<" constraints on "<<lFluxMask
168 0 : <<" free pixels";message.message(o);
169 0 : logSink().post(message);
170 : }
171 :
172 0 : if(4.0*Double(lMask)*Double(lFluxMask)>Double(HostInfo::memoryTotal(true))*1024.0) {
173 0 : ostringstream o;
174 0 : o << "Insufficient memory for PSF matrix: reduce the size of the masks";
175 0 : message.message(o);
176 0 : logSink().post(message);
177 0 : return false;
178 : }
179 0 : Matrix<Float> AMatrix(lMask, lFluxMask);
180 0 : Vector<Float> XVector(lFluxMask); // Unknown X vector
181 0 : Vector<Float> BVector(lMask); // Data Constraint Vector AX=B
182 :
183 0 : Int chan=0;
184 0 : for (imageStepli.reset(),imageli.reset(),psfli.reset();
185 0 : !imageStepli.atEnd();
186 0 : imageStepli++,imageli++,psfli++,chan++) {
187 : {
188 0 : ostringstream o; o<<"Processing channel "<<chan+1<<" of "
189 0 : <<nchan<<endl;message.message(o);
190 0 : logSink().post(message);
191 : }
192 :
193 0 : const Matrix<Float>& lpsf=psfli.matrixCursor();
194 : // Make IPositions and find position of peak of PSF
195 0 : IPosition psfposmax(lpsf.ndim());
196 0 : IPosition psfposmin(lpsf.ndim());
197 : Float psfmax;
198 : Float psfmin;
199 0 : minMax(psfmin, psfmax, psfposmin, psfposmax, lpsf);
200 :
201 : // Loop through masks setting AMatrix as required
202 0 : if(psfmax==0.0) {
203 0 : ostringstream o; o<<"No data for this channel: skipping";
204 0 : message.message(o);
205 0 : logSink().post(message);
206 : }
207 : else {
208 0 : int xvi=0;
209 0 : int bvi=0;
210 0 : for(fai.origin();!fai.pastEnd();fai.next()) {
211 0 : if(lfluxmask(fai.pos())>0.000001) {
212 0 : bvi=0;
213 0 : for(dai.origin();!dai.pastEnd();dai.next()) {
214 0 : if(lmask(dai.pos())>0.0) {
215 0 : AMatrix(bvi,xvi)=lpsf(psfposmax+fai.pos()-dai.pos());
216 0 : bvi++;
217 : }
218 : }
219 0 : xvi++;
220 : }
221 : }
222 :
223 : // Construct the X vector.
224 0 : xvi=0;
225 0 : for(fai.origin();!fai.pastEnd();fai.next()) {
226 0 : if(lfluxmask(fai.pos())>0.000001) {
227 0 : XVector(xvi)=limage(fai.pos());
228 0 : xvi++;
229 : }
230 : }
231 :
232 : // Construct the B vector
233 0 : bvi=0;
234 0 : for(dai.origin();!dai.pastEnd();dai.next()) {
235 0 : if(lmask(dai.pos())>0.0000001) {
236 0 : BVector(bvi)=limageStep(dai.pos());
237 0 : bvi++;
238 : }
239 : }
240 :
241 : // Set A matrix and B vector in the matrixSolver
242 0 : matrixSolver.setAB(AMatrix, BVector);
243 :
244 : // Call MatrixSolver
245 : {
246 0 : ostringstream o; o<<"Performing solution";message.message(o);
247 0 : logSink().post(message);
248 : }
249 0 : if(!matrixSolver.solve()) {
250 : {
251 0 : ostringstream o; o<<"matrixSolver nominally failed";
252 0 : message.message(o);
253 0 : logSink().post(message);
254 : }
255 : }
256 0 : XVector=matrixSolver.getSolution();
257 :
258 : // Fill in final image and residual image
259 0 : xvi=0;
260 0 : for(fai.origin();!fai.pastEnd();fai.next()) {
261 0 : if(lfluxmask(fai.pos())>0.0000001) {
262 0 : limage(fai.pos())=XVector(xvi);
263 0 : xvi++;
264 : }
265 : }
266 0 : bvi=0;
267 0 : for(dai.origin();!dai.pastEnd();dai.next()) {
268 0 : if(lmask(dai.pos())>0.0000001) {
269 0 : limageStep(dai.pos())=BVector(bvi);
270 0 : bvi++;
271 : }
272 : }
273 :
274 0 : if(npol>1) {
275 0 : imageli.rwCubeCursor().xyPlane(0)=limage;
276 : }
277 : else {
278 0 : imageli.woMatrixCursor()=limage;
279 : }
280 0 : if(npol>1) {
281 0 : imageStepli.rwCubeCursor().xyPlane(0)=limageStep;
282 : }
283 : else {
284 0 : imageStepli.woMatrixCursor()=limageStep;
285 : }
286 : }
287 : }
288 : }
289 0 : return(true);
290 : };
291 :
292 : // Zero Stokes I for masked region
293 0 : Bool NNLSImageSkyModel::maskedZeroI() {
294 :
295 0 : LogMessage message(LogOrigin("NNLSImageSkyModel","maskedZeroI"));
296 :
297 0 : if(isSolveable(0)) {
298 :
299 0 : Int nx=image(0).shape()(0);
300 0 : Int ny=image(0).shape()(1);
301 0 : Int npol=image(0).shape()(2);
302 :
303 0 : LatticeStepper ls(image(0).shape(), IPosition(4, nx, ny, npol, 1),
304 0 : IPosition(4, 0, 1, 2, 3));
305 0 : LatticeIterator<Float> imageli(image(0), ls);
306 :
307 : // Get the flux mask?
308 0 : ArrayPositionIterator fai(IPosition(4, nx, ny, npol, 1), 0);
309 0 : Matrix<Float> lfluxmask;
310 0 : if(hasFluxMask(0)) {
311 0 : LatticeStepper mls(fluxMask(0).shape(), IPosition(4, nx, ny, npol, 1),
312 0 : IPosition(4, 0, 1, 2, 3));
313 0 : RO_LatticeIterator<Float> fluxmaskli(fluxMask(0), mls);
314 0 : fluxmaskli.reset();
315 0 : if(npol>1) {
316 0 : lfluxmask=fluxmaskli.cubeCursor().xyPlane(0);
317 : }
318 : else {
319 0 : lfluxmask=fluxmaskli.matrixCursor();
320 : }
321 : }
322 : else {
323 0 : lfluxmask.resize(nx,ny);
324 0 : lfluxmask=1.0;
325 : }
326 :
327 0 : Int chan=0;
328 0 : for (imageli.reset();!imageli.atEnd();imageli++,chan++) {
329 :
330 0 : for(fai.origin();!fai.pastEnd();fai.next()) {
331 0 : if(lfluxmask(fai.pos())>0.000001) {
332 0 : if(npol>1) {
333 0 : imageli.rwCubeCursor().xyPlane(0)(fai.pos())=0.0;
334 : }
335 : else {
336 0 : imageli.rwMatrixCursor()(fai.pos())=0.0;
337 : }
338 : }
339 : }
340 : }
341 : }
342 :
343 0 : return(true);
344 : };
345 :
346 :
347 : } //# NAMESPACE CASA - END
348 :
|
