Repeatability issue with (py)RooFit in a notebook

Roofit and RooStats
, , ,
1

Hello Root experts,
I currently try to use pyROOT in a notebook, in order to fit a RooDataSet of mass with a 2 Crystal Ball (one right and one left sharing the same \mu and \sigma) + a Gaussian model. The RooDataset is made with the .from_numpy() method. When I run the cells containing the fits (eg one for the model definition and another one to perform the fit) for the first time my fit works (it’s fast). If I rerun just the cell that perform the fit, it is still working in the same way, but when I rerun the cell that define the model before (even if I don’t modify my model) and then reperform the fit => the cell run indefinitely and the fit is not done. The only solution to be able to rerun code is to restart the kernel …

As a zfit user, I give a name with an iterable to each parameter of the fit. My notebook runs in a Conda environment with ROOT V6.26.10 . Apparently when I try a simple exemple like fitting a normally distributed list with a Gaussian, I don’t have this issue.

Do you have any idea about how to solve my repeatability issue (if I want to play a bit with my model, I need to be able to rerun it easily) ?

2

Please find there the code of the three cells involved in this fit:

mass = ROOT.RooRealVar("mass", "mass", 4.6, 6.)
npsigmass=np.array(data[f"sig_{smearing}"][f"sel_plot"])
signal_mass = ROOT.RooDataSet.from_numpy(data={"mass": npsigmass}, variables=[mass])

low1=5.26
high1=5.3
low = 4.6
high = 6

rand+=1

#signal fit
#crystal ball parameters
mu=ROOT.RooRealVar(f"mu_{rand}", f"mu_{rand}", lp.B_0.mass/1000., low, high)
sigma = ROOT.RooRealVar(f"sigma_{rand}", f"sigma_{rand}", 0.04, 0.025, 0.2)
alphaL = ROOT.RooRealVar(f"alphaL_{rand}", f"alphaL_{rand}", 0.2, 0., 5.)
nL = ROOT.RooRealVar(f"nL_{rand}", f"nL_{rand}", 10., 0., 200.)
alphaR = ROOT.RooRealVar(f"alphaR_{rand}", f"alphaR_{rand}", -0.2, -5., 0.)
nR = ROOT.RooRealVar(f"nR_{rand}", f"nR_{rand}", 10., 0., 200.)
# Core gaussian parameters
muG = ROOT.RooRealVar(f"muG_{rand}", f"muG_{rand}", lp.B_0.mass/1000., low1, high1)
sigmaG = ROOT.RooRealVar(f"sigmaG_{rand}", f"sigmaG_{rand}", 0.015, 0., 1)
# Fractions
frac = ROOT.RooRealVar(f"frac_{rand}", f"frac_{rand}", 0.5, 0., 1.)
frac_CB = ROOT.RooRealVar(f"frac_CB_{rand}", f"frac_CB_{rand}", 0.5, 0., 1.)
#pdfs
CBR = ROOT.RooCrystalBall(f"CBR_{rand}",f"CBR_{rand}",mass,mu,sigma,alphaR, nR, doubleSided=False)
CBL = ROOT.RooCrystalBall(f"CBL_{rand}",f"CBL_{rand}",mass,mu,sigma,alphaL, nL, doubleSided=False)
Gau = ROOT.RooGaussian(f"gau_{rand}",f"gau_{rand}",mass,muG,sigmaG)
sig_a = ROOT.RooAddPdf( f"sig_a_{rand}",f"sig_a_{rand}",CBL,CBR,frac_CB)
sig_fit= ROOT.RooAddPdf( f"sig_fit_{rand}",f"sig_fit_{rand}",sig_a, Gau,frac)

sig_fit.fitTo(signal_mass)
c22=ROOT.TCanvas()
Mframe = mass.frame()
signal_mass.plotOn(Mframe)
sig_fit.plotOn(Mframe)
Mframe.Draw()
c22.Draw('mass')
3

Maybe @jonas can help

4

Hi @TristanM, welcome to the ROOT forum!

I don’t know what’s going on here yet. Unfortunately I’m not able to reproduce the issue, can you maybe help me with that that?

I have created this notebook where I copy-pasted your code into:
ReproducibilityIssue.ipynb.txt (5.8 KB)
(remove the .txt suffix, which I had to add because the forum doesn’t allow me to upload the file otherwise).

I used ROOT 6.26.10 to run the cells in the order you described, but I didn’t get the problem. Maybe you can create a reproducer for the problem that you send back to me so I can start debugging?

Thanks! Jonas

5

Hi @jonas, thanks for your answer !

I played a bit with your notebook, and it’s true, no problem occurs.
So I modify it by using as input the data I try to fit:
(here I would like to share my data.txt but this is not allowed to me …) .

Initially the fit as defined in your notebook run indefinitely on my data (and no output are print during execution, I have the feeling that ROOT print his output when the run of the cell is ended in a notebook), but by replacing the B0 mass you have fixed by the one from the particle package (e.g. +0.01) I manage to reproduce my issue with this notebook:
(there I would like to share an updated version of your notebook but this is not allowed too …) .

