Dev

.github/workflows/CI.yml

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

Project.toml

@@ -11,6 +11,7 @@ ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
 IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -21,6 +22,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"
@@ -34,14 +36,16 @@ ConstructionBase = "1.3"
 DensityInterface = "0.4"
 FillArrays = "0.12, 0.13"
 IfElse = "0.1"
-InverseFunctions = "0.1.7"
+IntervalSets = "0.7"
+InverseFunctions = "0.1.8"
 IrrationalConstants = "0.1"
 LogExpFunctions = "0.3"
 LogarithmicNumbers = "1"
 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

@@ -8,9 +8,16 @@ import Random: gentype
 using Statistics
 using LinearAlgebra
 
+import IntervalSets
+# This seems harder than it should be to get `IntervalSets.:(..)`
+@eval (using IntervalSets: $(Symbol(IntervalSets.:(..))))
+
+using IntervalSets: Interval, width
+
 import DensityInterface: logdensityof
 import DensityInterface: densityof
 import DensityInterface: DensityKind
+using DensityInterface: FuncDensity, LogFuncDensity
 using DensityInterface
 
 using InverseFunctions
@@ -19,6 +26,7 @@ using ChangesOfVariables
 import Base.iterate
 import ConstructionBase
 using ConstructionBase: constructorof
+using IntervalSets
 
 using PrettyPrinting
 const Pretty = PrettyPrinting
@@ -108,17 +116,18 @@ using Compat
 
 using IrrationalConstants
 
+include("smf.jl")
 include("getdof.jl")
 include("transport.jl")
 include("schema.jl")
 include("splat.jl")
 include("proxies.jl")
 include("kernel.jl")
 include("parameterized.jl")
-include("combinators/half.jl")
 include("domains.jl")
 include("primitive.jl")
 include("utils.jl")
+include("mass-interface.jl")
 # include("absolutecontinuity.jl")
 
 include("primitives/counting.jl")
@@ -144,9 +153,11 @@ include("standard/stdmeasure.jl")
 include("standard/stduniform.jl")
 include("standard/stdexponential.jl")
 include("standard/stdlogistic.jl")
-include("latent-joint.jl")
+include("standard/stdnormal.jl")
+include("combinators/half.jl")
 
 include("rand.jl")
+include("fixedrng.jl")
 
 include("density.jl")
 include("density-core.jl")

src/combinators/half.jl

@@ -26,4 +26,18 @@ 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)
+
+massof(μ::Half) = massof(unhalf(μ))
+
+function smf(μ::Half, x)
+    2 * smf(μ.parent, max(x, zero(x))) - 1
+end
+
+function invsmf(μ::Half, p)
+    @assert zero(p) ≤ p ≤ one(p)
+    invsmf(μ.parent, (p + 1) / 2)
+end
+
+transport_def(μ::Half, ::StdUniform, p) = invsmf(μ, p)
+transport_def(::StdUniform, μ::Half, x) = smf(μ, x)

src/combinators/likelihood.jl

@@ -202,4 +202,6 @@ more efficient than
 
     logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x)
 """
-likelihood_ratio(ℓ::Likelihood, p, q) = exp(ULogarithmic, logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x))
+function likelihood_ratio(ℓ::Likelihood, p, q)
+    exp(ULogarithmic, logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x))
+end

src/combinators/power.jl

@@ -128,3 +128,5 @@ end
 function checked_arg(μ::PowerMeasure, x::Any)
     throw(ArgumentError("Size of variate doesn't match size of power measure"))
 end
+
+massof(m::PowerMeasure) = massof(m.parent)^prod(m.axes)

src/combinators/product.jl

@@ -16,6 +16,8 @@ function Pretty.tile(μ::AbstractProductMeasure)
     result *= Pretty.literal(")")
 end
 
+massof(m::AbstractProductMeasure) = prod(massof, marginals(m))
+
 export marginals
 
 function Base.:(==)(a::AbstractProductMeasure, b::AbstractProductMeasure)
@@ -161,7 +163,9 @@ 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))
+function testvalue(::Type{T}, d::AbstractProductMeasure) where {T}
+    _map(m -> testvalue(T, m), marginals(d))
+end
 
 export ⊗
 
@@ -220,3 +224,16 @@ end
     end
     return true
 end
+
+getdof(d::AbstractProductMeasure) = mapreduce(getdof, +, marginals(d))
+
+function checked_arg(μ::ProductMeasure{<:NTuple{N,Any}}, x::NTuple{N,Any}) where {N}
+    map(checked_arg, marginals(μ), x)
+end
+
+function checked_arg(
+    μ::ProductMeasure{<:NamedTuple{names}},
+    x::NamedTuple{names},
+) where {names}
+    NamedTuple{names}(map(checked_arg, values(marginals(μ)), values(x)))
+end

src/combinators/smart-constructors.jl

@@ -85,7 +85,7 @@ superpose(nt::NamedTuple) = SuperpositionMeasure(nt)
 
 function superpose(μ::T, ν::T) where {T<:AbstractMeasure}
     if μ == ν
-        return weightedmeasure(logtwo, μ)
+        return weightedmeasure(static(float(logtwo)), μ)
     else
         return superpose((μ, ν))
     end

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

@@ -4,6 +4,8 @@ using LogExpFunctions
 
 export SuperpositionMeasure
 
+abstract type AbstractSuperpositionMeasure <: AbstractMeasure end
+
 @doc raw"""
     struct SuperpositionMeasure{NT} <: AbstractMeasure
         components :: NT
