Dev

Project.toml

@@ -1,7 +1,7 @@
 name = "MeasureBase"
 uuid = "fa1605e6-acd5-459c-a1e6-7e635759db14"
 authors = ["Chad Scherrer <[email protected]> and contributors"]
-version = "0.7.0"
+version = "0.8.0"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

src/MeasureBase.jl

@@ -13,6 +13,7 @@ import DensityInterface: densityof
 import DensityInterface: DensityKind
 using DensityInterface
 
+import Base.iterate
 import ConstructionBase
 using ConstructionBase: constructorof
 
@@ -31,7 +32,7 @@ export AbstractMeasure
 import IfElse: ifelse
 export logdensity_def
 export basemeasure
-export basekleisli
+export basekernel
 
 """
     inssupport(m, x)
@@ -97,23 +98,21 @@ using Compat
 include("schema.jl")
 include("splat.jl")
 include("proxies.jl")
-include("kleisli.jl")
+include("kernel.jl")
 include("parameterized.jl")
 include("combinators/half.jl")
 include("domains.jl")
 include("primitive.jl")
 include("utils.jl")
-include("absolutecontinuity.jl")
+# include("absolutecontinuity.jl")
 
 include("primitives/counting.jl")
 include("primitives/lebesgue.jl")
 include("primitives/dirac.jl")
 include("primitives/trivial.jl")
 
-include("combinators/conditional.jl")
 include("combinators/bind.jl")
 include("combinators/transformedmeasure.jl")
-include("combinators/factoredbase.jl")
 include("combinators/weighted.jl")
 include("combinators/superpose.jl")
 include("combinators/product.jl")
@@ -123,6 +122,8 @@ include("combinators/likelihood.jl")
 include("combinators/pointwise.jl")
 include("combinators/restricted.jl")
 include("combinators/smart-constructors.jl")
+include("combinators/powerweighted.jl")
+include("combinators/conditional.jl")
 
 include("rand.jl")
 

src/combinators/conditional.jl

@@ -25,4 +25,16 @@ upcoming changes.
 """
 Base.:|(μ::AbstractMeasure, constraint) = ConditionalMeasure(μ, constraint)
 
-@inline basemeasure(cm::ConditionalMeasure) = basemeasure(cm.parent) | cm.constraint
+@inline basemeasure(cm::ConditionalMeasure) = basemeasure(cm.parent) | cm.constraint
+
+# @generated function Base.:|(μ::ProductMeasure{NamedTuple{M,T}}, constraint::NamedTuple{N}) where {M,T,N}
+#     newkeys = tuple(setdiff(M, N)...)
+#     quote
+#         mar = marginals(μ)
+#         productmeasure(NamedTuple{$newkeys}(mar))
+#     end
+# end
+
+function Base.:|(μ::ProductMeasure{NamedTuple{M,T}}, constraint::NamedTuple{N}) where {M,T,N}
+    productmeasure(merge(marginals(μ),rmap(Dirac, constraint)))
+end

src/combinators/half.jl

@@ -11,13 +11,8 @@ end
 
 unhalf(μ::Half) = μ.parent
 
-isnonnegative(x) = x ≥ 0.0
-
 @inline function basemeasure(μ::Half)
-    constℓ = static(logtwo)
-    varℓ = Returns(0.0)
-    base = basemeasure(unhalf(μ))
-    FactoredBase(constℓ, varℓ, base)
+    weightedmeasure(static(logtwo), basemeasure(unhalf(μ)))
 end
 
 function Base.rand(rng::AbstractRNG, ::Type{T}, μ::Half) where {T}
@@ -27,7 +22,8 @@ end
 logdensity_def(μ::Half, x) = logdensity_def(unhalf(μ), x)
 
 @inline function insupport(d::Half, x)
-    ifelse(isnonnegative(x), insupport(unhalf(d), x), false)
+    x ≥ 0 || return false
+    insupport(unhalf(d), x)
 end
 
 testvalue(::Half) = 1.0

src/combinators/likelihood.jl

@@ -2,10 +2,16 @@ export AbstractLikelihood, Likelihood
 
 abstract type AbstractLikelihood end
 
