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: aips2-request@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 0 : SynthesisImager::SynthesisImager() : itsMappers(SIMapperCollection()), writeAccess_p(True),
104 0 : gridpars_p(), impars_p(), movingSource_p(""), doingCubeGridding_p(True)
105 : {
106 :
107 0 : imwgt_p=VisImagingWeight("natural");
108 0 : imageDefined_p=false;
109 0 : useScratch_p=false;
110 0 : readOnly_p=true;
111 :
112 0 : mss4vi_p.resize(0);
113 0 : wvi_p=0;
114 0 : rvi_p=0;
115 :
116 0 : facetsStore_p=-1;
117 0 : chanChunksStore_p=-1;
118 0 : unFacettedImStore_p=NULL;
119 0 : unChanChunkedImStore_p=NULL;
120 0 : dataSel_p.resize();
121 :
122 0 : nMajorCycles=0;
123 :
124 0 : itsDataLoopPerMapper=false;
125 :
126 0 : }
127 :
128 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
129 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
130 :
131 0 : SynthesisImager::~SynthesisImager()
132 : {
133 0 : LogIO os( LogOrigin("SynthesisImager","destructor",WHERE) );
134 0 : os << LogIO::DEBUG1 << "SynthesisImager destroyed" << LogIO::POST;
135 0 : cleanupTempFiles();
136 0 : if(rvi_p) delete rvi_p;
137 0 : rvi_p=NULL;
138 : // cerr << "IN DESTR"<< endl;
139 : // VisModelData::listModel(mss4vi_p[0]);
140 :
141 0 : SynthesisUtilMethods::getResource("End SynthesisImager");
142 0 : }
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 : }
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 : }
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 : }
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 : }
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 : }
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 : }
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 : }
501 0 : os << "selected spw " << blockSpw_p[msin] << " start channels " << blockStart_p[msin] << " nchannels "<< blockNChan_p[msin] << LogIO::POST;
502 :
503 : }
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 : }
545 0 : catch(AipsError &x)
546 : {
547 0 : unlockMSs();
548 0 : throw( AipsError("Error in selectData() : "+x.getMesg()) );
549 : }
550 :
551 0 : return true;
552 :
553 : }
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 : }
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 : }
737 0 : catch(AipsError &x)
738 : {
739 0 : os << "Error in building Coordinate System and Image Shape : " << x.getMesg() << LogIO::EXCEPTION;
740 : }
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 : }
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 : }
780 :
781 0 : return true;
782 : }
783 :
784 0 : Long SynthesisImager::estimateRAM(){
785 0 : Long mem=0;
786 0 : LogIO os( LogOrigin("SynthesisImager","estimateRAM",WHERE) );
787 0 : if(itsMappers.nMappers()==0)
788 0 : os << "SynthesisImage has not been setup to get an estimate of memory usage" << LogIO::WARN << LogIO::POST;
789 0 : mem=itsMappers.estimateRAM();
790 0 : return mem;
791 : }
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 : }
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 : }
894 0 : return dataSel_p;
895 :
896 : }
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 0 : Record SynthesisImager::executeMajorCycle(const Record& controlRecord)
960 : {
961 0 : LogIO os( LogOrigin("SynthesisImager","executeMajorCycle",WHERE) );
962 :
963 0 : nMajorCycles++;
964 0 : if(controlRecord.isDefined("nmajorcycles"))
965 0 : controlRecord.get("nmajorcycles", nMajorCycles);
966 0 : Record outRec=controlRecord;
967 0 : Bool lastcycle=false;
968 :
969 :
970 0 : if( controlRecord.isDefined("lastcycle") )
971 : {
972 0 : controlRecord.get( "lastcycle" , lastcycle );
973 : //cout << "lastcycle : " << lastcycle << endl;
974 : }
975 : //else {cout << "No lastcycle" << endl;}
976 :
977 0 : os << "----------------------------------------------------------- Run ";
978 : //if (lastcycle) os << "(Last) " ;
979 0 : os << "Major Cycle " << nMajorCycles << " -------------------------------------" << LogIO::POST;
980 :
981 : try
982 : {
983 0 : if( (itsMaxShape[3] > 1 || impars_p.mode.contains("cube"))&& doingCubeGridding_p ){/// and valid ftmachines
984 0 : runMajorCycleCube(false, controlRecord);
985 : }
986 : else{
987 0 : if( itsDataLoopPerMapper == false )
988 0 : { runMajorCycle(false, lastcycle);}
989 : else
990 0 : { runMajorCycle2(false, lastcycle);}
991 :
992 : }
993 0 : if(lastcycle){
994 0 : String mess="In Major Cycle";
995 0 : 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 0 : std::vector<String> tmpfiles=itsMappers.cleanupTempFiles(mess);
1000 0 : outRec.define("tempfilenames", Vector<String>(tmpfiles));
1001 : }
1002 0 : itsMappers.releaseImageLocks();
1003 :
1004 : }
1005 0 : catch(AipsError &x)
1006 : {
1007 0 : throw( AipsError("Error in running Major Cycle : "+x.getMesg()) );
1008 : }
1009 0 : return outRec;
1010 : }// end of executeMajorCycle
1011 : //////////////////////////////////////////////
1012 : /////////////////////////////////////////////
1013 0 : void SynthesisImager::cleanupTempFiles(){
1014 0 : LogIO os( LogOrigin("SynthesisImager","cleanupTempFiles",WHERE) );
1015 0 : 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 : }
1024 : }
1025 :
1026 :
1027 : }
1028 :
1029 :
1030 0 : }
1031 0 : Record SynthesisImager::makePSF()
1032 : {
1033 0 : Record outRec=Record();
1034 0 : LogIO os( LogOrigin("SynthesisImager","makePSF",WHERE) );
1035 :
1036 0 : os << "----------------------------------------------------------- Make PSF ---------------------------------------------" << LogIO::POST;
1037 :
1038 : try
1039 : {
1040 0 : if( (itsMaxShape[3] >1 || impars_p.mode.contains("cube")) && doingCubeGridding_p){///and valid ftmachines
1041 0 : runMajorCycleCube(true);
1042 : }
1043 : else{
1044 0 : if( itsDataLoopPerMapper == false )
1045 0 : {runMajorCycle(true, false);}
1046 : else
1047 0 : {runMajorCycle2(true, false);}
1048 : }
1049 : // makeImage();
1050 :
1051 0 : String mess("PSF wieghting scheme");
1052 0 : std::vector<String> tmpfiles=itsMappers.cleanupTempFiles(mess);
1053 0 : outRec.define("tempfilenames", Vector<String>(tmpfiles));
1054 0 : itsMappers.releaseImageLocks();
1055 :
1056 : }
1057 0 : catch(AipsError &x)
1058 : {
1059 0 : throw( AipsError("Error in making PSF : "+x.getMesg()) );
1060 : }
1061 0 : return outRec;
1062 : }
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 : }
1203 : else {
1204 : //this->unlock();
1205 : os << LogIO::SEVERE << "Unknown weighting " << type
1206 0 : << LogIO::EXCEPTION;
1207 0 : return false;
1208 : }
1209 : }
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 : }
1225 0 : catch(AipsError &x)
1226 : {
1227 0 : throw( AipsError("Error in Weighting : "+x.getMesg()) );
1228 : }
1229 :
1230 0 : return true;
1231 : }
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 : }
1263 0 : catch (AipsError &x)
1264 : {
1265 0 : throw(AipsError("In getWeightDensity : "+x.getMesg()));
1266 : }
1267 0 : return outname;
1268 : }
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 : }
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 : }
1298 0 : return true;
1299 : }
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 0 : CountedPtr<SIImageStore> SynthesisImager::createIMStore(String imageName,
1311 : CoordinateSystem& cSys,
1312 : IPosition imShape,
1313 : const Bool overwrite,
1314 : MSColumns& msc,
1315 : String mappertype,
1316 : uInt ntaylorterms,
1317 : Quantity distance,
1318 : const TcleanProcessingInfo& procInfo,
1319 : uInt facets,
1320 : Bool useweightimage,
1321 : const Vector<String> &startmodel)
1322 : {
1323 0 : LogIO os( LogOrigin("SynthesisImager","createIMStore",WHERE) );
1324 :
1325 0 : CountedPtr<SIImageStore> imstor;
1326 :
1327 : try
1328 : {
1329 : // Prepare miscellaneous image information
1330 0 : auto objectName = msc.field().name()(msc.fieldId()(0));
1331 : ///// misc info for ImageStore. This will go to the 'miscinfo' table keyword
1332 0 : Record miscInfo;
1333 0 : auto telescop=msc.observation().telescopeName()(0);
1334 0 : miscInfo.define("INSTRUME", telescop);
1335 0 : miscInfo.define("distance", distance.get("m").getValue());
1336 0 : miscInfo.define("mpiprocs", procInfo.mpiprocs);
1337 0 : miscInfo.define("chnchnks", procInfo.chnchnks);
1338 0 : miscInfo.define("memreq", procInfo.memreq);
1339 0 : miscInfo.define("memavail", procInfo.memavail);
1340 :
1341 0 : if( mappertype=="default" || mappertype=="imagemosaic" )
1342 : {
1343 0 : imstor = std::make_shared<SIImageStore>(imageName, cSys, imShape, objectName,
1344 : miscInfo, overwrite,
1345 0 : (useweightimage || (mappertype=="imagemosaic")
1346 0 : ));
1347 : }
1348 0 : else if (mappertype == "multiterm" ) // Currently does not support imagemosaic.
1349 : {
1350 : // upcast with shared_ptr and then assign to CountedPtr<SIImageStore>
1351 : std::shared_ptr<SIImageStore> multiTermStore =
1352 0 : std::make_shared<SIImageStoreMultiTerm>(imageName, cSys, imShape,
1353 : objectName, miscInfo, facets,
1354 0 : overwrite, ntaylorterms, useweightimage);
1355 0 : imstor = multiTermStore;
1356 : }
1357 : else
1358 : {
1359 0 : throw(AipsError("Internal Error : Invalid mapper type in SynthesisImager::createIMStore"));
1360 : }
1361 :
1362 : // Get polRep from 'msc' here, and send to imstore.
1363 0 : StokesImageUtil::PolRep polRep(StokesImageUtil::CIRCULAR);
1364 0 : Vector<String> polType=msc.feed().polarizationType()(0);
1365 0 : if (polType(0)!="X" && polType(0)!="Y" && polType(0)!="R" && polType(0)!="L") {
1366 : os << LogIO::WARN << "Unknown stokes types in feed table: ["
1367 0 : << polType(0) << ", " << polType(1) << "]" << endl
1368 0 : << "Results open to question!" << LogIO::POST;
1369 : }
1370 :
1371 0 : if (polType(0)=="X" || polType(0)=="Y") {
1372 0 : polRep=StokesImageUtil::LINEAR;
1373 0 : os << LogIO::DEBUG1 << "Preferred polarization representation is linear" << LogIO::POST;
1374 : }
1375 : else {
1376 0 : polRep=StokesImageUtil::CIRCULAR;
1377 0 : os << LogIO::DEBUG1 << "Preferred polarization representation is circular" << LogIO::POST;
1378 : }
1379 : /// end of reading polRep info
1380 :
1381 : ///////// Send this info into ImageStore.
1382 0 : imstor->setDataPolFrame(polRep);
1383 :
1384 : ///////// Set Starting model if it exists.
1385 : //cout << "In SI, set starting model to : " << startmodel << endl;
1386 0 : if( startmodel.nelements()>0 )
1387 : {
1388 0 : imstor->setModelImage( startmodel );
1389 : }
1390 :
1391 0 : imstor->releaseLocks();
1392 :
1393 : }
1394 0 : catch(AipsError &x)
1395 : {
1396 0 : throw(AipsError("Error in createImStore : " + x.getMesg() ) );
1397 : }
1398 :
1399 :
1400 0 : return imstor;
1401 : }
1402 :
1403 0 : CountedPtr<SIMapper> SynthesisImager::createSIMapper(String mappertype,
1404 : CountedPtr<SIImageStore> imagestore,
1405 : CountedPtr<FTMachine> ftmachine,
1406 : CountedPtr<FTMachine> iftmachine,
1407 : uInt /*ntaylorterms*/)
1408 : {
1409 0 : LogIO os( LogOrigin("SynthesisImager","createSIMapper",WHERE) );
1410 :
1411 0 : CountedPtr<SIMapper> localMapper;
1412 :
1413 : try
1414 : {
1415 :
1416 0 : if( mappertype == "default" || mappertype == "multiterm" )
1417 : {
1418 0 : localMapper = new SIMapper( imagestore, ftmachine, iftmachine );
1419 : }
1420 0 : else if( mappertype == "imagemosaic") // || mappertype == "mtimagemosaic" )
1421 : {
1422 0 : localMapper = new SIMapperImageMosaic( imagestore, ftmachine, iftmachine );
1423 : }
1424 : else
1425 : {
1426 0 : throw(AipsError("Unknown mapper type : " + mappertype));
1427 : }
1428 :
1429 : }
1430 0 : catch(AipsError &x) {
1431 0 : throw(AipsError("Error in createSIMapper : " + x.getMesg() ) );
1432 : }
1433 0 : return localMapper;
1434 : }
1435 :
1436 :
1437 0 : Block<CountedPtr<SIImageStore> > SynthesisImager::createFacetImageStoreList(
1438 : CountedPtr<SIImageStore> imagestore,
1439 : Int facets)
1440 : {
1441 0 : LogIO os(LogOrigin("SynthesisImager", "createFacetImageStoreList"));
1442 :
1443 0 : os << "Creating " << facets*facets << " facets in total " << LogIO::POST;
1444 :
1445 0 : Block<CountedPtr<SIImageStore> > facetList( facets*facets );
1446 :
1447 0 : if( facets==1 ) { facetList[0] = imagestore; return facetList; }
1448 :
1449 : // Remember, only the FIRST field in each run can have facets. So, check for this.
1450 0 : if( ! unFacettedImStore_p.null() ) {
1451 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. ") );
1452 : }
1453 :
1454 0 : unFacettedImStore_p = imagestore;
1455 0 : facetsStore_p = facets;
1456 :
1457 : // Note : facets : Number of facets on a side.
1458 : // Note : facet : index from range(0, facets*facets)
1459 0 : for (Int facet=0; facet< facets*facets; ++facet){
1460 0 : facetList[facet] = unFacettedImStore_p->getSubImageStore(facet, facets);
1461 : }
1462 :
1463 0 : return facetList;
1464 : }
1465 :
1466 : /*
1467 : void SynthesisImager::setPsfFromOneFacet()
1468 : {
1469 :
1470 : if( unFacettedImStore_p.null() ){
1471 : throw(AipsError("Internal Error in SynthesisImager : Setting PSF on Null unfacettedimage"));
1472 : }
1473 :
1474 : // Copy the PSF from one facet to the center of the full image, to use for the minor cycle
1475 : //
1476 : // This code segment will work for single and multi-term
1477 : // - for single term, the index will always be 0, and SIImageStore's access functions know this.
1478 : // - for multi-term, the index will be the taylor index and SIImageStoreMultiTerm knows this.
1479 : {
1480 : IPosition shape=(unFacettedImStore_p->psf(0))->shape();
1481 : IPosition blc(4, 0, 0, 0, 0);
1482 : IPosition trc=shape-1;
1483 : for(uInt tix=0; tix<2 * unFacettedImStore_p->getNTaylorTerms() - 1; tix++)
1484 : {
1485 : TempImage<Float> onepsf((itsMappers.imageStore(0)->psf(tix))->shape(),
1486 : (itsMappers.imageStore(0)->psf(tix))->coordinates());
1487 : onepsf.copyData(*(itsMappers.imageStore(0)->psf(tix)));
1488 : //now set the original to 0 as we have a copy of one facet psf
1489 : (unFacettedImStore_p->psf(tix))->set(0.0);
1490 : blc[0]=(shape[0]-(onepsf.shape()[0]))/2;
1491 : trc[0]=onepsf.shape()[0]+blc[0]-1;
1492 : blc[1]=(shape[1]-(onepsf.shape()[1]))/2;
1493 : trc[1]=onepsf.shape()[1]+blc[1]-1;
1494 : Slicer sl(blc, trc, Slicer::endIsLast);
1495 : SubImage<Float> sub(*(unFacettedImStore_p->psf(tix)), sl, true);
1496 : sub.copyData(onepsf);
1497 : }
1498 : }
1499 :
1500 : //cout << "In setPsfFromOneFacet : sumwt : " << unFacettedImStore_p->sumwt()->get() << endl;
1501 :
1502 : }
1503 : */
1504 :
1505 0 : Block<CountedPtr<SIImageStore> > SynthesisImager::createChanChunkImageStoreList(
1506 : CountedPtr<SIImageStore> imagestore,
1507 : Int chanchunks)
1508 : {
1509 0 : LogIO os(LogOrigin("SynthesisImager", "createChanChunkImageStoreList"));
1510 :
1511 0 : Int extrachunks=0;
1512 0 : Int chunksize=imagestore->getShape()[3]/chanchunks;
1513 0 : Int rem=imagestore->getShape()[3] % chanchunks;
1514 : ///Avoid an extra chunk with 1 channel as it cause bumps in linear interpolation
1515 : ///See CAS-12625
1516 0 : while((rem==1) && (chunksize >1)){
1517 0 : chanchunks +=1;
1518 0 : chunksize=imagestore->getShape()[3]/chanchunks;
1519 0 : rem=imagestore->getShape()[3] % chanchunks;
1520 : }
1521 :
1522 :
1523 :
1524 0 : if( rem>0 )
1525 : {
1526 : // os << LogIO::WARN << "chanchunks ["+String::toString(chanchunks)+"] is not a divisor of nchan ["+String::toString(imagestore->getShape()[3])+"].";
1527 : // 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 ;
1528 :
1529 0 : if( rem < chunksize ) extrachunks=1;
1530 : else
1531 : {
1532 0 : extrachunks=rem/chunksize;
1533 0 : if( rem % chunksize > 0 ) extrachunks += 1;
1534 : }
1535 : }
1536 :
1537 :
1538 0 : os << "Creating " << chanchunks +extrachunks << " reference subCubes (channel chunks) for gridding " << LogIO::POST;
1539 :
1540 0 : Block<CountedPtr<SIImageStore> > chunkList( chanchunks + extrachunks );
1541 :
1542 0 : if( chanchunks==1 ) { chunkList[0] = imagestore; return chunkList; }
1543 :
1544 : // Remember, only the FIRST field in each run can have chanchunks. So, check for this.
1545 0 : if( ! unChanChunkedImStore_p.null() ) {
1546 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. ") );
1547 : }
1548 :
1549 0 : unChanChunkedImStore_p = imagestore;
1550 0 : chanChunksStore_p = chanchunks;
1551 :
1552 0 : for (Int chunk=0; chunk< chanchunks+extrachunks; ++chunk){
1553 0 : chunkList[chunk] = unChanChunkedImStore_p->getSubImageStore(0,1,chunk, chanchunks);
1554 : }
1555 :
1556 0 : return chunkList;
1557 : }
1558 :
1559 :
1560 :
1561 :
1562 0 : void SynthesisImager::appendToMapperList(String imagename,
1563 : CoordinateSystem& csys,
1564 : IPosition imshape,
1565 : CountedPtr<FTMachine>& ftm,
1566 : CountedPtr<FTMachine>& iftm,
1567 : Quantity distance,
1568 : Int facets,
1569 : Int chanchunks,
1570 : const Bool overwrite,
1571 : String mappertype,
1572 : Float padding,
1573 : uInt ntaylorterms,
1574 : const Vector<String> &startmodel)
1575 : {
1576 0 : LogIO log_l(LogOrigin("SynthesisImager", "appendToMapperList(ftm)"));
1577 : //---------------------------------------------
1578 : // Some checks..
1579 0 : if(facets > 1 && itsMappers.nMappers() > 0)
1580 0 : log_l << "Facetted image has to be the first of multifields" << LogIO::EXCEPTION;
1581 :
1582 0 : TcleanProcessingInfo procInfo;
1583 0 : if(chanchunks<1)
1584 : {
1585 0 : log_l << "Automatically calculate chanchunks";
1586 0 : log_l << " using imshape : " << imshape << LogIO::POST;
1587 :
1588 : // Do calculation here.
1589 : // This runs once per image field (for multi-field imaging)
1590 : // This runs once per cube partition, and will see only its own partition's shape
1591 0 : chanchunks=1;
1592 :
1593 0 : CompositeNumber cn(uInt(imshape[0] * 2));
1594 : // heuristic factors multiplied to imshape based on gridder
1595 0 : size_t fudge_factor = 15;
1596 0 : if (ftm->name()=="MosaicFTNew") {
1597 0 : fudge_factor = 20;
1598 : }
1599 0 : else if (ftm->name()=="GridFT") {
1600 0 : fudge_factor = 9;
1601 : }
1602 :
1603 0 : Double required_mem = fudge_factor * sizeof(Float);
1604 0 : for (size_t i = 0; i < imshape.nelements(); i++) {
1605 : // gridding pads image and increases to composite number
1606 0 : if (i < 2) {
1607 0 : required_mem *= cn.nextLargerEven(Int(padding*Float(imshape[i])-0.5));
1608 : }
1609 : else {
1610 0 : required_mem *= imshape[i];
1611 : }
1612 : }
1613 :
1614 : // get number of tclean processes running on the same machine
1615 0 : size_t nlocal_procs = 1;
1616 0 : if (getenv("OMPI_COMM_WORLD_LOCAL_SIZE")) {
1617 0 : std::stringstream ss(getenv("OMPI_COMM_WORLD_LOCAL_SIZE"));
1618 0 : ss >> nlocal_procs;
1619 : }
1620 : // assumes all processes need the same amount of memory
1621 0 : required_mem *= nlocal_procs;
1622 : Double usr_memfrac, usr_mem;
1623 0 : AipsrcValue<Double>::find(usr_memfrac, "system.resources.memfrac", 80.);
1624 0 : AipsrcValue<Double>::find(usr_mem, "system.resources.memory", -1024.);
1625 : Double memory_avail;
1626 0 : if (usr_mem > 0.) {
1627 0 : memory_avail = usr_mem * 1024. * 1024.;
1628 : }
1629 : else {
1630 0 : memory_avail = Double(HostInfo::memoryFree()) * (usr_memfrac / 100.) * 1024.;
1631 : }
1632 :
1633 : // compute required chanchunks to fit into the available memory
1634 0 : chanchunks = (int)std::ceil((Double)required_mem / memory_avail);
1635 0 : if (imshape.nelements() == 4 && imshape[3] < chanchunks) {
1636 0 : chanchunks = imshape[3];
1637 : // TODO make chanchunks a divisor of nchannels?
1638 : }
1639 0 : chanchunks = chanchunks < 1 ? 1 : chanchunks;
1640 :
1641 0 : log_l << "Required memory " << required_mem / nlocal_procs / 1024. / 1024. / 1024.
1642 0 : << "\nAvailable memory " << memory_avail / 1024. / 1024 / 1024.
1643 : << " (rc: memory fraction " << usr_memfrac << "% memory " << usr_mem / 1024.
1644 : << ")\n" << nlocal_procs << " other processes on node\n"
1645 0 : << "Setting chanchunks to " << chanchunks << LogIO::POST;
1646 :
1647 0 : procInfo.mpiprocs = nlocal_procs;
1648 0 : procInfo.chnchnks = chanchunks;
1649 0 : const float toGB = 1024.0 * 1024.0 * 1024.0;
1650 0 : procInfo.memavail = memory_avail / toGB;
1651 0 : procInfo.memreq = required_mem / toGB;
1652 : }
1653 :
1654 0 : if( imshape.nelements()==4 && imshape[3]<chanchunks )
1655 : {
1656 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;
1657 : }
1658 :
1659 0 : if(chanchunks > 1 && itsMappers.nMappers() > 0)
1660 0 : log_l << "Channel chunking is currently not supported with multi(outlier)-fields. Please submit a feature request if needed." << LogIO::EXCEPTION;
1661 :
1662 0 : if(chanchunks > 1) itsDataLoopPerMapper=true;
1663 :
1664 0 : AlwaysAssert( ( ( ! (ftm->name()=="MosaicFTNew" && mappertype=="imagemosaic") ) &&
1665 : ( ! (ftm->name()=="AWProjectWBFTNew" && mappertype=="imagemosaic") )) ,
1666 : AipsError );
1667 : //---------------------------------------------
1668 :
1669 : // Create the ImageStore object
1670 0 : CountedPtr<SIImageStore> imstor;
1671 0 : MSColumns msc(mss4vi_p[0]);
1672 0 : imstor = createIMStore(imagename, csys, imshape, overwrite, msc, mappertype, ntaylorterms, distance, procInfo, facets, iftm->useWeightImage(), startmodel );
1673 :
1674 : // Create the Mappers
1675 0 : if( facets<2 && chanchunks<2) // One facet. Just add the above imagestore to the mapper list.
1676 : {
1677 0 : itsMappers.addMapper( createSIMapper( mappertype, imstor, ftm, iftm, ntaylorterms) );
1678 : }
1679 : else // This field is facetted. Make a list of reference imstores, and add all to the mapper list.
1680 : {
1681 :
1682 0 : if ( facets>1 && chanchunks==1 )
1683 : {
1684 : // Make and connect the list.
1685 0 : Block<CountedPtr<SIImageStore> > imstorList = createFacetImageStoreList( imstor, facets );
1686 0 : for( uInt facet=0; facet<imstorList.nelements(); facet++)
1687 : {
1688 0 : CountedPtr<FTMachine> new_ftm, new_iftm;
1689 0 : if(facet==0){ new_ftm = ftm; new_iftm = iftm; }
1690 0 : else{ new_ftm=ftm->cloneFTM(); new_iftm=iftm->cloneFTM(); }
1691 0 : itsMappers.addMapper(createSIMapper( mappertype, imstorList[facet], new_ftm, new_iftm, ntaylorterms));
1692 0 : }
1693 : }// facets
1694 0 : else if ( facets==1 && chanchunks>1 )
1695 : {
1696 : // Make and connect the list.
1697 0 : Block<CountedPtr<SIImageStore> > imstorList = createChanChunkImageStoreList( imstor, chanchunks );
1698 0 : for( uInt chunk=0; chunk<imstorList.nelements(); chunk++)
1699 : {
1700 0 : CountedPtr<FTMachine> new_ftm, new_iftm;
1701 0 : if(chunk==0){ new_ftm = ftm; new_iftm = iftm; }
1702 0 : else{ new_ftm=ftm->cloneFTM(); new_iftm=iftm->cloneFTM(); }
1703 0 : itsMappers.addMapper(createSIMapper( mappertype, imstorList[chunk], new_ftm, new_iftm, ntaylorterms));
1704 0 : }
1705 : }// chanchunks
1706 : else
1707 : {
1708 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. ") );
1709 : }
1710 :
1711 : }// facets or chunks
1712 :
1713 0 : }
1714 :
1715 : /////////////////////////
1716 : /*
1717 : Bool SynthesisImager::toUseWeightImage(CountedPtr<FTMachine>& ftm, String mappertype)
1718 : {
1719 : if( (ftm->name() == "GridFT" || ftm->name() == "WProjectFT")&&(mappertype!="imagemosaic") )
1720 : { return false; }
1721 : else
1722 : { return true; }
1723 : }
1724 : */
1725 :
1726 : // Make the FT-Machine and related objects (cfcache, etc.)
1727 0 : void SynthesisImager::createFTMachine(CountedPtr<FTMachine>& theFT,
1728 : CountedPtr<FTMachine>& theIFT,
1729 : const String& ftname,
1730 : const uInt nTaylorTerms,
1731 : const String mType,
1732 : const Int facets, //=1
1733 : //------------------------------
1734 : const Int wprojplane, //=1,
1735 : const Float padding, //=1.0,
1736 : const Bool useAutocorr, //=false,
1737 : const Bool useDoublePrec, //=true,
1738 : const String gridFunction, //=String("SF"),
1739 : //------------------------------
1740 : const Bool aTermOn, //= true,
1741 : const Bool psTermOn, //= true,
1742 : const Bool mTermOn, //= false,
1743 : const Bool wbAWP, //= true,
1744 : const String cfCache, //= "",
1745 : const Bool usePointing, //= false,
1746 : const Bool doPBCorr, //= true,
1747 : const Bool conjBeams, //= true,
1748 : const Float computePAStep, //=360.0
1749 : const Float rotatePAStep, //=5.0
1750 : const String interpolation, //="linear"
1751 : const Bool freqFrameValid, //=true
1752 : const Int cache, //=1000000000,
1753 : const Int tile, //=16
1754 : const String stokes, //=I
1755 : const String imageNamePrefix
1756 : )
1757 :
1758 : {
1759 0 : LogIO os( LogOrigin("SynthesisImager","createFTMachine",WHERE));
1760 :
1761 0 : if(ftname=="gridft"){
1762 0 : if(facets >1){
1763 0 : theFT=new GridFT(cache, tile, gridFunction, mLocation_p, phaseCenter_p, padding, useAutocorr, useDoublePrec);
1764 0 : theIFT=new GridFT(cache, tile, gridFunction, mLocation_p, phaseCenter_p, padding, useAutocorr, useDoublePrec);
1765 :
1766 : }
1767 : else{
1768 0 : theFT=new GridFT(cache, tile, gridFunction, mLocation_p, padding, useAutocorr, useDoublePrec);
1769 0 : theIFT=new GridFT(cache, tile, gridFunction, mLocation_p, padding, useAutocorr, useDoublePrec);
1770 : }
1771 : }
1772 0 : else if(ftname== "wprojectft"){
1773 0 : Double maxW=-1.0;
1774 0 : Double minW=-1.0;
1775 0 : Double rmsW=-1.0;
1776 0 : if(wprojplane <1)
1777 0 : WProjectFT::wStat(*rvi_p, minW, maxW, rmsW);
1778 0 : if(facets >1){
1779 0 : theFT=new WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
1780 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1781 0 : theIFT=new WProjectFT(wprojplane, phaseCenter_p, mLocation_p,
1782 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1783 : }
1784 : else{
1785 0 : theFT=new WProjectFT(wprojplane, mLocation_p,
1786 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1787 0 : theIFT=new WProjectFT(wprojplane, mLocation_p,
1788 0 : cache/2, tile, useAutocorr, padding, useDoublePrec, minW, maxW, rmsW);
1789 : }
1790 0 : CountedPtr<WPConvFunc> sharedconvFunc=static_cast<WProjectFT &>(*theFT).getConvFunc();
1791 : //static_cast<WProjectFT &>(*theFT).setConvFunc(sharedconvFunc);
1792 0 : static_cast<WProjectFT &>(*theIFT).setConvFunc(sharedconvFunc);
1793 : }
1794 0 : else if ((ftname == "awprojectft") || (ftname== "mawprojectft") || (ftname == "protoft")) {
1795 0 : createAWPFTMachine(theFT, theIFT, ftname, facets, wprojplane,
1796 : padding, useAutocorr, useDoublePrec, gridFunction,
1797 : aTermOn, psTermOn, mTermOn, wbAWP, cfCache,
1798 : usePointing, doPBCorr, conjBeams, computePAStep,
1799 : rotatePAStep, cache,tile,imageNamePrefix);
1800 : }
1801 0 : else if ( ftname == "mosaic" || ftname== "mosft" || ftname == "mosaicft" || ftname== "MosaicFT"){
1802 :
1803 0 : createMosFTMachine(theFT, theIFT, padding, useAutocorr, useDoublePrec, rotatePAStep, stokes, conjBeams);
1804 : }
1805 : else
1806 : {
1807 0 : throw( AipsError( "Invalid FTMachine name : " + ftname ) );
1808 : }
1809 : /* else if(ftname== "MosaicFT"){
1810 :
1811 : }*/
1812 :
1813 :
1814 :
1815 : ///////// Now, clone and pack the chosen FT into a MultiTermFT if needed.
1816 0 : if( mType=="multiterm" )
1817 : {
1818 0 : AlwaysAssert( nTaylorTerms>=1 , AipsError );
1819 :
1820 0 : CountedPtr<FTMachine> theMTFT = new MultiTermFTNew( theFT , nTaylorTerms, true/*forward*/ );
1821 0 : CountedPtr<FTMachine> theMTIFT = new MultiTermFTNew( theIFT , nTaylorTerms, false/*forward*/ );
1822 :
1823 0 : theFT = theMTFT;
1824 0 : theIFT = theMTIFT;
1825 : }
1826 :
1827 :
1828 :
1829 :
1830 : ////// Now, set the SkyJones if needed, and if not internally generated.
1831 0 : if( mType=="imagemosaic" &&
1832 0 : (ftname != "awprojectft" && ftname != "mawprojectft" && ftname != "proroft") )
1833 : {
1834 0 : CountedPtr<SkyJones> vp;
1835 0 : MSColumns msc(mss4vi_p[0]);
1836 0 : Quantity parang(0.0,"deg");
1837 0 : Quantity skyposthreshold(0.0,"deg");
1838 0 : vp = new VPSkyJones(msc, true, parang, BeamSquint::NONE,skyposthreshold);
1839 :
1840 0 : Vector<CountedPtr<SkyJones> > skyJonesList(1);
1841 0 : skyJonesList(0) = vp;
1842 0 : theFT->setSkyJones( skyJonesList );
1843 0 : theIFT->setSkyJones( skyJonesList );
1844 :
1845 : }
1846 :
1847 : //// For mode=cubedata, set the freq frame to invalid..
1848 : // get this info from buildCoordSystem
1849 0 : Vector<Int> tspws(0);
1850 : //theFT->setSpw( tspws, false );
1851 : //theIFT->setSpw( tspws, false );
1852 0 : theFT->setSpw( tspws, freqFrameValid );
1853 0 : theIFT->setSpw( tspws, freqFrameValid );
1854 :
1855 : //// Set interpolation mode
1856 0 : theFT->setFreqInterpolation( interpolation );
1857 0 : theIFT->setFreqInterpolation( interpolation );
1858 :
1859 : //channel selections from spw param
1860 0 : theFT->setSpwChanSelection(chanSel_p);
1861 0 : theIFT->setSpwChanSelection(chanSel_p);
1862 0 : }
1863 :
1864 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1865 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1866 :
1867 0 : void SynthesisImager::createAWPFTMachine(CountedPtr<FTMachine>& theFT, CountedPtr<FTMachine>& theIFT,
1868 : const String&,// ftmName,
1869 : const Int,// facets, //=1
1870 : //------------------------------
1871 : const Int wprojPlane, //=1,
1872 : const Float,// padding, //=1.0,
1873 : const Bool,// useAutocorr, //=false,
1874 : const Bool useDoublePrec, //=true,
1875 : const String,// gridFunction, //=String("SF"),
1876 : //------------------------------
1877 : const Bool aTermOn, //= true,
1878 : const Bool psTermOn, //= true,
1879 : const Bool mTermOn, //= false,
1880 : const Bool wbAWP, //= true,
1881 : const String cfCache, //= "",
1882 : const Bool usePointing, //= false,
1883 : const Bool doPBCorr, //= true,
1884 : const Bool conjBeams, //= true,
1885 : const Float computePAStep, //=360.0
1886 : const Float rotatePAStep, //=5.0
1887 : const Int cache, //=1000000000,
1888 : const Int tile, //=16
1889 : const String imageNamePrefix
1890 : )
1891 :
1892 : {
1893 0 : LogIO os( LogOrigin("SynthesisImager","createAWPFTMachine",WHERE));
1894 :
1895 0 : if (wprojPlane<=1)
1896 : {
1897 : os << LogIO::NORMAL
1898 : << "You are using wprojplanes=1. Doing co-planar imaging (no w-projection needed)"
1899 0 : << LogIO::POST;
1900 0 : os << LogIO::NORMAL << "Performing WBA-Projection" << LogIO::POST; // Loglevel PROGRESS
1901 : }
1902 : // if((wprojPlane>1)&&(wprojPlane<64))
1903 : // {
1904 : // os << LogIO::WARN
1905 : // << "No. of w-planes set too low for W projection - recommend at least 128"
1906 : // << LogIO::POST;
1907 : // os << LogIO::NORMAL << "Performing WBAW-Projection" << LogIO::POST; // Loglevel PROGRESS
1908 : // }
1909 :
1910 : // CountedPtr<ATerm> apertureFunction = createTelescopeATerm(mss4vi_p[0], aTermOn);
1911 : // CountedPtr<PSTerm> psTerm = new PSTerm();
1912 : // CountedPtr<WTerm> wTerm = new WTerm();
1913 :
1914 : // //
1915 : // // Selectively switch off CFTerms.
1916 : // //
1917 : // if (aTermOn == false) {apertureFunction->setOpCode(CFTerms::NOOP);}
1918 : // if (psTermOn == false) psTerm->setOpCode(CFTerms::NOOP);
1919 :
1920 : // //
1921 : // // Construct the CF object with appropriate CFTerms.
1922 : // //
1923 : // CountedPtr<ConvolutionFunction> tt;
1924 : // tt = AWProjectFT::makeCFObject(aTermOn, psTermOn, true, mTermOn, wbAWP);
1925 : // CountedPtr<ConvolutionFunction> awConvFunc;
1926 : // // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm, !wbAWP);
1927 : // if ((ftmName=="mawprojectft") || (mTermOn))
1928 : // awConvFunc = new AWConvFuncEPJones(apertureFunction,psTerm,wTerm,wbAWP);
1929 : // else
1930 : // awConvFunc = new AWConvFunc(apertureFunction,psTerm,wTerm,wbAWP);
1931 :
1932 0 : MSObservationColumns msoc(mss4vi_p[0].observation());
1933 0 : String telescopeName=msoc.telescopeName()(0);
1934 : CountedPtr<ConvolutionFunction> awConvFunc = AWProjectFT::makeCFObject(telescopeName,
1935 : aTermOn,
1936 : psTermOn, (wprojPlane > 1),
1937 0 : mTermOn, wbAWP, conjBeams);
1938 : //
1939 : // Construct the appropriate re-sampler.
1940 : //
1941 0 : CountedPtr<VisibilityResamplerBase> visResampler;
1942 : // if (ftmName=="protoft") visResampler = new ProtoVR();
1943 : //elsef
1944 0 : visResampler = new AWVisResampler();
1945 : // CountedPtr<VisibilityResamplerBase> visResampler = new VisibilityResampler();
1946 :
1947 : //
1948 : // Construct and initialize the CF cache object.
1949 : //
1950 :
1951 :
1952 : // CountedPtr<CFCache> cfCacheObj = new CFCache();
1953 : // cfCacheObj->setCacheDir(cfCache.data());
1954 : // // cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
1955 : // cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
1956 : // cfCacheObj->initCache2();
1957 :
1958 0 : CountedPtr<CFCache> cfCacheObj;
1959 :
1960 :
1961 : //
1962 : // Finally construct the FTMachine with the CFCache, ConvFunc and
1963 : // Re-sampler objects.
1964 : //
1965 0 : Float pbLimit_l=1e-2;
1966 0 : theFT = new AWProjectWBFTNew(wprojPlane, cache/2,
1967 : cfCacheObj, awConvFunc,
1968 : visResampler,
1969 : /*true */usePointing, doPBCorr,
1970 : tile, computePAStep, pbLimit_l, true,conjBeams,
1971 0 : useDoublePrec);
1972 :
1973 0 : cfCacheObj = new CFCache();
1974 0 : cfCacheObj->setCacheDir(cfCache.data());
1975 : // Get the LAZYFILL setting from the user configuration. If not
1976 : // found, default to False.
1977 : //
1978 : // With lazy fill ON, CFCache loads the required CFs on-demand
1979 : // from the disk. And periodically triggers garbage collection to
1980 : // release CFs that aren't required immediately.
1981 0 : cfCacheObj->setLazyFill(SynthesisUtils::getenv("CFCache.LAZYFILL",1)==1);
1982 :
1983 : // cerr << "Setting wtImagePrefix to " << imageNamePrefix.c_str() << endl;
1984 0 : cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
1985 0 : cfCacheObj->initCache2();
1986 :
1987 0 : theFT->setCFCache(cfCacheObj);
1988 :
1989 :
1990 0 : Quantity rotateOTF(rotatePAStep,"deg");
1991 0 : static_cast<AWProjectWBFTNew &>(*theFT).setObservatoryLocation(mLocation_p);
1992 0 : static_cast<AWProjectWBFTNew &>(*theFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
1993 :
1994 : // theIFT = new AWProjectWBFT(wprojPlane, cache/2,
1995 : // cfCacheObj, awConvFunc,
1996 : // visResampler,
1997 : // /*true */usePointing, doPBCorr,
1998 : // tile, computePAStep, pbLimit_l, true,conjBeams,
1999 : // useDoublePrec);
2000 :
2001 : // static_cast<AWProjectWBFT &>(*theIFT).setObservatoryLocation(mLocation_p);
2002 : // static_cast<AWProjectWBFT &>(*theIFT).setPAIncrement(Quantity(computePAStep,"deg"),rotateOTF);
2003 :
2004 0 : theIFT = new AWProjectWBFTNew(static_cast<AWProjectWBFTNew &>(*theFT));
2005 :
2006 : //// Send in Freq info.
2007 0 : os << "Sending frequency selection information " << mssFreqSel_p << " to AWP FTM." << LogIO::POST;
2008 0 : theFT->setSpwFreqSelection( mssFreqSel_p );
2009 0 : theIFT->setSpwFreqSelection( mssFreqSel_p );
2010 :
2011 0 : }
2012 :
2013 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2014 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2015 :
2016 0 : ATerm* SynthesisImager::createTelescopeATerm(const MeasurementSet& ms, const Bool& isATermOn)
2017 : {
2018 0 : LogIO os(LogOrigin("SynthesisImager", "createTelescopeATerm",WHERE));
2019 :
2020 0 : if (!isATermOn) return new NoOpATerm();
2021 :
2022 0 : MSObservationColumns msoc(ms.observation());
2023 0 : String ObsName=msoc.telescopeName()(0);
2024 0 : if ((ObsName == "EVLA") || (ObsName == "VLA"))
2025 0 : return new EVLAAperture();
2026 : else
2027 : {
2028 0 : os << "Telescope name ('"+
2029 0 : ObsName+"') in the MS not recognized to create the telescope specific ATerm"
2030 0 : << LogIO::WARN;
2031 : }
2032 :
2033 0 : return NULL;
2034 : }
2035 :
2036 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2037 :
2038 0 : void SynthesisImager::getVPRecord(Record &rec, PBMath::CommonPB &kpb, String telescop)
2039 : {
2040 0 : LogIO os(LogOrigin("SynthesisImager", "getVPRecord",WHERE));
2041 :
2042 0 : VPManager *vpman=VPManager::Instance();
2043 0 : if( itsVpTable != String("") )
2044 : {
2045 0 : os << "Loading Voltage Pattern information from " << itsVpTable << LogIO::POST;
2046 0 : vpman->loadfromtable(itsVpTable);
2047 : }
2048 : else
2049 : {
2050 : // os << "Using Voltage Pattern currently set in the VPManager" << LogIO::POST;
2051 0 : os << "Using default Voltage Patterns from the VPManager" << LogIO::POST;
2052 0 : vpman->reset();
2053 : }
2054 :
2055 :
2056 :
2057 :
2058 : // PBMath::CommonPB kpb;
2059 0 : PBMath::enumerateCommonPB(telescop, kpb);
2060 : // Record rec;
2061 0 : vpman->getvp(rec, telescop);
2062 :
2063 : //ostringstream oss; rec.print(oss);
2064 : //os << "Using VP model : " << oss.str() << LogIO::POST;
2065 :
2066 0 : if(rec.empty()){
2067 0 : if((telescop=="EVLA")){
2068 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;
2069 0 : telescop="VLA";
2070 0 : vpman->getvp(rec, telescop);
2071 0 : kpb=PBMath::VLA;
2072 : }
2073 : else{
2074 : //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;
2075 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;
2076 0 : kpb=PBMath::UNKNOWN;
2077 0 : rec.define("name","COMMONPB");
2078 0 : rec.define("commonpb","NONE");
2079 : }
2080 : }
2081 :
2082 0 : }// get VPRecord
2083 :
2084 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2085 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*/){
2086 :
2087 0 : LogIO os(LogOrigin("SynthesisImager", "createMosFTMachine",WHERE));
2088 :
2089 0 : rvi_p->originChunks();
2090 0 : MSColumns msc(rvi_p->ms());
2091 0 : String telescop=msc.observation().telescopeName()(0);
2092 : ///Multi ms with different telescop
2093 0 : Bool multiTel=False;
2094 0 : for(rvi_p->originChunks(); rvi_p->moreChunks(); rvi_p->nextChunk()){
2095 0 : if(rvi_p->newMS() && telescop != MSColumns(rvi_p->ms()).observation().telescopeName()(0))
2096 0 : multiTel=True;
2097 : }
2098 0 : rvi_p->originChunks();
2099 :
2100 : PBMath::CommonPB kpb;
2101 0 : Record rec;
2102 0 : getVPRecord( rec, kpb, telescop );
2103 :
2104 : /*
2105 : VPManager *vpman=VPManager::Instance();
2106 :
2107 : if( itsVpTable != String("") ) vpman->loadfromtable(itsVpTable);
2108 :
2109 : PBMath::CommonPB kpb;
2110 : PBMath::enumerateCommonPB(telescop, kpb);
2111 : Record rec;
2112 : vpman->getvp(rec, telescop);
2113 : if(rec.empty()){
2114 : if((telescop=="EVLA")){
2115 : os << LogIO::WARN << "vpmanager does not list EVLA. Using VLA beam parameters" << LogIO::POST;
2116 : telescop="VLA";
2117 : vpman->getvp(rec, telescop);
2118 : kpb=PBMath::VLA;
2119 : }
2120 : else{
2121 : os << LogIO::SEVERE << "vpmanager does not have a beam "+telescop <<"\n You can use VPManager to define one" << LogIO::POST;
2122 : }
2123 : }
2124 :
2125 : */
2126 :
2127 0 : if(rec.empty()){os << LogIO::SEVERE << "Cannot proceed with mosaicft gridder without a valid PB model" << LogIO::POST; }
2128 :
2129 :
2130 0 : VPSkyJones* vps=NULL;
2131 0 : if(rec.asString("name")=="COMMONPB" && kpb !=PBMath::UNKNOWN ){
2132 0 : vps= new VPSkyJones(msc, true, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
2133 : /////Don't know which parameter has pb threshold cutoff that the user want
2134 : ////leaving at default
2135 : ////vps.setThreshold(minPB);
2136 :
2137 : }
2138 : else{
2139 0 : PBMath myPB(rec);
2140 0 : String whichPBMath;
2141 0 : PBMathInterface::namePBClass(myPB.whichPBClass(), whichPBMath);
2142 0 : os << "Using the PB defined by " << whichPBMath << " for beam calculation for telescope " << telescop << LogIO::POST;
2143 0 : vps= new VPSkyJones(telescop, myPB, Quantity(rotatePAStep, "deg"), BeamSquint::GOFIGURE, Quantity(360.0, "deg"));
2144 0 : kpb=PBMath::DEFAULT;
2145 : }
2146 :
2147 :
2148 : //cerr << "SImager tel " << ((vps) ? vps->telescope(): "NOTEL " ) << " multiTel " << multiTel << endl;
2149 :
2150 0 : theFT = new MosaicFTNew(vps, mLocation_p, stokes, 1000000000, 16, useAutoCorr,
2151 0 : useDoublePrec);
2152 0 : PBMathInterface::PBClass pbtype=((kpb==PBMath::EVLA) || multiTel)? PBMathInterface::COMMONPB: PBMathInterface::AIRY;
2153 0 : if(rec.asString("name")=="IMAGE")
2154 0 : pbtype=PBMathInterface::IMAGE;
2155 : ///Use Heterogenous array mode for the following
2156 : ///Added EVLA in it to use different beam models for different frequencies
2157 0 : if((kpb == PBMath::UNKNOWN) || (kpb==PBMath::OVRO) || (kpb==PBMath::ACA)
2158 0 : || (kpb==PBMath::ALMA) || (kpb==PBMath::EVLA) || multiTel){
2159 0 : CountedPtr<SimplePBConvFunc> mospb=new HetArrayConvFunc(pbtype, "");
2160 0 : static_cast<MosaicFTNew &>(*theFT).setConvFunc(mospb);
2161 : }
2162 : ///////////////////make sure both FTMachine share the same conv functions.
2163 0 : theIFT= new MosaicFTNew(static_cast<MosaicFTNew &>(*theFT));
2164 :
2165 :
2166 0 : }
2167 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2168 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2169 :
2170 : // Do MS-Selection and set up vi/vb.
2171 : // Only this functions needs to know anything about the MS
2172 0 : void SynthesisImager::createVisSet(const Bool writeAccess)
2173 : {
2174 0 : LogIO os( LogOrigin("SynthesisImager","createVisSet",WHERE) );
2175 :
2176 : ////////////Temporary revert to vi/vb
2177 0 : Block<Int> sort(0);
2178 0 : Block<MeasurementSet> msblock(mss4vi_p.nelements());
2179 : //for (uInt k=0; k< msblock.nelements(); ++k){
2180 : // msblock[k]=*mss_p[k];
2181 : //}
2182 :
2183 : //vs_p= new VisSet(blockMSSel_p, sort, noChanSel, useModelCol_p);
2184 0 : if(!writeAccess){
2185 :
2186 0 : if(rvi_p) delete rvi_p;
2187 0 : rvi_p = NULL;
2188 0 : rvi_p=new ROVisibilityIterator(mss4vi_p, sort);
2189 :
2190 : }
2191 : else{
2192 : // if(wvi_p) delete wvi_p;
2193 : // wvi_p = NULL;
2194 0 : wvi_p=new VisibilityIterator(mss4vi_p, sort);
2195 0 : rvi_p=wvi_p;
2196 : }
2197 0 : Block<Vector<Int> > blockGroup(msblock.nelements());
2198 0 : for (uInt k=0; k < msblock.nelements(); ++k){
2199 0 : blockGroup[k].resize(blockSpw_p[k].nelements());
2200 0 : blockGroup[k].set(1);
2201 : // cerr << "start " << blockStart_p[k] << " nchan " << blockNChan_p[k] << " step " << blockStep_p[k] << " spw "<< blockSpw_p[k] <<endl;
2202 : }
2203 :
2204 0 : rvi_p->selectChannel(blockGroup, blockStart_p, blockNChan_p,
2205 0 : blockStep_p, blockSpw_p);
2206 0 : rvi_p->useImagingWeight(VisImagingWeight("natural"));
2207 : ////////////////////end of revert vi/vb
2208 0 : }// end of createVisSet
2209 :
2210 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2211 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2212 :
2213 :
2214 0 : void SynthesisImager::runMajorCycle(const Bool dopsf,
2215 : const Bool savemodel)
2216 : {
2217 0 : LogIO os( LogOrigin("SynthesisImager","runMajorCycle",WHERE) );
2218 :
2219 : // cout << "Savemodel : " << savemodel << " readonly : " << readOnly_p << " usescratch : " << useScratch_p << endl;
2220 :
2221 0 : Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
2222 0 : Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;
2223 :
2224 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
2225 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
2226 :
2227 0 : SynthesisUtilMethods::getResource("Start Major Cycle");
2228 :
2229 0 : itsMappers.checkOverlappingModels("blank");
2230 :
2231 : {
2232 0 : VisBufferAutoPtr vb(rvi_p);
2233 0 : rvi_p->originChunks();
2234 0 : rvi_p->origin();
2235 :
2236 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2237 0 : dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
2238 0 : Int cohDone=0;
2239 :
2240 :
2241 0 : if(!dopsf)itsMappers.initializeDegrid(*vb);
2242 0 : itsMappers.initializeGrid(*vb,dopsf);
2243 : //SynthesisUtilMethods::getResource("After initGrid for all mappers");
2244 :
2245 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2246 : {
2247 :
2248 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2249 : {
2250 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
2251 : // cerr << "nRows "<< vb->nRow() << " " << max(vb->visCube()) << endl;
2252 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2253 : {
2254 0 : if(!dopsf) {
2255 0 : { vb->setModelVisCube(Complex(0.0, 0.0)); }
2256 0 : itsMappers.degrid(*vb, savevirtualmodel );
2257 0 : if(savemodelcolumn && writeAccess_p )
2258 : {
2259 : // cout << "Vi: " << vb->modelVisCube() << endl;
2260 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2261 : }
2262 : }
2263 0 : itsMappers.grid(*vb, dopsf, datacol_p);
2264 0 : cohDone += vb->nRow();
2265 0 : pm.update(Double(cohDone));
2266 : }
2267 : }
2268 : }
2269 : //cerr << "IN SYNTHE_IMA" << endl;
2270 : //VisModelData::listModel(rvi_p->getMeasurementSet());
2271 0 : SynthesisUtilMethods::getResource("Before finalize for all mappers");
2272 0 : if(!dopsf) itsMappers.finalizeDegrid(*vb);
2273 0 : itsMappers.finalizeGrid(*vb, dopsf);
2274 :
2275 : }
2276 :
2277 0 : itsMappers.checkOverlappingModels("restore");
2278 :
2279 0 : unlockMSs();
2280 :
2281 0 : SynthesisUtilMethods::getResource("End Major Cycle");
2282 :
2283 0 : }// end runMajorCycle
2284 :
2285 :
2286 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2287 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2288 :
2289 : /// The mapper loop is outside the data iterator loop.
2290 : /// This is for cases where the image size is large compared to the RAM and
2291 : /// where data I/O is the relatively minor cost.
2292 0 : void SynthesisImager::runMajorCycle2(const Bool dopsf,
2293 : const Bool savemodel)
2294 : {
2295 0 : LogIO os( LogOrigin("SynthesisImager","runMajorCycle2",WHERE) );
2296 :
2297 : // cout << "Savemodel : " << savemodel << " readonly : " << readOnly_p << " usescratch : " << useScratch_p << endl;
2298 :
2299 0 : Bool savemodelcolumn = savemodel && !readOnly_p && useScratch_p;
2300 0 : Bool savevirtualmodel = savemodel && !readOnly_p && !useScratch_p;
2301 :
2302 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
2303 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
2304 :
2305 0 : itsMappers.checkOverlappingModels("blank");
2306 :
2307 0 : Bool resetModel=False;
2308 0 : if( savemodelcolumn && writeAccess_p)
2309 : {
2310 0 : resetModel=True;
2311 0 : os << "Iterating through the model column to reset it to zero" << LogIO::POST;
2312 0 : VisBufferAutoPtr vb(rvi_p);
2313 0 : rvi_p->originChunks();
2314 0 : rvi_p->origin();
2315 :
2316 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2317 0 : dopsf?"Seting model column to zero":"pre-Major Cycle", "","","",True);
2318 0 : Int cohDone=0;
2319 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2320 : {
2321 :
2322 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2323 : {
2324 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2325 : {
2326 0 : vb->setModelVisCube(Complex(0.0, 0.0));
2327 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2328 : }
2329 0 : cohDone += vb->nRow();
2330 0 : pm.update(Double(cohDone));
2331 : }
2332 : }
2333 : }// setting model to zero
2334 :
2335 0 : for(Int gmap=0;gmap<itsMappers.nMappers();gmap++)
2336 : {
2337 0 : os << "Running major cycle for chunk : " << gmap << LogIO::POST;
2338 :
2339 0 : SynthesisUtilMethods::getResource("Start Major Cycle for mapper"+String::toString(gmap));
2340 0 : VisBufferAutoPtr vb(rvi_p);
2341 0 : rvi_p->originChunks();
2342 0 : rvi_p->origin();
2343 :
2344 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2345 0 : dopsf?"Gridding Weights and PSF":"Major Cycle", "","","",true);
2346 0 : Int cohDone=0;
2347 :
2348 0 : if(!dopsf){
2349 0 : itsMappers.initializeDegrid(*vb, gmap);
2350 : //itsMappers.getMapper(gmap)->initializeDegrid(*vb);
2351 : }
2352 0 : itsMappers.initializeGrid(*vb,dopsf, gmap);
2353 : //itsMappers.getMapper(gmap)->initializeGrid(*vb,dopsf);
2354 :
2355 0 : SynthesisUtilMethods::getResource("After initialize for mapper"+String::toString(gmap));
2356 :
2357 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2358 : {
2359 :
2360 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2361 : {
2362 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; // No valid rows in this VB
2363 : // cerr << "nRows "<< vb->nRow() << " " << max(vb->visCube()) << endl;
2364 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2365 : {
2366 0 : if(!dopsf) {
2367 0 : if(resetModel==False) { vb->setModelVisCube(Complex(0.0, 0.0)); }
2368 0 : itsMappers.degrid(*vb, savevirtualmodel, gmap );
2369 : //itsMappers.getMapper(gmap)->degrid(*vb); //, savevirtualmodel );
2370 0 : if(savemodelcolumn && writeAccess_p )
2371 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2372 : }
2373 0 : itsMappers.grid(*vb, dopsf, datacol_p, gmap);
2374 : //itsMappers.getMapper(gmap)->grid(*vb, dopsf, datacol_p);
2375 0 : cohDone += vb->nRow();
2376 0 : pm.update(Double(cohDone));
2377 : }
2378 : }
2379 : }
2380 : //cerr << "IN SYNTHE_IMA" << endl;
2381 : //VisModelData::listModel(rvi_p->getMeasurementSet());
2382 :
2383 0 : SynthesisUtilMethods::getResource("Before finalize for mapper"+String::toString(gmap));
2384 :
2385 0 : if(!dopsf)
2386 : {
2387 0 : itsMappers.finalizeDegrid(*vb,gmap);
2388 : // auto back = itsMappers.getMapper(gmap)->imageStore()->forwardGrid();
2389 : // back->tempClose();
2390 : }
2391 0 : itsMappers.finalizeGrid(*vb, dopsf,gmap);
2392 :
2393 : //itsMappers.getMapper(gmap)->imageStore()->releaseLocks();
2394 0 : itsMappers.getMapper(gmap)->imageStore()->releaseComplexGrids();
2395 :
2396 :
2397 : // auto back = itsMappers.getMapper(gmap)->imageStore()->backwardGrid();
2398 : // back->tempClose();
2399 :
2400 0 : SynthesisUtilMethods::getResource("End Major Cycle for mapper"+String::toString(gmap));
2401 : }// end of mapper loop
2402 :
2403 0 : itsMappers.checkOverlappingModels("restore");
2404 :
2405 0 : unlockMSs();
2406 :
2407 0 : SynthesisUtilMethods::getResource("End Major Cycle");
2408 :
2409 0 : }// end runMajorCycle2
2410 :
2411 :
2412 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2413 0 : void SynthesisImager::predictModel(){
2414 0 : LogIO os( LogOrigin("SynthesisImager","predictModel ",WHERE) );
2415 :
2416 0 : os << "---------------------------------------------------- Predict Model ---------------------------------------------" << LogIO::POST;
2417 :
2418 0 : Bool savemodelcolumn = !readOnly_p && useScratch_p;
2419 0 : Bool savevirtualmodel = !readOnly_p && !useScratch_p;
2420 :
2421 0 : if( savemodelcolumn ) os << "Saving model column" << LogIO::POST;
2422 0 : if( savevirtualmodel ) os << "Saving virtual model" << LogIO::POST;
2423 :
2424 0 : itsMappers.checkOverlappingModels("blank");
2425 :
2426 :
2427 : {
2428 0 : VisBufferAutoPtr vb(rvi_p);
2429 0 : rvi_p->originChunks();
2430 0 : rvi_p->origin();
2431 :
2432 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2433 0 : "Predict Model", "","","",true);
2434 0 : Int cohDone=0;
2435 :
2436 0 : itsMappers.initializeDegrid(*vb);
2437 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2438 : {
2439 :
2440 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2441 : {
2442 : //if (SynthesisUtilMethods::validate(*vb)==SynthesisUtilMethods::NOVALIDROWS) break; //No valid rows in this MS
2443 : //if !usescratch ...just save
2444 0 : vb->setModelVisCube(Complex(0.0, 0.0));
2445 0 : itsMappers.degrid(*vb, savevirtualmodel);
2446 0 : if(savemodelcolumn && writeAccess_p )
2447 0 : wvi_p->setVis(vb->modelVisCube(),VisibilityIterator::Model);
2448 :
2449 : // cout << "nRows "<< vb->nRow() << " " << max(vb->modelVisCube()) << endl;
2450 0 : cohDone += vb->nRow();
2451 0 : pm.update(Double(cohDone));
2452 :
2453 : }
2454 : }
2455 0 : itsMappers.finalizeDegrid(*vb);
2456 : }
2457 :
2458 0 : itsMappers.checkOverlappingModels("restore");
2459 0 : unlockMSs();
2460 :
2461 0 : }// end of predictModel
2462 :
2463 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2464 0 : void SynthesisImager::makeSdImage(Bool dopsf)
2465 : {
2466 0 : LogIO os( LogOrigin("SynthesisImager","makeSdImage",WHERE) );
2467 :
2468 : // Bool dopsf=false;
2469 0 : if(datacol_p==FTMachine::PSF) dopsf=true;
2470 :
2471 : {
2472 0 : VisBufferAutoPtr vb(rvi_p);
2473 0 : rvi_p->originChunks();
2474 0 : rvi_p->origin();
2475 :
2476 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()),
2477 0 : String(datacol_p), "","","",true);
2478 0 : Int cohDone=0;
2479 :
2480 0 : itsMappers.initializeGrid(*vb,dopsf);
2481 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2482 : {
2483 :
2484 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2485 : {
2486 0 : itsMappers.grid(*vb, dopsf, datacol_p);
2487 0 : cohDone += vb->nRow();
2488 0 : pm.update(Double(cohDone));
2489 : }
2490 : }
2491 0 : itsMappers.finalizeGrid(*vb, dopsf);
2492 :
2493 : }
2494 :
2495 0 : unlockMSs();
2496 :
2497 0 : }// end makeSdImage
2498 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2499 :
2500 0 : void SynthesisImager::makeImage(String type, const String& imagename, const String& complexImage, const Int whichModel)
2501 : {
2502 0 : LogIO os( LogOrigin("SynthesisImager","makeImage",WHERE) );
2503 :
2504 :
2505 0 : FTMachine::Type seType(FTMachine::OBSERVED);
2506 0 : if(type=="observed") {
2507 0 : seType=FTMachine::OBSERVED;
2508 : os << LogIO::NORMAL // Loglevel INFO
2509 : << "Making dirty image from " << type << " data "
2510 0 : << LogIO::POST;
2511 : }
2512 0 : else if (type=="model") {
2513 0 : seType=FTMachine::MODEL;
2514 : os << LogIO::NORMAL // Loglevel INFO
2515 : << "Making dirty image from " << type << " data "
2516 0 : << LogIO::POST;
2517 : }
2518 0 : else if (type=="corrected") {
2519 0 : seType=FTMachine::CORRECTED;
2520 : os << LogIO::NORMAL // Loglevel INFO
2521 : << "Making dirty image from " << type << " data "
2522 0 : << LogIO::POST;
2523 : }
2524 0 : else if (type=="psf") {
2525 0 : seType=FTMachine::PSF;
2526 : os << "Making point spread function "
2527 0 : << LogIO::POST;
2528 : }
2529 0 : else if (type=="residual") {
2530 0 : seType=FTMachine::RESIDUAL;
2531 : os << LogIO::NORMAL // Loglevel INFO
2532 : << "Making dirty image from " << type << " data "
2533 0 : << LogIO::POST;
2534 : }
2535 0 : else if (type=="singledish-observed") {
2536 0 : seType=FTMachine::OBSERVED;
2537 : os << LogIO::NORMAL
2538 0 : << "Making single dish image from observed data" << LogIO::POST;
2539 : }
2540 0 : else if (type=="singledish") {
2541 0 : seType=FTMachine::CORRECTED;
2542 : os << LogIO::NORMAL
2543 0 : << "Making single dish image from corrected data" << LogIO::POST;
2544 : }
2545 0 : else if (type=="coverage") {
2546 0 : seType=FTMachine::COVERAGE;
2547 : os << LogIO::NORMAL
2548 : << "Making single dish coverage function "
2549 0 : << LogIO::POST;
2550 : }
2551 0 : else if (type=="holography") {
2552 0 : seType=FTMachine::CORRECTED;
2553 : os << LogIO::NORMAL
2554 : << "Making complex holographic image from corrected data "
2555 0 : << LogIO::POST;
2556 : }
2557 0 : else if (type=="holography-observed") {
2558 0 : seType=FTMachine::OBSERVED;
2559 : os << LogIO::NORMAL
2560 : << "Making complex holographic image from observed data "
2561 0 : << LogIO::POST;
2562 : }
2563 :
2564 :
2565 0 : String imageName=(itsMappers.imageStore(whichModel))->getName();
2566 0 : String cImageName(complexImage);
2567 0 : Bool keepComplexImage=(complexImage!="")||(type.contains("holography"));
2568 0 : if(complexImage=="") {
2569 0 : cImageName=imageName+".compleximage";
2570 : }
2571 0 : PagedImage<Float> theImage((itsMappers.imageStore(whichModel))->getShape(), (itsMappers.imageStore(whichModel))->getCSys(), imagename);
2572 0 : PagedImage<Complex> cImageImage(theImage.shape(),
2573 : theImage.coordinates(),
2574 0 : cImageName);
2575 0 : if(!keepComplexImage)
2576 0 : cImageImage.table().markForDelete();
2577 :
2578 :
2579 0 : Matrix<Float> weight;
2580 0 : (itsMappers.getFTM(whichModel))->makeImage(seType, *rvi_p, cImageImage, weight);
2581 :
2582 0 : if(seType==FTMachine::PSF){
2583 0 : StokesImageUtil::ToStokesPSF(theImage, cImageImage);
2584 0 : StokesImageUtil::normalizePSF(theImage);
2585 : }
2586 : else{
2587 0 : StokesImageUtil::To(theImage, cImageImage);
2588 : }
2589 :
2590 :
2591 0 : unlockMSs();
2592 :
2593 0 : }// end makeImage
2594 :
2595 :
2596 :
2597 :
2598 :
2599 :
2600 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2601 : // Method to run the AWProjectFT machine to seperate the CFCache
2602 : // construction from imaging. This is done by splitting the
2603 : // operation in two steps: (1) run the FTM in "dry" mode to create a
2604 : // "blank" CFCache, and (2) re-load the "blank" CFCache and "fill"
2605 : // it.
2606 : //
2607 : // If someone can get me (SB) out of the horrible statc_casts in the
2608 : // code below, I will be most grateful (we are out of it! :-)).
2609 : //
2610 0 : void SynthesisImager::dryGridding(const Vector<String>& cfList)
2611 : {
2612 0 : LogIO os( LogOrigin("SynthesisImager","dryGridding",WHERE) );
2613 0 : Int cohDone=0, whichFTM=0;
2614 : (void)cfList;
2615 : // If not an AWProject-class FTM, make this call a NoOp. Might be
2616 : // useful to extend it to other projection FTMs -- but later.
2617 0 : String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();
2618 :
2619 0 : if (!((itsMappers.getFTM(whichFTM,true))->isUsingCFCache())) return;
2620 :
2621 0 : os << "---------------------------------------------------- Dry Gridding ---------------------------------------------" << LogIO::POST;
2622 :
2623 : //
2624 : // Go through the entire MS in "dry" mode to set up a "blank"
2625 : // CFCache. This is done by setting the AWPWBFT in dryrun mode
2626 : // and gridding. The process of gridding emits CFCache, which
2627 : // will be "blank" in a dry run.
2628 : {
2629 0 : VisBufferAutoPtr vb(rvi_p);
2630 0 : rvi_p->originChunks();
2631 0 : rvi_p->origin();
2632 :
2633 0 : ProgressMeter pm(1.0, Double(vb->numberCoh()), "dryGridding", "","","",true);
2634 :
2635 0 : itsMappers.initializeGrid(*vb);
2636 :
2637 : // Set the gridder (iFTM) to run in dry-gridding mode
2638 0 : (itsMappers.getFTM(whichFTM,true))->setDryRun(true);
2639 :
2640 0 : Bool aTermIsOff=False;
2641 : {
2642 0 : CountedPtr<FTMachine> ftm=itsMappers.getFTM(whichFTM,True);
2643 0 : CountedPtr<ConvolutionFunction> cf=ftm->getAWConvFunc();
2644 : // Bool aTermIsOff = (((static_cast<AWConvFunc &> (*cf)).getTerm("ATerm")))->isNoOp();
2645 0 : aTermIsOff = cf->getTerm("ATerm")->isNoOp();
2646 : }
2647 :
2648 : os << "Making a \"blank\" CFCache"
2649 : << (aTermIsOff?" (without the A-Term)":"")
2650 0 : << LogIO::WARN << LogIO::POST;
2651 :
2652 : // Step through the MS. This triggers the logic in the Gridder
2653 : // to determine all the CFs that will be required. These empty
2654 : // CFs are written to the CFCache which can then be filled via
2655 : // a call to fillCFCache().
2656 0 : for (rvi_p->originChunks(); rvi_p->moreChunks();rvi_p->nextChunk())
2657 : {
2658 0 : for (rvi_p->origin(); rvi_p->more(); (*rvi_p)++)
2659 : {
2660 0 : if (SynthesisUtilMethods::validate(*vb)!=SynthesisUtilMethods::NOVALIDROWS)
2661 : {
2662 0 : itsMappers.grid(*vb, true, FTMachine::OBSERVED, whichFTM);
2663 0 : cohDone += vb->nRow();
2664 0 : pm.update(Double(cohDone));
2665 : // If there is no term that depends on time, don't iterate over the entire data base
2666 0 : if (aTermIsOff) break;
2667 : }
2668 : }
2669 : }
2670 : }
2671 0 : if (cohDone == 0) os << "No valid rows found in dryGridding." << LogIO::EXCEPTION << LogIO::POST;
2672 : // Unset the dry-gridding mode.
2673 0 : (itsMappers.getFTM(whichFTM,true))->setDryRun(false);
2674 :
2675 : //itsMappers.checkOverlappingModels("restore");
2676 0 : unlockMSs();
2677 : //fillCFCache(cfList);
2678 : }
2679 : //
2680 : // Re-load the CFCache from the disk using the supplied list of CFs
2681 : // (as cfList). Then extract the ConvFunc object (which was setup
2682 : // in the FTM) and call it's makeConvFunction2() to fill the CF.
2683 : // Finally, unset the dry-run mode of the FTM.
2684 : //
2685 0 : void SynthesisImager::fillCFCache(const Vector<String>& cfList,
2686 : const String& ftmName,
2687 : const String& cfcPath,
2688 : const Bool& psTermOn,
2689 : const Bool& aTermOn,
2690 : const Bool& conjBeams)
2691 : {
2692 0 : LogIO os( LogOrigin("SynthesisImager","fillCFCache",WHERE) );
2693 : // If not an AWProject-class FTM, make this call a NoOp. Might be
2694 : // useful to extend it to other projection FTMs -- but later.
2695 : // String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();
2696 :
2697 0 : if (!ftmName.contains("awproject") and
2698 0 : !ftmName.contains("multitermftnew")) return;
2699 :
2700 0 : os << "---------------------------------------------------- fillCFCache ---------------------------------------------" << LogIO::POST;
2701 :
2702 : //String cfcPath = itsMappers.getFTM(whichFTM)->getCacheDir();
2703 : //String imageNamePrefix=itsMappers.getFTM(whichFTM)->getCFCache()->getWtImagePrefix();
2704 :
2705 : //cerr << "Path = " << path << endl;
2706 :
2707 : // CountedPtr<AWProjectWBFTNew> tmpFT = new AWProjectWBFTNew(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM))));
2708 :
2709 :
2710 0 : Float dPA=360.0,selectedPA=2*360.0;
2711 0 : if (cfList.nelements() > 0)
2712 : {
2713 0 : CountedPtr<CFCache> cfCacheObj = new CFCache();
2714 : //Vector<String> wtCFList; wtCFList.resize(cfList.nelements());
2715 : //for (Int i=0; i<wtCFList.nelements(); i++) wtCFList[i] = "WT"+cfList[i];
2716 : //Directory dir(path);
2717 0 : Vector<String> cfList_p=cfList;//dir.find(Regex(Regex::fromPattern("CFS*")));
2718 0 : Vector<String> wtCFList_p;
2719 0 : wtCFList_p.resize(cfList_p.nelements());
2720 0 : for (Int i=0; i<(Int)wtCFList_p.nelements(); i++) wtCFList_p[i]="WT"+cfList_p[i];
2721 :
2722 : //cerr << cfList_p << endl;
2723 0 : cfCacheObj->setCacheDir(cfcPath.data());
2724 :
2725 0 : os << "Re-loading the \"blank\" CFCache for filling" << LogIO::WARN << LogIO::POST;
2726 :
2727 0 : cfCacheObj->initCacheFromList2(cfcPath, cfList_p, wtCFList_p,
2728 : selectedPA, dPA,1);
2729 : // tmpFT->setCFCache(cfCacheObj);
2730 0 : Vector<Double> uvScale, uvOffset;
2731 0 : Matrix<Double> vbFreqSelection;
2732 0 : CountedPtr<CFStore2> cfs2 = CountedPtr<CFStore2>(&cfCacheObj->memCache2_p[0],false);//new CFStore2;
2733 0 : CountedPtr<CFStore2> cfwts2 = CountedPtr<CFStore2>(&cfCacheObj->memCacheWt2_p[0],false);//new CFStore2;
2734 :
2735 : //
2736 : // Get whichFTM from itsMappers (SIMapperCollection) and
2737 : // cast it as AWProjectWBFTNew. Then get the ConvFunc from
2738 : // the FTM and cast it as AWConvFunc. Finally call
2739 : // AWConvFunc::makeConvFunction2().
2740 : //
2741 : // (static_cast<AWConvFunc &>
2742 : // (*(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getAWConvFunc())
2743 : // ).makeConvFunction2(String(path), uvScale, uvOffset, vbFreqSelection,
2744 : // *cfs2, *cfwts2);
2745 :
2746 : // This is a global methond in AWConvFunc. Does not require
2747 : // FTM to be constructed (which can be expensive in terms of
2748 : // memory footprint).
2749 0 : AWConvFunc::makeConvFunction2(String(cfcPath), uvScale, uvOffset, vbFreqSelection,
2750 0 : *cfs2, *cfwts2, psTermOn, aTermOn, conjBeams);
2751 : }
2752 : //cerr << "Mem used = " << itsMappers.getFTM(whichFTM)->getCFCache()->memCache2_p[0].memUsage() << endl;
2753 : //(static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).getCFCache()->initCache2();
2754 : }
2755 :
2756 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2757 0 : void SynthesisImager::reloadCFCache()
2758 : {
2759 0 : LogIO os( LogOrigin("SynthesisImager","reloadCFCache",WHERE) );
2760 0 : Int whichFTM=0;
2761 0 : String ftmName = ((*(itsMappers.getFTM(whichFTM)))).name();
2762 0 : if (!ftmName.contains("AWProject")) return;
2763 :
2764 0 : os << "-------------------------------------------- reloadCFCache ---------------------------------------------" << LogIO::POST;
2765 0 : String path = itsMappers.getFTM(whichFTM)->getCacheDir();
2766 0 : String imageNamePrefix=itsMappers.getFTM(whichFTM)->getCFCache()->getWtImagePrefix();
2767 :
2768 0 : CountedPtr<CFCache> cfCacheObj = new CFCache();
2769 0 : cfCacheObj->setCacheDir(path.data());
2770 0 : cfCacheObj->setWtImagePrefix(imageNamePrefix.c_str());
2771 0 : cfCacheObj->initCache2();
2772 :
2773 : // This assumes the itsMappers is always SIMapperCollection.
2774 0 : for (whichFTM = 0; whichFTM < itsMappers.nMappers(); whichFTM++)
2775 : {
2776 0 : (static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM)))).setCFCache(cfCacheObj,true); // Setup iFTM
2777 0 : (static_cast<AWProjectWBFTNew &> (*(itsMappers.getFTM(whichFTM,false)))).setCFCache(cfCacheObj,true); // Set FTM
2778 : }
2779 : }
2780 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2781 :
2782 : //Utility function to properly convert Double to String.
2783 : //With C++11 we can probably use STL to_string() function instead...
2784 0 : String SynthesisImager::doubleToString(const Double& df) {
2785 0 : std::ostringstream ss;
2786 0 : ss.precision(std::numeric_limits<double>::digits10+2);
2787 0 : ss << df;
2788 : //return ss.str();
2789 0 : return to_string(df);
2790 : }
2791 :
2792 : /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
2793 :
2794 0 : Bool SynthesisImager::makePB()
2795 : {
2796 0 : LogIO os( LogOrigin("SynthesisImager","makePB",WHERE) );
2797 :
2798 0 : if( itsMakeVP==False )
2799 : {
2800 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;
2801 : // Check that the .weight exists.. ?
2802 :
2803 0 : return False;
2804 : }
2805 : else
2806 : {
2807 0 : Bool doDefaultVP = itsVpTable.length()>0 ? False : True;
2808 :
2809 0 : CoordinateSystem coordsys=itsMappers.imageStore(0)->getCSys();
2810 0 : String telescope=coordsys.obsInfo().telescope();
2811 :
2812 0 : if (doDefaultVP) {
2813 :
2814 0 : MSAntennaColumns ac(mss4vi_p[0].antenna());
2815 0 : Double dishDiam=ac.dishDiameter()(0);
2816 0 : if(!allEQ(ac.dishDiameter().getColumn(), dishDiam))
2817 : os << LogIO::WARN
2818 : << "The MS has multiple antenna diameters ..PB could be wrong "
2819 0 : << LogIO::POST;
2820 0 : return makePBImage( telescope, False, dishDiam);
2821 : }
2822 : else{
2823 0 : return makePBImage(telescope );
2824 : }
2825 :
2826 : }
2827 : return False;
2828 : }
2829 :
2830 0 : bool SynthesisImager::makePBImage(const String& telescopeName,
2831 : Bool useSymmetricBeam, Double diam){
2832 :
2833 0 : LogIO os(LogOrigin("SynthesisImager", "makePBImage()", WHERE));
2834 :
2835 : // Check if this metadata info should come from each Mapper or if it's OK to be just the first one.
2836 : // Right now it assumes all mappers have the same FREQUENCY settings.
2837 :
2838 0 : CoordinateSystem imageCoord=itsMappers.imageStore(0)->getCSys();
2839 :
2840 0 : Int spectralIndex=imageCoord.findCoordinate(Coordinate::SPECTRAL);
2841 : SpectralCoordinate
2842 0 : spectralCoord=imageCoord.spectralCoordinate(spectralIndex);
2843 0 : Vector<String> units(1); units = "Hz";
2844 0 : spectralCoord.setWorldAxisUnits(units);
2845 0 : Vector<Double> spectralWorld(1);
2846 0 : Vector<Double> spectralPixel(1);
2847 0 : spectralPixel(0) = 0;
2848 0 : spectralCoord.toWorld(spectralWorld, spectralPixel);
2849 0 : Double freq = spectralWorld(0);
2850 0 : Quantity qFreq( freq, "Hz" );
2851 0 : String telName=telescopeName;
2852 0 : if(telName=="ALMA" && diam < 12.0)
2853 0 : telName="ACA";
2854 : //cerr << "Telescope Name is " << telName<< endl;
2855 : PBMath::CommonPB whichPB;
2856 0 : PBMath::enumerateCommonPB(telName, whichPB);
2857 0 : PBMath myPB;
2858 0 : if(whichPB!=PBMath::UNKNOWN && whichPB!=PBMath::NONE){
2859 :
2860 0 : myPB=PBMath(telName, useSymmetricBeam, qFreq);
2861 : }
2862 0 : else if(diam > 0.0){
2863 0 : myPB=PBMath(diam,useSymmetricBeam, qFreq);
2864 : }
2865 : else{
2866 : os << LogIO::WARN << "Telescope " << telName << " is not known\n "
2867 : << "Not making the PB image"
2868 0 : << LogIO::POST;
2869 0 : return false;
2870 : }
2871 0 : return makePrimaryBeam(myPB );
2872 : }
2873 :
2874 0 : Bool SynthesisImager::makePBImage(const String telescop){
2875 :
2876 : /*
2877 : ScalarColumn<TableRecord> recCol(vpTable, (String)"pbdescription");
2878 : PBMath myPB(recCol(0));
2879 : */
2880 0 : LogIO os( LogOrigin("SynthesisImager","makePBImage",WHERE) );
2881 :
2882 : PBMath::CommonPB kpb;
2883 0 : Record rec;
2884 0 : getVPRecord( rec, kpb, telescop );
2885 :
2886 0 : Bool ret=False;
2887 0 : if(rec.empty())
2888 0 : {os << LogIO::WARN << "Not making PB. Cannot use pbcor option later." << LogIO::POST; }
2889 : else
2890 : {
2891 0 : PBMath myPB( rec );
2892 0 : ret = makePrimaryBeam(myPB);
2893 : }
2894 0 : return ret;
2895 : }
2896 :
2897 :
2898 0 : Bool SynthesisImager::makePrimaryBeam(PBMath& pbMath)
2899 : {
2900 0 : LogIO os( LogOrigin("SynthesisImager","makePrimaryBeam",WHERE) );
2901 :
2902 0 : os << "Evaluating Primary Beam model onto image grid(s)" << LogIO::POST;
2903 :
2904 0 : itsMappers.initPB();
2905 :
2906 0 : VisBuffer vb(*rvi_p);
2907 0 : Int fieldCounter=0;
2908 0 : Vector<Int> fieldsDone;
2909 :
2910 0 : for(rvi_p->originChunks(); rvi_p->moreChunks(); rvi_p->nextChunk()){
2911 0 : Bool fieldDone=false;
2912 0 : for (uInt k=0; k < fieldsDone.nelements(); ++k){
2913 0 : fieldDone=fieldDone || (vb.fieldId()==fieldsDone(k));
2914 : }
2915 0 : if(!fieldDone){
2916 0 : ++fieldCounter;
2917 0 : fieldsDone.resize(fieldCounter, true);
2918 0 : fieldsDone(fieldCounter-1)=vb.fieldId();
2919 :
2920 0 : itsMappers.addPB(vb,pbMath);
2921 :
2922 : }
2923 : }
2924 :
2925 0 : itsMappers.releaseImageLocks();
2926 0 : unlockMSs();
2927 :
2928 0 : return True;
2929 : }// end makePB
2930 :
2931 : /////===========
2932 0 : void SynthesisImager::setMovingSource(const String& movingSource){
2933 0 : movingSource_p=movingSource;
2934 0 : }
2935 :
2936 0 : bool SynthesisImager::isSpectralCube(){
2937 0 : bool retval=False;
2938 0 : for (Int k=0; k < itsMappers.nMappers(); ++k){
2939 : //For some reason imstore sometime returns 0 shape
2940 : //trying to test with with psf size breaks parallel = true in some cases...go figure
2941 : //if((((itsMappers.imageStore(k))->psf())->shape()[3]) != ((itsMappers.imageStore(k))->getShape()[3])){
2942 : //cerr << "shapes " << ((itsMappers.imageStore(k))->psf())->shape() << " " << ((itsMappers.imageStore(k))->getShape()) << endl;
2943 : //throw(AipsError("images shape seem insistent "));
2944 0 : if((itsMappers.imageStore(k))->getShape()(3) ==0)
2945 0 : return True;
2946 : //}
2947 0 : if((itsMappers.imageStore(k))->getShape()(3) > 1)
2948 0 : retval=True;
2949 : }
2950 0 : return retval;
2951 :
2952 : }
2953 0 : void SynthesisImager::normalizerinfo(const Record& normpars){
2954 0 : normpars_p=normpars;
2955 0 : }
2956 : } //# NAMESPACE CASA - END
2957 :
|
