Line data Source code
1 : //# PClarkCleanImageSkyModel.cc: Implementation of PClarkCleanImageSkyModel.h
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: 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 <synthesis/MeasurementComponents/PClarkCleanImageSkyModel.h>
30 : #include <casacore/images/Images/PagedImage.h>
31 : #include <casacore/casa/OS/File.h>
32 : #include <casacore/lattices/Lattices/LatticeStepper.h>
33 : #include <casacore/lattices/Lattices/LatticeIterator.h>
34 : #include <casacore/lattices/LEL/LatticeExpr.h>
35 : #include <casacore/lattices/LEL/LatticeExprNode.h>
36 : #include <casacore/lattices/LRegions/LCBox.h>
37 : #include <casacore/lattices/Lattices/SubLattice.h>
38 : #include <synthesis/MeasurementEquations/SkyEquation.h>
39 : #include <casacore/casa/Exceptions/Error.h>
40 : #include <casacore/casa/BasicSL/String.h>
41 : #include <casacore/casa/Utilities/Assert.h>
42 :
43 : #include <sstream>
44 :
45 : #include <casacore/casa/Logging/LogMessage.h>
46 : #include <casacore/casa/Logging/LogIO.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 :
49 : #include <synthesis/MeasurementEquations/LatConvEquation.h>
50 : #include <synthesis/MeasurementEquations/ClarkCleanLatModel.h>
51 : #include <synthesis/MeasurementEquations/ClarkCleanProgress.h>
52 :
53 : #include <synthesis/MeasurementComponents/ClarkCleanAlgorithm.h>
54 : #include <synthesis/Parallel/Applicator.h>
55 :
56 : using namespace casacore;
57 : namespace casa { //# NAMESPACE CASA - BEGIN
58 :
59 : extern Applicator applicator;
60 :
61 : // Clean solver
62 0 : Bool PClarkCleanImageSkyModel::solve(SkyEquation& se) {
63 :
64 0 : LogIO os(LogOrigin("PClarkCleanImageSkyModel","solve",WHERE));
65 :
66 :
67 0 : if(numberOfModels()>1) {
68 0 : os << "Cannot process more than one field" << LogIO::EXCEPTION;
69 : }
70 :
71 : // Make the residual image
72 0 : makeNewtonRaphsonStep(se);
73 :
74 : //Make the PSF
75 0 : makeApproxPSFs(se);
76 :
77 0 : Int nx=image(0).shape()(0);
78 0 : Int ny=image(0).shape()(1);
79 0 : Int npol=image(0).shape()(2);
80 0 : Int nchan=image(0).shape()(3);
81 :
82 0 : AlwaysAssert((npol==1)||(npol==2)||(npol==4), AipsError);
83 :
84 : // Now read the mask and determine the bounding box
85 0 : Int xbeg=nx/4;
86 0 : Int ybeg=ny/4;
87 :
88 0 : Int xend=xbeg+nx/2-1;
89 0 : Int yend=ybeg+ny/2-1;
90 :
91 0 : if(hasMask(0)) {
92 0 : AlwaysAssert(mask(0).shape()(0)==nx, AipsError);
93 0 : AlwaysAssert(mask(0).shape()(1)==ny, AipsError);
94 :
95 0 : LatticeStepper mls(mask(0).shape(),
96 0 : IPosition(4, nx, ny, 1, 1),
97 0 : IPosition(4, 0, 1, 3, 2));
98 0 : RO_LatticeIterator<Float> maskli(mask(0), mls);
99 0 : maskli.reset();
100 0 : xbeg=nx-1;
101 0 : ybeg=ny-1;
102 0 : xend=0;
103 0 : yend=0;
104 0 : for (Int iy=0;iy<ny;iy++) {
105 0 : for (Int ix=0;ix<nx;ix++) {
106 0 : if(maskli.matrixCursor()(ix,iy)>0.000001) {
107 0 : xbeg=min(xbeg,ix);
108 0 : ybeg=min(ybeg,iy);
109 0 : xend=max(xend,ix);
110 0 : yend=max(yend,iy);
111 : }
112 : }
113 : }
114 : // Now have possible BLC. Make sure that we don't go over the
115 : // edge later
116 0 : if(xbeg>nx/2) xbeg=nx/2;
117 0 : if(ybeg>ny/2) ybeg=ny/2;
118 0 : }
119 0 : xend=min(xend,xbeg+nx/2-1);
120 0 : yend=min(yend,ybeg+ny/2-1);
121 :
122 : // Mask logic used here: one mask for all Stokes, for all Channels
123 0 : SubLattice<Float>* mask_sl_p = 0;
124 0 : if (hasMask(0) & (xend > xbeg) && (yend > ybeg) ) {
125 0 : LCBox maskbox (IPosition(4, xbeg, ybeg, 0, 0),
126 0 : IPosition(4, xend, yend, 0, 0), mask(0).shape());
127 0 : mask_sl_p = new SubLattice<Float> (mask(0), maskbox, false);
128 0 : }
129 :
130 : // Start of the parallelization (over channel no.)
131 0 : Bool rStat=true;
132 0 : Array<Float> maskTmp;
133 0 : if (mask_sl_p) {
134 0 : rStat = mask_sl_p->get(maskTmp);
135 : }
136 :
137 0 : os << "Begin parallel clean" << LogIO::NORMAL << LogIO::POST;
138 0 : ClarkCleanAlgorithm clarkClean;
139 :
140 : // Track the assignment of channel to process rank
141 0 : std::map<Int, Int> chanNo;
142 : Bool allDone, assigned;
143 : Int rank;
144 :
145 0 : for (Int ichan=0; ichan < nchan; ichan++) {
146 0 : LCBox imagebox(IPosition(4, xbeg, ybeg, 0, ichan),
147 0 : IPosition(4, xend, yend, npol-1, ichan),
148 0 : image(0).shape());
149 0 : LCBox psfbox(IPosition(4, 0, 0, 0, ichan),
150 0 : IPosition(4, nx-1, ny-1, 0, ichan),
151 0 : PSF(0).shape());
152 :
153 0 : SubLattice<Float> psf_sl (PSF(0), psfbox, false);
154 0 : SubLattice<Float> residual_sl (residual(0), imagebox, true);
155 0 : SubLattice<Float> model_sl (image(0), imagebox, true);
156 :
157 : // Assign the next available process to this algorithm
158 0 : assigned = applicator.nextAvailProcess(clarkClean, rank);
159 0 : while (!assigned) {
160 : // No free processes; wait for any Clark Clean worker
161 : // process to return.
162 0 : rank = applicator.nextProcessDone(clarkClean, allDone);
163 : // Get the resulting plane and insert in the correct place
164 0 : Array<Float> af;
165 0 : applicator.get(af);
166 0 : image(0).putSlice(af, IPosition(4, xbeg, ybeg, 0, chanNo.at(rank)));
167 : // Assign the next available process
168 0 : assigned = applicator.nextAvailProcess(clarkClean, rank);
169 0 : };
170 :
171 : // Send the current channel to the assigned worker process
172 0 : Array<Float> imageTmp;
173 0 : rStat = ((Lattice<Float> &)residual_sl).get(imageTmp);
174 0 : applicator.put(imageTmp);
175 :
176 0 : Array<Float> psfTmp;
177 0 : rStat = ((Lattice<Float> &)psf_sl).get(psfTmp);
178 0 : applicator.put(psfTmp);
179 0 : applicator.put(maskTmp);
180 :
181 0 : applicator.put(gain());
182 0 : applicator.put(threshold());
183 0 : applicator.put(numberIterations());
184 0 : applicator.put(ichan);
185 0 : applicator.put(nchan);
186 :
187 : // Record the assignment of channel to process rank
188 0 : chanNo.insert( std::pair<Int, Int>(rank, ichan) );
189 :
190 : // Execute the algorithm
191 0 : applicator.apply(clarkClean);
192 0 : };
193 :
194 : // Wait for all outstanding processes to return
195 0 : rank = applicator.nextProcessDone(clarkClean, allDone);
196 0 : while (!allDone) {
197 : // Get the resulting plane and insert in the correct place
198 0 : Array<Float> af;
199 0 : applicator.get(af);
200 0 : image(0).putSlice(af, IPosition(4, xbeg, ybeg, 0, chanNo.at(rank)));
201 : // Wait for the next process to complete
202 0 : rank = applicator.nextProcessDone(clarkClean, allDone);
203 0 : };
204 :
205 : (void)rStat; // avoid compiler warning
206 0 : if (mask_sl_p != 0 ) delete mask_sl_p;
207 0 : os << "End parallel clean" << LogIO::NORMAL << LogIO::POST;
208 0 : return(rStat);
209 0 : };
210 :
211 :
212 :
213 :
214 : } //# NAMESPACE CASA - END
215 :
|
