TF1::Integral and TF1::GetQuantiles problems

Hi,

Apparently I am having some issues computing integral and quantiles of simple functions (a quarter of a circle). The bug is related to a storage area which is overwritten. It appears from time to time and shows the following :

Fatal in <operator delete>: storage area overwritten
aborting
#0  0x00007fe6525cf4ca in __GI___waitpid (pid=9486, stat_loc=stat_loc
entry=0x7fffbf55f930, options=options
entry=0) at ../sysdeps/unix/sysv/linux/waitpid.c:29
#1  0x00007fe652548fbb in do_system (line=<optimised out>) at ../sysdeps/posix/system.c:148
#2  0x00007fe6587259bd in TUnixSystem::StackTrace() () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#3  0x00007fe65863ff6a in DefaultErrorHandler(int, bool, char const*, char const*) () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#4  0x00007fe65863f902 in ErrorHandler () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#5  0x00007fe65863fdff in Fatal(char const*, char const*, ...) () from /opt/cern/root/root_v6.10.02/lib/libCore.so
#6  0x00007fe6553db180 in operator delete(void*) () from /opt/cern/root/root_v6.10.02/lib/libNew.so
#7  0x00007fe65799410e in TF1::IntegralOneDim(double, double, double, double, double&) () from /opt/cern/root/root_v6.10.02/lib/libHist.so
#8  0x00007fe657987479 in TF1::Integral(double, double, double) () from /opt/cern/root/root_v6.10.02/lib/libHist.so
#9  0x0000000000402546 in main (argc=<optimised out>, argv=<optimised out>) at main.cpp:105
Aborted (core dumped)

Whereas when the output is correct, it gives the lines below where something looks wrong in the TF1::GetQuantiles function :

Curve length is : 9.42478 units
Probas:
0  0.125
1  0.25
2  0.375
3  0.5
4  0.625
5  0.75
6  0.875
7  1
Error in <TArrayD::operator[]>: index 10000 out of bounds (size: 10000, this: 0x7fffc91c0330)
Error in <TArrayD::operator[]>: index 10001 out of bounds (size: 10001, this: 0x7fffc91c0310)
Quantiles:
0  0.125  2.11529
1  0.25  2.45672
2  0.375  3.01118
3  0.5  3.75736
4  0.625  4.66658
5  0.75  5.7039
6  0.875  6.82946
7  1  8
ROfZ at quantiles:
0  2.11529  34.8295
1  2.45672  33.7039
2  3.01118  32.6666
3  3.75736  31.7574
4  4.66658  31.0112
5  5.7039  30.4567
6  6.82946  30.1153
7  8  30
Interval length:
1  2.45672  1.1781
2  3.01118  1.1781
3  3.75736  1.1781
4  4.66658  1.1781
5  5.7039  1.1781
6  6.82946  1.1781
7  8  1.1781

The simple and short program, compiled on Ubuntu 16.04 LTS using root 6.10/02 Built for linuxx8664gcc From tag v6-10-02, 6 July 2017 is the following :

//-----------------------------------------------------------------------------
//-----------------------------------------------------------------------------
#include <algorithm>
#include <cmath>
#include <ctype.h>
#include <float.h>      // required for LDBL_EPSILON, DBL_MAX
#include <fstream>
#include <iomanip>
#include <iostream>
#include <math.h>       // required for fabs(), fabsl(), sqrtl(), and M_PI_2
#include <sstream>
#include <stdexcept>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <utility>      // std::pair, std::make_pair
#include <vector>

#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_ellint.h>


// ROOT (http://root.cern.ch) specific files
#include "TROOT.h"
#include "Riostream.h"
#include "TArrayD.h"
#include "TCanvas.h"
#include "TColor.h"
#include "TF1.h"
#include "TGraph.h"
#include "TH1D.h"
#include "TNamed.h"
#include "TObject.h"
#include "TRint.h"
#include "TString.h"
#include "TMath.h"
#include "Math/IFunction.h"
#include "Math/SpecFuncMathMore.h"
#include "Math/WrappedTF1.h"
#include "Math/GSLIntegrator.h"


// Use functions from namespace std
using namespace std;


