Nested lambda functions overhead: ideas for algebra of functions

General Usage
, ,
1

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.

image
image709×470 27.3 KB

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

using AlgebraPDF # https://github.com/mmikhasenko/AlgebraPDF.j

# 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)
    end
    #
    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
end

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]),
        pdfb]
    #
    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
end

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 # https://github.com/mmikhasenko/HEPRecipes.jl
let 
    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")
end
1 Like
2

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)

3

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
       end

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}:
 0.0
 0.0
2 Likes
4

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

5

I am missing something in your example.

using Distributions
import Distributions: ContinuousUnivariateDistribution, pdf
struct MyD <: ContinuousUnivariateDistribution
end
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

6

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?

7

Hey,
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 https://arxiv.org/pdf/1904.03947.pdf

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
8

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?
9

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: https://github.com/JuliaStats/Distributions.jl/blob/e127e9d5ec0bdb83152e1af9f0a26a5a18bbbefe/src/truncate.jl#L53

1 Like
10

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.