@@ -24,17 +26,19 @@ Superposition measures satisfy
     \end{aligned}
 ```
 """
-struct SuperpositionMeasure{C} <: AbstractMeasure
+struct SuperpositionMeasure{C} <: AbstractSuperpositionMeasure
     components::C
 end
 
+massof(m::SuperpositionMeasure) = sum(massof, m.components)
+
 function Pretty.tile(d::SuperpositionMeasure)
     result = Pretty.literal("SuperpositionMeasure(")
     result *= Pretty.list_layout([Pretty.tile.(d.components)...])
     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

@@ -1,5 +1,7 @@
 # TODO: Compare with ChangesOfVariables.jl
 
+using InverseFunctions: FunctionWithInverse
+
 abstract type AbstractTransformedMeasure <: AbstractMeasure end
 
 abstract type AbstractPushforward <: AbstractTransformedMeasure end
@@ -17,16 +19,15 @@ function parent(::AbstractTransformedMeasure) end
 export PushforwardMeasure
 
 """
-    struct PushforwardMeasure{FF,IF,MU,VC<:TransformVolCorr} <: AbstractPushforward
-        f :: FF
-        inv_f :: IF
-        origin :: MU
-        volcorr :: VC
-    end
+    struct PushforwardMeasure{F,I,MU,VC<:TransformVolCorr} <:
+        AbstractPushforward f :: F finv :: I origin :: MU volcorr :: VC end
+
+    Users should not call `PushforwardMeasure` directly. Instead call or add
+    methods to `pushfwd`.
 """
-struct PushforwardMeasure{FF,IF,M,VC<:TransformVolCorr} <: AbstractPushforward
-    f::FF
-    inv_f::IF
+struct PushforwardMeasure{F,I,M,VC<:TransformVolCorr} <: AbstractPushforward
+    f::F
+    finv::I
     origin::M
     volcorr::VC
 end
@@ -35,14 +36,25 @@ gettransform(ν::PushforwardMeasure) = ν.f
 parent(ν::PushforwardMeasure) = ν.origin
 
 function Pretty.tile(ν::PushforwardMeasure)
-    Pretty.list_layout(Pretty.tile.([ν.f, ν.inv_f, ν.origin]); prefix = :PushforwardMeasure)
+    Pretty.list_layout(Pretty.tile.([ν.f, ν.origin]); prefix = :PushforwardMeasure)
 end
 
-@inline function logdensity_def(
-    ν::PushforwardMeasure{FF,IF,M,<:WithVolCorr},
-    y,
-) where {FF,IF,M}
-    x_orig, inv_ladj = with_logabsdet_jacobian(ν.inv_f, y)
+# TODO: THIS IS ALMOST CERTAINLY WRONG 
+# @inline function logdensity_rel(
+#     ν::PushforwardMeasure{FF1,IF1,M1,<:WithVolCorr},
+#     β::PushforwardMeasure{FF2,IF2,M2,<:WithVolCorr},
+#     y,
+# ) where {FF1,IF1,M1,FF2,IF2,M2}
+#     x = β.inv_f(y)
+#     f = ν.inv_f ∘ β.f
+#     inv_f = β.inv_f ∘ ν.f
+#     logdensity_rel(pushfwd(f, inv_f, ν.origin, WithVolCorr()), β.origin, x)
+# end
+
+@inline function logdensity_def(ν::PushforwardMeasure{F,I,M,<:WithVolCorr}, y) where {F,I,M}
+    f = ν.f
+    finv = ν.finv
+    x_orig, inv_ladj = with_logabsdet_jacobian(finv, y)
     logd_orig = logdensity_def(ν.origin, x_orig)
     logd = float(logd_orig + inv_ladj)
     neginf = oftype(logd, -Inf)
@@ -57,49 +69,93 @@ end
     )
 end
 
-@inline function logdensity_def(
-    ν::PushforwardMeasure{FF,IF,M,<:NoVolCorr},
-    y,
-) where {FF,IF,M}
-    x_orig = to_origin(ν, y)
-    return logdensity_def(ν.origin, x_orig)
+@inline function logdensity_def(ν::PushforwardMeasure{F,I,M,<:NoVolCorr}, y) where {F,I,M}
+    x = ν.finv(y)
+    return logdensity_def(ν.origin, x)
 end
 
-insupport(ν::PushforwardMeasure, y) = insupport(transport_origin(ν), to_origin(ν, y))
+insupport(ν::PushforwardMeasure, y) = insupport(ν.origin, ν.finv(y))
 
-testvalue(ν::PushforwardMeasure) = from_origin(ν, testvalue(transport_origin(ν)))
+function testvalue(::Type{T}, ν::PushforwardMeasure) where {T}
+    ν.f(testvalue(T, parent(ν)))
+end
 
 @inline function basemeasure(ν::PushforwardMeasure)
-    PushforwardMeasure(ν.f, ν.inv_f, basemeasure(transport_origin(ν)), NoVolCorr())
+    pushfwd(ν.f, basemeasure(parent(ν)), NoVolCorr())
 end
 
 _pushfwd_dof(::Type{MU}, ::Type, dof) where {MU} = NoDOF{MU}()
 _pushfwd_dof(::Type{MU}, ::Type{<:Tuple{Any,Real}}, dof) where {MU} = dof
 
-# Assume that DOF are preserved if with_logabsdet_jacobian is functional:
-@inline function getdof(ν::MU) where {MU<:PushforwardMeasure}
-    T = Core.Compiler.return_type(testvalue, Tuple{typeof(ν.origin)})
-    R = Core.Compiler.return_type(with_logabsdet_jacobian, Tuple{typeof(ν.f),T})
-    _pushfwd_dof(MU, R, getdof(ν.origin))
-end
+@inline getdof(ν::MU) where {MU<:PushforwardMeasure} = getdof(ν.origin)
 
 # Bypass `checked_arg`, would require potentially costly transformation:
 @inline checked_arg(::PushforwardMeasure, x) = x
 
 @inline transport_origin(ν::PushforwardMeasure) = ν.origin
 @inline from_origin(ν::PushforwardMeasure, x) = ν.f(x)
-@inline to_origin(ν::PushforwardMeasure, y) = ν.inv_f(y)
+@inline to_origin(ν::PushforwardMeasure, y) = ν.finv(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
 
+###############################################################################
+# pushfwd
+
 export pushfwd
 
+
+
 """
     pushfwd(f, μ, 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`.
