Line data Source code
1 : //# CubeMinorCycleAlgorithm.cc: implementation of class to deconvolve cube in parallel/serial
2 : //# Copyright (C) 2019
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 General Public License as published by
7 : //# the Free Software Foundation; either version 3 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 General Public
13 : //# License for more details.
14 : //#
15 : //# https://www.gnu.org/licenses/
16 : //#
17 : //# Queries concerning CASA should be submitted at
18 : //# https://help.nrao.edu
19 : //#
20 : //# Postal address: CASA Project Manager
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 : //#
26 : //# $Id$
27 : #include <casacore/lattices/Lattices/LatticeLocker.h>
28 : #include <synthesis/ImagerObjects/CubeMinorCycleAlgorithm.h>
29 : #include <synthesis/ImagerObjects/SynthesisDeconvolver.h>
30 : #include <casacore/casa/Containers/Record.h>
31 : #include <synthesis/ImagerObjects/SimpleSIImageStore.h>
32 : #include <imageanalysis/Utilities/SpectralImageUtil.h>
33 :
34 : using namespace casacore;
35 : namespace casa { //# NAMESPACE CASA - BEGIN
36 : extern Applicator applicator;
37 :
38 265 : CubeMinorCycleAlgorithm::CubeMinorCycleAlgorithm() : myName_p("CubeMinorCycleAlgorithm"), autoMaskOn_p(False),chanFlag_p(0), status_p(False){
39 :
40 265 : }
41 266 : CubeMinorCycleAlgorithm::~CubeMinorCycleAlgorithm() {
42 :
43 266 : }
44 :
45 264 : void CubeMinorCycleAlgorithm::get() {
46 : ///New instructions reset previous state
47 264 : reset();
48 : //cerr << "in get for child process " << applicator.isWorker() << endl;
49 528 : Record decParsRec;
50 :
51 : // get deconv params record #1
52 264 : applicator.get(decParsRec);
53 : // get iter control rec #2
54 264 : applicator.get(iterBotRec_p);
55 : // channel range to deconvolve #3
56 264 : applicator.get(chanRange_p);
57 : //get psf image name #4
58 264 : applicator.get(psfName_p);
59 : //get residual name #5
60 264 : applicator.get(residualName_p);
61 : //get model name #6
62 264 : applicator.get(modelName_p);
63 : //get mask name #7
64 264 : applicator.get(maskName_p);
65 : // get pb name #8
66 264 : applicator.get(pbName_p);
67 : //get beamsetrec #9
68 : //applicator.get(beamsetRec_p);
69 : //get psfsidelobelev #9
70 264 : applicator.get(psfSidelobeLevel_p);
71 : //get chanflag #10
72 528 : Record chanflagRec;
73 264 : applicator.get(chanflagRec);
74 264 : chanFlag_p.resize();
75 264 : chanflagRec.get("chanflag", chanFlag_p);
76 264 : statsRec_p=chanflagRec.asRecord("statsrec");
77 : //cerr <<"GET statsRec " << statsRec_p << endl;
78 264 : decPars_p.fromRecord(decParsRec);
79 :
80 :
81 :
82 264 : }
83 264 : void CubeMinorCycleAlgorithm::put() {
84 :
85 : ///# 1 chanrange processed
86 264 : applicator.put(chanRange_p);
87 : //cerr << "in put " << status_p << endl;
88 : //#2 chanflag
89 264 : chanFlagRec_p.define("chanflag", chanFlag_p);
90 : //cerr << "PUT statsRec "<< statsRec_p << endl;
91 264 : chanFlagRec_p.defineRecord("statsrec", statsRec_p);
92 264 : applicator.put(chanFlagRec_p);
93 : ///#3 return record of deconvolver
94 : // cerr << "nfield " << returnRec_p.nfields() << endl;
95 264 : applicator.put(returnRec_p);
96 :
97 264 : }
98 :
99 264 : void CubeMinorCycleAlgorithm::task(){
100 528 : LogIO logger(LogOrigin("CubeMinorCycleAlgorithm", "task", WHERE));
101 :
102 264 : status_p = False;
103 : try{
104 :
105 528 : SynthesisDeconvolver subDeconv;
106 264 : Bool writeBackAutomask=True;
107 264 : subDeconv.setupDeconvolution(decPars_p);
108 528 : std::shared_ptr<SIImageStore> subimstor=subImageStore();
109 : //ImageBeamSet bs=ImageBeamSet::fromRecord(beamsetRec_p);
110 528 : ImageBeamSet bs=(subimstor->psf()->imageInfo()).getBeamSet();
111 264 : subimstor->setBeamSet(bs);
112 264 : subimstor->setPSFSidelobeLevel(psfSidelobeLevel_p);
113 528 : LatticeLocker lock1 (*(subimstor->model()), FileLocker::Write, 30);
114 528 : Record prevresrec=subDeconv.initMinorCycle(subimstor);
115 :
116 : // unused
117 : // Float prevPeakRes=prevresrec.asFloat("peakresidual");
118 264 : Bool doDeconv=True;
119 264 : if(autoMaskOn_p){
120 15 : subDeconv.setChanFlag(chanFlag_p);
121 15 : subDeconv.setRobustStats(statsRec_p);
122 15 : Int automaskflag=iterBotRec_p.asInt("onlyautomask");
123 45 : LogIO os( LogOrigin("SynthesisDeconvolver",automaskflag==1 ? "excuteAutomask" : "executeMinorCycle",WHERE) );
124 15 : os << "Processing channels in range " << chanRange_p << LogIO::POST;
125 15 : if(automaskflag==1){
126 10 : doDeconv=False;
127 10 : if(iterBotRec_p.isDefined("cycleniter"))
128 5 : subDeconv.setMinorCycleControl(iterBotRec_p);
129 : }
130 : //cerr << "ITERDONE " << iterBotRec_p.asInt("iterdone")<< " itermask flag " << automaskflag << endl;
131 15 : subDeconv.setIterDone(iterBotRec_p.asInt("iterdone"));
132 15 : if(automaskflag !=0){
133 : //this is already sent in as part of subimstor
134 : //subDeconv.setPosMask(subimstor->tempworkimage());
135 :
136 10 : subDeconv.setAutoMask();
137 : /*
138 : Record resRec=subDeconv.initMinorCycle(subimstor);
139 : //cerr << "POST resrec " << resRec << endl;
140 : //cerr << "peakResidual " << resRec.asFloat("peakresidual") << " cyclethreshold " << iterBotRec_p.asFloat("cyclethreshold") << " as double " << iterBotRec_p.asDouble("cyclethreshold") << endl;
141 : if(resRec.isDefined("peakresidual") && iterBotRec_p.isDefined("cyclethreshold")){
142 :
143 : Float peakresidual=resRec.asFloat("peakresidual");
144 : Float peakresidualnomask=resRec.asFloat("peakresidualnomask");
145 : Float cyclethreshold=iterBotRec_p.asFloat("cyclethreshold");
146 : if(iterBotRec_p.isDefined("psffraction")){
147 : cyclethreshold=iterBotRec_p.asFloat("psffraction")*peakresidual;
148 : cyclethreshold= max(cyclethreshold, iterBotRec_p.asFloat("threshold"));
149 : if(automaskflag==-1){
150 :
151 : // cerr << "old cyclethreshold " <<iterBotRec_p.asFloat("cyclethreshold") << " new " << cyclethreshold << endl;
152 : iterBotRec_p.removeField("cyclethreshold");
153 : iterBotRec_p.define("cyclethreshold", cyclethreshold);
154 : }
155 : }
156 :
157 :
158 : //if(peakresidual < iterBotRec_p.asFloat("cyclethreshold"))
159 : if(peakresidual < cyclethreshold)
160 : writeBackAutomask=False;
161 :
162 : writeBackAutomask=False;
163 : //Its better to always write the automask
164 : //cerr << "chanRange " << chanRange_p << endl;
165 : //cerr << "WRITEBACK " << writeBackAutomask<< " peakres " << peakresidual << " prev " << prevPeakRes << " nomask " << peakresidualnomask << endl;
166 :
167 : }
168 : */
169 :
170 10 : writeBackAutomask=True;
171 : //Its better to always write the automask
172 :
173 : }
174 : else{
175 5 : writeBackAutomask=False;
176 : }
177 :
178 : } else {
179 747 : LogIO os( LogOrigin("CubeMinorCycleAlgorithm","task",WHERE) );
180 249 : os << "Processing channels in range " << chanRange_p << LogIO::POST;
181 : } // end if(autoMaskOn_p)
182 : //subDeconv.setupMask();
183 264 : if(doDeconv)
184 254 : returnRec_p=subDeconv.executeCoreMinorCycle(iterBotRec_p);
185 : else
186 10 : returnRec_p.define("doneautomask", True);
187 264 : chanFlag_p.resize();
188 264 : chanFlag_p=subDeconv.getChanFlag();
189 264 : statsRec_p=Record();
190 264 : statsRec_p=subDeconv.getRobustStats();
191 264 : if(doDeconv){
192 254 : writeBackToFullImage(modelName_p, chanRange_p[0], chanRange_p[1], (subimstor->model()));
193 254 : writeBackToFullImage(residualName_p, chanRange_p[0], chanRange_p[1], (subimstor->residual()));
194 :
195 : }
196 264 : if(autoMaskOn_p && writeBackAutomask){
197 10 : writeBackToFullImage(posMaskName_p, chanRange_p[0], chanRange_p[1], (subimstor->tempworkimage()));
198 10 : writeBackToFullImage(maskName_p, chanRange_p[0], chanRange_p[1], (subimstor->mask()));
199 : }
200 : }
201 0 : catch (AipsError x) {
202 0 : logger << LogIO::SEVERE << "Exception: " << x.getMesg() << LogIO::POST;
203 0 : returnRec_p=Record();
204 : }
205 0 : catch(std::exception& exc) {
206 0 : logger << LogIO::SEVERE << "Exception (std): " << exc.what() << LogIO::POST;
207 0 : returnRec_p=Record();
208 : }
209 0 : catch(...){
210 0 : logger << LogIO::SEVERE << "Unknown exception" << LogIO::POST;
211 0 : returnRec_p=Record();
212 : }
213 :
214 264 : status_p = True;
215 264 : }
216 1057 : String& CubeMinorCycleAlgorithm::name(){
217 1057 : return myName_p;
218 : }
219 :
220 264 : std::shared_ptr<SIImageStore> CubeMinorCycleAlgorithm::subImageStore(){
221 528 : std::shared_ptr<ImageInterface<Float> >subpsf=nullptr;
222 528 : std::shared_ptr<ImageInterface<Float> >subresid=nullptr;
223 528 : std::shared_ptr<ImageInterface<Float> >submodel=nullptr;
224 528 : std::shared_ptr<ImageInterface<Float> > submask=nullptr;
225 528 : std::shared_ptr<ImageInterface<Float> > subpb=nullptr;
226 528 : std::shared_ptr<ImageInterface<Float> > subposmask=nullptr;
227 264 : Int chanBeg=0;
228 264 : Int chanEnd=0;
229 264 : chanBeg=chanRange_p[0];
230 264 : chanEnd=chanRange_p[1];
231 528 : casacore::String imageName = decPars_p.imageName;
232 : //cerr << "chanBeg " << chanBeg << " chanEnd " << chanEnd << " imId " << imId << endl;
233 :
234 :
235 264 : makeTempImage(subpsf, psfName_p, chanBeg, chanEnd);
236 : ///Have to lock residual now that deconvolver task now will want to writeback residual
237 : ///tclean does not need writeback for residual so it is unnecessary step for now.
238 264 : makeTempImage(subresid, residualName_p, chanBeg, chanEnd, True);
239 264 : makeTempImage(submodel, modelName_p, chanBeg, chanEnd, True);
240 264 : makeTempImage(submask, maskName_p, chanBeg, chanEnd, True);
241 : // PagedImage<Float> model(modelName_p, TableLock::UserLocking);
242 : // submodel.reset(SpectralImageUtil::getChannel(model, chanBeg, chanEnd, true));
243 264 : if(!pbName_p.empty()){
244 264 : makeTempImage(subpb, pbName_p, chanBeg, chanEnd);
245 : }
246 264 : if(iterBotRec_p.isDefined("posmaskname") ){
247 15 : iterBotRec_p.get("posmaskname", posMaskName_p);
248 15 : if(Table::isReadable(posMaskName_p)){
249 15 : makeTempImage(subposmask, posMaskName_p, chanBeg, chanEnd, True);
250 15 : if(subposmask)
251 15 : autoMaskOn_p=True;
252 : }
253 : }
254 :
255 528 : std::shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imageName, submodel, subresid, subpsf, nullptr, nullptr, submask, nullptr, nullptr, subpb, nullptr, subposmask));
256 :
257 : //cerr << "subimagestor TYPE" << subimstor->getType() << endl;
258 528 : return subimstor;
259 : }
260 1335 : void CubeMinorCycleAlgorithm::makeTempImage(std::shared_ptr<ImageInterface<Float> >& outptr, const String& imagename, const Int chanBeg, const Int chanEnd, const Bool writelock){
261 : //For testing that locks are placed in the right places use userlocking
262 2670 : PagedImage<Float> im(imagename, writelock ? TableLock::UserLocking : TableLock::UserNoReadLocking);
263 : //PagedImage<Float> im(imagename, writelock ? TableLock::AutoLocking : TableLock::AutoNoReadLocking);
264 :
265 1335 : SubImage<Float> *tmpptr=nullptr;
266 : //LatticeLocker lockread (im, FileLocker::Read);
267 : //if(writelock)
268 : // im.lock(FileLocker::Write, 1000);
269 : ////TESTOO
270 : //outptr.reset(SpectralImageUtil::getChannel(im, chanBeg, chanEnd, writelock));
271 :
272 :
273 : ///END of TESTOO
274 1335 : if(writelock){
275 807 : im.lock(FileLocker::Write, 1000);
276 807 : tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, false);
277 : }
278 : else
279 528 : tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, false);
280 1335 : if(tmpptr){
281 2670 : IPosition tileshape=tmpptr->shape();
282 1335 : tileshape[2]=1; tileshape[3]=1;
283 2670 : TiledShape tshape(tmpptr->shape(),tileshape);
284 : ///TESTOO
285 : /*if(imagename.contains(".residual")){
286 : String newresid=File::newUniqueName("./", "newResid").absoluteName();
287 : outptr.reset(new PagedImage<Float>(tshape, tmpptr->coordinates(), newresid));
288 : }
289 : ////TESTOO
290 : else
291 : */
292 1335 : outptr.reset(new TempImage<Float>(tshape, tmpptr->coordinates()));
293 :
294 :
295 1335 : if(writelock){
296 1614 : LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
297 807 : outptr->copyData(*tmpptr);
298 : //cerr << "IMAGENAME " << imagename << " masked " << im.isMasked() << " tmptr " << tmpptr->isMasked() << endl;
299 807 : if(tmpptr->isMasked()){
300 246 : outptr->makeMask ("mask0", true, true, false, true);
301 246 : outptr->pixelMask().put(tmpptr->getMask());
302 : // cerr << "tempimage SUM of bit mask" << ntrue(tmpptr->pixelMask().get()) << " out " << ntrue(outptr->pixelMask().get()) << endl;
303 : }
304 : }
305 : else{
306 1056 : LatticeLocker lock1 (*(tmpptr), FileLocker::Read);
307 528 : outptr->copyData(*tmpptr);
308 : //cerr << "false IMAGENAME " << imagename << " masked " << im.isMasked() << " tmptr " << tmpptr->isMasked() << endl;
309 :
310 528 : if(tmpptr->isMasked()){
311 265 : outptr->makeMask ("mask0", true, true, false, true);
312 265 : outptr->pixelMask().put(tmpptr->getMask());
313 : //cerr << "tempimage SUM of bit mask" << ntrue(tmpptr->pixelMask().get()) << " out " << ntrue(outptr->pixelMask().get()) << endl;
314 : }
315 : }
316 2670 : ImageInfo iinfo=tmpptr->imageInfo();
317 1335 : outptr->setImageInfo(iinfo);
318 1335 : delete tmpptr;
319 : }
320 :
321 1335 : im.unlock();
322 1335 : }
323 528 : void CubeMinorCycleAlgorithm::writeBackToFullImage(const String imagename, const Int chanBeg, const Int chanEnd, std::shared_ptr<ImageInterface<Float> > subimptr){
324 1056 : PagedImage<Float> im(imagename, TableLock::UserLocking);
325 : //PagedImage<Float> im(imagename, TableLock::AutoLocking);
326 528 : SubImage<Float> *tmpptr=nullptr;
327 528 : tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, true);
328 : {
329 :
330 1056 : LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
331 528 : tmpptr->copyData(*(subimptr));
332 : }
333 528 : im.unlock();
334 528 : delete tmpptr;
335 :
336 528 : }
337 264 : void CubeMinorCycleAlgorithm::reset(){
338 :
339 264 : iterBotRec_p=Record();
340 264 : modelName_p="";
341 264 : residualName_p="";
342 264 : psfName_p="";
343 264 : maskName_p="";
344 264 : pbName_p="";
345 264 : posMaskName_p="";
346 264 : chanRange_p.resize();
347 264 : returnRec_p=Record();
348 264 : beamsetRec_p=Record();
349 : //psfSidelobeLevel_p;
350 264 : autoMaskOn_p=False;
351 264 : chanFlag_p.resize();
352 264 : statsRec_p=Record();
353 264 : chanFlagRec_p=Record();
354 264 : status_p=False;
355 :
356 :
357 :
358 264 : }
359 :
360 :
361 : } //# NAMESPACE CASA - END
|
