annotate armadillo-3.900.4/include/armadillo_bits/running_stat_vec_meat.hpp @ 84:55a047986812 tip

Update library URI so as not to be document-local
author Chris Cannam
date Wed, 22 Apr 2020 14:21:57 +0100
parents 1ec0e2823891
children
rev   line source
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 //! @}