Dev

.github/workflows/CI.yml

@@ -25,6 +25,7 @@ jobs:
         version:
           - '1.6'
           - '1.7'
+          - '1.8'
           - 'nightly'
         os:
           - ubuntu-latest

Project.toml

@@ -21,6 +21,7 @@ NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 PrettyPrinting = "54e16d92-306c-5ea0-a30b-337be88ac337"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -42,6 +43,7 @@ MappedArrays = "0.4"
 NaNMath = "0.3, 1"
 PrettyPrinting = "0.3, 0.4"
 Reexport = "1"
+SpecialFunctions = "2"
 Static = "0.5, 0.6"
 Tricks = "0.1"
 julia = "1.3"

src/MeasureBase.jl

@@ -144,9 +144,11 @@ include("standard/stdmeasure.jl")
 include("standard/stduniform.jl")
 include("standard/stdexponential.jl")
 include("standard/stdlogistic.jl")
+include("standard/stdnormal.jl")
 include("latent-joint.jl")
 
 include("rand.jl")
+include("fixedrng.jl")
 
 include("density.jl")
 

src/combinators/half.jl

@@ -26,4 +26,4 @@ logdensity_def(μ::Half, x) = logdensity_def(unhalf(μ), x)
     insupport(unhalf(d), x)
 end
 
-testvalue(::Half) = 1.0
+testvalue(::Type{T}, ::Half) where {T} = one(T)

src/combinators/product.jl

@@ -161,7 +161,7 @@ marginals(μ::ProductMeasure) = μ.marginals
 _map(f, args...) = map(f, args...)
 _map(f, x::MappedArrays.ReadonlyMappedArray) = mappedarray(f ∘ x.f, x.data)
 
-testvalue(d::AbstractProductMeasure) = _map(testvalue, marginals(d))
+testvalue(::Type{T}, d::AbstractProductMeasure) where {T} = _map(m -> testvalue(T, m), marginals(d))
 
 export ⊗
 

src/combinators/spikemixture.jl

@@ -37,6 +37,6 @@ function Base.rand(rng::AbstractRNG, T::Type, μ::SpikeMixture)
     return (rand(rng, T) < μ.w) * rand(rng, T, μ.m)
 end
 
-testvalue(μ::SpikeMixture) = testvalue(μ.m)
+testvalue(::Type{T}, μ::SpikeMixture) where {T} = zero(T)
 
 insupport(μ::SpikeMixture, x) = dynamic(insupport(μ.m, x)) || iszero(x)

src/combinators/superpose.jl

@@ -34,7 +34,7 @@ function Pretty.tile(d::SuperpositionMeasure)
     result *= Pretty.literal(")")
 end
 
-testvalue(μ::SuperpositionMeasure) = testvalue(first(μ.components))
+testvalue(::Type{T}, μ::SuperpositionMeasure) where {T} = testvalue(T, first(μ.components))
 
 # SuperpositionMeasure(ms :: AbstractMeasure...) = SuperpositionMeasure{X,length(ms)}(ms)
 

src/combinators/transformedmeasure.jl

@@ -31,9 +31,26 @@ struct PushforwardMeasure{FF,IF,M,VC<:TransformVolCorr} <: AbstractPushforward
     volcorr::VC
 end
 
+
 gettransform(ν::PushforwardMeasure) = ν.f
 parent(ν::PushforwardMeasure) = ν.origin
 
+function transport_def(ν::PushforwardMeasure{FF,IF,M}, μ::M, x) where {FF,IF,M}
+    if μ == parent(ν)
+        return ν.f(x)
+    else
+        invoke(transport_def, Tuple{Any, PushforwardMeasure, Any}, ν, μ, x)
+    end
+end
+
+function transport_def(μ::M, ν::PushforwardMeasure{FF,IF,M}, y) where {FF,IF,M}
+    if μ == parent(ν)
+        return ν.inv_f(y)
+    else
+        invoke(transport_def, Tuple{Any, PushforwardMeasure, Any}, μ, ν, y)
+    end
+end
+
 function Pretty.tile(ν::PushforwardMeasure)
     Pretty.list_layout(Pretty.tile.([ν.f, ν.inv_f, ν.origin]); prefix = :PushforwardMeasure)
 end
@@ -67,7 +84,7 @@ end
 
 insupport(ν::PushforwardMeasure, y) = insupport(transport_origin(ν), to_origin(ν, y))
 
