Decouple rand and eltype

[devmotion]

One of my take-aways from issues such as #1252, #1902, #894, #1783, #1041, and #1071 is that eltype is not only quite inconsistently implemented and unclearly documented currently but more generally it might be a bad design decision to use it for pre-allocating containers of samples: Setting it to a fixed type (historically Float64 for continuous and Int for discrete distributions) is too limiting, but trying to infer it from parameters is challenging and doomed to fail for more challenging scenarios that are not as clear as e.g. Normal.

This PR tries to decouple rand and eltype, to make it easier to possibly eventually deprecate and remove eltype.

Fixes #1041.
Fixes #1071.
Fixes #1082.
Fixes #1252.
Fixes #1783.
Fixes #1884.
Fixes #1902.
Fixes #1907.

[codecov-commenter]

Codecov Report

❌ Patch coverage is 83.25243% with 69 lines in your changes missing coverage. Please review.
✅ Project coverage is 85.95%. Comparing base (6b5aa0d) to head (053d7af).
⚠️ Report is 2 commits behind head on master.

Files with missing lines	Patch %	Lines
src/mixtures/unigmm.jl	27.77%	13 Missing ⚠️
src/product.jl	69.44%	11 Missing ⚠️
src/multivariate/product.jl	30.00%	7 Missing ⚠️
src/multivariate/mvlogitnormal.jl	68.75%	5 Missing ⚠️
src/multivariate/mvtdist.jl	79.16%	5 Missing ⚠️
src/common.jl	55.55%	4 Missing ⚠️
src/genericrand.jl	71.42%	4 Missing ⚠️
src/multivariate/mvlognormal.jl	60.00%	4 Missing ⚠️
src/univariate/continuous/chisq.jl	40.00%	3 Missing ⚠️
src/matrix/lkj.jl	95.12%	2 Missing ⚠️
... and 7 more

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1905      +/-   ##
==========================================
- Coverage   86.74%   85.95%   -0.80%     
==========================================
  Files         149      149              
  Lines        8882     8947      +65     
==========================================
- Hits         7705     7690      -15     
- Misses       1177     1257      +80     

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

[andreasnoack]

Given that, according to the docstring

Distributions.jl/src/common.jl

Lines 96 to 102 in a1010e4

	"""
	eltype(::Type{Sampleable})
	The default element type of a sample. This is the type of elements of the samples generated
	by the `rand` method. However, one can provide an array of different element types to
	store the samples using `rand!`.
	"""

the sole purpose of eltype is to return the type of rand, which I agree isn't really feasible, should we also deprecate the methods here as part of this PR?

[quildtide]

#1907 is related.

[devmotion]

#1907 is related.

17154a2 fixes it.

[quildtide]

My personal opinion:
rand(d::Distribution) should remain a valid method (I don't think anyone is trying to change this), and its return type should be predictable by knowing what d is. It would also be useful to have a method that tells you what type rand(d::Distribution) returns without running rand(d::Distribution). eltype(d::Distribution) is logical for this.

This is orthogonal to the fact that Distributions should probably support rand(d::Distribution, T::Type).

In an ideal situation, I think we would define rand(d::Distribution, T::Type) for distributions (to allow for potential type-specific optimizations and evade unnecessary typecasts) and then dispatch rand(d::Distribution) = rand(d, eltype(d)).

[devmotion]

its return type should be predictable by knowing what d is

I became convinced that in general this is not possible. It basically means trying to manually figure out type inference, a task that even the Julia compiler can fail at. Of course, for some number types and distributions it can be done (as also the compiler will succeed in many instances). But in general it's far from trivial. As in the initial stages of Distributions, one could be restrictive and limit samples to Int and Float but that's IMO not desirable either. Therefore I think it's best to remove eltype (which IMO has also been a weird name since I don't view distributions as containers or collections) from the API. IMO it's just too brittle and currently to inconsistent. Of course, that doesn't prevent users from trying to figure out the (element) type of samples with eg the help of the compiler, it's just not an official feature anymore.

[singularitti]

