Dev (#49)

* drop some unneeded methods

* update for DensityInterface

* updating densities to DensityInterface approach

* update domains

* moving things around

* updates to IntegerBounds

* more domain mucking

* fix some exports

* update deps

* update `using`

* working on tests

* Working toward tests passing

* some refactoring

* working on tests

* cleaning up

* get tests to pass

* tests passing

* update logdensity_def for pointwise product

* PrettyPrinting + tests

* speed up `rootmeasure`

* oops didn't mean to include that

* tile(::FactoredBase)

* Update Half and FactoredBase

* drop exp.jl

* simplify basemeasure

* drop old `include`

* drop redundant method

* update compat Returns

* add AbstractDensity

* Move Affine to MeasureTheory

* drop some old `For` code

* update counting measure

* add testvalue(::Type{T})

* bugfix

* some dispatch adjustments

* simplify show

* formatting

* transformed measures

* ZeroSets

* using LinearAlgebra, Statistics

* working on MeasureTheory tests

* updates

* updates

* typo

* update default to mimic Base

* fix tile(::Lebesgue)

* Add test_interface function

* small doc update

* DensityKind(::Likelihood)

* fixing show(::Likelihood)

* adding some docs

* update `rand` method

* drop old integration code

* add `rebase`

* export rebase

* law for ⊗

* typo

* Make Likelihood more flexible

* update kernel methods

* add Likelihood method (avoid stack overflow)

* refactoring

* compat

* Maybe Comat just works?

* refactoring

* some new stuff

* working on tests

* update powermeasure combinator

* bugfix

* comment out debugging lines

* more refactoring

* test @inferred basemeasure_depth(μ)

* drop `constructor` (just use ConstructionBase.constructorof`)

* debugging

* update help

* update interface

* make tests harder

* fixes

* Dirac bugfix

* formatting

* improve type inference

* working on type inference

* update interface

* udpates

* get test passing

* @test !isabstracttype(typejoin(...))

* work on show methods

* update CI

* remove old code

* update productmeasure

* prettyprinting stuff

* Drop te @constprop :aggressive stuff (maybe don't need it?)

* nerline

* dropping some old code

* update tbasemeasure_type(::PowerMeasure)

* moar tests

* update SpikeMixture

* update superpose type parameter name

* drop old tests

* func_string

* more updates

* getting closer

* almost there!

* generated function for type stability

* tests passing!

* newline

* more fixes

* exports and bugfix

* insupport(μ::Counting{T}, x) where {T<:Type}

* working on MeasureTheory tests

* MeasureTheory tests passing

* drop some old code

* inlining

* improve inference

* update `tile(::For)`

* tighten down infrerence

* update basemeasure(::For) for generators

* loosen type bound on instance_type

* drop debugging code

* small update for Likelihood, and a test

* fixing up likelihoods

* improve `basemeasure_depth` dispatch

* still some trouble with inferred basemeasure_depth

* clean up `For` dispatch

* simplify _logdensityof

* optimize for Returns{True} case

* rework basemeasure_depth

* aggressive tests passing!!

* drop type-level stuff

* drop help

* license

* affero

* copyright notice

* merge

* Drop Create Commons

* cleanup after merge

* update support computations

* insupport(d::SuperpositionMeasure, x)

* dorp ParamWeighted

* insupport(d::FactoredBase, x)

* export unsafe_logdensityof

* call promote_type instead of promote_rule

* logdensity_def for named tuple product measures

* type annotation for now

* debugging

* drop shows

* speed up mapped arrays

* throw an error for `Union{}` types

* MT tests passing

* updates

* get tests passing

* MIT license for MeasureBase

* bump version

* cleanup

* spacing

* Move ConditionalMeasure to MeasureBase

* add LogarithmicNumbers

* export basemeasure_sequence

* update superpose

* fix logdensity_rel

* remove FIXME (it's fixed!!)

* logdensityof(d::Density, x)

* simplify insupport(::Lebesgue, ::Real)

* clean up

* assume insupport yields Bool

* change logdensity_rel fall-through to warning and return NaN

* update logdensity_rel

* drop old code

* fix warning

* export logdensity_rel

* logdensity_def(μ::Dirac, ν::Dirac, x)

* logdensity_def methods

* drop `static`

* ]add StatsFuns

* Fixing up superposition

* [compat] entries

* trying to speed things up

* bugfixes

* logdensity_rel tests

* logdensity_rel tests

* drop qualifier, and add a test

* more tests

* type constraint in "logdensityof(μ::AbstractMeasure, x)" (was piracy, oops)

* add some docs

* docs

* docs

* typo

* moar speed

* don't export Test

* some more updates

* logdensity_rel for products

* `kleisli` docs

* update instance_type

* instance_type => Core.Typeof

* `powermeasure` bug fix

* fix logdensity_rel bug

* get `commonbase` to take x type into account

* test powers

* commonbase docstring

* deprecate instance_type

* avoid breakage

* switch || terms

* @ifelse macro

* simplify logdensity_rel

* give up on this @ifelse business

* bump version

* working on likelihoods

* update likelihood

* powerweightedmeasure

* powerweighted update

* more powerweighted methods

* bugfix

* dropFactoredBase

* drop FactoredBase

* (::ProductMeasure) | constraint

* update conditional measure

* update Dirac

* move conditional.jl down in the `include`s

* Kleisli => TransitionKernel

* simplify logdensity_def(::PowerMeasure, x)

* rename kleisli.jl to kernel.jl

* update Dirac tests

* update Half

* get tests passing

* update kernel

* Update Project.toml

* no call-site inlining

* restrict single-arg `kernel` to <:ParameterizedMeasure

* export log_likelihood_ratio

* Drop DensityKind(::Likelihood), at least for now

* isfinite(x) instead of x>-Inf

* add `condition` constructor

* EOF newline

* simplify logdensity_def for power measures

* finishing up

Author: cscherrer

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)
@@ -51,11 +52,9 @@ abstract type AbstractMeasure end
 
 using Static: @constprop
 
-function Pretty.tile(d::M) where {M<:AbstractMeasure}
+function Pretty.quoteof(d::M) where {M<:AbstractMeasure}
     the_names = fieldnames(typeof(d))
-    result = Pretty.literal(repr(M))
-    isempty(the_names) && return result * Pretty.literal("()")
-    Pretty.list_layout(Pretty.tile.([getfield(d, n) for n in the_names]); prefix=result)
+    :($M($([getfield(d, n) for n in the_names]...)))
 end
 
 @inline DensityKind(::AbstractMeasure) = HasDensity()
@@ -97,23 +96,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 +120,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

@@ -23,6 +23,24 @@ the notion of normalization may not even make sense.
 Because of this, this interface is not yet stable, and users should expect
 upcoming changes.
 """
-Base.:|(μ::AbstractMeasure, constraint) = ConditionalMeasure(μ, constraint)
+Base.:|(μ::AbstractMeasure, constraint) = condition(μ, constraint)
 
-@inline basemeasure(cm::ConditionalMeasure) = basemeasure(cm.parent) | cm.constraint
+condition(μ, constraint) = ConditionalMeasure(μ, 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
+
+function Pretty.tile(d::ConditionalMeasure)
+    Pretty.pair_layout(Pretty.tile(d.parent), Pretty.tile(d.constraint), sep=" | ")
+end

src/combinators/factoredbase.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
+testvalue(::Half) = 1.0

src/combinators/half.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 ℝ.
 
@@ -89,7 +95,8 @@ Finally, let's return to the expression for Bayes's Law,
 
 The product on the right side is computed pointwise. To work with this in
 MeasureBase, we have a "pointwise product" `⊙`, which takes a measure and a
-likelihood, and returns a new measure, that is, the unnormalized posterior that has density ``P(θ) P(x|θ)`` with respect to the base measure of the prior.
+likelihood, and returns a new measure, that is, the unnormalized posterior that
+has density ``P(θ) P(x|θ)`` with respect to the base measure of the prior.
 
 For example, say we have
 
@@ -109,14 +116,11 @@ 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
-@inline DensityKind(::Likelihood) = IsDensity()
-
 function Pretty.quoteof(ℓ::Likelihood)
     k = Pretty.quoteof(ℓ.k)
     x = Pretty.quoteof(ℓ.x)
@@ -128,10 +132,48 @@ 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::AbstractTransitionKernel, x; constraints...)
+    likelihood(k::AbstractTransitionKernel, x, constraints::NamedTuple)
+
+A likelihood is *not* a measure. Rather, a likelihood acts on a measure, through
+the "pointwise product" `⊙`, yielding another measure.
+"""
+function likelihood end
+
+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)
+
+export log_likelihood_ratio
+
+"""
+    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/likelihood.jl

@@ -1,40 +1,34 @@
 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
 
-function Base.show(io::IO, μ::PointwiseProductMeasure)
-    io = IOContext(io, :compact => true)
-    print(io, μ.measure, " ⊙ ", μ.likelihood)
-end
 
-function Base.show_unquoted(io::IO, μ::PointwiseProductMeasure, indent::Int, prec::Int)
-    io = IOContext(io, :compact => true)
-    if Base.operator_precedence(:*) ≤ prec
-        print(io, "(")
-        show(io, μ)
-        print(io, ")")
-    else
-        show(io, μ)
-    end
-    return nothing
+
+iterate(p::PointwiseProductMeasure, i=1) = iterate((p.prior, p.likelihood), i)
+
+function Pretty.tile(d::PointwiseProductMeasure)
+    Pretty.pair_layout(Pretty.tile(d.prior), Pretty.tile(d.likelihood), sep=" ⊙ ")
 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
 
-basemeasure(d::PointwiseProductMeasure) = @inbounds basemeasure(d.measure)
+@inbounds function insupport(d::PointwiseProductMeasure, p) 
+    μ, ℓ = d
+    insupport(μ, p) && insupport(ℓ.k(p), ℓ.x)
+end
+
+basemeasure(d::PointwiseProductMeasure, x) = d.prior
+
+basemeasure(d::PointwiseProductMeasure) = basemeasure(d.prior)

src/combinators/pointwise.jl

@@ -65,26 +65,16 @@ 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}
-    quote
-        $(Expr(:meta, :inline))
-        ℓ = 0.0
-        parent = d.parent
-        @inbounds @simd for j in 1:$N
-            Δℓ = logdensity_def(parent, x[j])
-            ℓ += Δℓ
-        end
-        ℓ
+@inline function logdensity_def(d::PowerMeasure{M, Tuple{Base.OneTo{StaticInt{N}}}}, x) where {M,N}
+    parent = d.parent
+    sum(1:N) do j
+        @inbounds logdensity_def(parent, x[j])
     end
 end
 

src/combinators/power.jl

@@ -0,0 +1,37 @@
+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
+
+function Pretty.tile(d::PowerWeightedMeasure)
+    Pretty.pair_layout(Pretty.tile(d.parent), Pretty.tile(d.exponent), sep=" ↑ ")
+end