Fix deprecations

src/genericrand.jl

@@ -26,6 +26,7 @@ rand(rng::AbstractRNG, s::Sampleable, dim1::Int, moredims::Int...) =
 
 # default fallback (redefined for univariate distributions)
 function rand(rng::AbstractRNG, s::Sampleable{<:ArrayLikeVariate})
+    Base.depwarn("Please implement `rand(rng::AbstractRNG, s::$(typeof(s)))`. The default fallback will be removed", :rand)
     return @inbounds rand!(rng, s, Array{eltype(s)}(undef, size(s)))
 end
 
@@ -45,7 +46,8 @@ end
 
 # this is a workaround for sampleables that incorrectly base `eltype` on the parameters
 function rand(rng::AbstractRNG, s::Sampleable{<:ArrayLikeVariate,Continuous})
-    return @inbounds rand!(rng, sampler(s), Array{float(eltype(s))}(undef, size(s)))
+    Base.depwarn("Please implement `rand(rng::AbstractRNG, s::$(typeof(s))`. The default fallback will be removed", :rand)
+    return @inbounds rand!(rng, s, Array{float(eltype(s))}(undef, size(s)))
 end
 
 """
@@ -63,16 +65,14 @@ form as specified above. The rules are summarized as below:
 """
 function rand! end
 Base.@propagate_inbounds rand!(s::Sampleable, X::AbstractArray) = rand!(default_rng(), s, X)
-Base.@propagate_inbounds function rand!(rng::AbstractRNG, s::Sampleable, X::AbstractArray)
-    return _rand!(rng, s, X)
-end
 
 # default definitions for arraylike variates
 @inline function rand!(
     rng::AbstractRNG,
     s::Sampleable{ArrayLikeVariate{N}},
     x::AbstractArray{<:Real,N},
 ) where {N}
+    Base.depwarn("Please implement `Random.rand!(rng::Random.AbstractRNG, s::$(typeof(s)), x::AbstractArray{<:Real,$N})`, the default fallback will be removed.", :rand!)
     @boundscheck begin
         size(x) == size(s) || throw(DimensionMismatch("inconsistent array dimensions"))
     end
@@ -93,7 +93,8 @@ end
             throw(DimensionMismatch("inconsistent array dimensions"))
     end
     # the function barrier fixes performance issues if `sampler(s)` is type unstable
