2D fitting with python functions in pyROOT

Hi,

I am trying to fit 2D histogram with 2D Gaussian function defined in python with pyROOT.
Using Eval, I am able to display the TF2 from the python function. However, I am unable to fit with the TF2.How can I do this?

My code:

import numpy as np
import ROOT

class Gausstwo:
    def __call__(self, arr, par):
        z = par[0]*np.exp(-(arr[0]-par[1])**2/(2*par[2]**2))*np.exp(-(arr[1]-par[3])**2/(2*par[4]**2))
        return z

def main():
    ROOT.EnableImplicitMT()
    ROOT.gROOT.SetBatch()
    c2 = ROOT.TCanvas("c", "c")

    h2 = ROOT.TH2F("hpxpy", "py vs px", 40, -4., 4., 40, -4., 4.)
    h2.SetStats(0)
    px, py = np.float32(), np.float32()
    for i in range(50000):
        ROOT.gRandom.Rannor(px, py)
        h2.Fill(px, py)

    pyg = Gausstwo()
    func = ROOT.TF2("func", pyg, -4., 4., -4., 4., 5)
    # func = ROOT.TF2("func", '[0]*exp(-pow((x-[1]),2)/(2*pow([2],2)))*exp(-pow((y-[3]),2)/(2*pow([4],2)))', -4., 4., -4., 4.)
    func.SetParameter(0, 320.)
    func.SetParameter(1, 0.001)
    func.SetParameter(2, 1)
    func.SetParameter(3, 0.001)
    func.SetParameter(4, 1.)
    
    # print(func.Eval(0., 2.,))
    h2.Fit(func, "V")

    ROOT.gStyle.SetOptFit()
    h2.Draw("surf0")
    func.Draw("same surf")
    c2.Draw()
    c2.Print("fittest.png")


if __name__ == '__main__':
    main()

The output is here.(USERNAME is my user name.):

 **********
 **    1 **SET PRINT           2
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 p0           3.20000e+02  9.60000e+01     no limits
     2 p1           1.00000e-03  3.00000e-04     no limits
     3 p2           1.00000e+00  3.00000e-01     no limits
     4 p3           1.00000e-03  3.00000e-04     no limits
     5 p4           1.00000e+00  3.00000e-01     no limits
 **********
 **    3 **SET ERR           1
 **********
 **********
 **    4 **SET PRINT           2
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD        1625        0.01
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
Assertion failed: (*I)->isCompleted() && "Nested transaction not completed!?", file C:\Users\USERNAME\root\root-6.26.02\interpreter\cling\lib\Interpreter\IncrementalParser.cpp, line 511
 *** Break *** Assertion failed: (*I)->isCompleted() && "Nested transaction not completed!?", file C:\Users\USERNAME\root\root-6.26.02\interpreter\cling\lib\Interpreter\IncrementalParser.cpp, line 511
abort

ROOT Version: 6.26/02
Platform: Windows 11
Compiler: Not Provided
using Visual Studio 2022, python 3.10.2 (32bit)

I tried the same code on Ubuntu-20.04 (WSL2).The terminal seems to stop working. There, the python version is 3.8.10 and the ROOT version is 6.26/02.
The output is here.:

 **********
 **    1 **SET PRINT           2
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 p0           3.20000e+02  9.60000e+01     no limits
     2 p1           1.00000e-03  3.00000e-04     no limits
     3 p2           1.00000e+00  3.00000e-01     no limits
     4 p3           1.00000e-03  3.00000e-04     no limits
     5 p4           1.00000e+00  3.00000e-01     no limits
 **********
 **    3 **SET ERR           1
 **********
 **********
 **    4 **SET PRINT           2
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD        1625        0.01
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.

Thanks


Hi @happyisland ,

sorry for the high latency, this must have slipped through the cracks. We need @moneta here, let’s ping him.

Cheers,
Enrico

Hi,
It seems copy get stacks calling the function. I have tried also defining the Python function as a free function and it is the same. It is a PyROOT issue, can you @etejedor have a look ?
As a workaround you should use the function defined as a TFormula, that is anyway faster than the Python function.

Cheers

Lorenzo

I think it might be an MT problem (deadlock?). If you remove the line ROOT.EnableImplicitMT() your reproducer finishes, but with that line it gets stuck. @moneta knows better how MT is handled here.

It could be that using Python callables for your function requires getting the GIL every time the function is called, so there is no real MT and a lot of contention.

Worker threads indeed take the GIL:

