Skip to content

Commit

Permalink
1. Add Ostap::MoreRooFit::Rational and `Ostap::MoreRooFit::Ratioal…
Browse files Browse the repository at this point in the history
…Bernstein`

  1.  Add `RationaFun` FUN
  • Loading branch information
VanyaBelyaev committed Sep 21, 2023
1 parent ba22586 commit 624ec49
Show file tree
Hide file tree
Showing 9 changed files with 504 additions and 21 deletions.
3 changes: 3 additions & 0 deletions ReleaseNotes/release_notes.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,8 @@
## New features:

1. Add `Ostap::MoreRooFit::Rational` and `Ostap::MoreRooFit::RatioalBernstein`
1. Add `RationaFun` FUN

## Backward incompatible:

## Bug fixes:
Expand Down
66 changes: 66 additions & 0 deletions ostap/fitting/roofuncs.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
'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)
Expand Down Expand Up @@ -538,6 +539,71 @@ def power ( self ) :
"""'power' : degree for BSpline object"""
return self.fun.degree()



# =============================================================================
## @class RationalFun
# A simple pole-free rational function at interval \f$ x_{min} \le x \le x_{max}\f$
# \f[ F(x) = \frac{p(x)}{q(x)} \f]
# Actually internally it uses
# the Floater-Hormann rational barycentric interpolant
# and parameters are the function values at Chebyshev's nodes
#
# @see Ostap::MoreRooFit::Rational
# @see Ostap::Math::Rational
# @see Ostap::Math::FloaterHormann
# @author Vanya BELYAEV [email protected]
# @date 2023-09-21
class RationalFun(FUN1,ParamsPoly) :
r"""A simple pole-free rational function at interval \f$ x_{min} \le x \le x_{max}\f$
>>> p1 = BernsteinPoly ( 'P1' , xvar = (0,1) , power = 3 )
>>> p2 = BernsteinPoly ( 'P2' , xvar = (0,1) , pars = ... )
- see Ostap.MoreRooFit.Bernstein
- see Ostap.Math.Bernstein
"""
def __init__ ( self , name , xvar ,
n = 3 ,
d = 1 , pars = None ) :

