From 271e81372d400bb329a20d0384eac92813acbe3c Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Fri, 17 Nov 2023 12:17:16 +0100 Subject: [PATCH] 1. Further optimisation in `Ostap::Math::ChebyshedSum` 1. add new test `ostap/math/tests/test_math.poly.py` --- ReleaseNotes/release_notes.md | 5 +- ostap/math/tests/test_math_poly.py | 184 ++++++++++++++++++++++++++ ostap/math/tests/test_math_poly_nD.py | 2 +- source/src/Polynomials.cpp | 72 +++++----- 4 files changed, 230 insertions(+), 33 deletions(-) create mode 100644 ostap/math/tests/test_math_poly.py diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index 12e1d372..d18f3e14 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -2,7 +2,8 @@ 1. Update `histo_compare` tests 1. Slight optimisation in `Ostap::Math::ChebyshedSum` - + 1. Further optimisation in `Ostap::Math::ChebyshedSum` + 1. add new test `ostap/math/tests/test_math.poly.py` ## Backward incompatible: @@ -23,7 +24,7 @@ ## Bug fixes: - 1. Fix a sad bug in `Ostap::Math::Bernstein` for incorrectuage of `elevate` + 1. Fix a sad bug in `Ostap::Math::Bernstein` for incorrect usage of `elevate` # v1.10.0.0 diff --git a/ostap/math/tests/test_math_poly.py b/ostap/math/tests/test_math_poly.py new file mode 100644 index 00000000..af3c0287 --- /dev/null +++ b/ostap/math/tests/test_math_poly.py @@ -0,0 +1,184 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +# ============================================================================= +# Copyright (c) Ostap developpers. +# ============================================================================= +## @file ostap/math/tests/test_math_poly.py +# Test script for polynomial functions +# @see Ostap::Math::LegendreSum +# @see Ostap::Math::HermiteSum +# @see Ostap::Math::ChebyshevSum +# @see Ostap::Math::Bernstein +# @see Ostap::Math::Polynomial +# ============================================================================= +""" Test script for 2&3-dimnsional polynomial functions +- see Ostap::Math::LegendreSum +- see Ostap::Math::HermiteSum +- see Ostap::Math::ChebyshevSum +- see Ostap::Math::Bernstein +- see Ostap::Math::Polynomial +""" +# ============================================================================= +from ostap.core.core import Ostap, SE +from ostap.logger.colorized import attention +from ostap.utils.progress_bar import progress_bar +import ostap.math.models +import ostap.math.integral as I +import ostap.math.derivative as D +import ostap.logger.table as T +import ROOT, random, math +# ============================================================================= +# logging +# ============================================================================= +from ostap.logger.logger import getLogger +if '__main__' == __name__ : logger = getLogger ( 'ostap.test_math_poly' ) +else : logger = getLogger ( __name__ ) +# ============================================================================= + +xmin = 0 +xmax = 10 +N = 5 +scale = 1.e+12 + +# ============================================================================= +## test derivative +def test_poly_derivative () : + """Test derivative + """ + logger = getLogger ( 'test_poly_derivative') + + poly = ( + Ostap.Math.ChebyshevSum ( N , xmin , xmax ) , + Ostap.Math.LegendreSum ( N , xmin , xmax ) , + Ostap.Math.HermiteSum ( N , xmin , xmax ) , + Ostap.Math.Bernstein ( N , xmin , xmax ) , + Ostap.Math.Polynomial ( N , xmin , xmax ) , + ) + + + for p in poly : + n = p.npars() + for i in range ( n ) : p [ i ] = random.uniform ( -10 , 10 ) + + ders = tuple ( [ p.derivative() for p in poly ] ) + + cnts = [ SE() for i in poly ] + + dx = 0.1 * ( xmax - xmin ) + for i in progress_bar ( range ( 1000 ) ) : + + x = random.uniform ( xmin + dx , xmax - dx ) + + for i, p in enumerate ( poly ) : + d1 = p.derivative ( x ) + d2 = ders [ i ] ( x ) + d3 = D.derivative ( p , x ) + + dd = max ( abs ( d1 - d2 ) / ( abs ( d1 ) + abs ( d2 ) ) , + abs ( d1 - d3 ) / ( abs ( d1 ) + abs ( d3 ) ) , + abs ( d2 - d3 ) / ( abs ( d2 ) + abs ( d3 ) ) ) + + ## dd = abs ( d2 - d3 ) / ( abs ( d2 ) + abs ( d3 ) ) + + + cnts [ i ] += 2 * dd * scale + + + names = 'Chebyshev' , 'Legendre' , 'Hermite' , 'Bernstein' , 'Polynomial' + + rows = [ ( 'Poly' , + 'mean [%.0g]' % ( 1/scale ) , + 'rms [%.0g]' % ( 1/scale ) , + 'max [%.0g]' % ( 1/scale ) ) ] + + for name, cnt in zip ( names , cnts ) : + + row = name , \ + '%.3g' % cnt.mean () , \ + '%.3g' % cnt.rms () , \ + '%.3g' % cnt.max () + + rows.append ( row ) + + title = 'Derivatives (analytical&numerical)' + table = T.table ( rows , title = title , prefix = '# ' , alignment = 'lccc' ) + logger.info ( '%s\n%s' % ( title , table ) ) + + +# ============================================================================= +## test integrals +def test_poly_integral () : + """Test derivative + """ + logger = getLogger ( 'test_poly_integral') + + poly = ( + Ostap.Math.ChebyshevSum ( N , xmin , xmax ) , + Ostap.Math.LegendreSum ( N , xmin , xmax ) , + Ostap.Math.HermiteSum ( N , xmin , xmax ) , + Ostap.Math.Bernstein ( N , xmin , xmax ) , + Ostap.Math.Polynomial ( N , xmin , xmax ) , + ) + + + for p in poly : + n = p.npars() + for i in range ( n ) : p [ i ] = random.uniform ( -10 , 10 ) + + ints = tuple ( [ p.indefinite_integral() for p in poly ] ) + + cnts = [ SE() for i in poly ] + + dx = 0.1 * ( xmax - xmin ) + for i in progress_bar ( range ( 1000 ) ) : + + x = random.uniform ( xmin , xmax ) + y = random.uniform ( xmin , xmax ) + + for i, p in enumerate ( poly ) : + d1 = p.integral ( x , y ) + d2 = ints [ i ] ( y ) - ints [ i ] ( x ) + d3 = I.integral ( p , x , y ) + + dd = max ( abs ( d1 - d2 ) / ( abs ( d1 ) + abs ( d2 ) ) , + abs ( d1 - d3 ) / ( abs ( d1 ) + abs ( d3 ) ) , + abs ( d2 - d3 ) / ( abs ( d2 ) + abs ( d3 ) ) ) + + ## dd = abs ( d2 - d3 ) / ( abs ( d2 ) + abs ( d3 ) ) + + + cnts [ i ] += 2 * dd * scale + + + names = 'Chebyshev' , 'Legendre' , 'Hermite' , 'Bernstein' , 'Polynomial' + + rows = [ ( 'Poly' , + 'mean [%.0g]' % ( 1/scale ) , + 'rms [%.0g]' % ( 1/scale ) , + 'max [%.0g]' % ( 1/scale ) ) ] + + for name, cnt in zip ( names , cnts ) : + + row = name , \ + '%.3g' % cnt.mean () , \ + '%.3g' % cnt.rms () , \ + '%.3g' % cnt.max () + + rows.append ( row ) + + title = 'Integrals (analytical&numerical)' + table = T.table ( rows , title = title , prefix = '# ' , alignment = 'lccc' ) + logger.info ( '%s\n%s' % ( title , table ) ) + + +# ============================================================================= + +if '__main__' == __name__ : + + test_poly_derivative () + test_poly_integral () + + +# ============================================================================= +## The END +# ============================================================================= diff --git a/ostap/math/tests/test_math_poly_nD.py b/ostap/math/tests/test_math_poly_nD.py index 8388eedc..d71ab730 100644 --- a/ostap/math/tests/test_math_poly_nD.py +++ b/ostap/math/tests/test_math_poly_nD.py @@ -11,7 +11,7 @@ # @see Ostap::Math::Bernstein3D # ============================================================================= """ Test script for 2&3-dimnsional polynomial functions -a- see Ostap::Math::LegendreSum2 +- see Ostap::Math::LegendreSum2 - see Ostap::Math::LegendreSum3 - see Ostap::Math::Bernstein2D - see Ostap::Math::Bernstein3D diff --git a/source/src/Polynomials.cpp b/source/src/Polynomials.cpp index 5e12ad13..c3ff1e29 100644 --- a/source/src/Polynomials.cpp +++ b/source/src/Polynomials.cpp @@ -1056,16 +1056,18 @@ Ostap::Math::Polynomial Ostap::Math::Polynomial::indefinite_integral ( const double C ) const { const double dx = 0.5 * ( m_xmax - m_xmin ) ; - std::vector npars ( m_pars.size() + 1 , 0 ) ; + // + Ostap::Math::Polynomial integ ( degree() + 1 , m_xmin , m_xmax ) ; + // for ( unsigned int i = 0 ; i < m_pars.size() ; ++i ) { if ( s_zero ( m_pars[i] ) ) { continue ; } // CONTINUE - npars [i+1] = m_pars[i] / ( i + 1 ) * dx ; + integ.m_pars [ i + 1 ] = m_pars [ i ] / ( i + 1 ) * dx ; } // - npars[0] = C ; + integ.m_pars[0] = C ; // - return Ostap::Math::Polynomial( npars , m_xmin , m_xmax ) ; + return integ ; } // ============================================================================ // get the derivative at point "x" @@ -1097,15 +1099,16 @@ Ostap::Math::Polynomial::derivative () const { return Ostap::Math::Polynomial( 0 , m_xmin , m_xmax ) ; } // const double dx = 2 / ( m_xmax - m_xmin ) ; - std::vector npars ( m_pars.size() - 1 , 0 ) ; - for ( unsigned int i = 0 ; i < npars.size() ; ++i ) + // + Ostap::Math::Polynomial deriv ( degree() - 1 , m_xmin , m_xmax ) ; + for ( unsigned int i = 0 ; i < deriv.npars () ; ++i ) { - const double p = m_pars[i+1] ; + const double p = m_pars [ i + 1 ] ; if ( s_zero ( p ) ) { continue ; } // CONTINUE - npars [i] = ( i + 1 ) * p * dx ; + deriv.m_pars [ i ] = ( i + 1 ) * p * dx ; } // - return Ostap::Math::Polynomial ( npars , m_xmin , m_xmax ) ; + return deriv ; } // ============================================================================ // simple manipulations with polynoms: shift it! @@ -1371,22 +1374,25 @@ Ostap::Math::ChebyshevSum::indefinite_integral ( const double C ) const { // const double dx = 0.5 * ( m_xmax - m_xmin ) ; - std::vector npars ( m_pars.size() + 1 , 0 ) ; + // + Ostap::Math::ChebyshevSum integ ( degree() + 1 , m_xmin , m_xmax ) ; + // for ( unsigned short i = 0 ; i < m_pars.size() ; ++i ) { if ( s_zero ( m_pars[i] ) ) { continue ; } // CONTINUE // - if ( 0 == i ) { npars [1] += m_pars[0] * dx ; } - else if ( 1 == i ) { npars [2] += 0.25 * m_pars[1] * dx ; } + if ( 0 == i ) { integ.m_pars [ 1 ] += m_pars [ 0 ] * dx ; } + else if ( 1 == i ) { integ.m_pars [ 2 ] += 0.25 * m_pars [ 1 ] * dx ; } else { - npars [ i + 1 ] += m_pars[i] * 0.5 / ( i + 1.0 ) * dx ; - npars [ i - 1 ] -= m_pars[i] * 0.5 / ( i - 1.0 ) * dx ; + integ.m_pars [ i + 1 ] += m_pars [ i ] * 0.5 / ( i + 1.0 ) * dx ; + integ.m_pars [ i - 1 ] -= m_pars [ i ] * 0.5 / ( i - 1.0 ) * dx ; } } // - npars[0] += C ; - return Ostap::Math::ChebyshevSum ( npars , m_xmin , m_xmax ) ; + integ.m_pars [ 0 ] += C ; + // + return integ ; } // ============================================================================ // get the derivative at point "x" @@ -1741,7 +1747,7 @@ double Ostap::Math::LegendreSum::evaluate ( const double x ) const // get the integral between xmin and xmax // ============================================================================ double Ostap::Math::LegendreSum::integral () const -{ return m_pars[0] * ( m_xmax - m_xmin ) ; } +{ return m_pars [ 0 ] * ( m_xmax - m_xmin ) ; } // ============================================================================ // get the integral between low and high // ============================================================================ @@ -1774,14 +1780,16 @@ double Ostap::Math::LegendreSum::integral // } // std::vector npars ( m_pars.size() + 1 , 0 ) ; - for ( unsigned int i = 0 ; i < m_pars.size() ; ++i ) + for ( unsigned int i = 1 ; i < m_pars.size() ; ++i ) { if ( s_zero ( m_pars[i] ) ) { continue ; } // CONTINUE // - npars [i+1] += m_pars[i] / ( 2*i + 1 ) ; - if ( 0 < i ) { npars [i-1] -= m_pars[i] / ( 2*i + 1 ) ; } + npars [ i + 1 ] += m_pars[i] / ( 2*i + 1 ) ; + npars [ i - 1 ] -= m_pars[i] / ( 2*i + 1 ) ; } // + npars [ 1 ] += m_pars [ 0 ] ; + // const double result = Ostap::Math::Clenshaw::legendre_sum ( npars.begin() , npars.end() , xh ) - Ostap::Math::Clenshaw::legendre_sum ( npars.begin() , npars.end() , xl ) ; @@ -1797,17 +1805,20 @@ Ostap::Math::LegendreSum::indefinite_integral ( const double C ) const { // const double dx = 0.5 * ( m_xmax - m_xmin ) ; - std::vector npars ( m_pars.size() + 1 , 0 ) ; - for ( unsigned int i = 0 ; i < m_pars.size() ; ++i ) + // + Ostap::Math::LegendreSum integ ( degree() + 1 , m_xmin , m_xmax ) ; + for ( unsigned int i = 1 ; i < m_pars.size() ; ++i ) { - if ( s_zero ( m_pars[i] ) ) { continue ; } // CONTINUE + if ( s_zero ( m_pars [ i ] ) ) { continue ; } // CONTINUE // - npars [i+1] += m_pars[i] / ( 2*i + 1 ) * dx ; - if ( 0 < i ) { npars [i-1] -= m_pars[i] / ( 2*i + 1 ) * dx ; } + integ.m_pars [ i + 1 ] += m_pars [ i ] / ( 2*i + 1 ) * dx ; + integ.m_pars [ i - 1 ] -= m_pars [ i ] / ( 2*i + 1 ) * dx ; } // - npars[0] += C ; - return Ostap::Math::LegendreSum ( npars , m_xmin , m_xmax ) ; + integ.m_pars [ 1 ] += m_pars [ 0 ] * dx ; + integ.m_pars [ 0 ] += C ; + // + return integ ; } // ============================================================================ // get the derivative at point "x" @@ -1840,16 +1851,17 @@ Ostap::Math::LegendreSum::derivative () const { return Ostap::Math::LegendreSum ( 0 , m_xmin , m_xmax ) ; } // const double dx = 2 / ( m_xmax - m_xmin ) ; - std::vector npars ( m_pars.size() - 1 , 0 ) ; + // + Ostap::Math::LegendreSum deriv ( degree() - 1 , m_xmin , m_xmax ) ; for ( unsigned short i = 1 ; i < m_pars.size() ; ++i ) { if ( s_zero ( m_pars[i] ) ) { continue ; } // CONTINUE // for ( int j = i-1 ; 0<=j ; j-=2 ) - { npars[j] += m_pars[i] * ( 2 * j + 1 ) * dx ; } + { deriv.m_pars [ j ] += m_pars[i] * ( 2 * j + 1 ) * dx ; } } // - return Ostap::Math::LegendreSum ( npars , m_xmin , m_xmax ) ; + return deriv ; } // ============================================================================ // simple manipulations with polynoms: shift it!