While @moneta is on vacation: could e.g. @vcroft have a look?
Thanks @Axel. Just checking if it is reproducible would already be great as I cannot rule out something wrong on my end.
Issue is reproducible in SWAN with JupyROOT 6.12/06, RooFit v3.60, and Python 3.
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.
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.
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.
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:
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()
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 andFCONV(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…)
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.
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. 
Your post just crossed an update I made to my earlier one with a plot that I hope is a better illustration.
This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.
- maybe 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
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. 
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.
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).
Hello,
After fixing some bugs in FFT, I come back to this one.
I have looked at the Landau case, and this is related to a problematic definition we have in the Landau distribution.
We have an implementation of the Landau(x) as defined in the original definition,
first equation of Landau distribution - Wikipedia and implemented from code coming from the CERNLIB.
See ROOT: Probability Density Functions (PDF)
In that definition there is no location and scale parameters defined and we artificially introduce a location parameter in the code as x' = (x-m)/s.
Now this is different than the correct definition of the family of Landau distributions described as p(x;mu,c) in Landau distribution - Wikipedia.
Using that correct definition is true that:
CONV(Landau(x,m,c), Landau(x,m,c)) = Landau(x, 2*m, 2*c)
but it is not true that
CONV( Landau( (x-m)/s), Landau((x-m)/s) ) = Landau ( (x - 2*m)/s )
This explains what is observed in the plot above running the given macro TwoLandau.py gives difference results for the PDF between twomip and onetheomip.
If the same macro is applied to another stable distribution, like the Gaussian, see TwoGaus.py (4.1 KB) ,
the result is what expected, confirming that there is no problem on the RooFit side when performing an FFT convolution.
Lorenzo