source: trunk/src/RowAccumulator.cpp@ 2127

Last change on this file since 2127 was 2127, checked in by WataruKawasaki, 14 years ago

New Development: No

JIRA Issue: Yes CAS-2776

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): SD

Description: a minor bugfix


File size: 5.5 KB
Line 
1//
2// C++ Implementation: RowAccumulator
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <Malte.Marquarding@csiro.au>, (C) 2005
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <casa/iomanip.h>
14#include <casa/Arrays/MaskArrMath.h>
15#include <casa/Arrays/MaskArrLogi.h>
16#include <casa/Arrays/ArrayMath.h>
17#include <casa/Arrays/ArrayLogical.h>
18#include "RowAccumulator.h"
19
20using namespace casa;
21using namespace asap;
22
23RowAccumulator::RowAccumulator(WeightType wt) : weightType_(wt), initialized_(False)
24{
25 reset();
26}
27
28RowAccumulator::~RowAccumulator()
29{
30}
31
32
33void RowAccumulator::reset(const uInt size, const uInt tsysSize)
34{
35 Vector<Bool> m(size, True);
36
37 spectrum_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
38 spectrumNoMask_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
39
40 n_.setData(Vector<uInt>(size, 0), Vector<Bool>(size, True));
41 nNoMask_.setData(Vector<uInt>(size, 0), Vector<Bool>(size, True));
42
43 weightSum_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
44 weightSumNoMask_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
45
46 tsysSum_.resize(tsysSize); tsysSum_=0.0;
47 tsysSumNoMask_.resize(tsysSize); tsysSumNoMask_=0.0;
48
49 intervalSum_ = 0.0;
50 intervalSumNoMask_ = 0.0;
51
52 timeSum_ = 0.0;
53 timeSumNoMask_ = 0.0;
54}
55
56void RowAccumulator::initialize(const uInt size, const uInt tsysSize)
57{
58 reset(size, tsysSize);
59 initialized_ = True;
60}
61
62void RowAccumulator::add(const Vector<Float>& v,
63 const Vector<Bool>& m,
64 const Vector<Float>& tsys,
65 const Double interval,
66 const Double time)
67{
68 uInt size = v.nelements();
69 //if (size != m.nelements()) raiseError;
70 if (!initialized_) initialize(size, tsys.nelements());
71
72 addSpectrum(v, m, tsys, interval, time);
73}
74
75void RowAccumulator::addSpectrum(const Vector<Float>& v,
76 const Vector<Bool>& m,
77 const Vector<Float>& tsys,
78 const Double interval,
79 const Double time)
80{
81 doAddSpectrum(v, m, tsys, interval, time, False);
82 doAddSpectrum(v, m, tsys, interval, time, True); // CAS-2776
83}
84
85void RowAccumulator::doAddSpectrum(const Vector<Float>& v,
86 const Vector<Bool>& m,
87 const Vector<Float>& tsys,
88 const Double interval,
89 const Double time,
90 const Bool ignoreMask)
91{
92 Vector<Float> vUse = v.copy();
93 Vector<Bool> mUse = m.copy();
94 if (ignoreMask) mUse = !mUse;
95
96 MaskedArray<Float> vadd(vUse, mUse);
97 Float totalWeight = getTotalWeight(vadd, tsys, interval, time, ignoreMask);
98 vadd *= totalWeight;
99 MaskedArray<Float> wadd(Vector<Float>(mUse.nelements(), totalWeight), mUse);
100 MaskedArray<uInt> inc(Vector<uInt>(mUse.nelements(), 1), mUse);
101
102 if (ignoreMask) {
103 spectrumNoMask_ += vadd;
104 weightSumNoMask_ += wadd;
105 nNoMask_ += inc;
106 } else {
107 spectrum_ += vadd;
108 weightSum_ += wadd;
109 n_ += inc;
110 }
111}
112
113Float RowAccumulator::getTotalWeight(const MaskedArray<Float>& data,
114 const Vector<Float>& tsys,
115 const Double interval,
116 const Double time,
117 const Bool ignoreMask)
118{
119 Float totalWeight = 1.0;
120
121 Vector<Bool> m = data.getMask();
122 if (!allEQ(m, False)) { // only add these if not everything masked
123 totalWeight *= addTsys(tsys, ignoreMask);
124 totalWeight *= addInterval(interval, ignoreMask);
125 addTime(time, ignoreMask);
126 }
127
128 if (weightType_ == W_VAR) {
129 Float fac = 1.0/variance(data);
130 if (!ignoreMask && (m.nelements() == userMask_.nelements()))
131 fac = 1.0/variance(data(userMask_));
132
133 totalWeight *= fac;
134 }
135
136 return totalWeight;
137}
138
139Float RowAccumulator::addTsys(const Vector<Float>& v, Bool ignoreMask)
140{
141 // @fixme this assume tsys is the same for all channels
142
143 Float w = 1.0;
144 if (ignoreMask) {
145 tsysSumNoMask_ += v[0];
146 } else {
147 tsysSum_ += v[0];
148 }
149 if ( weightType_ == W_TSYS || weightType_ == W_TINTSYS ) {
150 w /= (v[0]*v[0]);
151 }
152 return w;
153}
154
155void RowAccumulator::addTime(Double t, Bool ignoreMask)
156{
157 if (ignoreMask) {
158 timeSumNoMask_ += t;
159 } else {
160 timeSum_ += t;
161 }
162}
163
164Float RowAccumulator::addInterval(Double inter, Bool ignoreMask)
165{
166 Float w = 1.0;
167 if (ignoreMask) {
168 intervalSumNoMask_ += inter;
169 } else {
170 intervalSum_ += inter;
171 }
172 if (weightType_ == W_TINT || weightType_ == W_TINTSYS) {
173 w /= Float(inter);
174 }
175 return w;
176}
177
178Vector<Float> RowAccumulator::getSpectrum() const
179{
180 return (spectrum_/weightSum_).getArray();
181}
182
183Double RowAccumulator::getTime() const
184{
185 return timeSum_/Float(max(n_));
186}
187
188Double RowAccumulator::getInterval() const
189{
190 return intervalSum_;
191}
192
193Vector<Bool> RowAccumulator::getMask() const
194{
195 // Return the "total" mask - False where no points have been accumulated.
196 return (n_.getArray() > uInt(0));
197}
198
199Vector<Float> RowAccumulator::getTsys() const
200{
201 // @fixme this assumes tsys.nelements() == 1
202 return tsysSum_/Float(max(n_));
203}
204
205void RowAccumulator::setUserMask(const Vector<Bool>& m)
206{
207 userMask_.resize();
208 userMask_ = m;
209}
210
211// Added by TT check the state of RowAccumulator
212Bool RowAccumulator::state() const
213{
214 return initialized_;
215}
216
217void RowAccumulator::replaceNaN()
218{
219 Vector<Float> v = spectrum_.getArray();
220 Vector<Float> w = weightSum_.getArray();
221 Vector<Float> vRef = spectrumNoMask_.getArray();
222 Vector<Float> wRef = weightSumNoMask_.getArray();
223
224 for (uInt i = 0; i < v.size(); ++i) {
225 if (w[i] == 0.0) {
226 v[i] = vRef[i];
227 w[i] = wRef[i];
228 }
229 }
230
231 spectrum_.setData(v, Vector<Bool>(v.nelements(), True));
232 weightSum_.setData(w, Vector<Bool>(w.nelements(), True));
233}
Note: See TracBrowser for help on using the repository browser.