-testvalue(ν::PushforwardMeasure) = from_origin(ν, testvalue(transport_origin(ν)))
+testvalue(::Type{T}, ν::PushforwardMeasure) where {T} = from_origin(ν, testvalue(T, transport_origin(ν)))
 
 @inline function basemeasure(ν::PushforwardMeasure)
     PushforwardMeasure(ν.f, ν.inv_f, basemeasure(transport_origin(ν)), NoVolCorr())
@@ -91,15 +108,38 @@ end
 @inline to_origin(ν::PushforwardMeasure, y) = ν.inv_f(y)
 
 function Base.rand(rng::AbstractRNG, ::Type{T}, ν::PushforwardMeasure) where {T}
-    return from_origin(ν, rand(rng, T, transport_origin(ν)))
+    return ν.f(rand(rng, T, parent(ν)))
 end
 
 export pushfwd
 
 """
-    pushfwd(f, μ, volcorr = WithVolCorr())
+    pushfwd(f, [f_inverse,] μ, volcorr = WithVolCorr())
 
-Return the [pushforward measure](https://en.wikipedia.org/wiki/Pushforward_measure)
-from `μ` the [measurable function](https://en.wikipedia.org/wiki/Measurable_function) `f`.
+Return the [pushforward
+measure](https://en.wikipedia.org/wiki/Pushforward_measure) from `μ` the
+[measurable function](https://en.wikipedia.org/wiki/Measurable_function) `f`.
+
+If `f_inverse` is specified, it must be a valid inverse of the function given by
+restricting `f` to the support of `μ`.
 """
-pushfwd(f, μ, volcorr = WithVolCorr()) = PushforwardMeasure(f, inverse(f), μ, volcorr)
+pushfwd(f, μ, volcorr::TransformVolCorr = WithVolCorr()) = PushforwardMeasure(f, inverse(f), μ, volcorr)
+
+pushfwd(f, finv, μ, volcorr::TransformVolCorr = WithVolCorr()) = PushforwardMeasure(f, finv, μ, volcorr)
+
+@inline function pushfwd(f, μ::PushforwardMeasure, volcorr::TransformVolCorr = WithVolCorr())
+    _pushfwd(f, inverse(f), μ, volcorr, μ.volcorr)
+end
+
+
+@inline function pushfwd(f, finv, μ::PushforwardMeasure, volcorr::TransformVolCorr = WithVolCorr())
+    _pushfwd(f, finv, μ, volcorr, μ.volcorr)
+end
+
+@inline function _pushfwd(f, finv, μ::PushforwardMeasure, vf::V, vμ::V) where {V<:TransformVolCorr}
+    pushfwd(f ∘ μ.f, μ.inv_f ∘ finv, v)
+end
+
+@inline function _pushfwd(f, finv, μ::PushforwardMeasure, vf, vμ) 
+    PushforwardMeasure(f, finv, μ, vf)
+end

src/combinators/weighted.jl

@@ -20,6 +20,8 @@ function Base.rand(rng::AbstractRNG, ::Type{T}, μ::AbstractWeightedMeasure) whe
     rand(rng, T, basemeasure(μ))
 end
 
+testvalue(::Type{T}, μ::AbstractWeightedMeasure) where {T} = testvalue(T, basemeasure(μ))
+
 ###############################################################################
 
 struct WeightedMeasure{R,M} <: AbstractWeightedMeasure

src/domains.jl

@@ -32,9 +32,9 @@ Base.maximum(b::BoundedReals) = b.upper
 Base.show(io::IO, ::typeof(ℝ₊)) = print(io, "ℝ₊")
 Base.show(io::IO, ::typeof(𝕀)) = print(io, "𝕀")
 
-testvalue(::typeof(ℝ)) = 0.0
-testvalue(::typeof(ℝ₊)) = 1.0
-testvalue(::typeof(𝕀)) = 0.5
+testvalue(::Type{T}, ::typeof(ℝ)) where {T} = zero(T)
+testvalue(::Type{T}, ::typeof(ℝ₊)) where {T} = one(T)
+testvalue(::Type{T}, ::typeof(𝕀)) where {T} = one(T) / 2
 
 abstract type IntegerDomain <: AbstractDomain end
 

src/fixedrng.jl

