source: trunk/src/SDMath.cc@ 194

Last change on this file since 194 was 185, checked in by mar637, 20 years ago

cosmetics

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.5 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/Utilities/Assert.h>
45#include <casa/Exceptions.h>
46
47#include <scimath/Mathematics/VectorKernel.h>
48#include <scimath/Mathematics/Convolver.h>
49
50#include <tables/Tables/Table.h>
51#include <tables/Tables/ScalarColumn.h>
52#include <tables/Tables/ArrayColumn.h>
53
54#include <lattices/Lattices/LatticeUtilities.h>
55#include <lattices/Lattices/RebinLattice.h>
56#include <coordinates/Coordinates/SpectralCoordinate.h>
57#include <coordinates/Coordinates/CoordinateSystem.h>
58#include <coordinates/Coordinates/CoordinateUtil.h>
59
60#include "MathUtils.h"
61#include "SDContainer.h"
62#include "SDMemTable.h"
63
64#include "SDMath.h"
65
66using namespace casa;
67using namespace asap;
68
69
70SDMath::SDMath()
71{;}
72
73SDMath::SDMath(const SDMath& other)
74{
75
76// No state
77
78}
79
80SDMath& SDMath::operator=(const SDMath& other)
81{
82 if (this != &other) {
83// No state
84 }
85 return *this;
86}
87
88SDMath::~SDMath()
89{;}
90
91
92CountedPtr<SDMemTable> SDMath::average(const Block<CountedPtr<SDMemTable> >& in,
93 const Vector<Bool>& mask, Bool scanAv,
94 const std::string& weightStr)
95//
96// Weighted averaging of spectra from one or more Tables.
97//
98{
99
100// Convert weight type
101
102 WeightType wtType = NONE;
103 convertWeightString(wtType, weightStr);
104
105// Create output Table by cloning from the first table
106
107 SDMemTable* pTabOut = new SDMemTable(*in[0],True);
108
109// Setup
110
111 const uInt axis = 3; // Spectral axis
112 IPosition shp = in[0]->rowAsMaskedArray(0).shape(); // Must not change
113 Array<Float> arr(shp);
114 Array<Bool> barr(shp);
115 const Bool useMask = (mask.nelements() == shp(axis));
116
117// Columns from Tables
118
119 ROArrayColumn<Float> tSysCol;
120 ROScalarColumn<Double> mjdCol;
121 ROScalarColumn<String> srcNameCol;
122 ROScalarColumn<Double> intCol;
123 ROArrayColumn<uInt> fqIDCol;
124
125// Create accumulation MaskedArray. We accumulate for each channel,if,pol,beam
126// Note that the mask of the accumulation array will ALWAYS remain ALL True.
127// The MA is only used so that when data which is masked Bad is added to it,
128// that data does not contribute.
129
130 Array<Float> zero(shp);
131 zero=0.0;
132 Array<Bool> good(shp);
133 good = True;
134 MaskedArray<Float> sum(zero,good);
135
136// Counter arrays
137
138 Array<Float> nPts(shp); // Number of points
139 nPts = 0.0;
140 Array<Float> nInc(shp); // Increment
141 nInc = 1.0;
142
143// Create accumulation Array for variance. We accumulate for
144// each if,pol,beam, but average over channel. So we need
145// a shape with one less axis dropping channels.
146
147 const uInt nAxesSub = shp.nelements() - 1;
148 IPosition shp2(nAxesSub);
149 for (uInt i=0,j=0; i<(nAxesSub+1); i++) {
150 if (i!=axis) {
151 shp2(j) = shp(i);
152 j++;
153 }
154 }
155 Array<Float> sumSq(shp2);
156 sumSq = 0.0;
157 IPosition pos2(nAxesSub,0); // For indexing
158
159// Time-related accumulators
160
161 Double time;
162 Double timeSum = 0.0;
163 Double intSum = 0.0;
164 Double interval = 0.0;
165
166// To get the right shape for the Tsys accumulator we need to
167// access a column from the first table. The shape of this
168// array must not change
169
170 Array<Float> tSysSum;
171 {
172 const Table& tabIn = in[0]->table();
173 tSysCol.attach(tabIn,"TSYS");
174 tSysSum.resize(tSysCol.shape(0));
175 }
176 tSysSum =0.0;
177 Array<Float> tSys;
178
179// Scan and row tracking
180
181 Int oldScanID = 0;
182 Int outScanID = 0;
183 Int scanID = 0;
184 Int rowStart = 0;
185 Int nAccum = 0;
186 Int tableStart = 0;
187
188// Source and FreqID
189
190 String sourceName, oldSourceName, sourceNameStart;
191 Vector<uInt> freqID, freqIDStart, oldFreqID;
192
193// Loop over tables
194
195 Float fac = 1.0;
196 const uInt nTables = in.nelements();
197 for (uInt iTab=0; iTab<nTables; iTab++) {
198
199// Attach columns to Table
200
201 const Table& tabIn = in[iTab]->table();
202 tSysCol.attach(tabIn, "TSYS");
203 mjdCol.attach(tabIn, "TIME");
204 srcNameCol.attach(tabIn, "SRCNAME");
205 intCol.attach(tabIn, "INTERVAL");
206 fqIDCol.attach(tabIn, "FREQID");
207
208// Loop over rows in Table
209
210 const uInt nRows = in[iTab]->nRow();
211 for (uInt iRow=0; iRow<nRows; iRow++) {
212
213// Check conformance
214
215 IPosition shp2 = in[iTab]->rowAsMaskedArray(iRow).shape();
216 if (!shp.isEqual(shp2)) {
217 throw (AipsError("Shapes for all rows must be the same"));
218 }
219
220// If we are not doing scan averages, make checks for source and
221// frequency setup and warn if averaging across them
222
223// Get copy of Scan Container for this row
224
225 SDContainer sc = in[iTab]->getSDContainer(iRow);
226 scanID = sc.scanid;
227
228// Get quantities from columns
229
230 srcNameCol.getScalar(iRow, sourceName);
231 mjdCol.get(iRow, time);
232 tSysCol.get(iRow, tSys);
233 intCol.get(iRow, interval);
234 fqIDCol.get(iRow, freqID);
235
236// Initialize first source and freqID
237
238 if (iRow==0 && iTab==0) {
239 sourceNameStart = sourceName;
240 freqIDStart = freqID;
241 }
242
243// If we are doing scan averages, see if we are at the end of an
244// accumulation period (scan). We must check soutce names too,
245// since we might have two tables with one scan each but different
246// source names; we shouldn't average different sources together
247
248 if (scanAv && ( (scanID != oldScanID) ||
249 (iRow==0 && iTab>0 && sourceName!=oldSourceName))) {
250
251// Normalize data in 'sum' accumulation array according to weighting scheme
252
253 normalize(sum, sumSq, nPts, wtType, axis, nAxesSub);
254
255// Fill scan container. The source and freqID come from the
256// first row of the first table that went into this average (
257// should be the same for all rows in the scan average)
258
259 Float nR(nAccum);
260 fillSDC(sc, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
261 timeSum/nR, intSum, sourceNameStart, freqIDStart);
262
263// Write container out to Table
264
265 pTabOut->putSDContainer(sc);
266
267// Reset accumulators
268
269 sum = 0.0;
270 sumSq = 0.0;
271 nAccum = 0;
272//
273 tSysSum =0.0;
274 timeSum = 0.0;
275 intSum = 0.0;
276 nPts = 0.0; // reset this too!!!
277
278// Increment
279
280 rowStart = iRow; // First row for next accumulation
281 tableStart = iTab; // First table for next accumulation
282 sourceNameStart = sourceName; // First source name for next accumulation
283 freqIDStart = freqID; // First FreqID for next accumulation
284//
285 oldScanID = scanID;
286 outScanID += 1; // Scan ID for next accumulation period
287
288}
289
290// Accumulate
291
292 accumulate(timeSum, intSum, nAccum, sum, sumSq, nPts, tSysSum,
293 tSys, nInc, mask, time, interval, in, iTab, iRow, axis,
294 nAxesSub, useMask, wtType);
295//
296 oldSourceName = sourceName;
297 oldFreqID = freqID;
298 }
299 }
300
301// OK at this point we have accumulation data which is either
302// - accumulated from all tables into one row
303// or
304// - accumulated from the last scan average
305//
306// Normalize data in 'sum' accumulation array according to weighting scheme
307 normalize(sum, sumSq, nPts, wtType, axis, nAxesSub);
308
309// Create and fill container. The container we clone will be from
310// the last Table and the first row that went into the current
311// accumulation. It probably doesn't matter that much really...
312
313 Float nR(nAccum);
314 SDContainer sc = in[tableStart]->getSDContainer(rowStart);
315 fillSDC(sc, sum.getMask(), sum.getArray(), tSysSum/nR, outScanID,
316 timeSum/nR, intSum, sourceNameStart, freqIDStart);
317//
318 pTabOut->putSDContainer(sc);
319/*
320 cout << endl;
321 cout << "Last accumulation for output scan ID " << outScanID << endl;
322 cout << " The first row in this accumulation is " << rowStart << endl;
323 cout << " The number of rows accumulated is " << nAccum << endl;
324 cout << " The first table in this accumulation is " << tableStart << endl;
325 cout << " The first source in this accumulation is " << sourceNameStart << endl;
326 cout << " The first freqID in this accumulation is " << freqIDStart << endl;
327 cout << " Average time stamp = " << timeSum/nR << endl;
328 cout << " Integrated time = " << intSum << endl;
329*/
330 return CountedPtr<SDMemTable>(pTabOut);
331}
332
333
334
335CountedPtr<SDMemTable>
336SDMath::quotient(const CountedPtr<SDMemTable>& on,
337 const CountedPtr<SDMemTable>& off)
338{
339//
340// Compute quotient spectrum
341//
342 const uInt nRows = on->nRow();
343 if (off->nRow() != nRows) {
344 throw (AipsError("Input Scan Tables must have the same number of rows"));
345 }
346
347// Input Tables and columns
348
349 Table ton = on->table();
350 Table toff = off->table();
351 ROArrayColumn<Float> tsys(toff, "TSYS");
352 ROScalarColumn<Double> mjd(ton, "TIME");
353 ROScalarColumn<Double> integr(ton, "INTERVAL");
354 ROScalarColumn<String> srcn(ton, "SRCNAME");
355 ROArrayColumn<uInt> freqidc(ton, "FREQID");
356
357// Output Table cloned from input
358
359 SDMemTable* pTabOut = new SDMemTable(*on, True);
360
361// Loop over rows
362
363 for (uInt i=0; i<nRows; i++) {
364 MaskedArray<Float> mon(on->rowAsMaskedArray(i));
365 MaskedArray<Float> moff(off->rowAsMaskedArray(i));
366 IPosition ipon = mon.shape();
367 IPosition ipoff = moff.shape();
368//
369 Array<Float> tsarr;
370 tsys.get(i, tsarr);
371 if (ipon != ipoff && ipon != tsarr.shape()) {
372 throw(AipsError("on/off not conformant"));
373 }
374
375// Compute quotient
376
377 MaskedArray<Float> tmp = (mon-moff);
378 Array<Float> out(tmp.getArray());
379 out /= moff;
380 out *= tsarr;
381 Array<Bool> outflagsb = mon.getMask() && moff.getMask();
382
383// Fill container for this row
384
385 SDContainer sc = on->getSDContainer(i);
386//
387 putDataInSDC(sc, out, outflagsb);
388 sc.putTsys(tsarr);
389 sc.scanid = i;
390
391// Put new row in output Table
392
393 pTabOut->putSDContainer(sc);
394 }
395//
396 return CountedPtr<SDMemTable>(pTabOut);
397}
398
399
400
401std::vector<float> SDMath::statistic(const CountedPtr<SDMemTable>& in,
402 const std::vector<bool>& mask,
403 const String& which)
404//
405// Perhaps iteration over pol/beam/if should be in here
406// and inside the nrow iteration ?
407//
408{
409 const uInt nRow = in->nRow();
410 std::vector<float> result(nRow);
411 Vector<Bool> msk(mask);
412
413// Specify cursor location
414
415 IPosition start, end;
416 getCursorLocation(start, end, *in);
417
418// Loop over rows
419
420 const uInt nEl = msk.nelements();
421 for (uInt ii=0; ii < in->nRow(); ++ii) {
422
423// Get row and deconstruct
424
425 MaskedArray<Float> marr(in->rowAsMaskedArray(ii));
426 Array<Float> arr = marr.getArray();
427 Array<Bool> barr = marr.getMask();
428
429// Access desired piece of data
430
431 Array<Float> v((arr(start,end)).nonDegenerate());
432 Array<Bool> m((barr(start,end)).nonDegenerate());
433
434// Apply OTF mask
435
436 MaskedArray<Float> tmp;
437 if (m.nelements()==nEl) {
438 tmp.setData(v,m&&msk);
439 } else {
440 tmp.setData(v,m);
441 }
442
443// Get statistic
444
445 result[ii] = mathutil::statistics(which, tmp);
446 }
447//
448 return result;
449}
450
451
452SDMemTable* SDMath::bin(const SDMemTable& in, Int width)
453{
454 SDHeader sh = in.getSDHeader();
455 SDMemTable* pTabOut = new SDMemTable(in, True);
456
457// Bin up SpectralCoordinates
458
459 IPosition factors(1);
460 factors(0) = width;
461 for (uInt j=0; j<in.nCoordinates(); ++j) {
462 CoordinateSystem cSys;
463 cSys.addCoordinate(in.getCoordinate(j));
464 CoordinateSystem cSysBin =
465 CoordinateUtil::makeBinnedCoordinateSystem(factors, cSys, False);
466//
467 SpectralCoordinate sCBin = cSysBin.spectralCoordinate(0);
468 pTabOut->setCoordinate(sCBin, j);
469 }
470
471// Use RebinLattice to find shape
472
473 IPosition shapeIn(1,sh.nchan);
474 IPosition shapeOut = RebinLattice<Float>::rebinShape(shapeIn, factors);
475 sh.nchan = shapeOut(0);
476 pTabOut->putSDHeader(sh);
477
478
479// Loop over rows and bin along channel axis
480
481 const uInt axis = 3;
482 for (uInt i=0; i < in.nRow(); ++i) {
483 SDContainer sc = in.getSDContainer(i);
484//
485 Array<Float> tSys(sc.getTsys()); // Get it out before sc changes shape
486
487// Bin up spectrum
488
489 MaskedArray<Float> marr(in.rowAsMaskedArray(i));
490 MaskedArray<Float> marrout;
491 LatticeUtilities::bin(marrout, marr, axis, width);
492
493// Put back the binned data and flags
494
495 IPosition ip2 = marrout.shape();
496 sc.resize(ip2);
497//
498 putDataInSDC(sc, marrout.getArray(), marrout.getMask());
499
500// Bin up Tsys.
501
502 Array<Bool> allGood(tSys.shape(),True);
503 MaskedArray<Float> tSysIn(tSys, allGood, True);
504//
505 MaskedArray<Float> tSysOut;
506 LatticeUtilities::bin(tSysOut, tSysIn, axis, width);
507 sc.putTsys(tSysOut.getArray());
508//
509 pTabOut->putSDContainer(sc);
510 }
511 return pTabOut;
512}
513
514SDMemTable* SDMath::simpleOperate(const SDMemTable& in, Float val, Bool doAll,
515 uInt what)
516//
517// what = 0 Multiply
518// 1 Add
519{
520 SDMemTable* pOut = new SDMemTable(in,False);
521 const Table& tOut = pOut->table();
522 ArrayColumn<Float> spec(tOut,"SPECTRA");
523//
524 if (doAll) {
525 for (uInt i=0; i < tOut.nrow(); i++) {
526
527// Get
528
529 MaskedArray<Float> marr(pOut->rowAsMaskedArray(i));
530
531// Operate
532
533 if (what==0) {
534 marr *= val;
535 } else if (what==1) {
536 marr += val;
537 }
538
539// Put
540
541 spec.put(i, marr.getArray());
542 }
543 } else {
544
545// Get cursor location
546
547 IPosition start, end;
548 getCursorLocation(start, end, in);
549//
550 for (uInt i=0; i < tOut.nrow(); i++) {
551
552// Get
553
554 MaskedArray<Float> dataIn(pOut->rowAsMaskedArray(i));
555
556// Modify. More work than we would like to deal with the mask
557
558 Array<Float>& values = dataIn.getRWArray();
559 Array<Bool> mask(dataIn.getMask());
560//
561 Array<Float> values2 = values(start,end);
562 Array<Bool> mask2 = mask(start,end);
563 MaskedArray<Float> t(values2,mask2);
564 if (what==0) {
565 t *= val;
566 } else if (what==1) {
567 t += val;
568 }
569 values(start, end) = t.getArray(); // Write back into 'dataIn'
570
571// Put
572 spec.put(i, dataIn.getArray());
573 }
574 }
575//
576 return pOut;
577}
578
579
580
581SDMemTable* SDMath::averagePol(const SDMemTable& in, const Vector<Bool>& mask)
582//
583// Average all polarizations together, weighted by variance
584//
585{
586// WeightType wtType = NONE;
587// convertWeightString(wtType, weight);
588
589 const uInt nRows = in.nRow();
590 const uInt polAxis = 2; // Polarization axis
591 const uInt chanAxis = 3; // Spectrum axis
592
593// Create output Table and reshape number of polarizations
594
595 Bool clear=True;
596 SDMemTable* pTabOut = new SDMemTable(in, clear);
597 SDHeader header = pTabOut->getSDHeader();
598 header.npol = 1;
599 pTabOut->putSDHeader(header);
600
601// Shape of input and output data
602
603 const IPosition& shapeIn = in.rowAsMaskedArray(0u, False).shape();
604 IPosition shapeOut(shapeIn);
605 shapeOut(polAxis) = 1; // Average all polarizations
606//
607 const uInt nChan = shapeIn(chanAxis);
608 const IPosition vecShapeOut(4,1,1,1,nChan); // A multi-dim form of a Vector shape
609 IPosition start(4), end(4);
610
611// Output arrays
612
613 Array<Float> outData(shapeOut, 0.0);
614 Array<Bool> outMask(shapeOut, True);
615 const IPosition axes(2, 2, 3); // pol-channel plane
616//
617 const Bool useMask = (mask.nelements() == shapeIn(chanAxis));
618
619// Loop over rows
620
621 for (uInt iRow=0; iRow<nRows; iRow++) {
622
623// Get data for this row
624
625 MaskedArray<Float> marr(in.rowAsMaskedArray(iRow));
626 Array<Float>& arr = marr.getRWArray();
627 const Array<Bool>& barr = marr.getMask();
628
629// Make iterators to iterate by pol-channel planes
630
631 ReadOnlyArrayIterator<Float> itDataPlane(arr, axes);
632 ReadOnlyArrayIterator<Bool> itMaskPlane(barr, axes);
633
634// Accumulations
635
636 Float fac = 1.0;
637 Vector<Float> vecSum(nChan,0.0);
638
639// Iterate through data by pol-channel planes
640
641 while (!itDataPlane.pastEnd()) {
642
643// Iterate through plane by polarization and accumulate Vectors
644
645 Vector<Float> t1(nChan); t1 = 0.0;
646 Vector<Bool> t2(nChan); t2 = True;
647 MaskedArray<Float> vecSum(t1,t2);
648 Float varSum = 0.0;
649 {
650 ReadOnlyVectorIterator<Float> itDataVec(itDataPlane.array(), 1);
651 ReadOnlyVectorIterator<Bool> itMaskVec(itMaskPlane.array(), 1);
652 while (!itDataVec.pastEnd()) {
653
654// Create MA of data & mask (optionally including OTF mask) and get variance
655
656 if (useMask) {
657 const MaskedArray<Float> spec(itDataVec.vector(),mask&&itMaskVec.vector());
658 fac = 1.0 / variance(spec);
659 } else {
660 const MaskedArray<Float> spec(itDataVec.vector(),itMaskVec.vector());
661 fac = 1.0 / variance(spec);
662 }
663
664// Normalize spectrum (without OTF mask) and accumulate
665
666 const MaskedArray<Float> spec(fac*itDataVec.vector(), itMaskVec.vector());
667 vecSum += spec;
668 varSum += fac;
669
670// Next
671
672 itDataVec.next();
673 itMaskVec.next();
674 }
675 }
676
677// Normalize summed spectrum
678
679 vecSum /= varSum;
680
681// FInd position in input data array. We are iterating by pol-channel
682// plane so all that will change is beam and IF and that's what we want.
683
684 IPosition pos = itDataPlane.pos();
685
686// Write out data. This is a bit messy. We have to reform the Vector
687// accumulator into an Array of shape (1,1,1,nChan)
688
689 start = pos;
690 end = pos;
691 end(chanAxis) = nChan-1;
692 outData(start,end) = vecSum.getArray().reform(vecShapeOut);
693 outMask(start,end) = vecSum.getMask().reform(vecShapeOut);
694
695// Step to next beam/IF combination
696
697 itDataPlane.next();
698 itMaskPlane.next();
699 }
700
701// Generate output container and write it to output table
702
703 SDContainer sc = in.getSDContainer();
704 sc.resize(shapeOut);
705//
706 putDataInSDC(sc, outData, outMask);
707 pTabOut->putSDContainer(sc);
708 }
709//
710 return pTabOut;
711}
712
713
714SDMemTable* SDMath::smooth(const SDMemTable& in,
715 const casa::String& kernelType,
716 casa::Float width, Bool doAll)
717{
718
719// Number of channels
720
721 const uInt chanAxis = 3; // Spectral axis
722 SDHeader sh = in.getSDHeader();
723 const uInt nChan = sh.nchan;
724
725// Generate Kernel
726
727 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
728 Vector<Float> kernel = VectorKernel::make(type, width, nChan, True, False);
729
730// Generate Convolver
731
732 IPosition shape(1,nChan);
733 Convolver<Float> conv(kernel, shape);
734
735// New Table
736
737 SDMemTable* pTabOut = new SDMemTable(in,True);
738
739// Get cursor location
740
741 IPosition start, end;
742 getCursorLocation(start, end, in);
743//
744 IPosition shapeOut(4,1);
745
746// Output Vectors
747
748 Vector<Float> valuesOut(nChan);
749 Vector<Bool> maskOut(nChan);
750
751// Loop over rows in Table
752
753 for (uInt ri=0; ri < in.nRow(); ++ri) {
754
755// Get copy of data
756
757 const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(ri));
758 AlwaysAssert(dataIn.shape()(chanAxis)==nChan, AipsError);
759//
760 Array<Float> valuesIn = dataIn.getArray();
761 Array<Bool> maskIn = dataIn.getMask();
762
763// Branch depending on whether we smooth all locations or just
764// those pointed at by the current selection cursor
765
766 if (doAll) {
767 uInt axis = 3;
768 VectorIterator<Float> itValues(valuesIn, axis);
769 VectorIterator<Bool> itMask(maskIn, axis);
770 while (!itValues.pastEnd()) {
771
772// Smooth
773 if (kernelType==VectorKernel::HANNING) {
774 mathutil::hanning(valuesOut, maskOut, itValues.vector(), itMask.vector());
775 itMask.vector() = maskOut;
776 } else {
777 mathutil::replaceMaskByZero(itValues.vector(), itMask.vector());
778 conv.linearConv(valuesOut, itValues.vector());
779 }
780//
781 itValues.vector() = valuesOut;
782//
783 itValues.next();
784 itMask.next();
785 }
786 } else {
787
788// Set multi-dim Vector shape
789
790 shapeOut(chanAxis) = valuesIn.shape()(chanAxis);
791
792// Stuff about with shapes so that we don't have conformance run-time errors
793
794 Vector<Float> valuesIn2 = valuesIn(start,end).nonDegenerate();
795 Vector<Bool> maskIn2 = maskIn(start,end).nonDegenerate();
796
797// Smooth
798
799 if (kernelType==VectorKernel::HANNING) {
800 mathutil::hanning(valuesOut, maskOut, valuesIn2, maskIn2);
801 maskIn(start,end) = maskOut.reform(shapeOut);
802 } else {
803 mathutil::replaceMaskByZero(valuesIn2, maskIn2);
804 conv.linearConv(valuesOut, valuesIn2);
805 }
806//
807 valuesIn(start,end) = valuesOut.reform(shapeOut);
808 }
809
810// Create and put back
811
812 SDContainer sc = in.getSDContainer(ri);
813 putDataInSDC(sc, valuesIn, maskIn);
814//
815 pTabOut->putSDContainer(sc);
816 }
817//
818 return pTabOut;
819}
820
821
822
823
824
825// 'private' functions
826
827void SDMath::fillSDC(SDContainer& sc,
828 const Array<Bool>& mask,
829 const Array<Float>& data,
830 const Array<Float>& tSys,
831 Int scanID, Double timeStamp,
832 Double interval, const String& sourceName,
833 const Vector<uInt>& freqID)
834{
835// Data and mask
836
837 putDataInSDC(sc, data, mask);
838
839// TSys
840
841 sc.putTsys(tSys);
842
843// Time things
844
845 sc.timestamp = timeStamp;
846 sc.interval = interval;
847 sc.scanid = scanID;
848//
849 sc.sourcename = sourceName;
850 sc.putFreqMap(freqID);
851}
852
853void SDMath::normalize(MaskedArray<Float>& sum,
854 const Array<Float>& sumSq,
855 const Array<Float>& nPts,
856 WeightType wtType, Int axis,
857 Int nAxesSub)
858{
859 IPosition pos2(nAxesSub,0);
860//
861 if (wtType==NONE) {
862
863// We just average by the number of points accumulated.
864// We need to make a MA out of nPts so that no divide by
865// zeros occur
866
867 MaskedArray<Float> t(nPts, (nPts>Float(0.0)));
868 sum /= t;
869 } else if (wtType==VAR) {
870
871// Normalize each spectrum by sum(1/var) where the variance
872// is worked out for each spectrum
873
874 Array<Float>& data = sum.getRWArray();
875 VectorIterator<Float> itData(data, axis);
876 while (!itData.pastEnd()) {
877 pos2 = itData.pos().getFirst(nAxesSub);
878 itData.vector() /= sumSq(pos2);
879 itData.next();
880 }
881 } else if (wtType==TSYS) {
882 }
883}
884
885
886void SDMath::accumulate(Double& timeSum, Double& intSum, Int& nAccum,
887 MaskedArray<Float>& sum, Array<Float>& sumSq,
888 Array<Float>& nPts, Array<Float>& tSysSum,
889 const Array<Float>& tSys, const Array<Float>& nInc,
890 const Vector<Bool>& mask, Double time, Double interval,
891 const Block<CountedPtr<SDMemTable> >& in,
892 uInt iTab, uInt iRow, uInt axis,
893 uInt nAxesSub, Bool useMask,
894 WeightType wtType)
895{
896
897// Get data
898
899 MaskedArray<Float> dataIn(in[iTab]->rowAsMaskedArray(iRow));
900 Array<Float>& valuesIn = dataIn.getRWArray(); // writable reference
901 const Array<Bool>& maskIn = dataIn.getMask(); // RO reference
902//
903 if (wtType==NONE) {
904 const MaskedArray<Float> n(nInc,dataIn.getMask());
905 nPts += n; // Only accumulates where mask==T
906 } else if (wtType==VAR) {
907
908// We are going to average the data, weighted by the noise for each pol, beam and IF.
909// So therefore we need to iterate through by spectrum (axis 3)
910
911 VectorIterator<Float> itData(valuesIn, axis);
912 ReadOnlyVectorIterator<Bool> itMask(maskIn, axis);
913 Float fac = 1.0;
914 IPosition pos(nAxesSub,0);
915//
916 while (!itData.pastEnd()) {
917
918// Make MaskedArray of Vector, optionally apply OTF mask, and find scaling factor
919
920 if (useMask) {
921 MaskedArray<Float> tmp(itData.vector(),mask&&itMask.vector());
922 fac = 1.0/variance(tmp);
923 } else {
924 MaskedArray<Float> tmp(itData.vector(),itMask.vector());
925 fac = 1.0/variance(tmp);
926 }
927
928// Scale data
929
930 itData.vector() *= fac; // Writes back into 'dataIn'
931//
932// Accumulate variance per if/pol/beam averaged over spectrum
933// This method to get pos2 from itData.pos() is only valid
934// because the spectral axis is the last one (so we can just
935// copy the first nAXesSub positions out)
936
937 pos = itData.pos().getFirst(nAxesSub);
938 sumSq(pos) += fac;
939//
940 itData.next();
941 itMask.next();
942 }
943 } else if (wtType==TSYS) {
944 }
945
946// Accumulate sum of (possibly scaled) data
947
948 sum += dataIn;
949
950// Accumulate Tsys, time, and interval
951
952 tSysSum += tSys;
953 timeSum += time;
954 intSum += interval;
955 nAccum += 1;
956}
957
958
959
960
961void SDMath::getCursorLocation(IPosition& start, IPosition& end,
962 const SDMemTable& in)
963{
964 const uInt nDim = 4;
965 const uInt i = in.getBeam();
966 const uInt j = in.getIF();
967 const uInt k = in.getPol();
968 const uInt n = in.nChan();
969//
970 start.resize(nDim);
971 start(0) = i;
972 start(1) = j;
973 start(2) = k;
974 start(3) = 0;
975//
976 end.resize(nDim);
977 end(0) = i;
978 end(1) = j;
979 end(2) = k;
980 end(3) = n-1;
981}
982
983
984void SDMath::convertWeightString(WeightType& wtType, const std::string& weightStr)
985{
986 String tStr(weightStr);
987 tStr.upcase();
988 if (tStr.contains(String("NONE"))) {
989 wtType = NONE;
990 } else if (tStr.contains(String("VAR"))) {
991 wtType = VAR;
992 } else if (tStr.contains(String("TSYS"))) {
993 wtType = TSYS;
994 throw(AipsError("T_sys weighting not yet implemented"));
995 } else {
996 throw(AipsError("Unrecognized weighting type"));
997 }
998}
999
1000void SDMath::putDataInSDC(SDContainer& sc, const Array<Float>& data,
1001 const Array<Bool>& mask)
1002{
1003 sc.putSpectrum(data);
1004//
1005 Array<uChar> outflags(data.shape());
1006 convertArray(outflags,!mask);
1007 sc.putFlags(outflags);
1008}
Note: See TracBrowser for help on using the repository browser.