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 22 : 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 22 : }
44 23 : CubeMajorCycleAlgorithm::~CubeMajorCycleAlgorithm() {
45 :
46 23 : }
47 :
48 21 : void CubeMajorCycleAlgorithm::get() {
49 21 : reset();
50 : //cerr << "in get for child process " << applicator.isWorker() << endl;
51 21 : Record vecImParsRec;
52 21 : Record vecSelParsRec;
53 21 : Record vecGridParsRec;
54 : // get data sel params #1
55 21 : applicator.get(vecSelParsRec);
56 : // get image sel params #2
57 21 : applicator.get(vecImParsRec);
58 : // get gridders params #3
59 21 : applicator.get(vecGridParsRec);
60 : // get which channel to process #4
61 21 : applicator.get(chanRange_p);
62 : //cerr <<"GET chanRange " << chanRange_p << endl;
63 : // psf or residual CubeMajorCycleAlgorithm #5
64 21 : applicator.get(dopsf_p);
65 : // store modelvis and other controls #6
66 21 : applicator.get(controlRecord_p);
67 : // weight params
68 21 : applicator.get(weightParams_p);
69 : //Somewhere before this we have to make sure that vecSelPars has more than 0 fields)
70 21 : dataSel_p.resize(vecSelParsRec.nfields());
71 : /// Fill the private variables
72 42 : for (Int k=0; k < Int(dataSel_p.nelements()); ++k){
73 21 : (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 21 : Int nmajorcycles=0;
78 21 : if(controlRecord_p.isDefined("nmajorcycles"))
79 15 : controlRecord_p.get("nmajorcycles",nmajorcycles);
80 21 : imSel_p.resize(vecImParsRec.nfields());
81 21 : gridSel_p.resize(vecImParsRec.nfields());
82 21 : ftmRec_p.resize(vecImParsRec.nfields());
83 21 : iftmRec_p.resize(vecImParsRec.nfields());
84 21 : polRep_p.resize(vecImParsRec.nfields());
85 21 : polRep_p.set(-1);
86 21 : startmodel_p.resize(vecImParsRec.nfields());
87 21 : startmodel_p.set(Vector<String>(0));
88 : //TESTOO
89 : //Int CPUID;
90 : //MPI_Comm_rank(MPI_COMM_WORLD, &CPUID);
91 : //////////////////
92 42 : for (uInt k=0; k < imSel_p.nelements(); ++k){
93 42 : Record imSelRec=vecImParsRec.asRecord(String::toString(k));
94 : //cerr << k << " imsel " << imSelRec << endl;
95 21 : if(imSelRec.isDefined("polrep"))
96 21 : imSelRec.get("polrep", polRep_p[k]);
97 : ///Get moving source name
98 21 : if(imSelRec.isDefined("movingsource") && imSelRec.asString("movingsource") != ""){
99 0 : imSelRec.get("movingsource", movingSource_p);
100 : }
101 : //cerr << CPUID << " POLREP " << polRep_p << endl;
102 : //Only first major cycle we need to reset model
103 21 : if(nmajorcycles==1)
104 6 : 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 21 : imSelRec.define("startmodel", Vector<String>(0));
107 21 : (imSel_p[k]).fromRecord(imSelRec);
108 42 : Record recGridSel=vecGridParsRec.asRecord(String::toString(k));
109 21 : (gridSel_p[k]).fromRecord(recGridSel);
110 21 : if(!recGridSel.isDefined("ftmachine")){
111 21 : ftmRec_p.resize();
112 21 : iftmRec_p.resize();
113 : }
114 21 : if(ftmRec_p.nelements() >0){
115 0 : ftmRec_p[k]=recGridSel.asRecord("ftmachine");
116 0 : iftmRec_p[k]=recGridSel.asRecord("iftmachine");
117 : }
118 :
119 21 : }
120 :
121 :
122 21 : }
123 21 : 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 21 : applicator.put(returnRec_p);
132 21 : applicator.put(status_p);
133 :
134 21 : }
135 :
136 21 : void CubeMajorCycleAlgorithm::task(){
137 42 : LogIO logger(LogOrigin("CubeMajorCycleAlgorithm", "task", WHERE));
138 21 : 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 42 : for (uInt k=0; k < gridSel_p.nelements(); ++k)
146 21 : 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 21 : SynthesisImagerVi2 subImgr;
156 21 : subImgr.setCubeGridding(False);
157 42 : for (Int k=0; k < Int(dataSel_p.nelements()); ++k){
158 : //The original SynthesisImager would have cleared the model if it was requested
159 21 : dataSel_p[k].incrmodel=True;
160 21 : dataSel_p[k].freqbeg="";
161 21 : subImgr.selectData(dataSel_p[k]);
162 : }
163 : //cerr <<"***Time for select data " << tim.real() << endl;
164 : //tim.mark();
165 21 : subImgr.setMovingSource(movingSource_p);
166 21 : Vector<CountedPtr<SIImageStore> > subImStor(imSel_p.nelements());
167 : //copy as shared_ptr as we mix the usage of countedptr and sharedptr
168 21 : Vector<std::shared_ptr<SIImageStore> > subImStorShared(imSel_p.nelements());
169 21 : 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 21 : if(chanRange_p[0]==0){
173 42 : for (uInt k=0; k < imSel_p.nelements(); ++k){
174 21 : subImStorShared[k]=subImageStore(k);
175 21 : subImStor[k]=CountedPtr<SIImageStore>(subImStorShared[k]);
176 :
177 21 : if(ftmRec_p.nelements()>0){
178 0 : subImgr.defineImage(subImStor[k], ftmRec_p[k], iftmRec_p[k]);
179 : }else{
180 21 : 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 21 : 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 21 : if(!weightParams_p.isDefined("type") ){
198 0 : subImgr.weight("natural");
199 : }
200 : else{
201 21 : if(controlRecord_p.isDefined("weightdensity")){
202 0 : String densName=controlRecord_p.asString("weightdensity");
203 : //cerr << "Loading weightdensity " << densName << endl;
204 0 : if(Table::isReadable(densName))
205 0 : subImgr.setWeightDensity(densName);
206 0 : }
207 : else{
208 21 : subImgr.weight(weightParams_p);
209 : }
210 : }
211 : ///Now do the selection tuning if needed
212 21 : if(imSel_p[0].mode !="cubedata"){
213 : //cerr << "IN RETUNING " << endl;
214 21 : if(retuning_p)
215 21 : subImgr.tuneSelectData();
216 : }
217 :
218 : //cerr << "***Time for all other setting " << tim.real() << endl;
219 : //tim.mark();
220 21 : 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 15 : subImgr.loadMosaicSensitivity();
225 15 : Record outrec=subImgr.executeMajorCycle(controlRecord_p);
226 15 : if(outrec.isDefined("tempfilenames")){
227 12 : returnRec_p.define("tempfilenames", outrec.asArrayString("tempfilenames"));
228 : }
229 30 : for(uInt k=0; k < subImStor.nelements(); ++k){
230 15 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight")){
231 : {
232 15 : if(controlRecord_p.isDefined("normpars")){
233 15 : Record normpars=controlRecord_p.asRecord("normpars");
234 15 : SynthesisNormalizer norm;
235 15 : if(nterms_p[k] > 0){
236 0 : if(!normpars.isDefined("nterms"))
237 0 : normpars.define("nterms", uInt(nterms_p[k]));
238 0 : normpars.define("deconvolver", "mtmfs");
239 : }
240 :
241 15 : norm.setupNormalizer(normpars);
242 15 : norm.setImageStore(subImStorShared[k]);
243 15 : norm.divideResidualByWeight();
244 15 : }
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 15 : subImStor[k]->residual()->unlock();
253 15 : if(subImStor[k]->hasModel())
254 9 : subImStor[k]->model()->unlock();
255 : }
256 15 : Int chanBeg=0; Int chanEnd=0;
257 15 : if(k==0){
258 15 : chanBeg=chanRange_p[0];
259 15 : chanEnd=chanRange_p[1];
260 : }
261 : else{
262 0 : chanBeg=0;
263 0 : chanEnd=subImStor[k]->sumwt()->shape()[3]-1;
264 : }
265 15 : if(subImStor[k]->getType() != "multiterm"){
266 15 : writeBackToFullImage(residualNames_p[k], chanBeg, chanEnd, (subImStor[k]->residual()));
267 15 : writeBackToFullImage(sumwtNames_p[k], chanBeg, chanEnd, (subImStor[k]->sumwt()));
268 : }
269 :
270 : }
271 :
272 :
273 15 : }
274 : else{
275 6 : 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 6 : Record&& outrec=subImgr.makePSF();
286 6 : if(outrec.isDefined("tempfilenames")){
287 6 : 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 6 : subImgr.makePB();
293 :
294 12 : for(uInt k=0; k < subImStor.nelements(); ++k){
295 6 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight"))
296 : {
297 6 : if(controlRecord_p.isDefined("normpars")){
298 6 : Record normpars=controlRecord_p.asRecord("normpars");
299 6 : if(nterms_p[k] > 0){
300 0 : if(!normpars.isDefined("nterms"))
301 0 : normpars.define("nterms", uInt(nterms_p[k]));
302 0 : normpars.define("deconvolver", "mtmfs");
303 : }
304 : // cerr << k << " " << nterms_p[k] <<" NORMPARS " << normpars << endl;
305 6 : SynthesisNormalizer norm;
306 6 : norm.setupNormalizer(normpars);
307 6 : norm.setImageStore(subImStorShared[k]);
308 6 : norm.dividePSFByWeight();
309 6 : norm.divideWeightBySumWt();
310 6 : copyBeamSet(*(subImStorShared[k]->psf()), k);
311 6 : }
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 6 : Int chanBeg=0; Int chanEnd=0;
322 6 : if(k==0){
323 6 : chanBeg=chanRange_p[0];
324 6 : chanEnd=chanRange_p[1];
325 : }
326 : else{
327 0 : chanBeg=0;
328 0 : chanEnd=subImStor[k]->sumwt()->shape()[3]-1;
329 : }
330 6 : if(subImStor[k]->getType() != "multiterm"){
331 6 : writeBackToFullImage(psfNames_p[k], chanBeg, chanEnd, (subImStor[k]->psf()));
332 6 : writeBackToFullImage(sumwtNames_p[k], chanBeg, chanEnd, (subImStor[k]->sumwt()));
333 6 : if((subImStor[k]->hasSensitivity()) && Table::isWritable(weightNames_p[k]) && (gridSel_p[0]).ftmachine != "awphpg"){
334 :
335 3 : writeBackToFullImage(weightNames_p[k], chanBeg, chanEnd, (subImStor[k]->weight()));
336 :
337 : }
338 6 : if( doPB_p && Table::isWritable(pbNames_p[k])){
339 3 : writeBackToFullImage(pbNames_p[k], chanBeg, chanEnd, (subImStor[k]->pb()));
340 :
341 : }
342 :
343 : }
344 6 : subImStor[k]->psf()->unlock();
345 :
346 : }
347 6 : }
348 : //cerr << "***Time gridding/ffting " << tim.real() << endl;
349 21 : status_p = True;
350 21 : }
351 0 : catch (AipsError x) {
352 :
353 0 : LogIO os( LogOrigin("SynthesisImagerVi2","CubeMajorCycle",WHERE) );
354 0 : os << LogIO::WARN << "Exception for chan range " << chanRange_p << " --- "<< x.getMesg() << LogIO::POST;
355 0 : cerr << "##################################\n#############################\nException: " << x.getMesg() << endl;
356 :
357 0 : status_p=false;
358 0 : }
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 21 : }
379 85 : String& CubeMajorCycleAlgorithm::name(){
380 85 : return myName_p;
381 : }
382 :
383 21 : shared_ptr<SIImageStore> CubeMajorCycleAlgorithm::subImageStore(const int imId){
384 : //For some reason multiterm deconvolver is allowed with cubes !
385 21 : String isMTdeconv="";
386 21 : nterms_p.resize(imId+1, True);
387 21 : nterms_p[imId]=-1;
388 21 : if(imId==0 && imSel_p[imId].deconvolver=="mtmfs"){
389 0 : isMTdeconv=".tt0";
390 0 : nterms_p[0]=1;
391 : }
392 :
393 21 : String residname=imSel_p[imId].imageName+".residual"+isMTdeconv;
394 21 : String psfname=imSel_p[imId].imageName+".psf"+isMTdeconv;
395 21 : String sumwgtname=imSel_p[imId].imageName+".sumwt"+isMTdeconv;
396 21 : if(imId > 0 && !Table::isReadable(sumwgtname)){
397 0 : return multiTermImageStore(imId);
398 : }
399 21 : shared_ptr<ImageInterface<Float> >subpsf=nullptr;
400 21 : shared_ptr<ImageInterface<Float> >subresid=nullptr;
401 21 : shared_ptr<ImageInterface<Float> >submodel=nullptr;
402 21 : shared_ptr<ImageInterface<Float> > subweight=nullptr;
403 21 : shared_ptr<ImageInterface<Float> > subpb=nullptr;
404 : //cerr << imId << " sumwt name " << sumwgtname << endl;
405 21 : String workingdir="";
406 21 : if(controlRecord_p.isDefined("workingdirectory"))
407 21 : controlRecord_p.get("workingdirectory", workingdir);
408 : //cerr <<"WORKINGDIR " << workingdir << endl;
409 21 : if(Table::isReadable(workingdir+"/"+sumwgtname)){
410 21 : workingdir=workingdir+"/";
411 21 : residname=workingdir+residname;
412 21 : psfname=workingdir+psfname;
413 21 : 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 21 : PagedImage<Float> sumwt(sumwgtname, TableLock::UserLocking);
422 21 : sumwt.lock(FileLocker::Write, 200);
423 : //PagedImage<Float> sumwt(sumwgtname, TableLock::AutoNoReadLocking);
424 21 : if(sumwtNames_p.nelements() <= uInt(imId)){
425 21 : sumwtNames_p.resize(imId+1, True);
426 21 : psfNames_p.resize(imId+1, True);
427 21 : residualNames_p.resize(imId+1, True);
428 21 : sumwtNames_p[imId]=sumwgtname;
429 21 : residualNames_p[imId]=residname;
430 21 : psfNames_p[imId]=psfname;
431 : }
432 21 : Int nchannels=sumwt.shape()[3];
433 :
434 :
435 :
436 : //Should be partitioning for main image only
437 : //chanRange
438 21 : Int chanBeg=0;
439 21 : Int chanEnd=0;
440 21 : if(imId==0){
441 21 : chanBeg=chanRange_p[0];
442 21 : chanEnd=chanRange_p[1];
443 : }
444 : else{
445 0 : chanBeg=0;
446 0 : chanEnd=sumwt.shape()[3]-1;
447 : }
448 21 : 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 21 : Vector<String> weightnams(controlRecord_p.asArrayString("weightnames"));
458 21 : Vector<String> pbnams(controlRecord_p.asArrayString("pbnames"));
459 21 : pbNames_p.resize();
460 21 : pbNames_p=pbnams;
461 21 : weightNames_p.resize();
462 21 : weightNames_p=weightnams;
463 :
464 21 : if(imId >= int(weightNames_p.nelements()))
465 0 : throw(AipsError("Number of weight images does not match number of image fields defined"));
466 21 : if(Table::isReadable(workingdir+weightNames_p[imId])){
467 0 : weightNames_p[imId]=workingdir+weightNames_p[imId];
468 : }
469 21 : if(Table::isReadable(workingdir+pbNames_p[imId])){
470 0 : pbNames_p[imId]=workingdir+pbNames_p[imId];
471 : }
472 :
473 21 : if(dopsf_p){
474 : //PagedImage<Float> psf(psfname, TableLock::UserNoReadLocking);
475 : //subpsf.reset(SpectralImageUtil::getChannel(psf, chanBeg, chanEnd, true));
476 6 : getSubImage(subpsf, chanBeg, chanEnd, psfname, False);
477 : //cerr << "PBNAMES " << pbNames_p << endl;
478 6 : if(Table::isReadable(pbNames_p[imId])){
479 :
480 3 : doPB_p=True;
481 3 : 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 15 : getSubImage(subresid, chanBeg, chanEnd, residname, False);
489 :
490 15 : if(controlRecord_p.isDefined("modelnames")){
491 9 : Vector<String> modelnames(controlRecord_p.asArrayString("modelnames"));
492 9 : if(imId >= int(modelnames.nelements()))
493 0 : throw(AipsError("Number of model images does not match number of image fields defined"));
494 9 : if(Table::isReadable(workingdir+modelnames[imId])){
495 9 : modelnames[imId]=workingdir+modelnames[imId];
496 : }
497 9 : 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 9 : Int startmodchan=(chanBeg-2) >0 ? chanBeg-2 : 0;
502 9 : Int endmodchan=(chanEnd+2) < nchannels ? chanEnd+2 : nchannels-1 ;
503 9 : 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 9 : getSubImage(submodel, startmodchan, endmodchan, modelnames[imId], True);
516 : //
517 9 : if(Table::isReadable(weightNames_p[imId])){
518 6 : 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 9 : }
526 :
527 : }
528 :
529 21 : if(Table::isReadable(weightNames_p[imId])){
530 : //PagedImage<Float> weight(weightnames[imId], TableLock::UserNoReadLocking);
531 : //subweight.reset(SpectralImageUtil::getChannel(weight, chanBeg, chanEnd, true));
532 12 : getSubImage(subweight, chanBeg, chanEnd, weightNames_p[imId], True);
533 : }
534 21 : shared_ptr<ImageInterface<Float> >subsumwt=nullptr;
535 : // subsumwt.reset(SpectralImageUtil::getChannel(sumwt, chanBeg, chanEnd, true));
536 21 : getSubImage(subsumwt, chanBeg, chanEnd, sumwgtname, True);
537 21 : bool useweightimage=(subweight) ? true : false;
538 42 : shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imSel_p[imId].imageName, submodel, subresid, subpsf, subweight, nullptr, nullptr, subsumwt, nullptr, subpb, nullptr, nullptr, useweightimage));
539 21 : if(polRep_p[imId]< 0)
540 0 : throw(AipsError("data polarization type is not defined"));
541 21 : StokesImageUtil::PolRep polrep=(StokesImageUtil::PolRep)polRep_p[imId];
542 21 : subimstor->setDataPolFrame(polrep);
543 21 : if(startmodel_p[imId].nelements() >0){
544 0 : LatticeLocker lock1 (*(subimstor->model()), FileLocker::Write);
545 0 : subimstor->setModelImage(startmodel_p[imId]);
546 0 : }
547 : //cerr << "subimagestor TYPE" << subimstor->getType() << endl;
548 21 : return subimstor;
549 21 : }
550 :
551 :
552 0 : std::shared_ptr<SIImageStore> CubeMajorCycleAlgorithm::multiTermImageStore(const Int imId){
553 0 : uInt nterms=0;
554 0 : String sumwgtname=imSel_p[imId].imageName+".sumwt.tt"+String::toString(nterms);
555 0 : while (Table::isReadable(sumwgtname)){
556 0 : ++nterms;
557 0 : sumwgtname=imSel_p[imId].imageName+".sumwt.tt"+String::toString(nterms);
558 : }
559 0 : if(nterms==0){
560 0 : throw(AipsError("outlier "+String::toString(imId)+" field weight image is not defined"));
561 : }
562 0 : nterms=(nterms+1)/2;
563 0 : nterms_p[imId]=nterms;
564 0 : shared_ptr<SIImageStore> subimstor(new SIImageStoreMultiTerm(imSel_p[imId].imageName, nterms, True));
565 0 : if(polRep_p[imId]< 0)
566 0 : throw(AipsError("data polarization type is not defined"));
567 0 : StokesImageUtil::PolRep polrep=(StokesImageUtil::PolRep)polRep_p[imId];
568 0 : subimstor->setDataPolFrame(polrep);
569 0 : return subimstor;
570 0 : }
571 6 : void CubeMajorCycleAlgorithm::divideModelByWeight(casacore::String imageName, shared_ptr<ImageInterface<Float> >&submodel, const Int startchan, const Int endchan, const String weightname){
572 6 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight")){
573 6 : if(controlRecord_p.isDefined("normpars")){
574 6 : Record normpars=controlRecord_p.asRecord("normpars");
575 6 : SynthesisNormalizer norm;
576 6 : norm.setupNormalizer(normpars);
577 6 : shared_ptr<ImageInterface<Float> > subweight=nullptr;
578 6 : getSubImage(subweight, startchan, endchan, weightname, True);
579 6 : LatticeLocker lock1(*(subweight), FileLocker::Read);
580 6 : LatticeLocker lock2(*(submodel), FileLocker::Read);
581 12 : std::shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imageName, submodel, nullptr, nullptr, subweight, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, True));
582 6 : norm.setImageStore(subimstor);
583 6 : norm.divideModelByWeight();
584 :
585 6 : }
586 : }
587 6 : }
588 :
589 :
590 :
591 :
592 21 : void CubeMajorCycleAlgorithm::reset(){
593 :
594 21 : dataSel_p.resize();
595 21 : imSel_p.resize();
596 21 : gridSel_p.resize();
597 21 : ftmRec_p.resize();
598 21 : iftmRec_p.resize();
599 21 : polRep_p.resize();
600 21 : chanRange_p.resize();
601 21 : dopsf_p=False;
602 21 : controlRecord_p=Record();
603 21 : weightParams_p=Record();
604 21 : returnRec_p=Record();
605 21 : startmodel_p.resize();
606 21 : status_p=False;
607 21 : residualNames_p.resize();
608 21 : psfNames_p.resize();
609 21 : sumwtNames_p.resize();
610 21 : weightNames_p.resize();
611 21 : pbNames_p.resize();
612 21 : movingSource_p="";
613 21 : retuning_p=True;
614 21 : doPB_p=False;
615 21 : nterms_p.resize();
616 :
617 :
618 :
619 21 : }
620 :
621 72 : void CubeMajorCycleAlgorithm::getSubImage(std::shared_ptr<ImageInterface<Float> >& subimptr, const Int chanBeg, const Int chanEnd, const String imagename, const Bool copy){
622 144 : LogIO os( LogOrigin("CubeMajorCycle","getSubImage",WHERE) );
623 72 : CountedPtr<PagedImage<Float> > im=nullptr;
624 : {///work around
625 72 : uInt exceptCounter=0;
626 : while(true){
627 : try{
628 72 : im = new PagedImage<Float> (imagename, TableLock::UserLocking);
629 : //PagedImage<Float> im(imagename, TableLock::AutoNoReadLocking);
630 72 : im->lock(FileLocker::Read, 1000);
631 72 : SubImage<Float> *tmpptr=nullptr;
632 72 : tmpptr=SpectralImageUtil::getChannel(*im, chanBeg, chanEnd, false);
633 72 : subimptr.reset(new TempImage<Float>(TiledShape(tmpptr->shape(), tmpptr->niceCursorShape()), tmpptr->coordinates()));
634 :
635 72 : if(copy)
636 48 : subimptr->copyData(*tmpptr);
637 : //subimptr->flush();
638 72 : im->unlock();
639 72 : delete tmpptr;
640 72 : 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 72 : }
663 :
664 48 : void CubeMajorCycleAlgorithm::writeBackToFullImage(const String imagename, const Int chanBeg, const Int chanEnd, std::shared_ptr<ImageInterface<Float> > subimptr){
665 96 : LogIO os( LogOrigin("CubeMajorCycle","writebacktofullimage",WHERE) );
666 48 : CountedPtr<PagedImage<Float> > im=nullptr;
667 : {///work around
668 48 : uInt exceptCounter=0;
669 : while(true){
670 : try{
671 48 : im=new PagedImage<Float> (imagename, TableLock::UserLocking);
672 :
673 :
674 : //PagedImage<Float> im(imagename, TableLock::AutoLocking);
675 48 : SubImage<Float> *tmpptr=nullptr;
676 : {
677 48 : tmpptr=SpectralImageUtil::getChannel(*im, chanBeg, chanEnd, true);
678 48 : LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
679 48 : tmpptr->set(0.0);
680 48 : tmpptr->copyData(*(subimptr));
681 48 : }
682 48 : im->flush();
683 48 : im->unlock();
684 48 : delete tmpptr;
685 48 : 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 48 : }
709 6 : void CubeMajorCycleAlgorithm::copyBeamSet(ImageInterface<Float>& subpsf, const Int imageid){
710 : //For now lets do for the main image only
711 6 : if(imageid != 0)
712 0 : return;
713 6 : ImageBeamSet iibeamset=subpsf.imageInfo().getBeamSet();
714 6 : uInt nchan=chanRange_p[1]-chanRange_p[0]+1;
715 6 : if(nchan ==iibeamset.nchan()){
716 6 : Matrix<GaussianBeam> gbs=iibeamset.getBeams();
717 6 : Cube<Float> beams(nchan, iibeamset.nstokes(),3);
718 :
719 21 : for (uInt k=0; k < nchan; ++k){
720 30 : for (uInt j=0; j < iibeamset.nstokes(); ++j){
721 15 : beams(k,j, 0)=gbs(k, 0).getMajor("arcsec");
722 15 : beams(k,j, 1)=gbs(k, 0).getMinor("arcsec");
723 15 : beams(k,j, 2)=gbs(k, 0).getPA("deg", True);
724 : }
725 : }
726 6 : returnRec_p.define("imageid", imageid);
727 6 : returnRec_p.define("chanrange", chanRange_p);
728 6 : returnRec_p.define("beams", beams);
729 :
730 :
731 6 : }
732 :
733 :
734 :
735 6 : }
736 :
737 :
738 :
739 : } //# NAMESPACE CASA - END
|