@@ -0,0 +1,19 @@
+export FixedRNG
+struct FixedRNG <: AbstractRNG end
+
+Base.rand(::FixedRNG) = one(Float64) / 2
+Random.randn(::FixedRNG) = zero(Float64)
+Random.randexp(::FixedRNG) = one(Float64)
+
+Base.rand(::FixedRNG, ::Type{T}) where {T<:Real} = one(T) / 2
+Random.randn(::FixedRNG, ::Type{T}) where {T<:Real} = zero(T)
+Random.randexp(::FixedRNG, ::Type{T}) where {T<:Real} = one(T)
+
+# We need concrete type parameters to avoid amiguity for these cases
+for T in [Float16, Float32, Float64]
+    @eval begin
+        Base.rand(::FixedRNG, ::Type{$T}) = one($T) / 2
+        Random.randn(::FixedRNG, ::Type{$T}) = zero($T)
+        Random.randexp(::FixedRNG, ::Type{$T}) = one($T)
+    end
+end

src/interface.jl

@@ -54,15 +54,15 @@ function test_interface(μ::M) where {M}
             ###########################################################################
             # testvalue, logdensityof
 
-            x = @inferred testvalue(μ)
+            x = @inferred testvalue(Float64, μ)
             β = @inferred basemeasure(μ, x)
 
             ℓμ = @inferred logdensityof(μ, x)
             ℓβ = @inferred logdensityof(β, x)
 
             @test ℓμ ≈ logdensity_def(μ, x) + ℓβ
 
-            @test logdensity_def(μ, testvalue(μ)) isa Real
+            @test logdensity_def(μ, testvalue(Float64, μ)) isa Real
         end
     end
 end

src/primitives/counting.jl

@@ -1,10 +1,10 @@
 # Counting measure
 
-export Counting, CountingMeasure
+export Counting, CountingBase
 
-struct CountingMeasure <: PrimitiveMeasure end
+struct CountingBase <: PrimitiveMeasure end
 
-insupport(::CountingMeasure, x) = true
+insupport(::CountingBase, x) = true
 
 struct Counting{T} <: AbstractMeasure
     support::T
@@ -16,15 +16,15 @@ function logdensity_def(μ::Counting, x)
     insupport(μ, x) ? 0.0 : -Inf
 end
 
-basemeasure(::Counting) = CountingMeasure()
+basemeasure(::Counting) = CountingBase()
 
 Counting() = Counting(ℤ)
 
-testvalue(d::Counting) = testvalue(d.support)
+testvalue(::Type{T}, d::Counting) where {T} = testvalue(T, d.support)
 
-proxy(d::Counting) = restrict(in(d.support), CountingMeasure())
+proxy(d::Counting) = restrict(in(d.support), CountingBase())
 
-Base.:∘(::typeof(basemeasure), ::Type{Counting}) = CountingMeasure()
+Base.:∘(::typeof(basemeasure), ::Type{Counting}) = CountingBase()
 
 Base.show(io::IO, d::Counting) = print(io, "Counting(", d.support, ")")
 

src/primitives/dirac.jl

@@ -16,7 +16,7 @@ function (μ::Dirac{X})(s) where {X}
     return 0
 end
 
-basemeasure(d::Dirac) = CountingMeasure()
+basemeasure(d::Dirac) = CountingBase()
 
 logdensity_def(μ::Dirac, x) = 0.0
 
@@ -26,8 +26,6 @@ export dirac
 
 dirac(d::AbstractMeasure) = Dirac(rand(d))
 
-testvalue(d::Dirac) = d.x
-
 insupport(d::Dirac, x) = x == d.x
 
 @inline getdof(::Dirac) = static(0)

src/primitives/lebesgue.jl

@@ -2,13 +2,13 @@
 
 export Lebesgue
 
-struct LebesgueMeasure <: PrimitiveMeasure end
+struct LebesgueBase <: PrimitiveMeasure end
 
-testvalue(::LebesgueMeasure) = 0.0
+testvalue(::LebesgueBase) = 0.0
 
-insupport(::LebesgueMeasure, x) = true
+insupport(::LebesgueBase, x) = true
 
-insupport(::LebesgueMeasure) = Returns(true)
+insupport(::LebesgueBase) = Returns(true)
 
 struct Lebesgue{T} <: AbstractMeasure
     support::T
@@ -22,29 +22,29 @@ gentype(::Lebesgue) = Float64
 
 Lebesgue() = Lebesgue(ℝ)
 
-# basemeasure(::Lebesgue) = LebesgueMeasure()
+# basemeasure(::Lebesgue) = LebesgueBase()
 
-testvalue(d::Lebesgue) = testvalue(d.support)
+testvalue(::Type{T}, d::Lebesgue) where {T} = testvalue(T, d.support)::T
 
