[Internals] Switch to TaylorDiff.jl for generator derivatives. (#38)

lrnv · web-flow · commit 29216bfe8bcc · 2025-07-24T16:17:49.000+02:00
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.6'
+          - 'lts'
           - '1'
         os:
           - ubuntu-latest
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Copulas"
 uuid = "ae264745-0b69-425e-9d9d-cf662c5eec93"
+version = "0.1.30"
 authors = ["Oskar Laverny"]
-version = "0.1.29"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -19,7 +19,6 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"
 
 [compat]
@@ -41,10 +40,9 @@ SpecialFunctions = "2"
 StableRNGs = "1"
 StatsBase = "0.33, 0.34"
 StatsFuns = "0.9, 1.3"
-TaylorSeries = "0.15, 0.16, 0.17, 0.18, 0.19, 0.20"
 Test = "1"
 TestItemRunner = "1"
-WilliamsonTransforms = "0.1"
+WilliamsonTransforms = "0.2"
 julia = "1"
 
 [extras]
diff --git a/src/ArchimedeanCopula.jl b/src/ArchimedeanCopula.jl
@@ -52,7 +52,7 @@ end
 ϕ⁻¹(C::ArchimedeanCopula{d,TG},t)  where {d,TG} = ϕ⁻¹(C.G,t)
 ϕ⁽¹⁾(C::ArchimedeanCopula{d,TG},t) where {d,TG} = ϕ⁽¹⁾(C.G,t)
 ϕ⁽ᵏ⁾(C::ArchimedeanCopula{d,TG},t) where {d,TG} = ϕ⁽ᵏ⁾(C.G, Val(d), t)
-williamson_dist(C::ArchimedeanCopula{d,TG}) where {d,TG} = williamson_dist(C.G,d)
+williamson_dist(C::ArchimedeanCopula{d,TG}) where {d,TG} = williamson_dist(C.G,Val(d))
 
 
 function _cdf(C::CT,u) where {CT<:ArchimedeanCopula} 
@@ -66,7 +66,7 @@ function Distributions._logpdf(C::ArchimedeanCopula{d,TG}, u) where {d,TG}
     if !all(0 .<= u .<= 1)
         return eltype(u)(-Inf)
     end
-    T = promote_type(Float64, eltype(u)) # the FLoat64 here should be eltype(C) when copulas wil be type agnostic... 
+    T = promote_type(Float64, eltype(u)) # the Float64 here should be eltype(C) when copulas will be type agnostic... 
     logdenom = sum_ϕ⁻¹u = zero(T)
     for us in u
         ϕ⁻¹u = ϕ⁻¹(C,us)
diff --git a/src/Copulas.jl b/src/Copulas.jl
@@ -8,7 +8,6 @@ module Copulas
     import Distributions
     import StatsBase
     import StatsFuns
-    import TaylorSeries
     import ForwardDiff
     import HCubature
     import MvNormalCDF
diff --git a/src/Generator.jl b/src/Generator.jl
@@ -36,11 +36,11 @@ max_monotony(G::Generator) = throw("This generator does not have a defined max m
 ϕ⁻¹( G::Generator, x) = Roots.find_zero(t -> ϕ(G,t) - x, (0.0, Inf))
 ϕ⁽¹⁾(G::Generator, t) = ForwardDiff.derivative(x -> ϕ(G,x), t)
 function ϕ⁽ᵏ⁾(G::Generator, ::Val{k}, t) where k
-    coef = WilliamsonTransforms.taylor(x -> ϕ(G, x), t, k)[end]
+    coef = WilliamsonTransforms.taylor(x -> ϕ(G, x), t, Val(k))[end]
     der = coef * factorial(k)
     return der
 end
-williamson_dist(G::Generator, d) = WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t),d)
+williamson_dist(G::Generator, ::Val{d}) where d = WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t), Val(d))
 
 # τ(G::Generator) = @error("This generator has no kendall tau implemented.")
 # ρ(G::Generator) = @error ("This generator has no Spearman rho implemented.")
diff --git a/src/Generator/UnivariateGenerator/AMHGenerator.jl b/src/Generator/UnivariateGenerator/AMHGenerator.jl
@@ -38,7 +38,7 @@ max_monotony(::AMHGenerator) = Inf
 ϕ⁻¹(G::AMHGenerator, t) = log(G.θ + (1-G.θ)/t)
 # ϕ⁽¹⁾(G::AMHGenerator, t) =  First derivative of ϕ
 # ϕ⁽ᵏ⁾(G::AMHGenerator, k, t) = kth derivative of ϕ
