diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index d3122648..57edb9d9 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -1,5 +1,7 @@ ## New features + 1. Add estimators for haromins, geometic. power, Lehmer means and their weighted analogues + ## Backward incompatible ## Bug fixes diff --git a/ostap/stats/statvars.py b/ostap/stats/statvars.py index 82032a68..c938a99e 100644 --- a/ostap/stats/statvars.py +++ b/ostap/stats/statvars.py @@ -9,6 +9,7 @@ """Functions to collect statistics for trees and datasets - data_get_moment - calculate the moment - data_moment - get the moment (with uncertainty) +- data_get_stat - get the momentt-based statistics - data_the_moment - get the (central) moment - data_central_moment - get the central moment (with uncertainty) - data_mean - get the mean (with uncertainty) @@ -25,6 +26,10 @@ - data_quartiles - get three quartiles - data_quintiles - get four quintiles - data_deciles - get nine deciles +- data_harmonic_mean - get the (weighted) harmonic mean +- data_geometric_mean - get the geometric mean +- data_power_mean - get the power mean +- data_lehmer_mean - get Lehmer mean """ # ============================================================================= __version__ = "$Revision$" @@ -34,6 +39,7 @@ 'data_get_moment' , ## calculate the moment 'data_moment' , ## get the moment (with uncertainty) 'data_the_moment' , ## get the cental moment + 'data_get_stat' , ## get the momentt-based statistics 'data_central_moment' , ## get the central moment (with uncertainty) 'data_mean' , ## get the mean (with uncertainty) 'data_variance' , ## get the variance (with uncertainty) @@ -49,6 +55,10 @@ 'data_quartiles' , ## get three quartiles 'data_quintiles' , ## get four quintiles 'data_deciles' , ## get nine deciles + 'data_harmonic_mean' , ## get the (weighted) harmonic mean + 'data_geometric_mean' , ## get the geometric mean + 'data_power_mean' , ## get the power mean + 'data_lehmer_mean' , ## get Lehmer mean 'data_decorate' , ## technical function to decorate the class ) # ============================================================================= @@ -68,6 +78,7 @@ ## @var QEXACT # use it as threshold for exact/slow vs approximate/fast quantile calcualtion QEXACT = 10000 + # ============================================================================= ## get the moment of order 'order' relative to 'center' # @code @@ -112,8 +123,6 @@ def data_moment ( data , order , expression , cuts = '' , *args ) : order , expression , cuts , *args ) - - # ============================================================================= ## get the central moment (with uncertainty) of order 'order' @@ -143,6 +152,177 @@ def data_central_moment ( data , order , expression , cuts = '' , *args ) : expression , cuts , *args ) + +# ============================================================================= +## Get the (w)moments-base statsitocs from data +# @code +# statobj = Ostap.Math.MinValue() +# data = ... +# result = data.get_moment ( statobj , 'x+y' , 'pt>1' ) +# @encode +# @see Ostap::Math::Moment +# @see Ostap::Math::WMoment +# @see Ostap::statVar::the_moment +def data_get_stat ( data , statobj , expression , cuts = '' , *args ) : + """Get the (w)moments -based statistics + >>> data = ... + >>> stat = Ostap.Math.MinValue() + >>> result = data.get_thestat ( stat , 'x/y+z' , '00' ) +# @endcode +# @see Ostap::Math::Moment +# @see Ostap::Math::WMoment +# @see Ostap::Math::HarmonicMean +# @see Ostap::Math::WHarmonicMean +# @see Ostap::statVar::the_moment +def data_harmonic_mean ( data , expression , cuts = '' , *args ) : + """ Get harmonic mean over the data + >>> data = ... + >>> result = data_harmonic_mean ( data , 'pt' , 'eta>0' ) + - see `Ostap.Math::Moment` + - see `Ostap.Math::WMoment` + - see `Ostap.Math::HarmonicMean` + - see `Ostap.Math::WHarmonicMean` + - see `Ostap.statVar.the_moment` + """ + + if ( not cuts ) and isinstance ( data , Tree ) : + statobj = Ostap.Math.HarmonicMean() + return data_get_stat ( data , statobj , expression , '' , *args ) + + statobj = Ostap.Math.WHarmonicMean() + return data_get_stat ( data , statobj , expression , cuts , *args ) + +# ============================================================================= +## Get geometric mean over the data +# @code +# data = ... +# result = data_geometric_mean ( data , 'pt' , 'eta>0' ) +# @endcode +# @see Ostap::Math::Moment +# @see Ostap::Math::GeometricMean +# @see Ostap::Math::WMoment +# @see Ostap::Math::WGeometricMean +# @see Ostap::statVar::the_moment +def data_geometric_mean ( data , expression , cuts = '' , *args ) : + """ Get geometric mean over the data + >>> data = ... + >>> result = data_geometric_mean ( data , 'pt' , 'eta>0' ) + - see `Ostap.Math::Moment` + - see `Ostap.Math::HarmonicMean` + - see `Ostap.statVar.the_moment` + """ + + if ( not cuts ) and isinstance ( data , Tree ) : + statobj = Ostap.Math.GeometricMean() + return data_get_stat ( data , statobj , expression , '' , *args ) + + statobj = Ostap.Math.WGeometricMean() + return data_get_stat ( data , statobj , expression , cuts , *args ) + +# ============================================================================= +## Get power mean over the data +# @code +# data = ... +# result = data_power_mean ( data , 5 , 'pt' , 'eta>0' ) +# @endcode +# @see Ostap::Math::Moment +# @see Ostap::Math::WMoment +# @see Ostap::Math::PowerMean +# @see Ostap::Math::WPowerMean +# @see Ostap::statVar::the_moment +def data_power_mean ( data , p , expression , cuts = '' , *args ) : + """ Get power mean over the data + >>> data = ... + >>> result = data_power_mean ( data , 5 , 'pt' , 'eta>0' ) + - see `Ostap.Math::Moment` + - see `Ostap.Math::WMoment` + - see `Ostap.Math::PowerMean` + - see `Ostap.Math::WPowerMean` + - see `Ostap.statVar.the_moment` + """ + from ostap.core.ostap_types import num_types + from ostap.math.base import isequal, iszero + + assert isinstance ( p , num_types ) , 'Invalid p-parameter!' + + if ( not cuts ) and isinstance ( data , Tree ) : + if p == -1 or isequal ( p , -1. ) : statobj = Ostap.Math.HarmonicMean ( ) + elif p == 0 or iszero ( p ) : statobj = Ostap.Math.GeometricMean ( ) + elif p == 1 or isequa ( p , 1. ) : statobj = Ostap.Math.ArithmeticMean ( 0 ) + else : statobj = Ostap.Math.PowerMean ( p ) + return data_get_stat ( data , statobj , expression , '' , *args ) + + if p == -1 or isequal ( p , -1. ) : statobj = Ostap.Math.WHarmonicMean ( ) + elif p == 0 or iszero ( p ) : statobj = Ostap.Math.WGeometricMean ( ) + elif p == 1 or isequa ( p , 1. ) : statobj = Ostap.Math.WArithmeticMean ( ) + else : statobj = Ostap.Math.WPowerMean ( p ) + + return data_get_stat ( data , statobj , expression , cuts , *args ) + +# ============================================================================= +## Get Lehmer mean over the data +# @code +# data = ... +# result = data_lehmer_mean ( data , 5 , 'pt' , 'eta>0' ) +# @endcode +# @see Ostap::Math::Moment +# @see Ostap::Math::WMoment +# @see Ostap::Math::LehmerMean +# @see Ostap::Math::WLehmerMean +# @see Ostap::statVar::the_moment +def data_lehmer_mean ( data , p , expression , cuts = '' , *args ) : + """ Get Lehmer mean over the data + >>> data = ... + >>> result = data_lehmer_mean ( data , 5 , 'pt' , 'eta>0' ) + - see `Ostap.Math::Moment` + - see `Ostap.Math::WMoment` + - see `Ostap.Math::LehmerMean` + - see `Ostap.Math::WLehmerMean` + - see `Ostap.statVar.the_moment` + """ + from ostap.core.ostap_types import num_types + from ostap.math.base import isequal, iszero + + assert isinstance ( p , num_types ) , 'Invalid p-parameter!' + + if ( not cuts ) and isinstance ( data , Tree ) : + + if p == 0 or iszero ( p ) : statobj = Ostap.Math.HarmonicMean ( ) + elif 2 * p == 1 or isequal ( p , 0.5 ) : statobj = Ostap.Math.GeometricMean ( ) + elif p == 1 or isequal ( p , 1.0 ) : statobj = Ostap.Math.ArithmeticMean ( ) + else : statobj = Ostap.Math.LehmerMean ( p ) + + return data_get_stat ( data , statobj , expression , '' , *args ) + + if p == 0 or iszero ( p ) : statobj = Ostap.Math.WHarmonicMean ( ) + elif 2 * p == 1 or isequal ( p , 0.5 ) : statobj = Ostap.Math.WGeometricMean ( ) + elif p == 1 or isequal ( p , 1.0 ) : statobj = Ostap.Math.WArithmeticMean ( ) + else : statobj = Ostap.Math.WLehmerMean ( p ) + + return data_get_stat ( data , statobj , expression , cuts , *args ) + # ============================================================================= ## Get the moments or order order as Ostap::Math::(W)Moment_ # @code @@ -164,16 +344,17 @@ def data_the_moment ( data , order , expression , cuts = '' , *args ) : M = Ostap.Math. Moment_(order) moment = M () with rootException() : - sc = StatVar.the_moment ( data , moment , *args ) + sc = StatVar.the_moment ( data , moment , expression , *args ) assert sc.isSuccess() , 'Error %s from StatVar::the_moment' % sc return moment M = Ostap.Math.WMoment_(order) moment = M () with rootException() : - sc = StatVar.the_moment ( data , moment , cuts , *args ) + sc = StatVar.the_moment ( data , moment , expression , cuts , *args ) assert sc.isSuccess() , 'Error %s from StatVar::the_moment' % sc return moment + # ============================================================================= ## get the skewness (with uncertainty) @@ -561,6 +742,12 @@ def data_decorate ( klass ) : klass.quintiles = data_quintiles klass.deciles = data_deciles + klass.get_stats = data_get_stat + klass.harmonic_mean = data_harmonic_mean + klass.geometric_mean = data_geometric_mean + klass.power_mean = data_power_mean + klass.lehmer_mean = data_lehmer_mean + return ( klass.get_moment , klass.moment , klass.central_moment , @@ -577,7 +764,12 @@ def data_decorate ( klass ) : klass.terciles , klass.quartiles , klass.quintiles , - klass.deciles ) + klass.deciles , + klass.get_stats , + klass.harmonic_mean , + klass.geometric_mean , + klass.power_mean , + klass.lehmer_mean ) # ============================================================================= diff --git a/source/include/Ostap/Moments.h b/source/include/Ostap/Moments.h index f1ef9c34..b7291f39 100644 --- a/source/include/Ostap/Moments.h +++ b/source/include/Ostap/Moments.h @@ -726,7 +726,7 @@ namespace Ostap inline void swap ( Moment_& a , Moment_& b ) { a.swap( b ) ; } // ======================================================================== - + // ======================================================================== // Weighted moments // ======================================================================== @@ -1328,7 +1328,730 @@ namespace Ostap inline void swap ( WMoment_& a , WMoment_& b ) { a.swap( b ) ; } // ======================================================================== + + // ======================================================================== + // Other type of counters + // ======================================================================== + + // ======================================================================== + /** @class GeometricMean + * Calculate the geometric mean + * \f$ \left(x_1x_2...x_n\right)^{\frac{1}{n}} \f$ + * @see https://en.wikipedia.org/wiki/Geometric_mean + */ + class GeometricMean : public Moment + { + public: + // ====================================================================== + /// get the geometric mean + inline double mean () const { return std::pow ( 2 , m_log2.mean() ) ; } + // ====================================================================== + public: + // ====================================================================== + inline GeometricMean& operator+=( const double x ) { return add ( x ) ; } + inline GeometricMean& operator+=( const GeometricMean& x ) { return add ( x ) ; } + inline GeometricMean& operator*=( const double x ) { return add ( x ) ; } + inline GeometricMean& operator*=( const GeometricMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + inline GeometricMean& add ( const double x ) + { + if ( 0 < x ) { m_log2.add ( std::log2 ( x ) ) ; } + return *this ; + } + inline GeometricMean& add ( const GeometricMean& x ) + { + m_log2.add ( x.m_log2 ) ; + return *this ; + } + // ====================================================================== + /// add sequence of values + template + inline void add ( ITERATOR begin , ITERATOR end ) + { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a Moment + void update ( const double x ) override { add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_log2.size () ; } + /// empty ? + inline bool empty () const { return m_log2.empty () ; } + /// ok ? + inline bool ok () const { return m_log2.ok () ; } + // ====================================================================== + private: + // ====================================================================== + /// get the counter of log2(x) + Moment_<1> m_log2 {} ; + // ====================================================================== + }; + + // ======================================================================== + /** @class HarmonicMean + * Calcualet the harmonic mean + * \f$ \frac{n}{ \frac{1}{x_1} + ... + \frac{1}{x_n}} \f$ + * @see https://en.wikipedia.org/wiki/Harmonic_mean + */ + class HarmonicMean : public Moment + { + public: + // ====================================================================== + /// get the harmonic mean + inline double mean () const { return 1. / m_inv.mean() ; } + // ====================================================================== + public: + // ====================================================================== + inline HarmonicMean& operator+=( const double x ) { return add ( x ) ; } + inline HarmonicMean& operator+=( const HarmonicMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only non-zero entries + HarmonicMean& add ( const double x ) ; + /// add two counters togather + inline HarmonicMean& add ( const HarmonicMean& x ) + { + m_inv.add ( x.m_inv ) ; + return *this ; + } + // ====================================================================== + /// add sequence of values + template + inline void add ( ITERATOR begin , ITERATOR end ) + { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a Moment + void update ( const double x ) override { add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_inv.size () ; } + /// empty ? + inline bool empty () const { return m_inv.empty () ; } + /// ok ? + inline bool ok () const { return m_inv.ok () ; } + // ====================================================================== + private: + // ====================================================================== + /// get the counter of 1/x + Moment_<1> m_inv {} ; + // ====================================================================== + }; + // ======================================================================== + /** @class PowerMean + * Calculate the power mean + * \f$ \left(\frac{1}{n}\sum x_i^p \right)^{\frac{1}{p}}\f$ + * @see https://en.wikipedia.org/wiki/Power_mean + */ + class PowerMean : public Moment + { + public: + // ====================================================================== + // constructir from the power + PowerMean ( const double p = 1 ) ; + // ====================================================================== + public: + // ====================================================================== + /// get the power mean + inline double mean () const { return std::pow ( m_pow.mean() , 1 / m_p ) ; } + // ====================================================================== + public: + // ====================================================================== + inline PowerMean& operator+=( const double x ) { return add ( x ) ; } + inline PowerMean& operator+=( const PowerMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + inline PowerMean& add ( const double x ) + { + if ( 0 < x ) { m_pow.add ( std::pow ( x , m_p ) ) ; } + return *this ; + } + /// add two counters togather if p is common + PowerMean& add ( const PowerMean& x ) ; + // ====================================================================== + /// add sequence of values + template + inline void add ( ITERATOR begin , ITERATOR end ) + { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a Moment + void update ( const double x ) override { add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_pow.size () ; } + /// empty ? + inline bool empty () const { return m_pow.empty () ; } + /// ok ? + inline bool ok () const { return m_pow.ok () ; } + /// power + inline double p () const { return m_p ; } + // ====================================================================== + private: + // ====================================================================== + /// the power + double m_p {1} ; + /// get the counter of x^p + Moment_<1> m_pow {} ; + // ====================================================================== + }; + // ======================================================================== + /** @class LehmerMean + * Calculate the Lehmer mean + * \f$ \frac{ \sum_i x_i^p}{\sum_i x_i^{p-1}} \f$ + * - \f$ p \rigtharrow - \infty\f$ : minimal value + * - \f$ p = 0 \f$ : harmonic mean + * - \f$ p = 1/2 \f$ : geometric mean + * - \f$ p = 1 \f$ : arithmetic mean + * - \f$ p = 2 \f$ : contraharmonic mean + * - \f$ p \rigtharrow + \infty\f$ : maximal value + * @see https://en.wikipedia.org/wiki/Lehmer_mean + */ + class LehmerMean : public Moment + { + public: + // ====================================================================== + // constructir from the power + LehmerMean ( const double p = 1 ) ; + // ====================================================================== + public: + // ====================================================================== + /// get the power mean + inline double mean () const { return m_lp.mean() / m_lpm1.mean () ; } + // ====================================================================== + public: + // ====================================================================== + inline LehmerMean& operator+=( const double x ) { return add ( x ) ; } + inline LehmerMean& operator+=( const LehmerMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + inline LehmerMean& add ( const double x ) + { + if ( 0 < x ) + { + m_lp .add ( std::pow ( x , m_p ) ) ; + m_lpm1.add ( std::pow ( x , m_p - 1 ) ) ; + } + return *this ; + } + /// add two counters togather if p is common + LehmerMean& add ( const LehmerMean& x ) ; + // ====================================================================== + /// add sequence of values + template + inline void add ( ITERATOR begin , ITERATOR end ) + { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a Moment + void update ( const double x ) override { add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_lp.size () ; } + /// empty ? + inline bool empty () const { return m_lp.empty () ; } + /// ok ? + inline bool ok () const { return m_lp.ok () ; } + /// power + inline double p () const { return m_p ; } + // ====================================================================== + private: + // ====================================================================== + /// the power + double m_p {1} ; + /// get the counter of x^p + Moment_<1> m_lp {} ; + /// get the counter of x^(p-1) + Moment_<1> m_lpm1 {} ; + // ====================================================================== + }; + + // ======================================================================== + /** @class WGeometricMean + * Calculate the weighted geometric mean + * @see https://en.wikipedia.org/wiki/Geometric_mean + */ + class WGeometricMean : public WMoment + { + public: + // ====================================================================== + /// get the geometric mean + inline double mean () const { return std::pow ( 2 , m_log2.mean() ) ; } + // ====================================================================== + public: + // ====================================================================== + inline WGeometricMean& operator+=( const WGeometricMean& x ) { return add ( x ) ; } + inline WGeometricMean& operator*=( const WGeometricMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + inline WGeometricMean& add ( const double x , + const double w = 1 ) + { + if ( 0 < x ) { m_log2.add ( std::log2 ( x ) , w ) ; } + return *this ; + } + inline WGeometricMean& add ( const WGeometricMean& x ) + { + m_log2.add ( x.m_log2 ) ; + return *this ; + } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a WMoment + void update + ( const double x , + const double w ) override { add ( x , w ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_log2.size () ; } + /// number of effective entries + inline long double nEff () const { return m_log2.nEff () ; } + /// get sum of weighes \f$ \sum_i w_i \f$ + inline long double w () const { return m_log2.w () ; } + /// get sum of weights squared + inline long double w2 () const { return m_log2.w2 () ; } + /// empty ? + inline bool empty () const { return m_log2.empty () ; } + /// ok ? + inline bool ok () const { return m_log2.ok () ; } + // ====================================================================== + private: + // ====================================================================== + /// get the counter of log2(x) + WMoment_<1> m_log2 {} ; + // ====================================================================== + }; + // ======================================================================== + + // ======================================================================== + /** @class WHarmonicMean + * Calcualet the weighted harmonic mean + * @see https://en.wikipedia.org/wiki/Harmonic_mean + */ + class WHarmonicMean : public WMoment + { + public: + // ====================================================================== + /// get the harmonic mean + inline double mean () const { return 1. / m_inv.mean() ; } + // ====================================================================== + public: + // ====================================================================== + inline WHarmonicMean& operator+=( const WHarmonicMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only non-zero entries + WHarmonicMean& + add ( const double x , + const double w = 1 ) ; + /// add two counters togather + inline WHarmonicMean& add ( const WHarmonicMean& x ) + { + m_inv.add ( x.m_inv ) ; + return *this ; + } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a WMoment + void update + ( const double x , + const double w ) override { add ( x , w ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_inv.size () ; } + /// number of effective entries + inline long double nEff () const { return m_inv.nEff () ; } + /// get sum of weighes \f$ \sum_i w_i \f$ + inline long double w () const { return m_inv.w () ; } + /// get sum of weights squared + inline long double w2 () const { return m_inv.w2 () ; } + /// empty ? + inline bool empty () const { return m_inv.empty () ; } + /// ok ? + inline bool ok () const { return m_inv.ok () ; } + // ====================================================================== + private: + // ====================================================================== + /// get the counter of 1/x + WMoment_<1> m_inv {} ; + // ====================================================================== + }; + // ======================================================================== + /** @class WPowerMean + * Calculate the weighted power mean + * @see https://en.wikipedia.org/wiki/Power_mean + */ + class WPowerMean : public WMoment + { + public: + // ====================================================================== + // constructir from the power + WPowerMean ( const double p = 1 ) ; + // ====================================================================== + public: + // ====================================================================== + /// get the power mean + inline double mean () const { return std::pow ( m_pow.mean() , 1 / m_p ) ; } + // ====================================================================== + public: + // ====================================================================== + inline WPowerMean& operator+=( const WPowerMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + WPowerMean& add + ( const double x , + const double w = 1 ) ; + /// add two counters togather if p is common + WPowerMean& add ( const WPowerMean& x ) ; + // ====================================================================== + public: + // ====================================================================== + /// to use it as a Moment + void update + ( const double x , + const double w ) override { add ( x , w ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_pow.size () ; } + /// number of effective entries + inline long double nEff () const { return m_pow.nEff () ; } + /// get sum of weighes \f$ \sum_i w_i \f$ + inline long double w () const { return m_pow.w () ; } + /// get sum of weights squared + inline long double w2 () const { return m_pow.w2 () ; } + /// empty ? + inline bool empty () const { return m_pow.empty () ; } + /// ok ? + inline bool ok () const { return m_pow.ok () ; } + /// power + inline double p () const { return m_p ; } + // ====================================================================== + private: + // ====================================================================== + /// the power + double m_p {1} ; + /// get the counter of x^p + WMoment_<1> m_pow {} ; + // ====================================================================== + }; + + + // ======================================================================== + /** @class WLehmerMean + * Calculate the weighted Lehmer mean + * @see https://en.wikipedia.org/wiki/Lehmer_mean + */ + class WLehmerMean : public WMoment + { + public: + // ====================================================================== + // constructir from the power + WLehmerMean ( const double p = 1 ) ; + // ====================================================================== + public: + // ====================================================================== + /// get the power mean + inline double mean () const { return m_lp.mean() / m_lpm1.mean () ; } + // ====================================================================== + public: + // ====================================================================== + inline WLehmerMean& operator+=( const WLehmerMean& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries with non-zero weight + WLehmerMean& add + ( const double x , + const double w = 1 ) ; + /// add two counters togather if p is common + WLehmerMean& add ( const WLehmerMean& x ) ; + // ====================================================================== + public: + // ====================================================================== + /// to use it as a WMoment + void update + ( const double x , + const double w ) override { add ( x , w ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_lp.size () ; } + /// number of effective entries + inline long double nEff () const { return m_lp.nEff () ; } + /// get sum of weighes \f$ \sum_i w_i \f$ + inline long double w () const { return m_lp.w () ; } + /// get sum of weights squared + inline long double w2 () const { return m_lp.w2 () ; } + /// empty ? + inline bool empty () const { return m_lp.empty () ; } + /// ok ? + inline bool ok () const { return m_lp.ok () ; } + /// power + inline double p () const { return m_p ; } + // ====================================================================== + private: + // ====================================================================== + /// the power + double m_p {1} ; + /// get the counter of x^p + WMoment_<1> m_lp {} ; + /// get the counter of x^(p-1) + WMoment_<1> m_lpm1 {} ; + // ====================================================================== + }; + // ======================================================================== + /** @class ArithmeticMean + * Calculate the arithmetic mean + */ + class ArithmeticMean : public Moment_<0> {} ; + // ======================================================================== + /** @class WArithmeticMean + * Calculate the weighted arithmetic mean + */ + class WArithmeticMean : public WMoment_<0> {} ; + // ======================================================================== + /** @class MinMean + * Degenerate case of "mean" : minimal value + * \f$ min ( x_1 , .. , x_n ) \f$ + */ + class MinValue : public Moment + { + public: + // ====================================================================== + /// default constructor + MinValue () ; + // ====================================================================== + public: + // ====================================================================== + /// get the min mean + inline double mean () const { return m_min ; } + /// get the min value + inline double min () const { return m_min ; } + // ====================================================================== + public: + // ====================================================================== + inline MinValue& operator+=( const double x ) { return add ( x ) ; } + inline MinValue& operator+=( const MinValue& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + inline MinValue& add ( const double x ) + { + m_min = std::min ( m_min , x ) ; + m_cnt.add ( x ) ; + return *this ; + } + inline MinValue& add ( const MinValue& x ) + { + m_min = std::min ( m_min , x.m_min ) ; + m_cnt.add ( x.m_cnt ) ; + return *this ; + } + // ====================================================================== + /// add sequence of values + template + inline void add ( ITERATOR begin , ITERATOR end ) + { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a WMoment + void update ( const double x ) override { add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_cnt.size () ; } + /// empty ? + inline bool empty () const { return m_cnt.empty () ; } + /// ok ? + inline bool ok () const { return m_cnt.ok () ; } + // ====================================================================== + private: + // ====================================================================== + /// minimal value ; + double m_min ; + /// get the counter + Moment_<0> m_cnt {} ; + // ====================================================================== + } ; + // ======================================================================== + /** @class MaxValue + * Degenerate case of "mean" : maximal value + * \f$ max ( x_1 , .. , x_n ) \f$ + */ + class MaxValue : public Moment + { + public: + // ====================================================================== + /// default constructor + MaxValue () ; + // ====================================================================== + public: + // ====================================================================== + /// get the max value + inline double mean () const { return m_max ; } + /// get the max value + inline double max () const { return m_max ; } + // ====================================================================== + public: + // ====================================================================== + inline MaxValue& operator+=( const double x ) { return add ( x ) ; } + inline MaxValue& operator+=( const MaxValue& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + inline MaxValue& add ( const double x ) + { + m_max = std::max ( m_max , x ) ; + m_cnt.add ( x ) ; + return *this ; + } + inline MaxValue& add ( const MaxValue& x ) + { + m_max = std::max ( m_max , x.m_max ) ; + m_cnt.add ( x.m_cnt ) ; + return *this ; + } + // ====================================================================== + /// add sequence of values + template + inline void add ( ITERATOR begin , ITERATOR end ) + { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a WMoment + void update ( const double x ) override { add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_cnt.size () ; } + /// empty ? + inline bool empty () const { return m_cnt.empty () ; } + /// ok ? + inline bool ok () const { return m_cnt.ok () ; } + // ====================================================================== + private: + // ====================================================================== + /// maximal value ; + double m_max ; + /// get the counter + Moment_<0> m_cnt {} ; + // ====================================================================== + } ; + // ======================================================================== + /** @class MinMaxValue + */ + class MinMaxValue : public Moment + { + public: + // ====================================================================== + /// default constructor + MinMaxValue () ; + // ====================================================================== + public: + // ====================================================================== + /// get the minvalue + inline double min () const { return m_min ; } + /// get the max value + inline double max () const { return m_max ; } + // ====================================================================== + public: + // ====================================================================== + inline MinMaxValue& operator+=( const double x ) { return add ( x ) ; } + inline MinMaxValue& operator+=( const MinMaxValue& x ) { return add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// accumulate only positive entries + inline MinMaxValue& add ( const double x ) + { + m_min = std::min ( m_min , x ) ; + m_max = std::max ( m_max , x ) ; + m_cnt.add ( x ) ; + return *this ; + } + inline MinMaxValue& add ( const MinMaxValue& x ) + { + m_min = std::min ( m_min , x.m_min ) ; + m_max = std::max ( m_max , x.m_max ) ; + m_cnt.add ( x.m_cnt ) ; + return *this ; + } + // ====================================================================== + /// add sequence of values + template + inline void add ( ITERATOR begin , ITERATOR end ) + { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } + // ====================================================================== + public: + // ====================================================================== + /// to use it as a WMoment + void update ( const double x ) override { add ( x ) ; } + // ====================================================================== + public: + // ====================================================================== + /// number of entries + inline unsigned long long size () const { return m_cnt.size () ; } + /// empty ? + inline bool empty () const { return m_cnt.empty () ; } + /// ok ? + inline bool ok () const { return m_cnt.ok () ; } + // ====================================================================== + private: + // ====================================================================== + /// minimal value ; + double m_min ; + /// maximal value ; + double m_max ; + /// get the counter + Moment_<0> m_cnt {} ; + // ====================================================================== + } ; + + + + + // ======================================================================== /** @class Moments diff --git a/source/src/Moments.cpp b/source/src/Moments.cpp index 23b8e798..059f346b 100644 --- a/source/src/Moments.cpp +++ b/source/src/Moments.cpp @@ -9,6 +9,7 @@ // ============================================================================ // Local // ============================================================================ +#include "Exception.h" #include "local_math.h" // ============================================================================ /** @file @@ -33,6 +34,157 @@ Ostap::Math::WMoment::~WMoment(){} // =========================================================================== +// =========================================================================== +// HarmoincMean:: accumulate only non-zero entries +// =========================================================================== +Ostap::Math::HarmonicMean& +Ostap::Math::HarmonicMean::add ( const double x ) +{ + if ( !s_zero ( x ) ) { m_inv.add ( 1/x ) ; } + return *this ; +} +// =========================================================================== +// WHarmoincMean:: accumulate only non-zero entries +// =========================================================================== +Ostap::Math::WHarmonicMean& +Ostap::Math::WHarmonicMean::add +( const double x , + const double w ) +{ + if ( !s_zero ( x ) ) { m_inv.add ( 1/x , w ) ; } + return *this ; +} + + +// =========================================================================== +Ostap::Math::PowerMean::PowerMean +( const double p ) + : m_p ( p ) + , m_pow ( ) +{} +// =========================================================================== +// add two counters togather if p is common +// =========================================================================== +Ostap::Math::PowerMean& +Ostap::Math::PowerMean::add ( const Ostap::Math::PowerMean& x ) +{ + Ostap::Assert ( s_equal ( m_p , x.m_p ) , + "Cannot add counters with non-eual values of 'p'" , + "Ostap::Math::PowerMean" ) ; + m_pow.add ( x.m_pow ) ; + return *this ; +} + +// =========================================================================== +Ostap::Math::WPowerMean::WPowerMean +( const double p ) + : m_p ( p ) + , m_pow ( ) +{} +// =========================================================================== +// add two counters togather if p is common +// =========================================================================== +Ostap::Math::WPowerMean& +Ostap::Math::WPowerMean::add ( const Ostap::Math::WPowerMean& x ) +{ + Ostap::Assert ( s_equal ( m_p , x.m_p ) , + "Cannot add counters with non-eual values of 'p'" , + "Ostap::Math::WPowerMean" ) ; + m_pow.add ( x.m_pow ) ; + return *this ; +} +// =========================================================================== +Ostap::Math::WPowerMean& +Ostap::Math::WPowerMean::add +( const double x , + const double w ) +{ + if ( 0 < x && !s_zero ( w ) ) { m_pow.add ( std::pow ( x , m_p ) , w ) ; } + return *this ; +} + +// =========================================================================== +Ostap::Math::LehmerMean::LehmerMean +( const double p ) + : m_p ( p ) + , m_lp () + , m_lpm1 () +{} +// =========================================================================== +Ostap::Math::WLehmerMean::WLehmerMean +( const double p ) + : m_p ( p ) + , m_lp () + , m_lpm1 () +{} +// =========================================================================== +// add two counters togather if p is common +// =========================================================================== +Ostap::Math::LehmerMean& +Ostap::Math::LehmerMean::add ( const Ostap::Math::LehmerMean& x ) +{ + Ostap::Assert ( s_equal ( m_p , x.m_p ) , + "Cannot add counters with non-eual values of 'p'" , + "Ostap::Math::LehmerMean" ) ; + m_lp .add ( x.m_lp ) ; + m_lpm1 .add ( x.m_lpm1 ) ; + return *this ; +} +// =========================================================================== +Ostap::Math::WLehmerMean& +Ostap::Math::WLehmerMean::add +( const Ostap::Math::WLehmerMean& x ) +{ + Ostap::Assert ( s_equal ( m_p , x.m_p ) , + "Cannot add counters with non-eual values of 'p'" , + "Ostap::Math::WLehmerMean" ) ; + m_lp .add ( x.m_lp ) ; + m_lpm1 .add ( x.m_lpm1 ) ; + return *this ; +} +// =========================================================================== +// accumulate only positive entries with non-zero weight +// =========================================================================== +Ostap::Math::WLehmerMean& +Ostap::Math::WLehmerMean::add +( const double x , + const double w ) +{ + if ( ( 0 < x ) && !s_zero ( x ) && !s_zero ( w ) ) + { + m_lp .add ( std::pow ( x , m_p ) , w ) ; + m_lpm1.add ( std::pow ( x , m_p - 1 ) , w ) ; + } + return *this ; +} + + +// =========================================================================== +// default constructor +// =========================================================================== +Ostap::Math::MinValue::MinValue () + : m_min ( std::numeric_limits::max() ) + , m_cnt () +{} +// =========================================================================== +// default constructor +// =========================================================================== +Ostap::Math::MaxValue::MaxValue () + : m_max ( - std::numeric_limits::max() ) + , m_cnt () +{} +// =========================================================================== +// default constructor +// =========================================================================== +Ostap::Math::MinMaxValue::MinMaxValue () + : m_min ( std::numeric_limits::max() ) + , m_max ( - std::numeric_limits::max() ) + , m_cnt () +{} + +// ====================================================================== + + // int test() // { // Ostap::Math::Moment_<5> m1 ;