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})
+ @boundscheck 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})
+ @boundscheck length(s) == length(x)
  n = s.n
  k = length(s)
  if n^2 > k