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 : Record&& outrec=subImgr.makePSF();
276 0 : if(outrec.isDefined("tempfilenames")){
277 0 : returnRec_p.define("tempfilenames", outrec.asArrayString("tempfilenames"));
278 : }
279 : ////tclean expects a PB to be always there...
280 : //so for standard make it
281 0 : subImgr.makePB();
282 0 : for(uInt k=0; k < subImStor.nelements(); ++k){
283 0 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight"))
284 : {
285 0 : if(controlRecord_p.isDefined("normpars")){
286 0 : Record normpars=controlRecord_p.asRecord("normpars");
287 0 : if(nterms_p[k] > 0){
288 0 : if(!normpars.isDefined("nterms"))
289 0 : normpars.define("nterms", uInt(nterms_p[k]));
290 0 : normpars.define("deconvolver", "mtmfs");
291 : }
292 : // cerr << k << " " << nterms_p[k] <<" NORMPARS " << normpars << endl;
293 0 : SynthesisNormalizer norm;
294 0 : norm.setupNormalizer(normpars);
295 0 : norm.setImageStore(subImStorShared[k]);
296 0 : norm.dividePSFByWeight();
297 0 : copyBeamSet(*(subImStorShared[k]->psf()), k);
298 0 : }
299 : else{
300 0 : LatticeLocker lock1 (*(subImStor[k]->psf()), FileLocker::Write);
301 0 : LatticeLocker lock2 (*(subImStor[k]->sumwt()), FileLocker::Read);
302 0 : subImStor[k]->dividePSFByWeight();
303 0 : copyBeamSet(*(subImStor[k]->psf()), k);
304 : //subImStor[k]->psf()->flush();
305 0 : }
306 : }
307 0 : Int chanBeg=0; Int chanEnd=0;
308 0 : if(k==0){
309 0 : chanBeg=chanRange_p[0];
310 0 : chanEnd=chanRange_p[1];
311 : }
312 : else{
313 0 : chanBeg=0;
314 0 : chanEnd=subImStor[k]->sumwt()->shape()[3]-1;
315 : }
316 0 : if(subImStor[k]->getType() != "multiterm"){
317 0 : writeBackToFullImage(psfNames_p[k], chanBeg, chanEnd, (subImStor[k]->psf()));
318 0 : writeBackToFullImage(sumwtNames_p[k], chanBeg, chanEnd, (subImStor[k]->sumwt()));
319 0 : if((subImStor[k]->hasSensitivity()) && Table::isWritable(weightNames_p[k])){
320 0 : writeBackToFullImage(weightNames_p[k], chanBeg, chanEnd, (subImStor[k]->weight()));
321 :
322 : }
323 0 : if( doPB_p && Table::isWritable(pbNames_p[k])){
324 0 : writeBackToFullImage(pbNames_p[k], chanBeg, chanEnd, (subImStor[k]->pb()));
325 :
326 : }
327 :
328 : }
329 0 : subImStor[k]->psf()->unlock();
330 :
331 : }
332 0 : }
333 : //cerr << "***Time gridding/ffting " << tim.real() << endl;
334 0 : status_p = True;
335 0 : }
336 0 : catch (AipsError x) {
337 :
338 0 : LogIO os( LogOrigin("SynthesisImagerVi2","CubeMajorCycle",WHERE) );
339 0 : os << LogIO::WARN << "Exception for chan range " << chanRange_p << " --- "<< x.getMesg() << LogIO::POST;
340 0 : cerr << "##################################\n#############################\nException: " << x.getMesg() << endl;
341 :
342 0 : status_p=false;
343 0 : }
344 0 : catch(std::exception& exc) {
345 0 : logger << LogIO::SEVERE << "Exception (std): " << exc.what() << LogIO::POST;
346 0 : returnRec_p=Record();
347 0 : }
348 0 : catch(...){
349 0 : cerr << "###################################\n#######################3##\nUnknown exception " << endl;
350 0 : std::exception_ptr eptr=std::current_exception();
351 : try{
352 0 : if(eptr)
353 0 : std::rethrow_exception(eptr);
354 : }
355 0 : catch(const std::exception& a){
356 0 : LogIO os( LogOrigin("SynthesisImagerVi2","CubeMajorCycle",WHERE) );
357 0 : os << LogIO::WARN << "Unknown Exception for chan range " << chanRange_p << " --- "<< a.what() << LogIO::POST;
358 0 : }
359 :
360 0 : logger << LogIO::SEVERE << "Unknown exception " << LogIO::POST;
361 0 : status_p=False;
362 0 : }
363 0 : }
364 0 : String& CubeMajorCycleAlgorithm::name(){
365 0 : return myName_p;
366 : }
367 :
368 0 : shared_ptr<SIImageStore> CubeMajorCycleAlgorithm::subImageStore(const int imId){
369 : //For some reason multiterm deconvolver is allowed with cubes !
370 0 : String isMTdeconv="";
371 0 : nterms_p.resize(imId+1, True);
372 0 : nterms_p[imId]=-1;
373 0 : if(imId==0 && imSel_p[imId].deconvolver=="mtmfs"){
374 0 : isMTdeconv=".tt0";
375 0 : nterms_p[0]=1;
376 : }
377 :
378 0 : String residname=imSel_p[imId].imageName+".residual"+isMTdeconv;
379 0 : String psfname=imSel_p[imId].imageName+".psf"+isMTdeconv;
380 0 : String sumwgtname=imSel_p[imId].imageName+".sumwt"+isMTdeconv;
381 0 : if(imId > 0 && !Table::isReadable(sumwgtname)){
382 0 : return multiTermImageStore(imId);
383 : }
384 0 : shared_ptr<ImageInterface<Float> >subpsf=nullptr;
385 0 : shared_ptr<ImageInterface<Float> >subresid=nullptr;
386 0 : shared_ptr<ImageInterface<Float> >submodel=nullptr;
387 0 : shared_ptr<ImageInterface<Float> > subweight=nullptr;
388 0 : shared_ptr<ImageInterface<Float> > subpb=nullptr;
389 : //cerr << imId << " sumwt name " << sumwgtname << endl;
390 0 : String workingdir="";
391 0 : if(controlRecord_p.isDefined("workingdirectory"))
392 0 : controlRecord_p.get("workingdirectory", workingdir);
393 : //cerr <<"WORKINGDIR " << workingdir << endl;
394 0 : if(Table::isReadable(workingdir+"/"+sumwgtname)){
395 0 : workingdir=workingdir+"/";
396 0 : residname=workingdir+residname;
397 0 : psfname=workingdir+psfname;
398 0 : sumwgtname=workingdir+sumwgtname;
399 :
400 : }
401 0 : else if(!Table::isReadable(sumwgtname) )
402 0 : throw(AipsError("Programmer error: sumwt disk image is non existant"));
403 : else
404 0 : workingdir="";
405 : ///Use user locking to make sure locks are compliant
406 0 : PagedImage<Float> sumwt(sumwgtname, TableLock::UserLocking);
407 0 : sumwt.lock(FileLocker::Write, 200);
408 : //PagedImage<Float> sumwt(sumwgtname, TableLock::AutoNoReadLocking);
409 0 : if(sumwtNames_p.nelements() <= uInt(imId)){
410 0 : sumwtNames_p.resize(imId+1, True);
411 0 : psfNames_p.resize(imId+1, True);
412 0 : residualNames_p.resize(imId+1, True);
413 0 : sumwtNames_p[imId]=sumwgtname;
414 0 : residualNames_p[imId]=residname;
415 0 : psfNames_p[imId]=psfname;
416 : }
417 0 : Int nchannels=sumwt.shape()[3];
418 :
419 :
420 :
421 : //Should be partitioning for main image only
422 : //chanRange
423 0 : Int chanBeg=0;
424 0 : Int chanEnd=0;
425 0 : if(imId==0){
426 0 : chanBeg=chanRange_p[0];
427 0 : chanEnd=chanRange_p[1];
428 : }
429 : else{
430 0 : chanBeg=0;
431 0 : chanEnd=sumwt.shape()[3]-1;
432 : }
433 0 : sumwt.tempClose();
434 : ////For some small channel ms's retuning trigger a vi2/vb2 bug in nChannels
435 : ///avoid retuning for small images
436 : ////Skipping this here..
437 : //overloaded SynthesisImagerVi2::retune to do this check
438 : //if(nchannels < 30 && imId==0 && ((chanEnd-chanBeg) < 10)){
439 : // retuning_p=False;
440 : //}
441 : //cerr << "chanBeg " << chanBeg << " chanEnd " << chanEnd << " imId " << imId << endl;
442 0 : Vector<String> weightnams(controlRecord_p.asArrayString("weightnames"));
443 0 : Vector<String> pbnams(controlRecord_p.asArrayString("pbnames"));
444 0 : pbNames_p.resize();
445 0 : pbNames_p=pbnams;
446 0 : weightNames_p.resize();
447 0 : weightNames_p=weightnams;
448 :
449 0 : if(imId >= int(weightNames_p.nelements()))
450 0 : throw(AipsError("Number of weight images does not match number of image fields defined"));
451 0 : if(Table::isReadable(workingdir+weightNames_p[imId])){
452 0 : weightNames_p[imId]=workingdir+weightNames_p[imId];
453 : }
454 0 : if(Table::isReadable(workingdir+pbNames_p[imId])){
455 0 : pbNames_p[imId]=workingdir+pbNames_p[imId];
456 : }
457 :
458 0 : if(dopsf_p){
459 : //PagedImage<Float> psf(psfname, TableLock::UserNoReadLocking);
460 : //subpsf.reset(SpectralImageUtil::getChannel(psf, chanBeg, chanEnd, true));
461 0 : getSubImage(subpsf, chanBeg, chanEnd, psfname, False);
462 : //cerr << "PBNAMES " << pbNames_p << endl;
463 0 : if(Table::isReadable(pbNames_p[imId])){
464 :
465 0 : doPB_p=True;
466 0 : getSubImage(subpb, chanBeg, chanEnd, pbNames_p[imId], False);
467 : }
468 : }
469 : else{
470 : //need to loop over all fields somewhere
471 : //PagedImage<Float> resid(residname, TableLock::UserNoReadLocking);
472 : //subresid.reset(SpectralImageUtil::getChannel(resid, chanBeg, chanEnd, true));
473 0 : getSubImage(subresid, chanBeg, chanEnd, residname, False);
474 :
475 0 : if(controlRecord_p.isDefined("modelnames")){
476 0 : Vector<String> modelnames(controlRecord_p.asArrayString("modelnames"));
477 0 : if(imId >= int(modelnames.nelements()))
478 0 : throw(AipsError("Number of model images does not match number of image fields defined"));
479 0 : if(Table::isReadable(workingdir+modelnames[imId])){
480 0 : modelnames[imId]=workingdir+modelnames[imId];
481 : }
482 0 : if(Table::isReadable(modelnames[imId])){
483 :
484 :
485 : ///Pass some extra channels for interpolation while degridding..should match or be less than in SynthesisImager::tuneSelect
486 0 : Int startmodchan=(chanBeg-2) >0 ? chanBeg-2 : 0;
487 0 : Int endmodchan=(chanEnd+2) < nchannels ? chanEnd+2 : nchannels-1 ;
488 : //cerr << "START END mod " << startmodchan << " " << endmodchan << endl;
489 : //Darn has to lock it as writable because overlap in SIMapperCollection code
490 : //wants that...though we are not really modifying it here
491 : //Bool writeisneeded=(imSel_p.nelements()!=1 || startmodel_p[imId].nelements() >0);
492 : //PagedImage<Float> model(modelnames[imId], TableLock::UserNoReadLocking);
493 : //submodel.reset(SpectralImageUtil::getChannel(model, startmodchan, endmodchan, writeisneeded));
494 0 : getSubImage(submodel, startmodchan, endmodchan, modelnames[imId], True);
495 : //
496 0 : if(Table::isReadable(weightNames_p[imId])){
497 0 : divideModelByWeight(imSel_p[imId].imageName, submodel, startmodchan, endmodchan, weightNames_p[imId]);
498 : }
499 : //ImageInterface<Float>* modim=new PagedImage<Float>(modelnames[imId], TableLock::UserNoReadLocking);
500 : //submodel.reset(modim);
501 :
502 : }
503 :
504 0 : }
505 :
506 : }
507 :
508 0 : if(Table::isReadable(weightNames_p[imId])){
509 : //PagedImage<Float> weight(weightnames[imId], TableLock::UserNoReadLocking);
510 : //subweight.reset(SpectralImageUtil::getChannel(weight, chanBeg, chanEnd, true));
511 0 : getSubImage(subweight, chanBeg, chanEnd, weightNames_p[imId], True);
512 : }
513 0 : shared_ptr<ImageInterface<Float> >subsumwt=nullptr;
514 : // subsumwt.reset(SpectralImageUtil::getChannel(sumwt, chanBeg, chanEnd, true));
515 0 : getSubImage(subsumwt, chanBeg, chanEnd, sumwgtname, True);
516 0 : bool useweightimage=(subweight) ? true : false;
517 0 : shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imSel_p[imId].imageName, submodel, subresid, subpsf, subweight, nullptr, nullptr, subsumwt, nullptr, subpb, nullptr, nullptr, useweightimage));
518 0 : if(polRep_p[imId]< 0)
519 0 : throw(AipsError("data polarization type is not defined"));
520 0 : StokesImageUtil::PolRep polrep=(StokesImageUtil::PolRep)polRep_p[imId];
521 0 : subimstor->setDataPolFrame(polrep);
522 0 : if(startmodel_p[imId].nelements() >0){
523 0 : LatticeLocker lock1 (*(subimstor->model()), FileLocker::Write);
524 0 : subimstor->setModelImage(startmodel_p[imId]);
525 0 : }
526 : //cerr << "subimagestor TYPE" << subimstor->getType() << endl;
527 0 : return subimstor;
528 0 : }
529 :
530 :
531 0 : std::shared_ptr<SIImageStore> CubeMajorCycleAlgorithm::multiTermImageStore(const Int imId){
532 0 : uInt nterms=0;
533 0 : String sumwgtname=imSel_p[imId].imageName+".sumwt.tt"+String::toString(nterms);
534 0 : while (Table::isReadable(sumwgtname)){
535 0 : ++nterms;
536 0 : sumwgtname=imSel_p[imId].imageName+".sumwt.tt"+String::toString(nterms);
537 : }
538 0 : if(nterms==0){
539 0 : throw(AipsError("outlier "+String::toString(imId)+" field weight image is not defined"));
540 : }
541 0 : nterms=(nterms+1)/2;
542 0 : nterms_p[imId]=nterms;
543 0 : shared_ptr<SIImageStore> subimstor(new SIImageStoreMultiTerm(imSel_p[imId].imageName, nterms, True));
544 0 : if(polRep_p[imId]< 0)
545 0 : throw(AipsError("data polarization type is not defined"));
546 0 : StokesImageUtil::PolRep polrep=(StokesImageUtil::PolRep)polRep_p[imId];
547 0 : subimstor->setDataPolFrame(polrep);
548 0 : return subimstor;
549 0 : }
550 0 : void CubeMajorCycleAlgorithm::divideModelByWeight(casacore::String imageName, shared_ptr<ImageInterface<Float> >&submodel, const Int startchan, const Int endchan, const String weightname){
551 0 : if(controlRecord_p.isDefined("dividebyweight") && controlRecord_p.asBool("dividebyweight")){
552 0 : if(controlRecord_p.isDefined("normpars")){
553 0 : Record normpars=controlRecord_p.asRecord("normpars");
554 0 : SynthesisNormalizer norm;
555 0 : norm.setupNormalizer(normpars);
556 0 : shared_ptr<ImageInterface<Float> > subweight=nullptr;
557 0 : getSubImage(subweight, startchan, endchan, weightname, True);
558 0 : LatticeLocker lock1(*(subweight), FileLocker::Read);
559 0 : LatticeLocker lock2(*(submodel), FileLocker::Read);
560 0 : std::shared_ptr<SIImageStore> subimstor(new SimpleSIImageStore(imageName, submodel, nullptr, nullptr, subweight, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, nullptr, True));
561 0 : norm.setImageStore(subimstor);
562 0 : norm.divideModelByWeight();
563 :
564 0 : }
565 : }
566 0 : }
567 :
568 :
569 :
570 :
571 0 : void CubeMajorCycleAlgorithm::reset(){
572 :
573 0 : dataSel_p.resize();
574 0 : imSel_p.resize();
575 0 : gridSel_p.resize();
576 0 : ftmRec_p.resize();
577 0 : iftmRec_p.resize();
578 0 : polRep_p.resize();
579 0 : chanRange_p.resize();
580 0 : dopsf_p=False;
581 0 : controlRecord_p=Record();
582 0 : weightParams_p=Record();
583 0 : returnRec_p=Record();
584 0 : startmodel_p.resize();
585 0 : status_p=False;
586 0 : residualNames_p.resize();
587 0 : psfNames_p.resize();
588 0 : sumwtNames_p.resize();
589 0 : weightNames_p.resize();
590 0 : pbNames_p.resize();
591 0 : movingSource_p="";
592 0 : retuning_p=True;
593 0 : doPB_p=False;
594 0 : nterms_p.resize();
595 :
596 :
597 :
598 0 : }
599 :
600 0 : void CubeMajorCycleAlgorithm::getSubImage(std::shared_ptr<ImageInterface<Float> >& subimptr, const Int chanBeg, const Int chanEnd, const String imagename, const Bool copy){
601 0 : LogIO os( LogOrigin("CubeMajorCycle","getSubImage",WHERE) );
602 0 : CountedPtr<PagedImage<Float> > im=nullptr;
603 : {///work around
604 0 : uInt exceptCounter=0;
605 : while(true){
606 : try{
607 0 : im = new PagedImage<Float> (imagename, TableLock::UserLocking);
608 : //PagedImage<Float> im(imagename, TableLock::AutoNoReadLocking);
609 0 : im->lock(FileLocker::Read, 1000);
610 0 : SubImage<Float> *tmpptr=nullptr;
611 0 : tmpptr=SpectralImageUtil::getChannel(*im, chanBeg, chanEnd, false);
612 0 : subimptr.reset(new TempImage<Float>(TiledShape(tmpptr->shape(), tmpptr->niceCursorShape()), tmpptr->coordinates()));
613 :
614 0 : if(copy)
615 0 : subimptr->copyData(*tmpptr);
616 : //subimptr->flush();
617 0 : im->unlock();
618 0 : delete tmpptr;
619 0 : break; ///get out of while loop
620 : }
621 0 : catch(AipsError &x){
622 0 : im=nullptr;
623 0 : String mes=x.getMesg();
624 0 : if(mes.contains("FilebufIO::readBlock") ){
625 0 : std::this_thread::sleep_for(std::chrono::milliseconds(50));
626 0 : os << LogIO::WARN << "#####CATCHING a sleep because "<< mes<< LogIO::POST;
627 : }
628 : else
629 0 : throw(AipsError("Error in readimage "+imagename+" : "+mes));
630 :
631 0 : if(exceptCounter > 100){
632 0 : throw(AipsError("Error in readimage got 100 of this exeception: "+mes));
633 :
634 : }
635 :
636 0 : }
637 0 : ++exceptCounter;
638 0 : }
639 : }
640 :
641 0 : }
642 :
643 0 : void CubeMajorCycleAlgorithm::writeBackToFullImage(const String imagename, const Int chanBeg, const Int chanEnd, std::shared_ptr<ImageInterface<Float> > subimptr){
644 0 : LogIO os( LogOrigin("CubeMajorCycle","writebacktofullimage",WHERE) );
645 0 : CountedPtr<PagedImage<Float> > im=nullptr;
646 : {///work around
647 0 : uInt exceptCounter=0;
648 : while(true){
649 : try{
650 0 : im=new PagedImage<Float> (imagename, TableLock::UserLocking);
651 :
652 :
653 : //PagedImage<Float> im(imagename, TableLock::AutoLocking);
654 0 : SubImage<Float> *tmpptr=nullptr;
655 : {
656 0 : tmpptr=SpectralImageUtil::getChannel(*im, chanBeg, chanEnd, true);
657 0 : LatticeLocker lock1 (*(tmpptr), FileLocker::Write);
658 0 : tmpptr->set(0.0);
659 0 : tmpptr->copyData(*(subimptr));
660 0 : }
661 0 : im->flush();
662 0 : im->unlock();
663 0 : delete tmpptr;
664 0 : break; //get out of while loop
665 : }
666 0 : catch(AipsError &x){
667 :
668 0 : String mes=x.getMesg();
669 0 : if(mes.contains("FilebufIO::readBlock") ){
670 0 : std::this_thread::sleep_for(std::chrono::milliseconds(50));
671 0 : os << LogIO::WARN << "#####CATCHING a sleep because "<< mes<< LogIO::POST;
672 : }
673 : else
674 0 : throw(AipsError("Error in writing back image "+imagename+" : "+mes));
675 :
676 0 : if(exceptCounter > 100){
677 0 : throw(AipsError("Error in writeimage got 100 of this exeception: "+mes));
678 :
679 : }
680 :
681 0 : }
682 0 : ++exceptCounter;
683 0 : }
684 : }//End of work around for table disappearing bug
685 :
686 :
687 0 : }
688 0 : void CubeMajorCycleAlgorithm::copyBeamSet(ImageInterface<Float>& subpsf, const Int imageid){
689 : //For now lets do for the main image only
690 0 : if(imageid != 0)
691 0 : return;
692 0 : ImageBeamSet iibeamset=subpsf.imageInfo().getBeamSet();
693 0 : uInt nchan=chanRange_p[1]-chanRange_p[0]+1;
694 0 : if(nchan ==iibeamset.nchan()){
695 0 : Matrix<GaussianBeam> gbs=iibeamset.getBeams();
696 0 : Cube<Float> beams(nchan, iibeamset.nstokes(),3);
697 :
698 0 : for (uInt k=0; k < nchan; ++k){
699 0 : for (uInt j=0; j < iibeamset.nstokes(); ++j){
700 0 : beams(k,j, 0)=gbs(k, 0).getMajor("arcsec");
701 0 : beams(k,j, 1)=gbs(k, 0).getMinor("arcsec");
702 0 : beams(k,j, 2)=gbs(k, 0).getPA("deg", True);
703 : }
704 : }
705 0 : returnRec_p.define("imageid", imageid);
706 0 : returnRec_p.define("chanrange", chanRange_p);
707 0 : returnRec_p.define("beams", beams);
708 :
709 :
710 0 : }
711 :
712 :
713 :
714 0 : }
715 :
716 :
717 :
718 : } //# NAMESPACE CASA - END
|
