Dev (#64)

- `ScaledMeasure` is now `WeightedMeasure` to avoid confusion with scale parameters
- Added an `Exp` constructor for lazy exponentiation
- Poisson and Binomial distributions
- Added `ForArray` and `ForGenerator` type aliases
- `PowerMeasure` is now in terms of `ProductMeasure`, to avoid some performance hits due to compiler limitations
- Added `@half` macro to build measures like `HalfNormal` and `HalfCauchy`
- Added some tests

Project.toml

@@ -1,14 +1,13 @@
 name = "MeasureTheory"
 uuid = "eadaa1a4-d27c-401d-8699-e962e1bbc33b"
 authors = ["Chad Scherrer <[email protected]> and contributors"]
-version = "0.2.3"
+version = "0.3.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DistributionsAD = "ced4e74d-a319-5a8a-b0ac-84af2272839c"
-FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 KeywordDispatch = "5888135b-5456-5c80-a1b6-c91ef8180460"
@@ -26,18 +25,17 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
 
 [compat]
-ArrayInterface = "2"
+ArrayInterface = "2,3"
 Bijectors = "0.8"
 Distributions = "0.23, 0.24"
 DistributionsAD = "0.6"
-FillArrays = "0.11"
 InfiniteArrays = "0.9"
 IntervalSets = "0.5"
 KeywordDispatch = "0.3"
 LazyArrays = "0.20"
-MLStyle = "0.4"
 MacroTools = "0.5"
 MappedArrays = "0.3"
+MLStyle = "0.4"
 NamedTupleTools = "0.13"
 NestedTuples = "0.1"
 SpecialFunctions = "0.10, 1"

src/MeasureTheory.jl

@@ -13,6 +13,9 @@ export ≪
 export sampletype
 
 export AbstractMeasure
+using InfiniteArrays: ∞
+
+export ∞
 
 abstract type AbstractMeasure end
 
@@ -29,13 +32,14 @@ Methods for computing density relative to other measures will be
 """
 function logdensity end
 
+include("exp.jl")
 include("domains.jl")
 include("utils.jl")
 include("absolutecontinuity.jl")
 include("basemeasures.jl")
 include("parameterized.jl")
 include("macros.jl")
-include("combinators/scale.jl")
+include("combinators/weighted.jl")
 include("combinators/superpose.jl")
 include("combinators/product.jl")
 include("combinators/for.jl")
@@ -54,6 +58,8 @@ include("probability/exponential.jl")
 include("probability/mvnormal.jl")
 include("probability/inverse-gamma.jl")
 include("probability/bernoulli.jl")
+include("probability/poisson.jl")
+include("probability/binomial.jl")
 include("density.jl")
 include("pushforward.jl")
 include("kernel.jl")

src/basemeasures/counting.jl

@@ -1,7 +1,8 @@
-struct CountingMeasure{X} <: AbstractMeasure end
+export CountingMeasure
 
+struct CountingMeasure{X} <: AbstractMeasure end
 
-function Base.show(io::IO, μ::CountingMeasure{X}) where X
+function Base.show(io::IO, μ::CountingMeasure{X}) where {X}
     io = IOContext(io, :compact => true)
     print(io, "CountingMeasure(", X, ")")
 end
@@ -12,11 +13,13 @@ basemeasure(μ::CountingMeasure) = μ
 
 isprimitive(::CountingMeasure) = true
 
-sampletype(::CountingMeasure{ℝ}) = Float64
-sampletype(::CountingMeasure{ℝ₊}) = Float64
-sampletype(::CountingMeasure{𝕀}) = Float64
+# sampletype(::CountingMeasure{ℝ}) = Float64
+# sampletype(::CountingMeasure{ℝ₊}) = Float64
+# sampletype(::CountingMeasure{𝕀}) = Float64
+
+sampletype(::CountingMeasure{IntegerRange{lo,hi}}) where {lo, hi} = Int
+
 
-sampletype(::CountingMeasure{IntegerRange}) = Int
 
 logdensity(::CountingMeasure, x) = zero(float(x))
 

src/combinators/for.jl

@@ -3,15 +3,17 @@ export For
 using Random
 import Base
 
-For{D,N,T,F} = ProductMeasure{ReadonlyMappedArray{D, N, T, F}} 
+# ForArray
 
-function Base.show(io::IO, d::For{D,N,T,F}) where {D,N,T,F}
+ForArray{D,N,T,F} = ProductMeasure{ReadonlyMappedArray{D, N, T, F}} 
+
+function Base.show(io::IO, d::ForArray{D,N,T,F}) where {D,N,T,F}
     print(io, "For(")
     print(io, d.data.f, ", ")
     print(io, d.data.data, ")")
 end
 
-function Base.show(io::IO, d::For{D,N,T,F}) where {D,N,T <: CartesianIndices,F}
+function Base.show(io::IO, d::ForArray{D,N,T,F}) where {D,N,T <: CartesianIndices,F}
     print(io, "For(")
     print(io, d.data.f, ", ")
     join(io, size(d.data), ", ")
@@ -22,11 +24,11 @@ For(f, dims::AbstractArray...) = ProductMeasure(mappedarray(f, dims...))
 
 For(f, dims::Int...) = ProductMeasure(mappedarray(i -> f(Tuple(i)...), CartesianIndices(dims))) 
 
-function Base.eltype(::For{D,N,T,F}) where {D,N,T,F}
+function Base.eltype(::ForArray{D,N,T,F}) where {D,N,T,F}
     return eltype(D)
 end
 
-basemeasure(μ::For) = @inbounds basemeasure(μ.data[1])^size(μ.data)
+basemeasure(μ::ForArray) = @inbounds basemeasure(μ.data[1])^size(μ.data)
 
 # """
 #     indexstyle(a::AbstractArray, b::AbstractArray)
@@ -42,21 +44,43 @@ basemeasure(μ::For) = @inbounds basemeasure(μ.data[1])^size(μ.data)
 #     return IndexCartesian()
 # end
 
-# function Base.rand(rng::AbstractRNG, μ::For{D,N,T,F}) where {F,T<:AbstractArray,D,X}
+# function Base.rand(rng::AbstractRNG, μ::ForArray{D,N,T,F}) where {F,T<:AbstractArray,D,X}
 #     s = size(μ.θ)
 #     x = Array{X,length(s)}(undef, s...)
 #     rand!(rng, x, μ)
 # end
 
-# function logdensity(μ::For{D,N,T,F}, x)
+# function logdensity(μ::ForArray{D,N,T,F}, x)
 #     getℓ(θⱼ, xⱼ) = logdensity(μ.f(θⱼ), xⱼ)
 #     ℓ = mappedarray(getℓ, μ.θ, x)
 #     _logdensity(μ, x, indexstyle(μ.θ, x), result_type)
 # end
 
-# function _logdensity(μ::For{D,N,T,F}, x, ::IndexLinear, ::Type{R}) where {R<:AbstractFloat}
+# function _logdensity(μ::ForArray{D,N,T,F}, x, ::IndexLinear, ::Type{R}) where {R<:AbstractFloat}
 #     ℓ = zero(R)
 #     μ.f(μ.θ)
 # end
 
-# function basemeasure(μ::For{D,N,T,F}) where {F,T<:AbstractArray,D,X}
+# function basemeasure(μ::ForArray{D,N,T,F}) where {F,T<:AbstractArray,D,X}
+
+# ForGenerator
+
+ForGenerator{G} = ProductMeasure{G} where {G <: Base.Generator}
+
+For(f, dims::Base.Generator) = ProductMeasure(Base.Generator(f ∘ dims.f, dims.iter))
+
+sampletype(::ForGenerator) = Base.Generator
+
+function Base.rand(rng::AbstractRNG, T::Type{<:Base.Generator}, d::ForGenerator)
+    elT = sampletype(first(d.data))
+    f(x) = rand(rng, elT, x)
+    Base.Generator(f ∘ d.data.f, d.data.iter)
+end
+
+function Base.rand(rng::AbstractRNG, ::Type{<:Array{T}}, d::ForGenerator) where {T}
+    return collect(T, rand(rng, d))
+end
+
+function MeasureTheory.logdensity(d::ForGenerator, x)
+    sum((logdensity(dj, xj) for (dj, xj) in zip(d.data, x)))
+end

src/combinators/power.jl

@@ -1,42 +1,41 @@
 import Base
-import FillArrays
 
-"""
-A power measure is a product of a measure with itself. The number of elements in
-the product determines the dimensionality of the resulting support.
+# """
+# A power measure is a product of a measure with itself. The number of elements in
+# the product determines the dimensionality of the resulting support.
 
-Note that power measures are only well-defined for integer powers.
+# Note that power measures are only well-defined for integer powers.
 
-The nth power of a measure μ can be written μ^x.
-"""
-PowerMeasure{M,N,D} = ProductMeasure{Fill{M,N,D}}
+# The nth power of a measure μ can be written μ^x.
+# """
+# PowerMeasure{M,N,D} = ProductMeasure{Fill{M,N,D}}
 
-function Base.show(io::IO, μ::PowerMeasure)
-    io = IOContext(io, :compact => true)
-    print(io, μ.data.value, " ^ ", size(μ.data))
-end
+# function Base.show(io::IO, μ::PowerMeasure)
+#     io = IOContext(io, :compact => true)
+#     print(io, μ.data.value, " ^ ", size(μ.data))
+# end
 
-function Base.show_unquoted(io::IO, μ::PowerMeasure{M,N,D}, indent::Int, prec::Int) where {M,N,D}
-    io = IOContext(io, :compact => true)
-    if Base.operator_precedence(:^) ≤ prec
-        print(io, "(")
-        show(io, μ.data.value)
-        print(io, ")")
-    else
-        show(io, size(μ.data))
-    end
-    return nothing
-end
+# function Base.show_unquoted(io::IO, μ::PowerMeasure{M,N,D}, indent::Int, prec::Int) where {M,N,D}
+#     io = IOContext(io, :compact => true)
+#     if Base.operator_precedence(:^) ≤ prec
+#         print(io, "(")
+#         show(io, μ.data.value)
+#         print(io, ")")
+#     else
+#         show(io, size(μ.data))
+#     end
+#     return nothing
+# end
 
-Base.:^(μ::AbstractMeasure, n::Integer) = μ ^ (n,)
+Base.:^(μ::AbstractMeasure, n::Integer) = For(_ -> μ, 1:n)
 
-Base.:^(μ::AbstractMeasure, size::NTuple{N,I}) where {N, I <: Integer} = ProductMeasure(Fill(μ, size))
+Base.:^(μ::AbstractMeasure, size::NTuple{N,I}) where {N, I <: Integer} = For(_ -> μ, CartesianIndices(size))
 
-sampletype(d::PowerMeasure{M,N}) where {M,N} = @inbounds Array{sampletype(first(d.data)), N}
+# sampletype(d::PowerMeasure{M,N}) where {M,N} = @inbounds Array{sampletype(first(d.data)), N}
 
 function Base.:^(μ::WeightedMeasure, n::NTuple{N,Int}) where {N}
     k = prod(n) * μ.logweight
     return WeightedMeasure(k, μ.base^n)
 end
 
-basemeasure(μ::PowerMeasure) = @inbounds basemeasure(first(μ.data))^size(μ.data)
+# basemeasure(μ::PowerMeasure) = @inbounds basemeasure(first(μ.data))^size(μ.data)

src/combinators/product.jl

@@ -1,7 +1,6 @@
 export ProductMeasure
 
 using MappedArrays
-using FillArrays
 using Base: @propagate_inbounds
 
 struct ProductMeasure{T} <: AbstractMeasure

src/combinators/superpose.jl

@@ -10,36 +10,38 @@ measures need not be normalized) requires no scaling.
 
 The superposition of two measures μ and ν can be more concisely written as μ + ν.
 """
-struct SuperpositionMeasure{X,NT} <: AbstractMeasure
+struct SuperpositionMeasure{NT} <: AbstractMeasure
     components :: NT   
 end
 
-SuperpositionMeasure(ms :: AbstractMeasure...) where {X} = SuperpositionMeasure{X,length(ms)}(ms)
+# SuperpositionMeasure(ms :: AbstractMeasure...) = SuperpositionMeasure{X,length(ms)}(ms)
 
 # SuperpositionMeasure(m::NTuple{N, Measure{X}}) where {N,X} = SuperpositionMeasure(m...)
 
 # Base.length(::SuperpositionMeasure{X,N}) where {X,N} = N
 
-function Base.:+(μ::SuperpositionMeasure{X,N1}, ν::SuperpositionMeasure{X,N2}) where {X,N1,N2}
-    components = (μ.components..., ν.components...)
-    SuperpositionMeasure{X, N1+N2}(components)
-end
+# function Base.:+(μ::SuperpositionMeasure{N1}, ν::SuperpositionMeasure{N2}) where {N1,N2}
+#     components = (μ.components..., ν.components...)
+#     SuperpositionMeasure{X, N1+N2}(components)
+# end
 
-function Base.:+(μ::AbstractMeasure, ν::SuperpositionMeasure{X,N}) where {X,N}
-    components = (μ, ν.components...)
-    SuperpositionMeasure{X,N+1}(components)
-end
+# function Base.:+(μ::AbstractMeasure, ν::SuperpositionMeasure{X,N}) where {X,N}
+#     components = (μ, ν.components...)
+#     SuperpositionMeasure{X,N+1}(components)
+# end
 
-function Base.:+(μ::SuperpositionMeasure{X,N}, ν::AbstractMeasure) where {X,N}
-    components = (μ.components..., ν)
-    SuperpositionMeasure{X,N+1}(components)
-end
+# function Base.:+(μ::SuperpositionMeasure{X,N}, ν::AbstractMeasure) where {X,N}
+#     components = (μ.components..., ν)
+#     SuperpositionMeasure{X,N+1}(components)
+# end
 
-function Base.:+(μ::AbstractMeasure, ν::AbstractMeasure) where {X}
+function Base.:+(μ::AbstractMeasure, ν::AbstractMeasure)
     components = (μ, ν)
-    SuperpositionMeasure{X,2}(components)
+    SuperpositionMeasure(components)
 end
 
+logdensity(μ::SuperpositionMeasure, x) = logsumexp((logdensity(m,x) for m in μ.components))
+
 # TODO: Fix `rand` method (this one is wrong)
 # function Base.rand(μ::SuperpositionMeasure{X,N}) where {X,N}
 #     return rand(rand(μ.components))

src/combinators/scale.jl → src/combinators/weighted.jl

@@ -1,3 +1,5 @@
+export WeightedMeasure
+
 """
     struct WeightedMeasure{R,M} <: AbstractMeasure
         logweight :: R