-    return _rand!(rng, sampler(s), x)
+    _rand!(rng, sampler(s), x)
+    return x
 end
 
 function _rand!(

src/mixtures/mixturemodel.jl

@@ -470,16 +470,14 @@ end
 
 Base.length(s::MixtureSampler) = length(first(s.csamplers))
 
-rand(rng::AbstractRNG, s::MixtureSampler{Univariate}) =
+# mixture sampler
+rand(rng::AbstractRNG, s::MixtureSampler) =
     rand(rng, s.csamplers[rand(rng, s.psampler)])
-rand(rng::AbstractRNG, d::MixtureModel{Univariate}) =
-    rand(rng, component(d, rand(rng, d.prior)))
 
-# multivariate mixture sampler for a vector
-_rand!(rng::AbstractRNG, s::MixtureSampler{Multivariate}, x::AbstractVector{<:Real}) =
-    @inbounds rand!(rng, s.csamplers[rand(rng, s.psampler)], x)
 # if only a single sample is requested, no alias table is created
-_rand!(rng::AbstractRNG, d::MixtureModel{Multivariate}, x::AbstractVector{<:Real}) =
-    @inbounds rand!(rng, component(d, rand(rng, d.prior)), x)
+rand(rng::AbstractRNG, s::MixtureModel) = rand(rng, s.csamplers[rand(rng, s.psampler)])
+Base.@propagate_inbounds function rand!(rng::AbstractRNG, d::MixtureModel{ArrayLikeVariate{N}}, x::AbstractArray{<:Real,N}) where {N}
+    return rand!(rng, component(d, rand(rng, d.prior)), x)
+end
 
 sampler(d::MixtureModel) = MixtureSampler(d)

src/mixtures/unigmm.jl

@@ -25,10 +25,12 @@ probs(d::UnivariateGMM) = probs(d.prior)
 
 mean(d::UnivariateGMM) = dot(d.means, probs(d))
 
-rand(d::UnivariateGMM) = (k = rand(d.prior); d.means[k] + randn() * d.stds[k])
-
-rand(rng::AbstractRNG, d::UnivariateGMM) =
-    (k = rand(rng, d.prior); d.means[k] + randn(rng) * d.stds[k])
+function rand(rng::AbstractRNG, d::UnivariateGMM)
+    k = rand(rng, d.prior)
+    μ = d.means[k]
+    σ = d.std[k]
+    return muladd(randn(rng, float(Base.promote_typeof(μ, σ))), σ, μ)
+end
 
 params(d::UnivariateGMM) = (d.means, d.stds, d.prior)
 
@@ -38,6 +40,22 @@ struct UnivariateGMMSampler{VT1<:AbstractVector{<:Real},VT2<:AbstractVector{<:Re
     psampler::AliasTable
 end
 
-rand(rng::AbstractRNG, s::UnivariateGMMSampler) =
-    (k = rand(rng, s.psampler); s.means[k] + randn(rng) * s.stds[k])
+function rand(rng::AbstractRNG, s::UnivariateGMMSampler)
+    k = rand(rng, s.psampler)
+    μ = d.means[k]
+    σ = d.stds[k]
+    return muladd(randn(rng, float(Base.promote_typeof(μ, σ))), σ, μ)
+end
+function rand(rng::AbstractRNG, s::UnivariateGMMSampler, x::AbstractArray{<:Real})
+    psampler = s.psampler
+    means = s.means
+    stds = s.stds
+    randn!(rng, x)
+    for i in eachindex(x)
+        k = rand(rng, psampler)
+        x[i] = muladd(x[i], stds[k], means[k])
+    end
+    return x
+end
+
 sampler(d::UnivariateGMM) = UnivariateGMMSampler(d.means, d.stds, sampler(d.prior))

src/multivariate/dirichletmultinomial.jl

@@ -99,8 +99,10 @@ end
 # Sampling
 rand(rng::AbstractRNG, d::DirichletMultinomial) =
     multinom_rand(rng, ntrials(d), rand(rng, Dirichlet(d.α)))
-_rand!(rng::AbstractRNG, d::DirichletMultinomial, x::AbstractVector{<:Real}) =
+@inline function rand!(rng::AbstractRNG, d::DirichletMultinomial, x::AbstractVector{<:Real})
+    @checkbounds length(d) == length(x)
     multinom_rand!(rng, ntrials(d), rand(rng, Dirichlet(d.α)), x)
+end
 
 # Fit Model
 # Using https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2945396/pdf/nihms205488.pdf

src/multivariate/multinomial.jl

@@ -166,7 +166,7 @@ end
 
 # if only a single sample is requested, no alias table is created
 rand(rng::AbstractRNG, d::Multinomial) = multinom_rand(rng, ntrials(d), probs(d))
-_rand!(rng::AbstractRNG, d::Multinomial, x::AbstractVector{<:Real}) =
+rand!(rng::AbstractRNG, d::Multinomial, x::AbstractVector{<:Real}) =
     multinom_rand!(rng, ntrials(d), probs(d), x)
 
 sampler(d::Multinomial) = MultinomialSampler(ntrials(d), probs(d))

src/multivariate/mvlogitnormal.jl

@@ -99,11 +99,15 @@ function rand(rng::AbstractRNG, d::MvLogitNormal, n::Int)
     return x
 end
 
-function _rand!(rng::AbstractRNG, d::MvLogitNormal, x::AbstractVecOrMat{<:Real})
-    y = @views _drop1(x)
-    rand!(rng, d.normal, y)
-    _softmax1!(x, y)
-    return x
+for N in (1, 2)
+    @eval begin
+        Base.@propagate_inbounds function rand!(rng::AbstractRNG, d::MvLogitNormal, x::AbstractArray{<:Real, $N})
+            y = @views _drop1(x)
+            rand!(rng, d.normal, y)
+            _softmax1!(x, y)
+            return x
+        end
+    end
 end
 
 # Fitting

src/multivariate/mvlognormal.jl

@@ -241,10 +241,14 @@ function rand(rng::AbstractRNG, d::MvLogNormal, n::Int)
     return xs
 end
 
-function _rand!(rng::AbstractRNG, d::MvLogNormal, x::AbstractVecOrMat{<:Real})
-    _rand!(rng, d.normal, x)
-    map!(exp, x, x)
-    return x
+for N in (1, 2)
+    @eval begin
+        Base.@propagate_inbounds function rand!(rng::AbstractRNG, d::MvLogNormal, x::AbstractArray{<:Real,$N})
+            rand!(rng, d.normal, x)
+            map!(exp, x, x)
+            return x
+        end
+    end
 end
 
 _logpdf(d::MvLogNormal, x::AbstractVecOrMat{T}) where {T<:Real} = insupport(d, x) ? (_logpdf(d.normal, log.(x)) - sum(log.(x))) : -Inf

src/multivariate/mvnormal.jl

@@ -282,16 +282,16 @@ function rand(rng::AbstractRNG, d::MvNormal, n::Int)
     return x
 end
 
-function _rand!(rng::AbstractRNG, d::MvNormal, x::VecOrMat)
+Base.@propagate_inbounds function rand!(rng::AbstractRNG, d::MvNormal, x::VecOrMat{<:Real})
     unwhiten!(d.Σ, randn!(rng, x))
     x .+= d.μ
     return x
 end
 
 # Workaround: randn! only works for Array, but not generally for AbstractArray
-function _rand!(rng::AbstractRNG, d::MvNormal, x::AbstractVector)
+Base.@propagate_inbounds function rand!(rng::AbstractRNG, d::MvNormal, x::AbstractVector{<:Real})
     for i in eachindex(x)
-        @inbounds x[i] = randn(rng, eltype(x))
+        x[i] = randn(rng, eltype(x))
     end
     unwhiten!(d.Σ, x)
     x .+= d.μ

src/multivariate/product.jl

@@ -32,8 +32,11 @@ end
 
 length(d::Product) = length(d.v)
 
-_rand!(rng::AbstractRNG, d::Product, x::AbstractVector{<:Real}) =
-    map!(Base.Fix1(rand, rng), x, d.v)
+rand(rng::AbstractRNG, d::Product) = map(Base.Fix1(rand, rng), d.v)
+Base.@propagate_inbounds function rand!(rng::AbstractRNG, d::Product, x::AbstractVector{<:Real})
+    return map!(Base.Fix1(rand, rng), x, d.v)
+end
+
 function _logpdf(d::Product, x::AbstractVector{<:Real})
     dists = d.v
     if isempty(dists)

src/multivariate/vonmisesfisher.jl

@@ -77,11 +77,10 @@ _logpdf(d::VonMisesFisher, x::AbstractVector{T}) where {T<:Real} = d.logCκ + d.
 
 sampler(d::VonMisesFisher) = VonMisesFisherSampler(d.μ, d.κ)
 
-_rand!(rng::AbstractRNG, d::VonMisesFisher, x::AbstractVector) =
-    _rand!(rng, sampler(d), x)
-_rand!(rng::AbstractRNG, d::VonMisesFisher, x::AbstractMatrix) =
-    _rand!(rng, sampler(d), x)
-
+rand(rng::AbstractRNG, d::VonMisesFisher) = rand(rng, sampler(d))
+Base.@propagate_inbounds function rand!(rng::AbstractRNG, d::VonMisesFisher, x::AbstractVector)
+    return rand!(rng, sampler(d), x)
+end
 
 ### Estimation
 

src/product.jl

@@ -7,15 +7,15 @@ independent `M`-dimensional distributions by stacking them.
 Users should use [`product_distribution`](@ref) to construct a product distribution of
 independent distributions instead of constructing a `ProductDistribution` directly.
 """
-struct ProductDistribution{N,M,D,S<:ValueSupport,T} <: Distribution{ArrayLikeVariate{N},S}
+struct ProductDistribution{N,M,D,S<:ValueSupport} <: Distribution{ArrayLikeVariate{N},S}
     dists::D
     size::Dims{N}
 
     function ProductDistribution{N,M,D}(dists::D) where {N,M,D}
         if isempty(dists)
             throw(ArgumentError("a product distribution must consist of at least one distribution"))
         end
-        return new{N,M,D,_product_valuesupport(dists),_product_eltype(dists)}(
+        return new{N,M,D,_product_valuesupport(dists)}(
             dists,
             _product_size(dists),
         )
@@ -32,15 +32,11 @@ end
 
 # default definitions (type stable e.g. for arrays with concrete `eltype`)
 _product_valuesupport(dists) = mapreduce(value_support ∘ typeof, promote_type, dists)
-_product_eltype(dists) = mapreduce(eltype, promote_type, dists)
 
 # type-stable and faster implementations for tuples
 function _product_valuesupport(dists::NTuple{<:Any,Distribution})
     return __product_promote_type(value_support, typeof(dists))
 end
-function _product_eltype(dists::NTuple{<:Any,Distribution})
-    return __product_promote_type(eltype, typeof(dists))
-end
 
 __product_promote_type(f::F, ::Type{Tuple{D}}) where {F,D<:Distribution} = f(D)
 function __product_promote_type(f::F, ::Type{T}) where {F,T}
@@ -67,7 +63,7 @@ const VectorOfUnivariateDistribution{D,S<:ValueSupport,T} = ProductDistribution{
 const MatrixOfUnivariateDistribution{D,S<:ValueSupport,T} = ProductDistribution{2,0,D,S,T}
 const ArrayOfUnivariateDistribution{N,D,S<:ValueSupport,T} = ProductDistribution{N,0,D,S,T}
 
-const FillArrayOfUnivariateDistribution{N,D<:Fill{<:Any,N},S<:ValueSupport,T} = ProductDistribution{N,0,D,S,T}
+const FillArrayOfUnivariateDistribution{N,D<:FillArrays.AbstractFill{<:Any,N},S<:ValueSupport,T} = ProductDistribution{N,0,D,S,T}
 
 ## General definitions
 size(d::ProductDistribution) = d.size
@@ -104,8 +100,11 @@ length(d::VectorOfUnivariateDistribution) = length(d.dists)
 ## For matrix distributions
 cov(d::ProductDistribution{2}, ::Val{false}) = reshape(cov(d), size(d)..., size(d)...)
 
-# `_rand!` for arrays of univariate distributions
-function _rand!(
+# Arrays of univariate distributions
+function rand(rng::AbstractRNG, d::ArrayOfUnivariateDistribution)
+    return map(Base.Fix1(rand, rng), d.dists)
+end
+function rand!(
     rng::AbstractRNG,
     d::ArrayOfUnivariateDistribution{N},
     x::AbstractArray{<:Real,N},
@@ -129,8 +128,14 @@ function __logpdf(d::ArrayOfUnivariateDistribution, x::AbstractArray{<:Real,N})
     return sum(Broadcast.instantiate(broadcasted))
 end
 
-# more efficient implementation of `_rand!` for `Fill` array of univariate distributions
-function _rand!(
+# more efficient sampling for `Fill` array of univariate distributions
+function rand(
+    rng::AbstractRNG,
+    d::FillArrayOfUnivariateDistribution,
+)
+    return @inbounds rand(rng, sampler(first(d.dists)), size(d))
+end
+function rand!(
     rng::AbstractRNG,
     d::FillArrayOfUnivariateDistribution{N},
     x::AbstractArray{<:Real,N},
@@ -152,13 +157,19 @@ function __logpdf(
     return @inbounds loglikelihood(first(d.dists), x)
 end
 
-# `_rand! for arrays of distributions
-function _rand!(
+# sampling for arrays of distributions
+function rand(rng::AbstractRNG, d::ProductDistribution)
+    x = let rng = rng, d = d
+        mapreduce(di -> vec(rand(rng, di)), hcat, vec(d.dists))
+    end
+    return reshape(x, size(d))
+end
+Base.@propagate_inbounds function rand!(
     rng::AbstractRNG,
     d::ProductDistribution{N,M},
     A::AbstractArray{<:Real,N},
 ) where {N,M}
-    @inbounds for (di, Ai) in zip(d.dists, eachvariate(A, ArrayLikeVariate{M}))
+    for (di, Ai) in zip(d.dists, eachvariate(A, ArrayLikeVariate{M}))
         rand!(rng, di, Ai)
     end
     return A
@@ -180,32 +191,48 @@ function __logpdf(
     return sum(Broadcast.instantiate(broadcasted))
 end
 
-# more efficient implementation of `_rand!` for `Fill` arrays of distributions
+# more efficient sampling for `Fill` arrays of distributions
+function rand(rng::AbstractRNG, d::ProductDistribution{<:Any,<:Any,<:FillArrays.AbstractFill})
+    rand(rng, sampler(first(d.dists)), A)
+    x = let rng = rng, d = d
+        mapreduce(di -> vec(rand(rng, di)), hcat, vec(d.dists))
+    end
+    return reshape(x, size(d))
+end
+function _product_rand(rng::AbstractRNG, spl::Sampleable{ArrayLikeVariate{N}}, dims::Dims) where N
+    xi = rand(rng, spl)
+    x = Array{eltype(xi)}(undef, dims)
+    copyto!(x, xi)
+    vx = reshape(x, ntuple(i -> i <= N ? size(xi, i) : Colon(), N + 1))
+    @inbounds rand!(rng, spl, @view(vx[ntuple(i -> i <= N ? Colon() : 2:lastindex(vx, N + 1), N + 1)...]))
+    return x
+end
+
 function _rand!(
     rng::AbstractRNG,
-    d::ProductDistribution{N,M,<:Fill},
+    d::ProductDistribution{N,M,<:FillArrays.AbstractFill},
     A::AbstractArray{<:Real,N},
 ) where {N,M}
     @inbounds rand!(rng, sampler(first(d.dists)), A)
     return A
 end
 
-# more efficient implementation of `_logpdf` for `Fill` arrays of distributions
+# more efficient implementation of `_logpdf` for `AbstractFill` arrays of distributions
 # we have to fix a method ambiguity
 function _logpdf(
-    d::ProductDistribution{N,M,<:Fill},
+    d::ProductDistribution{N,M,<:FillArrays.AbstractFill},
     x::AbstractArray{<:Real,N},
 ) where {N,M}
     return __logpdf(d, x)
 end
 function _logpdf(
-    d::ProductDistribution{2,M,<:Fill},
+    d::ProductDistribution{2,M,<:FillArrays.AbstractFill},
     x::AbstractMatrix{<:Real},
 ) where {M}
     return __logpdf(d, x)
 end
 function __logpdf(
-    d::ProductDistribution{N,M,<:Fill},
+    d::ProductDistribution{N,M,<:FillArrays.AbstractFill},
     x::AbstractArray{<:Real,N},
 ) where {N,M}
     return @inbounds loglikelihood(first(d.dists), x)

src/samplers/multinomial.jl

@@ -53,10 +53,12 @@ function MultinomialSampler(n::Int, prob::Vector{<:Real})
 end
 
 function rand(rng::AbstractRNG, s::MultinomialSampler)
-    return _rand!(rng, s, Vector{Int}(undef, length(s.prob)))
+    x = Vector{Int}(undef, length(s.prob))
+    return rand!(rng, s, x)
 end
-function _rand!(rng::AbstractRNG, s::MultinomialSampler,
+@inline function rand!(rng::AbstractRNG, s::MultinomialSampler,
                 x::AbstractVector{<:Real})
+    @checkbounds length(s) == length(x)
     n = s.n
     k = length(s)
     if n^2 > k