From eda5bc26f44ee9a6f83dcf8c91f17296d7fc509d Mon Sep 17 00:00:00 2001 From: Nao Pross Date: Mon, 12 Feb 2024 14:52:43 +0100 Subject: Move into version control --- .../include/armadillo_bits/running_stat_meat.hpp | 463 +++++++++++++++++++++ 1 file changed, 463 insertions(+) create mode 100644 src/armadillo/include/armadillo_bits/running_stat_meat.hpp (limited to 'src/armadillo/include/armadillo_bits/running_stat_meat.hpp') diff --git a/src/armadillo/include/armadillo_bits/running_stat_meat.hpp b/src/armadillo/include/armadillo_bits/running_stat_meat.hpp new file mode 100644 index 0000000..35e6ba8 --- /dev/null +++ b/src/armadillo/include/armadillo_bits/running_stat_meat.hpp @@ -0,0 +1,463 @@ +// SPDX-License-Identifier: Apache-2.0 +// +// Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au) +// Copyright 2008-2016 National ICT Australia (NICTA) +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. +// ------------------------------------------------------------------------ + + +//! \addtogroup running_stat +//! @{ + + + +template +inline +arma_counter::~arma_counter() + { + arma_extra_debug_sigprint_this(this); + } + + + +template +inline +arma_counter::arma_counter() + : d_count( eT(0)) + , i_count(uword(0)) + { + arma_extra_debug_sigprint_this(this); + } + + + +template +inline +const arma_counter& +arma_counter::operator++() + { + if(i_count < ARMA_MAX_UWORD) + { + i_count++; + } + else + { + d_count += eT(ARMA_MAX_UWORD); + i_count = 1; + } + + return *this; + } + + + +template +inline +void +arma_counter::operator++(int) + { + operator++(); + } + + + +template +inline +void +arma_counter::reset() + { + d_count = eT(0); + i_count = uword(0); + } + + + +template +inline +eT +arma_counter::value() const + { + return d_count + eT(i_count); + } + + + +template +inline +eT +arma_counter::value_plus_1() const + { + if(i_count < ARMA_MAX_UWORD) + { + return d_count + eT(i_count + 1); + } + else + { + return d_count + eT(ARMA_MAX_UWORD) + eT(1); + } + } + + + +template +inline +eT +arma_counter::value_minus_1() const + { + if(i_count > 0) + { + return d_count + eT(i_count - 1); + } + else + { + return d_count - eT(1); + } + } + + + +// + + + +template +inline +running_stat::~running_stat() + { + arma_extra_debug_sigprint_this(this); + } + + + +template +inline +running_stat::running_stat() + : r_mean ( eT(0)) + , r_var (typename running_stat::T(0)) + , min_val ( eT(0)) + , max_val ( eT(0)) + , min_val_norm(typename running_stat::T(0)) + , max_val_norm(typename running_stat::T(0)) + { + arma_extra_debug_sigprint_this(this); + } + + + +//! update statistics to reflect new sample +template +inline +void +running_stat::operator() (const typename running_stat::T sample) + { + arma_extra_debug_sigprint(); + + if( arma_isfinite(sample) == false ) + { + arma_debug_warn_level(3, "running_stat: sample ignored as it is non-finite" ); + return; + } + + running_stat_aux::update_stats(*this, sample); + } + + + +//! update statistics to reflect new sample (version for complex numbers) +template +inline +void +running_stat::operator() (const std::complex< typename running_stat::T >& sample) + { + arma_extra_debug_sigprint(); + + if( arma_isfinite(sample) == false ) + { + arma_debug_warn_level(3, "running_stat: sample ignored as it is non-finite" ); + return; + } + + running_stat_aux::update_stats(*this, sample); + } + + + +//! set all statistics to zero +template +inline +void +running_stat::reset() + { + arma_extra_debug_sigprint(); + + // typedef typename running_stat::T T; + + counter.reset(); + + r_mean = eT(0); + r_var = T(0); + + min_val = eT(0); + max_val = eT(0); + + min_val_norm = T(0); + max_val_norm = T(0); + } + + + +//! mean or average value +template +inline +eT +running_stat::mean() const + { + arma_extra_debug_sigprint(); + + return r_mean; + } + + + +//! variance +template +inline +typename running_stat::T +running_stat::var(const uword norm_type) const + { + arma_extra_debug_sigprint(); + + const T N = counter.value(); + + if(N > T(1)) + { + if(norm_type == 0) + { + return r_var; + } + else + { + const T N_minus_1 = counter.value_minus_1(); + return (N_minus_1/N) * r_var; + } + } + else + { + return T(0); + } + } + + + +//! standard deviation +template +inline +typename running_stat::T +running_stat::stddev(const uword norm_type) const + { + arma_extra_debug_sigprint(); + + return std::sqrt( (*this).var(norm_type) ); + } + + + +//! minimum value +template +inline +eT +running_stat::min() const + { + arma_extra_debug_sigprint(); + + return min_val; + } + + + +//! maximum value +template +inline +eT +running_stat::max() const + { + arma_extra_debug_sigprint(); + + return max_val; + } + + + +template +inline +eT +running_stat::range() const + { + arma_extra_debug_sigprint(); + + return (max_val - min_val); + } + + + +//! number of samples so far +template +inline +typename get_pod_type::result +running_stat::count() const + { + arma_extra_debug_sigprint(); + + return counter.value(); + } + + + +//! update statistics to reflect new sample (version for non-complex numbers, non-complex sample) +template +inline +void +running_stat_aux::update_stats(running_stat& x, const eT sample, const typename arma_not_cx::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + typedef typename running_stat::T T; + + const T N = x.counter.value(); + + if(N > T(0)) + { + if(sample < x.min_val) + { + x.min_val = sample; + } + + if(sample > x.max_val) + { + x.max_val = sample; + } + + const T N_plus_1 = x.counter.value_plus_1(); + const T N_minus_1 = x.counter.value_minus_1(); + + // note: variance has to be updated before the mean + + const eT tmp = sample - x.r_mean; + + x.r_var = N_minus_1/N * x.r_var + (tmp*tmp)/N_plus_1; + + x.r_mean = x.r_mean + (sample - x.r_mean)/N_plus_1; + //x.r_mean = (N/N_plus_1)*x.r_mean + sample/N_plus_1; + //x.r_mean = (x.r_mean + sample/N) * N/N_plus_1; + } + else + { + x.r_mean = sample; + x.min_val = sample; + x.max_val = sample; + + // r_var is initialised to zero + // in the constructor and reset() + } + + x.counter++; + } + + + +//! update statistics to reflect new sample (version for non-complex numbers, complex sample) +template +inline +void +running_stat_aux::update_stats(running_stat& x, const std::complex& sample, const typename arma_not_cx::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + running_stat_aux::update_stats(x, std::real(sample)); + } + + + +//! update statistics to reflect new sample (version for complex numbers, non-complex sample) +template +inline +void +running_stat_aux::update_stats(running_stat& x, const typename eT::value_type sample, const typename arma_cx_only::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + typedef typename eT::value_type T; + + running_stat_aux::update_stats(x, std::complex(sample)); + } + + + +//! alter statistics to reflect new sample (version for complex numbers, complex sample) +template +inline +void +running_stat_aux::update_stats(running_stat& x, const eT& sample, const typename arma_cx_only::result* junk) + { + arma_extra_debug_sigprint(); + arma_ignore(junk); + + typedef typename eT::value_type T; + + const T sample_norm = std::norm(sample); + const T N = x.counter.value(); + + if(N > T(0)) + { + if(sample_norm < x.min_val_norm) + { + x.min_val_norm = sample_norm; + x.min_val = sample; + } + + if(sample_norm > x.max_val_norm) + { + x.max_val_norm = sample_norm; + x.max_val = sample; + } + + const T N_plus_1 = x.counter.value_plus_1(); + const T N_minus_1 = x.counter.value_minus_1(); + + x.r_var = N_minus_1/N * x.r_var + std::norm(sample - x.r_mean)/N_plus_1; + + x.r_mean = x.r_mean + (sample - x.r_mean)/N_plus_1; + //x.r_mean = (N/N_plus_1)*x.r_mean + sample/N_plus_1; + //x.r_mean = (x.r_mean + sample/N) * N/N_plus_1; + } + else + { + x.r_mean = sample; + x.min_val = sample; + x.max_val = sample; + x.min_val_norm = sample_norm; + x.max_val_norm = sample_norm; + + // r_var is initialised to zero + // in the constructor and reset() + } + + x.counter++; + } + + + +//! @} -- cgit v1.2.1