Line data Source code
1 : //# ImagerMultiMS.cc: Implementation of ImagerMultiMS.h
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : #include <cassert>
29 :
30 : #include <synthesis/MeasurementEquations/ImagerMultiMS.h>
31 : #include <synthesis/TransformMachines/VisModelData.h>
32 : #include <casacore/casa/BasicSL/String.h>
33 : #include <casacore/casa/Utilities/Assert.h>
34 : #include <casacore/casa/Logging.h>
35 : #include <casacore/casa/Logging/LogIO.h>
36 : #include <casacore/casa/Logging/LogMessage.h>
37 :
38 : #include <casacore/casa/Exceptions/Error.h>
39 : #include <iostream>
40 : #include <casacore/ms/MSSel/MSSelection.h>
41 : #include <casacore/ms/MSSel/MSDataDescIndex.h>
42 : #include <casacore/ms/MeasurementSets/MSHistoryHandler.h>
43 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
44 : #include <casacore/ms/MeasurementSets/MSDataDescColumns.h>
45 : #include <casacore/ms/MeasurementSets/MSColumns.h>
46 : #include <msvis/MSVis/SimpleSubMS.h>
47 : #include <msvis/MSVis/SubMS.h>
48 : #include <msvis/MSVis/VisSetUtil.h>
49 : #include <msvis/MSVis/VisibilityIterator.h>
50 : #include <casacore/casa/Arrays/Matrix.h>
51 : #include <casacore/casa/Arrays/ArrayMath.h>
52 :
53 : #include <casacore/tables/TaQL/ExprNode.h>
54 : #include <casacore/tables/TaQL/TableParse.h>
55 : #include <casacore/tables/Tables/SetupNewTab.h>
56 :
57 : #include <casacore/lattices/LEL/LatticeExpr.h>
58 :
59 : #include <casacore/casa/OS/File.h>
60 : #include <casacore/casa/OS/HostInfo.h>
61 : #include <casacore/casa/Containers/Record.h>
62 :
63 : #include <sstream>
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 :
68 1178 : ImagerMultiMS::ImagerMultiMS()
69 1178 : : Imager(), blockNChan_p(0), blockStart_p(0), blockStep_p(0), blockSpw_p(0),
70 2356 : blockMSSel_p(0), dataSet_p(false)
71 : {
72 :
73 1178 : lockCounter_p=0;
74 1178 : ms_p=0;
75 1178 : mssel_p=0;
76 1178 : se_p=0;
77 1178 : vs_p=0;
78 1178 : ft_p=0;
79 1178 : cft_p=0;
80 1178 : rvi_p=wvi_p=0;
81 1178 : numMS_p=0;
82 :
83 1178 : }
84 :
85 0 : Bool ImagerMultiMS::setDataToMemory(const String& msname,
86 : const String& /*mode*/,
87 : const Vector<Int>& nchan,
88 : const Vector<Int>& start,
89 : const Vector<Int>& step,
90 : const Vector<Int>& /*spectralwindowids*/,
91 : const Vector<Int>& /*fieldids*/,
92 : const String& msSelect,
93 : const String& timerng,
94 : const String& fieldnames,
95 : const Vector<Int>& /*antIndex*/,
96 : const String& antnames,
97 : const String& spwstring,
98 : const String& uvdist,
99 : const String& scan,
100 : const String& intent,
101 : const String& obs){
102 0 : useModelCol_p=false;
103 0 : LogIO os(LogOrigin("imager", "setDataToMemory()"), logSink_p);
104 0 : if(!Table::isReadable(msname)){
105 : os << LogIO::SEVERE << "MeasurementSet "
106 : << msname << " does not exist "
107 0 : << LogIO::POST;
108 0 : return false;
109 : }
110 0 : if(intent != ""){
111 : os << LogIO::WARN << "does not support INTENT selection here "
112 0 : << LogIO::POST;
113 : }
114 0 : MeasurementSet thisms(msname, TableLock(TableLock::AutoNoReadLocking),
115 0 : Table::Old);
116 0 : SimpleSubMS splitter(thisms);
117 : //SubMS splitter(thisms);
118 0 : splitter.setmsselect(spwstring, fieldnames, antnames, scan, obs, uvdist,
119 : msSelect, nchan, start, step, "");
120 0 : splitter.selectCorrelations("");
121 0 : splitter.selectTime(-1.0, timerng);
122 0 : MS::PredefinedColumns whichCol=MS::DATA;
123 0 : if(thisms.tableDesc().isColumn("CORRECTED_DATA"))
124 0 : whichCol=MS::CORRECTED_DATA;
125 0 : CountedPtr<MeasurementSet> subMS(splitter.makeMemSubMS(whichCol), true);
126 : //CountedPtr<MeasurementSet> subMS(splitter.makeScratchSubMS(Vector<MS::PredefinedColumns>(1,whichCol), true), true);
127 0 : return setDataOnThisMS(*subMS);
128 :
129 0 : }
130 :
131 :
132 59 : Bool ImagerMultiMS::setDataPerMS(const String& msname, const String& mode,
133 : const Vector<Int>& nchan,
134 : const Vector<Int>& start,
135 : const Vector<Int>& step,
136 : const Vector<Int>& spectralwindowids,
137 : const Vector<Int>& fieldids,
138 : const String& msSelect,
139 : const String& timerng,
140 : const String& fieldnames,
141 : const Vector<Int>& antIndex,
142 : const String& antnames,
143 : const String& spwstring,
144 : const String& uvdist,
145 : const String& scan,
146 : const String& intent,
147 : const String& obs,
148 : const Bool useModel,
149 : const Bool readonly){
150 59 : useModelCol_p=useModel;
151 : // Bool rd=readonly;
152 118 : LogIO os(LogOrigin("imager", "setDataPerMS()"), logSink_p);
153 : // if(useModel)
154 : // rd=true;
155 59 : if(!Table::isReadable(msname)){
156 : os << LogIO::SEVERE << "MeasurementSet "
157 : << msname << " does not exist "
158 0 : << LogIO::POST;
159 0 : return false;
160 : }
161 : else{
162 59 : MeasurementSet thisms;
163 59 : if(!readonly)
164 118 : thisms=MeasurementSet(msname, TableLock(TableLock::AutoNoReadLocking),
165 59 : Table::Update);
166 : else
167 0 : thisms=MeasurementSet(msname, TableLock(TableLock::AutoNoReadLocking),
168 0 : Table::Old);
169 : //thisms=MeasurementSet(msname, TableLock(),
170 : // Table::Old);
171 :
172 59 : thisms.setMemoryResidentSubtables (MrsEligibility::defaultEligible());
173 :
174 59 : return setDataOnThisMS(thisms, mode, nchan, start, step, spectralwindowids, fieldids, msSelect, timerng,
175 59 : fieldnames, antIndex, antnames, spwstring, uvdist, scan, intent, obs);
176 :
177 59 : }
178 :
179 59 : }
180 :
181 :
182 :
183 59 : Bool ImagerMultiMS::setDataOnThisMS(MeasurementSet& thisms, const String& mode,
184 : const Vector<Int>& nchan,
185 : const Vector<Int>& start,
186 : const Vector<Int>& step,
187 : const Vector<Int>& spectralwindowids,
188 : const Vector<Int>& fieldids,
189 : const String& msSelect,
190 : const String& timerng,
191 : const String& fieldnames,
192 : const Vector<Int>& antIndex,
193 : const String& antnames,
194 : const String& spwstring,
195 : const String& uvdist,
196 : const String& scan,
197 : const String& intent,
198 : const String& obs){
199 118 : LogIO os(LogOrigin("imager", "setDataOnThisMS()"), logSink_p);
200 :
201 :
202 59 : dataMode_p=mode;
203 59 : dataNchan_p.resize();
204 59 : dataStart_p.resize();
205 59 : dataStep_p.resize();
206 59 : dataNchan_p=nchan;
207 59 : dataStart_p=start;
208 59 : dataStep_p=step;
209 59 : dataspectralwindowids_p.resize(spectralwindowids.nelements());
210 59 : dataspectralwindowids_p=spectralwindowids;
211 59 : datafieldids_p.resize(fieldids.nelements());
212 59 : datafieldids_p=fieldids;
213 :
214 :
215 :
216 : //make ms
217 : // auto lock for now
218 : // make ms and visset
219 59 : ++numMS_p;
220 59 : blockMSSel_p.resize(numMS_p);
221 59 : blockNChan_p.resize(numMS_p);
222 59 : blockStart_p.resize(numMS_p);
223 59 : blockStep_p.resize(numMS_p);
224 59 : blockSpw_p.resize(numMS_p);
225 : //Using autolocking here
226 : //Will need to rethink this when in parallel mode with multi-cpu access
227 59 : blockMSSel_p[numMS_p-1]=thisms;
228 : //breaking reference
229 59 : if(ms_p.null())
230 35 : ms_p=new MeasurementSet();
231 : else
232 24 : (*ms_p)=MeasurementSet();
233 59 : (*ms_p)=thisms;
234 :
235 :
236 : try{
237 59 : openSubTables();
238 :
239 : //if spectralwindowids=-1 take all
240 118 : if(mode=="none" &&
241 59 : (spectralwindowids.nelements()==1 && spectralwindowids[0]<0)){
242 0 : Int nspw=thisms.spectralWindow().nrow();
243 0 : dataspectralwindowids_p.resize(nspw);
244 0 : indgen(dataspectralwindowids_p);
245 :
246 : }
247 :
248 : // If a selection has been made then close the current MS
249 : // and attach to a new selected MS. We do this on the original
250 : // MS.
251 : //I don't think i need this if statement
252 : // if(datafieldids_p.nelements()>0||datadescids_p.nelements()>0) {
253 59 : os << "Performing selection on MeasurementSet : " << thisms.tableName() << LogIO::POST;
254 : //if(vs_p) delete vs_p; vs_p=0;
255 : //if(mssel_p) delete mssel_p;
256 59 : mssel_p=0;
257 :
258 : // check that sorted table exists (it should), if not, make it now.
259 : //this->makeVisSet(thisms);
260 :
261 : //if you want to use scratch col...make sure they are there
262 59 : if(useModelCol_p){
263 0 : VisSetUtil::addScrCols(thisms, true, false, true, false);
264 0 : VisModelData::clearModel(thisms);
265 : }
266 : //MeasurementSet sorted=thisms.keywordSet().asTable("SORTED_TABLE");
267 :
268 :
269 : //Some MSSelection
270 59 : MSSelection thisSelection;
271 59 : if(datafieldids_p.nelements() > 0){
272 0 : thisSelection.setFieldExpr(MSSelection::indexExprStr(datafieldids_p));
273 0 : os << "Selecting on field ids : " << datafieldids_p << LogIO::POST;
274 : }
275 59 : if(fieldnames != ""){
276 7 : thisSelection.setFieldExpr(fieldnames);
277 7 : os << "Selecting on fields : " << fieldnames << LogIO::POST;
278 : }
279 59 : if(dataspectralwindowids_p.nelements() > 0){
280 0 : thisSelection.setSpwExpr(MSSelection::indexExprStr(dataspectralwindowids_p));
281 0 : os << "Selecting on spectral windows" << LogIO::POST;
282 : }
283 59 : else if(spwstring != ""){
284 46 : thisSelection.setSpwExpr(spwstring);
285 46 : os << "Selecting on spectral windows expression :"<< spwstring << LogIO::POST;
286 : }
287 59 : if(antIndex.nelements() >0){
288 0 : thisSelection.setAntennaExpr( MSSelection::indexExprStr(antIndex) );
289 0 : os << "Selecting on antenna ids : " << antIndex << LogIO::POST;
290 : }
291 59 : if(antnames != ""){
292 22 : Vector<String> antNames(1, antnames);
293 : // thisSelection.setAntennaExpr(MSSelection::nameExprStr( antNames));
294 22 : thisSelection.setAntennaExpr(antnames);
295 22 : os << "Selecting on antenna names : " << antnames << LogIO::POST;
296 :
297 22 : }
298 59 : if(timerng != ""){
299 0 : thisSelection.setTimeExpr(timerng);
300 0 : os << "Selecting on time range : " << timerng << LogIO::POST;
301 : }
302 59 : if(uvdist != ""){
303 0 : thisSelection.setUvDistExpr(uvdist);
304 0 : os << "Selecting on uvdist : " << uvdist << LogIO::POST;
305 : }
306 59 : if(scan != ""){
307 7 : thisSelection.setScanExpr(scan);
308 7 : os << "Selecting on scan : " << scan << LogIO::POST;
309 : }
310 59 : if(intent != ""){
311 36 : thisSelection.setStateExpr(intent);
312 36 : os << "Selecting on State(scan intent) Expr " << intent << LogIO::POST;
313 : }
314 59 : if(obs != ""){
315 0 : thisSelection.setObservationExpr(obs);
316 0 : os << "Selecting on Observation Expr : " << obs << LogIO::POST;
317 : }
318 59 : if(msSelect != ""){
319 0 : thisSelection.setTaQLExpr(msSelect);
320 0 : os << "Selecting via TaQL : " << msSelect << LogIO::POST;
321 : }
322 : //***************
323 :
324 59 : TableExprNode exprNode=thisSelection.toTableExprNode(&thisms);
325 : //if(exprNode.isNull())
326 : // throw(AipsError("Selection failed...review ms and selection parameters"));
327 59 : datafieldids_p.resize();
328 59 : datafieldids_p=thisSelection.getFieldList();
329 59 : if(datafieldids_p.nelements()==0){
330 52 : Int nf=ms_p->field().nrow();
331 52 : datafieldids_p.resize(nf);
332 52 : indgen(datafieldids_p);
333 : }
334 59 : if((numMS_p > 1) || datafieldids_p.nelements() > 1)
335 45 : multiFields_p= true;
336 : //Now lets see what was selected as spw and match it with datadesc
337 : //dataspectralwindowids_p.resize();
338 : //dataspectralwindowids_p=thisSelection.getSpwList();
339 59 : Matrix<Int> chansels=thisSelection.getChanList(NULL, 1, true);
340 59 : mssFreqSel_p.resize();
341 59 : mssFreqSel_p=thisSelection.getChanFreqList(NULL, true);
342 :
343 59 : uInt nms = numMS_p;
344 59 : uInt nrow = chansels.nrow();
345 59 : dataspectralwindowids_p.resize();
346 59 : const MSSpWindowColumns spwc(thisms.spectralWindow());
347 59 : uInt nspw = spwc.nrow();
348 59 : const ScalarColumn<Int> spwNchans(spwc.numChan());
349 59 : Vector<Int> nchanvec = spwNchans.getColumn();
350 : //cerr<<"SetDataOnThisMS::numMS_p="<<numMS_p<<" nchanvec="<<nchanvec<<endl;
351 59 : Int maxnchan = 0;
352 :
353 790 : for (uInt i=0;i<nchanvec.nelements();i++) {
354 731 : maxnchan=max(nchanvec[i],maxnchan);
355 : }
356 : //cout<<"spwchansels_p.shape()="<<spwchansels_p.shape()<<endl;
357 59 : uInt maxnspw = 0;
358 59 : if (numMS_p ==1) {
359 35 : maxnspw=nspw;
360 : }
361 : else {
362 50 : for (Int i=0;i<numMS_p-1;i++) {
363 26 : maxnspw=max(maxnspw,blockSpw_p[i].nelements());
364 : }
365 24 : maxnspw=max(nspw,maxnspw);
366 : }
367 59 : spwchansels_p.resize(nms,maxnspw,maxnchan,true);
368 : //cout<<"After resize: spwchansels_p.shape()="<<spwchansels_p.shape()<<endl;
369 59 : uInt nselspw=0;
370 59 : if (nrow==0) {
371 : //no channel selection, select all channels
372 13 : spwchansels_p.yzPlane(numMS_p-1)=1;
373 13 : dataspectralwindowids_p=thisSelection.getSpwList();
374 : }
375 : else {
376 46 : spwchansels_p.yzPlane(numMS_p-1)=0;
377 46 : Int prvspwid=-1;
378 46 : Vector<Int> selspw;
379 92 : for (uInt i=0;i<nrow;i++) {
380 46 : Vector<Int> sel = chansels.row(i);
381 46 : Int spwid = sel[0];
382 46 : if((sel[1] >= nchanvec[spwid]) || (sel[2] >=nchanvec[spwid]))
383 0 : throw(AipsError("Unexpected selection in spw selection of spwid "+String::toString(spwid)));
384 46 : if (spwid != prvspwid){
385 46 : nselspw++;
386 46 : selspw.resize(nselspw,true);
387 46 : selspw[nselspw-1]=spwid;
388 : }
389 46 : uInt minc= sel[1];
390 46 : uInt maxc = sel[2];
391 46 : uInt step = sel[3];
392 : // step as the same context as in im.selectvis
393 : // select channels
394 1115 : for (uInt k=minc;k<(maxc+1);k+=step) {
395 1069 : spwchansels_p(numMS_p-1,spwid,k)=1;
396 : }
397 46 : prvspwid=spwid;
398 46 : }
399 46 : dataspectralwindowids_p=selspw;
400 46 : }
401 :
402 : //cout<<"spwchansels_p(before shifting)="<<spwchansels_p<<endl;
403 59 : if(dataspectralwindowids_p.nelements()==0){
404 13 : Int nspwinms=thisms.spectralWindow().nrow();
405 13 : dataspectralwindowids_p.resize(nspwinms);
406 13 : indgen(dataspectralwindowids_p);
407 : }
408 :
409 : // old code
410 : /***
411 : if(dataspectralwindowids_p.nelements()==0){
412 : Int nspwinms=thisms.spectralWindow().nrow();
413 : dataspectralwindowids_p.resize(nspwinms);
414 : indgen(dataspectralwindowids_p);
415 : }
416 : ***/
417 :
418 : // Map the selected spectral window ids to data description ids
419 59 : MSDataDescIndex msDatIndex(thisms.dataDescription());
420 59 : datadescids_p.resize(0);
421 59 : datadescids_p=msDatIndex.matchSpwId(dataspectralwindowids_p);
422 :
423 59 : freqrange_p.resize(nms,2,true);
424 59 : if(mode=="none"){
425 : //check if we can find channel selection in the spw string
426 : //if(nselspw==dataspectralwindowids_p.nelements()){
427 59 : dataMode_p="channel";
428 59 : dataStep_p.resize(dataspectralwindowids_p.nelements());
429 59 : dataStart_p.resize(dataspectralwindowids_p.nelements());
430 59 : dataNchan_p.resize(dataspectralwindowids_p.nelements());
431 59 : Double fmin=C::dbl_max;
432 59 : Double fmax=-(C::dbl_max);
433 :
434 59 : Cube<Int> spwchansels_tmp=spwchansels_p;
435 121 : for (uInt k =0 ; k < dataspectralwindowids_p.nelements(); ++k){
436 62 : uInt curspwid=dataspectralwindowids_p[k];
437 62 : Vector<Double> chanFreq=spwc.chanFreq()(curspwid);
438 62 : Vector<Double> freqResolution = spwc.chanWidth()(curspwid);
439 :
440 :
441 : //dataStep_p[k]=chanselmat.row(k)(3);
442 62 : if (nrow==0) {
443 : //dataStep_p=step[0];
444 16 : dataStep_p[k]=step[0];
445 : }
446 : else {
447 46 : dataStep_p[k]=chansels.row(k)(3);
448 : }
449 : //if(dataStep_p[k] < 1)
450 : // dataStep_p[k]=1;
451 62 : dataStart_p[k]=0;
452 62 : dataNchan_p[k]=nchanvec(curspwid);
453 : //cout<<"SetDataOnThisMS: initial setting dataNchan_p["<<k<<"]="<<dataNchan_p[k]<<endl;
454 : //find start
455 62 : Bool first = true;
456 62 : uInt nchn = 0;
457 62 : uInt lastchan = 0;
458 7285 : for (Int j=0 ; j < nchanvec(curspwid); j++) {
459 7223 : if (spwchansels_p(numMS_p-1,curspwid,j)==1) {
460 7223 : if(first) {
461 62 : dataStart_p[k]=j;
462 62 : first = false;
463 : }
464 7223 : lastchan=j;
465 7223 : nchn++;
466 : }
467 : }
468 : //dataStart_p[k]=chanselmat.row(k)(1);
469 : //dataNchan_p[k]=Int(ceil(Double(chanselmat.row(k)(2)-dataStart_p[k]))/Double(dataStep_p[k]))+1;
470 62 : dataNchan_p[k]=Int(ceil(Double(lastchan-dataStart_p[k])/Double(dataStep_p[k])))+1;
471 : //cout<<"SetDataOnThisMS: after recalc. of nchan dataNchan_p["<<k<<"]="<<dataNchan_p[k]<<endl;
472 : //if(dataNchan_p[k]<1)
473 : // dataNchan_p[k]=1;
474 : //
475 : //Since msselet will be applied to the data before flags from spwchansels_p
476 : //are applied to the data in FTMachine, shift spwchansels_p by dataStart_p
477 7285 : for (Int j=0 ; j < nchanvec(curspwid); j++){
478 7223 : if (j<nchanvec(curspwid)-dataStart_p[k]) {
479 7223 : spwchansels_tmp(numMS_p-1,curspwid,j) = spwchansels_p(numMS_p-1,curspwid,j+dataStart_p[k]);
480 : }
481 : else {
482 0 : spwchansels_tmp(numMS_p-1,curspwid,j) = 0;
483 : }
484 : }
485 : //for mfs mode need to keep fmin,max info for later image setup
486 : //Int lastchan=dataStart_p[k]+ dataNchan_p[k]*dataStep_p[k];
487 : Int endchanused;
488 62 : if (nrow==0) {
489 : // default spw case
490 16 : endchanused=nchanvec(curspwid);
491 : }
492 : else {
493 46 : endchanused=lastchan;
494 : }
495 7239 : for(Int j=dataStart_p[k] ; j < endchanused ; j+=dataStep_p[k]){
496 7177 : fmin=min(fmin,chanFreq[j]-abs(freqResolution[j]*(dataStep_p[k]-0.5)));
497 7177 : fmax=max(fmax,chanFreq[j]+abs(freqResolution[j]*(dataStep_p[k]-0.5)));
498 : }
499 62 : }
500 : //cerr<<"numMS_p="<<numMS_p<<" fmin="<<fmin<<" fmax="<<fmax<<endl;
501 59 : freqrange_p(numMS_p-1,0)=fmin;
502 59 : freqrange_p(numMS_p-1,1)=fmax;
503 :
504 59 : spwchansels_p=spwchansels_tmp;
505 : //}
506 59 : }
507 :
508 :
509 :
510 : // Now remake the selected ms
511 59 : if(!(exprNode.isNull())){
512 47 : mssel_p = new MeasurementSet(thisms(exprNode));
513 : }
514 : else{
515 : // Null take all the ms ...setdata() blank means that
516 12 : mssel_p = new MeasurementSet(thisms);
517 : }
518 59 : AlwaysAssert(!mssel_p.null(), AipsError);
519 59 : if(mssel_p->nrow()==0) {
520 : //delete mssel_p;
521 0 : mssel_p=0;
522 :
523 : os << "Selection is empty: you may want to review this MSs selection"
524 0 : << LogIO::EXCEPTION;
525 : }
526 : else {
527 59 : mssel_p->flush();
528 : }
529 :
530 59 : if(mssel_p->nrow()!=thisms.nrow()) {
531 6 : os << "Selected " << mssel_p->nrow() << " out of "
532 3 : << thisms.nrow() << " rows." << LogIO::POST;
533 : }
534 : else {
535 56 : os << "Selected all " << mssel_p->nrow() << " rows" << LogIO::POST;
536 : }
537 : // }
538 :
539 : // Tell the user how many channels have been selected.
540 : // NOTE : This code is replicated in Imager.cc.
541 59 : Vector<Int> chancounts(dataspectralwindowids_p.nelements());
542 59 : chancounts=0;
543 : // if( spwstring == "" ) os << "Selected all spws and channels" << LogIO::POST;
544 : // else os << "Channel selection : " << spwstring << LogIO::POST;
545 59 : os << "Selected :";
546 121 : for(uInt k=0;k<dataspectralwindowids_p.nelements();k++)
547 : {
548 7285 : for(uInt ch=0;ch<uInt(nchanvec(dataspectralwindowids_p[k]));ch++)
549 7223 : {if(spwchansels_p(numMS_p-1,dataspectralwindowids_p[k],ch)) chancounts[k]++; }
550 : //if(chancounts[k]<1)
551 : // throw(AipsError("bad selection in spw "+ String::toString( dataspectralwindowids_p[k])));
552 62 : os << " [" << chancounts[k] << " chans in spw " << dataspectralwindowids_p[k] << "]";
553 : // os << "Selected " << chancounts[k] << " chans in spw "
554 : // << dataspectralwindowids_p[k] << LogIO::POST;
555 : }
556 59 : os << LogIO::POST;
557 :
558 :
559 59 : blockMSSel_p[numMS_p-1]=*mssel_p;
560 : //lets make the visSet now
561 59 : Block<Matrix<Int> > noChanSel;
562 59 : noChanSel.resize(numMS_p);
563 59 : Block<Int> sort(0);
564 : //if(vs_p) delete vs_p; vs_p=0;
565 59 : if(rvi_p) delete rvi_p;
566 59 : rvi_p=0;
567 59 : wvi_p=0;
568 :
569 : //vs_p= new VisSet(blockMSSel_p, sort, noChanSel, useModelCol_p);
570 59 : if(!(mssel_p->isWritable())){
571 0 : rvi_p=new ROVisibilityIterator(blockMSSel_p, sort);
572 :
573 : }
574 : else{
575 59 : wvi_p=new VisibilityIterator(blockMSSel_p, sort);
576 59 : rvi_p=wvi_p;
577 : }
578 :
579 59 : if(imwgt_p.getType()=="none")
580 35 : imwgt_p=VisImagingWeight("natural");
581 59 : rvi_p->useImagingWeight(imwgt_p);
582 : //rvi_p->slurp();
583 :
584 59 : selectDataChannel();
585 59 : dataSet_p=true;
586 :
587 59 : return dataSet_p;
588 59 : }
589 0 : catch(AipsError x){
590 : //Ayeee...lets back out of this one
591 0 : --numMS_p;
592 0 : blockMSSel_p.resize(numMS_p, true);
593 0 : blockNChan_p.resize(numMS_p, true);
594 0 : blockStart_p.resize(numMS_p, true);
595 0 : blockStep_p.resize(numMS_p, true);
596 0 : blockSpw_p.resize(numMS_p, true);
597 : //point it back to the previous ms
598 0 : if(numMS_p >0){
599 0 : mssel_p=new MeasurementSet(blockMSSel_p[numMS_p-1]);
600 : }
601 0 : throw(AipsError(x));
602 0 : }
603 :
604 59 : } //End of setDataPerMS
605 :
606 :
607 0 : Bool ImagerMultiMS::setimage(const Int nx, const Int ny,
608 : const Quantity& cellx, const Quantity& celly,
609 : const String& stokes,
610 : Bool doShift,
611 : const MDirection& phaseCenter,
612 : const Quantity& shiftx, const Quantity& shifty,
613 : const String& mode, const Int nchan,
614 : const Int start, const Int step,
615 : const MRadialVelocity& mStart,
616 : const MRadialVelocity& mStep,
617 : const Vector<Int>& spectralwindowids,
618 : const Int fieldid,
619 : const Int facets,
620 : const Quantity& distance)
621 : {
622 :
623 0 : if(!dataSet_p){
624 0 : LogIO os(LogOrigin("imagerMultiMS", "setimage()"), logSink_p);
625 : os << LogIO::SEVERE
626 : << "Please use setdata before setimage as imager need one ms at least "
627 0 : << LogIO::POST;
628 0 : return false;
629 :
630 0 : }
631 :
632 0 : Bool returnval=Imager::setimage(nx, ny, cellx, celly, stokes,doShift,
633 : phaseCenter, shiftx, shifty, mode, nchan,
634 : start, step, mStart, mStep,
635 : spectralwindowids,
636 : fieldid, facets, distance);
637 :
638 0 : return returnval;
639 : }
640 :
641 105 : Bool ImagerMultiMS::lock(){
642 : //Do nothing for now as its Autolocked..but a better scheme is necessary
643 : //for parallel access to same data for modification
644 :
645 105 : return true;
646 : }
647 :
648 105 : Bool ImagerMultiMS::unlock(){
649 :
650 325 : for (uInt k=0; k < blockMSSel_p.nelements(); ++k){
651 220 : blockMSSel_p[k].relinquishAutoLocks(true);
652 : }
653 105 : return true;
654 :
655 : }
656 :
657 :
658 59 : Bool ImagerMultiMS::selectDataChannel(){
659 :
660 118 : LogIO os(LogOrigin("imager", "selectDataChannel()"), logSink_p);
661 :
662 59 : (blockNChan_p[numMS_p-1]).resize(dataspectralwindowids_p.nelements());
663 59 : (blockStart_p[numMS_p-1]).resize(dataspectralwindowids_p.nelements());
664 59 : (blockStart_p[numMS_p-1]).set(0);
665 59 : (blockStep_p[numMS_p-1]).resize(dataspectralwindowids_p.nelements());
666 59 : (blockStep_p[numMS_p-1]).set(1);
667 59 : (blockSpw_p[numMS_p-1]).resize(dataspectralwindowids_p.nelements());
668 59 : (blockSpw_p[numMS_p-1]).resize();
669 59 : (blockSpw_p[numMS_p-1])=dataspectralwindowids_p;
670 : {
671 : //Fill default numChan for now
672 :
673 59 : MSSpWindowColumns msSpW(ms_p->spectralWindow());
674 59 : Vector<Int> numChan=msSpW.numChan().getColumn();
675 121 : for (uInt k=0; k < dataspectralwindowids_p.nelements(); ++k){
676 62 : blockNChan_p[numMS_p-1][k]=numChan[dataspectralwindowids_p[k]];
677 :
678 : }
679 59 : }
680 :
681 59 : if(dataMode_p=="channel") {
682 :
683 59 : if (dataNchan_p.nelements() != dataspectralwindowids_p.nelements()){
684 0 : if(dataNchan_p.nelements()==1){
685 0 : dataNchan_p.resize(dataspectralwindowids_p.nelements(), true);
686 0 : for(uInt k=1; k < dataspectralwindowids_p.nelements(); ++k){
687 0 : dataNchan_p[k]=dataNchan_p[0];
688 : }
689 : }
690 : else{
691 : os << LogIO::SEVERE
692 : << "Vector of nchan has to be of size 1 or be of the same shape as spw "
693 0 : << LogIO::POST;
694 0 : return false;
695 : }
696 : }
697 59 : if (dataStart_p.nelements() != dataspectralwindowids_p.nelements()){
698 0 : if(dataStart_p.nelements()==1){
699 0 : dataStart_p.resize(dataspectralwindowids_p.nelements(), true);
700 0 : for(uInt k=1; k < dataspectralwindowids_p.nelements(); ++k){
701 0 : dataStart_p[k]=dataStart_p[0];
702 : }
703 : }
704 : else{
705 : os << LogIO::SEVERE
706 : << "Vector of start has to be of size 1 or be of the same shape as spw "
707 0 : << LogIO::POST;
708 0 : return false;
709 : }
710 : }
711 59 : if (dataStep_p.nelements() != dataspectralwindowids_p.nelements()){
712 0 : if(dataStep_p.nelements()==1){
713 0 : dataStep_p.resize(dataspectralwindowids_p.nelements(), true);
714 0 : for(uInt k=1; k < dataspectralwindowids_p.nelements(); ++k){
715 0 : dataStep_p[k]=dataStep_p[0];
716 : }
717 : }
718 : else{
719 : os << LogIO::SEVERE
720 : << "Vector of step has to be of size 1 or be of the same shape as spw "
721 0 : << LogIO::POST;
722 0 : return false;
723 : }
724 : }
725 :
726 : {
727 59 : Int nch=0;
728 121 : for(uInt i=0;i<dataspectralwindowids_p.nelements();i++) {
729 62 : Int spwid=dataspectralwindowids_p(i);
730 62 : if(dataStart_p[i]<0) {
731 : os << LogIO::SEVERE << "Illegal start pixel = "
732 0 : << dataStart_p[i] << " for spw " << spwid
733 0 : << LogIO::POST;
734 0 : return false;
735 : }
736 :
737 62 : if(dataNchan_p[i]<=0){
738 0 : if(dataStep_p[i] <= 0)
739 0 : dataStep_p[i]=1;
740 0 : nch=(blockNChan_p[numMS_p-1](i)-dataStart_p[i])/Int(dataStep_p[i])+1;
741 :
742 : }
743 62 : else nch = dataNchan_p[i];
744 62 : while((nch*dataStep_p[i]+dataStart_p[i]) >
745 62 : blockNChan_p[numMS_p-1](i)){
746 0 : --nch;
747 : }
748 62 : Int end = Int(dataStart_p[i]) + Int(nch-1) * Int(dataStep_p[i]);
749 62 : if(end < 0 || end > (blockNChan_p[numMS_p-1](i)-1)) {
750 0 : os << LogIO::SEVERE << "Illegal step pixel = " << dataStep_p[i]
751 : << " for spw " << spwid
752 : << "\n end channel " << end << " out of range "
753 0 : << dataStart_p[i] << " to "
754 0 : << (blockNChan_p[numMS_p-1](i)-1)
755 :
756 0 : << LogIO::POST;
757 0 : return false;
758 : }
759 : os << LogIO::DEBUG1 << "Selecting "<< nch
760 : << " channels, starting at visibility channel "
761 124 : << dataStart_p[i] << " stepped by "
762 62 : << dataStep_p[i] << " for spw " << spwid << LogIO::POST;
763 62 : dataNchan_p[i]=nch;
764 62 : blockNChan_p[numMS_p-1][i]=nch;
765 62 : blockStep_p[numMS_p-1][i]=dataStep_p[i];
766 62 : blockStart_p[numMS_p-1][i]=dataStart_p[i];
767 : }
768 : }
769 :
770 : }
771 :
772 59 : Block<Vector<Int> > blockGroup(numMS_p);
773 144 : for (Int k=0; k < numMS_p; ++k){
774 85 : blockGroup[k].resize(blockSpw_p[k].nelements());
775 85 : blockGroup[k].set(1);
776 : }
777 :
778 59 : rvi_p->selectChannel(blockGroup, blockStart_p, blockNChan_p,
779 59 : blockStep_p, blockSpw_p);
780 59 : return true;
781 59 : }
782 :
783 59 : Bool ImagerMultiMS::openSubTables(){
784 :
785 : // Reopen disk-resident tables with the default table lock
786 :
787 59 : if((ms_p->tableType()) != Table::Memory){
788 :
789 59 : TableLock tableLock; // default table lock
790 :
791 59 : openSubTable (ms_p->antenna(), antab_p, tableLock);
792 59 : openSubTable (ms_p->dataDescription (), datadesctab_p, tableLock);
793 59 : openSubTable (ms_p->doppler(), dopplertab_p, tableLock);
794 59 : openSubTable (ms_p->feed(), feedtab_p, tableLock);
795 59 : openSubTable (ms_p->field(), fieldtab_p, tableLock);
796 59 : openSubTable (ms_p->flagCmd(), flagcmdtab_p, tableLock);
797 59 : openSubTable (ms_p->freqOffset(), freqoffsettab_p, tableLock);
798 59 : openSubTable (ms_p->observation(), obstab_p, tableLock);
799 59 : openSubTable (ms_p->pointing(), pointingtab_p, tableLock);
800 59 : openSubTable (ms_p->polarization(), poltab_p, tableLock);
801 59 : openSubTable (ms_p->processor(), proctab_p, tableLock);
802 59 : openSubTable (ms_p->source(), sourcetab_p, tableLock);
803 59 : openSubTable (ms_p->spectralWindow(), spwtab_p, tableLock);
804 59 : openSubTable (ms_p->state(), statetab_p, tableLock);
805 59 : openSubTable (ms_p->sysCal(), syscaltab_p, tableLock);
806 59 : openSubTable (ms_p->weather(), weathertab_p, tableLock);
807 :
808 : // Handle the history table specially
809 :
810 59 : if(ms_p->isWritable()){
811 :
812 59 : if(!(Table::isReadable(ms_p->historyTableName()))){
813 :
814 : // setup a new table in case its not there
815 0 : TableRecord &kws = ms_p->rwKeywordSet();
816 0 : SetupNewTable historySetup(ms_p->historyTableName(),
817 0 : MSHistory::requiredTableDesc(),Table::New);
818 0 : kws.defineTable(MS::keywordName(MS::HISTORY), Table(historySetup));
819 :
820 0 : }
821 :
822 59 : historytab_p=Table(ms_p->historyTableName(), TableLock(), Table::Update);
823 :
824 59 : hist_p= new MSHistoryHandler(*ms_p, "imager");
825 : }
826 : }
827 :
828 59 : return true;
829 : }
830 :
831 : #define INITIALIZE_DIRECTION_VECTOR(name) \
832 : do { \
833 : if (name.nelements() < 2) { \
834 : name.resize(2); \
835 : name = 0.0; \
836 : } \
837 : } while (false)
838 15 : Bool ImagerMultiMS::mapExtent(const String &referenceFrame, const String &movingSource,
839 : const String &pointingColumn, Vector<Double> ¢er, Vector<Double> &blc,
840 : Vector<Double> &trc, Vector<Double> &extent) {
841 : // initialize if necessary
842 15 : INITIALIZE_DIRECTION_VECTOR(center);
843 15 : INITIALIZE_DIRECTION_VECTOR(blc);
844 15 : INITIALIZE_DIRECTION_VECTOR(trc);
845 15 : INITIALIZE_DIRECTION_VECTOR(extent);
846 :
847 : try {
848 15 : Bool isValueSet = false;
849 32 : for (size_t i = 0; i < blockMSSel_p.nelements(); ++i) {
850 : //cout << "start MS " << i << endl;
851 17 : MeasurementSet ms = blockMSSel_p[i];
852 17 : Vector<Double> wcenter(2);
853 17 : Vector<Double> wblc(2);
854 17 : Vector<Double> wtrc(2);
855 17 : Vector<Double> wextent(2);
856 : //cout << "run getMapExtent" << endl;
857 17 : Bool status = getMapExtent(ms, referenceFrame, movingSource, pointingColumn,
858 : wcenter, wblc, wtrc, wextent);
859 : //cout << "done Imager::mapExtent" << endl;
860 17 : if (status) {
861 17 : if (!isValueSet) {
862 15 : center = wcenter;
863 15 : blc = wblc;
864 15 : trc = wtrc;
865 15 : extent = wextent;
866 15 : isValueSet = true;
867 : }
868 : else {
869 2 : blc[0] = min(blc[0], wblc[0]);
870 2 : blc[1] = min(blc[1], wblc[1]);
871 2 : trc[0] = max(trc[0], wtrc[0]);
872 2 : trc[1] = max(trc[1], wtrc[1]);
873 : }
874 : }
875 17 : }
876 :
877 15 : if (!isValueSet) {
878 0 : LogIO os(LogOrigin("ImagerMultiMS", "mapExtent", WHERE));
879 0 : os << LogIO::SEVERE << "No valid data found. Failed." << LogIO::POST;
880 :
881 0 : return false;
882 0 : }
883 :
884 : // extent is re-evaluated using true trc and blc
885 15 : extent = trc - blc;
886 :
887 : // declination correction factor
888 15 : Double declinationCorrection = cos(center[1]);
889 :
890 : // center needs to be updated according to specified column name
891 : // and moving source name
892 : MDirection::Types refType;
893 15 : Bool status = MDirection::getType(refType, movingSource);
894 32 : Bool doMovingSourceCorrection = (status == true &&
895 17 : MDirection::N_Types < refType &&
896 2 : refType < MDirection::N_Planets);
897 15 : Bool isOffsetColumn = (pointingColumn.contains("OFFSET")
898 15 : || pointingColumn.contains("Offset")
899 30 : || pointingColumn.contains("offset"));
900 15 : if (isOffsetColumn) {
901 : // center is always (0,0)
902 0 : center = 0.0;
903 0 : declinationCorrection = cos((blc[1] + trc[1]) / 2.0);
904 : }
905 15 : else if (!doMovingSourceCorrection) {
906 : // initial center value should be kept
907 : // if doMovingSourceCorrection is true.
908 : // otherwise, center should be re-evaluated
909 : // using true trc and blc
910 13 : center = (blc + trc) / 2.0;
911 13 : declinationCorrection = cos(center[1]);
912 : }
913 :
914 15 : assert(declinationCorrection != 0.0);
915 : //cout << "declinationCorrection = " << declinationCorrection << endl;
916 :
917 : // apply declinationCorrection to extent[0]
918 15 : extent[0] *= declinationCorrection;
919 : }
920 0 : catch (...) {
921 0 : LogIO os(LogOrigin("ImagerMultiMS", "mapExtent", WHERE));
922 0 : os << LogIO::SEVERE << "Failed due to unknown error" << LogIO::POST;
923 0 : return false;
924 0 : }
925 15 : return true;
926 : }
927 : } //# NAMESPACE CASA - END
|
