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 0 : CubeMinorCycleAlgorithm::CubeMinorCycleAlgorithm() : myName_p("CubeMinorCycleAlgorithm"), autoMaskOn_p(False),chanFlag_p(0), status_p(False){
39 :
40 0 : }
41 0 : CubeMinorCycleAlgorithm::~CubeMinorCycleAlgorithm() {
42 :
43 0 : }
44 :
45 0 : void CubeMinorCycleAlgorithm::get() {
46 : ///New instructions reset previous state
47 0 : reset();
48 : //cerr << "in get for child process " << applicator.isWorker() << endl;
49 0 : Record decParsRec;
50 :
51 : // get deconv params record #1
52 0 : applicator.get(decParsRec);
53 : // get iter control rec #2
54 0 : applicator.get(iterBotRec_p);
55 : // channel range to deconvolve #3
56 0 : applicator.get(chanRange_p);
57 : //get psf image name #4
58 0 : applicator.get(psfName_p);
59 : //get residual name #5
60 0 : applicator.get(residualName_p);
61 : //get model name #6
62 0 : applicator.get(modelName_p);
63 : //get mask name #7
64 0 : applicator.get(maskName_p);
65 : // get pb name #8
66 0 : applicator.get(pbName_p);
67 : //get beamsetrec #9
68 : //applicator.get(beamsetRec_p);
69 : //get psfsidelobelev #9
70 0 : applicator.get(psfSidelobeLevel_p);
71 : //get chanflag #10
72 0 : Record chanflagRec;
73 0 : applicator.get(chanflagRec);
74 0 : chanFlag_p.resize();
75 0 : chanflagRec.get("chanflag", chanFlag_p);
76 0 : statsRec_p=chanflagRec.asRecord("statsrec");
77 : //cerr <<"GET statsRec " << statsRec_p << endl;
78 0 : decPars_p.fromRecord(decParsRec);
79 :
80 :
81 :
82 0 : }
83 0 : void CubeMinorCycleAlgorithm::put() {
84 :
85 : ///# 1 chanrange processed
86 0 : applicator.put(chanRange_p);
87 : //cerr << "in put " << status_p << endl;
88 : //#2 chanflag
89 0 : chanFlagRec_p.define("chanflag", chanFlag_p);
90 : //cerr << "PUT statsRec "<< statsRec_p << endl;
91 0 : chanFlagRec_p.defineRecord("statsrec", statsRec_p);
92 0 : applicator.put(chanFlagRec_p);
93 : ///#3 return record of deconvolver
94 : // cerr << "nfield " << returnRec_p.nfields() << endl;
95 0 : applicator.put(returnRec_p);
96 :
97 0 : }
98 :
99 0 : void CubeMinorCycleAlgorithm::task(){
100 0 : LogIO logger(LogOrigin("CubeMinorCycleAlgorithm", "task", WHERE));
101 :
102 0 : status_p = False;
103 : try{
104 :
105 0 : SynthesisDeconvolver subDeconv;
106 0 : Bool writeBackAutomask=True;
107 0 : subDeconv.setupDeconvolution(decPars_p);
108 0 : std::shared_ptr<SIImageStore> subimstor=subImageStore();
109 : //ImageBeamSet bs=ImageBeamSet::fromRecord(beamsetRec_p);
110 0 : ImageBeamSet bs=(subimstor->psf()->imageInfo()).getBeamSet();
111 0 : subimstor->setBeamSet(bs);
112 0 : subimstor->setPSFSidelobeLevel(psfSidelobeLevel_p);
113 0 : LatticeLocker lock1 (*(subimstor->model()), FileLocker::Write, 30);
114 0 : Record prevresrec=subDeconv.initMinorCycle(subimstor);
115 :
116 : // unused
117 : // Float prevPeakRes=prevresrec.asFloat("peakresidual");
118 0 : Bool doDeconv=True;
119 0 : if(autoMaskOn_p){
120 0 : subDeconv.setChanFlag(chanFlag_p);
121 0 : subDeconv.setRobustStats(statsRec_p);
122 0 : Int automaskflag=iterBotRec_p.asInt("onlyautomask");
123 0 : LogIO os( LogOrigin("SynthesisDeconvolver",automaskflag==1 ? "excuteAutomask" : "executeMinorCycle",WHERE) );
124 0 : os << "Processing channels in range " << chanRange_p << LogIO::POST;
125 0 : if(automaskflag==1){
126 0 : doDeconv=False;
127 0 : if(iterBotRec_p.isDefined("cycleniter"))
128 0 : subDeconv.setMinorCycleControl(iterBotRec_p);
129 : }
130 : //cerr << "ITERDONE " << iterBotRec_p.asInt("iterdone")<< " itermask flag " << automaskflag << endl;
131 0 : subDeconv.setIterDone(iterBotRec_p.asInt("iterdone"));
132 0 : if(automaskflag !=0){
133 : //this is already sent in as part of subimstor
134 : //subDeconv.setPosMask(subimstor->tempworkimage());
135 :
136 0 : 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 0 : writeBackAutomask=True;
171 : //Its better to always write the automask
172 :
173 : }
174 : else{
175 0 : writeBackAutomask=False;
176 : }
177 :
178 : } else {
179 0 : LogIO os( LogOrigin("CubeMinorCycleAlgorithm","task",WHERE) );
180 0 : os << "Processing channels in range " << chanRange_p << LogIO::POST;
181 : } // end if(autoMaskOn_p)
182 : //subDeconv.setupMask();
183 0 : if(doDeconv)
184 0 : returnRec_p=subDeconv.executeCoreMinorCycle(iterBotRec_p);
185 : else
186 0 : returnRec_p.define("doneautomask", True);
187 0 : chanFlag_p.resize();
188 0 : chanFlag_p=subDeconv.getChanFlag();
189 0 : statsRec_p=Record();
190 0 : statsRec_p=subDeconv.getRobustStats();
191 0 : if(doDeconv){
192 0 : writeBackToFullImage(modelName_p, chanRange_p[0], chanRange_p[1], (subimstor->model()));
193 0 : writeBackToFullImage(residualName_p, chanRange_p[0], chanRange_p[1], (subimstor->residual()));
194 :
195 : }
196 0 : if(autoMaskOn_p && writeBackAutomask){
197 0 : writeBackToFullImage(posMaskName_p, chanRange_p[0], chanRange_p[1], (subimstor->tempworkimage()));
198 0 : 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 0 : status_p = True;
215 0 : }
216 0 : String& CubeMinorCycleAlgorithm::name(){
217 0 : return myName_p;
218 : }
219 :
220 0 : std::shared_ptr<SIImageStore> CubeMinorCycleAlgorithm::subImageStore(){
221 0 : std::shared_ptr<ImageInterface<Float> >subpsf=nullptr;
222 0 : std::shared_ptr<ImageInterface<Float> >subresid=nullptr;
223 0 : std::shared_ptr<ImageInterface<Float> >submodel=nullptr;
224 0 : std::shared_ptr<ImageInterface<Float> > submask=nullptr;
225 0 : std::shared_ptr<ImageInterface<Float> > subpb=nullptr;
226 0 : std::shared_ptr<ImageInterface<Float> > subposmask=nullptr;
227 0 : Int chanBeg=0;
228 0 : Int chanEnd=0;
229 0 : chanBeg=chanRange_p[0];
230 0 : chanEnd=chanRange_p[1];
231 0 : casacore::String imageName = decPars_p.imageName;
232 : //cerr << "chanBeg " << chanBeg << " chanEnd " << chanEnd << " imId " << imId << endl;
233 :
234 :
235 0 : 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 0 : makeTempImage(subresid, residualName_p, chanBeg, chanEnd, True);
239 0 : makeTempImage(submodel, modelName_p, chanBeg, chanEnd, True);
240 0 : 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 0 : if(!pbName_p.empty()){
244 0 : makeTempImage(subpb, pbName_p, chanBeg, chanEnd);
245 : }
246 0 : if(iterBotRec_p.isDefined("posmaskname") ){
247 0 : iterBotRec_p.get("posmaskname", posMaskName_p);
248 0 : if(Table::isReadable(posMaskName_p)){
249 0 : makeTempImage(subposmask, posMaskName_p, chanBeg, chanEnd, True);
250 0 : if(subposmask)
251 0 : autoMaskOn_p=True;
252 : }
253 : }
254 :
255 0 : 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 0 : return subimstor;
259 : }
260 0 : 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 0 : PagedImage<Float> im(imagename, writelock ? TableLock::UserLocking : TableLock::UserNoReadLocking);
263 : //PagedImage<Float> im(imagename, writelock ? TableLock::AutoLocking : TableLock::AutoNoReadLocking);
264 :
265 0 : 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 0 : if(writelock){
275 0 : im.lock(FileLocker::Write, 1000);
276 0 : tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, false);
277 : }
278 : else
279 0 : tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, false);
280 0 : if(tmpptr){
281 0 : IPosition tileshape=tmpptr->shape();
282 0 : tileshape[2]=1; tileshape[3]=1;
283 0 : 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 0 : outptr.reset(new TempImage<Float>(tshape, tmpptr->coordinates()));
293 :
294 :
295 0 : if(writelock){
296 0 : LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
297 0 : outptr->copyData(*tmpptr);
298 : //cerr << "IMAGENAME " << imagename << " masked " << im.isMasked() << " tmptr " << tmpptr->isMasked() << endl;
299 0 : if(tmpptr->isMasked()){
300 0 : outptr->makeMask ("mask0", true, true, false, true);
301 0 : 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 0 : LatticeLocker lock1 (*(tmpptr), FileLocker::Read);
307 0 : outptr->copyData(*tmpptr);
308 : //cerr << "false IMAGENAME " << imagename << " masked " << im.isMasked() << " tmptr " << tmpptr->isMasked() << endl;
309 :
310 0 : if(tmpptr->isMasked()){
311 0 : outptr->makeMask ("mask0", true, true, false, true);
312 0 : 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 0 : ImageInfo iinfo=tmpptr->imageInfo();
317 0 : outptr->setImageInfo(iinfo);
318 0 : delete tmpptr;
319 : }
320 :
321 0 : im.unlock();
322 0 : }
323 0 : void CubeMinorCycleAlgorithm::writeBackToFullImage(const String imagename, const Int chanBeg, const Int chanEnd, std::shared_ptr<ImageInterface<Float> > subimptr){
324 0 : PagedImage<Float> im(imagename, TableLock::UserLocking);
325 : //PagedImage<Float> im(imagename, TableLock::AutoLocking);
326 0 : SubImage<Float> *tmpptr=nullptr;
327 0 : tmpptr=SpectralImageUtil::getChannel(im, chanBeg, chanEnd, true);
328 : {
329 :
330 0 : LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
331 0 : tmpptr->copyData(*(subimptr));
332 : }
333 0 : im.unlock();
334 0 : delete tmpptr;
335 :
336 0 : }
337 0 : void CubeMinorCycleAlgorithm::reset(){
338 :
339 0 : iterBotRec_p=Record();
340 0 : modelName_p="";
341 0 : residualName_p="";
342 0 : psfName_p="";
343 0 : maskName_p="";
344 0 : pbName_p="";
345 0 : posMaskName_p="";
346 0 : chanRange_p.resize();
347 0 : returnRec_p=Record();
348 0 : beamsetRec_p=Record();
349 : //psfSidelobeLevel_p;
350 0 : autoMaskOn_p=False;
351 0 : chanFlag_p.resize();
352 0 : statsRec_p=Record();
353 0 : chanFlagRec_p=Record();
354 0 : status_p=False;
355 :
356 :
357 :
358 0 : }
359 :
360 :
361 : } //# NAMESPACE CASA - END
|
