From e0bbf7868b722c282c61c45d576a1eae229afced Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Wed, 4 Sep 2024 00:41:34 -0400
Subject: [PATCH 1/6] Implement ExpGammaIPSampler

Co-Authored-By: David Widmann <devmotion@users.noreply.github.com>
---
 src/samplers.jl          |  1 +
 src/samplers/expgamma.jl | 22 ++++++++++++++++++++++
 2 files changed, 23 insertions(+)
 create mode 100644 src/samplers/expgamma.jl

diff --git a/src/samplers.jl b/src/samplers.jl
index fa667fff2b..3b2c9cad87 100644
--- a/src/samplers.jl
+++ b/src/samplers.jl
@@ -20,6 +20,7 @@ for fname in ["aliastable.jl",
               "poisson.jl",
               "exponential.jl",
               "gamma.jl",
+              "expgamma.jl",
               "multinomial.jl",
               "vonmises.jl",
               "vonmisesfisher.jl",
diff --git a/src/samplers/expgamma.jl b/src/samplers/expgamma.jl
new file mode 100644
index 0000000000..8d6ca3bbd8
--- /dev/null
+++ b/src/samplers/expgamma.jl
@@ -0,0 +1,22 @@
+# These are used to bypass subnormals when sampling from
+
+# Inverse Power sampler
+# uses the x*u^(1/a) trick from Marsaglia and Tsang (2000) for when shape < 1
+struct ExpGammaIPSampler{S<:Sampleable{Univariate,Continuous},T<:Real} <: Sampleable{Univariate,Continuous}
+    s::S #sampler for Gamma(1+shape,scale)
+    nia::T #-1/scale
+end
+
+ExpGammaIPSampler(d::Gamma) = ExpGammaIPSampler(d, GammaMTSampler)
+function ExpGammaIPSampler(d::Gamma, ::Type{S}) where {S<:Sampleable}
+    shape_d = shape(d)
+    sampler = S(Gamma{partype(d)}(1 + shape_d, scale(d)))
+    return GammaIPSampler(sampler, -inv(shape_d))
+end
+
+function rand(rng::AbstractRNG, s::ExpGammaIPSampler)
+    x = log(rand(rng, s.s))
+    e = randexp(rng)
+    return muladd(s.nia, e, x)
+end
+

From 036e38c9a9a043c7e9c5913584caf4a9e06753be Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Wed, 4 Sep 2024 00:42:29 -0400
Subject: [PATCH 2/6] Implement ExpGammaSSSampler

Co-Authored-By: chelate <42802644+chelate@users.noreply.github.com>
---
 src/samplers/expgamma.jl | 36 ++++++++++++++++++++++++++++++++++++
 1 file changed, 36 insertions(+)

diff --git a/src/samplers/expgamma.jl b/src/samplers/expgamma.jl
index 8d6ca3bbd8..5d5b8d1637 100644
--- a/src/samplers/expgamma.jl
+++ b/src/samplers/expgamma.jl
@@ -20,3 +20,39 @@ function rand(rng::AbstractRNG, s::ExpGammaIPSampler)
     return muladd(s.nia, e, x)
 end
 
+
+# Small Shape sampler
+# From Liu, C., Martin, R., and Syring, N. (2015) for when shape < 0.3
+struct ExpGammaSSSampler{T<:Real} <: Sampleable{Univariate,Continuous}
+    α::T
+    θ::T
+    λ::T
+    ω::T
+    ωω::T
+    η::T
+end
+
+function ExpGammaSSSampler(d::Gamma)
+    α = shape(d)
+    ω = α / MathConstants.e / (1 - α)
+    return ExpGammaSSSampler(promote(
+        α,
+        scale(d),
+        inv(α) - 1,
+        ω,
+        inv(ω + 1)
+    ))
+end
+
+function rand(rng::AbstractRNG, s::ExpGammaSSSampler)
+    while true
+        U = rand(rng)
+        z = (U <= s.ωω) ? -log(U / s.ωω) : log(rand(rng)) / s.λ
+        h = exp(-z - exp(-z / α))
+        η = z >= 0 ? exp(-z) : s.ω * s.λ * exp(s.λ * z)
+        if h / η > rand(rng)
+            return z / α
+        end
+    end
+end
+

From a078362c9692db3cd8156c96d056496355282165 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Wed, 4 Sep 2024 02:44:14 -0400
Subject: [PATCH 3/6] Apply #1885 type change to ExpGammaIPSampler

---
 src/samplers/expgamma.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/samplers/expgamma.jl b/src/samplers/expgamma.jl
index 5d5b8d1637..aae725c634 100644
--- a/src/samplers/expgamma.jl
+++ b/src/samplers/expgamma.jl
@@ -16,7 +16,7 @@ end
 
 function rand(rng::AbstractRNG, s::ExpGammaIPSampler)
     x = log(rand(rng, s.s))
-    e = randexp(rng)
+    e = randexp(rng, typeof(x))
     return muladd(s.nia, e, x)
 end
 

From cb431410d0e806e2d0f6fc64f7c23db4f5630e9c Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Wed, 4 Sep 2024 02:33:01 -0400
Subject: [PATCH 4/6] Implement improved Dirichlet rand

---
 src/samplers/expgamma.jl | 46 ++++++++++++++++++++++++++++++++--------
 1 file changed, 37 insertions(+), 9 deletions(-)

diff --git a/src/samplers/expgamma.jl b/src/samplers/expgamma.jl
index aae725c634..3fefbb3a41 100644
--- a/src/samplers/expgamma.jl
+++ b/src/samplers/expgamma.jl
@@ -29,7 +29,6 @@ struct ExpGammaSSSampler{T<:Real} <: Sampleable{Univariate,Continuous}
     λ::T
     ω::T
     ωω::T
-    η::T
 end
 
 function ExpGammaSSSampler(d::Gamma)
@@ -41,18 +40,47 @@ function ExpGammaSSSampler(d::Gamma)
         inv(α) - 1,
         ω,
         inv(ω + 1)
-    ))
+    )...)
 end
 
