diff --git a/ostap/fitting/dataset.py b/ostap/fitting/dataset.py index c4bdd6cb..efa3e224 100644 --- a/ostap/fitting/dataset.py +++ b/ostap/fitting/dataset.py @@ -905,7 +905,6 @@ def target_reset ( t ) : logger.attention ("From ostap v1.10.1.9 variables are treted in natural order (no reverse!)") _to_print -= 1 - print ( 'VARIABLES/1', varlst, cuts ) nvars = len ( varlst ) assert ( 1 == dim and dim <= nvars ) or dim == nvars , \ 'Mismatch between the target/histo dimension %d and input variables %s' % ( dim , varlst ) @@ -934,7 +933,6 @@ def target_reset ( t ) : the_args = ( target , ) + varlst + tail ## very special case of projection several expressions into the same target if 1 == dim and dim < nvars : - print ( 'DS_PROJECT/1' , dim , nvars, varlst , tail ) ## very special case of projection several expressions into the same target htmp = target_copy ( target ) ## prepare temporary object for var in varlst : @@ -947,16 +945,13 @@ def target_reset ( t ) : target += htmp del htmp elif 1 == dim : - print ( 'DS_PROJECT/2' , dim , nvars, varlst , tail ) if progress : sc = Ostap.HistoProject.project ( dataset , progress_conf () , *the_args ) else : sc = Ostap.HistoProject.project ( dataset , *the_args ) if not sc.isSuccess() : logger.error ( "Error from Ostap.HistoProject.project %s" % sc ) elif 2 == dim : - print ( 'DS_PROJECT/3' , dim , nvars, varlst , tail ) if progress : sc = Ostap.HistoProject.project2 ( dataset , progress_conf () , *the_args ) else : sc = Ostap.HistoProject.project2 ( dataset , *the_args ) if not sc.isSuccess() : logger.error ( "Error from Ostap.HistoProject.project2 %s" % sc ) - print ( 'DS_PROJECT/4' , dim , nvars, varlst , tail ) elif 3 == dim : if progress : sc = Ostap.HistoProject.project3 ( dataset , progress_conf () , *the_args ) else : sc = Ostap.HistoProject.project3 ( dataset , *the_args ) diff --git a/ostap/math/minimize.py b/ostap/math/minimize.py index b9feeece..f5a0528f 100644 --- a/ostap/math/minimize.py +++ b/ostap/math/minimize.py @@ -40,7 +40,7 @@ from scipy.optimize import minimize_scalar scipy_OK = True except ImportError : - from ostap.math.local_minimize import scalar_minimuze as minimize_scalar + from ostap.math.local_minimize import scalar_minimize as minimize_scalar scipy_OK = False if scipy_OK : @@ -186,7 +186,6 @@ def sp_maximum_3D ( fun , 'sp_maximum_2D' , 'sp_maximum_3D' ) - # ============================================================================= if '__main__' == __name__ : diff --git a/ostap/math/models.py b/ostap/math/models.py index 6a11dfb8..d784c0aa 100644 --- a/ostap/math/models.py +++ b/ostap/math/models.py @@ -996,9 +996,13 @@ def _tf3_getattr_ ( self , attr ) : raise AttributeError("Can't get attribute: %s" % attr ) - -from ostap.math.minimize import sp_minimum_1D, sp_maximum_1D -from ostap.math.rootfinder import sp_solve +try : + from ostap.math.minimize import sp_minimum_1D, sp_maximum_1D + from ostap.math.rootfinder import sp_solve +except ImportError : + sp_minimum_1D = None + sp_maximum_1D = None + sp_solve = None # ============================================================================= ## decorate 1D-models/functions @@ -1365,7 +1369,12 @@ def _sp_print_ ( self , typ = 'BSpline' ) : Ostap.Math.Spline2D = Ostap.Math.PositiveSpline2D Ostap.Math.Spline2DSym = Ostap.Math.PositiveSpline2DSym -from ostap.math.minimize import sp_minimum_2D, sp_maximum_2D +try : + from ostap.math.minimize import sp_minimum_2D, sp_maximum_2D +except ImportError : + sp_minimum_2D = None + sp_maximum_2D = None + for model in ( Ostap.Math.BSpline2D , Ostap.Math.BSpline2DSym , @@ -1401,9 +1410,12 @@ def _sp_print_ ( self , typ = 'BSpline' ) : if sp_minimum_2D and not hasattr ( model , 'minimum' ) : model.minimum = sp_minimum_2D if sp_maximum_2D and not hasattr ( model , 'maximum' ) : model.maximum = sp_maximum_2D - -from ostap.math.minimize import sp_minimum_3D, sp_maximum_3D - +try : + from ostap.math.minimize import sp_minimum_3D, sp_maximum_3D +except ImportError : + sp_minimum_3D = None + sp_maximum_3D = None + # ============================================================================= ## Decorate 3D models # ============================================================================= diff --git a/ostap/math/reduce.py b/ostap/math/reduce.py index 2e21df8b..8835950b 100644 --- a/ostap/math/reduce.py +++ b/ostap/math/reduce.py @@ -27,11 +27,11 @@ def root_factory ( klass , *params ) : """ return klass ( *params ) + # ============================================================================= ## Simple (basic) polynomials # ============================================================================= - # ============================================================================= ## factory for deserialization of simple polynomians # @see Ostap.Math.Chebyshev diff --git a/ostap/stats/counters.py b/ostap/stats/counters.py index c441b73c..df3e8ce5 100755 --- a/ostap/stats/counters.py +++ b/ostap/stats/counters.py @@ -23,8 +23,8 @@ ) # ============================================================================= from ostap.math.ve import Ostap, VE -from ostap.math.base import isequal, isequalf -from ostap.core.ostap_types import dictlike_types, sequence_types +from ostap.math.base import isequal, isequalf, iszero +from ostap.core.ostap_types import dictlike_types, sequence_types import ROOT, cppyy # ============================================================================= _new_methods_ = [] @@ -51,7 +51,6 @@ _orig_mean = SE.mean SE._orig_mean = _orig_mean - SE. sum = lambda s : VE ( s._orig_sum () , s.sum2 () ) SE. minmax = lambda s : ( s.min () , s.max () ) SE. mean = lambda s : VE ( s._orig_mean () , s.meanErr()**2 ) @@ -289,7 +288,56 @@ def counters_table ( counters , prefix = '' , title = '' ) : counters_table ] - + +# ============================================================================== +## REDUCE +# ============================================================================== +from ostap.math.reduce import root_factory + + +# ============================================================================= +## Reduce Ostap::StatEntity +def _se_reduce_ ( cnt ) : + """Reduce `Ostap.StatEntity` object + """ + return root_factory , ( type ( cnt ) , + cnt.n () , + cnt.mu () , + cnt.mu2 () , + cnt.min () , + cnt.max () ) + +# ============================================================================= +## Reduce Ostap::WStatEntity +def _wse_reduce_ ( cnt ) : + """Reduce `Ostap.WStatEntity` object + """ + return root_factory , ( type ( cnt ) , + cnt.mu () , + cnt.mu2 () , + cnt.values () , + cnt.weights () ) +# ============================================================================= +## Reduce Ostap::NStatEntity +def _nse_reduce_ ( cnt ) : + """Reduce `Ostap.NStatEntity` object + """ + return root_factory , ( type ( cnt ) , + cnt.N () , + cnt.cnt1 () , + cnt.cnt2 () ) + +SE .__reduce__ = _se_reduce_ +WSE.__reduce__ = _wse_reduce_ +NSE.__reduce__ = _nse_reduce_ + +_new_methods_ += [ + SE .__reduce__ , + WSE.__reduce__ , + NSE.__reduce__ , +] + +# ============================================================================== _new_methods_ = tuple ( _new_methods_ ) # ============================================================================= diff --git a/ostap/stats/moment.py b/ostap/stats/moment.py index 17e4028f..e41b9ba2 100644 --- a/ostap/stats/moment.py +++ b/ostap/stats/moment.py @@ -33,7 +33,7 @@ __all__ = () # ============================================================================= from ostap.core.ostap_types import integer_types, num_types -from ostap.math.base import isfinite +from ostap.math.base import isfinite, isequal from ostap.core.core import Ostap, VE from ostap.core.meta_info import root_version_int, root_info import ROOT @@ -648,8 +648,8 @@ def _om_table ( obj , title = '' , prefix = '' , standard = False ) : M0 = Ostap.Math.Moment_(0) M1 = Ostap.Math.Moment_(1) -WM0 = Ostap.Math.Moment_(0) -WM1 = Ostap.Math.Moment_(1) +WM0 = Ostap.Math.WMoment_(0) +WM1 = Ostap.Math.WMoment_(1) for m in ( M0 , M1 , WM0 , WM1 ) : m.mean = _om_mean m.table = _om_table @@ -659,13 +659,236 @@ def _om_table ( obj , title = '' , prefix = '' , standard = False ) : if not hasattr ( WM0 , 'order' ) : WM0.order = 0 if not hasattr ( WM1 , 'order' ) : WM1.order = 1 + + +# ============================================================ +## equality for two counters +def _mom0_eq_ ( cnt , another ) : + """Equality for two Ostap.Math.Momemnt_<0> counters + """ + ## The same base type + if not isinstance ( another , Ostap.Math.Moment ) : return NotImplemented + ## order is defined + if not hasattr ( another , 'order' ) : return NotImplemented + ## the same order + if cnt.order != another.order : return NotImplemented + ### the same size + return cnt.size() == another.size () +# ============================================================ +## equality for two counters +def _mom1_eq_ ( cnt , another ) : + """Equality for two Ostap.Math.Momemnt_<1> counters + """ + ## The same base type + if not isinstance ( another , Ostap.Math.Moment ) : return NotImplemented + ## order is defined + if not hasattr ( another , 'order' ) : return NotImplemented + ## the same order + if cnt.order != another.order : return NotImplemented + return isequal ( cnt.mu() , another.mu() ) and cnt.previous() == another.previous() +# ============================================================ +## equality for two counters +def _momN_eq_ ( cnt , another ) : + """Equality for two Ostap.Math.Moment_ counters + """ + ## The same base type + if not isinstance ( another , Ostap.Math.Moment ) : return NotImplemented + ## order is defined + if not hasattr ( another , 'order' ) : return NotImplemented + ## the same order + if cnt.order != another.order : return NotImplemented + ## + return isequal ( cnt.M () , another.M () ) and cnt.previous() == another.previous() + +Ostap.Math.Moment . __eq__ = _momN_eq_ +M0 . __eq__ = _mom0_eq_ +M1 . __eq__ = _mom1_eq_ + +# ============================================================ +## equality for two counters +def _wmom0_eq_ ( cnt , another ) : + """Equality for two Ostap.Math.Momemnt_<0> counters + """ + ## The same base type + if not isinstance ( another , Ostap.Math.WMoment ) : return NotImplemented + ## order is defined + if not hasattr ( another , 'order' ) : return NotImplemented + ## the same order + if cnt.order != another.order : return NotImplemented + ### the same size + return \ + cnt.size() == another.size () and \ + isequal ( cnt.w () , another.w () ) and \ + isequal ( cnt.w2 () , another.w2 () ) + + +# ============================================================ +## equality for two counters +def _wmom1_eq_ ( cnt , another ) : + """Equality for two Ostap.Math.WMomemnt_<1> counters + """ + ## The same base type + if not isinstance ( another , Ostap.Math.WMoment ) : return NotImplemented + ## order is defined + if not hasattr ( another , 'order' ) : return NotImplemented + ## the same order + if cnt.order != another.order : return NotImplemented + return isequal ( cnt.mu() , another.mu() ) and cnt.previous() == another.previous() +# ============================================================ +## equality for two counters +def _wmomN_eq_ ( cnt , another ) : + """Equality for two Ostap.Math.WMoment_ counters + """ + ## The same base type + if not isinstance ( another , Ostap.Math.WMoment ) : return NotImplemented + ## order is defined + if not hasattr ( another , 'order' ) : return NotImplemented + ## the same order + if cnt.order != another.order : return NotImplemented + ## + return isequal ( cnt.M () , another.M () ) and cnt.previous() == another.previous() + + +Ostap.Math.WMoment . __eq__ = _wmomN_eq_ +WM0 . __eq__ = _wmom0_eq_ +WM1 . __eq__ = _wmom1_eq_ + + +# ========================================================== +## REDUCE +from ostap.math.reduce import root_factory +# ========================================================== +## Redude Ostap::Math::Moment_<0> +def _mom0_reduce_ ( cnt ) : + """Redude Ostap::Math::Moment_<0> + """ + return root_factory , ( type ( cnt ), cnt.size () ) +# ========================================================= +## Redude Ostap::Math::Moment_<1> +def _mom1_reduce_ ( cnt ) : + """Redude Ostap::Math::Moment_<1> + """ + return root_factory , ( type ( cnt ), cnt.mu () , cnt.previous () ) +# ======================================================== +## Redude Ostap::Math::Moment_ +def _momN_reduce_ ( cnt ) : + """Redude Ostap::Math::Moment_ + """ + return root_factory , ( type ( cnt ), cnt.moment () , cnt.previous () ) + +Ostap.Math.Moment . __reduce__ = _momN_reduce_ +M0 . __reduce__ = _mom0_reduce_ +M1 . __reduce__ = _mom1_reduce_ + +# ========================================================== +## Redude Ostap::Math::WMoment_<0> +def _wmom0_reduce_ ( cnt ) : + """Redude Ostap::Math::WMoment_<0> + """ + return root_factory , ( type ( cnt ), cnt.size () , cnt.w() , cnt.w2 () ) +# ========================================================= +## Redude Ostap::Math::WMoment_<1> +def _wmom1_reduce_ ( cnt ) : + """Redude Ostap::Math::WMoment_<1> + """ + return root_factory , ( type ( cnt ), cnt.mu () , cnt.previous () ) +# ======================================================== +## Redude Ostap::Math::WMoment_ +def _wmomN_reduce_ ( cnt ) : + """Redude Ostap::Math::WMoment_ + """ + return root_factory , ( type ( cnt ), cnt.moment () , cnt.previous () ) + +Ostap.Math.WMoment . __reduce__ = _wmomN_reduce_ +WM0 . __reduce__ = _wmom0_reduce_ +WM1 . __reduce__ = _wmom1_reduce_ + + + +# ============================================================================= +## serialization of Harmonic&Geometric means +# @see Ostap::Math::GeometricMean +# @see Ostap::Math::HarmonicMean +# @see Ostap::Math::WGeometricMean +# @see Ostap::Math::WHarmonicMean +def _mn_reduce_ ( cnt ) : + """Serialization of Harmonic&Geometric means + - see Ostap::Math::GeometricMean + - see Ostap::Math::HarmonicMean + - see Ostap::Math::WGeometricMean + - see Ostap::Math::WHarmonicMean + """ + return root_factory , ( type ( cnt ) , cnt.counter() ) + +Ostap.Math.GeometricMean . __reduce__ = _mn_reduce_ +Ostap.Math.HarmonicMean . __reduce__ = _mn_reduce_ +Ostap.Math.ArithmeticMean . __reduce__ = _mn_reduce_ +Ostap.Math.WGeometricMean . __reduce__ = _mn_reduce_ +Ostap.Math.WHarmonicMean . __reduce__ = _mn_reduce_ +Ostap.Math.WArithmeticMean . __reduce__ = _mn_reduce_ + +# ============================================================================= +## serialization of Power means +# @see Ostap::Math::PoweMean +# @see Ostap::Math::WPoweMean +def _pm_reduce_ ( cnt ) : + """Serialization of Power means + - see Ostap::Math::PowerMean + - see Ostap::Math::WPowerMean + """ + return root_factory , ( type ( cnt ) , cnt.p() , cnt.counter() ) + +Ostap.Math.PowerMean . __reduce__ = _pm_reduce_ +Ostap.Math.WPowerMean . __reduce__ = _pm_reduce_ + +# ============================================================================= +## serialization of Lehmermeans +# @see Ostap::Math::LehmerMean +# @see Ostap::Math::WLehmerMean +def _lm_reduce_ ( cnt ) : + """Serialization of Lehmermeans + - see Ostap::Math::LehmerMean + - see Ostap::Math::WLehmerMean + """ + return root_factory , ( type ( cnt ) , cnt.p() , cnt.counter1() , cnt.conuter2() ) + +Ostap.Math.LehmerMean . __reduce__ = _lm_reduce_ +Ostap.Math.WLehmerMean . __reduce__ = _lm_reduce_ + +# ============================================================================== +## equality for the mean +# @see Ostap::Math::GeometricMean +# @see Ostap::Math::HarmonicMean +# @see Ostap::Math::WGeometricMean +# @see Ostap::Math::WHarmonicMean +def _mn_eq_ ( cnt , another ) : + """Equality Harmonic&Geometric means + - see Ostap::Math::GeometricMean + - see Ostap::Math::HarmonicMean + - see Ostap::Math::WGeometricMean + - see Ostap::Math::WHarmonicMean + """ + if not type ( cnt ) is type ( another ) : return NotImplemented + return isequal ( cnt.value() , another.value() ) and \ + cnt.counter() == another.counter() + +Ostap.Math.GeometricMean . __eq__ = _mn_eq_ +Ostap.Math.HarmonicMean . __eq__ = _mn_eq_ +Ostap.Math.WGeometricMean . __eq__ = _mn_eq_ +Ostap.Math.WHarmonicMean . __eq__ = _mn_eq_ + _decorated_classes = ( - Ostap.Math.Moment , - Ostap.Math.WMoment , - Ostap.Math.Moment_(0) , - Ostap.Math.Moment_(1) , - Ostap.Math.Moment_(0) , - Ostap.Math.Moment_(1) , + Ostap.Math.Moment , + Ostap.Math.WMoment , + Ostap.Math.Moment_(0) , + Ostap.Math.Moment_(1) , + Ostap.Math.Moment_(0) , + Ostap.Math.Moment_(1) , + ## + Ostap.Math.GeometricMean , + Ostap.Math.HarmonicMean , + Ostap.Math.WGeometricMean , + Ostap.Math.WHarmonicMean , ) _new_methods_ = ( diff --git a/ostap/stats/moments.py b/ostap/stats/moments.py index b5dfc115..053aaf79 100644 --- a/ostap/stats/moments.py +++ b/ostap/stats/moments.py @@ -110,7 +110,7 @@ ) # ============================================================================= from ostap.core.ostap_types import integer_types, num_types -from ostap.core.core import Ostap +from ostap.core.core import Ostap # ============================================================================= # logging # ============================================================================= diff --git a/ostap/stats/statvars.py b/ostap/stats/statvars.py index e97b1429..8d7461b1 100644 --- a/ostap/stats/statvars.py +++ b/ostap/stats/statvars.py @@ -263,13 +263,9 @@ def data_statistics ( data , expressions , cuts = '' , *args ) : - see Ostap.StatVar.statVar """ - print ('DATA-STATISTICS1' , expressions , cuts , args ) - ## decode expressions & cuts var_lst, cuts, input_string = vars_and_cuts ( expressions , cuts ) - print ('DATA-STATISTICS2' , var_lst , cuts , input_string ) - ## only one name is specified as strnig if input_string : with rootException() : @@ -314,7 +310,6 @@ def data_minmax ( data , expressions , cuts = '' , *args ) : - see Ostap.Math.WStatEntity - see Ostap.StatVar.statVar """ - print ('DATA-MINMAX' , expressions , cuts , args ) results = data_statistics ( data , expressions , cuts , *args ) if isinstance ( results , dictlike_types ) : for k, r in loop_items ( results ) : @@ -345,7 +340,6 @@ def data_range ( data , >>> result = data_range ( data , 'sin(x)*100*y' , 'x<0' ) >>> results = data_range ( dataset , 'x,y,z,t,u,v' , 'x<0' ) ## as dictionary """ - print ('DATA-RANGE' , expressions , cuts , delta , args ) results = data_minmax ( data, expressions , cuts , *args ) if isinstance ( results , dictlike_types ) : for k , r in loop_items ( results ) : diff --git a/ostap/tools/tests/test_tools_reweight2.py b/ostap/tools/tests/test_tools_reweight2.py index 9c26c340..bda825b6 100755 --- a/ostap/tools/tests/test_tools_reweight2.py +++ b/ostap/tools/tests/test_tools_reweight2.py @@ -329,7 +329,6 @@ def prepare_data ( ) : tag = 'Reweighting iteration #%d' % iter logger.info ( allright ( tag ) ) - with timing ( tag + ': prepare MC-dataset:' , logger = logger ) : # ========================================================================= ## 0) The weighter object diff --git a/ostap/trees/trees.py b/ostap/trees/trees.py index 774341df..3db3d783 100755 --- a/ostap/trees/trees.py +++ b/ostap/trees/trees.py @@ -420,7 +420,6 @@ def tree_project ( tree , ## use frame if requested and if/when possible if use_frame and 0 == first and len ( tree ) < last : if ( input_histo and ( 6 , 19 ) <= root_info ) or ( 6,25 ) <= root_info : - print ( 'go to frame' ) import ostap.frames.frames as F frame = F.DataFrame ( tree ) if progress : frame , _ = F.frame_progress ( frame , len ( tree ) ) @@ -2804,7 +2803,6 @@ def get_slices ( self , first , last , chunk_size ) : ## yield lst [ i : min ( i + chunk_size , list_size ) ] yield slice ( first + i , first + min ( i + chunk_size , last -first ) ) - ## split the chain for several chains with at most chunk_size entries def split ( self , chunk_size = -1 , max_files = 10 ) : """Split the tree for several trees with chunk_size entries @@ -3040,7 +3038,6 @@ def __str__ ( self ) : return r + ")" __repr__ = __str__ - # ========================================================================= ## delegate all other attributes to the underlying chain object def __getattr__ ( self , attr ) : @@ -3052,10 +3049,8 @@ def tree ( self ) : """``tree'' : get the underlying tree/chain""" return self.chain - - # ============================================================================== -## get a full path of the TTree obnject +## get a full path of the TTree object # @code # rdir = ... # path = rdirt.full_path diff --git a/source/include/Ostap/Moments.h b/source/include/Ostap/Moments.h index 0f206426..6893a6eb 100644 --- a/source/include/Ostap/Moments.h +++ b/source/include/Ostap/Moments.h @@ -47,7 +47,7 @@ namespace Ostap public : // ====================================================================== virtual ~Moment() ; - // ====================================================================== + // ====================================================================== } ; // ======================================================================== /// forward declaration @@ -76,13 +76,25 @@ namespace Ostap // ====================================================================== enum _ { order = N } ; // ====================================================================== + public: + // ====================================================================== + /// default constructor + Moment_ () = default ; + /// constructor from the value and previous moment + Moment_ + ( const double mom , + const Moment_& prev ) + : m_prev ( prev ) + , m_M ( mom * prev.size() ) + {} + // ====================================================================== public: // ====================================================================== /** get the value of the Nth moment * \f[ \mu_N \equiv \frac{1}{N} \sum \left( x_i - \bar{x} \right)^N \f] * @return the value of the Nth central moment */ - inline double moment () const { return this->M ( N ) / this->size () ; } + inline double moment () const { return this->ok () ? this->M ( N ) / this->size () : 0 ; } // ====================================================================== /** get value of the kth moment for \f$ k \le N \f$ * \f[ \mu_k \equiv \frac{1}{N} \sum \left( x_i - \bar{x} \right)^k \f] @@ -91,16 +103,21 @@ namespace Ostap * @return the value of the kth central moment if \f$ 0 \le k \le N \f$, 0, otherwise */ inline double moment ( const unsigned short k ) const - { return N < k ? 0 : 0 == k ? 1 : 1 == k ? 0 : ( this->M ( k ) / this->size() ) ; } + { return + N < k ? 0 : + 0 == k ? 1 : + 1 == k ? 0 : + !this->ok () ? 0 : ( this->M ( k ) / this->size() ) ; } // ====================================================================== /// get value of the kth standartized moment for \f$ k \le N \f$ inline double std_moment ( const unsigned short k ) const { return - N < k ? 0 : - 0 == k ? 1 : - 1 == k ? 0 : - 2 == k ? 1 : this->moment ( k ) / std::pow ( this->moment ( 2 ) , 0.5 * k ) ; + N < k ? 0 : + 0 == k ? 1 : + 1 == k ? 0 : + 2 == k ? 1 : + !this->ok ? 0 : this->moment ( k ) / std::pow ( this->moment ( 2 ) , 0.5 * k ) ; } // ====================================================================== /// get number of entries @@ -159,7 +176,7 @@ namespace Ostap * @param k the order \f$ 0 \le k \le N \f$ * @return the value of \$ M_k \f$ if \f$ 0 \le k \le N \f$, 0 otherwise */ - inline long double M ( const unsigned short k ) const + inline long double M ( const unsigned short k = N ) const { return k > N ? 0.0L : k == N ? this->m_M : this->m_prev.M( k ) ; } // ====================================================================== public: @@ -197,7 +214,7 @@ namespace Ostap // ====================================================================== template ::type = 0 > inline long double moment_ () const - { return this->empty () ? 0 : this->m_M / this->size () ; } + { return !this->ok () ? 0 : this->m_M / this->size () ; } // ====================================================================== template K),int>::type = 0 > inline long double moment_ () const @@ -242,7 +259,8 @@ namespace Ostap template =K),int>::type = 0 > inline long double std_moment_ () const { - return this->empty () ? 0 : + return + !this->ok () ? 0 : this->template moment_ () / std::pow ( this->moment_<2>() , 0.5 * K ) ; } // ====================================================================== @@ -457,6 +475,14 @@ namespace Ostap // ====================================================================== enum _ { order = 0 } ; // ====================================================================== + public: + // ====================================================================== + /// (default) constructor + Moment_ + ( const unsigned long long size = 0 ) + : m_size ( size ) + {} + // ====================================================================== public: // ====================================================================== /** get the value of the Nth moment @@ -484,7 +510,7 @@ namespace Ostap /// empty ? inline bool empty () const { return 0 == m_size ; } /// ok ? - inline bool ok () const { return 0 != m_size ; } + inline bool ok () const { return m_size ; } // ====================================================================== public: // basic operations // ====================================================================== @@ -526,7 +552,7 @@ namespace Ostap * @param k the order \f$ 0 \le k \le N \f$ * @return the value of \$ M_k \f$ if \f$ 0 \le k \le \f$, 0 otherwise */ - inline long double M ( const unsigned short k ) const + inline long double M ( const unsigned short k = 0 ) const { return 0 < k ? 0 : m_size ; } // ====================================================================== public: @@ -552,7 +578,7 @@ namespace Ostap // ====================================================================== private: // ====================================================================== - unsigned long long m_size ; + unsigned long long m_size { 0 } ; // ====================================================================== private: // ====================================================================== @@ -570,6 +596,23 @@ namespace Ostap // ====================================================================== enum _ { order = 1 } ; // ====================================================================== + public: + // ====================================================================== + /// default constructor from mu and size + Moment_ + ( const double mu = 0 , + const unsigned long long size = 0 ) + : m_prev ( size ) + , m_mu ( mu ) + {} + /// constructor from mu and previous + Moment_ + ( const double mu , + const Moment_<0>& prev ) + : m_prev ( prev ) + , m_mu ( mu ) + {} + // ====================================================================== public: // ====================================================================== /** get the value of the Nth moment @@ -590,11 +633,11 @@ namespace Ostap { return 0 == k || 2 == k ? 1 : 0 ; } // ====================================================================== /// get number of entries - inline unsigned long long size () const { return m_prev.size () ; } + inline unsigned long long size () const { return m_prev.size () ; } /// get effective number of entries - inline unsigned long long nEff () const { return m_prev.nEff () ; } + inline unsigned long long nEff () const { return m_prev.nEff () ; } // get the mean value - inline long double mu () const { return m_mu ; } + inline long double mu () const { return m_mu ; } /// empty ? inline bool empty () const { return m_prev.empty () ; } /// ok ? @@ -658,7 +701,7 @@ namespace Ostap // ====================================================================== public: // ====================================================================== - inline long double M ( const unsigned short k ) const + inline long double M ( const unsigned short k = 1 ) const { return 1 < k ? 0 : 1 == k ? 0 : this->m_prev. M ( k ) ; } // ====================================================================== public: @@ -801,13 +844,25 @@ namespace Ostap // ====================================================================== enum _ { order = N } ; // ====================================================================== + public: + // ====================================================================== + /// default constructor + WMoment_ () = default ; + /// constructor from the value and previous moment + WMoment_ + ( const double mom , + const WMoment_& prev ) + : m_prev ( prev ) + , m_M ( mom * prev.w () ) + {} + // ====================================================================== public: // ====================================================================== /** get the value of the Nth weighted moment * \f[ \mu_N \equiv \frac{1}{\sum w_i} \sum w_i \left( x_i - \bar{x} \right)^N \f] * @return the value of the Nth central moment */ - inline double moment () const { return this->M ( N ) / this-> w () ; } + inline double moment () const { return this-> ok () ? this->M ( N ) / this-> w () : 0.0 ; } // ====================================================================== /** get value of the kth moment for \f$ k \le N \f$ * \f[ \mu_k \equiv \frac{1}{\sum w_i} \sum w_i \left( x_i - \bar{x} \right)^k \f] @@ -816,7 +871,11 @@ namespace Ostap * @return the value of the kth central moment if \f$ 0 \le k \le N \f$, 0, otherwise */ inline double moment ( const unsigned short k ) const - { return N < k ? 0 : 0 == k ? 1 : 1 == k ? 0 : ( this->M ( k ) / this-> w () ) ; } + { return + N < k ? 0 : + 0 == k ? 1 : + 1 == k ? 0 : + !this->ok () ? 0 : ( this->M ( k ) / this-> w () ) ; } // ====================================================================== /** get value of the kth standartized moment for \f$ k \le N \f$ * \f[ \mu_k \equiv \frac{1}{N} \sum \left( x_i - \bar{x} \right)^k \f] @@ -826,10 +885,11 @@ namespace Ostap */ inline double std_moment ( const unsigned short k ) const { return - N < k ? 0 : - 0 == k ? 1 : - 1 == k ? 0 : - 2 == k ? 1 : this->moment ( k ) / std::pow ( this->moment ( 2 ) , 0.5 * k ) ; + N < k ? 0 : + 0 == k ? 1 : + 1 == k ? 0 : + 2 == k ? 1 : + !this->ok () ? 0 : this->moment ( k ) / std::pow ( this->moment ( 2 ) , 0.5 * k ) ; } // ====================================================================== /// get number of entries @@ -865,7 +925,7 @@ namespace Ostap WMoment_ __radd__ ( const WMoment_& x ) const { return __add__ ( x ) ; } WMoment_& __iadd__ ( const WMoment_& x ) { return add ( x ) ; } // ====================================================================== - public: + public: // ====================================================================== /// update the counter void update ( const double x , const double w = 1 ) override { add ( x , w ) ; } @@ -882,7 +942,7 @@ namespace Ostap * @param k the order \f$ 0 \le k \le N \f$ * @return the value of \$ M_k \f$ if \f$ 0 \le k \le N \f$, 0 otherwise */ - inline long double M ( const unsigned short k ) const + inline long double M ( const unsigned short k = N ) const { return k > N ? 0.0L : k == N ? this->m_M : this->m_prev.M( k ) ; } // ====================================================================== public: @@ -1063,6 +1123,20 @@ namespace Ostap enum _ { order = 0 } ; // ====================================================================== public: + // ====================================================================== + /// default constructor + WMoment_ () = default ; + /** full constructor + * @param size number of entries + * @param sumw sum of weights + * @param sumw sum of squared weights + */ + WMoment_ + ( const unsigned long long size , + const double sumw , + const double sumw2 ) ; + // ====================================================================== + public : // ====================================================================== /** get the value of the 0th moment * \f[ \mu_N \equiv \frac{1}{\sum w_i } \sum w_i \left( x_i - \bar{x} \right)^0 \f] @@ -1093,7 +1167,7 @@ namespace Ostap /// empty ? inline bool empty () const { return 0 == m_size ; } /// ok ? - inline bool ok () const { return 0 != m_size && 0 != m_w && 0 < m_w2 ; } + inline bool ok () const { return m_size && m_w && ( 0 < m_w2 ) ; } // ====================================================================== public: // basic operations // ====================================================================== @@ -1143,7 +1217,7 @@ namespace Ostap * @param k the order \f$ 0 \le k \le N \f$ * @return the value of \$ M_k \f$ if \f$ 0 \le k \le \f$, 0 otherwise */ - inline long double M ( const unsigned short k ) const + inline long double M ( const unsigned short k = 0 ) const { return 0 < k ? 0 : m_w ; } // ====================================================================== public: @@ -1174,11 +1248,11 @@ namespace Ostap private: // ====================================================================== /// number of entries - unsigned long long m_size ; // number of entries + unsigned long long m_size { 0 } ; // number of entries /// sum of weights \f$ \sum w_i \f$ - long double m_w ; // sum of weights + long double m_w { 0 } ; // sum of weights /// sum of weights squared \f$ \sum w_i^2 \f$ - long double m_w2 ; // sum of weights squared + long double m_w2 { 0 } ; // sum of weights squared // ====================================================================== private: // ====================================================================== @@ -1196,6 +1270,18 @@ namespace Ostap // ====================================================================== enum _ { order = 1 } ; // ====================================================================== + public: + // ====================================================================== + /// default constructor + WMoment_ () = default ; + /// constructor from mu and previous moment + WMoment_ + ( const double mu , + const WMoment_<0>& prev ) + : m_prev ( prev ) + , m_mu ( mu ) + {} + // ====================================================================== public : // ====================================================================== /** get the value of the 1st moment @@ -1268,14 +1354,14 @@ namespace Ostap // return *this ; } - // ====================================================================== - public: // python operations + // ====================================================================== + public: // python operations // ====================================================================== WMoment_ __add__ ( const WMoment_& x ) const { WMoment_ r (*this) ; r += x ; return r ; } WMoment_ __radd__ ( const WMoment_& x ) const { return __add__ ( x ) ; } WMoment_& __iadd__ ( const WMoment_& x ) { return add ( x ) ; } // ====================================================================== - public: + public: // ====================================================================== /// update the counter void update ( const double x , const double w = 1 ) override { add ( x , w ) ; } @@ -1287,7 +1373,7 @@ namespace Ostap // ====================================================================== public: // ====================================================================== - inline long double M ( const unsigned short k ) const + inline long double M ( const unsigned short k = 1 ) const { return 1 < k ? 0 : 1 == k ? 0 : this->m_prev. M ( k ) ; } // ====================================================================== public: @@ -1335,9 +1421,9 @@ namespace Ostap // ======================================================================= private: // ====================================================================== - WMoment_<0> m_prev { } ; + WMoment_<0> m_prev { } ; /// mean value - long double m_mu {0} ; // mean value + long double m_mu { 0 } ; // mean value // ====================================================================== private: // ====================================================================== @@ -1370,6 +1456,18 @@ namespace Ostap */ class GeometricMean : public Statistic { + public: + // ====================================================================== + /// the actual underlying counter + typedef Moment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// default constructor + GeometricMean () = default ; + /// constructor for the counters + GeometricMean ( const Counter& cnt ) ; + // ====================================================================== public: // ====================================================================== /// get the geometric mean @@ -1383,6 +1481,11 @@ namespace Ostap inline GeometricMean& operator*=( const double x ) { return add ( x ) ; } inline GeometricMean& operator*=( const GeometricMean& x ) { return add ( x ) ; } // ====================================================================== + public: + // ====================================================================== + /// the counter of log2 values + const Counter& counter () const { return m_log ; } + // ====================================================================== public: // ====================================================================== /// accumulate only positive entries @@ -1419,7 +1522,7 @@ namespace Ostap private: // ====================================================================== /// get the counter of log2(x) - Moment_<1> m_log {} ; + Counter m_log {} ; // ====================================================================== }; // ======================================================================== @@ -1430,6 +1533,18 @@ namespace Ostap */ class HarmonicMean : public Statistic { + public: + // ====================================================================== + /// the actual underlying counter + typedef Moment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// default construictor + HarmonicMean () = default ; + /// constructor for the counters + HarmonicMean ( const Counter& cnt ) ; + // ====================================================================== public: // ====================================================================== /// get the harmonic mean @@ -1457,6 +1572,11 @@ namespace Ostap inline void add ( ITERATOR begin , ITERATOR end ) { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } // ====================================================================== + public: + // ====================================================================== + /// the counter of 1/x values + const Counter counter () const { return m_inv ; } + // ====================================================================== public: // ====================================================================== /// to use it as a Moment @@ -1474,7 +1594,7 @@ namespace Ostap private: // ====================================================================== /// get the counter of 1/x - Moment_<1> m_inv {} ; + Counter m_inv {} ; // ====================================================================== }; // ======================================================================== @@ -1487,8 +1607,15 @@ namespace Ostap { public: // ====================================================================== - // constructir from the power - PowerMean ( const double p = 1 ) ; + /// the actual underlying counter + typedef Moment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// defautl construictor + PowerMean ( const double = 1 ) ; + /// constructor for the counters + PowerMean ( const double p , const Counter& cnt ) ; // ====================================================================== public: // ====================================================================== @@ -1501,6 +1628,11 @@ namespace Ostap inline PowerMean& operator+=( const double x ) { return add ( x ) ; } inline PowerMean& operator+=( const PowerMean& x ) { return add ( x ) ; } // ====================================================================== + public: + // ====================================================================== + /// the counter of x**p values + const Counter counter () const { return m_pow ; } + // ====================================================================== public: // ====================================================================== /// accumulate only positive entries @@ -1536,9 +1668,9 @@ namespace Ostap private: // ====================================================================== /// the power - double m_p {1} ; + double m_p {1} ; /// get the counter of x^p - Moment_<1> m_pow {} ; + Counter m_pow {} ; // ====================================================================== }; // ======================================================================== @@ -1557,8 +1689,19 @@ namespace Ostap { public: // ====================================================================== - // constructir from the power - LehmerMean ( const double p = 1 ) ; + /// the actual underlying counter + typedef Moment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// defautl construictor + LehmerMean + ( const double = 1 ) ; + /// constructor for the counters + LehmerMean + ( const double p , + const Counter& cnt1 , + const Counter& cnt2 ) ; // ====================================================================== public: // ====================================================================== @@ -1571,6 +1714,13 @@ namespace Ostap inline LehmerMean& operator+=( const double x ) { return add ( x ) ; } inline LehmerMean& operator+=( const LehmerMean& x ) { return add ( x ) ; } // ====================================================================== + public: + // ====================================================================== + /// the counter of x**p values + const Counter counter1 () const { return m_lp ; } + /// the counter of x**(p-1) values + const Counter counter2 () const { return m_lpm1 ; } + // ====================================================================== public: // ====================================================================== /// accumulate only positive entries @@ -1624,6 +1774,18 @@ namespace Ostap */ class WGeometricMean : public WStatistic { + public: + // ====================================================================== + /// the actual underlying counter + typedef WMoment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// default constructor + WGeometricMean () = default ; + /// constructor for the counters + WGeometricMean ( const Counter& cnt ) ; + // ====================================================================== public: // ====================================================================== /// get the geometric mean @@ -1635,6 +1797,11 @@ namespace Ostap inline WGeometricMean& operator+=( const WGeometricMean& x ) { return add ( x ) ; } inline WGeometricMean& operator*=( const WGeometricMean& x ) { return add ( x ) ; } // ====================================================================== + public: + // ====================================================================== + /// the counter of log2 values + const Counter& counter () const { return m_log ; } + // ====================================================================== public: // ====================================================================== /// accumulate only positive entries @@ -1675,7 +1842,7 @@ namespace Ostap private: // ====================================================================== /// get the counter of log2(x) - WMoment_<1> m_log {} ; + Counter m_log {} ; // ====================================================================== }; @@ -1688,6 +1855,18 @@ namespace Ostap */ class WHarmonicMean : public WStatistic { + public: + // ====================================================================== + /// the actual underlying counter + typedef WMoment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// default construictor + WHarmonicMean () = default ; + /// constructor for the counters + WHarmonicMean ( const Counter& cnt ) ; + // ====================================================================== public: // ====================================================================== /// get the harmonic mean @@ -1698,6 +1877,11 @@ namespace Ostap // ====================================================================== inline WHarmonicMean& operator+=( const WHarmonicMean& x ) { return add ( x ) ; } // ====================================================================== + public: + // ====================================================================== + /// the counter of 1/x values + const Counter counter () const { return m_inv ; } + // ====================================================================== public: // ====================================================================== /// accumulate only non-zero entries @@ -1736,7 +1920,7 @@ namespace Ostap private: // ====================================================================== /// get the counter of 1/x - WMoment_<1> m_inv {} ; + Counter m_inv {} ; // ====================================================================== }; // ======================================================================== @@ -1748,8 +1932,15 @@ namespace Ostap { public: // ====================================================================== - // constructir from the power - WPowerMean ( const double p = 1 ) ; + /// the actual underlying counter + typedef WMoment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// defautl construictor + WPowerMean ( const double = 1 ) ; + /// constructor for the counters + WPowerMean ( const double p , const Counter& cnt ) ; // ====================================================================== public: // ====================================================================== @@ -1761,6 +1952,11 @@ namespace Ostap // ====================================================================== inline WPowerMean& operator+=( const WPowerMean& x ) { return add ( x ) ; } // ====================================================================== + public: + // ====================================================================== + /// the counter of x**p values + const Counter counter () const { return m_pow ; } + // ====================================================================== public: // ====================================================================== /// accumulate only positive entries @@ -1797,9 +1993,9 @@ namespace Ostap private: // ====================================================================== /// the power - double m_p {1} ; + double m_p {1} ; /// get the counter of x^p - WMoment_<1> m_pow {} ; + Counter m_pow {} ; // ====================================================================== }; // ======================================================================== @@ -1811,8 +2007,19 @@ namespace Ostap { public: // ====================================================================== - // constructir from the power - WLehmerMean ( const double p = 1 ) ; + /// the actual underlying counter + typedef WMoment_<1> Counter ; + // ====================================================================== + public: + // ====================================================================== + /// defautl construictor + WLehmerMean + ( const double = 1 ) ; + /// constructor for the counters + WLehmerMean + ( const double p , + const Counter& cnt1 , + const Counter& cnt2 ) ; // ====================================================================== public: // ====================================================================== @@ -1824,6 +2031,13 @@ namespace Ostap // ====================================================================== inline WLehmerMean& operator+=( const WLehmerMean& x ) { return add ( x ) ; } // ====================================================================== + public: + // ====================================================================== + /// the counter of x**p values + const Counter counter1 () const { return m_lp ; } + /// the counter of x**(p-1) values + const Counter counter2 () const { return m_lpm1 ; } + // ====================================================================== public: // ====================================================================== /// accumulate only positive entries with non-zero weight @@ -1860,11 +2074,11 @@ namespace Ostap private: // ====================================================================== /// the power - double m_p {1} ; + double m_p {1} ; /// get the counter of x^p - WMoment_<1> m_lp {} ; + Counter m_lp {} ; /// get the counter of x^(p-1) - WMoment_<1> m_lpm1 {} ; + Counter m_lpm1 {} ; // ====================================================================== }; // ======================================================================== @@ -1875,164 +2089,65 @@ namespace Ostap { public: // ====================================================================== - inline double value() const { return this->mean () ; } + typedef Moment_<1> Counter ; // ====================================================================== - } ; - // ======================================================================== - /** @class WArithmeticMean - * Calculate the weighted arithmetic mean - */ - class WArithmeticMean : public WMoment_<1> - { public: // ====================================================================== - inline double value() const { return this->mean () ; } + ArithmeticMean() = default ; + ArithmeticMean ( const Counter& cnt ) ; // ====================================================================== - } ; - // ======================================================================== - /** @class MinMean - * Degenerate case of "mean" : minimal value - * \f$ min ( x_1 , .. , x_n ) \f$ - */ - class MinValue : public Statistic - { public: // ====================================================================== - /// default constructor - MinValue () ; - // ====================================================================== - public: - // ====================================================================== - /// get the min mean - inline double value () 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 ) ; } } + inline double value() const { return this->mean () ; } // ====================================================================== 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 {} ; + Counter counter() const { return *this ; } // ====================================================================== } ; // ======================================================================== - /** @class MaxValue - * Degenerate case of "mean" : maximal value - * \f$ max ( x_1 , .. , x_n ) \f$ + /** @class WArithmeticMean + * Calculate the weighted arithmetic mean */ - class MaxValue : public Statistic + class WArithmeticMean : public WMoment_<1> { public: // ====================================================================== - /// default constructor - MaxValue () ; + typedef WMoment_<1> Counter ; // ====================================================================== public: // ====================================================================== - /// get the max value - inline double value () const { return m_max ; } - /// get the max value - inline double max () const { return m_max ; } + WArithmeticMean() = default ; + WArithmeticMean ( const Counter& cnt ) ; // ====================================================================== 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 ) ; } + inline double value() const { return this->mean () ; } // ====================================================================== 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 {} ; + Counter counter() const { return *this ; } // ====================================================================== - } ; + } ; // ======================================================================== /** @class MinMaxValue */ class MinMaxValue : public Statistic { // ====================================================================== + public: + // ====================================================================== + typedef Moment_<0> Counter ; + // ====================================================================== public: // ====================================================================== /// default constructor MinMaxValue () ; + /// constructor from min/max & counter + MinMaxValue + ( const double min , + const double max , + const Counter& cnt ) ; // ====================================================================== public: // ====================================================================== @@ -2072,6 +2187,10 @@ namespace Ostap inline void add ( ITERATOR begin , ITERATOR end ) { for ( ; begin != end ; ++begin ) { this->add ( *begin ) ; } } // ====================================================================== + public: + // ====================================================================== + const Counter& counter() const { return m_cnt ; } + // ====================================================================== public: // ====================================================================== /// to use it as a Moment @@ -2093,7 +2212,7 @@ namespace Ostap /// maximal value ; double m_max ; /// get the counter - Moment_<0> m_cnt {} ; + Counter m_cnt {} ; // ====================================================================== } ; // ======================================================================== @@ -2102,10 +2221,19 @@ namespace Ostap class WMinMaxValue : public WStatistic { // ====================================================================== + public: + // ====================================================================== + typedef WMoment_<0> Counter ; + // ====================================================================== public: // ====================================================================== /// default constructor WMinMaxValue () ; + /// constructor from min/max & counter + WMinMaxValue + ( const double min , + const double max , + const Counter& cnt ) ; // ====================================================================== public: // ====================================================================== @@ -2154,6 +2282,10 @@ namespace Ostap ( const double x , const double w = 1 ) override { add ( x , w ) ; } // ====================================================================== + public: + // ====================================================================== + const Counter& counter() const { return m_cnt ; } + // ====================================================================== public: // ====================================================================== /// number of entries @@ -2176,7 +2308,7 @@ namespace Ostap /// maximal value ; double m_max ; /// get the counter - WMoment_<0> m_cnt {} ; + Counter m_cnt {} ; // ====================================================================== } ; // ======================================================================== diff --git a/source/include/Ostap/NStatEntity.h b/source/include/Ostap/NStatEntity.h index dddc3d0c..eb4fa554 100644 --- a/source/include/Ostap/NStatEntity.h +++ b/source/include/Ostap/NStatEntity.h @@ -31,7 +31,14 @@ namespace Ostap public: // ====================================================================== /// constructor with N-parameter - NStatEntity ( const unsigned long N = 1000 ) ; + NStatEntity + ( const unsigned long N = 1000 ) ; + // ====================================================================== + /// constructor from two counters and sliding window + NStatEntity + ( const unsigned long N , + const StatEntity& cnt1 , + const StatEntity& cnt2 ) ; // ====================================================================== public: // ====================================================================== @@ -120,11 +127,11 @@ namespace Ostap private: // ====================================================================== /// the first counter - StatEntity m_cnt1 ; // the first counter + StatEntity m_cnt1 {} ; // the first counter /// the second counter - StatEntity m_cnt2 ; // the second counter + StatEntity m_cnt2 {} ; // the second counter /// the sliding window - unsigned long m_N ; // the sliding window + unsigned long m_N ; // the sliding window // ====================================================================== }; // ======================================================================== diff --git a/source/include/Ostap/WStatEntity.h b/source/include/Ostap/WStatEntity.h index be5bea97..bb28a693 100644 --- a/source/include/Ostap/WStatEntity.h +++ b/source/include/Ostap/WStatEntity.h @@ -24,12 +24,19 @@ namespace Ostap /// empty constructor WStatEntity () = default ; /// constructor from StatEntity of values - WStatEntity ( const StatEntity& values ) ; + WStatEntity \ + ( const StatEntity& values ) ; + /// full consrctor + WStatEntity + ( const double mu , + const double mu2 , + const StatEntity& values , + const StatEntity& weights ) ; // ====================================================================== public: // the basic getters // ====================================================================== /// total number of entries - unsigned long long n () const { return m_weights.n() ; } + unsigned long long n () const { return m_weights.n () ; } /// the first weighted moment/mean-value double mu () const { return m_mu ; } /// the second central weighted moment/dispersion/variance diff --git a/source/src/Moments.cpp b/source/src/Moments.cpp index e052c137..c03deba6 100644 --- a/source/src/Moments.cpp +++ b/source/src/Moments.cpp @@ -42,8 +42,158 @@ Ostap::Math::WMoment::~WMoment(){} // =========================================================================== + +// =========================================================================== +// constructor +// =========================================================================== +/* full constructor + * @param size number of entries + * @param sumw sum of weights + * @param sumw sum of squared weights + */ +// =========================================================================== +Ostap::Math::WMoment_<0>::WMoment_ +( const unsigned long long size , + const double sumw , + const double sumw2 ) + : m_size ( size ) + , m_w ( sumw ) + , m_w2 ( sumw2 ) +{ + if ( 0 == m_size ) + { + Ostap::Assert ( s_zero ( m_w ) && s_zero ( m_w2 ) , + "Non-zero sum of (squared) weights for empty counter!" , + "Ostap::Math::WMoment_<0>" ) ; + m_w = 0 ; + m_w2 = 0 ; + } + // + if ( s_zero ( m_w2 ) ) + { + Assert ( s_zero ( m_w ) , + "Non zero sum of weigth for zero sumw2 !" , + "Ostap::Math::WMoment_<0>" ) ; + m_w = 0 ; + m_w2 = 0 ; + } + // + if ( s_zero ( m_w ) ) { m_w = 0 ; } + // + Ostap::Assert ( 0 <= m_w2 , + "Negatoive sum of squared weights!" , + "Ostap::Math::WMoment_<0>" ) ; + +} +// =========================================================================== + + + +// =========================================================================== +Ostap::Math::GeometricMean::GeometricMean +( const Ostap::Math::GeometricMean::Counter& cnt ) + : m_log ( cnt ) +{} +// =========================================================================== +Ostap::Math::HarmonicMean::HarmonicMean +( const Ostap::Math::HarmonicMean::Counter& cnt ) + : m_inv ( cnt ) +{} +// =========================================================================== +Ostap::Math::ArithmeticMean::ArithmeticMean +( const Ostap::Math::ArithmeticMean::Counter& cnt ) + : Counter ( cnt ) +{} +// =========================================================================== +Ostap::Math::PowerMean::PowerMean +( const double p ) + : m_p ( p ) + , m_pow ( ) +{} +// =========================================================================== +Ostap::Math::PowerMean::PowerMean +( const double p , + const Ostap::Math::PowerMean::Counter& cnt ) + : m_p ( p ) + , m_pow ( cnt ) +{} +// =========================================================================== +Ostap::Math::LehmerMean::LehmerMean +( const double p ) + : m_p ( p ) + , m_lp ( ) + , m_lpm1 ( ) +{} +// =========================================================================== +Ostap::Math::LehmerMean::LehmerMean +( const double p , + const Ostap::Math::LehmerMean::Counter& cnt1 , + const Ostap::Math::LehmerMean::Counter& cnt2 ) + : m_p ( p ) + , m_lp ( cnt1 ) + , m_lpm1 ( cnt2 ) +{ + Ostap::Assert ( m_lp.size () == m_lpm1.size () , + "Inconsistent structure of two conuters!" , + "Ostap::Math::LehmerMean" ) ; + +} +// =========================================================================== +Ostap::Math::WGeometricMean::WGeometricMean +( const Ostap::Math::WGeometricMean::Counter& cnt ) + : m_log ( cnt ) +{} +// =========================================================================== +Ostap::Math::WHarmonicMean::WHarmonicMean +( const Ostap::Math::WHarmonicMean::Counter& cnt ) + : m_inv ( cnt ) +{} +// =========================================================================== +Ostap::Math::WArithmeticMean::WArithmeticMean +( const Ostap::Math::WArithmeticMean::Counter& cnt ) + : Counter ( cnt ) +{} +// =========================================================================== +Ostap::Math::WPowerMean::WPowerMean +( const double p ) + : m_p ( p ) + , m_pow ( ) +{} +// =========================================================================== +Ostap::Math::WPowerMean::WPowerMean +( const double p , + const Ostap::Math::WPowerMean::Counter& cnt ) + : m_p ( p ) + , m_pow ( cnt ) +{} +// =========================================================================== +Ostap::Math::WLehmerMean::WLehmerMean +( const double p ) + : m_p ( p ) + , m_lp ( ) + , m_lpm1 ( ) +{} +// =========================================================================== +Ostap::Math::WLehmerMean::WLehmerMean +( const double p , + const Ostap::Math::WLehmerMean::Counter& cnt1 , + const Ostap::Math::WLehmerMean::Counter& cnt2 ) + : m_p ( p ) + , m_lp ( cnt1 ) + , m_lpm1 ( cnt2 ) +{ + Ostap::Assert ( m_lp.size () == m_lpm1.size () , + "Inconsistent structure of two conuters!" , + "Ostap::Math::WLehmerMean" ) ; + Ostap::Assert ( s_equal ( m_lp.w () , m_lpm1.w () ) , + "Inconsistent structure of two conuters!" , + "Ostap::Math::WLehmerMean" ) ; + Ostap::Assert ( s_equal ( m_lp.w2 () , m_lpm1.w2 () ) , + "Inconsistent structure of two conuters!" , + "Ostap::Math::WLehmerMean" ) ; +} // =========================================================================== -// HarmoincMean:: accumulate only non-zero entries +// HarmoinicMean:: accumulate only non-zero entries // =========================================================================== Ostap::Math::HarmonicMean& Ostap::Math::HarmonicMean::add ( const double x ) @@ -64,12 +214,6 @@ Ostap::Math::WHarmonicMean::add } -// =========================================================================== -Ostap::Math::PowerMean::PowerMean -( const double p ) - : m_p ( p ) - , m_pow ( ) -{} // =========================================================================== // add two counters togather if p is common // =========================================================================== @@ -82,13 +226,6 @@ Ostap::Math::PowerMean::add ( const Ostap::Math::PowerMean& x ) 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 // =========================================================================== @@ -110,21 +247,6 @@ Ostap::Math::WPowerMean::add 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 // =========================================================================== @@ -165,22 +287,6 @@ Ostap::Math::WLehmerMean::add } 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 // =========================================================================== @@ -190,6 +296,24 @@ Ostap::Math::MinMaxValue::MinMaxValue () , m_cnt () {} // =========================================================================== +// constructor +// =========================================================================== +Ostap::Math::MinMaxValue::MinMaxValue +( const double minv , + const double maxv , + const Ostap::Math::MinMaxValue::Counter& cnt ) + : m_min ( minv ) + , m_max ( maxv ) + , m_cnt () +{ + Ostap::Assert ( ( empty() + && m_min == std::numeric_limits::max() + && m_max == - std::numeric_limits::max() ) + || ( !empty() && m_min <= m_max ) , + "Invalid min/max/empty structure!" , + "Ostap::Math::MinMaxValue!" ) ; +} +// =========================================================================== // default constructor // =========================================================================== Ostap::Math::WMinMaxValue::WMinMaxValue () @@ -198,31 +322,23 @@ Ostap::Math::WMinMaxValue::WMinMaxValue () , m_cnt () {} // ====================================================================== - - -// int test() -// { -// Ostap::Math::Moment_<5> m1 ; -// for ( unsigned i = 0 ; i < 100 ; ++i ) { m1 += i ; } -// // -// m1.M(0) ; -// m1.M(1) ; -// m1.M(2) ; -// Ostap::Math::Moment_<5> m2 ; -// for ( unsigned i = 0 ; i < 100 ; ++i ) { m2 += i ; } -// // -// m2.M(0) ; -// m2.M(1) ; -// m2.M(2) ; -// Ostap::Math::Moment_<2> m4 ; -// for ( unsigned i = 0 ; i < 100 ; ++i ) { m4 += i ; } -// const auto m3 = m1 + m2 ; -// const auto a1 = Ostap::Math::Moments::central_moment<2> ( m3 ) ; -// const auto a2 = Ostap::Math::Moments::central_moment<2> ( m4 ) ; -// return 0 ; -// } - - +// constructor +// =========================================================================== +Ostap::Math::WMinMaxValue::WMinMaxValue +( const double minv , + const double maxv , + const Ostap::Math::WMinMaxValue::Counter& cnt ) + : m_min ( minv ) + , m_max ( maxv ) + , m_cnt () +{ + Ostap::Assert ( ( empty() + && m_min == std::numeric_limits::max() + && m_max == - std::numeric_limits::max() ) + || ( !empty() && m_min <= m_max ) , + "Invalid min/max/empty structure!" , + "Ostap::Math::WMinMaxValue!" ) ; +} // ============================================================================ // The END // ============================================================================ diff --git a/source/src/NStatEntity.cpp b/source/src/NStatEntity.cpp index fba2d9f7..d8904ae7 100644 --- a/source/src/NStatEntity.cpp +++ b/source/src/NStatEntity.cpp @@ -13,6 +13,10 @@ #include "Ostap/StatEntity.h" #include "Ostap/NStatEntity.h" // ============================================================================ +// local +// ============================================================================ +#include "Exception.h" +// ============================================================================ /** @file * Implementation file for class Ostap::NStatEntity * @see Ostap::NStatEntity @@ -20,10 +24,6 @@ * @see Ostap::StatEntity * @author Vanya Belyaev Ivan.Belyaev@itep.ru * @date 2015-04-03 - * - * $Revision$ - * Last modification $Date$ - * by $Author$ */ // =========================================================================== namespace @@ -39,11 +39,41 @@ namespace // =========================================================================== // constructor with N-parameter // =========================================================================== -Ostap::NStatEntity::NStatEntity ( const unsigned long N ) - : m_cnt1 () - , m_cnt2 () - , m_N ( std::min ( std::max ( 1UL , N ) , s_max ) ) -{} +Ostap::NStatEntity::NStatEntity +( const unsigned long N ) + : m_cnt1 ( ) + , m_cnt2 ( ) + , m_N ( N ) +{ + Ostap::Assert ( 2 <= m_N , + "Invalid window size!" , + "Ostap::NstatEntity" ) ; + +} +// =========================================================================== +// constructor from two counters and sliding window +// =========================================================================== +Ostap::NStatEntity::NStatEntity +( const unsigned long N , + const StatEntity& cnt1 , + const StatEntity& cnt2 ) + : m_cnt1 ( cnt1 ) + , m_cnt2 ( cnt2 ) + , m_N ( N ) +{ + Ostap::Assert ( 2 <= m_N , + "Invalid window size!" , + "Ostap::NStatEntity" ) ; + Ostap::Assert ( m_cnt1.nEntries () <= 2 * m_N , + "Invalid first counter!" , + "Ostap::NstatEntity" ) ; + Ostap::Assert ( m_cnt2.nEntries () <= 2 * m_N , + "Invalid second counter!" , + "Ostap::NStatEntity" ) ; +} +// ====================================================================== + + // =========================================================================== /* printout to std::ostream * @param s the reference to the output stream diff --git a/source/src/StatEntity.cpp b/source/src/StatEntity.cpp index 89b75066..0ead191b 100644 --- a/source/src/StatEntity.cpp +++ b/source/src/StatEntity.cpp @@ -13,6 +13,7 @@ // ============================================================================ // Ostap // ============================================================================ +#include "Ostap/Math.h" #include "Ostap/StatEntity.h" // ============================================================================ // Local @@ -26,6 +27,13 @@ * @author Vanya BELYAEV Ivan.Belyaev@itep.ru */ // ============================================================================ +namespace +{ + // ========================================================================== + const Ostap::Math::Zero s_zero {} ; + // ========================================================================== +} +// ============================================================================ /* The full contructor from all important values * @see StatEntity::format * @param entries number of entries @@ -47,17 +55,34 @@ Ostap::StatEntity::StatEntity , m_min ( minv ) , m_max ( maxv ) { - // reset empty counter - if ( 0 == m_n ) { reset() ; } - else - { - Ostap::Assert ( m_min <= mu && mu <= m_max , - "Ostap::StatEntity: invalid minv/mu/maxv" , - "Ostap::StatEntity" ) ; - Ostap::Assert ( 0 <= m_mu2 , - "Ostap::StatEntity: invalid second moment", - "Ostap::StatEntity" ) ; - } + // empty counter (ignore min/max) + if ( 0 == m_n ) + { + Ostap::Assert ( s_zero ( m_mu ) && s_zero ( m_mu2 ) , + "Ostap::StatEntity: invalid mu/mu2 for empty counter!" , + "Ostap::StatEntity" ) ; + m_mu = 0 ; + m_mu2 = 0 ; + /// redefine/ignore min/max + m_min = std::numeric_limits::max() ; + m_max = - std::numeric_limits::max() ; + } + else + { + Ostap::Assert ( m_min <= mu && mu <= m_max , + "Ostap::StatEntity: invalid minv/mu/maxv" , + "Ostap::StatEntity" ) ; + } + // + if ( s_zero ( m_mu2 ) ) { m_mu2 = 0 ; } + // + Ostap::Assert ( ( empty() && !m_mu2 ) || ( !empty() && m_mu2 ) , + "Ostap::StatEntity: inconsistent mu2/empty!" , + "Ostap::StatEntity" ) ; + // + Ostap::Assert ( 0 <= m_mu2 , + "Ostap::StatEntity: invalid second moment", + "Ostap::StatEntity" ) ; } // ============================================================================ /* add a value : the main method diff --git a/source/src/WStatEntity.cpp b/source/src/WStatEntity.cpp index d269c411..58547f68 100644 --- a/source/src/WStatEntity.cpp +++ b/source/src/WStatEntity.cpp @@ -14,6 +14,7 @@ // local // ============================================================================ #include "format.h" +#include "Exception.h" // ============================================================================ /** @file * Implementation file for class Ostap::WStatEntity @@ -38,6 +39,42 @@ Ostap::WStatEntity::WStatEntity ( const Ostap::StatEntity& values ) , m_weights ( values.n () , 1 , 0 , 1 , 1 ) // weights are trivial {} // ============================================================================ +// full constructor +// ============================================================================ +Ostap::WStatEntity::WStatEntity +( const double mu , + const double mu2 , + const Ostap::StatEntity& values , + const Ostap::StatEntity& weights ) + : m_mu ( mu ) + , m_mu2 ( mu2 ) + , m_values ( values ) + , m_weights ( weights ) +{ + if ( empty () ) + { + Ostap::Assert ( s_zero ( m_mu ) && s_zero ( m_mu2 ) , + "Ostap::WStatWEntity: invalid mu/mu2 for empty counter!" , + "Ostap::WStatWEntity" ) ; + m_mu = 0 ; + m_mu2 = 0 ; + } + // + Ostap::Assert ( m_values.n () <= m_weights.n () , + "Ostap::WStatWEntity: inconsistent values/weights counters!" , + "Ostap::WStatWEntity" ) ; + // + if ( s_zero ( m_mu2 ) ) { m_mu2 = 0 ; } + // + Ostap::Assert ( ( empty() && !m_mu2 ) || ( !empty() && m_mu2 ) , + "Ostap::WStatEntity: inconsistent mu2/emptty!" , + "Ostap::WStatEntity" ) ; + // + Ostap::Assert ( 0 <= m_mu2 , + "Ostap::WStatEntity: invalid second moment", + "Ostap::WStatEntity" ) ; +} +// ============================================================================ // update statistics // ============================================================================ Ostap::WStatEntity& @@ -50,14 +87,14 @@ Ostap::WStatEntity::add if ( !std::isfinite ( value ) || !std::isfinite ( weight ) ) { return *this ; } // if ( 0 == n() ) - { - m_mu = value ; - m_mu2 = 0 ; - if ( weight ) { m_values += value ; } - m_weights += weight ; - // - return *this ; - } + { + m_mu = value ; + m_mu2 = 0 ; + if ( weight ) { m_values += value ; } + m_weights += weight ; + // + return *this ; + } // const long double wA = n () * m_weights.mu() ; const long double wB = weight ;