source: trunk/src/RowAccumulator.cpp@ 2693

Last change on this file since 2693 was 2580, checked in by ShinnosukeKawakami, 12 years ago

hpc33 merged asap-trunk

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