source: trunk/src/RowAccumulator.cpp@ 2997

Last change on this file since 2997 was 2986, checked in by Takeshi Nakazato, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6582

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

STMath::average is changed so that empty averaged result (data to be
averaged are all flagged) is kept with all channel flags of 128 and
row flag of 1.


File size: 8.4 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
23namespace {
24inline uInt nNominal(MaskedArray<uInt> nvalid, MaskedArray<uInt> ninvalid)
25{
26 return max((allEQ(nvalid, uInt(0))) ? ninvalid : nvalid);
27}
28} // anonymous namespace
29
30RowAccumulator::RowAccumulator(WeightType wt) : weightType_(wt), initialized_(False)
31{
32 reset();
33}
34
35RowAccumulator::~RowAccumulator()
36{
37}
38
39
40void RowAccumulator::reset(const uInt size, const uInt tsysSize)
41{
42 Vector<Bool> m(size, True);
43
44 spectrum_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
45 spectrumNoMask_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
46
47 n_.setData(Vector<uInt>(size, 0), Vector<Bool>(size, True));
48 nNoMask_.setData(Vector<uInt>(size, 0), Vector<Bool>(size, True));
49
50 weightSum_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
51 weightSumNoMask_.setData(Vector<Float>(size, 0.0), Vector<Bool>(size, True));
52
53 tsysSum_.resize(tsysSize); tsysSum_=0.0;
54 tsysSumNoMask_.resize(tsysSize); tsysSumNoMask_=0.0;
55
56 intervalSum_ = 0.0;
57 intervalSumNoMask_ = 0.0;
58
59 timeSum_ = 0.0;
60 timeSumNoMask_ = 0.0;
61
62 initialized_ = False;
63}
64
65void RowAccumulator::initialize(const uInt size, const uInt tsysSize)
66{
67 reset(size, tsysSize);
68 initialized_ = True;
69}
70
71void RowAccumulator::add(const Vector<Float>& v,
72 const Vector<Bool>& m,
73 const Vector<Float>& tsys,
74 const Double interval,
75 const Double time)
76{
77 uInt size = v.nelements();
78 //if (size != m.nelements()) raiseError;
79 if (!initialized_) initialize(size, tsys.nelements());
80
81 addSpectrum(v, m, tsys, interval, time);
82}
83
84void RowAccumulator::addSpectrum(const Vector<Float>& v,
85 const Vector<Bool>& m,
86 const Vector<Float>& tsys,
87 const Double interval,
88 const Double time)
89{
90// doAddSpectrum(v, m, tsys, interval, time, False);
91// doAddSpectrum(v, m, tsys, interval, time, True); // CAS-2776
92 doAddSpectrum2( v, m, tsys, interval, time ) ;
93}
94
95void RowAccumulator::doAddSpectrum(const Vector<Float>& v,
96 const Vector<Bool>& m,
97 const Vector<Float>& tsys,
98 const Double interval,
99 const Double time,
100 const Bool inverseMask)
101{
102 Vector<Float> vUse = v.copy();
103 Vector<Bool> mUse = m.copy();
104 if (inverseMask) mUse = !mUse;
105 MaskedArray<Float> vadd(vUse, mUse);
106 Float totalWeight = getTotalWeight(vadd, tsys, interval, time, inverseMask);
107 vadd *= totalWeight;
108 MaskedArray<Float> wadd(Vector<Float>(mUse.nelements(), totalWeight), mUse);
109 MaskedArray<uInt> inc(Vector<uInt>(mUse.nelements(), 1), mUse);
110
111 if (inverseMask) {
112 spectrumNoMask_ += vadd;
113 weightSumNoMask_ += wadd;
114 nNoMask_ += inc;
115 } else {
116 spectrum_ += vadd;
117 weightSum_ += wadd;
118 n_ += inc;
119 }
120}
121
122void RowAccumulator::doAddSpectrum2(const Vector<Float>& v,
123 const Vector<Bool>& m,
124 const Vector<Float>& tsys,
125 const Double interval,
126 const Double time)
127{
128 const MaskedArray<Float> vadd(v, m);
129 const MaskedArray<Float> vaddN(v, !m);
130 Float totalWeight = getTotalWeight( vadd, tsys, interval, time, False ) ;
131 Float totalWeightNoMask = getTotalWeight( vaddN, tsys, interval, time, True ) ;
132 // process spectrum with input mask and spectrum with inverted mask
133 // together
134 // prepare pointers
135 Bool vD, mD ;
136 const Float *v_p = v.getStorage( vD ) ;
137 const Float *v_wp = v_p ;
138 const Bool *m_p = m.getStorage( mD ) ;
139 const Bool *m_wp = m_p ;
140
141 Bool spD, spND, wgtD, wgtND, nD, nND ;
142 Float *sp_p = spectrum_.getRWArrayStorage( spD ) ;
143 Float *sp_wp = sp_p ;
144 Float *wgt_p = weightSum_.getRWArrayStorage( wgtD ) ;
145 Float *wgt_wp = wgt_p ;
146 Float *spN_p = spectrumNoMask_.getRWArrayStorage( spND ) ;
147 Float *spN_wp = spN_p ;
148 Float *wgtN_p = weightSumNoMask_.getRWArrayStorage( wgtND ) ;
149 Float *wgtN_wp = wgtN_p ;
150 uInt *n_p = n_.getRWArrayStorage( nD ) ;
151 uInt *n_wp = n_p ;
152 uInt *nN_p = nNoMask_.getRWArrayStorage( nND ) ;
153 uInt *nN_wp = nN_p ;
154
155 // actual processing
156 uInt len = m.nelements() ;
157 // loop over channels
158 for ( uInt i = 0 ; i < len ; i++ ) {
159 // masks for spectrum_ and spectrumNoMask_ are not needed since
160 // it is initialized as True for all channels and those are constant
161 if ( *m_wp ) {
162 // mask is True
163 // add v * totalWeight to spectrum_
164 // add totalWeight to weightSum_
165 // increment n_
166 *sp_wp += *v_wp * totalWeight ;
167 *wgt_wp += totalWeight ;
168 *n_wp += 1 ;
169 }
170 else {
171 // mask is False
172 // add v * totalWeightNoMask to spectrumNoMask_
173 // add totalWeightNoMask to weightSumNoMask_
174 // increment nNoMask_
175 *spN_wp += *v_wp * totalWeightNoMask ;
176 *wgtN_wp += totalWeightNoMask ;
177 *nN_wp += 1 ;
178 }
179 sp_wp++ ;
180 wgt_wp++ ;
181 n_wp++ ;
182 spN_wp++ ;
183 wgtN_wp++ ;
184 nN_wp++ ;
185 v_wp++ ;
186 m_wp++ ;
187 }
188
189 // free storage
190 spectrum_.putArrayStorage( sp_p, spD ) ;
191 weightSum_.putArrayStorage( wgt_p, wgtD ) ;
192 spectrumNoMask_.putArrayStorage( spN_p, spND ) ;
193 weightSumNoMask_.putArrayStorage( wgtN_p, wgtND ) ;
194 n_.putArrayStorage( n_p, nD ) ;
195 nNoMask_.putArrayStorage( nN_p, nND ) ;
196
197 v.freeStorage( v_p, vD ) ;
198 m.freeStorage( m_p, mD ) ;
199}
200
201Float RowAccumulator::getTotalWeight(const MaskedArray<Float>& data,
202 const Vector<Float>& tsys,
203 const Double interval,
204 const Double time,
205 const Bool inverseMask)
206{
207 Float totalWeight = 1.0;
208 Vector<Bool> m = data.getMask();
209 Float tsysWeight = addTsys(tsys, inverseMask);
210 Float intervalWeight = addInterval(interval, inverseMask);
211 addTime(time, inverseMask);
212 if (!allEQ(m, False)) { // only add these if not everything masked
213 totalWeight *= tsysWeight;
214 totalWeight *= intervalWeight;
215
216 if (weightType_ == W_VAR) {
217 Float fac = 1.0/variance(data);
218 if (!inverseMask && (m.nelements() == userMask_.nelements()))
219 fac = 1.0/variance(data(userMask_));
220
221 totalWeight *= fac;
222 }
223 }
224
225 return totalWeight;
226}
227
228Float RowAccumulator::addTsys(const Vector<Float>& v, Bool inverseMask)
229{
230 // @fixme this assume tsys is the same for all channels
231
232 Float w = 1.0;
233 if (inverseMask) {
234 tsysSumNoMask_ += v[0];
235 } else {
236 tsysSum_ += v[0];
237 }
238 if ( weightType_ == W_TSYS || weightType_ == W_TINTSYS ) {
239 w /= (v[0]*v[0]);
240 }
241 return w;
242}
243
244void RowAccumulator::addTime(Double t, Bool inverseMask)
245{
246 if (inverseMask) {
247 timeSumNoMask_ += t;
248 } else {
249 timeSum_ += t;
250 }
251}
252
253Float RowAccumulator::addInterval(Double inter, Bool inverseMask)
254{
255 Float w = 1.0;
256 if (inverseMask) {
257 intervalSumNoMask_ += inter;
258 } else {
259 intervalSum_ += inter;
260 }
261 if (weightType_ == W_TINT || weightType_ == W_TINTSYS) {
262 w *= Float(inter);
263 }
264 return w;
265}
266
267Vector<Float> RowAccumulator::getSpectrum() const
268{
269 return (spectrum_/weightSum_).getArray();
270}
271
272Double RowAccumulator::getTime() const
273{
274 return timeSum_/Double(nNominal(n_, nNoMask_));
275}
276
277Double RowAccumulator::getInterval() const
278{
279 return intervalSum_;
280}
281
282Vector<Bool> RowAccumulator::getMask() const
283{
284 // Return the "total" mask - False where no points have been accumulated.
285 return (n_.getArray() > uInt(0));
286}
287
288Vector<Float> RowAccumulator::getTsys() const
289{
290 // @fixme this assumes tsys.nelements() == 1
291 return tsysSum_/Float(nNominal(n_, nNoMask_));
292}
293
294void RowAccumulator::setUserMask(const Vector<Bool>& m)
295{
296 userMask_.resize();
297 userMask_ = m;
298}
299
300// Added by TT check the state of RowAccumulator
301Bool RowAccumulator::state() const
302{
303 return initialized_;
304}
305
306void RowAccumulator::replaceNaN()
307{
308 Vector<Float> v = spectrum_.getArray();
309 Vector<Float> w = weightSum_.getArray();
310 Vector<Float> vRef = spectrumNoMask_.getArray();
311 Vector<Float> wRef = weightSumNoMask_.getArray();
312
313 for (uInt i = 0; i < v.size(); ++i) {
314 if (w[i] == 0.0) {
315 v[i] = vRef[i];
316 w[i] = wRef[i];
317 }
318 }
319
320 spectrum_.setData(v, Vector<Bool>(v.nelements(), True));
321 weightSum_.setData(w, Vector<Bool>(w.nelements(), True));
322
323 tsysSum_ = tsysSumNoMask_;
324 intervalSum_ = intervalSumNoMask_;
325}
Note: See TracBrowser for help on using the repository browser.