Line data Source code
1 : //# SynthesisImagerVi2.cc: Implementation of SynthesisImager.h
2 : //# Copyright (C) 1997-2019
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //# This library is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU General Public License as published by
6 : //# the Free Software Foundation; either version 3 of the License, or (at your
7 : //# option) any later version.
8 : //#
9 : //# This library is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
12 : //# License for more details.
13 : //#
14 : //# https://www.gnu.org/licenses/
15 : //#
16 : //# You should have received a copy of the GNU General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Queries concerning CASA should be submitted at
21 : //# https://help.nrao.edu
22 : //#
23 : //# Postal address: CASA Project Manager
24 : //# National Radio Astronomy Observatory
25 : //# 520 Edgemont Road
26 : //# Charlottesville, VA 22903-2475 USA
27 : //#
28 : //#
29 : //# $Id$
30 :
31 : #define CFC_VERBOSE false /* Control the verbosity when building CFCache. */
32 :
33 : #include <casacore/casa/Exceptions/Error.h>
34 : #include <iostream>
35 : #include <sstream>
36 :
37 : #include <casacore/casa/Arrays/Matrix.h>
38 : #include <casacore/casa/Arrays/ArrayMath.h>
39 : #include <casacore/casa/Arrays/ArrayLogical.h>
40 :
41 :
42 : #include <casacore/casa/Logging.h>
43 : #include <casacore/casa/Logging/LogIO.h>
44 : #include <casacore/casa/Logging/LogMessage.h>
45 : #include <casacore/casa/Logging/LogSink.h>
46 : #include <casacore/casa/Logging/LogMessage.h>
47 : #include <casacore/casa/System/ProgressMeter.h>
48 :
49 : #include <casacore/casa/OS/DirectoryIterator.h>
50 : #include <casacore/casa/OS/File.h>
51 : #include <casacore/casa/OS/HostInfo.h>
52 : #include <casacore/casa/OS/Path.h>
53 : //#include <casa/OS/Memory.h>
54 :
55 : #include <casacore/lattices/LRegions/LCBox.h>
56 :
57 : #include <casacore/measures/Measures/MeasTable.h>
58 :
59 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
60 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
61 : #include <casacore/ms/MSSel/MSSelection.h>
62 :
63 :
64 : #include <synthesis/ImagerObjects/SynthesisImagerVi2.h>
65 :
66 : #include <synthesis/ImagerObjects/SynthesisUtilMethods.h>
67 : #include <synthesis/ImagerObjects/SIImageStore.h>
68 : #include <synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
69 : #include <synthesis/ImagerObjects/CubeMajorCycleAlgorithm.h>
70 : #include <synthesis/ImagerObjects/CubeMakeImageAlgorithm.h>
71 : #include <synthesis/MeasurementEquations/VPManager.h>
72 : #include <imageanalysis/Utilities/SpectralImageUtil.h>
73 : #include <msvis/MSVis/MSUtil.h>
74 : #include <msvis/MSVis/VisSetUtil.h>
75 : #include <msvis/MSVis/VisImagingWeight.h>
76 :
77 : #include <synthesis/TransformMachines2/GridFT.h>
78 : #include <synthesis/TransformMachines2/WPConvFunc.h>
79 : #include <synthesis/TransformMachines2/WProjectFT.h>
80 : #include <synthesis/TransformMachines2/VisModelData.h>
81 : #include <synthesis/TransformMachines2/AWProjectFT.h>
82 : #include <synthesis/TransformMachines2/HetArrayConvFunc.h>
83 : #include <synthesis/TransformMachines2/MosaicFTNew.h>
84 : #include <synthesis/TransformMachines2/AWPLPG.h>
85 : #include <synthesis/TransformMachines2/MultiTermFTNew.h>
86 : #include <synthesis/TransformMachines2/AWProjectWBFTNew.h>
87 : #include <synthesis/TransformMachines2/AWConvFunc.h>
88 : //#include <synthesis/TransformMachines2/AWConvFuncEPJones.h>
89 : #include <synthesis/TransformMachines2/NoOpATerm.h>
90 : #include <synthesis/TransformMachines2/SDGrid.h>
91 : #include <synthesis/TransformMachines/WProjectFT.h>
92 : #include <synthesis/TransformMachines2/BriggsCubeWeightor.h>
93 : #include <casacore/casa/Utilities/Regex.h>
94 : #include <casacore/casa/OS/Directory.h>
95 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
96 : //#include <casadbus/utilities/BusAccess.h>
97 : //#include <casadbus/session/DBusSession.h>
98 :
99 : #include <sys/types.h>
100 : #include <unistd.h>
101 : #include <iomanip>
102 : #include <thread>
103 : #include <synthesis/Parallel/Applicator.h>
104 :
105 : using namespace std;
106 :
107 : using namespace casacore;
108 :
109 : namespace casa { //# NAMESPACE CASA - BEGIN
110 : extern Applicator applicator;
111 1711 : SynthesisImagerVi2::SynthesisImagerVi2() : SynthesisImager(), vi_p(0), fselections_p(nullptr), imparsVec_p(0), gridparsVec_p(0) {
112 : /*cerr << "is applicator initialized " << applicator.initialized() << endl;
113 : if(!applicator.initialized()){
114 : int argc=1;
115 : char **argv;
116 : casa::applicator.init ( argc, argv );
117 : cerr << "controller ?" << applicator.isController() << " worker? " << applicator.isWorker() << " numprocs " << applicator.numProcs() << endl;
118 : }
119 : */
120 1711 : mss_p.resize(0);
121 1711 : }
122 :
123 2594 : SynthesisImagerVi2::~SynthesisImagerVi2(){
124 3500 : for (uInt k=0; k < mss_p.nelements(); ++k){
125 1789 : if(mss_p[k])
126 1789 : delete mss_p[k];
127 : }
128 1711 : SynthesisUtilMethods::getResource("End Run");
129 2594 : }
130 :
131 1804 : Bool SynthesisImagerVi2::selectData(const SynthesisParamsSelect& selpars){
132 3608 : LogIO os( LogOrigin("SynthesisImagerVi2","selectData",WHERE) );
133 1804 : Bool retval=True;
134 :
135 1804 : SynthesisUtilMethods::getResource("Start Run");
136 :
137 : try
138 : {
139 :
140 :
141 1804 : MeasurementSet thisms;
142 : { ///Table system seems to have a bug when running in multi-process as the source table disappears
143 : /// temporarily when other processes are updating it
144 1804 : uInt exceptCounter=0;
145 :
146 : while(true){
147 : try{
148 : //Respect the readonly flag...necessary for multi-process access
149 3608 : thisms=MeasurementSet(selpars.msname, TableLock(TableLock::UserNoReadLocking),
150 3608 : selpars.readonly ? Table::Old : Table::Update);
151 1804 : break;
152 : }
153 0 : catch(AipsError &x){
154 :
155 0 : String mes=x.getMesg();
156 0 : if(mes.contains("FilebufIO::readBlock") || mes.contains("SOURCE")){
157 0 : std::this_thread::sleep_for(std::chrono::milliseconds(50));
158 0 : os << LogIO::WARN << "#####CATCHING a sleep because "<< mes<< LogIO::POST;
159 : }
160 : else
161 0 : throw(AipsError("Error in selectdata: "+mes));
162 :
163 0 : if(exceptCounter > 100){
164 0 : throw(AipsError("Error in selectdata got 100 of this exeception: "+mes));
165 :
166 : }
167 :
168 0 : }
169 0 : ++exceptCounter;
170 0 : }
171 : }//End of work around for table disappearing bug
172 :
173 :
174 :
175 1804 : useScratch_p=selpars.usescratch;
176 1804 : readOnly_p = selpars.readonly;
177 1804 : lockMS(thisms);
178 1804 : thisms.setMemoryResidentSubtables (MrsEligibility::defaultEligible());
179 :
180 :
181 : // cout << "**************** usescr : " << useScratch_p << " readonly : " << readOnly_p << endl;
182 : //if you want to use scratch col...make sure they are there
183 : ///Need to clear this in first pass only...child processes or loops for cube ///should skip it
184 1804 : if(!(impars_p.mode.contains("cube")) || ((impars_p.mode.contains("cube")) && doingCubeGridding_p)){
185 1793 : if(selpars.usescratch && !selpars.readonly){
186 38 : VisSetUtil::addScrCols(thisms, true, false, true, false);
187 38 : refim::VisModelData::clearModel(thisms);
188 : }
189 : ////TESTOO
190 : //Int CPUID;
191 : //MPI_Comm_rank(MPI_COMM_WORLD, &CPUID);
192 : //cerr << " SELPARS " << selpars.toRecord() << endl;
193 1793 : if(!selpars.incrmodel && !selpars.usescratch && !selpars.readonly)
194 17 : refim::VisModelData::clearModel(thisms, selpars.field, selpars.spw);
195 : }//end of first pass if for cubes
196 : // unlockMSs();
197 :
198 :
199 1804 : os << "MS : " << selpars.msname << " | ";
200 :
201 : //Some MSSelection
202 : //If everything is empty (which is valid) it will throw an exception..below
203 : //So make sure the main defaults are not empy i.e field and spw
204 1804 : MSSelection thisSelection;
205 1804 : if(selpars.field != ""){
206 615 : thisSelection.setFieldExpr(selpars.field);
207 615 : os << "Selecting on fields : " << selpars.field << " | " ;//LogIO::POST;
208 : }else
209 1189 : thisSelection.setFieldExpr("*");
210 1804 : if(selpars.spw != ""){
211 665 : thisSelection.setSpwExpr(selpars.spw);
212 665 : os << "Selecting on spw :"<< selpars.spw << " | " ;//LogIO::POST;
213 : }else
214 1139 : thisSelection.setSpwExpr("*");
215 :
216 1804 : if(selpars.antenna != ""){
217 91 : Vector<String> antNames(1, selpars.antenna);
218 : // thisSelection.setAntennaExpr(MSSelection::nameExprStr( antNames));
219 91 : thisSelection.setAntennaExpr(selpars.antenna);
220 91 : os << "Selecting on antenna names : " << selpars.antenna << " | " ;//LogIO::POST;
221 :
222 91 : }
223 1804 : if(selpars.timestr != ""){
224 14 : thisSelection.setTimeExpr(selpars.timestr);
225 14 : os << "Selecting on time range : " << selpars.timestr << " | " ;//LogIO::POST;
226 : }
227 1804 : if(selpars.uvdist != ""){
228 0 : thisSelection.setUvDistExpr(selpars.uvdist);
229 0 : os << "Selecting on uvdist : " << selpars.uvdist << " | " ;//LogIO::POST;
230 : }
231 1804 : if(selpars.scan != ""){
232 14 : thisSelection.setScanExpr(selpars.scan);
233 14 : os << "Selecting on scan : " << selpars.scan << " | " ;//LogIO::POST;
234 : }
235 1804 : if(selpars.obs != ""){
236 0 : thisSelection.setObservationExpr(selpars.obs);
237 0 : os << "Selecting on Observation Expr : " << selpars.obs << " | " ;//LogIO::POST;
238 : }
239 1804 : if(selpars.state != ""){
240 72 : thisSelection.setStateExpr(selpars.state);
241 72 : os << "Selecting on Scan Intent/State : " << selpars.state << " | " ;//LogIO::POST;
242 : }
243 : // if(selpars.taql != ""){
244 : // thisSelection.setTaQLExpr(selpars.taql);
245 : // os << "Selecting via TaQL : " << selpars.taql << " | " ;//LogIO::POST;
246 : // }
247 1804 : os << "[Opened " << (readOnly_p?"in readonly mode":(useScratch_p?"with scratch model column":"with virtual model column")) << "]" << LogIO::POST;
248 1804 : TableExprNode exprNode=thisSelection.toTableExprNode(&thisms);
249 1801 : if(!(exprNode.isNull()))
250 : {
251 :
252 :
253 1801 : MeasurementSet thisMSSelected0 = MeasurementSet(thisms(exprNode));
254 1801 : mss_p.resize(mss_p.nelements()+1, false, true);
255 1801 : if(selpars.taql != "")
256 : {
257 0 : MSSelection mss0;
258 0 : mss0.setTaQLExpr(selpars.taql);
259 0 : os << "Selecting via TaQL : " << selpars.taql << " | " ;//LogIO::POST;
260 :
261 0 : TableExprNode tenWithTaQL=mss0.toTableExprNode(&thisMSSelected0);
262 0 : MeasurementSet thisMSSelected1 = MeasurementSet(thisMSSelected0(tenWithTaQL));
263 : //mss4vi_p[mss4vi_p.nelements()-1]=MeasurementSet(thisms(exprNode));
264 0 : mss_p[mss_p.nelements()-1]=new MeasurementSet(thisMSSelected1);
265 0 : }
266 : else
267 1801 : mss_p[mss_p.nelements()-1]=new MeasurementSet(thisMSSelected0);
268 :
269 1801 : os << " NRows selected : " << (mss_p[mss_p.nelements()-1])->nrow() << LogIO::POST;
270 1801 : unlockMSs();
271 1801 : }
272 : else{
273 0 : throw(AipsError("Selection for given MS "+selpars.msname+" is invalid"));
274 : }
275 1801 : if((mss_p[mss_p.nelements()-1])->nrow() ==0){
276 1 : delete mss_p[mss_p.nelements()-1];
277 1 : mss_p.resize(mss_p.nelements()-1, True, True);
278 1 : if(mss_p.nelements()==0)
279 1 : throw(AipsError("Data selection ended with 0 rows"));
280 : //Sill have some valid ms's so return false and do not proceed to do
281 : //channel selection
282 0 : unlockMSs();
283 0 : return False;
284 : }
285 :
286 :
287 :
288 : ///Channel selection
289 : {
290 1800 : if(!fselections_p) fselections_p=new FrequencySelections();
291 1800 : Matrix<Int> chanlist = thisSelection.getChanList(mss_p[mss_p.nelements()-1]);
292 1800 : Matrix<Double> freqList=thisSelection.getChanFreqList(mss_p[mss_p.nelements()-1]);
293 : //cerr << std::setprecision(12) << "FreqList " << freqList << endl;
294 1800 : IPosition shape = freqList.shape();
295 1800 : uInt nSelections = shape[0];
296 : ///temporary variable as we carry that for tunechunk...till we get rid of it
297 1800 : selFreqFrame_p=selpars.freqframe;
298 1800 : Bool ignoreframe=False;
299 1800 : MFrequency::Types freqFrame=MFrequency::castType(MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().measFreqRef()(Int(chanlist(0,0))));
300 :
301 1800 : if(selpars.freqframe == MFrequency::REST ||selpars.freqframe == MFrequency::Undefined){
302 21 : selFreqFrame_p=freqFrame;
303 21 : ignoreframe=True;
304 : }
305 1800 : if(selpars.freqbeg==""){
306 : // Going round the problem of CAS-8829
307 : /*vi::FrequencySelectionUsingChannels channelSelector;
308 :
309 : channelSelector.add(thisSelection, mss_p[mss_p.nelements()-1]);
310 :
311 : fselections_p.add(channelSelector);
312 : */
313 : ////////////////////////////
314 : Double lowfreq;
315 : Double topfreq;
316 1789 : Vector<Int> fieldList=thisSelection.getFieldList(mss_p[mss_p.nelements()-1]);
317 : // cerr << "chanlist " << chanlist.column(0) << "\n fieldList " << fieldList << endl;
318 :
319 : //cerr << "selpars.freqframe " << selpars.freqframe << endl;
320 1789 : vi::FrequencySelectionUsingFrame channelSelector(selFreqFrame_p);
321 : ///temporary variable as we carry that for tunechunk
322 :
323 1789 : Bool selectionValid=False;
324 4190 : for(uInt k=0; k < nSelections; ++k){
325 2401 : Bool thisSpwSelValid=False;
326 : //The getChanfreqList is wrong for beg and end..going round that too.
327 4802 : Vector<Double> freqies=MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().chanFreq()(Int(chanlist(k,0)));
328 4802 : Vector<Double> chanwidth=MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().chanWidth()(Int(chanlist(k,0)));
329 :
330 2401 : if(freqList(k,3) < 0.0){
331 94 : topfreq=freqies(chanlist(k,1));//-chanwidth(chanlist(k,1))/2.0;
332 94 : lowfreq=freqies(chanlist(k,2));//+chanwidth(chanlist(k,2))/2.0;
333 : //lowfreq=freqList(k,2); //+freqList(k,3)/2.0;
334 : //topfreq=freqList(k, 1); //-freqList(k,3)/2.0;
335 : }
336 : else{
337 2307 : lowfreq=freqies(chanlist(k,1));//-chanwidth(chanlist(k,1))/2.0;
338 2307 : topfreq=freqies(chanlist(k,2));//+chanwidth(chanlist(k,2))/2.0;
339 : //lowfreq=freqList(k,1);//-freqList(k,3)/2.0;
340 : //topfreq=freqList(k, 2);//+freqList(k,3)/2.0;
341 : }
342 :
343 2401 : if(!ignoreframe){
344 :
345 :
346 : //cerr << "begin " << lowfreq << " " << topfreq << endl;
347 : //vi::VisibilityIterator2 tmpvi(mss_p, vi::SortColumns(), false);
348 : //VisBufferUtil::getFreqRangeFromRange(lowfreq, topfreq, freqFrame, lowfreq, topfreq, tmpvi, selFreqFrame_p);
349 : // lockMS(*(const_cast<MeasurementSet*> (mss_p[mss_p.nelements()-1])));
350 2380 : if(MSUtil::getFreqRangeInSpw( lowfreq,
351 4760 : topfreq, Vector<Int>(1,chanlist(k,0)), Vector<Int>(1,chanlist(k,1)),
352 4760 : Vector<Int>(1, chanlist(k,2)-chanlist(k,1)+1),
353 2380 : *mss_p[mss_p.nelements()-1] ,
354 : selFreqFrame_p,
355 : fieldList, False))
356 : {
357 2362 : selectionValid=True;
358 2362 : thisSpwSelValid=True;
359 : }
360 : // unlockMSs();
361 :
362 : }
363 :
364 2401 : if(thisSpwSelValid || ignoreframe){
365 2383 : andFreqSelection(mss_p.nelements()-1, Int(freqList(k,0)), lowfreq, topfreq, selFreqFrame_p);
366 2383 : andChanSelection(mss_p.nelements()-1, Int(chanlist(k,0)), Int(chanlist(k,1)),Int(chanlist(k,2)));
367 : }
368 2401 : }
369 1789 : if(! (selectionValid && !ignoreframe)){
370 : //os << "Did not match spw selection in the selected ms " << LogIO::WARN << LogIO::POST;
371 21 : retval=False;
372 : }
373 : //fselections_p->add(channelSelector);
374 : //////////////////////////////////
375 1789 : }
376 : else{
377 :
378 : //////More workaroung CAS-8829
379 : //MFrequency::Types freqFrame=MFrequency::castType(MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().measFreqRef()(Int(freqList(0,0))));
380 11 : Quantity freq;
381 11 : Quantity::read(freq, selpars.freqbeg);
382 11 : Double lowfreq=freq.getValue("Hz");
383 11 : Quantity::read(freq, selpars.freqend);
384 11 : Double topfreq=freq.getValue("Hz");
385 :
386 : ////Work aroun CAS-8829
387 : // if(vi_p)
388 : //VisBufferUtil::getFreqRangeFromRange(lowfreq, topfreq, selpars.freqframe, lowfreq, topfreq, *vi_p, freqFrame);
389 : //cerr << "lowFreq "<< lowfreq << " topfreq " << topfreq << endl;
390 : //vi::FrequencySelectionUsingFrame channelSelector((vi_p ? freqFrame :selpars.freqframe));
391 : //vi::FrequencySelectionUsingFrame channelSelector(selpars.freqframe);
392 22 : for(uInt k=0; k < nSelections; ++k){
393 : //channelSelector.add(Int(freqList(k,0)), lowfreq, topfreq);
394 : //andFreqSelection((mss_p.nelements()-1), Int(freqList(k,0)), lowfreq, topfreq, vi_p ?freqFrame : selpars.freqframe);
395 11 : andFreqSelection((mss_p.nelements()-1), Int(freqList(k,0)), lowfreq, topfreq, selFreqFrame_p);
396 : }
397 : //fselections_p->add(channelSelector);
398 :
399 11 : }
400 1800 : }//End of channel selection
401 :
402 1800 : writeAccess_p=writeAccess_p && !selpars.readonly;
403 1800 : createVisSet(writeAccess_p);
404 :
405 : /////// Remove this when the new vi/vb is able to get the full freq range.
406 1800 : mssFreqSel_p.resize();
407 1800 : mssFreqSel_p = thisSelection.getChanFreqList(NULL,true);
408 :
409 : //// Set the data column on which to operate
410 : // TT: added checks for the requested data column existace
411 : // cout << "Using col : " << selpars.datacolumn << endl;
412 1800 : if( selpars.datacolumn.contains("data") || selpars.datacolumn.contains("obs") )
413 136 : { if( thisms.tableDesc().isColumn("DATA") ) { datacol_p = FTMachine::OBSERVED; }
414 0 : else { os << LogIO::SEVERE <<"DATA column does not exist" << LogIO::EXCEPTION;}
415 : }
416 1664 : else if( selpars.datacolumn.contains("corr") ) { datacol_p = FTMachine::CORRECTED; }
417 0 : else { os << LogIO::WARN << "Invalid data column : " << selpars.datacolumn << ". Using corrected (or observed if corrected doesn't exist)" << LogIO::POST; datacol_p = FTMachine::CORRECTED;}
418 :
419 :
420 1800 : dataSel_p.resize(dataSel_p.nelements()+1, true);
421 1800 : dataSel_p[dataSel_p.nelements()-1]=selpars;
422 : //cerr << "AT THE end of DATASEL " << selpars.toRecord() << endl;
423 :
424 1800 : unlockMSs();
425 1809 : }
426 4 : catch(AipsError &x)
427 : {
428 4 : unlockMSs();
429 4 : throw( AipsError("Error in selectData() : "+x.getMesg()) );
430 4 : }
431 :
432 1800 : return retval;
433 :
434 :
435 :
436 1804 : }
437 2383 : void SynthesisImagerVi2::andChanSelection(const Int msId, const Int spwId, const Int startchan, const Int endchan){
438 :
439 2383 : map<Int, Vector<Int> > spwsel;
440 2383 : auto it=channelSelections_p.find(msId);
441 2383 : if(it !=channelSelections_p.end())
442 594 : spwsel=it->second;
443 2383 : auto hasspw=spwsel.find(spwId);
444 2383 : Vector<Int>chansel(2,-1);
445 2383 : if(hasspw != spwsel.end()){
446 30 : chansel.resize();
447 30 : chansel=hasspw->second;
448 : }
449 2383 : Int nchan=endchan-startchan+1;
450 2383 : if(chansel(1)== -1)
451 2353 : chansel(1)=startchan;
452 2383 : if(chansel(1) >= startchan){
453 2354 : if(nchan > (chansel(1)-startchan+chansel(0))){
454 2353 : chansel(0)=nchan;
455 : }
456 : else{
457 1 : chansel(0)=chansel(1)-startchan+chansel(0);
458 : }
459 2354 : chansel(1)=startchan;
460 : }
461 : else{
462 29 : if((chansel(0) -(startchan - chansel(1)+1)) < nchan){
463 29 : chansel(0)=nchan+(startchan-chansel(1));
464 : }
465 : }
466 2383 : spwsel[spwId]=chansel;
467 2383 : channelSelections_p[msId]=spwsel;
468 : //cerr << "chansel "<< channelSelections_p << endl;
469 :
470 2383 : }
471 2394 : void SynthesisImagerVi2::andFreqSelection(const Int msId, const Int spwId, const Double freqBeg, const Double freqEnd, const MFrequency::Types frame){
472 :
473 :
474 2394 : Int key=msId;
475 :
476 2394 : Bool isDefined=False;
477 2394 : FrequencySelectionUsingFrame frameSel(frame);
478 3363 : for (uInt k =0; k<freqBegs_p.size(); ++k){
479 : //cerr <<freqBegs_p[k].first << " == " << key << " && " << freqSpws_p[k].second<< " == " << spwId << " && " << freqBeg << " < " << freqEnds_p[k].second<< " && " << freqEnd << " > " << freqBegs_p[k].second << endl;
480 969 : if((freqBegs_p[k].first == key || key <0 ) && (freqSpws_p[k].second==spwId || spwId <0) && (freqBeg < freqEnds_p[k].second) && (freqEnd > freqBegs_p[k].second)){
481 12 : isDefined=True;
482 : //cerr << k << " inside freqBegs " << freqBegs_p[k].second << " " << freqBeg << endl;
483 12 : if(freqBegs_p[k].second < freqBeg)
484 0 : freqBegs_p[k].second=freqBeg;
485 12 : if(freqEnds_p[k].second > freqEnd)
486 3 : freqEnds_p[k].second=freqEnd;
487 12 : if(msId < 0) key=freqBegs_p[k].first;
488 : //cerr << "modified " << freqBegs_p[k].second << " " << freqEnds_p[k].second << endl;
489 : }
490 : ///add only those that have the same msid
491 969 : if(freqBegs_p[k].first == key){
492 : //cerr << "added " << k << " freqBegs " << freqBegs_p[k].second << " " << freqEnds_p[k].second << endl;
493 885 : frameSel.add(freqSpws_p[k].second , freqBegs_p[k].second, freqEnds_p[k].second);
494 : }
495 : }
496 2394 : if(!isDefined && msId >=0){
497 : //cerr << "undefined " << key << " freqBegs " << freqBeg << " " << freqEnd << endl;
498 2382 : freqBegs_p.push_back(make_pair(key, freqBeg));
499 2382 : freqEnds_p.push_back(make_pair(key, freqEnd));
500 2382 : freqSpws_p.push_back(make_pair(key, spwId));
501 2382 : frameSel.add(spwId, freqBeg, freqEnd);
502 : }
503 2394 : CountedPtr<vi::FrequencySelections> copyFsels=fselections_p->clone();
504 2394 : uInt nMSs=copyFsels->size() <=msId ? msId+1 : copyFsels->size();
505 : //cerr << "nms " << nMSs << endl;
506 2394 : fselections_p=new FrequencySelections();
507 4872 : for (uInt k=0; k < nMSs ; ++k){
508 2478 : if(k==uInt(key)){
509 2394 : fselections_p->add(frameSel);
510 : //cerr <<"adding framesel " << frameSel.toString() << endl;
511 : }
512 : else{
513 84 : const FrequencySelectionUsingFrame& thissel= static_cast<const FrequencySelectionUsingFrame &> (copyFsels->get(k));
514 : //cerr <<"framesel orig " << thissel.toString() << endl;
515 84 : fselections_p->add(thissel);
516 :
517 : }
518 : }
519 :
520 :
521 :
522 2394 : }
523 :
524 0 : void SynthesisImagerVi2::tuneChunk(const Int gmap) {
525 :
526 0 : CoordinateSystem cs=itsMappers.imageStore(gmap)->getCSys();
527 0 : IPosition imshape=itsMappers.imageStore(gmap)->getShape();
528 : /////For some reason imagestore returns 0 channel image sometimes
529 : ////
530 0 : if(imshape(3) < 1) {
531 0 : return;
532 : }
533 :
534 0 : Double minFreq=SpectralImageUtil::worldFreq(cs, 0.0);
535 0 : Double maxFreq=SpectralImageUtil::worldFreq(cs,imshape(3)-1);
536 :
537 0 : if(maxFreq < minFreq){
538 0 : Double tmp=minFreq;
539 0 : minFreq=maxFreq;
540 0 : maxFreq=tmp;
541 : }
542 :
543 0 : Int spectralIndex=cs.findCoordinate(Coordinate::SPECTRAL);
544 0 : SpectralCoordinate spectralCoord=cs.spectralCoordinate(spectralIndex);
545 0 : maxFreq+=fabs(spectralCoord.increment()(0))/2.0;
546 0 : minFreq-=fabs(spectralCoord.increment()(0))/2.0;
547 0 : if(minFreq < 0.0) minFreq=0.0;
548 0 : MFrequency::Types intype=spectralCoord.frequencySystem(True);
549 :
550 0 : if(!VisBufferUtil::getFreqRangeFromRange(minFreq, maxFreq, intype, minFreq, maxFreq, *vi_p, selFreqFrame_p)){
551 : //Do not retune if conversion did not happen
552 0 : return;
553 : }
554 :
555 0 : maxFreq+=3.0*fabs(spectralCoord.increment()(0))/2.0;
556 0 : minFreq-=3.0*fabs(spectralCoord.increment()(0))/2.0;
557 0 : if(minFreq < 0.0) minFreq=0.0;
558 :
559 0 : auto copyFreqBegs=freqBegs_p;
560 0 : auto copyFreqEnds=freqEnds_p;
561 0 : auto copyFreqSpws= freqSpws_p;
562 : ///////////////TESTOO
563 : //cerr << std::setprecision(12) << "AFTER maxFreq " << maxFreq << " minFreq " << minFreq << endl;
564 : //for (Int k =0 ; k < (fselections_p->size()) ; ++k){
565 : // cerr << k << (fselections_p->get(k)).toString() << endl;
566 : // }
567 : ///////////////////////////////////////
568 : ///TESTOO
569 : // andFreqSelection(-1, -1, minFreq, maxFreq, MFrequency::TOPO);
570 0 : andFreqSelection(-1, -1, minFreq, maxFreq, selFreqFrame_p);
571 :
572 0 : vi_p->setFrequencySelection (*fselections_p);
573 :
574 0 : freqBegs_p=copyFreqBegs;
575 0 : freqEnds_p=copyFreqEnds;
576 0 : freqSpws_p=copyFreqSpws;
577 0 : }
578 :
579 :
580 892 : Bool SynthesisImagerVi2::defineImage(
581 : SynthesisParamsImage& impars,
582 : const SynthesisParamsGrid& gridpars)
583 : {
584 1784 : LogIO os( LogOrigin("SynthesisImagerVi2", "defineImage", WHERE) );
585 :
586 892 : if (mss_p.nelements() == 0)
587 0 : os << "SelectData has to be run before defineImage" << LogIO::EXCEPTION;
588 :
589 892 : CoordinateSystem csys;
590 892 : CountedPtr<refim::FTMachine> ftm, iftm;
591 892 : impars_p = impars;
592 892 : gridpars_p = gridpars;
593 :
594 : try {
595 892 : os << "Define image coordinates for [" << impars.imageName << "] : "
596 892 : << LogIO::POST;
597 892 : os << "Define image coordinates for [" << impars.imageName << "] : " << LogIO::POST;
598 : // cerr << "DEFIM " << gridpars_p.ftmachine << endl;
599 : // cerr << "###### gridpars compute " << gridpars.computePAStep << " " << gridpars_p.computePAStep << endl;
600 893 : csys = impars_p.buildCoordinateSystem( *vi_p, channelSelections_p, mss_p );
601 : //use the location defined for coordinates frame;
602 891 : mLocation_p=impars_p.obslocation;
603 891 : IPosition imshape = impars_p.shp();
604 :
605 :
606 890 : os << "Impars: start " << impars_p.start << LogIO::POST;
607 : os << "Shape: " << imshape
608 890 : << " Spectral: " << csys.spectralCoordinate().referenceValue()
609 1780 : << " at " << csys.spectralCoordinate().referencePixel()
610 1780 : << " with increment " << csys.spectralCoordinate().increment()
611 3560 : << LogIO::POST;
612 :
613 906 : if( (itsMappers.nMappers()==0) ||
614 16 : (impars_p.imsize[0]*impars_p.imsize[1] > itsMaxShape[0]*itsMaxShape[1]))
615 : {
616 876 : itsMaxShape=imshape;
617 876 : itsMaxCoordSys=csys;
618 : }
619 890 : itsNchan = imshape[3];
620 890 : itsCsysRec = impars_p.getcsys();
621 : /*
622 : os << "Define image [" << impars.imageName << "] : nchan : " << impars.nchan
623 : //<< ", freqstart:" << impars.freqStart.getValue() << impars.freqStart.getUnit()
624 : << ", start:" << impars.start
625 : << ", imsize:" << impars.imsize
626 : << ", cellsize: [" << impars.cellsize[0].getValue() << impars.cellsize[0].getUnit()
627 : << " , " << impars.cellsize[1].getValue() << impars.cellsize[1].getUnit()
628 : << LogIO::POST;
629 : */
630 : // phasecenter
631 890 : if (impars_p.phaseCenterFieldId == -1) {
632 : // user-specified
633 338 : phaseCenter_p = impars_p.phaseCenter;
634 552 : } else if (impars_p.phaseCenterFieldId >= 0) {
635 : // FIELD_ID
636 15 : auto const msobj = mss_p[0];
637 15 : MSFieldColumns msfield(msobj->field());
638 15 : phaseCenter_p=msfield.phaseDirMeas(impars_p.phaseCenterFieldId);
639 15 : } else {
640 : // use default FIELD_ID (0)
641 537 : auto const msobj = mss_p[0];
642 537 : MSFieldColumns msfield(msobj->field());
643 537 : phaseCenter_p=msfield.phaseDirMeas(0);
644 537 : }
645 :
646 :
647 906 : if ( (itsMappers.nMappers() == 0) or
648 16 : (impars_p.imsize[0]*impars_p.imsize[1] > itsMaxShape[0]*itsMaxShape[1])
649 : ) {
650 874 : itsMaxShape = imshape;
651 874 : itsMaxCoordSys = csys;
652 : }
653 890 : itsNchan = imshape[3];
654 890 : itsCsysRec = impars_p.getcsys();
655 :
656 : // os << "Define image [" << impars.imageName << "] : nchan : " << impars.nchan
657 : // //<< ", freqstart:" << impars.freqStart.getValue() << impars.freqStart.getUnit()
658 : // << ", start:" << impars.start
659 : // << ", imsize:" << impars.imsize
660 : // << ", cellsize: [" << impars.cellsize[0].getValue() << impars.cellsize[0].getUnit()
661 : // << " , " << impars.cellsize[1].getValue() << impars.cellsize[1].getUnit()
662 : // << LogIO::POST;
663 :
664 : // phasecenter
665 890 : if (impars_p.phaseCenterFieldId == -1) { // user-specified
666 338 : phaseCenter_p = impars_p.phaseCenter;
667 552 : } else if (impars_p.phaseCenterFieldId >= 0) { // FIELD_ID
668 15 : auto const msobj = mss_p[0];
669 15 : MSFieldColumns msfield(msobj->field());
670 15 : phaseCenter_p = msfield.phaseDirMeas(impars_p.phaseCenterFieldId);
671 15 : } else { // use default FIELD_ID (0)
672 537 : auto const msobj = mss_p[0];
673 537 : MSFieldColumns msfield(msobj->field());
674 537 : phaseCenter_p = msfield.phaseDirMeas(0);
675 537 : }
676 890 : }
677 2 : catch (AipsError &x) {
678 : os << "Error in building Coordinate System and Image Shape: "
679 : << x.getMesg()
680 2 : << LogIO::EXCEPTION;
681 2 : }
682 :
683 : try {
684 890 : os << "Set Gridding options for [" << impars_p.imageName << "]"
685 890 : << " with ftmachine: " << gridpars.ftmachine
686 890 : << LogIO::POST;
687 :
688 890 : itsVpTable = gridpars.vpTable;
689 1702 : itsMakeVP = ( gridpars.ftmachine.contains("mosaicft") or
690 1702 : gridpars.ftmachine.contains("awprojectft") ) ?
691 : False : True;
692 :
693 : //cerr << "DEFINEimage " << impars_p.toRecord() << endl;
694 :
695 1780 : createFTMachine(
696 890 : ftm, iftm, gridpars.ftmachine, impars_p.nTaylorTerms, gridpars.mType,
697 890 : gridpars.facets, gridpars.wprojplanes,
698 890 : gridpars.padding, gridpars.useAutoCorr, gridpars.useDoublePrec,
699 890 : gridpars.convFunc,
700 890 : gridpars.aTermOn, gridpars.psTermOn, gridpars.mTermOn,
701 1780 : gridpars.wbAWP, gridpars.cfCache, gridpars.usePointing, gridpars.pointingOffsetSigDev.tovector(),
702 890 : gridpars.doPBCorr, gridpars.conjBeams,
703 890 : gridpars.computePAStep, gridpars.rotatePAStep,
704 890 : gridpars.interpolation, impars_p.freqFrameValid, 1000000000, 16, impars_p.stokes,
705 890 : impars_p.imageName,
706 890 : gridpars.pointingDirCol, gridpars.convertFirst, gridpars.skyPosThreshold,
707 890 : gridpars.convSupport, gridpars.truncateSize, gridpars.gwidth, gridpars.jwidth,
708 890 : gridpars.minWeight, gridpars.clipMinMax, impars_p.pseudoi
709 : );
710 :
711 : }
712 0 : catch (AipsError &x) {
713 0 : os << "Error in setting up FTMachine(): " << x.getMesg() << LogIO::EXCEPTION;
714 0 : }
715 :
716 :
717 : try {
718 1780 : appendToMapperList(
719 1780 : impars_p.imageName, csys, impars_p.shp(),
720 : ftm, iftm,
721 890 : gridpars.distance, gridpars.facets, gridpars.chanchunks, impars_p.overwrite,
722 890 : gridpars.mType, gridpars.padding, impars_p.nTaylorTerms, impars_p.startModel
723 : );
724 :
725 890 : imageDefined_p = true;
726 : }
727 0 : catch(AipsError &x) {
728 0 : os << "Error in adding Mapper: " + x.getMesg() << LogIO::EXCEPTION;
729 0 : }
730 :
731 890 : imparsVec_p.resize(imparsVec_p.nelements()+1, true);
732 890 : imparsVec_p[imparsVec_p.nelements()-1] = impars_p;
733 : ///For now cannot deal with cube and mtmfs in C++ parallel mode
734 890 : if (imparsVec_p[0].deconvolver == "mtmfs") setCubeGridding(False);
735 : // cerr << "DECONV " << imparsVec_p[0].deconvolver
736 : // << " cube gridding " << doingCubeGridding_p << endl;
737 890 : gridparsVec_p.resize(gridparsVec_p.nelements()+1, true);
738 890 : gridparsVec_p[imparsVec_p.nelements()-1] = gridpars_p;
739 : // For now as awproject does not work with the c++ mpi cube gridding
740 : // make sure it works the old way as mfs
741 : // if ( gridparsVec_p[0].ftmachine.contains("awproject") )
742 : // setCubeGridding(False);
743 890 : itsMakeVP =
744 1702 : ( gridparsVec_p[0].ftmachine.contains("mosaicft") or
745 1702 : (gridparsVec_p[0].ftmachine.at(0,3) == "awp")
746 2592 : ) ? False : True;
747 890 : return true;
748 898 : }
749 :
750 842 : Bool SynthesisImagerVi2::defineImage(
751 : CountedPtr<SIImageStore> imstor,
752 : SynthesisParamsImage& impars,
753 : const SynthesisParamsGrid& gridpars)
754 : {
755 842 : gridpars_p=gridpars;
756 842 : Int id = itsMappers.nMappers();
757 842 : CoordinateSystem csys = imstor->getCSys();
758 842 : IPosition imshape = imstor->getShape();
759 842 : Int nx=imshape[0], ny=imshape[1];
760 842 : if ( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]) ) {
761 827 : itsMaxShape=imshape;
762 827 : itsMaxCoordSys=csys;
763 : }
764 842 : mLocation_p=impars.obslocation;
765 : // phasecenter
766 842 : if (impars.phaseCenterFieldId == -1) { // user-specified
767 395 : phaseCenter_p = impars.phaseCenter;
768 : }
769 447 : else if (impars.phaseCenterFieldId >= 0) { // FIELD_ID
770 9 : auto const msobj = mss_p[0];
771 9 : MSFieldColumns msfield(msobj->field());
772 9 : phaseCenter_p = msfield.phaseDirMeas(impars.phaseCenterFieldId);
773 9 : }
774 : else { // use default FIELD_ID (0)
775 438 : auto const msobj = mss_p[0];
776 438 : MSFieldColumns msfield(msobj->field());
777 438 : phaseCenter_p = msfield.phaseDirMeas(0);
778 438 : }
779 :
780 842 : itsVpTable = gridpars.vpTable;
781 712 : itsMakeVP = ( gridpars.ftmachine.contains("mosaicft") or
782 1554 : (gridpars.ftmachine.at(0,3) == "awp")
783 1554 : ) ? False : True;
784 842 : CountedPtr<refim::FTMachine> ftm, iftm;
785 :
786 1684 : createFTMachine(
787 842 : ftm, iftm, gridpars.ftmachine, impars.nTaylorTerms, gridpars.mType,
788 842 : gridpars.facets, gridpars.wprojplanes,
789 842 : gridpars.padding,gridpars.useAutoCorr,gridpars.useDoublePrec,
790 842 : gridpars.convFunc,
791 842 : gridpars.aTermOn,gridpars.psTermOn, gridpars.mTermOn,
792 842 : gridpars.wbAWP,gridpars.cfCache,gridpars.usePointing,
793 842 : gridpars.pointingOffsetSigDev.tovector(),
794 842 : gridpars.doPBCorr,gridpars.conjBeams,
795 842 : gridpars.computePAStep,gridpars.rotatePAStep,
796 842 : gridpars.interpolation, impars.freqFrameValid, 1000000000, 16, impars.stokes,
797 842 : impars.imageName,
798 842 : gridpars.pointingDirCol, gridpars.convertFirst, gridpars.skyPosThreshold,
799 842 : gridpars.convSupport, gridpars.truncateSize, gridpars.gwidth, gridpars.jwidth,
800 842 : gridpars.minWeight, gridpars.clipMinMax, impars.pseudoi);
801 :
802 842 : if (gridpars.facets >1) {
803 : // Make and connect the list.
804 : Block<CountedPtr<SIImageStore> > imstorList =
805 0 : createFacetImageStoreList( imstor, gridpars.facets );
806 0 : for( uInt facet=0; facet<imstorList.nelements(); facet++) {
807 0 : CountedPtr<refim::FTMachine> new_ftm, new_iftm;
808 0 : if (facet == 0) {
809 0 : new_ftm = ftm;
810 0 : new_iftm = iftm;
811 : }
812 : else {
813 0 : new_ftm = ftm->cloneFTM();
814 0 : new_iftm = iftm->cloneFTM();
815 : }
816 0 : itsMappers.addMapper(
817 0 : createSIMapper( gridpars.mType, imstorList[facet], new_ftm, new_iftm)
818 : );
819 0 : }
820 0 : }
821 : else {
822 842 : itsMappers.addMapper(
823 1684 : createSIMapper( gridpars.mType, imstor, ftm, iftm)
824 : );
825 : }
826 842 : impars_p = impars;
827 842 : gridpars_p = gridpars;
828 842 : imageDefined_p = true;
829 :
830 842 : imparsVec_p.resize(imparsVec_p.nelements()+1, true);
831 842 : imparsVec_p[imparsVec_p.nelements()-1] = impars_p;
832 :
833 842 : gridparsVec_p.resize(gridparsVec_p.nelements()+1, true);
834 842 : gridparsVec_p[gridparsVec_p.nelements()-1] = gridpars_p;
835 :
836 842 : return true;
837 842 : }
838 :
839 0 : Bool SynthesisImagerVi2::defineImage(CountedPtr<SIImageStore> imstor,
840 : const String& ftmachine)
841 : {
842 0 : CountedPtr<refim::FTMachine> ftm, iftm;
843 :
844 : // The following call to createFTMachine() uses the
845 : // following defaults
846 : //
847 : // facets=1, wprojplane=1, padding=1.0, useAutocorr=false,
848 : // useDoublePrec=true, gridFunction=String("SF")
849 : //
850 0 : createFTMachine(ftm, iftm, ftmachine);
851 :
852 0 : Int id=itsMappers.nMappers();
853 0 : CoordinateSystem csys =imstor->getCSys();
854 0 : IPosition imshape=imstor->getShape();
855 0 : Int nx=imshape[0], ny=imshape[1];
856 0 : if( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]))
857 : {
858 0 : itsMaxShape=imshape;
859 0 : itsMaxCoordSys=csys;
860 : }
861 :
862 0 : itsMappers.addMapper( createSIMapper( "default", imstor, ftm, iftm ) );
863 0 : imageDefined_p=true;
864 0 : return true;
865 0 : }
866 0 : Bool SynthesisImagerVi2::defineImage(CountedPtr<SIImageStore> imstor,
867 : const Record& ftmRec, const Record& iftmRec)
868 : {
869 0 : CountedPtr<refim::FTMachine> ftm, iftm;
870 0 : ftm=refim::VisModelData::NEW_FT(ftmRec);
871 0 : iftm=refim::VisModelData::NEW_FT(iftmRec);
872 :
873 0 : Int id=itsMappers.nMappers();
874 0 : CoordinateSystem csys =imstor->getCSys();
875 0 : IPosition imshape=imstor->getShape();
876 0 : Int nx=imshape[0], ny=imshape[1];
877 0 : if( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]))
878 : {
879 0 : itsMaxShape=imshape;
880 0 : itsMaxCoordSys=csys;
881 : }
882 :
883 0 : itsMappers.addMapper( createSIMapper( "default", imstor, ftm, iftm, id ) );
884 0 : imageDefined_p=true;
885 0 : return true;
886 0 : }
887 826 : Bool SynthesisImagerVi2::weight(const Record& inrec){
888 826 : String type, rmode, filtertype;
889 826 : Quantity noise, fieldofview,filterbmaj,filterbmin, filterbpa;
890 : Double robust, fracBW;
891 : Int npixels;
892 : Bool multiField, useCubeBriggs;
893 826 : SynthesisUtilMethods::getFromWeightRecord(type, rmode,noise, robust,fieldofview,npixels, multiField, useCubeBriggs,
894 : filtertype, filterbmaj,filterbmin, filterbpa, fracBW, inrec);
895 826 : return weight(type, rmode,noise, robust,fieldofview,npixels, multiField, useCubeBriggs,
896 1652 : filtertype, filterbmaj,filterbmin, filterbpa, fracBW );
897 :
898 :
899 826 : }
900 1545 : Bool SynthesisImagerVi2::weight(const String& type, const String& rmode,
901 : const Quantity& noise, const Double robust,
902 : const Quantity& fieldofview,
903 : const Int npixels, const Bool multiField, const Bool useCubeBriggs,
904 : const String& filtertype, const Quantity& filterbmaj,
905 : const Quantity& filterbmin, const Quantity& filterbpa, Double fracBW)
906 : {
907 3090 : LogIO os(LogOrigin("SynthesisImagerVi2", "weight()", WHERE));
908 :
909 1545 : if(rmode=="bwtaper") //See CAS-13021 for bwtaper algorithm details
910 : {
911 3 : if(fracBW == 0.0)
912 : {
913 1 : Double minFreq = 0.0;
914 1 : Double maxFreq = 0.0;
915 :
916 1 : if(itsMaxShape(3) < 1) {
917 0 : cout << "SynthesisImagerVi2::weight Only one channel in image " << endl;
918 : }
919 : else{
920 1 : minFreq=abs(SpectralImageUtil::worldFreq(itsMaxCoordSys, 0.0));
921 1 : maxFreq=abs(SpectralImageUtil::worldFreq(itsMaxCoordSys,itsMaxShape(3)-1));
922 :
923 1 : if(maxFreq < minFreq){
924 0 : Double tmp=minFreq;
925 0 : minFreq=maxFreq;
926 0 : maxFreq=tmp;
927 : }
928 :
929 1 : if((maxFreq != 0.0) || (minFreq != 0.0)) fracBW = 2*(maxFreq - minFreq)/(maxFreq + minFreq);
930 :
931 1 : os << LogIO::NORMAL << " Fractional bandwidth used by briggsbwtaper " << fracBW << endl; //<< LogIO::POST;
932 :
933 : }
934 : }
935 : }
936 :
937 3090 : weightParams_p=SynthesisUtilMethods::fillWeightRecord(type, rmode,noise, robust,fieldofview,
938 1545 : npixels, multiField, useCubeBriggs,filtertype, filterbmaj,filterbmin, filterbpa, fracBW);
939 :
940 : try {
941 : //Int nx=itsMaxShape[0];
942 : //Int ny=itsMaxShape[1];
943 :
944 :
945 : ///////////////////////
946 3090 : Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
947 3090 : Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
948 : os << LogIO::NORMAL // Loglevel INFO
949 1545 : << "Set imaging weights : " ; //<< LogIO::POST;
950 :
951 1545 : if (type=="natural") {
952 : os << LogIO::NORMAL // Loglevel INFO
953 1484 : << "Natural weighting" << LogIO::POST;
954 1484 : imwgt_p=VisImagingWeight("natural");
955 : }
956 61 : else if (type=="radial") {
957 1 : os << "Radial weighting" << LogIO::POST;
958 1 : imwgt_p=VisImagingWeight("radial");
959 : }
960 : else{
961 60 : vi_p->originChunks();
962 60 : vi_p->origin();
963 60 : if(!imageDefined_p)
964 0 : throw(AipsError("Need to define image"));
965 60 : Int nx=itsMaxShape[0];
966 60 : Int ny=itsMaxShape[1];
967 120 : Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
968 120 : Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
969 60 : if(type=="superuniform"){
970 1 : if(!imageDefined_p) throw(AipsError("Please define image first"));
971 1 : Int actualNpix=npixels;
972 1 : if(actualNpix <=0)
973 1 : actualNpix=3;
974 : os << LogIO::NORMAL // Loglevel INFO
975 : << "SuperUniform weighting over a square cell spanning ["
976 : << -actualNpix
977 1 : << ", " << actualNpix << "] in the uv plane" << LogIO::POST;
978 2 : imwgt_p=VisImagingWeight(*vi_p, rmode, noise, robust, nx,
979 : ny, cellx, celly, actualNpix,
980 1 : actualNpix, multiField);
981 : }
982 59 : else if ((type=="robust")||(type=="uniform")||(type=="briggs")) {
983 59 : if(!imageDefined_p) throw(AipsError("Please define image first"));
984 59 : Quantity actualFieldOfView_x(fieldofview), actualFieldOfView_y(fieldofview) ;
985 59 : Int actualNPixels_x(npixels),actualNPixels_y(npixels) ;
986 59 : String wtype;
987 59 : if(type=="briggs") {
988 46 : wtype = "Briggs";
989 : }
990 : else {
991 13 : wtype = "Uniform";
992 : }
993 59 : if(actualFieldOfView_x.get().getValue()==0.0&&actualNPixels_x==0) {
994 59 : actualNPixels_x=nx;
995 59 : actualFieldOfView_x=Quantity(actualNPixels_x*cellx.get("rad").getValue(),"rad");
996 59 : actualNPixels_y=ny;
997 59 : actualFieldOfView_y=Quantity(actualNPixels_y*celly.get("rad").getValue(),"rad");
998 : os << LogIO::NORMAL // Loglevel INFO
999 : << wtype
1000 : << " weighting: sidelobes will be suppressed over full image"
1001 59 : << LogIO::POST;
1002 : }
1003 0 : else if(actualFieldOfView_x.get().getValue()>0.0&&actualNPixels_x==0) {
1004 0 : actualNPixels_x=nx;
1005 0 : actualNPixels_y=ny;
1006 : os << LogIO::NORMAL // Loglevel INFO
1007 : << wtype
1008 : << " weighting: sidelobes will be suppressed over specified field of view: "
1009 0 : << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by "
1010 0 : << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
1011 : }
1012 0 : else if(actualFieldOfView_x.get().getValue()==0.0&&actualNPixels_x>0) {
1013 0 : actualFieldOfView_x=Quantity(actualNPixels_x*cellx.get("rad").getValue(),"rad");
1014 0 : actualFieldOfView_y=Quantity(actualNPixels_y*celly.get("rad").getValue(),"rad");
1015 : os << LogIO::NORMAL // Loglevel INFO
1016 : << wtype
1017 : << " weighting: sidelobes will be suppressed over full image field of view: "
1018 0 : << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by "
1019 0 : << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
1020 : }
1021 : else {
1022 : os << LogIO::NORMAL // Loglevel INFO
1023 : << wtype
1024 : << " weighting: sidelobes will be suppressed over specified field of view: "
1025 0 : << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by "
1026 0 : << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
1027 : }
1028 : os << LogIO::DEBUG1
1029 : << "Weighting used " << actualNPixels_x << " by " << actualNPixels_y << " uv pixels."
1030 59 : << LogIO::POST;
1031 59 : Quantity actualCellSize_x(actualFieldOfView_x.get("rad").getValue()/actualNPixels_x, "rad");
1032 59 : Quantity actualCellSize_y(actualFieldOfView_y.get("rad").getValue()/actualNPixels_y, "rad");
1033 :
1034 :
1035 : // cerr << "rmode " << rmode << " noise " << noise << " robust " << robust << " npixels " << actualNPixels << " cellsize " << actualCellSize << " multifield " << multiField << endl;
1036 : // Timer timer;
1037 : //timer.mark();
1038 : //Construct imwgt_p with old vi for now if old vi is in use as constructing with vi2 is slower
1039 : //Determine if any image is cube
1040 59 : if(useCubeBriggs){
1041 82 : String outstr=String("Doing spectral cube Briggs weighting formula -- " + rmode + (rmode=="abs" ? " with estimated noise "+ String::toString(noise.getValue())+noise.getUnit() : ""));
1042 41 : os << outstr << LogIO::POST;
1043 : //VisImagingWeight nat("natural");
1044 : //vi_p->useImagingWeight(nat);
1045 41 : if(rmode=="abs" && robust==0.0 && noise.getValue()==0.0)
1046 0 : throw(AipsError("Absolute Briggs formula does not allow for robust 0 and estimated noise per visibility 0"));
1047 :
1048 41 : CountedPtr<refim::BriggsCubeWeightor> bwgt=new refim::BriggsCubeWeightor(wtype=="Uniform" ? "none" : rmode, noise, robust,fracBW, npixels, multiField);
1049 82 : for (Int k=0; k < itsMappers.nMappers(); ++k){
1050 41 : itsMappers.getFTM2(k)->setBriggsCubeWeight(bwgt);
1051 : }
1052 :
1053 :
1054 41 : }
1055 : else
1056 : {
1057 36 : imwgt_p=VisImagingWeight(*vi_p, wtype=="Uniform" ? "none" : rmode, noise, robust,
1058 : actualNPixels_x, actualNPixels_y, actualCellSize_x,
1059 18 : actualCellSize_y, 0, 0, multiField);
1060 : }
1061 :
1062 : /*
1063 : if(rvi_p !=NULL){
1064 : imwgt_p=VisImagingWeight(*rvi_p, rmode, noise, robust,
1065 : actualNPixels, actualNPixels, actualCellSize,
1066 : actualCellSize, 0, 0, multiField);
1067 : }
1068 : else{
1069 : ////This is slower by orders of magnitude as of 2014/06/25
1070 : imwgt_p=VisImagingWeight(*vi_p, rmode, noise, robust,
1071 : actualNPixels, actualNPixels, actualCellSize,
1072 : actualCellSize, 0, 0, multiField);
1073 : }
1074 : */
1075 : //timer.show("After making visweight ");
1076 :
1077 59 : }
1078 : else {
1079 : //this->unlock();
1080 : os << LogIO::SEVERE << "Unknown weighting " << type
1081 0 : << LogIO::EXCEPTION;
1082 0 : return false;
1083 : }
1084 60 : }
1085 :
1086 : //// UV-Tapering
1087 : //cout << "Taper type : " << filtertype << " : " << (filtertype=="gaussian") << endl;
1088 1545 : if( filtertype == "gaussian" ) {
1089 : // os << "Setting uv-taper" << LogIO::POST;
1090 32 : imwgt_p.setFilter( filtertype, filterbmaj, filterbmin, filterbpa );
1091 : }
1092 1545 : vi_p->useImagingWeight(imwgt_p);
1093 : ///////////////////////////////
1094 :
1095 1545 : SynthesisUtilMethods::getResource("Set Weighting");
1096 :
1097 : /// return true;
1098 :
1099 1545 : }
1100 0 : catch(AipsError &x)
1101 : {
1102 0 : throw( AipsError("Error in Weighting : "+x.getMesg()) );
1103 0 : }
1104 :
1105 1545 : return true;
1106 1545 : }
1107 :
1108 890 : void SynthesisImagerVi2::appendToMapperList(String imagename,
1109 : CoordinateSystem& csys,
1110 : IPosition imshape,
1111 : CountedPtr<refim::FTMachine>& ftm,
1112 : CountedPtr<refim::FTMachine>& iftm,
1113 : Quantity distance,
1114 : Int facets,
1115 : Int chanchunks,
1116 : const Bool overwrite,
1117 : String mappertype,
1118 : Float padding,
1119 : uInt ntaylorterms,
1120 : const Vector<String> &startmodel)
1121 : {
1122 1780 : LogIO log_l(LogOrigin("SynthesisImagerVi2", "appendToMapperList(ftm)"));
1123 : //---------------------------------------------
1124 : // Some checks..
1125 :
1126 890 : if(facets > 1 && itsMappers.nMappers() > 0)
1127 0 : log_l << "Facetted image has to be the first of multifields" << LogIO::EXCEPTION;
1128 :
1129 890 : TcleanProcessingInfo procInfo;
1130 890 : CompositeNumber cn(uInt(imshape[0] * 2));
1131 : // heuristic factors multiplied to imshape based on gridder
1132 890 : size_t fudge_factor = 15;
1133 890 : if (ftm->name()=="MosaicFTNew") {
1134 54 : fudge_factor = 20;
1135 : }
1136 836 : else if (ftm->name()=="GridFT") {
1137 538 : fudge_factor = 9;
1138 : }
1139 890 : std::tie(procInfo, std::ignore, std::ignore) =
1140 1780 : nSubCubeFitInMemory(fudge_factor, imshape, padding);
1141 :
1142 : // chanchunks auto-calculation block, for now still here for awproject (CAS-12204)
1143 890 : if(chanchunks<1)
1144 : {
1145 0 : log_l << "Automatically calculated chanchunks";
1146 0 : log_l << " using imshape : " << imshape << LogIO::POST;
1147 :
1148 : // Do calculation here.
1149 : // This runs once per image field (for multi-field imaging)
1150 : // This runs once per cube partition, and will see only its own partition's shape
1151 : //chanchunks=1;
1152 :
1153 0 : chanchunks = procInfo.chnchnks;
1154 :
1155 : /*log_l << "Required memory " << required_mem / nlocal_procs / 1024. / 1024. / 1024.
1156 : << "\nAvailable memory " << memory_avail / 1024. / 1024 / 1024.
1157 : << " (rc: memory fraction " << usr_memfrac << "% rc memory " << usr_mem / 1024.
1158 : << ")\n" << nlocal_procs << " other processes on node\n"
1159 : << "Setting chanchunks to " << chanchunks << LogIO::POST;
1160 : */
1161 : }
1162 : //record this in gridpars_p
1163 890 : gridpars_p.chanchunks=chanchunks;
1164 890 : if( imshape.nelements()==4 && imshape[3]<chanchunks )
1165 : {
1166 0 : log_l << LogIO::WARN << "An image with " << imshape[3] << " channel(s) cannot be divided into " << chanchunks << " chunks. Please set chanchunks=1 or choose chanchunks<nchan." << LogIO::EXCEPTION;
1167 : }
1168 :
1169 890 : if(chanchunks > 1 && itsMappers.nMappers() > 0)
1170 0 : log_l << "Channel chunking is currently not supported with multi(outlier)-fields. Please submit a feature request if needed." << LogIO::EXCEPTION;
1171 :
1172 890 : if(chanchunks > 1) itsDataLoopPerMapper=true;
1173 :
1174 890 : AlwaysAssert( ( ( ! (ftm->name()=="MosaicFTNew" && mappertype=="imagemosaic") ) &&
1175 : ( ! (ftm->name()=="AWProjectWBFT" && mappertype=="imagemosaic") )) ,
1176 : AipsError );
1177 : //---------------------------------------------
1178 :
1179 : // Create the ImageStore object
1180 890 : CountedPtr<SIImageStore> imstor;
1181 890 : MSColumns msc(*(mss_p[0]));
1182 890 : imstor = createIMStore(imagename, csys, imshape, overwrite,msc, mappertype, ntaylorterms, distance, procInfo, facets, iftm->useWeightImage(), startmodel );
1183 :
1184 : // Create the Mappers
1185 890 : if( facets<2 && chanchunks<2) // One facet. Just add the above imagestore to the mapper list.
1186 : {
1187 886 : itsMappers.addMapper( createSIMapper( mappertype, imstor, ftm, iftm, ntaylorterms) );
1188 : }
1189 : else // This field is facetted. Make a list of reference imstores, and add all to the mapper list.
1190 : {
1191 :
1192 4 : if ( facets>1 && chanchunks==1 )
1193 : {
1194 : // Make and connect the list.
1195 4 : Block<CountedPtr<SIImageStore> > imstorList = createFacetImageStoreList( imstor, facets );
1196 44 : for( uInt facet=0; facet<imstorList.nelements(); facet++)
1197 : {
1198 40 : CountedPtr<refim::FTMachine> new_ftm, new_iftm;
1199 40 : if(facet==0){ new_ftm = ftm; new_iftm = iftm; }
1200 36 : else{ new_ftm=ftm->cloneFTM(); new_iftm=iftm->cloneFTM(); }
1201 : // imstorList[facet]->setDataPolFrame(imstor->getDataPolFrame());
1202 40 : itsMappers.addMapper(createSIMapper( mappertype, imstorList[facet], new_ftm, new_iftm, ntaylorterms));
1203 40 : }
1204 4 : }// facets
1205 0 : else if ( facets==1 && chanchunks>1 )
1206 : {
1207 : // Make and connect the list.
1208 0 : Block<CountedPtr<SIImageStore> > imstorList = createChanChunkImageStoreList( imstor, chanchunks );
1209 0 : for( uInt chunk=0; chunk<imstorList.nelements(); chunk++)
1210 : {
1211 :
1212 0 : CountedPtr<refim::FTMachine> new_ftm, new_iftm;
1213 0 : if(chunk==0){
1214 0 : new_ftm = ftm;
1215 0 : new_iftm = iftm; }
1216 : else{
1217 0 : new_ftm=ftm->cloneFTM();
1218 0 : new_iftm=iftm->cloneFTM(); }
1219 0 : imstorList[chunk]->setDataPolFrame(imstor->getDataPolFrame());
1220 0 : itsMappers.addMapper(createSIMapper( mappertype, imstorList[chunk], new_ftm, new_iftm, ntaylorterms));
1221 0 : }
1222 0 : }// chanchunks
1223 : else
1224 : {
1225 0 : throw( AipsError("Error in requesting "+String::toString(facets)+" facets on a side with " + String::toString(chanchunks) + " channel chunks. Support for faceting along with channel chunking is not yet available. Please submit a feature-request if you need multiple facets as well as chanchunks. ") );
1226 : }
1227 :
1228 : }// facets or chunks
1229 :
1230 890 : }
1231 :
1232 : /////////////////////////
1233 : /**
1234 : * Calculations of memory required / available -> nchunks .
1235 : *
1236 : * Returns a tuple with a TcleanProcessingInfo, vector of start channels per subchunk,
1237 : * vector of end channels.
1238 : */
1239 1718 : std::tuple<TcleanProcessingInfo, Vector<Int>, Vector<Int> > SynthesisImagerVi2::nSubCubeFitInMemory(const Int fudge_factor, const IPosition& imshape, const Float padding){
1240 3436 : LogIO log_l(LogOrigin("SynthesisImagerVi2", "nSubCubeFitInMemory"));
1241 :
1242 1718 : Double required_mem = fudge_factor * sizeof(Float);
1243 1718 : int nsubcube=1;
1244 1718 : CompositeNumber cn(uInt(imshape[0] * 2));
1245 8590 : for (size_t i = 0; i < imshape.nelements(); i++) {
1246 : // gridding pads image and increases to composite number
1247 6872 : if (i < 2) {
1248 3436 : required_mem *= cn.nextLargerEven(Int(padding*Float(imshape[i])-0.5));
1249 : }
1250 : else {
1251 3436 : required_mem *= imshape[i];
1252 : }
1253 : }
1254 :
1255 : // get number of tclean processes running on the same machine
1256 1718 : size_t nlocal_procs = 1;
1257 1718 : if (getenv("OMPI_COMM_WORLD_LOCAL_SIZE")) {
1258 0 : std::stringstream ss(getenv("OMPI_COMM_WORLD_LOCAL_SIZE"));
1259 0 : ss >> nlocal_procs;
1260 0 : }
1261 : //cerr << "NUM_PROC " << nlocal_procs << endl;
1262 : // assumes all processes need the same amount of memory
1263 1718 : required_mem *= nlocal_procs;
1264 : Double usr_memfrac, usr_mem;
1265 1718 : AipsrcValue<Double>::find(usr_memfrac, "system.resources.memfrac", 80.);
1266 1718 : AipsrcValue<Double>::find(usr_mem, "system.resources.memory", -1.);
1267 : Double memory_avail;
1268 1718 : if (usr_mem > 0.) {
1269 0 : memory_avail = usr_mem * 1024. * 1024.;
1270 : }
1271 : else {
1272 1718 : memory_avail = Double(HostInfo::memoryFree()) * (usr_memfrac / 100.) * 1024.;
1273 : }
1274 : // compute required chanchunks to fit into the available memory
1275 1718 : nsubcube = (int)std::ceil((Double)required_mem / memory_avail);
1276 1718 : Int nworkers= applicator.numProcs() < 2 ? 1 : applicator.numProcs()-1;
1277 1718 : if((nsubcube/nworkers) >1 && nworkers !=1){
1278 0 : nsubcube=(Int(std::floor(Float(nsubcube)/Float(nworkers)))+1)*nworkers;
1279 :
1280 : }
1281 1718 : if (imshape.nelements() == 4 && imshape[3] < nsubcube) {
1282 0 : nsubcube = imshape[3];
1283 : // TODO make chanchunks a divisor of nchannels?
1284 : }
1285 1718 : nsubcube = nsubcube < 1 ? 1 : nsubcube;
1286 1718 : if( (imshape[3] >= nworkers) && (nsubcube < nworkers)){
1287 0 : nsubcube=nworkers;
1288 : ///TESTOO
1289 : //if(imshape[3] > 2*nworkers)
1290 : // nsubcube=2*nworkers;
1291 :
1292 : }
1293 1718 : else if(imshape[3] < (applicator.numProcs()-1)){
1294 0 : nsubcube=imshape[3];
1295 : }
1296 1718 : Int chunksize=imshape[3]/nsubcube;
1297 1718 : Int rem=imshape[3] % nsubcube;
1298 : //case of nchan < numprocs
1299 1718 : if(chunksize==0 && rem > 0){
1300 0 : nsubcube=rem;
1301 0 : chunksize=1;
1302 0 : rem=0;
1303 : }
1304 : ///Avoid an extra chunk with 1 channel as it cause bumps in linear interpolation
1305 : ///See CAS-12625
1306 : /*
1307 : while((rem==1) && (chunksize >1)){
1308 : nsubcube +=1;
1309 : chunksize=imshape[3]/nsubcube;
1310 : rem=imshape[3] % nsubcube;
1311 : }
1312 : if(rem !=0) ++nsubcube;
1313 : . */
1314 1718 : Vector<Int> start(nsubcube,0);
1315 1718 : Vector<Int> end(nsubcube,chunksize-1);
1316 1718 : if(rem >0){
1317 0 : end(0)+=1;
1318 0 : --rem;
1319 : }
1320 1718 : for (Int k=1; k < nsubcube; ++k){
1321 0 : start(k)=end(k-1)+1;
1322 : // end(k)=((k !=nsubcube-1) || rem==0)? (start(k)+chunksize-1) : (start(k)+rem-1);
1323 0 : end(k)=(start(k)+chunksize-1);
1324 0 : if(rem > 0){
1325 0 : end(k)+=1;
1326 0 : --rem;
1327 : }
1328 : }
1329 :
1330 : // print mem related info to log
1331 1718 : const float toGB = 1024.0 * 1024.0 * 1024.0;
1332 1718 : std::ostringstream usr_mem_msg;
1333 1718 : if (usr_mem > 0.) {
1334 0 : usr_mem_msg << usr_mem / 1024.;
1335 : } else {
1336 1718 : usr_mem_msg << "-";
1337 : }
1338 1718 : std::ostringstream oss;
1339 1718 : oss << setprecision(4);
1340 1718 : oss << "Required memory: " << required_mem / toGB
1341 1718 : << " GB. Available mem.: " << memory_avail / toGB
1342 1718 : << " GB (rc, mem. fraction: " << usr_memfrac
1343 1718 : << "%, memory: " << usr_mem_msg.str()
1344 1718 : << ") => Subcubes: " << nsubcube
1345 1718 : << ". Processes on node: " << nlocal_procs << ".\n";
1346 1718 : log_l << oss.str() << LogIO::POST;
1347 :
1348 1718 : TcleanProcessingInfo procInfo;
1349 1718 : procInfo.mpiprocs = nlocal_procs;
1350 1718 : procInfo.chnchnks = nsubcube;
1351 1718 : procInfo.memavail = memory_avail / toGB;
1352 1718 : procInfo.memreq = required_mem / toGB;
1353 3436 : return make_tuple(procInfo, start, end);
1354 1718 : }
1355 :
1356 827 : void SynthesisImagerVi2::runMajorCycleCube( const Bool dopsf,
1357 : const Record inpcontrol) {
1358 1654 : LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycleCube",WHERE) );
1359 827 : if(dopsf){
1360 285 : runCubeGridding(True);
1361 : ///Store the beamsets in ImageInfo
1362 571 : for(Int k=0; k < itsMappers.nMappers(); ++k){
1363 :
1364 288 : (itsMappers.imageStore(k))->getBeamSet();
1365 : }
1366 : }
1367 : else
1368 543 : runCubeGridding(False, inpcontrol);
1369 :
1370 :
1371 :
1372 827 : }
1373 2117 : void SynthesisImagerVi2::runMajorCycle(const Bool dopsf,
1374 : const Bool savemodel)
1375 : {
1376 4234 : LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycle",WHERE) );
1377 :
1378 : // cout << "Savemodel : " << savemodel << " readonly : " << readOnly_p << " usescratch : " << useScratch_p << endl;
1379 2117 : Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
1380 2117 : Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;
1381 :
1382 2117 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
1383 2117 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
1384 :
1385 2117 : SynthesisUtilMethods::getResource("Start Major Cycle");
1386 :
1387 2117 : itsMappers.checkOverlappingModels("blank");
1388 :
1389 : {
1390 2117 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
1391 2117 : vi_p->originChunks();
1392 2117 : vi_p->origin();
1393 2117 : Double numcoh=0;
1394 4306 : for (uInt k=0; k< mss_p.nelements(); ++k)
1395 2189 : numcoh+=Double(mss_p[k]->nrow());
1396 : ProgressMeter pm(1.0, numcoh,
1397 4234 : dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
1398 2117 : rownr_t cohDone=0;
1399 :
1400 :
1401 2117 : if(!dopsf)itsMappers.initializeDegrid(*vb);
1402 2117 : itsMappers.initializeGrid(*vi_p,dopsf);
1403 2116 : SynthesisUtilMethods::getResource("After initGrid for all mappers");
1404 : ////Under some peculiar selection criterion and low channel ms vb2 seems to return more channels than in spw
1405 : {
1406 2116 : vi_p->originChunks();
1407 2116 : vi_p->origin();
1408 2116 : Int nchannow=vb->nChannels();
1409 2116 : Int spwnow=vb->spectralWindows()[0];
1410 2116 : Int nchaninms=MSColumns(vb->ms()).spectralWindow().numChan()(spwnow);
1411 : //cerr << "chans " << nchaninms << " " << nchannow << endl;
1412 :
1413 2116 : if (nchaninms < nchannow){
1414 0 : cerr << "NCHANS ms" << nchaninms << " now " << nchannow << " spw " << spwnow << " " << vb->spectralWindows() << endl;
1415 0 : throw(AipsError("A nasty Visbuffer2 error occured...wait for CNGI"));
1416 : }
1417 : }
1418 : //////
1419 6164 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
1420 : {
1421 :
1422 491821 : for (vi_p->origin(); vi_p->more(); vi_p->next())
1423 : {
1424 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
1425 : // cerr << "nRows "<< vb->nRow() << " " << max(vb->visCube()) << endl;
1426 487773 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
1427 : {
1428 487773 : if(!dopsf) {
1429 313921 : { Cube<Complex> mod(vb->nCorrelations(), vb->nChannels(), vb->nRows(), Complex(0.0));
1430 313921 : vb->setVisCubeModel(mod);
1431 313921 : }
1432 313921 : itsMappers.degrid(*vb, savevirtualmodel );
1433 313921 : if(savemodelcolumn && writeAccess_p ){
1434 9648 : const_cast<MeasurementSet& >((vi_p->ms())).lock(true);
1435 9648 : vi_p->writeVisModel(vb->visCubeModel());
1436 9648 : const_cast<MeasurementSet& >((vi_p->ms())).unlock();
1437 :
1438 : //static_cast<VisibilityIteratorImpl2 *> (vi_p->getImpl())->writeVisModel(vb->visCubeModel());
1439 :
1440 : // Cube<Complex> tt=vb->visCubeModel();
1441 : // tt = 20.0;
1442 : // cout << "Vis:" << tt << endl;
1443 : // static_cast<VisibilityIteratorImpl2 *> (vi_p->getImpl())->writeVisModel(tt);
1444 : }
1445 : }
1446 487773 : itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)datacol_p);
1447 :
1448 487772 : cohDone += vb->nRows();
1449 487772 : pm.update(Double(cohDone));
1450 : }
1451 : }
1452 : }
1453 :
1454 : // cerr << "VI2 data: " << cohDone << endl;
1455 : // exit(0);
1456 : //cerr << "IN SYNTHE_IMA" << endl;
1457 : //VisModelData::listModel(rvi_p->getMeasurementSet());
1458 2115 : SynthesisUtilMethods::getResource("Before finalize for all mappers");
1459 2115 : if(!dopsf) itsMappers.finalizeDegrid(*vb);
1460 2115 : itsMappers.finalizeGrid(*vb, dopsf);
1461 :
1462 2117 : }
1463 :
1464 2115 : itsMappers.checkOverlappingModels("restore");
1465 :
1466 2115 : unlockMSs();
1467 :
1468 2115 : SynthesisUtilMethods::getResource("End Major Cycle");
1469 :
1470 2117 : }// end runMajorCycle
1471 :
1472 :
1473 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1474 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1475 :
1476 : /// The mapper loop is outside the data iterator loop.
1477 : /// This is for cases where the image size is large compared to the RAM and
1478 : /// where data I/O is the relatively minor cost.
1479 0 : void SynthesisImagerVi2::runMajorCycle2(const Bool dopsf,
1480 : const Bool savemodel)
1481 : {
1482 0 : LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycle2",WHERE) );
1483 :
1484 : // cout << "Savemodel : " << savemodel << " readonly : " << readOnly_p << " usescratch : " << useScratch_p << endl;
1485 :
1486 0 : Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
1487 0 : Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;
1488 :
1489 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
1490 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
1491 :
1492 0 : itsMappers.checkOverlappingModels("blank");
1493 :
1494 0 : Bool resetModel=False;
1495 0 : if( savemodelcolumn && writeAccess_p)
1496 : {
1497 0 : resetModel=True;
1498 0 : os << "Iterating through the model column to reset it to zero" << LogIO::POST;
1499 0 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
1500 0 : vi_p->originChunks();
1501 0 : vi_p->origin();
1502 0 : Double numcoh=0;
1503 0 : for (uInt k=0; k< mss_p.nelements(); ++k)
1504 0 : numcoh+=Double(mss_p[k]->nrow());
1505 : ProgressMeter pm(1.0, numcoh,
1506 0 : dopsf?"Seting model column to zero":"pre-Major Cycle", "","","",True);
1507 0 : rownr_t cohDone=0;
1508 0 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
1509 : {
1510 :
1511 0 : for (vi_p->origin(); vi_p->more(); vi_p->next())
1512 : {
1513 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
1514 : {
1515 0 : { Cube<Complex> mod(vb->nCorrelations(), vb->nChannels(), vb->nRows(), Complex(0.0));
1516 0 : vb->setVisCubeModel(mod);
1517 0 : }
1518 0 : vi_p->writeVisModel(vb->visCubeModel());
1519 :
1520 : }
1521 0 : cohDone += vb->nRows();;
1522 0 : pm.update(Double(cohDone));
1523 : }
1524 : }
1525 0 : }// setting model to zero
1526 :
1527 : //Need to inialize the the forward ft machine to save the virtual model on first pass of each ms.
1528 0 : if(!dopsf && savevirtualmodel){
1529 0 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
1530 0 : vi_p->originChunks();
1531 0 : vi_p->origin();
1532 0 : itsMappers.initializeDegrid(*vb, -1);
1533 : }
1534 :
1535 0 : for(Int gmap=0;gmap<itsMappers.nMappers();gmap++)
1536 : {
1537 0 : os << "Running major cycle for chunk : " << gmap << LogIO::POST;
1538 :
1539 0 : SynthesisUtilMethods::getResource("Start Major Cycle for mapper"+String::toString(gmap));
1540 0 : CountedPtr<vi::FrequencySelections> copyFsels=fselections_p->clone();
1541 : ///CAS-12132 create a new visiter for each chunk
1542 0 : createVisSet(writeAccess_p);
1543 : ////////////////////////
1544 0 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
1545 : /// Careful where tunechunk
1546 0 : tuneChunk(gmap);
1547 :
1548 0 : vi_p->originChunks();
1549 0 : vi_p->origin();
1550 :
1551 0 : Double numcoh=0;
1552 0 : for (uInt k=0; k< mss_p.nelements(); ++k)
1553 0 : numcoh+=Double(mss_p[k]->nrow());
1554 :
1555 :
1556 : ProgressMeter pm(1.0, numcoh,
1557 0 : dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
1558 0 : rownr_t cohDone=0;
1559 :
1560 :
1561 0 : itsMappers.getFTM2(gmap, False)->reset();
1562 0 : itsMappers.getFTM2(gmap, True)->reset();
1563 :
1564 0 : if(!dopsf){
1565 0 : itsMappers.initializeDegrid(*vb, gmap);
1566 : //itsMappers.getMapper(gmap)->initializeDegrid(*vb);
1567 : }
1568 0 : itsMappers.initializeGrid(*vi_p,dopsf, gmap);
1569 : //itsMappers.getMapper(gmap)->initializeGrid(*vb,dopsf);
1570 :
1571 0 : SynthesisUtilMethods::getResource("After initialize for mapper"+String::toString(gmap));
1572 0 : Int iterNum=0;
1573 :
1574 :
1575 0 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
1576 : {
1577 :
1578 0 : for (vi_p->origin(); vi_p->more(); vi_p->next())
1579 : {
1580 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
1581 : // cerr << "nRows "<< vb->nRow() << " " << max(vb->visCube()) << endl;
1582 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
1583 : {
1584 :
1585 0 : if(!dopsf) {
1586 0 : if(resetModel==False)
1587 : {
1588 0 : Cube<Complex> mod(vb->nCorrelations(), vb->nChannels(), vb->nRows(), Complex(0.0));
1589 0 : vb->setVisCubeModel(mod);
1590 0 : }
1591 0 : itsMappers.degrid(*vb, savevirtualmodel, gmap );
1592 : //itsMappers.getMapper(gmap)->degrid(*vb); //, savevirtualmodel );
1593 0 : if(savemodelcolumn && writeAccess_p ){
1594 0 : vi_p->writeVisModel(vb->visCubeModel());
1595 : //vi_p->writeBackChanges(vb);
1596 : // static_cast<VisibilityIteratorImpl2 *> (vi_p->getImpl())->writeVisModel(vb->visCubeModel());
1597 : }
1598 :
1599 : }
1600 0 : itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)(datacol_p), gmap);
1601 : //itsMappers.getMapper(gmap)->grid(*vb, dopsf, datacol_p);
1602 0 : cohDone += vb->nRows();
1603 0 : ++iterNum;
1604 0 : pm.update(Double(cohDone));
1605 : }
1606 : }
1607 : }
1608 : //cerr << "IN SYNTHE_IMA" << endl;
1609 : //VisModelData::listModel(rvi_p->getMeasurementSet());
1610 :
1611 0 : SynthesisUtilMethods::getResource("Before finalize for mapper"+String::toString(gmap));
1612 :
1613 0 : if(!dopsf)
1614 : {
1615 0 : itsMappers.finalizeDegrid(*vb,gmap);
1616 : //itsMappers.getMapper(gmap)->finalizeDegrid();
1617 : }
1618 0 : itsMappers.finalizeGrid(*vb, dopsf,gmap);
1619 : //itsMappers.getMapper(gmap)->finalizeGrid(*vb, dopsf);
1620 :
1621 : // itsMappers.getMapper(gmap)->releaseImageLocks();
1622 0 : itsMappers.getMapper(gmap)->imageStore()->releaseComplexGrids();
1623 :
1624 0 : SynthesisUtilMethods::getResource("End Major Cycle for mapper"+String::toString(gmap));
1625 0 : fselections_p=copyFsels;
1626 0 : }// end of mapper loop
1627 : ///CAS-12132 create a new visiter for each chunk
1628 0 : createVisSet(writeAccess_p);
1629 : ////////////////////////
1630 : //////vi_p->setFrequencySelection(*fselections_p);
1631 :
1632 0 : itsMappers.checkOverlappingModels("restore");
1633 :
1634 0 : unlockMSs();
1635 :
1636 0 : SynthesisUtilMethods::getResource("End Major Cycle");
1637 :
1638 0 : }// end runMajorCycle2
1639 :
1640 : //////////////////////////////
1641 :
1642 : ////////////////////////////////////////////////////////////////////////////////////////////////
1643 :
1644 827 : bool SynthesisImagerVi2::runCubeGridding(Bool dopsf, const Record inpcontrol){
1645 :
1646 1654 : LogIO logger(LogOrigin("SynthesisImagerVi2", "runCubeGridding", WHERE));
1647 :
1648 : //dummy for now as init is overloaded on this signature
1649 827 : int argc=1;
1650 827 : char **argv=nullptr;
1651 827 : casa::applicator.init ( argc, argv );
1652 : //For serial or master only
1653 827 : if ( applicator.isController() ) {
1654 827 : CubeMajorCycleAlgorithm cmc;
1655 : //casa::applicator.defineAlgorithm(cmc);
1656 : //Initialize everything to get the setup as serial
1657 : {
1658 :
1659 827 : vi_p->originChunks();
1660 827 : vi_p->origin();
1661 :
1662 : }
1663 : ///Break things into chunks for parallelization or memory availabbility
1664 827 : Vector<Int> startchan;
1665 827 : Vector<Int> endchan;
1666 : Int numchunks;
1667 827 : Int fudge_factor = 15;
1668 827 : if ((itsMappers.getFTM2(0))->name()=="MosaicFTNew") {
1669 130 : fudge_factor = 20;
1670 : }
1671 697 : else if ((itsMappers.getFTM2(0))->name()=="GridFT") {
1672 619 : fudge_factor = 9;
1673 : }
1674 827 : TcleanProcessingInfo procInfo;
1675 827 : std::tie(procInfo, startchan, endchan)=nSubCubeFitInMemory(fudge_factor, itsMaxShape, gridpars_p.padding);
1676 827 : numchunks = procInfo.chnchnks;
1677 : ////TESTOO
1678 : //numchunks=2;
1679 : //startchan.resize(2);startchan[0]=0; startchan[1]=2;
1680 : //endchan.resize(2); endchan[0]=1; endchan[1]=2;
1681 :
1682 : /////END TESTOO
1683 : //cerr << "NUMCHUNKS " << numchunks << " start " << startchan << " end " << endchan << endl;
1684 827 : Record controlRecord=inpcontrol;
1685 : //For now just field 0 but should loop over all
1686 : ///This is to pass in explicit model, residual names etc
1687 827 : controlRecord.define("nfields", Int(imparsVec_p.nelements()));
1688 : //CountedPtr<SIImageStore> imstor = imageStore ( 0 );
1689 : // checking that psf, residual and sumwt is allDone
1690 : //cerr << "shapes " << imstor->residual()->shape() << " " << imstor->sumwt()->shape() << endl;
1691 827 : if(!dopsf){
1692 :
1693 : //controlRecord.define("lastcycle", savemodel);
1694 543 : controlRecord.define("nmajorcycles", nMajorCycles);
1695 543 : Vector<String> modelnames(Int(imparsVec_p.nelements()),"");
1696 1096 : for(uInt k=0; k < imparsVec_p.nelements(); ++k){
1697 553 : Int imageStoreId=k;
1698 553 : if(k>0){
1699 10 : if(gridparsVec_p[0].chanchunks >1 && uInt(itsMappers.nMappers()) >= (gridparsVec_p[0].chanchunks + gridparsVec_p.nelements()))
1700 0 : imageStoreId+=gridparsVec_p[0].chanchunks-1;
1701 10 : if(gridparsVec_p[0].facets >1)
1702 0 : imageStoreId+=gridparsVec_p[0].facets*gridparsVec_p[0].facets-1;
1703 : }
1704 553 : if((itsMappers.imageStore(imageStoreId))->hasModel()){
1705 289 : modelnames(k)=imparsVec_p[k].imageName+".model";
1706 289 : (itsMappers.imageStore(k))->model()->unlock();
1707 289 : controlRecord.define("modelnames", modelnames);
1708 : }
1709 : }
1710 :
1711 543 : }
1712 827 : Vector<String> weightnames(Int(imparsVec_p.nelements()),"");
1713 827 : Vector<String> pbnames(Int(imparsVec_p.nelements()),"");
1714 1669 : for(uInt k=0; k < imparsVec_p.nelements(); ++k){
1715 842 : Int imageStoreId=k;
1716 842 : if(k>0){
1717 15 : if(gridparsVec_p[0].chanchunks >1 && uInt(itsMappers.nMappers()) >= (gridparsVec_p[0].chanchunks + gridparsVec_p.nelements()))
1718 0 : imageStoreId+=gridparsVec_p[0].chanchunks-1;
1719 15 : if(gridparsVec_p[0].facets >1)
1720 0 : imageStoreId+=gridparsVec_p[0].facets*gridparsVec_p[0].facets-1;
1721 : }
1722 : //cerr << "FTMachine name " << (itsMappers.getFTM2(imageStoreId))->name() << endl;
1723 842 : if((itsMappers.getFTM2(imageStoreId))->useWeightImage()){
1724 : //cerr << "Mosaic weight image " << max(itsMappers.imageStore(imageStoreId)->weight(k)->get()) << endl;
1725 179 : weightnames(k)=itsMappers.imageStore(imageStoreId)->weight(k)->name();
1726 :
1727 :
1728 : }
1729 842 : if(itsMakeVP){
1730 663 : pbnames(k)=itsMappers.imageStore(imageStoreId)->pb(k)->name();
1731 663 : (itsMappers.imageStore(imageStoreId)->pb(k))->unlock();
1732 : }
1733 : }
1734 827 : controlRecord.define("weightnames", weightnames);
1735 827 : controlRecord.define("pbnames", pbnames);
1736 : // Tell the child processes not to do the dividebyweight process as this is done
1737 : // tell each child to do the normars stuff
1738 827 : controlRecord.define("dividebyweight", True);
1739 827 : controlRecord.defineRecord("normpars", normpars_p);
1740 : ///Let's see if no per chan weight density was chosen
1741 827 : String weightdensityimage=getWeightDensity();
1742 827 : if(weightdensityimage != "")
1743 2 : controlRecord.define("weightdensity", weightdensityimage);
1744 :
1745 827 : Record vecSelParsRec, vecImParsRec, vecGridParsRec;
1746 827 : Vector<Int>vecRange(2);
1747 1683 : for (uInt k = 0; k < dataSel_p.nelements(); ++k) {
1748 856 : Record selparsRec = dataSel_p[k].toRecord();
1749 856 : vecSelParsRec.defineRecord(String::toString(k), selparsRec);
1750 856 : }
1751 1669 : for (uInt k=0; k < imparsVec_p.nelements(); ++k){
1752 842 : Record imparsRec = imparsVec_p[k].toRecord();
1753 : //need to send polrep
1754 842 : imparsRec.define("polrep", Int((itsMappers.imageStore(k))->getDataPolFrame()));
1755 : //need to send movingSourceName if any
1756 842 : imparsRec.define("movingsource", movingSource_p);
1757 842 : Record gridparsRec = gridparsVec_p[k].toRecord();
1758 : /* Might need this to pass the state of the global ftmachines...test for parallel when needed
1759 : String err;
1760 : Record ftmrec, iftmrec;
1761 : if(!( (itsMappers.getFTM2(k)->toRecord(err, iftmrec,False)) && (itsMappers.getFTM2(k, false)->toRecord(err, ftmrec,False))))
1762 : throw(AipsError("FTMachines serialization failed"));
1763 : gridparsRec.defineRecord("ftmachine", ftmrec);
1764 : gridparsRec.defineRecord("iftmachine", iftmrec);
1765 : */
1766 842 : vecImParsRec.defineRecord(String::toString(k), imparsRec);
1767 842 : vecGridParsRec.defineRecord(String::toString(k), gridparsRec);
1768 842 : }
1769 827 : String workingdir="";
1770 : //Int numchan=(dopsf) ? imstor->psf()->shape()[3] : imstor->residual()->shape() [3];
1771 : //copy the imageinfo of main image here
1772 827 : if(dopsf)
1773 284 : cubePsfImageInfo_p=(itsMappers.imageStore(0)->psf())->imageInfo();
1774 1669 : for(Int k=0; k < itsMappers.nMappers(); ++k){
1775 :
1776 842 : if(dopsf){
1777 :
1778 580 : for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(true)); ++j){
1779 : ///TESTOO
1780 : //(itsMappers.imageStore(k))->psf(j)->set(0.0);
1781 : /////////
1782 291 : (itsMappers.imageStore(k))->psf(j)->unlock();
1783 291 : (itsMappers.imageStore(k))->pb()->unlock();
1784 : }
1785 : }
1786 : else{
1787 1108 : for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(false)); ++j){
1788 : /////////TESTOO
1789 : //(itsMappers.imageStore(k))->residual(j)->set(0.0);
1790 : ///////
1791 555 : (itsMappers.imageStore(k))->residual(j)->unlock();
1792 : }
1793 : }
1794 1690 : for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(true)); ++j){
1795 : //cerr << k << " type " << (itsMappers.imageStore(k))->sumwt(j)->imageType() << " name " << (itsMappers.imageStore(k))->sumwt(j)->name() << endl;
1796 :
1797 :
1798 1696 : Path namewgt( (itsMappers.imageStore(k))->sumwt(j)->name());
1799 848 : workingdir=namewgt.dirName();
1800 : ///TESTOO
1801 : //(itsMappers.imageStore(k))->sumwt(j)->set(0.0);
1802 : ////
1803 848 : (itsMappers.imageStore(k))->sumwt(j)->unlock();
1804 : //(itsMappers.imageStore(k))->releaseLocks();
1805 848 : }
1806 842 : (itsMappers.imageStore(k))->releaseLocks();
1807 : }
1808 : //Send the working directory as the child and master may be at different places
1809 :
1810 827 : controlRecord.define("workingdirectory", workingdir);
1811 : // For now this contains lastcycle if necessary in the future this
1812 : // should come from the master control record
1813 :
1814 : //Int numprocs = applicator.numProcs();
1815 : //cerr << "Number of procs: " << numprocs << endl;
1816 827 : Int rank ( 0 );
1817 : Bool assigned; //(casa::casa::applicator.nextAvailProcess(pwrite, rank));
1818 827 : Bool allDone ( false );
1819 827 : Vector<Bool> retvals(numchunks, False);
1820 827 : Int indexofretval=0;
1821 1654 : for ( Int k=0; k < numchunks; ++k ) {
1822 827 : assigned=casa::applicator.nextAvailProcess ( cmc, rank );
1823 : //cerr << "assigned "<< assigned << endl;
1824 827 : while ( !assigned ) {
1825 0 : rank = casa::applicator.nextProcessDone ( cmc, allDone );
1826 : //cerr << "while rank " << rank << endl;
1827 : Bool status;
1828 0 : Record returnRec;
1829 0 : casa::applicator.get(returnRec);
1830 0 : casa::applicator.get ( status );
1831 0 : retvals(indexofretval)=status;
1832 0 : if(dopsf)
1833 0 : updateImageBeamSet(returnRec);
1834 0 : if(returnRec.isDefined("tempfilenames")){
1835 0 : std::vector<String> b=returnRec.asArrayString("tempfilenames").tovector();
1836 0 : tempFileNames_p.insert(std::end(tempFileNames_p), std::begin(b), std::end(b));
1837 0 : }
1838 :
1839 0 : ++indexofretval;
1840 0 : if ( status )
1841 : //cerr << k << " rank " << rank << " successful " << endl;
1842 0 : cerr << "" ;
1843 : else
1844 0 : logger << LogIO::SEVERE << k << " rank " << rank << " failed " << LogIO::POST;
1845 0 : assigned = casa::applicator.nextAvailProcess ( cmc, rank );
1846 :
1847 0 : }
1848 :
1849 : ///send process info
1850 : // put data sel params #1
1851 827 : applicator.put ( vecSelParsRec );
1852 : // put image sel params #2
1853 827 : applicator.put ( vecImParsRec );
1854 : // put gridders params #3
1855 827 : applicator.put ( vecGridParsRec );
1856 : // put which channel to process #4
1857 827 : vecRange(0)=startchan(k);
1858 827 : vecRange(1)=endchan(k);
1859 827 : applicator.put ( vecRange );
1860 : // psf or residual CubeMajorCycleAlgorithm #5
1861 827 : applicator.put ( dopsf );
1862 : // store modelvis and other controls #6
1863 827 : applicator.put ( controlRecord );
1864 : // weighting scheme #7
1865 827 : applicator.put( weightParams_p);
1866 : /// Tell worker to process it
1867 827 : applicator.apply ( cmc );
1868 :
1869 : }
1870 : // Wait for all outstanding processes to return
1871 827 : rank = casa::applicator.nextProcessDone ( cmc, allDone );
1872 1654 : while ( !allDone ) {
1873 : Bool status;
1874 827 : Record returnRec;
1875 827 : casa::applicator.get(returnRec);
1876 827 : casa::applicator.get ( status );
1877 827 : if(dopsf)
1878 284 : updateImageBeamSet(returnRec);
1879 827 : if(returnRec.isDefined("tempfilenames")){
1880 807 : std::vector<String> b=returnRec.asArrayString("tempfilenames").tovector();
1881 807 : tempFileNames_p.insert(std::end(tempFileNames_p), std::begin(b), std::end(b));
1882 807 : }
1883 827 : retvals(indexofretval)=status;
1884 827 : ++indexofretval;
1885 827 : if ( status )
1886 : //cerr << "remainder rank " << rank << " successful " << endl;
1887 826 : cerr << "";
1888 : else
1889 1 : logger << LogIO::SEVERE << "remainder rank " << rank << " failed " << LogIO::POST;
1890 :
1891 827 : rank = casa::applicator.nextProcessDone ( cmc, allDone );
1892 827 : if(casa::applicator.isSerial())
1893 827 : allDone=true;
1894 827 : }
1895 827 : if(anyEQ(retvals, False)){
1896 : //cerr << retvals << endl;
1897 1 : ostringstream oss;
1898 : oss << "One or more of the cube section failed in de/gridding. Return values for "
1899 1 : "the sections: " << retvals;
1900 1 : throw(AipsError(oss));
1901 1 : }
1902 826 : if(!dopsf && normpars_p.isDefined("pblimit") && (normpars_p.asFloat("pblimit") > 0.0) ){
1903 : try{
1904 490 : SIImageStore::copyMask(itsMappers.imageStore(0)->pb(), itsMappers.imageStore(0)->residual());
1905 469 : (itsMappers.imageStore(0))->residual()->unlock();
1906 : //(itsMappers.imageStore(0)->pb())->pixelMask().unlock();
1907 469 : (itsMappers.imageStore(0))->pb()->unlock();
1908 : }
1909 3 : catch(AipsError &x) {
1910 3 : if(!String(x.getMesg()).contains("T/F"))
1911 0 : throw(AipsError(x.getMesg()));
1912 : else{
1913 3 : logger << LogIO::WARN << "Error : " << x.getMesg() << LogIO::POST;
1914 : //cout << "x.getMesg() " << endl;
1915 : }
1916 : ///ignore copy mask error and proceed as this happens with interactive
1917 3 : }
1918 : }
1919 : else{
1920 708 : LatticeLocker lock1 (*(itsMappers.imageStore(0)->psf()), FileLocker::Write);
1921 354 : itsMappers.imageStore(0)->psf()->setImageInfo(cubePsfImageInfo_p);
1922 354 : itsMappers.imageStore(0)->psf()->unlock();
1923 354 : (itsMappers.imageStore(0))->pb()->unlock();
1924 354 : }
1925 :
1926 839 : }
1927 :
1928 :
1929 826 : return true;
1930 :
1931 827 : }
1932 : /////////////////////////
1933 20 : void SynthesisImagerVi2::predictModel(){
1934 40 : LogIO os( LogOrigin("SynthesisImagerVi2","predictModel ",WHERE) );
1935 :
1936 20 : os << "---------------------------------------------------- Predict Model ---------------------------------------------" << LogIO::POST;
1937 :
1938 20 : Bool savemodelcolumn = !readOnly_p && useScratch_p;
1939 20 : Bool savevirtualmodel = !readOnly_p && !useScratch_p;
1940 : //os<<"PREDICTMODEL: readOnly_p=="<<readOnly_p<<" useScratch_p=="<<useScratch_p<<LogIO::POST;
1941 20 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
1942 20 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
1943 :
1944 20 : itsMappers.checkOverlappingModels("blank");
1945 :
1946 :
1947 : {
1948 20 : vi::VisBuffer2* vb = vi_p->getVisBuffer();;
1949 20 : vi_p->originChunks();
1950 20 : vi_p->origin();
1951 20 : Double numberCoh=0;
1952 42 : for (uInt k=0; k< mss_p.nelements(); ++k)
1953 22 : numberCoh+=Double(mss_p[k]->nrow());
1954 :
1955 40 : ProgressMeter pm(1.0, numberCoh, "Predict Model", "","","",true);
1956 20 : rownr_t cohDone=0;
1957 :
1958 20 : itsMappers.initializeDegrid(*vb);
1959 42 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
1960 : {
1961 :
1962 8750 : for (vi_p->origin(); vi_p->more(); vi_p->next())
1963 : {
1964 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; //No valid rows in this MS
1965 : //if !usescratch ...just save
1966 8728 : vb->setVisCubeModel(Complex(0.0, 0.0));
1967 8728 : itsMappers.degrid(*vb, savevirtualmodel);
1968 :
1969 8728 : if(savemodelcolumn && writeAccess_p ){
1970 :
1971 4888 : const_cast<MeasurementSet& >((vi_p->ms())).lock(true);
1972 :
1973 4888 : vi_p->writeVisModel(vb->visCubeModel());
1974 :
1975 : //cerr << "nRows "<< vb->nRows() << " " << max(vb->visCubeModel()) << endl;
1976 4888 : const_cast<MeasurementSet& >((vi_p->ms())).unlock();
1977 : }
1978 :
1979 8728 : cohDone += vb->nRows();
1980 8728 : pm.update(Double(cohDone));
1981 :
1982 : }
1983 : }
1984 20 : itsMappers.finalizeDegrid(*vb);
1985 20 : }
1986 :
1987 20 : itsMappers.checkOverlappingModels("restore");
1988 20 : itsMappers.releaseImageLocks();
1989 20 : unlockMSs();
1990 :
1991 20 : }// end of predictModel
1992 :
1993 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1994 144 : void SynthesisImagerVi2::makeSdImage(Bool dopsf)
1995 : {
1996 288 : LogIO os( LogOrigin("SynthesisImagerVi2","makeSdImage",WHERE) );
1997 :
1998 : // Bool dopsf=false;
1999 144 : if(datacol_p==FTMachine::PSF) dopsf=true;
2000 :
2001 : {
2002 144 : vi::VisBuffer2* vb = vi_p->getVisBuffer();;
2003 144 : vi_p->originChunks();
2004 144 : vi_p->origin();
2005 :
2006 144 : Double numberCoh=0;
2007 314 : for (uInt k=0; k< mss_p.nelements(); ++k)
2008 170 : numberCoh+=Double(mss_p[k]->nrow());
2009 :
2010 288 : ProgressMeter pm(1.0, numberCoh, "Predict Model", "","","",true);
2011 144 : rownr_t cohDone=0;
2012 :
2013 144 : itsMappers.initializeGrid(*vi_p,dopsf);
2014 766 : for (vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
2015 : {
2016 622 : if (vi_p->getImpl()->isNewMs()) {
2017 170 : itsMappers.handleNewMs(vi_p->ms());
2018 : }
2019 204441 : for (vi_p->origin(); vi_p->more(); vi_p->next())
2020 : {
2021 203819 : itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)datacol_p);
2022 203819 : cohDone += vb->nRows();
2023 203819 : pm.update(Double(cohDone));
2024 : }
2025 : }
2026 144 : itsMappers.finalizeGrid(*vb, dopsf);
2027 :
2028 144 : }
2029 :
2030 144 : unlockMSs();
2031 :
2032 144 : }// end makeSDImage
2033 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2034 2 : void SynthesisImagerVi2::makeImage(String type, const String& imagename, const String& complexImage, const Int whichModel)
2035 : {
2036 4 : LogIO os( LogOrigin("SynthesisImager","makeImage",WHERE) );
2037 :
2038 :
2039 2 : refim::FTMachine::Type seType(refim::FTMachine::OBSERVED);
2040 2 : if(type=="observed") {
2041 0 : seType=refim::FTMachine::OBSERVED;
2042 : os << LogIO::NORMAL // Loglevel INFO
2043 : << "Making dirty image from " << type << " data "
2044 0 : << LogIO::POST;
2045 : }
2046 2 : else if (type=="model") {
2047 0 : seType=refim::FTMachine::MODEL;
2048 : os << LogIO::NORMAL // Loglevel INFO
2049 : << "Making dirty image from " << type << " data "
2050 0 : << LogIO::POST;
2051 : }
2052 2 : else if (type=="corrected") {
2053 0 : seType=refim::FTMachine::CORRECTED;
2054 : os << LogIO::NORMAL // Loglevel INFO
2055 : << "Making dirty image from " << type << " data "
2056 0 : << LogIO::POST;
2057 : }
2058 2 : else if (type=="psf") {
2059 2 : seType=refim::FTMachine::PSF;
2060 : os << "Making point spread function "
2061 2 : << LogIO::POST;
2062 : }
2063 0 : else if (type=="residual") {
2064 0 : seType=refim::FTMachine::RESIDUAL;
2065 : os << LogIO::NORMAL // Loglevel INFO
2066 : << "Making dirty image from " << type << " data "
2067 0 : << LogIO::POST;
2068 : }
2069 0 : else if (type=="singledish-observed") {
2070 0 : seType=refim::FTMachine::OBSERVED;
2071 : os << LogIO::NORMAL
2072 0 : << "Making single dish image from observed data" << LogIO::POST;
2073 : }
2074 0 : else if (type=="singledish") {
2075 0 : seType=refim::FTMachine::CORRECTED;
2076 : os << LogIO::NORMAL
2077 0 : << "Making single dish image from corrected data" << LogIO::POST;
2078 : }
2079 0 : else if (type=="coverage") {
2080 0 : seType=refim::FTMachine::COVERAGE;
2081 : os << LogIO::NORMAL
2082 : << "Making single dish coverage function "
2083 0 : << LogIO::POST;
2084 : }
2085 0 : else if (type=="holography") {
2086 0 : seType=refim::FTMachine::CORRECTED;
2087 : os << LogIO::NORMAL
2088 : << "Making complex holographic image from corrected data "
2089 0 : << LogIO::POST;
2090 : }
2091 0 : else if (type=="holography-observed") {
2092 0 : seType=refim::FTMachine::OBSERVED;
2093 : os << LogIO::NORMAL
2094 : << "Making complex holographic image from observed data "
2095 0 : << LogIO::POST;
2096 : }
2097 :
2098 2 : String imageName=(itsMappers.imageStore(whichModel))->getName();
2099 2 : String cImageName(complexImage);
2100 2 : if(complexImage=="") {
2101 2 : cImageName=imageName+".compleximage";
2102 : }
2103 2 : Bool keepComplexImage=(complexImage!="")||(type.contains("holography"));
2104 : //cerr << "name " << (itsMappers.getFTM2(whichModel))->name() << endl;
2105 2 : if(((itsMappers.getFTM2(whichModel))->name())!="MultiTermFTNew"){
2106 : ////Non multiterm case
2107 : //cerr << "whichModel " << whichModel << " itsMappers num " << itsMappers.nMappers() << " shape " << (itsMappers.imageStore(whichModel))->getShape() << endl;
2108 : ///the SIImageStore sometimes return 0 shape...
2109 2 : CoordinateSystem csys=itsMaxCoordSys;
2110 2 : IPosition shp=itsMaxShape;
2111 2 : if((itsMappers.imageStore(whichModel))->getShape().product()!=0 ){
2112 2 : csys= (itsMappers.imageStore(whichModel))->getCSys();
2113 2 : shp=(itsMappers.imageStore(whichModel))->getShape();
2114 : }
2115 2 : itsMappers.releaseImageLocks();
2116 2 : PagedImage<Float> theImage( shp, csys, imagename);
2117 0 : PagedImage<Complex> cImageImage(theImage.shape(),
2118 : theImage.coordinates(),
2119 4 : cImageName);
2120 2 : if(!keepComplexImage)
2121 2 : cImageImage.table().markForDelete();
2122 :
2123 :
2124 2 : Matrix<Float> weight;
2125 2 : if(cImageImage.shape()[3] > 1){
2126 1 : cImageImage.set(Complex(0.0));
2127 1 : cImageImage.tempClose();
2128 1 : makeComplexCubeImage(cImageName, seType, whichModel);
2129 : }
2130 : else{
2131 1 : (itsMappers.getFTM2(whichModel))->makeImage(seType, *vi_p, cImageImage, weight);
2132 : }
2133 2 : if(seType==refim::FTMachine::PSF){
2134 2 : StokesImageUtil::ToStokesPSF(theImage, cImageImage);
2135 2 : StokesImageUtil::normalizePSF(theImage);
2136 : }
2137 : else{
2138 0 : StokesImageUtil::To(theImage, cImageImage);
2139 : }
2140 2 : }
2141 : else{
2142 : ///Multiterm
2143 : //refim::MultiTermFTNew *theft=static_cast<refim::MultiTermFTNew *>( (itsMappers.getFTM2(whichModel))->cloneFTM());
2144 0 : refim::MultiTermFTNew *theft=static_cast<refim::MultiTermFTNew *>( (itsMappers.getFTM2(whichModel)).get());
2145 0 : Int ntaylor= seType==refim::FTMachine::PSF ? theft->psfNTerms() : theft->nTerms();
2146 0 : if(ntaylor<2)
2147 0 : throw(AipsError("some issue with muti term setting "));
2148 0 : Vector<CountedPtr<ImageInterface<Float> > >theImage(ntaylor);
2149 0 : Vector<CountedPtr<ImageInterface<Complex> > >cImageImage(ntaylor);
2150 0 : Vector<CountedPtr<Matrix<Float> > >weight(ntaylor);
2151 0 : for (Int taylor=0; taylor < ntaylor; ++taylor){
2152 0 : theImage[taylor]=new PagedImage<Float>((itsMappers.imageStore(whichModel))->getShape(), (itsMappers.imageStore(whichModel))->getCSys(), imagename+".tt"+String::toString(taylor));
2153 0 : cImageImage[taylor]=new PagedImage<Complex> (theImage[taylor]->shape(),
2154 0 : theImage[taylor]->coordinates(),
2155 0 : cImageName+".tt"+String::toString(taylor));
2156 0 : if(!keepComplexImage)
2157 0 : static_cast<PagedImage<Complex> *> (cImageImage[taylor].get())->table().markForDelete();
2158 0 : weight[taylor]=new Matrix<Float>();
2159 :
2160 : }
2161 0 : theft->makeMTImages(seType, *vi_p, cImageImage, weight);
2162 0 : Float maxpsf=1.0;
2163 0 : for (Int taylor=0; taylor < ntaylor; ++taylor){
2164 0 : if(seType==refim::FTMachine::PSF){
2165 0 : StokesImageUtil::ToStokesPSF(*(theImage[taylor]), *(cImageImage[taylor]));
2166 0 : if(taylor==0){
2167 0 : maxpsf=StokesImageUtil::normalizePSF(*theImage[taylor]);
2168 : //cerr << "maxpsf " << maxpsf << endl;
2169 : }
2170 : else{
2171 : ///divide by max;
2172 0 : (*theImage[taylor]).copyData((LatticeExpr<Float>)((*theImage[taylor])/maxpsf));
2173 : }
2174 : }
2175 : else{
2176 0 : StokesImageUtil::To(*(theImage[taylor]), *(cImageImage[taylor]));
2177 : }
2178 : }
2179 : //delete theft;
2180 :
2181 0 : }
2182 2 : unlockMSs();
2183 :
2184 2 : }// end makeImage
2185 : /////////////////////////////////////////////////////
2186 :
2187 1 : void SynthesisImagerVi2::makeComplexCubeImage(const String& cimage, const refim::FTMachine::Type imtype, const Int whichmodel){
2188 1 : CubeMakeImageAlgorithm *cmi=new CubeMakeImageAlgorithm();
2189 : // Dummies for using the right signature for init
2190 1 : Path cimpath(cimage);
2191 1 : String cimname=cimpath.absoluteName();
2192 : //cerr << "ABSOLUTE path " << cimname << endl;
2193 1 : Int argc = 1;
2194 1 : char **argv = nullptr;
2195 1 : casa::applicator.init(argc, argv);
2196 1 : if(applicator.isController())
2197 : {
2198 1 : Record vecSelParsRec;
2199 2 : for (uInt k = 0; k < dataSel_p.nelements(); ++k) {
2200 1 : Record selparsRec = dataSel_p[k].toRecord();
2201 1 : vecSelParsRec.defineRecord(String::toString(k), selparsRec);
2202 1 : }
2203 1 : Record imparsRec = imparsVec_p[whichmodel].toRecord();
2204 : //cerr << "which model " << whichmodel << " record " << imparsRec << endl;
2205 1 : Record gridparsRec = gridparsVec_p[whichmodel].toRecord();
2206 : ///Break things into chunks for parallelization or memory availabbility
2207 1 : Vector<Int> startchan;
2208 1 : Vector<Int> endchan;
2209 : Int numchunks;
2210 1 : Int fudge_factor = 15;
2211 1 : if ((itsMappers.getFTM2(0))->name()=="MosaicFTNew") {
2212 1 : fudge_factor = 20;
2213 : }
2214 0 : else if ((itsMappers.getFTM2(0))->name()=="GridFT") {
2215 0 : fudge_factor = 9;
2216 : }
2217 1 : TcleanProcessingInfo procInfo;
2218 1 : std::tie(procInfo, startchan, endchan)=nSubCubeFitInMemory(fudge_factor, itsMaxShape, gridpars_p.padding);
2219 1 : numchunks = procInfo.chnchnks;
2220 :
2221 1 : Int imageType=Int(imtype);
2222 1 : Int rank(0);
2223 : Bool assigned;
2224 1 : Bool allDone(false);
2225 1 : Vector<Int> chanRange(2);
2226 2 : for (Int k=0; k < numchunks; ++k) {
2227 1 : assigned=casa::applicator.nextAvailProcess ( *cmi, rank );
2228 : //cerr << "assigned "<< assigned << endl;
2229 1 : while ( !assigned ) {
2230 0 : rank = casa::applicator.nextProcessDone ( *cmi, allDone );
2231 : //cerr << "while rank " << rank << endl;
2232 : Bool status;
2233 0 : casa::applicator.get(status);
2234 : /*if ( status )
2235 : cerr << k << " rank " << rank << " successful " << endl;
2236 : else
2237 : cerr << k << " rank " << rank << " failed " << endl;
2238 : */
2239 0 : assigned = casa::applicator.nextAvailProcess ( *cmi, rank );
2240 :
2241 : }
2242 1 : applicator.put(vecSelParsRec);
2243 : // put image sel params #2
2244 1 : applicator.put(imparsRec);
2245 : // put gridder params #3
2246 1 : applicator.put(gridparsRec);
2247 : // put which channel to process #4
2248 1 : chanRange(0)=startchan(k);
2249 1 : chanRange(1)=endchan(k);
2250 1 : applicator.put(chanRange);
2251 : //Type of image #5
2252 1 : applicator.put(imageType);
2253 : // weighting parameters #6
2254 1 : applicator.put( weightParams_p);
2255 : // complex imagename #7
2256 1 : applicator.put(cimname);
2257 1 : applicator.apply(*cmi);
2258 : }
2259 : // Wait for all outstanding processes to return
2260 1 : rank = casa::applicator.nextProcessDone(*cmi, allDone);
2261 2 : while (!allDone) {
2262 : Bool status;
2263 1 : casa::applicator.get(status);
2264 : /*
2265 : if(status)
2266 : cerr << "remainder rank " << rank << " successful " << endl;
2267 : else
2268 : cerr << "remainder rank " << rank << " failed " << endl;
2269 : */
2270 1 : rank = casa::applicator.nextProcessDone(*cmi, allDone);
2271 1 : if(casa::applicator.isSerial())
2272 1 : allDone=true;
2273 : }
2274 1 : }//applicator controller section
2275 1 : }
2276 :
2277 :
2278 :
2279 1768 : CountedPtr<SIMapper> SynthesisImagerVi2::createSIMapper(String mappertype,
2280 : CountedPtr<SIImageStore> imagestore,
2281 : CountedPtr<refim::FTMachine> ftmachine,
2282 : CountedPtr<refim::FTMachine> iftmachine,
2283 : uInt /*ntaylorterms*/)
2284 : {
2285 3536 : LogIO os( LogOrigin("SynthesisImagerVi2","createSIMapper",WHERE) );
2286 :
2287 1768 : CountedPtr<SIMapper> localMapper;
2288 :
2289 : try
2290 : {
2291 :
2292 1768 : if( mappertype == "default" || mappertype == "multiterm" )
2293 : {
2294 1768 : localMapper = new SIMapper( imagestore, ftmachine, iftmachine );
2295 : }
2296 0 : else if( mappertype == "imagemosaic") // || mappertype == "mtimagemosaic" )
2297 : {
2298 0 : localMapper = new SIMapperImageMosaic( imagestore, ftmachine, iftmachine );
2299 : }
2300 : else
2301 : {
2302 0 : throw(AipsError("Unknown mapper type : " + mappertype));
2303 : }
2304 :
2305 : }
2306 0 : catch(AipsError &x) {
2307 0 : throw(AipsError("Error in createSIMapper : " + x.getMesg() ) );
2308 0 : }
2309 3536 : return localMapper;
2310 1768 : }
2311 :
2312 3688 : void SynthesisImagerVi2::lockMS(MeasurementSet& thisms){
2313 3688 : thisms.lock(!readOnly_p);
2314 3688 : thisms.antenna().lock(false);
2315 3688 : thisms.dataDescription().lock(false);
2316 3688 : thisms.feed().lock(false);
2317 3688 : thisms.field().lock(false);
2318 3688 : thisms.observation().lock(false);
2319 3688 : thisms.polarization().lock(false);
2320 3688 : thisms.processor().lock(false);
2321 3688 : thisms.spectralWindow().lock(false);
2322 3688 : thisms.state().lock(false);
2323 3688 : if(!thisms.doppler().isNull())
2324 150 : thisms.doppler().lock(false);
2325 3688 : if(!thisms.flagCmd().isNull())
2326 3688 : thisms.flagCmd().lock(false);
2327 3688 : if(!thisms.freqOffset().isNull())
2328 150 : thisms.freqOffset().lock(false);
2329 : //True here as we can write in that
2330 3688 : if(!thisms.history().isNull())
2331 : // thisms.history().lock(!readOnly_p);
2332 3688 : thisms.history().lock(false);
2333 3688 : if(!thisms.pointing().isNull())
2334 3688 : thisms.pointing().lock(false);
2335 : //we write virtual model here
2336 3688 : if(!thisms.source().isNull())
2337 3688 : thisms.source().lock(!readOnly_p);
2338 3688 : if(!thisms.sysCal().isNull())
2339 315 : thisms.sysCal().lock(false);
2340 3688 : if(!thisms.weather().isNull())
2341 315 : thisms.weather().lock(false);
2342 :
2343 3688 : }
2344 : ///////////////////////
2345 :
2346 821 : Vector<SynthesisParamsSelect> SynthesisImagerVi2::tuneSelectData(){
2347 821 : vi_p->originChunks();
2348 821 : vi_p->origin();
2349 821 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
2350 821 : Int spwnow=vb->spectralWindows()[0];
2351 821 : Int nchaninms=MSColumns(vb->ms()).spectralWindow().numChan()(spwnow);
2352 : //For some small number of channels in the ms vi/vb2 retuning the selection
2353 : //will sometimes return more channels than what is in the ms...this is a
2354 : //kludge here to bypass it...mostly seen in test_tclean
2355 : /// write to the test !! till someboody fixes this is vi2 or wait for cngi
2356 : //if savescratch column we have tune...otherwise some channel may be 0
2357 : // when chunking or in parallel
2358 : //cerr << "nchanims " << nchaninms << endl;
2359 821 : if(nchaninms <30 && !(!readOnly_p && useScratch_p))
2360 810 : return dataSel_p;
2361 :
2362 11 : return SynthesisImager::tuneSelectData();
2363 : }
2364 : //////////////////////
2365 1800 : void SynthesisImagerVi2::lockMSs(){
2366 3684 : for(uInt i=0;i<mss_p.nelements();i++)
2367 : {
2368 :
2369 1884 : MeasurementSet *ms_l = const_cast<MeasurementSet* >(mss_p[i]);
2370 1884 : lockMS(*ms_l);
2371 : }
2372 :
2373 1800 : }
2374 8272 : void SynthesisImagerVi2::unlockMSs()
2375 : {
2376 16544 : LogIO os( LogOrigin("SynthesisImagerVi2","unlockMSs",WHERE) );
2377 16915 : for(uInt i=0;i<mss_p.nelements();i++)
2378 : {
2379 8643 : os << LogIO::NORMAL2 << "Unlocking : " << (mss_p[i])->tableName() << LogIO::POST;
2380 8643 : MeasurementSet *ms_l = const_cast<MeasurementSet* >(mss_p[i]);
2381 8643 : ms_l->unlock();
2382 8643 : ms_l->antenna().unlock();
2383 8643 : ms_l->dataDescription().unlock();
2384 8643 : ms_l->feed().unlock();
2385 8643 : ms_l->field().unlock();
2386 8643 : ms_l->observation().unlock();
2387 8643 : ms_l->polarization().unlock();
2388 8643 : ms_l->processor().unlock();
2389 8643 : ms_l->spectralWindow().unlock();
2390 8643 : ms_l->state().unlock();
2391 : //
2392 : // Unlock the optional sub-tables as well, if they are present
2393 : //
2394 8643 : if(!(ms_l->source().isNull())) ms_l->source().unlock();
2395 8643 : if(!(ms_l->doppler().isNull())) ms_l->doppler().unlock();
2396 8643 : if(!(ms_l->flagCmd().isNull())) ms_l->flagCmd().unlock();
2397 8643 : if(!(ms_l->freqOffset().isNull())) ms_l->freqOffset().unlock();
2398 8643 : if(!(ms_l->history().isNull())) ms_l->history().unlock();
2399 8643 : if(!(ms_l->pointing().isNull())) ms_l->pointing().unlock();
2400 8643 : if(!(ms_l->sysCal().isNull())) ms_l->sysCal().unlock();
2401 8643 : if(!(ms_l->weather().isNull())) ms_l->weather().unlock();
2402 : }
2403 8272 : }
2404 1732 : void SynthesisImagerVi2::createFTMachine(
2405 : CountedPtr<refim::FTMachine>& theFT,
2406 : CountedPtr<refim::FTMachine>& theIFT,
2407 : const String& ftname,
2408 : const uInt nTaylorTerms,
2409 : const String mType,
2410 : const Int facets, //=1
2411 : //------------------------------
2412 : const Int wprojplane, //=1,
2413 : const Float padding, //=1.0,
2414 : const Bool useAutocorr, //=false,
2415 : const Bool useDoublePrec, //=true,
2416 : const String gridFunction, //=String("SF"),
2417 : //------------------------------
2418 : const Bool aTermOn, //= true,
2419 : const Bool psTermOn, //= true,
2420 : const Bool mTermOn, //= false,
2421 : const Bool wbAWP, //= true,
2422 : const String cfCache, //= "",
2423 : const Bool usePointing, //= false,
2424 : // const Vector<Float> pointingOffsetSigDev, //= 10.0,
2425 : const vector<float> pointingOffsetSigDev, // = {10,10}
2426 : const Bool doPBCorr, //= true,
2427 : const Bool conjBeams, //= true,
2428 : const Float computePAStep, //=360.0
2429 : const Float rotatePAStep, //=5.0
2430 : const String interpolation, //="linear"
2431 : const Bool freqFrameValid, //=true
2432 : const Int cache, //=1000000000,
2433 : const Int tile, //=16
2434 : const String stokes, //=I
2435 : const String imageNamePrefix,
2436 : //---------------------------
2437 : const String &pointingDirCol,
2438 : const String &convertFirst,
2439 : const Float skyPosThreshold,
2440 : const Int convSupport,
2441 : const Quantity &truncateSize,
2442 : const Quantity &gwidth,
2443 : const Quantity &jwidth,
2444 : const Float minWeight,
2445 : const Bool clipMinMax,
2446 : const Bool pseudoI
2447 : )
2448 :
2449 : {
2450 3464 : LogIO os( LogOrigin("SynthesisImagerVi2","createFTMachine",WHERE));
2451 :
2452 :
2453 1732 : if (ftname == "gridft") {
2454 1242 : if (facets >1) {
2455 : theFT = new refim::GridFT(
2456 : cache, tile,
2457 6 : gridFunction, mLocation_p, phaseCenter_p, padding,
2458 : useAutocorr, useDoublePrec
2459 3 : );
2460 : theIFT = new refim::GridFT(
2461 : cache, tile,
2462 6 : gridFunction, mLocation_p, phaseCenter_p, padding,
2463 : useAutocorr, useDoublePrec
2464 3 : );
2465 : }
2466 : else {
2467 : theFT = new refim::GridFT(
2468 : cache, tile,
2469 1239 : gridFunction, mLocation_p, padding,
2470 : useAutocorr, useDoublePrec
2471 1239 : );
2472 : theIFT = new refim::GridFT(
2473 : cache, tile,
2474 1239 : gridFunction, mLocation_p, padding,
2475 : useAutocorr, useDoublePrec
2476 1239 : );
2477 : }
2478 : }
2479 490 : else if (ftname == "wprojectft") {
2480 45 : Double maxW = -1.0;
2481 45 : Double minW = -1.0;
2482 45 : Double rmsW = -1.0;
2483 45 : if (wprojplane < 1)
2484 7 : casa::refim::WProjectFT::wStat(*vi_p, minW, maxW, rmsW);
2485 45 : if (facets > 1) {
2486 2 : theFT = new refim::WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
2487 1 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2488 2 : theIFT = new refim::WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
2489 1 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2490 : }
2491 : else {
2492 44 : theFT = new refim::WProjectFT(wprojplane, mLocation_p,
2493 44 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2494 44 : theIFT = new refim::WProjectFT(wprojplane, mLocation_p,
2495 44 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2496 : }
2497 : CountedPtr<refim::WPConvFunc> sharedconvFunc =
2498 45 : static_cast<refim::WProjectFT &>(*theFT).getConvFunc();
2499 : //static_cast<WProjectFT &>(*theFT).setConvFunc(sharedconvFunc);
2500 45 : static_cast<refim::WProjectFT &>(*theIFT).setConvFunc(sharedconvFunc);
2501 45 : }
2502 :
2503 445 : else if ( ftname == "mosaic" || ftname== "mosft" || ftname == "mosaicft" || ftname== "MosaicFT" || ftname == "awp2"){
2504 :
2505 208 : createMosFTMachine(theFT, theIFT, padding, useAutocorr, useDoublePrec, rotatePAStep, stokes, conjBeams);
2506 : }
2507 237 : else if ((ftname.at(0,3)=="awp") || (ftname== "mawprojectft") || (ftname == "protoft")) {
2508 93 : createAWPFTMachine(theFT, theIFT, ftname, facets, wprojplane,
2509 : padding, useAutocorr, useDoublePrec, gridFunction,
2510 : aTermOn, psTermOn, mTermOn, wbAWP, cfCache,
2511 : usePointing, pointingOffsetSigDev, doPBCorr, conjBeams, computePAStep,
2512 : rotatePAStep, cache,tile,imageNamePrefix);
2513 : }
2514 :
2515 288 : else if ( ftname == "mosaic" or
2516 288 : ftname == "mosft" or
2517 432 : ftname == "mosaicft" or
2518 144 : ftname == "MosaicFT" ) {
2519 0 : createMosFTMachine(
2520 : theFT, theIFT,
2521 : padding, useAutocorr, useDoublePrec, rotatePAStep, stokes, conjBeams
2522 : );
2523 : }
2524 144 : else if ( ftname == "sd" ) {
2525 144 : createSDFTMachine(
2526 : theFT, theIFT,
2527 : pointingDirCol, convertFirst,
2528 : skyPosThreshold, doPBCorr, rotatePAStep,
2529 : gridFunction, convSupport, truncateSize, gwidth, jwidth,
2530 : minWeight, clipMinMax, cache, tile, stokes
2531 : );
2532 : }
2533 : else {
2534 0 : throw( AipsError( "Invalid FTMachine name : " + ftname ) );
2535 : }
2536 :
2537 : // Now, clone and pack the chosen FT into a MultiTermFT if needed.
2538 1732 : if( mType == "multiterm" ) {
2539 119 : AlwaysAssert( nTaylorTerms>=1 , AipsError );
2540 :
2541 : CountedPtr<refim::FTMachine> theMTFT =
2542 119 : new refim::MultiTermFTNew( theFT, nTaylorTerms, true /*forward*/ );
2543 : CountedPtr<refim::FTMachine> theMTIFT =
2544 119 : new refim::MultiTermFTNew( theIFT, nTaylorTerms, false /*forward*/ );
2545 :
2546 119 : theFT = theMTFT;
2547 119 : theIFT = theMTIFT;
2548 119 : }
2549 :
2550 : // Now, set the SkyJones if needed, and if not internally generated.
2551 1732 : if ( mType == "imagemosaic" and
2552 0 : ( ftname != "awprojectft" and
2553 0 : ftname != "mawprojectft" and
2554 0 : ftname != "proroft" ) ) {
2555 0 : CountedPtr<refim::SkyJones> vp;
2556 0 : MSColumns msc(*(mss_p[0]));
2557 0 : Quantity parang(0.0,"deg");
2558 0 : Quantity skyposthreshold(0.0,"deg");
2559 0 : vp = new refim::VPSkyJones(
2560 : msc, true, parang, BeamSquint::NONE, skyposthreshold
2561 0 : );
2562 :
2563 0 : Vector<CountedPtr<refim::SkyJones> > skyJonesList(1);
2564 0 : skyJonesList(0) = vp;
2565 0 : theFT->setSkyJones( skyJonesList );
2566 0 : theIFT->setSkyJones( skyJonesList );
2567 0 : }
2568 :
2569 : // For mode=cubedata, set the freq frame to invalid.
2570 : // get this info from buildCoordSystem
2571 : // theFT->setSpw( tspws, false );
2572 : // theIFT->setSpw( tspws, false );
2573 1732 : theFT->setFrameValidity( freqFrameValid );
2574 1732 : theIFT->setFrameValidity( freqFrameValid );
2575 :
2576 : // Set interpolation mode
2577 1732 : theFT->setFreqInterpolation( interpolation );
2578 1732 : theIFT->setFreqInterpolation( interpolation );
2579 :
2580 : // Set tracking of moving source if any
2581 1732 : if (movingSource_p != "") {
2582 22 : theFT->setMovingSource(movingSource_p);
2583 22 : theIFT->setMovingSource(movingSource_p);
2584 : }
2585 : /* vi_p has chanselection now
2586 : //channel selections from spw param
2587 : theFT->setSpwChanSelection(chanSel_p);
2588 : theIFT->setSpwChanSelection(chanSel_p);
2589 : */
2590 :
2591 : // Set pseudo-I if requested.
2592 1732 : if (pseudoI == true) {
2593 1 : os << "Turning on Pseudo-I gridding" << LogIO::POST;
2594 1 : theFT->setPseudoIStokes(true);
2595 1 : theIFT->setPseudoIStokes(true);
2596 : }
2597 :
2598 1732 : }
2599 :
2600 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2601 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2602 93 : void SynthesisImagerVi2::createAWPFTMachine(CountedPtr<refim::FTMachine>& theFT, CountedPtr<refim::FTMachine>& theIFT,
2603 : const String& ftmName,
2604 : const Int,// facets, //=1
2605 : //------------------------------
2606 : const Int wprojPlane, //=1,
2607 : const Float,// padding, //=1.0,
2608 : const Bool,// useAutocorr, //=false,
2609 : const Bool useDoublePrec, //=true,
2610 : const String,// gridFunction, //=String("SF"),
2611 : //------------------------------
2612 : const Bool aTermOn, //= true,
2613 : const Bool psTermOn, //= true,
2614 : const Bool mTermOn, //= false,
2615 : const Bool wbAWP, //= true,
2616 : const String cfCache, //= "",
2617 : const Bool usePointing, //= false,
2618 : const vector<Float> pointingOffsetSigDev,//={10,10},
2619 : const Bool doPBCorr, //= true,
2620 : const Bool conjBeams, //= true,
2621 : const Float computePAStep, //=360.0
2622 : const Float rotatePAStep, //=5.0
2623 : const Int cache, //=1000000000,
2624 : const Int tile, //=16
2625 : const String imageNamePrefix
2626 : )
2627 :
2628 : {
2629 186 : LogIO os( LogOrigin("SynthesisImagerVi2","createAWPFTMachine",WHERE));
2630 :
2631 93 : if (wprojPlane<=1)
2632 : {
2633 : os << LogIO::NORMAL
2634 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
2635 92 : << LogIO::POST;
2636 92 : os << LogIO::NORMAL << "Performing WBA-Projection" << LogIO::POST; // Loglevel PROGRESS
2637 : }
2638 : // if((wprojPlane>1)&&(wprojPlane<64))
2639 : // {
2640 : // os << LogIO::WARN
2641 : // << "No. of w-planes set too low for W projection - recommend at least 128"
2642 : // << LogIO::POST;
2643 : // os << LogIO::NORMAL << "Performing WBAW-Projection" << LogIO::POST; // Loglevel PROGRESS
2644 : // }
2645 :
2646 : // CountedPtr<ATerm> apertureFunction = createTelescopeATerm(mss4vi_p[0], aTermOn);
2647 : // CountedPtr<PSTerm> psTerm = new PSTerm();
2648 : // CountedPtr<WTerm> wTerm = new WTerm();
2649 :
2650 : // //
2651 : // // Selectively switch off CFTerms.
2652 : // //
2653 : // if (aTermOn == false) {apertureFunction->setOpCode(CFTerms::NOOP);}
2654 : // if (psTermOn == false) psTerm->setOpCode(CFTerms::NOOP);
2655 :
2656 : // //
2657 : // // Construct the CF object with appropriate CFTerms.
2658 : // //
2659 : // CountedPtr<ConvolutionFunction> tt;
2660 : // tt = AWProjectFT::makeCFObject(aTermOn, psTermOn, true, mTermOn, wbAWP);
2661 : // CountedPtr<ConvolutionFunction> awConvFunc;
2662 : // // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm, !wbAWP);
2663 : // if ((ftmName=="mawprojectft") || (mTermOn))
2664 : // awConvFunc = new AWConvFuncEPJones(apertureFunction,psTerm,wTerm,wbAWP);
2665 : // else
2666 : // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm,wbAWP);
2667 :
2668 93 : MSObservationColumns msoc((mss_p[0])->observation());
2669 93 : String telescopeName=msoc.telescopeName()(0);
2670 : CountedPtr<refim::ConvolutionFunction> awConvFunc = refim::AWProjectFT::makeCFObject(telescopeName,
2671 : aTermOn,
2672 : psTermOn, (wprojPlane > 1),
2673 93 : mTermOn, wbAWP, conjBeams);
2674 :
2675 : //
2676 : // Construct the appropriate re-sampler.
2677 : //
2678 93 : CountedPtr<refim::VisibilityResamplerBase> visResampler;
2679 : // if (ftmName=="protoft") visResampler = new ProtoVR();
2680 : //elsef
2681 93 : visResampler = new refim::AWVisResampler();
2682 : // CountedPtr<VisibilityResamplerBase> visResampler = new VisibilityResampler();
2683 :
2684 : //
2685 : // Construct and initialize the CF cache object.
2686 : //
2687 :
2688 :
2689 : // CountedPtr<CFCache> cfCacheObj = new CFCache();
2690 : // cfCacheObj->setCacheDir(cfCache.data());
2691 : // // cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
2692 : // cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
2693 : // cfCacheObj->initCache2();
2694 :
2695 93 : CountedPtr<refim::CFCache> cfCacheObj;
2696 :
2697 :
2698 : //
2699 : // Finally construct the FTMachine with the CFCache, ConvFunc and
2700 : // Re-sampler objects.
2701 : //
2702 93 : Float pbLimit_l=1e-3;
2703 :
2704 93 : theFT = new refim::AWProjectWBFT(wprojPlane, cache/2,
2705 : cfCacheObj, awConvFunc,
2706 : visResampler,
2707 : /*true */usePointing, pointingOffsetSigDev ,doPBCorr,
2708 : tile, computePAStep, pbLimit_l, true,conjBeams,
2709 93 : useDoublePrec);
2710 :
2711 93 : cfCacheObj = new refim::CFCache();
2712 93 : cfCacheObj->setCacheDir(cfCache.data());
2713 : // Get the LAZYFILL setting from the user configuration. If not
2714 : // found, default to False.
2715 : //
2716 : // With lazy fill ON, CFCache loads the required CFs on-demand
2717 : // from the disk. And periodically triggers garbage collection to
2718 : // release CFs that aren't required immediately.
2719 93 : if(impars_p.mode.contains("cube")){
2720 19 : cfCacheObj->setLazyFill(False);
2721 : }
2722 : else
2723 74 : cfCacheObj->setLazyFill(refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1);
2724 : // cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
2725 93 : cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
2726 93 : cfCacheObj->initCache2(CFC_VERBOSE);
2727 :
2728 93 : theFT->setCFCache(cfCacheObj);
2729 :
2730 :
2731 93 : Quantity rotateOTF(rotatePAStep,"deg");
2732 93 : static_cast<refim::AWProjectWBFT &>(*theFT).setObservatoryLocation(mLocation_p);
2733 93 : static_cast<refim::AWProjectWBFT &>(*theFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
2734 :
2735 : // theIFT = new AWProjectWBFT(wprojPlane, cache/2,
2736 : // cfCacheObj, awConvFunc,
2737 : // visResampler,
2738 : // /*true */usePointing, doPBCorr,
2739 : // tile, computePAStep, pbLimit_l, true,conjBeams,
2740 : // useDoublePrec);
2741 :
2742 : // static_cast<AWProjectWBFT &>(*theIFT).setObservatoryLocation(mLocation_p);
2743 : // static_cast<AWProjectWBFT &>(*theIFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
2744 :
2745 93 : theIFT = new refim::AWProjectWBFT(static_cast<refim::AWProjectWBFT &>(*theFT));
2746 :
2747 93 : os << "Sending frequency selection information " << mssFreqSel_p << " to AWP FTM." << LogIO::POST;
2748 93 : theFT->setSpwFreqSelection( mssFreqSel_p );
2749 93 : theIFT->setSpwFreqSelection( mssFreqSel_p );
2750 :
2751 :
2752 93 : }
2753 :
2754 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2755 :
2756 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2757 :
2758 :
2759 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2760 :
2761 208 : void SynthesisImagerVi2:: createMosFTMachine(CountedPtr<refim::FTMachine>& theFT,CountedPtr<refim::FTMachine>& theIFT, const Float /*padding*/, const Bool useAutoCorr, const Bool useDoublePrec, const Float rotatePAStep, const String stokes, const Bool doConjBeams){
2762 :
2763 416 : LogIO os(LogOrigin("SynthesisImagerVi2", "createMosFTMachine",WHERE));
2764 :
2765 208 : MSColumns msc(vi_p->ms());
2766 208 : String telescop=msc.observation().telescopeName()(0);
2767 208 : Bool multiTel=False;
2768 208 : Int msid=0;
2769 1737 : for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk()){
2770 1529 : if(((vi_p->getVisBuffer())->msId() != msid) && telescop != MSColumns(vi_p->ms()).observation().telescopeName()(0)){
2771 0 : msid=(vi_p->getVisBuffer())->msId();
2772 0 : multiTel=True;
2773 : }
2774 : }
2775 208 : vi_p->originChunks();
2776 :
2777 :
2778 :
2779 : PBMath::CommonPB kpb;
2780 208 : Record rec;
2781 208 : getVPRecord( rec, kpb, telescop );
2782 :
2783 :
2784 208 : if(rec.empty()){os << LogIO::SEVERE << "Cannot proceed with mosaicft gridder without a valid PB model" << LogIO::POST; }
2785 :
2786 : /*
2787 : VPManager *vpman=VPManager::Instance();
2788 : PBMath::CommonPB kpb;
2789 : PBMath::enumerateCommonPB(telescop, kpb);
2790 : Record rec;
2791 : vpman->getvp(rec, telescop);
2792 : */
2793 :
2794 208 : refim::VPSkyJones* vps= nullptr;
2795 : //cerr << "rec " << rec << " kpb " << kpb << endl;
2796 : //cerr << "createMOs ftname " << gridpars_p.ftmachine << endl;
2797 208 : if (!gridpars_p.ftmachine.contains("mos")) {
2798 0 : cerr << "PASTERP " << rotatePAStep << " " << gridpars_p.computePAStep << endl;
2799 0 : bool dosquint = (gridpars_p.computePAStep < 180); //anything beneath 180 deg ...you are not serious about squint correction
2800 : // TESTOO
2801 0 : dosquint = False;
2802 : ///////
2803 :
2804 0 : cerr << "Doing AWPLPG" << " wprojplanes " << gridpars_p.wprojplanes << endl;
2805 0 : theFT = new refim::AWPLPG(vps , gridpars_p.wprojplanes, dosquint, rotatePAStep*(C::pi)/180.0, mLocation_p, stokes, useAutoCorr, useDoublePrec, gridpars_p.usePointing);
2806 0 : theIFT = new refim::AWPLPG(vps , gridpars_p.wprojplanes, dosquint, rotatePAStep*(C::pi)/180.0, mLocation_p, stokes, useAutoCorr, useDoublePrec, gridpars_p.usePointing);
2807 0 : CountedPtr<refim::SimplePBConvFunc> mospb=new refim::HetArrayConvFunc();
2808 0 : static_cast<refim::AWPLPG &>(*theFT).setConvFunc(mospb);
2809 0 : static_cast<refim::AWPLPG &>(*theIFT).setConvFunc(mospb);
2810 :
2811 0 : }
2812 : else{
2813 208 : if(rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2814 191 : vps= new refim::VPSkyJones(msc, true, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
2815 : /////Don't know which parameter has pb threshold cutoff that the user want
2816 : ////leaving at default
2817 : ////vps.setThreshold(minPB);
2818 :
2819 : }
2820 : else{
2821 17 : PBMath myPB(rec);
2822 17 : String whichPBMath;
2823 17 : PBMathInterface::namePBClass(myPB.whichPBClass(), whichPBMath);
2824 17 : os << "Using the PB defined by " << whichPBMath << " for beam calculation for telescope " << telescop << LogIO::POST;
2825 17 : vps= new refim::VPSkyJones(telescop, myPB, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
2826 : //kpb=PBMath::DEFAULT;
2827 17 : }
2828 :
2829 208 : theFT = new refim::MosaicFTNew(vps, mLocation_p, stokes, 1000000000, 16, useAutoCorr,
2830 208 : useDoublePrec, doConjBeams, gridpars_p.usePointing);
2831 208 : PBMathInterface::PBClass pbtype=((kpb==PBMath::EVLA) || multiTel)? PBMathInterface::COMMONPB: PBMathInterface::AIRY;
2832 208 : if(rec.asString("name")=="IMAGE")
2833 0 : pbtype=PBMathInterface::IMAGE;
2834 : ///Use Heterogenous array mode for the following
2835 208 : if((kpb == PBMath::UNKNOWN) || (kpb==PBMath::OVRO) || (kpb==PBMath::ACA)
2836 195 : || (kpb==PBMath::ALMA) || (kpb==PBMath::EVLA) || multiTel){
2837 182 : CountedPtr<refim::SimplePBConvFunc> mospb=new refim::HetArrayConvFunc(pbtype, itsVpTable);
2838 182 : static_cast<refim::MosaicFTNew &>(*theFT).setConvFunc(mospb);
2839 182 : }
2840 : ///////////////////make sure both FTMachine share the same conv functions.
2841 208 : theIFT= new refim::MosaicFTNew(static_cast<refim::MosaicFTNew &>(*theFT));
2842 : }
2843 208 : }
2844 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2845 :
2846 144 : void SynthesisImagerVi2::createSDFTMachine(CountedPtr<refim::FTMachine>& theFT,
2847 : CountedPtr<refim::FTMachine>& theIFT,
2848 : const String &pointingDirCol,
2849 : const String &convertFirst,
2850 : const Float skyPosThreshold,
2851 : const Bool /*doPBCorr*/,
2852 : const Float rotatePAStep,
2853 : const String& gridFunction,
2854 : const Int convSupport,
2855 : const Quantity& truncateSize,
2856 : const Quantity& gwidth,
2857 : const Quantity& jwidth,
2858 : const Float minWeight,
2859 : const Bool clipMinMax,
2860 : const Int cache,
2861 : const Int tile,
2862 : const String &stokes,
2863 : const Bool pseudoI) {
2864 : // // member variable itsVPTable is VP table name
2865 288 : LogIO os(LogOrigin("SynthesisImagerVi2", "createSDFTMachine", WHERE));
2866 : os << LogIO::NORMAL // Loglevel INFO
2867 144 : << "Performing single dish gridding..." << LogIO::POST;
2868 : os << LogIO::NORMAL1 // gridFunction is too cryptic for most users.
2869 144 : << "with convolution function " << gridFunction << LogIO::POST;
2870 :
2871 : // Now make the Single Dish Gridding
2872 : os << LogIO::NORMAL // Loglevel INFO
2873 144 : << "Gridding will use specified common tangent point:" << LogIO::POST;
2874 144 : os << LogIO::NORMAL << SynthesisUtilMethods::asComprehensibleDirectionString(phaseCenter_p)
2875 144 : << LogIO::POST; // Loglevel INFO
2876 144 : String const gridfunclower = downcase(gridFunction);
2877 144 : if(gridfunclower=="pb") {
2878 22 : refim::SkyJones *vp = nullptr;
2879 22 : MSColumns msc(*mss_p[0]);
2880 : // todo: NONE is OK?
2881 22 : BeamSquint::SquintType squintType = BeamSquint::NONE;
2882 22 : Quantity skyPosThresholdQuant((Double)skyPosThreshold+180.0, "deg");
2883 22 : Quantity parAngleStepQuant((Double)rotatePAStep, "deg");
2884 22 : if (itsVpTable.empty()) {
2885 : os << LogIO::NORMAL // Loglevel INFO
2886 22 : << "Using defaults for primary beams used in gridding" << LogIO::POST;
2887 22 : vp=new refim::VPSkyJones(msc, true, parAngleStepQuant, squintType,
2888 22 : skyPosThresholdQuant);
2889 : } else {
2890 : os << LogIO::NORMAL // Loglevel INFO
2891 0 : << "Using VP as defined in " << itsVpTable << LogIO::POST;
2892 0 : Table vpTable( itsVpTable );
2893 0 : vp=new refim::VPSkyJones(msc, vpTable, parAngleStepQuant, squintType,
2894 0 : skyPosThresholdQuant);
2895 0 : }
2896 22 : theFT = new refim::SDGrid(mLocation_p, *vp, cache/2, tile, gridfunclower,
2897 22 : convSupport, minWeight, clipMinMax);
2898 22 : theIFT = new refim::SDGrid(mLocation_p, *vp, cache/2, tile, gridfunclower,
2899 22 : convSupport, minWeight, clipMinMax);
2900 144 : } else if (gridfunclower=="gauss" || gridfunclower=="gjinc") {
2901 3 : DirectionCoordinate dirCoord = itsMaxCoordSys.directionCoordinate();
2902 3 : Vector<String> units = dirCoord.worldAxisUnits();
2903 3 : Vector<Double> increments = dirCoord.increment();
2904 3 : Quantity cellx(increments[0], units[0]);
2905 3 : Quantity celly(increments[1], units[1]);
2906 6 : if (cellx != celly &&
2907 3 : ((!truncateSize.getUnit().empty()||truncateSize.getUnit()=="pixel")
2908 3 : || (!gwidth.getUnit().empty()||gwidth.getUnit()=="pixel")
2909 3 : || (!jwidth.getUnit().empty()||jwidth.getUnit()=="pixel"))) {
2910 : os << LogIO::WARN
2911 : << "The " << gridFunction << " gridding doesn't support non-square grid." << endl
2912 0 : << "Result may be wrong." << LogIO::POST;
2913 : }
2914 : Float truncateValue, gwidthValue, jwidthValue;
2915 3 : if (truncateSize.getUnit().empty() || truncateSize.getUnit()=="pixel")
2916 3 : truncateValue = truncateSize.getValue();
2917 : else
2918 0 : truncateValue = truncateSize.getValue("rad")/celly.getValue("rad");
2919 3 : if (gwidth.getUnit().empty() || gwidth.getUnit()=="pixel")
2920 3 : gwidthValue = gwidth.getValue();
2921 : else
2922 0 : gwidthValue = gwidth.getValue("rad")/celly.getValue("rad");
2923 3 : if (jwidth.getUnit().empty() || jwidth.getUnit()=="pixel")
2924 3 : jwidthValue = jwidth.getValue();
2925 : else
2926 0 : jwidthValue = jwidth.getValue("rad")/celly.getValue("rad");
2927 3 : theFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2928 3 : truncateValue, gwidthValue, jwidthValue, minWeight, clipMinMax);
2929 3 : theIFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2930 3 : truncateValue, gwidthValue, jwidthValue, minWeight, clipMinMax);
2931 3 : }
2932 : else {
2933 119 : theFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2934 119 : convSupport, minWeight, clipMinMax);
2935 119 : theIFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2936 119 : convSupport, minWeight, clipMinMax);
2937 : }
2938 144 : theFT->setPointingDirColumn(pointingDirCol);
2939 144 : static_cast<refim::SDGrid*>(theFT.get())->setConvertFirst(convertFirst);
2940 144 : static_cast<refim::SDGrid*>(theIFT.get())->setConvertFirst(convertFirst);
2941 :
2942 : // turn on Pseudo Stokes mode if necessary
2943 144 : if (pseudoI || stokes == "XX" || stokes == "YY" || stokes == "XXYY"
2944 288 : || stokes == "RR" || stokes == "LL" || stokes == "RRLL") {
2945 4 : theFT->setPseudoIStokes(True);
2946 4 : theIFT->setPseudoIStokes(True);
2947 : }
2948 144 : }
2949 :
2950 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2951 : //// Get/Set Weight Grid.... write to disk and read
2952 :
2953 : /// todo : do for full mapper list, and taylor terms.
2954 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2955 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2956 2 : Bool SynthesisImagerVi2::setWeightDensity( const String& weightimagename)
2957 : {
2958 4 : LogIO os(LogOrigin("SynthesisImagerVi2", "setWeightDensity()", WHERE));
2959 : try
2960 : {
2961 2 : if(weightimagename.size() !=0){
2962 2 : Table::isReadable(weightimagename, True);
2963 2 : PagedImage<Float> im(weightimagename);
2964 2 : imwgt_p=VisImagingWeight(im);
2965 2 : im.unlock();
2966 2 : }
2967 : else{
2968 : ////In memory weight densities is being deprecated...we should get rid of this bit
2969 0 : Int ndensities=1;
2970 0 : if((itsMappers.imageStore(0)->gridwt()->nelements())==5)
2971 0 : ndensities=(itsMappers.imageStore(0)->gridwt())->shape()[4];
2972 0 : Int nx=(itsMappers.imageStore(0)->gridwt())->shape()[0];
2973 0 : Int ny=(itsMappers.imageStore(0)->gridwt())->shape()[1];
2974 0 : Block<Matrix<Float> > densitymatrices(ndensities);
2975 0 : if(((itsMappers.imageStore(0)->gridwt())->shape().nelements())==5){
2976 0 : IPosition blc(Vector<Int>(5,0));
2977 0 : for (uInt fid=0;fid<densitymatrices.nelements();fid++)
2978 : {
2979 0 : densitymatrices[fid].resize();
2980 0 : Array<Float> lala;
2981 0 : blc[4]=fid;
2982 0 : itsMappers.imageStore(0)->gridwt()->getSlice(lala, blc, IPosition(5, nx, ny,1,1,1), True);
2983 0 : densitymatrices[fid].reference( lala.reform(IPosition(2, nx, ny)));
2984 0 : }
2985 0 : }
2986 : else{
2987 0 : Array<Float> lala;
2988 0 : itsMappers.imageStore(0)->gridwt()->get(lala, True);
2989 0 : densitymatrices[0].reference( lala.reform(IPosition(2, nx, ny)));
2990 0 : }
2991 0 : imwgt_p.setWeightDensity( densitymatrices );
2992 0 : }
2993 :
2994 2 : vi_p->useImagingWeight(imwgt_p);
2995 2 : itsMappers.releaseImageLocks();
2996 :
2997 : }
2998 0 : catch (AipsError &x)
2999 : {
3000 0 : throw(AipsError("In setWeightDensity : "+x.getMesg()));
3001 0 : }
3002 2 : return true;
3003 2 : }
3004 :
3005 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3006 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3007 :
3008 :
3009 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3010 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3011 1800 : void SynthesisImagerVi2::createVisSet(const Bool /*writeAccess*/)
3012 : {
3013 3600 : LogIO os( LogOrigin("SynthesisImagerVi2","createVisSet",WHERE) );
3014 : //cerr << "mss_p num" << mss_p.nelements() << " sel " << fselections_p->size() << endl;
3015 1800 : lockMSs();
3016 1800 : if(mss_p.nelements() > uInt(fselections_p->size()) && (fselections_p->size() !=0)){
3017 0 : throw(AipsError("Discrepancy between Number of MSs and Frequency selections"));
3018 : }
3019 1800 : vi_p=new vi::VisibilityIterator2(mss_p, vi::SortColumns(), true); //writeAccess);
3020 :
3021 1800 : if(fselections_p->size() !=0){
3022 1800 : CountedPtr<vi::FrequencySelections> tmpfselections=new FrequencySelections();
3023 : //Temporary fix till we get rid of old vi and we can get rid of tuneSelect
3024 1800 : if(uInt(fselections_p->size()) > mss_p.nelements()){
3025 0 : for(uInt k=0 ; k < mss_p.nelements(); ++k){
3026 0 : tmpfselections->add(fselections_p->get(k));
3027 : }
3028 : }
3029 : else{
3030 1800 : tmpfselections=fselections_p;
3031 : }
3032 : ////end of fix for tuneSelectdata
3033 1800 : vi_p->setFrequencySelection (*tmpfselections);
3034 :
3035 1800 : }
3036 : //
3037 1800 : vi_p->originChunks();
3038 1800 : vi_p->origin();
3039 : ////make sure to use the latest imaging weight scheme
3040 1800 : vi_p->useImagingWeight(imwgt_p);
3041 1800 : unlockMSs();
3042 1800 : }// end of createVisSet
3043 :
3044 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3045 : // Method to run the AWProjectFT machine to seperate the CFCache
3046 : // construction from imaging. This is done by splitting the
3047 : // operation in two steps: (1) run the FTM in "dry" mode to create a
3048 : // "blank" CFCache, and (2) re-load the "blank" CFCache and "fill"
3049 : // it.
3050 : //
3051 : // If someone can get me (SB) out of the horrible statc_casts in the
3052 : // code below, I will be most grateful (we are out of it! :-)).
3053 : //
3054 7 : void SynthesisImagerVi2::dryGridding(const Vector<String>& cfList)
3055 : {
3056 14 : LogIO os( LogOrigin("SynthesisImagerVi2","dryGridding",WHERE) );
3057 :
3058 7 : rownr_t cohDone=0;
3059 7 : Int whichFTM=0;
3060 : (void)cfList;
3061 : // If not an AWProject-class FTM, make this call a NoOp. Might be
3062 : // useful to extend it to other projection FTMs -- but later.
3063 7 : String ftmName = ((*(itsMappers.getFTM2(whichFTM)))).name();
3064 :
3065 7 : if (!((itsMappers.getFTM2(whichFTM,true))->isUsingCFCache())) return;
3066 :
3067 7 : os << "---------------------------------------------------- Dry Gridding ---------------------------------------------" << LogIO::POST;
3068 :
3069 : //
3070 : // Go through the entire MS in "dry" mode to set up a "blank"
3071 : // CFCache. This is done by setting the AWPWBFT in dryrun mode
3072 : // and gridding. The process of gridding emits CFCache, which
3073 : // will be "blank" in a dry run.
3074 : {
3075 7 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
3076 7 : vi_p->originChunks();
3077 7 : vi_p->origin();
3078 7 : Double numberCoh=0;
3079 14 : for (uInt k=0; k< mss_p.nelements(); ++k)
3080 7 : numberCoh+=Double(mss_p[k]->nrow());
3081 :
3082 14 : ProgressMeter pm(1.0, numberCoh, "dryGridding", "","","",true);
3083 :
3084 7 : itsMappers.initializeGrid(*vi_p);
3085 :
3086 : // Set the gridder (iFTM) to run in dry-gridding mode
3087 7 : (itsMappers.getFTM2(whichFTM,true))->setDryRun(true);
3088 :
3089 7 : Bool aTermIsOff=False;
3090 : {
3091 7 : CountedPtr<refim::FTMachine> ftm=itsMappers.getFTM2(whichFTM,True);
3092 7 : CountedPtr<refim::ConvolutionFunction> cf=ftm->getAWConvFunc();
3093 7 : aTermIsOff = cf->getTerm("ATerm")->isNoOp();
3094 7 : }
3095 :
3096 : os << "Making a \"blank\" CFCache"
3097 : << (aTermIsOff?" (without the A-Term)":"")
3098 7 : << LogIO::WARN << LogIO::POST;
3099 :
3100 : // Step through the MS. This triggers the logic in the Gridder
3101 : // to determine all the CFs that will be required. These empty
3102 : // CFs are written to the CFCache which can then be filled via
3103 : // a call to fillCFCache().
3104 31 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
3105 : {
3106 732 : for (vi_p->origin(); vi_p->more(); vi_p->next())
3107 : {
3108 708 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
3109 : {
3110 708 : itsMappers.grid(*vb, true, refim::FTMachine::OBSERVED, whichFTM);
3111 708 : cohDone += vb->nRows();
3112 708 : pm.update(Double(cohDone));
3113 : // If there is no term that depends on time, don't iterate over the entire data base
3114 708 : if (aTermIsOff) break;
3115 : }
3116 : }
3117 : }
3118 7 : }
3119 7 : if (cohDone == 0) os << "No valid rows found in dryGridding." << LogIO::EXCEPTION << LogIO::POST;
3120 : // Unset the dry-gridding mode.
3121 7 : (itsMappers.getFTM2(whichFTM,true))->setDryRun(false);
3122 :
3123 : //itsMappers.checkOverlappingModels("restore");
3124 7 : unlockMSs();
3125 : //fillCFCache(cfList);
3126 7 : }
3127 : //
3128 : // Re-load the CFCache from the disk using the supplied list of CFs
3129 : // (as cfList). Then extract the ConvFunc object (which was setup
3130 : // in the FTM) and call it's makeConvFunction2() to fill the CF.
3131 : // Finally, unset the dry-run mode of the FTM.
3132 : //
3133 7 : void SynthesisImagerVi2::fillCFCache(const Vector<String>& cfList,
3134 : const String& ftmName,
3135 : const String& cfcPath,
3136 : const Bool& psTermOn,
3137 : const Bool& aTermOn,
3138 : const Bool& conjBeams)
3139 : {
3140 14 : LogIO os( LogOrigin("SynthesisImagerVi2","fillCFCache",WHERE) );
3141 : // If not an AWProject-class FTM, make this call a NoOp. Might be
3142 : // useful to extend it to other projection FTMs -- but later.
3143 : // String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();
3144 :
3145 7 : if (!ftmName.contains("awproject") and
3146 0 : !ftmName.contains("multitermftnew")) return;
3147 : //if (!ftmName.contains("awproject")) return;
3148 :
3149 7 : os << "---------------------------------------------------- fillCFCache ---------------------------------------------" << LogIO::POST;
3150 :
3151 : //String cfcPath = itsMappers.getFTM(whichFTM)->getCacheDir();
3152 : //String imageNamePrefix=itsMappers.getFTM(whichFTM)->getCFCache()->getWtImagePrefix();
3153 :
3154 : //cerr << "Path = " << path << endl;
3155 :
3156 : // CountedPtr<AWProjectWBFTNew> tmpFT = new AWProjectWBFTNew(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM))));
3157 :
3158 :
3159 7 : Float dPA=360.0,selectedPA=2*360.0;
3160 7 : if (cfList.nelements() > 0)
3161 : {
3162 7 : CountedPtr<refim::CFCache> cfCacheObj = new refim::CFCache();
3163 : //Vector<String> wtCFList; wtCFList.resize(cfList.nelements());
3164 : //for (Int i=0; i<wtCFList.nelements(); i++) wtCFList[i] = "WT"+cfList[i];
3165 : //Directory dir(path);
3166 7 : Vector<String> cfList_p=cfList;//dir.find(Regex(Regex::fromPattern("CFS*")));
3167 7 : Vector<String> wtCFList_p;
3168 7 : wtCFList_p.resize(cfList_p.nelements());
3169 49 : for (Int i=0; i<(Int)wtCFList_p.nelements(); i++) wtCFList_p[i]="WT"+cfList_p[i];
3170 :
3171 : //cerr << cfList_p << endl;
3172 7 : cfCacheObj->setCacheDir(cfcPath.data());
3173 :
3174 7 : os << "Re-loading the \"blank\" CFCache for filling" << LogIO::WARN << LogIO::POST;
3175 :
3176 7 : cfCacheObj->initCacheFromList2(cfcPath, cfList_p, wtCFList_p,
3177 : selectedPA, dPA,CFC_VERBOSE);
3178 :
3179 : // tmpFT->setCFCache(cfCacheObj);
3180 7 : Vector<Double> uvScale, uvOffset;
3181 7 : Matrix<Double> vbFreqSelection;
3182 7 : CountedPtr<refim::CFStore2> cfs2 = CountedPtr<refim::CFStore2>(&cfCacheObj->memCache2_p[0],false);//new CFStore2;
3183 7 : CountedPtr<refim::CFStore2> cfwts2 = CountedPtr<refim::CFStore2>(&cfCacheObj->memCacheWt2_p[0],false);//new CFStore2;
3184 :
3185 : //
3186 : // Get whichFTM from itsMappers (SIMapperCollection) and
3187 : // cast it as AWProjectWBFTNew. Then get the ConvFunc from
3188 : // the FTM and cast it as AWConvFunc. Finally call
3189 : // AWConvFunc::makeConvFunction2().
3190 : //
3191 : // (static_cast<AWConvFunc &>
3192 : // (*(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getAWConvFunc())
3193 : // ).makeConvFunction2(String(path), uvScale, uvOffset, vbFreqSelection,
3194 : // *cfs2, *cfwts2);
3195 :
3196 : // This is a global methond in AWConvFunc. Does not require
3197 : // FTM to be constructed (which can be expensive in terms of
3198 : // memory footprint).
3199 7 : refim::AWConvFunc::makeConvFunction2(String(cfcPath), uvScale, uvOffset, vbFreqSelection,
3200 7 : *cfs2, *cfwts2, psTermOn, aTermOn, conjBeams);
3201 7 : }
3202 : //cerr << "Mem used = " << itsMappers.getFTM(whichFTM)->getCFCache()->memCache2_p[0].memUsage() << endl;
3203 : //(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getCFCache()->initCache2();
3204 7 : }
3205 :
3206 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3207 7 : void SynthesisImagerVi2::reloadCFCache()
3208 : {
3209 14 : LogIO os( LogOrigin("SynthesisImagerVi2","reloadCFCache",WHERE) );
3210 7 : Int whichFTM=0;
3211 7 : CountedPtr<refim::FTMachine> ftm=itsMappers.getFTM2(whichFTM,true);
3212 :
3213 : // Proceed only if FMTs uses the CFCache mechanism. The first FTM
3214 : // in the Mapper is used to make this decision. Not sure if the
3215 : // framework pipes allow other FTMs in SIMapper to be
3216 : // fundamentally different. If it does, and if that is
3217 : // triggered, the current decision may be insufficient.
3218 7 : if (!(ftm->isUsingCFCache())) return; // Better check than checking against FTM name
3219 :
3220 7 : os << "-------------------------------------------- Re-load CFCache ---------------------------------------------" << LogIO::POST;
3221 :
3222 : // Following code that distinguishes between MultiTermFTM and
3223 : // all others should ideally be replaced with a polymorphic
3224 : // solution. I.e. all FTMs should have a working getFTM2(bool)
3225 : // method. This is required since MultiTermFTM is a container
3226 : // FTM and it's getFTM2() returns the internal (per-MTMFS term)
3227 : // FTMs. Non-container FTMs must return a pointer to
3228 : // themselves. The if-else below is because attempts to make
3229 : // AWProjectFT::getFTM2() work have failed.
3230 : //
3231 : // Control reaches this stage only if the isUsingCFCache() test
3232 : // above return True. The only FTMs what will pass that test
3233 : // for now are AWProjectFT (and its derivatives) and
3234 : // MultiTermFTM if it is constructed with AWP.
3235 : //
3236 7 : CountedPtr<refim::CFCache> cfc;
3237 7 : if (ftm->name().contains("MultiTerm")) cfc = ftm->getFTM2(true)->getCFCache();
3238 7 : else cfc = ftm->getCFCache();
3239 7 : cfc->setLazyFill((refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1));
3240 7 : cfc->initCache2();
3241 :
3242 :
3243 : // String path,imageNamePrefix;
3244 : // if (ftm->name().contains("MultiTerm"))
3245 : // {
3246 : // path = ftm->getFTM2(true)->getCacheDir();
3247 : // imageNamePrefix = ftm->getFTM2(true)->getCFCache()->getWtImagePrefix();
3248 : // }
3249 : // else
3250 : // {
3251 : // path = ftm->getCacheDir();
3252 : // imageNamePrefix = ftm->getCFCache()->getWtImagePrefix();
3253 : // }
3254 :
3255 :
3256 : // CountedPtr<refim::CFCache> cfCacheObj = new refim::CFCache();
3257 : // cfCacheObj->setCacheDir(path.c_str());
3258 : // cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
3259 : // cfCacheObj->setLazyFill((refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1));
3260 : // cfCacheObj->initCache2();
3261 :
3262 : // // This assumes the itsMappers is always SIMapperCollection.
3263 : // for (whichFTM = 0; whichFTM < itsMappers.nMappers(); whichFTM++)
3264 : // {
3265 : // CountedPtr<refim::FTMachine> ifftm=itsMappers.getFTM2(whichFTM,true),
3266 : // fftm=itsMappers.getFTM2(whichFTM,false);
3267 :
3268 : // ifftm->setCFCache(cfCacheObj,true);
3269 : // fftm->setCFCache(cfCacheObj,true);
3270 : // }
3271 7 : }
3272 : //////////////////
3273 0 : bool SynthesisImagerVi2::makeMosaicSensitivity(){
3274 : ///We will bother with the first image. As A projection style gridding
3275 : ///usually is done on that first image.
3276 : /// if necessary in the future we will need to migrate this to SIMapper to
3277 : /// do it for all fields if multiple fields are A-projected.
3278 0 : if(!itsMappers.getFTM2(0))
3279 0 : return False;
3280 : /////////////////
3281 0 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
3282 0 : vi_p->originChunks();
3283 0 : vi_p->origin();
3284 0 : Double numcoh=0;
3285 0 : for (uInt k=0; k< mss_p.nelements(); ++k)
3286 0 : numcoh+=Double(mss_p[k]->nrow());
3287 : ProgressMeter pm(1.0, numcoh,
3288 0 : "Gridding Weights for PB", "","","",true);
3289 0 : rownr_t cohDone=0;
3290 :
3291 :
3292 : ///This will initialize weight grid too.
3293 0 : itsMappers.initializeGrid(*vi_p,True);
3294 0 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
3295 : {
3296 :
3297 0 : for (vi_p->origin(); vi_p->more(); vi_p->next())
3298 : {
3299 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
3300 : {
3301 0 : itsMappers.getFTM2(0)->gridImgWeights(*vb);
3302 0 : cohDone += vb->nRows();
3303 0 : pm.update(Double(cohDone));
3304 : }
3305 : }
3306 : }
3307 : //now load the images in weight and sumwt
3308 0 : itsMappers.getFTM2(0)-> finalizeToWeightImage(*vb, imageStore(0));
3309 : //cerr << "@@@@@@@MAKING PB " << endl;
3310 0 : return True;
3311 :
3312 :
3313 0 : }
3314 :
3315 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3316 676 : Bool SynthesisImagerVi2::loadMosaicSensitivity(){
3317 676 : String ftmname=itsMappers.getFTM2(0)->name();
3318 :
3319 676 : if(ftmname.contains("Mosaic") || ftmname.contains("AWProjectWB")){
3320 : //sumwt has been calcuated
3321 252 : Bool donesumwt=(max(itsMappers.imageStore(0)->sumwt()->get()) > 0.0);
3322 : //cerr << "Done sumwght " << donesumwt << max(itsMappers.imageStore(0)->sumwt()->get()) << endl;
3323 252 : if(donesumwt){
3324 504 : IPosition shp=itsMappers.imageStore(0)->weight()->shape();
3325 504 : CoordinateSystem cs=itsMappers.imageStore(0)->weight()->coordinates();
3326 252 : CountedPtr<TempImage<Float> > wgtim=new TempImage<Float>(shp, cs);
3327 252 : wgtim->copyData(*(itsMappers.imageStore(0)->weight()));
3328 252 : (static_cast<refim::FTMachine &>( *(itsMappers.getFTM2(0,False)))).setWeightImage(*wgtim);
3329 252 : static_cast<refim::FTMachine &>( *(itsMappers.getFTM2(0,True))).setWeightImage(*wgtim);
3330 252 : return true;
3331 252 : }
3332 :
3333 :
3334 : }
3335 424 : return false;
3336 676 : }
3337 : /////////////////////////////////////////////////
3338 0 : Record SynthesisImagerVi2::apparentSensitivity()
3339 : {
3340 0 : LogIO os(LogOrigin("SynthesisImagerVi2", "apparentSensitivity()", WHERE));
3341 :
3342 0 : Record outrec;
3343 : try {
3344 :
3345 : os << LogIO::NORMAL // Loglevel INFO
3346 : << "Calculating apparent sensitivity from MS weights, as modified by gridding weight function"
3347 0 : << LogIO::POST;
3348 : os << LogIO::NORMAL // Loglevel INFO
3349 : << "(assuming that MS weights have correct scale and units)"
3350 0 : << LogIO::POST;
3351 :
3352 0 : Double sumNatWt=0.0;
3353 0 : Double sumGridWt=0.0;
3354 0 : Double sumGridWt2OverNatWt=0.0;
3355 :
3356 0 : Float iNatWt(0.0),iGridWt(0.0);
3357 :
3358 0 : vi::VisBuffer2* vb = vi_p->getVisBuffer();
3359 0 : vi_p->originChunks();
3360 0 : vi_p->origin();
3361 :
3362 : // Discover if weightSpectrum non-trivially available
3363 0 : Bool doWtSp=vi_p->weightSpectrumExists();
3364 :
3365 : //////
3366 0 : for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
3367 : {
3368 0 : for (vi_p->origin(); vi_p->more(); vi_p->next())
3369 : {
3370 0 : auto nRow=vb->nRows();
3371 0 : const Vector<Bool>& rowFlags(vb->flagRow());
3372 :
3373 0 : const Vector<Int>& a1(vb->antenna1()), a2(vb->antenna2());
3374 :
3375 : // Extract weights correctly (use WEIGHT_SPECTRUM, if available)
3376 0 : Int nCorr=vb->nCorrelations();
3377 0 : Matrix<Float> wtm;
3378 0 : Cube<Float> wtc;
3379 0 : if (doWtSp)
3380 : // WS available [nCorr,nChan,nRow]
3381 0 : wtc.reference(vb->weightSpectrum());
3382 : else {
3383 : // WS UNavailable weight()[nCorr,nRow] --> [nCorr,nChan,nRow]
3384 0 : wtc.reference(vb->weight().reform(IPosition(3,nCorr,1,nRow))); // unchan'd weight as single-chan
3385 : }
3386 0 : Int nChanWt=wtc.shape()(1); // Might be 1 (no WtSp)
3387 :
3388 0 : Cube<Bool> flagCube(vb->flagCube());
3389 0 : for (rownr_t row=0; row<nRow; row++) {
3390 0 : if (!rowFlags(row) && a1(row)!=a2(row)) { // exclude ACs
3391 :
3392 0 : for (Int ich=0;ich<vb->nChannels();++ich) {
3393 0 : if( !flagCube(0,ich,row) && !flagCube(nCorr-1,ich,row)) { // p-hands unflagged
3394 :
3395 : // Accumulate relevant info
3396 :
3397 : // Simple sum of p-hand for now
3398 0 : iNatWt=wtc(0,ich%nChanWt,row)+wtc(nCorr-1,ich%nChanWt,row);
3399 :
3400 0 : iGridWt=2.0f*vb->imagingWeight()(ich,row);
3401 :
3402 0 : if (iGridWt>0.0 && iNatWt>0.0) {
3403 0 : sumNatWt+=(iNatWt);
3404 0 : sumGridWt+=(iGridWt);
3405 0 : sumGridWt2OverNatWt+=(iGridWt*iGridWt/iNatWt);
3406 : }
3407 : }
3408 : }
3409 : }
3410 : } // row
3411 0 : } // vb
3412 : } // chunks
3413 :
3414 0 : if (sumNatWt==0.0) {
3415 0 : os << "Cannot calculate sensitivity: sum of selected natural weights is zero" << LogIO::EXCEPTION;
3416 : }
3417 0 : if (sumGridWt==0.0) {
3418 0 : os << "Cannot calculate sensitivity: sum of gridded weights is zero" << LogIO::EXCEPTION;
3419 : }
3420 :
3421 0 : Double effSensitivity = sqrt(sumGridWt2OverNatWt)/sumGridWt;
3422 :
3423 0 : Double natSensitivity = 1.0/sqrt(sumNatWt);
3424 0 : Double relToNat=effSensitivity/natSensitivity;
3425 :
3426 : os << LogIO::NORMAL << "RMS Point source sensitivity : " // Loglevel INFO
3427 : << effSensitivity // << " Jy/beam" // actually, units are arbitrary
3428 0 : << LogIO::POST;
3429 : os << LogIO::NORMAL // Loglevel INFO
3430 0 : << "Relative to natural weighting : " << relToNat << LogIO::POST;
3431 :
3432 : // Fill output Record
3433 0 : outrec.define("relToNat",relToNat);
3434 0 : outrec.define("effSens",effSensitivity);
3435 :
3436 0 : } catch (AipsError x) {
3437 0 : throw(x);
3438 : return outrec;
3439 0 : }
3440 0 : return outrec;
3441 :
3442 0 : }
3443 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3444 :
3445 752 : Bool SynthesisImagerVi2::makePB()
3446 : {
3447 1504 : LogIO os( LogOrigin("SynthesisImagerVi2","makePB",WHERE) );
3448 :
3449 752 : if( itsMakeVP==False )
3450 : {
3451 173 : if( ((itsMappers.getFTM2(0))->name())!="MultiTermFTNew")
3452 133 : if(!loadMosaicSensitivity()){
3453 0 : if(!makeMosaicSensitivity())
3454 0 : throw(AipsError("Problem with making/loading sensitivity image for A -projection gridder"));
3455 : }
3456 : }
3457 : else
3458 : {
3459 579 : Bool doDefaultVP = itsVpTable.length()>0 ? False : True;
3460 :
3461 579 : CoordinateSystem coordsys=itsMappers.imageStore(0)->getCSys();
3462 579 : String telescope=coordsys.obsInfo().telescope();
3463 :
3464 579 : if (doDefaultVP) {
3465 :
3466 573 : MSAntennaColumns ac(mss_p[0]->antenna());
3467 573 : Double dishDiam=ac.dishDiameter()(0);
3468 573 : if(!allEQ(ac.dishDiameter().getColumn(), dishDiam))
3469 : os << LogIO::WARN
3470 : << "The MS has multiple antenna diameters ..PB could be wrong "
3471 0 : << LogIO::POST;
3472 573 : return makePBImage( telescope, False, dishDiam);
3473 573 : }
3474 : else{
3475 6 : return makePBImage(telescope );
3476 : }
3477 :
3478 579 : }
3479 :
3480 173 : return False;
3481 752 : }
3482 :
3483 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3484 :
3485 579 : Bool SynthesisImagerVi2::makePrimaryBeam(PBMath& pbMath)
3486 : {
3487 1158 : LogIO os( LogOrigin("SynthesisImagerVi2","makePrimaryBeam",WHERE) );
3488 :
3489 579 : os << "vi2 : Evaluating Primary Beam model onto image grid(s)" << LogIO::POST;
3490 :
3491 579 : itsMappers.initPB();
3492 :
3493 579 : vi::VisBuffer2* vb = vi_p->getVisBuffer();
3494 579 : vi_p->originChunks();
3495 579 : vi_p->origin();
3496 579 : std::map<Int, std::set<Int>> fieldsDone;
3497 579 : VisBufferUtil vbU(*vb);
3498 : ///////if tracking a moving source
3499 579 : MDirection origMovingDir;
3500 579 : MDirection newPhaseCenter;
3501 579 : Bool trackBeam=getMovingDirection(*vb, origMovingDir, True);
3502 : //////
3503 1307 : for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
3504 : {
3505 154098 : for (vi_p->origin(); vi_p->more(); vi_p->next())
3506 : {
3507 153370 : Bool fieldDone=False;
3508 153370 : if(fieldsDone.count(vb->msId() >0)){
3509 152768 : fieldDone=fieldDone || (fieldsDone[vb->msId()].count(vb->fieldId()(0)) > 0);
3510 : }
3511 : else{
3512 602 : fieldsDone[vb->msId()]=std::set<int>();
3513 : }
3514 153370 : if(!fieldDone){
3515 613 : fieldsDone[vb->msId()].insert(vb->fieldId()(0));
3516 613 : if(trackBeam){
3517 5 : MDirection newMovingDir;
3518 5 : getMovingDirection(*vb, newMovingDir);
3519 : //newPhaseCenter=vb->phaseCenter();
3520 5 : newPhaseCenter=vbU.getPhaseCenter(*vb);
3521 5 : newPhaseCenter.shift(MVDirection(-newMovingDir.getAngle()+origMovingDir.getAngle()), False);
3522 5 : }
3523 613 : itsMappers.addPB(*vb,pbMath, newPhaseCenter, trackBeam);
3524 :
3525 : }
3526 : }
3527 : }
3528 579 : itsMappers.releaseImageLocks();
3529 579 : unlockMSs();
3530 :
3531 579 : return True;
3532 579 : }// end makePB
3533 :
3534 584 : Bool SynthesisImagerVi2::getMovingDirection(const vi::VisBuffer2& vb, MDirection& outDir, const Bool useImageEpoch){
3535 584 : MDirection movingDir;
3536 584 : Bool trackBeam=False;
3537 :
3538 1168 : MeasFrame mFrame(MEpoch(Quantity(vb.time()(0), "s"), MSColumns(vb.ms()).timeMeas()(0).getRef()), mLocation_p);
3539 584 : if(useImageEpoch){
3540 579 : mFrame.resetEpoch((itsMappers.imageStore(0))->getCSys().obsInfo().obsDate());
3541 :
3542 : }
3543 584 : if(movingSource_p != ""){
3544 : MDirection::Types refType;
3545 10 : trackBeam=True;
3546 10 : if(Table::isReadable(movingSource_p, False)){
3547 : //seems to be a table so assuming ephemerides table
3548 6 : Table laTable(movingSource_p);
3549 6 : Path leSentier(movingSource_p);
3550 6 : MeasComet laComet(laTable, leSentier.absoluteName());
3551 6 : movingDir.setRefString("COMET");
3552 6 : mFrame.set(laComet);
3553 6 : }
3554 : ///if not a table
3555 4 : else if(casacore::MDirection::getType(refType, movingSource_p)){
3556 0 : if(refType > casacore::MDirection::N_Types && refType < casacore::MDirection:: N_Planets ){
3557 : ///A known planet
3558 0 : movingDir.setRefString(movingSource_p);
3559 : }
3560 : }
3561 4 : else if(upcase(movingSource_p)=="TRACKFIELD"){
3562 4 : VisBufferUtil vbU(vb);
3563 4 : movingDir=vbU.getEphemDir(vb, MEpoch(mFrame.epoch()).get("s").getValue());
3564 4 : }
3565 : else{
3566 0 : throw(AipsError("Erroneous tracking direction set to make pb"));
3567 : }
3568 10 : MDirection::Ref outref1(MDirection::AZEL, mFrame);
3569 10 : MDirection tmphazel=MDirection::Convert(movingDir, outref1)();
3570 10 : MDirection::Ref outref(vb.phaseCenter().getRef().getType(), mFrame);
3571 10 : outDir=MDirection::Convert(tmphazel, outref)();
3572 10 : }
3573 : else{
3574 574 : outDir=vb.phaseCenter();
3575 574 : trackBeam=False;
3576 : }
3577 584 : return trackBeam;
3578 :
3579 :
3580 584 : }
3581 1 : CountedPtr<vi::VisibilityIterator2> SynthesisImagerVi2::getVi(){
3582 1 : return vi_p;
3583 : }
3584 1 : CountedPtr<refim::FTMachine> SynthesisImagerVi2::getFTM(const Int fid, Bool ift){
3585 1 : if(itsMappers.nMappers() <=fid)
3586 0 : throw(AipsError("Mappers have not been set for field "+String::toString(fid)));
3587 1 : return (itsMappers.getFTM2(fid, ift));
3588 :
3589 : }
3590 284 : void SynthesisImagerVi2::updateImageBeamSet(Record& returnRec){
3591 284 : if(returnRec.isDefined("imageid") && returnRec.asInt("imageid")==0){
3592 : //ImageInfo ii=(itsMappers.imageStore(0)->psf())->imageInfo();
3593 283 : Vector<Int> chanRange(0);
3594 283 : if(returnRec.isDefined("chanrange"))
3595 283 : chanRange=returnRec.asArrayInt("chanrange");
3596 283 : Int npol=(itsMappers.imageStore(0)->getShape())(2);
3597 283 : Int nchan=(itsMappers.imageStore(0)->getShape())(3);
3598 283 : if(chanRange.nelements() ==2 && chanRange(1) < nchan){
3599 :
3600 283 : ImageBeamSet iibeamset=cubePsfImageInfo_p.getBeamSet();
3601 283 : Matrix<GaussianBeam> mbeams=iibeamset.getBeams();
3602 283 : if(mbeams.nelements()==0){
3603 262 : mbeams.resize(itsMappers.imageStore(0)->getShape()(3), itsMappers.imageStore(0)->getShape()(2));
3604 262 : mbeams.set(GaussianBeam(Vector<Quantity>(3, Quantity(1e-12, "arcsec"))));
3605 : }
3606 283 : Cube<Float> recbeams(returnRec.asArrayFloat("beams"));
3607 2549 : for(Int c=chanRange[0]; c <= chanRange[1]; ++c){
3608 4703 : for(Int p=0; p < npol; ++p){
3609 2437 : mbeams(c,p)=GaussianBeam(Quantity(recbeams(c-chanRange[0], p, 0),"arcsec"), Quantity(recbeams(c-chanRange[0], p, 1),"arcsec"), Quantity(recbeams(c-chanRange[0], p, 2),"deg"));
3610 :
3611 : }
3612 : }
3613 283 : cubePsfImageInfo_p.setBeams(ImageBeamSet(mbeams));
3614 :
3615 283 : }
3616 : //itsMappers.imageStore(0)->psf()->setImageInfo(ii);
3617 : //itsMappers.imageStore(0)->psf()->unlock();
3618 :
3619 :
3620 283 : }
3621 284 : }
3622 :
3623 : } //# NAMESPACE CASA - END
3624 :
|
