# A mutable struct that does not get mutated ? or does it?

**URL:** <https://discourse.julialang.org/t/a-mutable-struct-that-does-not-get-mutated-or-does-it/50537>\
**Category:** Statistics\
**Created:** [November 21, 2020, 2:49pm UTC](https://discourse.julialang.org/t/a-mutable-struct-that-does-not-get-mutated-or-does-it/50537 "2020-11-21T14:49:11Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 21, 2020, 2:49pm UTC](https://discourse.julialang.org/t/a-mutable-struct-that-does-not-get-mutated-or-does-it/50537/1 "2020-11-21T14:49:12Z")

</div>

Hi,

I’m trying to code an extension do Distributions.jl for a certain distribution (Some gamma convolutions), and implementing the Moschopoulos algorithm.

I have an issue: To compute the pdf, i need to cache certain values, that i compute on demand. When running `pdf.(dist,1:10)`, the caching works correctly and each new call to `pdf` has values cached from the previous one.

On the other hand, when looking at the object `dist` **after** calling the function, cached values are gone.

Could someone explain to me what is going on ? Here is the code :

```julia

using Distributions
# Deal wich moshopoulos parameters :
mutable struct MoschopoulosParameters{T}
        θ₁::T
        C::T
        to_power::Array{T,1}
        γ::Array{T,1}
        ρ::T
        δ::Array{T,1}
end
function MoschopoulosParameters(α,θ)
    T = Base.promote_eltype(α,θ,[1.0])
    α = T.(α)
    θ = T.(θ)
    δ = T.([1])
    θ₁ = Base.minimum(θ)
    C = prod((θ ./ θ₁) .^ α)
    ρ = sum(α)
    to_power = (-θ₁ ./ θ) .+1
    γ = [sum(α .* to_power)] # gamma1
    return MoschopoulosParameters(θ₁,C,to_power,γ,ρ,δ)
end
eltype(P::MoschopoulosParameters) = typeof(P.θ₁)

struct UnivariateGammaConvolution{T<:Real} <: Distributions.ContinuousUnivariateDistribution where T
    α::Array{T}
    θ::Array{T}
    P::MoschopoulosParameters{T}
end

# Constructors :
UnivariateGammaConvolution(α::Real,θ::Real) = Distributions.Gamma(α,θ)
function UnivariateGammaConvolution(α::Array{T1,1},θ::Array{T2,1}) where {T1 <: Real, T2 <: Real}

    # Start by promoving alpha and theta to the same type :
    T = Base.promote_eltype(α,θ,[1.0]) # At least float.
    α = T.(α)
    θ = T.(θ)

    # Remove alphas and theta that are non-positives
    if any(α .== T(0))
        are_pos = α .* θ .!= T(0)
        α = α[are_pos]
        θ = θ[are_pos]
    end

    # regroup theta that are the sames :
    n = length(θ)
    for i in 1:(n-1)
        if i >= length(θ)
            break
        end
        for j in (i+1):length(θ)
            if θ[i] == θ[j]
                α[i] += α[j]
                deleteat!(α,j)
                deleteat!(θ,j)
                j = j-1
            end
        end
    end

    if length(α) == 1
        return Distributions.Gamma(α,θ)
    end
    P = MoschopoulosParameters(α,θ)
    tt = eltype(P)
    α = tt.(α)
    θ = tt.(θ)
    return UnivariateGammaConvolution(α,θ,MoschopoulosParameters(α,θ))
end

#### Support
Distributions.@distr_support UnivariateGammaConvolution 0.0 Inf

#### Conversions
eltype(d::UnivariateGammaConvolution) = typeof(d.α[1])
Base.convert(::Type{UnivariateGammaConvolution{T}}, d::UnivariateGammaConvolution{S}) where {T <: Real, S <: Real} = UnivariateGammaConvolution(T.(d.α), T.(d.θ))

#### parameters
shapes(d::UnivariateGammaConvolution) = d.α
scales(d::UnivariateGammaConvolution) = d.θ
rates(d::UnivariateGammaConvolution) = 1 / d.θ

#### Statistics
mean(d::UnivariateGammaConvolution) = sum(d.α .* d.θ)
var(d::UnivariateGammaConvolution) = sum(d.α .* d.θ .^ 2)

#### Sampling
function Base.rand(rng::Distributions.AbstractRNG,d::UnivariateGammaConvolution)
    sum(rand.(rng,Distributions.Gamma.(d.α,d.θ)))
end

#### Characteristic functions: pdf, cdf, mgf, cf

# Moshopoulos algorithm for pdf and cdf.
function MoschopoulosAlgo!(d::UnivariateGammaConvolution,x::Real, which, tolerence)

    @assert(which in ["pdf","cdf"], "which should be etiher pdf or cdf")
    print(length(d.P.γ),"at beginning\n")
    T = eltype(d)
    x = T(x)
    if x < T(0)
        return 0
    end
    k = 0
    out = T(0)
    while true
        if length(d.P.δ) < k+1
            # then compute the new deltas:
            n = length(d.P.δ)
            for k in n:(k+1)
                pushfirst!(d.P.γ, sum(d.α .* d.P.to_power .^ (k+1)))
                push!(d.P.δ,sum(d.P.γ[2:end] .* d.P.δ)/(k))
            end
        end
        dist = Distributions.Gamma(d.P.ρ + k,d.P.θ₁)
        if which == "pdf"
            step = d.P.δ[k+1] * Distributions.pdf(dist,x)
        elseif which == "cdf"
            step = d.P.δ[k+1] * Distributions.cdf(dist,x)
        end

        @assert(isfinite(step),"inf or nan append, the algorithm did not converge")
        out += step
        if isapprox(step,T(0),atol=tolerence,rtol=T(0))
            break
        end
        k = k+1
    end
    out *= d.P.C
    print(length(d.P.γ),"at end\n")
    return out
end

pdf(d::UnivariateGammaConvolution,x) = MoschopoulosAlgo!(d,x,"pdf", eps(eltype(d)))
cdf(d::UnivariateGammaConvolution, x) = MoschopoulosAlgo!(d,x,"cdf", eps(eltype(d)))
logpdf(d::UnivariateGammaConvolution, x) = log(pdf(d,x))
mgf(d::UnivariateGammaConvolution, t::Real) = prod((1 - t .* d.θ) .^ (-d.α))
cf(d::UnivariateGammaConvolution, t::Real) = prod((1 - (im * t) .* d.θ) .^ (-d.α))

# Let's now test that :
dist = UnivariateGammaConvolution([1,0.5],[4,2])
sample = zeros(Float64,100)
import Random
Random.rand!(dist,sample)
x = pdf.(dist,1:10)
print(length(dist.P.γ)) # Returns 1 !!

```

Edit : Looks like `pdf(dist,x)` caches the values correctly, while `pdf.(dist,1:100)` does not. I do not get why.

---

<div class="post-metadata">

**Author:** ![jebej](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jebej/32/1784_2.png) [@jebej](https://discourse.julialang.org/u/jebej)\
**Post date:** [November 21, 2020, 5:03pm UTC](https://discourse.julialang.org/t/a-mutable-struct-that-does-not-get-mutated-or-does-it/50537/2 "2020-11-21T17:03:46Z")

</div>

It probably has to do with the fact that you are broadcasting over `dist`.

By default, a dot call will try to iterate over all its arguments. Try wrapping `dist` in a `Ref` or a 1-element tuple: `pdf.(Ref(dist),1:100)` to prevent this from happening.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 21, 2020, 5:22pm UTC](https://discourse.julialang.org/t/a-mutable-struct-that-does-not-get-mutated-or-does-it/50537/3 "2020-11-21T17:22:30Z")

</div>

> [@jebej](#):
>
> t probably has to do with the fact that you are broadcasting over `dist` .
> 
> By default, a dot call will try to iterate over all its arguments. Try wrapping `dist` in a `Ref` or a 1-element tuple: `pdf.(Ref(dist),1:100)` to prevent this from happening.

It did not work, sadly. I dont get why a broadcasted call will not mutate it’s objects as a normal call, though.

Edit: The `Ref()` did not work, but the one ellement tuple did !

---

<div class="post-metadata">

**Author:** ![jebej](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jebej/32/1784_2.png) [@jebej](https://discourse.julialang.org/u/jebej)\
**Post date:** [November 21, 2020, 5:38pm UTC](https://discourse.julialang.org/t/a-mutable-struct-that-does-not-get-mutated-or-does-it/50537/4 "2020-11-21T17:38:21Z")

</div>

> [@lrnv](#):
>
> It did not work, sadly. I dont get why a broadcasted call will not mutate it’s objects as a normal call, though.

I’m not sure exactly why, but since the broadcast call didn’t just entirely error out, I assume that your type `UnivariateGammaConvolution` must have a broadcasting rule defined implicitly for it by the `Distributions.jl` package because it is a subtype of `Distributions.ContinuousUnivariateDistribution`.

It might be worth opening an issue there to see if that is an intended behavior or not. I am also surprised that the `Ref` trick didn’t work, maybe mention that too.