From fc26c7e1b9bc0a92a610512da4fe7cdba3dfa4ba Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Thu, 21 Sep 2023 16:02:51 +0200 Subject: [PATCH] 1. Add `RationalFun` & `RationalBernsteinFun` FUN --- ReleaseNotes/release_notes.md | 4 +- ostap/fitting/roofuncs.py | 118 ++++++++++++++++--- ostap/fitting/tests/test_fitting_roofuncs.py | 30 ++--- source/include/Ostap/MoreVars.h | 18 +++ source/src/MoreVars.cpp | 13 +- source/src/Rational.cpp | 1 - 6 files changed, 143 insertions(+), 41 deletions(-) diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index 5d740306..a1d40d42 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -1,7 +1,7 @@ ## New features: - 1. Add `Ostap::MoreRooFit::Rational` and `Ostap::MoreRooFit::RatioalBernstein` - 1. Add `RationaFun` FUN + 1. Add `Ostap::MoreRooFit::Rational` and `Ostap::MoreRooFit::RationalBernstein` + 1. Add `RationalFun` & `RationalBernsteinFun` FUN ## Backward incompatible: diff --git a/ostap/fitting/roofuncs.py b/ostap/fitting/roofuncs.py index 7aae29ac..5f95da19 100644 --- a/ostap/fitting/roofuncs.py +++ b/ostap/fitting/roofuncs.py @@ -14,17 +14,18 @@ __date__ = "2020-03-08" # ============================================================================= __all__ = ( - 'BernsteinPoly' , ## generic polynomial in Bernstein form (RooAbsReal) - 'MonotonicPoly' , ## monotonic polynomial (RooAbsReal) - 'ConvexPoly' , ## monotonic convex/concave polynomial (RooAbsReal) - 'ConvexOnlyPoly' , ## convex/concave polynomial (RooAbsReal) - 'RooPoly' , ## simple wrapper for RooPolyVar (RooAbsReal) - 'ScaleAndShift' , ## scale and shift (RooAbsReal) - 'BSplineFun' , ## BSpline (RooAbsReal) - 'RationalFun' , ## Ratioal function (RooAbsReal) - 'Shape1D_fun' , ## arbitrary fixed shape (RooAbsReal) - 'Histo1D_fun' , ## fixed shap form historgam (RooAbsReal) - 'Histo1DErr_fun' , ## fixed shap from historgam errors (RooAbsReal) + 'BernsteinPoly' , ## generic polynomial in Bernstein form (RooAbsReal) + 'MonotonicPoly' , ## monotonic polynomial (RooAbsReal) + 'ConvexPoly' , ## monotonic convex/concave polynomial (RooAbsReal) + 'ConvexOnlyPoly' , ## convex/concave polynomial (RooAbsReal) + 'RooPoly' , ## simple wrapper for RooPolyVar (RooAbsReal) + 'ScaleAndShift' , ## scale and shift (RooAbsReal) + 'BSplineFun' , ## BSpline (RooAbsReal) + 'RationalFun' , ## Ratioal function (RooAbsReal) + 'RationalBernsteinFun' , ## Ratioal function (RooAbsReal) + 'Shape1D_fun' , ## arbitrary fixed shape (RooAbsReal) + 'Histo1D_fun' , ## fixed shap form historgam (RooAbsReal) + 'Histo1DErr_fun' , ## fixed shap from historgam errors (RooAbsReal) ## ## 'var_sum' , ## sum for RooAbsReal objects ## 'var_mul' , ## product for RooAbsReal objects @@ -59,7 +60,7 @@ ) # ============================================================================= from ostap.core.core import Ostap, VE -from ostap.core.ostap_types import num_types +from ostap.core.ostap_types import num_types, integer_types from ostap.fitting.fithelpers import ParamsPoly , ShiftScalePoly import ostap.fitting.variables from ostap.fitting.funbasic import FUN1, Fun1D, Fun2D, Fun3D @@ -540,7 +541,6 @@ def power ( self ) : return self.fun.degree() - # ============================================================================= ## @class RationalFun # A simple pole-free rational function at interval \f$ x_{min} \le x \le x_{max}\f$ @@ -561,9 +561,12 @@ class RationalFun(FUN1,ParamsPoly) : - see Ostap.MoreRooFit.Bernstein - see Ostap.Math.Bernstein """ - def __init__ ( self , name , xvar , - n = 3 , - d = 1 , pars = None ) : + def __init__ ( self , + name , + xvar , + n = 3 , + d = 1 , + pars = None ) : ## initialize the base class FUN1 .__init__ ( self , name , xvar = xvar ) @@ -571,8 +574,8 @@ def __init__ ( self , name , xvar , npars = n + 1 , pars = pars ) - assert isinstance ( n , int ) and 1 <= n , 'Invalid parameter n' - assert isinstance ( d , int ) and 0 <= d <= n , 'Invalid parameter d' + assert isinstance ( n , integer_types ) and 0 <= n , 'Invalid parameter n' + assert isinstance ( d , integer_types ) and 0 <= d <= n , 'Invalid parameter d' xmin , xmax = self.xminmax () @@ -603,7 +606,84 @@ def n ( self ) : def d ( self ) : """'d' : d-parameter of rational function""" return self.fun.d () - + + + +# ============================================================================= +## @class RationalBernsteinFun +# A simple rational function at interval \f$ x_{min} \le x \le x_{max}\f$ as ratio +# of Bernstein and positive Bernstein polynomials +# \f[ F(x) = \frac{p(x)}{q(x)} \f] +# @see Ostap::MoreRooFit::RationalBernstein +# @see Ostap::Math::RationalBernstein +# @see Ostap::Math::Bernstein +# @see Ostap::Math::Positive +# @author Vanya BELYAEV Ivan.Belyaev@cern.ch +# @date 2023-09-21 +class RationalBernsteinFun(FUN1,ParamsPoly) : + """ A simple rational function at interval \f$ x_{min} \le x \le x_{max}\f$ as ratio + of Bernstein and positive Bernstein polynomials + - see Ostap.MoreRooFit.RationalBernstein + - see Ostap.Math.RationalBernstein + - see Ostap.Math.Bernstein + - see Ostap.Math.Positive + """ + def __init__ ( self , + name , + xvar , + p = 3 , + q = 3 , + pars = None ) : + + assert isinstance ( p , integer_types ) and 0 <= p , 'Invalid parameter p' + assert isinstance ( q , integer_types ) and 0 <= q , 'Invalid parameter q' + + ## initialize the base class + FUN1 .__init__ ( self , name , xvar = xvar ) + ParamsPoly.__init__ ( self , + npars = p + q + 1 , + pars = pars ) + + assert self.pars , 'Invalid number of parameters!' + p = min ( p , self.npars - 1 ) + + if not pars : + for i, v in enumerate ( self.pars ) : + if i < p + 1 : v.setVal ( 1 ) + else : + v.setMin ( -5 * math.pi ) + v.setMax ( +5 * math.pi ) + + xmin , xmax = self.xminmax () + + ## create the function + self.fun = Ostap.MoreRooFit.RationalBernstein ( + self.roo_name ( 'rbf_' ) , + 'RationalBernstein %s' % self.name , + self.xvar , + self.pars_lst , + p , + xmin , + xmax ) + + ## self.tricks = True + self.config = { + 'name' : self.name , + 'xvar' : self.xvar , + 'p' : self.fun.p () , + 'q' : self.fun.q () , + 'pars' : self.pars , + } + + @property + def p ( self ) : + """'p' : p-parameter of rational function""" + return self.fun.p () + @property + def q ( self ) : + """'q' : q-parameter of rational function""" + return self.fun.q () + # ============================================================================= ## Generic 1D-shape from C++ callable diff --git a/ostap/fitting/tests/test_fitting_roofuncs.py b/ostap/fitting/tests/test_fitting_roofuncs.py index e112bbe8..fb514ece 100644 --- a/ostap/fitting/tests/test_fitting_roofuncs.py +++ b/ostap/fitting/tests/test_fitting_roofuncs.py @@ -42,7 +42,7 @@ # ============================================================================= ## test basic operationsw with functions def test_roofuncs_1 () : - """Test basic operationsw with functions + """Test basic operations with functions """ logger = getLogger ( 'test_roofuncs_1' ) @@ -50,23 +50,23 @@ def test_roofuncs_1 () : power = 6 - B = RF.BernsteinPoly ( 'Bernstein' , xvar = x , power = power ) - I = RF.MonotonicPoly ( 'BernsteinI' , xvar = x , power = power , increasing = True ) - D = RF.MonotonicPoly ( 'DernsteinD' , xvar = x , power = power , increasing = False ) - IX = RF.ConvexPoly ( 'ConvexIX' , xvar = x , power = power , increasing = True , convex = True ) - IV = RF.ConvexPoly ( 'ConvexIV' , xvar = x , power = power , increasing = True , convex = False ) - DX = RF.ConvexPoly ( 'ConvexDX' , xvar = x , power = power , increasing = False , convex = True ) - DV = RF.ConvexPoly ( 'ConvexDV' , xvar = x , power = power , increasing = False , convex = False ) - CX = RF.ConvexOnlyPoly ( 'ConvexOnlyX' , xvar = x , power = power , convex = True ) - CV = RF.ConvexOnlyPoly ( 'ConvexOnlyV' , xvar = x , power = power , convex = False ) + B = RF.BernsteinPoly ( 'Bernstein' , xvar = x , power = power ) + I = RF.MonotonicPoly ( 'BernsteinI' , xvar = x , power = power , increasing = True ) + D = RF.MonotonicPoly ( 'DernsteinD' , xvar = x , power = power , increasing = False ) + IX = RF.ConvexPoly ( 'ConvexIX' , xvar = x , power = power , increasing = True , convex = True ) + IV = RF.ConvexPoly ( 'ConvexIV' , xvar = x , power = power , increasing = True , convex = False ) + DX = RF.ConvexPoly ( 'ConvexDX' , xvar = x , power = power , increasing = False , convex = True ) + DV = RF.ConvexPoly ( 'ConvexDV' , xvar = x , power = power , increasing = False , convex = False ) + CX = RF.ConvexOnlyPoly ( 'ConvexOnlyX' , xvar = x , power = power , convex = True ) + CV = RF.ConvexOnlyPoly ( 'ConvexOnlyV' , xvar = x , power = power , convex = False ) + RT = RF.RationalFun ( 'Rational' , xvar = x , n = power , d = 2 ) + RB = RF.RationalBernsteinFun ( 'RationalBernstein' , xvar = x , p = 3 , q = 3 ) - for p in ( B , I , D , IX , IV , DX , DV , CX , CV ) : - p.pars = [ random.uniform ( -5 , 5 ) for i in range ( power + 2 ) ] - - with wait ( 1 ), use_canvas ( 'test_roofuncs:%s' % p.name ) : p.draw () + for p in ( B , I , D , IX , IV , DX , DV , CX , CV , RT , RB ) : + p.pars = [ random.uniform ( -5 , 5 ) for i in range ( power + 2 ) ] + with use_canvas ( 'test_roofuncs:%s' % p.name , wait = 1 ) : p.draw () funs.add ( p ) - # ============================================================================= ## check that everything is serializable # ============================================================================= diff --git a/source/include/Ostap/MoreVars.h b/source/include/Ostap/MoreVars.h index 7a1bf101..7c8ab8fe 100644 --- a/source/include/Ostap/MoreVars.h +++ b/source/include/Ostap/MoreVars.h @@ -88,6 +88,8 @@ namespace Ostap public: // ====================================================================== /// get the varibale + const RooAbsReal& x () const { return m_xvar.arg() ; } + /// get the varibale const RooAbsReal& xvar () const { return m_xvar.arg() ; } /// get parameters const RooArgList& pars () const { return m_pars ; } @@ -186,6 +188,8 @@ namespace Ostap public: // ====================================================================== /// get the varibale + const RooAbsReal& x () const { return m_xvar.arg() ; } + /// get the varibale const RooAbsReal& xvar () const { return m_xvar.arg() ; } /// get parameters const RooArgList& pars () const { return m_pars ; } @@ -283,6 +287,8 @@ namespace Ostap public: // ====================================================================== /// get the varibale + const RooAbsReal& x () const { return m_xvar.arg() ; } + /// get the varibale const RooAbsReal& xvar () const { return m_xvar.arg() ; } /// get parameters const RooArgList& pars () const { return m_pars ; } @@ -386,6 +392,8 @@ namespace Ostap public: // ====================================================================== /// get the variable + const RooAbsReal& x () const { return m_xvar.arg() ; } + /// get the variable const RooAbsReal& xvar () const { return m_xvar.arg() ; } /// get parameters const RooArgList& pars () const { return m_pars ; } @@ -462,6 +470,8 @@ namespace Ostap public: // ====================================================================== /// get the varibale + const RooAbsReal& x () const { return m_xvar.arg() ; } + /// get the varibale const RooAbsReal& xvar () const { return m_xvar.arg() ; } /// get parameters const RooArgList& pars () const { return m_pars ; } @@ -545,6 +555,8 @@ namespace Ostap public: // ====================================================================== /// get the variable/observable + const RooAbsReal& x () const { return m_xvar.arg() ; } + /// get the variable/observable const RooAbsReal& xvar () const { return m_xvar.arg() ; } /// get parameters const RooArgList& pars () const { return m_pars ; } @@ -552,6 +564,8 @@ namespace Ostap unsigned short n () const { return m_rational.n () ; } /// get d unsigned short d () const { return m_rational.d () ; } + /// get p (==d) + unsigned short p () const { return m_rational.d () ; } /// xmin for Rational double xmin () const { return m_rational.xmin () ; } /// xmax for Rational @@ -634,11 +648,15 @@ namespace Ostap public: // ====================================================================== /// get the variable/observable + const RooAbsReal& x () const { return m_xvar.arg() ; } + /// get the variable/observable const RooAbsReal& xvar () const { return m_xvar.arg() ; } /// get parameters const RooArgList& pars () const { return m_pars ; } /// get p unsigned short p () const { return m_rational.pdegree () ; } + /// get q + unsigned short q () const { return m_rational.qdegree () ; } /// xmin for Rational double xmin () const { return m_rational.xmin () ; } /// xmax for Rational diff --git a/source/src/MoreVars.cpp b/source/src/MoreVars.cpp index 450bfed3..768ec44b 100644 --- a/source/src/MoreVars.cpp +++ b/source/src/MoreVars.cpp @@ -801,13 +801,18 @@ Ostap::MoreRooFit::RationalBernstein::RationalBernstein : RooAbsReal ( name.c_str () , title.c_str () ) , m_xvar ( "!x" , "Observable" , this , xvar ) , m_pars ( "!pars" , "Parameters" , this ) - , m_rational ( ::size ( pars ) , p , xmin , xmax ) + , m_rational ( 0 , 0 , xmin , xmax ) // TMP!! { // Ostap::Assert ( 1 <= ::size ( pars ) , s_EMPTYPARS , s_v7 , 510 ) ; // ::copy_real ( pars , m_pars , s_INVALIDPAR , s_v7 ) ; // + const unsigned short np = ::size ( m_pars ) ; + const unsigned short pp = p + 1 <= np ? p : np - 1 ; + const unsigned short qq = ( np - 1 ) - pp ; + m_rational = Ostap::Math::RationalBernstein ( pp , qq , xmin , xmax ) ; + // Ostap::Assert ( m_rational.npars() == ::size ( m_pars ), s_INVALIDPARS , s_v7 , 512 ) ; // } @@ -818,17 +823,17 @@ Ostap::MoreRooFit::RationalBernstein::RationalBernstein ( const std::string& name , const std::string& title , RooAbsReal& xvar , - const RooArgList& p , + const RooArgList& p , const RooArgList& q , const double xmin , const double xmax ) : RooAbsReal ( name.c_str () , title.c_str () ) , m_xvar ( "!x" , "Observable" , this , xvar ) , m_pars ( "!pars" , "Parameters" , this ) - , m_rational ( ::size ( p ) , ::size ( q ) , xmin , xmax ) + , m_rational ( 1 <= ::size ( p ) ? ::size ( p ) - 1 : 0 , ::size ( q ) , xmin , xmax ) { // - Ostap::Assert ( 1 <= ::size ( p ) + ::size ( q ) , s_EMPTYPARS , s_v7 , 510 ) ; + Ostap::Assert ( 1 <= ::size ( p ) , s_EMPTYPARS , s_v7 , 510 ) ; // ::copy_real ( p , m_pars , s_INVALIDPAR , s_v7 ) ; ::copy_real ( q , m_pars , s_INVALIDPAR , s_v7 ) ; diff --git a/source/src/Rational.cpp b/source/src/Rational.cpp index 1b3dfe82..9b1cba26 100644 --- a/source/src/Rational.cpp +++ b/source/src/Rational.cpp @@ -149,7 +149,6 @@ std::size_t Ostap::Math::Rational::tag () const - // ============================================================================ Ostap::Math::RationalBernstein::RationalBernstein ( const unsigned short p ,