Line data Source code
1 : //# CubeMajorCycleAlgorithm.cc: implementation of class to grid and degrid (and write model vis when necessary) 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/CubeMajorCycleAlgorithm.h>
29 : #include <synthesis/ImagerObjects/SynthesisImagerVi2.h>
30 : #include <synthesis/ImagerObjects/SynthesisNormalizer.h>
31 : #include <casacore/casa/Containers/Record.h>
32 : #include <synthesis/ImagerObjects/SimpleSIImageStore.h>
33 : #include <imageanalysis/Utilities/SpectralImageUtil.h>
34 : #include <casacore/casa/OS/Timer.h>
35 : #include <chrono>
36 : #include <thread>
37 : using namespace casacore;
38 : namespace casa { //# NAMESPACE CASA - BEGIN
39 : extern Applicator applicator;
40 :
41 887 : CubeMajorCycleAlgorithm::CubeMajorCycleAlgorithm() : myName_p("CubeMajorCycleAlgorithm"), ftmRec_p(0), iftmRec_p(0), polRep_p(0),startmodel_p(0), residualNames_p(0), psfNames_p(0), sumwtNames_p(0), weightNames_p(0), pbNames_p(0), movingSource_p(""),status_p(False), retuning_p(True), doPB_p(False), nterms_p(0){
42 :
43 887 : }
44 891 : CubeMajorCycleAlgorithm::~CubeMajorCycleAlgorithm() {
45 :
46 891 : }
47 :
48 883 : void CubeMajorCycleAlgorithm::get() {
49 883 : reset();
50 : //cerr << "in get for child process " << applicator.isWorker() << endl;
51 883 : Record vecImParsRec;
52 883 : Record vecSelParsRec;
53 883 : Record vecGridParsRec;
54 : // get data sel params #1
55 883 : applicator.get(vecSelParsRec);
56 : // get image sel params #2
57 883 : applicator.get(vecImParsRec);
58 : // get gridders params #3
59 883 : applicator.get(vecGridParsRec);
60 : // get which channel to process #4
61 883 : applicator.get(chanRange_p);
62 : //cerr <<"GET chanRange " << chanRange_p << endl;
63 : // psf or residual CubeMajorCycleAlgorithm #5
64 883 : applicator.get(dopsf_p);
65 : // store modelvis and other controls #6
66 883 : applicator.get(controlRecord_p);
67 : // weight params
68 883 : applicator.get(weightParams_p);
69 : //Somewhere before this we have to make sure that vecSelPars has more than 0 fields)
70 883 : dataSel_p.resize(vecSelParsRec.nfields());
71 : /// Fill the private variables
72 1795 : for (Int k=0; k < Int(dataSel_p.nelements()); ++k){
73 912 : (dataSel_p[k]).fromRecord(vecSelParsRec.asRecord(String::toString(k)));
74 : // cerr << k << "datasel " << vecSelParsRec.asRecord(String::toString(k)) << endl;
75 : }
76 : //imsel and gridsel should be the same numbers (number of image fields)
77 883 : Int nmajorcycles=0;
78 883 : if(controlRecord_p.isDefined("nmajorcycles"))
79 583 : controlRecord_p.get("nmajorcycles",nmajorcycles);
80 883 : imSel_p.resize(vecImParsRec.nfields());
81 883 : gridSel_p.resize(vecImParsRec.nfields());
82 883 : ftmRec_p.resize(vecImParsRec.nfields());
83 883 : iftmRec_p.resize(vecImParsRec.nfields());
84 883 : polRep_p.resize(vecImParsRec.nfields());
85 883 : polRep_p.set(-1);
86 883 : startmodel_p.resize(vecImParsRec.nfields());
87 883 : startmodel_p.set(Vector<String>(0));
88 : //TESTOO
89 : //Int CPUID;
90 : //MPI_Comm_rank(MPI_COMM_WORLD, &CPUID);
91 : //////////////////
92 1781 : for (uInt k=0; k < imSel_p.nelements(); ++k){
93 1796 : Record imSelRec=vecImParsRec.asRecord(String::toString(k));
94 : //cerr << k << " imsel " << imSelRec << endl;
95 898 : if(imSelRec.isDefined("polrep"))
96 898 : imSelRec.get("polrep", polRep_p[k]);
97 : ///Get moving source name
98 898 : if(imSelRec.isDefined("movingsource") && imSelRec.asString("movingsource") != ""){
99 8 : imSelRec.get("movingsource", movingSource_p);
100 : }
101 : //cerr << CPUID << " POLREP " << polRep_p << endl;
102 : //Only first major cycle we need to reset model
103 898 : if(nmajorcycles==1)
104 308 : imSelRec.get("startmodel", startmodel_p[k]);
105 : //have to do that as fromRecord check does not like to have both model and startmodel on disk !
106 898 : imSelRec.define("startmodel", Vector<String>(0));
107 898 : (imSel_p[k]).fromRecord(imSelRec);
108 1796 : Record recGridSel=vecGridParsRec.asRecord(String::toString(k));
109 898 : (gridSel_p[k]).fromRecord(recGridSel);
110 898 : if(!recGridSel.isDefined("ftmachine")){
111 898 : ftmRec_p.resize();
112 898 : iftmRec_p.resize();
113 : }
114 898 : if(ftmRec_p.nelements() >0){
115 0 : ftmRec_p[k]=recGridSel.asRecord("ftmachine");
116 0 : iftmRec_p[k]=recGridSel.asRecord("iftmachine");
117 : }
118 :
119 898 : }
120 :
121 :
122 883 : }
123 883 : void CubeMajorCycleAlgorithm::put() {
124 :
125 : //if(applicator.isSerial()){
126 : // serialBug_p=Applicator::DONE;
127 : //applicator.put(serialBug_p);
128 :
129 : // }
130 : //cerr << "in put " << status_p << endl;
131 883 : applicator.put(returnRec_p);
132 883 : applicator.put(status_p);
133 :
134 883 : }
135 :
136 883 : void CubeMajorCycleAlgorithm::task(){
137 1766 : LogIO logger(LogOrigin("CubeMajorCycleAlgorithm", "task", WHERE));
138 883 : status_p = True;
139 : try{
140 : //Timer tim;
141 : //tim.mark();
142 : //SynthesisImagerVi2 imgr;
143 : //imgr.selectData(dataSel_p);
144 : // We do not use chanchunking in this model
145 1781 : for (uInt k=0; k < gridSel_p.nelements(); ++k)
146 898 : gridSel_p[k].chanchunks = 1;
147 :
148 : //imgr.defineImage(imSel_p,gridSel_p);
149 : // need to find how many subfields/outliers have been set
150 : //CountedPtr<SIImageStore> imstor =imgr.imageStore(0);
151 : //CountedPtr<ImageInterface<Float> > resid=imstor->residual();
152 : //Int nchan = resid->shape()(3);
153 : //std::shared_ptr<SIImageStore> subImStor=imstor->getSubImageStore(0, 1, chanId_p, nchan, 0,1);
154 :
155 883 : SynthesisImagerVi2 subImgr;
156 883 : subImgr.setCubeGridding(False);
157 1795 : for (Int k=0; k < Int(dataSel_p.nelements()); ++k){
158 : //The original SynthesisImager would have cleared the model if it was requested
159 912 : dataSel_p[k].incrmodel=True;
160 912 : dataSel_p[k].freqbeg="";
161 912 : subImgr.selectData(dataSel_p[k]);
162 : }
163 : //cerr <<"***Time for select data " << tim.real() << endl;
164 : //tim.mark();
165 883 : subImgr.setMovingSource(movingSource_p);
166 883 : Vector<CountedPtr<SIImageStore> > subImStor(imSel_p.nelements());
167 : //copy as shared_ptr as we mix the usage of countedptr and sharedptr
168 883 : Vector<std::shared_ptr<SIImageStore> > subImStorShared(imSel_p.nelements());
169 883 : Vector<SIImageStore *> p(imSel_p.nelements());
170 : //Do multifield in one process only for now
171 : //TODO if all fields have same nchan then partition all on all subcalls
172 883 : if(chanRange_p[0]==0){
173 1781 : for (uInt k=0; k < imSel_p.nelements(); ++k){
174 898 : subImStorShared[k]=subImageStore(k);
175 898 : subImStor[k]=CountedPtr<SIImageStore>(subImStorShared[k]);
176 :
177 898 : if(ftmRec_p.nelements()>0){
178 0 : subImgr.defineImage(subImStor[k], ftmRec_p[k], iftmRec_p[k]);
179 : }else{
180 898 : subImgr.defineImage(subImStor[k], imSel_p[k], gridSel_p[k]);
181 : }
182 : }
183 : }else{
184 0 : subImStor.resize(1);
185 0 : subImStorShared[0]=subImageStore(0);
186 0 : subImStor[0]=CountedPtr<SIImageStore>(subImStorShared[0]);
187 :
188 0 : if(ftmRec_p.nelements()>0){
189 0 : subImgr.defineImage(subImStor[0], ftmRec_p[0], iftmRec_p[0]);
190 : }else{
191 0 : subImgr.defineImage(subImStor[0], imSel_p[0], gridSel_p[0]);
192 : }
193 : }
194 883 : subImgr.setCubeGridding(False);
195 : // TO DO get weight param and set weight
196 : // Modified to fix CAS-13660 bug for natural weight setting with uvtaper
197 883 : if(!weightParams_p.isDefined("type") ){
198 0 : subImgr.weight("natural");
199 : }
200 : else{
201 883 : if(controlRecord_p.isDefined("weightdensity")){
202 15 : String densName=controlRecord_p.asString("weightdensity");
203 : //cerr << "Loading weightdensity " << densName << endl;
204 15 : if(Table::isReadable(densName))
205 15 : subImgr.setWeightDensity(densName);
206 15 : }
207 : else{
208 868 : subImgr.weight(weightParams_p);
209 : }
210 : }
211 : ///Now do the selection tuning if needed
212 883 : if(imSel_p[0].mode !="cubedata"){
213 : //cerr << "IN RETUNING " << endl;
214 877 : if(retuning_p)
215 877 : subImgr.tuneSelectData();
216 : }
217 :
218 : //cerr << "***Time for all other setting " << tim.real() << endl;
219 : //tim.mark();
220 883 : if (!dopsf_p){
221 : ///In case weightimages for mosaicft is done load it...we can get rid of this if we are using fromrecord ftm
222 : ///doing it for non-psf only ...psf divides it by sumwt for some reason in
223 : ///SIImageStore ..so restart with psf creates havoc
224 583 : subImgr.loadMosaicSensitivity();
225 583 : Record outrec=subImgr.executeMajorCycle(controlRecord_p);
226 583 : if(outrec.isDefined("tempfilenames")){
227 562 : returnRec_p.define("tempfilenames", outrec.asArrayString("tempfilenames"));
228 : }
229 1176 : for(uInt k=0; k < subImStor.nelements(); ++k){
230 593 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight")){
231 : {
232 593 : if(controlRecord_p.isDefined("normpars")){
233 593 : Record normpars=controlRecord_p.asRecord("normpars");
234 593 : SynthesisNormalizer norm;
235 593 : if(nterms_p[k] > 0){
236 2 : if(!normpars.isDefined("nterms"))
237 0 : normpars.define("nterms", uInt(nterms_p[k]));
238 2 : normpars.define("deconvolver", "mtmfs");
239 : }
240 :
241 593 : norm.setupNormalizer(normpars);
242 593 : norm.setImageStore(subImStorShared[k]);
243 593 : norm.divideResidualByWeight();
244 593 : }
245 : else{
246 0 : LatticeLocker lock1 (*(subImStor[k]->residual()), FileLocker::Write);
247 0 : LatticeLocker lock2 (*(subImStor[k]->sumwt()), FileLocker::Read);
248 0 : subImStor[k]->divideResidualByWeight();
249 : //subImStor[k]->residual()->flush();
250 0 : }
251 : }
252 593 : subImStor[k]->residual()->unlock();
253 593 : if(subImStor[k]->hasModel())
254 313 : subImStor[k]->model()->unlock();
255 : }
256 593 : Int chanBeg=0; Int chanEnd=0;
257 593 : if(k==0){
258 583 : chanBeg=chanRange_p[0];
259 583 : chanEnd=chanRange_p[1];
260 : }
261 : else{
262 10 : chanBeg=0;
263 10 : chanEnd=subImStor[k]->sumwt()->shape()[3]-1;
264 : }
265 593 : if(subImStor[k]->getType() != "multiterm"){
266 591 : writeBackToFullImage(residualNames_p[k], chanBeg, chanEnd, (subImStor[k]->residual()));
267 591 : writeBackToFullImage(sumwtNames_p[k], chanBeg, chanEnd, (subImStor[k]->sumwt()));
268 : }
269 :
270 : }
271 :
272 :
273 583 : }
274 : else{
275 300 : if((gridSel_p[0]).ftmachine=="awphpg"){
276 0 : if(!subImgr.loadMosaicSensitivity()){
277 0 : subImgr.makeMosaicSensitivity();
278 0 : writeBackToFullImage(weightNames_p[0], chanRange_p[0],chanRange_p[1], (subImStor[0]->weight()));
279 0 : writeBackToFullImage(sumwtNames_p[0], chanRange_p[0], chanRange_p[1], (subImStor[0]->sumwt()));
280 :
281 :
282 : }
283 0 : subImgr.loadMosaicSensitivity();
284 : }
285 300 : Record&& outrec=subImgr.makePSF();
286 299 : if(outrec.isDefined("tempfilenames")){
287 299 : returnRec_p.define("tempfilenames", outrec.asArrayString("tempfilenames"));
288 : }
289 : ////tclean expects a PB to be always there...
290 : //so for standard make it
291 :
292 299 : subImgr.makePB();
293 :
294 603 : for(uInt k=0; k < subImStor.nelements(); ++k){
295 304 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight"))
296 : {
297 304 : if(controlRecord_p.isDefined("normpars")){
298 304 : Record normpars=controlRecord_p.asRecord("normpars");
299 304 : if(nterms_p[k] > 0){
300 1 : if(!normpars.isDefined("nterms"))
301 0 : normpars.define("nterms", uInt(nterms_p[k]));
302 1 : normpars.define("deconvolver", "mtmfs");
303 : }
304 : // cerr << k << " " << nterms_p[k] <<" NORMPARS " << normpars << endl;
305 304 : SynthesisNormalizer norm;
306 304 : norm.setupNormalizer(normpars);
307 304 : norm.setImageStore(subImStorShared[k]);
308 304 : norm.dividePSFByWeight();
309 304 : norm.divideWeightBySumWt();
310 304 : copyBeamSet(*(subImStorShared[k]->psf()), k);
311 304 : }
312 : else{
313 0 : LatticeLocker lock1 (*(subImStor[k]->psf()), FileLocker::Write);
314 0 : LatticeLocker lock2 (*(subImStor[k]->sumwt()), FileLocker::Read);
315 0 : subImStor[k]->dividePSFByWeight();
316 0 : subImStor[k]->divideWeightBySumwt();
317 0 : copyBeamSet(*(subImStor[k]->psf()), k);
318 : //subImStor[k]->psf()->flush();
319 0 : }
320 : }
321 304 : Int chanBeg=0; Int chanEnd=0;
322 304 : if(k==0){
323 299 : chanBeg=chanRange_p[0];
324 299 : chanEnd=chanRange_p[1];
325 : }
326 : else{
327 5 : chanBeg=0;
328 5 : chanEnd=subImStor[k]->sumwt()->shape()[3]-1;
329 : }
330 304 : if(subImStor[k]->getType() != "multiterm"){
331 303 : writeBackToFullImage(psfNames_p[k], chanBeg, chanEnd, (subImStor[k]->psf()));
332 303 : writeBackToFullImage(sumwtNames_p[k], chanBeg, chanEnd, (subImStor[k]->sumwt()));
333 303 : if((subImStor[k]->hasSensitivity()) && Table::isWritable(weightNames_p[k]) && (gridSel_p[0]).ftmachine != "awphpg"){
334 :
335 72 : writeBackToFullImage(weightNames_p[k], chanBeg, chanEnd, (subImStor[k]->weight()));
336 :
337 : }
338 303 : if( doPB_p && Table::isWritable(pbNames_p[k])){
339 231 : writeBackToFullImage(pbNames_p[k], chanBeg, chanEnd, (subImStor[k]->pb()));
340 :
341 : }
342 :
343 : }
344 304 : subImStor[k]->psf()->unlock();
345 :
346 : }
347 299 : }
348 : //cerr << "***Time gridding/ffting " << tim.real() << endl;
349 882 : status_p = True;
350 886 : }
351 1 : catch (AipsError x) {
352 :
353 2 : LogIO os( LogOrigin("SynthesisImagerVi2","CubeMajorCycle",WHERE) );
354 1 : os << LogIO::WARN << "Exception for chan range " << chanRange_p << " --- "<< x.getMesg() << LogIO::POST;
355 1 : cerr << "##################################\n#############################\nException: " << x.getMesg() << endl;
356 :
357 1 : status_p=false;
358 1 : }
359 0 : catch(std::exception& exc) {
360 0 : logger << LogIO::SEVERE << "Exception (std): " << exc.what() << LogIO::POST;
361 0 : returnRec_p=Record();
362 0 : }
363 0 : catch(...){
364 0 : cerr << "###################################\n#######################3##\nUnknown exception " << endl;
365 0 : std::exception_ptr eptr=std::current_exception();
366 : try{
367 0 : if(eptr)
368 0 : std::rethrow_exception(eptr);
369 : }
370 0 : catch(const std::exception& a){
371 0 : LogIO os( LogOrigin("SynthesisImagerVi2","CubeMajorCycle",WHERE) );
372 0 : os << LogIO::WARN << "Unknown Exception for chan range " << chanRange_p << " --- "<< a.what() << LogIO::POST;
373 0 : }
374 :
375 0 : logger << LogIO::SEVERE << "Unknown exception " << LogIO::POST;
376 0 : status_p=False;
377 0 : }
378 883 : }
379 3536 : String& CubeMajorCycleAlgorithm::name(){
380 3536 : return myName_p;
381 : }
382 :
383 898 : shared_ptr<SIImageStore> CubeMajorCycleAlgorithm::subImageStore(const int imId){
384 : //For some reason multiterm deconvolver is allowed with cubes !
385 898 : String isMTdeconv="";
386 898 : nterms_p.resize(imId+1, True);
387 898 : nterms_p[imId]=-1;
388 898 : if(imId==0 && imSel_p[imId].deconvolver=="mtmfs"){
389 0 : isMTdeconv=".tt0";
390 0 : nterms_p[0]=1;
391 : }
392 :
393 898 : String residname=imSel_p[imId].imageName+".residual"+isMTdeconv;
394 898 : String psfname=imSel_p[imId].imageName+".psf"+isMTdeconv;
395 898 : String sumwgtname=imSel_p[imId].imageName+".sumwt"+isMTdeconv;
396 898 : if(imId > 0 && !Table::isReadable(sumwgtname)){
397 3 : return multiTermImageStore(imId);
398 : }
399 895 : shared_ptr<ImageInterface<Float> >subpsf=nullptr;
400 895 : shared_ptr<ImageInterface<Float> >subresid=nullptr;
401 895 : shared_ptr<ImageInterface<Float> >submodel=nullptr;
402 895 : shared_ptr<ImageInterface<Float> > subweight=nullptr;
403 895 : shared_ptr<ImageInterface<Float> > subpb=nullptr;
404 : //cerr << imId << " sumwt name " << sumwgtname << endl;
405 895 : String workingdir="";
406 895 : if(controlRecord_p.isDefined("workingdirectory"))
407 895 : controlRecord_p.get("workingdirectory", workingdir);
408 : //cerr <<"WORKINGDIR " << workingdir << endl;
409 895 : if(Table::isReadable(workingdir+"/"+sumwgtname)){
410 895 : workingdir=workingdir+"/";
411 895 : residname=workingdir+residname;
412 895 : psfname=workingdir+psfname;
413 895 : sumwgtname=workingdir+sumwgtname;
414 :
415 : }
416 0 : else if(!Table::isReadable(sumwgtname) )
417 0 : throw(AipsError("Programmer error: sumwt disk image is non existant"));
418 : else
419 0 : workingdir="";
420 : ///Use user locking to make sure locks are compliant
421 895 : PagedImage<Float> sumwt(sumwgtname, TableLock::UserLocking);
422 895 : sumwt.lock(FileLocker::Write, 200);
423 : //PagedImage<Float> sumwt(sumwgtname, TableLock::AutoNoReadLocking);
424 895 : if(sumwtNames_p.nelements() <= uInt(imId)){
425 895 : sumwtNames_p.resize(imId+1, True);
426 895 : psfNames_p.resize(imId+1, True);
427 895 : residualNames_p.resize(imId+1, True);
428 895 : sumwtNames_p[imId]=sumwgtname;
429 895 : residualNames_p[imId]=residname;
430 895 : psfNames_p[imId]=psfname;
431 : }
432 895 : Int nchannels=sumwt.shape()[3];
433 :
434 :
435 :
436 : //Should be partitioning for main image only
437 : //chanRange
438 895 : Int chanBeg=0;
439 895 : Int chanEnd=0;
440 895 : if(imId==0){
441 883 : chanBeg=chanRange_p[0];
442 883 : chanEnd=chanRange_p[1];
443 : }
444 : else{
445 12 : chanBeg=0;
446 12 : chanEnd=sumwt.shape()[3]-1;
447 : }
448 895 : sumwt.tempClose();
449 : ////For some small channel ms's retuning trigger a vi2/vb2 bug in nChannels
450 : ///avoid retuning for small images
451 : ////Skipping this here..
452 : //overloaded SynthesisImagerVi2::retune to do this check
453 : //if(nchannels < 30 && imId==0 && ((chanEnd-chanBeg) < 10)){
454 : // retuning_p=False;
455 : //}
456 : //cerr << "chanBeg " << chanBeg << " chanEnd " << chanEnd << " imId " << imId << endl;
457 895 : Vector<String> weightnams(controlRecord_p.asArrayString("weightnames"));
458 895 : Vector<String> pbnams(controlRecord_p.asArrayString("pbnames"));
459 895 : pbNames_p.resize();
460 895 : pbNames_p=pbnams;
461 895 : weightNames_p.resize();
462 895 : weightNames_p=weightnams;
463 :
464 895 : if(imId >= int(weightNames_p.nelements()))
465 0 : throw(AipsError("Number of weight images does not match number of image fields defined"));
466 895 : if(Table::isReadable(workingdir+weightNames_p[imId])){
467 0 : weightNames_p[imId]=workingdir+weightNames_p[imId];
468 : }
469 895 : if(Table::isReadable(workingdir+pbNames_p[imId])){
470 0 : pbNames_p[imId]=workingdir+pbNames_p[imId];
471 : }
472 :
473 895 : if(dopsf_p){
474 : //PagedImage<Float> psf(psfname, TableLock::UserNoReadLocking);
475 : //subpsf.reset(SpectralImageUtil::getChannel(psf, chanBeg, chanEnd, true));
476 304 : getSubImage(subpsf, chanBeg, chanEnd, psfname, False);
477 : //cerr << "PBNAMES " << pbNames_p << endl;
478 304 : if(Table::isReadable(pbNames_p[imId])){
479 :
480 232 : doPB_p=True;
481 232 : getSubImage(subpb, chanBeg, chanEnd, pbNames_p[imId], False);
482 : }
483 : }
484 : else{
485 : //need to loop over all fields somewhere
486 : //PagedImage<Float> resid(residname, TableLock::UserNoReadLocking);
487 : //subresid.reset(SpectralImageUtil::getChannel(resid, chanBeg, chanEnd, true));
488 591 : getSubImage(subresid, chanBeg, chanEnd, residname, False);
489 :
490 591 : if(controlRecord_p.isDefined("modelnames")){
491 312 : Vector<String> modelnames(controlRecord_p.asArrayString("modelnames"));
492 312 : if(imId >= int(modelnames.nelements()))
493 0 : throw(AipsError("Number of model images does not match number of image fields defined"));
494 312 : if(Table::isReadable(workingdir+modelnames[imId])){
495 312 : modelnames[imId]=workingdir+modelnames[imId];
496 : }
497 312 : if(Table::isReadable(modelnames[imId])){
498 :
499 :
500 : ///Pass some extra channels for interpolation while degridding..should match or be less than in SynthesisImager::tuneSelect
501 312 : Int startmodchan=(chanBeg-2) >0 ? chanBeg-2 : 0;
502 312 : Int endmodchan=(chanEnd+2) < nchannels ? chanEnd+2 : nchannels-1 ;
503 312 : if((gridSel_p[0]).ftmachine=="awphpg"){
504 : ////Only one channel is used
505 0 : startmodchan=chanBeg;
506 0 : endmodchan=chanEnd;
507 : }
508 :
509 : //cerr << "START END mod " << startmodchan << " " << endmodchan << endl;
510 : //Darn has to lock it as writable because overlap in SIMapperCollection code
511 : //wants that...though we are not really modifying it here
512 : //Bool writeisneeded=(imSel_p.nelements()!=1 || startmodel_p[imId].nelements() >0);
513 : //PagedImage<Float> model(modelnames[imId], TableLock::UserNoReadLocking);
514 : //submodel.reset(SpectralImageUtil::getChannel(model, startmodchan, endmodchan, writeisneeded));
515 312 : getSubImage(submodel, startmodchan, endmodchan, modelnames[imId], True);
516 : //
517 312 : if(Table::isReadable(weightNames_p[imId])){
518 92 : divideModelByWeight(imSel_p[imId].imageName, submodel, startmodchan, endmodchan, weightNames_p[imId]);
519 : }
520 : //ImageInterface<Float>* modim=new PagedImage<Float>(modelnames[imId], TableLock::UserNoReadLocking);
521 : //submodel.reset(modim);
522 :
523 : }
524 :
525 312 : }
526 :
527 : }
528 :
529 895 : if(Table::isReadable(weightNames_p[imId])){
530 : //PagedImage<Float> weight(weightnames[imId], TableLock::UserNoReadLocking);
531 : //subweight.reset(SpectralImageUtil::getChannel(weight, chanBeg, chanEnd, true));
532 224 : getSubImage(subweight, chanBeg, chanEnd, weightNames_p[imId], True);
533 : }
534 895 : shared_ptr<ImageInterface<Float> >subsumwt=nullptr;
535 : // subsumwt.reset(SpectralImageUtil::getChannel(sumwt, chanBeg, chanEnd, true));
536 895 : getSubImage(subsumwt, chanBeg, chanEnd, sumwgtname, True);
537 895 : bool useweightimage=(subweight) ? true : false;
538 1790 : shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imSel_p[imId].imageName, submodel, subresid, subpsf, subweight, nullptr, nullptr, subsumwt, nullptr, subpb, nullptr, nullptr, useweightimage));
539 895 : if(polRep_p[imId]< 0)
540 0 : throw(AipsError("data polarization type is not defined"));
541 895 : StokesImageUtil::PolRep polrep=(StokesImageUtil::PolRep)polRep_p[imId];
542 895 : subimstor->setDataPolFrame(polrep);
543 895 : if(startmodel_p[imId].nelements() >0){
544 3 : LatticeLocker lock1 (*(subimstor->model()), FileLocker::Write);
545 3 : subimstor->setModelImage(startmodel_p[imId]);
546 3 : }
547 : //cerr << "subimagestor TYPE" << subimstor->getType() << endl;
548 895 : return subimstor;
549 898 : }
550 :
551 :
552 3 : std::shared_ptr<SIImageStore> CubeMajorCycleAlgorithm::multiTermImageStore(const Int imId){
553 3 : uInt nterms=0;
554 6 : String sumwgtname=imSel_p[imId].imageName+".sumwt.tt"+String::toString(nterms);
555 12 : while (Table::isReadable(sumwgtname)){
556 9 : ++nterms;
557 9 : sumwgtname=imSel_p[imId].imageName+".sumwt.tt"+String::toString(nterms);
558 : }
559 3 : if(nterms==0){
560 0 : throw(AipsError("outlier "+String::toString(imId)+" field weight image is not defined"));
561 : }
562 3 : nterms=(nterms+1)/2;
563 3 : nterms_p[imId]=nterms;
564 3 : shared_ptr<SIImageStore> subimstor(new SIImageStoreMultiTerm(imSel_p[imId].imageName, nterms, True));
565 3 : if(polRep_p[imId]< 0)
566 0 : throw(AipsError("data polarization type is not defined"));
567 3 : StokesImageUtil::PolRep polrep=(StokesImageUtil::PolRep)polRep_p[imId];
568 3 : subimstor->setDataPolFrame(polrep);
569 6 : return subimstor;
570 3 : }
571 92 : void CubeMajorCycleAlgorithm::divideModelByWeight(casacore::String imageName, shared_ptr<ImageInterface<Float> >&submodel, const Int startchan, const Int endchan, const String weightname){
572 92 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight")){
573 92 : if(controlRecord_p.isDefined("normpars")){
574 92 : Record normpars=controlRecord_p.asRecord("normpars");
575 92 : SynthesisNormalizer norm;
576 92 : norm.setupNormalizer(normpars);
577 92 : shared_ptr<ImageInterface<Float> > subweight=nullptr;
578 92 : getSubImage(subweight, startchan, endchan, weightname, True);
579 92 : LatticeLocker lock1(*(subweight), FileLocker::Read);
580 92 : LatticeLocker lock2(*(submodel), FileLocker::Read);
581 184 : std::shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imageName, submodel, nullptr, nullptr, subweight, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, True));
582 92 : norm.setImageStore(subimstor);
583 92 : norm.divideModelByWeight();
584 :
585 92 : }
586 : }
587 92 : }
588 :
589 :
590 :
591 :
592 883 : void CubeMajorCycleAlgorithm::reset(){
593 :
594 883 : dataSel_p.resize();
595 883 : imSel_p.resize();
596 883 : gridSel_p.resize();
597 883 : ftmRec_p.resize();
598 883 : iftmRec_p.resize();
599 883 : polRep_p.resize();
600 883 : chanRange_p.resize();
601 883 : dopsf_p=False;
602 883 : controlRecord_p=Record();
603 883 : weightParams_p=Record();
604 883 : returnRec_p=Record();
605 883 : startmodel_p.resize();
606 883 : status_p=False;
607 883 : residualNames_p.resize();
608 883 : psfNames_p.resize();
609 883 : sumwtNames_p.resize();
610 883 : weightNames_p.resize();
611 883 : pbNames_p.resize();
612 883 : movingSource_p="";
613 883 : retuning_p=True;
614 883 : doPB_p=False;
615 883 : nterms_p.resize();
616 :
617 :
618 :
619 883 : }
620 :
621 2650 : void CubeMajorCycleAlgorithm::getSubImage(std::shared_ptr<ImageInterface<Float> >& subimptr, const Int chanBeg, const Int chanEnd, const String imagename, const Bool copy){
622 5300 : LogIO os( LogOrigin("CubeMajorCycle","getSubImage",WHERE) );
623 2650 : CountedPtr<PagedImage<Float> > im=nullptr;
624 : {///work around
625 2650 : uInt exceptCounter=0;
626 : while(true){
627 : try{
628 2650 : im = new PagedImage<Float> (imagename, TableLock::UserLocking);
629 : //PagedImage<Float> im(imagename, TableLock::AutoNoReadLocking);
630 2650 : im->lock(FileLocker::Read, 1000);
631 2650 : SubImage<Float> *tmpptr=nullptr;
632 2650 : tmpptr=SpectralImageUtil::getChannel(*im, chanBeg, chanEnd, false);
633 2650 : subimptr.reset(new TempImage<Float>(TiledShape(tmpptr->shape(), tmpptr->niceCursorShape()), tmpptr->coordinates()));
634 :
635 2650 : if(copy)
636 1523 : subimptr->copyData(*tmpptr);
637 : //subimptr->flush();
638 2650 : im->unlock();
639 2650 : delete tmpptr;
640 2650 : break; ///get out of while loop
641 : }
642 0 : catch(AipsError &x){
643 0 : im=nullptr;
644 0 : String mes=x.getMesg();
645 0 : if(mes.contains("FilebufIO::readBlock") ){
646 0 : std::this_thread::sleep_for(std::chrono::milliseconds(50));
647 0 : os << LogIO::WARN << "#####CATCHING a sleep because "<< mes<< LogIO::POST;
648 : }
649 : else
650 0 : throw(AipsError("Error in readimage "+imagename+" : "+mes));
651 :
652 0 : if(exceptCounter > 100){
653 0 : throw(AipsError("Error in readimage got 100 of this exeception: "+mes));
654 :
655 : }
656 :
657 0 : }
658 0 : ++exceptCounter;
659 0 : }
660 : }
661 :
662 2650 : }
663 :
664 2091 : void CubeMajorCycleAlgorithm::writeBackToFullImage(const String imagename, const Int chanBeg, const Int chanEnd, std::shared_ptr<ImageInterface<Float> > subimptr){
665 4182 : LogIO os( LogOrigin("CubeMajorCycle","writebacktofullimage",WHERE) );
666 2091 : CountedPtr<PagedImage<Float> > im=nullptr;
667 : {///work around
668 2091 : uInt exceptCounter=0;
669 : while(true){
670 : try{
671 2091 : im=new PagedImage<Float> (imagename, TableLock::UserLocking);
672 :
673 :
674 : //PagedImage<Float> im(imagename, TableLock::AutoLocking);
675 2091 : SubImage<Float> *tmpptr=nullptr;
676 : {
677 2091 : tmpptr=SpectralImageUtil::getChannel(*im, chanBeg, chanEnd, true);
678 2091 : LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
679 2091 : tmpptr->set(0.0);
680 2091 : tmpptr->copyData(*(subimptr));
681 2091 : }
682 2091 : im->flush();
683 2091 : im->unlock();
684 2091 : delete tmpptr;
685 2091 : break; //get out of while loop
686 : }
687 0 : catch(AipsError &x){
688 :
689 0 : String mes=x.getMesg();
690 0 : if(mes.contains("FilebufIO::readBlock") ){
691 0 : std::this_thread::sleep_for(std::chrono::milliseconds(50));
692 0 : os << LogIO::WARN << "#####CATCHING a sleep because "<< mes<< LogIO::POST;
693 : }
694 : else
695 0 : throw(AipsError("Error in writing back image "+imagename+" : "+mes));
696 :
697 0 : if(exceptCounter > 100){
698 0 : throw(AipsError("Error in writeimage got 100 of this exeception: "+mes));
699 :
700 : }
701 :
702 0 : }
703 0 : ++exceptCounter;
704 0 : }
705 : }//End of work around for table disappearing bug
706 :
707 :
708 2091 : }
709 304 : void CubeMajorCycleAlgorithm::copyBeamSet(ImageInterface<Float>& subpsf, const Int imageid){
710 : //For now lets do for the main image only
711 304 : if(imageid != 0)
712 5 : return;
713 299 : ImageBeamSet iibeamset=subpsf.imageInfo().getBeamSet();
714 299 : uInt nchan=chanRange_p[1]-chanRange_p[0]+1;
715 299 : if(nchan ==iibeamset.nchan()){
716 299 : Matrix<GaussianBeam> gbs=iibeamset.getBeams();
717 299 : Cube<Float> beams(nchan, iibeamset.nstokes(),3);
718 :
719 2614 : for (uInt k=0; k < nchan; ++k){
720 4849 : for (uInt j=0; j < iibeamset.nstokes(); ++j){
721 2534 : beams(k,j, 0)=gbs(k, 0).getMajor("arcsec");
722 2534 : beams(k,j, 1)=gbs(k, 0).getMinor("arcsec");
723 2534 : beams(k,j, 2)=gbs(k, 0).getPA("deg", True);
724 : }
725 : }
726 299 : returnRec_p.define("imageid", imageid);
727 299 : returnRec_p.define("chanrange", chanRange_p);
728 299 : returnRec_p.define("beams", beams);
729 :
730 :
731 299 : }
732 :
733 :
734 :
735 299 : }
736 :
737 :
738 :
739 : } //# NAMESPACE CASA - END
|
