Home | News | Documentation | Download

Chained PDF definitions


#1

Hi,

I’m trying to define PDFs as a function of others in a chain, using JupyROOT 6.14/00, RooFit v3.60, and Pyhton 3.6, on OS X 10.13.6:

import ROOT as r

w = r.RooWorkspace()
onemip = w.factory('Landau::onemip(x[0,500], mpv[40,0,100], sigma_landau[5,0,20])')
twomip = w.factory('FCONV::twomip(x, onemip, onemip )')
threemip = w.factory('FCONV::threemip(x, twomip, onemip )')
threemip.getVal()

mips = w.factory('SUM::mips(fone[0.1,0,1]*onemip, ftwo[0.3,0,1]*twomip, threemip)')
mips.getVal()
noise = w.factory('Gaussian::noise(x, mu[0,-1,1], sigma_noise[2,0,20])')
noise.getVal()
signal = w.factory('FCONV::signal(x, noise, mips)')
signal.getVal()

The last getVal() seems to trigger something bad:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-13-c4996a9ad9b2> in <module>()
----> 1 signal.getVal()

TypeError: none of the 2 overloaded methods succeeded. Full details:
  double RooAbsReal::getVal(const RooArgSet* set = 0) =>
    unhandled, unknown C++ exception
  double RooAbsReal::getVal(const RooArgSet& set) =>
    takes at least 1 arguments (0 given)[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'refrange_fft_signal' created with bounds [0,500]
[#0] ERROR:InputArguments -- RooAbsDataStore::initialize(twomip_fft_onemip_fft_CONV_onemip_fft_CACHEHIST_Obs[x_shifted_FFTBuffer2]_BufFrac0.1_BufStrat0): Data set cannot contain non-fundamental types, ignoring x_shifted_FFTBuffer2
[#0] ERROR:InputArguments -- RooHistPdf::ctor(onemip_fft_CONV_onemip_fft_CACHE_Obs[x]) ERROR histogram variable list and RooDataHist must contain the same variables.
[#0] ERROR:InputArguments -- RooAbsDataStore::initialize(twomip_fft_onemip_fft_CONV_onemip_fft_CACHEHIST_Obs[x_shifted_FFTBuffer2]_BufFrac0.1_BufStrat0): Data set cannot contain non-fundamental types, ignoring x_shifted_FFTBuffer2
[#0] ERROR:InputArguments -- RooHistPdf::ctor(onemip_fft_CONV_onemip_fft_CACHE_Obs[x]) ERROR histogram variable list and RooDataHist must contain the same variables.

Am I doing something wrong?


#2

While @moneta is on vacation: could e.g. @vcroft have a look?


#3

Thanks @Axel. Just checking if it is reproducible would already be great as I cannot rule out something wrong on my end.


#4

Issue is reproducible in SWAN with JupyROOT 6.12/06, RooFit v3.60, and Python 3.


#5

While trying a different strategy, I’m hitting the same issue:

import ROOT as r

w = r.RooWorkspace()

noise = w.factory('Gaussian::noise(x[0,500], mu[0], sigma_noise[2,0,20])')

onemip = w.factory('Landau::onemip(x, mpv[40,0,100], sigma_landau[5,0,20])')
onerealmip = w.factory('FCONV::onerealmip(x, noise, onemip )')
onerealmip.getVal()

twomip = w.factory('FCONV::twomip(x, onemip, onemip )')
tworealmip = w.factory('FCONV::tworealmip(x, noise, twomip )')
tworealmip.getVal()
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-9-49e4da1dd4f8> in <module>()
----> 1 tworealmip.getVal()

TypeError: none of the 2 overloaded methods succeeded. Full details:
  double RooAbsReal::getVal(const RooArgSet* set = 0) =>
    unhandled, unknown C++ exception
  double RooAbsReal::getVal(const RooArgSet& set) =>
    takes at least 1 arguments (0 given)[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'refrange_fft_tworealmip' created with bounds [0,500]
[#0] ERROR:InputArguments -- RooAbsDataStore::initialize(twomip_fft_onemip_fft_CONV_onemip_fft_CACHEHIST_Obs[x,x_shifted_FFTBuffer2]_BufFrac0.1_BufStrat0): Data set cannot contain non-fundamental types, ignoring x_shifted_FFTBuffer2
[#0] ERROR:InputArguments -- RooHistPdf::ctor(onemip_fft_CONV_onemip_fft_CACHE_Obs[x]) ERROR histogram variable list and RooDataHist must contain the same variables.
[#0] ERROR:InputArguments -- RooAbsDataStore::initialize(twomip_fft_onemip_fft_CONV_onemip_fft_CACHEHIST_Obs[x,x_shifted_FFTBuffer2]_BufFrac0.1_BufStrat0): Data set cannot contain non-fundamental types, ignoring x_shifted_FFTBuffer2
[#0] ERROR:InputArguments -- RooHistPdf::ctor(onemip_fft_CONV_onemip_fft_CACHE_Obs[x]) ERROR histogram variable list and RooDataHist must contain the same variables.

#6

Hi all I can reproduce. I’ll have a think, but at first glance it seems to be an issue with automatic object renaming. I’m flying today but I’ll get back to you when I’ve had a play with it.


#7

Thanks.

FWIW, and against what would be expected from a stable distribution like the Landau, the convolution of Landau(mpv, sigma) with Landau(mpv, sigma) is not yielding the same result as Landau(2mpv, 2sigma). I wonder if this is related to numerical issues with the Landau implementation in ROOT/RooFit or with something else.


#8

Hi there adavid,
I’m still puzzling over this. Especially working out what you’re trying to do. I got the example from the comments working (by swapping noise and twomip in the tworealmip definition) is this what you wanted/expected?

For example mipone with the convolution looks like this:


but the second convolution twomip:


#9

for me this works…

onemip = w.factory('Landau::onemip(x[0,500], mpv[40,0,100], sigma_landau[5,0,20])')
twomip = w.factory('FCONV::twomip(x, onemip, onemip )')
threemip = w.factory('FCONV::threemip(x, twomip, onemip )')
threemip.getVal()

mips = w.factory('SUM::mips(fone[0.1,0,1]*onemip, ftwo[0.3,0,1]*twomip, threemip)')
mips.getVal()
noise = w.factory('Gaussian::noise(x, mu[0,-1,1], sigma_noise[2,0,20])')
noise.getVal()
signal = w.factory('FCONV::signal(x, mips, noise)')
signal.getVal()

#10

Hi and thanks for looking into this:

  • I confirm that changing the order of the arguments seems to work. OTOH it’s rather unexpected (and frustrating) that FCONV(mips, noise) works and FCONV(noise, mips) doesn’t. I don’t think anyone can expect to not rely on convolution commutativity, though I understand computer science is not (just) mathematics.

  • As for what I am trying to do it’s pretty straightforward: modelling the signal of particles passing through thin silicon detector (Landaus) with noise (Gaussian).

  • Finally, I am still puzzling as to why in ROOT/RooFit the convolution of Landau(mpv, sigma) with Landau(mpv, sigma) does not yield the same result as Landau(2mpv, 2sigma) (see also the comment above).
    I looked a little more into this, tossed N pseudo-data from a Landau(mpv, sigma) and then summed pairs (or triplets) of values (data points). These are not at all described by Landau(2mpv, 2sigma) (or Landau(3mpv, 3sigma)) (dotted curves). OTOH, the convolutions of Landau(mpv, sigma) (full lines) seem to match the pseudo-data as expected (modulo some very likely numerical issues at the lowest values):

So, could the problem with dotted vs full be a numerical problem with the underlying Landau implementation? (The issues with the full lines at the lowest x values are certainly numerical…)


Landau distribution properties in RooFit (ROOT?)
#11

Just to note that this while this evaluates to something, it also yields an ERROR:

0.005754568929568781[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'refrange_fft_threemip' created with bounds [0,500]
[#0] ERROR:InputArguments -- RooArgSet::checkForDup: ERROR argument with name onemip is already in this set
[#1] INFO:Eval -- RooRealVar::setRange(x) new range named 'refrange_fft_twomip' created with bounds [0,500]
[#1] INFO:NumericIntegration -- RooRealIntegral::init(onemip_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(onemip_fft_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(twomip) creating new cache 0x7fbb711232d0 with pdf onemip_CONV_onemip_CACHE_Obs[x] for nset (x) with code 0
[#1] INFO:NumericIntegration -- RooRealIntegral::init(onemip_CONV_onemip_CACHE_Obs[x]_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(onemip_fft_Int[x]) using numeric integrator RooIntegrator1D to calculate Int(x)
[#1] INFO:Caching -- RooAbsCachedPdf::getCache(threemip) creating new cache 0x7fbb75340130 with pdf twomip_CONV_onemip_CACHE_Obs[x] for nset () with code 0

Every subsequent operation with threemip yields this pesky (useless?) ERROR.


#12

Hi again.
Indeed it’s a little unexpected but I expect it’s just a lack of function overloading for combinations that aren’t used so often in the factory (which is what the original error said - and still it took me forever to notice). Also the Errors are to do with the fact that the factory isn’t often used for more in depth studies such as these and therefore the errors are because the factory presumes that all of the distributions created will be new.

Whilst the factory is great for generating things quickly most people prefer to work with the pdf objects directly as can be seen in examples such as https://root.cern.ch/root/html/tutorials/roofit/rf708_bphysics.C.html that uses RooGaussian, RooFormulaVar, RooTruthModel etc for their convolutions. I’m not sure it helps but the factory can be a bit pernickety at times.

If I get a chance tomorrow I’ll look some more in to the underlying differences between the two landau implementations and see if I can dig up anything on the numerical issues. But indeed maybe we need Lorenzo’s expertise here. :wink:


#13

Your post just crossed an update I made to my earlier one with a plot that I hope is a better illustration.


#14

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


#15

#16

With @moneta being on vacation - again :wink: - maybe @Hageboeck can have a look?


#17

Hello,

alright - FFTs are a difficult topic. We likely are seeing numerical and logical (as in: roofit has a bug) issues here.

Bins

The first thing I want to recommend is to use a high number of bins for x. The factory uses 500, as I found out with

onemip = w.factory('Landau::onemip(x[0,500], mpv[40,0,100], sigma_landau[5,0,20])')
x = w.obj("x")
x.getBinning().Print()

Now, what you don’t see, unfortunately (because Doxygen for many RooFit classes is misconfigured improvements coming soon, at least in the master documentation), the numerical FFT needs a lot of bins for decent accuracy. See the comments at the top of:
https://root.cern.ch/doc/master/RooFFTConvPdf_8cxx_source.html

The factory doesn’t seem to have a way to communicate the number of bins, so you unfortunately have to

x = w.obj("x")
x.setBins(2000)

and please do this before evaluating anything, because this allocates the caches. I hope 2000 will do a good job.

Periodicity

The next issue might be the cyclical nature of FFTs. Since they are analysing frequencies, they see your input data as a grid with cyclical boundary conditions. After all, frequencies only make sense if something oscillates …
If you had only 3 bins, the FFT sees the bins as follows:

... 0 1 2 0 1 2 ...

Now, if the Landau has a tail to the right (because the range of x is not large enough), the FFT sees periodic jumps. This leads to artifacts in Fourier space and might well be the cause of the non-linearity of the FFT. You could try to choose the range of x large such that the dropoff at the high end of x’s value range is small. Also, try a generous buffer zone setBufferFraction(0.2) (default is 0.1), which extends the grid, on which the FFTs run. What you get is an x with an extended value range, this is FFT-ed, and artifacts due to boundary conditions should mostly affect the buffer zone. When evaluating the PDFs, the buffer zone is stripped away, so you are left with a better-behaved PDF.

Expert-level FFT

To avoid problems with the boundary conditions, one could multiply a window function on top of the Landaus before transforming them. A proper window would pull the Landaus to zero at the boundaries, removing all problems with the periodicity.
The poor-man’s solution is the generous buffer zone. Play around with it.

Whether 2000 bins and 0.2 will do the job, isn’t clear to me. I guess that these are good starting points, maybe overdoing it a bit. You might tweak them up or down a bit.

Speed

One more thing: FFTs like powers of 2. It would be better to use 2048 bins. However, since the buffer zone comes on top of that, you will be left with something like 245x bins. If you want to squeeze out the last bits of speed, you can play around with bin numbers around 1706.6 (2048/1.2). You might hit the magic number and get a bit more speed out of it.

Commutativity:

That must be a bug. I will try to investigate.

Cheers,
Stephan


#18

Hi Stephan,

Thanks for the hints on setBins and setBufferFraction. With those in, I confirm that the bad behaviour at zero goes away in my setup.

Then of course, the other issues are as you say: odd. :slight_smile: