Skip to content

Commit

Permalink
?
Browse files Browse the repository at this point in the history
  • Loading branch information
VanyaBelyaev committed Apr 5, 2024
1 parent f55fd32 commit 1f5d07a
Show file tree
Hide file tree
Showing 5 changed files with 179 additions and 34 deletions.
2 changes: 1 addition & 1 deletion .aux/test_with_lcg
Original file line number Diff line number Diff line change
Expand Up @@ -4,5 +4,5 @@ CMTCONFIG=$2
source /cvmfs/sft.cern.ch/lcg/views/${LCG}/${CMTCONFIG}/setup.sh
source build/INSTALL/thisostap.sh
cd build
ctest -N && cmake .. -DCMAKE_INSTALL_PREFIX=./INSTALL/ && ctest -j6 -R '(linalg|pypdf|morph|fill|eff)' --output-on-failure
ctest -N && cmake .. -DCMAKE_INSTALL_PREFIX=./INSTALL/ && ctest -j6 -R '(linalg|pypdf|morph|fill|eff|model)' --output-on-failure

113 changes: 90 additions & 23 deletions ostap/fitting/roofitresult.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,12 @@
from ostap.core.meta_info import root_info
from ostap.core.core import Ostap, VE, valid_pointer, iszero, isequal
from ostap.core.ostap_types import string_types , integer_types
import ostap.math.linalg
import ostap.fitting.variables
import ostap.fitting.printable
import ostap.math.linalg as LA
from ostap.logger.colorized import allright, attention
from ostap.logger.utils import pretty_float, pretty_ve, pretty_2ve
import ROOT, math, sys
import ROOT, math, sys, ctypes
# =============================================================================
from ostap.logger.logger import getLogger
if '__main__' == __name__ : logger = getLogger ( 'ostap.fitting.roofitresult' )
Expand Down Expand Up @@ -205,16 +205,19 @@ def _rfr_cov_matrix_ ( self , var1 , var2 , *vars ) :
## get the covariance matrix
# @author Vanya BELYAEV [email protected]
# @date 2011-06-07
def _rfr_covmatrix_ ( self ) :
def _rfr_covmatrix_ ( self , tmatrix = False ) :
"""Get covariance ( matrix
>>> result = ...
>>> cov = results.covmatrix()
>>> print corr
"""

cm = self.covarianceMatrix ()
N = cm.GetNrows()

if tmatrix : return cm ## RETURN NATIVE ROOT TMAtrixSym<double

N = cm.GetNrows()

import ostap.math.linalg
m = Ostap.Math.SymMatrix ( N )()

Expand Down Expand Up @@ -621,7 +624,7 @@ def _rfr_table_ ( rr , title = '' , prefix = '' , more_vars = {} ) :
else : n = ''

rows.append ( ( 'Estimated distance to minimum' , n , ' ' + s , '' ) )

cq = r.covQual()
cn = ''
if -1 == cq :
Expand Down Expand Up @@ -652,9 +655,9 @@ def _rfr_table_ ( rr , title = '' , prefix = '' , more_vars = {} ) :

with_globcorr = not ( (6,24) <= root_info < (6,28) )
with_globcorr = True or not ( (6,24) <= root_info < (6,28) )
with_globcorr = True
with_globcorr = not ( (6,24) <= root_info < (6,26) )

with_globcorr = True

if with_globcorr : rows = [ ( '', 'Unit', 'Value' , 'Global/max correlation [%]') ] + rows
else : rows = [ ( '', 'Unit', 'Value' , 'Max correlation [%]') ] + rows
Expand Down Expand Up @@ -778,6 +781,47 @@ def _rfr_print_ ( self , opts = 'v' ) :
result = result.encode ('utf-8')
return result

# =============================================================================
## Get global correlation coeffcient
# @code
# fit_reuslt = ///
# coeff = fit_results.global_cc ( 3 )
# @endcode
# @see Ostap::Utils::global_cc
def _rfr_global_cc_ ( rfr , par ) :
"""Get global correlation coeffcient
- see `Ostap.Utils.global_cc `
>>> fit_reuslt = ///
>>> coeff = fit_results.global_cc ( 3 )
"""
if par is None :

## get all values at once
values = Ostap.Utils.global_cc ( rfr )
if len ( values ) != len ( rfr ) :
raise ValueError( "Error from Ostap::Utils::global_cc" )
return tuple ( v for v in values )

elif isinstance ( par , string_types ) :

index = rfr.floatParsFinal().index ( par )
if index < 0 : raise KeyError ( 'Invalid parameter %s' % par )
return _rfr_global_cc_ ( rfr , index )

elif isinstance ( par , ROOT.RooAbsArg ) :
return _rfr_global_cc_ ( rfr , par.name )