## initialize the base class
FUN1 .__init__ ( self , name , xvar = xvar )
ParamsPoly.__init__ ( self ,
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'

xmin , xmax = self.xminmax ()

## create the function
self.fun = Ostap.MoreRooFit.Rational (
self.roo_name ( 'rfun_' ) ,
'Rational %s' % self.name ,
self.xvar ,
self.pars_lst ,
d ,
xmin ,
xmax )

## self.tricks = True
self.config = {
'name' : self.name ,
'xvar' : self.xvar ,
'n' : self.fun.n () ,
'd' : self.fun.d () ,
'pars' : self.pars ,
}

@property
def n ( self ) :
"""'n' : n-parameter of rational function"""
return self.fun.n ()
@property
def d ( self ) :
"""'d' : d-parameter of rational function"""
return self.fun.d ()


# =============================================================================
## Generic 1D-shape from C++ callable
Expand Down
8 changes: 6 additions & 2 deletions ostap/fitting/rooreduce.py
Original file line number Diff line number Diff line change
Expand Up @@ -2668,7 +2668,8 @@ def _rgauss3d_reduce_ ( pdf ):
# =============================================================================
## reduce Rational
def _rational_reduce_ ( pdf ):
"""Reduce Ostap::Models::Rational"""
"""Reduce Ostap::Models::Rational & Ostap::MoreRooFit::Rational
"""
return root_store_factory , ( type ( pdf ) ,
pdf.name ,
pdf.title ,
Expand All @@ -2678,7 +2679,10 @@ def _rational_reduce_ ( pdf ):
pdf.xmin () ,
pdf.xmax () )

Ostap.Models.Rational.__reduce__ = _rational_reduce_
Ostap.Models.Rational .__reduce__ = _rational_reduce_
Ostap.MoreRooFit.Rational .__reduce__ = _rational_reduce_
Ostap.MoreRooFit.RationalBernstein .__reduce__ = _rational_reduce_



# =============================================================================
Expand Down
20 changes: 12 additions & 8 deletions ostap/fitting/tests/test_fitting_transform.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
__author__ = "Ostap developers"
__all__ = () ## nothing to import
# =============================================================================
import ostap.fitting.roofit
from ostap.core.meta_info import root_info
from ostap.core.core import cpp, VE, dsID
import ostap.fitting.models as Models
from ostap.fitting.transform_pdf import TrPDF
Expand All @@ -23,6 +23,7 @@
from ostap.utils.timing import timing
from ostap.plotting.canvas import use_canvas
from ostap.utils.utils import wait
import ostap.fitting.roofit
import ROOT, random
# =============================================================================
# logging
Expand All @@ -40,11 +41,14 @@ def test_transform () :

logger = getLogger ( 'test_transform' )

if not 62000 <= ROOT.gROOT.GetVersionInt() :
logger.info ("Not for this version of ROOT")
if not (6,20) <= root_info :
logger.info ("Not for this version of ROOT!")
return

lmin = 0
lmax = 5

x = ROOT.RooRealVar('x','',1,100000)
x = ROOT.RooRealVar('x','', 10.**lmin , 10.**lmax )
mean = 1000
sigma = 1000

Expand All @@ -58,15 +62,15 @@ def test_transform () :
dataset.add ( varset )

dataset.add_var ( 'lx' , 'log10(x)' )
dataset.lx.setMin ( 0 )
dataset.lx.setMax ( 5 )
dataset.lx.setMin ( lmin )
dataset.lx.setMax ( lmax )

## original PDF
gauss = Models.Gauss_pdf ( 'G' , xvar = x ,
mean = ( mean , mean / 2 , mean * 2 ) ,
sigma = ( sigma , sigma / 2 , sigma * 2 ) )

with use_canvas ( 'test_transform' ) , wait ( 1 ) :
with use_canvas ( 'test_transform/x' ) , wait ( 1 ) :
r1 , f1 = gauss.fitTo ( dataset , draw = True , silent = True )
logger.info ( 'Fit x:\n%s' % r1.table() )

Expand All @@ -77,7 +81,7 @@ def test_transform () :
## transformed PDF
tgauss = TrPDF ( pdf = gauss , new_var = NX )

with use_canvas ( 'test_transform' ) , wait ( 1 ):
with use_canvas ( 'test_transform/log10' ) , wait ( 1 ):
r2 , f2 = tgauss.fitTo ( dataset , draw = True , silent = True )
logger.info ( 'Fit log10(x):\n%s' % r2.table() )

Expand Down
22 changes: 22 additions & 0 deletions source/include/Ostap/MoreRooFit.h
Original file line number Diff line number Diff line change
Expand Up @@ -2729,6 +2729,28 @@ namespace Ostap
/// clone
AddDeps* clone ( const char* newname ) const override ;
// ======================================================================
public: // delegation
// ======================================================================
Double_t analyticalIntegral
( Int_t code ,
const char* range = nullptr ) const override ;
//
Double_t analyticalIntegralWN
( Int_t code ,
const RooArgSet* normset ,
const char* range = nullptr ) const override ;
//
Int_t getAnalyticalIntegral
( RooArgSet& allVars ,
RooArgSet& analVars ,
const char* range = nullptr ) const override ;
//
Int_t getAnalyticalIntegralWN
( RooArgSet& allVars ,
RooArgSet& analVars ,
const RooArgSet* normset ,
const char* range = nullptr ) const override ;
// ======================================================================
public:
// ======================================================================
/// the actual evaluation of the result
Expand Down
173 changes: 172 additions & 1 deletion source/include/Ostap/MoreVars.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include "Ostap/Bernstein.h"
#include "Ostap/Bernstein1D.h"
#include "Ostap/BSpline.h"
#include "Ostap/Rational.h"
#include "Ostap/HistoInterpolators.h"
// ============================================================================
/// forward declarations
Expand Down Expand Up @@ -494,7 +495,177 @@ namespace Ostap
// ======================================================================
} ; // The end of class Ostap::MoreRooFit::BSpline
// ========================================================================
/** @class Shape
/** @class Rational
* A simple pole-free rational function at interval \f$ x_{min} \le x \le x_{max}\f$
* \f[ F(x) = \frac{p(x)}{q(x)} \f]
* Actually internally it uses
* the Floater-Hormann rational barycentric interpolant
* and parameters are the function valeus at Chebyshev's nodes
*
* @see Ostap::Math::Rational
* @see Ostap::Math::FloaterHormann
* @author Vanya BELYAEV [email protected]
* @date 2023-09-21
*/
class Rational final : public RooAbsReal
{
public:
// ======================================================================
ClassDefOverride(Ostap::MoreRooFit::Rational, 1) ;
// ======================================================================
public:
// ======================================================================
Rational
( const std::string& name ,
const std::string& title ,
RooAbsReal& xvar ,
const RooArgList& pars ,
const unsigned short d ,
const double xmin ,
const double xmax ) ;
/// copy constructor
Rational ( const Rational& right , const char* name = nullptr ) ;
/// default
Rational () ;
/// clone method
Rational* clone ( const char* name ) const override ;
/// virtual destructor
virtual ~Rational() ;
// ======================================================================
public:
// ======================================================================
Int_t getAnalyticalIntegral
( RooArgSet& allVars ,
RooArgSet& analVars ,
const char* rangeName = nullptr ) const override ;
Double_t analyticalIntegral
( Int_t code ,
const char* rangeName = nullptr ) const override ;
// ======================================================================
public:
// ======================================================================
/// get the variable/observable
const RooAbsReal& xvar () const { return m_xvar.arg() ; }
/// get parameters
const RooArgList& pars () const { return m_pars ; }
/// get n
unsigned short n () const { return m_rational.n () ; }
/// get d
unsigned short d () const { return m_rational.d () ; }
/// xmin for Rational
double xmin () const { return m_rational.xmin () ; }
/// xmax for Rational
double xmax () const { return m_rational.xmax () ; }
// ======================================================================
public:
// ======================================================================
void setPars () const ;
// ======================================================================
public:
// ======================================================================
const Ostap::Math::Rational& function () const { return m_rational ; }
const Ostap::Math::Rational& rational () const { return m_rational ; }
// ======================================================================
public:
// ======================================================================
Double_t evaluate () const override ;
// ======================================================================
private:
// ======================================================================
RooRealProxy m_xvar {} ;
RooListProxy m_pars {} ;
mutable Ostap::Math::Rational m_rational { 3 , 1 } ;
// ======================================================================
} ;
// ========================================================================
/** @class RationalBernstein
* Rational fnuction as ratio of Bernstein polyhomial and
* positive Bernstein polynomial
* \f[ R ( x ) = \frac{B(x)}{P(x)\f]
* @see Ostap::Math::RationalBernstein
* @see Ostap::Math::Bernstein
* @see Ostap::Math::Positive
* @author Vanya BELYAEV [email protected]
* @date 2023-09-21
*/
class RationalBernstein final : public RooAbsReal
{
public:
// ======================================================================
ClassDefOverride(Ostap::MoreRooFit::RationalBernstein, 1) ;
// ======================================================================
public:
// ======================================================================
RationalBernstein
( const std::string& name ,
const std::string& title ,
RooAbsReal& xvar ,
const RooArgList& p , // numerator
const RooArgList& q , // denominator
const double xmin ,
const double xmax ) ;
RationalBernstein
( const std::string& name ,
const std::string& title ,
RooAbsReal& xvar ,
const RooArgList& pars , // all pars
const unsigned short p , // degree of numerator
const double xmin ,
const double xmax ) ;
/// copy constructor
RationalBernstein ( const RationalBernstein& right , const char* name = nullptr ) ;
/// default
RationalBernstein () ;
/// clone method
RationalBernstein* clone ( const char* name ) const override ;
/// virtual destructor
virtual ~RationalBernstein() ;
// ======================================================================
public:
// ======================================================================
Int_t getAnalyticalIntegral
( RooArgSet& allVars ,
RooArgSet& analVars ,
const char* rangeName = nullptr ) const override ;
Double_t analyticalIntegral
( Int_t code ,
const char* rangeName = nullptr ) const override ;
// ======================================================================
public:
// ======================================================================
/// 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 () ; }
/// xmin for Rational
double xmin () const { return m_rational.xmin () ; }
/// xmax for Rational
double xmax () const { return m_rational.xmax () ; }
// ======================================================================
public:
// ======================================================================
void setPars () const ;
// ======================================================================
public:
// ======================================================================
const Ostap::Math::RationalBernstein& function () const { return m_rational ; }
const Ostap::Math::RationalBernstein& rational () const { return m_rational ; }
// ======================================================================
public:
// ======================================================================
Double_t evaluate () const override ;
// ======================================================================
private:
// ======================================================================
RooRealProxy m_xvar {} ;
RooListProxy m_pars {} ;
mutable Ostap::Math::RationalBernstein m_rational { 3 , 1 } ;
// ======================================================================
} ;
// ========================================================================
/** @class Shape1D
* The generic "fixed shape" function
* @author Vanya BELYAEV [email protected]
* @date 2023-01-30
Expand Down
2 changes: 1 addition & 1 deletion source/include/Ostap/Rational.h
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ namespace Ostap
* the Floater-Hormann rational barycentric interpolant
* and parameters are the function valeus at Chebyshev's nodes
*
# @see Ostap::Math::FloaterHormann
* @see Ostap::Math::FloaterHormann
* @author Vanya BELYAEV [email protected]
* @date 2023-09-14
*/
Expand Down
Loading

0 comments on commit 624ec49

Please sign in to comment.