Line data Source code
1 : //# MSContinuumSubtractor.cc: Subtract continuum from spectral line data
2 : //# Copyright (C) 2004
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library 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 Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 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 <msvis/MSVis/MSContinuumSubtractor.h>
29 :
30 : #include <casacore/casa/Quanta/QuantumHolder.h>
31 : #include <casacore/casa/Containers/RecordFieldId.h>
32 : #include <casacore/measures/Measures/Stokes.h>
33 : #include <casacore/ms/MeasurementSets/MSColumns.h>
34 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
35 : #include <casacore/ms/MSSel/MSSelector.h>
36 : #include <casacore/ms/MSSel/MSSelection.h>
37 : #include <casacore/ms/MSSel/MSSelectionTools.h>
38 : #include <msvis/MSVis/VisSet.h>
39 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
40 : #include <msvis/MSVis/ViFrequencySelection.h>
41 : #include <msvis/MSVis/IteratingParameters.h>
42 : #include <msvis/MSVis/AveragingVi2Factory.h>
43 : #include <msvis/MSVis/LayeredVi2Factory.h>
44 : #include <casacore/scimath/Fitting/LinearFit.h>
45 : #include <casacore/scimath/Functionals/Polynomial.h>
46 : #include <casacore/casa/Arrays/ArrayLogical.h>
47 : #include <casacore/casa/Arrays/ArrayMath.h>
48 : #include <casacore/casa/Arrays/ArrayUtil.h>
49 : #include <casacore/casa/Containers/Record.h>
50 : #include <casacore/casa/Exceptions/Error.h>
51 : #include <casacore/casa/Logging/LogIO.h>
52 : #include <casacore/tables/TaQL/TableParse.h>
53 : #include <iostream>
54 :
55 : using namespace casacore;
56 : namespace casa { //# NAMESPACE CASA - BEGIN
57 :
58 : //
59 : // Constructor assigns pointer (if MS goes out of scope you will get rubbish)
60 0 : MSContinuumSubtractor::MSContinuumSubtractor (MeasurementSet& ms)
61 0 : : ms_p(&ms),itsSolInt(0.0),itsOrder(0),itsMode("subtract")
62 : {
63 :
64 : // Make a VisSet so scratch columns are created
65 0 : Block<Int> nosort(0);
66 0 : Matrix<Int> noselection;
67 0 : Double timeInterval=0.0;
68 0 : Bool compress(False);
69 0 : VisSet vs(ms,nosort,noselection,timeInterval,compress);
70 :
71 0 : nSpw_= vs.numberSpw();
72 :
73 0 : }
74 :
75 :
76 : //
77 : // Assignment operator
78 : //
79 0 : MSContinuumSubtractor& MSContinuumSubtractor::operator=(MSContinuumSubtractor& other)
80 : {
81 0 : if (this==&other) return *this;
82 0 : ms_p = other.ms_p;
83 0 : return *this;
84 : }
85 :
86 :
87 : //
88 : // Destructor does nothing
89 : //
90 0 : MSContinuumSubtractor::~MSContinuumSubtractor()
91 0 : {}
92 :
93 : // Set the required field Ids
94 0 : void MSContinuumSubtractor::setField(const String& field)
95 : {
96 :
97 0 : MSSelection mssel;
98 0 : mssel.setFieldExpr(field);
99 0 : Vector<Int> fldlist;
100 0 : fldlist=mssel.getFieldList(ms_p);
101 :
102 0 : setFields(fldlist);
103 :
104 0 : }
105 :
106 0 : void MSContinuumSubtractor::setFields(const Vector<Int>& fieldIds)
107 : {
108 0 : itsFieldIds = fieldIds;
109 0 : }
110 :
111 :
112 : // Set the channels to use in the fit
113 0 : void MSContinuumSubtractor::setFitSpw(const String& fitspw)
114 : {
115 : // NB: this method assumes spwids == ddids!
116 :
117 : // Using MSSelection
118 0 : MSSelection mssel;
119 0 : mssel.setSpwExpr(fitspw);
120 0 : Vector<Int> spwlist;
121 0 : spwlist=mssel.getSpwList(ms_p);
122 :
123 0 : if (spwlist.nelements()==0) {
124 0 : spwlist.resize(nSpw_);
125 0 : indgen(spwlist);
126 : }
127 :
128 0 : setDataDescriptionIds(spwlist);
129 :
130 0 : itsFitChans= mssel.getChanList(ms_p);
131 :
132 0 : }
133 :
134 : // Set the channels from which the fit should be subtracted
135 0 : void MSContinuumSubtractor::setSubSpw(const String& subspw)
136 : {
137 : // Using MSSelection
138 0 : MSSelection mssel;
139 0 : mssel.setSpwExpr(subspw);
140 0 : itsSubChans= mssel.getChanList(ms_p);
141 :
142 0 : }
143 :
144 0 : void MSContinuumSubtractor::setDataDescriptionIds(const Vector<Int>& ddIds)
145 : {
146 0 : itsDDIds = ddIds;
147 0 : }
148 :
149 : // Set the solution interval in seconds, the value zero implies scan averaging
150 0 : void MSContinuumSubtractor::setSolutionInterval(Float solInt)
151 : {
152 0 : itsSolInt = solInt;
153 0 : }
154 :
155 : // Set the solution interval in seconds, the value zero implies scan averaging
156 0 : void MSContinuumSubtractor::setSolutionInterval(String solInt)
157 : {
158 :
159 0 : LogIO os(LogOrigin("MSContinuumSubtractor","setSolutionInterval"));
160 :
161 0 : os << LogIO::NORMAL << "Fitting continuum on ";
162 :
163 0 : if (upcase(solInt).contains("INF")) {
164 : // ~Infinite (this actually means per-scan in UV contsub)
165 0 : solInt="inf";
166 0 : itsSolInt=0.0;
167 0 : os <<"per-scan ";
168 : }
169 0 : else if (upcase(solInt).contains("INT")) {
170 : // Per integration
171 0 : itsSolInt=-1.0;
172 0 : os << "per-integration ";
173 : }
174 : else {
175 : // User-selected timescale
176 0 : QuantumHolder qhsolint;
177 0 : String error;
178 0 : Quantity qsolint;
179 0 : qhsolint.fromString(error,solInt);
180 0 : if (error.length()!=0)
181 0 : throw(AipsError("Unrecognized units for solint."));
182 :
183 0 : qsolint=qhsolint.asQuantumDouble();
184 :
185 0 : if (qsolint.isConform("s"))
186 0 : itsSolInt=qsolint.get("s").getValue();
187 : else {
188 : // assume seconds
189 0 : itsSolInt=qsolint.getValue();
190 : }
191 :
192 0 : os << itsSolInt <<"-second (per scan) ";
193 : }
194 :
195 0 : os << "timescale." << LogIO::POST;
196 :
197 0 : }
198 :
199 :
200 :
201 : // Set the order of the fit (1=linear)
202 0 : void MSContinuumSubtractor::setOrder(Int order)
203 : {
204 0 : itsOrder = order;
205 0 : }
206 :
207 : // Set the processing mode: subtract, model or replace
208 0 : void MSContinuumSubtractor::setMode(const String& mode)
209 : {
210 0 : itsMode = mode;
211 0 : }
212 :
213 : // Do the subtraction (or save the model)
214 0 : void MSContinuumSubtractor::subtract()
215 : {
216 :
217 0 : LogIO os(LogOrigin("MSContinuumSubtractor","subtract"));
218 : os << LogIO::NORMAL<< "MSContinuumSubtractor::subtract() - parameters:"
219 : << "ddIds="<<itsDDIds<<", fieldIds="<<itsFieldIds
220 : << ", order="<<itsOrder
221 0 : << ", mode="<<itsMode<<LogIO::POST;
222 :
223 0 : if (itsFitChans.nelements()>0) {
224 0 : os<<"fit channels: " << LogIO::POST;
225 0 : for (uInt i=0;i<itsFitChans.nrow();++i)
226 0 : os<<" spw="<<itsFitChans.row(i)(0)
227 0 : <<": "<<itsFitChans.row(i)(1)<<"~"<<itsFitChans.row(i)(2)
228 0 : << LogIO::POST;
229 : }
230 :
231 0 : ostringstream select;
232 0 : select <<"select from $1 where ANTENNA1!=ANTENNA2";
233 0 : if (itsFieldIds.nelements()>0) {
234 0 : select<<" && FIELD_ID IN ["<<itsFieldIds(0);
235 0 : for (uInt j=1; j<itsFieldIds.nelements(); j++) select<<", "<<itsFieldIds(j);
236 0 : select<<"]";
237 : }
238 0 : if (itsDDIds.nelements()>0) {
239 0 : select<<" && DATA_DESC_ID IN ["<<itsDDIds(0);
240 0 : for (uInt j=1; j<itsDDIds.nelements(); j++) select<<", "<<itsDDIds(j);
241 0 : select<<"]";
242 : }
243 : //os <<"Selection string: "<<select.str()<<LogIO::POST;
244 : //os <<" nrow="<<ms_p->nrow()<<LogIO::POST;
245 0 : MeasurementSet selectedMS(tableCommand(select,*ms_p).table());
246 0 : MSSelector msSel(selectedMS);
247 : //os <<" nrow="<<msSel.nrow()<<LogIO::POST;
248 0 : MSColumns msc(selectedMS);
249 0 : if (itsDDIds.nelements()>1) {
250 0 : cout<<"Processing "<<itsDDIds.nelements()<<" spectral windows"<<endl;
251 : }
252 0 : for (uInt iDD=0; iDD<itsDDIds.nelements(); iDD++) {
253 0 : Vector<Int> ddIDs(1,itsDDIds(iDD));
254 0 : msSel.initSelection(ddIDs);
255 : //os <<" nrow="<<msSel.nrow()<<LogIO::POST;
256 0 : if (msSel.nrow()==0) continue;
257 0 : Int nChan=msc.spectralWindow()
258 0 : .numChan()(msc.dataDescription().spectralWindowId()(ddIDs(0)));
259 0 : Vector<Int> corrTypes=msc.polarization().
260 0 : corrType()(msc.dataDescription().polarizationId()(ddIDs(0)));
261 :
262 :
263 : // default to fit all channels
264 0 : Vector<Bool> fitChanMask(nChan,True);
265 :
266 : // Handle non-trivial channel selection:
267 :
268 0 : if (itsFitChans.nelements()>0 && anyEQ(itsFitChans.column(0),itsDDIds(iDD))) {
269 : // If subset of channels selected, set mask all False...
270 0 : fitChanMask=False;
271 :
272 0 : IPosition blc(1,0);
273 0 : IPosition trc(1,0);
274 :
275 : // ... and set only selected channels True:
276 0 : for (uInt i=0;i<itsFitChans.nrow();++i) {
277 0 : Vector<Int> chansel(itsFitChans.row(i));
278 :
279 : // match current spwId/DDid
280 0 : if (chansel(0)==itsDDIds(iDD)) {
281 0 : blc(0)=chansel(1);
282 0 : trc(0)=chansel(2);
283 0 : fitChanMask(blc,trc)=True;
284 : }
285 : }
286 : }
287 :
288 : // cout << "fitChanMask = " << fitChanMask << endl;
289 :
290 : // default to subtract from all channels
291 0 : Vector<Bool> subChanMask(nChan,True);
292 0 : if (itsSubChans.nelements()>0 && anyEQ(itsSubChans.column(0),itsDDIds(iDD))) {
293 : // If subset of channels selected, set sub mask all False...
294 0 : subChanMask=False;
295 :
296 0 : IPosition blc(1,0);
297 0 : IPosition trc(1,0);
298 :
299 : // ... and set only selected sub channels True:
300 0 : for (uInt i=0;i<itsSubChans.nrow();++i) {
301 0 : Vector<Int> chansel(itsSubChans.row(i));
302 :
303 : // match current spwId/DDid
304 0 : if (chansel(0)==itsDDIds(iDD)) {
305 0 : blc(0)=chansel(1);
306 0 : trc(0)=chansel(2);
307 0 : subChanMask(blc,trc)=True;
308 : }
309 : }
310 : }
311 :
312 : // select parallel hand polarizations
313 0 : Vector<String> polSel(corrTypes.nelements());
314 0 : Int nPol = 0;
315 0 : for (uInt j=0; j<corrTypes.nelements(); j++) {
316 0 : if (corrTypes(j)==Stokes::XX||corrTypes(j)==Stokes::YY||
317 0 : corrTypes(j)==Stokes::RR||corrTypes(j)==Stokes::LL) {
318 0 : polSel(nPol++)=Stokes::name(Stokes::type(corrTypes(j)));
319 : }
320 : }
321 0 : polSel.resize(nPol,True);
322 0 : msSel.selectPolarization(polSel);
323 :
324 0 : msSel.iterInit(
325 0 : stringToVector("ARRAY_ID,DATA_DESC_ID,SCAN_NUMBER,FIELD_ID,TIME"),
326 0 : itsSolInt,0,False);
327 0 : msSel.iterOrigin();
328 0 : Int nIter=1;
329 0 : while (msSel.iterNext()) nIter++;
330 0 : os<<"Processing "<<nIter<<" slots."<< LogIO::POST;
331 0 : msSel.iterOrigin();
332 0 : do {
333 0 : Record avRec = msSel.getData(stringToVector("corrected_data"),True,0,1,True);
334 0 : Record dataRec = msSel.getData(stringToVector("model_data,corrected_data"),
335 0 : True,0,1);
336 0 : Array<Complex> avCorData(avRec.asArrayComplex("corrected_data"));
337 0 : Array<Complex> modelData(dataRec.asArrayComplex("model_data"));
338 0 : Array<Complex> correctedData(dataRec.asArrayComplex("corrected_data"));
339 0 : Int nTime=modelData.shape()[3];
340 0 : Int nIfr=modelData.shape()[2];
341 :
342 : // fit
343 0 : Vector<Float> x(nChan);
344 0 : for (Int i=0; i<nChan; i++) x(i)=i;
345 0 : LinearFit<Float> fitter;
346 0 : Polynomial<AutoDiff<Float> > apoly(itsOrder);
347 0 : fitter.setFunction(apoly);
348 0 : Polynomial<Float> poly(itsOrder);
349 :
350 0 : Vector<Float> y1(nChan),y2(nChan);
351 0 : Vector<Float> sol(itsOrder+1);
352 0 : Vector<Complex> tmp(nChan);
353 0 : IPosition start(3,0,0,0),end(3,0,nChan-1,0);
354 0 : for (; start[2]<nIfr; start[2]++,end[2]++) {
355 0 : for (start[0]=end[0]=0; start[0]<nPol; start[0]++,end[0]++) {
356 0 : Vector<Complex> c(avCorData(start,end).nonDegenerate());
357 0 : tmp=c; // copy into contiguous storage
358 0 : c.set(Complex(0.0)); // zero the input
359 0 : real(y1,tmp);
360 0 : imag(y2,tmp);
361 0 : sol = fitter.fit(x,y1,&fitChanMask);
362 0 : poly.setCoefficients(sol);
363 0 : y1=x;
364 0 : y1.apply(poly);
365 0 : sol = fitter.fit(x,y2,&fitChanMask);
366 0 : poly.setCoefficients(sol);
367 0 : y2=x;
368 0 : y2.apply(poly);
369 0 : for (Int chn=0; chn<nChan; chn++)
370 0 : if (subChanMask(chn))
371 0 : c(chn)=Complex(y1(chn),y2(chn));
372 : }
373 : }
374 :
375 0 : IPosition start4(4,0,0,0,0),end4(4,nPol-1,nChan-1,nIfr-1,0);
376 0 : for (Int iTime=0; iTime<nTime; iTime++,start4[3]++,end4[3]++) {
377 0 : if (itsMode!="replace") {
378 0 : Array<Complex> model(modelData(start4,end4).nonDegenerate(3));
379 0 : model=avCorData;
380 : }
381 0 : if (itsMode!="model") {
382 0 : Array<Complex> corr(correctedData(start4,end4).nonDegenerate(3));
383 0 : if (itsMode=="replace") corr=avCorData;
384 0 : if (itsMode=="subtract" || itsMode=="sub") corr-=avCorData;
385 : }
386 : }
387 :
388 0 : Record newDataRec;
389 0 : if (itsMode=="model"||itsMode=="subtract"||itsMode=="sub") {
390 0 : newDataRec.define("model_data",modelData);
391 : }
392 0 : if (itsMode=="replace"||itsMode=="subtract"||itsMode=="sub") {
393 0 : newDataRec.define("corrected_data",correctedData);
394 : }
395 0 : msSel.putData(newDataRec);
396 :
397 0 : } while (msSel.iterNext());
398 :
399 : }
400 0 : }
401 :
402 : // Do the subtraction (or save the model)
403 0 : void MSContinuumSubtractor::subtract2()
404 : {
405 : // Log selections for subtraction
406 0 : LogIO os(LogOrigin("MSContinuumSubtractor","subtract2"));
407 : os << LogIO::NORMAL << "parameters:"
408 0 : << LogIO::POST;
409 : os << LogIO::NORMAL << " ddIds=" << itsDDIds
410 : << ", fieldIds="<< itsFieldIds
411 : << ", order="<< itsOrder
412 0 : << ", mode="<< itsMode
413 0 : << LogIO::POST;
414 0 : if (itsFitChans.nelements()>0) {
415 0 : os << "fit channels: " << LogIO::POST;
416 0 : for (uInt i=0; i<itsFitChans.nrow(); ++i)
417 0 : os << " spw=" << itsFitChans.row(i)(0)
418 0 : << ": " << itsFitChans.row(i)(1) << "~" << itsFitChans.row(i)(2)
419 0 : << LogIO::POST;
420 : }
421 :
422 : // Use taQL for selection (no DDiD sel in MSSelection)
423 0 : ostringstream select;
424 0 : select << "select from $1 where ANTENNA1!=ANTENNA2";
425 0 : if (itsFieldIds.nelements() > 0) {
426 0 : select << " && FIELD_ID IN [" << itsFieldIds(0);
427 0 : for (uInt j=1; j<itsFieldIds.nelements(); j++)
428 0 : select << ", " << itsFieldIds(j);
429 0 : select << "]";
430 : }
431 0 : if (itsDDIds.nelements() > 0) {
432 0 : select << " && DATA_DESC_ID IN [" << itsDDIds(0);
433 0 : for (uInt j=1; j<itsDDIds.nelements(); j++)
434 0 : select << ", " << itsDDIds(j);
435 0 : select << "]";
436 : }
437 0 : MeasurementSet selectedMS(tableCommand(select, *ms_p).table());
438 0 : MSColumns msc(selectedMS);
439 0 : MSSelection mssel; // for channel and poln selection
440 :
441 0 : for (uInt iDD=0; iDD<itsDDIds.nelements(); iDD++) {
442 0 : Int thisDDId(itsDDIds(iDD));
443 0 : Int nChan = msc.spectralWindow().numChan()
444 0 : (msc.dataDescription().spectralWindowId()(thisDDId));
445 0 : Vector<Int> corrTypes = msc.polarization().corrType()
446 0 : (msc.dataDescription().polarizationId()(thisDDId));
447 :
448 : // default to fit all channels
449 0 : Vector<Bool> fitChanMask(nChan,True);
450 : // Handle non-trivial channel selection:
451 0 : if (itsFitChans.nelements()>0 && anyEQ(itsFitChans.column(0),thisDDId)) {
452 : // If subset of channels selected, set fit mask all False...
453 0 : fitChanMask = False;
454 0 : IPosition blc(1,0);
455 0 : IPosition trc(1,0);
456 : // ... and set only selected fit channels True:
457 0 : for (uInt i=0; i<itsFitChans.nrow(); ++i) {
458 0 : Vector<Int> chansel(itsFitChans.row(i)); // [spw, startChan, endChan]
459 : // match current spwId/DDid
460 0 : if (chansel(0) == thisDDId) {
461 0 : blc(0) = chansel(1);
462 0 : trc(0) = chansel(2);
463 0 : fitChanMask(blc,trc) = True;
464 : }
465 : }
466 : }
467 :
468 : // default to subtract from all channels
469 0 : Vector<Bool> subChanMask(nChan,True);
470 : // Handle non-trivial channel selection:
471 0 : if (itsSubChans.nelements()>0 && anyEQ(itsSubChans.column(0),thisDDId)) {
472 : // If subset of channels selected, set sub mask all False...
473 0 : subChanMask = False;
474 0 : IPosition blc(1,0);
475 0 : IPosition trc(1,0);
476 : // ... and set only selected sub channels True:
477 0 : for (uInt i=0; i<itsSubChans.nrow(); ++i) {
478 0 : Vector<Int> chansel(itsSubChans.row(i));
479 : // match current spwId/DDid
480 0 : if (chansel(0) == thisDDId) {
481 0 : blc(0) = chansel(1);
482 0 : trc(0) = chansel(2);
483 0 : subChanMask(blc,trc) = True;
484 : }
485 : }
486 : }
487 :
488 : // select parallel hand polarizations with MSSelection
489 0 : Vector<String> polSel(corrTypes.nelements());
490 0 : Int nPol = 0;
491 0 : for (uInt j=0; j<corrTypes.nelements(); j++) {
492 0 : if (corrTypes(j)==Stokes::XX || corrTypes(j)==Stokes::YY ||
493 0 : corrTypes(j)==Stokes::RR || corrTypes(j)==Stokes::LL) {
494 0 : polSel(nPol++) = Stokes::name(Stokes::type(corrTypes(j)));
495 : }
496 : }
497 0 : polSel.resize(nPol,True);
498 0 : String polnExpr("");
499 0 : if (nPol > 0) {
500 0 : polnExpr = polSel[0];
501 0 : for (Int np=1; np<nPol; ++np) {
502 0 : polnExpr += "," + polSel[np];
503 : }
504 : }
505 : // apply selection to get correlation slices
506 0 : String spwExpr = String::toString(thisDDId);
507 0 : MeasurementSet spwSelMS(selectedMS);
508 0 : Vector<Vector<Slice> > chansel;
509 0 : Vector<Vector<Slice> > corrsel;
510 0 : mssel.clear();
511 0 : mssSetData(selectedMS, spwSelMS, chansel, corrsel,
512 : "", "", "", "", spwExpr, "", "",
513 : polnExpr, "", "", "", "", 1, &mssel);
514 :
515 0 : casacore::Float chunkInt(itsSolInt);
516 0 : if (chunkInt < 0.0) {
517 0 : MSMainColumns msmain(spwSelMS);
518 0 : chunkInt = msmain.interval().getColumn()(0);
519 : }
520 :
521 : // Set up FrequencySelection with chansel and corrsel for VisibilityIterator2
522 0 : vi::FrequencySelectionUsingChannels freqSel = vi::FrequencySelectionUsingChannels();
523 0 : freqSel.add(mssel, &spwSelMS);
524 0 : freqSel.addCorrelationSlices(corrsel);
525 :
526 : // Sort Columns
527 0 : Block<Int> columns(5);
528 0 : Int i(0);
529 0 : columns[i++] = MS::ARRAY_ID;
530 0 : columns[i++] = MS::DATA_DESC_ID;
531 0 : columns[i++] = MS::SCAN_NUMBER;
532 0 : columns[i++] = MS::FIELD_ID;
533 0 : columns[i++] = MS::TIME;
534 0 : vi::SortColumns sortcols(columns, false);
535 :
536 : // Set up averaged non-calibrating vi
537 0 : vi::IteratingParameters iterpar(chunkInt, sortcols);
538 : vi::AveragingOptions avgopt(vi::AveragingOptions::AverageCorrected |
539 0 : vi::AveragingOptions::CorrectedWeightAvgFromWEIGHT);
540 0 : vi::AveragingParameters avgpar(1e8, 0.0, sortcols, avgopt);
541 0 : vi::LayeredVi2Factory factory(&spwSelMS, &iterpar, &avgpar);
542 0 : vi::VisibilityIterator2* vi2 = new vi::VisibilityIterator2(factory);
543 0 : vi2->setFrequencySelection(freqSel);
544 0 : vi::VisBuffer2* vb2 = vi2->getVisBuffer();
545 0 : vi2->originChunks();
546 0 : vi2->origin();
547 :
548 : // Get continuum model
549 : // fit
550 0 : Vector<Float> x(nChan);
551 0 : for (Int i=0; i<nChan; i++) x(i)=i;
552 0 : LinearFit<Float> fitter;
553 0 : Polynomial<AutoDiff<Float> > apoly(itsOrder);
554 0 : fitter.setFunction(apoly);
555 0 : Polynomial<Float> poly(itsOrder);
556 0 : Vector<Float> y1(nChan),y2(nChan);
557 0 : Vector<Float> sol(itsOrder+1);
558 0 : Vector<Complex> tmp(nChan);
559 0 : Cube<Complex> avgCorrData;
560 : // map is <scan, fit model>
561 0 : map<Int, Cube<Complex>> fitCorrData;
562 0 : map<Int, Vector<Int>> fitCorrIfr;
563 0 : Int nIter(0);
564 0 : IPosition start(3,0,0,0), end(3,0,nChan-1,0);
565 0 : for (vi2->originChunks(); vi2->moreChunks(); vi2->nextChunk()){
566 0 : for (vi2->origin(); vi2->more(); vi2->next()){
567 0 : avgCorrData.resize();
568 0 : avgCorrData = vb2->visCubeCorrected();
569 0 : for (start[2]=end[2]=0; start[2]<avgCorrData.shape()[2];
570 0 : start[2]++,end[2]++) {
571 0 : for (start[0]=end[0]=0; start[0]<nPol; start[0]++,end[0]++) {
572 0 : Vector<Complex> c(avgCorrData(start,end).nonDegenerate());
573 0 : tmp=c; // copy into contiguous storage
574 0 : c.set(Complex(0.0)); // zero the input
575 0 : real(y1,tmp);
576 0 : imag(y2,tmp);
577 0 : sol = fitter.fit(x,y1,&fitChanMask);
578 0 : poly.setCoefficients(sol);
579 0 : y1=x;
580 0 : y1.apply(poly);
581 0 : sol = fitter.fit(x,y2,&fitChanMask);
582 0 : poly.setCoefficients(sol);
583 0 : y2=x;
584 0 : y2.apply(poly);
585 0 : for (Int chn=0; chn<nChan; chn++) {
586 0 : if (subChanMask(chn))
587 0 : c(chn)=Complex(y1(chn),y2(chn));
588 : }
589 : }
590 : }
591 0 : fitCorrData[nIter] = avgCorrData;
592 0 : Vector<Int> ant1 = vb2->antenna1();
593 0 : Vector<Int> ant2 = vb2->antenna2();
594 0 : fitCorrIfr[nIter] = ant1*1000 + ant2;
595 0 : nIter++;
596 : }
597 : }
598 0 : os << "Processing " << nIter << " slots for ddId " << iDD << LogIO::POST;
599 0 : delete vi2;
600 :
601 : // unaveraged non-calibrating vi
602 0 : vi::LayeredVi2Factory factory2(&spwSelMS, &iterpar);
603 0 : vi2 = new vi::VisibilityIterator2(factory2);
604 0 : vi2->setFrequencySelection(freqSel);
605 0 : vb2 = vi2->getVisBuffer();
606 0 : vi2->originChunks();
607 0 : vi2->origin();
608 0 : Cube<Complex> modelData, correctedData, fitModel;
609 0 : nIter=0;
610 0 : for (vi2->originChunks(); vi2->moreChunks(); vi2->nextChunk()){
611 0 : for (vi2->origin(); vi2->more(); vi2->next()){
612 : // get model and corrected data with poln selection applied
613 0 : modelData.resize();
614 0 : modelData = vb2->visCubeModel();
615 0 : correctedData.resize();
616 0 : correctedData = vb2->visCubeCorrected();
617 0 : IPosition datashape = correctedData.shape();
618 0 : fitModel.resize();
619 0 : fitModel = fitCorrData[nIter];
620 :
621 : // Unaveraged data may have different shape
622 0 : if (fitModel.shape() != datashape) {
623 : // only use baselines in this chunk
624 0 : Cube<Complex> newFitModel(datashape);
625 0 : Vector<Int> ant1 = vb2->antenna1();
626 0 : Vector<Int> ant2 = vb2->antenna2();
627 0 : Vector<Int> ifrs = ant1*1000 + ant2;
628 0 : Vector<Int> fitIfr = fitCorrIfr[nIter];
629 0 : IPosition start(3,0,0,0), start2(3,0,0,0);
630 0 : Int nPol(datashape[0]), nChan(datashape[1]), nIfr(datashape[2]);
631 0 : IPosition end(3, nPol-1, nChan-1, nIfr-1);
632 0 : IPosition end2(3, nPol-1, nChan-1, nIfr-1);
633 0 : for (uInt i=0; i<ifrs.size(); ++i) {
634 0 : start[2] = end[2] = i;
635 0 : Int ifr = ifrs[i];
636 : // now find ifr in fitIfr
637 0 : for (uInt j=0; j<fitIfr.size(); ++j) {
638 0 : if (fitIfr[j]==ifr) {
639 0 : start2[2] = end2[2] = j;
640 0 : newFitModel(start,end) = fitModel(start2,end2);
641 0 : break;
642 : }
643 : }
644 : }
645 0 : fitModel.resize(datashape);
646 0 : fitModel = newFitModel;
647 : }
648 :
649 : // write out cubes to MS
650 0 : if (itsMode=="subtract" || itsMode=="model" || itsMode=="sub") {
651 0 : modelData = fitModel;
652 0 : vi2->getImpl()->writeVisModel(modelData);
653 : }
654 0 : if (itsMode=="replace") {
655 0 : correctedData = fitModel;
656 0 : vi2->getImpl()->writeVisCorrected(correctedData);
657 : }
658 0 : if (itsMode=="subtract" || itsMode=="sub") {
659 0 : correctedData -= fitModel;
660 0 : vi2->getImpl()->writeVisCorrected(correctedData);
661 : }
662 : }
663 0 : nIter++;
664 : }
665 0 : delete vi2;
666 : } // end DDId loop
667 0 : }
668 :
669 :
670 : } //# NAMESPACE CASA - END
671 :
|
