# Nested lambda functions overhead: ideas for algebra of functions

**URL:** <https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335>\
**Category:** General Usage\
**Tags:** question, statistics, optimization\
**Created:** [September 26, 2020, 7:29pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335 "2020-09-26T19:29:17Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [September 26, 2020, 7:29pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/1 "2020-09-26T19:29:18Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/c/8/c843cea56148d32a937bf561f4e49722ed39f2cd.png)

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

```julia
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.

```julia
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 ☹

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [September 26, 2020, 7:33pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/2 "2020-09-26T19:33:43Z")

</div>

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](https://juliastats.org/Distributions.jl/stable/extends/#Univariate-Distribution-1) looked tedious)

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [September 26, 2020, 7:58pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/3 "2020-09-26T19:58:00Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [September 26, 2020, 9:37pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/4 "2020-09-26T21:37:16Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [September 26, 2020, 9:46pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/5 "2020-09-26T21:46:08Z")

</div>

I am missing something in your example.

```julia
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

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

```

EDIT: resolved

---

<div class="post-metadata">

**Author:** ![anon92994695](https://avatars.discourse-cdn.com/v4/letter/a/ce7236/32.png) [@anon92994695](https://discourse.julialang.org/u/anon92994695)\
**Post date:** [September 26, 2020, 10:10pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/6 "2020-09-26T22:10:51Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [September 26, 2020, 10:21pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/7 "2020-09-26T22:21:21Z")

</div>

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](https://arxiv.org/pdf/1904.03947.pdf)

The peaks are described by the [Breit-Wigner functions](https://en.wikipedia.org/wiki/Relativistic_Breit%E2%80%93Wigner_distribution), often convoluted with gaussian resolution.  
What makes it non-standard always is the phase space factor that multiplied the decay probability function.

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [September 26, 2020, 10:30pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/8 "2020-09-26T22:30:08Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [September 27, 2020, 11:27am UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/9 "2020-09-27T11:27:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![misha\_mikhasenko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/misha_mikhasenko/32/5060_2.png) [@misha\_mikhasenko](https://discourse.julialang.org/u/misha_mikhasenko)\
**Post date:** [September 27, 2020, 2:49pm UTC](https://discourse.julialang.org/t/nested-lambda-functions-overhead-ideas-for-algebra-of-functions/47335/10 "2020-09-27T14:49:35Z")

</div>

> [@johnmyleswhite](#):
>
> 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

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.