source: trunk/src/SDMath.cc@ 256

Last change on this file since 256 was 248, checked in by kil064, 20 years ago

remove function 'quotient' and put its functionality into
function 'binaryOperate'

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 38.9 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDMath.cc: A collection of single dish mathematical operations
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2004
5//# ATNF
6//#
7//# This program is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU General Public License as published by the Free
9//# Software Foundation; either version 2 of the License, or (at your option)
10//# any later version.
11//#
12//# This program is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of
14//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15//# Public License for more details.
16//#
17//# You should have received a copy of the GNU General Public License along
18//# with this program; if not, write to the Free Software Foundation, Inc.,
19//# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: Malte.Marquarding@csiro.au
23//# Postal address: Malte Marquarding,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
29//# $Id:
30//#---------------------------------------------------------------------------
31#include <vector>
32
33#include <casa/aips.h>
34#include <casa/BasicSL/String.h>
35#include <casa/Arrays/IPosition.h>
36#include <casa/Arrays/Array.h>
37#include <casa/Arrays/ArrayIter.h>
38#include <casa/Arrays/VectorIter.h>
39#include <casa/Arrays/ArrayMath.h>
40#include <casa/Arrays/ArrayLogical.h>
41#include <casa/Arrays/MaskedArray.h>
42#include <casa/Arrays/MaskArrMath.h>
43#include <casa/Arrays/MaskArrLogi.h>
44#include <casa/BasicMath/Math.h>
45#include <casa/Containers/Block.h>
46#include <casa/Quanta/QC.h>
47#include <casa/Utilities/Assert.h>
48#include <casa/Exceptions.h>
49
50#include <scimath/Mathematics/VectorKernel.h>
51#include <scimath/Mathematics/Convolver.h>
52#include <scimath/Mathematics/InterpolateArray1D.h>
53#include <scimath/Functionals/Polynomial.h>
54
55#include <tables/Tables/Table.h>
56#include <tables/Tables/ScalarColumn.h>
57#include <tables/Tables/ArrayColumn.h>
58#include <tables/Tables/ReadAsciiTable.h>
59
60#include <lattices/Lattices/LatticeUtilities.h>
61#include <lattices/Lattices/RebinLattice.h>
62#include <coordinates/Coordinates/SpectralCoordinate.h>
63#include <coordinates/Coordinates/CoordinateSystem.h>
64#include <coordinates/Coordinates/CoordinateUtil.h>
65#include <coordinates/Coordinates/VelocityAligner.h>
66
67#include "MathUtils.h"
68#include "SDDefs.h"
69#include "SDContainer.h"
70#include "SDMemTable.h"
71
72#include "SDMath.h"
73
74using namespace casa;
75using namespace asap;
76
77
78SDMath::SDMath()
79{;}
80
81SDMath::SDMath(const SDMath& other)
82{
83
84// No state
85
86}
87
88SDMath& SDMath::operator=(const SDMath& other)
89{
90 if (this != &other) {
91// No state
92 }
93 return *this;
94}
95
96SDMath::~SDMath()
97{;}
98
99
100CountedPtr<SDMemTable> SDMath::average(const Block<CountedPtr<SDMemTable> >& in,
101 const Vector<Bool>& mask, Bool scanAv,
102 const String& weightStr, Bool alignVelocity) const
103//
104// Weighted averaging of spectra from one or more Tables.
105//
106{
107
108// Convert weight type
109
110 WeightType wtType = NONE;
111 convertWeightString(wtType, weightStr);
112
113// Create output Table by cloning from the first table
114
115 SDMemTable* pTabOut = new SDMemTable(*in[0],True);
116
117// Setup
118
119 IPosition shp = in[0]->rowAsMaskedArray(0).shape(); // Must not change
120 Array<Float> arr(shp);
121 Array<Bool> barr(shp);
122 const Bool useMask = (mask.nelements() == shp(asap::ChanAxis));
123
124// Columns from Tables
125
126 ROArrayColumn<Float> tSysCol;
127 ROScalarColumn<Double> mjdCol;
128 ROScalarColumn<String> srcNameCol;
129 ROScalarColumn<Double> intCol;
130 ROArrayColumn<uInt> fqIDCol;
131
132// Create accumulation MaskedArray. We accumulate for each channel,if,pol,beam
133// Note that the mask of the accumulation array will ALWAYS remain ALL True.
134// The MA is only used so that when data which is masked Bad is added to it,
135// that data does not contribute.
136
137 Array<Float> zero(shp);
138 zero=0.0;
139 Array<Bool> good(shp);
140 good = True;
141 MaskedArray<Float> sum(zero,good);
142
143// Counter arrays
144
145 Array<Float> nPts(shp); // Number of points
146 nPts = 0.0;
147 Array<Float> nInc(shp); // Increment
148 nInc = 1.0;
149
150// Create accumulation Array for variance. We accumulate for
151// each if,pol,beam, but average over channel. So we need
152// a shape with one less axis dropping channels.
153
154 const uInt nAxesSub = shp.nelements() - 1;
155 IPosition shp2(nAxesSub);
156 for (uInt i=0,j=0; i<(nAxesSub+1); i++) {
157 if (i!=asap::ChanAxis) {
158 shp2(j) = shp(i);
159 j++;
160 }
161 }
162 Array<Float> sumSq(shp2);
163 sumSq = 0.0;
164 IPosition pos2(nAxesSub,0); // For indexing
165
166// Time-related accumulators
167
168 Double time;
169 Double timeSum = 0.0;
170 Double intSum = 0.0;
171 Double interval = 0.0;
172
173// To get the right shape for the Tsys accumulator we need to
174// access a column from the first table. The shape of this
175// array must not change
176
177 Array<Float> tSysSum;
178 {
179 const Table& tabIn = in[0]->table();
180 tSysCol.attach(tabIn,"TSYS");
181 tSysSum.resize(tSysCol.shape(0));
182 }
183 tSysSum =0.0;
184 Array<Float> tSys;
185
186// Scan and row tracking
187
188 Int oldScanID = 0;
189 Int outScanID = 0;
190 Int scanID = 0;
191 Int rowStart = 0;
192 Int nAccum = 0;
193 Int tableStart = 0;
194
195// Source and FreqID
196
197 String sourceName, oldSourceName, sourceNameStart;
198 Vector<uInt> freqID, freqIDStart, oldFreqID;
199
200// Velocity Aligner. We need an aligner for each Direction and FreqID
201// combination. I don't think there is anyway to know how many
202// directions there are.
203// For now, assume all Tables have the same Frequency Table
204
205/*
206 {
207 MEpoch::Ref timeRef(MEpoch::UTC); // Should be in header
208 MDirection::Types dirRef(MDirection::J2000); // Should be in header
209//
210 SDHeader sh = in[0].getSDHeader();
211 const uInt nChan = sh.nchan;
212//
213 const SDFrequencyTable freqTab = in[0]->getSDFreqTable();
214 const uInt nFreqID = freqTab.length();
215 PtrBlock<const VelocityAligner<Float>* > vA(nFreqID);
216
217// Get first time from first table
218
219 const Table& tabIn0 = in[0]->table();
220 mjdCol.attach(tabIn0, "TIME");
221 Double dTmp;
222 mjdCol.get(0, dTmp);
223 MVEpoch tmp2(Quantum<Double>(dTmp, Unit(String("d"))));
224 MEpoch epoch(tmp2, timeRef);
225//
226 for (uInt freqID=0; freqID<nFreqID; freqID++) {
227 SpectralCoordinate sC = in[0]->getCoordinate(freqID);
228 vA[freqID] = new VelocityAligner<Float>(sC, nChan, epoch, const MDirection& dir,
229 const MPosition& pos, const String& velUnit,
230 MDoppler::Types velType, MFrequency::Types velFreqSystem)
231 }
232 }
233*/
234
235// Loop over tables
236
237 Float fac = 1.0;
238 const uInt nTables = in.nelements();
239 for (uInt iTab=0; iTab<nTables; iTab++) {
240
241// Should check that the frequency tables don't change if doing VelocityAlignment
242
243// Attach columns to Table
244
245 const Table& tabIn = in[iTab]->table();
246 tSysCol.attach(tabIn, "TSYS");
247 mjdCol.attach(tabIn, "TIME");
248 srcNameCol.attach(tabIn, "SRCNAME");
249 intCol.attach(tabIn, "INTERVAL");
250 fqIDCol.attach(tabIn, "FREQID");
251
252// Loop over rows in Table
253
254 const uInt nRows = in[iTab]->nRow();
255 for (uInt iRow=0; iRow<nRows; iRow++) {
256
257// Check conformance
258
259 IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
260 if (!shp.isEqual(shp2)) {
261 throw (AipsError("Shapes for all rows must be the same"));
262 }
263
264// If we are not doing scan averages, make checks for source and
265// frequency setup and warn if averaging across them
266
267// Get copy of Scan Container for this row
268
269 SDContainer sc = in[iTab]->getSDContainer(iRow);
270 scanID = sc.scanid;
271
272// Get quantities from columns
273
274 srcNameCol.getScalar(iRow, sourceName);
275 mjdCol.get(iRow, time);
276 tSysCol.get(iRow, tSys);
277 intCol.get(iRow, interval);
278 fqIDCol.get(iRow, freqID);
279
280// Initialize first source and freqID
281
282 if (iRow==0 && iTab==0) {
283 sourceNameStart = sourceName;
284 freqIDStart = freqID;
285 }
286
287// If we are doing scan averages, see if we are at the end of an
288// accumulation period (scan). We must check soutce names too,
289// since we might have two tables with one scan each but different
290// source names; we shouldn't average different sources together
291
292 if (scanAv && ( (scanID != oldScanID) ||
293 (iRow==0 && iTab>0 && sourceName!=oldSourceName))) {
294
295// Normalize data in 'sum' accumulation array according to weighting scheme
296
297 normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
298
299// Fill scan container. The source and freqID come from the
300// first row of the first table that went into this average (
301// should be the same for all rows in the scan average)
302
303 Float nR(nAccum);
304 fillSDC(sc, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
305 timeSum/nR, intSum, sourceNameStart, freqIDStart);
306
307// Write container out to Table
308
309 pTabOut->putSDContainer(sc);
310
311// Reset accumulators
312
313 sum = 0.0;
314 sumSq = 0.0;
315 nAccum = 0;
316//
317 tSysSum =0.0;
318 timeSum = 0.0;
319 intSum = 0.0;
320 nPts = 0.0;
321
322// Increment
323
324 rowStart = iRow; // First row for next accumulation
325 tableStart = iTab; // First table for next accumulation
326 sourceNameStart = sourceName; // First source name for next accumulation
327 freqIDStart = freqID; // First FreqID for next accumulation
328//
329 oldScanID = scanID;
330 outScanID += 1; // Scan ID for next accumulation period
331 }
332
333// Accumulate
334
335 accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum,
336 tSys, nInc, mask, time, interval, in, iTab, iRow, asap::ChanAxis,
337 nAxesSub, useMask, wtType);
338//
339 oldSourceName = sourceName;
340 oldFreqID = freqID;
341 }
342 }
343
344// OK at this point we have accumulation data which is either
345// - accumulated from all tables into one row
346// or
347// - accumulated from the last scan average
348//
349// Normalize data in 'sum' accumulation array according to weighting scheme
350 normalize(sum, sumSq, nPts, wtType, asap::ChanAxis, nAxesSub);
351
352// Create and fill container. The container we clone will be from
353// the last Table and the first row that went into the current
354// accumulation. It probably doesn't matter that much really...
355
356 Float nR(nAccum);
357 SDContainer sc = in[tableStart]->getSDContainer(rowStart);
358 fillSDC(sc, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
359 timeSum/nR, intSum, sourceNameStart, freqIDStart);
360 pTabOut->putSDContainer(sc);
361//
362 return CountedPtr<SDMemTable>(pTabOut);
363}
364
365
366
367CountedPtr<SDMemTable> SDMath::binaryOperate (const CountedPtr<SDMemTable>& left,
368 const CountedPtr<SDMemTable>& right,
369 const String& op, Bool preserve) const
370{
371
372// Check operator
373
374 String op2(op);
375 op2.upcase();
376 uInt what = 0;
377 if (op2=="ADD") {
378 what = 0;
379 } else if (op2=="SUB") {
380 what = 1;
381 } else if (op2=="MUL") {
382 what = 2;
383 } else if (op2=="DIV") {
384 what = 3;
385 } else if (op2=="QUOTIENT") {
386 what = 4;
387 } else {
388 throw( AipsError("Unrecognized operation"));
389 }
390
391// Check rows
392
393 const uInt nRowLeft = left->nRow();
394 const uInt nRowRight = right->nRow();
395 Bool ok = (nRowRight==1&&nRowLeft>0) ||
396 (nRowRight>=1&&nRowLeft==nRowRight);
397 if (!ok) {
398 throw (AipsError("The right Scan Table can have one row or the same number of rows as the left Scan Table"));
399 }
400
401// Input Tables
402
403 const Table& tLeft = left->table();
404 const Table& tRight = right->table();
405
406// TSys columns
407
408 ROArrayColumn<Float> tSysLeft(tLeft, "TSYS");
409 ROArrayColumn<Float> tSysRight(tRight, "TSYS");
410
411// First row for right
412
413 Array<Float> tSysLeftArr, tSysRightArr;
414 tSysRight.get(0, tSysRightArr);
415 MaskedArray<Float>* pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(0));
416 IPosition shpRight = pMRight->shape();
417
418// Output Table cloned from left
419
420 SDMemTable* pTabOut = new SDMemTable(*left, True);
421
422// Loop over rows
423
424 for (uInt i=0; i<nRowLeft; i++) {
425
426// Get data
427
428 MaskedArray<Float> mLeft(left->rowAsMaskedArray(i));
429 IPosition shpLeft = mLeft.shape();
430 tSysLeft.get(i, tSysLeftArr);
431//
432 if (nRowRight>1) {
433 delete pMRight;
434 pMRight = new MaskedArray<Float>(right->rowAsMaskedArray(i));
435 shpRight = pMRight->shape();
436 tSysRight.get(i, tSysRightArr);
437 }
438//
439 if (!shpRight.isEqual(shpLeft)) {
440 throw(AipsError("left and right scan tables are not conformant"));
441 }
442 if (!tSysRightArr.shape().isEqual(tSysRightArr.shape())) {
443 throw(AipsError("left and right Tsys data are not conformant"));
444 }
445 if (!shpRight.isEqual(tSysRightArr.shape())) {
446 throw(AipsError("left and right scan tables are not conformant"));
447 }
448
449// Make container
450
451 SDContainer sc = left->getSDContainer(i);
452
453// Operate on data and TSys
454
455 if (what==0) {
456 MaskedArray<Float> tmp = mLeft + *pMRight;
457 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
458 sc.putTsys(tSysLeftArr+tSysRightArr);
459 } else if (what==1) {
460 MaskedArray<Float> tmp = mLeft - *pMRight;
461 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
462 sc.putTsys(tSysLeftArr-tSysRightArr);
463 } else if (what==2) {
464 MaskedArray<Float> tmp = mLeft * *pMRight;
465 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
466 sc.putTsys(tSysLeftArr*tSysRightArr);
467 } else if (what==3) {
468 MaskedArray<Float> tmp = mLeft / *pMRight;
469 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
470 sc.putTsys(tSysLeftArr/tSysRightArr);
471 } else if (what==4) {
472 if (preserve) {
473 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - tSysRightArr;
474 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
475 } else {
476 MaskedArray<Float> tmp = (tSysRightArr * mLeft / *pMRight) - tSysLeftArr;
477 putDataInSDC(sc, tmp.getArray(), tmp.getMask());
478 }
479 sc.putTsys(tSysRightArr);
480 }
481
482// Put new row in output Table
483
484 pTabOut->putSDContainer(sc);
485 }
486 if (pMRight) delete pMRight;
487//
488 return CountedPtr<SDMemTable>(pTabOut);
489}
490
491
492
493std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in,
494 const Vector<Bool>& mask,
495 const String& which, Int row) const
496//
497// Perhaps iteration over pol/beam/if should be in here
498// and inside the nrow iteration ?
499//
500{
501 const uInt nRow = in->nRow();
502
503// Specify cursor location
504
505 IPosition start, end;
506 getCursorLocation(start, end, *in);
507
508// Loop over rows
509
510 const uInt nEl = mask.nelements();
511 uInt iStart = 0;
512 uInt iEnd = in->nRow()-1;
513//
514 if (row>=0) {
515 iStart = row;
516 iEnd = row;
517 }
518//
519 std::vector<float> result(iEnd-iStart+1);
520 for (uInt ii=iStart; ii <= iEnd; ++ii) {
521
522// Get row and deconstruct
523
524 MaskedArray<Float> marr(in->rowAsMaskedArray(ii));
525 Array<Float> arr = marr.getArray();
526 Array<Bool> barr = marr.getMask();
527
528// Access desired piece of data
529
530 Array<Float> v((arr(start,end)).nonDegenerate());
531 Array<Bool> m((barr(start,end)).nonDegenerate());
532
533// Apply OTF mask
534
535 MaskedArray<Float> tmp;
536 if (m.nelements()==nEl) {
537 tmp.setData(v,m&&mask);
538 } else {
539 tmp.setData(v,m);
540 }
541
542// Get statistic
543
544 result[ii-iStart] = mathutil::statistics(which, tmp);
545 }
546//
547 return result;
548}
549
550
551SDMemTable* SDMath::bin(const SDMemTable& in, Int width) const
552{
553 SDHeader sh = in.getSDHeader();
554 SDMemTable* pTabOut = new SDMemTable(in, True);
555
556// Bin up SpectralCoordinates
557
558 IPosition factors(1);
559 factors(0) = width;
560 for (uInt j=0; j<in.nCoordinates(); ++j) {
561 CoordinateSystem cSys;
562 cSys.addCoordinate(in.getCoordinate(j));
563 CoordinateSystem cSysBin =
564 CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False);
565//
566 SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
567 pTabOut->setCoordinate(sCBin, j);
568 }
569
570// Use RebinLattice to find shape
571
572 IPosition shapeIn(1,sh.nchan);
573 IPosition shapeOut = RebinLattice<Float>::rebinShape(shapeIn, factors);
574 sh.nchan = shapeOut(0);
575 pTabOut->putSDHeader(sh);
576
577
578// Loop over rows and bin along channel axis
579
580 for (uInt i=0; i < in.nRow(); ++i) {
581 SDContainer sc = in.getSDContainer(i);
582//
583 Array<Float> tSys(sc.getTsys()); // Get it out before sc changes shape
584
585// Bin up spectrum
586
587 MaskedArray<Float> marr(in.rowAsMaskedArray(i));
588 MaskedArray<Float> marrout;
589 LatticeUtilities::bin(marrout, marr, asap::ChanAxis, width);
590
591// Put back the binned data and flags
592
593 IPosition ip2 = marrout.shape();
594 sc.resize(ip2);
595//
596 putDataInSDC(sc, marrout.getArray(), marrout.getMask());
597
598// Bin up Tsys.
599
600 Array<Bool> allGood(tSys.shape(),True);
601 MaskedArray<Float> tSysIn(tSys, allGood, True);
602//
603 MaskedArray<Float> tSysOut;
604 LatticeUtilities::bin(tSysOut, tSysIn, asap::ChanAxis, width);
605 sc.putTsys(tSysOut.getArray());
606//
607 pTabOut->putSDContainer(sc);
608 }
609 return pTabOut;
610}
611
612SDMemTable* SDMath::unaryOperate(const SDMemTable& in, Float val, Bool doAll,
613 uInt what) const
614//
615// what = 0 Multiply
616// 1 Add
617{
618 SDMemTable* pOut = new SDMemTable(in,False);
619 const Table& tOut = pOut->table();
620 ArrayColumn<Float> spec(tOut,"SPECTRA");
621//
622 if (doAll) {
623 for (uInt i=0; i < tOut.nrow(); i++) {
624
625// Get
626
627 MaskedArray<Float> marr(pOut->rowAsMaskedArray(i));
628
629// Operate
630
631 if (what==0) {
632 marr *= val;
633 } else if (what==1) {
634 marr += val;
635 }
636
637// Put
638
639 spec.put(i, marr.getArray());
640 }
641 } else {
642
643// Get cursor location
644
645 IPosition start, end;
646 getCursorLocation(start, end, in);
647//
648 for (uInt i=0; i < tOut.nrow(); i++) {
649
650// Get
651
652 MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
653
654// Modify. More work than we would like to deal with the mask
655
656 Array<Float>& values = dataIn.getRWArray();
657 Array<Bool> mask(dataIn.getMask());
658//
659 Array<Float> values2 = values(start,end);
660 Array<Bool> mask2 = mask(start,end);
661 MaskedArray<Float> t(values2,mask2);
662 if (what==0) {
663 t *= val;
664 } else if (what==1) {
665 t += val;
666 }
667 values(start, end) = t.getArray(); // Write back into 'dataIn'
668
669// Put
670 spec.put(i, dataIn.getArray());
671 }
672 }
673//
674 return pOut;
675}
676
677
678
679SDMemTable* SDMath::averagePol(const SDMemTable& in, const Vector<Bool>& mask) const
680//
681// Average all polarizations together, weighted by variance
682//
683{
684// WeightType wtType = NONE;
685// convertWeightString(wtType, weight);
686
687 const uInt nRows = in.nRow();
688 const uInt polAxis = asap::PolAxis; // Polarization axis
689 const uInt chanAxis = asap::ChanAxis; // Spectrum axis
690
691// Create output Table and reshape number of polarizations
692
693 Bool clear=True;
694 SDMemTable* pTabOut = new SDMemTable(in, clear);
695 SDHeader header = pTabOut->getSDHeader();
696 header.npol = 1;
697 pTabOut->putSDHeader(header);
698
699// Shape of input and output data
700
701 const IPosition& shapeIn = in.rowAsMaskedArray(0u, False).shape();
702 IPosition shapeOut(shapeIn);
703 shapeOut(polAxis) = 1; // Average all polarizations
704//
705 const uInt nChan = shapeIn(chanAxis);
706 const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape
707 IPosition start(4), end(4);
708
709// Output arrays
710
711 Array<Float> outData(shapeOut, 0.0);
712 Array<Bool> outMask(shapeOut, True);
713 const IPosition axes(2, 2, 3); // pol-channel plane
714//
715 const Bool useMask = (mask.nelements() == shapeIn(chanAxis));
716
717// Loop over rows
718
719 for (uInt iRow=0; iRow<nRows; iRow++) {
720
721// Get data for this row
722
723 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
724 Array<Float>& arr = marr.getRWArray();
725 const Array<Bool>& barr = marr.getMask();
726
727// Make iterators to iterate by pol-channel planes
728
729 ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
730 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
731
732// Accumulations
733
734 Float fac = 1.0;
735 Vector<Float> vecSum(nChan,0.0);
736
737// Iterate through data by pol-channel planes
738
739 while (!itDataPlane.pastEnd()) {
740
741// Iterate through plane by polarization and accumulate Vectors
742
743 Vector<Float> t1(nChan); t1 = 0.0;
744 Vector<Bool> t2(nChan); t2 = True;
745 MaskedArray<Float> vecSum(t1,t2);
746 Float varSum = 0.0;
747 {
748 ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
749 ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
750 while (!itDataVec.pastEnd()) {
751
752// Create MA of data & mask (optionally including OTF mask) and get variance
753
754 if (useMask) {
755 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
756 fac = 1.0 / variance(spec);
757 } else {
758 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
759 fac = 1.0 / variance(spec);
760 }
761
762// Normalize spectrum (without OTF mask) and accumulate
763
764 const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
765 vecSum += spec;
766 varSum += fac;
767
768// Next
769
770 itDataVec.next();
771 itMaskVec.next();
772 }
773 }
774
775// Normalize summed spectrum
776
777 vecSum /= varSum;
778
779// FInd position in input data array. We are iterating by pol-channel
780// plane so all that will change is beam and IF and that's what we want.
781
782 IPosition pos = itDataPlane.pos();
783
784// Write out data. This is a bit messy. We have to reform the Vector
785// accumulator into an Array of shape (1,1,1,nChan)
786
787 start = pos;
788 end = pos;
789 end(chanAxis) = nChan-1;
790 outData(start,end) = vecSum.getArray().reform(vecShapeOut);
791 outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
792
793// Step to next beam/IF combination
794
795 itDataPlane.next();
796 itMaskPlane.next();
797 }
798
799// Generate output container and write it to output table
800
801 SDContainer sc = in.getSDContainer();
802 sc.resize(shapeOut);
803//
804 putDataInSDC(sc, outData, outMask);
805 pTabOut->putSDContainer(sc);
806 }
807//
808 return pTabOut;
809}
810
811
812SDMemTable* SDMath::smooth(const SDMemTable& in,
813 const casa::String& kernelType,
814 casa::Float width, Bool doAll) const
815{
816
817// Number of channels
818
819 const uInt chanAxis = asap::ChanAxis; // Spectral axis
820 SDHeader sh = in.getSDHeader();
821 const uInt nChan = sh.nchan;
822
823// Generate Kernel
824
825 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
826 Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
827
828// Generate Convolver
829
830 IPosition shape(1,nChan);
831 Convolver<Float> conv(kernel, shape);
832
833// New Table
834
835 SDMemTable* pTabOut = new SDMemTable(in,True);
836
837// Get cursor location
838
839 IPosition start, end;
840 getCursorLocation(start, end, in);
841//
842 IPosition shapeOut(4,1);
843
844// Output Vectors
845
846 Vector<Float> valuesOut(nChan);
847 Vector<Bool> maskOut(nChan);
848
849// Loop over rows in Table
850
851 for (uInt ri=0; ri < in.nRow(); ++ri) {
852
853// Get copy of data
854
855 const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(ri));
856 AlwaysAssert(dataIn.shape()(chanAxis)==nChan, AipsError);
857//
858 Array<Float> valuesIn = dataIn.getArray();
859 Array<Bool> maskIn = dataIn.getMask();
860
861// Branch depending on whether we smooth all locations or just
862// those pointed at by the current selection cursor
863
864 if (doAll) {
865 uInt axis = asap::ChanAxis;
866 VectorIterator<Float> itValues(valuesIn, axis);
867 VectorIterator<Bool> itMask(maskIn, axis);
868 while (!itValues.pastEnd()) {
869
870// Smooth
871 if (kernelType==VectorKernel::HANNING) {
872 mathutil::hanning(valuesOut, maskOut, itValues.vector(), itMask.vector());
873 itMask.vector() = maskOut;
874 } else {
875 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
876 conv.linearConv(valuesOut, itValues.vector());
877 }
878//
879 itValues.vector() = valuesOut;
880//
881 itValues.next();
882 itMask.next();
883 }
884 } else {
885
886// Set multi-dim Vector shape
887
888 shapeOut(chanAxis) = valuesIn.shape()(chanAxis);
889
890// Stuff about with shapes so that we don't have conformance run-time errors
891
892 Vector<Float> valuesIn2 = valuesIn(start,end).nonDegenerate();
893 Vector<Bool> maskIn2 = maskIn(start,end).nonDegenerate();
894
895// Smooth
896
897 if (kernelType==VectorKernel::HANNING) {
898 mathutil::hanning(valuesOut, maskOut, valuesIn2, maskIn2);
899 maskIn(start,end) = maskOut.reform(shapeOut);
900 } else {
901 mathutil::replaceMaskByZero(valuesIn2, maskIn2);
902 conv.linearConv(valuesOut, valuesIn2);
903 }
904//
905 valuesIn(start,end) = valuesOut.reform(shapeOut);
906 }
907
908// Create and put back
909
910 SDContainer sc = in.getSDContainer(ri);
911 putDataInSDC(sc, valuesIn, maskIn);
912//
913 pTabOut->putSDContainer(sc);
914 }
915//
916 return pTabOut;
917}
918
919
920SDMemTable* SDMath::convertFlux (const SDMemTable& in, Float a, Float eta, Bool doAll) const
921//
922// As it is, this function could be implemented with 'simpleOperate'
923// However, I anticipate that eventually we will look the conversion
924// values up in a Table and apply them in a frequency dependent way,
925// so I have implemented it fully here
926//
927{
928 SDHeader sh = in.getSDHeader();
929 SDMemTable* pTabOut = new SDMemTable(in, True);
930
931// FInd out how to convert values into Jy and K (e.g. units might be mJy or mK)
932// Also automatically find out what we are converting to according to the
933// flux unit
934
935 Unit fluxUnit(sh.fluxunit);
936 Unit K(String("K"));
937 Unit JY(String("Jy"));
938//
939 Bool toKelvin = True;
940 Double inFac = 1.0;
941 if (fluxUnit==JY) {
942 cerr << "Converting to K" << endl;
943//
944 Quantum<Double> t(1.0,fluxUnit);
945 Quantum<Double> t2 = t.get(JY);
946 inFac = (t2 / t).getValue();
947//
948 toKelvin = True;
949 sh.fluxunit = "K";
950 } else if (fluxUnit==K) {
951 cerr << "Converting to Jy" << endl;
952//
953 Quantum<Double> t(1.0,fluxUnit);
954 Quantum<Double> t2 = t.get(K);
955 inFac = (t2 / t).getValue();
956//
957 toKelvin = False;
958 sh.fluxunit = "Jy";
959 } else {
960 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
961 }
962 pTabOut->putSDHeader(sh);
963
964// Compute conversion factor. 'a' and 'eta' are really frequency, time and
965// telescope dependent and should be looked// up in a table
966
967 Float factor = 2.0 * inFac * 1.0e-7 * 1.0e26 *
968 QC::k.getValue(Unit(String("erg/K"))) / a / eta;
969 if (toKelvin) {
970 factor = 1.0 / factor;
971 }
972 cerr << "Applying conversion factor = " << factor << endl;
973
974// For operations only on specified cursor location
975
976 IPosition start, end;
977 getCursorLocation(start, end, in);
978
979// Loop over rows and apply factor to spectra
980
981 const uInt axis = asap::ChanAxis;
982 for (uInt i=0; i < in.nRow(); ++i) {
983
984// Get data
985
986 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
987 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
988 const Array<Bool>& maskIn = dataIn.getMask();
989
990// Need to apply correct conversion factor (frequency and time dependent)
991// which should be sourced from a Table. For now we just apply the given
992// factor to everything
993
994 if (doAll) {
995 VectorIterator<Float> itValues(valuesIn, asap::ChanAxis);
996 while (!itValues.pastEnd()) {
997 itValues.vector() *= factor; // Writes back into dataIn
998//
999 itValues.next();
1000 }
1001 } else {
1002 Array<Float> valuesIn2 = valuesIn(start,end);
1003 valuesIn2 *= factor;
1004 valuesIn(start,end) = valuesIn2;
1005 }
1006
1007// Write out
1008
1009 SDContainer sc = in.getSDContainer(i);
1010 putDataInSDC(sc, valuesIn, maskIn);
1011//
1012 pTabOut->putSDContainer(sc);
1013 }
1014 return pTabOut;
1015}
1016
1017
1018
1019SDMemTable* SDMath::gainElevation (const SDMemTable& in, const Vector<Float>& coeffs,
1020 const String& fileName,
1021 const String& methodStr, Bool doAll) const
1022{
1023
1024// Get header and clone output table
1025
1026 SDHeader sh = in.getSDHeader();
1027 SDMemTable* pTabOut = new SDMemTable(in, True);
1028
1029// Get elevation data from SDMemTable and convert to degrees
1030
1031 const Table& tab = in.table();
1032 ROScalarColumn<Float> elev(tab, "ELEVATION");
1033 Vector<Float> x = elev.getColumn();
1034 x *= Float(180 / C::pi);
1035//
1036 const uInt nC = coeffs.nelements();
1037 if (fileName.length()>0 && nC>0) {
1038 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1039 }
1040
1041// Correct
1042
1043 if (nC>0 || fileName.length()==0) {
1044
1045// Find instrument
1046
1047 Bool throwIt = True;
1048 Instrument inst = SDMemTable::convertInstrument (sh.antennaname, throwIt);
1049
1050// Set polynomial
1051
1052 Polynomial<Float>* pPoly = 0;
1053 Vector<Float> coeff;
1054 String msg;
1055 if (nC>0) {
1056 pPoly = new Polynomial<Float>(nC);
1057 coeff = coeffs;
1058 msg = String("user");
1059 } else {
1060 if (inst==PKSMULTIBEAM) {
1061 } else if (inst==PKSSINGLEBEAM) {
1062 } else if (inst==TIDBINBILLA) {
1063 pPoly = new Polynomial<Float>(3);
1064 coeff.resize(3);
1065 coeff(0) = 3.58788e-1;
1066 coeff(1) = 2.87243e-2;
1067 coeff(2) = -3.219093e-4;
1068 } else if (inst==MOPRA) {
1069 }
1070 msg = String("built in");
1071 }
1072//
1073 if (coeff.nelements()>0) {
1074 pPoly->setCoefficients(coeff);
1075 } else {
1076 throw(AipsError("There is no known gain-el polynomial known for this instrument"));
1077 }
1078//
1079 cerr << "Making polynomial correction with " << msg << " coefficients" << endl;
1080 const uInt nRow = in.nRow();
1081 Vector<Float> factor(nRow);
1082 for (uInt i=0; i<nRow; i++) {
1083 factor[i] = (*pPoly)(x[i]);
1084 }
1085 delete pPoly;
1086//
1087 correctFromVector (pTabOut, in, doAll, factor);
1088 } else {
1089
1090// Indicate which columns to read from ascii file
1091
1092 String col0("ELEVATION");
1093 String col1("FACTOR");
1094
1095// Read and correct
1096
1097 cerr << "Making correction from ascii Table" << endl;
1098 correctFromAsciiTable (pTabOut, in, fileName, col0, col1,
1099 methodStr, doAll, x);
1100 }
1101//
1102 return pTabOut;
1103}
1104
1105
1106
1107SDMemTable* SDMath::opacity (const SDMemTable& in, Float tau, Bool doAll) const
1108{
1109
1110// Get header and clone output table
1111
1112 SDHeader sh = in.getSDHeader();
1113 SDMemTable* pTabOut = new SDMemTable(in, True);
1114
1115// Get elevation data from SDMemTable and convert to degrees
1116
1117 const Table& tab = in.table();
1118 ROScalarColumn<Float> elev(tab, "ELEVATION");
1119 Vector<Float> zDist = elev.getColumn();
1120 zDist = Float(C::pi_2) - zDist;
1121
1122// Generate correction factor
1123
1124 const uInt nRow = in.nRow();
1125 Vector<Float> factor(nRow);
1126 Vector<Float> factor2(nRow);
1127 for (uInt i=0; i<nRow; i++) {
1128 factor[i] = exp(tau)/cos(zDist[i]);
1129 }
1130
1131// Correct
1132
1133 correctFromVector (pTabOut, in, doAll, factor);
1134//
1135 return pTabOut;
1136}
1137
1138
1139
1140
1141// 'private' functions
1142
1143void SDMath::fillSDC(SDContainer& sc,
1144 const Array<Bool>& mask,
1145 const Array<Float>& data,
1146 const Array<Float>& tSys,
1147 Int scanID, Double timeStamp,
1148 Double interval, const String& sourceName,
1149 const Vector<uInt>& freqID) const
1150{
1151// Data and mask
1152
1153 putDataInSDC(sc, data, mask);
1154
1155// TSys
1156
1157 sc.putTsys(tSys);
1158
1159// Time things
1160
1161 sc.timestamp = timeStamp;
1162 sc.interval = interval;
1163 sc.scanid = scanID;
1164//
1165 sc.sourcename = sourceName;
1166 sc.putFreqMap(freqID);
1167}
1168
1169void SDMath::normalize(MaskedArray<Float>& sum,
1170 const Array<Float>& sumSq,
1171 const Array<Float>& nPts,
1172 WeightType wtType, Int axis,
1173 Int nAxesSub) const
1174{
1175 IPosition pos2(nAxesSub,0);
1176//
1177 if (wtType==NONE) {
1178
1179// We just average by the number of points accumulated.
1180// We need to make a MA out of nPts so that no divide by
1181// zeros occur
1182
1183 MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
1184 sum /= t;
1185 } else if (wtType==VAR) {
1186
1187// Normalize each spectrum by sum(1/var) where the variance
1188// is worked out for each spectrum
1189
1190 Array<Float>& data = sum.getRWArray();
1191 VectorIterator<Float> itData(data, axis);
1192 while (!itData.pastEnd()) {
1193 pos2 = itData.pos().getFirst(nAxesSub);
1194 itData.vector() /= sumSq(pos2);
1195 itData.next();
1196 }
1197 } else if (wtType==TSYS) {
1198 }
1199}
1200
1201
1202void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
1203 MaskedArray<Float>& sum, Array<Float>& sumSq,
1204 Array<Float>& nPts, Array<Float>& tSysSum,
1205 const Array<Float>& tSys, const Array<Float>& nInc,
1206 const Vector<Bool>& mask, Double time, Double interval,
1207 const Block<CountedPtr<SDMemTable> >& in,
1208 uInt iTab, uInt iRow, uInt axis,
1209 uInt nAxesSub, Bool useMask,
1210 WeightType wtType) const
1211{
1212
1213// Get data
1214
1215 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
1216 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
1217 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
1218//
1219 if (wtType==NONE) {
1220 const MaskedArray<Float> n(nInc,dataIn.getMask());
1221 nPts += n; // Only accumulates where mask==T
1222 } else if (wtType==VAR) {
1223
1224// We are going to average the data, weighted by the noise for each pol, beam and IF.
1225// So therefore we need to iterate through by spectrum (axis 3)
1226
1227 VectorIterator<Float> itData(valuesIn, axis);
1228 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
1229 Float fac = 1.0;
1230 IPosition pos(nAxesSub,0);
1231//
1232 while (!itData.pastEnd()) {
1233
1234// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
1235
1236 if (useMask) {
1237 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
1238 fac = 1.0/variance(tmp);
1239 } else {
1240 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
1241 fac = 1.0/variance(tmp);
1242 }
1243
1244// Scale data
1245
1246 itData.vector() *= fac; // Writes back into 'dataIn'
1247//
1248// Accumulate variance per if/pol/beam averaged over spectrum
1249// This method to get pos2 from itData.pos() is only valid
1250// because the spectral axis is the last one (so we can just
1251// copy the first nAXesSub positions out)
1252
1253 pos = itData.pos().getFirst(nAxesSub);
1254 sumSq(pos) += fac;
1255//
1256 itData.next();
1257 itMask.next();
1258 }
1259 } else if (wtType==TSYS) {
1260 }
1261
1262// Accumulate sum of (possibly scaled) data
1263
1264 sum += dataIn;
1265
1266// Accumulate Tsys, time, and interval
1267
1268 tSysSum += tSys;
1269 timeSum += time;
1270 intSum += interval;
1271 nAccum += 1;
1272}
1273
1274
1275
1276
1277void SDMath::getCursorLocation(IPosition& start, IPosition& end,
1278 const SDMemTable& in) const
1279{
1280 const uInt nDim = 4;
1281 const uInt i = in.getBeam();
1282 const uInt j = in.getIF();
1283 const uInt k = in.getPol();
1284 const uInt n = in.nChan();
1285//
1286 start.resize(nDim);
1287 start(0) = i;
1288 start(1) = j;
1289 start(2) = k;
1290 start(3) = 0;
1291//
1292 end.resize(nDim);
1293 end(0) = i;
1294 end(1) = j;
1295 end(2) = k;
1296 end(3) = n-1;
1297}
1298
1299
1300void SDMath::convertWeightString(WeightType& wtType, const String& weightStr) const
1301{
1302 String tStr(weightStr);
1303 tStr.upcase();
1304 if (tStr.contains(String("NONE"))) {
1305 wtType = NONE;
1306 } else if (tStr.contains(String("VAR"))) {
1307 wtType = VAR;
1308 } else if (tStr.contains(String("TSYS"))) {
1309 wtType = TSYS;
1310 throw(AipsError("T_sys weighting not yet implemented"));
1311 } else {
1312 throw(AipsError("Unrecognized weighting type"));
1313 }
1314}
1315
1316void SDMath::convertInterpString(Int& type, const String& interp) const
1317{
1318 String tStr(interp);
1319 tStr.upcase();
1320 if (tStr.contains(String("NEAR"))) {
1321 type = InterpolateArray1D<Float,Float>::nearestNeighbour;
1322 } else if (tStr.contains(String("LIN"))) {
1323 type = InterpolateArray1D<Float,Float>::linear;
1324 } else if (tStr.contains(String("CUB"))) {
1325 type = InterpolateArray1D<Float,Float>::cubic;
1326 } else if (tStr.contains(String("SPL"))) {
1327 type = InterpolateArray1D<Float,Float>::spline;
1328 } else {
1329 throw(AipsError("Unrecognized interpolation type"));
1330 }
1331}
1332
1333void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
1334 const Array<Bool>& mask) const
1335{
1336 sc.putSpectrum(data);
1337//
1338 Array<uChar> outflags(data.shape());
1339 convertArray(outflags,!mask);
1340 sc.putFlags(outflags);
1341}
1342
1343Table SDMath::readAsciiFile (const String& fileName) const
1344{
1345 String formatString;
1346 Table tbl = readAsciiTable (formatString, Table::Memory, fileName, "", "", False);
1347 return tbl;
1348}
1349
1350
1351
1352void SDMath::correctFromAsciiTable(SDMemTable* pTabOut,
1353 const SDMemTable& in, const String& fileName,
1354 const String& col0, const String& col1,
1355 const String& methodStr, Bool doAll,
1356 const Vector<Float>& xOut) const
1357{
1358
1359// Read gain-elevation ascii file data into a Table.
1360
1361 Table geTable = readAsciiFile (fileName);
1362//
1363 correctFromTable (pTabOut, in, geTable, col0, col1, methodStr, doAll, xOut);
1364}
1365
1366void SDMath::correctFromTable(SDMemTable* pTabOut, const SDMemTable& in,
1367 const Table& tTable, const String& col0,
1368 const String& col1,
1369 const String& methodStr, Bool doAll,
1370 const Vector<Float>& xOut) const
1371{
1372
1373// Get data from Table
1374
1375 ROScalarColumn<Float> geElCol(tTable, col0);
1376 ROScalarColumn<Float> geFacCol(tTable, col1);
1377 Vector<Float> xIn = geElCol.getColumn();
1378 Vector<Float> yIn = geFacCol.getColumn();
1379 Vector<Bool> maskIn(xIn.nelements(),True);
1380
1381// Interpolate (and extrapolate) with desired method
1382
1383 Int method = 0;
1384 convertInterpString(method, methodStr);
1385//
1386 Vector<Float> yOut;
1387 Vector<Bool> maskOut;
1388 InterpolateArray1D<Float,Float>::interpolate(yOut, maskOut, xOut,
1389 xIn, yIn, maskIn, method,
1390 True, True);
1391// Apply
1392
1393 correctFromVector (pTabOut, in, doAll, yOut);
1394}
1395
1396
1397void SDMath::correctFromVector (SDMemTable* pTabOut, const SDMemTable& in,
1398 Bool doAll, const Vector<Float>& factor) const
1399{
1400// For operations only on specified cursor location
1401
1402 IPosition start, end;
1403 getCursorLocation(start, end, in);
1404
1405// Loop over rows and interpolate correction factor
1406
1407 const uInt axis = asap::ChanAxis;
1408 for (uInt i=0; i < in.nRow(); ++i) {
1409
1410// Get data
1411
1412 MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
1413 Array<Float>& valuesIn = dataIn.getRWArray();
1414 const Array<Bool>& maskIn = dataIn.getMask();
1415
1416// Apply factor
1417
1418 if (doAll) {
1419 VectorIterator<Float> itValues(valuesIn, asap::ChanAxis);
1420 while (!itValues.pastEnd()) {
1421 itValues.vector() *= factor(i);
1422 itValues.next();
1423 }
1424 } else {
1425 Array<Float> valuesIn2 = valuesIn(start,end);
1426 valuesIn2 *= factor(i);
1427 valuesIn(start,end) = valuesIn2;
1428 }
1429
1430// Write out
1431
1432 SDContainer sc = in.getSDContainer(i);
1433 putDataInSDC(sc, valuesIn, maskIn);
1434//
1435 pTabOut->putSDContainer(sc);
1436 }
1437}
1438
1439
Note: See TracBrowser for help on using the repository browser.