Skip to content

Commit

Permalink
1. Further optimisation in Ostap::Math::ChebyshedSum
Browse files Browse the repository at this point in the history
  1. add new test `ostap/math/tests/test_math.poly.py`
  • Loading branch information
VanyaBelyaev committed Nov 17, 2023
1 parent 6518046 commit 271e813
Show file tree
Hide file tree
Showing 4 changed files with 230 additions and 33 deletions.
5 changes: 3 additions & 2 deletions ReleaseNotes/release_notes.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:

Expand All @@ -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

Expand Down
184 changes: 184 additions & 0 deletions ostap/math/tests/test_math_poly.py
Original file line number Diff line number Diff line change
@@ -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
# =============================================================================
2 changes: 1 addition & 1 deletion ostap/math/tests/test_math_poly_nD.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
72 changes: 42 additions & 30 deletions source/src/Polynomials.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<double> 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"
Expand Down Expand Up @@ -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<double> 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!
Expand Down Expand Up @@ -1371,22 +1374,25 @@ Ostap::Math::ChebyshevSum::indefinite_integral ( const double C ) const
{
//
const double dx = 0.5 * ( m_xmax - m_xmin ) ;
std::vector<double> 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"
Expand Down Expand Up @@ -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
// ============================================================================
Expand Down Expand Up @@ -1774,14 +1780,16 @@ double Ostap::Math::LegendreSum::integral
// }
//
std::vector<double> 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 ) ;
Expand All @@ -1797,17 +1805,20 @@ Ostap::Math::LegendreSum::indefinite_integral ( const double C ) const
{
//
const double dx = 0.5 * ( m_xmax - m_xmin ) ;
std::vector<double> 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"
Expand Down Expand Up @@ -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<double> 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!
Expand Down

0 comments on commit 271e813

Please sign in to comment.