Line data Source code
1 : //# SIMapperCollection.cc: Implementation of Imager.h
2 : //# Copyright (C) 1997-2008
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 : #include <casacore/casa/Logging.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 : #include <casacore/casa/Logging/LogMessage.h>
41 :
42 : #include <casacore/casa/OS/DirectoryIterator.h>
43 : #include <casacore/casa/OS/File.h>
44 : #include <casacore/casa/OS/Path.h>
45 :
46 : #include <casacore/casa/OS/HostInfo.h>
47 :
48 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
49 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
50 :
51 : #include <synthesis/ImagerObjects/SIMapperCollection.h>
52 : #include <synthesis/ImagerObjects/SynthesisUtilMethods.h>
53 :
54 : #include <synthesis/TransformMachines/VisModelData.h>
55 : #include <casacore/images/Regions/WCBox.h>
56 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
57 : #include <msvis/MSVis/VisBufferComponents2.h>
58 :
59 : #include <sys/types.h>
60 : #include <unistd.h>
61 : using namespace std;
62 :
63 : using namespace casacore;
64 : namespace casa { //# NAMESPACE CASA - BEGIN
65 :
66 : //////////////////////////////////////////////////////////////////////////////////////////////////////
67 : //////////////////////////////////////////////////////////////////////////////////////////////////////
68 :
69 1703 : SIMapperCollection::SIMapperCollection()
70 : {
71 3406 : LogIO os( LogOrigin("SIMapperCollection","Construct a mapperCollection",WHERE) );
72 :
73 1703 : itsMappers.resize(0);
74 1703 : oldMsId_p=-1;
75 1703 : itsIsNonZeroModel=false;
76 :
77 1703 : }
78 :
79 : //////////////////////////////////////////////////////////////////////////////////////////////////////
80 : //////////////////////////////////////////////////////////////////////////////////////////////////////
81 :
82 1703 : SIMapperCollection::~SIMapperCollection()
83 : {
84 1703 : }
85 :
86 : //////////////////////////////////////////////////////////////////////////////////////////////////////
87 : //////////////////////////////////////////////////////////////////////////////////////////////////////
88 :
89 3544 : Bool SIMapperCollection::releaseImageLocks()
90 : {
91 3544 : Bool validflag=true;
92 7282 : for(Int mapperid=0;mapperid<nMappers();mapperid++ )
93 : {
94 3738 : validflag &= itsMappers[mapperid]->releaseImageLocks();
95 : }
96 3544 : return validflag;
97 : }
98 :
99 : //////////////////////////////////////////////////////////////////////////////////////////////////////
100 : //////////////////////////////////////////////////////////////////////////////////////////////////////
101 :
102 2642 : std::vector<String> SIMapperCollection::cleanupTempFiles(const String& mess)
103 : {
104 2642 : std::vector<String> outstr;
105 5556 : auto appvectors = [](std::vector<String>& a, const std::vector<String> b) { a.insert(std::end(a), std::begin(b), std::end(b));};
106 5420 : for(Int mapperid=0;mapperid<nMappers();mapperid++ )
107 : {
108 2778 : if((itsMappers[mapperid]->getFTM2())){
109 2778 : appvectors( outstr, (itsMappers[mapperid]->getFTM2(True)->cleanupTempFiles(mess)).tovector());
110 2778 : appvectors(outstr, (itsMappers[mapperid]->getFTM2(False)->cleanupTempFiles(mess)).tovector());
111 : }
112 : }
113 5284 : return outstr;
114 0 : }
115 :
116 : //////////////////////////////////////////////////////////////////////////////////////////////////////
117 : //////////////////////////////////////////////////////////////////////////////////////////////////////
118 :
119 : /*
120 : // Allocate Memory and open images.
121 : void SIMapperCollection::addMapper( String mappertype,
122 : CountedPtr<SIImageStore> imagestore,
123 : CountedPtr<FTMachine> ftmachine,
124 : CountedPtr<FTMachine> iftmachine)
125 : {
126 :
127 : LogIO os( LogOrigin("SIMapperCollection","addMapper",WHERE) );
128 :
129 : CountedPtr<SIMapper> localMapper=NULL;
130 : Int nMappers = itsMappers.nelements();
131 : // Check 'mappertype' for valid types....
132 : if( mappertype == "default" )
133 : {
134 : localMapper = new SIMapperSingle( imagestore, ftmachine, iftmachine, nMappers );
135 : }
136 :
137 : else if( mappertype == "multiterm" )
138 : {
139 : localMapper = new SIMapperMultiTerm( imagestore, ftmachine, iftmachine,nMappers, ntaylorterms );
140 : }
141 :
142 : else
143 : {
144 : throw ( AipsError("Internal Error : Unrecognized Mapper Type in MapperCollection.addMapper") );
145 : }
146 :
147 : // If all is well, add to the list.
148 : itsMappers.resize(nMappers+1, true);
149 : itsMappers[nMappers] = localMapper;
150 :
151 : }
152 : */
153 :
154 : //////////////////////////////////////////////////////////////////////////////////////////////////////
155 :
156 :
157 : //////////////////////////////////////////////////////////////////////////////////////////////////////
158 1760 : void SIMapperCollection::addMapper( CountedPtr<SIMapper> map){
159 1760 : Int nMappers = itsMappers.nelements();
160 1760 : itsMappers.resize(nMappers+1, true);
161 1760 : itsMappers[nMappers]=map;
162 1760 : }
163 :
164 25695 : Int SIMapperCollection::nMappers()
165 : {
166 25695 : return itsMappers.nelements();
167 : }
168 :
169 : //////////////////////////////////////////////////////////////////////////////////////////////////////
170 : //////////////////////////////////////////////////////////////////////////////////////////////////////
171 :
172 0 : Vector<String> SIMapperCollection::getImageNames()
173 : {
174 0 : Vector<String> names( nMappers() );
175 :
176 0 : for(Int mapperid=0;mapperid<nMappers();mapperid++ )
177 : {
178 0 : names[mapperid] = itsMappers[mapperid]->getImageName();
179 : }
180 :
181 0 : return names;
182 0 : }
183 :
184 : //////////////////////////////////////////////////////////////////////////////////////////////////////
185 : /////////////////// Start VB dependent code /////////////////////////////////////////////
186 2259 : void SIMapperCollection::initializeGrid(vi::VisBuffer2& vb, Bool dopsf, const Int mapperid)
187 : {
188 2259 : if(mapperid<0)
189 : {
190 4644 : for (uInt k=0; k < itsMappers.nelements(); ++k)
191 : {
192 2386 : (itsMappers[k])->initializeGrid(vb,dopsf,true);
193 : }
194 : }
195 : else
196 : {
197 0 : if (mapperid > (Int)itsMappers.nelements())
198 0 : throw ( AipsError("Internal Error : SIMapperCollection::initializeGrid(): mapperid out of range") );
199 0 : else itsMappers[mapperid]->initializeGrid(vb, dopsf, true);
200 : }
201 2258 : }
202 :
203 2259 : void SIMapperCollection::initializeGrid(vi::VisibilityIterator2& vi, Bool dopsf, const Int mapperid)
204 : {
205 :
206 2259 : vi::VisBuffer2 *vb=vi.getVisBuffer();
207 2259 : initializeGrid(*vb, dopsf, mapperid);
208 2258 : if (mapperid < 0) {
209 4643 : for (uInt k=0; k < itsMappers.nelements(); ++k) {
210 2385 : ((itsMappers[k])->getFTM2())->initBriggsWeightor(vi);
211 : }
212 : }
213 0 : else if (mapperid > (Int)itsMappers.nelements()) {
214 0 : throw (
215 : AipsError("Internal Error : "
216 : "SIMapperCollection::initializeGrid(): mapperid out of range")
217 0 : );
218 : }
219 : else {
220 0 : (itsMappers[mapperid]->getFTM2())->initBriggsWeightor(vi);
221 : }
222 2258 : }
223 :
224 : ////////////////////////////////////////////////////////////////////////////////////
225 : /////////////////////////////////OLD vi/vb //////////////////////////////////////////////
226 0 : void SIMapperCollection::initializeGrid(VisBuffer& vb, Bool dopsf, const Int mapperid)
227 : {
228 0 : if (mapperid < 0) {
229 0 : for (uInt k=0; k < itsMappers.nelements(); ++k) {
230 0 : (itsMappers[k])->initializeGrid(vb,dopsf,true);
231 : }
232 : }
233 0 : else if (mapperid > (Int)itsMappers.nelements()) {
234 0 : throw (
235 : AipsError("Internal Error : "
236 : "SIMapperCollection::initializeGrid(): mapperid out of range")
237 0 : );
238 : }
239 : else {
240 0 : itsMappers[mapperid]->initializeGrid(vb, dopsf, true);
241 : }
242 0 : }
243 :
244 : //////////////////////////////////////////////////////////////////////////////
245 :
246 160 : void SIMapperCollection::handleNewMs(const MeasurementSet &ms, const Int mapperid)
247 : {
248 160 : if (mapperid < 0) {
249 320 : for (uInt k=0; k < itsMappers.nelements(); ++k) {
250 160 : (itsMappers[k])->handleNewMs(ms);
251 : }
252 : }
253 0 : else if (mapperid > (Int)itsMappers.nelements()) {
254 0 : throw (
255 : AipsError("Internal Error : "
256 : "SIMapperCollection::handleNewMs(): mapperid out of range")
257 0 : );
258 : }
259 : else {
260 0 : itsMappers[mapperid]->handleNewMs(ms);
261 : }
262 160 : }
263 :
264 675482 : void SIMapperCollection::grid(vi::VisBuffer2& vb, Bool dopsf, refim::FTMachine::Type col,
265 : Int mapperid)
266 : {
267 675482 : if( itsIsNonZeroModel == true ) // Try to subtract model visibilities only if a model exists.
268 : {
269 151500 : if(col==refim::FTMachine::CORRECTED){
270 : //Dang i thought the new vb will return Data or FloatData if correctedData was
271 : //not there
272 146388 : if(!vb.existsColumn(VisBufferComponent2::VisibilityCorrected)){
273 62468 : col=refim::FTMachine::OBSERVED;
274 : // cerr << "Max of visCube" << max(vb.visCube()) << " model " << max(vb.modelVisCube())<< endl;
275 62468 : vb.setVisCube(vb.visCube()-vb.visCubeModel());
276 : }
277 : else{
278 83920 : vb.setVisCubeCorrected(vb.visCubeCorrected()-vb.visCubeModel());
279 : }
280 : }
281 5112 : else if (col==refim::FTMachine::OBSERVED) {
282 5112 : vb.setVisCube(vb.visCube()-vb.visCubeModel());
283 : }
284 : }// if non zero model
285 :
286 675482 : if(col==refim::FTMachine::CORRECTED && !vb.existsColumn(VisBufferComponent2::VisibilityCorrected)) {
287 : //cout << "Corrected column isn't there, using data instead" << endl;
288 300971 : col=refim::FTMachine::OBSERVED;
289 : }
290 :
291 675482 : if (mapperid < 0)
292 : {
293 : //cout << "Using column : " << col << endl;
294 1364187 : for (uInt k=0; k < itsMappers.nelements(); ++k)
295 : {
296 689414 : (itsMappers[k])->grid(vb, dopsf, col);
297 : }
298 : }
299 : else
300 : {
301 708 : if (mapperid > (Int)itsMappers.nelements())
302 0 : throw ( AipsError("Internal Error : SIMapperCollection::grid(): mapperid out of range") );
303 708 : else itsMappers[mapperid]->grid(vb, dopsf, col);
304 : }
305 675481 : }
306 :
307 : ///////////////////
308 : ///////////////////////////////////////OLD VI/VB ///////////////////////////////////
309 0 : void SIMapperCollection::grid(VisBuffer& vb, Bool dopsf, FTMachine::Type col,
310 : Int mapperid)
311 : {
312 0 : if( itsIsNonZeroModel == true ) // Try to subtract model visibilities only if a model exists.
313 : {
314 0 : if(col==FTMachine::CORRECTED){
315 0 : if(vb.msColumns().correctedData().isNull()){
316 0 : col=FTMachine::OBSERVED;
317 : // cerr << "Max of visCube" << max(vb.visCube()) << " model " << max(vb.modelVisCube())<< endl;
318 0 : vb.visCube()-=vb.modelVisCube();
319 : }
320 : else{
321 0 : vb.correctedVisCube()-=vb.modelVisCube();
322 : }
323 : }
324 0 : else if (col==FTMachine::OBSERVED) {
325 0 : vb.visCube()-=vb.modelVisCube();
326 : }
327 : }// if non zero model
328 :
329 0 : if(col==FTMachine::CORRECTED && vb.msColumns().correctedData().isNull())
330 0 : { col=FTMachine::OBSERVED;}
331 :
332 0 : if (mapperid < 0)
333 : {
334 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
335 : {
336 0 : (itsMappers[k])->grid(vb, dopsf, col);
337 :
338 : }
339 : }
340 : else
341 : {
342 0 : if (mapperid > (Int)itsMappers.nelements())
343 0 : throw ( AipsError("Internal Error : SIMapperCollection::grid(): mapperid out of range") );
344 0 : else itsMappers[mapperid]->grid(vb, dopsf, col);
345 : }
346 0 : }
347 : ///////////////////////////////
348 : ////////////////////////////////
349 2250 : void SIMapperCollection::finalizeGrid(vi::VisBuffer2& vb, Bool dopsf,const Int mapperid)
350 : {
351 2250 : if(mapperid<0)
352 : {
353 4627 : for (uInt k=0; k < itsMappers.nelements(); ++k)
354 : {
355 2377 : (itsMappers[k])->finalizeGrid(vb, dopsf);
356 : }
357 : }
358 : else
359 : {
360 0 : if (mapperid > (Int)itsMappers.nelements())
361 0 : throw ( AipsError("Internal Error : SIMapperCollection::finalizeGrid(): mapperid out of range") );
362 0 : else itsMappers[mapperid]->finalizeGrid(vb, dopsf);
363 : }
364 2250 : }
365 :
366 :
367 : ////////////////////////////////////////////////////////////////////////////////////////////////////////
368 : /////////////////////////////////////////////OLD VI/VB////////////////////////////////////////////////
369 0 : void SIMapperCollection::finalizeGrid(VisBuffer& vb, Bool dopsf,const Int mapperid)
370 : {
371 0 : if(mapperid<0)
372 : {
373 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
374 : {
375 0 : (itsMappers[k])->finalizeGrid(vb, dopsf);
376 : }
377 : }
378 : else
379 : {
380 0 : if (mapperid > (Int)itsMappers.nelements())
381 0 : throw ( AipsError("Internal Error : SIMapperCollection::finalizeGrid(): mapperid out of range") );
382 0 : else itsMappers[mapperid]->finalizeGrid(vb, dopsf);
383 : }
384 0 : }
385 :
386 :
387 : //////////////////////////////////////////////////////////////////////////////////////////////////////
388 1454 : void SIMapperCollection::initializeDegrid(vi::VisBuffer2& vb, const Int mapperid)
389 : {
390 :
391 1454 : itsIsNonZeroModel = anyNonZeroModels();
392 :
393 1454 : if( itsIsNonZeroModel == true )
394 : {
395 : // vb.setModelVisCube( Complex(0.0,0.0) );
396 :
397 806 : if(mapperid<0)
398 : {
399 1635 : for (uInt k=0; k < itsMappers.nelements(); ++k)
400 : {
401 829 : (itsMappers[k])->initializeDegrid(vb);
402 : }
403 : }
404 : else
405 : {
406 0 : if (mapperid > (Int)itsMappers.nelements())
407 0 : throw ( AipsError("Internal Error : SIMapperCollection::initializeDegrid(): mapperid out of range") );
408 0 : else itsMappers[mapperid]->initializeDegrid(vb);
409 : }
410 :
411 : }// if non zero model
412 1454 : }
413 :
414 :
415 :
416 : ///////////////////////////////////////OLD VI/VB ///////////////////////////////////////////////////
417 0 : void SIMapperCollection::initializeDegrid(VisBuffer& vb, const Int mapperid)
418 : {
419 :
420 0 : itsIsNonZeroModel = anyNonZeroModels();
421 :
422 0 : if( itsIsNonZeroModel == true )
423 : {
424 : // vb.setModelVisCube( Complex(0.0,0.0) );
425 :
426 0 : if(mapperid<0)
427 : {
428 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
429 : {
430 0 : (itsMappers[k])->initializeDegrid(vb);
431 : }
432 : }
433 : else
434 : {
435 0 : if (mapperid > (Int)itsMappers.nelements())
436 0 : throw ( AipsError("Internal Error : SIMapperCollection::initializeDegrid(): mapperid out of range") );
437 0 : else itsMappers[mapperid]->initializeDegrid(vb);
438 : }
439 :
440 : }// if non zero model
441 0 : }
442 :
443 : ////////////////////////////////////////////////////////////////////////////////////////////////
444 321929 : void SIMapperCollection::degrid(vi::VisBuffer2& vb, Bool saveVirtualMod, const Int mapperid)
445 : {
446 321929 : if( itsIsNonZeroModel == true )
447 : {
448 160228 : if(mapperid<0)
449 : {
450 324596 : for (uInt k=0; k < itsMappers.nelements(); ++k)
451 : {
452 164368 : (itsMappers[k])->degrid(vb);
453 : }
454 : }
455 : else
456 : {
457 0 : if (mapperid > (Int)itsMappers.nelements())
458 0 : throw ( AipsError("Internal Error : SIMapperCollection::degrid(): mapperid out of range") );
459 0 : else itsMappers[mapperid]->degrid(vb);
460 : }
461 :
462 160228 : if(saveVirtualMod){
463 8160 : saveVirtualModel(vb);
464 : }
465 : }// if non zero model
466 321929 : }
467 :
468 613 : void SIMapperCollection::addPB(vi::VisBuffer2& vb, PBMath& pbMath, const MDirection& altDir, const Bool useAltDir)
469 : {
470 1278 : for (uInt k=0; k < itsMappers.nelements(); ++k)
471 : {
472 665 : (itsMappers[k])->addPB(vb,pbMath, altDir, useAltDir);
473 :
474 : }
475 613 : }
476 :
477 :
478 : /////////////////////////////////////OLD VI/VB ////////////////////////////////////////////////////
479 0 : void SIMapperCollection::degrid(VisBuffer& vb, Bool saveVirtualMod, const Int mapperid)
480 : {
481 0 : if( itsIsNonZeroModel == true )
482 : {
483 0 : if(mapperid<0)
484 : {
485 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
486 : {
487 0 : (itsMappers[k])->degrid(vb);
488 : }
489 : }
490 : else
491 : {
492 0 : if (mapperid > (Int)itsMappers.nelements())
493 0 : throw ( AipsError("Internal Error : SIMapperCollection::degrid(): mapperid out of range") );
494 0 : else itsMappers[mapperid]->degrid(vb);
495 : }
496 :
497 0 : if(saveVirtualMod){
498 0 : saveVirtualModel(vb);
499 : }
500 : }// if non zero model
501 0 : }
502 :
503 :
504 : /////
505 : /////////////
506 0 : void SIMapperCollection::saveVirtualModel(VisBuffer& vb){
507 :
508 :
509 :
510 0 : if(vb.msId() != oldMsId_p){
511 0 : oldMsId_p=vb.msId();
512 : /*Block< Vector<Int> > blockNGroup;
513 : Block< Vector<Int> > blockStart;
514 : Block< Vector<Int> > blockWidth;
515 : Block< Vector<Int> > blockIncr;
516 : Block< Vector<Int> > blockSpw;
517 : vb.getChannelSelection(blockNGroup, blockStart, blockWidth, blockIncr, blockSpw);
518 : Vector<Int> fields = vb.msColumns().fieldId().getColumn();
519 : const Int option = Sort::HeapSort | Sort::NoDuplicates;
520 : const Sort::Order order = Sort::Ascending;
521 : Int nfields = GenSort<Int>::sort (fields, order, option);
522 :
523 : // Make sure we have the right size
524 :
525 : fields.resize(nfields, true);
526 : */
527 : //Int msid = vb.msId();
528 0 : ROVisibilityIterator *viloc=vb.getVisibilityIterator();
529 0 : for (uInt k=0; k < itsMappers.nelements(); ++k){
530 0 : Record rec;
531 0 : String modImage=viloc->ms().getPartNames()[0];
532 0 : if(!((viloc->ms()).source().isNull()))
533 0 : modImage=(viloc->ms()).source().tableName();
534 0 : modImage=File::newUniqueName(modImage, "FT_MODEL").absoluteName();
535 0 : Bool iscomp=itsMappers[k]->getCLRecord(rec);
536 0 : if(iscomp || itsMappers[k]->getFTMRecord(rec, modImage)){
537 :
538 : //VisModelData::putModel(vb.getVisibilityIterator()->ms(), rec, fields, blockSpw[msid], blockStart[msid],
539 : // blockWidth[msid], blockIncr[msid],
540 : // iscomp, true);
541 0 : VisibilityIterator * elvi=(dynamic_cast<VisibilityIterator* >(vb.getVisibilityIterator()));
542 0 : if(elvi)
543 0 : elvi->putModel(rec, iscomp, true);
544 : // VisModelData::listModel(vb.getVisibilityIterator()->ms());
545 : }
546 :
547 0 : }
548 :
549 :
550 :
551 : }
552 :
553 :
554 :
555 0 : }
556 :
557 8160 : void SIMapperCollection::saveVirtualModel(vi::VisBuffer2& vb){
558 :
559 :
560 :
561 8160 : if(vb.msId() != oldMsId_p){
562 17 : oldMsId_p=vb.msId();
563 :
564 34 : for (uInt k=0; k < itsMappers.nelements(); ++k){
565 17 : Record rec;
566 17 : String modImage=vb.ms().getPartNames()[0];
567 17 : if(!((vb.ms()).source().isNull()))
568 17 : modImage=(vb.ms()).source().tableName();
569 17 : modImage=File::newUniqueName(modImage, "FT_MODEL").absoluteName();
570 17 : Bool iscomp=itsMappers[k]->getCLRecord(rec);
571 17 : if(iscomp || itsMappers[k]->getFTMRecord(rec, modImage)){
572 :
573 : ////Darn not implemented
574 : //static_cast<VisibilityIteratorImpl2 *>(viloc->getImpl())->writeModel(rec, //iscomp, true);
575 :
576 17 : if(!iscomp && Table::isReadable(modImage)){
577 : //make sure complex image is of compliant size/shape
578 15 : (itsMappers[k]->imageStore())->intersectComplexImage(modImage);
579 :
580 : }
581 17 : VisibilityIterator2* vi=const_cast<VisibilityIterator2*>(vb.getVi());
582 17 : const_cast<MeasurementSet& >(vi->ms()).lock();
583 : /////TESTOO
584 : //Int CPUID;
585 : //MPI_Comm_rank(MPI_COMM_WORLD, &CPUID);
586 : //cerr << CPUID << " writing " << modImage << endl;
587 : /////////////////
588 17 : vi->writeModel(rec, iscomp, true);
589 17 : const_cast<MeasurementSet& >(vi->ms()).unlock();
590 : // VisModelData::listModel(vb.getVisibilityIterator()->ms());
591 : }
592 :
593 17 : }
594 :
595 :
596 :
597 : }
598 :
599 :
600 :
601 8160 : }
602 :
603 :
604 : /////////////////////////////////////////////////////////
605 1454 : void SIMapperCollection::finalizeDegrid(vi::VisBuffer2& /*vb*/, const Int mapperid)
606 : {
607 1454 : if( itsIsNonZeroModel == true )
608 : {
609 806 : if(mapperid<0)
610 : {
611 1635 : for (uInt k=0; k < itsMappers.nelements(); ++k)
612 : {
613 829 : (itsMappers[k])->finalizeDegrid();
614 :
615 : }
616 : }
617 : else
618 : {
619 0 : if (mapperid > (Int)itsMappers.nelements())
620 0 : throw ( AipsError("Internal Error : SIMapperCollection::finalizeDegrid(): mapperid out of range") );
621 0 : else itsMappers[mapperid]->finalizeDegrid();
622 : }
623 : }// if non zero model
624 1454 : }
625 :
626 : ////////////////////////////////////////////////////////
627 0 : void SIMapperCollection::finalizeDegrid(VisBuffer& /*vb*/, const Int mapperid)
628 : {
629 0 : if( itsIsNonZeroModel == true )
630 : {
631 0 : if(mapperid<0)
632 : {
633 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
634 : {
635 0 : (itsMappers[k])->finalizeDegrid();
636 :
637 : }
638 : }
639 : else
640 : {
641 0 : if (mapperid > (Int)itsMappers.nelements())
642 0 : throw ( AipsError("Internal Error : SIMapperCollection::finalizeDegrid(): mapperid out of range") );
643 0 : else itsMappers[mapperid]->finalizeDegrid();
644 : }
645 : }// if non zero model
646 0 : }
647 :
648 :
649 579 : void SIMapperCollection::initPB()
650 : {
651 1210 : for (uInt k=0; k < itsMappers.nelements(); ++k)
652 : {
653 631 : (itsMappers[k])->initPB();
654 : }
655 579 : }
656 :
657 0 : void SIMapperCollection::addPB(VisBuffer& vb, PBMath& pbMath)
658 : {
659 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
660 : {
661 0 : (itsMappers[k])->addPB(vb,pbMath);
662 :
663 : }
664 0 : }
665 :
666 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
667 : //////////// End of VB dependent code.
668 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
669 :
670 :
671 :
672 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
673 : /////////////////////
674 18045 : CountedPtr<SIImageStore> SIMapperCollection::imageStore(Int id)
675 : {
676 18045 : if(uInt(id) < itsMappers.nelements())
677 : {
678 18045 : return (itsMappers[id])->imageStore();
679 : }
680 0 : return make_shared<SIImageStore>( );
681 : }
682 : /////////////////////////////////////////////////////////////////////////////////////////////
683 : /////////////////////////////////////////////////////////////////////////////////////////////
684 0 : Record SIMapperCollection::getFTMRecord(Int mapperid)
685 : {
686 0 : AlwaysAssert( mapperid >=0 && mapperid < nMappers() , AipsError );
687 : //return itsMappers[mapperid]->getFTMRecord();
688 0 : Record rec;
689 0 : return rec;
690 : }
691 :
692 : //////////////////////////////////////////////////////////////////////////////////////////////////////
693 : //////////////////////////////////////////////////////////////////////////////////////////////////////
694 :
695 4268 : void SIMapperCollection::checkOverlappingModels(String action)
696 : {
697 : // If nothing that overlaps, don't check.
698 4268 : if(nMappers()==1) return;
699 :
700 104 : Int nmodels = nMappers();
701 :
702 : // If there is no model image (i.e. first major cycle with no starting model), don't check.
703 104 : Bool hasmodel=true;
704 358 : for (Int model=0;model<(nmodels-1); ++model)
705 254 : { hasmodel = hasmodel && ((itsMappers[model])->imageStore())->hasModel(); }
706 104 : if( hasmodel==false ) {
707 : //cout << "No model images to check overlap for." << endl;
708 76 : return;
709 : }
710 :
711 : // Internal check
712 28 : AlwaysAssert( action=="blank" || action=="restore" , AipsError );
713 :
714 74 : for (Int model=0;model<(nmodels-1); ++model)
715 : {
716 : // Connect to one image for aux info.
717 92 : LatticeLocker lockmodel (*(((itsMappers[model])->imageStore())->model()), FileLocker::Write);
718 92 : SubImage<Float> modelimage( *(((itsMappers[model])->imageStore())->model()), true );
719 :
720 46 : uInt nTaylor0 = ((itsMappers[model])->imageStore())->getNTaylorTerms();
721 :
722 46 : CoordinateSystem cs0=modelimage.coordinates();
723 46 : IPosition iblc0(modelimage.shape().nelements(),0);
724 46 : IPosition itrc0(modelimage.shape());
725 46 : itrc0=itrc0-Int(1);
726 46 : LCBox lbox0(iblc0, itrc0, modelimage.shape());
727 46 : ImageRegion imagreg0(WCBox(lbox0, cs0));
728 :
729 118 : for (Int nextmodel=model+1; nextmodel < nmodels; ++nextmodel)
730 : {
731 144 : LatticeLocker nextlockmodel (*(((itsMappers[nextmodel])->imageStore())->model()), FileLocker::Write);
732 144 : SubImage<Float> nextmodelimage( *(((itsMappers[nextmodel])->imageStore())->model()) , true);
733 :
734 72 : uInt nTaylor1 = ((itsMappers[nextmodel])->imageStore())->getNTaylorTerms();
735 :
736 72 : CoordinateSystem cs=nextmodelimage.coordinates();
737 72 : IPosition iblc(nextmodelimage.shape().nelements(),0);
738 72 : IPosition itrc(nextmodelimage.shape());
739 72 : itrc=itrc-Int(1);
740 72 : LCBox lbox(iblc, itrc, nextmodelimage.shape());
741 72 : ImageRegion imagreg(WCBox(lbox, cs));
742 :
743 : try{
744 :
745 72 : if( action.matches("blank") )
746 : {
747 :
748 : //cerr << "blank MODEL image shape " << modelimage.shape() << " " << nextmodelimage.shape() << endl;
749 :
750 69 : LatticeRegion latReg=imagreg.toLatticeRegion(modelimage.coordinates(), modelimage.shape());
751 :
752 7 : for(uInt taylor=0;taylor<min(nTaylor0,nTaylor1);taylor++)
753 : { // loop for taylor term
754 8 : SubImage<Float> modelim( *(((itsMappers[model])->imageStore())->model(taylor)), true );
755 4 : SubImage<Float> partToMask(modelim, imagreg, true);
756 4 : LatticeLocker lock1 (partToMask, FileLocker::Write);
757 4 : ArrayLattice<Bool> pixmask(latReg.get());
758 8 : LatticeExpr<Float> myexpr(iif(pixmask, 0.0, partToMask) );
759 4 : partToMask.copyData(myexpr);
760 4 : }
761 :
762 :
763 3 : }
764 : else // "restore"
765 : {
766 : //cerr << "rsetore MODEL image shape " << modelimage.shape() << " " << nextmodelimage.shape() << endl;
767 40 : LatticeRegion latReg0=imagreg0.toLatticeRegion(nextmodelimage.coordinates(), nextmodelimage.shape());
768 61 : LatticeRegion latReg=imagreg.toLatticeRegion(modelimage.coordinates(), modelimage.shape());
769 3 : ArrayLattice<Bool> pixmask(latReg.get());
770 :
771 :
772 7 : for(uInt taylor=0;taylor<min(nTaylor0,nTaylor1);taylor++)
773 : {// loop for taylor term
774 8 : SubImage<Float> modelim( *(((itsMappers[model])->imageStore())->model(taylor)), true );
775 8 : SubImage<Float> nextmodelim( *(((itsMappers[nextmodel])->imageStore())->model(taylor)), true );
776 :
777 4 : SubImage<Float> partToMerge(nextmodelim, imagreg0, true);
778 4 : SubImage<Float> partToUnmask(modelim, imagreg, true);
779 4 : LatticeLocker lock1 (partToUnmask, FileLocker::Write);
780 4 : LatticeLocker lock2 (partToMerge, FileLocker::Write);
781 8 : LatticeExpr<Float> myexpr0(iif(pixmask,partToMerge,partToUnmask));
782 4 : partToUnmask.copyData(myexpr0);
783 4 : }
784 :
785 32 : }
786 : }
787 66 : catch(AipsError &x){
788 : // cout << "Hmm.... in here : "<< x.getMesg() << endl;
789 : //no overlap you think ?
790 : /*
791 : os << LogIO::WARN
792 : << "no overlap or failure of copying the clean components"
793 : << x.getMesg()
794 : << LogIO::POST;
795 : */
796 66 : continue;
797 66 : }
798 468 : }
799 46 : }
800 : }
801 :
802 : ////////////////////////////////////////////////////////////////////////////////////////////////////
803 : ///////////////////////////////////////////////////////////////////////////////////////////////////
804 694 : Long SIMapperCollection::estimateRAM(){
805 694 : Long mem=0;
806 1440 : for (uInt k=0; k < itsMappers.nelements(); ++k)
807 : {
808 746 : if(! ((itsMappers[k])->getFTM2()))
809 0 : throw(AipsError("No VI/VB2 FTMachine set"));
810 : ///IFT
811 746 : if(((itsMappers[k])->getFTM2())->estimateRAM(((itsMappers[k])->imageStore()))> 0){
812 746 : mem+=((itsMappers[k])->getFTM2())->estimateRAM(((itsMappers[k])->imageStore()));
813 : //FT
814 746 : mem+=((itsMappers[k])->getFTM2(False))->estimateRAM(((itsMappers[k])->imageStore()));
815 : }
816 : else{
817 : //Assuming double precision...ignoring padding etc.
818 0 : mem+=6*sizeof(float)*(((itsMappers[k])->imageStore())->getShape().product());
819 : }
820 : //Imagestorages
821 746 : mem+=((itsMappers[k])->imageStore())->estimateRAM();
822 : }
823 694 : return mem;
824 : }
825 :
826 : //////////////////////////////////////////////////////////////////////////////////////////////////////
827 : //////////////////////////////////////////////////////////////////////////////////////////////////////
828 :
829 1454 : Bool SIMapperCollection::anyNonZeroModels()
830 : {
831 1454 : Bool validmodel=false;
832 : // If any one Mapper has a valid and nonzero model, return true.
833 2983 : for (Int model=0;model<nMappers(); ++model)
834 : {
835 1529 : validmodel |= (! ( ((itsMappers[model])->imageStore())->isModelEmpty() ));
836 : }
837 1454 : return validmodel;
838 : }
839 : //////////////////////////////////////////////////////////////////////////////////////////////////////
840 :
841 : } //# NAMESPACE CASA - END
842 :
|
