Make boundary constructors type stable

[lrnv]

Enhance type stability for boundary constructors by implementing reduced tail and generator functions. This change improves performance and reliability in handling various tail types.

Fixes #333

[github-actions]

🟢 Tachometer — no performance regressions detected (1 improvement 🟢)

25 benchmarks compared · a5babd5 → 6724d35 · Julia 1.12.7 · 5% tolerance · 3× runs

	Benchmark	Time	Memory
🟢	fitting/sklar_ifm	53.6 ms → 16.6 ms (−69%)	37.9 MiB → 987 KiB (−97%)

Full results (25 benchmarks)

	Benchmark	Time	Memory
⚪️	cdf/bb1_d2	1.59 ms → 1.61 ms (+1%)	—
⚪️	cdf/galambos_d2	1.12 ms → 1.14 ms (+2%)	—
⚪️	conditioning/cdf_bb1_d2	340 µs → 336 µs (−1%)	—
⚪️	conditioning/construct_bb1_d2_1000	344 µs → 343 µs (−0%)	—
⚪️	conditioning/inverse_rosenblatt_gaussian_d5	263 µs → 264 µs (+0%)	—
⚪️	conditioning/quantile_galambos_d2	11.5 ms → 11.5 ms (+0%)	—
⚪️	conditioning/rosenblatt_gaussian_d5	296 µs → 297 µs (+0%)	—
⚪️	data/beta_logpdf	10.3 ms → 10.4 ms (+0%)	—
⚪️	data/checkerboard_cdf	54.8 ms → 53.8 ms (−2%)	—
⚪️	data/empirical_cdf	15.2 ms → 15.2 ms (−0%)	—
⚪️	data/pseudos_5x10000	1.96 ms → 1.96 ms (−0%)	—
⚪️	density/bb1_d2	1.86 ms → 1.88 ms (+1%)	—
⚪️	density/galambos_d2	2.53 ms → 2.57 ms (+2%)	—
⚪️	density/gaussian_d10	8.5 ms → 8.55 ms (+1%)	—
⚪️	density/gumbel_d5	9.09 ms → 9.13 ms (+0%)	—
⚪️	density/nested_d6	10.5 ms → 10.5 ms (−0%)	—
⚪️	fitting/gaussian_mle	55 µs → 54.8 µs (−0%)	—
⚪️	fitting/gumbel_itau	141 µs → 142 µs (+1%)	—
🟢	fitting/sklar_ifm	53.6 ms → 16.6 ms (−69%)	37.9 MiB → 987 KiB (−97%)
⚪️	sampling/archimax_d2	16.3 ms → 16.4 ms (+0%)	—
⚪️	sampling/clayton_d5	2.3 ms → 2.3 ms (−0%)	—
⚪️	sampling/galambos_d2	80.6 ms → 80.8 ms (+0%)	—
⚪️	sampling/gaussian_d10	1.73 ms → 1.74 ms (+0%)	—
⚪️	sampling/nested_d6	122 ms → 124 ms (+1%)	—
⚪️	sampling/student_d2	4.42 ms → 4.49 ms (+1%)	—

_{baseline a5babd5 · target 6724d35 · Julia 1.12.7 · ARM Neoverse-N2 · min estimator · tol 5%/5% · floor 1 µs · 3× runs · run · Tachometer.jl}

[lrnv]

I don't think the current explosion in the number of methods is actually necessary.

Details

Some additional boundary logic is unavoidable once constructors become type-stable, but I think #437 currently puts too much of that logic too high in the stack.

The clearest symptom is ArchimaxCopula.jl: the PR adds separate specializations for copula_measure_style, _cdf, _logpdf, _rand!, DistortionFromCop, etc., then has to resolve intersections such as IndependentGenerator × NoTail/MTail, and BB4/BB5 add yet another layer of operation-specific delegation.

That looks like a sign that we are re-encoding the constructor-reduction graph inside each public operation.

What is actually unavoidable

If we want

ClaytonCopula{2}(θ)::ClaytonCopula{2}

for every admissible θ, and θ is a runtime value, then Julia dispatch can no longer distinguish

θ = 1.5
θ = 0.0
θ = Inf

from the type of C.

So for a limit at which the internal representation genuinely degenerates, there has to be a runtime branch somewhere. There is no dispatch trick that removes that requirement without making the return type value-dependent again.

But that branch does not have to be duplicated across 8 operations × 20 families.

I would push boundary handling down to the lowest mathematical layer possible

I think the reduction graph actually contains three very different kinds of boundaries.

1. Boundaries that are genuine extensions of the representation

For example:

Clayton(0)   -> Independent
Gumbel(1)    -> Independent
Frank(0)     -> Independent
Galambos(0)  -> Independent
Log(1)       -> Independent

I would not treat those as operational reductions at all.

For instance, ClaytonGenerator(0) should simply satisfy the correct limiting identities:

