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