Chris@49
|
1 // Copyright (C) 2009-2011 NICTA (www.nicta.com.au)
|
Chris@49
|
2 // Copyright (C) 2009-2011 Conrad Sanderson
|
Chris@49
|
3 //
|
Chris@49
|
4 // This Source Code Form is subject to the terms of the Mozilla Public
|
Chris@49
|
5 // License, v. 2.0. If a copy of the MPL was not distributed with this
|
Chris@49
|
6 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
Chris@49
|
7
|
Chris@49
|
8
|
Chris@49
|
9 //! \addtogroup running_stat_vec
|
Chris@49
|
10 //! @{
|
Chris@49
|
11
|
Chris@49
|
12
|
Chris@49
|
13
|
Chris@49
|
14 template<typename eT>
|
Chris@49
|
15 running_stat_vec<eT>::~running_stat_vec()
|
Chris@49
|
16 {
|
Chris@49
|
17 arma_extra_debug_sigprint_this(this);
|
Chris@49
|
18 }
|
Chris@49
|
19
|
Chris@49
|
20
|
Chris@49
|
21
|
Chris@49
|
22 template<typename eT>
|
Chris@49
|
23 running_stat_vec<eT>::running_stat_vec(const bool in_calc_cov)
|
Chris@49
|
24 : calc_cov(in_calc_cov)
|
Chris@49
|
25 {
|
Chris@49
|
26 arma_extra_debug_sigprint_this(this);
|
Chris@49
|
27 }
|
Chris@49
|
28
|
Chris@49
|
29
|
Chris@49
|
30
|
Chris@49
|
31 template<typename eT>
|
Chris@49
|
32 running_stat_vec<eT>::running_stat_vec(const running_stat_vec<eT>& in_rsv)
|
Chris@49
|
33 : calc_cov (in_rsv.calc_cov)
|
Chris@49
|
34 , counter (in_rsv.counter)
|
Chris@49
|
35 , r_mean (in_rsv.r_mean)
|
Chris@49
|
36 , r_var (in_rsv.r_var)
|
Chris@49
|
37 , r_cov (in_rsv.r_cov)
|
Chris@49
|
38 , min_val (in_rsv.min_val)
|
Chris@49
|
39 , max_val (in_rsv.max_val)
|
Chris@49
|
40 , min_val_norm(in_rsv.min_val_norm)
|
Chris@49
|
41 , max_val_norm(in_rsv.max_val_norm)
|
Chris@49
|
42 {
|
Chris@49
|
43 arma_extra_debug_sigprint_this(this);
|
Chris@49
|
44 }
|
Chris@49
|
45
|
Chris@49
|
46
|
Chris@49
|
47
|
Chris@49
|
48 template<typename eT>
|
Chris@49
|
49 const running_stat_vec<eT>&
|
Chris@49
|
50 running_stat_vec<eT>::operator=(const running_stat_vec<eT>& in_rsv)
|
Chris@49
|
51 {
|
Chris@49
|
52 arma_extra_debug_sigprint();
|
Chris@49
|
53
|
Chris@49
|
54 access::rw(calc_cov) = in_rsv.calc_cov;
|
Chris@49
|
55
|
Chris@49
|
56 counter = in_rsv.counter;
|
Chris@49
|
57 r_mean = in_rsv.r_mean;
|
Chris@49
|
58 r_var = in_rsv.r_var;
|
Chris@49
|
59 r_cov = in_rsv.r_cov;
|
Chris@49
|
60 min_val = in_rsv.min_val;
|
Chris@49
|
61 max_val = in_rsv.max_val;
|
Chris@49
|
62 min_val_norm = in_rsv.min_val_norm;
|
Chris@49
|
63 max_val_norm = in_rsv.max_val_norm;
|
Chris@49
|
64
|
Chris@49
|
65 return *this;
|
Chris@49
|
66 }
|
Chris@49
|
67
|
Chris@49
|
68
|
Chris@49
|
69
|
Chris@49
|
70 //! update statistics to reflect new sample
|
Chris@49
|
71 template<typename eT>
|
Chris@49
|
72 template<typename T1>
|
Chris@49
|
73 arma_hot
|
Chris@49
|
74 inline
|
Chris@49
|
75 void
|
Chris@49
|
76 running_stat_vec<eT>::operator() (const Base<typename get_pod_type<eT>::result, T1>& X)
|
Chris@49
|
77 {
|
Chris@49
|
78 arma_extra_debug_sigprint();
|
Chris@49
|
79
|
Chris@49
|
80 //typedef typename get_pod_type<eT>::result T;
|
Chris@49
|
81
|
Chris@49
|
82 const unwrap<T1> tmp(X.get_ref());
|
Chris@49
|
83 const Mat<eT>& sample = tmp.M;
|
Chris@49
|
84
|
Chris@49
|
85 if( sample.is_empty() )
|
Chris@49
|
86 {
|
Chris@49
|
87 return;
|
Chris@49
|
88 }
|
Chris@49
|
89
|
Chris@49
|
90 if( sample.is_finite() == false )
|
Chris@49
|
91 {
|
Chris@49
|
92 arma_warn(true, "running_stat_vec: sample ignored as it has non-finite elements");
|
Chris@49
|
93 return;
|
Chris@49
|
94 }
|
Chris@49
|
95
|
Chris@49
|
96 running_stat_vec_aux::update_stats(*this, sample);
|
Chris@49
|
97 }
|
Chris@49
|
98
|
Chris@49
|
99
|
Chris@49
|
100
|
Chris@49
|
101 //! update statistics to reflect new sample (version for complex numbers)
|
Chris@49
|
102 template<typename eT>
|
Chris@49
|
103 template<typename T1>
|
Chris@49
|
104 arma_hot
|
Chris@49
|
105 inline
|
Chris@49
|
106 void
|
Chris@49
|
107 running_stat_vec<eT>::operator() (const Base<std::complex<typename get_pod_type<eT>::result>, T1>& X)
|
Chris@49
|
108 {
|
Chris@49
|
109 arma_extra_debug_sigprint();
|
Chris@49
|
110
|
Chris@49
|
111 //typedef typename std::complex<typename get_pod_type<eT>::result> eT;
|
Chris@49
|
112
|
Chris@49
|
113 const unwrap<T1> tmp(X.get_ref());
|
Chris@49
|
114 const Mat<eT>& sample = tmp.M;
|
Chris@49
|
115
|
Chris@49
|
116 if( sample.is_empty() )
|
Chris@49
|
117 {
|
Chris@49
|
118 return;
|
Chris@49
|
119 }
|
Chris@49
|
120
|
Chris@49
|
121 if( sample.is_finite() == false )
|
Chris@49
|
122 {
|
Chris@49
|
123 arma_warn(true, "running_stat_vec: sample ignored as it has non-finite elements");
|
Chris@49
|
124 return;
|
Chris@49
|
125 }
|
Chris@49
|
126
|
Chris@49
|
127 running_stat_vec_aux::update_stats(*this, sample);
|
Chris@49
|
128 }
|
Chris@49
|
129
|
Chris@49
|
130
|
Chris@49
|
131
|
Chris@49
|
132 //! set all statistics to zero
|
Chris@49
|
133 template<typename eT>
|
Chris@49
|
134 inline
|
Chris@49
|
135 void
|
Chris@49
|
136 running_stat_vec<eT>::reset()
|
Chris@49
|
137 {
|
Chris@49
|
138 arma_extra_debug_sigprint();
|
Chris@49
|
139
|
Chris@49
|
140 counter.reset();
|
Chris@49
|
141
|
Chris@49
|
142 r_mean.reset();
|
Chris@49
|
143 r_var.reset();
|
Chris@49
|
144 r_cov.reset();
|
Chris@49
|
145
|
Chris@49
|
146 min_val.reset();
|
Chris@49
|
147 max_val.reset();
|
Chris@49
|
148
|
Chris@49
|
149 min_val_norm.reset();
|
Chris@49
|
150 max_val_norm.reset();
|
Chris@49
|
151
|
Chris@49
|
152 r_var_dummy.reset();
|
Chris@49
|
153 r_cov_dummy.reset();
|
Chris@49
|
154
|
Chris@49
|
155 tmp1.reset();
|
Chris@49
|
156 tmp2.reset();
|
Chris@49
|
157 }
|
Chris@49
|
158
|
Chris@49
|
159
|
Chris@49
|
160
|
Chris@49
|
161 //! mean or average value
|
Chris@49
|
162 template<typename eT>
|
Chris@49
|
163 inline
|
Chris@49
|
164 const Mat<eT>&
|
Chris@49
|
165 running_stat_vec<eT>::mean() const
|
Chris@49
|
166 {
|
Chris@49
|
167 arma_extra_debug_sigprint();
|
Chris@49
|
168
|
Chris@49
|
169 return r_mean;
|
Chris@49
|
170 }
|
Chris@49
|
171
|
Chris@49
|
172
|
Chris@49
|
173
|
Chris@49
|
174 //! variance
|
Chris@49
|
175 template<typename eT>
|
Chris@49
|
176 inline
|
Chris@49
|
177 const Mat<typename get_pod_type<eT>::result>&
|
Chris@49
|
178 running_stat_vec<eT>::var(const uword norm_type)
|
Chris@49
|
179 {
|
Chris@49
|
180 arma_extra_debug_sigprint();
|
Chris@49
|
181
|
Chris@49
|
182 const T N = counter.value();
|
Chris@49
|
183
|
Chris@49
|
184 if(N > T(1))
|
Chris@49
|
185 {
|
Chris@49
|
186 if(norm_type == 0)
|
Chris@49
|
187 {
|
Chris@49
|
188 return r_var;
|
Chris@49
|
189 }
|
Chris@49
|
190 else
|
Chris@49
|
191 {
|
Chris@49
|
192 const T N_minus_1 = counter.value_minus_1();
|
Chris@49
|
193
|
Chris@49
|
194 r_var_dummy = (N_minus_1/N) * r_var;
|
Chris@49
|
195
|
Chris@49
|
196 return r_var_dummy;
|
Chris@49
|
197 }
|
Chris@49
|
198 }
|
Chris@49
|
199 else
|
Chris@49
|
200 {
|
Chris@49
|
201 r_var_dummy.zeros(r_mean.n_rows, r_mean.n_cols);
|
Chris@49
|
202
|
Chris@49
|
203 return r_var_dummy;
|
Chris@49
|
204 }
|
Chris@49
|
205
|
Chris@49
|
206 }
|
Chris@49
|
207
|
Chris@49
|
208
|
Chris@49
|
209
|
Chris@49
|
210 //! standard deviation
|
Chris@49
|
211 template<typename eT>
|
Chris@49
|
212 inline
|
Chris@49
|
213 Mat<typename get_pod_type<eT>::result>
|
Chris@49
|
214 running_stat_vec<eT>::stddev(const uword norm_type) const
|
Chris@49
|
215 {
|
Chris@49
|
216 arma_extra_debug_sigprint();
|
Chris@49
|
217
|
Chris@49
|
218 const T N = counter.value();
|
Chris@49
|
219
|
Chris@49
|
220 if(N > T(1))
|
Chris@49
|
221 {
|
Chris@49
|
222 if(norm_type == 0)
|
Chris@49
|
223 {
|
Chris@49
|
224 return sqrt(r_var);
|
Chris@49
|
225 }
|
Chris@49
|
226 else
|
Chris@49
|
227 {
|
Chris@49
|
228 const T N_minus_1 = counter.value_minus_1();
|
Chris@49
|
229
|
Chris@49
|
230 return sqrt( (N_minus_1/N) * r_var );
|
Chris@49
|
231 }
|
Chris@49
|
232 }
|
Chris@49
|
233 else
|
Chris@49
|
234 {
|
Chris@49
|
235 return Mat<T>();
|
Chris@49
|
236 }
|
Chris@49
|
237 }
|
Chris@49
|
238
|
Chris@49
|
239
|
Chris@49
|
240
|
Chris@49
|
241 //! covariance
|
Chris@49
|
242 template<typename eT>
|
Chris@49
|
243 inline
|
Chris@49
|
244 const Mat<eT>&
|
Chris@49
|
245 running_stat_vec<eT>::cov(const uword norm_type)
|
Chris@49
|
246 {
|
Chris@49
|
247 arma_extra_debug_sigprint();
|
Chris@49
|
248
|
Chris@49
|
249 if(calc_cov == true)
|
Chris@49
|
250 {
|
Chris@49
|
251 const T N = counter.value();
|
Chris@49
|
252
|
Chris@49
|
253 if(N > T(1))
|
Chris@49
|
254 {
|
Chris@49
|
255 if(norm_type == 0)
|
Chris@49
|
256 {
|
Chris@49
|
257 return r_cov;
|
Chris@49
|
258 }
|
Chris@49
|
259 else
|
Chris@49
|
260 {
|
Chris@49
|
261 const T N_minus_1 = counter.value_minus_1();
|
Chris@49
|
262
|
Chris@49
|
263 r_cov_dummy = (N_minus_1/N) * r_cov;
|
Chris@49
|
264
|
Chris@49
|
265 return r_cov_dummy;
|
Chris@49
|
266 }
|
Chris@49
|
267 }
|
Chris@49
|
268 else
|
Chris@49
|
269 {
|
Chris@49
|
270 r_cov_dummy.zeros(r_mean.n_rows, r_mean.n_cols);
|
Chris@49
|
271
|
Chris@49
|
272 return r_cov_dummy;
|
Chris@49
|
273 }
|
Chris@49
|
274 }
|
Chris@49
|
275 else
|
Chris@49
|
276 {
|
Chris@49
|
277 r_cov_dummy.reset();
|
Chris@49
|
278
|
Chris@49
|
279 return r_cov_dummy;
|
Chris@49
|
280 }
|
Chris@49
|
281
|
Chris@49
|
282 }
|
Chris@49
|
283
|
Chris@49
|
284
|
Chris@49
|
285
|
Chris@49
|
286 //! vector with minimum values
|
Chris@49
|
287 template<typename eT>
|
Chris@49
|
288 inline
|
Chris@49
|
289 const Mat<eT>&
|
Chris@49
|
290 running_stat_vec<eT>::min() const
|
Chris@49
|
291 {
|
Chris@49
|
292 arma_extra_debug_sigprint();
|
Chris@49
|
293
|
Chris@49
|
294 return min_val;
|
Chris@49
|
295 }
|
Chris@49
|
296
|
Chris@49
|
297
|
Chris@49
|
298
|
Chris@49
|
299 //! vector with maximum values
|
Chris@49
|
300 template<typename eT>
|
Chris@49
|
301 inline
|
Chris@49
|
302 const Mat<eT>&
|
Chris@49
|
303 running_stat_vec<eT>::max() const
|
Chris@49
|
304 {
|
Chris@49
|
305 arma_extra_debug_sigprint();
|
Chris@49
|
306
|
Chris@49
|
307 return max_val;
|
Chris@49
|
308 }
|
Chris@49
|
309
|
Chris@49
|
310
|
Chris@49
|
311
|
Chris@49
|
312 //! number of samples so far
|
Chris@49
|
313 template<typename eT>
|
Chris@49
|
314 inline
|
Chris@49
|
315 typename get_pod_type<eT>::result
|
Chris@49
|
316 running_stat_vec<eT>::count() const
|
Chris@49
|
317 {
|
Chris@49
|
318 arma_extra_debug_sigprint();
|
Chris@49
|
319
|
Chris@49
|
320 return counter.value();
|
Chris@49
|
321 }
|
Chris@49
|
322
|
Chris@49
|
323
|
Chris@49
|
324
|
Chris@49
|
325 //
|
Chris@49
|
326
|
Chris@49
|
327
|
Chris@49
|
328
|
Chris@49
|
329 //! update statistics to reflect new sample
|
Chris@49
|
330 template<typename eT>
|
Chris@49
|
331 inline
|
Chris@49
|
332 void
|
Chris@49
|
333 running_stat_vec_aux::update_stats(running_stat_vec<eT>& x, const Mat<eT>& sample)
|
Chris@49
|
334 {
|
Chris@49
|
335 arma_extra_debug_sigprint();
|
Chris@49
|
336
|
Chris@49
|
337 typedef typename running_stat_vec<eT>::T T;
|
Chris@49
|
338
|
Chris@49
|
339 const T N = x.counter.value();
|
Chris@49
|
340
|
Chris@49
|
341 if(N > T(0))
|
Chris@49
|
342 {
|
Chris@49
|
343 arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch");
|
Chris@49
|
344
|
Chris@49
|
345 const uword n_elem = sample.n_elem;
|
Chris@49
|
346 const eT* sample_mem = sample.memptr();
|
Chris@49
|
347 eT* r_mean_mem = x.r_mean.memptr();
|
Chris@49
|
348 T* r_var_mem = x.r_var.memptr();
|
Chris@49
|
349 eT* min_val_mem = x.min_val.memptr();
|
Chris@49
|
350 eT* max_val_mem = x.max_val.memptr();
|
Chris@49
|
351
|
Chris@49
|
352 const T N_plus_1 = x.counter.value_plus_1();
|
Chris@49
|
353 const T N_minus_1 = x.counter.value_minus_1();
|
Chris@49
|
354
|
Chris@49
|
355 if(x.calc_cov == true)
|
Chris@49
|
356 {
|
Chris@49
|
357 Mat<eT>& tmp1 = x.tmp1;
|
Chris@49
|
358 Mat<eT>& tmp2 = x.tmp2;
|
Chris@49
|
359
|
Chris@49
|
360 tmp1 = sample - x.r_mean;
|
Chris@49
|
361
|
Chris@49
|
362 if(sample.n_cols == 1)
|
Chris@49
|
363 {
|
Chris@49
|
364 tmp2 = tmp1*trans(tmp1);
|
Chris@49
|
365 }
|
Chris@49
|
366 else
|
Chris@49
|
367 {
|
Chris@49
|
368 tmp2 = trans(tmp1)*tmp1;
|
Chris@49
|
369 }
|
Chris@49
|
370
|
Chris@49
|
371 x.r_cov *= (N_minus_1/N);
|
Chris@49
|
372 x.r_cov += tmp2 / N_plus_1;
|
Chris@49
|
373 }
|
Chris@49
|
374
|
Chris@49
|
375
|
Chris@49
|
376 for(uword i=0; i<n_elem; ++i)
|
Chris@49
|
377 {
|
Chris@49
|
378 const eT val = sample_mem[i];
|
Chris@49
|
379
|
Chris@49
|
380 if(val < min_val_mem[i])
|
Chris@49
|
381 {
|
Chris@49
|
382 min_val_mem[i] = val;
|
Chris@49
|
383 }
|
Chris@49
|
384
|
Chris@49
|
385 if(val > max_val_mem[i])
|
Chris@49
|
386 {
|
Chris@49
|
387 max_val_mem[i] = val;
|
Chris@49
|
388 }
|
Chris@49
|
389
|
Chris@49
|
390 const eT r_mean_val = r_mean_mem[i];
|
Chris@49
|
391 const eT tmp = val - r_mean_val;
|
Chris@49
|
392
|
Chris@49
|
393 r_var_mem[i] = N_minus_1/N * r_var_mem[i] + (tmp*tmp)/N_plus_1;
|
Chris@49
|
394
|
Chris@49
|
395 r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1;
|
Chris@49
|
396 }
|
Chris@49
|
397 }
|
Chris@49
|
398 else
|
Chris@49
|
399 {
|
Chris@49
|
400 arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector");
|
Chris@49
|
401
|
Chris@49
|
402 x.r_mean.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
403
|
Chris@49
|
404 x.r_var.zeros(sample.n_rows, sample.n_cols);
|
Chris@49
|
405
|
Chris@49
|
406 if(x.calc_cov == true)
|
Chris@49
|
407 {
|
Chris@49
|
408 x.r_cov.zeros(sample.n_elem, sample.n_elem);
|
Chris@49
|
409 }
|
Chris@49
|
410
|
Chris@49
|
411 x.min_val.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
412 x.max_val.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
413
|
Chris@49
|
414
|
Chris@49
|
415 const uword n_elem = sample.n_elem;
|
Chris@49
|
416 const eT* sample_mem = sample.memptr();
|
Chris@49
|
417 eT* r_mean_mem = x.r_mean.memptr();
|
Chris@49
|
418 eT* min_val_mem = x.min_val.memptr();
|
Chris@49
|
419 eT* max_val_mem = x.max_val.memptr();
|
Chris@49
|
420
|
Chris@49
|
421
|
Chris@49
|
422 for(uword i=0; i<n_elem; ++i)
|
Chris@49
|
423 {
|
Chris@49
|
424 const eT val = sample_mem[i];
|
Chris@49
|
425
|
Chris@49
|
426 r_mean_mem[i] = val;
|
Chris@49
|
427 min_val_mem[i] = val;
|
Chris@49
|
428 max_val_mem[i] = val;
|
Chris@49
|
429 }
|
Chris@49
|
430 }
|
Chris@49
|
431
|
Chris@49
|
432 x.counter++;
|
Chris@49
|
433 }
|
Chris@49
|
434
|
Chris@49
|
435
|
Chris@49
|
436
|
Chris@49
|
437 //! update statistics to reflect new sample (version for complex numbers)
|
Chris@49
|
438 template<typename T>
|
Chris@49
|
439 inline
|
Chris@49
|
440 void
|
Chris@49
|
441 running_stat_vec_aux::update_stats(running_stat_vec< std::complex<T> >& x, const Mat<T>& sample)
|
Chris@49
|
442 {
|
Chris@49
|
443 arma_extra_debug_sigprint();
|
Chris@49
|
444
|
Chris@49
|
445 const Mat< std::complex<T> > tmp = conv_to< Mat< std::complex<T> > >::from(sample);
|
Chris@49
|
446
|
Chris@49
|
447 running_stat_vec_aux::update_stats(x, tmp);
|
Chris@49
|
448 }
|
Chris@49
|
449
|
Chris@49
|
450
|
Chris@49
|
451
|
Chris@49
|
452 //! alter statistics to reflect new sample (version for complex numbers)
|
Chris@49
|
453 template<typename T>
|
Chris@49
|
454 inline
|
Chris@49
|
455 void
|
Chris@49
|
456 running_stat_vec_aux::update_stats(running_stat_vec< std::complex<T> >& x, const Mat< std::complex<T> >& sample)
|
Chris@49
|
457 {
|
Chris@49
|
458 arma_extra_debug_sigprint();
|
Chris@49
|
459
|
Chris@49
|
460 typedef typename std::complex<T> eT;
|
Chris@49
|
461
|
Chris@49
|
462 const T N = x.counter.value();
|
Chris@49
|
463
|
Chris@49
|
464 if(N > T(0))
|
Chris@49
|
465 {
|
Chris@49
|
466 arma_debug_assert_same_size(x.r_mean, sample, "running_stat_vec(): dimensionality mismatch");
|
Chris@49
|
467
|
Chris@49
|
468 const uword n_elem = sample.n_elem;
|
Chris@49
|
469 const eT* sample_mem = sample.memptr();
|
Chris@49
|
470 eT* r_mean_mem = x.r_mean.memptr();
|
Chris@49
|
471 T* r_var_mem = x.r_var.memptr();
|
Chris@49
|
472 eT* min_val_mem = x.min_val.memptr();
|
Chris@49
|
473 eT* max_val_mem = x.max_val.memptr();
|
Chris@49
|
474 T* min_val_norm_mem = x.min_val_norm.memptr();
|
Chris@49
|
475 T* max_val_norm_mem = x.max_val_norm.memptr();
|
Chris@49
|
476
|
Chris@49
|
477 const T N_plus_1 = x.counter.value_plus_1();
|
Chris@49
|
478 const T N_minus_1 = x.counter.value_minus_1();
|
Chris@49
|
479
|
Chris@49
|
480 if(x.calc_cov == true)
|
Chris@49
|
481 {
|
Chris@49
|
482 Mat<eT>& tmp1 = x.tmp1;
|
Chris@49
|
483 Mat<eT>& tmp2 = x.tmp2;
|
Chris@49
|
484
|
Chris@49
|
485 tmp1 = sample - x.r_mean;
|
Chris@49
|
486
|
Chris@49
|
487 if(sample.n_cols == 1)
|
Chris@49
|
488 {
|
Chris@49
|
489 tmp2 = arma::conj(tmp1)*strans(tmp1);
|
Chris@49
|
490 }
|
Chris@49
|
491 else
|
Chris@49
|
492 {
|
Chris@49
|
493 tmp2 = trans(tmp1)*tmp1; //tmp2 = strans(conj(tmp1))*tmp1;
|
Chris@49
|
494 }
|
Chris@49
|
495
|
Chris@49
|
496 x.r_cov *= (N_minus_1/N);
|
Chris@49
|
497 x.r_cov += tmp2 / N_plus_1;
|
Chris@49
|
498 }
|
Chris@49
|
499
|
Chris@49
|
500
|
Chris@49
|
501 for(uword i=0; i<n_elem; ++i)
|
Chris@49
|
502 {
|
Chris@49
|
503 const eT& val = sample_mem[i];
|
Chris@49
|
504 const T val_norm = std::norm(val);
|
Chris@49
|
505
|
Chris@49
|
506 if(val_norm < min_val_norm_mem[i])
|
Chris@49
|
507 {
|
Chris@49
|
508 min_val_norm_mem[i] = val_norm;
|
Chris@49
|
509 min_val_mem[i] = val;
|
Chris@49
|
510 }
|
Chris@49
|
511
|
Chris@49
|
512 if(val_norm > max_val_norm_mem[i])
|
Chris@49
|
513 {
|
Chris@49
|
514 max_val_norm_mem[i] = val_norm;
|
Chris@49
|
515 max_val_mem[i] = val;
|
Chris@49
|
516 }
|
Chris@49
|
517
|
Chris@49
|
518 const eT& r_mean_val = r_mean_mem[i];
|
Chris@49
|
519
|
Chris@49
|
520 r_var_mem[i] = N_minus_1/N * r_var_mem[i] + std::norm(val - r_mean_val)/N_plus_1;
|
Chris@49
|
521
|
Chris@49
|
522 r_mean_mem[i] = r_mean_val + (val - r_mean_val)/N_plus_1;
|
Chris@49
|
523 }
|
Chris@49
|
524
|
Chris@49
|
525 }
|
Chris@49
|
526 else
|
Chris@49
|
527 {
|
Chris@49
|
528 arma_debug_check( (sample.is_vec() == false), "running_stat_vec(): given sample is not a vector");
|
Chris@49
|
529
|
Chris@49
|
530 x.r_mean.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
531
|
Chris@49
|
532 x.r_var.zeros(sample.n_rows, sample.n_cols);
|
Chris@49
|
533
|
Chris@49
|
534 if(x.calc_cov == true)
|
Chris@49
|
535 {
|
Chris@49
|
536 x.r_cov.zeros(sample.n_elem, sample.n_elem);
|
Chris@49
|
537 }
|
Chris@49
|
538
|
Chris@49
|
539 x.min_val.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
540 x.max_val.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
541
|
Chris@49
|
542 x.min_val_norm.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
543 x.max_val_norm.set_size(sample.n_rows, sample.n_cols);
|
Chris@49
|
544
|
Chris@49
|
545
|
Chris@49
|
546 const uword n_elem = sample.n_elem;
|
Chris@49
|
547 const eT* sample_mem = sample.memptr();
|
Chris@49
|
548 eT* r_mean_mem = x.r_mean.memptr();
|
Chris@49
|
549 eT* min_val_mem = x.min_val.memptr();
|
Chris@49
|
550 eT* max_val_mem = x.max_val.memptr();
|
Chris@49
|
551 T* min_val_norm_mem = x.min_val_norm.memptr();
|
Chris@49
|
552 T* max_val_norm_mem = x.max_val_norm.memptr();
|
Chris@49
|
553
|
Chris@49
|
554 for(uword i=0; i<n_elem; ++i)
|
Chris@49
|
555 {
|
Chris@49
|
556 const eT& val = sample_mem[i];
|
Chris@49
|
557 const T val_norm = std::norm(val);
|
Chris@49
|
558
|
Chris@49
|
559 r_mean_mem[i] = val;
|
Chris@49
|
560 min_val_mem[i] = val;
|
Chris@49
|
561 max_val_mem[i] = val;
|
Chris@49
|
562
|
Chris@49
|
563 min_val_norm_mem[i] = val_norm;
|
Chris@49
|
564 max_val_norm_mem[i] = val_norm;
|
Chris@49
|
565 }
|
Chris@49
|
566 }
|
Chris@49
|
567
|
Chris@49
|
568 x.counter++;
|
Chris@49
|
569 }
|
Chris@49
|
570
|
Chris@49
|
571
|
Chris@49
|
572
|
Chris@49
|
573 //! @}
|