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 0 : 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 0 : mss_p.resize(0);
121 0 : }
122 :
123 0 : SynthesisImagerVi2::~SynthesisImagerVi2(){
124 0 : for (uInt k=0; k < mss_p.nelements(); ++k){
125 0 : if(mss_p[k])
126 0 : delete mss_p[k];
127 : }
128 0 : SynthesisUtilMethods::getResource("End Run");
129 0 : }
130 :
131 0 : Bool SynthesisImagerVi2::selectData(const SynthesisParamsSelect& selpars){
132 0 : LogIO os( LogOrigin("SynthesisImagerVi2","selectData",WHERE) );
133 0 : Bool retval=True;
134 :
135 0 : SynthesisUtilMethods::getResource("Start Run");
136 :
137 : try
138 : {
139 :
140 :
141 0 : 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 0 : uInt exceptCounter=0;
145 :
146 : while(true){
147 : try{
148 : //Respect the readonly flag...necessary for multi-process access
149 0 : thisms=MeasurementSet(selpars.msname, TableLock(TableLock::UserNoReadLocking),
150 0 : selpars.readonly ? Table::Old : Table::Update);
151 0 : 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 0 : useScratch_p=selpars.usescratch;
176 0 : readOnly_p = selpars.readonly;
177 0 : lockMS(thisms);
178 0 : 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 0 : if(!(impars_p.mode.contains("cube")) || ((impars_p.mode.contains("cube")) && doingCubeGridding_p)){
185 0 : if(selpars.usescratch && !selpars.readonly){
186 0 : VisSetUtil::addScrCols(thisms, true, false, true, false);
187 0 : 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 0 : if(!selpars.incrmodel && !selpars.usescratch && !selpars.readonly)
194 0 : refim::VisModelData::clearModel(thisms, selpars.field, selpars.spw);
195 : }//end of first pass if for cubes
196 : // unlockMSs();
197 :
198 :
199 0 : 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 0 : MSSelection thisSelection;
205 0 : if(selpars.field != ""){
206 0 : thisSelection.setFieldExpr(selpars.field);
207 0 : os << "Selecting on fields : " << selpars.field << " | " ;//LogIO::POST;
208 : }else
209 0 : thisSelection.setFieldExpr("*");
210 0 : if(selpars.spw != ""){
211 0 : thisSelection.setSpwExpr(selpars.spw);
212 0 : os << "Selecting on spw :"<< selpars.spw << " | " ;//LogIO::POST;
213 : }else
214 0 : thisSelection.setSpwExpr("*");
215 :
216 0 : if(selpars.antenna != ""){
217 0 : Vector<String> antNames(1, selpars.antenna);
218 : // thisSelection.setAntennaExpr(MSSelection::nameExprStr( antNames));
219 0 : thisSelection.setAntennaExpr(selpars.antenna);
220 0 : os << "Selecting on antenna names : " << selpars.antenna << " | " ;//LogIO::POST;
221 :
222 0 : }
223 0 : if(selpars.timestr != ""){
224 0 : thisSelection.setTimeExpr(selpars.timestr);
225 0 : os << "Selecting on time range : " << selpars.timestr << " | " ;//LogIO::POST;
226 : }
227 0 : if(selpars.uvdist != ""){
228 0 : thisSelection.setUvDistExpr(selpars.uvdist);
229 0 : os << "Selecting on uvdist : " << selpars.uvdist << " | " ;//LogIO::POST;
230 : }
231 0 : if(selpars.scan != ""){
232 0 : thisSelection.setScanExpr(selpars.scan);
233 0 : os << "Selecting on scan : " << selpars.scan << " | " ;//LogIO::POST;
234 : }
235 0 : if(selpars.obs != ""){
236 0 : thisSelection.setObservationExpr(selpars.obs);
237 0 : os << "Selecting on Observation Expr : " << selpars.obs << " | " ;//LogIO::POST;
238 : }
239 0 : if(selpars.state != ""){
240 0 : thisSelection.setStateExpr(selpars.state);
241 0 : 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 0 : os << "[Opened " << (readOnly_p?"in readonly mode":(useScratch_p?"with scratch model column":"with virtual model column")) << "]" << LogIO::POST;
248 0 : TableExprNode exprNode=thisSelection.toTableExprNode(&thisms);
249 0 : if(!(exprNode.isNull()))
250 : {
251 :
252 :
253 0 : MeasurementSet thisMSSelected0 = MeasurementSet(thisms(exprNode));
254 0 : mss_p.resize(mss_p.nelements()+1, false, true);
255 0 : 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 0 : mss_p[mss_p.nelements()-1]=new MeasurementSet(thisMSSelected0);
268 :
269 0 : os << " NRows selected : " << (mss_p[mss_p.nelements()-1])->nrow() << LogIO::POST;
270 0 : unlockMSs();
271 0 : }
272 : else{
273 0 : throw(AipsError("Selection for given MS "+selpars.msname+" is invalid"));
274 : }
275 0 : if((mss_p[mss_p.nelements()-1])->nrow() ==0){
276 0 : delete mss_p[mss_p.nelements()-1];
277 0 : mss_p.resize(mss_p.nelements()-1, True, True);
278 0 : if(mss_p.nelements()==0)
279 0 : 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 0 : if(!fselections_p) fselections_p=new FrequencySelections();
291 0 : Matrix<Int> chanlist = thisSelection.getChanList(mss_p[mss_p.nelements()-1]);
292 0 : Matrix<Double> freqList=thisSelection.getChanFreqList(mss_p[mss_p.nelements()-1]);
293 : //cerr << std::setprecision(12) << "FreqList " << freqList << endl;
294 0 : IPosition shape = freqList.shape();
295 0 : uInt nSelections = shape[0];
296 : ///temporary variable as we carry that for tunechunk...till we get rid of it
297 0 : selFreqFrame_p=selpars.freqframe;
298 0 : Bool ignoreframe=False;
299 0 : MFrequency::Types freqFrame=MFrequency::castType(MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().measFreqRef()(Int(chanlist(0,0))));
300 :
301 0 : if(selpars.freqframe == MFrequency::REST ||selpars.freqframe == MFrequency::Undefined){
302 0 : selFreqFrame_p=freqFrame;
303 0 : ignoreframe=True;
304 : }
305 0 : 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 0 : 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 0 : vi::FrequencySelectionUsingFrame channelSelector(selFreqFrame_p);
321 : ///temporary variable as we carry that for tunechunk
322 :
323 0 : Bool selectionValid=False;
324 0 : for(uInt k=0; k < nSelections; ++k){
325 0 : Bool thisSpwSelValid=False;
326 : //The getChanfreqList is wrong for beg and end..going round that too.
327 0 : Vector<Double> freqies=MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().chanFreq()(Int(chanlist(k,0)));
328 0 : Vector<Double> chanwidth=MSColumns(*mss_p[mss_p.nelements()-1]).spectralWindow().chanWidth()(Int(chanlist(k,0)));
329 :
330 0 : if(freqList(k,3) < 0.0){
331 0 : topfreq=freqies(chanlist(k,1));//-chanwidth(chanlist(k,1))/2.0;
332 0 : 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 0 : lowfreq=freqies(chanlist(k,1));//-chanwidth(chanlist(k,1))/2.0;
338 0 : 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 0 : 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 0 : if(MSUtil::getFreqRangeInSpw( lowfreq,
351 0 : topfreq, Vector<Int>(1,chanlist(k,0)), Vector<Int>(1,chanlist(k,1)),
352 0 : Vector<Int>(1, chanlist(k,2)-chanlist(k,1)+1),
353 0 : *mss_p[mss_p.nelements()-1] ,
354 : selFreqFrame_p,
355 : fieldList, False))
356 : {
357 0 : selectionValid=True;
358 0 : thisSpwSelValid=True;
359 : }
360 : // unlockMSs();
361 :
362 : }
363 :
364 0 : if(thisSpwSelValid || ignoreframe){
365 0 : andFreqSelection(mss_p.nelements()-1, Int(freqList(k,0)), lowfreq, topfreq, selFreqFrame_p);
366 0 : andChanSelection(mss_p.nelements()-1, Int(chanlist(k,0)), Int(chanlist(k,1)),Int(chanlist(k,2)));
367 : }
368 0 : }
369 0 : if(! (selectionValid && !ignoreframe)){
370 : //os << "Did not match spw selection in the selected ms " << LogIO::WARN << LogIO::POST;
371 0 : retval=False;
372 : }
373 : //fselections_p->add(channelSelector);
374 : //////////////////////////////////
375 0 : }
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 0 : Quantity freq;
381 0 : Quantity::read(freq, selpars.freqbeg);
382 0 : Double lowfreq=freq.getValue("Hz");
383 0 : Quantity::read(freq, selpars.freqend);
384 0 : 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 0 : 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 0 : andFreqSelection((mss_p.nelements()-1), Int(freqList(k,0)), lowfreq, topfreq, selFreqFrame_p);
396 : }
397 : //fselections_p->add(channelSelector);
398 :
399 0 : }
400 0 : }//End of channel selection
401 :
402 0 : writeAccess_p=writeAccess_p && !selpars.readonly;
403 0 : createVisSet(writeAccess_p);
404 :
405 : /////// Remove this when the new vi/vb is able to get the full freq range.
406 0 : mssFreqSel_p.resize();
407 0 : 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 0 : if( selpars.datacolumn.contains("data") || selpars.datacolumn.contains("obs") )
413 0 : { if( thisms.tableDesc().isColumn("DATA") ) { datacol_p = FTMachine::OBSERVED; }
414 0 : else { os << LogIO::SEVERE <<"DATA column does not exist" << LogIO::EXCEPTION;}
415 : }
416 0 : 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 0 : dataSel_p.resize(dataSel_p.nelements()+1, true);
421 0 : dataSel_p[dataSel_p.nelements()-1]=selpars;
422 : //cerr << "AT THE end of DATASEL " << selpars.toRecord() << endl;
423 :
424 0 : unlockMSs();
425 0 : }
426 0 : catch(AipsError &x)
427 : {
428 0 : unlockMSs();
429 0 : throw( AipsError("Error in selectData() : "+x.getMesg()) );
430 0 : }
431 :
432 0 : return retval;
433 :
434 :
435 :
436 0 : }
437 0 : void SynthesisImagerVi2::andChanSelection(const Int msId, const Int spwId, const Int startchan, const Int endchan){
438 :
439 0 : map<Int, Vector<Int> > spwsel;
440 0 : auto it=channelSelections_p.find(msId);
441 0 : if(it !=channelSelections_p.end())
442 0 : spwsel=it->second;
443 0 : auto hasspw=spwsel.find(spwId);
444 0 : Vector<Int>chansel(2,-1);
445 0 : if(hasspw != spwsel.end()){
446 0 : chansel.resize();
447 0 : chansel=hasspw->second;
448 : }
449 0 : Int nchan=endchan-startchan+1;
450 0 : if(chansel(1)== -1)
451 0 : chansel(1)=startchan;
452 0 : if(chansel(1) >= startchan){
453 0 : if(nchan > (chansel(1)-startchan+chansel(0))){
454 0 : chansel(0)=nchan;
455 : }
456 : else{
457 0 : chansel(0)=chansel(1)-startchan+chansel(0);
458 : }
459 0 : chansel(1)=startchan;
460 : }
461 : else{
462 0 : if((chansel(0) -(startchan - chansel(1)+1)) < nchan){
463 0 : chansel(0)=nchan+(startchan-chansel(1));
464 : }
465 : }
466 0 : spwsel[spwId]=chansel;
467 0 : channelSelections_p[msId]=spwsel;
468 : //cerr << "chansel "<< channelSelections_p << endl;
469 :
470 0 : }
471 0 : void SynthesisImagerVi2::andFreqSelection(const Int msId, const Int spwId, const Double freqBeg, const Double freqEnd, const MFrequency::Types frame){
472 :
473 :
474 0 : Int key=msId;
475 :
476 0 : Bool isDefined=False;
477 0 : FrequencySelectionUsingFrame frameSel(frame);
478 0 : 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 0 : 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 0 : isDefined=True;
482 : //cerr << k << " inside freqBegs " << freqBegs_p[k].second << " " << freqBeg << endl;
483 0 : if(freqBegs_p[k].second < freqBeg)
484 0 : freqBegs_p[k].second=freqBeg;
485 0 : if(freqEnds_p[k].second > freqEnd)
486 0 : freqEnds_p[k].second=freqEnd;
487 0 : 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 0 : if(freqBegs_p[k].first == key){
492 : //cerr << "added " << k << " freqBegs " << freqBegs_p[k].second << " " << freqEnds_p[k].second << endl;
493 0 : frameSel.add(freqSpws_p[k].second , freqBegs_p[k].second, freqEnds_p[k].second);
494 : }
495 : }
496 0 : if(!isDefined && msId >=0){
497 : //cerr << "undefined " << key << " freqBegs " << freqBeg << " " << freqEnd << endl;
498 0 : freqBegs_p.push_back(make_pair(key, freqBeg));
499 0 : freqEnds_p.push_back(make_pair(key, freqEnd));
500 0 : freqSpws_p.push_back(make_pair(key, spwId));
501 0 : frameSel.add(spwId, freqBeg, freqEnd);
502 : }
503 0 : CountedPtr<vi::FrequencySelections> copyFsels=fselections_p->clone();
504 0 : uInt nMSs=copyFsels->size() <=msId ? msId+1 : copyFsels->size();
505 : //cerr << "nms " << nMSs << endl;
506 0 : fselections_p=new FrequencySelections();
507 0 : for (uInt k=0; k < nMSs ; ++k){
508 0 : if(k==uInt(key)){
509 0 : fselections_p->add(frameSel);
510 : //cerr <<"adding framesel " << frameSel.toString() << endl;
511 : }
512 : else{
513 0 : const FrequencySelectionUsingFrame& thissel= static_cast<const FrequencySelectionUsingFrame &> (copyFsels->get(k));
514 : //cerr <<"framesel orig " << thissel.toString() << endl;
515 0 : fselections_p->add(thissel);
516 :
517 : }
518 : }
519 :
520 :
521 :
522 0 : }
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 0 : Bool SynthesisImagerVi2::defineImage(
581 : SynthesisParamsImage& impars,
582 : const SynthesisParamsGrid& gridpars)
583 : {
584 0 : LogIO os( LogOrigin("SynthesisImagerVi2", "defineImage", WHERE) );
585 :
586 0 : if (mss_p.nelements() == 0)
587 0 : os << "SelectData has to be run before defineImage" << LogIO::EXCEPTION;
588 :
589 0 : CoordinateSystem csys;
590 0 : CountedPtr<refim::FTMachine> ftm, iftm;
591 0 : impars_p = impars;
592 0 : gridpars_p = gridpars;
593 :
594 : try {
595 0 : os << "Define image coordinates for [" << impars.imageName << "] : "
596 0 : << LogIO::POST;
597 0 : 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 0 : csys = impars_p.buildCoordinateSystem( *vi_p, channelSelections_p, mss_p );
601 : //use the location defined for coordinates frame;
602 0 : mLocation_p=impars_p.obslocation;
603 0 : IPosition imshape = impars_p.shp();
604 :
605 :
606 0 : os << "Impars: start " << impars_p.start << LogIO::POST;
607 : os << "Shape: " << imshape
608 0 : << " Spectral: " << csys.spectralCoordinate().referenceValue()
609 0 : << " at " << csys.spectralCoordinate().referencePixel()
610 0 : << " with increment " << csys.spectralCoordinate().increment()
611 0 : << LogIO::POST;
612 :
613 0 : if( (itsMappers.nMappers()==0) ||
614 0 : (impars_p.imsize[0]*impars_p.imsize[1] > itsMaxShape[0]*itsMaxShape[1]))
615 : {
616 0 : itsMaxShape=imshape;
617 0 : itsMaxCoordSys=csys;
618 : }
619 0 : itsNchan = imshape[3];
620 0 : 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 0 : if (impars_p.phaseCenterFieldId == -1) {
632 : // user-specified
633 0 : phaseCenter_p = impars_p.phaseCenter;
634 0 : } else if (impars_p.phaseCenterFieldId >= 0) {
635 : // FIELD_ID
636 0 : auto const msobj = mss_p[0];
637 0 : MSFieldColumns msfield(msobj->field());
638 0 : phaseCenter_p=msfield.phaseDirMeas(impars_p.phaseCenterFieldId);
639 0 : } else {
640 : // use default FIELD_ID (0)
641 0 : auto const msobj = mss_p[0];
642 0 : MSFieldColumns msfield(msobj->field());
643 0 : phaseCenter_p=msfield.phaseDirMeas(0);
644 0 : }
645 :
646 :
647 0 : if ( (itsMappers.nMappers() == 0) or
648 0 : (impars_p.imsize[0]*impars_p.imsize[1] > itsMaxShape[0]*itsMaxShape[1])
649 : ) {
650 0 : itsMaxShape = imshape;
651 0 : itsMaxCoordSys = csys;
652 : }
653 0 : itsNchan = imshape[3];
654 0 : 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 0 : if (impars_p.phaseCenterFieldId == -1) { // user-specified
666 0 : phaseCenter_p = impars_p.phaseCenter;
667 0 : } else if (impars_p.phaseCenterFieldId >= 0) { // FIELD_ID
668 0 : auto const msobj = mss_p[0];
669 0 : MSFieldColumns msfield(msobj->field());
670 0 : phaseCenter_p = msfield.phaseDirMeas(impars_p.phaseCenterFieldId);
671 0 : } else { // use default FIELD_ID (0)
672 0 : auto const msobj = mss_p[0];
673 0 : MSFieldColumns msfield(msobj->field());
674 0 : phaseCenter_p = msfield.phaseDirMeas(0);
675 0 : }
676 0 : }
677 0 : catch (AipsError &x) {
678 : os << "Error in building Coordinate System and Image Shape: "
679 : << x.getMesg()
680 0 : << LogIO::EXCEPTION;
681 0 : }
682 :
683 : try {
684 0 : os << "Set Gridding options for [" << impars_p.imageName << "]"
685 0 : << " with ftmachine: " << gridpars.ftmachine
686 0 : << LogIO::POST;
687 :
688 0 : itsVpTable = gridpars.vpTable;
689 0 : itsMakeVP = ( gridpars.ftmachine.contains("mosaicft") or
690 0 : gridpars.ftmachine.contains("awprojectft") ) ?
691 : False : True;
692 :
693 : //cerr << "DEFINEimage " << impars_p.toRecord() << endl;
694 :
695 0 : createFTMachine(
696 0 : ftm, iftm, gridpars.ftmachine, impars_p.nTaylorTerms, gridpars.mType,
697 0 : gridpars.facets, gridpars.wprojplanes,
698 0 : gridpars.padding, gridpars.useAutoCorr, gridpars.useDoublePrec,
699 0 : gridpars.convFunc,
700 0 : gridpars.aTermOn, gridpars.psTermOn, gridpars.mTermOn,
701 0 : gridpars.wbAWP, gridpars.cfCache, gridpars.usePointing, gridpars.pointingOffsetSigDev.tovector(),
702 0 : gridpars.doPBCorr, gridpars.conjBeams,
703 0 : gridpars.computePAStep, gridpars.rotatePAStep,
704 0 : gridpars.interpolation, impars_p.freqFrameValid, 1000000000, 16, impars_p.stokes,
705 0 : impars_p.imageName,
706 0 : gridpars.pointingDirCol, gridpars.convertFirst, gridpars.skyPosThreshold,
707 0 : gridpars.convSupport, gridpars.truncateSize, gridpars.gwidth, gridpars.jwidth,
708 0 : 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 0 : appendToMapperList(
719 0 : impars_p.imageName, csys, impars_p.shp(),
720 : ftm, iftm,
721 0 : gridpars.distance, gridpars.facets, gridpars.chanchunks, impars_p.overwrite,
722 0 : gridpars.mType, gridpars.padding, impars_p.nTaylorTerms, impars_p.startModel
723 : );
724 :
725 0 : imageDefined_p = true;
726 : }
727 0 : catch(AipsError &x) {
728 0 : os << "Error in adding Mapper: " + x.getMesg() << LogIO::EXCEPTION;
729 0 : }
730 :
731 0 : imparsVec_p.resize(imparsVec_p.nelements()+1, true);
732 0 : imparsVec_p[imparsVec_p.nelements()-1] = impars_p;
733 : ///For now cannot deal with cube and mtmfs in C++ parallel mode
734 0 : if (imparsVec_p[0].deconvolver == "mtmfs") setCubeGridding(False);
735 : // cerr << "DECONV " << imparsVec_p[0].deconvolver
736 : // << " cube gridding " << doingCubeGridding_p << endl;
737 0 : gridparsVec_p.resize(gridparsVec_p.nelements()+1, true);
738 0 : 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 0 : itsMakeVP =
744 0 : ( gridparsVec_p[0].ftmachine.contains("mosaicft") or
745 0 : (gridparsVec_p[0].ftmachine.at(0,3) == "awp")
746 0 : ) ? False : True;
747 0 : return true;
748 0 : }
749 :
750 0 : Bool SynthesisImagerVi2::defineImage(
751 : CountedPtr<SIImageStore> imstor,
752 : SynthesisParamsImage& impars,
753 : const SynthesisParamsGrid& gridpars)
754 : {
755 0 : gridpars_p=gridpars;
756 0 : Int id = itsMappers.nMappers();
757 0 : CoordinateSystem csys = imstor->getCSys();
758 0 : IPosition imshape = imstor->getShape();
759 0 : Int nx=imshape[0], ny=imshape[1];
760 0 : if ( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]) ) {
761 0 : itsMaxShape=imshape;
762 0 : itsMaxCoordSys=csys;
763 : }
764 0 : mLocation_p=impars.obslocation;
765 : // phasecenter
766 0 : if (impars.phaseCenterFieldId == -1) { // user-specified
767 0 : phaseCenter_p = impars.phaseCenter;
768 : }
769 0 : else if (impars.phaseCenterFieldId >= 0) { // FIELD_ID
770 0 : auto const msobj = mss_p[0];
771 0 : MSFieldColumns msfield(msobj->field());
772 0 : phaseCenter_p = msfield.phaseDirMeas(impars.phaseCenterFieldId);
773 0 : }
774 : else { // use default FIELD_ID (0)
775 0 : auto const msobj = mss_p[0];
776 0 : MSFieldColumns msfield(msobj->field());
777 0 : phaseCenter_p = msfield.phaseDirMeas(0);
778 0 : }
779 :
780 0 : itsVpTable = gridpars.vpTable;
781 0 : itsMakeVP = ( gridpars.ftmachine.contains("mosaicft") or
782 0 : (gridpars.ftmachine.at(0,3) == "awp")
783 0 : ) ? False : True;
784 0 : CountedPtr<refim::FTMachine> ftm, iftm;
785 :
786 0 : createFTMachine(
787 0 : ftm, iftm, gridpars.ftmachine, impars.nTaylorTerms, gridpars.mType,
788 0 : gridpars.facets, gridpars.wprojplanes,
789 0 : gridpars.padding,gridpars.useAutoCorr,gridpars.useDoublePrec,
790 0 : gridpars.convFunc,
791 0 : gridpars.aTermOn,gridpars.psTermOn, gridpars.mTermOn,
792 0 : gridpars.wbAWP,gridpars.cfCache,gridpars.usePointing,
793 0 : gridpars.pointingOffsetSigDev.tovector(),
794 0 : gridpars.doPBCorr,gridpars.conjBeams,
795 0 : gridpars.computePAStep,gridpars.rotatePAStep,
796 0 : gridpars.interpolation, impars.freqFrameValid, 1000000000, 16, impars.stokes,
797 0 : impars.imageName,
798 0 : gridpars.pointingDirCol, gridpars.convertFirst, gridpars.skyPosThreshold,
799 0 : gridpars.convSupport, gridpars.truncateSize, gridpars.gwidth, gridpars.jwidth,
800 0 : gridpars.minWeight, gridpars.clipMinMax, impars.pseudoi);
801 :
802 0 : 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 0 : itsMappers.addMapper(
823 0 : createSIMapper( gridpars.mType, imstor, ftm, iftm)
824 : );
825 : }
826 0 : impars_p = impars;
827 0 : gridpars_p = gridpars;
828 0 : imageDefined_p = true;
829 :
830 0 : imparsVec_p.resize(imparsVec_p.nelements()+1, true);
831 0 : imparsVec_p[imparsVec_p.nelements()-1] = impars_p;
832 :
833 0 : gridparsVec_p.resize(gridparsVec_p.nelements()+1, true);
834 0 : gridparsVec_p[gridparsVec_p.nelements()-1] = gridpars_p;
835 :
836 0 : return true;
837 0 : }
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 0 : Bool SynthesisImagerVi2::weight(const Record& inrec){
888 0 : String type, rmode, filtertype;
889 0 : Quantity noise, fieldofview,filterbmaj,filterbmin, filterbpa;
890 : Double robust, fracBW;
891 : Int npixels;
892 : Bool multiField, useCubeBriggs;
893 0 : SynthesisUtilMethods::getFromWeightRecord(type, rmode,noise, robust,fieldofview,npixels, multiField, useCubeBriggs,
894 : filtertype, filterbmaj,filterbmin, filterbpa, fracBW, inrec);
895 0 : return weight(type, rmode,noise, robust,fieldofview,npixels, multiField, useCubeBriggs,
896 0 : filtertype, filterbmaj,filterbmin, filterbpa, fracBW );
897 :
898 :
899 0 : }
900 0 : 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 0 : LogIO os(LogOrigin("SynthesisImagerVi2", "weight()", WHERE));
908 :
909 0 : if(rmode=="bwtaper") //See CAS-13021 for bwtaper algorithm details
910 : {
911 0 : if(fracBW == 0.0)
912 : {
913 0 : Double minFreq = 0.0;
914 0 : Double maxFreq = 0.0;
915 :
916 0 : if(itsMaxShape(3) < 1) {
917 0 : cout << "SynthesisImagerVi2::weight Only one channel in image " << endl;
918 : }
919 : else{
920 0 : minFreq=abs(SpectralImageUtil::worldFreq(itsMaxCoordSys, 0.0));
921 0 : maxFreq=abs(SpectralImageUtil::worldFreq(itsMaxCoordSys,itsMaxShape(3)-1));
922 :
923 0 : if(maxFreq < minFreq){
924 0 : Double tmp=minFreq;
925 0 : minFreq=maxFreq;
926 0 : maxFreq=tmp;
927 : }
928 :
929 0 : if((maxFreq != 0.0) || (minFreq != 0.0)) fracBW = 2*(maxFreq - minFreq)/(maxFreq + minFreq);
930 :
931 0 : os << LogIO::NORMAL << " Fractional bandwidth used by briggsbwtaper " << fracBW << endl; //<< LogIO::POST;
932 :
933 : }
934 : }
935 : }
936 :
937 0 : weightParams_p=SynthesisUtilMethods::fillWeightRecord(type, rmode,noise, robust,fieldofview,
938 0 : 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 0 : Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
947 0 : Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
948 : os << LogIO::NORMAL // Loglevel INFO
949 0 : << "Set imaging weights : " ; //<< LogIO::POST;
950 :
951 0 : if (type=="natural") {
952 : os << LogIO::NORMAL // Loglevel INFO
953 0 : << "Natural weighting" << LogIO::POST;
954 0 : imwgt_p=VisImagingWeight("natural");
955 : }
956 0 : else if (type=="radial") {
957 0 : os << "Radial weighting" << LogIO::POST;
958 0 : imwgt_p=VisImagingWeight("radial");
959 : }
960 : else{
961 0 : vi_p->originChunks();
962 0 : vi_p->origin();
963 0 : if(!imageDefined_p)
964 0 : throw(AipsError("Need to define image"));
965 0 : Int nx=itsMaxShape[0];
966 0 : Int ny=itsMaxShape[1];
967 0 : Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
968 0 : Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
969 0 : if(type=="superuniform"){
970 0 : if(!imageDefined_p) throw(AipsError("Please define image first"));
971 0 : Int actualNpix=npixels;
972 0 : if(actualNpix <=0)
973 0 : actualNpix=3;
974 : os << LogIO::NORMAL // Loglevel INFO
975 : << "SuperUniform weighting over a square cell spanning ["
976 : << -actualNpix
977 0 : << ", " << actualNpix << "] in the uv plane" << LogIO::POST;
978 0 : imwgt_p=VisImagingWeight(*vi_p, rmode, noise, robust, nx,
979 : ny, cellx, celly, actualNpix,
980 0 : actualNpix, multiField);
981 : }
982 0 : else if ((type=="robust")||(type=="uniform")||(type=="briggs")) {
983 0 : if(!imageDefined_p) throw(AipsError("Please define image first"));
984 0 : Quantity actualFieldOfView_x(fieldofview), actualFieldOfView_y(fieldofview) ;
985 0 : Int actualNPixels_x(npixels),actualNPixels_y(npixels) ;
986 0 : String wtype;
987 0 : if(type=="briggs") {
988 0 : wtype = "Briggs";
989 : }
990 : else {
991 0 : wtype = "Uniform";
992 : }
993 0 : if(actualFieldOfView_x.get().getValue()==0.0&&actualNPixels_x==0) {
994 0 : actualNPixels_x=nx;
995 0 : actualFieldOfView_x=Quantity(actualNPixels_x*cellx.get("rad").getValue(),"rad");
996 0 : actualNPixels_y=ny;
997 0 : 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 0 : << 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 0 : << LogIO::POST;
1031 0 : Quantity actualCellSize_x(actualFieldOfView_x.get("rad").getValue()/actualNPixels_x, "rad");
1032 0 : 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 0 : if(useCubeBriggs){
1041 0 : String outstr=String("Doing spectral cube Briggs weighting formula -- " + rmode + (rmode=="abs" ? " with estimated noise "+ String::toString(noise.getValue())+noise.getUnit() : ""));
1042 0 : os << outstr << LogIO::POST;
1043 : //VisImagingWeight nat("natural");
1044 : //vi_p->useImagingWeight(nat);
1045 0 : 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 0 : CountedPtr<refim::BriggsCubeWeightor> bwgt=new refim::BriggsCubeWeightor(wtype=="Uniform" ? "none" : rmode, noise, robust,fracBW, npixels, multiField);
1049 0 : for (Int k=0; k < itsMappers.nMappers(); ++k){
1050 0 : itsMappers.getFTM2(k)->setBriggsCubeWeight(bwgt);
1051 : }
1052 :
1053 :
1054 0 : }
1055 : else
1056 : {
1057 0 : imwgt_p=VisImagingWeight(*vi_p, wtype=="Uniform" ? "none" : rmode, noise, robust,
1058 : actualNPixels_x, actualNPixels_y, actualCellSize_x,
1059 0 : 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 0 : }
1078 : else {
1079 : //this->unlock();
1080 : os << LogIO::SEVERE << "Unknown weighting " << type
1081 0 : << LogIO::EXCEPTION;
1082 0 : return false;
1083 : }
1084 0 : }
1085 :
1086 : //// UV-Tapering
1087 : //cout << "Taper type : " << filtertype << " : " << (filtertype=="gaussian") << endl;
1088 0 : if( filtertype == "gaussian" ) {
1089 : // os << "Setting uv-taper" << LogIO::POST;
1090 0 : imwgt_p.setFilter( filtertype, filterbmaj, filterbmin, filterbpa );
1091 : }
1092 0 : vi_p->useImagingWeight(imwgt_p);
1093 : ///////////////////////////////
1094 :
1095 0 : SynthesisUtilMethods::getResource("Set Weighting");
1096 :
1097 : /// return true;
1098 :
1099 0 : }
1100 0 : catch(AipsError &x)
1101 : {
1102 0 : throw( AipsError("Error in Weighting : "+x.getMesg()) );
1103 0 : }
1104 :
1105 0 : return true;
1106 0 : }
1107 :
1108 0 : 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 0 : LogIO log_l(LogOrigin("SynthesisImagerVi2", "appendToMapperList(ftm)"));
1123 : //---------------------------------------------
1124 : // Some checks..
1125 :
1126 0 : if(facets > 1 && itsMappers.nMappers() > 0)
1127 0 : log_l << "Facetted image has to be the first of multifields" << LogIO::EXCEPTION;
1128 :
1129 0 : TcleanProcessingInfo procInfo;
1130 0 : CompositeNumber cn(uInt(imshape[0] * 2));
1131 : // heuristic factors multiplied to imshape based on gridder
1132 0 : size_t fudge_factor = 15;
1133 0 : if (ftm->name()=="MosaicFTNew") {
1134 0 : fudge_factor = 20;
1135 : }
1136 0 : else if (ftm->name()=="GridFT") {
1137 0 : fudge_factor = 9;
1138 : }
1139 0 : std::tie(procInfo, std::ignore, std::ignore) =
1140 0 : nSubCubeFitInMemory(fudge_factor, imshape, padding);
1141 :
1142 : // chanchunks auto-calculation block, for now still here for awproject (CAS-12204)
1143 0 : 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 0 : gridpars_p.chanchunks=chanchunks;
1164 0 : 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 0 : 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 0 : if(chanchunks > 1) itsDataLoopPerMapper=true;
1173 :
1174 0 : AlwaysAssert( ( ( ! (ftm->name()=="MosaicFTNew" && mappertype=="imagemosaic") ) &&
1175 : ( ! (ftm->name()=="AWProjectWBFT" && mappertype=="imagemosaic") )) ,
1176 : AipsError );
1177 : //---------------------------------------------
1178 :
1179 : // Create the ImageStore object
1180 0 : CountedPtr<SIImageStore> imstor;
1181 0 : MSColumns msc(*(mss_p[0]));
1182 0 : imstor = createIMStore(imagename, csys, imshape, overwrite,msc, mappertype, ntaylorterms, distance, procInfo, facets, iftm->useWeightImage(), startmodel );
1183 :
1184 : // Create the Mappers
1185 0 : if( facets<2 && chanchunks<2) // One facet. Just add the above imagestore to the mapper list.
1186 : {
1187 0 : 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 0 : if ( facets>1 && chanchunks==1 )
1193 : {
1194 : // Make and connect the list.
1195 0 : Block<CountedPtr<SIImageStore> > imstorList = createFacetImageStoreList( imstor, facets );
1196 0 : for( uInt facet=0; facet<imstorList.nelements(); facet++)
1197 : {
1198 0 : CountedPtr<refim::FTMachine> new_ftm, new_iftm;
1199 0 : if(facet==0){ new_ftm = ftm; new_iftm = iftm; }
1200 0 : else{ new_ftm=ftm->cloneFTM(); new_iftm=iftm->cloneFTM(); }
1201 : // imstorList[facet]->setDataPolFrame(imstor->getDataPolFrame());
1202 0 : itsMappers.addMapper(createSIMapper( mappertype, imstorList[facet], new_ftm, new_iftm, ntaylorterms));
1203 0 : }
1204 0 : }// 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 0 : }
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 0 : std::tuple<TcleanProcessingInfo, Vector<Int>, Vector<Int> > SynthesisImagerVi2::nSubCubeFitInMemory(const Int fudge_factor, const IPosition& imshape, const Float padding){
1240 0 : LogIO log_l(LogOrigin("SynthesisImagerVi2", "nSubCubeFitInMemory"));
1241 :
1242 0 : Double required_mem = fudge_factor * sizeof(Float);
1243 0 : int nsubcube=1;
1244 0 : CompositeNumber cn(uInt(imshape[0] * 2));
1245 0 : for (size_t i = 0; i < imshape.nelements(); i++) {
1246 : // gridding pads image and increases to composite number
1247 0 : if (i < 2) {
1248 0 : required_mem *= cn.nextLargerEven(Int(padding*Float(imshape[i])-0.5));
1249 : }
1250 : else {
1251 0 : required_mem *= imshape[i];
1252 : }
1253 : }
1254 :
1255 : // get number of tclean processes running on the same machine
1256 0 : size_t nlocal_procs = 1;
1257 0 : 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 0 : required_mem *= nlocal_procs;
1264 : Double usr_memfrac, usr_mem;
1265 0 : AipsrcValue<Double>::find(usr_memfrac, "system.resources.memfrac", 80.);
1266 0 : AipsrcValue<Double>::find(usr_mem, "system.resources.memory", -1.);
1267 : Double memory_avail;
1268 0 : if (usr_mem > 0.) {
1269 0 : memory_avail = usr_mem * 1024. * 1024.;
1270 : }
1271 : else {
1272 0 : memory_avail = Double(HostInfo::memoryFree()) * (usr_memfrac / 100.) * 1024.;
1273 : }
1274 : // compute required chanchunks to fit into the available memory
1275 0 : nsubcube = (int)std::ceil((Double)required_mem / memory_avail);
1276 0 : Int nworkers= applicator.numProcs() < 2 ? 1 : applicator.numProcs()-1;
1277 0 : if((nsubcube/nworkers) >1 && nworkers !=1){
1278 0 : nsubcube=(Int(std::floor(Float(nsubcube)/Float(nworkers)))+1)*nworkers;
1279 :
1280 : }
1281 0 : if (imshape.nelements() == 4 && imshape[3] < nsubcube) {
1282 0 : nsubcube = imshape[3];
1283 : // TODO make chanchunks a divisor of nchannels?
1284 : }
1285 0 : nsubcube = nsubcube < 1 ? 1 : nsubcube;
1286 0 : 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 0 : else if(imshape[3] < (applicator.numProcs()-1)){
1294 0 : nsubcube=imshape[3];
1295 : }
1296 0 : Int chunksize=imshape[3]/nsubcube;
1297 0 : Int rem=imshape[3] % nsubcube;
1298 : //case of nchan < numprocs
1299 0 : 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 0 : Vector<Int> start(nsubcube,0);
1315 0 : Vector<Int> end(nsubcube,chunksize-1);
1316 0 : if(rem >0){
1317 0 : end(0)+=1;
1318 0 : --rem;
1319 : }
1320 0 : 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 0 : const float toGB = 1024.0 * 1024.0 * 1024.0;
1332 0 : std::ostringstream usr_mem_msg;
1333 0 : if (usr_mem > 0.) {
1334 0 : usr_mem_msg << usr_mem / 1024.;
1335 : } else {
1336 0 : usr_mem_msg << "-";
1337 : }
1338 0 : std::ostringstream oss;
1339 0 : oss << setprecision(4);
1340 0 : oss << "Required memory: " << required_mem / toGB
1341 0 : << " GB. Available mem.: " << memory_avail / toGB
1342 0 : << " GB (rc, mem. fraction: " << usr_memfrac
1343 0 : << "%, memory: " << usr_mem_msg.str()
1344 0 : << ") => Subcubes: " << nsubcube
1345 0 : << ". Processes on node: " << nlocal_procs << ".\n";
1346 0 : log_l << oss.str() << LogIO::POST;
1347 :
1348 0 : TcleanProcessingInfo procInfo;
1349 0 : procInfo.mpiprocs = nlocal_procs;
1350 0 : procInfo.chnchnks = nsubcube;
1351 0 : procInfo.memavail = memory_avail / toGB;
1352 0 : procInfo.memreq = required_mem / toGB;
1353 0 : return make_tuple(procInfo, start, end);
1354 0 : }
1355 :
1356 0 : void SynthesisImagerVi2::runMajorCycleCube( const Bool dopsf,
1357 : const Record inpcontrol) {
1358 0 : LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycleCube",WHERE) );
1359 0 : if(dopsf){
1360 0 : runCubeGridding(True);
1361 : ///Store the beamsets in ImageInfo
1362 0 : for(Int k=0; k < itsMappers.nMappers(); ++k){
1363 :
1364 0 : (itsMappers.imageStore(k))->getBeamSet();
1365 : }
1366 : }
1367 : else
1368 0 : runCubeGridding(False, inpcontrol);
1369 :
1370 :
1371 :
1372 0 : }
1373 0 : void SynthesisImagerVi2::runMajorCycle(const Bool dopsf,
1374 : const Bool savemodel)
1375 : {
1376 0 : LogIO os( LogOrigin("SynthesisImagerVi2","runMajorCycle",WHERE) );
1377 :
1378 : // cout << "Savemodel : " << savemodel << " readonly : " << readOnly_p << " usescratch : " << useScratch_p << endl;
1379 0 : Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
1380 0 : Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;
1381 :
1382 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
1383 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
1384 :
1385 0 : SynthesisUtilMethods::getResource("Start Major Cycle");
1386 :
1387 0 : itsMappers.checkOverlappingModels("blank");
1388 :
1389 : {
1390 0 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
1391 0 : vi_p->originChunks();
1392 0 : vi_p->origin();
1393 0 : Double numcoh=0;
1394 0 : for (uInt k=0; k< mss_p.nelements(); ++k)
1395 0 : numcoh+=Double(mss_p[k]->nrow());
1396 : ProgressMeter pm(1.0, numcoh,
1397 0 : dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
1398 0 : rownr_t cohDone=0;
1399 :
1400 :
1401 0 : if(!dopsf)itsMappers.initializeDegrid(*vb);
1402 0 : itsMappers.initializeGrid(*vi_p,dopsf);
1403 0 : 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 0 : vi_p->originChunks();
1407 0 : vi_p->origin();
1408 0 : Int nchannow=vb->nChannels();
1409 0 : Int spwnow=vb->spectralWindows()[0];
1410 0 : Int nchaninms=MSColumns(vb->ms()).spectralWindow().numChan()(spwnow);
1411 : //cerr << "chans " << nchaninms << " " << nchannow << endl;
1412 :
1413 0 : 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 0 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
1420 : {
1421 :
1422 0 : 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 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
1427 : {
1428 0 : if(!dopsf) {
1429 0 : { Cube<Complex> mod(vb->nCorrelations(), vb->nChannels(), vb->nRows(), Complex(0.0));
1430 0 : vb->setVisCubeModel(mod);
1431 0 : }
1432 0 : itsMappers.degrid(*vb, savevirtualmodel );
1433 0 : if(savemodelcolumn && writeAccess_p ){
1434 0 : const_cast<MeasurementSet& >((vi_p->ms())).lock(true);
1435 0 : vi_p->writeVisModel(vb->visCubeModel());
1436 0 : 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 0 : itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)datacol_p);
1447 :
1448 0 : cohDone += vb->nRows();
1449 0 : 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 0 : SynthesisUtilMethods::getResource("Before finalize for all mappers");
1459 0 : if(!dopsf) itsMappers.finalizeDegrid(*vb);
1460 0 : itsMappers.finalizeGrid(*vb, dopsf);
1461 :
1462 0 : }
1463 :
1464 0 : itsMappers.checkOverlappingModels("restore");
1465 :
1466 0 : unlockMSs();
1467 :
1468 0 : SynthesisUtilMethods::getResource("End Major Cycle");
1469 :
1470 0 : }// 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 0 : bool SynthesisImagerVi2::runCubeGridding(Bool dopsf, const Record inpcontrol){
1645 :
1646 0 : LogIO logger(LogOrigin("SynthesisImagerVi2", "runCubeGridding", WHERE));
1647 :
1648 : //dummy for now as init is overloaded on this signature
1649 0 : int argc=1;
1650 0 : char **argv=nullptr;
1651 0 : casa::applicator.init ( argc, argv );
1652 : //For serial or master only
1653 0 : if ( applicator.isController() ) {
1654 0 : CubeMajorCycleAlgorithm cmc;
1655 : //casa::applicator.defineAlgorithm(cmc);
1656 : //Initialize everything to get the setup as serial
1657 : {
1658 :
1659 0 : vi_p->originChunks();
1660 0 : vi_p->origin();
1661 :
1662 : }
1663 : ///Break things into chunks for parallelization or memory availabbility
1664 0 : Vector<Int> startchan;
1665 0 : Vector<Int> endchan;
1666 : Int numchunks;
1667 0 : Int fudge_factor = 15;
1668 0 : if ((itsMappers.getFTM2(0))->name()=="MosaicFTNew") {
1669 0 : fudge_factor = 20;
1670 : }
1671 0 : else if ((itsMappers.getFTM2(0))->name()=="GridFT") {
1672 0 : fudge_factor = 9;
1673 : }
1674 0 : TcleanProcessingInfo procInfo;
1675 0 : std::tie(procInfo, startchan, endchan)=nSubCubeFitInMemory(fudge_factor, itsMaxShape, gridpars_p.padding);
1676 0 : 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 0 : 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 0 : 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 0 : if(!dopsf){
1692 :
1693 : //controlRecord.define("lastcycle", savemodel);
1694 0 : controlRecord.define("nmajorcycles", nMajorCycles);
1695 0 : Vector<String> modelnames(Int(imparsVec_p.nelements()),"");
1696 0 : for(uInt k=0; k < imparsVec_p.nelements(); ++k){
1697 0 : Int imageStoreId=k;
1698 0 : if(k>0){
1699 0 : 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 0 : if(gridparsVec_p[0].facets >1)
1702 0 : imageStoreId+=gridparsVec_p[0].facets*gridparsVec_p[0].facets-1;
1703 : }
1704 0 : if((itsMappers.imageStore(imageStoreId))->hasModel()){
1705 0 : modelnames(k)=imparsVec_p[k].imageName+".model";
1706 0 : (itsMappers.imageStore(k))->model()->unlock();
1707 0 : controlRecord.define("modelnames", modelnames);
1708 : }
1709 : }
1710 :
1711 0 : }
1712 0 : Vector<String> weightnames(Int(imparsVec_p.nelements()),"");
1713 0 : Vector<String> pbnames(Int(imparsVec_p.nelements()),"");
1714 0 : for(uInt k=0; k < imparsVec_p.nelements(); ++k){
1715 0 : Int imageStoreId=k;
1716 0 : if(k>0){
1717 0 : 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 0 : 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 0 : if((itsMappers.getFTM2(imageStoreId))->useWeightImage()){
1724 : //cerr << "Mosaic weight image " << max(itsMappers.imageStore(imageStoreId)->weight(k)->get()) << endl;
1725 0 : weightnames(k)=itsMappers.imageStore(imageStoreId)->weight(k)->name();
1726 :
1727 :
1728 : }
1729 0 : if(itsMakeVP){
1730 0 : pbnames(k)=itsMappers.imageStore(imageStoreId)->pb(k)->name();
1731 0 : (itsMappers.imageStore(imageStoreId)->pb(k))->unlock();
1732 : }
1733 : }
1734 0 : controlRecord.define("weightnames", weightnames);
1735 0 : 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 0 : controlRecord.define("dividebyweight", True);
1739 0 : controlRecord.defineRecord("normpars", normpars_p);
1740 : ///Let's see if no per chan weight density was chosen
1741 0 : String weightdensityimage=getWeightDensity();
1742 0 : if(weightdensityimage != "")
1743 0 : controlRecord.define("weightdensity", weightdensityimage);
1744 :
1745 0 : Record vecSelParsRec, vecImParsRec, vecGridParsRec;
1746 0 : Vector<Int>vecRange(2);
1747 0 : for (uInt k = 0; k < dataSel_p.nelements(); ++k) {
1748 0 : Record selparsRec = dataSel_p[k].toRecord();
1749 0 : vecSelParsRec.defineRecord(String::toString(k), selparsRec);
1750 0 : }
1751 0 : for (uInt k=0; k < imparsVec_p.nelements(); ++k){
1752 0 : Record imparsRec = imparsVec_p[k].toRecord();
1753 : //need to send polrep
1754 0 : imparsRec.define("polrep", Int((itsMappers.imageStore(k))->getDataPolFrame()));
1755 : //need to send movingSourceName if any
1756 0 : imparsRec.define("movingsource", movingSource_p);
1757 0 : 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 0 : vecImParsRec.defineRecord(String::toString(k), imparsRec);
1767 0 : vecGridParsRec.defineRecord(String::toString(k), gridparsRec);
1768 0 : }
1769 0 : String workingdir="";
1770 : //Int numchan=(dopsf) ? imstor->psf()->shape()[3] : imstor->residual()->shape() [3];
1771 : //copy the imageinfo of main image here
1772 0 : if(dopsf)
1773 0 : cubePsfImageInfo_p=(itsMappers.imageStore(0)->psf())->imageInfo();
1774 0 : for(Int k=0; k < itsMappers.nMappers(); ++k){
1775 :
1776 0 : if(dopsf){
1777 :
1778 0 : for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(true)); ++j){
1779 : ///TESTOO
1780 : //(itsMappers.imageStore(k))->psf(j)->set(0.0);
1781 : /////////
1782 0 : (itsMappers.imageStore(k))->psf(j)->unlock();
1783 0 : (itsMappers.imageStore(k))->pb()->unlock();
1784 : }
1785 : }
1786 : else{
1787 0 : for(uInt j =0; j <(itsMappers.imageStore(k)->getNTaylorTerms(false)); ++j){
1788 : /////////TESTOO
1789 : //(itsMappers.imageStore(k))->residual(j)->set(0.0);
1790 : ///////
1791 0 : (itsMappers.imageStore(k))->residual(j)->unlock();
1792 : }
1793 : }
1794 0 : 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 0 : Path namewgt( (itsMappers.imageStore(k))->sumwt(j)->name());
1799 0 : workingdir=namewgt.dirName();
1800 : ///TESTOO
1801 : //(itsMappers.imageStore(k))->sumwt(j)->set(0.0);
1802 : ////
1803 0 : (itsMappers.imageStore(k))->sumwt(j)->unlock();
1804 : //(itsMappers.imageStore(k))->releaseLocks();
1805 0 : }
1806 0 : (itsMappers.imageStore(k))->releaseLocks();
1807 : }
1808 : //Send the working directory as the child and master may be at different places
1809 :
1810 0 : 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 0 : Int rank ( 0 );
1817 : Bool assigned; //(casa::casa::applicator.nextAvailProcess(pwrite, rank));
1818 0 : Bool allDone ( false );
1819 0 : Vector<Bool> retvals(numchunks, False);
1820 0 : Int indexofretval=0;
1821 0 : for ( Int k=0; k < numchunks; ++k ) {
1822 0 : assigned=casa::applicator.nextAvailProcess ( cmc, rank );
1823 : //cerr << "assigned "<< assigned << endl;
1824 0 : 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 0 : applicator.put ( vecSelParsRec );
1852 : // put image sel params #2
1853 0 : applicator.put ( vecImParsRec );
1854 : // put gridders params #3
1855 0 : applicator.put ( vecGridParsRec );
1856 : // put which channel to process #4
1857 0 : vecRange(0)=startchan(k);
1858 0 : vecRange(1)=endchan(k);
1859 0 : applicator.put ( vecRange );
1860 : // psf or residual CubeMajorCycleAlgorithm #5
1861 0 : applicator.put ( dopsf );
1862 : // store modelvis and other controls #6
1863 0 : applicator.put ( controlRecord );
1864 : // weighting scheme #7
1865 0 : applicator.put( weightParams_p);
1866 : /// Tell worker to process it
1867 0 : applicator.apply ( cmc );
1868 :
1869 : }
1870 : // Wait for all outstanding processes to return
1871 0 : rank = casa::applicator.nextProcessDone ( cmc, allDone );
1872 0 : while ( !allDone ) {
1873 : Bool status;
1874 0 : Record returnRec;
1875 0 : casa::applicator.get(returnRec);
1876 0 : casa::applicator.get ( status );
1877 0 : if(dopsf)
1878 0 : updateImageBeamSet(returnRec);
1879 0 : if(returnRec.isDefined("tempfilenames")){
1880 0 : std::vector<String> b=returnRec.asArrayString("tempfilenames").tovector();
1881 0 : tempFileNames_p.insert(std::end(tempFileNames_p), std::begin(b), std::end(b));
1882 0 : }
1883 0 : retvals(indexofretval)=status;
1884 0 : ++indexofretval;
1885 0 : if ( status )
1886 : //cerr << "remainder rank " << rank << " successful " << endl;
1887 0 : cerr << "";
1888 : else
1889 0 : logger << LogIO::SEVERE << "remainder rank " << rank << " failed " << LogIO::POST;
1890 :
1891 0 : rank = casa::applicator.nextProcessDone ( cmc, allDone );
1892 0 : if(casa::applicator.isSerial())
1893 0 : allDone=true;
1894 0 : }
1895 0 : if(anyEQ(retvals, False)){
1896 : //cerr << retvals << endl;
1897 0 : ostringstream oss;
1898 : oss << "One or more of the cube section failed in de/gridding. Return values for "
1899 0 : "the sections: " << retvals;
1900 0 : throw(AipsError(oss));
1901 0 : }
1902 0 : if(!dopsf && normpars_p.isDefined("pblimit") && (normpars_p.asFloat("pblimit") > 0.0) ){
1903 : try{
1904 0 : SIImageStore::copyMask(itsMappers.imageStore(0)->pb(), itsMappers.imageStore(0)->residual());
1905 0 : (itsMappers.imageStore(0))->residual()->unlock();
1906 : //(itsMappers.imageStore(0)->pb())->pixelMask().unlock();
1907 0 : (itsMappers.imageStore(0))->pb()->unlock();
1908 : }
1909 0 : catch(AipsError &x) {
1910 0 : if(!String(x.getMesg()).contains("T/F"))
1911 0 : throw(AipsError(x.getMesg()));
1912 : else{
1913 0 : 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 0 : }
1918 : }
1919 : else{
1920 0 : LatticeLocker lock1 (*(itsMappers.imageStore(0)->psf()), FileLocker::Write);
1921 0 : itsMappers.imageStore(0)->psf()->setImageInfo(cubePsfImageInfo_p);
1922 0 : itsMappers.imageStore(0)->psf()->unlock();
1923 0 : (itsMappers.imageStore(0))->pb()->unlock();
1924 0 : }
1925 :
1926 0 : }
1927 :
1928 :
1929 0 : return true;
1930 :
1931 0 : }
1932 : /////////////////////////
1933 0 : void SynthesisImagerVi2::predictModel(){
1934 0 : LogIO os( LogOrigin("SynthesisImagerVi2","predictModel ",WHERE) );
1935 :
1936 0 : os << "---------------------------------------------------- Predict Model ---------------------------------------------" << LogIO::POST;
1937 :
1938 0 : Bool savemodelcolumn = !readOnly_p && useScratch_p;
1939 0 : Bool savevirtualmodel = !readOnly_p && !useScratch_p;
1940 : //os<<"PREDICTMODEL: readOnly_p=="<<readOnly_p<<" useScratch_p=="<<useScratch_p<<LogIO::POST;
1941 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
1942 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
1943 :
1944 0 : itsMappers.checkOverlappingModels("blank");
1945 :
1946 :
1947 : {
1948 0 : vi::VisBuffer2* vb = vi_p->getVisBuffer();;
1949 0 : vi_p->originChunks();
1950 0 : vi_p->origin();
1951 0 : Double numberCoh=0;
1952 0 : for (uInt k=0; k< mss_p.nelements(); ++k)
1953 0 : numberCoh+=Double(mss_p[k]->nrow());
1954 :
1955 0 : ProgressMeter pm(1.0, numberCoh, "Predict Model", "","","",true);
1956 0 : rownr_t cohDone=0;
1957 :
1958 0 : itsMappers.initializeDegrid(*vb);
1959 0 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
1960 : {
1961 :
1962 0 : 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 0 : vb->setVisCubeModel(Complex(0.0, 0.0));
1967 0 : itsMappers.degrid(*vb, savevirtualmodel);
1968 :
1969 0 : if(savemodelcolumn && writeAccess_p ){
1970 :
1971 0 : const_cast<MeasurementSet& >((vi_p->ms())).lock(true);
1972 :
1973 0 : vi_p->writeVisModel(vb->visCubeModel());
1974 :
1975 : //cerr << "nRows "<< vb->nRows() << " " << max(vb->visCubeModel()) << endl;
1976 0 : const_cast<MeasurementSet& >((vi_p->ms())).unlock();
1977 : }
1978 :
1979 0 : cohDone += vb->nRows();
1980 0 : pm.update(Double(cohDone));
1981 :
1982 : }
1983 : }
1984 0 : itsMappers.finalizeDegrid(*vb);
1985 0 : }
1986 :
1987 0 : itsMappers.checkOverlappingModels("restore");
1988 0 : itsMappers.releaseImageLocks();
1989 0 : unlockMSs();
1990 :
1991 0 : }// end of predictModel
1992 :
1993 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1994 0 : void SynthesisImagerVi2::makeSdImage(Bool dopsf)
1995 : {
1996 0 : LogIO os( LogOrigin("SynthesisImagerVi2","makeSdImage",WHERE) );
1997 :
1998 : // Bool dopsf=false;
1999 0 : if(datacol_p==FTMachine::PSF) dopsf=true;
2000 :
2001 : {
2002 0 : vi::VisBuffer2* vb = vi_p->getVisBuffer();;
2003 0 : vi_p->originChunks();
2004 0 : vi_p->origin();
2005 :
2006 0 : Double numberCoh=0;
2007 0 : for (uInt k=0; k< mss_p.nelements(); ++k)
2008 0 : numberCoh+=Double(mss_p[k]->nrow());
2009 :
2010 0 : ProgressMeter pm(1.0, numberCoh, "Predict Model", "","","",true);
2011 0 : rownr_t cohDone=0;
2012 :
2013 0 : itsMappers.initializeGrid(*vi_p,dopsf);
2014 0 : for (vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
2015 : {
2016 0 : if (vi_p->getImpl()->isNewMs()) {
2017 0 : itsMappers.handleNewMs(vi_p->ms());
2018 : }
2019 0 : for (vi_p->origin(); vi_p->more(); vi_p->next())
2020 : {
2021 0 : itsMappers.grid(*vb, dopsf, (refim::FTMachine::Type)datacol_p);
2022 0 : cohDone += vb->nRows();
2023 0 : pm.update(Double(cohDone));
2024 : }
2025 : }
2026 0 : itsMappers.finalizeGrid(*vb, dopsf);
2027 :
2028 0 : }
2029 :
2030 0 : unlockMSs();
2031 :
2032 0 : }// end makeSDImage
2033 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2034 0 : void SynthesisImagerVi2::makeImage(String type, const String& imagename, const String& complexImage, const Int whichModel)
2035 : {
2036 0 : LogIO os( LogOrigin("SynthesisImager","makeImage",WHERE) );
2037 :
2038 :
2039 0 : refim::FTMachine::Type seType(refim::FTMachine::OBSERVED);
2040 0 : 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 0 : 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 0 : 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 0 : else if (type=="psf") {
2059 0 : seType=refim::FTMachine::PSF;
2060 : os << "Making point spread function "
2061 0 : << 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 0 : String imageName=(itsMappers.imageStore(whichModel))->getName();
2099 0 : String cImageName(complexImage);
2100 0 : if(complexImage=="") {
2101 0 : cImageName=imageName+".compleximage";
2102 : }
2103 0 : Bool keepComplexImage=(complexImage!="")||(type.contains("holography"));
2104 : //cerr << "name " << (itsMappers.getFTM2(whichModel))->name() << endl;
2105 0 : 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 0 : CoordinateSystem csys=itsMaxCoordSys;
2110 0 : IPosition shp=itsMaxShape;
2111 0 : if((itsMappers.imageStore(whichModel))->getShape().product()!=0 ){
2112 0 : csys= (itsMappers.imageStore(whichModel))->getCSys();
2113 0 : shp=(itsMappers.imageStore(whichModel))->getShape();
2114 : }
2115 0 : itsMappers.releaseImageLocks();
2116 0 : PagedImage<Float> theImage( shp, csys, imagename);
2117 0 : PagedImage<Complex> cImageImage(theImage.shape(),
2118 : theImage.coordinates(),
2119 0 : cImageName);
2120 0 : if(!keepComplexImage)
2121 0 : cImageImage.table().markForDelete();
2122 :
2123 :
2124 0 : Matrix<Float> weight;
2125 0 : if(cImageImage.shape()[3] > 1){
2126 0 : cImageImage.set(Complex(0.0));
2127 0 : cImageImage.tempClose();
2128 0 : makeComplexCubeImage(cImageName, seType, whichModel);
2129 : }
2130 : else{
2131 0 : (itsMappers.getFTM2(whichModel))->makeImage(seType, *vi_p, cImageImage, weight);
2132 : }
2133 0 : if(seType==refim::FTMachine::PSF){
2134 0 : StokesImageUtil::ToStokesPSF(theImage, cImageImage);
2135 0 : StokesImageUtil::normalizePSF(theImage);
2136 : }
2137 : else{
2138 0 : StokesImageUtil::To(theImage, cImageImage);
2139 : }
2140 0 : }
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 0 : unlockMSs();
2183 :
2184 0 : }// end makeImage
2185 : /////////////////////////////////////////////////////
2186 :
2187 0 : void SynthesisImagerVi2::makeComplexCubeImage(const String& cimage, const refim::FTMachine::Type imtype, const Int whichmodel){
2188 0 : CubeMakeImageAlgorithm *cmi=new CubeMakeImageAlgorithm();
2189 : // Dummies for using the right signature for init
2190 0 : Path cimpath(cimage);
2191 0 : String cimname=cimpath.absoluteName();
2192 : //cerr << "ABSOLUTE path " << cimname << endl;
2193 0 : Int argc = 1;
2194 0 : char **argv = nullptr;
2195 0 : casa::applicator.init(argc, argv);
2196 0 : if(applicator.isController())
2197 : {
2198 0 : Record vecSelParsRec;
2199 0 : for (uInt k = 0; k < dataSel_p.nelements(); ++k) {
2200 0 : Record selparsRec = dataSel_p[k].toRecord();
2201 0 : vecSelParsRec.defineRecord(String::toString(k), selparsRec);
2202 0 : }
2203 0 : Record imparsRec = imparsVec_p[whichmodel].toRecord();
2204 : //cerr << "which model " << whichmodel << " record " << imparsRec << endl;
2205 0 : Record gridparsRec = gridparsVec_p[whichmodel].toRecord();
2206 : ///Break things into chunks for parallelization or memory availabbility
2207 0 : Vector<Int> startchan;
2208 0 : Vector<Int> endchan;
2209 : Int numchunks;
2210 0 : Int fudge_factor = 15;
2211 0 : if ((itsMappers.getFTM2(0))->name()=="MosaicFTNew") {
2212 0 : fudge_factor = 20;
2213 : }
2214 0 : else if ((itsMappers.getFTM2(0))->name()=="GridFT") {
2215 0 : fudge_factor = 9;
2216 : }
2217 0 : TcleanProcessingInfo procInfo;
2218 0 : std::tie(procInfo, startchan, endchan)=nSubCubeFitInMemory(fudge_factor, itsMaxShape, gridpars_p.padding);
2219 0 : numchunks = procInfo.chnchnks;
2220 :
2221 0 : Int imageType=Int(imtype);
2222 0 : Int rank(0);
2223 : Bool assigned;
2224 0 : Bool allDone(false);
2225 0 : Vector<Int> chanRange(2);
2226 0 : for (Int k=0; k < numchunks; ++k) {
2227 0 : assigned=casa::applicator.nextAvailProcess ( *cmi, rank );
2228 : //cerr << "assigned "<< assigned << endl;
2229 0 : 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 0 : applicator.put(vecSelParsRec);
2243 : // put image sel params #2
2244 0 : applicator.put(imparsRec);
2245 : // put gridder params #3
2246 0 : applicator.put(gridparsRec);
2247 : // put which channel to process #4
2248 0 : chanRange(0)=startchan(k);
2249 0 : chanRange(1)=endchan(k);
2250 0 : applicator.put(chanRange);
2251 : //Type of image #5
2252 0 : applicator.put(imageType);
2253 : // weighting parameters #6
2254 0 : applicator.put( weightParams_p);
2255 : // complex imagename #7
2256 0 : applicator.put(cimname);
2257 0 : applicator.apply(*cmi);
2258 : }
2259 : // Wait for all outstanding processes to return
2260 0 : rank = casa::applicator.nextProcessDone(*cmi, allDone);
2261 0 : while (!allDone) {
2262 : Bool status;
2263 0 : 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 0 : rank = casa::applicator.nextProcessDone(*cmi, allDone);
2271 0 : if(casa::applicator.isSerial())
2272 0 : allDone=true;
2273 : }
2274 0 : }//applicator controller section
2275 0 : }
2276 :
2277 :
2278 :
2279 0 : 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 0 : LogIO os( LogOrigin("SynthesisImagerVi2","createSIMapper",WHERE) );
2286 :
2287 0 : CountedPtr<SIMapper> localMapper;
2288 :
2289 : try
2290 : {
2291 :
2292 0 : if( mappertype == "default" || mappertype == "multiterm" )
2293 : {
2294 0 : 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 0 : return localMapper;
2310 0 : }
2311 :
2312 0 : void SynthesisImagerVi2::lockMS(MeasurementSet& thisms){
2313 0 : thisms.lock(!readOnly_p);
2314 0 : thisms.antenna().lock(false);
2315 0 : thisms.dataDescription().lock(false);
2316 0 : thisms.feed().lock(false);
2317 0 : thisms.field().lock(false);
2318 0 : thisms.observation().lock(false);
2319 0 : thisms.polarization().lock(false);
2320 0 : thisms.processor().lock(false);
2321 0 : thisms.spectralWindow().lock(false);
2322 0 : thisms.state().lock(false);
2323 0 : if(!thisms.doppler().isNull())
2324 0 : thisms.doppler().lock(false);
2325 0 : if(!thisms.flagCmd().isNull())
2326 0 : thisms.flagCmd().lock(false);
2327 0 : if(!thisms.freqOffset().isNull())
2328 0 : thisms.freqOffset().lock(false);
2329 : //True here as we can write in that
2330 0 : if(!thisms.history().isNull())
2331 : // thisms.history().lock(!readOnly_p);
2332 0 : thisms.history().lock(false);
2333 0 : if(!thisms.pointing().isNull())
2334 0 : thisms.pointing().lock(false);
2335 : //we write virtual model here
2336 0 : if(!thisms.source().isNull())
2337 0 : thisms.source().lock(!readOnly_p);
2338 0 : if(!thisms.sysCal().isNull())
2339 0 : thisms.sysCal().lock(false);
2340 0 : if(!thisms.weather().isNull())
2341 0 : thisms.weather().lock(false);
2342 :
2343 0 : }
2344 : ///////////////////////
2345 :
2346 0 : Vector<SynthesisParamsSelect> SynthesisImagerVi2::tuneSelectData(){
2347 0 : vi_p->originChunks();
2348 0 : vi_p->origin();
2349 0 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
2350 0 : Int spwnow=vb->spectralWindows()[0];
2351 0 : 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 0 : if(nchaninms <30 && !(!readOnly_p && useScratch_p))
2360 0 : return dataSel_p;
2361 :
2362 0 : return SynthesisImager::tuneSelectData();
2363 : }
2364 : //////////////////////
2365 0 : void SynthesisImagerVi2::lockMSs(){
2366 0 : for(uInt i=0;i<mss_p.nelements();i++)
2367 : {
2368 :
2369 0 : MeasurementSet *ms_l = const_cast<MeasurementSet* >(mss_p[i]);
2370 0 : lockMS(*ms_l);
2371 : }
2372 :
2373 0 : }
2374 0 : void SynthesisImagerVi2::unlockMSs()
2375 : {
2376 0 : LogIO os( LogOrigin("SynthesisImagerVi2","unlockMSs",WHERE) );
2377 0 : for(uInt i=0;i<mss_p.nelements();i++)
2378 : {
2379 0 : os << LogIO::NORMAL2 << "Unlocking : " << (mss_p[i])->tableName() << LogIO::POST;
2380 0 : MeasurementSet *ms_l = const_cast<MeasurementSet* >(mss_p[i]);
2381 0 : ms_l->unlock();
2382 0 : ms_l->antenna().unlock();
2383 0 : ms_l->dataDescription().unlock();
2384 0 : ms_l->feed().unlock();
2385 0 : ms_l->field().unlock();
2386 0 : ms_l->observation().unlock();
2387 0 : ms_l->polarization().unlock();
2388 0 : ms_l->processor().unlock();
2389 0 : ms_l->spectralWindow().unlock();
2390 0 : ms_l->state().unlock();
2391 : //
2392 : // Unlock the optional sub-tables as well, if they are present
2393 : //
2394 0 : if(!(ms_l->source().isNull())) ms_l->source().unlock();
2395 0 : if(!(ms_l->doppler().isNull())) ms_l->doppler().unlock();
2396 0 : if(!(ms_l->flagCmd().isNull())) ms_l->flagCmd().unlock();
2397 0 : if(!(ms_l->freqOffset().isNull())) ms_l->freqOffset().unlock();
2398 0 : if(!(ms_l->history().isNull())) ms_l->history().unlock();
2399 0 : if(!(ms_l->pointing().isNull())) ms_l->pointing().unlock();
2400 0 : if(!(ms_l->sysCal().isNull())) ms_l->sysCal().unlock();
2401 0 : if(!(ms_l->weather().isNull())) ms_l->weather().unlock();
2402 : }
2403 0 : }
2404 0 : 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 0 : LogIO os( LogOrigin("SynthesisImagerVi2","createFTMachine",WHERE));
2451 :
2452 :
2453 0 : if (ftname == "gridft") {
2454 0 : if (facets >1) {
2455 : theFT = new refim::GridFT(
2456 : cache, tile,
2457 0 : gridFunction, mLocation_p, phaseCenter_p, padding,
2458 : useAutocorr, useDoublePrec
2459 0 : );
2460 : theIFT = new refim::GridFT(
2461 : cache, tile,
2462 0 : gridFunction, mLocation_p, phaseCenter_p, padding,
2463 : useAutocorr, useDoublePrec
2464 0 : );
2465 : }
2466 : else {
2467 : theFT = new refim::GridFT(
2468 : cache, tile,
2469 0 : gridFunction, mLocation_p, padding,
2470 : useAutocorr, useDoublePrec
2471 0 : );
2472 : theIFT = new refim::GridFT(
2473 : cache, tile,
2474 0 : gridFunction, mLocation_p, padding,
2475 : useAutocorr, useDoublePrec
2476 0 : );
2477 : }
2478 : }
2479 0 : else if (ftname == "wprojectft") {
2480 0 : Double maxW = -1.0;
2481 0 : Double minW = -1.0;
2482 0 : Double rmsW = -1.0;
2483 0 : if (wprojplane < 1)
2484 0 : casa::refim::WProjectFT::wStat(*vi_p, minW, maxW, rmsW);
2485 0 : if (facets > 1) {
2486 0 : theFT = new refim::WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
2487 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2488 0 : theIFT = new refim::WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
2489 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2490 : }
2491 : else {
2492 0 : theFT = new refim::WProjectFT(wprojplane, mLocation_p,
2493 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2494 0 : theIFT = new refim::WProjectFT(wprojplane, mLocation_p,
2495 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
2496 : }
2497 : CountedPtr<refim::WPConvFunc> sharedconvFunc =
2498 0 : static_cast<refim::WProjectFT &>(*theFT).getConvFunc();
2499 : //static_cast<WProjectFT &>(*theFT).setConvFunc(sharedconvFunc);
2500 0 : static_cast<refim::WProjectFT &>(*theIFT).setConvFunc(sharedconvFunc);
2501 0 : }
2502 :
2503 0 : else if ( ftname == "mosaic" || ftname== "mosft" || ftname == "mosaicft" || ftname== "MosaicFT" || ftname == "awp2"){
2504 :
2505 0 : createMosFTMachine(theFT, theIFT, padding, useAutocorr, useDoublePrec, rotatePAStep, stokes, conjBeams);
2506 : }
2507 0 : else if ((ftname.at(0,3)=="awp") || (ftname== "mawprojectft") || (ftname == "protoft")) {
2508 0 : 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 0 : else if ( ftname == "mosaic" or
2516 0 : ftname == "mosft" or
2517 0 : ftname == "mosaicft" or
2518 0 : ftname == "MosaicFT" ) {
2519 0 : createMosFTMachine(
2520 : theFT, theIFT,
2521 : padding, useAutocorr, useDoublePrec, rotatePAStep, stokes, conjBeams
2522 : );
2523 : }
2524 0 : else if ( ftname == "sd" ) {
2525 0 : 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 0 : if( mType == "multiterm" ) {
2539 0 : AlwaysAssert( nTaylorTerms>=1 , AipsError );
2540 :
2541 : CountedPtr<refim::FTMachine> theMTFT =
2542 0 : new refim::MultiTermFTNew( theFT, nTaylorTerms, true /*forward*/ );
2543 : CountedPtr<refim::FTMachine> theMTIFT =
2544 0 : new refim::MultiTermFTNew( theIFT, nTaylorTerms, false /*forward*/ );
2545 :
2546 0 : theFT = theMTFT;
2547 0 : theIFT = theMTIFT;
2548 0 : }
2549 :
2550 : // Now, set the SkyJones if needed, and if not internally generated.
2551 0 : 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 0 : theFT->setFrameValidity( freqFrameValid );
2574 0 : theIFT->setFrameValidity( freqFrameValid );
2575 :
2576 : // Set interpolation mode
2577 0 : theFT->setFreqInterpolation( interpolation );
2578 0 : theIFT->setFreqInterpolation( interpolation );
2579 :
2580 : // Set tracking of moving source if any
2581 0 : if (movingSource_p != "") {
2582 0 : theFT->setMovingSource(movingSource_p);
2583 0 : 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 0 : if (pseudoI == true) {
2593 0 : os << "Turning on Pseudo-I gridding" << LogIO::POST;
2594 0 : theFT->setPseudoIStokes(true);
2595 0 : theIFT->setPseudoIStokes(true);
2596 : }
2597 :
2598 0 : }
2599 :
2600 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2601 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2602 0 : 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 0 : LogIO os( LogOrigin("SynthesisImagerVi2","createAWPFTMachine",WHERE));
2630 :
2631 0 : if (wprojPlane<=1)
2632 : {
2633 : os << LogIO::NORMAL
2634 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
2635 0 : << LogIO::POST;
2636 0 : 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 0 : MSObservationColumns msoc((mss_p[0])->observation());
2669 0 : String telescopeName=msoc.telescopeName()(0);
2670 : CountedPtr<refim::ConvolutionFunction> awConvFunc = refim::AWProjectFT::makeCFObject(telescopeName,
2671 : aTermOn,
2672 : psTermOn, (wprojPlane > 1),
2673 0 : mTermOn, wbAWP, conjBeams);
2674 :
2675 : //
2676 : // Construct the appropriate re-sampler.
2677 : //
2678 0 : CountedPtr<refim::VisibilityResamplerBase> visResampler;
2679 : // if (ftmName=="protoft") visResampler = new ProtoVR();
2680 : //elsef
2681 0 : 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 0 : CountedPtr<refim::CFCache> cfCacheObj;
2696 :
2697 :
2698 : //
2699 : // Finally construct the FTMachine with the CFCache, ConvFunc and
2700 : // Re-sampler objects.
2701 : //
2702 0 : Float pbLimit_l=1e-3;
2703 :
2704 0 : 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 0 : useDoublePrec);
2710 :
2711 0 : cfCacheObj = new refim::CFCache();
2712 0 : 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 0 : if(impars_p.mode.contains("cube")){
2720 0 : cfCacheObj->setLazyFill(False);
2721 : }
2722 : else
2723 0 : cfCacheObj->setLazyFill(refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1);
2724 : // cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
2725 0 : cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
2726 0 : cfCacheObj->initCache2(CFC_VERBOSE);
2727 :
2728 0 : theFT->setCFCache(cfCacheObj);
2729 :
2730 :
2731 0 : Quantity rotateOTF(rotatePAStep,"deg");
2732 0 : static_cast<refim::AWProjectWBFT &>(*theFT).setObservatoryLocation(mLocation_p);
2733 0 : 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 0 : theIFT = new refim::AWProjectWBFT(static_cast<refim::AWProjectWBFT &>(*theFT));
2746 :
2747 0 : os << "Sending frequency selection information " << mssFreqSel_p << " to AWP FTM." << LogIO::POST;
2748 0 : theFT->setSpwFreqSelection( mssFreqSel_p );
2749 0 : theIFT->setSpwFreqSelection( mssFreqSel_p );
2750 :
2751 :
2752 0 : }
2753 :
2754 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2755 :
2756 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2757 :
2758 :
2759 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2760 :
2761 0 : 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 0 : LogIO os(LogOrigin("SynthesisImagerVi2", "createMosFTMachine",WHERE));
2764 :
2765 0 : MSColumns msc(vi_p->ms());
2766 0 : String telescop=msc.observation().telescopeName()(0);
2767 0 : Bool multiTel=False;
2768 0 : Int msid=0;
2769 0 : for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk()){
2770 0 : 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 0 : vi_p->originChunks();
2776 :
2777 :
2778 :
2779 : PBMath::CommonPB kpb;
2780 0 : Record rec;
2781 0 : getVPRecord( rec, kpb, telescop );
2782 :
2783 :
2784 0 : 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 0 : refim::VPSkyJones* vps= nullptr;
2795 : //cerr << "rec " << rec << " kpb " << kpb << endl;
2796 : //cerr << "createMOs ftname " << gridpars_p.ftmachine << endl;
2797 0 : 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 0 : if(rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2814 0 : 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 0 : PBMath myPB(rec);
2822 0 : String whichPBMath;
2823 0 : PBMathInterface::namePBClass(myPB.whichPBClass(), whichPBMath);
2824 0 : os << "Using the PB defined by " << whichPBMath << " for beam calculation for telescope " << telescop << LogIO::POST;
2825 0 : vps= new refim::VPSkyJones(telescop, myPB, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
2826 : //kpb=PBMath::DEFAULT;
2827 0 : }
2828 :
2829 0 : theFT = new refim::MosaicFTNew(vps, mLocation_p, stokes, 1000000000, 16, useAutoCorr,
2830 0 : useDoublePrec, doConjBeams, gridpars_p.usePointing);
2831 0 : PBMathInterface::PBClass pbtype=((kpb==PBMath::EVLA) || multiTel)? PBMathInterface::COMMONPB: PBMathInterface::AIRY;
2832 0 : if(rec.asString("name")=="IMAGE")
2833 0 : pbtype=PBMathInterface::IMAGE;
2834 : ///Use Heterogenous array mode for the following
2835 0 : if((kpb == PBMath::UNKNOWN) || (kpb==PBMath::OVRO) || (kpb==PBMath::ACA)
2836 0 : || (kpb==PBMath::ALMA) || (kpb==PBMath::EVLA) || multiTel){
2837 0 : CountedPtr<refim::SimplePBConvFunc> mospb=new refim::HetArrayConvFunc(pbtype, itsVpTable);
2838 0 : static_cast<refim::MosaicFTNew &>(*theFT).setConvFunc(mospb);
2839 0 : }
2840 : ///////////////////make sure both FTMachine share the same conv functions.
2841 0 : theIFT= new refim::MosaicFTNew(static_cast<refim::MosaicFTNew &>(*theFT));
2842 : }
2843 0 : }
2844 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2845 :
2846 0 : 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 0 : LogIO os(LogOrigin("SynthesisImagerVi2", "createSDFTMachine", WHERE));
2866 : os << LogIO::NORMAL // Loglevel INFO
2867 0 : << "Performing single dish gridding..." << LogIO::POST;
2868 : os << LogIO::NORMAL1 // gridFunction is too cryptic for most users.
2869 0 : << "with convolution function " << gridFunction << LogIO::POST;
2870 :
2871 : // Now make the Single Dish Gridding
2872 : os << LogIO::NORMAL // Loglevel INFO
2873 0 : << "Gridding will use specified common tangent point:" << LogIO::POST;
2874 0 : os << LogIO::NORMAL << SynthesisUtilMethods::asComprehensibleDirectionString(phaseCenter_p)
2875 0 : << LogIO::POST; // Loglevel INFO
2876 0 : String const gridfunclower = downcase(gridFunction);
2877 0 : if(gridfunclower=="pb") {
2878 0 : refim::SkyJones *vp = nullptr;
2879 0 : MSColumns msc(*mss_p[0]);
2880 : // todo: NONE is OK?
2881 0 : BeamSquint::SquintType squintType = BeamSquint::NONE;
2882 0 : Quantity skyPosThresholdQuant((Double)skyPosThreshold+180.0, "deg");
2883 0 : Quantity parAngleStepQuant((Double)rotatePAStep, "deg");
2884 0 : if (itsVpTable.empty()) {
2885 : os << LogIO::NORMAL // Loglevel INFO
2886 0 : << "Using defaults for primary beams used in gridding" << LogIO::POST;
2887 0 : vp=new refim::VPSkyJones(msc, true, parAngleStepQuant, squintType,
2888 0 : 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 0 : theFT = new refim::SDGrid(mLocation_p, *vp, cache/2, tile, gridfunclower,
2897 0 : convSupport, minWeight, clipMinMax);
2898 0 : theIFT = new refim::SDGrid(mLocation_p, *vp, cache/2, tile, gridfunclower,
2899 0 : convSupport, minWeight, clipMinMax);
2900 0 : } else if (gridfunclower=="gauss" || gridfunclower=="gjinc") {
2901 0 : DirectionCoordinate dirCoord = itsMaxCoordSys.directionCoordinate();
2902 0 : Vector<String> units = dirCoord.worldAxisUnits();
2903 0 : Vector<Double> increments = dirCoord.increment();
2904 0 : Quantity cellx(increments[0], units[0]);
2905 0 : Quantity celly(increments[1], units[1]);
2906 0 : if (cellx != celly &&
2907 0 : ((!truncateSize.getUnit().empty()||truncateSize.getUnit()=="pixel")
2908 0 : || (!gwidth.getUnit().empty()||gwidth.getUnit()=="pixel")
2909 0 : || (!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 0 : if (truncateSize.getUnit().empty() || truncateSize.getUnit()=="pixel")
2916 0 : truncateValue = truncateSize.getValue();
2917 : else
2918 0 : truncateValue = truncateSize.getValue("rad")/celly.getValue("rad");
2919 0 : if (gwidth.getUnit().empty() || gwidth.getUnit()=="pixel")
2920 0 : gwidthValue = gwidth.getValue();
2921 : else
2922 0 : gwidthValue = gwidth.getValue("rad")/celly.getValue("rad");
2923 0 : if (jwidth.getUnit().empty() || jwidth.getUnit()=="pixel")
2924 0 : jwidthValue = jwidth.getValue();
2925 : else
2926 0 : jwidthValue = jwidth.getValue("rad")/celly.getValue("rad");
2927 0 : theFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2928 0 : truncateValue, gwidthValue, jwidthValue, minWeight, clipMinMax);
2929 0 : theIFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2930 0 : truncateValue, gwidthValue, jwidthValue, minWeight, clipMinMax);
2931 0 : }
2932 : else {
2933 0 : theFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2934 0 : convSupport, minWeight, clipMinMax);
2935 0 : theIFT = new refim::SDGrid(mLocation_p, cache/2, tile, gridfunclower,
2936 0 : convSupport, minWeight, clipMinMax);
2937 : }
2938 0 : theFT->setPointingDirColumn(pointingDirCol);
2939 0 : static_cast<refim::SDGrid*>(theFT.get())->setConvertFirst(convertFirst);
2940 0 : static_cast<refim::SDGrid*>(theIFT.get())->setConvertFirst(convertFirst);
2941 :
2942 : // turn on Pseudo Stokes mode if necessary
2943 0 : if (pseudoI || stokes == "XX" || stokes == "YY" || stokes == "XXYY"
2944 0 : || stokes == "RR" || stokes == "LL" || stokes == "RRLL") {
2945 0 : theFT->setPseudoIStokes(True);
2946 0 : theIFT->setPseudoIStokes(True);
2947 : }
2948 0 : }
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 0 : Bool SynthesisImagerVi2::setWeightDensity( const String& weightimagename)
2957 : {
2958 0 : LogIO os(LogOrigin("SynthesisImagerVi2", "setWeightDensity()", WHERE));
2959 : try
2960 : {
2961 0 : if(weightimagename.size() !=0){
2962 0 : Table::isReadable(weightimagename, True);
2963 0 : PagedImage<Float> im(weightimagename);
2964 0 : imwgt_p=VisImagingWeight(im);
2965 0 : im.unlock();
2966 0 : }
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 0 : vi_p->useImagingWeight(imwgt_p);
2995 0 : itsMappers.releaseImageLocks();
2996 :
2997 : }
2998 0 : catch (AipsError &x)
2999 : {
3000 0 : throw(AipsError("In setWeightDensity : "+x.getMesg()));
3001 0 : }
3002 0 : return true;
3003 0 : }
3004 :
3005 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3006 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3007 :
3008 :
3009 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3010 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3011 0 : void SynthesisImagerVi2::createVisSet(const Bool /*writeAccess*/)
3012 : {
3013 0 : LogIO os( LogOrigin("SynthesisImagerVi2","createVisSet",WHERE) );
3014 : //cerr << "mss_p num" << mss_p.nelements() << " sel " << fselections_p->size() << endl;
3015 0 : lockMSs();
3016 0 : 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 0 : vi_p=new vi::VisibilityIterator2(mss_p, vi::SortColumns(), true); //writeAccess);
3020 :
3021 0 : if(fselections_p->size() !=0){
3022 0 : CountedPtr<vi::FrequencySelections> tmpfselections=new FrequencySelections();
3023 : //Temporary fix till we get rid of old vi and we can get rid of tuneSelect
3024 0 : 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 0 : tmpfselections=fselections_p;
3031 : }
3032 : ////end of fix for tuneSelectdata
3033 0 : vi_p->setFrequencySelection (*tmpfselections);
3034 :
3035 0 : }
3036 : //
3037 0 : vi_p->originChunks();
3038 0 : vi_p->origin();
3039 : ////make sure to use the latest imaging weight scheme
3040 0 : vi_p->useImagingWeight(imwgt_p);
3041 0 : unlockMSs();
3042 0 : }// 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 0 : void SynthesisImagerVi2::dryGridding(const Vector<String>& cfList)
3055 : {
3056 0 : LogIO os( LogOrigin("SynthesisImagerVi2","dryGridding",WHERE) );
3057 :
3058 0 : rownr_t cohDone=0;
3059 0 : 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 0 : String ftmName = ((*(itsMappers.getFTM2(whichFTM)))).name();
3064 :
3065 0 : if (!((itsMappers.getFTM2(whichFTM,true))->isUsingCFCache())) return;
3066 :
3067 0 : 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 0 : vi::VisBuffer2* vb=vi_p->getVisBuffer();
3076 0 : vi_p->originChunks();
3077 0 : vi_p->origin();
3078 0 : Double numberCoh=0;
3079 0 : for (uInt k=0; k< mss_p.nelements(); ++k)
3080 0 : numberCoh+=Double(mss_p[k]->nrow());
3081 :
3082 0 : ProgressMeter pm(1.0, numberCoh, "dryGridding", "","","",true);
3083 :
3084 0 : itsMappers.initializeGrid(*vi_p);
3085 :
3086 : // Set the gridder (iFTM) to run in dry-gridding mode
3087 0 : (itsMappers.getFTM2(whichFTM,true))->setDryRun(true);
3088 :
3089 0 : Bool aTermIsOff=False;
3090 : {
3091 0 : CountedPtr<refim::FTMachine> ftm=itsMappers.getFTM2(whichFTM,True);
3092 0 : CountedPtr<refim::ConvolutionFunction> cf=ftm->getAWConvFunc();
3093 0 : aTermIsOff = cf->getTerm("ATerm")->isNoOp();
3094 0 : }
3095 :
3096 : os << "Making a \"blank\" CFCache"
3097 : << (aTermIsOff?" (without the A-Term)":"")
3098 0 : << 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 0 : for (vi_p->originChunks(); vi_p->moreChunks();vi_p->nextChunk())
3105 : {
3106 0 : for (vi_p->origin(); vi_p->more(); vi_p->next())
3107 : {
3108 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
3109 : {
3110 0 : itsMappers.grid(*vb, true, refim::FTMachine::OBSERVED, whichFTM);
3111 0 : cohDone += vb->nRows();
3112 0 : pm.update(Double(cohDone));
3113 : // If there is no term that depends on time, don't iterate over the entire data base
3114 0 : if (aTermIsOff) break;
3115 : }
3116 : }
3117 : }
3118 0 : }
3119 0 : if (cohDone == 0) os << "No valid rows found in dryGridding." << LogIO::EXCEPTION << LogIO::POST;
3120 : // Unset the dry-gridding mode.
3121 0 : (itsMappers.getFTM2(whichFTM,true))->setDryRun(false);
3122 :
3123 : //itsMappers.checkOverlappingModels("restore");
3124 0 : unlockMSs();
3125 : //fillCFCache(cfList);
3126 0 : }
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 0 : 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 0 : 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 0 : if (!ftmName.contains("awproject") and
3146 0 : !ftmName.contains("multitermftnew")) return;
3147 : //if (!ftmName.contains("awproject")) return;
3148 :
3149 0 : 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 0 : Float dPA=360.0,selectedPA=2*360.0;
3160 0 : if (cfList.nelements() > 0)
3161 : {
3162 0 : 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 0 : Vector<String> cfList_p=cfList;//dir.find(Regex(Regex::fromPattern("CFS*")));
3167 0 : Vector<String> wtCFList_p;
3168 0 : wtCFList_p.resize(cfList_p.nelements());
3169 0 : for (Int i=0; i<(Int)wtCFList_p.nelements(); i++) wtCFList_p[i]="WT"+cfList_p[i];
3170 :
3171 : //cerr << cfList_p << endl;
3172 0 : cfCacheObj->setCacheDir(cfcPath.data());
3173 :
3174 0 : os << "Re-loading the \"blank\" CFCache for filling" << LogIO::WARN << LogIO::POST;
3175 :
3176 0 : cfCacheObj->initCacheFromList2(cfcPath, cfList_p, wtCFList_p,
3177 : selectedPA, dPA,CFC_VERBOSE);
3178 :
3179 : // tmpFT->setCFCache(cfCacheObj);
3180 0 : Vector<Double> uvScale, uvOffset;
3181 0 : Matrix<Double> vbFreqSelection;
3182 0 : CountedPtr<refim::CFStore2> cfs2 = CountedPtr<refim::CFStore2>(&cfCacheObj->memCache2_p[0],false);//new CFStore2;
3183 0 : 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 0 : refim::AWConvFunc::makeConvFunction2(String(cfcPath), uvScale, uvOffset, vbFreqSelection,
3200 0 : *cfs2, *cfwts2, psTermOn, aTermOn, conjBeams);
3201 0 : }
3202 : //cerr << "Mem used = " << itsMappers.getFTM(whichFTM)->getCFCache()->memCache2_p[0].memUsage() << endl;
3203 : //(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getCFCache()->initCache2();
3204 0 : }
3205 :
3206 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3207 0 : void SynthesisImagerVi2::reloadCFCache()
3208 : {
3209 0 : LogIO os( LogOrigin("SynthesisImagerVi2","reloadCFCache",WHERE) );
3210 0 : Int whichFTM=0;
3211 0 : 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 0 : if (!(ftm->isUsingCFCache())) return; // Better check than checking against FTM name
3219 :
3220 0 : 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 0 : CountedPtr<refim::CFCache> cfc;
3237 0 : if (ftm->name().contains("MultiTerm")) cfc = ftm->getFTM2(true)->getCFCache();
3238 0 : else cfc = ftm->getCFCache();
3239 0 : cfc->setLazyFill((refim::SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1));
3240 0 : 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 0 : }
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 0 : Bool SynthesisImagerVi2::loadMosaicSensitivity(){
3317 0 : String ftmname=itsMappers.getFTM2(0)->name();
3318 :
3319 0 : if(ftmname.contains("Mosaic") || ftmname.contains("AWProjectWB")){
3320 : //sumwt has been calcuated
3321 0 : Bool donesumwt=(max(itsMappers.imageStore(0)->sumwt()->get()) > 0.0);
3322 : //cerr << "Done sumwght " << donesumwt << max(itsMappers.imageStore(0)->sumwt()->get()) << endl;
3323 0 : if(donesumwt){
3324 0 : IPosition shp=itsMappers.imageStore(0)->weight()->shape();
3325 0 : CoordinateSystem cs=itsMappers.imageStore(0)->weight()->coordinates();
3326 0 : CountedPtr<TempImage<Float> > wgtim=new TempImage<Float>(shp, cs);
3327 0 : wgtim->copyData(*(itsMappers.imageStore(0)->weight()));
3328 0 : (static_cast<refim::FTMachine &>( *(itsMappers.getFTM2(0,False)))).setWeightImage(*wgtim);
3329 0 : static_cast<refim::FTMachine &>( *(itsMappers.getFTM2(0,True))).setWeightImage(*wgtim);
3330 0 : return true;
3331 0 : }
3332 :
3333 :
3334 : }
3335 0 : return false;
3336 0 : }
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 0 : Bool SynthesisImagerVi2::makePB()
3446 : {
3447 0 : LogIO os( LogOrigin("SynthesisImagerVi2","makePB",WHERE) );
3448 :
3449 0 : if( itsMakeVP==False )
3450 : {
3451 0 : if( ((itsMappers.getFTM2(0))->name())!="MultiTermFTNew")
3452 0 : 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 0 : Bool doDefaultVP = itsVpTable.length()>0 ? False : True;
3460 :
3461 0 : CoordinateSystem coordsys=itsMappers.imageStore(0)->getCSys();
3462 0 : String telescope=coordsys.obsInfo().telescope();
3463 :
3464 0 : if (doDefaultVP) {
3465 :
3466 0 : MSAntennaColumns ac(mss_p[0]->antenna());
3467 0 : Double dishDiam=ac.dishDiameter()(0);
3468 0 : 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 0 : return makePBImage( telescope, False, dishDiam);
3473 0 : }
3474 : else{
3475 0 : return makePBImage(telescope );
3476 : }
3477 :
3478 0 : }
3479 :
3480 0 : return False;
3481 0 : }
3482 :
3483 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
3484 :
3485 0 : Bool SynthesisImagerVi2::makePrimaryBeam(PBMath& pbMath)
3486 : {
3487 0 : LogIO os( LogOrigin("SynthesisImagerVi2","makePrimaryBeam",WHERE) );
3488 :
3489 0 : os << "vi2 : Evaluating Primary Beam model onto image grid(s)" << LogIO::POST;
3490 :
3491 0 : itsMappers.initPB();
3492 :
3493 0 : vi::VisBuffer2* vb = vi_p->getVisBuffer();
3494 0 : vi_p->originChunks();
3495 0 : vi_p->origin();
3496 0 : std::map<Int, std::set<Int>> fieldsDone;
3497 0 : VisBufferUtil vbU(*vb);
3498 : ///////if tracking a moving source
3499 0 : MDirection origMovingDir;
3500 0 : MDirection newPhaseCenter;
3501 0 : Bool trackBeam=getMovingDirection(*vb, origMovingDir, True);
3502 : //////
3503 0 : for(vi_p->originChunks(); vi_p->moreChunks(); vi_p->nextChunk())
3504 : {
3505 0 : for (vi_p->origin(); vi_p->more(); vi_p->next())
3506 : {
3507 0 : Bool fieldDone=False;
3508 0 : if(fieldsDone.count(vb->msId() >0)){
3509 0 : fieldDone=fieldDone || (fieldsDone[vb->msId()].count(vb->fieldId()(0)) > 0);
3510 : }
3511 : else{
3512 0 : fieldsDone[vb->msId()]=std::set<int>();
3513 : }
3514 0 : if(!fieldDone){
3515 0 : fieldsDone[vb->msId()].insert(vb->fieldId()(0));
3516 0 : if(trackBeam){
3517 0 : MDirection newMovingDir;
3518 0 : getMovingDirection(*vb, newMovingDir);
3519 : //newPhaseCenter=vb->phaseCenter();
3520 0 : newPhaseCenter=vbU.getPhaseCenter(*vb);
3521 0 : newPhaseCenter.shift(MVDirection(-newMovingDir.getAngle()+origMovingDir.getAngle()), False);
3522 0 : }
3523 0 : itsMappers.addPB(*vb,pbMath, newPhaseCenter, trackBeam);
3524 :
3525 : }
3526 : }
3527 : }
3528 0 : itsMappers.releaseImageLocks();
3529 0 : unlockMSs();
3530 :
3531 0 : return True;
3532 0 : }// end makePB
3533 :
3534 0 : Bool SynthesisImagerVi2::getMovingDirection(const vi::VisBuffer2& vb, MDirection& outDir, const Bool useImageEpoch){
3535 0 : MDirection movingDir;
3536 0 : Bool trackBeam=False;
3537 :
3538 0 : MeasFrame mFrame(MEpoch(Quantity(vb.time()(0), "s"), MSColumns(vb.ms()).timeMeas()(0).getRef()), mLocation_p);
3539 0 : if(useImageEpoch){
3540 0 : mFrame.resetEpoch((itsMappers.imageStore(0))->getCSys().obsInfo().obsDate());
3541 :
3542 : }
3543 0 : if(movingSource_p != ""){
3544 : MDirection::Types refType;
3545 0 : trackBeam=True;
3546 0 : if(Table::isReadable(movingSource_p, False)){
3547 : //seems to be a table so assuming ephemerides table
3548 0 : Table laTable(movingSource_p);
3549 0 : Path leSentier(movingSource_p);
3550 0 : MeasComet laComet(laTable, leSentier.absoluteName());
3551 0 : movingDir.setRefString("COMET");
3552 0 : mFrame.set(laComet);
3553 0 : }
3554 : ///if not a table
3555 0 : 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 0 : else if(upcase(movingSource_p)=="TRACKFIELD"){
3562 0 : VisBufferUtil vbU(vb);
3563 0 : movingDir=vbU.getEphemDir(vb, MEpoch(mFrame.epoch()).get("s").getValue());
3564 0 : }
3565 : else{
3566 0 : throw(AipsError("Erroneous tracking direction set to make pb"));
3567 : }
3568 0 : MDirection::Ref outref1(MDirection::AZEL, mFrame);
3569 0 : MDirection tmphazel=MDirection::Convert(movingDir, outref1)();
3570 0 : MDirection::Ref outref(vb.phaseCenter().getRef().getType(), mFrame);
3571 0 : outDir=MDirection::Convert(tmphazel, outref)();
3572 0 : }
3573 : else{
3574 0 : outDir=vb.phaseCenter();
3575 0 : trackBeam=False;
3576 : }
3577 0 : return trackBeam;
3578 :
3579 :
3580 0 : }
3581 0 : CountedPtr<vi::VisibilityIterator2> SynthesisImagerVi2::getVi(){
3582 0 : return vi_p;
3583 : }
3584 0 : CountedPtr<refim::FTMachine> SynthesisImagerVi2::getFTM(const Int fid, Bool ift){
3585 0 : if(itsMappers.nMappers() <=fid)
3586 0 : throw(AipsError("Mappers have not been set for field "+String::toString(fid)));
3587 0 : return (itsMappers.getFTM2(fid, ift));
3588 :
3589 : }
3590 0 : void SynthesisImagerVi2::updateImageBeamSet(Record& returnRec){
3591 0 : if(returnRec.isDefined("imageid") && returnRec.asInt("imageid")==0){
3592 : //ImageInfo ii=(itsMappers.imageStore(0)->psf())->imageInfo();
3593 0 : Vector<Int> chanRange(0);
3594 0 : if(returnRec.isDefined("chanrange"))
3595 0 : chanRange=returnRec.asArrayInt("chanrange");
3596 0 : Int npol=(itsMappers.imageStore(0)->getShape())(2);
3597 0 : Int nchan=(itsMappers.imageStore(0)->getShape())(3);
3598 0 : if(chanRange.nelements() ==2 && chanRange(1) < nchan){
3599 :
3600 0 : ImageBeamSet iibeamset=cubePsfImageInfo_p.getBeamSet();
3601 0 : Matrix<GaussianBeam> mbeams=iibeamset.getBeams();
3602 0 : if(mbeams.nelements()==0){
3603 0 : mbeams.resize(itsMappers.imageStore(0)->getShape()(3), itsMappers.imageStore(0)->getShape()(2));
3604 0 : mbeams.set(GaussianBeam(Vector<Quantity>(3, Quantity(1e-12, "arcsec"))));
3605 : }
3606 0 : Cube<Float> recbeams(returnRec.asArrayFloat("beams"));
3607 0 : for(Int c=chanRange[0]; c <= chanRange[1]; ++c){
3608 0 : for(Int p=0; p < npol; ++p){
3609 0 : 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 0 : cubePsfImageInfo_p.setBeams(ImageBeamSet(mbeams));
3614 :
3615 0 : }
3616 : //itsMappers.imageStore(0)->psf()->setImageInfo(ii);
3617 : //itsMappers.imageStore(0)->psf()->unlock();
3618 :
3619 :
3620 0 : }
3621 0 : }
3622 :
3623 : } //# NAMESPACE CASA - END
3624 :
|
