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