Dev

Project.toml

@@ -38,7 +38,7 @@ TupleVectors = "615932cf-77b6-4358-adcd-5b7eba981d7e"
 
 [compat]
 Accessors = "0.1"
-ArrayInterface = "4, 5, 6"
+ArrayInterface = "5, 6"
 DataStructures = "0.18"
 DensityInterface = "0.4"
 DiffResults = "1"
@@ -47,7 +47,7 @@ JuliaVariables = "0.2"
 MLStyle = "0.3,0.4"
 MacroTools = "0.5"
 MappedArrays = "0.3, 0.4"
-MeasureBase = "0.9"
+MeasureBase = "0.12"
 MeasureTheory = "0.16"
 NamedTupleTools = "0.12, 0.13, 0.14"
 NestedTuples = "0.3"
@@ -63,7 +63,7 @@ StatsFuns = "0.9, 1"
 TransformVariables = "0.5, 0.6"
 Tricks = "0.1"
 TupleVectors = "0.1"
-julia = "1.5"
+julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"

benchmarks/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
+OnlineStats = "a15396b6-48d5-5d58-9928-6d29437db91e"
 Pathfinder = "b1d3bc72-d0e7-4279-b92f-7fa5d6d2d454"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

benchmarks/bouncy.jl

@@ -79,12 +79,9 @@ M = pf_result.fit_distribution.Σ;
 x0 = pf_result.fit_distribution.μ;
 v0 = PDMats.unwhiten(M, randn(length(x0)));
 
-
-
-
-
 MAP = pf_result.optim_solution; # MAP, could be useful for control variates
 