-williamson_dist(G::AMHGenerator, d) = G.θ >= 0 ? WilliamsonFromFrailty(1 + Distributions.Geometric(1-G.θ),d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t),d)
+williamson_dist(G::AMHGenerator, ::Val{d}) where d = G.θ >= 0 ? WilliamsonFromFrailty(1 + Distributions.Geometric(1-G.θ), Val(d)) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t), Val(d))
 
 function τ(G::AMHGenerator)
     θ = G.θ
diff --git a/src/Generator/UnivariateGenerator/ClaytonGenerator.jl b/src/Generator/UnivariateGenerator/ClaytonGenerator.jl
@@ -46,4 +46,4 @@ max_monotony(G::ClaytonGenerator) = G.θ >= 0 ? Inf : Int(floor(1 - 1/G.θ))
 # ϕ⁽ᵏ⁾(G::ClaytonGenerator, k, t) = kth derivative of ϕ
 τ(G::ClaytonGenerator) = ifelse(isfinite(G.θ), G.θ/(G.θ+2), 1)
 τ⁻¹(::Type{T},τ) where T<:ClaytonGenerator = ifelse(τ == 1,Inf,2τ/(1-τ))
-williamson_dist(G::ClaytonGenerator, d) = G.θ >= 0 ? WilliamsonFromFrailty(Distributions.Gamma(1/G.θ,G.θ),d) : ClaytonWilliamsonDistribution(G.θ,d)
+williamson_dist(G::ClaytonGenerator, ::Val{d}) where d = G.θ >= 0 ? WilliamsonFromFrailty(Distributions.Gamma(1/G.θ,G.θ), Val(d)) : ClaytonWilliamsonDistribution(G.θ,d)
diff --git a/src/Generator/UnivariateGenerator/FrankGenerator.jl b/src/Generator/UnivariateGenerator/FrankGenerator.jl
@@ -39,11 +39,11 @@ struct FrankGenerator{T} <: UnivariateGenerator
 end
 max_monotony(G::FrankGenerator) = G.θ < 0 ? 2 : Inf
 ϕ(  G::FrankGenerator, t) = G.θ > 0 ? -LogExpFunctions.log1mexp(LogExpFunctions.log1mexp(-G.θ)-t)/G.θ : -log1p(exp(-t) * expm1(-G.θ))/G.θ
-ϕ(  G::FrankGenerator, t::TaylorSeries.Taylor1) = G.θ > 0 ? -log(-expm1(LogExpFunctions.log1mexp(-G.θ)-t))/G.θ : -log1p(exp(-t) * expm1(-G.θ))/G.θ
+# ϕ(  G::FrankGenerator, t::TaylorSeries.Taylor1) = G.θ > 0 ? -log(-expm1(LogExpFunctions.log1mexp(-G.θ)-t))/G.θ : -log1p(exp(-t) * expm1(-G.θ))/G.θ
 ϕ⁻¹(G::FrankGenerator, t) = G.θ > 0 ? LogExpFunctions.log1mexp(-G.θ) - LogExpFunctions.log1mexp(-t*G.θ) : -log(expm1(-t*G.θ)/expm1(-G.θ))
 ϕ⁽¹⁾(G::FrankGenerator, t) = (one(t) - one(t) / (one(t) + exp(-t)*expm1(-G.θ))) / G.θ  # First derivative of ϕ
 # ϕ⁽ᵏ⁾(G::FrankGenerator, k, t) = kth derivative of ϕ
-williamson_dist(G::FrankGenerator, d) = G.θ > 0 ?  WilliamsonFromFrailty(Logarithmic(-G.θ), d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t),d)
+williamson_dist(G::FrankGenerator, ::Val{d}) where d = G.θ > 0 ?  WilliamsonFromFrailty(Logarithmic(-G.θ), Val(d)) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t), Val(d))
 
 Debye(x, k::Int=1) = k / x^k * QuadGK.quadgk(t -> t^k/expm1(t), 0, x)[1]
 function τ(G::FrankGenerator)
