From 589b6bacdf217b585e32c448592c0a774b873dcc Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Sun, 1 Sep 2024 17:34:14 +0200 Subject: [PATCH] 1.improve all toys machinery: toys, jackknife, bootstrapping and significance 1. extend all tests for toys --- ReleaseNotes/release_notes.md | 6 +- ostap/fitting/tests/test_fitting_toys.py | 398 +++++++++++++++++---- ostap/fitting/toys.py | 396 ++++++++++++++++---- ostap/parallel/parallel_toys.py | 74 +++- ostap/parallel/tests/test_parallel_toys.py | 391 ++++++++++++++++---- 5 files changed, 1047 insertions(+), 218 deletions(-) diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index 85b09d6d..aff26078 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -1,5 +1,9 @@ ## New features - + + 1.improve all toys machinery: toys, jackknife, bootstrapping and significance + 1. extend all tests for toys + + ## Backward incompatible ## Bug fixes diff --git a/ostap/fitting/tests/test_fitting_toys.py b/ostap/fitting/tests/test_fitting_toys.py index 6b12c6be..5ae572fd 100755 --- a/ostap/fitting/tests/test_fitting_toys.py +++ b/ostap/fitting/tests/test_fitting_toys.py @@ -16,14 +16,18 @@ __author__ = "Ostap developers" __all__ = () ## nothing to import # ============================================================================= -from ostap.core.pyrouts import hID +from ostap.core.pyrouts import hID, VE import ostap.fitting.models as Models import ostap.fitting.toys as Toys import ostap.histos.histos -from ostap.utils.timing import timing +from ostap.utils.timing import timing +from ostap.utils.memory import memory from ostap.plotting.canvas import use_canvas from ostap.utils.utils import wait from ostap.fitting.toys import pull_var +from ostap.math.base import num_range +from ostap.utils.progress_bar import progress_bar +import ostap.logger.table as T import ROOT, time, random # ============================================================================= # logging @@ -34,16 +38,36 @@ else : logger = getLogger ( __name__ ) # ============================================================================= -nominal_mean = 0.4 -nominal_sigma = 0.1 +nominal_mean = 0.4 +nominal_sigma = 0.1 +nominal_S = 1000 +nominal_B = 100 +nominal_F = float ( nominal_S ) / ( nominal_S + nominal_B ) -mass = ROOT.RooRealVar ( 'mass' , '', 0 , 1 ) -gen_gauss = Models.Gauss_pdf ( 'GG' , xvar = mass ) -fit_gauss = Models.Gauss_pdf ( 'FG' , xvar = mass ) +mass = ROOT.RooRealVar ( 'mass' , '', 0 , 1 ) +gen_gauss = Models.Gauss_pdf ( 'GG' , xvar = mass ) +fit_gauss = Models.Gauss_pdf ( 'FG' , xvar = mass ) gen_gauss.mean = nominal_mean gen_gauss.sigma = nominal_sigma -model = Models.Fit1D ( signal = gen_gauss , background = 'flat' ) +model = Models.Fit1D ( signal = gen_gauss , background = 'flat' ) +model_NE = Models.Fit1D ( signal = gen_gauss , background = 'flat' , extended = False ) +model.S = nominal_S +model.B = nominal_B +model_NE.F = nominal_F + +fit_config = { 'silent' : True , 'refit' : 5 } + +def reset_funcs () : + + gen_gauss.mean = nominal_mean + gen_gauss.sigma = nominal_sigma + fit_gauss.mean = nominal_mean + fit_gauss.sigma = nominal_sigma + model.S = nominal_S + model.B = nominal_B + +toy_results = {} # ============================================================================== ## Perform toy-study for possible fit bias and correct uncertainty evaluation @@ -54,7 +78,7 @@ # - fill distributions for fit results # - fill distribution of pulls def test_toys ( ) : - """Perform toys-study for possible fit bias and correct uncertainty evaluation + """ Perform toys-study for possible fit bias and correct uncertainty evaluation - generate `nToys` pseudoexperiments with some PDF `pdf` - fit teach experiment with the same PDF - store fit results @@ -67,14 +91,19 @@ def test_toys ( ) : more_vars = { 'pull:mean_GG' : pull_var ( 'mean_GG' , nominal_mean ) , 'pull:sigma_GG' : pull_var ( 'sigma_GG' , nominal_sigma ) } +\ + ## reset all functions + reset_funcs () + + pdf = model with timing ( 'Toys analysis' , logger = logger ) : results , stats = Toys.make_toys ( - pdf = gen_gauss , + pdf = pdf , nToys = 1000 , data = [ mass ] , - gen_config = { 'nEvents' : 200 , 'sample' : True } , - fit_config = { 'silent' : True , 'refit' : 5 } , + gen_config = { 'nEvents' : 1100 , 'sample' : True } , + fit_config = fit_config , init_pars = { 'mean_GG' : nominal_mean , 'sigma_GG' : nominal_sigma } , more_vars = more_vars , @@ -91,6 +120,8 @@ def test_toys ( ) : for h in ( h_mean , h_sigma ) : with use_canvas ( 'test_toys %s' % h.title , wait = 1 ) : h.draw() + + toy_results [ 'test_toys' ] = results , stats # ============================================================================= ## Perform toy-study for possible fit bias and correct uncertainty evaluation @@ -111,6 +142,9 @@ def test_toys2 ( ) : more_vars = { 'pull:mean_FG' : pull_var ( 'mean_FG' , nominal_mean ) , 'pull:sigma_FG' : pull_var ( 'sigma_FG', nominal_sigma ) } + ## reset all functions + reset_funcs () + with timing ( 'Toys2 analysis' , logger = logger ) : results , stats = Toys.make_toys2 ( gen_pdf = gen_gauss , @@ -118,7 +152,7 @@ def test_toys2 ( ) : nToys = 1000 , data = [ mass ] , gen_config = { 'nEvents' : 200 , 'sample' : True } , - fit_config = { 'silent' : True } , + fit_config = fit_config , gen_pars = { 'mean_GG' : nominal_mean , 'sigma_GG' : nominal_sigma } , fit_pars = { 'mean_FG' : nominal_mean , 'sigma_FG' : nominal_sigma } , more_vars = more_vars , @@ -137,49 +171,197 @@ def test_toys2 ( ) : for h in ( h_mean , h_sigma ) : with use_canvas ( 'test_toys2 %s' % h.title , wait = 1 ) : h.draw() + toy_results [ 'test_toys2' ] = results , stats # ============================================================================== -## Perform toy-study for Jackknife -def test_jackknife ( ) : - """Perform toys-study for Jackknife +## Perform toy-study for Jackknife (non-extended model) +def test_jackknife_NE1 ( ) : + """ Perform toys-study for Jackknife (non-extended model) """ - logger = getLogger ( 'test_parallel_jackknife' ) + logger = getLogger ( 'test_parallel_jackknife_NE1' ) + logger.info ( 'Jackknife analysis for non-extended model' ) + + pdf = model_NE - with timing ( 'Jackknife analysis' , logger = logger ) : - model.S = 1000 - model.B = 100 + with timing ( 'Jackknife analysis (non-extended)' , logger = logger ) : - data = model.generate ( 1100 ) + ## reset all functions + reset_funcs () - Toys.make_jackknife ( pdf = model , - data = data , - fit_config = { 'silent' : True } , - progress = True , - silent = True ) - + data = pdf.generate ( nominal_S ) + with use_canvas ( "Jackknife analysis (non-extended)" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + ## reset all functions + reset_funcs () + + results, stats = Toys.make_jackknife ( pdf = pdf , + data = data , + fit_config = fit_config , + progress = True , + silent = True ) + + toy_results [ 'test_jackknife_NE1' ] = results , stats + # ============================================================================== -## Perform toy-study for Bootstrap -def test_bootstrap ( ) : - """Perform toys-study for Bootstrap +## Perform toy-study for Jackknife (non-extended model) +def test_jackknife_NE2 ( ) : + """ Perform toys-study for Jackknife (non-extended model) """ - logger = getLogger ( 'test_parallel_bootstrap' ) - - with timing ( 'Bootstrap analysis' , logger = logger ) : + logger = getLogger ( 'test_parallel_jackknife_NE2' ) + logger.info ( 'Jackknife analysis for non-extended model' ) + + pdf = fit_gauss + + with timing ( 'Jackknife analysis (non-extended)' , logger = logger ) : + + ## reset all functions + reset_funcs () + + data = pdf.generate ( nominal_S ) + with use_canvas ( "Jackknife analysis (non-extended)" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + ## reset all functions + reset_funcs () + + results, stats = Toys.make_jackknife ( pdf = pdf , + data = data , + fit_config = fit_config , + progress = True , + silent = True ) + + toy_results [ 'test_jackknife_NE2' ] = results , stats + +# ============================================================================== +## Perform toy-study for Jackknife (extended model) +def test_jackknife_EXT ( ) : + """ Perform toys-study for Jackknife (extended model) + """ - model.S = 1000 - model.B = 100 + logger = getLogger ( 'test_jackknife_EXT' ) + logger.info ( 'Jackknife analysis for extended model' ) + + pdf = model + + with timing ( 'Jackknife analysis (extended)' , logger = logger ) : - data = model.generate ( 1100 ) + ## reset all functions + reset_funcs () - Toys.make_bootstrap ( pdf = model , - size = 100 , - data = data , - fit_config = { 'silent' : True } , - progress = True , - silent = True ) + data = pdf.generate ( nominal_S ) + with use_canvas ( "Jackknife analysis (extended)" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + ## reset all functions + reset_funcs () + + results, stats = Toys.make_jackknife ( pdf = pdf , + data = data , + fit_config = fit_config , + progress = True , + silent = True ) + + toy_results [ 'test_jackknife_EXT' ] = results , stats + +# ============================================================================== +## Perform toy-study for Bootstrap (non-extended) +def test_bootstrap_NE1 ( ) : + """ Perform toys-study for Bootstrap (non-extended) + """ + + logger = getLogger ( 'test_bootstrap_NE1' ) + logger.info ( 'Bootstrap analysis for non-extended model' ) + + pdf = model_NE + with timing ( 'Bootstrap analysis (non-extended)' , logger = logger ) : + + ## reset all functions + reset_funcs () + + data = pdf.generate ( nominal_S + nominal_B ) + with use_canvas ( "Booststrap analysis" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + ## reset all functions + reset_funcs () + results, stats = Toys.make_bootstrap ( pdf = pdf , + size = 400 , + data = data , + fit_config = fit_config , + extended = False , + progress = True , + silent = True ) + + toy_results [ 'test_bootstrap_NE1' ] = results , stats + +# ============================================================================== +## Perform toy-study for Bootstrap (non-extended) +def test_bootstrap_NE2 ( ) : + """ Perform toys-study for Bootstrap (non-extended) + """ + + logger = getLogger ( 'test_parallel_bootstrap_NE2' ) + logger.info ( 'Bootstrap analysis for non-extended model' ) + + pdf = fit_gauss + with timing ( 'Bootstrap analysis (non-extended)' , logger = logger ) : + + ## reset all functions + reset_funcs () + + data = pdf.generate ( nominal_S + nominal_B ) + with use_canvas ( "Booststrap analysis" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + ## reset all functions + reset_funcs () + + results, stats = Toys.make_bootstrap ( pdf = pdf , + size = 400 , + data = data , + fit_config = fit_config , + extended = False , + progress = True , + silent = True ) + + toy_results [ 'test_bootstrap_NE2' ] = results , stats + +# ============================================================================== +## Perform toy-study for Bootstrap (extended) +def test_bootstrap_EXT ( ) : + """ Perform toys-study for Bootstrap (extended) + """ + + logger = getLogger ( 'test_parallel_bootstrap_EXT' ) + logger.info ( 'Bootstrap analysis for extended model' ) + + pdf = model + with timing ( 'Bootstrap analysis (extended)' , logger = logger ) : + + ## reset all functions + reset_funcs () + + data = pdf.generate ( nominal_S + nominal_B ) + with use_canvas ( "(Extended) Booststrap analysis" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + ## reset all functions + reset_funcs () + + results, stats = Toys.make_bootstrap ( pdf = pdf , + size = 400 , + data = data , + fit_config = fit_config , + extended = True , + progress = True , + silent = True ) + + toy_results [ 'test_bootstrap_EXT' ] = results , stats + # ============================================================================= ## Perform toy-study for significance of the signal # - generate nToys pseudoexperiments using background-only hypothesis @@ -187,7 +369,7 @@ def test_bootstrap ( ) : # - store fit results # - fill distributions for fit results def test_significance ( ) : - """Perform toy-study for significance of the signal + """ Perform toy-study for significance of the signal - generate `nToys` pseudoexperiments using background-only hypothesis - fit each experiment with signal+background hypothesis - store fit results @@ -196,6 +378,9 @@ def test_significance ( ) : logger = getLogger ( 'test_significance' ) + ## reset all functions + reset_funcs () + with timing ( 'Significance analysis' , logger = logger ) : ## only background hypothesis @@ -209,35 +394,120 @@ def test_significance ( ) : model.background.tau.fix ( 0 ) results , stats = Toys.make_toys2 ( - gen_pdf = bkg_only , - fit_pdf = model , - nToys = 1000 , - data = [ mass ] , - gen_config = { 'nEvents' : 100 , 'sample' : True } , - fit_config = { 'silent' : True } , + gen_pdf = bkg_only , + fit_pdf = model , + nToys = 2000 , + data = [ mass ] , + gen_config = { 'nEvents' : 110 , 'sample' : True } , + fit_config = fit_config , gen_pars = { 'tau_BKG' : 0. } , ## initial values for generation - fit_pars = { 'B' : 100 , - 'S' : 10 , - 'phi0_Bkg_S': 0.0 } , ## initial fit values for parameters + fit_pars = { 'B' : 100 , + 'S' : 10 , + 'phi0_Bkg_S' : 0.0 } , ## initial fit values for parameters silent = True , progress = True ) - - h_S = ROOT.TH1F ( hID() , '#S' , 60 , 0 , 60 ) - - for r in results ['S' ] : h_S .Fill ( r ) + + + # ========================================================================= + ## yields + h_Y = ROOT.TH1F ( hID() , '#S' , 140 , 0 , 70 ) + for r in results [ 'S' ] : h_Y .Fill ( r ) + + ## get p-value and significance histograms: + h_P , h_S = h_Y.significance ( cut_low = True ) + + h_S.red () + h_P.blue () + + minv = 1.e+100 + for i,_,y in h_P.items() : + yv = y.value() + if 0 < yv and yv < minv : minv = yv + minv , _ = num_range ( 0.75 * minv ) + h_P.SetMinimum ( minv ) + h_S.SetMinimum ( 0 ) - for h in ( h_S , ) : - with use_canvas ( 'test_significance %s' % h.title , wait = 1 ) : h.draw() - + with use_canvas ( 'test_significance: yields' , wait = 1 ) : h_Y.draw ( ) + with use_canvas ( 'test_significance: p-value' , wait = 1 ) : h_P.draw ( logy = True ) + with use_canvas ( 'test_significance: significance' , wait = 3 ) : h_S.draw () + + toy_results [ 'test_significance' ] = results , stats +# ============================================================================= +## Save resuls of toys to DBASE +def test_db () : + """ Save resuls of toys to DBASE + """ + logger = getLogger ( 'test_db' ) + logger.info ( 'Saving all toys results into DBASE' ) + import ostap.io.zipshelve as DBASE + from ostap.utils.timing import timing + with timing( 'Save everything to DBASE', logger ), DBASE.tmpdb() as db : + for key in progress_bar ( toy_results ) : + r , s = toy_results [ key ] + key1 = '%s:results' % key + key2 = '%s:stats' % key + db [ key1 ] = r + db [ key2 ] = s + db.ls() + while toy_results : toy_results.popitem() +# ============================================================================= + + # ============================================================================= if '__main__' == __name__ : + + """ + with memory ( 'test_toys' , logger = logger ) as mtt : test_toys () + with memory ( 'test_toys2' , logger = logger ) as mt2 : test_toys () + with memory ( 'test_jackknife_NE1' , logger = logger ) as mj1 : test_jackknife_NE1 () + with memory ( 'test_jackknife_NE2' , logger = logger ) as mj2 : test_jackknife_NE2 () + with memory ( 'test_jackknife_EXT' , logger = logger ) as mje : test_jackknife_EXT () + with memory ( 'test_bootstrap_NE1' , logger = logger ) as mb1 : test_bootstrap_NE1 () + with memory ( 'test_bootstrap_NE2' , logger = logger ) as mb2 : test_bootstrap_NE2 () + with memory ( 'test_bootstrap_EXT' , logger = logger ) as mbe : test_bootstrap_EXT () + """ + + with memory ( 'test_significance' , logger = logger ) as mts : test_significance () + with memory ( 'test_db' , logger = logger ) as mdb : test_db () + + rows = [ ( 'Test' , 'Memory [MB]' ) ] + + """ + row = 'Toys' , '%.1f' % mtt.delta + rows.append ( row ) + + row = 'Toys2' , '%.1f' % mt2.delta + rows.append ( row ) + + row = 'Jackknife_NE1' , '%.1f' % mj1.delta + rows.append ( row ) + + row = 'Jackknife_NE2' , '%.1f' % mj2.delta + rows.append ( row ) + + row = 'Jackknife_EXT' , '%.1f' % mje.delta + rows.append ( row ) + + row = 'Bootstrap_NE1' , '%.1f' % mb1.delta + rows.append ( row ) + + row = 'Bootstrap_NE2' , '%.1f' % mb2.delta + rows.append ( row ) + + row = 'Bootstrap_EXT' , '%.1f' % mbe.delta + rows.append ( row ) + """ + + row = 'Significance' , '%.1f' % mts.delta + rows.append ( row ) + + row = 'test_db' , '%.1f' % mdb.delta + rows.append ( row ) - test_toys ( ) - test_toys2 ( ) - test_jackknife ( ) - test_bootstrap ( ) - test_significance ( ) + title = 'Memory usage for various toys' + table = T.table ( rows , title = title , prefix = '# ' , alignment = 'lc' ) + logger.info ( '%s:\n%s' % ( title , table ) ) # ============================================================================= diff --git a/ostap/fitting/toys.py b/ostap/fitting/toys.py index d0401dec..07c52cd1 100644 --- a/ostap/fitting/toys.py +++ b/ostap/fitting/toys.py @@ -26,6 +26,8 @@ from builtins import range from ostap.core.ostap_types import string_types, integer_types from ostap.core.core import VE, SE +from ostap.logger.pretty import pretty_ve, pretty_float, fmt_pretty_float +from ostap.logger.colorized import attention import ROOT # ============================================================================= # logging @@ -34,11 +36,36 @@ if '__main__' == __name__ : logger = getLogger( 'ostap.fitting.toys' ) else : logger = getLogger( __name__ ) # ============================================================================= -logger.debug ( 'Utilities to run fitting toys, Jackknife, bootstrap, ... ') +logger.debug ( 'Utilities to run toys, Jackknife, bootstrap, ... ') # ============================================================================== +extmodes = { + ROOT.RooAbsPdf.MustBeExtended : 'MustBeExtended' , + ROOT.RooAbsPdf.CanBeExtended : 'CanBeExtended' , + ROOT.RooAbsPdf.CanNotBeExtended : 'CanBeExtended' , +} +# ============================================================================== +## Check the consistency of the PDF and `extended` flag +def check_extended ( pdf , extended ) : + """ Check the consistency of the PDF and `extended` flag + """ + if isinstance ( pdf , ROOT.RooAbsPdf ) : pass + else : pdf = pdf.pdf + ## + extmode = pdf.extendMode() + ## + if pdf.canBeExtended () : return True + elif pdf.mustBeExtended () and extended : return True + elif ROOT.RooAbsPdf.MustBeExtended == extmode and extended : return True + elif ROOT.RooAbsPdf.CanNotBeExtended == extmode and not extended : return True + elif ROOT.RooAbsPdf.CanBeExtended == extmode and extended : return True + ## + logger.warning ( "Mismatch of `ExtendMode`:`%s` and `extended`:%s" % ( extmodes [ extmode ] , extended ) ) + ## + return False +# ============================================================================= ## Technical transformation to the dictionary : { 'name' : float_value } def vars_transform ( vars ) : - """Technical transformation to the dictionary : `{ 'name' : float_value }` + """ Technical transformation to the dictionary : `{ 'name' : float_value }` """ from ostap.core.ostap_types import dictlike_types @@ -76,7 +103,7 @@ def pull_var ( name , value ) : return PoolVar ( name , value ) # statistics = make_stats ( results ) # @endcode def make_stats ( results , fits = {} , covs = {} , accept = SE () ) : - """Prepare statistics from results of toys/jackknifes/boostrap studies + """ Prepare statistics from results of toys/jackknifes/boostrap studies >>> results = ... >>> statistics = make_stats ( results ) """ @@ -98,25 +125,55 @@ def make_stats ( results , fits = {} , covs = {} , accept = SE () ) : for k in covs : stats ['- CovQual %s' % k ] = covs [ k ] - return stats # ============================================================================= -## print statistics of pseudoexperiments +## Print statistics of pseudoexperiments def print_stats ( stats , ntoys = '' , logger = logger ) : - """print statistics of pseudoexperiments + """ Print statistics of pseudoexperiments """ - - table = [ ( 'Parameter' , '#', 'mean' , 'rms' , '%13s / %-13s' % ( 'min' , 'max' ) ) ] + ## 0 1 2 3 4 5 6 7 + table = [ ( 'Parameter' , '#', 'mean', 'x[..]' , 'rms' , 'x[..]', 'min / max' , 'x[..]' ) ] - def make_row ( c ) : - n = "{:^11}".format ( c.nEntries() ) + has_e1 = [] + has_e2 = [] + has_e3 = [] + + def make_row ( c ) : + + n = '%d' % c.nEntries() mean = c.mean () - mean = "%+13.6g +/- %-13.6g" % ( mean.value() , mean.error () ) - rms = "%13.6g" % c.rms () - minmax = "%+13.6g / %-+13.6g" % ( c.min() , c.max () ) - return p , n , mean , rms , minmax + rms = c.rms () + minmax = c.min() , c.max() + + mean , expo1 = pretty_ve ( mean , precision = 6 , width = 8 , parentheses = False ) + rms , expo2 = pretty_float ( rms , precision = 4 , width = 6 ) + fmt , expo3 = fmt_pretty_float ( max ( abs ( c.min() ) , abs ( c.max() ) ) , precision = 3 , width = 5 ) + fmt = '%s / %s ' % ( fmt , fmt ) + + if expo3 : scale = 10**expo3 + else : scale = 1 + minmax = fmt % ( c.min() / scale , c.max() / scale ) + + ## mean = "%+13.6g +/- %-13.6g" % ( mean.value() , mean.error () ) + ## rms = "%13.6g" % c.rms () + ## minmax = "%+13.6g / %-+13.6g" % ( c.min() , c.max () ) + + if expo1 : expo1 = '10^%+d' % expo1 + else : expo1 = '' + if expo2 : expo2 = '10^%+d' % expo2 + else : expo2 = '' + + if expo3 : expo3 = '10^%+d' % expo3 + else : expo3 = '' + + if expo1 : has_e1.append ( '' ) + if expo2 : has_e2.append ( '' ) + if expo3 : has_e3.append ( '' ) + + return p , n , mean, expo1, rms, expo2, minmax, expo3 + for p in sorted ( stats ) : if p.startswith('pull:') : continue c = stats [ p ] @@ -127,16 +184,45 @@ def make_row ( c ) : c = stats [ p ] table.append ( make_row ( c ) ) + ## skip empty column + if not has_e3 : + rows = [] + for row in table : + r = list ( row ) + del r [ 7 ] + rows.append ( r ) + table = rows + + ## skip empty column + if not has_e2 : + rows = [] + for row in table : + r = list ( row ) + del r [ 5 ] + rows.append ( r ) + table = rows + + ## skip empty column + if not has_e1 : + rows = [] + for row in table : + r = list ( row ) + del r [ 3 ] + rows.append ( r ) + table = rows + + if not ntoys : ntoys = 0 for s in stats : ntoys = max ( ntoys , stats[s].nEntries() ) - + import ostap.logger.table as Table table = Table.table ( table , title = "Results of %s toys" % ntoys , - alignment = 'lcccc' , + alignment = 'lccccccc' , prefix = "# " ) + logger.info ( 'Results of %s toys:\n%s' % ( ntoys , table ) ) @@ -148,11 +234,12 @@ def make_row ( c ) : # jackknife , theta_jack = jackknife_esiimator ( statistics , value ) # @endcode def jackknife_statistics ( statistics , theta = None ) : - """Jackknife estimator from jackknife statistic + """ Jackknife estimator from jackknife statistic >>> statistics = .... >>> jacknife = jackknife_estimator ( statistics ) >>> jacknife , theta_corr , bias = jackknife_esiimator ( statistics , value ) """ + assert isinstance ( theta , VE ) or theta is None ,\ "jackknife_statistics: invalid type of 'value': %s" % type ( value ) @@ -167,24 +254,35 @@ def jackknife_statistics ( statistics , theta = None ) : bias = ( N - 1 ) * ( theta_dot.value() - theta.value() ) - ## theta corrected for the bias and with Jackknife variance + ## theta corrected for the bias and with Jackknife variance + theta_jack = VE ( theta.value() - bias , jvar ) ## corrected value return jackknife , theta_jack # ============================================================================= ## print Jackknife statistics -def print_jackknife ( fitresult , +def print_jackknife ( fitresult , ## restl of the fit of th e total datasample stats , morevars = {} , logger = logger , title = '' ) : - """print Jackknife statistics + """ Print Jackknife statistics """ + - header = ( 'Parameter' , 'theta' , 'theta_(.)' , 'theta_jack' , 'bias/sigma [%]' , 'error [%]' ) - table = [] + header = ( 'Parameter' , + 'theta' , 'x[...]' , + 'theta_(.)' , 'x[...]' , + 'bias' , 'x[...]' , + 'bias [%]' , 's/s' ) + + table = [ header ] + has_e1 = False + has_e2 = False + has_e3 = False + N = 0 for name in sorted ( stats ) : @@ -206,20 +304,45 @@ def print_jackknife ( fitresult , N = max ( N , statistics.nEntries() ) + ## jackknife estimates jackknife , theta_jack = jackknife_statistics ( statistics , theta ) - bias = theta_jack.value () - theta .value () + bias = ( N - 1 ) * ( jackknife - theta ).value() + + ## bias = theta_jack.value () - theta .value () scale = theta .error () / theta_jack.error () - row = ( name , - "%+13.6g +/- %-13.6g" % ( theta . value () , theta .error () ) , - "%+13.6g +/- %-13.6g" % ( jackknife . value () , jackknife .error () ) , - "%+13.6g +/- %-13.6g" % ( theta_jack . value () , theta_jack .error () ) , - '%+6.2f' % ( bias / theta.error() * 100 ) , - '%+6.2f' % ( scale * 100 - 100 ) ) + th1 , expo1 = pretty_ve ( theta , precision = 4 , width = 6 , parentheses = False ) + th2 , expo2 = pretty_ve ( jackknife , precision = 4 , width = 6 , parentheses = False ) + ## th3 , expo3 = pretty_ve ( theta_jack , precision = 4 , width = 6 , parentheses = False ) + th3 , expo3 = pretty_float ( bias , precision = 4 , width = 6 ) + + if expo1 : expo1 = '10^%+d' % expo1 + else : expo1 = '' + + if expo2 : expo2 = '10^%+d' % expo2 + else : expo2 = '' + + if expo3 : expo3 = '10^%+d' % expo3 + else : expo3 = '' + + bias = bias / theta.error() * 100 + + fbias = '%+.1f' % bias + errs = '%+.2f' % scale + + if expo1 : has_e1 = True + if expo2 : has_e2 = True + if expo3 : has_e3 = True + + if 50 < abs ( bias ) : fbias = attention ( fbias ) + if 0.5 < abs ( scale - 1 ) : errs = attention ( errs ) + + row = name , th1, expo1 , th2 , expo2 , th3 , expo3 , fbias , errs + table.append ( row ) @@ -230,13 +353,44 @@ def print_jackknife ( fitresult , statistics = stats [ name ] jackknife = jackknife_statistics ( statistics ) + + th2, expo2 = pretty_ve ( jackknife , precision = 4 , width = 6 , parentheses = False ) + + if expo2 : expo2 = '10^%+d' % expo2 + else : expo2 = '' + + if expo2 : has_e2 = True - row = name , '' , "%+13.6g +/- %-13.6g" % ( jackknife . value () , jackknife .error () ) , '' , '' , '' + row = name , '' , '' , th2 , expo2 , '' , '' , '' , '' table.append ( row ) - - table = [ header ] + table - + ## skip empty column + if not has_e3 : + rows = [] + for row in table : + r = list ( row ) + del r [ 6 ] + rows.append ( r ) + table = rows + + ## skip empty column + if not has_e2 : + rows = [] + for row in table : + r = list ( row ) + del r [ 4 ] + rows.append ( r ) + table = rows + + ## skip empty column + if not has_e1 : + rows = [] + for row in table : + r = list ( row ) + del r [ 2 ] + rows.append ( r ) + table = rows + title = title if title else "Jackknife results (N=%d)" % N import ostap.logger.table as Table @@ -254,14 +408,17 @@ def print_bootstrap ( fitresult , morevars = {} , logger = logger , title = '' ) : - """print Bootstrap statistics + """ Print Bootstrap statistics """ - header = ( 'Parameter' , 'theta' , 'theta_boot' , 'bias/sigma [%]' , 'error [%]' ) + header = ( 'Parameter' , 'theta' , 'x[...]' , 'theta_boot' , 'x[...]', 'bias[%]' , 's/s' ) table = [] n = 0 + has_e1 = False + has_e2 = False + for name in sorted ( stats ) : if name in fitresult : @@ -284,15 +441,30 @@ def print_bootstrap ( fitresult , theta_boot = VE ( statistics.mean().value() , statistics.mu2() ) - bias = theta_boot.value () - theta .value () - scale = theta .error () / theta_boot.error () + scale = theta .error () / theta_boot.error () + + th1 , expo1 = pretty_ve ( theta , precision = 4 , width = 6 , parentheses = False ) + th2 , expo2 = pretty_ve ( theta_boot , precision = 4 , width = 6 , parentheses = False ) + + if expo1 : expo1 = '10^%+d' % expo1 + else : expo1 = '' + + if expo2 : expo2 = '10^%+d' % expo2 + else : expo2 = '' + + bias = bias / theta.error() * 100 + + fbias = '%+.1f' % bias + errs = '%+.2f' % scale - row = ( name , - "%+13.6g +/- %-13.6g" % ( theta . value () , theta .error () ) , - "%+13.6g +/- %-13.6g" % ( theta_boot . value () , theta_boot .error () ) , - '%+6.2f' % ( bias / theta.error() * 100 ) , - '%+6.2f' % ( scale * 100 - 100 ) ) + if 50 < abs ( bias ) : fbias = attention ( fbias ) + if 0.5 < abs ( scale - 1 ) : errs = attention ( errs ) + + row = name , th1, expo1 , th2 , expo2 , fbias , errs + + if expo1 : has_e1 = True + if expo2 : has_e2 = True table.append ( row ) @@ -304,19 +476,42 @@ def print_bootstrap ( fitresult , statistics = stats [ name ] theta_boot = VE ( statistics.mean().value() , statistics.mu2() ) - row = name , '', "%+13.6g +/- %-13.6g" % ( theta_boot . value () , theta_boot .error () ) , '' , '' + th2 , expo2 = pretty_ve ( theta_boot , precision = 4 , width = 6 , parentheses = False ) + + if expo2 : expo2 = '10^%+d' % expo2 + else : expo2 = '' + + if expo2 : has_e2 = True + row = name , '','' , th2 , expo2 , '' , '' table.append ( row ) - + table = [ header ] + table - table = [ header ] + table + ## skip empty column + if not has_e2 : + rows = [] + for row in table : + r = list ( row ) + del r [ 4 ] + rows.append ( r ) + table = rows + + ## skip empty column + if not has_e1 : + rows = [] + for row in table : + r = list ( row ) + del r [ 2 ] + rows.append ( r ) + table = rows + title = title if title else "Bootstrapping with #%d samples" % n import ostap.logger.table as Table table = Table.table ( table , title = title , - alignment = 'lcccc' , + alignment = 'lccccccc' , prefix = "# " ) logger.info ( '%s:\n%s' % ( title , table ) ) @@ -325,7 +520,7 @@ def print_bootstrap ( fitresult , ## Default function to generate the data # - simple call for PDF.generate def generate_data ( pdf , varset , **config ) : - """Default function to generate the data + """ Default function to generate the data - simple call for `PDF.generate` """ return pdf.generate ( varset = varset , **config ) @@ -353,7 +548,7 @@ def make_fit ( pdf , dataset , **config ) : # @param dataset pdf # def accept_fit ( result , pdf = None , dataset = None ) : - """Accept the fit result? + """ Accept the fit result? - valid fit result - fit status is 0 (SUCCESS) - covariance matrix quality is either 0 or -1 @@ -431,7 +626,7 @@ def make_toys ( pdf , progress = True , logger = logger , frequency = 500 ) : ## - """Make `nToys` pseudoexperiments + """ Make `nToys` pseudoexperiments - Schematically: >>> for toy in range ( nToys ) : @@ -575,11 +770,13 @@ def make_toys ( pdf , ## 3. fit it! r = fit_fun ( pdf , dataset , **fitcnf ) - ## fit status - fits [ r.status () ] += 1 + ## fit status + st = r.status() + if 0 != st : fits [ st ] += 1 ## covariance matrix quality - covs [ r.covQual () ] += 1 + cq = r.covQual() + if not cq in ( -1 , 3 ) : covs [ cq ] += 1 ## ok ? ok = accept_fun ( r , pdf , dataset ) @@ -839,10 +1036,12 @@ def make_toys2 ( gen_pdf , ## pdf to generate toys r = fit_fun ( fit_pdf , dataset , **fitcnf ) ## fit status - fits [ r.status () ] += 1 + st = r.status() + if 0 != st : fits [ st ] += 1 ## covariance matrix quality - covs [ r.covQual () ] += 1 + cq = r.covQual() + if not cq in ( -1 , 3 ) : covs [ cq ] += 1 ## ok ? ok = accept_fun ( r , fit_pdf , dataset ) @@ -925,14 +1124,14 @@ def make_jackknife ( pdf , fit_config = {} , ## parameters for pdf.fitTo fit_pars = {} , ## fit-parameters to reset/use more_vars = {} , ## additional results to be calculated - fit_fun = None , ## fit function ( pdf , dataset , **fit_config ) - accept_fun = None , ## accept function ( fit-result, pdf, dataset ) + fit_fun = None , ## - fit function ( pdf , dataset , **fit_config ) + accept_fun = None , ## - accept function ( fit-result, pdf, dataset ) event_range = () , ## event range for jackknife silent = True , progress = True , logger = logger , frequency = 500 ) : - """Run Jackknife analysis, useful for evaluaton of fit biased and uncertainty estimates + """ Run Jackknife analysis, useful for evaluaton of fit biased and uncertainty estimates For each i remove event with index i from the dataset, and refit it. >>> dataset = ... >>> model = ... @@ -970,7 +1169,11 @@ def make_jackknife ( pdf , assert 0 <= begin and begin < end and begin < N , 'make_jackknife: invalid event range (%s,%s)/%d' % ( begin , end , N ) ## adjust the end end = min ( end , N ) - + + ## check if pdf and `extended` flag are in agreemrnt, and print warning message otherwise + ok1 = check_extended ( pdf , False ) + if not ok1 and not silent : logger.warning ( "Jackknife: estimates for `yield` parameters could be wrong!") + ## 1. fitting function? if fit_fun is None : if not silent : logger.info ( "make_jackknife: use default 'make_fit' function!") @@ -1024,10 +1227,12 @@ def make_jackknife ( pdf , r = fit_fun ( pdf , ds , **fitcnf ) ## 4. fit status - fits [ r.status () ] += 1 + st = r.status() + if 0 != st : fits [ st ] += 1 ## 5. covariance matrix quality - covs [ r.covQual () ] += 1 + cq = r.covQual() + if not cq in ( -1 , 3 ) : covs [ cq ] += 1 ## ok ? ok = accept_fun ( r , pdf , ds ) @@ -1050,8 +1255,10 @@ def make_jackknife ( pdf , results [ '#' ] .append ( len ( ds ) ) results [ '#sumw' ] .append ( ds.sumVar ( '1' ) ) + ## reset/remove/delete dataset ds.clear() - + del ds + if progress or not silent : if 0 < frequency and 1 <= i and 0 == i % frequency : stats = make_stats ( results , fits , covs , accept ) @@ -1066,19 +1273,27 @@ def make_jackknife ( pdf , r_tot = fit_fun ( pdf , data , **fitcnf ) r_tot = fit_fun ( pdf , data , **fitcnf ) - ## 10. the final table + ## 10. the final table + title = '' + if not ok1 : + logger.warning ( "Jackknife: estimates for `yield` parameters are likely wrong!") + NN = 0 + for k in stats : NN = max ( NN , stats[k].nEntries() ) + title = "Jackknife results (N=%d).[Estimates for `yield` are likely wrong!]" % NN print_jackknife ( r_tot , stats , morevars = dict ( ( k , more_vars [ k ]( r_tot , pdf ) ) for k in more_vars ) , - logger = logger ) - + logger = logger , + title = title ) + + if not ok1 and not silent : logger.warning ( "Jackknife: estimates for `yield` parameters are likely wrong!") return results , stats # ============================================================================= ## Run Bootstrap analysis, useful for evaluaton of fit biases and uncertainty estimates # -# In total size datasets are sampled (with replacement) from the original dataste +# In total size datasets are sampled (with replacement) from the original dataset # data and each sampled dataset is fit # @code # dataset = ... @@ -1102,6 +1317,7 @@ def make_jackknife ( pdf , # @param fit_config configuration of pdf.FitTo( data , ... ) # @param fit_pars redefine these parameters before each fit # @param more_vars calculate more variables from the fit-results +# @param extended use extended bootstrap? # @param fit_fun specific fitting action (if needed) # @param accept_fun specific accept action (if needed) # @param silent silent processing @@ -1115,14 +1331,15 @@ def make_bootstrap ( pdf , fit_config = {} , ## parameters for pdf.fitTo fit_pars = {} , ## fit-parameters to reset/use more_vars = {} , ## additional results to be calculated - fit_fun = None , ## fit function ( pdf , dataset , **fit_config ) + extended = True , ## use extended/non-extended bootstrtap + fit_fun = None , ## fit function ( pdf , dataset , **fit_config ) accept_fun = None , ## accept function ( fit-result, pdf, dataset ) silent = True , ## silent processing? progress = True , ## shpow progress bar? - logger = logger , ## use this logger + logger = logger , ## use this logger frequency = 500 ) : - """Run Bootstrap analysis, useful for evaluaton of fit biased and uncertainty estimates + """ Run Bootstrap analysis, useful for evaluaton of fit biased and uncertainty estimates In total `size` datasets are sampled (with replacement) from the original dataste `data` and each sampled dataset is fit >>> dataset = ... @@ -1137,6 +1354,7 @@ def make_bootstrap ( pdf , - `fit_config` : configuration of `pdf.FitTo( data , ... )` - `fit_pars` : redefine these parameters before each fit - `more_vars` : calculate more variables from the fit-results + - `extended` : use extended bootstrap? - `fit_fun` : specific fitting acion (if needed) - `accept_fun` : specific accept action (if needed) - `silent` : silent processing? @@ -1152,6 +1370,10 @@ def make_bootstrap ( pdf , assert isinstance ( size , integer_types ) and 1 <= size ,\ 'Invalid "size" argument %s/%s' % ( size , type ( size ) ) + ## check if pdf and `extended` flag are in agreemrnt, and print warning message otherwise + ok1 = check_extended ( pdf , extended ) + if not ok1 and not silent : logger.warning ( "Bootstrap: estimates for `yield` parameters could be wrong!") + ## 1. fitting function? if fit_fun is None : if not silent : logger.info ( "make_bootstrap: use default 'make_fit' function!") @@ -1193,10 +1415,14 @@ def make_bootstrap ( pdf , pdf.load_params ( params = fix_fit_pars , silent = True ) ## silent = silent ) r_tot = fit_fun ( pdf , data , **fitcnf ) + nn = 0 + from ostap.utils.progress_bar import progress_bar - ## run jackknife bootstrapping - for i , ds in progress_bar ( enumerate ( data.bootstrap ( size ) ) , max_value = size , silent = not progress ) : - + ## run bootstrapping + for i , ds in progress_bar ( enumerate ( data.bootstrap ( size , extended = extended ) ) , + max_value = size , + silent = not progress ) : + ## 2. reset parameters of fit_pdf pdf.load_params ( params = fix_fit_init , silent = True ) ## silent = silent ) pdf.load_params ( params = fix_fit_pars , silent = True ) ## silent = silent ) @@ -1205,10 +1431,12 @@ def make_bootstrap ( pdf , r = fit_fun ( pdf , ds , **fitcnf ) ## 4. fit status - fits [ r.status () ] += 1 + st = r.status() + if 0 != st : fits [ st ] += 1 ## 5. covariance matrix quality - covs [ r.covQual () ] += 1 + cq = r.covQual() + if not cq in ( -1 , 3 ) : covs [ cq ] += 1 ## ok ? ok = accept_fun ( r , pdf , ds ) @@ -1230,18 +1458,24 @@ def make_bootstrap ( pdf , results [ '#' ] .append ( len ( ds ) ) results [ '#sumw' ] .append ( ds.sumVar ( '1' ) ) + nn += 1 + ## reset/remove/delete dataset ds.clear() + del ds if progress or not silent : if 0 < frequency and 1 <= i and 0 == i % frequency : stats = make_stats ( results , fits , covs , accept ) - ## print_stats ( stats , i + 1 , logger = logger ) + title = "Bootstrapping with #%d samples" % nn + if extended : title = '(Extended) %s' % title + if not ok1 : title += '[Estimates for `yield` are likely wrong!]' print_bootstrap ( r_tot , stats , morevars = dict ( ( k , more_vars [ k ] ( r_tot , pdf ) ) for k in more_vars ), - logger = logger ) - + logger = logger , + title = title ) + if not ok1 and not silent : logger.warning ( "Bootstrap: estimates for `yield` parameters are likely wrong!") ## 8. make a final statistics stats = make_stats ( results , fits , covs , accept ) @@ -1253,11 +1487,19 @@ def make_bootstrap ( pdf , r_tot = fit_fun ( pdf , data , **fitcnf ) ## 10. the final table + title = "Bootstrapping with #%d samples" % nn + if extended : title = '(Extended) %s' % title + if not ok1 : + title += '[Estimates for `yield` are likely wrong!]' + logger.warning ( "Bootstrap: estimates for `yield` parameters are likely wrong!") print_bootstrap ( r_tot , stats , morevars = dict ( ( k , more_vars [ k ]( r_tot , pdf ) ) for k in more_vars ), - logger = logger ) + logger = logger , + title = title ) + if not ok1 and not silent : logger.warning ( "Bootstrap: estimates for `yield` parameters are likely wrong!") + return results , stats diff --git a/ostap/parallel/parallel_toys.py b/ostap/parallel/parallel_toys.py index f825bbbc..fd506d94 100644 --- a/ostap/parallel/parallel_toys.py +++ b/ostap/parallel/parallel_toys.py @@ -302,7 +302,8 @@ def __init__ ( self , silent = True , progress = True , frequency = 100 ) : - + + TheTask_.__init__ ( self , data = data , fit_config = fit_config , @@ -354,9 +355,38 @@ def process ( self , jobid , event_range ) : # @author Vanya BELYAEV Ivan.Belyaev@itep.ru # @date 2023-06-22 class BootstrapTask(JackknifeTask) : - """The simple task object for parallel Jackknife + """ The simple task object for parallel Jackknife - see ostap.fitting.toys.make_jackknife """ + + def __init__ ( self , + pdf , + data , + fit_config = {} , ## parameters for pdf.fitTo + fit_pars = {} , ## fit-parameters to reset/use + more_vars = {} , ## additional results to be calculated + extended = True , ## use extended bootstrap + fit_fun = None , ## fit function ( pdf , dataset , **fit_config ) + accept_fun = None , ## accept function ( fit-result, pdf, dataset ) + silent = True , + progress = True , + frequency = 100 ) : + + JackknifeTask.__init__ ( self , + pdf = pdf , + data = data , + fit_config = fit_config , ## parameters for pdf.fitTo + fit_pars = fit_pars , ## fit-parameters to reset/use + more_vars = more_vars , ## additional results to be calculated + fit_fun = fit_fun , ## fit function ( pdf , dataset , **fit_config ) + accept_fun = accept_fun , ## accept function ( fit-result, pdf, dataset ) + silent = silent , + progress = progress , + frequency = frequency ) + + self.extended = True if extended else False + + ## the actual processing def process ( self , jobid , size ) : @@ -370,16 +400,17 @@ def process ( self , jobid , size ) : import ostap.fitting.variables from ostap.core.ostap_types import integer_types - assert isinstance ( size , integer_types ) and 0 < size ,\ + assert isinstance ( size , integer_types ) and 1 <= size ,\ 'Jobid %s: Invalid "size" argument %s/%s' % ( jobid , size , type ( size ) ) import ostap.fitting.toys as Toys return Toys.make_bootstrap ( pdf = self.pdf , data = self.data , - size = size , ## ATTENTION + size = size , ## ATTENTION !!! fit_config = self.fit_config , fit_pars = self.fit_pars , more_vars = self.more_vars , + extended = self.extended , fit_fun = self.fit_fun , accept_fun = self.accept_fun , silent = self.silent , @@ -867,6 +898,10 @@ def parallel_jackknife ( pdf , import ostap.fitting.toys as Toys if nSplit < 2 : return Toys.make_jackknife ( progress = progress , **config ) + + ## check if pdf and `extended` flag are in agreemrnt, and print warning message otherwise + ok1 = Toys.check_extended ( pdf , False ) + if not ok1 : logger.warning ( "Jackknife: estimates for `yield` parameters could be wrong!") ## create the task task = JackknifeTask ( progress = progress and not silent , **config ) @@ -885,8 +920,6 @@ def parallel_jackknife ( pdf , results , stats = task.results () if progress or not silent : - ## Toys.print_stats ( stats ) - fitcnf = {} fitcnf.update ( fit_config ) if not 'silent' in fitcnf : fitcnf [ 'silent' ] = silent @@ -896,14 +929,21 @@ def parallel_jackknife ( pdf , ## 9. fit total dataset (twice) r_tot = fit_fun ( pdf , data , **fitcnf ) r_tot = fit_fun ( pdf , data , **fitcnf ) - + + title = '' + if not ok1 : + logger.warning ( "Jackknife: estimates for `yield` parameters are likely wrong!") + NN = 0 + for k in stats : NN = max ( NN , stats[k].nEntries() ) + title = "Jackknife results (N=%d).[Estimates for `yield` are likely wrong!]" % NN ## the final table Toys.print_jackknife ( r_tot , stats , - morevars = dict ( ( k , more_vars [ k ]( r_tot , pdf ) ) for k in more_vars ) ) - + morevars = dict ( ( k , more_vars [ k ]( r_tot , pdf ) ) for k in more_vars ) , + title = title ) + if not ok1 : logger.warning ( "Jackknife: estimates for `yield` parameters are likely wrong!") return results, stats @@ -948,6 +988,7 @@ def parallel_bootstrap ( pdf , fit_config = {} , ## parameters for pdf.fitTo fit_pars = {} , ## fit-parameters to reset/use more_vars = {} , ## additional results to be calculated + extended = True , ## use extended bootstrap fit_fun = None , ## fit function ( pdf , dataset , **fit_config ) accept_fun = None , ## accept function ( fit-result, pdf, dataset ) silent = True , ## silent processing? @@ -973,7 +1014,10 @@ def parallel_bootstrap ( pdf , import ostap.fitting.toys as Toys if nSplit < 2 : return Toys.make_boostrap ( size = size , progress = progress , **config ) - + + ## check if pdf and `extended` flag are in agreemrnt, and print warning message otherwise + ok1 = Toys.check_extended ( pdf , extended ) + if not ok1 : logger.warning ( "Bootstrap: estimates for `yield` parameters could be wrong!") ## create teh task task = BootstrapTask ( progress = progress and not silent , **config ) @@ -1002,13 +1046,19 @@ def parallel_bootstrap ( pdf , r_tot = fit_fun ( pdf , data , **fitcnf ) r_tot = fit_fun ( pdf , data , **fitcnf ) - ## the final table + ## the final table + title = '' + if not ok1 : + logger.warning ( "Bootstrap: estimates for `yield` parameters are likely wrong!") + NN = 0 + for k in stats : NN = max ( NN , stats[k].nEntries() ) + title = "Bootstrap results (N=%d).[Estimates for `yield` are likely wrong!]" % NN Toys.print_bootstrap ( r_tot , stats , morevars = dict ( ( k , more_vars [ k ]( r_tot , pdf ) ) for k in more_vars ) ) - + if not ok1 : logger.warning ( "Bootstrap: estimates for `yield` parameters are likely wrong!") return results, stats diff --git a/ostap/parallel/tests/test_parallel_toys.py b/ostap/parallel/tests/test_parallel_toys.py index e5ece3ea..a89f9e6c 100755 --- a/ostap/parallel/tests/test_parallel_toys.py +++ b/ostap/parallel/tests/test_parallel_toys.py @@ -18,12 +18,16 @@ from ostap.core.meta_info import root_info import ostap.fitting.roofit import ostap.fitting.models as Models -from ostap.utils.timing import timing +from ostap.utils.timing import timing +from ostap.utils.memory import memory from ostap.core.core import cpp, VE, dsID, hID, rooSilent import ostap.fitting.models as Models from ostap.plotting.canvas import use_canvas from ostap.utils.utils import wait from ostap.fitting.toys import pull_var +from ostap.math.base import num_range +from ostap.utils.progress_bar import progress_bar +import ostap.logger.table as T import ROOT, sys, os, random, time # ============================================================================= # logging @@ -47,16 +51,27 @@ nominal_mean = 0.4 nominal_sigma = 0.1 +nominal_S = 1000 +nominal_B = 100 +nominal_F = float ( nominal_S ) / ( nominal_S + nominal_B ) -mass = ROOT.RooRealVar ( 'mass' , '', 0 , 1 ) -gen_gauss = Models.Gauss_pdf ( 'GG' , xvar = mass ) -fit_gauss = Models.Gauss_pdf ( 'FG' , xvar = mass ) + +mass = ROOT.RooRealVar ( 'mass' , '', 0 , 1 ) +gen_gauss = Models.Gauss_pdf ( 'GG' , xvar = mass ) +fit_gauss = Models.Gauss_pdf ( 'FG' , xvar = mass ) gen_gauss.mean = nominal_mean gen_gauss.sigma = nominal_sigma -model = Models.Fit1D ( signal = gen_gauss , - background = 'flat' ) +model = Models.Fit1D ( signal = gen_gauss , background = 'flat' ) +model_NE = Models.Fit1D ( signal = gen_gauss , background = 'flat' , extended = False ) + +model.S = nominal_S +model.B = nominal_B +model_NE.F = nominal_F + +toy_results = {} +fit_config = { 'silent' : True , 'refit' : 5 } # ============================================================================== ## Perform toy-study for possible fit bias and correct uncertainty evaluation @@ -67,7 +82,7 @@ # - fill distributions for fit results # - fill distribution of pulls def test_parallel_toys ( ) : - """Perform toys-study for possible fit bias and correct uncertainty evaluation + """ Perform toys-study for possible fit bias and correct uncertainty evaluation - generate `nToys` pseudoexperiments with some PDF `pdf` - fit teach experiment with the same PDF - store fit results @@ -84,11 +99,11 @@ def test_parallel_toys ( ) : with timing ( 'Toys analysis' , logger = logger ) : results , stats = Toys.parallel_toys ( pdf = gen_gauss , - nToys = 1000 , + nToys = 2000 , nSplit = 20 , data = [ mass ] , gen_config = { 'nEvents' : 200 , 'sample' : True } , - fit_config = { 'silent' : True } , + fit_config = fit_config , init_pars = { 'mean_GG' : nominal_mean , 'sigma_GG' : nominal_sigma } , more_vars = more_vars , silent = True , @@ -105,6 +120,8 @@ def test_parallel_toys ( ) : for h in ( h_mean , h_sigma ) : with use_canvas ( 'test_paralllel_toys %s' % h.title , wait = 1 ) : h.draw() + toy_results['test_parallel_toys'] = results , stats + # ============================================================================= ## Perform toy-study for possible fit bias and correct uncertainty evaluation # - generate nToys pseudoexperiments with some PDF gen_pdf @@ -112,7 +129,7 @@ def test_parallel_toys ( ) : # - store fit results # - fill distributions for fit results def test_parallel_toys2 ( ) : - """Perform toys-study for possible fit bias and correct uncertainty evaluation + """ Perform toys-study for possible fit bias and correct uncertainty evaluation - generate `nToys` pseudoexperiments with some PDF `gen_pdf` - fit teach experiment with the PDF `fit_pdf` - store fit results @@ -128,11 +145,11 @@ def test_parallel_toys2 ( ) : results , stats = Toys.parallel_toys2 ( gen_pdf = gen_gauss , fit_pdf = fit_gauss , - nToys = 1000 , + nToys = 2000 , nSplit = 20 , data = [ mass ] , gen_config = { 'nEvents' : 200 , 'sample' : True } , - fit_config = { 'silent' : True } , + fit_config = fit_config , gen_pars = { 'mean_GG' : nominal_mean , 'sigma_GG' : nominal_sigma } , fit_pars = { 'mean_FG' : nominal_mean , 'sigma_FG' : nominal_sigma } , more_vars = more_vars , @@ -150,58 +167,221 @@ def test_parallel_toys2 ( ) : for h in ( h_mean , h_sigma ) : with use_canvas ( 'test_parallel_toys2 %s' % h.title , wait = 1 ) : h.draw() + toy_results['test_parallel_toys2'] = results , stats -# ============================================================================== -## Perform toy-study for Jackknife -def test_parallel_jackknife ( ) : - """Perform toys-study for Jackknife - """ - logger = getLogger ( 'test_parallel_jackknife' ) - if root_info < (6,18) : - logger.warning ( "(parallel) Jackknife is temporarily disabled for this version of ROOT" ) - return +# ============================================================================== +## Perform toy-study for Jackknife (non-extended model) +def test_parallel_jackknife_NE1 ( ) : + """ Perform toys-study for Jackknife (non-extended model) + """ - with timing ( 'Jackknife analysis' , logger = logger ) : - - model.S = 1000 - model.B = 100 - - dataset = model.generate ( 1100 ) - - Toys.parallel_jackknife ( pdf = model , - nSplit = 6 , - data = dataset , - fit_config = { 'silent' : True } , - progress = True , - silent = True ) - + logger = getLogger ( 'test_parallel_jackknife_NE1' ) + logger.info ( 'Jackknife analysis for non-extended model' ) + + pdf = model_NE + + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.F = nominal_F + + with timing ( 'Jackknife analysis (non-extended)' , logger = logger ) : + data = pdf.generate ( nominal_S ) + with use_canvas ( "Jackknife analysis (non-extended)" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.F = nominal_F + + results, stats = Toys.parallel_jackknife ( pdf = pdf , + nSplit = 10 , + data = data , + fit_config = fit_config , + progress = True , + silent = True ) + + toy_results [ 'test_parallel_jackknife_NE1' ] = results , stats + +# ============================================================================== +## Perform toy-study for Jackknife (non-extended model) +def test_parallel_jackknife_NE2 ( ) : + """ Perform toys-study for Jackknife (non-extended model) + """ + + logger = getLogger ( 'test_parallel_jackknife_NE2' ) + logger.info ( 'Jackknife analysis for non-extended model' ) + + pdf = fit_gauss + + pdf.mean = nominal_mean + pdf.sigma = nominal_sigma + + with timing ( 'Jackknife analysis (non-extended)' , logger = logger ) : + data = pdf.generate ( nominal_S ) + with use_canvas ( "Jackknife analysis (non-extended)" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + pdf.mean = nominal_mean + pdf.sigma = nominal_sigma + + results, stats = Toys.parallel_jackknife ( pdf = pdf , + nSplit = 10 , + data = data , + fit_config = fit_config , + progress = True , + silent = True ) + + toy_results [ 'test_parallel_jackknife_NE2' ] = results , stats + # ============================================================================== -## Perform toy-study for Bootstrap -def test_parallel_bootstrap ( ) : - """Perform toys-study for Bootstrap +## Perform toy-study for Jackknife (extended model) +def test_parallel_jackknife_EXT ( ) : + """ Perform toys-study for Jackknife (extended model) """ - logger = getLogger ( 'test_parallel_bootstrap' ) + logger = getLogger ( 'test_parallel_jackknife_EXT' ) + logger.info ( 'Jackknife analysis for extended model' ) + + pdf = model - with timing ( 'Bootstrap analysis' , logger = logger ) : - - model.S = 1000 - model.B = 100 + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.S = nominal_S + pdf.B = nominal_B + + with timing ( 'Jackknife analysis (extended)' , logger = logger ) : - dataset = model.generate ( 1100 ) + data = pdf.generate ( nominal_S ) + with use_canvas ( "Jackknife analysis (extended)" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.S = nominal_S + pdf.B = nominal_B + + results, stats = Toys.parallel_jackknife ( pdf = pdf , + nSplit = 10 , + data = data , + fit_config = fit_config , + progress = True , + silent = True ) + + toy_results [ 'test_parallel_jackknife_EXT' ] = results , stats - Toys.parallel_bootstrap ( pdf = model , - size = 1000 , - nSplit = 6 , - data = dataset , - fit_config = { 'silent' : True } , - progress = True , - silent = True ) +# ============================================================================== +## Perform toy-study for Bootstrap (non-extended) +def test_parallel_bootstrap_NE1 ( ) : + """ Perform toys-study for Bootstrap (non-extended) + """ + + logger = getLogger ( 'test_parallel_bootstrap_NE1' ) + logger.info ( 'Bootstrap analysis for non-extended model' ) + + pdf = model_NE + + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.F = nominal_F + + with timing ( 'Bootstrap analysis (non-extended)' , logger = logger ) : + + data = pdf.generate ( nominal_S + nominal_B ) + with use_canvas ( "Booststrap analysis" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.F = nominal_F + + results, stats = Toys.parallel_bootstrap ( pdf = pdf , + nSplit = 10 , + size = 400 , + data = data , + fit_config = fit_config , + extended = False , + progress = True , + silent = True ) + + toy_results [ 'test_parallel_bootstrap_NE1' ] = results , stats + +# ============================================================================== +## Perform toy-study for Bootstrap (non-extended) +def test_parallel_bootstrap_NE2 ( ) : + """ Perform toys-study for Bootstrap (non-extended) + """ + + logger = getLogger ( 'test_parallel_bootstrap_NE2' ) + logger.info ( 'Bootstrap analysis for non-extended model' ) + + pdf = fit_gauss + + pdf.mean = nominal_mean + pdf.sigma = nominal_sigma + + with timing ( 'Bootstrap analysis (non-extended)' , logger = logger ) : + + data = pdf.generate ( nominal_S + nominal_B ) + with use_canvas ( "Booststrap analysis" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + pdf.mean = nominal_mean + pdf.sigma = nominal_sigma + + results, stats = Toys.parallel_bootstrap ( pdf = pdf , + nSplit = 10 , + size = 400 , + data = data , + fit_config = fit_config , + extended = False , + progress = True , + silent = True ) + + toy_results [ 'test_parallel_bootstrap_NE2' ] = results , stats + +# ============================================================================== +## Perform toy-study for Bootstrap (extended) +def test_parallel_bootstrap_EXT ( ) : + """ Perform toys-study for Bootstrap (extended) + """ + + logger = getLogger ( 'test_parallel_bootstrap_EXT' ) + logger.info ( 'Bootstrap analysis for extended model' ) + + pdf = model + + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.S = nominal_S + pdf.B = nominal_B + + with timing ( 'Bootstrap analysis (extended)' , logger = logger ) : + + data = pdf.generate ( nominal_S + nominal_B ) + with use_canvas ( "(Extended) Booststrap analysis" , wait = 3 ) : + r , f = pdf.fitTo ( data , draw = True , nbins = 100 , **fit_config ) + + pdf.signal.mean = nominal_mean + pdf.signal.sigma = nominal_sigma + pdf.S = nominal_S + pdf.B = nominal_B + + results, stats = Toys.parallel_bootstrap ( pdf = pdf , + nSplit = 10 , + size = 400 , + data = data , + fit_config = fit_config , + extended = True , + progress = True , + silent = True ) + + toy_results [ 'test_parallel_bootstrap_EXT' ] = results , stats + # ============================================================================= ## Perform toy-study for significance of the signal # - generate nToys pseudoexperiments using background-only hypothesis @@ -209,7 +389,7 @@ def test_parallel_bootstrap ( ) : # - store fit results # - fill distributions for fit results def test_parallel_significance ( ) : - """Perform toy-study for significance of the signal + """ Perform toy-study for significance of the signal - generate `nToys` pseudoexperiments using background-only hypothesis - fit each experiment with signal+background hypothesis - store fit results @@ -236,34 +416,117 @@ def test_parallel_significance ( ) : results , stats = Toys.parallel_toys2 ( gen_pdf = bkg_only , fit_pdf = model , - nToys = 1000 , + nToys = 10000 , nSplit = 20 , data = [ mass ] , gen_config = { 'nEvents' : 100 , 'sample' : True } , - fit_config = { 'silent' : True } , + fit_config = fit_config , gen_pars = { 'tau_BKG' : 0. } , ## initial values for generation fit_pars = { 'B' : 100 , 'S' : 10 , 'phi0_Bkg_S': 0.0 } , ## initial fit values for parameters silent = True , progress = True ) - - h_S = ROOT.TH1F ( hID() , '#S' , 60 , 0 , 60 ) - for r in results ['S' ] : h_S .Fill ( r ) + + + # ========================================================================= + ## yields + h_Y = ROOT.TH1F ( hID() , '#S' , 140 , 0 , 70 ) + for r in results [ 'S' ] : h_Y .Fill ( r ) + + ## get p-value and significance histograms: + h_P , h_S = h_Y.significance ( cut_low = True ) + + h_S.red () + h_P.blue () + + minv = 1.e+100 + for i,_,y in h_P.items() : + yv = y.value() + if 0 < yv and yv < minv : minv = yv + minv , _ = num_range ( 0.75 * minv ) + h_P.SetMinimum ( minv ) + h_S.SetMinimum ( 0 ) + + with use_canvas ( 'test_significance: yields' , wait = 1 ) : h_Y.draw ( ) + with use_canvas ( 'test_significance: p-value' , wait = 1 ) : h_P.draw ( logy = True ) + with use_canvas ( 'test_significance: significance' , wait = 3 ) : h_S.draw () - for h in ( h_S , ) : - with use_canvas ( 'test_parallel_significance %s' % h.title , wait = 1 ) : h.draw() + toy_results [ 'test_parallel_significance' ] = results , stats + +# ============================================================================= +## Save resuls of toys to DBASE +def test_db () : + """ Save resuls of toys to DBASE + """ + logger = getLogger ( 'test_db' ) + logger.info ( 'Saving all toys results into DBASE' ) + import ostap.io.zipshelve as DBASE + from ostap.utils.timing import timing + with timing( 'Save everything to DBASE', logger ), DBASE.tmpdb() as db : + for key in progress_bar ( toy_results ) : + r , s = toy_results [ key ] + key1 = '%s:results' % key + key2 = '%s:stats' % key + db [ key1 ] = r + db [ key2 ] = s + db.ls() + while toy_results : toy_results.popitem() +# ============================================================================= + + # ============================================================================= if '__main__' == __name__ : + with memory ( 'test_toys' , logger = logger ) as mtt : test_parallel_toys () + with memory ( 'test_toys2' , logger = logger ) as mt2 : test_parallel_toys2 () + with memory ( 'test_jackknife_NE1' , logger = logger ) as mj1 : test_parallel_jackknife_NE1 () + with memory ( 'test_jackknife_NE2' , logger = logger ) as mj2 : test_parallel_jackknife_NE2 () + with memory ( 'test_jackknife_EXT' , logger = logger ) as mje : test_parallel_jackknife_EXT () + with memory ( 'test_bootstrap_NE1' , logger = logger ) as mb1 : test_parallel_bootstrap_NE1 () + with memory ( 'test_bootstrap_NE2' , logger = logger ) as mb2 : test_parallel_bootstrap_NE2 () + with memory ( 'test_bootstrap_EXT' , logger = logger ) as mbe : test_parallel_bootstrap_EXT () + + with memory ( 'test_significance' , logger = logger ) as mts : test_parallel_significance () + with memory ( 'test_db' , logger = logger ) as mdb : test_db () + + rows = [ ( 'Test' , 'Memory [MB]' ) ] + + row = 'Toys' , '%.1f' % mtt.delta + rows.append ( row ) - test_parallel_toys ( ) - test_parallel_toys2 ( ) - test_parallel_significance ( ) + row = 'Toys2' , '%.1f' % mt2.delta + rows.append ( row ) + + row = 'Jackknife_NE1' , '%.1f' % mj1.delta + rows.append ( row ) + + row = 'Jackknife_NE2' , '%.1f' % mj2.delta + rows.append ( row ) + + row = 'Jackknife_EXT' , '%.1f' % mje.delta + rows.append ( row ) + + row = 'Bootstrap_NE1' , '%.1f' % mb1.delta + rows.append ( row ) + + row = 'Bootstrap_NE2' , '%.1f' % mb2.delta + rows.append ( row ) + + row = 'Bootstrap_EXT' , '%.1f' % mbe.delta + rows.append ( row ) + + row = 'Significance' , '%.1f' % mts.delta + rows.append ( row ) + + row = 'test_db' , '%.1f' % mdb.delta + rows.append ( row ) + + title = 'Memory usage for various toys' + table = T.table ( rows , title = title , prefix = '# ' , alignment = 'lc' ) + logger.info ( '%s:\n%s' % ( title , table ) ) - test_parallel_jackknife ( ) - test_parallel_bootstrap ( ) # =============================================================================