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

@@ -17,6 +17,29 @@ using Pathfinder.PDMats
 
 Random.seed!(1)
 
+function make_grads(post)    
+    as_post = as(post)
+    obj(θ) = -Tilde.unsafe_logdensityof(post, transform(as_post, θ))
+    ℓ(θ) = -obj(θ)
+    @inline function dneglogp(t, x, v) # two directional derivatives
+        f(t) = obj(x + t*v)
+        u = ForwardDiff.derivative(f, Dual{:hSrkahPmmC}(0.0, 1.0))
+        u.value, u.partials[]
+    end
+
+    gconfig = ForwardDiff.GradientConfig(obj, rand(25), ForwardDiff.Chunk{25}())
+    function ∇neglogp!(y, t, x)
+        ForwardDiff.gradient!(y, obj, x, gconfig)
+        y
+    end
+    ℓ, dneglogp, ∇neglogp!
+end
+
+
+# ↑ general purpose
+############################################################
+# ↓ problem-specific
+
 # read data
 function readlrdata()
     fname = joinpath("lr.data")
@@ -26,45 +49,36 @@ function readlrdata()
     y = z[:,end] .- 1
     return A, y
 end
-A, y = readlrdata();
-At = collect(A');
+
 
 model_lr = @model (At, y, σ) begin
     d,n = size(At)
     θ ~ Normal(σ=σ)^d
     for j in 1:n
-        logitp = dot(view(At,:,j), θ)
+        logitp = view(At,:,j)' * θ
         y[j] ~ Bernoulli(logitp = logitp)
     end
 end
+
+A, y = readlrdata();
+At = collect(A');
 σ = 100.0
 
-function make_grads(model_lr, At, y, σ)    
-    post = model_lr(At, y, σ) | (;y)
-    as_post = as(post)
-    obj(θ) = -Tilde.unsafe_logdensityof(post, transform(as_post, θ))
-    ℓ(θ) = -obj(θ)
-    @inline function dneglogp(t, x, v) # two directional derivatives
-        f(t) = obj(x + t*v)
-        u = ForwardDiff.derivative(f, Dual{:hSrkahPmmC}(0.0, 1.0))
-        u.value, u.partials[]
-    end
-
-    gconfig = ForwardDiff.GradientConfig(obj, rand(25), ForwardDiff.Chunk{25}())
-    function ∇neglogp!(y, t, x)
-        ForwardDiff.gradient!(y, obj, x, gconfig)
-        return
-    end
-    post, ℓ, dneglogp, ∇neglogp!
-end
 
-post, ℓ, dneglogp, ∇neglogp! = make_grads(model_lr, At, y, σ)  
-# Try things out
-dneglogp(2.4, randn(25), randn(25));
-∇neglogp!(randn(25), 2.1, randn(25));
+post = model_lr(At, y, σ) | (;y)
+
+d = as(post).dimension
+
+let
+    ℓ, dneglogp, ∇neglogp! = make_grads(post)  
+    # Try things out
+    @show dneglogp(2.4, randn(d), randn(d));
+    y = Vector{Float64}(undef, d)
+    @show ∇neglogp!(y, 2.1, randn(d));
+    nothing
+end
 
 
-d = 25 # number of parameters 
 t0 = 0.0;
 x0 = zeros(d); # starting point sampler
 # estimated posterior mean (n=100000, 797s)
@@ -79,12 +93,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 +114,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 +185,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")
@@ -103,6 +102,7 @@ include("primitives/testparams.jl")
 include("primitives/weightedsampling.jl")
 include("primitives/measures.jl")
 include("primitives/basemeasure.jl")
+include("primitives/predict.jl")
 
 include("transforms/utils.jl")
 

src/core/models/model.jl

@@ -1,3 +1,47 @@
+struct Model{A,B,M<:GG.TypeLevel} <: AbstractModel{A,B,M}
+    args :: Vector{Symbol}
+    body :: Expr
+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}(args, body)
+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}}) where {A,B,M}
+    args = Symbol[fieldnames(A)...]
+    body = from_type(B)
+    Model{A,B,M}(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)