Chained PDF definitions

adavid · July 19, 2018, 1:15am

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?

Axel · July 19, 2018, 10:10am

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

adavid · July 19, 2018, 10:33am

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

adavid · July 19, 2018, 11:56am

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

adavid · July 19, 2018, 12:20pm

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.

vcroft · July 20, 2018, 1:01pm

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.

adavid · July 20, 2018, 2:59pm

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.

vcroft · July 26, 2018, 3:54pm

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:

vcroft · July 26, 2018, 4:03pm

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()

adavid · July 26, 2018, 9:00pm

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):

image.png1252×848 125 KB

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…)

adavid · July 26, 2018, 9:15pm

vcroft:

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()

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.

vcroft · July 26, 2018, 10:44pm

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.

adavid · July 26, 2018, 11:01pm

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

system · August 9, 2018, 11:01pm

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

Axel · October 18, 2018, 12:52pm

With @moneta being on vacation - again - maybe @StephanH can have a look?

StephanH · October 19, 2018, 9:13am

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

adavid · October 19, 2018, 1:10pm

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.

StephanH · March 20, 2019, 6:20pm

Hello @adavid,

I’m coming back to this now. I found how to fix a bug with FFT convolutions, and now I would like to test if this has an effect on what we were discussing here. Could you send me the code that you used to plot these three coloured examples post #10?

The fix could potentially solve the Conv(x,y) != Conv(2x, 2y) mystery, and maybe the commutativity could be fixed as well.

adavid · March 20, 2019, 7:54pm

Good news, thanks!
I had to check if the code was still running and with ROOT 6.16.00 not everything is pretty, but one should get the gist of it from TwoLandau.py (4.0 KB).