You can simply run first all the code and then rerun the two last.

Thanks in advance.

6

Is there a way to share file when we are new ?

7

I have the ability now, let me repost my previous post @jonas :

I played a bit with your notebook, and it’s true, no problem occurs.
So I modify it by using as input the data I try to fit:
data.txt (895,7 Ko) .

Initially the fit as defined in your notebook run indefinitely on my data (and no output are print during execution, I have the feeling that ROOT print his output when the run of the cell is ended in a notebook), but by replacing the B0 mass you have fixed by the one from the particle package (e.g. +0.01) I manage to reproduce my issue with this notebook:
ReproducibilityIssue.ipynb.txt (5,1 Ko) .

You can simply run first all the code and then rerun the two last cells.

Thanks in advance.

8

Hi @TristanM, thanks for sharing the reproducer code and data, this is very much appreciated!

I can reproduce the problem now, at least a similar problem. Indeed, when I run the notebook the first time, the last cell with the fit takes 6.5 seconds. If I rerun the cell with the model definition and then do the fit again, the fit takes 21.6 seconds! So I suppose you don’t need to restart the kernel, just wait longer…

Anyway, it’s weird that the same thing takes more than three times as long when you rerun the cell. I will investigate the problem next week.

9

I see that it has nothing to do with Python. I converted the reproducer to standalone C++ code:

const double b0_mass = 5279.66;

int rand = 0;

RooRealVar mass{"mass", "mass", 4.6, 6.0};
std::unique_ptr<RooDataSet> signal_mass{RooDataSet::read("data.txt", mass)};

double low1 = 5.26;
double high1 = 5.3;
double low = 4.6;
double high = 6;

for (std::size_t i = 0; i < 3; ++i) {

   rand += 1;

   auto suffix = [&](std::string const &s) { return s + "_" + std::to_string(rand); };

   // signal fit
   // crystal ball parameters
   RooRealVar mu{suffix("mu").c_str(), "", b0_mass / 1000.0, low, high};
   RooRealVar sigma{suffix("sigma").c_str(), "", 0.04, 0.025, 0.2};
   RooRealVar alphaL{suffix("alphaL").c_str(), "", 0.2, 0.0, 5.0};
   RooRealVar nL{suffix("nL").c_str(), "", 10.0, 0.0, 200.0};
   RooRealVar alphaR{suffix("alphaR").c_str(), "", -0.2, -5.0, 0.0};
   RooRealVar nR{suffix("nR").c_str(), "", 10.0, 0.0, 200.0};
   // Core gaussian parameters
   RooRealVar muG{suffix("muG").c_str(), "", b0_mass / 1000.0, low1, high1};
   RooRealVar sigmaG{suffix("sigmaG").c_str(), "", 0.015, 0.0, 1};
   // Fractions
   RooRealVar frac{suffix("frac").c_str(), "", 0.5, 0.0, 1.0};
   RooRealVar frac_CB{suffix("frac_CB").c_str(), "", 0.5, 0.0, 1.0};

   RooCrystalBall CBR{suffix("CBR").c_str(), "", mass, mu, sigma, alphaR, nR, false};
   RooCrystalBall CBL{suffix("CBL").c_str(), "", mass, mu, sigma, alphaL, nL, false};
   RooGaussian Gau{suffix("gau").c_str(), "", mass, muG, sigmaG};
   RooAddPdf sig_a{suffix("sig_a").c_str(), "", CBL, CBR, frac_CB};
   RooAddPdf sig_fit{suffix("sig_fit").c_str(), "", sig_a, Gau, frac};

   using namespace RooFit;
   RooAbsReal::clearEvalErrorLog();
   std::unique_ptr<RooFitResult> res{
      sig_fit.fitTo(*signal_mass, PrintLevel(-1), PrintEvalErrors(-1), Save(), Timer(true), Optimize(0))};

   res->Print();
}

I see that the fit result is different the first time compared to when you are calling fitTo again. Probably this relates to some global state of the minimizer that is not correctly reset. I will continue to investigate that.

But probably, the problem would not occur if your fit converges well. Are you aware that RooFit reports the fit as unsuccessful? Here is the output of my C++ code:

[#1] INFO:DataHandling -- RooDataSet::read: reading file data.txt
[#1] INFO:DataHandling -- RooDataSet::read: read 36687 events (ignored 0 out of range events)
[#0] WARNING:InputArguments -- The parameter 'nR_1' with range [0, 200] of the RooCrystalBall 'CBR_1' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'nL_1' with range [0, 200] of the RooCrystalBall 'CBL_1' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'sigmaG_1' with range [0, 1] of the RooGaussian 'gau_1' exceeds the safe range of (0, inf). Advise to limit its range.
[#1] INFO:Minimization -- Command timer: Real time 0:00:07, CP time 7.230
[#1] INFO:Minimization -- Session timer: Real time 0:00:07, CP time 7.230
[#1] INFO:Minimization -- Command timer: Real time 0:00:00, CP time 0.750
[#1] INFO:Minimization -- Session timer: Real time 0:00:08, CP time 7.980, 2 slices

  RooFitResult: minimized FCN value: -11813.7, estimated distance to minimum: 1.22137
                covariance matrix quality: Full matrix, but forced positive-definite
                Status : MINIMIZE=-1 HESSE=4

    Floating Parameter    FinalValue +/-  Error
  --------------------  --------------------------
              alphaL_1    2.2833e-01 +/-  3.54e-03
              alphaR_1   -3.2005e-01 +/-  4.15e-04
                frac_1    9.9244e-01 +/-  7.78e-03
             frac_CB_1    4.5868e-01 +/-  3.12e-04
                 muG_1    5.2600e+00 +/-  7.40e-05
                  mu_1    5.2793e+00 +/-  4.18e-05
                  nL_1    1.1021e+02 +/-  4.08e-02
                  nR_1    2.0154e+00 +/-  3.75e-03
              sigmaG_1    9.5086e-02 +/-  2.74e-02
               sigma_1    3.5422e-02 +/-  2.83e-05

[#0] WARNING:InputArguments -- The parameter 'nR_2' with range [0, 200] of the RooCrystalBall 'CBR_2' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'nL_2' with range [0, 200] of the RooCrystalBall 'CBL_2' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'sigmaG_2' with range [0, 1] of the RooGaussian 'gau_2' exceeds the safe range of (0, inf). Advise to limit its range.
[#1] INFO:Minimization -- Command timer: Real time 0:00:25, CP time 25.160
[#1] INFO:Minimization -- Session timer: Real time 0:00:25, CP time 25.160
[#1] INFO:Minimization -- Command timer: Real time 0:00:00, CP time 0.820
[#1] INFO:Minimization -- Session timer: Real time 0:00:26, CP time 25.980, 2 slices

  RooFitResult: minimized FCN value: -11823.9, estimated distance to minimum: 5.39854
                covariance matrix quality: Full matrix, but forced positive-definite
                Status : MINIMIZE=-1 HESSE=4

    Floating Parameter    FinalValue +/-  Error
  --------------------  --------------------------
              alphaL_2    2.0316e-01 +/-  1.04e-05
              alphaR_2   -2.7724e-01 +/-  1.14e-02
                frac_2    9.4489e-01 +/-  1.09e-02
             frac_CB_2    4.6144e-01 +/-  5.14e-03
                 muG_2    5.2986e+00 +/-  2.58e-02
                  mu_2    5.2780e+00 +/-  1.10e-03
                  nL_2    1.1180e+02 +/-  3.00e-03
                  nR_2    1.9739e+00 +/-  4.43e-02
              sigmaG_2    7.4775e-02 +/-  5.63e-02
               sigma_2    3.1501e-02 +/-  2.58e-04

[#0] WARNING:InputArguments -- The parameter 'nR_3' with range [0, 200] of the RooCrystalBall 'CBR_3' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'nL_3' with range [0, 200] of the RooCrystalBall 'CBL_3' exceeds the safe range of (0, inf). Advise to limit its range.
[#0] WARNING:InputArguments -- The parameter 'sigmaG_3' with range [0, 1] of the RooGaussian 'gau_3' exceeds the safe range of (0, inf). Advise to limit its range.
[#1] INFO:Minimization -- Command timer: Real time 0:00:28, CP time 28.800
[#1] INFO:Minimization -- Session timer: Real time 0:00:28, CP time 28.800
[#1] INFO:Minimization -- Command timer: Real time 0:00:00, CP time 0.760
[#1] INFO:Minimization -- Session timer: Real time 0:00:29, CP time 29.560, 2 slices

  RooFitResult: minimized FCN value: -11823.9, estimated distance to minimum: 5.39854
                covariance matrix quality: Full matrix, but forced positive-definite
                Status : MINIMIZE=-1 HESSE=4

    Floating Parameter    FinalValue +/-  Error
  --------------------  --------------------------
              alphaL_3    2.0316e-01 +/-  1.04e-05
              alphaR_3   -2.7724e-01 +/-  1.14e-02
                frac_3    9.4489e-01 +/-  1.09e-02
             frac_CB_3    4.6144e-01 +/-  5.14e-03
                 muG_3    5.2986e+00 +/-  2.58e-02
                  mu_3    5.2780e+00 +/-  1.10e-03
                  nL_3    1.1180e+02 +/-  3.00e-03
                  nR_3    1.9739e+00 +/-  4.43e-02
              sigmaG_3    7.4775e-02 +/-  5.63e-02
               sigma_3    3.1501e-02 +/-  2.58e-04

If always says MINIMIZE=-1 HESSE=4, which means there was a problem. If things go well, the status codes will be zero. Probably the issue only happens if the fit doesn’t converge, and that’s why if was depending of the input data whether it can be reproduced.

10

Thanks, ok I see this strange behavior seems to come first from the fact that my model is not good enough to fit my data.

11

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.