Dev (#41)

* some refactoring

* pretty printing

* some updates for resettablerng

* deepcopy on iterate

* Move ResettableRNG to MeasureTheory

* Make Pretty a const

* drop extra spaces

* formatting

* formatting

* formatting

* add DensityInterface dependency

* faster `rootmeasure`

* updates

* primitive measure docs

* get tests passing

* cleaning up

* formatting

* 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

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.6.1"
+version = "0.7.1"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
@@ -13,12 +13,12 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 LogarithmicNumbers = "aa2f6b4e-9042-5d33-9679-40d3a6b85899"
 MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
+NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 PrettyPrinting = "54e16d92-306c-5ea0-a30b-337be88ac337"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Tricks = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775"
 
@@ -31,10 +31,10 @@ IfElse = "0.1"
 LogExpFunctions = "0.3"
 LogarithmicNumbers = "1"
 MappedArrays = "0.4"
+NaNMath = "0.3, 1"
 PrettyPrinting = "0.3, 0.4"
 Reexport = "1"
 Static = "0.5, 0.6"
-StatsFuns = "0.9"
 Tricks = "0.1"
 julia = "1.3"
 

src/MeasureBase.jl

@@ -64,8 +64,11 @@ gentype(μ::AbstractMeasure) = typeof(testvalue(μ))
 
 # gentype(μ::AbstractMeasure) = gentype(basemeasure(μ))
 
+using NaNMath
 using LogExpFunctions: logsumexp
 
+@deprecate instance_type(x) Core.Typeof(x) false
+
 """
 `logdensity_def` is the standard way to define a log-density for a new measure.
 Note that this definition does not include checking for membership in the

src/combinators/product.jl

@@ -65,6 +65,10 @@ struct ProductMeasure{M} <: AbstractProductMeasure
     marginals::M
 end
 
+@inline function logdensity_rel(μ::ProductMeasure, ν::ProductMeasure, x)
+    mapreduce(logdensity_rel, +, marginals(μ), marginals(ν), x)
+end
+
 function Pretty.tile(d::ProductMeasure{T}) where {T<:Tuple}
     Pretty.list_layout(Pretty.tile.([marginals(d)...]), sep=" ⊗ ")
 end

src/combinators/smart-constructors.jl

@@ -14,11 +14,16 @@ end
 ###############################################################################
 # PowerMeaure
 
-function powermeasure(μ::WeightedMeasure, dims::NTuple{N,I}) where {N,I}
+function powermeasure(μ::WeightedMeasure, dims::NTuple{N,I}) where {N,I<:AbstractArray}
     k = mapreduce(length, *, dims) * μ.logweight
     return weightedmeasure(k, μ.base^dims)
 end
 
+function powermeasure(μ::WeightedMeasure, dims::NTuple{N,I}) where {N,I}
+    k = prod(dims) * μ.logweight
+    return weightedmeasure(k, μ.base^dims)
+end
+
 ###############################################################################
 # ProductMeasure
 
@@ -84,7 +89,7 @@ function superpose(μ::T, ν::T) where {T}
     end
 end
 
-function superpose(μ::AbstractMeasure, ν::AbstractMeasure)
+function superpose(μ, ν)
     components = (μ, ν)
     superpose(components)
 end

src/combinators/superpose.jl

@@ -54,8 +54,7 @@ testvalue(μ::SuperpositionMeasure) = testvalue(first(μ.components))
 # end
 
 function Base.:+(μ::AbstractMeasure, ν::AbstractMeasure)
-    components = (μ, ν)
-    superpose(components)
+    superpose(μ, ν)
 end
 
 using LogarithmicNumbers
@@ -64,8 +63,10 @@ oneplus(x::ULogarithmic) = exp(ULogarithmic, log1pexp(x.log))
 
 @inline function density_def(s::SuperpositionMeasure{Tuple{A,B}}, x) where {A,B}
     (μ, ν) = s.components
+
     insupport(μ, x) || return exp(ULogarithmic, logdensity_def(ν, x))
     insupport(ν, x) || return exp(ULogarithmic, logdensity_def(μ, x))
+  
     α = basemeasure(μ)
     β = basemeasure(ν)
     dμ_dα = exp(ULogarithmic, logdensity_def(μ, x))
@@ -75,7 +76,7 @@ oneplus(x::ULogarithmic) = exp(ULogarithmic, log1pexp(x.log))
     return dμ_dα / oneplus(dβ_dα) + dν_dβ / oneplus(dα_dβ)
 end
 
-using StatsFuns
+using LogExpFunctions
 
 @inline function logdensity_def(μ::T, ν::T, x::Any) where T<:(SuperpositionMeasure{Tuple{A, B}} where {A, B})
     if μ === ν
@@ -87,6 +88,7 @@ end
 
 @inline function logdensity_def(s::SuperpositionMeasure{Tuple{A,B}}, β, x) where {A,B}
     (μ, ν) = s.components
+
     insupport(μ, x) || return logdensity_rel(ν, β, x)
     insupport(ν, x) || return logdensity_rel(μ, β, x)
     return logaddexp(logdensity_rel(μ, β, x), logdensity_rel(ν, β, x))
@@ -105,7 +107,7 @@ end
 
 @inline logdensity_def(s::SuperpositionMeasure, x) = log(density_def(s, x))
 
-basemeasure(μ::SuperpositionMeasure) = superpose(map(basemeasure, μ.components))
+basemeasure(μ::SuperpositionMeasure) = superpose(map(basemeasure, μ.components)...)
 
 # TODO: Fix `rand` method (this one is wrong)
 # function Base.rand(μ::SuperpositionMeasure{X,N}) where {X,N}

src/density.jl

@@ -174,12 +174,11 @@ known to be in the support of both, it can be more efficient to call
 `unsafe_logdensity_rel`. 
 """
 @inline function logdensity_rel(μ::M, ν::N, x::X) where {M,N,X}
-    T = unstatic(float(promote_type(return_type(logdensity_def, (μ, x)), return_type(logdensity_def, (ν, x)))))
-    insupport(μ, x) || begin
-        insupport(ν, x) || return convert(T, NaN)
-        return convert(T, -Inf)
-    end
-    insupport(ν, x) || return convert(T, Inf)
+    T = unstatic(promote_type(return_type(logdensity_def, (μ, x)), return_type(logdensity_def, (ν, x))))
+    inμ = insupport(μ, x)
+    inν = insupport(ν, x)
+    inμ || return convert(T, ifelse(inν, -Inf, NaN))
+    inν || return convert(T, Inf)
 
     return unsafe_logdensity_rel(μ, ν, x)
 end
@@ -198,7 +197,24 @@ See also `logdensity_rel`.
     end
     μs = basemeasure_sequence(μ)
     νs = basemeasure_sequence(ν)
-    return _logdensity_rel(μs, νs, x)
+    cb = commonbase(μs, νs, X)
+    # _logdensity_rel(μ, ν)
+    isnothing(cb) && begin
+        μ = μs[end]
+        ν = νs[end]
+        @warn """
+        No common base measure for
+            $μ
+        and
+            $ν
+
+        Returning a relative log-density of NaN. If this is incorrect, add a
+        three-argument method
+            logdensity_def($μ, $ν, x)
+        """
+        return NaN
+    end
+    return _logdensity_rel(μs, νs, cb, x)
 end
 
 # Note that this method assumes `μ` and `ν` to have the same type
@@ -210,44 +226,29 @@ function logdensity_def(μ::T, ν::T, x) where {T}
     end
 end
 
-@generated function _logdensity_rel(μs::Tμ, νs::Tν, x::X)  where {Tμ, Tν, X}
+@generated function _logdensity_rel(μs::Tμ, νs::Tν, ::Tuple{StaticInt{M},StaticInt{N}}, x::X)  where {Tμ, Tν,M,N,X}
     sμ = schema(Tμ)
     sν = schema(Tν)
 
     q = quote 
         $(Expr(:meta, :inline))
+        ℓ = logdensity_def(μs[$M], νs[$N], x)
     end
 
-    for it in Iterators.product(enumerate(sμ), enumerate(sν))
-        ((nμ, μtype), (nν, νtype)) = it
-        if static_hasmethod(logdensity_def, Tuple{μtype, νtype, X})
-            push!(q.args, :(ℓ = logdensity_def(μs[$nμ], νs[$nν], x)))
-            for i in 1:nμ-1
-                push!(q.args, :(ℓ += logdensity_def(μs[$i], x)))
-            end
-            for j in 1:nν-1
-                push!(q.args, :(ℓ -= logdensity_def(νs[$j], x)))
-            end
-
-            return q
-        end
+    for i in 1:M-1
+        push!(q.args, :(Δℓ = logdensity_def(μs[$i], x)))
+        # push!(q.args, :(println("Adding", Δℓ)))
+        push!(q.args, :(ℓ += Δℓ))
     end
 
-    return quote
-        μ = μs[end]
-        ν = νs[end]
-        @warn """
-        No common base measure for
-            $μ
-        and
-            $ν
-
-        Returning a relative log-density of NaN. If this is incorrect, add a
-        three-argument method
-            logdensity_def(μ, ν, x)
-        """
-        NaN
+    for j in 1:N-1
+        push!(q.args, :(Δℓ = logdensity_def(νs[$j], x)))
+        # push!(q.args, :(println("Subtracting", Δℓ)))
+        push!(q.args, :(ℓ -= Δℓ))
     end
+
+    push!(q.args, :(return ℓ))
+    return q
 end
 
 export densityof

src/interface.jl

@@ -5,12 +5,14 @@ using Reexport
 @reexport using MeasureBase
 
 using MeasureBase:basemeasure_depth, proxy
-using MeasureBase: insupport
+using MeasureBase: insupport, basemeasure_sequence, commonbase
 
 export test_interface
 export basemeasure_depth
 export proxy
 export insupport
+export basemeasure_sequence
+export commonbase
 
 using Test
 

src/kleisli.jl

@@ -15,17 +15,32 @@ end
     kleisli(f, M)
     kleisli((f1, f2, ...), M)
 
-A kleisli `κ = kleisli(f, m)` returns a wrapper around
-a function `f` giving the parameters for a measure of type `M`,
-such that `κ(x) = M(f(x)...)`
-respective `κ(x) = M(f1(x), f2(x), ...)`
+`kleisli` was originally called `kernel`, as in a *Markov kernel*. Such a kernel
+can be considered to map each value in its domain to a probability measure.
+
+In the context of MeasureTheory, the codomain is not required to be a
+*probability* measure; any measure will do. This makes "Markov" not really fit,
+since the map need not be Markovian.
+
+This leaves us with "kernel", which can mean too wide a range of things to be
+useful in such a general context as measure theory. See for example
+https://github.com/JuliaGaussianProcesses/KernelFunctions.jl for one common use
+of this term.
+
+We solve this problem by changing to a term from an even more general context.
+In category theory, a *Kleisli arrow* is a function taking monadic values.
+Since measures comprise a monad, our use is a special case of this.
+
+A kleisli `κ = kleisli(f, m)` returns a wrapper around a function `f` giving the
+parameters for a measure of type `M`, such that `κ(x) = M(f(x)...)` respective
+`κ(x) = M(f1(x), f2(x), ...)`
 
 If the argument is a named tuple `(;a=f1, b=f1)`, `κ(x)` is defined as
 `M(;a=f(x),b=g(x))`.
 
 # Reference
 
-* https://en.wikipedia.org/wiki/Markov_kleisli
+* https://en.wikipedia.org/wiki/Markov_kernel
 """
 function kleisli end
 

src/primitives/counting.jl

@@ -9,7 +9,7 @@ insupport(::CountingMeasure, x) = true
 struct Counting{T} <: AbstractMeasure
     support::T
 
-    Counting(supp) = new{instance_type(supp)}(supp)
+    Counting(supp) = new{Core.Typeof(supp)}(supp)
 end
 
 function logdensity_def(μ::Counting, x)

src/utils.jl

@@ -52,11 +52,6 @@ end
     return getfield(T, :instance)::T
 end
 
-# See https://github.com/cscherrer/KeywordCalls.jl/issues/22
-@inline instance_type(f::F) where {F} = F
-@inline instance_type(T::UnionAll) = Type{T}
-@inline instance_type(T::DataType) = Type{T}
-
 export basemeasure_depth
 
 @inline function basemeasure_depth(μ::M) where {M}
@@ -70,9 +65,12 @@ export basemeasure_depth
     return static(10)
 end
 
+"""
+    basemeasure_sequence(m)
 
-export basemeasure_sequence
-
+Construct the longest `Tuple` starting with `m` having each term as the base
+measure of the previous term, and with no repeated entries.
+"""
 @inline function basemeasure_sequence(μ::M) where {M}
     b_1 = μ
     done = false
@@ -86,21 +84,32 @@ export basemeasure_sequence
     return filter(!isnothing, Base.Cartesian.@ntuple 10 b)
 end
 
-# @inline function basemeasure_depth(μ::M) where {M}
-#     return basemeasure_depth(μ, basemeasure(μ), static(0))
-# end
-
-# @inline function basemeasure_depth(μ::M, β::M, s::StaticInt{N}) where {M,N}
-#     s
-# end
-
-# @generated function basemeasure_depth(μ::M, β::B, ::StaticInt{N}) where {M,B,N}
-#     s = Expr(:call, Expr(:curly, :StaticInt, N + 1))
-#     quote
-#         $(Expr(:meta, :inline))
-#         basemeasure_depth(β, basemeasure(β), $s)
-#     end
-# end
+commonbase(μ, ν) = commonbase(μ, ν, Any) 
+
+"""
+    commonbase(μ, ν, T) -> Tuple{StaticInt{i}, StaticInt{j}}
+
+Find minimal (with respect to their sum) `i` and `j` such that there is a method
+
+    logdensity_def(basemeasure_sequence(μ)[i], basemeasure_sequence(ν)[j], ::T)
+
+This is used in `logdensity_rel` to help make that function efficient.
+"""
+@inline function commonbase(μ, ν, ::Type{T}) where {T}
+    return commonbase(basemeasure_sequence(μ), basemeasure_sequence(ν), T)
+end
+
+@generated function commonbase(μ::M, ν::N, ::Type{T}) where {M<:Tuple,N<:Tuple,T}
+    m = schema(M)
+    n = schema(N)
+
+    sols = Iterators.filter(((i,j),) ->  static_hasmethod(logdensity_def, Tuple{m[i], n[j], T}), Iterators.product(1:length(m), 1:length(n))) 
+    isempty(sols) && return :(nothing)
+    minsol = static.(argmin(((i,j),) -> i+j, sols))
+    quote
+        $minsol
+    end
+end
 
 mymap(f, gen::Base.Generator) = mymap(f ∘ gen.f, gen.iter)
 mymap(f, inds...) = Iterators.map(f, inds...)