Nested lambda functions overhead: ideas for algebra of functions

I have an implementation of algebraic operations (addition, multiplication) with use-defined 1d functions for the purpose of the LLH fitting.

The modular construction of PDFs is nice and handy, but I noticed a significant overhead (x10)
due to nested lambdas.

Are there any suggestions for possible improvement?

Here is a toy example on a spectrum I would like to fit.

The PDF is implemented below: it consists of mostly non-standard PDF components.

using AlgebraPDF #

# initial parameters of the model
const αv = 0.37
const βv = -0.002
const σev = [0.58, 1.19, 1.35, 1.58] .* 1e-3;
const mth = 0.96;
const Mv = [0.99, 1.02, 1.06, 1.11]
const Γv = [0.005, 0.001, 0.004, 0.010]
const fitlims = (0, 0.22);

pq(s) = s< mth^2 ? 0.0 : sqrt(s-mth^2)/sqrt(s)
BW(s,m,Γ) = 1/(m^2-s-1im*m*Γ)
BWmth(e,ΔM,Γ) = (e+mth)*Γ*BW((e+mth)^2, ΔM+mth, Γ)
pqmth(e) = pq((e+mth)^2)
bgd(e,α,β) = e^α*exp(-e*β)

# construction of PDFs
function construct_pdf_directly_integral_conv()
    function four_peaks_and_background(e; p)
               (conv_with_gauss.(e, y->abs2(BWmth(y, p.m1, p.Γ1)), σev[1]) .+
        p.f2.*conv_with_gauss.(e, y->abs2(BWmth(y, p.m2, p.Γ2)), σev[2]) .+
        p.f3.*conv_with_gauss.(e, y->abs2(BWmth(y, p.m3, p.Γ3)), σev[3]) .+
        p.f4.*conv_with_gauss.(e, y->abs2(BWmth(y, p.m4, p.Γ4)), σev[4]) .+
        p.fb.*bgd.(e, p.α, p.β)).*pqmth.(e)
    pdf_direct = 
        pdf(four_peaks_and_background; lims=fitlims, p0=(
            m1=Mv[1]-mth, Γ1=Γv[1],
            m2=Mv[2]-mth, Γ2=Γv[2],
            m3=Mv[3]-mth, Γ3=Γv[3],
            m4=Mv[4]-mth, Γ4=Γv[4],
            f2=5.0, f3=0.2, f4=0.1, fb=0.1,
            α = αv, β = βv
    return pdf_direct

I need modular implementation for easy calculation of integrals over the signal and background.

function construct_pdf_modular_integral_conv()
    pdf1 = pdf(@. (e;p)->abs2(BWmth(e, p.m1, p.Γ1));  p0 = (m1=Mv[1]-mth, Γ1=Γv[1]), lims = fitlims)
    pdf2 = pdf(@. (e;p)->abs2(BWmth(e, p.m2, p.Γ2));  p0 = (m2=Mv[2]-mth, Γ2=Γv[2]), lims = fitlims)
    pdf3 = pdf(@. (e;p)->abs2(BWmth(e, p.m3, p.Γ3));  p0 = (m3=Mv[3]-mth, Γ3=Γv[3]), lims = fitlims)
    pdf4 = pdf(@. (e;p)->abs2(BWmth(e, p.m4, p.Γ4));  p0 = (m4=Mv[4]-mth, Γ4=Γv[4]), lims = fitlims)
    pdfb = pdf(@. (e;p)->bgd(e, αv, βv) * pqmth(e);  p0 = NamedTuple(), lims = fitlims)
    pdfS = [
        conv_with_gauss(pdf1, σev[1]),
        conv_with_gauss(pdf2, σev[2]),
        conv_with_gauss(pdf3, σev[3]),
        conv_with_gauss(pdf4, σev[4]),
    phsp = pdf(@. (e;p)->pqmth(e);  p0 = NamedTuple(), lims = fitlims)
    [d *= phsp for d in pdfS]
    pdfS[2] *= (f2=5.0,)
    pdfS[3] *= (f3=0.2,)
    pdfS[4] *= (f4=0.1,)
    pdfS[5] *= (fb=0.1,)
    pdf_modular = sum(pdfS);
    return pdf_modular

Here is the benchmarks: the modular implementation is x10 slower :frowning:

using BenchmarkTools
const xv = range(fitlims..., length=10)
pdf_direct = construct_pdf_directly_integral_conv()
@btime pdf_direct(0.1) # 19 ms
@btime pdf_direct(xv)  # 20 ms
pdf_modular = construct_pdf_modular_integral_conv()
@btime pdf_modular(0.1) # 221 ms
@btime pdf_modular(xv) # 212 ms 

Plotting of the picture above (for completeness):

using Plots
using HEPRecipes #
    Nev = 1000
    data = generate(Nev, pdf_direct; Nbins=500)
    bins = range(fitlims..., length=100)
    N = Nev / (length(bins)-1) * (fitlims[2]-fitlims[1])
    fine = range(fitlims..., length=300)
    plot(fine, N .* pdf_direct(fine), l=(2, :red), lab="model")
    wh = weightedHistogram(data, bins=bins)
    plot!(wh, m=(3,:black), l=(2, :red), lab="data")
1 Like

Foreseeing a question on the usage of Distributions.jl: I would love to, I just could not find an easy way to implement customary density functions (always codding 10 interface function looked tedious)

Not an answer to your real question, but this claim does not seem correct to me: you can easily create a new distribution without defining those functions if you never try to call those functions. I think the assertion that you need to implement them reflects the bar expected of a distribution that would be contributed to the library itself, not what you can do with your own local code.

julia> import Distributions: ContinuousUnivariateDistribution, pdf

julia> struct MyD <: ContinuousUnivariateDistribution

julia> pdf(d::MyD, x::Real) = x < 0 ? 0.0 : x > 1 ? 0.0 : 1.0
pdf (generic function with 81 methods)

julia> pdf(d, [1.1, 2.2])
2-element Array{Float64,1}:

great! What is about parameters?
can I make the pdf dependent on parameters to be optimized by the llh fit?

I am missing something in your example.

using Distributions
import Distributions: ContinuousUnivariateDistribution, pdf
struct MyD <: ContinuousUnivariateDistribution
pdf(d::MyD, x::Real) = x < 0 ? 0.0 : x > 2 ? 0.0 : 1.0
d = MyD()
pdf(d, [1.1, 2.2])

It tells

MethodError: no method matching pdf(::MyD, ::Array{Float64,1})

EDIT: resolved

Mostly unrelated question, what kind of spectra are these? Almost looks like a small chunk from some XPS spectra. Some appear more lorentzian than guassian. Have you considered using voigt profiles?

I work on high energy physics. These are the peaks of resonances (particles) in the energy spectrum of invariant masses of the decay products. E.g. see Fig.5 in

The peaks are described by the Breit-Wigner functions, often convoluted with gaussian resolution.
What makes it non-standard always is the phase space factor that multiplied the decay probability function.

1 Like

Sorry, that works! (originally had a conflict in the session with existing pdf)

I also find how you deal with parameters: they are in the structure. What I would be still missing is automatic computation of normalization integral.

  • Are there plans to define a pdf type that would automatically normalize itself using QuadGK for example?

I can’t speak for the maintainers, but it seems likely you’d be better off doing that computation once in the constructor and caching it rather than recomputing it per-evaluation of the PDF. You see a pattern along those lines in the existing Truncated type:

1 Like

Re-evaluation of normalization looks unavoidable when parameters change. To evaluate PDF on a set with the normalization computed once I just run myPDF(dataset), i.e. taking the vector as an argument.

Concerning Truncated, is it just to limit pdf to the range (l,r)?
Can you give a reference to an example in Distribution.jl that implements an example close to my problem? I would be very curious to code my pdf one day and check the speed.