diff --git a/src/Generator/UnivariateGenerator/GumbelGenerator.jl b/src/Generator/UnivariateGenerator/GumbelGenerator.jl
@@ -59,4 +59,4 @@ function τ⁻¹(::Type{T},τ) where T<:GumbelGenerator
         return θ
     end
 end
-williamson_dist(G::GumbelGenerator, d) = WilliamsonFromFrailty(AlphaStable(α = 1/G.θ, β = 1,scale = cos(π/(2G.θ))^G.θ, location = (G.θ == 1 ? 1 : 0)), d)
+williamson_dist(G::GumbelGenerator, ::Val{d}) where d = WilliamsonFromFrailty(AlphaStable(α = 1/G.θ, β = 1,scale = cos(π/(2G.θ))^G.θ, location = (G.θ == 1 ? 1 : 0)), Val(d))
diff --git a/src/Generator/UnivariateGenerator/InvGaussianGenerator.jl b/src/Generator/UnivariateGenerator/InvGaussianGenerator.jl
@@ -73,4 +73,4 @@ function τ⁻¹(::Type{T}, tau) where T<:InvGaussianGenerator
     end
     return Roots.find_zero(x -> τ(InvGaussianGenerator(x)) - tau, (sqrt(eps(tau)), Inf))
 end
-williamson_dist(G::InvGaussianGenerator, d) = WilliamsonFromFrailty(Distributions.InverseGaussian(G.θ,1),d)
+williamson_dist(G::InvGaussianGenerator, ::Val{d}) where d = WilliamsonFromFrailty(Distributions.InverseGaussian(G.θ,1), Val(d))
diff --git a/src/Generator/UnivariateGenerator/JoeGenerator.jl b/src/Generator/UnivariateGenerator/JoeGenerator.jl
@@ -58,4 +58,4 @@ function τ⁻¹(::Type{T},tau) where T<:JoeGenerator
         return Roots.find_zero(θ -> τ(JoeGenerator(θ)) - tau, (one(tau),tau*Inf))
     end
 end
-williamson_dist(G::JoeGenerator, d) = WilliamsonFromFrailty(Sibuya(1/G.θ), d)
+williamson_dist(G::JoeGenerator, ::Val{d}) where d = WilliamsonFromFrailty(Sibuya(1/G.θ), Val(d))
diff --git a/src/Generator/WilliamsonGenerator.jl b/src/Generator/WilliamsonGenerator.jl
@@ -47,14 +47,14 @@ struct WilliamsonGenerator{TX} <: Generator
 end
 const i𝒲 = WilliamsonGenerator
 max_monotony(G::WilliamsonGenerator) = G.d
-function williamson_dist(G::WilliamsonGenerator, d)
+function williamson_dist(G::WilliamsonGenerator, ::Val{d}) where d
     if d == G.d 
         return G.X
     end
     # what about d < G.d ? Mayeb we can do some frailty stuff ? 
-    return WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t),d)
+    return WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t), Val(d))
 end
-ϕ(G::WilliamsonGenerator, t) = WilliamsonTransforms.𝒲(G.X,G.d)(t)
+ϕ(G::WilliamsonGenerator, t) = WilliamsonTransforms.𝒲(G.X, Val(G.d))(t)
 
 # McNeil & Neshelova 2009
-τ(G::WilliamsonGenerator) = 4*Distributions.expectation(x -> Copulas.ϕ(G,x), Copulas.williamson_dist(G,2))-1
+τ(G::WilliamsonGenerator) = 4*Distributions.expectation(x -> Copulas.ϕ(G,x), Copulas.williamson_dist(G, Val(2)))-1
diff --git a/src/UnivariateDistribution/WilliamsonFromFrailty.jl b/src/UnivariateDistribution/WilliamsonFromFrailty.jl
@@ -1,8 +1,9 @@
 struct WilliamsonFromFrailty{TF,d} <: Distributions.ContinuousUnivariateDistribution
     frailty_dist::TF
-    function WilliamsonFromFrailty(frailty_dist,d)
+    function WilliamsonFromFrailty(frailty_dist,::Val{d}) where d
         return new{typeof(frailty_dist),d}(frailty_dist)
     end
+    WilliamsonFromFrailty(frailty_dist,d) = WilliamsonFromFrailty(frailty_dist,Val(d))
 end
 function Distributions.rand(rng::Distributions.AbstractRNG, D::WilliamsonFromFrailty{TF,d}) where {TF,d}
     f = rand(rng,D.frailty_dist)