//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
//-------------------------------------------------------------------------------
TROOT root("Rint", "The ROOT Interactive Interface");
int main(int argc, char **argv)
{
  
  TRint theApp("Electrobem", NULL, NULL);
  gROOT->Reset();

  Int_t icell = 0;

  Double_t zmin = 2;
  Double_t zmax = 8;

  const Int_t nz = 8;
  Double_t dz[nz+1], PointAlongZ[nz+1];


  TCanvas *cnv = new TCanvas("cnv", "cnv", 700., 700.);


  TF1 * ROfZ = new TF1("ROfZ","36-sqrt(36-(x-8)*(x-8))", zmin, zmax);
  ROfZ->SetNpx(2000);
  ROfZ->SetLineColor(kRed);
  ROfZ->Draw();


  TF1 * SOfZ = new TF1("SOfZ","sqrt(1+(x-8)*(x-8)/(36-(x-8)*(x-8)))/[0]", zmin, zmax);
  SOfZ->SetParameter(0, 1.);
  SOfZ->SetNpx(10000);// useful for integral !!

  
  Double_t CurveL = SOfZ->Integral(zmin, zmax, 1.e-07);
  cout << "Curve length is : " << CurveL << " units" << endl;
  SOfZ->SetParameter(0, CurveL); // PDF normalization
  SOfZ->SetNpx(10000);


  
  // Number of desired equidistant point along the curve
  Double_t q[nz], pr[nz], foq[nz];

  
  // Array of probabilities in which to divide Y axis before taking
  // reciproqual of the distribution called quantiles
  cout << "Probas:" << endl;
  for (Int_t ipr = 0; ipr < nz; ++ipr){
    pr[ipr] = (ipr + 1) / ((double)nz);
    cout << ipr << "  " << pr[ipr] << endl;
  }

  // Quantiles determination
  SOfZ->GetQuantiles(nz, q, pr); // does not like the quantile for xq=1
  q[nz-1] = zmax;

  cout << "Quantiles:" << endl;
  for (Int_t ipr = 0; ipr<nz; ipr++){
    cout << ipr << "  " << pr[ipr] << "  " << q[ipr] << endl;
  }

  cout << "ROfZ at quantiles:" << endl;
  for (Int_t ipr = 0; ipr < nz; ++ipr){
    foq[ipr] = ROfZ->Eval(q[ipr]);
    cout << ipr << "  " << q[ipr] << "  " << foq[ipr] << endl;
  }

  cout << "Interval length:" << endl;
  SOfZ->SetParameter(0, 1.);
  for (Int_t ipr = 1; ipr < nz; ++ipr){
    cout << ipr << "  " << q[ipr] << "  " << SOfZ->Integral(q[ipr-1], q[ipr]) << endl;
  }
    
  TGraph *gr = new TGraph(nz, q, foq);
  gr->SetLineColor(1);
  gr->SetLineWidth(1);
  gr->SetMarkerColor(1);
  gr->SetMarkerSize(1);
  gr->SetMarkerStyle(20);
  gr->Draw("PLsame");

  
  theApp.Run();
  return 0;
}

I have not been able to find a bug on my side, i.e. in the above lines, therefore I am wondering if something is wrong in the functions TF1::Integral() and TF1::GetQuantiles().
Any hint would be helpful.
Thanks a lot,
Gilles

Hi Gilles,
Looking at your code seems that the problem is due to something in the GetQuantile method of TF1.
The error is for for index equal to 10000 which is the number of points of your function.

Maybe if you change this line with

with this one

SOfZ->GetQuantiles(nz- 1, q, pr);

In this way you avoid the calculation of the quantile for the last point, that in any case you do by hand.

Cheers
Stefano

Hi Stefano,

Thanks for your answer, now only remains the TF1::Integral crash from time to time.
Your trick to change nz in nz-1 solves the TF1::GetQuantiles problem.
Many thanks,

Gilles

Hi Gilles,
maybe the integral crash is due to the behaviour of the function close to 2.
In 2 the derivative of the function goes to infinite and this can create problem to the integral.

Maybe if your zmin instead of 2 is 2 + 0.001 or something similar you should achieve similar results without the Integral crash.

Stefano

Hi Stefano,

Thanks for this suggestion, but unfortunately this does not solve the problem.
Even when setting zmin to 3 and zmax to 5 the crash still happens from time to time.

Gilles

Hi,

The analysis of Stefano is correct (Thanks Stefano!).
There is a bug in TF1::GetQuantile for probability=1. I will fix this now. However, I am not sure this caused the Fatal error you are getting, since the bounds are checked in TArrayD.

Lorenzo

1 Like

Hi Lorenzo,

Thanks for correcting this bug and for Stefano’s suggestion !
The 2 encountered problems are indeed not correlated at all. I have only performed some simple integral and get this fatal error. Perhaps it is due to some wrong compilation of ROOT on our systems ? In order to test that I will try to find some other ROOT versions in the lab.

Gilles

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