+
 # define BouncyParticle sampler (has two relevant parameters) 
 Z = BouncyParticle(missing, # graphical structure 
     MAP, # MAP estimate, unused
@@ -103,6 +100,7 @@ sampler = ZZB.NotFactSampler(Z, (dneglogp, ∇neglogp!), ZZB.LocalBound(c), t0 =
 using TupleVectors: chainvec
 using Tilde.MeasureTheory: transform
 
+# @time first(Iterators.drop(tvs,1000))
 
 function collect_sampler(t, sampler, n; progress=true, progress_stops=20)
     if progress
@@ -173,4 +171,5 @@ ylabel!(plt, "DynamicHMC");
 plt_bounds = collect(extrema(ess_hmc));
 lineplot!(plt, plt_bounds, plt_bounds);
 plt
-@info "For each coordinate, a point (x,y) shows the effective sample size per second for BPS (x) and HMC (y) . In blue is the diagonal x=y"
+@info "For each coordinate, a point (x,y) shows the effective sample size per second for BPS (x) and HMC (y) . In blue is the diagonal x=y"
+

src/Tilde.jl

@@ -83,7 +83,6 @@ end
 include("optics.jl")
 include("maybe.jl")
 include("core/models/abstractmodel.jl")
-include("core/models/astmodel/astmodel.jl")
 include("core/models/model.jl")
 include("core/dependencies.jl")
 include("core/utils.jl")

src/core/models/abstractmodel.jl

@@ -3,6 +3,12 @@
 # struct MixedSupport <: ValueSupport end
 # struct MixedVariate <: VariateForm end
 
+abstract type ModelSupport end
+
+struct LatentSupport <: ModelSupport end
+struct PushforwardSupport <: ModelSupport end
+struct JointSupport <: ModelSupport end
+
 """
     AbstractModel{A,B}
 
@@ -15,7 +21,7 @@ N gives the Names of arguments (each a Symbol)
 B gives the Body, as an Expr
 M gives the Module where the model is defined
 """
-abstract type AbstractModel{A,B,M} <: AbstractTransitionKernel end
+abstract type AbstractModel{A,B,M,F} <: AbstractTransitionKernel end
 
 abstract type AbstractConditionalModel{M, Args, Obs} <: AbstractMeasure end
 

src/core/models/model.jl

@@ -1,3 +1,49 @@
+struct Model{A,B,M<:GG.TypeLevel,F} <: AbstractModel{A,B,M,F}
+    args :: Vector{Symbol}
+    body :: Expr
+    f    :: F
+end
+
+function Model(theModule::Module, args::Vector{Symbol}, body::Expr)
+    A = NamedTuple{Tuple(args)}
+
+    B = to_type(body)
+    M = to_type(theModule)
+    return Model{A,B,M,typeof(last)}(args, body, last)
+end
+
+model(m::Model) = m
+
+function Base.convert(::Type{Expr}, m::Model)
+    numArgs = length(m.args)
+    args = if numArgs == 1
+       m.args[1]
+    elseif numArgs > 1
+        Expr(:tuple, [x for x in m.args]...)
+    end
+
+    body = m.body
+
+    q = if numArgs == 0
+        @q begin
+            @model $body
+        end
+    else
+        @q begin
+            @model $(args) $body
+        end
+    end
+
+    striplines(q).args[1] 
+end
+
+Base.show(io::IO, m :: Model) = println(io, convert(Expr, m))
+
+function type2model(::Type{Model{A,B,M,F}}) where {A,B,M,S}
+    args = [fieldnames(A)...]
+    body = from_type(B)
+    Model{A,B,M,F}(from_type(M), convert(Vector{Symbol},args), body)
+end
 
 toargs(vs :: Vector{Symbol}) = Tuple(vs)
 toargs(vs :: NTuple{N,Symbol} where {N}) = vs

src/core/utils.jl

@@ -110,15 +110,6 @@ end
 import MacroTools: striplines, @q
 
 
-
-
-# function arguments(model::DAGModel)
-#     model.args
-# end
-
-
-
-
 allequal(xs) = all(xs[1] .== xs)
 
 

src/primitives/insupport.jl

@@ -2,5 +2,5 @@ import MeasureBase: insupport
 export insupport
 
 @inline function insupport(m::AbstractConditionalModel, x::NamedTuple)
-    mapreduce(insupport, (a,b) -> a && b, measures!(m,x), x)
+    mapreduce(insupport, (a,b) -> a && b, measures(m), x)
 end

src/primitives/interpret.jl

@@ -45,13 +45,21 @@ function make_body(M, f, ast::Expr, retfun, argsT, obsT, parsT)
                 # X = to_type(unsolved_lhs)
                 # M = to_type(unsolve(rhs))
 
-                inargs = inkeys(sx, argsT)
+                # inargs = inkeys(sx, argsT)
                 inobs = inkeys(sx, obsT)
-                inpars = inkeys(sx, parsT)
+                # inpars = inkeys(sx, parsT)
                 rhs = unsolve(rhs)
 
-                xval = inobs ? :($Observed($x)) : (x ∈ knownvars ? :($Unobserved($x)) : :($Unobserved(missing)))
+                xval = if inkeys(sx, obsT)
+                    :($Observed($x))
+                elseif x ∈ knownvars 
+                    :($Unobserved($x)) 
+                else 
+                    :($Unobserved(missing))
+                end
+
                 st = :(($x, _ctx, _retn) = $tilde($f, $l, $sx, $xval, $rhs, _cfg, _ctx))
+
                 qst = QuoteNode(st)
                 q = quote
                     # println($qst)

src/primitives/rand.jl

@@ -4,48 +4,73 @@ using TupleVectors: chainvec
 export rand
 EmptyNTtype = NamedTuple{(),Tuple{}} where T<:Tuple
 
-@inline function Base.rand(rng::AbstractRNG, d::AbstractConditionalModel, N::Int)
-    r = chainvec(rand(rng, d), N)
+@inline function Base.rand(m::ModelClosure, args...; kwargs...) 
+    rand(GLOBAL_RNG, Float64, m, args...; kwargs...)
+end
+
+@inline function Base.rand(rng::AbstractRNG, m::ModelClosure, args...; kwargs...) 
+    rand(rng, Float64, m, args...; kwargs...)
+end
+
+@inline function Base.rand(::Type{T_rng}, m::ModelClosure, args...; kwargs...) where {T_rng}
+    rand(GLOBAL_RNG, T_rng, m, args...; kwargs...)
+end
+
+@inline function Base.rand(m::ModelClosure, d::Integer, dims::Integer...; kwargs...) 
+    rand(GLOBAL_RNG, Float64, m, d, dims...; kwargs...)
+end
+
+@inline function Base.rand(rng::AbstractRNG, m::ModelClosure, d::Integer, dims::Integer...; kwargs...) 
+    rand(rng, Float64, m, d, dims...; kwargs...)
+end
+
+@inline function Base.rand(::Type{T_rng}, m::ModelClosure, d::Integer, dims::Integer...; kwargs...) where {T_rng}
+    rand(GLOBAL_RNG, T_rng, m, d, dims...; kwargs...)
+end
+
+@inline function Base.rand(rng::AbstractRNG, ::Type{T_rng}, d::ModelClosure, N::Integer, v::Vararg{Integer}) where {T_rng}
+    @assert isempty(v)
+    r = chainvec(rand(rng, T_rng, d), N)
     for j in 2:N
-        @inbounds r[j] = rand(rng, d)
+        @inbounds r[j] = rand(rng, T_rng, d)
     end
     return r
 end
 
-@inline Base.rand(d::AbstractConditionalModel, N::Int) = rand(GLOBAL_RNG, d, N)
+@inline Base.rand(d::ModelClosure, N::Int) = rand(GLOBAL_RNG, d, N)
 
-@inline function Base.rand(m::AbstractConditionalModel; kwargs...) 
+@inline function Base.rand(m::ModelClosure; kwargs...) 
     rand(GLOBAL_RNG, m; kwargs...)
 end
 
-@inline function Base.rand(rng::AbstractRNG, m::AbstractConditionalModel; ctx=NamedTuple(), retfun = (r, ctx) -> r)
+@inline function Base.rand(rng::AbstractRNG, m::ModelClosure; ctx=NamedTuple(), retfun = (r, ctx) -> r)
     cfg = (rng=rng,)
+
+@inline function Base.rand(rng::AbstractRNG, ::Type{T_rng}, m::ModelClosure; ctx=NamedTuple(), retfun = (r, ctx) -> r) where {T_rng}
+    cfg = (rng=rng, T_rng=T_rng)
     gg_call(rand, m, NamedTuple(), cfg, ctx, retfun)
 end
 
 ###############################################################################
 # ctx::NamedTuple
-@inline function tilde(::typeof(Base.rand), lens, xname, x::Unobserved, d, cfg, ctx::NamedTuple)
+@inline function tilde(::typeof(Base.rand), lens, xname, x, d, cfg, ctx::NamedTuple)
     xnew = set(x.value, Lens!!(lens), rand(cfg.rng, d))
     ctx′ = merge(ctx, NamedTuple{(dynamic(xname),)}((xnew,)))
     (xnew, ctx′, nothing)
 end
 
-@inline function tilde(::typeof(Base.rand), lens, xname, x::Observed, d, cfg, ctx::NamedTuple)
-    (x.value, ctx, nothing)
-end
 
 
 ###############################################################################
 # ctx::Dict
 
 @inline function tilde(::typeof(Base.rand), lens::typeof(identity), xname, x, d, cfg, ctx::Dict)
-    x = rand(cfg.rng, d)
+    x = rand(cfg.rng, cfg.T_rng, d)
     ctx[dynamic(xname)] = x 
     (x, ctx, nothing)
 end
 
-@inline function tilde(::typeof(Base.rand), lens, xname, x, m::AbstractConditionalModel, cfg, ctx::Dict)
+@inline function tilde(::typeof(Base.rand), lens, xname, x, m::ModelClosure, cfg, ctx::Dict)
     args = get(cfg.args, dynamic(xname), Dict())
     cfg = merge(cfg, (args = args,))
     tilde(rand, lens, xname, x, m(cfg.args), cfg, ctx)

test/transforms.jl

@@ -1,12 +1,3 @@
-# Check for Model equality up to reorderings of a few fields
-function ≊(m1::DAGModel,m2::DAGModel)
-    function eq_tuples(nt1::NamedTuple,nt2::NamedTuple)
-        return length(nt1)==length(nt2) && all(nt1[k]==nt2[k] for k in keys(nt1))
-    end
-    return Set(arguments(m1))==Set(arguments(m2)) && m1.retn==m2.retn && eq_tuples(m1.dists,m2.dists) && eq_tuples(m1.vals,m2.vals)
-end
-
-
 m = @model (n,α,β) begin
     p ~ Beta(α, β)
     x ~ Binomial(n, p)