Line data Source code
1 : //# VisSetUtil.cc: VisSet Utilities
2 : //# Copyright (C) 1996,1997,1998,1999,2001,2002
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 adressed 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 : //#
27 : //# $Id$
28 :
29 : #include <casacore/casa/aips.h>
30 :
31 : #include <casacore/casa/Arrays/Matrix.h>
32 : #include <casacore/casa/Arrays/MatrixMath.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/Cube.h>
35 : #include <casacore/casa/BasicMath/Math.h>
36 : #include <casacore/casa/BasicSL/Complex.h>
37 : #include <casacore/casa/BasicSL/Constants.h>
38 : #include <casacore/casa/Utilities/Assert.h>
39 :
40 : #include <casacore/tables/Tables/ArrColDesc.h>
41 : #include <casacore/tables/Tables/ScaColDesc.h>
42 : #include <casacore/tables/Tables/TableCopy.h>
43 : #include <casacore/tables/DataMan/TiledDataStMan.h>
44 : #include <casacore/tables/DataMan/TiledShapeStMan.h>
45 : #include <casacore/tables/DataMan/StandardStMan.h>
46 : #include <casacore/tables/DataMan/TiledDataStManAccessor.h>
47 : #include <casacore/tables/DataMan/CompressComplex.h>
48 : #include <casacore/tables/DataMan/CompressFloat.h>
49 :
50 :
51 : #include <casacore/ms/MeasurementSets/MSColumns.h>
52 : #include <msvis/MSVis/VisModelDataI.h>
53 :
54 : #include <msvis/MSVis/VisSet.h>
55 : #include <msvis/MSVis/VisBuffer.h>
56 : #include <msvis/MSVis/VisSetUtil.h>
57 :
58 : #include <casacore/casa/Quanta/UnitMap.h>
59 : #include <casacore/casa/Quanta/UnitVal.h>
60 : #include <casacore/measures/Measures/Stokes.h>
61 : #include <casacore/casa/Quanta/MVAngle.h>
62 :
63 : #include <casacore/casa/Logging/LogIO.h>
64 :
65 : #include <iostream>
66 :
67 : using std::vector;
68 : using std::pair;
69 :
70 : using namespace casacore;
71 : namespace casa { //# NAMESPACE CASA - BEGIN
72 :
73 : // <summary>
74 : // </summary>
75 :
76 : // <reviewed reviewer="" date="" tests="tMEGI" demos="">
77 :
78 : // <prerequisite>
79 : // </prerequisite>
80 : //
81 : // <etymology>
82 : // </etymology>
83 : //
84 : // <synopsis>
85 : // </synopsis>
86 : //
87 : // <example>
88 : // <srcblock>
89 : // </srcblock>
90 : // </example>
91 : //
92 : // <motivation>
93 : // </motivation>
94 : //
95 : // <todo asof="">
96 : // </todo>
97 : // Calculate sensitivity
98 0 : void VisSetUtil::Sensitivity(VisSet &vs, Matrix<Double>& mssFreqSel,
99 : Matrix<Int>& mssChanSel,
100 : Quantity& pointsourcesens,
101 : Double& relativesens,
102 : Double& sumwt,
103 : Double& effectiveBandwidth,
104 : Double& effectiveIntegration,
105 : Int& nBaselines,
106 : Vector<Vector<Int> >& nData,
107 : Vector<Vector<Double> >& sumwtChan,
108 : Vector<Vector<Double> >& sumwtsqChan,
109 : Vector<Vector<Double> >& sumInverseVarianceChan)
110 : {
111 0 : ROVisIter& vi(vs.iter());
112 0 : VisSetUtil::Sensitivity(vi, mssFreqSel, mssChanSel, pointsourcesens, relativesens,
113 : sumwt, effectiveBandwidth, effectiveIntegration, nBaselines,
114 : nData, sumwtChan, sumwtsqChan, sumInverseVarianceChan);
115 :
116 0 : }
117 0 : void VisSetUtil::Sensitivity(ROVisIter &vi, Matrix<Double>& /*mssFreqSel*/,
118 : Matrix<Int>& mssChanSel,
119 : Quantity& pointsourcesens,
120 : Double& relativesens,
121 : Double& sumwt,
122 : Double& effectiveBandwidth,
123 : Double& effectiveIntegration,
124 : Int& nBaselines,
125 : Vector<Vector<Int> >& nData,
126 : Vector<Vector<Double> >& sumwtChan,
127 : Vector<Vector<Double> >& sumwtsqChan,
128 : Vector<Vector<Double> >& sumInverseVarianceChan)
129 : {
130 0 : LogIO os(LogOrigin("VisSetUtil", "Sensitivity()", WHERE));
131 :
132 0 : sumwt=0.0;
133 0 : Vector<Double> timeInterval, chanWidth;
134 0 : Double sumwtsq=0.0;
135 0 : Double sumInverseVariance=0.0;
136 : // Vector<Vector<Double> > sumwtChan, sumwtsqChan, sumInverseVarianceChan;
137 0 : VisBuffer vb(vi);
138 0 : Int nd=0,totalRows=0,spw=0;
139 0 : nBaselines=0;
140 0 : effectiveBandwidth=effectiveIntegration=0.0;
141 0 : Vector<Double> bwList;
142 0 : bwList.resize(mssChanSel.shape()(0)); bwList=0.0;
143 :
144 0 : sumwtChan.resize(mssChanSel.shape()(0));
145 0 : sumwtsqChan.resize(mssChanSel.shape()(0));
146 0 : sumInverseVarianceChan.resize(mssChanSel.shape()(0));
147 0 : nData.resize(mssChanSel.shape()(0));
148 0 : for (spw=0; spw<mssChanSel.shape()(0); spw++)
149 : {
150 0 : Int nc=mssChanSel(spw,2)-mssChanSel(spw,1)+1;
151 0 : for (int c=0; c<nc; c++)
152 : {
153 0 : sumwtChan(spw).resize(nc); sumwtChan(spw)=0.0;
154 0 : sumwtsqChan(spw).resize(nc); sumwtsqChan(spw) = 0.0;
155 0 : sumInverseVarianceChan(spw).resize(nc); sumInverseVarianceChan(spw)=0.0;
156 0 : nData(spw).resize(nc); nData(spw)=0;
157 : }
158 : }
159 : // sumwtChan=0.0;
160 : // sumwtsqChan=0.0;
161 : // sumInverseVarianceChan=0.0;
162 :
163 : // Now iterate through the data
164 0 : Int nant=vi.msColumns().antenna().nrow();
165 0 : Matrix<Int> baselines(nant,nant); baselines=0;
166 0 : Vector<Int> a1,a2;
167 0 : Vector<Double> t0, spwIntegTime;
168 0 : t0.resize(mssChanSel.shape()(0));
169 0 : spwIntegTime.resize(mssChanSel.shape()(0));
170 0 : t0 = spwIntegTime = 0.0;
171 :
172 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk())
173 : {
174 0 : for (vi.origin();vi.more();vi++)
175 : {
176 0 : Int nRow=vb.nRow();
177 0 : vb.nChannel();
178 : Int spwIndex;
179 0 : spw=vb.spectralWindow();
180 0 : spwIndex=-1;
181 0 : for (Int i=0;i<mssChanSel.shape()(0); i++)
182 0 : if (mssChanSel(i,0) == spw) {spwIndex=i; break;}
183 0 : if (spwIndex == -1)
184 : os << "Internal error: Could not locate the SPW index in the list of selected SPWs."
185 0 : << LogIO::EXCEPTION;
186 :
187 0 : timeInterval.assign(vb.timeInterval());
188 0 : Vector<Bool> rowFlags=vb.flagRow();
189 0 : Matrix<Bool> flag = vb.flag();
190 : // cerr << "SPW shape = " << vb.msColumns().spectralWindow().chanWidth().shape(spw) << " ";
191 0 : chanWidth.assign(vb.msColumns().spectralWindow().chanWidth().get(spw));
192 0 : a1.assign(vb.antenna1()); a2.assign(vb.antenna2());
193 :
194 0 : for (Int row=0; row<nRow; row++)
195 : {
196 : // TBD: Should probably use weight() here, which updates with calibration
197 0 : if (!rowFlags(row))
198 : {
199 : // Double variance=square(vb.sigma()(row));
200 0 : Double variance=1.0/(vb.weight()(row));
201 0 : if (abs(vb.time()(row) - t0(spwIndex)) > timeInterval(row))
202 : {
203 0 : t0(spwIndex)=vb.time()(row);
204 0 : spwIntegTime(spwIndex) += timeInterval(row);
205 : }
206 0 : totalRows++;
207 0 : baselines(a1(row), a2(row))++;
208 0 : for (Int chn=mssChanSel(spwIndex,1); chn<mssChanSel(spwIndex,2); chn+=mssChanSel(spwIndex,3))
209 : {
210 0 : if(!flag(chn,row)&&(variance>0.0))
211 : {
212 0 : sumwt+=vb.imagingWeight()(chn,row)*variance;
213 0 : sumwtsq+=square(vb.imagingWeight()(chn,row)*variance);
214 : // sumwtsq+=square(vb.imagingWeight()(chn,row));
215 0 : sumInverseVariance+=1.0/variance;
216 :
217 0 : sumwtChan(spwIndex)(chn) += vb.imagingWeight()(chn,row);
218 : //sumwtsqChan(spwIndex)(chn)+=square(vb.imagingWeight()(chn,row))*variance;
219 0 : sumwtsqChan(spwIndex)(chn)+=square(vb.imagingWeight()(chn,row));
220 :
221 0 : bwList(spwIndex) += abs(chanWidth(chn));
222 0 : (nData(spwIndex)(chn))++;
223 0 : nd++;
224 : }
225 : }
226 : }
227 : }
228 : }
229 : }
230 :
231 0 : if (totalRows == 0)
232 0 : os << "Cannot calculate sensitivty: No unflagged rows found" << LogIO::EXCEPTION;
233 :
234 0 : for (uInt spwndx=0; spwndx<bwList.nelements(); spwndx++)
235 : {
236 : // cerr << spwndx << " " << bwList(spwndx) << " " << nData(spwndx) << endl;
237 : // cerr << spwndx << " " << spwIntegTime(spwndx) << " " << endl;
238 0 : spw=mssChanSel(spwndx,0); // Get the SPW ID;
239 0 : chanWidth.assign(vi.msColumns().spectralWindow().chanWidth().get(spw)); // Extract the list of chan widths
240 :
241 0 : Int nchan=nData(spwndx).nelements();
242 : // If a channel from the current SPW was used, uses its width
243 : // for effective bandwidth calculation.
244 0 : for (Int j=0; j<nchan; j++)
245 0 : if (nData(spwndx)(j) > 0)
246 0 : effectiveBandwidth += abs(chanWidth(j));
247 0 : effectiveIntegration += spwIntegTime(spwndx)/bwList.nelements();
248 : }
249 :
250 0 : for (Int i=0; i<nant; i++)
251 0 : for (Int j=0; j<nant; j++)
252 0 : if (baselines(i,j) > 0)
253 0 : nBaselines++;
254 :
255 :
256 0 : if(sumwt==0.0) {
257 : os << "Cannot calculate sensitivity: sum of weights is zero" << endl
258 0 : << "Perhaps you need to weight the data" << LogIO::EXCEPTION;
259 : }
260 0 : if(sumInverseVariance==0.0) {
261 : os << "Cannot calculate sensitivity: sum of inverse variances is zero" << endl
262 0 : << "Perhaps you need to weight the data" << LogIO::EXCEPTION;
263 : }
264 :
265 : // Double naturalsens=1.0/sqrt(sumInverseVariance);
266 : // pointsourcesens=Quantity(sqrt(sumwtsq)/sumwt, "Jy");
267 : // relativesens=sqrt(sumwtsq)/sumwt/naturalsens;
268 :
269 : // cerr << "sumwt, sumwtsq, nd: " << sumwt << " " << sumwtsq << " " << nd << endl;
270 0 : relativesens=sqrt(nd*sumwtsq)/sumwt;
271 : //Double naturalsens=1.0/sqrt(sumInverseVariance);
272 0 : Double KB=1.3806488e-23;
273 0 : pointsourcesens=Quantity((10e26*2*KB*relativesens/sqrt(nBaselines*effectiveIntegration*effectiveBandwidth)), "Jy m^2/K");
274 :
275 : // cerr << "Pt src. = " << pointsourcesens << " " << integTime << " " << bandWidth/totalRows << endl;
276 : // cerr << baselines << endl;
277 0 : }
278 2 : void VisSetUtil::HanningSmooth(VisSet &vs, const String& dataCol, const Bool& doFlagAndWeight)
279 : {
280 2 : VisIter& vi(vs.iter());
281 2 : VisSetUtil::HanningSmooth(vi, dataCol, doFlagAndWeight);
282 2 : }
283 2 : void VisSetUtil::HanningSmooth(VisIter &vi, const String& dataCol, const Bool& doFlagAndWeight)
284 : {
285 : using casacore::operator*;
286 :
287 6 : LogIO os(LogOrigin("VisSetUtil", "HanningSmooth()"));
288 :
289 4 : VisBuffer vb(vi);
290 : Int row, chn, pol;
291 :
292 8 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
293 6 : if (vi.existsWeightSpectrum()) {
294 126 : for (vi.origin();vi.more();vi++) {
295 :
296 120 : Cube<Complex>& vc = ( dataCol=="data" ? vb.visCube() : vb.correctedVisCube());
297 :
298 120 : Cube<Bool>& fc= vb.flagCube();
299 120 : Cube<Float>& wc= vb.weightSpectrum();
300 :
301 120 : Int nRow=vb.nRow();
302 120 : Int nChan=vb.nChannel();
303 120 : if (nChan < 3) break;
304 120 : Int nPol=vi.visibilityShape()(0);
305 240 : Cube<Complex> smoothedData(nPol,nChan,nRow);
306 240 : Cube<Bool> newFlag(nPol,nChan,nRow);
307 240 : Cube<Float> newWeight(nPol,nChan,nRow);
308 45426 : for (row=0; row<nRow; row++) {
309 135918 : for (pol=0; pol<nPol; pol++) {
310 : ///Handle first channel and flag it
311 90612 : smoothedData(pol,0,row) = vc(pol,0,row)*0.5 + vc(pol,1,row)*0.5;
312 90612 : newWeight(pol,0,row) = 0.0;
313 90612 : newFlag(pol,0,row) = true;
314 5617944 : for (chn=1; chn<nChan-1; chn++) {
315 5527332 : smoothedData(pol,chn,row) =
316 5527332 : vc(pol,chn-1,row)*0.25 + vc(pol,chn,row)*0.50 +
317 11054664 : vc(pol,chn+1,row)*0.25;
318 11054664 : if (wc(pol,chn-1,row) != 0 && wc(pol,chn,row) != 0
319 11054664 : && wc(pol,chn+1,row) != 0) {
320 5527332 : newWeight(pol,chn,row) = 1.0 /
321 5527332 : (1.0/(wc(pol,chn-1,row)*16.0) + 1.0/(wc(pol,chn,row)*4.0)
322 5527332 : + 1.0/(wc(pol,chn+1,row)*16.0));
323 : } else {
324 0 : newWeight(pol,chn,row) = 0.0;
325 : }
326 5527332 : newFlag(pol,chn,row) =
327 5527332 : fc(pol,chn-1,row)||fc(pol,chn,row)||fc(pol,chn+1,row);
328 : }
329 : //Handle last channel and flag it
330 90612 : smoothedData(pol,nChan-1,row) =
331 181224 : vc(pol,nChan-2,row)*0.5+vc(pol,nChan-1,row)*0.5;
332 90612 : newWeight(pol,nChan-1,row) = 0.0;
333 90612 : newFlag(pol,nChan-1,row) = true; // flag last channel
334 : }
335 : }
336 :
337 120 : if(dataCol=="data"){
338 60 : if(doFlagAndWeight){
339 60 : vi.setVisAndFlag(smoothedData,newFlag,VisibilityIterator::Observed);
340 60 : vi.setWeightSpectrum(newWeight);
341 : }
342 : else{
343 0 : vi.setVis(smoothedData,VisibilityIterator::Observed);
344 : }
345 : }
346 : else{
347 60 : if(doFlagAndWeight){
348 0 : vi.setVisAndFlag(smoothedData,newFlag,VisibilityIterator::Corrected);
349 0 : vi.setWeightSpectrum(newWeight);
350 : }
351 : else{
352 60 : vi.setVis(smoothedData,VisibilityIterator::Corrected);
353 : }
354 : }
355 :
356 : }
357 : } else {
358 0 : for (vi.origin();vi.more();vi++) {
359 :
360 0 : Cube<Complex>& vc = (dataCol=="data" ? vb.visCube() : vb.correctedVisCube());
361 :
362 0 : Cube<Bool>& fc= vb.flagCube();
363 0 : Matrix<Float>& wm = vb.weightMat();
364 :
365 0 : Int nRow=vb.nRow();
366 0 : Int nChan=vb.nChannel();
367 0 : if (nChan < 3) break;
368 0 : Int nPol=vi.visibilityShape()(0);
369 0 : Cube<Complex> smoothedData(nPol,nChan,nRow);
370 0 : Cube<Bool> newFlag(nPol,nChan,nRow);
371 0 : Matrix<Float> newWeight(nPol, nRow);
372 0 : for (row=0; row<nRow; row++) {
373 0 : for (pol=0; pol<nPol; pol++) {
374 : ///Handle first channel and flag it
375 0 : smoothedData(pol,0,row) = vc(pol,0,row)*0.5 + vc(pol,1,row)*0.5;
376 0 : newFlag(pol,0,row) = true;
377 : ///Handle chan-independent weights
378 0 : newWeight(pol, row) = 8.0*wm(pol, row)/3.0;
379 0 : for (chn=1; chn<nChan-1; chn++) {
380 0 : smoothedData(pol,chn,row) =
381 0 : vc(pol,chn-1,row)*0.25 + vc(pol,chn,row)*0.50 +
382 0 : vc(pol,chn+1,row)*0.25;
383 0 : newFlag(pol,chn,row) =
384 0 : fc(pol,chn-1,row)||fc(pol,chn,row)||fc(pol,chn+1,row);
385 : }
386 : //Handle last channel and flag it
387 0 : smoothedData(pol,nChan-1,row) =
388 0 : vc(pol,nChan-2,row)*0.5+vc(pol,nChan-1,row)*0.5;
389 0 : newFlag(pol,nChan-1,row) = true; // flag last channel
390 : }
391 : }
392 :
393 0 : if(dataCol=="data"){
394 0 : if(doFlagAndWeight){
395 0 : vi.setVisAndFlag(smoothedData,newFlag,VisibilityIterator::Observed);
396 0 : vi.setWeightMat(newWeight);
397 : }
398 : else{
399 0 : vi.setVis(smoothedData,VisibilityIterator::Observed);
400 : }
401 : }
402 : else{
403 0 : if(doFlagAndWeight){
404 0 : vi.setVisAndFlag(smoothedData,newFlag,VisibilityIterator::Corrected);
405 0 : vi.setWeightMat(newWeight);
406 : }
407 : else{
408 0 : vi.setVis(smoothedData,VisibilityIterator::Corrected);
409 : }
410 : }
411 : }
412 : }
413 : }
414 2 : }
415 :
416 :
417 393 : void VisSetUtil::addScrCols(MeasurementSet& ms, Bool addModel, Bool addCorr,
418 : Bool alsoinit, Bool compress) {
419 :
420 : // Add but DO NOT INITIALIZE calibration set (comprising a set of CORRECTED_DATA
421 : // and MODEL_DATA columns) to the MeasurementSet.
422 :
423 : // Sense if anything is to be done
424 393 : addModel=addModel && !(ms.tableDesc().isColumn("MODEL_DATA"));
425 393 : addCorr=addCorr && !(ms.tableDesc().isColumn("CORRECTED_DATA"));
426 :
427 : // Escape (silently) without doing anything if nothing
428 : // to be done
429 393 : if (!addModel && !addCorr)
430 220 : return;
431 : else {
432 : // Remove SORTED_TABLE and continue
433 : // This is needed because old SORTED_TABLE won't see
434 : // the added column(s)
435 173 : if (ms.keywordSet().isDefined("SORT_COLUMNS"))
436 53 : ms.rwKeywordSet().removeField("SORT_COLUMNS");
437 173 : if (ms.keywordSet().isDefined("SORTED_TABLE"))
438 45 : ms.rwKeywordSet().removeField("SORTED_TABLE");
439 : }
440 :
441 : // If we are adding the MODEL_DATA column, make it
442 : // the exclusive origin for model visibilities by
443 : // deleting any OTF model keywords
444 173 : if (addModel)
445 102 : VisSetUtil::remOTFModel(ms);
446 :
447 :
448 : // Form log message
449 346 : String addMessage("Adding ");
450 173 : if (addModel) {
451 102 : addMessage+="MODEL_DATA ";
452 102 : if (addCorr) addMessage+="and ";
453 : }
454 173 : if (addCorr) addMessage+="CORRECTED_DATA ";
455 173 : addMessage+="column(s).";
456 346 : LogSink logSink;
457 519 : LogMessage message(addMessage,LogOrigin("VisSetUtil","addScrCols"));
458 173 : logSink.post(message);
459 :
460 : {
461 : // Define a column accessor to the observed data
462 : TableColumn* data;
463 173 : const bool data_is_float = ms.tableDesc().isColumn(MS::columnName(MS::FLOAT_DATA));
464 173 : if (data_is_float) {
465 12 : data = new ArrayColumn<Float> (ms, MS::columnName(MS::FLOAT_DATA));
466 : } else {
467 161 : data = new ArrayColumn<Complex> (ms, MS::columnName(MS::DATA));
468 : };
469 :
470 : // Check if the data column is tiled and, if so, the
471 : // smallest tile shape used.
472 346 : TableDesc td = ms.actualTableDesc();
473 173 : const ColumnDesc& cdesc = td[data->columnDesc().name()];
474 346 : String dataManType = cdesc.dataManagerType();
475 346 : String dataManGroup = cdesc.dataManagerGroup();
476 :
477 346 : IPosition dataTileShape;
478 173 : Bool tiled = (dataManType.contains("Tiled"));
479 173 : Bool simpleTiling = false;
480 :
481 173 : if (!tiled || !simpleTiling) {
482 : // Untiled, or tiled at a higher than expected dimensionality
483 : // Use a canonical tile shape of 1 MB size
484 :
485 173 : MSSpWindowColumns msspwcol(ms.spectralWindow());
486 173 : Int maxNchan = max (msspwcol.numChan().getColumn());
487 173 : Int tileSize = maxNchan/10 + 1;
488 173 : Int nCorr = data->shape(0)(0);
489 173 : dataTileShape = IPosition(3, nCorr, tileSize, 131072/nCorr/tileSize + 1);
490 : };
491 :
492 173 : delete data;
493 :
494 173 : if(addModel){
495 : // Add the MODEL_DATA column
496 204 : TableDesc tdModel, tdModelComp, tdModelScale;
497 102 : CompressComplex* ccModel=NULL;
498 204 : String colModel=MS::columnName(MS::MODEL_DATA);
499 :
500 102 : tdModel.addColumn(ArrayColumnDesc<Complex>(colModel,"model data", 2));
501 102 : td.addColumn(ArrayColumnDesc<Complex>(colModel,"model data", 2));
502 204 : IPosition modelTileShape = dataTileShape;
503 102 : if (compress) {
504 0 : tdModelComp.addColumn(ArrayColumnDesc<Int>(colModel+"_COMPRESSED",
505 0 : "model data compressed",2));
506 0 : tdModelScale.addColumn(ScalarColumnDesc<Float>(colModel+"_SCALE"));
507 0 : tdModelScale.addColumn(ScalarColumnDesc<Float>(colModel+"_OFFSET"));
508 0 : ccModel = new CompressComplex(colModel, colModel+"_COMPRESSED",
509 0 : colModel+"_SCALE", colModel+"_OFFSET", true);
510 :
511 0 : StandardStMan modelScaleStMan("ModelScaleOffset");
512 0 : ms.addColumn(tdModelScale, modelScaleStMan);
513 :
514 :
515 0 : TiledShapeStMan modelCompStMan("ModelCompTiled", modelTileShape);
516 0 : ms.addColumn(tdModelComp, modelCompStMan);
517 0 : ms.addColumn(tdModel, *ccModel);
518 :
519 : } else {
520 :
521 : // Make a clone of the appropriate data column. The clone can differ by
522 : // by having a different element datatype.
523 :
524 : try {
525 10 : String dataColumnName = data_is_float ? MS::columnName (MS::FLOAT_DATA)
526 113 : : MS::columnName (MS::DATA);
527 :
528 104 : TableCopy::cloneColumnTyped<Complex> (ms, dataColumnName, ms, MS::columnName(MS::MODEL_DATA),
529 : "TiledModel");
530 :
531 : // Now fill the model column with complex zeroes. Using the function below
532 : // to do the filling preserves the tiling of the hypercubes making up
533 : // the column. The later call to initScrCols will set these appropriately
534 : // since depending on the way the correlations are represented, only one
535 : // of the correlations will be set to one. Unfortunately, this results in
536 : // two writes of the scratch column where one would do.
537 :
538 101 : TableCopy::fillColumnData (ms, MS::columnName (MS::MODEL_DATA),
539 101 : Complex (0.0, 0.0),
540 : ms, dataColumnName,
541 101 : true); // last arg preserves tile shape
542 2 : } catch (AipsError & e){
543 :
544 : // Column cloning failed (presumably). Try it the old way.
545 :
546 1 : MeasurementSet::addColumnToDesc(tdModel, MeasurementSet::MODEL_DATA, 2);
547 2 : TiledShapeStMan modelStMan("ModelTiled", modelTileShape);
548 1 : ms.addColumn(tdModel, modelStMan);
549 : }
550 : }
551 :
552 102 : if (ccModel) delete ccModel;
553 : }
554 173 : if (addCorr) {
555 : // Add the CORRECTED_DATA column
556 206 : TableDesc tdCorr, tdCorrComp, tdCorrScale;
557 103 : CompressComplex* ccCorr=NULL;
558 206 : String colCorr=MS::columnName(MS::CORRECTED_DATA);
559 :
560 103 : tdCorr.addColumn(ArrayColumnDesc<Complex>(colCorr,"corrected data", 2));
561 206 : IPosition corrTileShape = dataTileShape;
562 103 : if (compress) {
563 0 : tdCorrComp.addColumn(ArrayColumnDesc<Int>(colCorr+"_COMPRESSED",
564 0 : "corrected data compressed",2));
565 0 : tdCorrScale.addColumn(ScalarColumnDesc<Float>(colCorr+"_SCALE"));
566 0 : tdCorrScale.addColumn(ScalarColumnDesc<Float>(colCorr+"_OFFSET"));
567 0 : ccCorr = new CompressComplex(colCorr, colCorr+"_COMPRESSED",
568 0 : colCorr+"_SCALE", colCorr+"_OFFSET", true);
569 :
570 0 : StandardStMan corrScaleStMan("CorrScaleOffset");
571 0 : ms.addColumn(tdCorrScale, corrScaleStMan);
572 :
573 0 : TiledShapeStMan corrCompStMan("CorrectedCompTiled", corrTileShape);
574 0 : ms.addColumn(tdCorrComp, corrCompStMan);
575 0 : ms.addColumn(tdCorr, *ccCorr);
576 :
577 : } else {
578 :
579 : try {
580 :
581 : // Make a clone of the appropriate data column, except the new column
582 : // will have datatype complex.
583 :
584 3 : String dataColumnName = data_is_float ? MS::columnName (MS::FLOAT_DATA)
585 107 : : MS::columnName (MS::DATA);
586 :
587 105 : TableCopy::cloneColumnTyped<Complex> (ms, dataColumnName,
588 : ms, colCorr,
589 : "TiledCorrected");
590 :
591 : // Copy the appropriate data column into the new one, preserving the tiling
592 : // of the column's hypercubes.
593 :
594 102 : TableCopy::copyColumnData (ms, dataColumnName,
595 : ms, colCorr,
596 : true); // last arg preserves tile shape
597 :
598 102 : addCorr = false; // We've already done the copying
599 :
600 2 : } catch (AipsError & e){
601 :
602 : // Table clone failed (presumably): try it the old way
603 :
604 2 : TiledShapeStMan corrStMan("CorrectedTiled", corrTileShape);
605 1 : ms.addColumn(tdCorr, corrStMan);
606 : }
607 : // Copy the column keywords
608 3 : String dataColumnName = data_is_float ? MS::columnName (MS::FLOAT_DATA)
609 209 : : MS::columnName (MS::DATA);
610 103 : const TableRecord& dataColKeyWord = ms.tableDesc().columnDesc(dataColumnName).keywordSet();
611 103 : if (dataColKeyWord.nfields() > 0) {
612 12 : LogMessage message("Start copying column keyword(s) of "+colCorr+" from "+dataColumnName,LogOrigin("VisSetUtil","addScrCols"));
613 4 : logSink.post(message);
614 8 : ArrayColumn<Complex> corrCol(ms, colCorr);
615 4 : TableRecord &corrColKeyWord = corrCol.rwKeywordSet();
616 4 : if (corrColKeyWord.nfields()==0) {
617 4 : corrColKeyWord.assign(dataColKeyWord);
618 : } else {
619 0 : corrColKeyWord.merge(dataColKeyWord, RecordInterface::OverwriteDuplicates);
620 : }
621 : }
622 : }
623 :
624 103 : if (ccCorr) delete ccCorr;
625 : }
626 173 : ms.flush();
627 :
628 : }
629 :
630 173 : if (alsoinit)
631 : // Initialize only what we added
632 173 : VisSetUtil::initScrCols(ms,addModel,addCorr);
633 :
634 173 : return;
635 :
636 :
637 : }
638 :
639 0 : void VisSetUtil::remScrCols(MeasurementSet& ms, Bool remModel, Bool remCorr) {
640 :
641 0 : Vector<String> colNames(2);
642 0 : colNames(0)=MS::columnName(MS::MODEL_DATA);
643 0 : colNames(1)=MS::columnName(MS::CORRECTED_DATA);
644 :
645 0 : Vector<Bool> doCol(2);
646 0 : doCol(0)=remModel;
647 0 : doCol(1)=remCorr;
648 :
649 :
650 0 : for (uInt j=0; j<colNames.nelements(); j++) {
651 0 : if (doCol(j)) {
652 0 : if (ms.tableDesc().isColumn(colNames(j))) {
653 0 : ms.removeColumn(colNames(j));
654 : };
655 0 : if (ms.tableDesc().isColumn(colNames(j)+"_COMPRESSED")) {
656 0 : ms.removeColumn(colNames(j)+"_COMPRESSED");
657 : };
658 0 : if (ms.tableDesc().isColumn(colNames(j)+"_SCALE")) {
659 0 : ms.removeColumn(colNames(j)+"_SCALE");
660 : };
661 0 : if (ms.tableDesc().isColumn(colNames(j)+"_OFFSET")) {
662 0 : ms.removeColumn(colNames(j)+"_OFFSET");
663 : };
664 : };
665 : };
666 :
667 :
668 0 : }
669 :
670 102 : void VisSetUtil::remOTFModel(MeasurementSet& ms) {
671 :
672 204 : CountedPtr<VisModelDataI> visModelData = VisModelDataI::create();
673 102 : visModelData->clearModelI(ms);
674 102 : }
675 :
676 173 : void VisSetUtil::initScrCols(MeasurementSet& ms, Bool initModel, Bool initCorr) {
677 :
678 : // Sense if anything is to be done (relevant scr cols must be present)
679 173 : initModel=initModel && ms.tableDesc().isColumn("MODEL_DATA");
680 173 : initCorr=initCorr && ms.tableDesc().isColumn("CORRECTED_DATA");
681 :
682 : // Do nothing?
683 173 : if (!initModel && !initCorr)
684 71 : return;
685 :
686 : /* Reconsider this trap?
687 : // Trap missing columns:
688 : if (initModel && !ms.tableDesc().isColumn("MODEL_DATA"))
689 : throw(AipsError("Cannot initialize MODEL_DATA if column is absent."));
690 : if (initCorr && !ms.tableDesc().isColumn("CORRECTED_DATA"))
691 : throw(AipsError("Cannot initialize CORRECTED_DATA if column is absent."));
692 : */
693 :
694 : // Create ordinary (un-row-selected) VisibilityIterator from the MS
695 204 : VisibilityIterator vi(ms,Block<Int>(),0.0);
696 :
697 : // Pass to VisibilityIterator-oriented method
698 102 : VisSetUtil::initScrCols(vi,initModel,initCorr);
699 : }
700 :
701 :
702 102 : void VisSetUtil::initScrCols(VisibilityIterator& vi, Bool initModel, Bool initCorr) {
703 :
704 : // Sense if anything is to be done (relevant scr cols must be present)
705 102 : initModel=initModel && vi.ms().tableDesc().isColumn("MODEL_DATA");
706 102 : initCorr=initCorr && vi.ms().tableDesc().isColumn("CORRECTED_DATA");
707 :
708 : // Do nothing?
709 102 : if (!initModel && !initCorr)
710 0 : return;
711 :
712 : /* Reconsider this trap?
713 : // Trap missing columns:
714 : if (initModel && !vi.ms().tableDesc().isColumn("MODEL_DATA"))
715 : throw(AipsError("Cannot initialize MODEL_DATA if column is absent."));
716 : if (initCorr && !vi.ms().tableDesc().isColumn("CORRECTED_DATA"))
717 : throw(AipsError("Cannot initialize CORRECTED_DATA if column is absent."));
718 : */
719 :
720 : // Form log message
721 204 : String initMessage("Initializing ");
722 102 : if (initModel) {
723 102 : initMessage+="MODEL_DATA to (unity)";
724 102 : if (initCorr) initMessage+=" and ";
725 : }
726 102 : if (initCorr) initMessage+="CORRECTED_DATA (to DATA)";
727 102 : initMessage+=".";
728 204 : LogSink logSink;
729 306 : LogMessage message(initMessage,LogOrigin("VisSetUtil","initScrCols"));
730 102 : logSink.post(message);
731 :
732 204 : Vector<Int> lastCorrType;
733 204 : Vector<Bool> zero;
734 102 : Int nRows(0);
735 102 : vi.setRowBlocking(10000);
736 820 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
737 :
738 1436 : Vector<Int> corrType; vi.corrType(corrType);
739 718 : uInt nCorr = corrType.nelements();
740 718 : if (initModel) {
741 : // figure out which correlations to set to 1. and 0. for the model.
742 1334 : if (nCorr != lastCorrType.nelements() ||
743 616 : !allEQ(corrType, lastCorrType)) {
744 :
745 103 : lastCorrType.resize(nCorr);
746 103 : lastCorrType=corrType;
747 103 : zero.resize(nCorr);
748 :
749 317 : for (uInt i=0; i<nCorr; i++)
750 : {
751 630 : zero[i]=(corrType[i]==Stokes::RL || corrType[i]==Stokes::LR ||
752 416 : corrType[i]==Stokes::XY || corrType[i]==Stokes::YX);
753 : }
754 : }
755 : }
756 1980 : for (vi.origin(); vi.more(); vi++) {
757 1262 : nRows+=vi.nRow();
758 :
759 2524 : Cube<Complex> data;
760 1262 : vi.visibility(data, VisibilityIterator::Observed);
761 1262 : if (initCorr) {
762 : // Read DATA and set CORRECTED_DATA
763 16 : vi.setVis(data,VisibilityIterator::Corrected);
764 : }
765 1262 : if (initModel) {
766 : // Set MODEL_DATA
767 1262 : data.set(1.0);
768 1262 : if (ntrue(zero)>0) {
769 455 : for (uInt i=0; i < nCorr; i++) {
770 364 : if (zero[i]) data(Slice(i), Slice(), Slice()) = Complex(0.0,0.0);
771 : }
772 : }
773 1262 : vi.setVis(data,VisibilityIterator::Model);
774 : }
775 : }
776 : }
777 102 : vi.ms().relinquishAutoLocks();
778 :
779 102 : vi.originChunks();
780 102 : vi.setRowBlocking(0);
781 :
782 204 : ostringstream os;
783 102 : os << "Initialized " << nRows << " rows.";
784 102 : message.message(os.str());
785 102 : logSink.post(message);
786 :
787 : }
788 :
789 0 : void VisSetUtil::removeCalSet(MeasurementSet& ms, Bool removeModel) {
790 : // Remove an existing calibration set (comprising a set of CORRECTED_DATA
791 : // and MODEL_DATA columns) from the MeasurementSet.
792 :
793 : //Remove model in header
794 0 : if(removeModel)
795 0 : VisSetUtil::remOTFModel(ms);
796 :
797 0 : VisSetUtil::remScrCols(ms,true,true);
798 :
799 0 : }
800 :
801 :
802 2 : void VisSetUtil::UVSub(VisSet &vs, Bool reverse)
803 : {
804 2 : VisIter& vi(vs.iter());
805 2 : VisSetUtil::UVSub(vi, reverse);
806 2 : }
807 2 : void VisSetUtil::UVSub(VisIter &vi, Bool reverse)
808 : {
809 6 : LogIO os(LogOrigin("VisSetUtil", "UVSub()"));
810 :
811 :
812 4 : VisBuffer vb(vi);
813 :
814 : Int row, chn, pol;
815 :
816 4 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
817 362 : for (vi.origin();vi.more();vi++) {
818 :
819 360 : Cube<Complex>& vc= vb.correctedVisCube();
820 360 : Cube<Complex>& mc= vb.modelVisCube();
821 :
822 360 : Int nRow=vb.nRow();
823 360 : Int nChan=vb.nChannel();
824 360 : Int nPol=vi.visibilityShape()(0);
825 720 : Cube<Complex> residualData(nPol,nChan,nRow);
826 360 : if (reverse) {
827 0 : for (row=0; row<nRow; row++) {
828 0 : for (pol=0; pol<nPol; pol++) {
829 0 : for (chn=0; chn<nChan; chn++) {
830 0 : residualData(pol,chn,row) = vc(pol,chn,row)+mc(pol,chn ,row);
831 : }
832 : }
833 : }
834 : } else {
835 126720 : for (row=0; row<nRow; row++) {
836 379080 : for (pol=0; pol<nPol; pol++) {
837 758160 : for (chn=0; chn<nChan; chn++) {
838 505440 : residualData(pol,chn,row) = vc(pol,chn,row)-mc(pol,chn ,row);
839 : }
840 : }
841 : }
842 : }
843 360 : vi.setVis(residualData,VisibilityIterator::Corrected);
844 : }
845 : }
846 2 : }
847 :
848 : } //# NAMESPACE CASA - END
849 :
|