A few distributions still give me Float64 after this PR, while others work fine:

julia> dist = Semicircle(50.0f0)
Semicircle{Float32}(r=50.0f0)

julia> rand(dist, 5)
5-element Array{Float64, 1}:
  36.80671953487414
 -18.355635129900335
 -11.701855436648922
 -21.444118928985656
  -5.80120463505302

julia> dist = JohnsonSU(0.0f0, 1.0f0, 0.0f0, 1.0f0)
JohnsonSU{Float32}(ξ=0.0f0, λ=1.0f0, γ=0.0f0, δ=1.0f0)

julia> rand(dist, 5)
5-element Array{Float64, 1}:
  0.5311696707298299
  1.632313034117999
  0.04951771555318912
  0.4721610259428258
 -3.052321854866766

julia> dist = Chisq(5.0f0)
Chisq{Float32}(ν=5.0f0)

julia> rand(dist, 5)
5-element Array{Float32, 1}:
 15.465032
 1.888659
 7.013455
 4.258529
 3.9611576

Can this be fixed?

[devmotion]

The reason is that for quite a few of the older, probably less used and definitely less frequently updated distributions the rand implementation contains hardcoded Float64, either as literals or as calls of rand(rng) or randn(rng). I fixed the SemicircleandJohnsonSU`, but I'm not sure if it's possible (or even desirable) to address all of these in a single PR (it's already a quite massive undertaking). Maybe it would be best to open issues for others and address them in follow-up PRs (most (all?) other such changes in this PR are included to fix open issues).

[dom-linkevicius]

eltype issues described here are also quite painful when using parameters of Dual type in a distribution d. This makes rand(d, 10) throw due to pre-allocation based on eltype (which is often Float64), as well as eltype just being plain incorrect (doing rand(d)) since Duals actually get returned.

[mschauer]

Sorry, late to the party, but with iterators we have a "HasEltype" attribute to address exactly this issue, how is this design related?

[sethaxen]

Here's a review of censored, Cholesky-, and matrix-variate distributions.

[nilsbecker]

I also believe that it is reasonable that quantile(d, x) should return eltype(d); it is analogous to getindex to some extent. When would it make sense for rand and quantile to return different types?

Sometimes quantiles of integer-valued distributions are taken to interpolate between the integers flanking the 50% value of the cumulative. E.g. the median of an unbiased Bernoulli distribution would be 0.5 in this convention. Just saying.

[sethaxen]

Here's a review of the mixture models and multivariate distributions.

[sethaxen]

According to the docs for map!:

Behavior can be unexpected when any mutated argument shares memory with any other argument.

Maybe we could broadcast here instead?

Suggested change

	map!(exp, xs, xs)
	xs .= exp.(xs)

[devmotion]

I'm not sure. Generally IMO it's better to avoid broadcasting unless it's necessary (ie unless you're actually broadcasting objects with different dimensions). Broadcasting is slow, both at compile and runtime.

[sethaxen]

Suggested change

	map!(exp, x, x)
	x .= exp.(x)

[sethaxen]

Suggested change

	map!(exp, x, x)
	x .= exp.(x)

[sethaxen]

Why not broadcast the rand call to avoid allocating an intermediate container?

Suggested change

	y = similar(x, (1, cols))
	unwhiten!(d.Σ, randn!(rng, x))
	rand!(rng, chisqd, y)
	x .= x ./ sqrt.(y ./ d.df) .+ d.μ
	unwhiten!(d.Σ, randn!(rng, x))
	rand!(rng, chisqd, y)
	x .= x ./ sqrt.(rand.(rng, chisqd) ./ d.df) .+ d.μ

[devmotion]

We could, but I chose this version on purpose: Generally, vectorized sampling tends to be much faster.

[bgroenks96]

Hey guys, what's the status of this? I just ran into #1089 and #1907 when trying to do rand(truncated(Normal(0f0, 1f0), 0f0, Inf32) which fails due to the mismatched types. Given that lower precision floats are becoming increasingly widespread, this kind of seems like a major issue?