Dev (#166)

* Dirichlet(k::Integer, α) = Dirichlet(Fill(α, k))

* export TransformVariables as TV

* drop redundant import

* 0.0 => zero(Float64)

* drop outdated Dists.logpdf

* update StudentT

* drop redundant import

* update Uniform

* bump MeasureBase version

* reworking beta

* small update to StudentT

* basemeasure for discrete Distributions

* using LogExpFunctions => import LogExpFunctions

* quoteof(::Chain)

* prettyprinting and chain-mucking

* Some refactoring for Markov chains

* import MeasureBase: ≪

* version bound for PrettyPrinting

* copy(rng) might change its type (e.g. GLOBAL_RNG)

* tests pass

* cleaning up

* more cleanup

* bump versions

* change import..as so it's 1.5 friendly

* bump version

Author: cscherrer

Project.toml

@@ -1,7 +1,7 @@
 name = "MeasureTheory"
 uuid = "eadaa1a4-d27c-401d-8699-e962e1bbc33b"
 authors = ["Chad Scherrer <[email protected]> and contributors"]
-version = "0.13.0"
+version = "0.13.1"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -11,6 +11,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DynamicIterators = "6c76993d-992e-5bf1-9e63-34920a5a5a38"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 KeywordCalls = "4d827475-d3e4-43d6-abe3-9688362ede9f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
@@ -21,6 +22,7 @@ MeasureBase = "fa1605e6-acd5-459c-a1e6-7e635759db14"
 NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
 NestedTuples = "a734d2a7-8d68-409b-9419-626914d4061d"
 PositiveFactorizations = "85a6dd25-e78a-55b7-8502-1745935b8125"
+PrettyPrinting = "54e16d92-306c-5ea0-a30b-337be88ac337"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -35,16 +37,17 @@ ConstructionBase = "1.3"
 Distributions = "0.25"
 DynamicIterators = "0.4"
 FillArrays = "0.12"
-InfiniteArrays = "0.11"
+InfiniteArrays = "0.11, 0.12"
 KeywordCalls = "0.2"
 LogExpFunctions = "0.3.3"
 MLStyle = "0.4"
 MacroTools = "0.5"
 MappedArrays = "0.4"
-MeasureBase = "0.4"
+MeasureBase = "0.5"
 NamedTupleTools = "0.13"
 NestedTuples = "0.3"
 PositiveFactorizations = "0.2"
+PrettyPrinting = "0.3"
 Reexport = "1"
 SpecialFunctions = "1"
 StatsFuns = "0.9"

src/MeasureTheory.jl

@@ -13,7 +13,7 @@ import Base
 import Distributions
 const Dists = Distributions
 
-# export TV
+export TV
 export ≪
 export sampletype
 export For
@@ -38,12 +38,17 @@ using ConstructionBase
 using Accessors
 using StatsFuns
 using SpecialFunctions
-using LogExpFunctions 
 
+import LogExpFunctions 
 import NamedTupleTools
 
 import MeasureBase: testvalue, logdensity, density, basemeasure, kernel, params, ∫
-import MeasureBase: affine, supportdim
+import MeasureBase: affine, supportdim, ≪
+using MeasureBase: constructor
+
+import PrettyPrinting
+
+const Pretty = PrettyPrinting
 
 import Base: rand
 
@@ -65,9 +70,6 @@ sampletype(μ::AbstractMeasure) = typeof(testvalue(μ))
 
 # sampletype(μ::AbstractMeasure) = sampletype(basemeasure(μ))
 
-
-
-
 import Distributions: logpdf, pdf
 
 export pdf, logpdf
@@ -94,6 +96,12 @@ const AFFINEPARS = [
     (:μ,)
 ]
 
+xlogy(x::Number, y::Number) = LogExpFunctions.xlogy(x, y)
+xlogy(x, y) = x * log(y)
+
+xlog1py(x::Number, y::Number) = LogExpFunctions.xlog1py(x, y)
+xlog1py(x, y) = x * log1p(y)
+
 
 include("const.jl")
 # include("traits.jl")
@@ -104,6 +112,9 @@ include("combinators/affine.jl")
 include("combinators/weighted.jl")
 include("combinators/product.jl")
 include("combinators/transforms.jl")
+
+include("resettable-rng.jl")
+include("realized.jl")
 include("combinators/chain.jl")
 
 include("distributions.jl")

src/combinators/chain.jl

@@ -1,18 +1,23 @@
 using ConcreteStructs
 using DynamicIterators
-using DynamicIterators: dub
+import DynamicIterators: dub, dyniterate, evolve
 using Base.Iterators: SizeUnknown, IsInfinite
-import DynamicIterators: dyniterate, evolve
 
 import MeasureBase: For
 
+
+
 export Chain
 
 @concrete terse struct Chain{K,M} <: AbstractMeasure
     κ::K
     μ::M
 end
 
+Pretty.quoteof(c::Chain) = :(Chain($(Pretty.quoteof(c.κ)), $(Pretty.quoteof(c.μ))))
+
+Base.length(::Chain) = ∞
+
 function basemeasure(mc::Chain)
     Chain(basemeasure ∘ mc.κ, basemeasure(mc.μ))
 end
@@ -39,73 +44,19 @@ dyniterate(E::Chain, ::Nothing) = dub(evolve(E))
 Base.iterate(E::Chain) = dyniterate(E, nothing)
 Base.iterate(E::Chain, value) = dyniterate(E, value)
 
-function DynamicIterators.dyniterate(r::Chain, (u,rng)::Sample)
-    μ = r.κ(u) 
-    u = rand(rng, μ)
-    return u, Sample(u, rng)
+function DynamicIterators.dyniterate(r::Chain, (x,rng)::Sample)
+    μ = r.κ(x) 
+    y = rand(rng, μ)
+    return (μ ↝ y), Sample(y, rng)
 end
+
 Base.IteratorSize(::Chain) = IsInfinite()
 Base.IteratorSize(::Type{Chain}) = IsInfinite()
 
 
-@concrete terse struct Realized{R,S,T} <: DynamicIterators.DynamicIterator
-    rng::ResettableRNG{R,S}
-    iter::T
-end
-
-Base.IteratorEltype(mc::Realized) = Base.HasEltype()
-
-function Base.eltype(::Type{Rz}) where {R,S,T,Rz <: Realized{R,S,T}}
-    eltype(T)
-end
-
-Base.length(r::Realized) = length(r.iter)
-
-Base.size(r::Realized) = size(r.iter)
-
-Base.IteratorSize(::Type{Rz}) where {R,S,T, Rz <: Realized{R,S,T}} = Base.IteratorSize(T)
-Base.IteratorSize(r::Rz) where {R,S,T, Rz <: Realized{R,S,T}} = Base.IteratorSize(r.iter)
-
-
-function Base.iterate(rv::Realized{R,S,T}) where {R,S,T}
-    if static_hasmethod(evolve, Tuple{T})
-        dyniterate(rv, nothing)
-    else
-        !isnothing(rv.rng.seed) && reset!(rv.rng)
-        μ,s = iterate(rv.iter)
-        x = rand(rv.rng, μ)
-        x,s
-    end
-end
-
-
-function Base.iterate(rv::Realized{R,S,T}, s) where {R,S,T}
-    if static_hasmethod(evolve, Tuple{T})
-        dyniterate(rv, s)
-    else
-        μs = iterate(rv.iter, s)
-        isnothing(μs) && return nothing
-        (μ,s) = μs
-        x = rand(rv.rng, μ)
-        return x,s
-    end
-end
-
-
-function dyniterate(rv::Realized, ::Nothing)
-    !isnothing(rv.rng.seed) && reset!(rv.rng)
-    μ = evolve(rv.iter)
-    x = rand(rv.rng, μ)
-    x, Sample(x, rv.rng)
-end
-function dyniterate(rv::Realized, u::Sample)
-    dyniterate(rv.iter, u)
-end
-
 function Base.rand(rng::AbstractRNG, T::Type, chain::Chain)
-    seed = rand(rng, UInt)
-    r = ResettableRNG(rng, seed)
-    return Realized(r, chain)
+    r = ResettableRNG(rng)
+    return RealizedSamples(r, chain)
 end
 
 ###############################################################################
@@ -115,38 +66,53 @@ end
 
 @concrete terse struct DynamicFor{T,K,S} <: AbstractMeasure
     κ ::K
-    sampler :: S        
+    iter :: S        
 end
 
-function DynamicFor(κ::K,sampler::S) where {K,S}
-    T = typeof(κ(first(sampler)))
-    DynamicFor{T,K,S}(κ,sampler)
+Pretty.quoteof(r::DynamicFor) = :(DynamicFor($(Pretty.quoteof(r.κ)), $(Pretty.quoteof(r.iter))))
+
+function DynamicFor(κ::K,iter::S) where {K,S}
+    T = typeof(κ(first(iter)))
+    DynamicFor{T,K,S}(κ,iter)
+end
+
+function Base.rand(rng::AbstractRNG, T::Type, df::DynamicFor)
+    r = ResettableRNG(rng)
+    return RealizedSamples(r, df)
+end
+
+function logdensity(df::DynamicFor, y)
+    ℓ = 0.0
+    for (xj, yj) in zip(df.iter, y)
+        ℓ += logdensity(df.κ(xj), yj)
+    end
+    return ℓ
 end
 
 Base.eltype(::Type{D}) where {T,D<:DynamicFor{T}} = eltype(T)
 
 Base.IteratorEltype(d::DynamicFor) = Base.HasEltype()
 
-Base.IteratorSize(d::DynamicFor) = Base.IteratorSize(d.sampler)
+Base.IteratorSize(d::DynamicFor) = Base.IteratorSize(d.iter)
 
 function Base.iterate(d::DynamicFor)
-    (x,s) = iterate(d.sampler)
+    (x,s) = iterate(d.iter)
     (d.κ(x), s)
 end
 
 function Base.iterate(d::DynamicFor, s)
-    (x,s) = iterate(d.sampler, s)
+    (x,s) = iterate(d.iter, s)
     (d.κ(x), s)
 end
 
-Base.length(d::DynamicFor) = length(d.sampler)
+Base.length(d::DynamicFor) = length(d.iter)
 
 
-For(f, r::Realized) = DynamicFor(f,r)
+For(f, r::Realized) = DynamicFor(f, copy(r))
 
 function Base.rand(rng::AbstractRNG, dfor::DynamicFor)
-    seed = rand(rng, UInt)
-    return Realized(seed, copy(rng), dfor)
+    r = ResettableRNG(rng)
+    return RealizedSamples(r, dfor)
 end
 
 function dyniterate(df::DynamicFor, st, args...)
@@ -158,11 +124,11 @@ For(f, it::DynamicIterator) = DynamicFor(f, it)
 
 For(f, it::DynamicFor) = DynamicFor(f, it)
 
-function dyniterate(fr::DynamicFor, state)
-      ϕ = dyniterate(fr.iter, state)
+function dyniterate(df::DynamicFor, state)
+      ϕ = dyniterate(df.iter, state)
       ϕ === nothing && return nothing
       u, state = ϕ
-      fr.f(u), state
+      df.f(u), state
 end
 
 function Base.collect(r::Realized)

src/distributions.jl

@@ -33,6 +33,10 @@ function basemeasure(μ::Dists.Distribution{Dists.Univariate,Dists.Continuous})
     return Lebesgue(ℝ)
 end
 
+function basemeasure(μ::Dists.Distribution{Dists.Univariate,Dists.Discrete})
+    return CountingMeasure(ℤ)
+end
+
 ∫(::typeof(identity), ::Dists.Distribution) = 1.0
 
 logdensity(μ::Dists.Distribution, x) = Dists.logpdf(μ,x)

src/parameterized/beta.jl

@@ -2,7 +2,7 @@
 
 export Beta
 
-@parameterized Beta(α,β) ≃ Lebesgue(𝕀)
+@parameterized Beta(α,β)
 
 @kwstruct Beta(α, β)
 
@@ -15,8 +15,20 @@ export Beta
 
 TV.as(::Beta) = as𝕀
 
-function logdensity(d::Beta{(:α, :β)}, x)
-    return xlogy(d.α - 1, x) + xlog1py(d.β - 1, -x) - logbeta(d.α, d.β)
+function logdensity(d::Beta{(:α, :β), Tuple{A,B}}, x::X) where {A,B,X}
+    if static_hasmethod(xlogy, Tuple{A,X}) && static_hasmethod(xlog1py, Tuple{B,X})
+        return xlogy(d.α - 1, x) + xlog1py(d.β - 1, -x) 
+    else
+        return (d.α - 1) * log(x) + (d.β - 1) * log1p(-x)
+    end
+end
+
+function basemeasure(d::Beta{(:α,:β)})
+    inbounds(x) = 0 < x < 1
+    constℓ = 0.0
+    varℓ() = - logbeta(d.α, d.β)
+    base = Lebesgue(ℝ)
+    FactoredBase(inbounds, constℓ, varℓ, base)
 end
 
 Base.rand(rng::AbstractRNG, T::Type, μ::Beta) = rand(rng, Dists.Beta(μ.α, μ.β))

src/parameterized/studentt.jl

@@ -44,8 +44,8 @@ function logdensity(d::StudentT{(:ν,)}, x)
     return  (ν + 1) / (-2) * log1p(x^2 / ν)
 end
 
-function basemeasure(d::StudentT{(:ν,)})
-    inbounds(x) = true
+@inline function basemeasure(d::StudentT{(:ν,)})
+    inbounds = Returns(true)
     constℓ = 0.0
     varℓ() = loggamma((d.ν+1)/2) - loggamma(d.ν/2) - log(π * d.ν) / 2
     base = Lebesgue(ℝ)