Line data Source code
1 : //# SynthesisImager.cc: Implementation of SynthesisImager.h
2 : //# Copyright (C) 1997-2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <iostream>
30 : #include <sstream>
31 :
32 : #include <casacore/casa/Arrays/Matrix.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 :
36 :
37 : #include <casacore/casa/Logging.h>
38 : #include <casacore/casa/Logging/LogIO.h>
39 : #include <casacore/casa/Logging/LogMessage.h>
40 : #include <casacore/casa/Logging/LogSink.h>
41 : #include <casacore/casa/Logging/LogMessage.h>
42 : #include <casacore/casa/System/ProgressMeter.h>
43 :
44 : #include <casacore/casa/OS/DirectoryIterator.h>
45 : #include <casacore/casa/OS/File.h>
46 : #include <casacore/casa/OS/HostInfo.h>
47 : #include <casacore/casa/OS/Path.h>
48 : //#include <casa/OS/Memory.h>
49 :
50 : #include <casacore/lattices/LRegions/LCBox.h>
51 :
52 : #include <casacore/measures/Measures/MeasTable.h>
53 :
54 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
55 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
56 : #include <casacore/ms/MSSel/MSSelection.h>
57 :
58 : #include <casacore/tables/Tables/TableUtil.h>
59 :
60 : #include <synthesis/ImagerObjects/SynthesisImager.h>
61 :
62 : #include <synthesis/ImagerObjects/SynthesisUtilMethods.h>
63 : #include <synthesis/ImagerObjects/SIImageStore.h>
64 : #include <synthesis/ImagerObjects/SIImageStoreMultiTerm.h>
65 :
66 : #include <synthesis/MeasurementEquations/ImagerMultiMS.h>
67 : #include <synthesis/MeasurementEquations/VPManager.h>
68 : #include <msvis/MSVis/MSUtil.h>
69 : #include <msvis/MSVis/VisSetUtil.h>
70 : #include <msvis/MSVis/VisImagingWeight.h>
71 :
72 : #include <synthesis/TransformMachines/GridFT.h>
73 : #include <synthesis/TransformMachines/WPConvFunc.h>
74 : #include <synthesis/TransformMachines/WProjectFT.h>
75 : #include <synthesis/TransformMachines/VisModelData.h>
76 : #include <synthesis/TransformMachines/AWProjectFT.h>
77 : #include <synthesis/TransformMachines/HetArrayConvFunc.h>
78 : #include <synthesis/TransformMachines/MosaicFTNew.h>
79 : #include <synthesis/TransformMachines/MultiTermFTNew.h>
80 : #include <synthesis/TransformMachines/AWProjectWBFTNew.h>
81 : #include <synthesis/TransformMachines/AWConvFunc.h>
82 : #include <synthesis/TransformMachines/AWConvFuncEPJones.h>
83 : #include <synthesis/TransformMachines/NoOpATerm.h>
84 :
85 : #include <casacore/casa/Utilities/Regex.h>
86 : #include <casacore/casa/OS/Directory.h>
87 :
88 : //#include <casadbus/utilities/BusAccess.h>
89 : //#include <casadbus/session/DBusSession.h>
90 :
91 : #include <sys/types.h>
92 : #include <unistd.h>
93 :
94 :
95 : using namespace std;
96 :
97 : using namespace casacore;
98 : namespace casa { //# NAMESPACE CASA - BEGIN
99 :
100 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
101 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
102 :
103 34 : SynthesisImager::SynthesisImager() : itsMappers(SIMapperCollection()), writeAccess_p(True),
104 17 : gridpars_p(), impars_p(), movingSource_p(""), doingCubeGridding_p(True)
105 : {
106 :
107 17 : imwgt_p=VisImagingWeight("natural");
108 17 : imageDefined_p=false;
109 17 : useScratch_p=false;
110 17 : readOnly_p=true;
111 :
112 17 : mss4vi_p.resize(0);
113 17 : wvi_p=0;
114 17 : rvi_p=0;
115 :
116 17 : facetsStore_p=-1;
117 17 : chanChunksStore_p=-1;
118 17 : unFacettedImStore_p=NULL;
119 17 : unChanChunkedImStore_p=NULL;
120 17 : dataSel_p.resize();
121 :
122 17 : nMajorCycles=0;
123 :
124 17 : itsDataLoopPerMapper=false;
125 :
126 17 : }
127 :
128 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
129 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
130 :
131 17 : SynthesisImager::~SynthesisImager()
132 : {
133 34 : LogIO os( LogOrigin("SynthesisImager","destructor",WHERE) );
134 17 : os << LogIO::DEBUG1 << "SynthesisImager destroyed" << LogIO::POST;
135 17 : cleanupTempFiles();
136 17 : if(rvi_p) delete rvi_p;
137 17 : rvi_p=NULL;
138 : // cerr << "IN DESTR"<< endl;
139 : // VisModelData::listModel(mss4vi_p[0]);
140 :
141 17 : SynthesisUtilMethods::getResource("End SynthesisImager");
142 17 : }
143 :
144 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
145 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
146 :
147 0 : Bool SynthesisImager::selectData(const String& msname,
148 : const String& spw,
149 : const String& freqBeg,
150 : const String& freqEnd,
151 : const MFrequency::Types freqframe,
152 : const String& field,
153 : const String& antenna,
154 : const String& timestr,
155 : const String& scan,
156 : const String& obs,
157 : const String& state,
158 : const String& uvdist,
159 : const String& taql,
160 : const Bool usescratch,
161 : const Bool readonly,
162 : const Bool incrModel)
163 : {
164 0 : SynthesisParamsSelect pars;
165 0 : pars.msname=msname;
166 0 : pars.spw=spw;
167 0 : pars.freqbeg=freqBeg;
168 0 : pars.freqend=freqEnd;
169 0 : pars.freqframe=freqframe;
170 0 : pars.field=field;
171 0 : pars.antenna=antenna;
172 0 : pars.timestr=timestr;
173 0 : pars.scan=scan;
174 0 : pars.obs=obs;
175 0 : pars.state=state;
176 0 : pars.uvdist=uvdist;
177 0 : pars.taql=taql;
178 0 : pars.usescratch=usescratch;
179 0 : pars.readonly=readonly;
180 0 : pars.incrmodel=incrModel;
181 :
182 :
183 0 : String err = pars.verify();
184 :
185 0 : if( err.length()>0 ) throw(AipsError("Invalid Selection parameters : " + err));
186 :
187 0 : selectData( pars );
188 :
189 0 : return true;
190 0 : }
191 :
192 0 : Bool SynthesisImager::selectData(const SynthesisParamsSelect& selpars)
193 : {
194 0 : LogIO os( LogOrigin("SynthesisImager","selectData",WHERE) );
195 :
196 0 : SynthesisUtilMethods::getResource("Start SelectData");
197 :
198 : try
199 : {
200 :
201 : //Respect the readonly flag...necessary for multi-process access
202 0 : MeasurementSet thisms(selpars.msname, TableLock(TableLock::AutoNoReadLocking),
203 0 : selpars.readonly ? Table::Old : Table::Update);
204 0 : thisms.setMemoryResidentSubtables (MrsEligibility::defaultEligible());
205 :
206 0 : useScratch_p=selpars.usescratch;
207 0 : readOnly_p = selpars.readonly;
208 : // cout << "**************** usescr : " << useScratch_p << " readonly : " << readOnly_p << endl;
209 : //if you want to use scratch col...make sure they are there
210 0 : if(selpars.usescratch && !selpars.readonly){
211 0 : VisSetUtil::addScrCols(thisms, true, false, true, false);
212 0 : VisModelData::clearModel(thisms);
213 : }
214 0 : if(!selpars.incrmodel && !selpars.usescratch && !selpars.readonly)
215 0 : VisModelData::clearModel(thisms, selpars.field, selpars.spw);
216 :
217 0 : os << "MS : " << selpars.msname << " | ";
218 :
219 : //Some MSSelection
220 : //If everything is empty (which is valid) it will throw an exception..below
221 : //So make sure the main defaults are not empy i.e field and spw
222 0 : MSSelection thisSelection;
223 0 : if(selpars.field != ""){
224 0 : thisSelection.setFieldExpr(selpars.field);
225 0 : os << "Selecting on fields : " << selpars.field << " | " ;//LogIO::POST;
226 : }else
227 0 : thisSelection.setFieldExpr("*");
228 0 : if(selpars.spw != ""){
229 0 : thisSelection.setSpwExpr(selpars.spw);
230 0 : os << "Selecting on spw :"<< selpars.spw << " | " ;//LogIO::POST;
231 : }else
232 0 : thisSelection.setSpwExpr("*");
233 :
234 0 : if(selpars.antenna != ""){
235 0 : Vector<String> antNames(1, selpars.antenna);
236 : // thisSelection.setAntennaExpr(MSSelection::nameExprStr( antNames));
237 0 : thisSelection.setAntennaExpr(selpars.antenna);
238 0 : os << "Selecting on antenna names : " << selpars.antenna << " | " ;//LogIO::POST;
239 :
240 0 : }
241 0 : if(selpars.timestr != ""){
242 0 : thisSelection.setTimeExpr(selpars.timestr);
243 0 : os << "Selecting on time range : " << selpars.timestr << " | " ;//LogIO::POST;
244 : }
245 0 : if(selpars.uvdist != ""){
246 0 : thisSelection.setUvDistExpr(selpars.uvdist);
247 0 : os << "Selecting on uvdist : " << selpars.uvdist << " | " ;//LogIO::POST;
248 : }
249 0 : if(selpars.scan != ""){
250 0 : thisSelection.setScanExpr(selpars.scan);
251 0 : os << "Selecting on scan : " << selpars.scan << " | " ;//LogIO::POST;
252 : }
253 0 : if(selpars.obs != ""){
254 0 : thisSelection.setObservationExpr(selpars.obs);
255 0 : os << "Selecting on Observation Expr : " << selpars.obs << " | " ;//LogIO::POST;
256 : }
257 0 : if(selpars.state != ""){
258 0 : thisSelection.setStateExpr(selpars.state);
259 0 : os << "Selecting on Scan Intent/State : " << selpars.state << " | " ;//LogIO::POST;
260 : }
261 : // if(selpars.taql != ""){
262 : // thisSelection.setTaQLExpr(selpars.taql);
263 : // os << "Selecting via TaQL : " << selpars.taql << " | " ;//LogIO::POST;
264 : // }
265 0 : os << "[Opened " << (readOnly_p?"in readonly mode":(useScratch_p?"with scratch model column":"with virtual model column")) << "]" << LogIO::POST;
266 0 : TableExprNode exprNode=thisSelection.toTableExprNode(&thisms);
267 0 : if(!(exprNode.isNull()))
268 : {
269 0 : mss4vi_p.resize(mss4vi_p.nelements()+1, false, true);
270 0 : MeasurementSet thisMSSelected0 = MeasurementSet(thisms(exprNode));
271 :
272 0 : if(selpars.taql != "")
273 : {
274 0 : MSSelection mss0;
275 0 : mss0.setTaQLExpr(selpars.taql);
276 0 : os << "Selecting via TaQL : " << selpars.taql << " | " ;//LogIO::POST;
277 :
278 0 : TableExprNode tenWithTaQL=mss0.toTableExprNode(&thisMSSelected0);
279 0 : MeasurementSet thisMSSelected1 = MeasurementSet(thisMSSelected0(tenWithTaQL));
280 : //mss4vi_p[mss4vi_p.nelements()-1]=MeasurementSet(thisms(exprNode));
281 0 : mss4vi_p[mss4vi_p.nelements()-1]=thisMSSelected1;
282 0 : }
283 : else
284 0 : mss4vi_p[mss4vi_p.nelements()-1]=thisMSSelected0;
285 :
286 0 : os << " NRows selected : " << (mss4vi_p[mss4vi_p.nelements()-1]).nrow() << LogIO::POST;
287 0 : thisms.unlock();
288 0 : }
289 : else{
290 0 : throw(AipsError("Selection for given MS "+selpars.msname+" is invalid"));
291 : }
292 : //We should do the select channel here for the VI construction later
293 : //Need a cross check between channel selection and ms
294 : // replace below if/when viFrquencySelectionUsingChannels takes in a MSSelection
295 : // rather than the following gymnastics
296 : {
297 0 : Matrix<Int> chanlist = thisSelection.getChanList( & mss4vi_p[mss4vi_p.nelements()-1] );
298 :
299 0 : IPosition shape = chanlist.shape();
300 0 : uInt nSelections = shape[0];
301 0 : Int spw=0,chanStart=0,chanEnd=0,chanStep=1;
302 :
303 : ///////////////Temporary revert to using Vi/vb
304 0 : Int msin=mss4vi_p.nelements();
305 0 : blockNChan_p.resize(msin, false, true);
306 0 : blockStart_p.resize(msin, false, true);
307 0 : blockStep_p.resize(msin, false, true);
308 0 : blockSpw_p.resize(msin, false, true);
309 0 : msin-=1;
310 0 : blockNChan_p[msin].resize(nSelections);
311 0 : blockStart_p[msin].resize(nSelections);
312 0 : blockStep_p[msin].resize(nSelections);
313 0 : blockSpw_p[msin].resize(nSelections);
314 : ///////////////////////
315 :
316 : ////Channel selection 'flags' need for when using old VI/VB
317 : //set up Cube for storing the 'flags' for all MSes
318 : //find max no. channels from the current ms
319 0 : const MSSpWindowColumns spwc(thisms.spectralWindow());
320 0 : uInt nspw = spwc.nrow();
321 0 : const ScalarColumn<Int> spwNchans(spwc.numChan());
322 0 : Vector<Int> nchanvec = spwNchans.getColumn();
323 0 : Int maxnchan = 0;
324 0 : for (uInt i=0;i<nchanvec.nelements();i++) {
325 0 : maxnchan=max(nchanvec[i],maxnchan);
326 : }
327 0 : uInt maxnspw = 0;
328 : // msin is now zero based index
329 0 : for (Int i=0;i<msin+1;i++) {
330 0 : maxnspw=max(maxnspw,nspw);
331 : }
332 : //maxnspw=max(nspw,maxnspw);
333 : // update the channel selections and initialize to 1's (all selected)
334 0 : chanSel_p.resize(msin+1,nspw,maxnchan,true);
335 0 : chanSel_p.yzPlane(msin)=1;
336 :
337 0 : if(selpars.freqbeg==""){
338 : //re-initialize to false
339 0 : chanSel_p.yzPlane(msin)=0;
340 0 : Vector<Int> loChans(nspw, -1);
341 : /////So this gymnastic is needed
342 0 : Int nUsedSpw = 0;
343 0 : Vector<Int> chanStepPerSpw(nspw,-1);
344 0 : Vector<Int> maxChanEnd(nspw,0);
345 0 : Vector<Int> usedSpw(nspw,-1);
346 0 : for(uInt k=0; k < nSelections; ++k)
347 : {
348 :
349 0 : spw = chanlist(k,0);
350 :
351 : // channel selection
352 0 : chanStart = chanlist(k,1);
353 0 : chanEnd = chanlist(k,2);
354 0 : chanStep = chanlist(k,3);
355 : //nchan = chanEnd-chanStart+1;
356 : //nchan = Int(ceil(Double(chanEnd-chanStart+1)/Double(chanStep)));
357 0 : maxChanEnd(spw) = max(maxChanEnd(spw), chanEnd);
358 0 : chanStepPerSpw(spw) = chanStep;
359 : // find lowest selected channel for each spw
360 0 : if(loChans(spw) == -1) {
361 0 : loChans(spw) = chanStart;
362 0 : nUsedSpw++;
363 0 : usedSpw.resize(nUsedSpw,true);
364 0 : usedSpw(nUsedSpw-1) = spw;
365 : }
366 : else {
367 0 : if ( loChans(spw) > chanStart) loChans(spw) = chanStart;
368 : }
369 :
370 : //channelSelector.add (spw, chanStart, nchan,chanStep);
371 : ///////////////Temporary revert to using Vi/vb
372 : // will not work with multi-selections within a spw
373 : //blockNChan_p[msin][k]=nchan;
374 : //blockStart_p[msin][k]=chanStart;
375 : //blockStep_p[msin][k]=chanStep;
376 : //blockSpw_p[msin][k]=spw;
377 : }
378 :
379 : // vi.selectChannel() does not handle multiple chan sels within a spw???
380 : // So store a single channel sel(range) per spw... + chanflags...
381 0 : blockNChan_p[msin].resize(nUsedSpw);
382 0 : blockStart_p[msin].resize(nUsedSpw);
383 0 : blockStep_p[msin].resize(nUsedSpw);
384 0 : blockSpw_p[msin].resize(nUsedSpw);
385 0 : for(uInt j=0; j < (uInt) nUsedSpw; ++j)
386 : {
387 0 : Int ispw = usedSpw(j);
388 0 : if(loChans(ispw)!=-1)
389 : {
390 0 : blockNChan_p[msin][j] = Int(ceil(Double(maxChanEnd(ispw)-loChans(ispw)+1)/Double(chanStepPerSpw(ispw))));
391 0 : blockStart_p[msin][j] = loChans(ispw);
392 0 : blockStep_p[msin][j] = chanStepPerSpw(ispw);
393 0 : blockSpw_p[msin][j] = ispw;
394 : }
395 : }
396 :
397 0 : Int relStart=0;
398 0 : for(uInt k=0; k < nSelections; ++k)
399 : {
400 0 : spw = chanlist(k,0);
401 0 : chanStart = chanlist(k,1);
402 0 : chanEnd = chanlist(k,2);
403 :
404 : // for 'channel flags' (for old VI/VB)
405 : // MSSelect will be applied before the channel flags in FTMachine so
406 : // chanSel_p will be relative to the start chan
407 : //for (uInt k=chanStart;k<(chanEnd+1);k+=chanStep) {
408 0 : if (loChans(spw) == chanStart ) {
409 0 : relStart = 0;
410 : }
411 : else {
412 0 : relStart = chanStart - loChans(spw);
413 : }
414 0 : for (uInt k=relStart;k<(uInt)(chanEnd-loChans(spw)+1);k+=chanStep) {
415 0 : chanSel_p(msin,spw,k)=1;
416 : //cerr<<"chanselflag spw="<<spw<<" chan="<<k<<endl;
417 : }
418 : /////////////////////////////////////////
419 :
420 : }
421 :
422 0 : }
423 : else{
424 0 : Quantity freq;
425 0 : Quantity::read(freq, selpars.freqbeg);
426 0 : Double lowfreq=freq.getValue("Hz");
427 0 : Quantity::read(freq, selpars.freqend);
428 0 : Double topfreq=freq.getValue("Hz");
429 : //////////OLD VI/VB
430 : //ImagerMultiMS im;
431 0 : Vector<Int> elspw, elstart, elnchan;
432 : // Intent can be used to select a field id so check
433 : // field ids actually present in the selection-applied MS
434 0 : Vector<Int>fields=thisSelection.getFieldList( & mss4vi_p[msin]);
435 0 : MSColumns tmpmsc(mss4vi_p[msin]);
436 0 : Vector<Int> fldidv=tmpmsc.fieldId().getColumn();
437 0 : std::set<Int> ufldids(fldidv.begin(),fldidv.end());
438 0 : std::vector<Int> tmpv(ufldids.begin(), ufldids.end());
439 0 : fields.resize(tmpv.size());
440 0 : uInt count=0;
441 0 : for (std::vector<int>::const_iterator it=tmpv.begin();it != tmpv.end(); it++)
442 : {
443 0 : fields(count) = *it;
444 0 : count++;
445 : }
446 :
447 0 : Int fieldid=fields.nelements() ==0 ? 0: fields[0];
448 0 : Vector<Int> chanSteps = chanlist.column(3);
449 0 : Vector<Int> spws = chanlist.column(0);
450 0 : Int nspw = genSort( spws, Sort::Ascending, (Sort::QuickSort | Sort::NoDuplicates ));
451 0 : Vector<Int> chanStepPerSpw(nspw);
452 0 : uInt icount = 0;
453 0 : Int prevspw = -1;
454 0 : for (uInt j=0; j < spws.nelements(); j++ ) {
455 0 : if (spws[j] != prevspw) {
456 0 : chanStepPerSpw[icount] = chanSteps[j];
457 0 : prevspw = spws[j];
458 0 : icount++;
459 : }
460 : }
461 : //getSpwInFreqRange seems to assume the freqs to be the center of a channel while freqbeg and freqend appear to be
462 : //frequencies at the channel range edges. So set freqStep = 0.0 .
463 : //MSUtil::getSpwInFreqRange(elspw, elstart,elnchan,mss4vi_p[msin],lowfreq, topfreq,1.0, selpars.freqframe, fieldid);
464 0 : MSUtil::getSpwInFreqRange(elspw, elstart,elnchan,mss4vi_p[msin],lowfreq, topfreq,0.0, selpars.freqframe, fieldid);
465 : //im.adviseChanSelex(lowfreq, topfreq, 1.0, selpars.freqframe, elspw, elstart, elnchan, selpars.msname, fieldid, false);
466 : // Refine elspw if there is spw selection
467 0 : if ( elspw.nelements()==0 || elstart.nelements()==0 || elnchan.nelements()==0 )
468 0 : throw(AipsError("Failed to select vis. data by frequency range"));
469 0 : Vector<Int> tempspwlist, tempnchanlist, tempstartlist;
470 0 : uInt nselspw=0;
471 0 : if (selpars.spw != "") {
472 0 : Vector<Int> spwids = thisSelection.getSpwList( & mss4vi_p[msin]);
473 0 : for (uInt k=0; k < elspw.nelements(); k++ ) {
474 0 : if (std::find(spwids.begin(), spwids.end(), elspw[k]) != spwids.end()) {
475 0 : tempspwlist.resize(tempspwlist.nelements()+1,true);
476 0 : tempnchanlist.resize(tempnchanlist.nelements()+1,true);
477 0 : tempstartlist.resize(tempstartlist.nelements()+1,true);
478 0 : tempspwlist[nselspw] = elspw[k];
479 0 : tempnchanlist[nselspw] = elnchan[k];
480 0 : tempstartlist[nselspw] = elstart[k];
481 0 : nselspw++;
482 : }
483 : }
484 0 : elspw.resize(tempspwlist.nelements());
485 0 : elnchan.resize(tempnchanlist.nelements());
486 0 : elstart.resize(tempstartlist.nelements());
487 0 : elspw = tempspwlist;
488 0 : elnchan = tempnchanlist;
489 0 : elstart = tempstartlist;
490 0 : }
491 0 : blockNChan_p[msin].resize(elspw.nelements());
492 0 : blockStart_p[msin].resize(elspw.nelements());
493 0 : blockStep_p[msin].resize(elspw.nelements());
494 0 : blockSpw_p[msin].resize(elspw.nelements());
495 0 : blockNChan_p[msin]=elnchan;
496 0 : blockStart_p[msin]=elstart;
497 0 : blockStep_p[msin].set(1);
498 0 : blockSpw_p[msin]=elspw;
499 : //////////////////////
500 0 : }
501 0 : os << "selected spw " << blockSpw_p[msin] << " start channels " << blockStart_p[msin] << " nchannels "<< blockNChan_p[msin] << LogIO::POST;
502 :
503 0 : }
504 0 : writeAccess_p=writeAccess_p && !selpars.readonly;
505 0 : createVisSet(writeAccess_p);
506 :
507 : /////// Remove this when the new vi/vb is able to get the full freq range.
508 0 : mssFreqSel_p.resize();
509 0 : mssFreqSel_p = thisSelection.getChanFreqList(NULL,true);
510 :
511 : //// Set the data column on which to operate
512 : // TT: added checks for the requested data column existace
513 : // cout << "Using col : " << selpars.datacolumn << endl;
514 0 : if( selpars.datacolumn.contains("data") || selpars.datacolumn.contains("obs") )
515 0 : { if( thisms.tableDesc().isColumn("DATA") ) { datacol_p = FTMachine::OBSERVED; }
516 0 : else { os << LogIO::SEVERE <<"DATA column does not exist" << LogIO::EXCEPTION;}
517 : }
518 0 : else if( selpars.datacolumn.contains("corr") ) { datacol_p = FTMachine::CORRECTED; }
519 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;}
520 :
521 : /*
522 : else if( selpars.datacolumn.contains("corr") ) {
523 : if( thisms.tableDesc().isColumn("CORRECTED_DATA") ) { datacol_p = FTMachine::CORRECTED; }
524 : else
525 : {
526 : if( thisms.tableDesc().isColumn("DATA") ) {
527 : datacol_p = FTMachine::OBSERVED;
528 : os << "CORRECTED_DATA column does not exist. Using DATA column instead" << LogIO::POST;
529 : }
530 : else {
531 : os << LogIO::SEVERE <<"Neither CORRECTED_DATA nor DATA columns exist" << LogIO::EXCEPTION;
532 : }
533 : }
534 :
535 : }
536 : else { os << LogIO::WARN << "Invalid data column : " << datacol_p << ". Using corrected (or observed if corrected doesn't exist)" << LogIO::POST; datacol_p = thisms.tableDesc().isColumn("CORRECTED_DATA") ? FTMachine::CORRECTED : FTMachine::OBSERVED; }
537 : */
538 0 : dataSel_p.resize(dataSel_p.nelements()+1, true);
539 :
540 0 : dataSel_p[dataSel_p.nelements()-1]=selpars;
541 :
542 :
543 0 : unlockMSs();
544 0 : }
545 0 : catch(AipsError &x)
546 : {
547 0 : unlockMSs();
548 0 : throw( AipsError("Error in selectData() : "+x.getMesg()) );
549 0 : }
550 :
551 0 : return true;
552 :
553 0 : }
554 :
555 :
556 0 : void SynthesisImager::unlockMSs()
557 : {
558 0 : LogIO os( LogOrigin("SynthesisImager","unlockMSs",WHERE) );
559 0 : for(uInt i=0;i<mss4vi_p.nelements();i++)
560 : {
561 0 : os << LogIO::NORMAL2 << "Unlocking : " << mss4vi_p[i].tableName() << LogIO::POST;
562 0 : MeasurementSet &ms_l = const_cast<MeasurementSet& >(mss4vi_p[i]);
563 0 : ms_l.unlock();
564 0 : ms_l.antenna().unlock();
565 0 : ms_l.dataDescription().unlock();
566 0 : ms_l.feed().unlock();
567 0 : ms_l.field().unlock();
568 0 : ms_l.observation().unlock();
569 0 : ms_l.polarization().unlock();
570 0 : ms_l.processor().unlock();
571 0 : ms_l.spectralWindow().unlock();
572 0 : ms_l.state().unlock();
573 : //
574 : // Unlock the optional sub-tables as well, if they are present
575 : //
576 0 : if(!(ms_l.source().isNull())) ms_l.source().unlock();
577 0 : if(!(ms_l.doppler().isNull())) ms_l.doppler().unlock();
578 0 : if(!(ms_l.flagCmd().isNull())) ms_l.flagCmd().unlock();
579 0 : if(!(ms_l.freqOffset().isNull())) ms_l.freqOffset().unlock();
580 0 : if(!(ms_l.history().isNull())) ms_l.history().unlock();
581 0 : if(!(ms_l.pointing().isNull())) ms_l.pointing().unlock();
582 0 : if(!(ms_l.sysCal().isNull())) ms_l.sysCal().unlock();
583 0 : if(!(ms_l.weather().isNull())) ms_l.weather().unlock();
584 : }
585 0 : }
586 : ///////////////////////////////////
587 0 : bool SynthesisImager::unlockImages()
588 : {
589 :
590 0 : return itsMappers.releaseImageLocks();
591 :
592 : }
593 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
594 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
595 0 : Bool SynthesisImager::defineImage(const String& imagename, const Int nx, const Int ny,
596 : const Quantity& cellx, const Quantity& celly,
597 : const String& stokes,
598 : const MDirection& phaseCenter,
599 : const Int nchan,
600 : const Quantity&freqStart,
601 : const Quantity& freqStep,
602 : const Vector<Quantity>& restFreq,
603 : const Int facets,
604 : //s const Int chanchunks,
605 : const String ftmachine,
606 : const Int nTaylorTerms,
607 : const Quantity& refFreq,
608 : const Projection& projection,
609 : const Quantity& distance,
610 : const MFrequency::Types& freqFrame,
611 : const Bool trackSource,
612 : const MDirection& trackDir,
613 : const Bool overwrite,
614 : const Float padding,
615 : const Bool useAutocorr,
616 : const Bool useDoublePrec,
617 : const Int wprojplanes,
618 : const String convFunc,
619 : // const String vpTable,
620 : const String startmodel,
621 : // The extra params for WB-AWP
622 : const Bool aTermOn,// = true,
623 : const Bool psTermOn,// = true,
624 : const Bool mTermOn,// = false,
625 : const Bool wbAWP,// = true,
626 : const String cfCache,// = "",
627 : const Bool usePointing,// = false,
628 : const Bool doPBCorr,// = true,
629 : const Bool conjBeams,// = true,
630 : const Float computePAStep, //=360.0
631 : const Float rotatePAStep //=5.0
632 : )
633 : {
634 0 : String err("");
635 0 : LogIO os( LogOrigin("SynthesisImager","defineImage",WHERE) );
636 0 : SynthesisParamsImage impars;
637 0 : impars.imageName=imagename;
638 0 : Vector<Int> ims(2);ims[0]=nx; ims[1]=ny;
639 0 : impars.imsize=ims;
640 0 : Vector<Quantity> cells(2); cells[0]=cellx, cells[1]=celly;
641 0 : impars.cellsize=cells;
642 0 : impars.stokes=stokes;
643 0 : impars.phaseCenter=phaseCenter;
644 0 : impars.nchan=nchan;
645 0 : impars.mode= nchan < 0 ? "mfs" : "cube";
646 0 : impars.freqStart=freqStart;
647 0 : impars.freqStep=freqStep;
648 0 : impars.restFreq=restFreq;
649 0 : impars.nTaylorTerms=nTaylorTerms;
650 0 : impars.refFreq=refFreq;
651 0 : impars.projection=projection;
652 0 : impars.freqFrame=freqFrame;
653 0 : impars.overwrite=overwrite;
654 0 : impars.startModel.resize(1); impars.startModel[0]=startmodel;
655 :
656 0 : err += impars.verify();
657 :
658 0 : if(err.length() >0)
659 0 : os << LogIO::WARN << "impars verify error " << err << LogIO::POST;
660 0 : err="";
661 0 : SynthesisParamsGrid gridpars;
662 0 : gridpars.ftmachine=ftmachine;
663 0 : gridpars.distance=distance;
664 0 : gridpars.trackSource=trackSource;
665 0 : gridpars.trackDir=trackDir;
666 0 : gridpars.padding=padding;
667 0 : gridpars.facets=facets;
668 : // gridpars.chanchunks=chanchunks;
669 0 : gridpars.useAutoCorr=useAutocorr;
670 0 : gridpars.useDoublePrec=useDoublePrec;
671 0 : gridpars.wprojplanes=wprojplanes;
672 0 : gridpars.convFunc=convFunc;
673 : // gridpars.vpTable=vpTable;
674 0 : gridpars.aTermOn=aTermOn;
675 0 : gridpars.psTermOn=psTermOn;
676 0 : gridpars.mTermOn=mTermOn;
677 0 : gridpars.wbAWP=wbAWP;
678 0 : gridpars.cfCache=cfCache;
679 0 : gridpars.usePointing=usePointing;
680 0 : gridpars.doPBCorr=doPBCorr;
681 0 : gridpars.conjBeams=conjBeams;
682 0 : gridpars.computePAStep=computePAStep;
683 0 : gridpars.rotatePAStep=rotatePAStep;
684 0 : gridpars.imageName="YAKYAK"; // no clue why gridpars has imageName..verify needs it
685 0 : err += gridpars.verify();
686 0 : if(err.length() >0)
687 0 : os << LogIO::WARN << "gridpars verify error " << err << LogIO::POST;
688 : //if( err.length()>0 ) throw(AipsError("Invalid Image/Gridding parameters : " + err));
689 :
690 0 : defineImage( impars, gridpars );
691 :
692 0 : return true;
693 0 : }
694 :
695 0 : Bool SynthesisImager::defineImage(SynthesisParamsImage& impars,
696 : const SynthesisParamsGrid& gridpars)
697 : {
698 :
699 0 : LogIO os( LogOrigin("SynthesisImager","defineImage",WHERE) );
700 0 : if(mss4vi_p.nelements() ==0)
701 0 : os << "SelectData has to be run before defineImage" << LogIO::EXCEPTION;
702 :
703 0 : CoordinateSystem csys;
704 0 : CountedPtr<FTMachine> ftm, iftm;
705 :
706 0 : impars_p = impars;
707 0 : gridpars_p = gridpars;
708 : try
709 : {
710 :
711 0 : os << "Define image coordinates for [" << impars.imageName << "] : " << LogIO::POST;
712 :
713 0 : csys = impars.buildCoordinateSystem( rvi_p );
714 0 : IPosition imshape = impars.shp();
715 :
716 0 : os << "Impars : start " << impars.start << LogIO::POST;
717 0 : os << "Shape : " << imshape << "Spectral : " << csys.spectralCoordinate().referenceValue() << " at " << csys.spectralCoordinate().referencePixel() << " with increment " << csys.spectralCoordinate().increment() << LogIO::POST;
718 :
719 0 : if( (itsMappers.nMappers()==0) ||
720 0 : (impars.imsize[0]*impars.imsize[1] > itsMaxShape[0]*itsMaxShape[1]))
721 : {
722 0 : itsMaxShape=imshape;
723 0 : itsMaxCoordSys=csys;
724 : }
725 0 : itsNchan = imshape[3];
726 0 : itsCsysRec = impars.getcsys();
727 : /*
728 : os << "Define image [" << impars.imageName << "] : nchan : " << impars.nchan
729 : //<< ", freqstart:" << impars.freqStart.getValue() << impars.freqStart.getUnit()
730 : << ", start:" << impars.start
731 : << ", imsize:" << impars.imsize
732 : << ", cellsize: [" << impars.cellsize[0].getValue() << impars.cellsize[0].getUnit()
733 : << " , " << impars.cellsize[1].getValue() << impars.cellsize[1].getUnit()
734 : << LogIO::POST;
735 : */
736 0 : }
737 0 : catch(AipsError &x)
738 : {
739 0 : os << "Error in building Coordinate System and Image Shape : " << x.getMesg() << LogIO::EXCEPTION;
740 0 : }
741 :
742 :
743 : try
744 : {
745 0 : os << "Set Gridding options for [" << impars.imageName << "] with ftmachine : " << gridpars.ftmachine << LogIO::POST;
746 :
747 0 : itsVpTable=gridpars.vpTable;
748 0 : itsMakeVP= ( gridpars.ftmachine.contains("mosaicft") ||
749 0 : gridpars.ftmachine.contains("awprojectft") )?False:True;
750 :
751 0 : createFTMachine(ftm, iftm, gridpars.ftmachine, impars.nTaylorTerms, gridpars.mType,
752 0 : gridpars.facets, gridpars.wprojplanes,
753 0 : gridpars.padding,gridpars.useAutoCorr,gridpars.useDoublePrec,
754 0 : gridpars.convFunc,
755 0 : gridpars.aTermOn,gridpars.psTermOn, gridpars.mTermOn,
756 0 : gridpars.wbAWP,gridpars.cfCache,gridpars.usePointing,
757 0 : gridpars.doPBCorr,gridpars.conjBeams,
758 0 : gridpars.computePAStep,gridpars.rotatePAStep,
759 0 : gridpars.interpolation, impars.freqFrameValid, 1000000000, 16, impars.stokes,
760 0 : impars.imageName);
761 :
762 : }
763 0 : catch(AipsError &x)
764 : {
765 0 : os << "Error in setting up FTMachine() : " << x.getMesg() << LogIO::EXCEPTION;
766 0 : }
767 :
768 : try
769 : {
770 0 : appendToMapperList(impars.imageName, csys, impars.shp(),
771 : ftm, iftm,
772 0 : gridpars.distance, gridpars.facets, gridpars.chanchunks,impars.overwrite,
773 0 : gridpars.mType, gridpars.padding, impars.nTaylorTerms, impars.startModel);
774 0 : imageDefined_p=true;
775 : }
776 0 : catch(AipsError &x)
777 : {
778 0 : os << "Error in adding Mapper : "+x.getMesg() << LogIO::EXCEPTION;
779 0 : }
780 :
781 0 : return true;
782 0 : }
783 :
784 17 : Long SynthesisImager::estimateRAM(){
785 17 : Long mem=0;
786 34 : LogIO os( LogOrigin("SynthesisImager","estimateRAM",WHERE) );
787 17 : if(itsMappers.nMappers()==0)
788 0 : os << "SynthesisImage has not been setup to get an estimate of memory usage" << LogIO::WARN << LogIO::POST;
789 17 : mem=itsMappers.estimateRAM();
790 17 : return mem;
791 17 : }
792 :
793 0 : Bool SynthesisImager::defineImage(CountedPtr<SIImageStore> imstor,
794 : const String& ftmachine)
795 : {
796 0 : CountedPtr<FTMachine> ftm, iftm;
797 :
798 : // The following call to createFTMachine() uses the
799 : // following defaults
800 : //
801 : // facets=1, wprojplane=1, padding=1.0, useAutocorr=false,
802 : // useDoublePrec=true, gridFunction=String("SF")
803 : //
804 0 : createFTMachine(ftm, iftm, ftmachine);
805 :
806 0 : Int id=itsMappers.nMappers();
807 0 : CoordinateSystem csys =imstor->residual()->coordinates();
808 0 : IPosition imshape=imstor->residual()->shape();
809 0 : Int nx=imshape[0], ny=imshape[1];
810 0 : if( (id==0) || (nx*ny > itsMaxShape[0]*itsMaxShape[1]))
811 : {
812 0 : itsMaxShape=imshape;
813 0 : itsMaxCoordSys=csys;
814 : }
815 :
816 0 : itsMappers.addMapper( createSIMapper( "default", imstor, ftm, iftm, id ) );
817 :
818 0 : return true;
819 0 : }
820 :
821 : //////////////////////////////////////////////////////////////////////////////////////////////////////
822 : /////////////////////////////////////////////////////////////////////////////////////////////////////
823 : // Record SynthesisImage::runGridParamsHeuristics(SynthesisParamsGrid& gridPars)
824 : // {
825 : // return gridPars;
826 : // }
827 : //////////////////////////////////////////////////////////////////////////////////////////////////////
828 : /////////////////////////////////////////////////////////////////////////////////////////////////////
829 :
830 0 : Vector<SynthesisParamsSelect> SynthesisImager::tuneSelectData(){
831 0 : LogIO os( LogOrigin("SynthesisImager","tuneSelectData",WHERE) );
832 0 : if(itsMappers.nMappers() < 1)
833 0 : ThrowCc("defineimage has to be run before tuneSelectData");
834 0 : if(impars_p.mode=="cubesource" || impars_p.mode=="cubedata")
835 0 : return dataSel_p;
836 0 : os << "Tuning frequency data selection to match image spectral coordinates" << LogIO::POST;
837 :
838 0 : Vector<SynthesisParamsSelect> origDatSel(dataSel_p.nelements());
839 0 : origDatSel=dataSel_p;
840 : /*Record selpars;
841 : for(uInt k=0; k < origDatSel.nelements(); ++k){
842 : Record inSelRec=origDatSel[k].toRecord();
843 : selpars.defineRecord("ms"+String::toString(k), inSelRec);
844 : }
845 : */
846 0 : Int nchannel=itsMaxShape[3];
847 0 : CoordinateSystem cs=itsMaxCoordSys;
848 0 : MFrequency::Types freqframe=cs.spectralCoordinate(cs.findCoordinate(Coordinate::SPECTRAL)).frequencySystem(False);
849 0 : if(freqframe != MFrequency::REST && freqframe != MFrequency::Undefined)
850 0 : cs.setSpectralConversion("LSRK");
851 0 : Vector<Double> pix(4);
852 : //increase edge part a bit vi/vb2 is making lots of grief for selecting properly if the channel are negative
853 : // specially...
854 0 : pix[0]=0; pix[1]=0; pix[2]=0; pix[3]=-1.5;
855 0 : Double freq1=cs.toWorld(pix)[3];
856 0 : pix[3]=Double(nchannel)+0.5;
857 0 : Double freq2=cs.toWorld(pix)[3];
858 0 : String units=cs.worldAxisUnits()[3];
859 0 : if(freq2 < freq1){
860 0 : Double tmp=freq1;
861 0 : freq1=freq2;
862 0 : freq2=tmp;
863 : }
864 : //String freqbeg=String::toString(freq1)+units;
865 : //String freqend=String::toString(freq2)+units;
866 : // String::toString drops the digits for double precision
867 0 : String freqbeg=doubleToString(freq1)+units;
868 0 : String freqend=doubleToString(freq2)+units;
869 : //Record outRec=SynthesisUtilMethods::cubeDataPartition(selpars, 1, freq1, freq2);
870 : //Record partRec=outRec.asRecord("0");
871 :
872 0 : if(mss_p.nelements() >0){
873 0 : for (uInt k=0; k < mss_p.nelements(); ++k){
874 0 : if(mss_p[k])
875 0 : delete mss_p[k];
876 : }
877 0 : mss_p.resize(0, true, false);
878 : }
879 : ///resetting the block ms
880 0 : mss4vi_p.resize(0,true, false);
881 : //resetting data selection stored
882 0 : dataSel_p.resize();
883 : // reset rvi_p
884 0 : if(rvi_p) delete rvi_p;
885 0 : rvi_p=NULL;
886 :
887 0 : for(uInt k=0; k< origDatSel.nelements(); ++k){
888 0 : SynthesisParamsSelect outsel=origDatSel[k];
889 0 : outsel.freqbeg=freqbeg;
890 0 : outsel.freqend=freqend;
891 0 : outsel.freqframe=MFrequency::LSRK;
892 0 : selectData(outsel);
893 0 : }
894 0 : return dataSel_p;
895 :
896 0 : }
897 :
898 :
899 : ///////////////////////////////////////////////////////////////////////////////////////////
900 : ///////////////////////////////////////////////////////////////////////////////////////////
901 0 : void SynthesisImager::setComponentList(const ComponentList& cl, Bool sdgrid){
902 0 : String cft="SimpleComponentFTMachine";
903 0 : if(sdgrid)
904 0 : cft="SimpCompGridFTMachine";
905 0 : CountedPtr<SIMapper> sm=new SIMapper(cl, cft);
906 0 : itsMappers.addMapper(sm);
907 : ////itsMappers.addMapper( createSIMapper( mappertype, imstor, ftm, iftm, id) );
908 :
909 0 : }
910 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
911 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
912 :
913 : //////////////////////Reset the Mapper
914 : ////////////////////
915 0 : void SynthesisImager::resetMappers(){
916 : ////reset code
917 0 : itsMappers=SIMapperCollection();
918 0 : unFacettedImStore_p=NULL;
919 0 : unChanChunkedImStore_p=NULL;
920 0 : }
921 : //////////////////////////////////////////////////////////////////
922 : /////////////////////////////////////////////
923 0 : CountedPtr<SIImageStore> SynthesisImager::imageStore(const Int id)
924 : {
925 0 : AlwaysAssert( ! ( facetsStore_p>1 && chanChunksStore_p>1 ) , AipsError);
926 :
927 0 : if(facetsStore_p >1)
928 : {
929 0 : if(id==0)
930 : {
931 0 : return unFacettedImStore_p;
932 : }
933 : else
934 : {
935 0 : return itsMappers.imageStore(facetsStore_p*facetsStore_p+id-1);
936 : }
937 : }
938 :
939 0 : if(chanChunksStore_p >1)
940 : {
941 0 : if(id==0)
942 : {
943 0 : return unChanChunkedImStore_p;
944 : }
945 : else
946 : {
947 0 : return itsMappers.imageStore(chanChunksStore_p+id-1);
948 : }
949 : }
950 :
951 :
952 0 : return itsMappers.imageStore(id);
953 : }
954 :
955 :
956 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
957 : ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
958 :
959 17 : Record SynthesisImager::executeMajorCycle(const Record& controlRecord)
960 : {
961 34 : LogIO os( LogOrigin("SynthesisImager","executeMajorCycle",WHERE) );
962 :
963 17 : nMajorCycles++;
964 17 : if(controlRecord.isDefined("nmajorcycles"))
965 0 : controlRecord.get("nmajorcycles", nMajorCycles);
966 17 : Record outRec=controlRecord;
967 17 : Bool lastcycle=false;
968 :
969 :
970 17 : if( controlRecord.isDefined("lastcycle") )
971 : {
972 17 : controlRecord.get( "lastcycle" , lastcycle );
973 : //cout << "lastcycle : " << lastcycle << endl;
974 : }
975 : //else {cout << "No lastcycle" << endl;}
976 :
977 17 : os << "----------------------------------------------------------- Run ";
978 : //if (lastcycle) os << "(Last) " ;
979 17 : os << "Major Cycle " << nMajorCycles << " -------------------------------------" << LogIO::POST;
980 :
981 : try
982 : {
983 17 : if( (itsMaxShape[3] > 1 || impars_p.mode.contains("cube"))&& doingCubeGridding_p ){/// and valid ftmachines
984 0 : runMajorCycleCube(false, controlRecord);
985 : }
986 : else{
987 17 : if( itsDataLoopPerMapper == false )
988 17 : { runMajorCycle(false, lastcycle);}
989 : else
990 0 : { runMajorCycle2(false, lastcycle);}
991 :
992 : }
993 17 : if(lastcycle){
994 17 : String mess="In Major Cycle";
995 17 : if(controlRecord.isDefined("usemask") && controlRecord.asString("usemask").contains("auto")){
996 0 : mess="\nFor Automasking most major cycles may appear wrongly as the last one ";
997 : }
998 :
999 17 : std::vector<String> tmpfiles=itsMappers.cleanupTempFiles(mess);
1000 17 : outRec.define("tempfilenames", Vector<String>(tmpfiles));
1001 17 : }
1002 17 : itsMappers.releaseImageLocks();
1003 :
1004 : }
1005 0 : catch(AipsError &x)
1006 : {
1007 0 : throw( AipsError("Error in running Major Cycle : "+x.getMesg()) );
1008 0 : }
1009 34 : return outRec;
1010 17 : }// end of executeMajorCycle
1011 : //////////////////////////////////////////////
1012 : /////////////////////////////////////////////
1013 17 : void SynthesisImager::cleanupTempFiles(){
1014 34 : LogIO os( LogOrigin("SynthesisImager","cleanupTempFiles",WHERE) );
1015 17 : for (auto it = tempFileNames_p.begin(); it != tempFileNames_p.end(); ++it) {
1016 : //cerr <<"Trying to cleanup " << (*it) << endl;
1017 0 : if(Table::isReadable(*it)){
1018 : try{
1019 0 : TableUtil::deleteTable(*it);
1020 : }
1021 0 : catch(AipsError &x){
1022 0 : os << LogIO::WARN<< "YOU may have to delete the temporary file " << *it << " because " << x.getMesg() << LogIO::POST;
1023 0 : }
1024 : }
1025 :
1026 :
1027 : }
1028 :
1029 :
1030 17 : }
1031 17 : Record SynthesisImager::makePSF()
1032 : {
1033 17 : Record outRec=Record();
1034 34 : LogIO os( LogOrigin("SynthesisImager","makePSF",WHERE) );
1035 :
1036 17 : os << "----------------------------------------------------------- Make PSF ---------------------------------------------" << LogIO::POST;
1037 :
1038 : try
1039 : {
1040 17 : if( (itsMaxShape[3] >1 || impars_p.mode.contains("cube")) && doingCubeGridding_p){///and valid ftmachines
1041 0 : runMajorCycleCube(true);
1042 : }
1043 : else{
1044 17 : if( itsDataLoopPerMapper == false )
1045 17 : {runMajorCycle(true, false);}
1046 : else
1047 0 : {runMajorCycle2(true, false);}
1048 : }
1049 : // makeImage();
1050 :
1051 17 : String mess("PSF wieghting scheme");
1052 17 : std::vector<String> tmpfiles=itsMappers.cleanupTempFiles(mess);
1053 17 : outRec.define("tempfilenames", Vector<String>(tmpfiles));
1054 17 : itsMappers.releaseImageLocks();
1055 :
1056 17 : }
1057 0 : catch(AipsError &x)
1058 : {
1059 0 : throw( AipsError("Error in making PSF : "+x.getMesg()) );
1060 0 : }
1061 34 : return outRec;
1062 17 : }
1063 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1064 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1065 :
1066 0 : Record SynthesisImager::apparentSensitivity()
1067 : {
1068 :
1069 0 : throw( AipsError("apparentSensitivity calculation not supported in SynthesisImager (VI1)") );
1070 : return Record();
1071 :
1072 : }
1073 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1074 :
1075 :
1076 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1077 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1078 0 : Bool SynthesisImager::weight(const String& type, const String& rmode,
1079 : const Quantity& noise, const Double robust,
1080 : const Quantity& fieldofview,
1081 : const Int npixels, const Bool multiField,
1082 : const Bool /*useCubeBriggs*/,
1083 : const String& filtertype, const Quantity& filterbmaj,
1084 : const Quantity& filterbmin, const Quantity& filterbpa, Double /*fracBW*/ )
1085 : {
1086 0 : LogIO os(LogOrigin("SynthesisImager", "weight()", WHERE));
1087 :
1088 : try {
1089 : //Int nx=itsMaxShape[0];
1090 : //Int ny=itsMaxShape[1];
1091 0 : Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
1092 0 : Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
1093 : os << LogIO::NORMAL // Loglevel INFO
1094 0 : << "Set imaging weights : " ; //<< LogIO::POST;
1095 :
1096 0 : if (type=="natural") {
1097 : os << LogIO::NORMAL // Loglevel INFO
1098 0 : << "Natural weighting" << LogIO::POST;
1099 0 : imwgt_p=VisImagingWeight("natural");
1100 : }
1101 0 : else if (type=="radial") {
1102 0 : os << "Radial weighting" << LogIO::POST;
1103 0 : imwgt_p=VisImagingWeight("radial");
1104 : }
1105 : else{
1106 0 : if(!imageDefined_p)
1107 0 : throw(AipsError("Need to define image"));
1108 0 : Int nx=itsMaxShape[0];
1109 0 : Int ny=itsMaxShape[1];
1110 0 : Quantity cellx=Quantity(itsMaxCoordSys.increment()[0], itsMaxCoordSys.worldAxisUnits()[0]);
1111 0 : Quantity celly=Quantity(itsMaxCoordSys.increment()[1], itsMaxCoordSys.worldAxisUnits()[1]);
1112 0 : if(type=="superuniform"){
1113 0 : if(!imageDefined_p) throw(AipsError("Please define image first"));
1114 0 : Int actualNpix=npixels;
1115 0 : if(actualNpix <=0)
1116 0 : actualNpix=3;
1117 : os << LogIO::NORMAL // Loglevel INFO
1118 : << "SuperUniform weighting over a square cell spanning ["
1119 : << -actualNpix
1120 0 : << ", " << actualNpix << "] in the uv plane" << LogIO::POST;
1121 0 : imwgt_p=VisImagingWeight(*rvi_p, rmode, noise, robust, nx,
1122 : ny, cellx, celly, actualNpix,
1123 0 : actualNpix, multiField);
1124 : }
1125 0 : else if ((type=="robust")||(type=="uniform")||(type=="briggs")) {
1126 0 : if(!imageDefined_p) throw(AipsError("Please define image first"));
1127 0 : Quantity actualFieldOfView_x(fieldofview), actualFieldOfView_y(fieldofview) ;
1128 0 : Int actualNPixels_x(npixels),actualNPixels_y(npixels) ;
1129 0 : String wtype;
1130 0 : if(type=="briggs") {
1131 0 : wtype = "Briggs";
1132 : }
1133 : else {
1134 0 : wtype = "Uniform";
1135 : }
1136 0 : if(actualFieldOfView_x.get().getValue()==0.0&&actualNPixels_x==0) {
1137 0 : actualNPixels_x=nx;
1138 0 : actualFieldOfView_x=Quantity(actualNPixels_x*cellx.get("rad").getValue(),"rad");
1139 0 : actualNPixels_y=ny;
1140 0 : actualFieldOfView_y=Quantity(actualNPixels_y*celly.get("rad").getValue(),"rad");
1141 : os << LogIO::NORMAL // Loglevel INFO
1142 : << wtype
1143 : << " weighting: sidelobes will be suppressed over full image"
1144 0 : << LogIO::POST;
1145 : }
1146 0 : else if(actualFieldOfView_x.get().getValue()>0.0&&actualNPixels_x==0) {
1147 0 : actualNPixels_x=nx;
1148 0 : actualNPixels_y=ny;
1149 : os << LogIO::NORMAL // Loglevel INFO
1150 : << wtype
1151 : << " weighting: sidelobes will be suppressed over specified field of view: "
1152 0 : << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by "
1153 0 : << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
1154 : }
1155 0 : else if(actualFieldOfView_x.get().getValue()==0.0&&actualNPixels_x>0) {
1156 0 : actualFieldOfView_x=Quantity(actualNPixels_x*cellx.get("rad").getValue(),"rad");
1157 0 : actualFieldOfView_y=Quantity(actualNPixels_y*celly.get("rad").getValue(),"rad");
1158 : os << LogIO::NORMAL // Loglevel INFO
1159 : << wtype
1160 : << " weighting: sidelobes will be suppressed over full image field of view: "
1161 0 : << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by "
1162 0 : << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
1163 : }
1164 : else {
1165 : os << LogIO::NORMAL // Loglevel INFO
1166 : << wtype
1167 : << " weighting: sidelobes will be suppressed over specified field of view: "
1168 0 : << actualFieldOfView_x.get("arcsec").getValue() << " arcsec by "
1169 0 : << actualFieldOfView_y.get("arcsec").getValue() << " arcsec" << LogIO::POST;
1170 : }
1171 : os << LogIO::DEBUG1
1172 : << "Weighting used " << actualNPixels_x << " by " << actualNPixels_y << " uv pixels."
1173 0 : << LogIO::POST;
1174 0 : Quantity actualCellSize_x(actualFieldOfView_x.get("rad").getValue()/actualNPixels_x, "rad");
1175 0 : Quantity actualCellSize_y(actualFieldOfView_y.get("rad").getValue()/actualNPixels_y, "rad");
1176 :
1177 : // cerr << "rmode " << rmode << " noise " << noise << " robust " << robust << " npixels " << actualNPixels << " cellsize " << actualCellSize << " multifield " << multiField << endl;
1178 : // Timer timer;
1179 : //timer.mark();
1180 : //Construct imwgt_p with old vi for now if old vi is in use as constructing with vi2 is slower
1181 :
1182 :
1183 0 : imwgt_p=VisImagingWeight(*rvi_p, wtype=="Uniform" ? "none" : rmode, noise, robust,
1184 : actualNPixels_x, actualNPixels_y, actualCellSize_x,
1185 0 : actualCellSize_y, 0, 0, multiField);
1186 :
1187 : /*
1188 : if(rvi_p !=NULL){
1189 : imwgt_p=VisImagingWeight(*rvi_p, rmode, noise, robust,
1190 : actualNPixels, actualNPixels, actualCellSize,
1191 : actualCellSize, 0, 0, multiField);
1192 : }
1193 : else{
1194 : ////This is slower by orders of magnitude as of 2014/06/25
1195 : imwgt_p=VisImagingWeight(*vi_p, rmode, noise, robust,
1196 : actualNPixels, actualNPixels, actualCellSize,
1197 : actualCellSize, 0, 0, multiField);
1198 : }
1199 : */
1200 : //timer.show("After making visweight ");
1201 :
1202 0 : }
1203 : else {
1204 : //this->unlock();
1205 : os << LogIO::SEVERE << "Unknown weighting " << type
1206 0 : << LogIO::EXCEPTION;
1207 0 : return false;
1208 : }
1209 0 : }
1210 :
1211 : //// UV-Tapering
1212 : //cout << "Taper type : " << filtertype << " : " << (filtertype=="gaussian") << endl;
1213 0 : if( filtertype == "gaussian" ) {
1214 : // os << "Setting uv-taper" << LogIO::POST;
1215 0 : imwgt_p.setFilter( filtertype, filterbmaj, filterbmin, filterbpa );
1216 : }
1217 0 : rvi_p->useImagingWeight(imwgt_p);
1218 : ///////////////////////////////
1219 :
1220 0 : SynthesisUtilMethods::getResource("Set Weighting");
1221 :
1222 : /// return true;
1223 :
1224 0 : }
1225 0 : catch(AipsError &x)
1226 : {
1227 0 : throw( AipsError("Error in Weighting : "+x.getMesg()) );
1228 0 : }
1229 :
1230 0 : return true;
1231 0 : }
1232 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1233 :
1234 : //// Get/Set Weight Grid.... write to disk and read
1235 :
1236 : /// todo : do for full mapper list, and taylor terms.
1237 0 : String SynthesisImager::getWeightDensity( )
1238 : {
1239 0 : String outname("");
1240 0 : LogIO os(LogOrigin("SynthesisImager", "getWeightDensity()", WHERE));
1241 : try
1242 : {
1243 : //if
1244 0 : if(imwgt_p.getType() != "uniform"){
1245 0 : return outname;
1246 : }
1247 0 : IPosition newshape;
1248 0 : Vector<Int> shpOfGrid=imwgt_p.shapeOfdensityGrid();
1249 0 : if(shpOfGrid(2) > 1){
1250 0 : newshape=IPosition(5,shpOfGrid[0], shpOfGrid[1],1,1,shpOfGrid[2]);
1251 : }
1252 : else{
1253 0 : newshape=IPosition(4,shpOfGrid[0], shpOfGrid[1],1,1);
1254 : }
1255 0 : IPosition where= IPosition(Vector<Int>((itsMappers.imageStore(0)->gridwt(newshape))->shape().nelements(),0));
1256 0 : if ( shpOfGrid[2] > 0 )
1257 : {
1258 0 : imwgt_p.toImageInterface(*(itsMappers.imageStore(0)->gridwt(newshape))); }
1259 0 : outname=itsMappers.imageStore(0)->gridwt()->name();
1260 0 : itsMappers.releaseImageLocks();
1261 :
1262 0 : }
1263 0 : catch (AipsError &x)
1264 : {
1265 0 : throw(AipsError("In getWeightDensity : "+x.getMesg()));
1266 0 : }
1267 0 : return outname;
1268 0 : }
1269 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1270 : /// todo : do for full mapper list, and taylor terms.
1271 :
1272 0 : Bool SynthesisImager::setWeightDensity(const String& weightimagename)
1273 : {
1274 0 : LogIO os(LogOrigin("SynthesisImager", "setWeightDensity()", WHERE));
1275 : try
1276 : {
1277 :
1278 : //Array<Float> arr;
1279 : ///Use image 0 for weight density shape
1280 : //itsMappers.imageStore(0)->gridwt()->get(arr,true);
1281 0 : if(weightimagename.size() !=0){
1282 0 : Table::isReadable(weightimagename, True);
1283 0 : PagedImage<Float> im(weightimagename);
1284 0 : imwgt_p=VisImagingWeight(im);
1285 0 : }
1286 : else{
1287 0 : imwgt_p=VisImagingWeight(*(itsMappers.imageStore(0)->gridwt()));
1288 :
1289 : }
1290 0 : rvi_p->useImagingWeight(imwgt_p);
1291 0 : itsMappers.releaseImageLocks();
1292 :
1293 : }
1294 0 : catch (AipsError &x)
1295 : {
1296 0 : throw(AipsError("In setWeightDensity : "+x.getMesg()));
1297 0 : }
1298 0 : return true;
1299 0 : }
1300 :
1301 : //////////////////////////////////////////////////////////////////////////////
1302 : //////////////////////////////////////////////////////////////////////////////
1303 :
1304 :
1305 : //////////////////////////////////////////////////////////////////////////////
1306 : //// Internal Functions start here. These are not visible to the tool layer.
1307 : //////////////////////////////////////////////////////////////////////////////
1308 :
1309 : // This should be called at each defineimage
1310 17 : CountedPtr<SIImageStore> SynthesisImager::createIMStore(
1311 : String imageName,
1312 : CoordinateSystem& cSys,
1313 : IPosition imShape,
1314 : const Bool overwrite,
1315 : MSColumns& msc,
1316 : String mappertype,
1317 : uInt ntaylorterms,
1318 : Quantity distance,
1319 : const TcleanProcessingInfo& procInfo,
1320 : uInt facets,
1321 : Bool useweightimage,
1322 : const Vector<String> &startmodel,
1323 : const Bool makeSingleDishStore)
1324 : {
1325 34 : LogIO os( LogOrigin("SynthesisImager","createIMStore",WHERE) );
1326 :
1327 17 : CountedPtr<SIImageStore> imstor;
1328 :
1329 : try {
1330 : // Prepare miscellaneous image information
1331 17 : auto objectName = msc.field().name()(msc.fieldId()(0));
1332 : // Misc info for ImageStore. This will go to the 'miscinfo' table keyword
1333 17 : Record miscInfo;
1334 17 : auto telescop=msc.observation().telescopeName()(0);
1335 17 : miscInfo.define("INSTRUME", telescop);
1336 17 : miscInfo.define("distance", distance.get("m").getValue());
1337 17 : miscInfo.define("mpiprocs", procInfo.mpiprocs);
1338 17 : miscInfo.define("chnchnks", procInfo.chnchnks);
1339 17 : miscInfo.define("memreq", procInfo.memreq);
1340 17 : miscInfo.define("memavail", procInfo.memavail);
1341 :
1342 17 : if (mappertype == "default" or mappertype == "imagemosaic") {
1343 51 : imstor = std::make_shared<SIImageStore>(
1344 : imageName, cSys, imShape, objectName,
1345 : miscInfo, overwrite,
1346 34 : (useweightimage or mappertype == "imagemosaic"),
1347 : makeSingleDishStore
1348 17 : );
1349 : }
1350 0 : else if (mappertype == "multiterm") {
1351 : // Upcast with shared_ptr and then ...
1352 : std::shared_ptr<SIImageStore>
1353 0 : multiTermStore = std::make_shared<SIImageStoreMultiTerm>(
1354 : imageName, cSys, imShape,
1355 : objectName, miscInfo, facets,
1356 : overwrite, ntaylorterms, useweightimage
1357 0 : );
1358 : // ... assign to CountedPtr<SIImageStore>
1359 0 : imstor = multiTermStore;
1360 0 : }
1361 : else { // imagemosaic currently not supported
1362 0 : throw(
1363 : AipsError(
1364 : "Internal Error : Invalid mapper type in"
1365 : " SynthesisImager::createIMStore"
1366 : )
1367 0 : );
1368 : }
1369 :
1370 : // Get polRep from 'msc' here, and send to imstore.
1371 17 : StokesImageUtil::PolRep polRep(StokesImageUtil::CIRCULAR);
1372 17 : Vector<String> polType=msc.feed().polarizationType()(0);
1373 51 : if (polType(0) != "X" and polType(0) != "Y" and
1374 51 : polType(0) != "R" and polType(0) != "L") {
1375 : os << LogIO::WARN << "Unknown stokes types in feed table: ["
1376 0 : << polType(0) << ", " << polType(1) << "]" << endl
1377 0 : << "Results open to question!" << LogIO::POST;
1378 : }
1379 :
1380 17 : if (polType(0) == "X" or polType(0) == "Y") {
1381 0 : polRep=StokesImageUtil::LINEAR;
1382 : os << LogIO::DEBUG1
1383 : << "Preferred polarization representation is linear"
1384 0 : << LogIO::POST;
1385 : }
1386 : else {
1387 17 : polRep=StokesImageUtil::CIRCULAR;
1388 : os << LogIO::DEBUG1
1389 : << "Preferred polarization representation is circular"
1390 17 : << LogIO::POST;
1391 : }
1392 : // End of reading polRep info
1393 :
1394 : // Send this info into ImageStore.
1395 17 : imstor->setDataPolFrame(polRep);
1396 :
1397 : // Set Starting model if it exists.
1398 : // cout << "In SI, set starting model to : " << startmodel << endl;
1399 17 : if (startmodel.nelements() > 0) {
1400 0 : imstor->setModelImage(startmodel);
1401 : }
1402 :
1403 17 : imstor->releaseLocks();
1404 :
1405 17 : }
1406 0 : catch (AipsError &x) {
1407 0 : throw(
1408 : AipsError(
1409 0 : "Error in createImStore : " + x.getMesg()
1410 : )
1411 0 : );
1412 0 : }
1413 :
1414 34 : return imstor;
1415 17 : }
1416 :
1417 0 : CountedPtr<SIMapper> SynthesisImager::createSIMapper(String mappertype,
1418 : CountedPtr<SIImageStore> imagestore,
1419 : CountedPtr<FTMachine> ftmachine,
1420 : CountedPtr<FTMachine> iftmachine,
1421 : uInt /*ntaylorterms*/)
1422 : {
1423 0 : LogIO os( LogOrigin("SynthesisImager","createSIMapper",WHERE) );
1424 :
1425 0 : CountedPtr<SIMapper> localMapper;
1426 :
1427 : try
1428 : {
1429 :
1430 0 : if( mappertype == "default" || mappertype == "multiterm" )
1431 : {
1432 0 : localMapper = new SIMapper( imagestore, ftmachine, iftmachine );
1433 : }
1434 0 : else if( mappertype == "imagemosaic") // || mappertype == "mtimagemosaic" )
1435 : {
1436 0 : localMapper = new SIMapperImageMosaic( imagestore, ftmachine, iftmachine );
1437 : }
1438 : else
1439 : {
1440 0 : throw(AipsError("Unknown mapper type : " + mappertype));
1441 : }
1442 :
1443 : }
1444 0 : catch(AipsError &x) {
1445 0 : throw(AipsError("Error in createSIMapper : " + x.getMesg() ) );
1446 0 : }
1447 0 : return localMapper;
1448 0 : }
1449 :
1450 :
1451 0 : Block<CountedPtr<SIImageStore> > SynthesisImager::createFacetImageStoreList(
1452 : CountedPtr<SIImageStore> imagestore,
1453 : Int facets)
1454 : {
1455 0 : LogIO os(LogOrigin("SynthesisImager", "createFacetImageStoreList"));
1456 :
1457 0 : os << "Creating " << facets*facets << " facets in total " << LogIO::POST;
1458 :
1459 0 : Block<CountedPtr<SIImageStore> > facetList( facets*facets );
1460 :
1461 0 : if( facets==1 ) { facetList[0] = imagestore; return facetList; }
1462 :
1463 : // Remember, only the FIRST field in each run can have facets. So, check for this.
1464 0 : if( ! unFacettedImStore_p.null() ) {
1465 0 : throw( AipsError("A facetted image has already been set. Facets are supported only for the main (first) field. Please submit a feature-request if you need multiple facets for outlier fields as well. ") );
1466 : }
1467 :
1468 0 : unFacettedImStore_p = imagestore;
1469 0 : facetsStore_p = facets;
1470 :
1471 : // Note : facets : Number of facets on a side.
1472 : // Note : facet : index from range(0, facets*facets)
1473 0 : for (Int facet=0; facet< facets*facets; ++facet){
1474 0 : facetList[facet] = unFacettedImStore_p->getSubImageStore(facet, facets);
1475 : }
1476 :
1477 0 : return facetList;
1478 0 : }
1479 :
1480 : /*
1481 : void SynthesisImager::setPsfFromOneFacet()
1482 : {
1483 :
1484 : if( unFacettedImStore_p.null() ){
1485 : throw(AipsError("Internal Error in SynthesisImager : Setting PSF on Null unfacettedimage"));
1486 : }
1487 :
1488 : // Copy the PSF from one facet to the center of the full image, to use for the minor cycle
1489 : //
1490 : // This code segment will work for single and multi-term
1491 : // - for single term, the index will always be 0, and SIImageStore's access functions know this.
1492 : // - for multi-term, the index will be the taylor index and SIImageStoreMultiTerm knows this.
1493 : {
1494 : IPosition shape=(unFacettedImStore_p->psf(0))->shape();
1495 : IPosition blc(4, 0, 0, 0, 0);
1496 : IPosition trc=shape-1;
1497 : for(uInt tix=0; tix<2 * unFacettedImStore_p->getNTaylorTerms() - 1; tix++)
1498 : {
1499 : TempImage<Float> onepsf((itsMappers.imageStore(0)->psf(tix))->shape(),
1500 : (itsMappers.imageStore(0)->psf(tix))->coordinates());
1501 : onepsf.copyData(*(itsMappers.imageStore(0)->psf(tix)));
1502 : //now set the original to 0 as we have a copy of one facet psf
1503 : (unFacettedImStore_p->psf(tix))->set(0.0);
1504 : blc[0]=(shape[0]-(onepsf.shape()[0]))/2;
1505 : trc[0]=onepsf.shape()[0]+blc[0]-1;
1506 : blc[1]=(shape[1]-(onepsf.shape()[1]))/2;
1507 : trc[1]=onepsf.shape()[1]+blc[1]-1;
1508 : Slicer sl(blc, trc, Slicer::endIsLast);
1509 : SubImage<Float> sub(*(unFacettedImStore_p->psf(tix)), sl, true);
1510 : sub.copyData(onepsf);
1511 : }
1512 : }
1513 :
1514 : //cout << "In setPsfFromOneFacet : sumwt : " << unFacettedImStore_p->sumwt()->get() << endl;
1515 :
1516 : }
1517 : */
1518 :
1519 0 : Block<CountedPtr<SIImageStore> > SynthesisImager::createChanChunkImageStoreList(
1520 : CountedPtr<SIImageStore> imagestore,
1521 : Int chanchunks)
1522 : {
1523 0 : LogIO os(LogOrigin("SynthesisImager", "createChanChunkImageStoreList"));
1524 :
1525 0 : Int extrachunks=0;
1526 0 : Int chunksize=imagestore->getShape()[3]/chanchunks;
1527 0 : Int rem=imagestore->getShape()[3] % chanchunks;
1528 : ///Avoid an extra chunk with 1 channel as it cause bumps in linear interpolation
1529 : ///See CAS-12625
1530 0 : while((rem==1) && (chunksize >1)){
1531 0 : chanchunks +=1;
1532 0 : chunksize=imagestore->getShape()[3]/chanchunks;
1533 0 : rem=imagestore->getShape()[3] % chanchunks;
1534 : }
1535 :
1536 :
1537 :
1538 0 : if( rem>0 )
1539 : {
1540 : // os << LogIO::WARN << "chanchunks ["+String::toString(chanchunks)+"] is not a divisor of nchan ["+String::toString(imagestore->getShape()[3])+"].";
1541 : // os << "Therefore, "+String::toString(imagestore->getShape()[3] % chanchunks)+" channel(s) at the end of the cube will be treated as an extra chunk." << LogIO::POST ;
1542 :
1543 0 : if( rem < chunksize ) extrachunks=1;
1544 : else
1545 : {
1546 0 : extrachunks=rem/chunksize;
1547 0 : if( rem % chunksize > 0 ) extrachunks += 1;
1548 : }
1549 : }
1550 :
1551 :
1552 0 : os << "Creating " << chanchunks +extrachunks << " reference subCubes (channel chunks) for gridding " << LogIO::POST;
1553 :
1554 0 : Block<CountedPtr<SIImageStore> > chunkList( chanchunks + extrachunks );
1555 :
1556 0 : if( chanchunks==1 ) { chunkList[0] = imagestore; return chunkList; }
1557 :
1558 : // Remember, only the FIRST field in each run can have chanchunks. So, check for this.
1559 0 : if( ! unChanChunkedImStore_p.null() ) {
1560 0 : throw( AipsError("A channel chunked image has already been set. Chanchunks are supported only for the main (first) field. Please submit a feature-request if you need multiple chanchunks for outlier fields as well. ") );
1561 : }
1562 :
1563 0 : unChanChunkedImStore_p = imagestore;
1564 0 : chanChunksStore_p = chanchunks;
1565 :
1566 0 : for (Int chunk=0; chunk< chanchunks+extrachunks; ++chunk){
1567 0 : chunkList[chunk] = unChanChunkedImStore_p->getSubImageStore(0,1,chunk, chanchunks);
1568 : }
1569 :
1570 0 : return chunkList;
1571 0 : }
1572 :
1573 :
1574 :
1575 :
1576 0 : void SynthesisImager::appendToMapperList(String imagename,
1577 : CoordinateSystem& csys,
1578 : IPosition imshape,
1579 : CountedPtr<FTMachine>& ftm,
1580 : CountedPtr<FTMachine>& iftm,
1581 : Quantity distance,
1582 : Int facets,
1583 : Int chanchunks,
1584 : const Bool overwrite,
1585 : String mappertype,
1586 : Float padding,
1587 : uInt ntaylorterms,
1588 : const Vector<String> &startmodel)
1589 : {
1590 0 : LogIO log_l(LogOrigin("SynthesisImager", "appendToMapperList(ftm)"));
1591 : //---------------------------------------------
1592 : // Some checks..
1593 0 : if(facets > 1 && itsMappers.nMappers() > 0)
1594 0 : log_l << "Facetted image has to be the first of multifields" << LogIO::EXCEPTION;
1595 :
1596 0 : TcleanProcessingInfo procInfo;
1597 0 : if(chanchunks<1)
1598 : {
1599 0 : log_l << "Automatically calculate chanchunks";
1600 0 : log_l << " using imshape : " << imshape << LogIO::POST;
1601 :
1602 : // Do calculation here.
1603 : // This runs once per image field (for multi-field imaging)
1604 : // This runs once per cube partition, and will see only its own partition's shape
1605 0 : chanchunks=1;
1606 :
1607 0 : CompositeNumber cn(uInt(imshape[0] * 2));
1608 : // heuristic factors multiplied to imshape based on gridder
1609 0 : size_t fudge_factor = 15;
1610 0 : if (ftm->name()=="MosaicFTNew") {
1611 0 : fudge_factor = 20;
1612 : }
1613 0 : else if (ftm->name()=="GridFT") {
1614 0 : fudge_factor = 9;
1615 : }
1616 :
1617 0 : Double required_mem = fudge_factor * sizeof(Float);
1618 0 : for (size_t i = 0; i < imshape.nelements(); i++) {
1619 : // gridding pads image and increases to composite number
1620 0 : if (i < 2) {
1621 0 : required_mem *= cn.nextLargerEven(Int(padding*Float(imshape[i])-0.5));
1622 : }
1623 : else {
1624 0 : required_mem *= imshape[i];
1625 : }
1626 : }
1627 :
1628 : // get number of tclean processes running on the same machine
1629 0 : size_t nlocal_procs = 1;
1630 0 : if (getenv("OMPI_COMM_WORLD_LOCAL_SIZE")) {
1631 0 : std::stringstream ss(getenv("OMPI_COMM_WORLD_LOCAL_SIZE"));
1632 0 : ss >> nlocal_procs;
1633 0 : }
1634 : // assumes all processes need the same amount of memory
1635 0 : required_mem *= nlocal_procs;
1636 : Double usr_memfrac, usr_mem;
1637 0 : AipsrcValue<Double>::find(usr_memfrac, "system.resources.memfrac", 80.);
1638 0 : AipsrcValue<Double>::find(usr_mem, "system.resources.memory", -1024.);
1639 : Double memory_avail;
1640 0 : if (usr_mem > 0.) {
1641 0 : memory_avail = usr_mem * 1024. * 1024.;
1642 : }
1643 : else {
1644 0 : memory_avail = Double(HostInfo::memoryFree()) * (usr_memfrac / 100.) * 1024.;
1645 : }
1646 :
1647 : // compute required chanchunks to fit into the available memory
1648 0 : chanchunks = (int)std::ceil((Double)required_mem / memory_avail);
1649 0 : if (imshape.nelements() == 4 && imshape[3] < chanchunks) {
1650 0 : chanchunks = imshape[3];
1651 : // TODO make chanchunks a divisor of nchannels?
1652 : }
1653 0 : chanchunks = chanchunks < 1 ? 1 : chanchunks;
1654 :
1655 0 : log_l << "Required memory " << required_mem / nlocal_procs / 1024. / 1024. / 1024.
1656 0 : << "\nAvailable memory " << memory_avail / 1024. / 1024 / 1024.
1657 : << " (rc: memory fraction " << usr_memfrac << "% memory " << usr_mem / 1024.
1658 : << ")\n" << nlocal_procs << " other processes on node\n"
1659 0 : << "Setting chanchunks to " << chanchunks << LogIO::POST;
1660 :
1661 0 : procInfo.mpiprocs = nlocal_procs;
1662 0 : procInfo.chnchnks = chanchunks;
1663 0 : const float toGB = 1024.0 * 1024.0 * 1024.0;
1664 0 : procInfo.memavail = memory_avail / toGB;
1665 0 : procInfo.memreq = required_mem / toGB;
1666 0 : }
1667 :
1668 0 : if( imshape.nelements()==4 && imshape[3]<chanchunks )
1669 : {
1670 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;
1671 : }
1672 :
1673 0 : if(chanchunks > 1 && itsMappers.nMappers() > 0)
1674 0 : log_l << "Channel chunking is currently not supported with multi(outlier)-fields. Please submit a feature request if needed." << LogIO::EXCEPTION;
1675 :
1676 0 : if(chanchunks > 1) itsDataLoopPerMapper=true;
1677 :
1678 0 : AlwaysAssert( ( ( ! (ftm->name()=="MosaicFTNew" && mappertype=="imagemosaic") ) &&
1679 : ( ! (ftm->name()=="AWProjectWBFTNew" && mappertype=="imagemosaic") )) ,
1680 : AipsError );
1681 : //---------------------------------------------
1682 :
1683 : // Create the ImageStore object
1684 0 : CountedPtr<SIImageStore> imstor;
1685 0 : MSColumns msc(mss4vi_p[0]);
1686 0 : imstor = createIMStore(imagename, csys, imshape, overwrite, msc, mappertype, ntaylorterms, distance, procInfo, facets, iftm->useWeightImage(), startmodel );
1687 :
1688 : // Create the Mappers
1689 0 : if( facets<2 && chanchunks<2) // One facet. Just add the above imagestore to the mapper list.
1690 : {
1691 0 : itsMappers.addMapper( createSIMapper( mappertype, imstor, ftm, iftm, ntaylorterms) );
1692 : }
1693 : else // This field is facetted. Make a list of reference imstores, and add all to the mapper list.
1694 : {
1695 :
1696 0 : if ( facets>1 && chanchunks==1 )
1697 : {
1698 : // Make and connect the list.
1699 0 : Block<CountedPtr<SIImageStore> > imstorList = createFacetImageStoreList( imstor, facets );
1700 0 : for( uInt facet=0; facet<imstorList.nelements(); facet++)
1701 : {
1702 0 : CountedPtr<FTMachine> new_ftm, new_iftm;
1703 0 : if(facet==0){ new_ftm = ftm; new_iftm = iftm; }
1704 0 : else{ new_ftm=ftm->cloneFTM(); new_iftm=iftm->cloneFTM(); }
1705 0 : itsMappers.addMapper(createSIMapper( mappertype, imstorList[facet], new_ftm, new_iftm, ntaylorterms));
1706 0 : }
1707 0 : }// facets
1708 0 : else if ( facets==1 && chanchunks>1 )
1709 : {
1710 : // Make and connect the list.
1711 0 : Block<CountedPtr<SIImageStore> > imstorList = createChanChunkImageStoreList( imstor, chanchunks );
1712 0 : for( uInt chunk=0; chunk<imstorList.nelements(); chunk++)
1713 : {
1714 0 : CountedPtr<FTMachine> new_ftm, new_iftm;
1715 0 : if(chunk==0){ new_ftm = ftm; new_iftm = iftm; }
1716 0 : else{ new_ftm=ftm->cloneFTM(); new_iftm=iftm->cloneFTM(); }
1717 0 : itsMappers.addMapper(createSIMapper( mappertype, imstorList[chunk], new_ftm, new_iftm, ntaylorterms));
1718 0 : }
1719 0 : }// chanchunks
1720 : else
1721 : {
1722 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. ") );
1723 : }
1724 :
1725 : }// facets or chunks
1726 :
1727 0 : }
1728 :
1729 : /////////////////////////
1730 : /*
1731 : Bool SynthesisImager::toUseWeightImage(CountedPtr<FTMachine>& ftm, String mappertype)
1732 : {
1733 : if( (ftm->name() == "GridFT" || ftm->name() == "WProjectFT")&&(mappertype!="imagemosaic") )
1734 : { return false; }
1735 : else
1736 : { return true; }
1737 : }
1738 : */
1739 :
1740 : // Make the FT-Machine and related objects (cfcache, etc.)
1741 0 : void SynthesisImager::createFTMachine(CountedPtr<FTMachine>& theFT,
1742 : CountedPtr<FTMachine>& theIFT,
1743 : const String& ftname,
1744 : const uInt nTaylorTerms,
1745 : const String mType,
1746 : const Int facets, //=1
1747 : //------------------------------
1748 : const Int wprojplane, //=1,
1749 : const Float padding, //=1.0,
1750 : const Bool useAutocorr, //=false,
1751 : const Bool useDoublePrec, //=true,
1752 : const String gridFunction, //=String("SF"),
1753 : //------------------------------
1754 : const Bool aTermOn, //= true,
1755 : const Bool psTermOn, //= true,
1756 : const Bool mTermOn, //= false,
1757 : const Bool wbAWP, //= true,
1758 : const String cfCache, //= "",
1759 : const Bool usePointing, //= false,
1760 : const Bool doPBCorr, //= true,
1761 : const Bool conjBeams, //= true,
1762 : const Float computePAStep, //=360.0
1763 : const Float rotatePAStep, //=5.0
1764 : const String interpolation, //="linear"
1765 : const Bool freqFrameValid, //=true
1766 : const Int cache, //=1000000000,
1767 : const Int tile, //=16
1768 : const String stokes, //=I
1769 : const String imageNamePrefix
1770 : )
1771 :
1772 : {
1773 0 : LogIO os( LogOrigin("SynthesisImager","createFTMachine",WHERE));
1774 :
1775 0 : if(ftname=="gridft"){
1776 0 : if(facets >1){
1777 0 : theFT=new GridFT(cache, tile, gridFunction, mLocation_p, phaseCenter_p, padding, useAutocorr, useDoublePrec);
1778 0 : theIFT=new GridFT(cache, tile, gridFunction, mLocation_p, phaseCenter_p, padding, useAutocorr, useDoublePrec);
1779 :
1780 : }
1781 : else{
1782 0 : theFT=new GridFT(cache, tile, gridFunction, mLocation_p, padding, useAutocorr, useDoublePrec);
1783 0 : theIFT=new GridFT(cache, tile, gridFunction, mLocation_p, padding, useAutocorr, useDoublePrec);
1784 : }
1785 : }
1786 0 : else if(ftname== "wprojectft"){
1787 0 : Double maxW=-1.0;
1788 0 : Double minW=-1.0;
1789 0 : Double rmsW=-1.0;
1790 0 : if(wprojplane <1)
1791 0 : WProjectFT::wStat(*rvi_p, minW, maxW, rmsW);
1792 0 : if(facets >1){
1793 0 : theFT=new WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
1794 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1795 0 : theIFT=new WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
1796 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1797 : }
1798 : else{
1799 0 : theFT=new WProjectFT(wprojplane, mLocation_p,
1800 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1801 0 : theIFT=new WProjectFT(wprojplane, mLocation_p,
1802 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1803 : }
1804 0 : CountedPtr<WPConvFunc> sharedconvFunc=static_cast<WProjectFT &>(*theFT).getConvFunc();
1805 : //static_cast<WProjectFT &>(*theFT).setConvFunc(sharedconvFunc);
1806 0 : static_cast<WProjectFT &>(*theIFT).setConvFunc(sharedconvFunc);
1807 0 : }
1808 0 : else if ((ftname == "awprojectft") || (ftname== "mawprojectft") || (ftname == "protoft")) {
1809 0 : createAWPFTMachine(theFT, theIFT, ftname, facets, wprojplane,
1810 : padding, useAutocorr, useDoublePrec, gridFunction,
1811 : aTermOn, psTermOn, mTermOn, wbAWP, cfCache,
1812 : usePointing, doPBCorr, conjBeams, computePAStep,
1813 : rotatePAStep, cache,tile,imageNamePrefix);
1814 : }
1815 0 : else if ( ftname == "mosaic" || ftname== "mosft" || ftname == "mosaicft" || ftname== "MosaicFT"){
1816 :
1817 0 : createMosFTMachine(theFT, theIFT, padding, useAutocorr, useDoublePrec, rotatePAStep, stokes, conjBeams);
1818 : }
1819 : else
1820 : {
1821 0 : throw( AipsError( "Invalid FTMachine name : " + ftname ) );
1822 : }
1823 : /* else if(ftname== "MosaicFT"){
1824 :
1825 : }*/
1826 :
1827 :
1828 :
1829 : ///////// Now, clone and pack the chosen FT into a MultiTermFT if needed.
1830 0 : if( mType=="multiterm" )
1831 : {
1832 0 : AlwaysAssert( nTaylorTerms>=1 , AipsError );
1833 :
1834 0 : CountedPtr<FTMachine> theMTFT = new MultiTermFTNew( theFT , nTaylorTerms, true/*forward*/ );
1835 0 : CountedPtr<FTMachine> theMTIFT = new MultiTermFTNew( theIFT , nTaylorTerms, false/*forward*/ );
1836 :
1837 0 : theFT = theMTFT;
1838 0 : theIFT = theMTIFT;
1839 0 : }
1840 :
1841 :
1842 :
1843 :
1844 : ////// Now, set the SkyJones if needed, and if not internally generated.
1845 0 : if( mType=="imagemosaic" &&
1846 0 : (ftname != "awprojectft" && ftname != "mawprojectft" && ftname != "proroft") )
1847 : {
1848 0 : CountedPtr<SkyJones> vp;
1849 0 : MSColumns msc(mss4vi_p[0]);
1850 0 : Quantity parang(0.0,"deg");
1851 0 : Quantity skyposthreshold(0.0,"deg");
1852 0 : vp = new VPSkyJones(msc, true, parang, BeamSquint::NONE,skyposthreshold);
1853 :
1854 0 : Vector<CountedPtr<SkyJones> > skyJonesList(1);
1855 0 : skyJonesList(0) = vp;
1856 0 : theFT->setSkyJones( skyJonesList );
1857 0 : theIFT->setSkyJones( skyJonesList );
1858 :
1859 0 : }
1860 :
1861 : //// For mode=cubedata, set the freq frame to invalid..
1862 : // get this info from buildCoordSystem
1863 0 : Vector<Int> tspws(0);
1864 : //theFT->setSpw( tspws, false );
1865 : //theIFT->setSpw( tspws, false );
1866 0 : theFT->setSpw( tspws, freqFrameValid );
1867 0 : theIFT->setSpw( tspws, freqFrameValid );
1868 :
1869 : //// Set interpolation mode
1870 0 : theFT->setFreqInterpolation( interpolation );
1871 0 : theIFT->setFreqInterpolation( interpolation );
1872 :
1873 : //channel selections from spw param
1874 0 : theFT->setSpwChanSelection(chanSel_p);
1875 0 : theIFT->setSpwChanSelection(chanSel_p);
1876 0 : }
1877 :
1878 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1879 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1880 :
1881 0 : void SynthesisImager::createAWPFTMachine(CountedPtr<FTMachine>& theFT, CountedPtr<FTMachine>& theIFT,
1882 : const String&,// ftmName,
1883 : const Int,// facets, //=1
1884 : //------------------------------
1885 : const Int wprojPlane, //=1,
1886 : const Float,// padding, //=1.0,
1887 : const Bool,// useAutocorr, //=false,
1888 : const Bool useDoublePrec, //=true,
1889 : const String,// gridFunction, //=String("SF"),
1890 : //------------------------------
1891 : const Bool aTermOn, //= true,
1892 : const Bool psTermOn, //= true,
1893 : const Bool mTermOn, //= false,
1894 : const Bool wbAWP, //= true,
1895 : const String cfCache, //= "",
1896 : const Bool usePointing, //= false,
1897 : const Bool doPBCorr, //= true,
1898 : const Bool conjBeams, //= true,
1899 : const Float computePAStep, //=360.0
1900 : const Float rotatePAStep, //=5.0
1901 : const Int cache, //=1000000000,
1902 : const Int tile, //=16
1903 : const String imageNamePrefix
1904 : )
1905 :
1906 : {
1907 0 : LogIO os( LogOrigin("SynthesisImager","createAWPFTMachine",WHERE));
1908 :
1909 0 : if (wprojPlane<=1)
1910 : {
1911 : os << LogIO::NORMAL
1912 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
1913 0 : << LogIO::POST;
1914 0 : os << LogIO::NORMAL << "Performing WBA-Projection" << LogIO::POST; // Loglevel PROGRESS
1915 : }
1916 : // if((wprojPlane>1)&&(wprojPlane<64))
1917 : // {
1918 : // os << LogIO::WARN
1919 : // << "No. of w-planes set too low for W projection - recommend at least 128"
1920 : // << LogIO::POST;
1921 : // os << LogIO::NORMAL << "Performing WBAW-Projection" << LogIO::POST; // Loglevel PROGRESS
1922 : // }
1923 :
1924 : // CountedPtr<ATerm> apertureFunction = createTelescopeATerm(mss4vi_p[0], aTermOn);
1925 : // CountedPtr<PSTerm> psTerm = new PSTerm();
1926 : // CountedPtr<WTerm> wTerm = new WTerm();
1927 :
1928 : // //
1929 : // // Selectively switch off CFTerms.
1930 : // //
1931 : // if (aTermOn == false) {apertureFunction->setOpCode(CFTerms::NOOP);}
1932 : // if (psTermOn == false) psTerm->setOpCode(CFTerms::NOOP);
1933 :
1934 : // //
1935 : // // Construct the CF object with appropriate CFTerms.
1936 : // //
1937 : // CountedPtr<ConvolutionFunction> tt;
1938 : // tt = AWProjectFT::makeCFObject(aTermOn, psTermOn, true, mTermOn, wbAWP);
1939 : // CountedPtr<ConvolutionFunction> awConvFunc;
1940 : // // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm, !wbAWP);
1941 : // if ((ftmName=="mawprojectft") || (mTermOn))
1942 : // awConvFunc = new AWConvFuncEPJones(apertureFunction,psTerm,wTerm,wbAWP);
1943 : // else
1944 : // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm,wbAWP);
1945 :
1946 0 : MSObservationColumns msoc(mss4vi_p[0].observation());
1947 0 : String telescopeName=msoc.telescopeName()(0);
1948 : CountedPtr<ConvolutionFunction> awConvFunc = AWProjectFT::makeCFObject(telescopeName,
1949 : aTermOn,
1950 : psTermOn, (wprojPlane > 1),
1951 0 : mTermOn, wbAWP, conjBeams);
1952 : //
1953 : // Construct the appropriate re-sampler.
1954 : //
1955 0 : CountedPtr<VisibilityResamplerBase> visResampler;
1956 : // if (ftmName=="protoft") visResampler = new ProtoVR();
1957 : //elsef
1958 0 : visResampler = new AWVisResampler();
1959 : // CountedPtr<VisibilityResamplerBase> visResampler = new VisibilityResampler();
1960 :
1961 : //
1962 : // Construct and initialize the CF cache object.
1963 : //
1964 :
1965 :
1966 : // CountedPtr<CFCache> cfCacheObj = new CFCache();
1967 : // cfCacheObj->setCacheDir(cfCache.data());
1968 : // // cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
1969 : // cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
1970 : // cfCacheObj->initCache2();
1971 :
1972 0 : CountedPtr<CFCache> cfCacheObj;
1973 :
1974 :
1975 : //
1976 : // Finally construct the FTMachine with the CFCache, ConvFunc and
1977 : // Re-sampler objects.
1978 : //
1979 0 : Float pbLimit_l=1e-2;
1980 0 : theFT = new AWProjectWBFTNew(wprojPlane, cache/2,
1981 : cfCacheObj, awConvFunc,
1982 : visResampler,
1983 : /*true */usePointing, doPBCorr,
1984 : tile, computePAStep, pbLimit_l, true,conjBeams,
1985 0 : useDoublePrec);
1986 :
1987 0 : cfCacheObj = new CFCache();
1988 0 : cfCacheObj->setCacheDir(cfCache.data());
1989 : // Get the LAZYFILL setting from the user configuration. If not
1990 : // found, default to False.
1991 : //
1992 : // With lazy fill ON, CFCache loads the required CFs on-demand
1993 : // from the disk. And periodically triggers garbage collection to
1994 : // release CFs that aren't required immediately.
1995 0 : cfCacheObj->setLazyFill(SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1);
1996 :
1997 : // cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
1998 0 : cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
1999 0 : cfCacheObj->initCache2();
2000 :
2001 0 : theFT->setCFCache(cfCacheObj);
2002 :
2003 :
2004 0 : Quantity rotateOTF(rotatePAStep,"deg");
2005 0 : static_cast<AWProjectWBFTNew &>(*theFT).setObservatoryLocation(mLocation_p);
2006 0 : static_cast<AWProjectWBFTNew &>(*theFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
2007 :
2008 : // theIFT = new AWProjectWBFT(wprojPlane, cache/2,
2009 : // cfCacheObj, awConvFunc,
2010 : // visResampler,
2011 : // /*true */usePointing, doPBCorr,
2012 : // tile, computePAStep, pbLimit_l, true,conjBeams,
2013 : // useDoublePrec);
2014 :
2015 : // static_cast<AWProjectWBFT &>(*theIFT).setObservatoryLocation(mLocation_p);
2016 : // static_cast<AWProjectWBFT &>(*theIFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
2017 :
2018 0 : theIFT = new AWProjectWBFTNew(static_cast<AWProjectWBFTNew &>(*theFT));
2019 :
2020 : //// Send in Freq info.
2021 0 : os << "Sending frequency selection information " << mssFreqSel_p << " to AWP FTM." << LogIO::POST;
2022 0 : theFT->setSpwFreqSelection( mssFreqSel_p );
2023 0 : theIFT->setSpwFreqSelection( mssFreqSel_p );
2024 :
2025 0 : }
2026 :
2027 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2028 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2029 :
2030 0 : ATerm* SynthesisImager::createTelescopeATerm(const MeasurementSet& ms, const Bool& isATermOn)
2031 : {
2032 0 : LogIO os(LogOrigin("SynthesisImager", "createTelescopeATerm",WHERE));
2033 :
2034 0 : if (!isATermOn) return new NoOpATerm();
2035 :
2036 0 : MSObservationColumns msoc(ms.observation());
2037 0 : String ObsName=msoc.telescopeName()(0);
2038 0 : if ((ObsName == "EVLA") || (ObsName == "VLA"))
2039 0 : return new EVLAAperture();
2040 : else
2041 : {
2042 0 : os << "Telescope name ('"+
2043 0 : ObsName+"') in the MS not recognized to create the telescope specific ATerm"
2044 0 : << LogIO::WARN;
2045 : }
2046 :
2047 0 : return NULL;
2048 0 : }
2049 :
2050 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2051 :
2052 0 : void SynthesisImager::getVPRecord(Record &rec, PBMath::CommonPB &kpb, String telescop)
2053 : {
2054 0 : LogIO os(LogOrigin("SynthesisImager", "getVPRecord",WHERE));
2055 :
2056 0 : VPManager *vpman=VPManager::Instance();
2057 0 : if( itsVpTable != String("") )
2058 : {
2059 0 : os << "Loading Voltage Pattern information from " << itsVpTable << LogIO::POST;
2060 0 : vpman->loadfromtable(itsVpTable);
2061 : }
2062 : else
2063 : {
2064 : // os << "Using Voltage Pattern currently set in the VPManager" << LogIO::POST;
2065 0 : os << "Using default Voltage Patterns from the VPManager" << LogIO::POST;
2066 0 : vpman->reset();
2067 : }
2068 :
2069 :
2070 :
2071 :
2072 : // PBMath::CommonPB kpb;
2073 0 : PBMath::enumerateCommonPB(telescop, kpb);
2074 : // Record rec;
2075 0 : vpman->getvp(rec, telescop);
2076 :
2077 : //ostringstream oss; rec.print(oss);
2078 : //os << "Using VP model : " << oss.str() << LogIO::POST;
2079 :
2080 0 : if(rec.empty()){
2081 0 : if((telescop=="EVLA")){
2082 0 : os << LogIO::WARN << "vpmanager does not list EVLA. Using VLA beam parameters. To use EVLA beams, please use the vpmanager and set the vptable parameter in tclean (see inline help for vptable)." << LogIO::POST;
2083 0 : telescop="VLA";
2084 0 : vpman->getvp(rec, telescop);
2085 0 : kpb=PBMath::VLA;
2086 : }
2087 : else{
2088 : //os << LogIO::WARN << "vpmanager does not have a beam for antenna : "+telescop <<".\n Please use the vpanager to define one (and optionally save its state as a table and supply its name via the vptable parameter.)" << LogIO::POST;
2089 0 : os << LogIO::WARN << "vpmanager does not have a beam for antenna : "+telescop <<".\n If needed, please use the vpanager to define one (and optionally save its state as a table and supply its name via the vptable parameter.). \n For now, we will proceed by reading dish diameters from the ANTENNA subtable, and form Airy disk beams." << LogIO::POST;
2090 0 : kpb=PBMath::UNKNOWN;
2091 0 : rec.define("name","COMMONPB");
2092 0 : rec.define("commonpb","NONE");
2093 : }
2094 : }
2095 :
2096 0 : }// get VPRecord
2097 :
2098 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2099 0 : void SynthesisImager:: createMosFTMachine(CountedPtr<FTMachine>& theFT,CountedPtr<FTMachine>& theIFT, const Float /*padding*/, const Bool useAutoCorr, const Bool useDoublePrec, const Float rotatePAStep, const String stokes, const Bool /*doConjBeam*/){
2100 :
2101 0 : LogIO os(LogOrigin("SynthesisImager", "createMosFTMachine",WHERE));
2102 :
2103 0 : rvi_p->originChunks();
2104 0 : MSColumns msc(rvi_p->ms());
2105 0 : String telescop=msc.observation().telescopeName()(0);
2106 : ///Multi ms with different telescop
2107 0 : Bool multiTel=False;
2108 0 : for(rvi_p->originChunks(); rvi_p->moreChunks(); rvi_p->nextChunk()){
2109 0 : if(rvi_p->newMS() && telescop != MSColumns(rvi_p->ms()).observation().telescopeName()(0))
2110 0 : multiTel=True;
2111 : }
2112 0 : rvi_p->originChunks();
2113 :
2114 : PBMath::CommonPB kpb;
2115 0 : Record rec;
2116 0 : getVPRecord( rec, kpb, telescop );
2117 :
2118 : /*
2119 : VPManager *vpman=VPManager::Instance();
2120 :
2121 : if( itsVpTable != String("") ) vpman->loadfromtable(itsVpTable);
2122 :
2123 : PBMath::CommonPB kpb;
2124 : PBMath::enumerateCommonPB(telescop, kpb);
2125 : Record rec;
2126 : vpman->getvp(rec, telescop);
2127 : if(rec.empty()){
2128 : if((telescop=="EVLA")){
2129 : os << LogIO::WARN << "vpmanager does not list EVLA. Using VLA beam parameters" << LogIO::POST;
2130 : telescop="VLA";
2131 : vpman->getvp(rec, telescop);
2132 : kpb=PBMath::VLA;
2133 : }
2134 : else{
2135 : os << LogIO::SEVERE << "vpmanager does not have a beam "+telescop <<"\n You can use VPManager to define one" << LogIO::POST;
2136 : }
2137 : }
2138 :
2139 : */
2140 :
2141 0 : if(rec.empty()){os << LogIO::SEVERE << "Cannot proceed with mosaicft gridder without a valid PB model" << LogIO::POST; }
2142 :
2143 :
2144 0 : VPSkyJones* vps=NULL;
2145 0 : if(rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2146 0 : vps= new VPSkyJones(msc, true, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
2147 : /////Don't know which parameter has pb threshold cutoff that the user want
2148 : ////leaving at default
2149 : ////vps.setThreshold(minPB);
2150 :
2151 : }
2152 : else{
2153 0 : PBMath myPB(rec);
2154 0 : String whichPBMath;
2155 0 : PBMathInterface::namePBClass(myPB.whichPBClass(), whichPBMath);
2156 0 : os << "Using the PB defined by " << whichPBMath << " for beam calculation for telescope " << telescop << LogIO::POST;
2157 0 : vps= new VPSkyJones(telescop, myPB, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
2158 0 : kpb=PBMath::DEFAULT;
2159 0 : }
2160 :
2161 :
2162 : //cerr << "SImager tel " << ((vps) ? vps->telescope(): "NOTEL " ) << " multiTel " << multiTel << endl;
2163 :
2164 0 : theFT = new MosaicFTNew(vps, mLocation_p, stokes, 1000000000, 16, useAutoCorr,
2165 0 : useDoublePrec);
2166 0 : PBMathInterface::PBClass pbtype=((kpb==PBMath::EVLA) || multiTel)? PBMathInterface::COMMONPB: PBMathInterface::AIRY;
2167 0 : if(rec.asString("name")=="IMAGE")
2168 0 : pbtype=PBMathInterface::IMAGE;
2169 : ///Use Heterogenous array mode for the following
2170 : ///Added EVLA in it to use different beam models for different frequencies
2171 0 : if((kpb == PBMath::UNKNOWN) || (kpb==PBMath::OVRO) || (kpb==PBMath::ACA)
2172 0 : || (kpb==PBMath::ALMA) || (kpb==PBMath::EVLA) || multiTel){
2173 0 : CountedPtr<SimplePBConvFunc> mospb=new HetArrayConvFunc(pbtype, "");
2174 0 : static_cast<MosaicFTNew &>(*theFT).setConvFunc(mospb);
2175 0 : }
2176 : ///////////////////make sure both FTMachine share the same conv functions.
2177 0 : theIFT= new MosaicFTNew(static_cast<MosaicFTNew &>(*theFT));
2178 :
2179 :
2180 0 : }
2181 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2182 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2183 :
2184 : // Do MS-Selection and set up vi/vb.
2185 : // Only this functions needs to know anything about the MS
2186 0 : void SynthesisImager::createVisSet(const Bool writeAccess)
2187 : {
2188 0 : LogIO os( LogOrigin("SynthesisImager","createVisSet",WHERE) );
2189 :
2190 : ////////////Temporary revert to vi/vb
2191 0 : Block<Int> sort(0);
2192 0 : Block<MeasurementSet> msblock(mss4vi_p.nelements());
2193 : //for (uInt k=0; k< msblock.nelements(); ++k){
2194 : // msblock[k]=*mss_p[k];
2195 : //}
2196 :
2197 : //vs_p= new VisSet(blockMSSel_p, sort, noChanSel, useModelCol_p);
2198 0 : if(!writeAccess){
2199 :
2200 0 : if(rvi_p) delete rvi_p;
2201 0 : rvi_p = NULL;
2202 0 : rvi_p=new ROVisibilityIterator(mss4vi_p, sort);
2203 :
2204 : }
2205 : else{
2206 : // if(wvi_p) delete wvi_p;
2207 : // wvi_p = NULL;
2208 0 : wvi_p=new VisibilityIterator(mss4vi_p, sort);
2209 0 : rvi_p=wvi_p;
2210 : }
2211 0 : Block<Vector<Int> > blockGroup(msblock.nelements());
2212 0 : for (uInt k=0; k < msblock.nelements(); ++k){
2213 0 : blockGroup[k].resize(blockSpw_p[k].nelements());
2214 0 : blockGroup[k].set(1);
2215 : // cerr << "start " << blockStart_p[k] << " nchan " << blockNChan_p[k] << " step " << blockStep_p[k] << " spw "<< blockSpw_p[k] <<endl;
2216 : }
2217 :
2218 0 : rvi_p->selectChannel(blockGroup, blockStart_p, blockNChan_p,
2219 0 : blockStep_p, blockSpw_p);
2220 0 : rvi_p->useImagingWeight(VisImagingWeight("natural"));
2221 : ////////////////////end of revert vi/vb
2222 0 : }// end of createVisSet
2223 :
2224 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2225 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2226 :
2227 :
2228 0 : void SynthesisImager::runMajorCycle(const Bool dopsf,
2229 : const Bool savemodel)
2230 : {
2231 0 : LogIO os( LogOrigin("SynthesisImager","runMajorCycle",WHERE) );
2232 :
2233 : // cout << "Savemodel : " << savemodel << " readonly : " << readOnly_p << " usescratch : " << useScratch_p << endl;
2234 :
2235 0 : Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
2236 0 : Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;
2237 :
2238 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
2239 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
2240 :
2241 0 : SynthesisUtilMethods::getResource("Start Major Cycle");
2242 :
2243 0 : itsMappers.checkOverlappingModels("blank");
2244 :
2245 : {
2246 0 : VisBufferAutoPtr vb(rvi_p);
2247 0 : rvi_p->originChunks();
2248 0 : rvi_p->origin();
2249 :
2250 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2251 0 : dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
2252 0 : Int cohDone=0;
2253 :
2254 :
2255 0 : if(!dopsf)itsMappers.initializeDegrid(*vb);
2256 0 : itsMappers.initializeGrid(*vb,dopsf);
2257 : //SynthesisUtilMethods::getResource("After initGrid for all mappers");
2258 :
2259 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2260 : {
2261 :
2262 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2263 : {
2264 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
2265 : // cerr << "nRows "<< vb->nRow() << " " << max(vb->visCube()) << endl;
2266 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2267 : {
2268 0 : if(!dopsf) {
2269 0 : { vb->setModelVisCube(Complex(0.0, 0.0)); }
2270 0 : itsMappers.degrid(*vb, savevirtualmodel );
2271 0 : if(savemodelcolumn && writeAccess_p )
2272 : {
2273 : // cout << "Vi: " << vb->modelVisCube() << endl;
2274 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2275 : }
2276 : }
2277 0 : itsMappers.grid(*vb, dopsf, datacol_p);
2278 0 : cohDone += vb->nRow();
2279 0 : pm.update(Double(cohDone));
2280 : }
2281 : }
2282 : }
2283 : //cerr << "IN SYNTHE_IMA" << endl;
2284 : //VisModelData::listModel(rvi_p->getMeasurementSet());
2285 0 : SynthesisUtilMethods::getResource("Before finalize for all mappers");
2286 0 : if(!dopsf) itsMappers.finalizeDegrid(*vb);
2287 0 : itsMappers.finalizeGrid(*vb, dopsf);
2288 :
2289 0 : }
2290 :
2291 0 : itsMappers.checkOverlappingModels("restore");
2292 :
2293 0 : unlockMSs();
2294 :
2295 0 : SynthesisUtilMethods::getResource("End Major Cycle");
2296 :
2297 0 : }// end runMajorCycle
2298 :
2299 :
2300 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2301 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2302 :
2303 : /// The mapper loop is outside the data iterator loop.
2304 : /// This is for cases where the image size is large compared to the RAM and
2305 : /// where data I/O is the relatively minor cost.
2306 0 : void SynthesisImager::runMajorCycle2(const Bool dopsf,
2307 : const Bool savemodel)
2308 : {
2309 0 : LogIO os( LogOrigin("SynthesisImager","runMajorCycle2",WHERE) );
2310 :
2311 : // cout << "Savemodel : " << savemodel << " readonly : " << readOnly_p << " usescratch : " << useScratch_p << endl;
2312 :
2313 0 : Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
2314 0 : Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;
2315 :
2316 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
2317 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
2318 :
2319 0 : itsMappers.checkOverlappingModels("blank");
2320 :
2321 0 : Bool resetModel=False;
2322 0 : if( savemodelcolumn && writeAccess_p)
2323 : {
2324 0 : resetModel=True;
2325 0 : os << "Iterating through the model column to reset it to zero" << LogIO::POST;
2326 0 : VisBufferAutoPtr vb(rvi_p);
2327 0 : rvi_p->originChunks();
2328 0 : rvi_p->origin();
2329 :
2330 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2331 0 : dopsf?"Seting model column to zero":"pre-Major Cycle", "","","",True);
2332 0 : Int cohDone=0;
2333 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2334 : {
2335 :
2336 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2337 : {
2338 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2339 : {
2340 0 : vb->setModelVisCube(Complex(0.0, 0.0));
2341 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2342 : }
2343 0 : cohDone += vb->nRow();
2344 0 : pm.update(Double(cohDone));
2345 : }
2346 : }
2347 0 : }// setting model to zero
2348 :
2349 0 : for(Int gmap=0;gmap<itsMappers.nMappers();gmap++)
2350 : {
2351 0 : os << "Running major cycle for chunk : " << gmap << LogIO::POST;
2352 :
2353 0 : SynthesisUtilMethods::getResource("Start Major Cycle for mapper"+String::toString(gmap));
2354 0 : VisBufferAutoPtr vb(rvi_p);
2355 0 : rvi_p->originChunks();
2356 0 : rvi_p->origin();
2357 :
2358 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2359 0 : dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
2360 0 : Int cohDone=0;
2361 :
2362 0 : if(!dopsf){
2363 0 : itsMappers.initializeDegrid(*vb, gmap);
2364 : //itsMappers.getMapper(gmap)->initializeDegrid(*vb);
2365 : }
2366 0 : itsMappers.initializeGrid(*vb,dopsf, gmap);
2367 : //itsMappers.getMapper(gmap)->initializeGrid(*vb,dopsf);
2368 :
2369 0 : SynthesisUtilMethods::getResource("After initialize for mapper"+String::toString(gmap));
2370 :
2371 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2372 : {
2373 :
2374 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2375 : {
2376 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
2377 : // cerr << "nRows "<< vb->nRow() << " " << max(vb->visCube()) << endl;
2378 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2379 : {
2380 0 : if(!dopsf) {
2381 0 : if(resetModel==False) { vb->setModelVisCube(Complex(0.0, 0.0)); }
2382 0 : itsMappers.degrid(*vb, savevirtualmodel, gmap );
2383 : //itsMappers.getMapper(gmap)->degrid(*vb); //, savevirtualmodel );
2384 0 : if(savemodelcolumn && writeAccess_p )
2385 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2386 : }
2387 0 : itsMappers.grid(*vb, dopsf, datacol_p, gmap);
2388 : //itsMappers.getMapper(gmap)->grid(*vb, dopsf, datacol_p);
2389 0 : cohDone += vb->nRow();
2390 0 : pm.update(Double(cohDone));
2391 : }
2392 : }
2393 : }
2394 : //cerr << "IN SYNTHE_IMA" << endl;
2395 : //VisModelData::listModel(rvi_p->getMeasurementSet());
2396 :
2397 0 : SynthesisUtilMethods::getResource("Before finalize for mapper"+String::toString(gmap));
2398 :
2399 0 : if(!dopsf)
2400 : {
2401 0 : itsMappers.finalizeDegrid(*vb,gmap);
2402 : // auto back = itsMappers.getMapper(gmap)->imageStore()->forwardGrid();
2403 : // back->tempClose();
2404 : }
2405 0 : itsMappers.finalizeGrid(*vb, dopsf,gmap);
2406 :
2407 : //itsMappers.getMapper(gmap)->imageStore()->releaseLocks();
2408 0 : itsMappers.getMapper(gmap)->imageStore()->releaseComplexGrids();
2409 :
2410 :
2411 : // auto back = itsMappers.getMapper(gmap)->imageStore()->backwardGrid();
2412 : // back->tempClose();
2413 :
2414 0 : SynthesisUtilMethods::getResource("End Major Cycle for mapper"+String::toString(gmap));
2415 0 : }// end of mapper loop
2416 :
2417 0 : itsMappers.checkOverlappingModels("restore");
2418 :
2419 0 : unlockMSs();
2420 :
2421 0 : SynthesisUtilMethods::getResource("End Major Cycle");
2422 :
2423 0 : }// end runMajorCycle2
2424 :
2425 :
2426 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2427 0 : void SynthesisImager::predictModel(){
2428 0 : LogIO os( LogOrigin("SynthesisImager","predictModel ",WHERE) );
2429 :
2430 0 : os << "---------------------------------------------------- Predict Model ---------------------------------------------" << LogIO::POST;
2431 :
2432 0 : Bool savemodelcolumn = !readOnly_p && useScratch_p;
2433 0 : Bool savevirtualmodel = !readOnly_p && !useScratch_p;
2434 :
2435 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
2436 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
2437 :
2438 0 : itsMappers.checkOverlappingModels("blank");
2439 :
2440 :
2441 : {
2442 0 : VisBufferAutoPtr vb(rvi_p);
2443 0 : rvi_p->originChunks();
2444 0 : rvi_p->origin();
2445 :
2446 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2447 0 : "Predict Model", "","","",true);
2448 0 : Int cohDone=0;
2449 :
2450 0 : itsMappers.initializeDegrid(*vb);
2451 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2452 : {
2453 :
2454 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2455 : {
2456 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; //No valid rows in this MS
2457 : //if !usescratch ...just save
2458 0 : vb->setModelVisCube(Complex(0.0, 0.0));
2459 0 : itsMappers.degrid(*vb, savevirtualmodel);
2460 0 : if(savemodelcolumn && writeAccess_p )
2461 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2462 :
2463 : // cout << "nRows "<< vb->nRow() << " " << max(vb->modelVisCube()) << endl;
2464 0 : cohDone += vb->nRow();
2465 0 : pm.update(Double(cohDone));
2466 :
2467 : }
2468 : }
2469 0 : itsMappers.finalizeDegrid(*vb);
2470 0 : }
2471 :
2472 0 : itsMappers.checkOverlappingModels("restore");
2473 0 : unlockMSs();
2474 :
2475 0 : }// end of predictModel
2476 :
2477 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2478 0 : void SynthesisImager::makeSdImage(Bool dopsf)
2479 : {
2480 0 : LogIO os( LogOrigin("SynthesisImager","makeSdImage",WHERE) );
2481 :
2482 : // Bool dopsf=false;
2483 0 : if(datacol_p==FTMachine::PSF) dopsf=true;
2484 :
2485 : {
2486 0 : VisBufferAutoPtr vb(rvi_p);
2487 0 : rvi_p->originChunks();
2488 0 : rvi_p->origin();
2489 :
2490 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2491 0 : String(datacol_p), "","","",true);
2492 0 : Int cohDone=0;
2493 :
2494 0 : itsMappers.initializeGrid(*vb,dopsf);
2495 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2496 : {
2497 :
2498 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2499 : {
2500 0 : itsMappers.grid(*vb, dopsf, datacol_p);
2501 0 : cohDone += vb->nRow();
2502 0 : pm.update(Double(cohDone));
2503 : }
2504 : }
2505 0 : itsMappers.finalizeGrid(*vb, dopsf);
2506 :
2507 0 : }
2508 :
2509 0 : unlockMSs();
2510 :
2511 0 : }// end makeSdImage
2512 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2513 :
2514 0 : void SynthesisImager::makeImage(String type, const String& imagename, const String& complexImage, const Int whichModel)
2515 : {
2516 0 : LogIO os( LogOrigin("SynthesisImager","makeImage",WHERE) );
2517 :
2518 :
2519 0 : FTMachine::Type seType(FTMachine::OBSERVED);
2520 0 : if(type=="observed") {
2521 0 : seType=FTMachine::OBSERVED;
2522 : os << LogIO::NORMAL // Loglevel INFO
2523 : << "Making dirty image from " << type << " data "
2524 0 : << LogIO::POST;
2525 : }
2526 0 : else if (type=="model") {
2527 0 : seType=FTMachine::MODEL;
2528 : os << LogIO::NORMAL // Loglevel INFO
2529 : << "Making dirty image from " << type << " data "
2530 0 : << LogIO::POST;
2531 : }
2532 0 : else if (type=="corrected") {
2533 0 : seType=FTMachine::CORRECTED;
2534 : os << LogIO::NORMAL // Loglevel INFO
2535 : << "Making dirty image from " << type << " data "
2536 0 : << LogIO::POST;
2537 : }
2538 0 : else if (type=="psf") {
2539 0 : seType=FTMachine::PSF;
2540 : os << "Making point spread function "
2541 0 : << LogIO::POST;
2542 : }
2543 0 : else if (type=="residual") {
2544 0 : seType=FTMachine::RESIDUAL;
2545 : os << LogIO::NORMAL // Loglevel INFO
2546 : << "Making dirty image from " << type << " data "
2547 0 : << LogIO::POST;
2548 : }
2549 0 : else if (type=="singledish-observed") {
2550 0 : seType=FTMachine::OBSERVED;
2551 : os << LogIO::NORMAL
2552 0 : << "Making single dish image from observed data" << LogIO::POST;
2553 : }
2554 0 : else if (type=="singledish") {
2555 0 : seType=FTMachine::CORRECTED;
2556 : os << LogIO::NORMAL
2557 0 : << "Making single dish image from corrected data" << LogIO::POST;
2558 : }
2559 0 : else if (type=="coverage") {
2560 0 : seType=FTMachine::COVERAGE;
2561 : os << LogIO::NORMAL
2562 : << "Making single dish coverage function "
2563 0 : << LogIO::POST;
2564 : }
2565 0 : else if (type=="holography") {
2566 0 : seType=FTMachine::CORRECTED;
2567 : os << LogIO::NORMAL
2568 : << "Making complex holographic image from corrected data "
2569 0 : << LogIO::POST;
2570 : }
2571 0 : else if (type=="holography-observed") {
2572 0 : seType=FTMachine::OBSERVED;
2573 : os << LogIO::NORMAL
2574 : << "Making complex holographic image from observed data "
2575 0 : << LogIO::POST;
2576 : }
2577 :
2578 :
2579 0 : String imageName=(itsMappers.imageStore(whichModel))->getName();
2580 0 : String cImageName(complexImage);
2581 0 : Bool keepComplexImage=(complexImage!="")||(type.contains("holography"));
2582 0 : if(complexImage=="") {
2583 0 : cImageName=imageName+".compleximage";
2584 : }
2585 0 : PagedImage<Float> theImage((itsMappers.imageStore(whichModel))->getShape(), (itsMappers.imageStore(whichModel))->getCSys(), imagename);
2586 0 : PagedImage<Complex> cImageImage(theImage.shape(),
2587 : theImage.coordinates(),
2588 0 : cImageName);
2589 0 : if(!keepComplexImage)
2590 0 : cImageImage.table().markForDelete();
2591 :
2592 :
2593 0 : Matrix<Float> weight;
2594 0 : (itsMappers.getFTM(whichModel))->makeImage(seType, *rvi_p, cImageImage, weight);
2595 :
2596 0 : if(seType==FTMachine::PSF){
2597 0 : StokesImageUtil::ToStokesPSF(theImage, cImageImage);
2598 0 : StokesImageUtil::normalizePSF(theImage);
2599 : }
2600 : else{
2601 0 : StokesImageUtil::To(theImage, cImageImage);
2602 : }
2603 :
2604 :
2605 0 : unlockMSs();
2606 :
2607 0 : }// end makeImage
2608 :
2609 :
2610 :
2611 :
2612 :
2613 :
2614 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2615 : // Method to run the AWProjectFT machine to seperate the CFCache
2616 : // construction from imaging. This is done by splitting the
2617 : // operation in two steps: (1) run the FTM in "dry" mode to create a
2618 : // "blank" CFCache, and (2) re-load the "blank" CFCache and "fill"
2619 : // it.
2620 : //
2621 : // If someone can get me (SB) out of the horrible statc_casts in the
2622 : // code below, I will be most grateful (we are out of it! :-)).
2623 : //
2624 0 : void SynthesisImager::dryGridding(const Vector<String>& cfList)
2625 : {
2626 0 : LogIO os( LogOrigin("SynthesisImager","dryGridding",WHERE) );
2627 0 : Int cohDone=0, whichFTM=0;
2628 : (void)cfList;
2629 : // If not an AWProject-class FTM, make this call a NoOp. Might be
2630 : // useful to extend it to other projection FTMs -- but later.
2631 0 : String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();
2632 :
2633 0 : if (!((itsMappers.getFTM(whichFTM,true))->isUsingCFCache())) return;
2634 :
2635 0 : os << "---------------------------------------------------- Dry Gridding ---------------------------------------------" << LogIO::POST;
2636 :
2637 : //
2638 : // Go through the entire MS in "dry" mode to set up a "blank"
2639 : // CFCache. This is done by setting the AWPWBFT in dryrun mode
2640 : // and gridding. The process of gridding emits CFCache, which
2641 : // will be "blank" in a dry run.
2642 : {
2643 0 : VisBufferAutoPtr vb(rvi_p);
2644 0 : rvi_p->originChunks();
2645 0 : rvi_p->origin();
2646 :
2647 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()), "dryGridding", "","","",true);
2648 :
2649 0 : itsMappers.initializeGrid(*vb);
2650 :
2651 : // Set the gridder (iFTM) to run in dry-gridding mode
2652 0 : (itsMappers.getFTM(whichFTM,true))->setDryRun(true);
2653 :
2654 0 : Bool aTermIsOff=False;
2655 : {
2656 0 : CountedPtr<FTMachine> ftm=itsMappers.getFTM(whichFTM,True);
2657 0 : CountedPtr<ConvolutionFunction> cf=ftm->getAWConvFunc();
2658 : // Bool aTermIsOff = (((static_cast<AWConvFunc &> (*cf)).getTerm("ATerm")))->isNoOp();
2659 0 : aTermIsOff = cf->getTerm("ATerm")->isNoOp();
2660 0 : }
2661 :
2662 : os << "Making a \"blank\" CFCache"
2663 : << (aTermIsOff?" (without the A-Term)":"")
2664 0 : << LogIO::WARN << LogIO::POST;
2665 :
2666 : // Step through the MS. This triggers the logic in the Gridder
2667 : // to determine all the CFs that will be required. These empty
2668 : // CFs are written to the CFCache which can then be filled via
2669 : // a call to fillCFCache().
2670 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2671 : {
2672 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2673 : {
2674 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2675 : {
2676 0 : itsMappers.grid(*vb, true, FTMachine::OBSERVED, whichFTM);
2677 0 : cohDone += vb->nRow();
2678 0 : pm.update(Double(cohDone));
2679 : // If there is no term that depends on time, don't iterate over the entire data base
2680 0 : if (aTermIsOff) break;
2681 : }
2682 : }
2683 : }
2684 0 : }
2685 0 : if (cohDone == 0) os << "No valid rows found in dryGridding." << LogIO::EXCEPTION << LogIO::POST;
2686 : // Unset the dry-gridding mode.
2687 0 : (itsMappers.getFTM(whichFTM,true))->setDryRun(false);
2688 :
2689 : //itsMappers.checkOverlappingModels("restore");
2690 0 : unlockMSs();
2691 : //fillCFCache(cfList);
2692 0 : }
2693 : //
2694 : // Re-load the CFCache from the disk using the supplied list of CFs
2695 : // (as cfList). Then extract the ConvFunc object (which was setup
2696 : // in the FTM) and call it's makeConvFunction2() to fill the CF.
2697 : // Finally, unset the dry-run mode of the FTM.
2698 : //
2699 0 : void SynthesisImager::fillCFCache(const Vector<String>& cfList,
2700 : const String& ftmName,
2701 : const String& cfcPath,
2702 : const Bool& psTermOn,
2703 : const Bool& aTermOn,
2704 : const Bool& conjBeams)
2705 : {
2706 0 : LogIO os( LogOrigin("SynthesisImager","fillCFCache",WHERE) );
2707 : // If not an AWProject-class FTM, make this call a NoOp. Might be
2708 : // useful to extend it to other projection FTMs -- but later.
2709 : // String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();
2710 :
2711 0 : if (!ftmName.contains("awproject") and
2712 0 : !ftmName.contains("multitermftnew")) return;
2713 :
2714 0 : os << "---------------------------------------------------- fillCFCache ---------------------------------------------" << LogIO::POST;
2715 :
2716 : //String cfcPath = itsMappers.getFTM(whichFTM)->getCacheDir();
2717 : //String imageNamePrefix=itsMappers.getFTM(whichFTM)->getCFCache()->getWtImagePrefix();
2718 :
2719 : //cerr << "Path = " << path << endl;
2720 :
2721 : // CountedPtr<AWProjectWBFTNew> tmpFT = new AWProjectWBFTNew(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM))));
2722 :
2723 :
2724 0 : Float dPA=360.0,selectedPA=2*360.0;
2725 0 : if (cfList.nelements() > 0)
2726 : {
2727 0 : CountedPtr<CFCache> cfCacheObj = new CFCache();
2728 : //Vector<String> wtCFList; wtCFList.resize(cfList.nelements());
2729 : //for (Int i=0; i<wtCFList.nelements(); i++) wtCFList[i] = "WT"+cfList[i];
2730 : //Directory dir(path);
2731 0 : Vector<String> cfList_p=cfList;//dir.find(Regex(Regex::fromPattern("CFS*")));
2732 0 : Vector<String> wtCFList_p;
2733 0 : wtCFList_p.resize(cfList_p.nelements());
2734 0 : for (Int i=0; i<(Int)wtCFList_p.nelements(); i++) wtCFList_p[i]="WT"+cfList_p[i];
2735 :
2736 : //cerr << cfList_p << endl;
2737 0 : cfCacheObj->setCacheDir(cfcPath.data());
2738 :
2739 0 : os << "Re-loading the \"blank\" CFCache for filling" << LogIO::WARN << LogIO::POST;
2740 :
2741 0 : cfCacheObj->initCacheFromList2(cfcPath, cfList_p, wtCFList_p,
2742 : selectedPA, dPA,1);
2743 : // tmpFT->setCFCache(cfCacheObj);
2744 0 : Vector<Double> uvScale, uvOffset;
2745 0 : Matrix<Double> vbFreqSelection;
2746 0 : CountedPtr<CFStore2> cfs2 = CountedPtr<CFStore2>(&cfCacheObj->memCache2_p[0],false);//new CFStore2;
2747 0 : CountedPtr<CFStore2> cfwts2 = CountedPtr<CFStore2>(&cfCacheObj->memCacheWt2_p[0],false);//new CFStore2;
2748 :
2749 : //
2750 : // Get whichFTM from itsMappers (SIMapperCollection) and
2751 : // cast it as AWProjectWBFTNew. Then get the ConvFunc from
2752 : // the FTM and cast it as AWConvFunc. Finally call
2753 : // AWConvFunc::makeConvFunction2().
2754 : //
2755 : // (static_cast<AWConvFunc &>
2756 : // (*(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getAWConvFunc())
2757 : // ).makeConvFunction2(String(path), uvScale, uvOffset, vbFreqSelection,
2758 : // *cfs2, *cfwts2);
2759 :
2760 : // This is a global methond in AWConvFunc. Does not require
2761 : // FTM to be constructed (which can be expensive in terms of
2762 : // memory footprint).
2763 0 : AWConvFunc::makeConvFunction2(String(cfcPath), uvScale, uvOffset, vbFreqSelection,
2764 0 : *cfs2, *cfwts2, psTermOn, aTermOn, conjBeams);
2765 0 : }
2766 : //cerr << "Mem used = " << itsMappers.getFTM(whichFTM)->getCFCache()->memCache2_p[0].memUsage() << endl;
2767 : //(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getCFCache()->initCache2();
2768 0 : }
2769 :
2770 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2771 0 : void SynthesisImager::reloadCFCache()
2772 : {
2773 0 : LogIO os( LogOrigin("SynthesisImager","reloadCFCache",WHERE) );
2774 0 : Int whichFTM=0;
2775 0 : String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();
2776 0 : if (!ftmName.contains("AWProject")) return;
2777 :
2778 0 : os << "-------------------------------------------- reloadCFCache ---------------------------------------------" << LogIO::POST;
2779 0 : String path = itsMappers.getFTM(whichFTM)->getCacheDir();
2780 0 : String imageNamePrefix=itsMappers.getFTM(whichFTM)->getCFCache()->getWtImagePrefix();
2781 :
2782 0 : CountedPtr<CFCache> cfCacheObj = new CFCache();
2783 0 : cfCacheObj->setCacheDir(path.data());
2784 0 : cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
2785 0 : cfCacheObj->initCache2();
2786 :
2787 : // This assumes the itsMappers is always SIMapperCollection.
2788 0 : for (whichFTM = 0; whichFTM < itsMappers.nMappers(); whichFTM++)
2789 : {
2790 0 : (static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).setCFCache(cfCacheObj,true); // Setup iFTM
2791 0 : (static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM,false)))).setCFCache(cfCacheObj,true); // Set FTM
2792 : }
2793 0 : }
2794 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2795 :
2796 : //Utility function to properly convert Double to String.
2797 : //With C++11 we can probably use STL to_string() function instead...
2798 0 : String SynthesisImager::doubleToString(const Double& df) {
2799 0 : std::ostringstream ss;
2800 0 : ss.precision(std::numeric_limits<double>::digits10+2);
2801 0 : ss << df;
2802 : //return ss.str();
2803 0 : return to_string(df);
2804 0 : }
2805 :
2806 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2807 :
2808 0 : Bool SynthesisImager::makePB()
2809 : {
2810 0 : LogIO os( LogOrigin("SynthesisImager","makePB",WHERE) );
2811 :
2812 0 : if( itsMakeVP==False )
2813 : {
2814 0 : os << LogIO::NORMAL1 << "Not making .pb by direct evaluation. The gridder will make a .weight and a .pb will be computed from it." << LogIO::POST;
2815 : // Check that the .weight exists.. ?
2816 :
2817 0 : return False;
2818 : }
2819 : else
2820 : {
2821 0 : Bool doDefaultVP = itsVpTable.length()>0 ? False : True;
2822 :
2823 0 : CoordinateSystem coordsys=itsMappers.imageStore(0)->getCSys();
2824 0 : String telescope=coordsys.obsInfo().telescope();
2825 :
2826 0 : if (doDefaultVP) {
2827 :
2828 0 : MSAntennaColumns ac(mss4vi_p[0].antenna());
2829 0 : Double dishDiam=ac.dishDiameter()(0);
2830 0 : if(!allEQ(ac.dishDiameter().getColumn(), dishDiam))
2831 : os << LogIO::WARN
2832 : << "The MS has multiple antenna diameters ..PB could be wrong "
2833 0 : << LogIO::POST;
2834 0 : return makePBImage( telescope, False, dishDiam);
2835 0 : }
2836 : else{
2837 0 : return makePBImage(telescope );
2838 : }
2839 :
2840 0 : }
2841 : return False;
2842 0 : }
2843 :
2844 17 : bool SynthesisImager::makePBImage(const String& telescopeName,
2845 : Bool useSymmetricBeam, Double diam){
2846 :
2847 34 : LogIO os(LogOrigin("SynthesisImager", "makePBImage()", WHERE));
2848 :
2849 : // Check if this metadata info should come from each Mapper or if it's OK to be just the first one.
2850 : // Right now it assumes all mappers have the same FREQUENCY settings.
2851 :
2852 17 : CoordinateSystem imageCoord=itsMappers.imageStore(0)->getCSys();
2853 :
2854 17 : Int spectralIndex=imageCoord.findCoordinate(Coordinate::SPECTRAL);
2855 : SpectralCoordinate
2856 17 : spectralCoord=imageCoord.spectralCoordinate(spectralIndex);
2857 17 : Vector<String> units(1); units = "Hz";
2858 17 : spectralCoord.setWorldAxisUnits(units);
2859 17 : Vector<Double> spectralWorld(1);
2860 17 : Vector<Double> spectralPixel(1);
2861 17 : spectralPixel(0) = 0;
2862 17 : spectralCoord.toWorld(spectralWorld, spectralPixel);
2863 17 : Double freq = spectralWorld(0);
2864 17 : Quantity qFreq( freq, "Hz" );
2865 17 : String telName=telescopeName;
2866 17 : if(telName=="ALMA" && diam < 12.0)
2867 0 : telName="ACA";
2868 : //cerr << "Telescope Name is " << telName<< endl;
2869 : PBMath::CommonPB whichPB;
2870 17 : PBMath::enumerateCommonPB(telName, whichPB);
2871 17 : PBMath myPB;
2872 17 : if(whichPB!=PBMath::UNKNOWN && whichPB!=PBMath::NONE){
2873 :
2874 17 : myPB=PBMath(telName, useSymmetricBeam, qFreq);
2875 : }
2876 0 : else if(diam > 0.0){
2877 0 : myPB=PBMath(diam,useSymmetricBeam, qFreq);
2878 : }
2879 : else{
2880 : os << LogIO::WARN << "Telescope " << telName << " is not known\n "
2881 : << "Not making the PB image"
2882 0 : << LogIO::POST;
2883 0 : return false;
2884 : }
2885 17 : return makePrimaryBeam(myPB );
2886 17 : }
2887 :
2888 0 : Bool SynthesisImager::makePBImage(const String telescop){
2889 :
2890 : /*
2891 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
2892 : PBMath myPB(recCol(0));
2893 : */
2894 0 : LogIO os( LogOrigin("SynthesisImager","makePBImage",WHERE) );
2895 :
2896 : PBMath::CommonPB kpb;
2897 0 : Record rec;
2898 0 : getVPRecord( rec, kpb, telescop );
2899 :
2900 0 : Bool ret=False;
2901 0 : if(rec.empty())
2902 0 : {os << LogIO::WARN << "Not making PB. Cannot use pbcor option later." << LogIO::POST; }
2903 : else
2904 : {
2905 0 : PBMath myPB( rec );
2906 0 : ret = makePrimaryBeam(myPB);
2907 0 : }
2908 0 : return ret;
2909 0 : }
2910 :
2911 :
2912 0 : Bool SynthesisImager::makePrimaryBeam(PBMath& pbMath)
2913 : {
2914 0 : LogIO os( LogOrigin("SynthesisImager","makePrimaryBeam",WHERE) );
2915 :
2916 0 : os << "Evaluating Primary Beam model onto image grid(s)" << LogIO::POST;
2917 :
2918 0 : itsMappers.initPB();
2919 :
2920 0 : VisBuffer vb(*rvi_p);
2921 0 : Int fieldCounter=0;
2922 0 : Vector<Int> fieldsDone;
2923 :
2924 0 : for(rvi_p->originChunks(); rvi_p->moreChunks(); rvi_p->nextChunk()){
2925 0 : Bool fieldDone=false;
2926 0 : for (uInt k=0; k < fieldsDone.nelements(); ++k){
2927 0 : fieldDone=fieldDone || (vb.fieldId()==fieldsDone(k));
2928 : }
2929 0 : if(!fieldDone){
2930 0 : ++fieldCounter;
2931 0 : fieldsDone.resize(fieldCounter, true);
2932 0 : fieldsDone(fieldCounter-1)=vb.fieldId();
2933 :
2934 0 : itsMappers.addPB(vb,pbMath);
2935 :
2936 : }
2937 : }
2938 :
2939 0 : itsMappers.releaseImageLocks();
2940 0 : unlockMSs();
2941 :
2942 0 : return True;
2943 0 : }// end makePB
2944 :
2945 : /////===========
2946 0 : void SynthesisImager::setMovingSource(const String& movingSource){
2947 0 : movingSource_p=movingSource;
2948 0 : }
2949 :
2950 0 : bool SynthesisImager::isSpectralCube(){
2951 0 : bool retval=False;
2952 0 : for (Int k=0; k < itsMappers.nMappers(); ++k){
2953 : //For some reason imstore sometime returns 0 shape
2954 : //trying to test with with psf size breaks parallel = true in some cases...go figure
2955 : //if((((itsMappers.imageStore(k))->psf())->shape()[3]) != ((itsMappers.imageStore(k))->getShape()[3])){
2956 : //cerr << "shapes " << ((itsMappers.imageStore(k))->psf())->shape() << " " << ((itsMappers.imageStore(k))->getShape()) << endl;
2957 : //throw(AipsError("images shape seem insistent "));
2958 0 : if((itsMappers.imageStore(k))->getShape()(3) ==0)
2959 0 : return True;
2960 : //}
2961 0 : if((itsMappers.imageStore(k))->getShape()(3) > 1)
2962 0 : retval=True;
2963 : }
2964 0 : return retval;
2965 :
2966 : }
2967 17 : void SynthesisImager::normalizerinfo(const Record& normpars){
2968 17 : normpars_p=normpars;
2969 17 : }
2970 : } //# NAMESPACE CASA - END
2971 :
|