index = par
assert isinstance ( index , integer_types ) and 0 <= index , \
'Invald index %s/%s' % ( par , index )

if not 0 <= index < len ( rfr ) :
raise IndexError("Invalid index %s/%s" % ( par , index ) )

return Ostap.Utils.global_cc ( rfr , index )

# =============================================================================
## Get global correlation coefficient for the parameter
# \f$ \rho_k = \sqrt{ 1 - \left[ C_{kk} V_{kk}\right]^{-1} } \f$
Expand All @@ -789,7 +833,7 @@ def _rfr_print_ ( self , opts = 'v' ) :
# It should be accessible via <code>RooFitResult::globalCorr</code> method
# but often it results in segfault
# @see RooFitResult::globalCorr
def _rfr_global_corr_ ( self , par ) :
def _rfr_global_corr_ ( rfr , par ) :
""" Get global correlation coefficient for the parameter
rho_k = sqrt ( 1 - 1/( C_{kk} V_{kk}) )
- where C is covariance matrix and V is inverse
Expand All @@ -799,31 +843,52 @@ def _rfr_global_corr_ ( self , par ) :
- see `RooFitResult.globalCorr`
"""
if isinstance ( par , string_types ) :
index = self.floatParsFinal().index ( par )
if index < 0 : return -1
index = rfr.floatParsFinal().index ( par )
if index < 0 : return -1
elif isinstance ( par , integer_types ) :
if not 0 <= par < len ( self.floatParsFinal () ) : return -1
if not 0 <= par < len ( rfr.floatParsFinal () ) : return -1
index = par
elif isinstance ( par , ROOT.RooAbsArg ) :
return _rfr_global_corr_ ( self , par.name )
return _rfr_global_corr_ ( rfr , par.name )
else : return -1

np = len ( rfr )

## get the covariance matrix and invert it!
v = self.covmatrix()
if not v.InvertChol () : return -1
return _rfr_global_cc_ (rfr , par )

## get the covariance matrix
c = self.covmatrix()
## ## ostap.math.linalg2 machinery fails for large matrices
## if 10 >= np or ( np , 'double' ) in LA.LinAlgT.known_ssymmatrices :

## ## get the covariance matrix and invert it
## mc = self.covmatrix ()
## if not mc.InvertChol () : return -1

## ## get the covariance matrix
## cm = self.covmatrix()

## else : ## use TMatrixTSym


## cm = self.covarianceMatrix()
## mc = type ( cm ) ( cm )
## det = ctypes.c_double(0)
## mc.Invert( det )
## det = det.value
## if iszero ( det ) : return -1


## cv = cm ( index , index ) * mc ( index , index )
## rho2 = 1.0 - 1.0 / cv
## if rho2 < 0 : return -1
## return math.sqrt ( rho2 )

cv = c ( index , index ) * v ( index , index )

rho2 = 1.0 - 1.0 / cv
if rho2 < 0 : return -1.0

return math.sqrt ( rho2 )

# =============================================================================
##
## The size of the problem == number of floating parameters
def _rfr_len_ ( rfr ) :
""" The size of the problem == number of floating parameters
"""
return len ( rfr.floatParsFinal() )

