# Optimization question

**URL:** <https://discourse.julialang.org/t/optimization-question/120707>\
**Category:** New to Julia\
**Created:** [September 30, 2024, 4:22am UTC](https://discourse.julialang.org/t/optimization-question/120707 "2024-09-30T04:22:20Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [September 30, 2024, 4:22am UTC](https://discourse.julialang.org/t/optimization-question/120707/1 "2024-09-30T04:22:20Z")

</div>

I’m trying to optimize the code block below:

```julia
using Distributions
order = 2
σ_Priors = Array{Distribution,2}(undef,order,order)
ϵ_Priors = Array{Distribution,2}(undef,order,order)
σ_Priors[:,:] .= Gamma(10,0.1)
ϵ_Priors[:,:] .= Gamma(10,0.1)
std_Prior = Gamma(10,0.3)

function totalEnergy(crystal::Crystal,LJ::LJ)
    totalEnergy = 0
    for i in eachindex(LJ.ϵ)
        totalEnergy += -LJ.ϵ[i] * LJ.σ[i] ^6 * crystal.r6[i] + LJ.ϵ[i] * LJ.σ[i] .^12 * crystal.r12[i]
    end
# @time "total Energy" totalEnergy = sum(-LJ.ϵ .* LJ.σ .^6 .* crystal.r6 .+ LJ.ϵ .* LJ.σ .^12 .* crystal.r12 ) 
end

function logNormal(data::MatSim.DataSet,LJ::MatSim.LJ,σ::Float64)
    thesum = 0
    for i = 1:data.nData
        thesum += (data.crystals[i].energyFP -MatSim.totalEnergy(data.crystals[i],LJ))^2
    end
    thesum *= -data.nData/2 *log(σ^2) - 1/(2 * σ^2)
# @time "thesum" thesum = -data.nData/2 *log(σ^2) - 1/(2 * σ^2) * sum([(i.energyFP - totalEnergy(i,LJ))^2 for i in data.crystals])
    return thesum

end

function logPost(data,model,σ)
    # return 3.0
    thesum = logNormal(data,model,σ) + logpdf(std_Prior,σ)
    for i in eachindex(model.σ)
        thesum += logpdf(σ_Priors[i],model.σ[i])
        thesum += logpdf(ϵ_Priors[i],model.ϵ[i])
    end
    return thesum
end

@time logPost(trainingSet,LJ,0.3)

```

I recognize that I’m not giving you all of the code but I’m hoping it’s enough to go on. Currently a single call to logPost requires 34 allocations. Doesn’t seem too bad except that this function is being called hundreds of thousands of times. I’m convinced that this function is the performance bottleneck. Can anyone see any glaring problems here? Perhaps using globals for `σ_Priors` instead of passing them into the function?

Thanks in advance.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [September 30, 2024, 4:30am UTC](https://discourse.julialang.org/t/optimization-question/120707/2 "2024-09-30T04:30:08Z")

</div>

> [@lancejnelson](#):
>
> Perhaps using globals for `σ_Priors` instead of passing them into the function?

Yea, mark that const is a good start

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [September 30, 2024, 4:31am UTC](https://discourse.julialang.org/t/optimization-question/120707/3 "2024-09-30T04:31:24Z")

</div>

Or if I just pass them into the function that would work too right?

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [September 30, 2024, 4:52am UTC](https://discourse.julialang.org/t/optimization-question/120707/4 "2024-09-30T04:52:54Z")

</div>

`Distribution` is an abstract type, hence

```julia
logpdf(σ_Priors[i],model.σ[i])

```

is a dynamic (and thus expensive) dispatch. If you only have gamma distributions, do this instead

```julia
σ_Priors = fill(Gamma(10,0.1), order, order)

```

(same for `ϵ_Priors`)

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [September 30, 2024, 11:52am UTC](https://discourse.julialang.org/t/optimization-question/120707/5 "2024-09-30T11:52:24Z")

</div>

For example, following the previous tips of the other thread:

With a parametric struct:

```julia-repl
julia> struct ModelA{T}
           σ::Vector{Float64}
           ϵ::Vector{Float64}
           σ_Priors::Matrix{T}
           ϵ_Priors::Matrix{T}
       end

julia> function logPost(model)
           # return 3.0
           thesum = 0.0 # simplified version
           for i in eachindex(model.σ)
               thesum += logpdf(model.σ_Priors[i],model.σ[i])
               thesum += logpdf(model.ϵ_Priors[i],model.ϵ[i])
           end
           return thesum
       end
logPost (generic function with 2 methods)

julia> m = ModelA(rand(100), rand(100), fill(Gamma(10,0.1), 10, 10), fill(Gamma(10,0.1), 10, 10));

julia> typeof(m)
ModelA{Gamma{Float64}}

julia> @btime logPost($m)
  3.100 μs (0 allocations: 0 bytes)
-675.8052525561492

```

With the abstract field, instead:

```julia-repl
julia> struct ModelB
           σ::Vector{Float64}
           ϵ::Vector{Float64}
           σ_Priors::Matrix{Distribution}
           ϵ_Priors::Matrix{Distribution}
       end

julia> mB = ModelB(m.σ, m.ϵ, m.σ_Priors, m.ϵ_Priors);

julia> typeof(mB)
ModelB

julia> @btime logPost($mB)
  6.455 μs (600 allocations: 9.38 KiB)
-675.8052525561492

```

FWIW, I think this blog post might be helpful: [New blog post about Julia parametric types and constructors](https://discourse.julialang.org/t/new-blog-post-about-julia-parametric-types-and-constructors/120717)

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [September 30, 2024, 12:42pm UTC](https://discourse.julialang.org/t/optimization-question/120707/6 "2024-09-30T12:42:58Z")

</div>

> [@lancejnelson](#):
>
> `LJ::LJ`

I think perhaps this is not causing any problems right here, but I don’t think it’s good style to shadow the typename like this.

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [September 30, 2024, 2:54pm UTC](https://discourse.julialang.org/t/optimization-question/120707/7 "2024-09-30T14:54:57Z")

</div>

> [@lmiq](#):
>
> ```julia
> struct ModelA{T}
> σ::Vector{Float64}
> ϵ::Vector{Float64}
> σ_Priors::Matrix{T}
> ϵ_Priors::Matrix{T}
> end
> 
> julia> function logPost(model)
> # return 3.0
> thesum = 0.0 # simplified version
> for i in eachindex(model.σ)
> thesum += logpdf(model.σ_Priors[i],model.σ[i])
> thesum += logpdf(model.ϵ_Priors[i],model.ϵ[i])
> end
> return thesum
> end
> logPost (generic function with 2 methods)
> 
> julia> m = ModelA(rand(100), rand(100), fill(Gamma(10,0.1), 10, 10), fill(Gamma(10,0.1), 10, 10));
> 
> julia> typeof(m)
> 
> isconcretetype(m) #false
> 
> @time logPost(m)
> 
> ```

I ran this exact code and got

```julia

0.000010 seconds (1 allocation: 16 bytes)

```

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [September 30, 2024, 3:13pm UTC](https://discourse.julialang.org/t/optimization-question/120707/8 "2024-09-30T15:13:44Z")

</div>

Note that it’s `@btime` not `@time` (from BenchmarkTools.jl).

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [September 30, 2024, 4:24pm UTC](https://discourse.julialang.org/t/optimization-question/120707/9 "2024-09-30T16:24:16Z")

</div>

I see what you’re saying. This declaration:

```julia
σ_Priors = Array{Distribution,2}(undef,order,order)

```

is problematic. It turns out that not all of the distributions are identical. They may all be Gammas, but with different parameters. How could I avoid this problem in this case?

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [September 30, 2024, 4:29pm UTC](https://discourse.julialang.org/t/optimization-question/120707/10 "2024-09-30T16:29:42Z")

</div>

If you know all possible Distributions then you can just pack them into a single type by using DynamicalSumTypes.jl

---

<div class="post-metadata">

**Author:** ![DanielVandH](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielvandh/32/31134_2.png) [@DanielVandH](https://discourse.julialang.org/u/DanielVandH)\
**Post date:** [September 30, 2024, 4:34pm UTC](https://discourse.julialang.org/t/optimization-question/120707/11 "2024-09-30T16:34:09Z")

</div>

Does the gamma distribution use different types for different parameters? If they are all the same type, you could just use `typeof(Gamma(10, 0.1))` for example. If you have a few different distribution types, wrapping them in a `Union` type could be fine. Alternatively, use LightSumTypes.jl as mentioned above

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [September 30, 2024, 11:00pm UTC](https://discourse.julialang.org/t/optimization-question/120707/12 "2024-09-30T23:00:01Z")

</div>

I’m showing that the function `logPost` is type unstable. Since `model.σ_Priors[i]` is a distribution the type of `logpdf` is `Any`.

Here’s the snippet from `@code_warntype`:

```julia
%138 = thesum::Any
%139 = Base.getproperty(metrop, :ϵ_Priors)::Matrix{Distribution}
%140 = i@_29::Int64
%141 = Base.getindex(%139, %140, j)::Distribution
%142 = Base.getproperty(model, :ϵ)::LinearAlgebra.UpperTriangular{Float64, Matrix{Float64}}
%143 = i@_29::Int64
%144 = Base.getindex(%142, %143, j)::Float64
%145 = MatSim.logpdf(%141, %144)::Any
%146 = (%138 + %145)::Any

```

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [September 30, 2024, 11:05pm UTC](https://discourse.julialang.org/t/optimization-question/120707/13 "2024-09-30T23:05:38Z")

</div>

My `model.σ` is an upper triangular matrix:

```julia

σ:: UpperTriangular{Float64, Matrix{Float64}}
ϵ:: UpperTriangular{Float64, Matrix{Float64}}

```

Could that be part of the problem?

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [September 30, 2024, 11:40pm UTC](https://discourse.julialang.org/t/optimization-question/120707/14 "2024-09-30T23:40:13Z")

</div>

I don’t think so. Where that `UpperTriangular` type is coming from?

The main reason is probably the abstract `Distribution` type. If you only need distributions of type `Gamma` as elements of of the priors (or other of the same type with varying parameters), you should parameterize the type.

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [October 1, 2024, 3:07pm UTC](https://discourse.julialang.org/t/optimization-question/120707/15 "2024-10-01T15:07:36Z")

</div>

Let’s assume the matrix was going to hold some Gammas and some Uniforms, would the following type definition work?

```julia
struct ModelA{T<:Union{Gamma,Uniform}}
           σ::Vector{Float64}
           ϵ::Vector{Float64}
           σ_Priors::Matrix{T}
           ϵ_Priors::Matrix{T}
       end

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [October 1, 2024, 3:53pm UTC](https://discourse.julialang.org/t/optimization-question/120707/17 "2024-10-01T15:53:01Z")

</div>

You should possibly parameterize `Gamma{Float64}` and `Uniform{Float64}`, so the compiler can predict the type of the sampled data.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [October 1, 2024, 4:04pm UTC](https://discourse.julialang.org/t/optimization-question/120707/18 "2024-10-01T16:04:27Z")

</div>

Actually that can be specified on the structure creation, not on the type. With the same generic type, you get, for example:

```julia-repl
julia> using Distributions, BenchmarkTools

julia> @kwdef struct ModelA{T}
           σ::Vector{Float64}
           ϵ::Vector{Float64}
           σ_Priors::Matrix{T}
           ϵ_Priors::Matrix{T}
       end

julia> function logPost(model)
           # return 3.0
           thesum = 0.0 # simplified version
           for i in eachindex(model.σ)
               thesum += logpdf(model.σ_Priors[i],model.σ[i])
               thesum += logpdf(model.ϵ_Priors[i],model.ϵ[i])
           end
           return thesum
       end
logPost (generic function with 1 method)

julia> T = Union{Gamma{Float64},Uniform{Float64}}

julia> m = ModelA(
           σ=rand(100), 
           ϵ=rand(100), 
           σ_Priors=T[rand(Bool) ? Gamma(10,0.1) : Uniform(0.1,10) for _ in 1:10, _ in 1:10],
           ϵ_Priors=T[rand(Bool) ? Gamma(10,0.1) : Uniform(0.1,10) for _ in 1:10, _ in 1:10],
       );

julia> typeof(m)
ModelA{Union{Gamma{Float64}, Uniform{Float64}}}

julia> @btime logPost($m)
  1.911 μs (0 allocations: 0 bytes)
-Inf

```

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [October 1, 2024, 4:23pm UTC](https://discourse.julialang.org/t/optimization-question/120707/19 "2024-10-01T16:23:17Z")

</div>

Interesting. Thank you.

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [October 1, 2024, 4:30pm UTC](https://discourse.julialang.org/t/optimization-question/120707/20 "2024-10-01T16:30:48Z")

</div>

Can I ask one more kinda related question and then I’ll be quiet.

1. One “type” of array is Array{Float64,2}, meaning a 2 dimensional array filled with floats. Are users allowed to define types with integers in their type definition? If so, what does that look like?

```julia
struct myType{T,N}
    x::T
    ...?

```

1. For best performance, user-defined types should be concrete when instantiated right? (Not sure I’m even asking this question right. I’m just wanting to build intuition for how to use types correctly and most efficiently. ). For example, I notice that `isconcretetype(Matrix{Distribution})` is true, but I wouldn’t have thought it was concrete because Distribution isn’t concrete. I guess I’m asking when is a composite type deemed concrete and is that something that should be strived for.

---

<div class="post-metadata">

**Author:** ![lancejnelson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lancejnelson/32/21898_2.png) [@lancejnelson](https://discourse.julialang.org/u/lancejnelson)\
**Post date:** [October 1, 2024, 6:19pm UTC](https://discourse.julialang.org/t/optimization-question/120707/21 "2024-10-01T18:19:39Z")

</div>

> [@lmiq](#):
>
> ```julia
> julia> @kwdef struct ModelA{T}
> σ::Vector{Float64}
> ϵ::Vector{Float64}
> σ_Priors::Matrix{T}
> ϵ_Priors::Matrix{T}
> end
> 
> ```

Why wouldn’t you also parameterize the Float64?

```julia
struct ModelA{F,T}
     σ::Vector{F}
     ϵ::Vector{F}
     σ_Priors::Matrix{T}
     ϵ_Priors::Matrix{T}
end

```

I understand that parameterizing with a concrete doesn’t create a family of types because you’re at the bottom of the type tree, but is there no performance gain by doing this?

[Next page](https://discourse.julialang.org/t/optimization-question/120707.md?page=2)