GammaIPSampler rand: avoid recursion

[nsajko]

The recursion here seems to prevent Julia from inferring that rand(::FixedRNG, ::Gamma) terminates, preventing constant folding.

Eliminating it is as simple as short-circuiting to GammaGDSampler. This is correct because the shape parameter here is one more than with the original Gamma distribution, so it must be greater than one, assuming it was positive to begin with.

Note that I'm not an expert and in fact I'm having much trouble with passing Probability&Statistics, so there might be something I'm not taking into account.

[nsajko]

Inspired by this Discourse discussion: https://discourse.julialang.org/t/constant-propagation-irrationals-doing-much-better-than-float64/87142/

Note that the pull request doesn't actually accomplish the goal of making rand constant-foldable, but it does manage to convince Julia that it's constant-foldable according to its effects-analysis:

Before my patches:

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.9.0-DEV.1323 (2022-09-10)
 _/ |\__'_|_|_|\__'_|  |  Commit d4f6fd24138 (2 days old master)
|__/                   |

julia> using Distributions, Random

julia> struct FixedRNG <: AbstractRNG end

julia> Base.rand(::FixedRNG) = one(Float64) / 2

julia> Random.randn(::FixedRNG) = zero(Float64)

julia> Random.randexp(::FixedRNG) = one(Float64)

julia> Base.rand(::FixedRNG, ::Type{T}) where {T<:Real} = one(T) / 2

julia> Random.randn(::FixedRNG, ::Type{T}) where {T<:Real} = zero(T)

julia> Random.randexp(::FixedRNG, ::Type{T}) where {T<:Real} = one(T)

julia> # We need concrete type parameters to avoid ambiguity for these cases
       for T in [Float16, Float32, Float64]
           @eval begin
               Base.rand(::FixedRNG, ::Type{$T}) = one($T) / 2
               Random.randn(::FixedRNG, ::Type{$T}) = zero($T)
               Random.randexp(::FixedRNG, ::Type{$T}) = one($T)
           end
       end

julia> rand(FixedRNG(), TDist(3.1))
0.0

julia> Base.infer_effects(rand, (FixedRNG, Gamma{Float64}))
(+c,+e,!n,!t,+s,+m)

julia> Base.infer_effects(rand, (FixedRNG, Chisq{Float64}))
(+c,+e,!n,!t,+s,+m)

julia> Base.infer_effects(rand, (FixedRNG, TDist{Float64}))
(+c,+e,!n,!t,+s,+m)

After (now we have "+t"):

julia> Base.infer_effects(rand, (FixedRNG, Gamma{Float64}))
(+c,+e,!n,+t,+s,+m)

julia> Base.infer_effects(rand, (FixedRNG, Chisq{Float64}))
(+c,+e,!n,+t,+s,+m)

julia> Base.infer_effects(rand, (FixedRNG, TDist{Float64}))
(+c,+e,!n,+t,+s,+m)

[codecov-commenter]

Codecov Report

Base: 85.94% // Head: 85.94% // No change to project coverage 👍

Coverage data is based on head (1a98fd7) compared to base (e67d6b3).
Patch coverage: 100.00% of modified lines in pull request are covered.

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1615   +/-   ##
=======================================
  Coverage   85.94%   85.94%           
=======================================
  Files         129      129           
  Lines        8080     8080           
=======================================
  Hits         6944     6944           
  Misses       1136     1136           

Impacted Files	Coverage Δ
src/samplers/gamma.jl	100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

[devmotion]

The recursion here seems to prevent Julia from inferring that rand(::FixedRNG, ::Gamma) terminates, preventing constant folding.

What's FixedRNG? Is there any problem with standard RNGs (from Random and/or RandomNumbers) that's fixed by this PR?

[nsajko]

What's FixedRNG?

It's just an example static RNG from the linked Discourse post, also defined in my previous message.

Is there any problem with standard RNGs (from Random and/or RandomNumbers) that's fixed by this PR?

I don't know.

When making the commit message, I didn't realize that I refer to FixedRNG, which is not defined in the commit message. If You want I can expand the commit message with information from my second message here to fix this.

[devmotion]

I think the motivation for this change is a bit unclear. However, there's actually a bug here but IMO it should be fixed differently (see my comment).

[devmotion]

Actually, I think this change is undesirable - but the current implementation is as well. The motivation for samplers is to pre-compute and pre-allocate everything required for sampling since it should pay off when you want to generate many samples. Hence by default, distributions only use dedicated samplers if multiple samples are requested (similar to the concept in Random).

So s.s should already be a sampler at this point here. However, it seems that's currently not the case - even though it seems that was the intention behing

Distributions.jl/src/samplers/gamma.jl

Lines 202 to 205 in e67d6b3

	function GammaIPSampler(d::Gamma,::Type{S}) where S<:Sampleable
	GammaIPSampler(Gamma(1.0 + shape(d), scale(d)), -1.0 / shape(d))
	end
	GammaIPSampler(d::Gamma) = GammaIPSampler(d,GammaMTSampler)

. The bug is that it does not use S though (and that it should be designed a bit differently probably).

[nsajko]

Yeah, I noticed that S and GammaMTSampler are unused. I think I found some issues on this repo from like eight years ago that mention that GammaMTSampler was "disabled".

[devmotion]

That comment does not refer to these lines, in any case. This bug was already present when GammaIPSampler was added: #313

[devmotion]

There was a bug in the implementaiton of GammaMTSampler, fixing it improves performance significantly: #1617 (comment)

[devmotion]

I close this PR in favour of #1617. Generally, we don't want to construct the sampler in the rand call but at construction time.