source: trunk/src/SDMath.cc@ 157

Last change on this file since 157 was 155, checked in by kil064, 21 years ago

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