-@inline logdensityof(ℓ::AbstractLikelihood, par) = logdensity_def(ℓ, par)
+# @inline function logdensityof(ℓ::AbstractLikelihood, p)
+#     t() = dynamic(unsafe_logdensityof(ℓ, p))
+#     f() = -Inf
+#     ifelse(insupport(ℓ, p), t, f)()
+# end
+
+# insupport(ℓ::AbstractLikelihood, p) = insupport(ℓ.k(p), ℓ.x)
 
 @doc raw"""
-    Likelihood(k::AbstractKleisli, x)
+    Likelihood(k::AbstractTransitionKernel, x)
 
 "Observe" a value `x`, yielding a function from the parameters to ℝ.
 
@@ -109,9 +115,9 @@ struct Likelihood{K,X} <: AbstractLikelihood
     k::K
     x::X
 
-    Likelihood(k::K, x::X) where {K<:AbstractKleisli,X} = new{K,X}(k,x)
+    Likelihood(k::K, x::X) where {K<:AbstractTransitionKernel,X} = new{K,X}(k,x)
     Likelihood(k::K, x::X) where {K<:Function,X} = new{K,X}(k,x)
-    Likelihood(μ, x) = Likelihood(kleisli(μ), x)
+    Likelihood(μ, x) = Likelihood(kernel(μ), x)
 end
 
 # Not really a density, but this makes the code work
@@ -128,10 +134,37 @@ function Base.show(io::IO, ℓ::Likelihood)
     Pretty.pprint(io, ℓ)
 end
 
-@inline function logdensity_def(ℓ::Likelihood, p::Tuple)
-    return logdensity_def(ℓ.k(p), ℓ.x)
-end
+# @inline function logdensity_def(ℓ::Likelihood, p)
+#     return logdensity_def(ℓ.k(p), ℓ.x)
+# end
 
-@inline function logdensity_def(ℓ::Likelihood, p)
-    return logdensity_def(ℓ.k((p,)), ℓ.x)
-end
+# basemeasure(ℓ::Likelihood, p) = basemeasure(ℓ.k(p), ℓ.x)
+
+# basemeasure(ℓ::Likelihood) = @error "Likelihood requires local base measure"
+
+export likelihood
+
+likelihood(k, x, ::NamedTuple{()}) = Likelihood(k, x)
+
+likelihood(k, x; kwargs...) = likelihood(k, x, NamedTuple(kwargs))
+
+likelihood(k, x, pars::NamedTuple) = likelihood(kernel(k, pars), x)
+
+likelihood(k::AbstractTransitionKernel, x) = Likelihood(k, x)
+
+"""
+    log_likelihood_ratio(ℓ::Likelihood, p, q)
+
+Compute the log of the likelihood ratio, in order to compare two choices for
+parameters. This is computed as
+
+    logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x)
+
+Since `logdensity_rel` can leave common base measure unevaluated, this can be
+more efficient than
+
+    logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x)
+"""
+log_likelihood_ratio(ℓ::Likelihood, p, q) = logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x)
+
+# likelihood(k, x; kwargs...) = likelihood(k, x, NamedTuple(kwargs))

src/combinators/pointwise.jl

@@ -1,18 +1,17 @@
 export ⊙
 
-struct PointwiseProductMeasure{M,L} <: AbstractMeasure
-    measure::M
+struct PointwiseProductMeasure{P,L} <: AbstractMeasure
+    prior::P
     likelihood::L
-
-    function PointwiseProductMeasure(μ::M, ℓ::L) where {M,L}
-        @assert static_hasmethod(logdensity_def, Tuple{L, gentype(μ)})
-        return new{M,L}(μ, ℓ)
-    end
 end
 
+
+
+iterate(p::PointwiseProductMeasure, i=1) = iterate((p.prior, p.likelihood), i)
+
 function Base.show(io::IO, μ::PointwiseProductMeasure)
     io = IOContext(io, :compact => true)
-    print(io, μ.measure, " ⊙ ", μ.likelihood)
+    print(io, μ.prior, " ⊙ ", μ.likelihood)
 end
 
 function Base.show_unquoted(io::IO, μ::PointwiseProductMeasure, indent::Int, prec::Int)
@@ -29,12 +28,20 @@ end
 
 ⊙(μ, ℓ) = pointwiseproduct(μ, ℓ)
 
-@inline function logdensity_def(d::PointwiseProductMeasure, x)
-    logdensity_def(d.measure, x) + logdensity_def(d.likelihood, x)
+@inline function logdensity_def(d::PointwiseProductMeasure, p)
+    μ, ℓ = d
+    logdensityof(ℓ.k(p), ℓ.x)
 end
 
 function gentype(d::PointwiseProductMeasure)
-    @inbounds gentype(d.measure)
+    gentype(d.prior)
+end
+
+@inbounds function insupport(d::PointwiseProductMeasure, p) 
+    μ, ℓ = d
+    insupport(μ, p) && insupport(ℓ.k(p), ℓ.x)
 end
 
-basemeasure(d::PointwiseProductMeasure) = @inbounds basemeasure(d.measure)
+basemeasure(d::PointwiseProductMeasure, x) = d.prior
+
+basemeasure(d::PointwiseProductMeasure) = basemeasure(d.prior)

src/combinators/power.jl

@@ -65,14 +65,10 @@ params(d::PowerMeasure) = params(first(marginals(d)))
 end
 
 @inline function logdensity_def(d::PowerMeasure{M}, x) where {M}
-    ℓ = 0.0
-    # ℓ = zero(typeintersect(AbstractFloat,Core.Compiler.return_type(logdensity_def, Tuple{M,T})))
     parent = d.parent
-    @simd for xj in x
-        Δℓ = logdensity_def(parent, xj)
-        ℓ += Δℓ
+    sum(x) do xj
+        logdensity_def(parent, xj)
     end
-    ℓ
 end
 
 @generated function logdensity_def(d::PowerMeasure{M, Tuple{Base.OneTo{StaticInt{N}}}}, x) where {M,N}

src/combinators/powerweighted.jl

@@ -0,0 +1,33 @@
+export ↑
+
+struct PowerWeightedMeasure{M,A} <: AbstractMeasure
+    parent::M
+    exponent::A
+end
+
+logdensity_def(d::PowerWeightedMeasure, x) = d.exponent * logdensity_def(d.parent, x)
+
+basemeasure(d::PowerWeightedMeasure, x) = basemeasure(d.parent, x) ↑ d.exponent
+
+basemeasure(d::PowerWeightedMeasure) = basemeasure(d.parent) ↑ d.exponent
+
+function powerweightedmeasure(d, α)
+    isone(α) && return d
+    PowerWeightedMeasure(d, α)
+end
+
+(d::AbstractMeasure) ↑ α = powerweightedmeasure(d, α)
+
+insupport(d::PowerWeightedMeasure, x) = insupport(d.parent, x)
+
+function Base.show(io::IO, d::PowerWeightedMeasure)
+    print(io, d.parent, " ↑ ", d.exponent)
+end
+
+function powerweightedmeasure(d::PowerWeightedMeasure, α)
+    powerweightedmeasure(d.parent, α * d.exponent)
+end
+
+function powerweightedmeasure(d::WeightedMeasure, α)
+    weightedmeasure(α*d.logweight, powerweightedmeasure(d.base, α))
+end

src/combinators/product.jl

@@ -94,19 +94,19 @@ end
     return q
 end
 
-@generated function basemeasure(d::ProductMeasure{NamedTuple{N,T}}, x) where {N,T}
-    q = quote
-        m = marginals(d)
-    end
-    for k in N
-        qk = QuoteNode(k)
-        push!(q.args, :($k = basemeasure(getproperty(m, $qk))))
-    end
-
-    vals = map(x -> Expr(:(=), x,x), N)
-    push!(q.args, Expr(:tuple, vals...))
-    return q
-end
+# @generated function basemeasure(d::ProductMeasure{NamedTuple{N,T}}, x) where {N,T}
+#     q = quote
+#         m = marginals(d)
+#     end
+#     for k in N
+#         qk = QuoteNode(k)
+#         push!(q.args, :($k = basemeasure(getproperty(m, $qk))))
+#     end
+
+#     vals = map(x -> Expr(:(=), x,x), N)
+#     push!(q.args, Expr(:tuple, vals...))
+#     return q
+# end
 
 function basemeasure(μ::ProductMeasure{Base.Generator{I,F}}) where {I,F}
     mar = marginals(μ)
@@ -133,7 +133,7 @@ end
 
 function _basemeasure(μ::ProductMeasure{Base.Generator{I,F}}, ::Type{B}, ::False) where {I,F,B}
     mar = marginals(μ)
-    productmeasure(Base.Generator(basekleisli(mar.f), mar.iter))
+    productmeasure(Base.Generator(basekernel(mar.f), mar.iter))
 end
 
 marginals(μ::ProductMeasure) = μ.marginals
@@ -174,7 +174,7 @@ using Random: rand!, GLOBAL_RNG, AbstractRNG
 end
 
 export rand!
-using Random: rand!, GLOBAL_RNG, AbstractRNG
+using Random: rand!, GLOBAL_RNG
 
 function _rand(rng::AbstractRNG, ::Type{T}, d::ProductMeasure, mar::AbstractArray) where {T}
     elT = typeof(rand(rng, T, first(mar)))
@@ -187,8 +187,8 @@ end
 
 @inline function insupport(d::AbstractProductMeasure, x::AbstractArray)
     mar = marginals(d)
-    for j in eachindex(x)
-        @inbounds dynamic(insupport(mar[j], x[j])) || return false
+    for (j,mj) in enumerate(mar)
+        dynamic(insupport(mj, x[j])) || return false
     end
     return true
 end