+
+To manually specify an inverse, call 
+`pushfwd(InverseFunctions.setinverse(f, finv), μ, volcorr)`.
 """
-pushfwd(f, μ, volcorr = WithVolCorr()) = PushforwardMeasure(f, inverse(f), μ, volcorr)
+function pushfwd(f, μ, volcorr::TransformVolCorr = WithVolCorr())
+    PushforwardMeasure(f, inverse(f), μ, volcorr)
+end
+
+function pushfwd(f, μ::PushforwardMeasure, volcorr::TransformVolCorr = WithVolCorr())
+    _pushfwd_of_pushfwd(f, μ, μ.volcorr, volcorr)
+end
+
+function pushfwd(f, μ::PushforwardMeasure, ::WithVolCorr)
+    _pushfwd_of_pushfwd(f, μ, μ.volcorr, WithVolCorr())
+end
+
+# Either both WithVolCorr or both NoVolCorr, so we can merge them
+function _pushfwd_of_pushfwd(f, μ::PushforwardMeasure, ::V, v::V) where {V}
+    pushfwd(f ∘ μ.f, μ.origin, v)
+end
+
+function _pushfwd_of_pushfwd(f, μ::PushforwardMeasure, _, v)
+    PushforwardMeasure(f, inverse(f), μ, v)
+end
+
+###############################################################################
+# pullback
+
+"""
+    pullback(f, μ, volcorr = WithVolCorr())
+
+A _pullback_ is a dual concept to a _pushforward_. While a pushforward needs a
+map _from_ the support of a measure, a pullback requires a map _into_ the
+support of a measure. The log-density is then computed through function
+composition, together with a volume correction as needed.
+
+This can be useful, since the log-density of a `PushforwardMeasure` is computing
+in terms of the inverse function; the "forward" function is not used at all. In
+some cases, we may be focusing on log-density (and not, for example, sampling).
+
+To manually specify an inverse, call 
+`pullback(InverseFunctions.setinverse(f, finv), μ, volcorr)`.
+"""
+function pullback(f, μ, volcorr::TransformVolCorr = WithVolCorr())
+    pushfwd(setinverse(inverse(f), f), μ, volcorr)
+end