source: trunk/src/RowAccumulator.cpp@ 3030

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

New Development: No

JIRA Issue: Yes CAS-6582

Ready for Test: Yes

Interface Changes: Yes/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...

Fix for asap unit test failure.


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 if (!allEQ(m, False)) { // only add these if not everything masked
210 totalWeight *= addTsys(tsys, inverseMask);
211 totalWeight *= addInterval(interval, inverseMask);
212 addTime(time, inverseMask);
213
214 if (weightType_ == W_VAR) {
215 Float fac = 1.0/variance(data);
216 if (!inverseMask && (m.nelements() == userMask_.nelements()))
217 fac = 1.0/variance(data(userMask_));
218
219 totalWeight *= fac;
220 }
221 }
222
223 return totalWeight;
224}
225
226Float RowAccumulator::addTsys(const Vector<Float>& v, Bool inverseMask)
227{
228 // @fixme this assume tsys is the same for all channels
229
230 Float w = 1.0;
231 if (inverseMask) {
232 tsysSumNoMask_ += v[0];
233 } else {
234 tsysSum_ += v[0];
235 }
236 if ( weightType_ == W_TSYS || weightType_ == W_TINTSYS ) {
237 w /= (v[0]*v[0]);
238 }
239 return w;
240}
241
242void RowAccumulator::addTime(Double t, Bool inverseMask)
243{
244 if (inverseMask) {
245 timeSumNoMask_ += t;
246 } else {
247 timeSum_ += t;
248 }
249}
250
251Float RowAccumulator::addInterval(Double inter, Bool inverseMask)
252{
253 Float w = 1.0;
254 if (inverseMask) {
255 intervalSumNoMask_ += inter;
256 } else {
257 intervalSum_ += inter;
258 }
259 if (weightType_ == W_TINT || weightType_ == W_TINTSYS) {
260 w *= Float(inter);
261 }
262 return w;
263}
264
265Vector<Float> RowAccumulator::getSpectrum() const
266{
267 return (spectrum_/weightSum_).getArray();
268}
269
270Double RowAccumulator::getTime() const
271{
272 return timeSum_/Double(nNominal(n_, nNoMask_));
273}
274
275Double RowAccumulator::getInterval() const
276{
277 return intervalSum_;
278}
279
280Vector<Bool> RowAccumulator::getMask() const
281{
282 // Return the "total" mask - False where no points have been accumulated.
283 return (n_.getArray() > uInt(0));
284}
285
286Vector<Float> RowAccumulator::getTsys() const
287{
288 // @fixme this assumes tsys.nelements() == 1
289 return tsysSum_/Float(nNominal(n_, nNoMask_));
290}
291
292void RowAccumulator::setUserMask(const Vector<Bool>& m)
293{
294 userMask_.resize();
295 userMask_ = m;
296}
297
298// Added by TT check the state of RowAccumulator
299Bool RowAccumulator::state() const
300{
301 return initialized_;
302}
303
304void RowAccumulator::replaceNaN()
305{
306 Vector<Float> v = spectrum_.getArray();
307 Vector<Float> w = weightSum_.getArray();
308 Vector<Float> vRef = spectrumNoMask_.getArray();
309 Vector<Float> wRef = weightSumNoMask_.getArray();
310
311 for (uInt i = 0; i < v.size(); ++i) {
312 if (w[i] == 0.0) {
313 v[i] = vRef[i];
314 w[i] = wRef[i];
315 }
316 }
317
318 spectrum_.setData(v, Vector<Bool>(v.nelements(), True));
319 weightSum_.setData(w, Vector<Bool>(w.nelements(), True));
320
321 if (allEQ(tsysSum_, 0.0f))
322 tsysSum_ = tsysSumNoMask_;
323 if (intervalSum_ == 0.0)
324 intervalSum_ = intervalSumNoMask_;
325}
Note: See TracBrowser for help on using the repository browser.