annotate armadillo-2.4.4/include/armadillo_bits/running_stat_vec_meat.hpp @ 0:8b6102e2a9b0

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