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