ϕ(G, t)   == exp(-t)
ϕ⁻¹(G, u) == -log(u)
𝒲₋₁(G, d) == Gamma(d, 1)

and then the generic ArchimedeanCopula implementation should keep working.

Likewise for EV tails: if GalambosTail(0) represents independence, then its STDF/Pickands primitives and derivatives should have the correct extension at θ = 0. ExtremeValueCopula should ideally not even need to know that this is a boundary value.

This is especially important for the goal of the PR: Clayton(0) should really behave as a Clayton copula, not merely be an object that delegates all of its operations to IndependentCopula.

2. Compositional reductions

For example:

BB1(θ, 1)         -> Clayton(θ)
BB6(...)          -> Joe / Gumbel
AsymLog(...)      -> Log
AsymGalambos(...) -> Galambos
BB4               -> Clayton / Galambos
BB5               -> Gumbel / Galambos

Here again, I would strongly avoid turning the reduction graph into

constructor graph
      ↓
cdf graph
logpdf graph
rand graph
conditioning graph
...

which is what the PR is starting to do.

For example, BB4 now effectively does things like

iszero(δ) && return cdf(ClaytonCopula{2}(θ), u)
iszero(θ) && return cdf(GalambosCopula{2}(δ), u)

and then repeats the same idea in logpdf, _rand!, etc. BB5 follows the same pattern.

That feels like a code smell.

BB4 is Archimax(ClaytonGenerator, GalambosTail). If the Clayton generator and Galambos tail are both correctly defined at their boundaries, then ideally the generic Archimax machinery should produce the right result by itself.

In other words:

An identity between composite models should preferably emerge from the primitives of their components, rather than being reimplemented operation by operation.

3. Genuine degeneracies: Inf -> M, sometimes -> W

This is where I think a small central runtime mechanism is justified.

Examples:

Gumbel(Inf)  -> M
Joe(Inf)     -> M
Clayton(Inf) -> M
Frank(Inf)   -> M
Frank(-Inf)  -> W

At Inf, there is not necessarily a clean generator/frailty/Williamson representation that can be extended in a meaningful way.

Here I would introduce a tiny concrete runtime trait, but importantly not something like the first implementation in #437 that returned

nothing | NoTail | MTail

and then reconstructed another copula and redispatched.

Something along these lines:

@enum LimitKind::UInt8 begin
    NO_LIMIT
    M_LIMIT
    W_LIMIT
end

@inline limit_kind(::Generator, ::Val) = NO_LIMIT

@inline function limit_kind(G::ClaytonGenerator, ::Val{d}) where {d}
    isinf(G.θ) && return M_LIMIT
    d == 2 && G.θ == -1 && return W_LIMIT
    return NO_LIMIT
end

@inline function limit_kind(G::FrankGenerator, ::Val{2})
    G.θ == Inf  && return M_LIMIT
    G.θ == -Inf && return W_LIMIT
    return NO_LIMIT
end

@inline limit_kind(G::GumbelGenerator, ::Val) =
    isinf(G.θ) ? M_LIMIT : NO_LIMIT

@inline limit_kind(G::JoeGenerator, ::Val) =
    isinf(G.θ) ? M_LIMIT : NO_LIMIT

The important part is that limit_kind always returns exactly LimitKind.

Not:

Union{Nothing, IndependentCopula, MCopula, ...}

So there is no reconstruction of another model and no dynamic redispatch over a union of unrelated types.

Then, at the Archimedean level, this can be handled once:

@inline function _cdf(C::ArchimedeanCopula{d}, u) where {d}
    kind = limit_kind(C.gen, Val(d))

    kind === M_LIMIT && return minimum(u)

    if kind === W_LIMIT
        @assert d == 2
        return max(u[1] + u[2] - 1, zero(eltype(u)))
    end

    return _archimedean_cdf(C, u)
end

And the same idea can be used for sampling:

function Distributions._rand!(rng, C::ArchimedeanCopula{d}, A) where {d}
    kind = limit_kind(C.gen, Val(d))

    kind === M_LIMIT &&
        return _rand_M!(rng, A)

    kind === W_LIMIT &&
        return _rand_W!(rng, A)

    return _rand_archimedean!(rng, C, A)
end

This gives us roughly one branch per structural Archimedean operation, instead of one specialization per family × boundary × operation.

A UInt8-backed enum plus a predictable branch is also very different from the earlier strategy that returned different object types and redispatched. It still needs to be benchmarked, of course, but it is much friendlier to inference.

The same principle applies to EV families

For example:

@inline limit_kind(::Tail) = NO_LIMIT

@inline limit_kind(T::GalambosTail) =
    isinf(T.θ) ? M_LIMIT : NO_LIMIT

@inline limit_kind(T::LogTail) =
    isinf(T.θ) ? M_LIMIT : NO_LIMIT