Thread 20 (Thread 0x7fffb5921640 (LWP 43426) "python"):
#0  __futex_abstimed_wait_common64 (private=0, cancel=true, abstime=0x7fffb59204c0, op=137, expected=0, futex_word=0x7ffff7dafa48 <_PyRuntime+424>) at futex-internal.c:57
#1  __futex_abstimed_wait_common (futex_word=futex_word@entry=0x7ffff7dafa48 <_PyRuntime+424>, expected=expected@entry=0, clockid=clockid@entry=1, abstime=abstime@entry=0x7fffb59204c0, private=private@entry=0, cancel=cancel@entry=true) at futex-internal.c:87
#2  0x00007ffff768917f in __GI___futex_abstimed_wait_cancelable64 (futex_word=futex_word@entry=0x7ffff7dafa48 <_PyRuntime+424>, expected=expected@entry=0, clockid=clockid@entry=1, abstime=abstime@entry=0x7fffb59204c0, private=private@entry=0) at futex-internal.c:139
#3  0x00007ffff768bc24 in __pthread_cond_wait_common (abstime=0x7fffb59204c0, clockid=1, mutex=0x7ffff7dafa50 <_PyRuntime+432>, cond=0x7ffff7dafa20 <_PyRuntime+384>) at pthread_cond_wait.c:503
#4  ___pthread_cond_timedwait64 (cond=0x7ffff7dafa20 <_PyRuntime+384>, mutex=0x7ffff7dafa50 <_PyRuntime+432>, abstime=0x7fffb59204c0) at pthread_cond_wait.c:643
#5  0x00007ffff7b234ec in ?? () from /usr/lib/libpython3.10.so.1.0
#6  0x00007ffff7b43976 in PyEval_RestoreThread () from /usr/lib/libpython3.10.so.1.0
#7  0x00007ffff7c175fd in PyGILState_Ensure () from /usr/lib/libpython3.10.so.1.0
#8  0x00007fffce275220 in __cppyy_internal::fptr_wrap1 (arg0=<optimized out>, arg1=<optimized out>) at input_line_47:11
#9  0x00007fffb677a3c5 in ROOT::Math::ParamFunctorHandler<ROOT::Math::ParamFunctorTempl<double>, double (*)(double*, double*)>::FuncEvaluator<double (*)(double*, double*), double>::Eval (f=0x7fffce275000 <__cppyy_internal::fptr_wrap1(double*, double*)>, x=0x7fffac001340, p=0x5555668b7140) at ../math/mathcore/inc/Math/ParamFunctor.h:120
#10 0x00007fffb6779f07 in ROOT::Math::ParamFunctorHandler<ROOT::Math::ParamFunctorTempl<double>, double (*)(double*, double*)>::operator() (this=0x55555c8ef0f0, x=0x7fffac001340, p=0x5555668b7140) at ../math/mathcore/inc/Math/ParamFunctor.h:90
#11 0x00007fffb6776d89 in ROOT::Math::ParamFunctorTempl<double>::operator() (this=0x55555dd2e688, x=0x7fffac001340, p=0x5555668b7140) at ../math/mathcore/inc/Math/ParamFunctor.h:362
#12 0x00007fffb6766b1b in TF1::EvalPar (this=0x555565705710, x=0x7fffac001340, params=0x0) at ../hist/hist/src/TF1.cxx:1500
#13 0x00007fffb672cb78 in ROOT::Math::WrappedMultiTF1Templ<double>::DoEval (this=0x5555670f9be0, x=0x7fffac001340) at ../hist/hist/inc/Math/WrappedMultiTF1.h:165
#14 0x00007fffb7646946 in ROOT::Math::IParametricFunctionMultiDimTempl<double>::operator() (this=0x5555670f9be0, x=0x7fffac001340) at ../math/mathcore/inc/Math/IParamFunction.h:127
#15 0x00007fffb763673c in operator() (__closure=0x7fffffffbf90, i=502) at ../math/mathcore/src/FitUtil.cxx:318
#16 0x00007fffb763e0d7 in operator() (__closure=0x5555688932a0, i=498) at ../core/imt/inc/ROOT/TThreadExecutor.hxx:389

EDIT: at a quick glance I can’t tell whether the thread is deadlocked or not, but this definitely deserves a look

Thank you for your replies.

I removed ROOT.EnableImplicitMT() from my code and it work well. I don’t know much about GIL but it seems MT needs to be disabled when using the Python function.

Thank you Enric, we should then disable MT when using Python functions. If MT is enables one can still use the Serial fitting option.

@happyisland, a workaround, in case you need to have enabled MT for other reasons, is to do the fit as following:

h2.Fit(func, "SERIAL V")

Lorenzo