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

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 :

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

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.

1 Like

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 !

1 Like

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.