But I would not put Galambos(0) into limit_kind if we can make ℓ(GalambosTail(0), x) naturally equal to sum(x).

The trait should mean:

the generic mathematical representation genuinely stops being usable here

not:

the test graph happens to know that this boundary is equivalent to another family.

That distinction matters a lot.

This would also clarify the role of CONSTRUCTOR_REDUCTIONS

The graph contains a lot of identities: simple Archimedean families, BB families, EV families, asymmetric EV families, Archimax, Liouville, Nested, and miscellaneous constructors.

But I don't think that graph should become a specification of internal dispatch.

It should remain only a correctness oracle:

source and target must have the same semantics at this boundary.

So:

test_boundary_equivalence(source, target)

can tell us that rand fails at some boundary.

Then we inspect why.

If it is Clayton(0), the fix should be something like:

fix 𝒲₋₁(ClaytonGenerator(0), d).

Not:

add _rand!(::ClaytonCopula{d}) that detects zero and constructs an IndependentCopula.

If it is Gumbel(Inf), then the representation really does degenerate, so an architectural M shortcut is justified.

I think following that discipline would remove a large amount of the method growth we currently see.

Archimax is probably the best place to test this approach first

The current PR adds variants corresponding to things like:

Archimax{Generator, NoTail}             -> Archimedean
Archimax{IndependentGenerator, T}       -> ExtremeValue
Archimax{IndependentGenerator, NoTail}  -> Independent
Archimax{..., MTail}                    -> M

for measure_style, cdf, logpdf, rand!, conditioning, etc., plus intersection cases.

I think there are two opportunities here.

First, IndependentGenerator could become a real mathematical generator rather than only a reduction marker:

ϕ(::IndependentGenerator, t) = exp(-t)
ϕ⁻¹(::IndependentGenerator, u) = -log(u)
# derivatives, Williamson transform, etc.

Then Archimax(IndependentGenerator(), tail) should ideally become an EV model through the Archimax formulas themselves, rather than because cdf, rand!, conditioning, and other operations all explicitly convert it to ExtremeValueCopula.

Likewise, NoTail is a genuine independence STDF. It should be complete enough that the generic formulas naturally reduce to the Archimedean case.

MTail is different: it may require a few singular paths, but those paths should live at the Archimax/EV engine level, not be repeated in BB4 and BB5.

What I would keep from #437

I would keep:

1. type-stable constructors
2. the new documented constructor contract
3. CONSTRUCTOR_REDUCTIONS
4. test_boundary_equivalence
5. mathematical extensions of Generator/Tail primitives at their boundaries

Then I would revisit the operation-specific shortcuts that have accumulated in the PR.

The rule I would use is:

Handle a boundary at the lowest mathematical layer that can express the limit correctly.

And I would debug failures in this order:

Generator / Tail primitive
        ↓ if impossible
Generic Archimedean / EV / Archimax engine
        ↓ if genuinely family-specific
Family operation

rather than:

reduction graph
        ↓
special method for every public operation

I suspect that would let us remove a substantial fraction of the methods currently added by #437, especially the BB4/BB5 forwarding methods and a significant part of the Archimax-specific duplication.

The nice thing is that, once test_boundary_equivalence is reasonably complete, we can test this mechanically: remove the current shortcuts one by one and see which ones are actually required. That should give us a much cleaner way to simplify the PR than trying to predict all necessary special cases up front.

[codecov]

Codecov Report

❌ Patch coverage is 90.50555% with 77 lines in your changes missing coverage. Please review.
✅ Project coverage is 85.34%. Comparing base (a5babd5) to head (6724d35).
⚠️ Report is 5 commits behind head on main.

Files with missing lines	Patch %	Lines
src/Generator/AMHGenerator.jl	64.28%	35 Missing ⚠️
src/Generator.jl	75.00%	10 Missing ⚠️
src/Tail/MTail.jl	0.00%	5 Missing ⚠️
src/ArchimedeanCopula.jl	89.47%	4 Missing ⚠️
src/ExtremeValueCopula.jl	95.77%	3 Missing ⚠️
src/MiscellaneousCopulas/SurvivalCopula.jl	94.44%	3 Missing ⚠️
src/Generator/ClaytonGenerator.jl	95.23%	2 Missing ⚠️
src/MiscellaneousCopulas/FGMCopula.jl	85.71%	2 Missing ⚠️
src/Tail/NoTail.jl	75.00%	2 Missing ⚠️
...rc/UnivariateDistribution/Frailties/AlphaStable.jl	0.00%	2 Missing ⚠️
... and 8 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #437      +/-   ##
==========================================
+ Coverage   85.32%   85.34%   +0.02%     
==========================================
  Files          90       90              
  Lines        7984     8310     +326     
==========================================
+ Hits         6812     7092     +280     
- Misses       1172     1218      +46     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.