From 638f0284e9221b24728b77d2cb402e5cc5d45c35 Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Tue, 27 Aug 2024 07:32:52 +0200 Subject: [PATCH] 1. improve pretty-print for matrices 1. add `pretty_array` function for nice print of arrays/sequences 1. add some auxillary methods for matrices: `(min/max/minabs/maxabs)_element` and `(min/max)_diagonal` and `(min/max/minabs/maxabs)_element_index` --- CMakeLists.txt | 2 +- ReleaseNotes/release_notes.md | 11 +- ReleaseNotes/v1.13.0.0.md | 20 -- ostap/core/ostap_types.py | 72 +++--- ostap/logger/pretty.py | 87 ++++++- ostap/math/base.py | 183 ++++++++------- ostap/math/linalg2.py | 338 ++++++++++++++++++++++------ ostap/math/tests/test_math_base.py | 68 ++++++ source/include/Ostap/MatrixUtils.h | 140 ++++++------ source/include/Ostap/MatrixUtilsT.h | 40 +++- 10 files changed, 678 insertions(+), 283 deletions(-) delete mode 100644 ReleaseNotes/v1.13.0.0.md create mode 100644 ostap/math/tests/test_math_base.py diff --git a/CMakeLists.txt b/CMakeLists.txt index bdedf60f..052422bb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -23,7 +23,7 @@ include(CTest) set(OSTAP_VERSION_MAJOR 1) set(OSTAP_VERSION_MINOR 13) set(OSTAP_VERSION_PATCH 0) -set(OSTAP_VERSION_TWEAK 0) +set(OSTAP_VERSION_TWEAK 1) set(OSTAP_VERSION ${OSTAP_VERSION_MAJOR}.${OSTAP_VERSION_MINOR}.${OSTAP_VERSION_PATCH}.${OSTAP_VERSION_TWEAK}) diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index cddd85d2..0711aa66 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -1,3 +1,13 @@ +## New features + + 1. improve prettry-print for matrices + 1. add `pretty_array` function for nice print of arrays/sequences + 1. add some auxillary methods for matrices: `(min/max/minabs/maxabs)_element` and `(min/max)_diagonal` and `(min/max/minabs/maxabs)_element_index` + +## Backward incompatible + +## Bug fixes + # v1.13.0.0 ## New features @@ -19,7 +29,6 @@ 1. fix a minor bug in `bernstein.solve` 1. fix couple of recent bugs in `histos.graphs` ` - # v1.12.0.0 ## New features diff --git a/ReleaseNotes/v1.13.0.0.md b/ReleaseNotes/v1.13.0.0.md deleted file mode 100644 index bb19da2c..00000000 --- a/ReleaseNotes/v1.13.0.0.md +++ /dev/null @@ -1,20 +0,0 @@ -# v1.13.0.0 - -## New features - - 1. sight improvements for `bernstein` - 1. disable Thiele rational interpolation from the tests. Sometnjug wrong with the code. - 1. extend a bit functionality of asymmetric errors (needed for graphs&plots) - 1. collect pretty-print functions into new module `ostap.logger.pretty` - 1. extend and improve pretty-print stuff - 1. make the dbase.Item a bit lighter - 1. For shelve-like databases: allow to specify the preferred backend as `dbtype = ...` - 1. reimplement `SqliteShelf` as `ZipShelf` with `dbtype='sqlite'` - 1. eliminate messy stuff with extensions for `XXXShelf``-classes - -## Backward incompatible - -## Bug fixes - - 1. fix a minor bug in `bernstein.solve` - 1. fix couple of recent bugs in `histos.graphs` ` diff --git a/ostap/core/ostap_types.py b/ostap/core/ostap_types.py index 5569e6a7..02c2f2c9 100644 --- a/ostap/core/ostap_types.py +++ b/ostap/core/ostap_types.py @@ -43,8 +43,9 @@ 'all_strings' , ## all argumets of string types? ) # ============================================================================= -import math, os -from sys import version_info as python_version +import math, os, array +from sys import version_info as python_version +# ============================================================================= if ( 3 , 0 ) <= python_version : long = int string_types = bytes , str @@ -56,17 +57,17 @@ integer_types = int , long long_type = long python_version = 2 - import collections as C + import collections as C # ============================================================================= if ( 3 , 5 ) <= python_version : - from collections.abc import Collection, Sequence, Iterable, Sized, Generator + from collections.abc import Collection, Sequence, Iterable, Sized, Generator elif ( 3 , 3 ) <= python_version : - from collections.abc import Collection, Sequence, Iterable, Sized - from types import GeneratorType as Generator + from collections.abc import Collection, Sequence, Iterable, Sized + from types import GeneratorType as Generator else : - from collections import Sequence , Iterable , Sized - from collections import Container as Collection - from types import GeneratorType as Generator + from collections import Sequence , Iterable , Sized + from collections import Container as Collection + from types import GeneratorType as Generator # ============================================================================= # logging @@ -81,7 +82,6 @@ num_types = integer_types + ( float , ) str_types = str, list_types = list , tuple -import array listlike_types = list_types + ( set , C.Sequence , array.array ) # ============================================================================= try : @@ -95,10 +95,10 @@ sequence_types = listlike_types + ( Sequence , Collection , Iterable , Generator ) sized_types = Sized , path_types = string_types -if (3,6) <= python_version : +if ( 3 , 6 ) <= python_version : path_types = string_types + ( os.PathLike , ) # ============================================================================= -## sometimewe we need to ensure that dictionary is ordered +## sometimes we need to ensure that dictionary is ordered ordered_dict = dict if python_version < ( 3 , 7 ) : from collections import OrderedDict as ordered_dict @@ -108,26 +108,35 @@ # ============================================================================= ## Is this number of a proper integer? def is_integer ( v ) : - """Is this number of a proper integer?""" + """ Is this number of a proper integer?""" return isinstance ( v , integer_types ) # ============================================================================= ## Is this number of a proper numeric type def is_number ( v ) : - """Is this number of a proper numeric?""" + """ Is this number of a proper numeric?""" return isinstance ( v , num_types ) # ============================================================================= -## good numeric value -def is_good_number ( v ) : - """Is numeric type and good value?""" - return isinstance ( v , num_types ) and \ - ( not math.isinf ( v ) ) and ( not math.isnan ( v ) ) +if ( 3 , 2 ) <= python_version : + # ========================================================================= + ## good numeric value + def is_good_number ( v ) : + """ Is numeric type and good value?""" + return isinstance ( v , num_types ) and math.isfinite ( v ) + # ========================================================================= +else : + # ========================================================================= + ## good numeric value + def is_good_number ( v ) : + """ Is numeric type and good value?""" + return isinstance ( v , num_types ) and \ + ( not math.isinf ( v ) ) and ( not math.isnan ( v ) ) # ============================================================================= ## is value of str-type? def is_string ( v ) : - """Is value of str-type?""" + """ Is value of str-type?""" return isinstance ( v , str_types ) # ============================================================================= @@ -139,42 +148,35 @@ def is_string_like ( v ) : # ============================================================================= ## is list type? def is_list ( v ) : - """Is value of list type (list ot tuple)""" + """ Is value of list type (list ot tuple)""" return isinstance ( v , list_type ) # ============================================================================= ## is list-like type? def is_list_like ( v ) : - """Is value of list-like type (list but not a string-like!)""" + """ Is value of list-like type (list but not a string-like!)""" return isinstance ( v , listlike_type ) and not isinstance ( v , string_types ) # ============================================================================= ## are all values of integer-line types? def all_integers ( *args ) : - """Are all value of integer-line types? + """ Are all value of integer-line types? """ - for a in args : - if not isinstance ( a , integer_types ) : return False - return True + return all ( isinstance ( a , integer_types ) for a in args ) # ============================================================================= ## are all values of numeric types? def all_numerics ( *args ) : - """Are all values of numeric types? + """ Are all values of numeric types? """ - for a in args : - if not isinstance ( a , num_types ) : return False - return True + return all ( isinstance ( a , num_types ) for a in args ) # ============================================================================= ## are all values of string types? def all_strings ( *args ) : - """Are all values of string types? + """ Are all values of string types? """ - for a in args : - if not isinstance ( a , stgring_types ) : return False - return True - + return all ( isinstance ( a , string_types ) for a in args ) # ============================================================================= if '__main__' == __name__ : diff --git a/ostap/logger/pretty.py b/ostap/logger/pretty.py index 3f1a3a5f..f5157f40 100644 --- a/ostap/logger/pretty.py +++ b/ostap/logger/pretty.py @@ -35,12 +35,14 @@ 'add_expo' , ## add an exponetaion factor for sting representaion of the object ) # ============================================================================= -from ostap.logger.utils import ( pretty_float , - pretty_err1 , - pretty_err2 , - pretty_errs , - add_expo , - fmt_pretty_err1 ) +from ostap.core.ostap_types import num_types, sequence_types +from ostap.logger.utils import ( pretty_float , + pretty_err1 , + pretty_err2 , + pretty_errs , + add_expo , + fmt_pretty_err1 , + fmt_pretty_float ) # ============================================================================= # logging # ============================================================================= @@ -95,6 +97,79 @@ def pretty_ve ( value , ## delegate return pretty_err1 ( v , e , width = width , precision = precision ) +# ======================================================================= +## nice printout of complex number ( string + exponent) +# @code +# ae = complex ( ... ) +# s , expo = pretty_complex ( ae ) +# @endcode +# @return nice string and the separate exponent +def pretty_complex ( value , + width = 6 , + precision = 4 , + parentheses = True ) : + """ Nice printout of complex number ( string + exponent) + - return nice stirng and the separate exponent + >>> s , expo = pretty_complex( number ) + """ + v = complex ( value ) + vv = max ( abs ( v.real() ) , abs ( v.imag ) ) + ## + fmtv , expo = fmt_pretty_float ( value = vv , + width = width , + precision = precision , + parentheses = False ) + ## + fmt = '%s%sj' % ( fmtv , fmtv ) + if parentheses : fmt = '(' + fmt + ')' + ## + if expo : + scale = 10 ** expo + v /= scale + return fmt % ( v.real , v.imag ) , expo + +# ======================================================================= +## Nice pritout of arrays of numerical values +# @code +# array = ... +# result, expo = pretty_array ( array ) +# @endcode +# @return nice string and the separate exponent +def pretty_array ( value , + width = 6 , + precision = 4 , + parentheses = True ) : + """ Nice pritout of arrays of numerical values + - return nice stirng and the separate exponent + >>> array = ... + >>> result, expo = pretty_array ( array ) + """ + assert isinstance ( value , sequence_types ) , \ + "Invalid type of `value':%s" % type ( value ) + + if not value : return '' if not parentheses else [] + + assert all ( isinstance ( v , num_types ) for v in value ) , \ + "Invalid content of `value': %s" % str ( value ) + + value1 = value [ 0 ] + value2 = max ( abs ( v ) for v in value ) + vv = max ( abs ( value1 ) , abs ( value2 ) ) + # + fmtv , expo = fmt_pretty_float ( value = vv , + width = width , + precision = precision ) + + if expo : scale = 10 ** expo + else : scale = 1 + ## + result = [] + for v in value : result.append ( fmtv % ( v / scale ) ) + result = ', '.join ( result ) + ## + if parentheses : result = ' [ ' + result + ' ]' + return result, expo + # ======================================================================= ## nice printout of asymmetric errors ( string + exponent) # @code diff --git a/ostap/math/base.py b/ostap/math/base.py index 3abda84e..8c1ce96c 100755 --- a/ostap/math/base.py +++ b/ostap/math/base.py @@ -34,7 +34,7 @@ # @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl # @date 2009-09-12 # ============================================================================= -"""Simple file to provide 'easy' access in python for the basic ROOT::Math classes +""" Simple file to provide 'easy' access in python for the basic ROOT::Math classes see Ostap/Point3DTypes.h see Ostap/Vector3DTypes.h @@ -112,7 +112,7 @@ 'neg_infinity' , ## negative infinity ) # ============================================================================= -from ostap.core.meta_info import root_version_int +from ostap.core.meta_info import root_info import ROOT, cppyy, sys, math # ============================================================================= # logging @@ -128,14 +128,18 @@ import math if not hasattr ( math , 'inf' ) : math.inf = float('inf' ) if not hasattr ( math , 'nan' ) : math.nan = float('nan' ) - + +# ============================================================================= ## get global C++ namespace -cpp = cppyy.gbl -## C++ namespace Gaudi -std = cpp.std +cpp = cppyy.gbl +# ============================================================================= +## C++ namespace std +std = cpp.std +# ============================================================================= ## C++ namespace Ostap Ostap = cpp.Ostap +# ============================================================================= ## positive and negative infinities pos_infinity = float('+inf') neg_infinity = float('-inf') @@ -144,12 +148,12 @@ ## Get the type name # @code # obj = ... -# print ( 'Obnject type name is %s' % typename ( obj ) ) +# print ( 'Object type name is %s' % typename ( obj ) ) # @endcode def typename ( o ) : """Get the type name >>> obj = ... - >>> print ( 'Obnject type name is %s' % typename ( obj ) ) + >>> print ( 'Object type name is %s' % typename ( obj ) ) """ return type ( o ) .__name__ @@ -212,7 +216,7 @@ def __exit__ ( self , *_ ) : vLongs = std.vector ( 'long' ) # ============================================================================= -if (3,2) <= sys.version_info : +if ( 3 , 2 ) <= sys.version_info : ## local version of isfinite isfinite = math.isfinite else : @@ -224,17 +228,33 @@ def isfinite ( x ) : return ( not math.isinf ( y ) ) and ( not math.isnan ( y ) ) # ============================================================================= -if (3,5) <= sys.version_info : +if ( 3 , 5 ) <= sys.version_info : + # ========================================================================= ## local version of isclose isclose = math.isclose - ## def isclose ( a , b , rel_tol = 1.e-9 , abs_tol = 0.0 ) : - ## """Local version of `isclose`""" - ## return math.isclose ( a , b , rel_tol = rel_tol, abs_tol = abs_tol ) else : # ========================================================================= - ## local version of isclose + ## local version of math.isclose def isclose ( a , b , rel_tol = 1.e-9 , abs_tol = 0.0 ) : - """Local version of `isclose`""" + """ Local version of `math.isclose` + Determine whether two floating point numbers are close in value. + + - rel_tol + maximum difference for being considered "close", relative to the + magnitude of the input values + - abs_tol + maximum difference for being considered "close", regardless of the + magnitude of the input values + + Return True if a is close in value to b, and False otherwise. + + For the values to be considered close, the difference between them + must be smaller than at least one of the tolerances. + + -inf, inf and NaN behave similarly to the IEEE 754 Standard. That + is, NaN is not close to anything, even itself. inf and -inf are + only close to themselves. + """ return abs ( a - b ) <= max ( rel_tol * max ( abs ( a ) , abs ( b ) ) , abs_tol ) # ============================================================================= @@ -259,6 +279,7 @@ def samesign ( a , b ) : # ============================================================================= ## natural number ? natural_number = Ostap.Math.natural_number +# ============================================================================= ## natural etry in histo-bin ? natural_entry = Ostap.Math.natural_entry @@ -270,7 +291,7 @@ def samesign ( a , b ) : # print 'Is %s<=%s<=%s ? %s' % ( a , x , b , inrange ( a , x , b ) ) # @endcode def inrange ( a , x , b ) : - """Is float value in range? + """ Is float value in range? >>> x = 1.1 >>> a,b = 1,2 >>> print 'Is x between a and b? %s' % inrange ( a , x , b ) @@ -296,7 +317,7 @@ def inrange ( a , x , b ) : v = std.vector( t ) v.asList = lambda s : [ i for i in s ] ## convert vector into list v.toList = v.asList - v.asTuple = lambda s : tuple ( s.asList() ) + v.asTuple = lambda s : tuple ( i for i in s ) v.toTuple = v.asTuple v.__repr__ = lambda s : str ( [ i for i in s ] ) ## print it ! v.__str__ = lambda s : str ( [ i for i in s ] ) ## print it ! @@ -310,7 +331,7 @@ def _v_iter_ ( v ) : # ============================================================================= ## Convert vector of floating types to string def float_vct_str ( vct , format = '%.5g' ) : - """Convert vector of floating types to string""" + """ Convert vector of floating types to string""" try : return '[ ' + ', '.join ( [ format % v for v in vct ] ) + ' ]' except TypeError : @@ -320,7 +341,7 @@ def float_vct_str ( vct , format = '%.5g' ) : # ============================================================================= ## Convert vector of complex types to string def complex_vct_str ( vct , format = '%.5g%-+.5gj' ) : - """Convert vector of complex types to string""" + """ Convert vector of complex types to string""" try : lst = [] for c in vct : @@ -345,9 +366,9 @@ def complex_vct_str ( vct , format = '%.5g%-+.5gj' ) : v.__str__ = complex_vct_str # ============================================================================= -## self-printout of TMaxtrix +## self-printout of TMatrix def _tmg_str_ ( self , fmt = ' %+11.4g') : - """Self-printout of TMatrix + """ Self-printout of `TMatrix` """ _rows = self.GetNrows() _cols = self.GetNcols() @@ -364,11 +385,10 @@ def _tmg_str_ ( self , fmt = ' %+11.4g') : ROOT.TMatrix.__repr__ = _tmg_str_ ROOT.TMatrix.__str__ = _tmg_str_ - # ============================================================================= -## add something to std::vector +## add something to std::vector def _add_to ( vct , cnv , arg1 , *args ) : - """Add something to std::vector + """ Add something to `std.vector` """ from ostap.core.ostap_types import sequence_types @@ -431,7 +451,7 @@ def make_vector ( TYPE , cnv , *args ) : # ============================================================================= ## construct std::vector from the arguments def doubles ( arg1 , *args ) : - """Construct the std::vector from the arguments + """ Construct the std::vector from the arguments >>> v1 = doubles ( 1.01 ) >>> v2 = doubles ( 1.01 , 1.02 , 1.03 ) >>> v3 = doubles ( [ 1.01 , 1.02 , 1.03 ] ) @@ -441,7 +461,7 @@ def doubles ( arg1 , *args ) : # ============================================================================= ## construct std::vector from the arguments def ints ( arg1 , *args ) : - """Construct the std::vector from the arguments + """ Construct the std::vector from the arguments >>> v1 = ints ( 1 ) >>> v2 = ints ( 1 , 1 , 1 ) >>> v3 = ints ( [ 1 , 2 , 3 ] ) @@ -451,7 +471,7 @@ def ints ( arg1 , *args ) : # ============================================================================= ## construct std::vector from the arguments def uints ( arg1 , *args ) : - """Construct the std::vector from the arguments + """ Construct the std::vector from the arguments >>> v1 = uints ( 1 ) >>> v2 = uints ( 1 , 1 , 1 ) >>> v3 = uints ( [ 1 , 2 , 3 ] ) @@ -486,26 +506,22 @@ def strings ( *args ) : >>> v2 = strings( ['a','b','c'] ) >>> v3 = strings( ['a','b'],'c',['d','e'] ) """ - return make_vector ( 'std::string' , std.string , *args ) - ## VS = std.vector(std.string) - ## vs = VS() - ## if not args : return vs - ## return _add_to ( vs , std.string , args[0] , *args[1:] ) - -SPD = std.pair('double','double') +# ============================================================================= +SPD = std.pair ( 'double' , 'double' ) SPD.asTuple = lambda s : ( s.first , s.second ) SPD.__str__ = lambda s : str( ( s.first , s.second ) ) SPD.__repr__ = SPD.__str__ SPD.__len__ = lambda s : 2 def _spd_getitem_ ( s , i ) : - """Get item form the pair: + """Get item from the pair: >>> p = ... >>> p[0], p[1] """ - if 0 == i : return s.first - elif 1 == i : return s.second + if 0 == i : return s.first ## first + elif 1 == i : return s.second ## second + elif -1 == i : return s.second ## last raise IndexError('Invalid index %s' % i ) def _spd_setitem_ ( s , i , value ) : """Set item from the pair: @@ -513,8 +529,9 @@ def _spd_setitem_ ( s , i , value ) : >>> p[0] = 1 >>> p[1] = 2 """ - if 0 == i : s.first = value - elif 1 == i : s.second = value + if 0 == i : s.first = value ## first + elif 1 == i : s.second = value ## second + elif -1 == i : s.second = value ## last raise IndexError('Invalid index %s' % i ) SPD.__getitem__ = _spd_getitem_ SPD.__setitem__ = _spd_setitem_ @@ -522,9 +539,9 @@ def _spd_setitem_ ( s , i , value ) : # ============================================================================= # Improve operations with std.complex # ============================================================================= -COMPLEX = cpp.std.complex ( 'double' ) -COMPLEXf = cpp.std.complex ( 'float' ) -COMPLEXl = cpp.std.complex ( 'long double' ) +COMPLEX = cpp.std.complex ( 'double' ) +COMPLEXf = cpp.std.complex ( 'float' ) +COMPLEXl = cpp.std.complex ( 'long double' ) VCOMPLEX = cpp.std.vector ( COMPLEX ) VDOUBLE = cpp.std.vector ('double' ) @@ -533,25 +550,22 @@ def _spd_setitem_ ( s , i , value ) : ## list of complex types: native and X++ complex_types = complex , COMPLEX, COMPLEXf, COMPLEXl # ============================================================================= -if root_version_int < 62200 : +if root_info < ( 6 , 22 ) : def _real_ ( s ) : return s.real () def _imag_ ( s ) : return s.imag () def _cmplx_to_complex_ ( s ) : - """Convert C++ complex to Python's complex""" + """ Convert C++ complex to Python's complex""" return complex ( s.real () , s.imag () ) else : def _real_ ( s ) : return s.real def _imag_ ( s ) : return s.imag def _cmplx_to_complex_ ( s ) : - """Convert C++ complex to Python's complex""" + """ Convert C++ complex to Python's complex""" return complex ( s.real , s.imag ) -# ============================================================================= - - # ============================================================================= def _cmplx_negate_ ( s ) : - """Negation: + """ Negation: >>> v = ... >>> v1 = -v """ @@ -559,7 +573,7 @@ def _cmplx_negate_ ( s ) : # ============================================================================= def _cmplx_abs_ ( s ) : - """Absolute value + """ Absolute value >>> print abs(v) """ import math @@ -567,7 +581,7 @@ def _cmplx_abs_ ( s ) : # ============================================================================= def _cmplx_norm_ ( s ) : - """Norm (squared absolute value) + """ Norm ( squared absolute value ) >>> print v.norm() """ sr = _real_ ( s ) @@ -583,49 +597,49 @@ def _cmplx_conjugate_ ( s ) : # ============================================================================= def _cmplx_add_ ( s , o ) : - """add complex values + """ Add complex values >>> r = v + other """ return o + complex ( s ) # ============================================================================= def _cmplx_mul_ ( s , o ) : - """multiply complex values + """ Multiply complex values >>> r = v * other """ return o * complex ( s ) # ============================================================================= def _cmplx_div_ ( s , o ) : - """divide complex values + """ Divide complex values >>> r = v / other """ return ( 1.0 / o ) * complex ( s ) # ============================================================================= def _cmplx_rdiv_ ( s , o ) : - """divide complex values + """ Divide complex values >>> r = other / v """ return o * ( 1.0 / complex ( s ) ) # ============================================================================= def _cmplx_sub_ ( s , o ) : - """subtract complex values + """ Subtract complex values >>> r = v - other """ return (-o ) + complex ( s ) # ============================================================================= def _cmplx_rsub_ ( s , o ) : - """subtract complex values + """ Subtract complex values >>> r = other - v """ return o - complex ( s ) # ============================================================================= def _cmplx_pow_ ( s , o ) : - """power function + """ Power function >>> r = v ** other """ if isinstance ( o , COMPLEX ) : o = complex ( o ) @@ -633,15 +647,14 @@ def _cmplx_pow_ ( s , o ) : # ============================================================================= def _cmplx_rpow_ ( s , o ) : - """power function - >>> r = other **v + """ Power function + >>> r = other ** v """ return o ** complex ( s ) - # ============================================================================= def _cmplx_eq_ ( s , o ) : - """equality: + """ Equality: >>> r = v == other """ if isinstance ( o, COMPLEX ) : @@ -650,28 +663,27 @@ def _cmplx_eq_ ( s , o ) : # ============================================================================= def _cmplx_ne_ ( s , o ) : - """non-equality: + """ Non-equality: >>> r = v != other """ if isinstance ( o, COMPLEX ) : return _real_ ( s ) != _real_ ( o ) or _imag_ ( s ) != _imag_ ( o ) return complex ( s ) != o - # ============================================================================== -## deserialize the complex numbers +## Deserialize the complex numbers def _cmplx_factory_ ( cmplxt , re , im ) : - """Deserialize the complex numbers + """ Deserialize the complex numbers """ return cmplxt ( re , im ) # ============================================================================== ## reduce complex numbers def _cmplx_reduce_ ( c ) : - """Reduce complex numbers""" + """ Reduce complex numbers""" return _cmplx_factory_ , ( type ( c ) , c.real , c.imag ) # ============================================================================= -if root_version_int < 62200 : +if root_info < ( 6 , 22 ) : # ========================================================================= def _cmplx_iadd_ ( s , o ) : x = s + o @@ -727,7 +739,6 @@ def _cmplx_idiv_ ( s , o ) : s.__assign__ ( t ) return s # ========================================================================= - # ============================================================================= for CMPLX in ( COMPLEX , COMPLEXf , COMPLEXl ) : @@ -744,11 +755,11 @@ def _cmplx_new_init_ ( s , a = 0 , *b ) : elif 1 != len( b ) : raise TypeError("Can't create std::complex!") else : - b = b[0] + b = b [ 0 ] return s._old_init_ ( a , b ) - CMPLX.__init__ = _cmplx_new_init_ + CMPLX.__init__ = _cmplx_new_init_ CMPLX.__complex__ = _cmplx_to_complex_ @@ -835,7 +846,7 @@ def is_complex ( value ) : # @author Vanya BELYAEV Ivan.Belyaev@itep.ru # @date 2015-07-20 def frexp10 ( value ) : - """Get the mantissa (0.1<=m<1) and exponent for radix10 + """ Get the mantissa (0.1<=m<1) and exponent for radix10 (similar for frexp, but use radix=10) >>> a,b = frexp10 ( value ) @@ -850,7 +861,7 @@ def frexp10 ( value ) : q = math.floor ( math.log10 ( float ( xv ) ) ) if 0 < q : xv /= 10**q - elif 0 > q : xv *= 10**abs(q) + elif 0 > q : xv *= 10**abs ( q ) if 1 <= xv : xv /= 10 @@ -866,13 +877,16 @@ def frexp10 ( value ) : # mn, mx = num_range ( value ) # @endcode def num_range ( value , N = 1 ) : - """ Define some "range" for the given value: + """ Define some `range' for the given value: >>> value = ... >>> mn, mx = num_range ( value ) """ - + if iszero ( value ) : return ( -0.5 , 0.5 ) + assert isinstance ( N , int ) and 1 <= N , \ + "Inavalid `N' parameter:%s" % N + a , b = frexp10 ( value ) b -= N @@ -888,7 +902,6 @@ def num_range ( value , N = 1 ) : xmax = ac * ( 10 ** b ) * 0.5 return xmin , xmax - # ============================================================================= ## Find suitable range for histogram axis @@ -904,12 +917,10 @@ def axis_range ( xmin , xmax , delta = 0.02 , log = False ) : ## 2) special case if isequal ( xmn , xmx ) : return num_range ( 0.5 * ( xmn + xmx ) ) - ## 3) special case if islong ( xmn - 0.5 ) and islong ( xmx + 0.5 ) : return math.floor ( xmn - 0.1 ) , math.ceil ( xmx + 0.1 ) - d = xmx - xmn delta = abs ( delta ) @@ -952,22 +963,24 @@ def axis_range ( xmin , xmax , delta = 0.02 , log = False ) : return xmin, xmax # ============================================================================= -if (3,9) <= sys.version_info : +if ( 3 , 9 ) <= sys.version_info : # ========================================================================= ## Least Common Multiple. lcm = math.lcm # ========================================================================= ## Greatest Common Divisor - gcd = math.gcd -elif (3,5) <= sys.version_info : + gcd = math.gcd + # ========================================================================= +elif ( 3 , 5 ) <= sys.version_info : # ========================================================================= ## Greatest Common Divisor gcd = math.gcd ## Least Common Multiple. # ========================================================================= def lcm ( a , b ) : - """ Least Common Multiple""" + """ Least Common Multiple """ return abs ( a , b ) // gcd ( a , b ) + # ========================================================================= else : # ========================================================================= ## Greatest Common Divisor @@ -976,6 +989,7 @@ def lcm ( a , b ) : def lcm ( a , b ) : """ Least Common Multiple""" return ( a * b ) // gcd ( a , b ) + # ========================================================================= # ============================================================================= ## imports at the end of the module to avoid the circular dependency @@ -984,7 +998,6 @@ def lcm ( a , b ) : import ostap.math.reduce import ostap.math.polynomials - # ============================================================================= if '__main__' == __name__ : @@ -992,12 +1005,12 @@ def lcm ( a , b ) : docme ( __name__ , logger = logger ) _v = [ l for l in dir(Ostap ) if 0 != l.find('__') ] - logger.info (' dir(Ostap) : ') + logger.info ('dir(Ostap) : ') _v.sort() for v in _v : logger.info ( v ) _v = [ l for l in dir(Ostap.Math) if 0 != l.find('__') ] - logger.info (' dir(Ostap.Math) : ') + logger.info ('dir(Ostap.Math) : ') _v.sort() for v in _v : logger.info ( v ) diff --git a/ostap/math/linalg2.py b/ostap/math/linalg2.py index 732d15ae..8fcdd445 100644 --- a/ostap/math/linalg2.py +++ b/ostap/math/linalg2.py @@ -6,7 +6,7 @@ # @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl # @date 2009-09-12 # ============================================================================= -"""Few utilities to simplify linear algebra manipulations +""" Few utilities to simplify linear algebra manipulations """ # ============================================================================= from __future__ import print_function @@ -21,8 +21,14 @@ 'correlation' , ## get i,j-correlation coeffiecient from matrix-like object ) # ============================================================================= -from sys import version_info as python_version -from builtins import range +from sys import version_info as python_version +from builtins import range +from ostap.math.base import isequal , iszero, std , Ostap, typename +from ostap.core.ostap_types import num_types , integer_types +from ostap.utils.clsgetter import classgetter +from ostap.logger.pretty import pretty_array, fmt_pretty_float +import ostap.logger.table as T +from ostap.logger.colorized import infostr import ROOT, math, re, ctypes, array, random # ============================================================================= # logging @@ -31,11 +37,6 @@ if '__main__' == __name__ : logger = getLogger ( 'ostap.math.linalg2' ) else : logger = getLogger ( __name__ ) # ============================================================================= -from ostap.math.base import isequal , iszero, std , Ostap, typename -from ostap.core.ostap_types import num_types , integer_types -from ostap.utils.clsgetter import classgetter -import ostap.logger.table as T -# ============================================================================= try : import numpy as np except ImportError : @@ -991,7 +992,6 @@ def SIMT ( a , b ) : raise NotImplementedError ( "No SIMT for %s/%s and %s/%s" % ( a , typename ( a ) , b , typename ( b ) ) ) - # ========================================================================= ## power function for square matrices # @code @@ -1073,7 +1073,192 @@ def M_ASYM ( a ) : return result raise NotImplementedError ( "Cannot anti-symmetrise %s/%s" % ( a , typename ( a ) ) ) + + + # ============================================================================= + ## Get minimal element of the matrix + # @code + # mtrx = ... + # mtrx.min_element() + # @endcode + # @see Ostap::Math::min_element + @staticmethod + def M_MINELEMENT ( mtrx ) : + """ Get minimal element of the matrix + >>> mtrx = ... + >>> mtrx.min_element() + """ + return Ostap.Math.min_element ( mtrx ) + + # ============================================================================= + ## Get maximal element of the matrix + # @code + # mtrx = ... + # mtrx.max_element() + # @endcode + # @see Ostap::Math::max_element + @staticmethod + def M_MAXELEMENT ( mtrx ) : + """ Get maximal element of the matrix + >>> mtrx = ... + >>> mtrx.max_element() + """ + return Ostap.Math.max_element ( mtrx ) + + # ============================================================================= + ## Get element with the minimal absolute value of the matrix + # @code + # mtrx = ... + # mtrx.minabs_element() + # @endcode + # @see Ostap::Math::minabs_element + @staticmethod + def M_MINABSELEMENT ( mtrx ) : + """ Get element with the minimal absolute value of the matrix + >>> mtrx = ... + >>> mtrx.minabs_element() + """ + return Ostap.Math.minabs_element ( mtrx ) + + # ============================================================================= + ## Get the element with the maximal absolute value of the matrix + # @code + # mtrx = ... + # mtrx.maxabs_element() + # @endcode + # @see Ostap::Math::maxabs_element + @staticmethod + def M_MAXABSELEMENT ( mtrx ) : + """ Get element with the maximal absoluye value of the matrix + >>> mtrx = ... + >>> mtrx.maxabs_element() + """ + return Ostap.Math.maxabs_element ( mtrx ) + # ============================================================================= + ## Get index of the minimal element of the matrix + # @code + # mtrx = ... + # mtrx.min_element_index () + # @endcode + # @see Ostap::Math::ind_min_element + @staticmethod + def M_MININDEX ( mtrx ) : + """ Get index of the minimal element of the matrix + >>> mtrx = ... + >>> mtrx.min_element_index () + """ + i , j = Ostap.Math.ind_min_element ( mtrx ) + return i , j + + # ============================================================================= + ## Get idnex of maximal element of the matrix + # @code + # mtrx = ... + # mtrx.max_element_index () + # @endcode + # @see Ostap::Math::ind_max_element + @staticmethod + def M_MAXINDEX ( mtrx ) : + """ Get maximal element of the matrix + >>> mtrx = ... + >>> mtrx.max_element_index () + """ + i, j = Ostap.Math.ind_max_element ( mtrx ) + return i, j + + # ============================================================================= + ## Get index of the element with minial absolute value of the matrix + # @code + # mtrx = ... + # mtrx.minabs_element_index () + # @endcode + # @see Ostap::Math::ind_minabs_element + @staticmethod + def M_MINABSINDEX ( mtrx ) : + """ Get index of the element with minimal absoluet value of the matrix + >>> mtrx = ... + >>> mtrx.minabs_element_index () + """ + i , j = Ostap.Math.ind_minabs_element ( mtrx ) + return i , j + + # ============================================================================= + ## Get index of element with maximal absolute value of the matrix + # @code + # mtrx = ... + # mtrx.maxabs_element_index () + # @endcode + # @see Ostap::Math::ind_maxabs_element + @staticmethod + def M_MAXABSINDEX ( mtrx ) : + """ Get index of element with maximal absolute value of the matrix + >>> mtrx = ... + >>> mtrx.maxabs_element_index () + """ + i, j = Ostap.Math.ind_maxabs_element ( mtrx ) + return i, j + + # ============================================================================= + ## Get minimal diagonal element of the matrix + # @code + # mtrx = ... + # mtrx.min_diagnal() + # @endcode + # @see Ostap::Math::min_diagonal + @staticmethod + def M_MINDIAGONAL ( mtrx ) : + """ Get minimal diagonal element of the matrix + >>> mtrx = ... + >>> mtrx.min_diagonal () + """ + return Ostap.Math.min_diagonal( mtrx ) + + # ============================================================================= + ## Get maximal diagonal element of the matrix + # @code + # mtrx = ... + # mtrx.max_diagnal() + # @endcode + # @see Ostap::Math::max_diagonal + @staticmethod + def M_MAXDIAGONAL ( mtrx ) : + """ Get maximal diagonal element of the matrix + >>> mtrx = ... + >>> mtrx.max_diagonal () + """ + return Ostap.Math.max_diagonal( mtrx ) + + # ============================================================================= + ## Get diagonal element with minamal absoluite value of the matrix + # @code + # mtrx = ... + # mtrx.minabs_diagnal() + # @endcode + # @see Ostap::Math::minabs_diagonal + @staticmethod + def M_MINABSDIAGONAL ( mtrx ) : + """ Get diagonal element with mininal diagoinal value of the matrix + >>> mtrx = ... + >>> mtrx.minabs_diagonal () + """ + return Ostap.Math.minabs_diagonal( mtrx ) + + # ============================================================================= + ## Get diagonal element with maximal absolute value of the matrix + # @code + # mtrx = ... + # mtrx.maxabs_diagnal() + # @endcode + # @see Ostap::Math::maxabs_diagonal + @staticmethod + def M_MAXABSDIAGONAL ( mtrx ) : + """ Get diagonal element with maximal diagonal value of the matrix + >>> mtrx = ... + >>> mtrx.maxabs_diagonal () + """ + return Ostap.Math.maxabs_diagonal( mtrx ) + # ============================================================================= ## get matrix shape @@ -1124,7 +1309,7 @@ def V_ITER ( vct ) : # @endcode @staticmethod def V_KEYS ( vct ) : - """Iterator for SVector + """ Iterator for SVector >>> vct = ... >>> for i in vct.keys() : print i """ @@ -1140,7 +1325,7 @@ def V_KEYS ( vct ) : # @endcode @staticmethod def V_ITEMS ( vct ) : - """Iterator for SVector + """ Iterator for SVector >>> vct = ... >>> for i,v in vct.items () : print i,v >>> for i,v in vct.iteritems() : print i,v ## ditto @@ -1155,9 +1340,8 @@ def V_ITEMS ( vct ) : # @date 2009-09-12 @staticmethod def V_STR ( vct , fmt = ' %g' ) : - """Self-printout of SVectors: (...) + """ Self-printout of SVectors: (...) """ - result = '(' for i , v in enumerate ( vct ) : if 0 != i : result += ', ' @@ -1170,13 +1354,11 @@ def V_STR ( vct , fmt = ' %g' ) : # @date 2020-05-15 @staticmethod def VE_STR ( vct , fmt = '' , prefix = '' , title = '' , correlations = False ) : - """Self-printout of SVectorWithError: (...) + """ Self-printout of SVectorWithError: (...) """ - v = vct.value () - + v = vct.value () c = vct.covariance () - if correlations : c = c.correlations () fmt = '%+.3f' @@ -1357,35 +1539,39 @@ def M_ITEMS ( mtrx ) : # print matrix # @endcode @staticmethod - def M_STR ( mtrx , fmt = '' , prefix = '' , title = '' ) : - """Self-printout of matrices + def M_STR ( mtrx , fmt = '' , prefix = '' , title = '' , width = 6 , precision = 4 ) : + """ Self-printout of matrices >>> matrix = ... - >>> print matrix + >>> print (matrix) """ rows = mtrx.kRows cols = mtrx.kCols - if fmt : pass - elif cols <= 9 : fmt = '%+.4g' - else : fmt = ' %+11.4g' - - if cols <= 9 : - header = tuple ( [ '\\' ] + [ '%d' % i for i in range ( cols ) ] ) - table = [ header ] - for i in range ( rows ) : - row = [ '%d' % i ] + [ fmt % mtrx ( i , j ) for j in range ( cols ) ] - table.append ( row ) - return T.table ( table , alignment = 'r'+cols*'l' , prefix = prefix , title = title ) + mae = abs ( Ostap.Math.maxabs_element ( mtrx ) ) + fmtv , expo = fmt_pretty_float ( mae , width = width , precision = precision ) + + zeros = fmtv % ( +0.0 ) , fmtv % ( -0.0 ) + if expo : + scale = 10 ** expo + title = title + ( '[x10^%+d]' % expo ) + else : + scale = 1 - line = '' + table = [ tuple ( [ '\\' ] + [ '%d' % i for i in range ( cols ) ] ) ] + for i in range ( rows ) : - line += ' |' - for j in range ( cols ) : - line += fmt % mtrx ( i , j ) - line += ' |' - if ( rows - 1 ) != i : line += '\n' - return line - + row = [ infostr ( '%d' % i ) ] + for j in range ( cols ) : + value = mtrx ( i , j ) / scale + if iszero ( value ) : item = '0' + else : + item = fmtv % value + if item in zeros : item = '0' + row.append ( item ) + table.append ( row ) + + return T.table ( table , alignment = 'r'+cols*'l' , prefix = prefix , title = title ) + # ============================================================================= ## Self-printout of symmetric matrices # @code @@ -1393,38 +1579,39 @@ def M_STR ( mtrx , fmt = '' , prefix = '' , title = '' ) : # print matrix # @endcode @staticmethod - def MS_STR ( mtrx , fmt = '' , width = 12 , prefix = '' , title = '' ) : - """Self-printout of symmetric matrices + def MS_STR ( mtrx , fmt = '' , prefix = '' , title = '' , width = 6 , precision = 4 ) : + """ Self-printout of symmetric matrices >>> matrix = ... - >>> print matrix + >>> print(matrix) """ rows = mtrx.kRows cols = mtrx.kCols - if fmt : pass - elif cols <= 9 : fmt = '%+.4g' - else : fmt = ' %+11.4g' + mae = abs ( Ostap.Math.maxabs_element ( mtrx ) ) + fmtv , expo = fmt_pretty_float ( mae , width = width , precision = precision ) - if cols <= 9 : - header = tuple ( [ '\\' ] + [ '%d' % i for i in range ( cols ) ] ) - table = [ header ] - for i in range ( rows ) : - row = [ '%d' % i ] - for j in range ( cols ) : - if j < i : row.append ( '' ) - else : row.append ( fmt % mtrx ( i , j ) ) - table.append ( row ) - return T.table ( table , alignment = 'r'+cols*'l' , prefix = prefix , title = title ) + zeros = fmtv % ( +0.0 ) , fmtv % ( -0.0 ) + if expo : + scale = 10 ** expo + title = title + ( '[x10^%+d]' % expo ) + else : + scale = 1 - line = '' + table = [ tuple ( [ '\\' ] + [ '%d' % i for i in range ( cols ) ] ) ] for i in range ( rows ) : - line += ' |' - for j in range ( cols ) : - if j < i : line += width*' ' - else : line += fmt % mtrx ( i , j ) - line += ' |' - if ( rows - 1 ) != i : line += '\n' - return line + row = [ infostr ( '%d' % i ) ] + for j in range ( cols ) : + if j < i : row.append ( '' ) + else : + value = mtrx ( i , j ) / scale + if iszero ( value ) : item = '0' + else : + item = fmtv % value + if item in zeros : item = '0' + row.append ( item ) + table.append ( row ) + + return T.table ( table , alignment = 'r'+cols*'l' , prefix = prefix , title = title ) # ========================================================================= ## get the correlation matrix @@ -1802,7 +1989,7 @@ def deco_vector ( t ) : ## Decorate SMatrix @staticmethod def deco_matrix ( m ) : - """Decorate SMatrix + """ Decorate SMatrix """ if m in LinAlg.decorated_matrices : return m @@ -1857,6 +2044,21 @@ def deco_matrix ( m ) : m.tmatrix = LinAlg.M_TM + m.min_element = LinAlg.M_MINELEMENT + m.max_element = LinAlg.M_MAXELEMENT + m.minabs_element = LinAlg.M_MINABSELEMENT + m.maxabs_element = LinAlg.M_MAXABSELEMENT + + m.min_diagonal = LinAlg.M_MINDIAGONAL + m.max_diagonal = LinAlg.M_MAXDIAGONAL + m.minabs_diagonal = LinAlg.M_MINABSDIAGONAL + m.maxabs_diagonal = LinAlg.M_MAXABSDIAGONAL + + m.min_element_index = LinAlg.M_MININDEX + m.max_element_index = LinAlg.M_MAXINDEX + m.minabs_element_index = LinAlg.M_MINABSINDEX + m.maxabs_element_index = LinAlg.M_MAXABSINDEX + if m.kRows == m.kCols : m.inverse = LinAlg.M_INVERSE @@ -1871,8 +2073,7 @@ def deco_matrix ( m ) : m.sym = LinAlg.M_SYM m.asym = LinAlg.M_ASYM m.skew = LinAlg.M_ASYM - - + m.__reduce__ = LinAlg.M_REDUCE m.rep_size = classgetter ( lambda cls : cls.rep_type.kSize ) @@ -1948,7 +2149,8 @@ def deco_vectore ( t ) : t.__reduce__ = LinAlg.VE_REDUCE t.random = LinAlg.VE_RANDOM - + + t.pretty_print = pretty_array return t @@ -2379,8 +2581,6 @@ def checkops ( a , b , logger = logger ) : ( 'ssym' , LinAlg.methods_ASYM ) , ) - - import ostap.logger.table as T title = 'Allowed binary operations' table = T.table ( rows , title = title , prefix = '# ' , alignment = 'lcccl' ) diff --git a/ostap/math/tests/test_math_base.py b/ostap/math/tests/test_math_base.py new file mode 100644 index 00000000..96c149a8 --- /dev/null +++ b/ostap/math/tests/test_math_base.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# ============================================================================= +# Copyright (c) Ostap developpers. +# ============================================================================= +## @file ostap/math/tests/test_math_base.py +# Test module for the file ostap/math/base.py +# ============================================================================= +""" Test module for ostap/math/base.py +""" +# ============================================================================= +import random +import ostap.math.base as MB +import ostap.logger.table as T +from ostap.logger.pretty import pretty_float +from ostap.logger.colorized import allright, attention +# ============================================================================= +# logging +# ============================================================================= +from ostap.logger.logger import getLogger +if '__main__' == __name__ : logger = getLogger ( 'test_math_base' ) +else : logger = getLogger ( __name__ ) +# ============================================================================= + +# ============================================================================= +## test frexp10 +def test_frexp10 () : + + logger = getLogger ( 'test_frexp10' ) + + rows = [ ( 'Value' , 'exp' , 'm/py' , 'exp/py' , 'm/C++' , 'exp/C++' , 'OK?' ) ] + + ok = True + for i in range ( 100 ) : + + mexp = random.randrange ( -100 , 101 ) + val = random.uniform ( -1 , 1 ) + value = val * 10**mexp + m1, e1 = MB.frexp10 ( value ) + m2, e2 = MB.cpp_frexp10 ( value ) + + v, n = pretty_float ( value ) + + if MB.isequal ( m1 , m2 ) and MB.isequal ( e1 , e2 ) : + good = '' + else : + ok = False + good = attention ( 'No!' ) + + row = v , '%+3d' % n if n else '' , \ + '%+.4f' % ( 10 * m1 ) , '%+3d' % ( e1 - 1 ) if e1 != 1 else '' , \ + '%+.4f' % ( 10 * m2 ) , '%+3d' % ( e2 - 1 ) if e2 != 1 else '' , good + + rows.append ( row ) + + title = 'Test for frexp10' + table = T.table ( rows , title = title , prefix = '# ' , alignment = 'lrlrlrc' ) + if ok : logger.info ( '%s:\n%s' % ( title , table ) ) + else : logger.error ( '%s:\n%s' % ( title , table ) ) + +# ============================================================================= +if '__main__' == __name__ : + + test_frexp10() + +# ============================================================================= +## The END +# ============================================================================= diff --git a/source/include/Ostap/MatrixUtils.h b/source/include/Ostap/MatrixUtils.h index 73a2d3a8..9b504332 100755 --- a/source/include/Ostap/MatrixUtils.h +++ b/source/include/Ostap/MatrixUtils.h @@ -190,7 +190,7 @@ namespace Ostap * * @param m (input/output) matrix to be modified * @param value (input) new value for all matrix elements - * @return number of modified matrix elemenets + * @return number of modified matrix elements * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2006-05-24 */ @@ -221,15 +221,16 @@ namespace Ostap * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2006-05-24 */ - template + template inline std::size_t setToUnit - ( ROOT::Math::SMatrix& m , const T& value = T(1) ) + ( ROOT::Math::SMatrix& m , const T& value = T ( 1 ) ) { - /// nullify the matrix: - std::fill ( m.begin() , m.end() , T(0.0) ) ; - /// set diagonal elements - for ( unsigned int i = 0 ; i < D ; ++i ) { m(i,i) = value ; } + /// nullify the whole matrix: + std::fill ( m.begin() , m.end() , T ( 0.0 ) ) ; + /// set diagonal elements + const unsigned int D = std::min ( D1 , D2 ) ; + for ( unsigned int i = 0 ; i < D ; ++i ) { m ( i , i ) = value ; } return m.end() - m.begin() ; } // ======================================================================== @@ -693,12 +694,13 @@ namespace Ostap * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2006-05-24 */ - template + template inline T trace - ( const ROOT::Math::SMatrix& m ) + ( const ROOT::Math::SMatrix& m ) { T result = m ( 0 , 0 ) ; + const unsigned int D = std::min ( D1 , D2 ) ; for ( unsigned int i = 1 ; i < D ; ++i ) { result += m ( i , i ) ; } return result ; } @@ -709,12 +711,13 @@ namespace Ostap * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2006-05-24 */ - template + template inline T trace - ( const ROOT::Math::Expr& m ) + ( const ROOT::Math::Expr& m ) { T result = m ( 0 , 0 ) ; + const unsigned int D = std::min ( D1 , D2 ) ; for ( unsigned int i = 1 ; i < D ; ++i ) { result += m ( i , i ) ; } return result ; } @@ -724,36 +727,38 @@ namespace Ostap * @param cmp comparison criteria * @return "min" diagonal element (in the sense of comparison criteria) */ - template + template inline T min_diagonal - ( const ROOT::Math::SMatrix& m , CMP cmp ) + ( const ROOT::Math::SMatrix& m , CMP cmp ) { T result = m(0,0); + const unsigned int D = std::min ( D1 , D2 ) ; for ( unsigned int i = 1 ; i < D ; ++i ) - { - const T value = m(i,i) ; - if ( cmp ( value , result ) ) { result = value ; } - } + { + const T value = m ( i , i ) ; + if ( cmp ( value , result ) ) { result = value ; } + } return result ; } // ======================================================================== /** find the maximal diagonal element - * @param m (input) square matrix to be studied + * @param m (input) matrix to be studied * @param cmp comparison criteria * @return "max" diagonal element (in the sense of comparison criteria) */ - template + template inline T max_diagonal - ( const ROOT::Math::SMatrix& m , CMP cmp ) + ( const ROOT::Math::SMatrix& m , CMP cmp ) { T result = m(0,0); + const unsigned int D = std::min ( D1 , D2 ) ; for ( unsigned int i = 1 ; i < D ; ++i ) - { - const T value = m(i,i) ; - if ( cmp ( result , value ) ) { result = value ; } - } + { + const T value = m ( i , i ) ; + if ( cmp ( result , value ) ) { result = value ; } + } return result ; } // ======================================================================== @@ -777,10 +782,10 @@ namespace Ostap * @author Vanya BELYAEV ibelyaev@physics.cyr.edu * @date 2006-05-24 */ - template + template inline T max_diagonal - ( const ROOT::Math::SMatrix& m ) + ( const ROOT::Math::SMatrix& m ) { return max_diagonal( m , std::less() ) ; } // ======================================================================== /** find the maximal diagonal element of the square matrix @@ -807,10 +812,10 @@ namespace Ostap * @author Vanya BELYAEV ibelyaev@physics.cyr.edu * @date 2006-05-24 */ - template + template inline T min_diagonal - ( const ROOT::Math::SMatrix& m ) + ( const ROOT::Math::SMatrix& m ) { return min_diagonal( m , std::less() ) ; } // ======================================================================== /** find the diagonal element of square matrix with maximal absolute value @@ -833,10 +838,10 @@ namespace Ostap * @author Vanya BELYAEV ibelyaev@physics.cyr.edu * @date 2006-05-24 */ - template + template inline T maxabs_diagonal - ( const ROOT::Math::SMatrix& m ) + ( const ROOT::Math::SMatrix& m ) { return max_diagonal( m , _AbsCompare() ) ; } // ======================================================================== /** find the diagonal element of the square matrix with @@ -863,12 +868,12 @@ namespace Ostap * @author Vanya BELYAEV ibelyaev@physics.cyr.edu * @date 2006-05-24 */ - template + template inline T minabs_diagonal - ( const ROOT::Math::SMatrix& m ) + ( const ROOT::Math::SMatrix& m ) { return min_diagonal( m , _AbsCompare() ) ; } - // ======================================================================== + // ======================================================================== /** count the number of elements in matrix, which satisfy the certain criteria * * @code @@ -946,13 +951,13 @@ namespace Ostap { std::size_t result = 0 ; for ( unsigned int i = 0 ; i < D ; ++i ) - { - if ( pred ( m ( i , i ) ) ) { result += 1 ; } - for ( unsigned int j = i + 1 ; j < D ; ++j ) { - if ( pred ( m ( i , j ) ) ) { result +=2 ; } // ATTENTION! + if ( pred ( m ( i , i ) ) ) { result += 1 ; } + for ( unsigned int j = i + 1 ; j < D ; ++j ) + { + if ( pred ( m ( i , j ) ) ) { result += 2 ; } // ATTENTION! + } } - } return result ; } // ======================================================================== @@ -975,18 +980,19 @@ namespace Ostap * * @endcode * - * @param m (input) square matrix to be studied + * @param m (input) matrix to be studied * @param pred (input) predicate to be tested * @return number of diagonal elements for which the predicate is valid * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2006-04-24 */ - template + template inline std::size_t count_diagonal - ( const ROOT::Math::SMatrix& m , P pred ) + ( const ROOT::Math::SMatrix& m , P pred ) { std::size_t result = 0 ; + const unsigned int D = std::min ( D1 , D2 ) ; for ( unsigned int i = 0 ; i < D ; ++i ) { if ( pred ( m ( i , i ) ) ) { result += 1 ; } } return result ; @@ -1052,13 +1058,14 @@ namespace Ostap * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2006-04-24 */ - template + template inline bool check_diagonal - ( const ROOT::Math::SMatrix& m , P pred ) - { + ( const ROOT::Math::SMatrix& m , P pred ) + { + const unsigned int D = std::min ( D1 , D2 ) ; for ( unsigned int i = 0 ; i < D ; ++i ) - { if ( pred ( m ( i , i ) ) ) { return true ; } } + { if ( pred ( m ( i , i ) ) ) { return true ; } } return false ; } // ======================================================================== @@ -1176,10 +1183,11 @@ namespace Ostap const double scale ) { for ( unsigned int i = 0 ; i < D ; ++i ) - { - for ( unsigned int j = i ; j < D ; ++j ) - { left ( i , j ) += scale * vect(i) * vect(j) ; } - } + { + const double vi = vect ( i ) ; + for ( unsigned int j = i ; j < D ; ++j ) + { left ( i , j ) += scale * vi * vect(j) ; } + } } // ========================================================================= /** update the symmetric matrix according to the rule m += s*v*v^T @@ -1197,10 +1205,11 @@ namespace Ostap const double scale ) { for ( unsigned int i = 0 ; i < D ; ++i ) - { - for ( unsigned int j = i ; j < D ; ++j ) - { left ( i , j ) += scale * vect(i) * vect(j) ; } - } + { + const double vi = vect ( i ) ; + for ( unsigned int j = i ; j < D ; ++j ) + { left ( i , j ) += scale * vi * vect(j) ; } + } } // ========================================================================= /** update the matrix according to the rule m += s*v1*v2^T @@ -1220,10 +1229,11 @@ namespace Ostap const double scale = 1.0 ) { for ( unsigned int i = 0 ; i < D1 ; ++i ) - { - for ( unsigned int j = 0 ; j < D2 ; ++j ) - { left ( i , j ) += scale * vct1(i) * vct2(j) ; } - } + { + const double vi = vct1 ( i ) ; + for ( unsigned int j = 0 ; j < D2 ; ++j ) + { left ( i , j ) += scale * vi * vct2(j) ; } + } } // ========================================================================= /** useful shortcut for product of vector, matrix and vector (v1^T*M*v2) @@ -1258,10 +1268,10 @@ namespace Ostap const double scale = 1.0 ) { for ( unsigned int i = 0 ; i < D ; ++i ) - { - for ( unsigned int j = i ; j < D ; ++j ) - { left ( i , j ) += scale * ( right ( i , j ) + right ( j , i ) ) ; } - } + { + for ( unsigned int j = i ; j < D ; ++j ) + { left ( i , j ) += scale * ( right ( i , j ) + right ( j , i ) ) ; } + } } // ========================================================================= /** update the symmetric matrix according to the rule m += scale * ( m + m^T ) @@ -1279,10 +1289,10 @@ namespace Ostap const double scale = 1.0 ) { for ( unsigned int i = 0 ; i < D ; ++i ) - { - for ( unsigned int j = i ; j < D ; ++j ) - { left ( i , j ) += scale * ( right ( i , j ) + right ( j , i ) ) ; } - } + { + for ( unsigned int j = i ; j < D ; ++j ) + { left ( i , j ) += scale * ( right ( i , j ) + right ( j , i ) ) ; } + } } // ======================================================================== /** inversion of symmetric positively defined matrices diff --git a/source/include/Ostap/MatrixUtilsT.h b/source/include/Ostap/MatrixUtilsT.h index a7b10f0d..d2310b13 100644 --- a/source/include/Ostap/MatrixUtilsT.h +++ b/source/include/Ostap/MatrixUtilsT.h @@ -317,7 +317,45 @@ namespace Ostap // ====================================================================== } ; // ======================================================================== - + template + inline T + maxabs_element + ( const TMatrixT& m ) + { + if ( !m.IsValid() ) { return 0 ; } + T result = m ( 0 , 0 ) ; + const unsigned int rows = m.GetNrows () ; + const unsigned int cols = m.GetNrows () ; + for ( unsigned int i = 0 ; i < rows ; ++i ) + { + for ( unsigned int j = 0 ; j < cols ; ++j ) + { + const double value = m ( i , j ) ; + if ( std::abs ( result ) < std::abs ( value ) ) { result = value ; } + } + } + return result ; + } + // ======================================================================== + template + inline T + maxabs_element + ( const TMatrixTSym& m ) + { + if ( !m.IsValid() ) { return 0 ; } + T result = m ( 0 , 0 ) ; + const unsigned int rows = m.GetNrows () ; + const unsigned int cols = m.GetNrows () ; + for ( unsigned int i = 0 ; i < rows ; ++i ) + { + for ( unsigned int j = i ; j < cols ; ++j ) + { + const double value = m ( i , j ) ; + if ( std::abs ( result ) < std::abs ( value ) ) { result = value ; } + } + } + return result ; + } // ======================================================================== namespace Ops {