From e8a9fa07cd318e2b4600031f6dd6562168b42249 Mon Sep 17 00:00:00 2001 From: Kyle Daruwalla Date: Mon, 9 Dec 2019 08:10:19 -0600 Subject: [PATCH 01/11] Allow probability vectors other than Float64 vector for Multinomial distributions. Tests still pass. --- src/samplers/multinomial.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/samplers/multinomial.jl b/src/samplers/multinomial.jl index e8464a1b4f..89c0cef354 100644 --- a/src/samplers/multinomial.jl +++ b/src/samplers/multinomial.jl @@ -1,5 +1,5 @@ -function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{Float64}, - x::AbstractVector{<:Real}) +function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{PT}, + x::AbstractVector{T}) where {T<:Real, PT<:Real} k = length(p) length(x) == k || throw(DimensionMismatch("Invalid argument dimension.")) From f5f8d5f024383269dc735d382ebee43a3b2e3e81 Mon Sep 17 00:00:00 2001 From: Kyle Daruwalla Date: Mon, 9 Dec 2019 10:50:39 -0600 Subject: [PATCH 02/11] Remove explicit concrete types from multinom_rand --- src/samplers/multinomial.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/samplers/multinomial.jl b/src/samplers/multinomial.jl index 89c0cef354..73e2287516 100644 --- a/src/samplers/multinomial.jl +++ b/src/samplers/multinomial.jl @@ -1,5 +1,5 @@ -function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{PT}, - x::AbstractVector{T}) where {T<:Real, PT<:Real} +function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{<:Real}, + x::AbstractVector{<:Real}) k = length(p) length(x) == k || throw(DimensionMismatch("Invalid argument dimension.")) From 654a293374e8e0cfd7e11e5c1043690c77700255 Mon Sep 17 00:00:00 2001 From: Kyle Daruwalla Date: Wed, 29 Jan 2020 09:47:42 -0600 Subject: [PATCH 03/11] Added test to make sure rand works for real types other than Float64 --- test/multivariate/multinomial.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/multivariate/multinomial.jl b/test/multivariate/multinomial.jl index d944835067..9e0a1ef707 100644 --- a/test/multivariate/multinomial.jl +++ b/test/multivariate/multinomial.jl @@ -1,6 +1,6 @@ # Tests for Multinomial -using Distributions, Random, StaticArrays +using Distributions, Random, StaticArrays, ForwardDiff using Test @@ -210,3 +210,9 @@ p_v = [0.1, 0.4, 0.3, 0.8] @test_throws DomainError Multinomial(10, p_v) @test_throws DomainError Multinomial(10, p_v; check_args=true) Multinomial(10, p_v; check_args=false) # should not warn + +# check different prob vector types +p = [0.2, 0.4, 0.3, 0.1] +@test (rand(Multinomial(10, p)); true) +@test (rand(Multinomial(10, convert.(Float32, p))); true) +# @test (rand(Multinomial(10, convert.(ForwardDiff.Dual{Float64}, p))); true) From ce0895ee8f2fc8ddaa6f3f8092446093f7028250 Mon Sep 17 00:00:00 2001 From: Kyle Daruwalla Date: Mon, 3 Feb 2020 10:56:24 -0600 Subject: [PATCH 04/11] Removed extra test --- test/multivariate/multinomial.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/multivariate/multinomial.jl b/test/multivariate/multinomial.jl index 9e0a1ef707..488ade18da 100644 --- a/test/multivariate/multinomial.jl +++ b/test/multivariate/multinomial.jl @@ -215,4 +215,3 @@ Multinomial(10, p_v; check_args=false) # should not warn p = [0.2, 0.4, 0.3, 0.1] @test (rand(Multinomial(10, p)); true) @test (rand(Multinomial(10, convert.(Float32, p))); true) -# @test (rand(Multinomial(10, convert.(ForwardDiff.Dual{Float64}, p))); true) From d7a2284795d693786f5e588b85715aa8c4a0c752 Mon Sep 17 00:00:00 2001 From: Lior Blech Date: Mon, 19 Jun 2023 19:08:54 +0300 Subject: [PATCH 05/11] updated tests for type stability --- test/multivariate/multinomial.jl | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/test/multivariate/multinomial.jl b/test/multivariate/multinomial.jl index 488ade18da..9e91b32208 100644 --- a/test/multivariate/multinomial.jl +++ b/test/multivariate/multinomial.jl @@ -1,7 +1,8 @@ # Tests for Multinomial -using Distributions, Random, StaticArrays, ForwardDiff +using Distributions, Random, StaticArrays using Test +import Distributions: multinom_rand! p = [0.2, 0.5, 0.3] @@ -215,3 +216,13 @@ Multinomial(10, p_v; check_args=false) # should not warn p = [0.2, 0.4, 0.3, 0.1] @test (rand(Multinomial(10, p)); true) @test (rand(Multinomial(10, convert.(Float32, p))); true) + +# check type stability +@inferred rand(Multinomial(10, convert.(Float32, p))) +@inferred rand(Multinomial(10, convert.(Float64, p))) +@inferred rand(Multinomial(10, convert.(Float16, p))) + +x = zeros(Float64, size(p)...) +@inferred multinom_rand!(rng, 10, convert.(Float32, p), convert.(Float32, x)) +@inferred multinom_rand!(rng, 10, convert.(Float64, p), convert.(Float64, x)) +@inferred multinom_rand!(rng, 10, convert.(Float16, p), convert.(Float16, x)) From 0721c7540b20385dadd6efe45009a0ea24f6216f Mon Sep 17 00:00:00 2001 From: Lior Blech Date: Mon, 19 Jun 2023 19:50:14 +0300 Subject: [PATCH 06/11] same test but not using internal function --- test/multivariate/multinomial.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/multivariate/multinomial.jl b/test/multivariate/multinomial.jl index 9e91b32208..dd410b0fc3 100644 --- a/test/multivariate/multinomial.jl +++ b/test/multivariate/multinomial.jl @@ -2,7 +2,6 @@ using Distributions, Random, StaticArrays using Test -import Distributions: multinom_rand! p = [0.2, 0.5, 0.3] @@ -223,6 +222,6 @@ p = [0.2, 0.4, 0.3, 0.1] @inferred rand(Multinomial(10, convert.(Float16, p))) x = zeros(Float64, size(p)...) -@inferred multinom_rand!(rng, 10, convert.(Float32, p), convert.(Float32, x)) -@inferred multinom_rand!(rng, 10, convert.(Float64, p), convert.(Float64, x)) -@inferred multinom_rand!(rng, 10, convert.(Float16, p), convert.(Float16, x)) +@inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float32, x)) +@inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float64, x)) +@inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float16, x)) From eed61f757c3e895533c7748d414fb9ebac5a3b61 Mon Sep 17 00:00:00 2001 From: Lior Blech Date: Wed, 21 Jun 2023 00:59:33 +0300 Subject: [PATCH 07/11] repeat tests under different types instead --- test/multivariate/multinomial.jl | 415 ++++++++++++++++--------------- 1 file changed, 209 insertions(+), 206 deletions(-) diff --git a/test/multivariate/multinomial.jl b/test/multivariate/multinomial.jl index dd410b0fc3..9c5e19d792 100644 --- a/test/multivariate/multinomial.jl +++ b/test/multivariate/multinomial.jl @@ -4,224 +4,227 @@ using Distributions, Random, StaticArrays using Test -p = [0.2, 0.5, 0.3] -nt = 10 -d = Multinomial(nt, p) -rng = MersenneTwister(123) - -@testset "Testing Multinomial with $key" for (key, func) in - Dict("rand(...)" => [rand, rand], - "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)]) - -# Basics - -@test length(d) == 3 -@test d.n == nt -@test mean(d) ≈ [2., 5., 3.] -@test var(d) ≈ [1.6, 2.5, 2.1] -@test cov(d) ≈ [1.6 -1.0 -0.6; -1.0 2.5 -1.5; -0.6 -1.5 2.1] - -@test insupport(d, [1, 6, 3]) -@test !insupport(d, [2, 6, 3]) -@test partype(d) == Float64 -@test partype(Multinomial(nt, Vector{Float32}(p))) == Float32 - -# Conversion -@test typeof(d) == Multinomial{Float64, Vector{Float64}} -@test typeof(Multinomial(nt, Vector{Float32}(p))) == Multinomial{Float32, Vector{Float32}} -@test typeof(convert(Multinomial{Float32}, d)) == Multinomial{Float32, Vector{Float32}} -@test typeof(convert(Multinomial{Float32, Vector{Float32}}, d)) == Multinomial{Float32, Vector{Float32}} -@test typeof(convert(Multinomial{Float32}, params(d)...)) == Multinomial{Float32, Vector{Float32}} -@test typeof(convert(Multinomial{Float32, Vector{Float32}}, params(d)...)) == Multinomial{Float32, Vector{Float32}} -@test convert(Multinomial{Float64}, d) === d -@test convert(Multinomial{Float64, Vector{Float64}}, d) === d - -# random sampling - -x = func[1](d) -@test isa(x, Vector{Int}) -@test sum(x) == nt -@test insupport(d, x) -@test size(x) == size(d) -@test length(x) == length(d) -@test d == typeof(d)(params(d)...) -@test d == deepcopy(d) - -x = func[2](d, 100) -@test isa(x, Matrix{Int}) -@test all(sum(x, dims=1) .== nt) -@test all(insupport(d, x)) - -x = func[1](sampler(d)) -@test isa(x, Vector{Int}) -@test sum(x) == nt -@test insupport(d, x) - -# logpdf - -x1 = [1, 6, 3] - -@test isapprox(pdf(d, x1), 0.070875, atol=1.0e-8) -@test logpdf(d, x1) ≈ log(pdf(d, x1)) - -x = func[2](d, 100) -pv = pdf(d, x) -lp = logpdf(d, x) -for i in 1 : size(x, 2) - @test pv[i] ≈ pdf(d, x[:,i]) - @test lp[i] ≈ logpdf(d, x[:,i]) -end - -# test type stability of logpdf -@test typeof(logpdf(convert(Multinomial{Float32}, d), x1)) == Float32 - -# test degenerate cases of logpdf -d1 = Multinomial(1, [0.5, 0.5, 0.0]) -d2 = Multinomial(0, [0.5, 0.5, 0.0]) -x2 = [1, 0, 0] -x3 = [0, 0, 1] -x4 = [1, 0, 1] - -@test logpdf(d1, x2) ≈ log(0.5) -@test logpdf(d2, x2) == -Inf -@test logpdf(d1, x3) == -Inf -@test logpdf(d2, x3) == -Inf - -# suffstats - -d0 = d -n0 = 100 -x = func[2](d0, n0) -w = func[1](n0) - -ss = suffstats(Multinomial, x) -@test isa(ss, Distributions.MultinomialStats) -@test ss.n == nt -@test ss.scnts == vec(sum(Float64[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)], dims=2)) -@test ss.tw == n0 - -ss = suffstats(Multinomial, x, w) -@test isa(ss, Distributions.MultinomialStats) -@test ss.n == nt -@test ss.scnts ≈ Float64[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)] * w -@test ss.tw ≈ sum(w) - -# fit - -x = func[2](d0, 10^5) -@test size(x) == (length(d0), 10^5) -@test all(sum(x, dims=1) .== nt) - -r = fit(Multinomial, x) -@test r.n == nt -@test length(r) == length(p) -@test isapprox(probs(r), p, atol=0.02) -r = fit(Multinomial{Float64}, x) -@test r.n == nt -@test length(r) == length(p) -@test isapprox(probs(r), p, atol=0.02) - -r = fit_mle(Multinomial, x, fill(2.0, size(x,2))) -@test r.n == nt -@test length(r) == length(p) -@test isapprox(probs(r), p, atol=0.02) - -# behavior for n = 0 -d0 = Multinomial(0, p) -@test func[1](d0) == [0, 0, 0] -@test pdf(d0, [0, 0, 0]) == 1 -@test pdf(d0, [0, 1, 0]) == 0 -@test mean(d0) == [0, 0, 0] -@test var(d0) == [0, 0, 0] -@test cov(d0) == zeros(3, 3) -@test entropy(d0) == 0 -@test insupport(d0, [0, 0, 0]) == true -@test insupport(d0, [0, 0, 4]) == false -@test length(d0) == 3 -@test size(d0) == (3,) - -# Abstract vector p - -@test typeof(Multinomial(nt, SVector{length(p), Float64}(p))) == Multinomial{Float64, SVector{3, Float64}} - -end - -@testset "Testing Multinomial with $key" for (key, func) in - Dict("rand!(...)" => (dist, X) -> rand!(dist, X), - "rand!(rng, ...)" => (dist, X) -> rand!(rng, dist, X)) +@testset "testing Multinomial with type $type" for (type, tol) in zip((Float16, Float32, Float64), (1e-3, 1e-6, 1e-8)) + p = [0.2, 0.5, 0.3] + nt = 10 + rng = MersenneTwister(123) + d = Multinomial(nt, convert.(type, p)) + + @testset "Testing Multinomial with $key" for (key, func) in + Dict("rand(...)" => [rand, rand], + "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)]) + + # Basics + + @test length(d) == 3 + @test d.n == nt + @test mean(d) ≈ [2., 5., 3.] + @test var(d) ≈ [1.6, 2.5, 2.1] + @test cov(d) ≈ [1.6 -1.0 -0.6; -1.0 2.5 -1.5; -0.6 -1.5 2.1] + + @test insupport(d, [1, 6, 3]) + @test !insupport(d, [2, 6, 3]) + @test partype(d) == type #Float64 + @test partype(Multinomial(nt, Vector{Float32}(p))) == Float32 + + # Conversion + @test typeof(d) == Multinomial{type, Vector{type}} + @test typeof(Multinomial(nt, Vector{Float32}(p))) == Multinomial{Float32, Vector{Float32}} + @test typeof(convert(Multinomial{Float32}, d)) == Multinomial{Float32, Vector{Float32}} + @test typeof(convert(Multinomial{Float32, Vector{Float32}}, d)) == Multinomial{Float32, Vector{Float32}} + @test typeof(convert(Multinomial{Float32}, params(d)...)) == Multinomial{Float32, Vector{Float32}} + @test typeof(convert(Multinomial{Float32, Vector{Float32}}, params(d)...)) == Multinomial{Float32, Vector{Float32}} + @test convert(Multinomial{type}, d) === d + @test convert(Multinomial{type, Vector{type}}, d) === d + # random sampling - X = Matrix{Int}(undef, length(p), 100) - x = func(d, X) - @test x ≡ X + + x = func[1](d) + @test isa(x, Vector{Int}) + @test sum(x) == nt + @test insupport(d, x) + @test size(x) == size(d) + @test length(x) == length(d) + @test d == typeof(d)(params(d)...) + @test d == deepcopy(d) + + x = func[2](d, 100) @test isa(x, Matrix{Int}) @test all(sum(x, dims=1) .== nt) @test all(insupport(d, x)) + + x = func[1](sampler(d)) + @test isa(x, Vector{Int}) + @test sum(x) == nt + @test insupport(d, x) + + # logpdf + + x1 = [1, 6, 3] + + @test isapprox(pdf(d, x1), 0.070875, atol=tol) + @test logpdf(d, x1) ≈ log(pdf(d, x1)) + + x = func[2](d, 100) pv = pdf(d, x) lp = logpdf(d, x) for i in 1 : size(x, 2) @test pv[i] ≈ pdf(d, x[:,i]) @test lp[i] ≈ logpdf(d, x[:,i]) end + + # test type stability of logpdf + @test typeof(logpdf(convert(Multinomial{Float32}, d), x1)) == Float32 + + # test degenerate cases of logpdf + d1 = Multinomial(1, [0.5, 0.5, 0.0]) + d2 = Multinomial(0, [0.5, 0.5, 0.0]) + x2 = [1, 0, 0] + x3 = [0, 0, 1] + x4 = [1, 0, 1] + + @test logpdf(d1, x2) ≈ log(0.5) + @test logpdf(d2, x2) == -Inf + @test logpdf(d1, x3) == -Inf + @test logpdf(d2, x3) == -Inf + + # suffstats + + d0 = d + n0 = 100 + x = func[2](d0, n0) + w = func[1](n0) + + ss = suffstats(Multinomial, x) + @test isa(ss, Distributions.MultinomialStats) + @test ss.n == nt + @test ss.scnts == vec(sum(Float64[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)], dims=2)) + @test ss.tw == n0 + + ss = suffstats(Multinomial, x, w) + @test isa(ss, Distributions.MultinomialStats) + @test ss.n == nt + @test ss.scnts ≈ Float64[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)] * w + @test ss.tw ≈ sum(w) + + # fit + + x = func[2](d0, 10^5) + @test size(x) == (length(d0), 10^5) + @test all(sum(x, dims=1) .== nt) + + r = fit(Multinomial, x) + @test r.n == nt + @test length(r) == length(p) + @test isapprox(probs(r), p, atol=0.02) + r = fit(Multinomial{Float64}, x) + @test r.n == nt + @test length(r) == length(p) + @test isapprox(probs(r), p, atol=0.02) + + r = fit_mle(Multinomial, x, fill(2.0, size(x,2))) + @test r.n == nt + @test length(r) == length(p) + @test isapprox(probs(r), p, atol=0.02) + + # behavior for n = 0 + d0 = Multinomial(0, p) + @test func[1](d0) == [0, 0, 0] + @test pdf(d0, [0, 0, 0]) == 1 + @test pdf(d0, [0, 1, 0]) == 0 + @test mean(d0) == [0, 0, 0] + @test var(d0) == [0, 0, 0] + @test cov(d0) == zeros(3, 3) + @test entropy(d0) == 0 + @test insupport(d0, [0, 0, 0]) == true + @test insupport(d0, [0, 0, 4]) == false + @test length(d0) == 3 + @test size(d0) == (3,) + + # Abstract vector p + + @test typeof(Multinomial(nt, SVector{length(p), Float64}(p))) == Multinomial{Float64, SVector{3, Float64}} + + end + + @testset "Testing Multinomial with $key" for (key, func) in + Dict("rand!(...)" => (dist, X) -> rand!(dist, X), + "rand!(rng, ...)" => (dist, X) -> rand!(rng, dist, X)) + # random sampling + X = Matrix{Int}(undef, length(p), 100) + x = func(d, X) + @test x ≡ X + @test isa(x, Matrix{Int}) + @test all(sum(x, dims=1) .== nt) + @test all(insupport(d, x)) + pv = pdf(d, x) + lp = logpdf(d, x) + for i in 1 : size(x, 2) + @test pv[i] ≈ pdf(d, x[:,i]) + @test lp[i] ≈ logpdf(d, x[:,i]) + end + end + + @testset "Testing Multinomial with $key" for (key, func) in + Dict("rand!(..., true)" => (dist, X) -> rand!(dist, X, true), + "rand!(rng, ..., true)" => (dist, X) -> rand!(rng, dist, X, true)) + # random sampling + X = Vector{Vector{Int}}(undef, 100) + x = func(d, X) + @test x ≡ X + @test all(sum.(x) .== nt) + @test all(insupport(d, a) for a in x) + end + + @testset "Testing Multinomial with $key" for (key, func) in + Dict("rand!(..., false)" => (dist, X) -> rand!(dist, X, false), + "rand!(rng, ..., false)" => (dist, X) -> rand!(rng, dist, X, false)) + # random sampling + X = [Vector{Int}(undef, length(p)) for _ in Base.OneTo(100)] + x1 = X[1] + x = func(d, X) + @test x1 ≡ X[1] + @test all(sum.(x) .== nt) + @test all(insupport(d, a) for a in x) + end + + repeats = 10 + m = Vector{Vector{partype(d)}}(undef, repeats) + rand!(d, m) + @test isassigned(m, 1) + m1=m[1] + rand!(d, m) + @test m1 ≡ m[1] + rand!(d, m, true) + @test m1 ≢ m[1] + m1 = m[1] + rand!(d, m, false) + @test m1 ≡ m[1] + + p = [0.2, 0.4, 0.3, 0.1] + nt = 10 + d = Multinomial(nt, p) + @test_throws DimensionMismatch rand!(d, m, false) + @test_nowarn rand!(d, m) + + p_v = [0.1, 0.4, 0.3, 0.8] + @test_throws DomainError Multinomial(10, p_v) + @test_throws DomainError Multinomial(10, p_v; check_args=true) + Multinomial(10, p_v; check_args=false) # should not warn end + +@testset "Type stability" begin -@testset "Testing Multinomial with $key" for (key, func) in - Dict("rand!(..., true)" => (dist, X) -> rand!(dist, X, true), - "rand!(rng, ..., true)" => (dist, X) -> rand!(rng, dist, X, true)) - # random sampling - X = Vector{Vector{Int}}(undef, 100) - x = func(d, X) - @test x ≡ X - @test all(sum.(x) .== nt) - @test all(insupport(d, a) for a in x) -end - -@testset "Testing Multinomial with $key" for (key, func) in - Dict("rand!(..., false)" => (dist, X) -> rand!(dist, X, false), - "rand!(rng, ..., false)" => (dist, X) -> rand!(rng, dist, X, false)) - # random sampling - X = [Vector{Int}(undef, length(p)) for _ in Base.OneTo(100)] - x1 = X[1] - x = func(d, X) - @test x1 ≡ X[1] - @test all(sum.(x) .== nt) - @test all(insupport(d, a) for a in x) -end - -repeats = 10 -m = Vector{Vector{partype(d)}}(undef, repeats) -rand!(d, m) -@test isassigned(m, 1) -m1=m[1] -rand!(d, m) -@test m1 ≡ m[1] -rand!(d, m, true) -@test m1 ≢ m[1] -m1 = m[1] -rand!(d, m, false) -@test m1 ≡ m[1] - -p = [0.2, 0.4, 0.3, 0.1] -nt = 10 -d = Multinomial(nt, p) -@test_throws DimensionMismatch rand!(d, m, false) -@test_nowarn rand!(d, m) + p = [0.2, 0.5, 0.3] -p_v = [0.1, 0.4, 0.3, 0.8] -@test_throws DomainError Multinomial(10, p_v) -@test_throws DomainError Multinomial(10, p_v; check_args=true) -Multinomial(10, p_v; check_args=false) # should not warn + @inferred rand(Multinomial(10, convert.(Float32, p))) + @inferred rand(Multinomial(10, convert.(Float64, p))) + @inferred rand(Multinomial(10, convert.(Float16, p))) -# check different prob vector types -p = [0.2, 0.4, 0.3, 0.1] -@test (rand(Multinomial(10, p)); true) -@test (rand(Multinomial(10, convert.(Float32, p))); true) - -# check type stability -@inferred rand(Multinomial(10, convert.(Float32, p))) -@inferred rand(Multinomial(10, convert.(Float64, p))) -@inferred rand(Multinomial(10, convert.(Float16, p))) + x = zeros(Float64, size(p)...) + @inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float32, x)) + @inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float64, x)) + @inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float16, x)) +end -x = zeros(Float64, size(p)...) -@inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float32, x)) -@inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float64, x)) -@inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float16, x)) + \ No newline at end of file From 0eabf8e45885e3bb3ebe9624804237fc24c4b8fd Mon Sep 17 00:00:00 2001 From: Lior Blech Date: Wed, 21 Jun 2023 01:22:42 +0300 Subject: [PATCH 08/11] type stability within multinomial. may want to extend into Binomial --- src/samplers/multinomial.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/samplers/multinomial.jl b/src/samplers/multinomial.jl index 73e2287516..b478915566 100644 --- a/src/samplers/multinomial.jl +++ b/src/samplers/multinomial.jl @@ -1,9 +1,9 @@ -function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{<:Real}, - x::AbstractVector{<:Real}) +function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{Tp}, + x::AbstractVector{Tx}) where {Tx <: Real, Tp <: Real} k = length(p) length(x) == k || throw(DimensionMismatch("Invalid argument dimension.")) - rp = 1.0 # remaining total probability + rp = Tp(1.0) # remaining total probability i = 0 km1 = k - 1 @@ -11,23 +11,23 @@ function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{<:Real}, i += 1 @inbounds pi = p[i] if pi < rp - xi = rand(rng, Binomial(n, pi / rp)) + xi = Tx(rand(rng, Binomial(n, Float64(pi / rp)))) @inbounds x[i] = xi - n -= xi + n -= Int64(xi) rp -= pi else # In this case, we don't even have to sample # from Binomial. Just assign remaining counts # to xi. - @inbounds x[i] = n + @inbounds x[i] = Tx(n) n = 0 # rp = 0.0 (no need for this, as rp is no longer needed) end end if i == km1 - @inbounds x[k] = n + @inbounds x[k] = Tx(n) else # n must have been zero z = zero(eltype(x)) for j = i+1 : k From 19aa408df56365d84d8a685e221188b4b42e6f49 Mon Sep 17 00:00:00 2001 From: Lior Blech Date: Wed, 21 Jun 2023 15:49:43 +0300 Subject: [PATCH 09/11] Apply suggestions from code review Type conversions implicitly at runtime covers more systems better Co-authored-by: David Widmann --- src/samplers/multinomial.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/samplers/multinomial.jl b/src/samplers/multinomial.jl index b478915566..864064ec93 100644 --- a/src/samplers/multinomial.jl +++ b/src/samplers/multinomial.jl @@ -1,9 +1,10 @@ -function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{Tp}, - x::AbstractVector{Tx}) where {Tx <: Real, Tp <: Real} +function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{<:Real}, + x::AbstractVector{<:Real}) k = length(p) length(x) == k || throw(DimensionMismatch("Invalid argument dimension.")) - rp = Tp(1.0) # remaining total probability + z = zero(eltype(p)) + rp = oftype(z + z, 1) # remaining total probability (widens type if needed) i = 0 km1 = k - 1 @@ -13,21 +14,21 @@ function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{Tp}, if pi < rp xi = Tx(rand(rng, Binomial(n, Float64(pi / rp)))) @inbounds x[i] = xi - n -= Int64(xi) + n -= xi rp -= pi else # In this case, we don't even have to sample # from Binomial. Just assign remaining counts # to xi. - @inbounds x[i] = Tx(n) + @inbounds x[i] = n n = 0 # rp = 0.0 (no need for this, as rp is no longer needed) end end if i == km1 - @inbounds x[k] = Tx(n) + @inbounds x[k] = n else # n must have been zero z = zero(eltype(x)) for j = i+1 : k From 1529164cc24684046134bd966d32fec42adba051 Mon Sep 17 00:00:00 2001 From: Lior Blech Date: Wed, 21 Jun 2023 15:52:26 +0300 Subject: [PATCH 10/11] delete leftover type variable --- src/samplers/multinomial.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/samplers/multinomial.jl b/src/samplers/multinomial.jl index 864064ec93..7ce8b89730 100644 --- a/src/samplers/multinomial.jl +++ b/src/samplers/multinomial.jl @@ -12,7 +12,7 @@ function multinom_rand!(rng::AbstractRNG, n::Int, p::AbstractVector{<:Real}, i += 1 @inbounds pi = p[i] if pi < rp - xi = Tx(rand(rng, Binomial(n, Float64(pi / rp)))) + xi = rand(rng, Binomial(n, Float64(pi / rp))) @inbounds x[i] = xi n -= xi rp -= pi From 23d69d44c21e0687af31ca8bdae93a0ca8fb352e Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 30 Jun 2023 16:54:20 +0200 Subject: [PATCH 11/11] Simplify tests --- test/multivariate/multinomial.jl | 101 +++++++++++++------------------ 1 file changed, 43 insertions(+), 58 deletions(-) diff --git a/test/multivariate/multinomial.jl b/test/multivariate/multinomial.jl index 9c5e19d792..6a9ebcbb19 100644 --- a/test/multivariate/multinomial.jl +++ b/test/multivariate/multinomial.jl @@ -4,11 +4,11 @@ using Distributions, Random, StaticArrays using Test -@testset "testing Multinomial with type $type" for (type, tol) in zip((Float16, Float32, Float64), (1e-3, 1e-6, 1e-8)) - p = [0.2, 0.5, 0.3] +@testset "testing Multinomial with type $T" for T in (Float16, Float32, Float64) + p = T[0.2, 0.5, 0.3] nt = 10 rng = MersenneTwister(123) - d = Multinomial(nt, convert.(type, p)) + d = Multinomial(nt, p) @testset "Testing Multinomial with $key" for (key, func) in Dict("rand(...)" => [rand, rand], @@ -18,28 +18,29 @@ using Test @test length(d) == 3 @test d.n == nt - @test mean(d) ≈ [2., 5., 3.] - @test var(d) ≈ [1.6, 2.5, 2.1] - @test cov(d) ≈ [1.6 -1.0 -0.6; -1.0 2.5 -1.5; -0.6 -1.5 2.1] + @test mean(d) ≈ T[2., 5., 3.] + @test var(d) ≈ T[1.6, 2.5, 2.1] + @test cov(d) ≈ T[1.6 -1.0 -0.6; -1.0 2.5 -1.5; -0.6 -1.5 2.1] @test insupport(d, [1, 6, 3]) @test !insupport(d, [2, 6, 3]) - @test partype(d) == type #Float64 - @test partype(Multinomial(nt, Vector{Float32}(p))) == Float32 + @test partype(d) === T # Conversion - @test typeof(d) == Multinomial{type, Vector{type}} - @test typeof(Multinomial(nt, Vector{Float32}(p))) == Multinomial{Float32, Vector{Float32}} - @test typeof(convert(Multinomial{Float32}, d)) == Multinomial{Float32, Vector{Float32}} - @test typeof(convert(Multinomial{Float32, Vector{Float32}}, d)) == Multinomial{Float32, Vector{Float32}} - @test typeof(convert(Multinomial{Float32}, params(d)...)) == Multinomial{Float32, Vector{Float32}} - @test typeof(convert(Multinomial{Float32, Vector{Float32}}, params(d)...)) == Multinomial{Float32, Vector{Float32}} - @test convert(Multinomial{type}, d) === d - @test convert(Multinomial{type, Vector{type}}, d) === d + @test typeof(d) === Multinomial{T, Vector{T}} + for S in (Float16, Float32, Float64) + S === T && continue + @test typeof(convert(Multinomial{S}, d)) === Multinomial{S, Vector{S}} + @test typeof(convert(Multinomial{S, Vector{S}}, d)) === Multinomial{S, Vector{S}} + @test typeof(convert(Multinomial{S}, params(d)...)) === Multinomial{S, Vector{S}} + @test typeof(convert(Multinomial{S, Vector{S}}, params(d)...)) == Multinomial{S, Vector{S}} + end + @test convert(Multinomial{T}, d) === d + @test convert(Multinomial{T, Vector{T}}, d) === d # random sampling - x = func[1](d) + x = @inferred(func[1](d)) @test isa(x, Vector{Int}) @test sum(x) == nt @test insupport(d, x) @@ -48,12 +49,12 @@ using Test @test d == typeof(d)(params(d)...) @test d == deepcopy(d) - x = func[2](d, 100) + x = @inferred(func[2](d, 100)) @test isa(x, Matrix{Int}) @test all(sum(x, dims=1) .== nt) @test all(insupport(d, x)) - x = func[1](sampler(d)) + x = @inferred(func[1](sampler(d))) @test isa(x, Vector{Int}) @test sum(x) == nt @test insupport(d, x) @@ -62,8 +63,8 @@ using Test x1 = [1, 6, 3] - @test isapprox(pdf(d, x1), 0.070875, atol=tol) - @test logpdf(d, x1) ≈ log(pdf(d, x1)) + @test @inferred(pdf(d, x1)) ≈ T(0.070875) + @test @inferred(logpdf(d, x1)) ≈ log(pdf(d, x1)) x = func[2](d, 100) pv = pdf(d, x) @@ -73,17 +74,17 @@ using Test @test lp[i] ≈ logpdf(d, x[:,i]) end - # test type stability of logpdf - @test typeof(logpdf(convert(Multinomial{Float32}, d), x1)) == Float32 + # test result type of logpdf + @test typeof(logpdf(d, x1)) === T # test degenerate cases of logpdf - d1 = Multinomial(1, [0.5, 0.5, 0.0]) - d2 = Multinomial(0, [0.5, 0.5, 0.0]) + d1 = Multinomial(1, T[0.5, 0.5, 0.0]) + d2 = Multinomial(0, T[0.5, 0.5, 0.0]) x2 = [1, 0, 0] x3 = [0, 0, 1] x4 = [1, 0, 1] - @test logpdf(d1, x2) ≈ log(0.5) + @test logpdf(d1, x2) ≈ T(log(0.5)) @test logpdf(d2, x2) == -Inf @test logpdf(d1, x3) == -Inf @test logpdf(d2, x3) == -Inf @@ -98,13 +99,13 @@ using Test ss = suffstats(Multinomial, x) @test isa(ss, Distributions.MultinomialStats) @test ss.n == nt - @test ss.scnts == vec(sum(Float64[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)], dims=2)) + @test ss.scnts == vec(sum(T[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)], dims=2)) @test ss.tw == n0 ss = suffstats(Multinomial, x, w) @test isa(ss, Distributions.MultinomialStats) @test ss.n == nt - @test ss.scnts ≈ Float64[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)] * w + @test ss.scnts ≈ T[x[i,j] for i = 1:size(x,1), j = 1:size(x,2)] * w @test ss.tw ≈ sum(w) # fit @@ -116,16 +117,16 @@ using Test r = fit(Multinomial, x) @test r.n == nt @test length(r) == length(p) - @test isapprox(probs(r), p, atol=0.02) + @test probs(r) ≈ p atol=0.02 r = fit(Multinomial{Float64}, x) @test r.n == nt @test length(r) == length(p) - @test isapprox(probs(r), p, atol=0.02) + @test probs(r) ≈ p atol=0.02 r = fit_mle(Multinomial, x, fill(2.0, size(x,2))) @test r.n == nt @test length(r) == length(p) - @test isapprox(probs(r), p, atol=0.02) + @test probs(r) ≈ p atol=0.02 # behavior for n = 0 d0 = Multinomial(0, p) @@ -143,7 +144,7 @@ using Test # Abstract vector p - @test typeof(Multinomial(nt, SVector{length(p), Float64}(p))) == Multinomial{Float64, SVector{3, Float64}} + @test typeof(Multinomial(nt, SVector{length(p), T}(p))) == Multinomial{T, SVector{3, T}} end @@ -152,7 +153,7 @@ using Test "rand!(rng, ...)" => (dist, X) -> rand!(rng, dist, X)) # random sampling X = Matrix{Int}(undef, length(p), 100) - x = func(d, X) + x = @inferred(func(d, X)) @test x ≡ X @test isa(x, Matrix{Int}) @test all(sum(x, dims=1) .== nt) @@ -170,7 +171,7 @@ using Test "rand!(rng, ..., true)" => (dist, X) -> rand!(rng, dist, X, true)) # random sampling X = Vector{Vector{Int}}(undef, 100) - x = func(d, X) + x = @inferred(func(d, X)) @test x ≡ X @test all(sum.(x) .== nt) @test all(insupport(d, a) for a in x) @@ -182,49 +183,33 @@ using Test # random sampling X = [Vector{Int}(undef, length(p)) for _ in Base.OneTo(100)] x1 = X[1] - x = func(d, X) + x = @inferred(func(d, X)) @test x1 ≡ X[1] @test all(sum.(x) .== nt) @test all(insupport(d, a) for a in x) end repeats = 10 - m = Vector{Vector{partype(d)}}(undef, repeats) - rand!(d, m) + m = Vector{Vector{T}}(undef, repeats) + @inferred(rand!(d, m)) @test isassigned(m, 1) m1=m[1] - rand!(d, m) + @inferred(rand!(d, m)) @test m1 ≡ m[1] - rand!(d, m, true) + @inferred(rand!(d, m, true)) @test m1 ≢ m[1] m1 = m[1] - rand!(d, m, false) + @inferred(rand!(d, m, false)) @test m1 ≡ m[1] - p = [0.2, 0.4, 0.3, 0.1] + p = T[0.2, 0.4, 0.3, 0.1] nt = 10 d = Multinomial(nt, p) @test_throws DimensionMismatch rand!(d, m, false) @test_nowarn rand!(d, m) - p_v = [0.1, 0.4, 0.3, 0.8] + p_v = T[0.1, 0.4, 0.3, 0.8] @test_throws DomainError Multinomial(10, p_v) @test_throws DomainError Multinomial(10, p_v; check_args=true) Multinomial(10, p_v; check_args=false) # should not warn end - -@testset "Type stability" begin - - p = [0.2, 0.5, 0.3] - - @inferred rand(Multinomial(10, convert.(Float32, p))) - @inferred rand(Multinomial(10, convert.(Float64, p))) - @inferred rand(Multinomial(10, convert.(Float16, p))) - - x = zeros(Float64, size(p)...) - @inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float32, x)) - @inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float64, x)) - @inferred rand!(Multinomial(10, convert.(Float32, p)), convert.(Float16, x)) -end - - \ No newline at end of file