# =============================================================================
## Symmetrized Kullback-Leibler divergency between two RooFitResult objects
Expand Down Expand Up @@ -1051,6 +1116,8 @@ def _rm_contour_ ( self ,
ROOT.RooFitResult . global_corr = _rfr_global_corr_
ROOT.RooFitResult . kullback = _rfr_kullback_
ROOT.RooFitResult . kullback_leibler = _rfr_kullback_
ROOT.RooFitResult . __len__ = _rfr_len_
ROOT.RooFitResult . global_cc = _rfr_global_cc_

_new_methods_ += [
ROOT.RooFitResult . __repr__ ,
Expand Down
7 changes: 4 additions & 3 deletions ostap/math/linalg.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,8 @@
# =============================================================================
__all__ = (
'mgetter' , ## get (i,j) element from matrix-like object
'checkops' , ## check the allowed operations
'checkops' , ## check the allowed operations
'LinAlgT' , ## LinAlgenra type&decorator store
)
# =============================================================================
import ROOT
Expand All @@ -26,8 +27,8 @@
if '__main__' == __name__ : logger = getLogger ( 'ostap.math.linalg2' )
else : logger = getLogger ( __name__ )
# =============================================================================
from ostap.math.linalg2 import mgetter, checkops
import ostap.math.linalgt
from ostap.math.linalg2 import mgetter, checkops
from ostap.math.linalgt import LinAlgT

# =============================================================================
if '__main__' == __name__ :
Expand Down
22 changes: 22 additions & 0 deletions source/include/Ostap/FitResult.h
Original file line number Diff line number Diff line change
Expand Up @@ -97,6 +97,28 @@ namespace Ostap
// ======================================================================
}; // The end of the class Ostap::FitResults
// ========================================================================
/** calculate the global correlation coefficient
* \f$ \rho_k = \sqrt{ 1 - \left[ C_{kk} V_{kk}\right]^{-1} } \f$
* where \f$ C \f$ is covarinace matrix and \f$ V = C^{-1}\$ is inverse
* @code
* @param index parameter index
* @param r fit result
* @return global correlation coefficients (or -1 in case of failure)
*/
double global_cc
( const RooFitResult& r ,
const unsigned short index ) ;
// ========================================================================
/** calculate the global correlation coefficients
* \f$ \rho_k = \sqrt{ 1 - \left[ C_{kk} V_{kk}\right]^{-1} } \f$
* where \f$ C \f$ is covarinace matrix and \f$ V = C^{-1}\$ is inverse
* @code
* @param r fit result
* @return vector global correlation coefficients (empty in case of failure)
*/
std::vector<double> global_cc
( const RooFitResult& r ) ;
// ========================================================================
} // The end of namespace Ostap::Utils
// ==========================================================================
} // The end of namespace Ostap
Expand Down
69 changes: 62 additions & 7 deletions source/src/FitResult.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
// ============================================================================
// Ostap
// ============================================================================
#include "Ostap/Math.h"
#include "Ostap/FitResult.h"
// ============================================================================
/** @file
Expand Down Expand Up @@ -119,10 +120,15 @@ Ostap::Utils::FitResults::clone () const { return Clone() ; }
std::vector<double>
Ostap::Utils::FitResults::global_cc () const
{
if ( nullptr == _GC ) { return std::vector<double>() ; }
return std::vector<double>( _GC->GetMatrixArray() ,
_GC->GetMatrixArray() +
_GC->GetNrows() ) ;
// use (pre)calcualted values
if ( nullptr != _GC )
{
return std::vector<double>( _GC->GetMatrixArray () ,
_GC->GetMatrixArray () +
_GC->GetNrows () ) ;
}
// recalculate:
return Ostap::Utils::global_cc ( *this ) ;
}
// ============================================================================
// add label/status pair to the history
Expand All @@ -132,9 +138,58 @@ void Ostap::Utils::FitResults::add_to_history
const int status )
{ _statusHistory.push_back ( std::make_pair ( label , status ) ) ; }
// ============================================================================



/* calculate the global correlation coefficients
* \f$ \rho_k = \sqrt{ 1 - \left[ C_{kk} V_{kk}\right]^{-1} } \f$
* where \f$ C \f$ is covarinace matrix and \f$ V = C^{-1}\$ is inverse
* @code
* @param r fit result
* @return vector global correlation coefficients (empty in case of failure)
*/
// ============================================================================
std::vector<double>
Ostap::Utils::global_cc ( const RooFitResult& r )
{
/// zero for doubles
const static Ostap::Math::Zero<double> s_zero {} ; // zero for doubles
///
/// 1) try to invert covariance matrix
TMatrixTSym<double> cinv { r.covarianceMatrix() } ;
double det = 0 ;
cinv.Invert ( &det ) ;
if ( s_zero ( det ) ) { return std::vector<double>() ; }
///
const int ncols = cinv.GetNcols() ;
std::vector<double> result ( ncols, 0.0 ) ;
//
const TMatrixTSym<double>& cm = r.covarianceMatrix() ;
for ( int index = 0 ; index < ncols ; ++index )
{
const double cv = cinv ( index , index ) * cm ( index , index ) ;
if ( s_zero ( cv ) ) { return std::vector<double>() ; }
const double rho2 = 1.0 - 1.0 / cv ;
if ( rho2 < 0 ) { return std::vector<double>() ; }
result [ index ] = std::sqrt ( rho2 ) ;
}
//
return result ;
}
// ============================================================================
/* calculate the global correlation coefficient
* \f$ \rho_k = \sqrt{ 1 - \left[ C_{kk} V_{kk}\right]^{-1} } \f$
* where \f$ C \f$ is covarinace matrix and \f$ V = C^{-1}\$ is inverse
* @code
* @param r fit result
* @return global correlation coefficients (or -1 in case of failure)
*/
// ============================================================================
double Ostap::Utils::global_cc
( const RooFitResult& r ,
const unsigned short index )
{
const std::vector<double> result { Ostap::Utils::global_cc ( r ) } ;
if ( result.size() <= index ) { return -1 ; }
return result [ index ] ;
}
// ============================================================================
ClassImp(Ostap::Utils::FitResults )
// ============================================================================
Expand Down

0 comments on commit 1f5d07a

Please sign in to comment.