Running out of memory to calculate multi-dimensional integral

ROOT
1

I am calculating the elements of a matrix, and each element involves the integration of multi-dimensional functions. I have attached a part of my code to calculate row number 6 using ROOT::Math::AdaptiveIntegratorMultiDim. I would like to have an accuracy better than 1%, however, when I choose the relative tolerance to 1.0e-2, I run out of memory. Is there another integration method in ROOT that can give me good accuracy without running out of memory? I am aware of a ML method published in Nature Scientific Reports 11, 18965 (2021), but as far as I know it was not yet implemented in ROOT.

double f(double x) {
    string filet="xyscan-ExpData.txt";
    TString flt(filet);
        
    TGraph *gr1 = new TGraph(flt);
    gr1->SetMarkerStyle(27);
    TF1 *tf1 = new TF1("tf1","([0]**2*[1])/([1]-[0])**2*(exp(-[0]*x)*(([1]-[0])*x-1)+exp(-[1]*x))*[2]+[3]",0.0,100.0); //[0] is lambda and [1] stands for beta.
    tf1->SetParameters(0.154,0.0722,100.0,0.0001,0.0);
    gr1->Fit("tf1","QR","",0,100);
    Double_t top=tf1->Eval(tf1->GetMaximumX())+100.0;
    Double_t T=dt(1);
    gr1->GetYaxis()->SetRangeUser(0,top);
    gr1->GetXaxis()->SetRangeUser(0,T);
    TF1 *tf2 = new TF1("tf2","tf1-[3]",0.0,100.0);
    TF1 *tf4 = new TF1("tf4","tf2/([0])",0.0,T);
    tf4->SetParameter(0,tf2->Integral(0.0,T));

    Double_t p0=tf1->GetParameter(0), p1=tf1->GetParameter(1), p2=tf1->GetParameter(2), p3=tf2->Integral(0.0,T);
    
    return ((pow(p0,2)*p1/pow(p1-p0,2))*(exp(-p0*x)*((p1-p0)*x-1.0)+exp(-p1*x))*p2)/p3;
}

//____________________________

double fint1(double t)
{
    ROOT::Math::Functor1D wf(&f);
    ROOT::Math::Integrator ig(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR);
    ig.SetFunction(wf);
    Double_t tau=dt(0);
    Double_t fintval = ig.Integral(t,t+tau);
    
    return fintval;
}
//_______________________________i=6-------------------------------------------------------

double ff6(const double *x, const double *par)
{
    const int i=6;
    Double_t tau=0.9;
    Double_t T=60.0;
    Double_t t2min=x[0]+tau;
    Double_t t2max=T-(i-2.0)*tau;
    Double_t t2=(t2max-t2min)*x[1]+t2min;
    Double_t dt2=t2max-t2min;
    Double_t t3min=t2+tau;
    Double_t t3max=T-(i-3.0)*tau;
    Double_t t3=(t3max-t3min)*x[2]+t3min;
    Double_t dt3=t3max-t3min;
    Double_t t4min=t3+tau;
    Double_t t4max=T-(i-4.0)*tau;
    Double_t t4=(t4max-t4min)*x[3]+t4min;
    Double_t dt4=t4max-t4min;
    Double_t t5min=t4+tau;
    Double_t t5max=T-(i-5.0)*tau;
    Double_t t5=(t5max-t5min)*x[4]+t5min;
    Double_t dt5=t5max-t5min;
    Double_t t6min=t5+tau;
    Double_t t6max=T-(i-6.0)*tau;
    Double_t t6=(t6max-t6min)*x[5]+t6min;
    Double_t dt6=t6max-t6min;
    Double_t ffval = (TMath::Factorial(par[0])/TMath::Factorial(par[0]-i)) *f(x[0]) *f(t2) *f(t3) *f(t4) *f(t5) *f(t6) *dt2 *dt3 *dt4 *dt5 *dt6 *pow(fint1(x[0]) +fint1(t2) +fint1(t3) +fint1(t4) +fint1(t5) +fint1(t6),par[0]-i);
    return ffval;
}

//____________________________

double S6j(double j)
{
    const int i=6;
    Double_t tau=0.9;
    Double_t T=60.0;
    Double_t epsrel=1e-2;
    ROOT::Math::WrappedParamFunction<> wf( &ff6, i, 1);
    Double_t p[1] = {j};
    wf.SetParameters(p);
    // Create the Integrator
    ROOT::Math::AdaptiveIntegratorMultiDim ig;
    // Set parameters of the integration
    ig.SetFunction(wf);
    ig.SetRelTolerance(epsrel);
    Double_t xmin[i] = {0.0,0.0,0.0,0.0,0.0,0.0};
    Double_t xmax[i] = {T-(i-1.0)*tau,1.0,1.0,1.0,1.0,1.0};
    Double_t Sjvalue =ig.Integral(xmin, xmax);
    return Sjvalue;
}
void spencerIntegrals() {
    ROOT::EnableImplicitMT(0);
    cout << "Matrix element S66: " << S6j(6) << endl;
    cout << "Matrix element S67: " << S6j(7) << endl;
}

ROOT Version: v6-28-06
Platform: Red Hat Enterprise Linux 8.8; Memory 125.1 GiB; Intel Core i7-9800x CPU 3.8 GHz x16.
Compiler: C++ (GCC) 8.5.0 20210514 (Red Hat 8.5.0-18)

2

Hi @garciaduarte1,

thank you for the question. Maybe @moneta could share some insights here?

Cheers,
Marta

3

Hi,

The adaptive integrator multi-dim is designed to work for low dimensions (e.g. 2,3,4), and not really for higher like 6. You can try to use the MC integrator based for example on VEGAS, it should use less memory and work better in this case. See the example code here:

However, the fact that is running out of memory is due to a massive memory leak in your code. You are
creating for each function evaluation TGraph and TF1 objects without deleting them. I would recommend to create these objects outside the function evaluation code, and re-use them every time you need to compute the function f(x).

Best regards

Lorenzo

1 Like
4

I appreciate your help @mczurylo.

Thank you for your answer @moneta.

I was able to fix the memory problem. However, the MC integrator with VEGAS takes forever to converge even for a relative tolerance of 1.0 or 0.1 in the case of a 5-dimensional integration. I still think the Adaptive Integrator works better in my case. From the example @moneta attached, it seems I just need to replace

ROOT::Math::AdaptiveIntegratorMultiDim ig;

by

ROOT::Math::GSLMCIntegrator ig(ROOT::Math::MCIntegration::kVEGAS);

While the Adaptive Integrator takes about an hour to converge with a tolerance of 0.01, the MC integrator with VEGAS was still running after 3 hours. Maybe part of the problem is that the Adaptive Integrator runs well in multi-thread and apparently the MC Integrator is not made to run in multi-thread. Do you have other thoughts on this @moneta?

Thank you and best regards,
Jeremias

5

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