-proxy(d::Lebesgue) = restrict(in(d.support), LebesgueMeasure())
+proxy(d::Lebesgue) = restrict(in(d.support), LebesgueBase())
 @useproxy Lebesgue
 
-Base.:∘(::typeof(basemeasure), ::Type{Lebesgue}) = LebesgueMeasure()
+Base.:∘(::typeof(basemeasure), ::Type{Lebesgue}) = LebesgueBase()
 
 Base.show(io::IO, d::Lebesgue) = print(io, "Lebesgue(", d.support, ")")
 
 insupport(μ::Lebesgue, x) = x ∈ μ.support
 
 insupport(::Lebesgue{RealNumbers}, ::Real) = true
 
-logdensity_def(::LebesgueMeasure, ::CountingMeasure, x) = -Inf
+logdensity_def(::LebesgueBase, ::CountingBase, x) = -Inf
 
-logdensity_def(::CountingMeasure, ::LebesgueMeasure, x) = Inf
+logdensity_def(::CountingBase, ::LebesgueBase, x) = Inf
 
-@inline getdof(::Lebesgue) = static(1)
+@inline getdof(::LebesgueBase) = static(1)
 
-@inline checked_arg(::Lebesgue, x::Real) = x
+@inline checked_arg(::LebesgueBase, x::Real) = x
 
-@propagate_inbounds function checked_arg(::Lebesgue, x::Any)
+@propagate_inbounds function checked_arg(::LebesgueBase, x::Any)
     @boundscheck throw(ArgumentError("Invalid variate type for measure"))
 end

src/standard/stdexponential.jl

@@ -5,7 +5,7 @@ export StdExponential
 insupport(d::StdExponential, x) = x ≥ zero(x)
 
 @inline logdensity_def(::StdExponential, x) = -x
-@inline basemeasure(::StdExponential) = Lebesgue()
+@inline basemeasure(::StdExponential) = LebesgueBase()
 
 @inline transport_def(::StdUniform, μ::StdExponential, x) = -expm1(-x)
 @inline transport_def(::StdExponential, μ::StdUniform, x) = -log1p(-x)

src/standard/stdlogistic.jl

@@ -5,7 +5,7 @@ export StdLogistic
 @inline insupport(d::StdLogistic, x) = true
 
 @inline logdensity_def(::StdLogistic, x) = (u = -abs(x); u - 2 * log1pexp(u))
-@inline basemeasure(::StdLogistic) = Lebesgue()
+@inline basemeasure(::StdLogistic) = LebesgueBase()
 
 @inline transport_def(::StdUniform, μ::StdLogistic, x) = logistic(x)
 @inline transport_def(::StdLogistic, μ::StdUniform, x) = logit(x)

src/standard/stdmeasure.jl

@@ -2,6 +2,7 @@ abstract type StdMeasure <: AbstractMeasure end
 
 StdMeasure(::typeof(rand)) = StdUniform()
 StdMeasure(::typeof(randexp)) = StdExponential()
+StdMeasure(::typeof(randn)) = StdNormal()
 
 @inline check_dof(::StdMeasure, ::StdMeasure) = nothing
 

src/standard/stdnormal.jl

@@ -0,0 +1,20 @@
+using SpecialFunctions: erfc, erfcinv
+using IrrationalConstants: invsqrt2
+
+struct StdNormal <: StdMeasure end
+
+export StdNormal
+
+@inline insupport(d::StdNormal, x) = true
+
+@inline logdensity_def(::StdNormal, x) = -x^2 / 2
+@inline basemeasure(::StdNormal) = WeightedMeasure(static(-0.5 * log2π), LebesgueBase())
+
+@inline getdof(::StdNormal) = static(1)
+
+@inline transport_def(::StdUniform, μ::StdNormal, x) = erfc(-x * invsqrt2) / 2
+@inline transport_def(::StdNormal, μ::StdUniform, p) = -erfcinv(2*p) * sqrt2
+
+@inline Base.rand(rng::Random.AbstractRNG, ::Type{T}, ::StdNormal) where {T} = randn(rng, T)
+
+

src/standard/stduniform.jl

@@ -5,6 +5,6 @@ export StdUniform
 insupport(d::StdUniform, x) = zero(x) ≤ x ≤ one(x)
 
 @inline logdensity_def(::StdUniform, x) = zero(x)
-@inline basemeasure(::StdUniform) = Lebesgue()
+@inline basemeasure(::StdUniform) = LebesgueBase()
 
 Base.rand(rng::Random.AbstractRNG, ::Type{T}, ::StdUniform) where {T} = rand(rng, T)