-function rand(rng::AbstractRNG, s::ExpGammaSSSampler)
+function rand(rng::AbstractRNG, s::ExpGammaSSSampler{T})::Float64 where T
+    flT = float(T)
     while true
-        U = rand(rng)
-        z = (U <= s.ωω) ? -log(U / s.ωω) : log(rand(rng)) / s.λ
-        h = exp(-z - exp(-z / α))
-        η = z >= 0 ? exp(-z) : s.ω * s.λ * exp(s.λ * z)
-        if h / η > rand(rng)
-            return z / α
+        U = rand(rng, flT)
+        z = (U <= s.ωω) ? -log(U / s.ωω) : log(rand(rng, flT)) / s.λ
+        h = exp(-z - exp(-z / s.α))
+        η = z >= zero(T) ? exp(-z) : s.ω * s.λ * exp(s.λ * z)
+        if h / η > rand(rng, flT)
+            return s.θ - z / s.α
         end
     end
 end
 
+
+function _logsampler(d::Gamma)
+    if shape(d) < 0.3
+        return ExpGammaSSSampler(d)
+    else
+        return ExpGammaIPSampler(d)
+    end
+end
+
+function _logrand(rng::AbstractRNG, d::Gamma)
+    if shape(d) < 0.3
+        return rand(rng, ExpGammaSSSampler(d))
+    else
+        return rand(rng, ExpGammaIPSampler(d))
+    end
+end
+
+function _logrand!(rng::AbstractRNG, d::Gamma, A::AbstractArray{<:Real})
+    if shape(d) < 0.3
+        @inbounds for i in eachindex(A)
+            A[i] = rand(rng, ExpGammaSSSampler(d))
+        end
+    else
+        @inbounds for i in eachindex(A)
+            A[i] = rand(rng, ExpGammaIPSampler(d))
+        end
+    end
+end

From 799f90190525993ee39f45f4fe0dff10b638a37f Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Fri, 28 Mar 2025 17:23:43 -0400
Subject: [PATCH 5/6] Revert  _logsampler and _logrand additions

---
 src/samplers/expgamma.jl | 29 -----------------------------
 1 file changed, 29 deletions(-)

diff --git a/src/samplers/expgamma.jl b/src/samplers/expgamma.jl
index 3fefbb3a41..ed05dadece 100644
--- a/src/samplers/expgamma.jl
+++ b/src/samplers/expgamma.jl
@@ -55,32 +55,3 @@ function rand(rng::AbstractRNG, s::ExpGammaSSSampler{T})::Float64 where T
         end
     end
 end
-
-
-function _logsampler(d::Gamma)
-    if shape(d) < 0.3
-        return ExpGammaSSSampler(d)
-    else
-        return ExpGammaIPSampler(d)
-    end
-end
-
-function _logrand(rng::AbstractRNG, d::Gamma)
-    if shape(d) < 0.3
-        return rand(rng, ExpGammaSSSampler(d))
-    else
-        return rand(rng, ExpGammaIPSampler(d))
-    end
-end
-
-function _logrand!(rng::AbstractRNG, d::Gamma, A::AbstractArray{<:Real})
-    if shape(d) < 0.3
-        @inbounds for i in eachindex(A)
-            A[i] = rand(rng, ExpGammaSSSampler(d))
-        end
-    else
-        @inbounds for i in eachindex(A)
-            A[i] = rand(rng, ExpGammaIPSampler(d))
-        end
-    end
-end

From d67890488623f714790b85ee48a05b7cb1dd0a5e Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Fri, 28 Mar 2025 20:59:18 -0400
Subject: [PATCH 6/6] Fix wrong return in ExpGammaIPSampler

---
 src/samplers/expgamma.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/samplers/expgamma.jl b/src/samplers/expgamma.jl
index ed05dadece..875a98aa45 100644
--- a/src/samplers/expgamma.jl
+++ b/src/samplers/expgamma.jl
@@ -11,7 +11,7 @@ ExpGammaIPSampler(d::Gamma) = ExpGammaIPSampler(d, GammaMTSampler)
 function ExpGammaIPSampler(d::Gamma, ::Type{S}) where {S<:Sampleable}
     shape_d = shape(d)
     sampler = S(Gamma{partype(d)}(1 + shape_d, scale(d)))
-    return GammaIPSampler(sampler, -inv(shape_d))
+    return ExpGammaIPSampler(sampler, -inv(shape_d))
 end
 
 function rand(rng::AbstractRNG, s::ExpGammaIPSampler)