@@ -38,7 +40,8 @@ function Base.:*(k::T, m::AbstractMeasure) where {T <: Number}
     if hasmethod(iszero, (T,)) && iszero(k)
         return TrivialMeasure
     else
-        return WeightedMeasure{typeof(k), typeof(m)}(log(k),m)
+        logk = log(k)
+        return WeightedMeasure{typeof(logk), typeof(m)}(logk,m)
     end
 end
 

src/density.jl

@@ -54,6 +54,11 @@ function Base.show(io::IO, μ::DensityMeasure{F,B,Val{L}}) where {F,B,L}
     print(io, "; log = ", L, ")")
 end
 
+function Base.rand(rng::AbstractRNG, T::Type, d::DensityMeasure)
+    x = rand(d.base)
+    WeightedMeasure(d.f(x), Dirac(x))
+end
+
 basemeasure(μ::DensityMeasure) = μ.base
 
 logdensity(μ::DensityMeasure{F,B,Val{true}}, x) where {F,B} = μ.f(x)
@@ -68,6 +73,8 @@ Define a new measure in terms of a density `f` over some measure `base`. If
 """
 ∫(f, base::AbstractMeasure; log=true) = DensityMeasure(f, base, Val(log))
 
+∫(μ::AbstractMeasure, base::AbstractMeasure; log=true) = ∫(𝒹(μ,base), base; log=log)
+
 # TODO: `density` and `logdensity` functions for `DensityMeasure`
 
 function logdensity(μ::AbstractMeasure, ν::AbstractMeasure, x)
@@ -80,3 +87,9 @@ function logdensity(μ::AbstractMeasure, ν::AbstractMeasure, x)
 end
 
 _logdensity(::Lebesgue{ℝ}, ::Lebesgue{ℝ}, x) = zero(float(x))
+
+export density
+
+density(μ::AbstractMeasure, ν::AbstractMeasure, x) = Exp(logdensity(μ, ν, x))
+
+density(μ::AbstractMeasure, x) = Exp(logdensity(μ, x))