diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index d312d9b8..16c0a523 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -12,6 +12,7 @@ 1. Add `Ostap::Math::RationalPositive` : rational function as ratio of two positve Bernstein polynomials 1. Add `Ostap::Models::Rational` : rational PDF as ratio of two positve Bernstein polynomials 1. Add `Rational_pdf` : rational PDF as ratio of two positve Bernstein polynomials + 1. Add 1D-histogram parameterisations in terms of rational functions: `rational_fun`, `rational` and `brational` ## Backward incompatible: diff --git a/ostap/histos/param.py b/ostap/histos/param.py index 674538db..8253849c 100755 --- a/ostap/histos/param.py +++ b/ostap/histos/param.py @@ -13,16 +13,16 @@ # - as Fourier sum (if numpy/scipy available) # - as Cosine/Fourier sum (if numpy/scipy available) # - as plain monomal sum -# - as positive Bernstein/Bezier sum -# - as positive even Bernstein/Bezier sum +# - as positive Bernstein/Bezier sum +# - as positive even Bernstein/Bezier sum # - as positive monotonic Bernstein/Bezier sum -# - as positive convex/concave Bernstein/Bezier sum +# - as positive convex/concave Bernstein/Bezier sum # - as positive monotonic convex/concave Bernstein/Bezier sum -# - as generic spline (b-spline) -# - as positive spline (p-spline) +# - as generic spline (b-spline) +# - as positive spline (p-spline) # - as positive monotonic spline (m-spline) # - as positive monotonic convex/concave spline (c-spline) -# - as positive convex/concave spline (convex/concave-spline) +# - as positive convex/concave spline (convex/concave-spline) # # where possible, fit starts from reasonable approximatuion, taken from (1) # @@ -43,15 +43,15 @@ # ## (2) using RooFit pdfs: # -# - as positive Bernstein/Bezier sum -# - as positive even Bernstein/Bezier sum +# - as positive Bernstein/Bezier sum +# - as positive even Bernstein/Bezier sum # - as positive monotonic Bernstein/Bezier sum -# - as positive convex/concave Bernstein/Bezier sum +# - as positive convex/concave Bernstein/Bezier sum # - as positive monotonic convex/concave Bernstein/Bezier sum -# - as positive spline (p-spline) +# - as positive spline (p-spline) # - as positive monotonic spline (m-spline) # - as positive monotonic convex/concave spline (c-spline) -# - as positive convex/concave spline (convex/concave-spline) +# - as positive convex/concave spline (convex/concave-spline) # # Typical usage: # @code @@ -79,14 +79,14 @@ - as plain monomal sum - as positive Bernstein/Bezier sum - as positive even Bernstein/Bezier sum -- as positive monotonic Bernstein/Bezier sum +- as positive monotonic Bernstein/Bezier sum - as positive convex/concave Bernstein/Bezier sum -- as positive monotonic convex/concave Bernstein/Bezier sum -- as generic spline (b-spline) -- as positive spline (p-spline) +- as positive monotonic convex/concave Bernstein/Bezier sum +- as generic spline (b-spline) +- as positive spline (p-spline) - as positive monotonic spline (m-spline) - as positive monotonic convex/concave spline (c-spline) -- as positive convex/concave spline (convex/concave-spline) +- as positive convex/concave spline (convex/concave-spline) where possible, fit starts from reasonable approximatuion, taken from (1) @@ -109,13 +109,13 @@ - as positive Bernstein/Bezier sum - as positive even Bernstein/Bezier sum -- as positive monotonic Bernstein/Bezier sum +- as positive monotonic Bernstein/Bezier sum - as positive convex/concave Bernstein/Bezier sum -- as positive monotonic convex/concave Bernstein/Bezier sum -- as positive spline (p-spline) +- as positive monotonic convex/concave Bernstein/Bezier sum +- as positive spline (p-spline) - as positive monotonic spline (m-spline) - as positive monotonic convex/concave spline (c-spline) -- as positive convex/concave spline (convex/concave-spline) +- as positive convex/concave spline (convex/concave-spline) Typical usage: @@ -139,8 +139,10 @@ bezier_sum , bernstein_sum , beziereven_sum , - bernsteineven_sum ) -from ostap.core.ostap_types import integer_types, long_type, num_types + bernsteineven_sum , + rational_fun ) +from ostap.core.ostap_types import integer_types, long_type, num_types +from ostap.utils.utils import vrange from collections import namedtuple import ROOT, math # ============================================================================= @@ -224,8 +226,7 @@ def _h1_param_sum_ ( h1 , for i in range ( 0, b.npars() ) : fun.SetParameter ( i , 0 ) - - + if not opts : opts = 'S' if not 'S' in opts.upper() : opts += 'S' @@ -595,6 +596,159 @@ def _h1_legendre_fast_ ( h1 , degree , xmin = inf_neg , xmax = inf_pos ) : return func +# ============================================================================= +## represent 1D-histo as rational function. based on Floater-Hormann's +# rational barycentric interpolant +# +# @code +# h = ... ## the histogram +# b = h.rational ( 5, 2 ) ## make a fit... +# +# tf1 = b[0] ## TF1 object +# obj = b[1] ## helper object +# fun = b[2] ## underlying normalzed C++ object +# fit_result = b[3] ## fit result & status +# norm = b[4] ## normalization +# +# x = ... +# print 'TF1(%s) = %s' % ( x , tf1 ( x ) ) +# print 'fun(%s) = %s' % ( x , fun ( x ) * norm ) +# @endcode +def _h1_rational_ ( h1 , + n , ## degree of numerator + d , + opts = 'SQ0' , + xmin = inf_neg , + xmax = inf_pos , + fixes = () , + params = () , + limits = () , + refit = 1 ) : + """Represent histo as ratioal fnuction, based on Floater-Hormann's + rational barycentric interpolant + + >>> h = ... # the histogram + >>> b = h.rational ( 5 , 2 ) ## make a fit... + + >>> tf1 = b[0] ## TF1 object + >>> obj = b[1] ## helper object + >>> fun = b[2] ## underlying normalzed C++ object + >>> fit_result = b[3] ## fit result & status + >>> norm = b[4] ## normalization + + >>> x = ... + >>> print 'TF1(%s) = %s' % ( x , tf1 ( x ) ) + >>> print 'FUN(%s) = %s' % ( x , norm * fun ( x ) ) + """ + xmin = max ( xmin , h1.xmin() ) + xmax = min ( xmax , h1.xmax() ) + # make reasonable approximation + func = rational_fun ( h1 , n , d , xmin , xmax ) + ## make a fit + if not params : params = tuple ( [ p for p in func.pars() ] ) + from ostap.fitting.param import H_fit + return _h1_param_sum_ ( h1 , + func , + H_fit , + opts = opts , + xmin = xmin , + xmax = xmax , + fixes = fixes , + params = params , + limits = limits , + refit = refit ) + + +# ============================================================================= +## represent 1D-histo as rational function - ration of Bernstein and positive Bernstein polynomials +# +# @code +# h = ... ## the histogram +# b = h.brational ( 3 , 3 ) ## make a fit... +# +# tf1 = b[0] ## TF1 object +# obj = b[1] ## helper object +# fun = b[2] ## underlying normalzed C++ object +# fit_result = b[3] ## fit result & status +# norm = b[4] ## normalization +# +# x = ... +# print 'TF1(%s) = %s' % ( x , tf1 ( x ) ) +# print 'fun(%s) = %s' % ( x , fun ( x ) * norm ) +# @endcode +def _h1_brational_ ( h1 , + n , ## degree of numerator + d , ## degree of denominator + opts = 'SQ0' , + xmin = inf_neg , + xmax = inf_pos , + fixes = () , + params = () , + limits = () , + refit = 1 ) : + """Represent 1D-histo as rational function - ratio + of Bernstein and positive Bernstein polynomials + + >>> h = ... # the histogram + >>> b = h.brational ( 3 , 3 ) ## make a fit... + + >>> tf1 = b[0] ## TF1 object + >>> obj = b[1] ## helper object + >>> fun = b[2] ## underlying normalzed C++ object + >>> fit_result = b[3] ## fit result & status + >>> norm = b[4] ## normalization + + >>> x = ... + >>> print 'TF1(%s) = %s' % ( x , tf1 ( x ) ) + >>> print 'FUN(%s) = %s' % ( x , norm * fun ( x ) ) + """ + xmin = max ( xmin , h1.xmin() ) + xmax = min ( xmax , h1.xmax() ) + + ii = 1.0 / ( xmax - xmin ) + + # make reasonable approximation + func = Ostap.Math.RationalBernstein ( n, d , xmin , xmax ) + + first = not params and not fixes and not limits + + if not params and not fixes : + np = func.pnpars() + for i , x in enumerate ( vrange ( xmin , xmax , n ) ) : + func.setPar ( i , float ( h1 ( x ) ) * ii / max ( 1 , np ) ) + + from ostap.fitting.param import H_fit + + if not fixes : + + nopts = opts.lower() + if not 'q' in nopts : nopts += 'q' + if not '0' in nopts : nopts += '0' + + ## fix denominator + p0 = tuple ( [ func.par ( i ) for p in range ( func.npars() ) ] ) + f0 = [ ( i , func.par ( i ) ) for i in range ( np , func.npars() ) ] + r0 = _h1_param_sum_ ( h1 , + func , + H_fit , + opts = nopts , + xmin = xmin , + xmax = xmax , + fixes = f0 , + params = p0 ) + + ## make a fit + if not params : params = tuple ( [ p for p in func.pars() ] ) + return _h1_param_sum_ ( h1 , + func , + H_fit , + opts = opts , + xmin = xmin , + xmax = xmax , + fixes = fixes , + params = params , + limits = limits , + refit = refit ) try : @@ -1780,9 +1934,11 @@ def _h1_karlinstudden_ ( h1 , t.chebyshev = _h1_chebyshev_ t.legendre = _h1_legendre_ t.polynomial = _h1_polinomial_ + t.rational = _h1_rational_ + t.brational = _h1_brational_ t.positive = _h1_positive_ t.positiveeven = _h1_positiveeven_ - t.monotonic = _h1_monotonic_ + t.monotonic = _h1_monotonic_ t.convex = _h1_convex_ t.convexpoly = _h1_convexpoly_ t.concavepoly = _h1_concavepoly_ @@ -1818,7 +1974,9 @@ def _h1_karlinstudden_ ( h1 , _h1_concavespline_ , _h1_legendre_fast_ , _h1_karlinshapley_ , - _h1_karlinstudden_ + _h1_karlinstudden_ , + _h1_rational_ , + _h1_brational_ , ] # ============================================================================= @@ -2692,6 +2850,41 @@ def _h1_beziereven_sum_ ( h1 , N , **kwargs ) : return beziereven_sum ( h1 , N , xmin , xmax ) + +# ============================================================================= +## make a histogram representation in terms of Bezier(Bernstein) sum +# (sum over Bernstein polynomials) using even polynomials: +# \f$ f( \frac{x_{min}+x_{max}}{2} - x ) = f( \frac{x_{min}+x_{max}}{2} + x ) \f$ +# @code +# histo = ... +# fsum = histo.beziereven_sum ( 2 ) +# fsum = histo.bernsteineven_sum ( 2 ) ## ditto +# print fsum +# x = ... +# print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) +# @endcode +# @see Ostap::Math::BernsteinEven +# @param h1 the histogram +# @param N the degree of polynomial +# @author Vanya Belyaev Ivan.Belyaev@itep.ru +# @date 2015-07-26 +def _h1_rational_fun_ ( h1 , n , d , **kwargs ) : + """Make a histogram representation in terms of Rational function + inspired by Floated-Hormann's rational bary centric interpolant + >>> histo = ... + >>> rat = rational_fun ( 2 ) ## ditto + >>> print ( rat ) + >>> x = ... + >>> print 'FUN(%s) = %s ' % ( x , rat ( x ) ) + """ + ## + xmin = max ( kwargs.get ( 'xmin' , h1.xmin() ) , h1.xmin () ) + xmax = min ( kwargs.get ( 'xmax' , h1.xmax() ) , h1.xmax () ) + ## + return rational_fun ( h1 , n , d , xmin , xmax ) + + + for h in ( ROOT.TH1F , ROOT.TH1D ) : @@ -2711,6 +2904,8 @@ def _h1_beziereven_sum_ ( h1 , N , **kwargs ) : h.bezier_sum_orig = _h1_bezier_sum_orig_ h.bernstein_sum_orig = _h1_bezier_sum_orig_ + h.rational_fun = _h1_rational_fun_ + h.legendre_sum = _h1_legendre_sum_orig_ ## ATTENTION h.chebyshev_sum = _h1_chebyshev_sum_orig_ ## ATTENTION h.bezier_sum = _h1_bezier_sum_orig_ ## ATTENTION @@ -2733,7 +2928,9 @@ def _h1_beziereven_sum_ ( h1 , N , **kwargs ) : _new_methods_ .append ( h.beziereven_sum ) _new_methods_ .append ( h.bernsteineven_sum ) - + + _new_methods_ .append ( h.rational_fun ) + # ============================================================================== ## (relatively) fast parameterization of 2D histogram as sum of @@ -3038,6 +3235,7 @@ def _h3_bernstein_fast_ ( h3 , NX , NY , NZ , _new_methods_ .append ( h.bezier ) _new_methods_ .append ( h.bezier_fast ) _new_methods_ .append ( h.bezier_fill ) + _new_methods_ .append ( _h3_legendre_fast_ ) # ============================================================================= diff --git a/ostap/histos/tests/test_histos_parameterisation.py b/ostap/histos/tests/test_histos_parameterisation.py index 7f2b76b2..2b276cc4 100755 --- a/ostap/histos/tests/test_histos_parameterisation.py +++ b/ostap/histos/tests/test_histos_parameterisation.py @@ -328,7 +328,21 @@ def test_chebyshev_sum() : f.draw('same') logger.info ( "%-25s : difference %s" % ( h.title , diff1 ( f , h ) ) ) - +# ============================================================================= +def test_rational_fun () : + + logger = getLogger("test_rational_fun") + + n = 8 + with timing ( 'Rational_fun[%s]' % n , logger ) : + for h in histos : + with use_canvas ( 'test_rational_fun: ' + h.GetTitle() , wait = 1 ) : + h.draw() + for d in range ( 0 , n + 1 ) : + f = h.rational_fun ( n , d ) + f.draw ( 'same' , color = d + 1 ) + logger.info ( "%-25s : [%d] difference %s" % ( h.title , d , diff1 ( f , h ) ) ) + # ============================================================================= def test_fourier_sum() : @@ -337,7 +351,6 @@ def test_fourier_sum() : logger.warning("No numpy/scipy is avilable, skip the test test") return - hh = [ h for h in histos if hasattr ( h , 'fourier_sum' ) ] with timing ( 'Fourier-sum[10]' , logger ) : params = [ h.fourier_sum ( 10 ) for h in hh ] @@ -642,6 +655,8 @@ def rdif ( x , y ) : test_chebyshev_sum_fill () test_chebyshev_sum () + test_rational_fun () + test_fourier_sum () test_cosine_sum () diff --git a/ostap/histos/tests/test_histos_parameterisation2.py b/ostap/histos/tests/test_histos_parameterisation2.py index 196efbf4..5edd082b 100755 --- a/ostap/histos/tests/test_histos_parameterisation2.py +++ b/ostap/histos/tests/test_histos_parameterisation2.py @@ -287,6 +287,39 @@ def test_convex_poly () : f.tf1.draw('same') logger.info ( "%-25s : difference %s" % ( h.title , diff2 ( f , h ) ) ) + +# ============================================================================= +def test_rational () : + + n = 6 + funs = set() + logger = getLogger("test_rational") + with timing ( 'Rational [%s]' % n , logger ) : + for h in ( h1 , h2 , h3 , h4 ) : + with use_canvas ( 'test_rational %s' % h.GetTitle() , wait = 2 ) : + h.draw() + for d in range ( 0 , n + 1 ) : + f = h.rational ( n = n , d = d ) + f.tf1.draw('same', color = d + 1 ) + funs.add ( f ) + logger.info ( "%-25s : [%d] difference %s" % ( h.title , d , diff2 ( f , h ) ) ) + +# ============================================================================= +def test_brational () : + + n = 3 + funs = set() + logger = getLogger("test_brational") + with timing ( 'BRational [%s]' % n , logger ) : + for h in ( h1 , h2 , h3 , h4 ) : + with use_canvas ( 'test_brational %s' % h.GetTitle() , wait = 2 ) : + h.draw() + for d in range ( 2 , n + 2 ) : + f = h.brational ( n = n , d = d ) + f.tf1.draw('same', color = d + 1 ) + funs.add ( f ) + logger.info ( "%-25s : [%d] difference %s" % ( h.title , d , diff2 ( f , h ) ) ) + # ============================================================================= def test_fourier () : @@ -389,7 +422,6 @@ def test_convex_only_spline () : f.tf1.draw('same') logger.info ( "%-25s : difference %s" % ( h.title , diff2 ( f , h ) ) ) - # ============================================================================ ## Karlin-Shapley # ============================================================================ @@ -500,6 +532,9 @@ def test_legendre3_fast () : test_convex () test_convex_poly () + test_rational () + test_brational () + test_fourier () test_cosine () diff --git a/ostap/math/param.py b/ostap/math/param.py index 5e60792a..fb88dec3 100644 --- a/ostap/math/param.py +++ b/ostap/math/param.py @@ -12,7 +12,8 @@ # - as Cosine/Fourier sum (relies on scipy.fftpack) # - as Bernstein/Bezier sum # - as even Bernstein/Bezier sum -# +# - as rational function inspired by Floater-Hormann's rationa barycentric interpolant + # Typical usage: # # @code @@ -35,7 +36,7 @@ # @author Vanya BELYAEV Ivan.Belyaev@itep.ru # @date 2011-06-07 # ============================================================================= -"""Module with utilities for parameterization of histograms +"""Module with utilities for parameterization of functions/histograms ## (1) using histogram data only: @@ -74,13 +75,14 @@ 'bernstein_sum' , ## - ditto - 'beziereven_sum' , ## even Bezier/Bernstein sum for the given function/object 'bernsteineven_sum' , ## - ditto - + 'rational_fun' , ## Rational (a'la Floater-Hormann's interpolant) parameterisaton ) # ============================================================================= -import ROOT, ctypes -# ============================================================================= from ostap.core.core import cpp, Ostap from ostap.core.ostap_types import is_integer, num_types +from ostap.utils.utils import crange import ostap.math.models +import ROOT, ctypes # ============================================================================= # logging # ============================================================================= @@ -138,7 +140,7 @@ def _get_xminmax_ ( func , xmin , xmax , name = 'get_xminmax') : # x = ... # print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) # @endcode -# @see Gaudi::Math::LegendreSum +# @see Ostap::Math::LegendreSum # @author Vanya Belyaev Ivan.Belyaev@iter.ru # It is not very CPU efficient, but stable enough... # @date 2015-07-26 @@ -150,7 +152,7 @@ def legendre_sum ( func , N , xmin , xmax , **kwargs ) : >>> print fsum >>> x = ... >>> print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) - see Gaudi::Math::LegendreSum + see Ostap.Math.LegendreSum """ from copy import deepcopy @@ -203,7 +205,7 @@ def chebyshev_sum ( func , N , xmin , xmax ) : >>> print fsum >>> x = ... >>> print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) - see Gaudi::Math::ChebyshevSum + see Ostap::Math::ChebyshevSum """ assert is_integer ( N ) and 0 <= N, "chebyshev_sum: invalid N %s " % N @@ -238,6 +240,9 @@ def chebyshev_sum ( func , N , xmin , xmax ) : return csum + + + # ============================================================================= try : # ========================================================================= @@ -251,7 +256,7 @@ def chebyshev_sum ( func , N , xmin , xmax ) : # x = ... # print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) # @endcode - # @see Gaudi::Math::FourierSum + # @see Ostap::Math::FourierSum # @author Vanya Belyaev Ivan.Belyaev@itep.ru # @date 2015-07-26 def fourier_sum ( func , N , xmin , xmax , fejer = False ) : @@ -330,7 +335,7 @@ def fourier_sum ( func , N , xmin , xmax , fejer = False ) : # x = ... # print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) # @endcode - # @see Gaudi::Math::CosineSum + # @see Ostap::Math::CosineSum # @author Vanya Belyaev Ivan.Belyaev@itep.ru # @date 2015-07-26 def cosine_sum ( func , N , xmin , xmax , fejer = False ) : @@ -381,7 +386,7 @@ def cosine_sum ( func , N , xmin , xmax , fejer = False ) : # x = ... # print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) # @endcode -# @see Gaudi::Math::Bernstein +# @see Ostap::Math::Bernstein # @param func (INPUT) the function # @param N (INPUT) the polynomial degree # @param xmin (INPUT) the low edge of the domain @@ -452,7 +457,7 @@ def bezier_sum ( func , N , xmin , xmax , **kwargs ) : # x = ... # print 'FUN(%s) = %s ' % ( x , fsum ( x ) ) # @endcode -# @see Gaudi::Math::BernsteinEven +# @see Ostap::Math::BernsteinEven # @param func (INPUT) the function # @param N (INPUT) the degree of even polynomial # @param xmin (INPUT) the low edge of the domain @@ -527,6 +532,46 @@ def _sym_func_ ( x ) : return bsum + +# ============================================================================= +## make a function representation in terms of rational function +# based on the Floater-Hormann's rational barycentric interpolant +# @code +# func = lambda x : x * x +# rat = rational( func , 3 , 3 , xmin = -1 , xmax = 1 ) +# print ( rat ) +# x = ... +# print 'FUN(%s) = %s ' % ( x , rat ( x ) ) +# @endcode +# @see Ostap::Math::Rational +# @see Ostap::Math::FloaterHormann +# @author Vanya Belyaev Ivan.Belyaev@iter.ru +# @date 2023-09-18 +def rational_fun ( func , n , d , xmin , xmax , **kwargs ) : + """ Make a function representation in terms of rational function + based on the Floater-Hormann's rational barycentric interpolant + >>> func = lambda x : x * x + >>> rat = rational ( func , N = 4 , d = 1 , xmin = -1 , xmax = 1 ) + >>> print ( rat ) + >>> x = ... + >>> print 'FUN(%s) = %s ' % ( x , rat ( x ) ) + - see Ostap.Math.Rational + - see Ostap.Math.FloaterHormann + """ + from copy import deepcopy + + assert is_integer ( n ) and 0 <= n , "rational: invalid n %s " % N + assert is_integer ( d ) and 0 <= d <= n , "rational: invalid d %s " % d + + xmin , xmax = _get_xminmax_ ( func , xmin , xmax , 'rational' ) + + ## prepare the result + rational = Ostap.Math.Rational( n + 1 , d , xmin , xmax ) + for i, r in enumerate ( crange ( xmin , xmax , n + 1 ) ) : + rational [ i ] = float ( func ( r ) ) + + return rational + # ============================================================================= ## make a function representation in terms of even Bezier sum (sum over Bernstein polynomials) bernsteineven_sum = beziereven_sum diff --git a/source/src/Interpolation.cpp b/source/src/Interpolation.cpp index f50da9a5..c1735bc7 100644 --- a/source/src/Interpolation.cpp +++ b/source/src/Interpolation.cpp @@ -557,7 +557,7 @@ Ostap::Math::Berrut2nd::Berrut2nd Ostap::Math::FloaterHormann::Weights::Weights ( const Ostap::Math::Interpolation::Abscissas& a , const unsigned short d ) - : m_d ( d <= a.size() ? d : a.size() ) + : m_d ( a.size() ? ( d < a.size() ? d : a.size() - 1 ) : 0 ) , m_weights ( a.size() ) { // @@ -565,7 +565,7 @@ Ostap::Math::FloaterHormann::Weights::Weights const unsigned int n = 1 <= nn ? nn - 1 : 0 ; // for ( unsigned int i = 0 ; i <= n ; ++i ) - { + { const long double xi = a.x ( i ) ; // long double ib = 0 ; diff --git a/source/src/Rational.cpp b/source/src/Rational.cpp index 2db40b7e..f687f425 100644 --- a/source/src/Rational.cpp +++ b/source/src/Rational.cpp @@ -27,7 +27,6 @@ */ // ============================================================================ - // ============================================================================ /* constructor * @param n degree of numerator @@ -192,9 +191,9 @@ double Ostap::Math::RationalBernstein::evaluate ( const double x ) const std::vector Ostap::Math::RationalBernstein::pars () const { - std::vector result ( npars() , 0 ) ; - std::copy ( m_p.pars().begin() , m_p.pars().end() , result.begin() ) ; - std::copy ( m_q.pars().begin() , m_q.pars().end() , result.begin() + m_p.npars() ) ; + const unsigned short np = npars() ; + std::vector result ( np , 0.0 ) ; + for ( unsigned short i = 0 ; i < np ; ++i ) { result [ i ] = par ( i ) ; } return result ; } // ============================================================================