source: trunk/src/SDMath.cc @ 163

Last change on this file since 163 was 163, checked in by kil064, 19 years ago

consolidate code in 'private' functions
rerwork function 'avergae_pol' to reshape the output
to 1 pol rather than replciating in all pols

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