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: aips2-request@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 0 : SIMapperCollection::SIMapperCollection()
70 : {
71 0 : LogIO os( LogOrigin("SIMapperCollection","Construct a mapperCollection",WHERE) );
72 :
73 0 : itsMappers.resize(0);
74 0 : oldMsId_p=-1;
75 0 : itsIsNonZeroModel=false;
76 :
77 0 : }
78 :
79 : //////////////////////////////////////////////////////////////////////////////////////////////////////
80 : //////////////////////////////////////////////////////////////////////////////////////////////////////
81 :
82 0 : SIMapperCollection::~SIMapperCollection()
83 : {
84 0 : }
85 :
86 : //////////////////////////////////////////////////////////////////////////////////////////////////////
87 : //////////////////////////////////////////////////////////////////////////////////////////////////////
88 :
89 0 : Bool SIMapperCollection::releaseImageLocks()
90 : {
91 0 : Bool validflag=true;
92 0 : for(Int mapperid=0;mapperid<nMappers();mapperid++ )
93 : {
94 0 : validflag &= itsMappers[mapperid]->releaseImageLocks();
95 : }
96 0 : return validflag;
97 : }
98 :
99 : //////////////////////////////////////////////////////////////////////////////////////////////////////
100 : //////////////////////////////////////////////////////////////////////////////////////////////////////
101 :
102 0 : std::vector<String> SIMapperCollection::cleanupTempFiles(const String& mess)
103 : {
104 0 : std::vector<String> outstr;
105 0 : auto appvectors = [](std::vector<String>& a, const std::vector<String> b) { a.insert(std::end(a), std::begin(b), std::end(b));};
106 0 : for(Int mapperid=0;mapperid<nMappers();mapperid++ )
107 : {
108 0 : if((itsMappers[mapperid]->getFTM2())){
109 0 : appvectors( outstr, (itsMappers[mapperid]->getFTM2(True)->cleanupTempFiles(mess)).tovector());
110 0 : appvectors(outstr, (itsMappers[mapperid]->getFTM2(False)->cleanupTempFiles(mess)).tovector());
111 : }
112 : }
113 0 : return outstr;
114 : }
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 0 : void SIMapperCollection::addMapper( CountedPtr<SIMapper> map){
159 0 : Int nMappers = itsMappers.nelements();
160 0 : itsMappers.resize(nMappers+1, true);
161 0 : itsMappers[nMappers]=map;
162 0 : }
163 :
164 0 : Int SIMapperCollection::nMappers()
165 : {
166 0 : 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 : }
183 :
184 : //////////////////////////////////////////////////////////////////////////////////////////////////////
185 : /////////////////// Start VB dependent code /////////////////////////////////////////////
186 0 : void SIMapperCollection::initializeGrid(vi::VisBuffer2& vb, Bool dopsf, const Int mapperid)
187 : {
188 0 : if(mapperid<0)
189 : {
190 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
191 : {
192 0 : (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 0 : }
202 :
203 0 : void SIMapperCollection::initializeGrid(vi::VisibilityIterator2& vi, Bool dopsf, const Int mapperid)
204 : {
205 :
206 0 : vi::VisBuffer2 *vb=vi.getVisBuffer();
207 0 : initializeGrid(*vb, dopsf, mapperid);
208 0 : if(mapperid<0)
209 : {
210 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
211 : {
212 0 : ((itsMappers[k])->getFTM2())->initBriggsWeightor(vi);
213 : }
214 : }
215 : else
216 : {
217 0 : if (mapperid > (Int)itsMappers.nelements())
218 0 : throw ( AipsError("Internal Error : SIMapperCollection::initializeGrid(): mapperid out of range") );
219 0 : else (itsMappers[mapperid]->getFTM2())->initBriggsWeightor(vi);
220 : }
221 0 : }
222 :
223 : ////////////////////////////////////////////////////////////////////////////////////
224 : /////////////////////////////////OLD vi/vb //////////////////////////////////////////////
225 0 : void SIMapperCollection::initializeGrid(VisBuffer& vb, Bool dopsf, const Int mapperid)
226 : {
227 0 : if(mapperid<0)
228 : {
229 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
230 : {
231 0 : (itsMappers[k])->initializeGrid(vb,dopsf,true);
232 : }
233 : }
234 : else
235 : {
236 0 : if (mapperid > (Int)itsMappers.nelements())
237 0 : throw ( AipsError("Internal Error : SIMapperCollection::initializeGrid(): mapperid out of range") );
238 0 : else itsMappers[mapperid]->initializeGrid(vb, dopsf, true);
239 : }
240 0 : }
241 :
242 : //////////////////////////////////////////////////////////////////////////////
243 :
244 0 : void SIMapperCollection::grid(vi::VisBuffer2& vb, Bool dopsf, refim::FTMachine::Type col,
245 : Int mapperid)
246 : {
247 0 : if( itsIsNonZeroModel == true ) // Try to subtract model visibilities only if a model exists.
248 : {
249 0 : if(col==refim::FTMachine::CORRECTED){
250 : //Dang i thought the new vb will return Data or FloatData if correctedData was
251 : //not there
252 0 : if(!vb.existsColumn(VisBufferComponent2::VisibilityCorrected)){
253 0 : col=refim::FTMachine::OBSERVED;
254 : // cerr << "Max of visCube" << max(vb.visCube()) << " model " << max(vb.modelVisCube())<< endl;
255 0 : vb.setVisCube(vb.visCube()-vb.visCubeModel());
256 : }
257 : else{
258 0 : vb.setVisCubeCorrected(vb.visCubeCorrected()-vb.visCubeModel());
259 : }
260 : }
261 0 : else if (col==refim::FTMachine::OBSERVED) {
262 0 : vb.setVisCube(vb.visCube()-vb.visCubeModel());
263 : }
264 : }// if non zero model
265 :
266 0 : if(col==refim::FTMachine::CORRECTED && !vb.existsColumn(VisBufferComponent2::VisibilityCorrected)) {
267 : //cout << "Corrected column isn't there, using data instead" << endl;
268 0 : col=refim::FTMachine::OBSERVED;
269 : }
270 :
271 0 : if (mapperid < 0)
272 : {
273 : //cout << "Using column : " << col << endl;
274 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
275 : {
276 0 : (itsMappers[k])->grid(vb, dopsf, col);
277 : }
278 : }
279 : else
280 : {
281 0 : if (mapperid > (Int)itsMappers.nelements())
282 0 : throw ( AipsError("Internal Error : SIMapperCollection::grid(): mapperid out of range") );
283 0 : else itsMappers[mapperid]->grid(vb, dopsf, col);
284 : }
285 0 : }
286 :
287 : ///////////////////
288 : ///////////////////////////////////////OLD VI/VB ///////////////////////////////////
289 0 : void SIMapperCollection::grid(VisBuffer& vb, Bool dopsf, FTMachine::Type col,
290 : Int mapperid)
291 : {
292 0 : if( itsIsNonZeroModel == true ) // Try to subtract model visibilities only if a model exists.
293 : {
294 0 : if(col==FTMachine::CORRECTED){
295 0 : if(vb.msColumns().correctedData().isNull()){
296 0 : col=FTMachine::OBSERVED;
297 : // cerr << "Max of visCube" << max(vb.visCube()) << " model " << max(vb.modelVisCube())<< endl;
298 0 : vb.visCube()-=vb.modelVisCube();
299 : }
300 : else{
301 0 : vb.correctedVisCube()-=vb.modelVisCube();
302 : }
303 : }
304 0 : else if (col==FTMachine::OBSERVED) {
305 0 : vb.visCube()-=vb.modelVisCube();
306 : }
307 : }// if non zero model
308 :
309 0 : if(col==FTMachine::CORRECTED && vb.msColumns().correctedData().isNull())
310 0 : { col=FTMachine::OBSERVED;}
311 :
312 0 : if (mapperid < 0)
313 : {
314 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
315 : {
316 0 : (itsMappers[k])->grid(vb, dopsf, col);
317 :
318 : }
319 : }
320 : else
321 : {
322 0 : if (mapperid > (Int)itsMappers.nelements())
323 0 : throw ( AipsError("Internal Error : SIMapperCollection::grid(): mapperid out of range") );
324 0 : else itsMappers[mapperid]->grid(vb, dopsf, col);
325 : }
326 0 : }
327 : ///////////////////////////////
328 : ////////////////////////////////
329 0 : void SIMapperCollection::finalizeGrid(vi::VisBuffer2& vb, Bool dopsf,const Int mapperid)
330 : {
331 0 : if(mapperid<0)
332 : {
333 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
334 : {
335 0 : (itsMappers[k])->finalizeGrid(vb, dopsf);
336 : }
337 : }
338 : else
339 : {
340 0 : if (mapperid > (Int)itsMappers.nelements())
341 0 : throw ( AipsError("Internal Error : SIMapperCollection::finalizeGrid(): mapperid out of range") );
342 0 : else itsMappers[mapperid]->finalizeGrid(vb, dopsf);
343 : }
344 0 : }
345 :
346 :
347 : ////////////////////////////////////////////////////////////////////////////////////////////////////////
348 : /////////////////////////////////////////////OLD VI/VB////////////////////////////////////////////////
349 0 : void SIMapperCollection::finalizeGrid(VisBuffer& vb, Bool dopsf,const Int mapperid)
350 : {
351 0 : if(mapperid<0)
352 : {
353 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
354 : {
355 0 : (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 0 : }
365 :
366 :
367 : //////////////////////////////////////////////////////////////////////////////////////////////////////
368 0 : void SIMapperCollection::initializeDegrid(vi::VisBuffer2& vb, const Int mapperid)
369 : {
370 :
371 0 : itsIsNonZeroModel = anyNonZeroModels();
372 :
373 0 : if( itsIsNonZeroModel == true )
374 : {
375 : // vb.setModelVisCube( Complex(0.0,0.0) );
376 :
377 0 : if(mapperid<0)
378 : {
379 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
380 : {
381 0 : (itsMappers[k])->initializeDegrid(vb);
382 : }
383 : }
384 : else
385 : {
386 0 : if (mapperid > (Int)itsMappers.nelements())
387 0 : throw ( AipsError("Internal Error : SIMapperCollection::initializeDegrid(): mapperid out of range") );
388 0 : else itsMappers[mapperid]->initializeDegrid(vb);
389 : }
390 :
391 : }// if non zero model
392 0 : }
393 :
394 :
395 :
396 : ///////////////////////////////////////OLD VI/VB ///////////////////////////////////////////////////
397 0 : void SIMapperCollection::initializeDegrid(VisBuffer& vb, const Int mapperid)
398 : {
399 :
400 0 : itsIsNonZeroModel = anyNonZeroModels();
401 :
402 0 : if( itsIsNonZeroModel == true )
403 : {
404 : // vb.setModelVisCube( Complex(0.0,0.0) );
405 :
406 0 : if(mapperid<0)
407 : {
408 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
409 : {
410 0 : (itsMappers[k])->initializeDegrid(vb);
411 : }
412 : }
413 : else
414 : {
415 0 : if (mapperid > (Int)itsMappers.nelements())
416 0 : throw ( AipsError("Internal Error : SIMapperCollection::initializeDegrid(): mapperid out of range") );
417 0 : else itsMappers[mapperid]->initializeDegrid(vb);
418 : }
419 :
420 : }// if non zero model
421 0 : }
422 :
423 : ////////////////////////////////////////////////////////////////////////////////////////////////
424 0 : void SIMapperCollection::degrid(vi::VisBuffer2& vb, Bool saveVirtualMod, const Int mapperid)
425 : {
426 0 : if( itsIsNonZeroModel == true )
427 : {
428 0 : if(mapperid<0)
429 : {
430 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
431 : {
432 0 : (itsMappers[k])->degrid(vb);
433 : }
434 : }
435 : else
436 : {
437 0 : if (mapperid > (Int)itsMappers.nelements())
438 0 : throw ( AipsError("Internal Error : SIMapperCollection::degrid(): mapperid out of range") );
439 0 : else itsMappers[mapperid]->degrid(vb);
440 : }
441 :
442 0 : if(saveVirtualMod){
443 0 : saveVirtualModel(vb);
444 : }
445 : }// if non zero model
446 0 : }
447 :
448 0 : void SIMapperCollection::addPB(vi::VisBuffer2& vb, PBMath& pbMath, const MDirection& altDir, const Bool useAltDir)
449 : {
450 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
451 : {
452 0 : (itsMappers[k])->addPB(vb,pbMath, altDir, useAltDir);
453 :
454 : }
455 0 : }
456 :
457 :
458 : /////////////////////////////////////OLD VI/VB ////////////////////////////////////////////////////
459 0 : void SIMapperCollection::degrid(VisBuffer& vb, Bool saveVirtualMod, const Int mapperid)
460 : {
461 0 : if( itsIsNonZeroModel == true )
462 : {
463 0 : if(mapperid<0)
464 : {
465 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
466 : {
467 0 : (itsMappers[k])->degrid(vb);
468 : }
469 : }
470 : else
471 : {
472 0 : if (mapperid > (Int)itsMappers.nelements())
473 0 : throw ( AipsError("Internal Error : SIMapperCollection::degrid(): mapperid out of range") );
474 0 : else itsMappers[mapperid]->degrid(vb);
475 : }
476 :
477 0 : if(saveVirtualMod){
478 0 : saveVirtualModel(vb);
479 : }
480 : }// if non zero model
481 0 : }
482 :
483 :
484 : /////
485 : /////////////
486 0 : void SIMapperCollection::saveVirtualModel(VisBuffer& vb){
487 :
488 :
489 :
490 0 : if(vb.msId() != oldMsId_p){
491 0 : oldMsId_p=vb.msId();
492 : /*Block< Vector<Int> > blockNGroup;
493 : Block< Vector<Int> > blockStart;
494 : Block< Vector<Int> > blockWidth;
495 : Block< Vector<Int> > blockIncr;
496 : Block< Vector<Int> > blockSpw;
497 : vb.getChannelSelection(blockNGroup, blockStart, blockWidth, blockIncr, blockSpw);
498 : Vector<Int> fields = vb.msColumns().fieldId().getColumn();
499 : const Int option = Sort::HeapSort | Sort::NoDuplicates;
500 : const Sort::Order order = Sort::Ascending;
501 : Int nfields = GenSort<Int>::sort (fields, order, option);
502 :
503 : // Make sure we have the right size
504 :
505 : fields.resize(nfields, true);
506 : */
507 : //Int msid = vb.msId();
508 0 : ROVisibilityIterator *viloc=vb.getVisibilityIterator();
509 0 : for (uInt k=0; k < itsMappers.nelements(); ++k){
510 0 : Record rec;
511 0 : String modImage=viloc->ms().getPartNames()[0];
512 0 : if(!((viloc->ms()).source().isNull()))
513 0 : modImage=(viloc->ms()).source().tableName();
514 0 : modImage=File::newUniqueName(modImage, "FT_MODEL").absoluteName();
515 0 : Bool iscomp=itsMappers[k]->getCLRecord(rec);
516 0 : if(iscomp || itsMappers[k]->getFTMRecord(rec, modImage)){
517 :
518 : //VisModelData::putModel(vb.getVisibilityIterator()->ms(), rec, fields, blockSpw[msid], blockStart[msid],
519 : // blockWidth[msid], blockIncr[msid],
520 : // iscomp, true);
521 0 : VisibilityIterator * elvi=(dynamic_cast<VisibilityIterator* >(vb.getVisibilityIterator()));
522 0 : if(elvi)
523 0 : elvi->putModel(rec, iscomp, true);
524 : // VisModelData::listModel(vb.getVisibilityIterator()->ms());
525 : }
526 :
527 : }
528 :
529 :
530 :
531 : }
532 :
533 :
534 :
535 0 : }
536 :
537 0 : void SIMapperCollection::saveVirtualModel(vi::VisBuffer2& vb){
538 :
539 :
540 :
541 0 : if(vb.msId() != oldMsId_p){
542 0 : oldMsId_p=vb.msId();
543 :
544 0 : for (uInt k=0; k < itsMappers.nelements(); ++k){
545 0 : Record rec;
546 0 : String modImage=vb.ms().getPartNames()[0];
547 0 : if(!((vb.ms()).source().isNull()))
548 0 : modImage=(vb.ms()).source().tableName();
549 0 : modImage=File::newUniqueName(modImage, "FT_MODEL").absoluteName();
550 0 : Bool iscomp=itsMappers[k]->getCLRecord(rec);
551 0 : if(iscomp || itsMappers[k]->getFTMRecord(rec, modImage)){
552 :
553 : ////Darn not implemented
554 : //static_cast<VisibilityIteratorImpl2 *>(viloc->getImpl())->writeModel(rec, //iscomp, true);
555 :
556 0 : if(!iscomp && Table::isReadable(modImage)){
557 : //make sure complex image is of compliant size/shape
558 0 : (itsMappers[k]->imageStore())->intersectComplexImage(modImage);
559 :
560 : }
561 0 : VisibilityIterator2* vi=const_cast<VisibilityIterator2*>(vb.getVi());
562 0 : const_cast<MeasurementSet& >(vi->ms()).lock();
563 : /////TESTOO
564 : //Int CPUID;
565 : //MPI_Comm_rank(MPI_COMM_WORLD, &CPUID);
566 : //cerr << CPUID << " writing " << modImage << endl;
567 : /////////////////
568 0 : vi->writeModel(rec, iscomp, true);
569 0 : const_cast<MeasurementSet& >(vi->ms()).unlock();
570 : // VisModelData::listModel(vb.getVisibilityIterator()->ms());
571 : }
572 :
573 : }
574 :
575 :
576 :
577 : }
578 :
579 :
580 :
581 0 : }
582 :
583 :
584 : /////////////////////////////////////////////////////////
585 0 : void SIMapperCollection::finalizeDegrid(vi::VisBuffer2& /*vb*/, const Int mapperid)
586 : {
587 0 : if( itsIsNonZeroModel == true )
588 : {
589 0 : if(mapperid<0)
590 : {
591 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
592 : {
593 0 : (itsMappers[k])->finalizeDegrid();
594 :
595 : }
596 : }
597 : else
598 : {
599 0 : if (mapperid > (Int)itsMappers.nelements())
600 0 : throw ( AipsError("Internal Error : SIMapperCollection::finalizeDegrid(): mapperid out of range") );
601 0 : else itsMappers[mapperid]->finalizeDegrid();
602 : }
603 : }// if non zero model
604 0 : }
605 :
606 : ////////////////////////////////////////////////////////
607 0 : void SIMapperCollection::finalizeDegrid(VisBuffer& /*vb*/, const Int mapperid)
608 : {
609 0 : if( itsIsNonZeroModel == true )
610 : {
611 0 : if(mapperid<0)
612 : {
613 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
614 : {
615 0 : (itsMappers[k])->finalizeDegrid();
616 :
617 : }
618 : }
619 : else
620 : {
621 0 : if (mapperid > (Int)itsMappers.nelements())
622 0 : throw ( AipsError("Internal Error : SIMapperCollection::finalizeDegrid(): mapperid out of range") );
623 0 : else itsMappers[mapperid]->finalizeDegrid();
624 : }
625 : }// if non zero model
626 0 : }
627 :
628 :
629 0 : void SIMapperCollection::initPB()
630 : {
631 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
632 : {
633 0 : (itsMappers[k])->initPB();
634 : }
635 0 : }
636 :
637 0 : void SIMapperCollection::addPB(VisBuffer& vb, PBMath& pbMath)
638 : {
639 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
640 : {
641 0 : (itsMappers[k])->addPB(vb,pbMath);
642 :
643 : }
644 0 : }
645 :
646 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
647 : //////////// End of VB dependent code.
648 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
649 :
650 :
651 :
652 : ///////////////////////////////////////////////////////////////////////////////////////////////////////
653 : /////////////////////
654 0 : CountedPtr<SIImageStore> SIMapperCollection::imageStore(Int id)
655 : {
656 0 : if(uInt(id) < itsMappers.nelements())
657 : {
658 0 : return (itsMappers[id])->imageStore();
659 : }
660 0 : return make_shared<SIImageStore>( );
661 : }
662 : /////////////////////////////////////////////////////////////////////////////////////////////
663 : /////////////////////////////////////////////////////////////////////////////////////////////
664 0 : Record SIMapperCollection::getFTMRecord(Int mapperid)
665 : {
666 0 : AlwaysAssert( mapperid >=0 && mapperid < nMappers() , AipsError );
667 : //return itsMappers[mapperid]->getFTMRecord();
668 0 : Record rec;
669 0 : return rec;
670 : }
671 :
672 : //////////////////////////////////////////////////////////////////////////////////////////////////////
673 : //////////////////////////////////////////////////////////////////////////////////////////////////////
674 :
675 0 : void SIMapperCollection::checkOverlappingModels(String action)
676 : {
677 : // If nothing that overlaps, don't check.
678 0 : if(nMappers()==1) return;
679 :
680 0 : Int nmodels = nMappers();
681 :
682 : // If there is no model image (i.e. first major cycle with no starting model), don't check.
683 0 : Bool hasmodel=true;
684 0 : for (Int model=0;model<(nmodels-1); ++model)
685 0 : { hasmodel = hasmodel && ((itsMappers[model])->imageStore())->hasModel(); }
686 0 : if( hasmodel==false ) {
687 : //cout << "No model images to check overlap for." << endl;
688 0 : return;
689 : }
690 :
691 : // Internal check
692 0 : AlwaysAssert( action=="blank" || action=="restore" , AipsError );
693 :
694 0 : for (Int model=0;model<(nmodels-1); ++model)
695 : {
696 : // Connect to one image for aux info.
697 0 : LatticeLocker lockmodel (*(((itsMappers[model])->imageStore())->model()), FileLocker::Write);
698 0 : SubImage<Float> modelimage( *(((itsMappers[model])->imageStore())->model()), true );
699 :
700 0 : uInt nTaylor0 = ((itsMappers[model])->imageStore())->getNTaylorTerms();
701 :
702 0 : CoordinateSystem cs0=modelimage.coordinates();
703 0 : IPosition iblc0(modelimage.shape().nelements(),0);
704 0 : IPosition itrc0(modelimage.shape());
705 0 : itrc0=itrc0-Int(1);
706 0 : LCBox lbox0(iblc0, itrc0, modelimage.shape());
707 0 : ImageRegion imagreg0(WCBox(lbox0, cs0));
708 :
709 0 : for (Int nextmodel=model+1; nextmodel < nmodels; ++nextmodel)
710 : {
711 0 : LatticeLocker nextlockmodel (*(((itsMappers[nextmodel])->imageStore())->model()), FileLocker::Write);
712 0 : SubImage<Float> nextmodelimage( *(((itsMappers[nextmodel])->imageStore())->model()) , true);
713 :
714 0 : uInt nTaylor1 = ((itsMappers[nextmodel])->imageStore())->getNTaylorTerms();
715 :
716 0 : CoordinateSystem cs=nextmodelimage.coordinates();
717 0 : IPosition iblc(nextmodelimage.shape().nelements(),0);
718 0 : IPosition itrc(nextmodelimage.shape());
719 0 : itrc=itrc-Int(1);
720 0 : LCBox lbox(iblc, itrc, nextmodelimage.shape());
721 0 : ImageRegion imagreg(WCBox(lbox, cs));
722 :
723 : try{
724 :
725 0 : if( action.matches("blank") )
726 : {
727 :
728 : //cerr << "blank MODEL image shape " << modelimage.shape() << " " << nextmodelimage.shape() << endl;
729 :
730 0 : LatticeRegion latReg=imagreg.toLatticeRegion(modelimage.coordinates(), modelimage.shape());
731 :
732 0 : for(uInt taylor=0;taylor<min(nTaylor0,nTaylor1);taylor++)
733 : { // loop for taylor term
734 0 : SubImage<Float> modelim( *(((itsMappers[model])->imageStore())->model(taylor)), true );
735 0 : SubImage<Float> partToMask(modelim, imagreg, true);
736 0 : LatticeLocker lock1 (partToMask, FileLocker::Write);
737 0 : ArrayLattice<Bool> pixmask(latReg.get());
738 0 : LatticeExpr<Float> myexpr(iif(pixmask, 0.0, partToMask) );
739 0 : partToMask.copyData(myexpr);
740 : }
741 :
742 :
743 : }
744 : else // "restore"
745 : {
746 : //cerr << "rsetore MODEL image shape " << modelimage.shape() << " " << nextmodelimage.shape() << endl;
747 0 : LatticeRegion latReg0=imagreg0.toLatticeRegion(nextmodelimage.coordinates(), nextmodelimage.shape());
748 0 : LatticeRegion latReg=imagreg.toLatticeRegion(modelimage.coordinates(), modelimage.shape());
749 0 : ArrayLattice<Bool> pixmask(latReg.get());
750 :
751 :
752 0 : for(uInt taylor=0;taylor<min(nTaylor0,nTaylor1);taylor++)
753 : {// loop for taylor term
754 0 : SubImage<Float> modelim( *(((itsMappers[model])->imageStore())->model(taylor)), true );
755 0 : SubImage<Float> nextmodelim( *(((itsMappers[nextmodel])->imageStore())->model(taylor)), true );
756 :
757 0 : SubImage<Float> partToMerge(nextmodelim, imagreg0, true);
758 0 : SubImage<Float> partToUnmask(modelim, imagreg, true);
759 0 : LatticeLocker lock1 (partToUnmask, FileLocker::Write);
760 0 : LatticeLocker lock2 (partToMerge, FileLocker::Write);
761 0 : LatticeExpr<Float> myexpr0(iif(pixmask,partToMerge,partToUnmask));
762 0 : partToUnmask.copyData(myexpr0);
763 : }
764 :
765 : }
766 : }
767 0 : catch(AipsError &x){
768 : // cout << "Hmm.... in here : "<< x.getMesg() << endl;
769 : //no overlap you think ?
770 : /*
771 : os << LogIO::WARN
772 : << "no overlap or failure of copying the clean components"
773 : << x.getMesg()
774 : << LogIO::POST;
775 : */
776 0 : continue;
777 : }
778 : }
779 : }
780 : }
781 :
782 : ////////////////////////////////////////////////////////////////////////////////////////////////////
783 : ///////////////////////////////////////////////////////////////////////////////////////////////////
784 0 : Long SIMapperCollection::estimateRAM(){
785 0 : Long mem=0;
786 0 : for (uInt k=0; k < itsMappers.nelements(); ++k)
787 : {
788 0 : if(! ((itsMappers[k])->getFTM2()))
789 0 : throw(AipsError("No VI/VB2 FTMachine set"));
790 : ///IFT
791 0 : if(((itsMappers[k])->getFTM2())->estimateRAM(((itsMappers[k])->imageStore()))> 0){
792 0 : mem+=((itsMappers[k])->getFTM2())->estimateRAM(((itsMappers[k])->imageStore()));
793 : //FT
794 0 : mem+=((itsMappers[k])->getFTM2(False))->estimateRAM(((itsMappers[k])->imageStore()));
795 : }
796 : else{
797 : //Assuming double precision...ignoring padding etc.
798 0 : mem+=6*sizeof(float)*(((itsMappers[k])->imageStore())->getShape().product());
799 : }
800 : //Imagestorages
801 0 : mem+=((itsMappers[k])->imageStore())->estimateRAM();
802 : }
803 0 : return mem;
804 : }
805 :
806 : //////////////////////////////////////////////////////////////////////////////////////////////////////
807 : //////////////////////////////////////////////////////////////////////////////////////////////////////
808 :
809 0 : Bool SIMapperCollection::anyNonZeroModels()
810 : {
811 0 : Bool validmodel=false;
812 : // If any one Mapper has a valid and nonzero model, return true.
813 0 : for (Int model=0;model<nMappers(); ++model)
814 : {
815 0 : validmodel |= (! ( ((itsMappers[model])->imageStore())->isModelEmpty() ));
816 : }
817 0 : return validmodel;
818 : }
819 : //////////////////////////////////////////////////////////////////////////////////////////////////////
820 :
821 : } //# NAMESPACE CASA - END
822 :
|
