diff --git a/Project.toml b/Project.toml
index c97e3df7..b80b3b80 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,7 +5,6 @@ version = "0.8.0"
 
 [deps]
 AdvancedHMC = "0bf59076-c3b1-5ca4-86bd-e02cd72cde3d"
-Bijectors = "76274a88-744f-5084-9051-94815aaf08c4"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DynamicHMC = "bbc10e6e-7c05-544b-b16e-64fede858acb"
@@ -13,10 +12,8 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeneralizedGenerated = "6b9d7cbe-bcb9-11e9-073f-15a7a543e2eb"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-IRTools = "7869d1d1-7146-5819-86e3-90919afe41df"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogDensityProblems = "6fdf6af0-433a-55f7-b3ed-c6c6e0b8df7c"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
@@ -34,14 +31,11 @@ SimplePartitions = "ec83eff0-a5b5-5643-ae32-5cbf6eedec9d"
 SimplePosets = "b2aef97b-4721-5af9-b440-0bad754dc5ba"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-Stheno = "8188c328-b5d6-583d-959b-9690869a5511"
 SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 TransformVariables = "84d833dd-6860-57f9-a1a7-6da5db126cff"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
 AdvancedHMC = "0.2"
-Bijectors = "0.4"
 DiffResults = "0.0.4"
 Distributions = "0.21"
 DynamicHMC = "2.1"
@@ -49,7 +43,6 @@ FillArrays = "0.8"
 ForwardDiff = "0.10"
 GeneralizedGenerated = "0.2"
 Graphs = "0.10"
-IRTools = "0.3"
 IterTools = "1.3"
 LazyArrays = "0.14"
 LogDensityProblems = "0.9"
@@ -66,10 +59,8 @@ SimpleGraphs = "0.3"
 SimplePartitions = "0.2"
 SimplePosets = "0.0"
 StatsFuns = "0.9"
-Stheno = "0.3, 0.4, 0.5"
 SymPy = "1.0"
 TransformVariables = "0.3"
-Zygote = "0.4"
 julia = "^1.3"
 
 [extras]
diff --git a/docs/src/probprog2020/probprog2020.jmd b/docs/src/probprog2020/probprog2020.jmd
new file mode 100644
index 00000000..539d572f
--- /dev/null
+++ b/docs/src/probprog2020/probprog2020.jmd
@@ -0,0 +1,90 @@
+
+```julia
+using Revise, Soss
+import SymPy
+```
+
+
+A `Model` is a 
+
+```julia
+m = @model x begin
+    α ~ Normal()
+    β ~ Normal()
+    σ ~ HalfNormal()
+    yhat = α .+ β .* x
+    n = length(x)
+    y ~ For(n) do j
+        Normal(yhat[j], σ)
+    end
+end;
+```
+
+
+## Symbolic Simplification
+
+A common bottleneck in probabilistic programming comes from evaluation of the log-posterior density, often required in an inner loop. Because of this, optimizing this computation can be a significant benefit to overall computational cost.
+
+## Code Generation
+
+Given values for model arguments and all stochastic variables, ...
+
+- 
+
+```julia
+Soss.sourceSymlogpdf()(m)
+```
+
+
+## Model Transformations
+
+### Markov Blankets
+
+In a DAG-based model, a _Markov blanket_ at a node $v$ is given by the following "moralizing" procedure:
+
+*TODO reference*
+
+1. Connect $v$ to parents of its children
+2. Replace directed edges with undirected edges
+3. Return the nodes connected to $v$ (including $v$)
+
+Markov blankets are commonly used in Gibbs sampling, where a given node can be updated by sampling from the posterior distribution for that node. This is complicated by the presence of deterministic nodes, which can introduce constraints that disallow any update. we address this problem be extending the blanket to ancestors, stopping only when a stochastic node is reached. For example, we have
+
+```julia
+markovBlanket(m,:α)
+```
+
+Note that the arguments of this model consist of the parents of a `:α`, together with its "stochastic childrens'" parents. Equivalently (and more simply), we must include any referenced but undefined variable as an argument.
+
+For deterministic 
+
+```julia
+markovBlanket(m, :yhat)
+```
+
+# Internals
+
+When designing a PPL, it's important to consider the form of the intermediate representation. 
+
+In Soss, a `Model` is a `struct` of the form
+
+```julia
+struct Model{A,B,M} 
+    args  :: Vector{Symbol}
+    vals  :: NamedTuple
+    dists :: NamedTuple
+    retn  :: Union{Nothing, Symbol, Expr}
+end
+```
+
+
+
+
+
+```julia
+symlogpdf(m)
+```
+
+```julia
+symlogpdf(m) |> expandSums |> foldConstants
+```
diff --git a/examples/2019-11-07-demo.jmd b/examples/2019-11-07-demo.jmd
index 601508f1..f05b4f0f 100644
--- a/examples/2019-11-07-demo.jmd
+++ b/examples/2019-11-07-demo.jmd
@@ -166,7 +166,7 @@ scatter!(truth.x,truth.y, legend=false, c=1)
 ```
 
 ```julia
-pred2 = predictive(m2, setdiff(stochastic(m2), [:y])...)
+pred2 = predictive(m2, setdiff(sampled(m2), [:y])...)
 ```
 
 ```julia
diff --git a/src/Soss.jl b/src/Soss.jl
index 0541fe8e..db1d1021 100644
--- a/src/Soss.jl
+++ b/src/Soss.jl
@@ -11,7 +11,6 @@ import MacroTools: prewalk, postwalk, @q, striplines, replace, flatten, @capture
 import MLStyle
 @reexport using MonteCarloMeasurements
 
-using LinearAlgebra
 using LazyArrays
 using FillArrays
 
diff --git a/src/core/model.jl b/src/core/model.jl
index 605b6c45..9a79725c 100644
--- a/src/core/model.jl
+++ b/src/core/model.jl
@@ -13,11 +13,11 @@ struct Model{A,B,M}
     retn  :: Union{Nothing, Symbol, Expr}
 end
 
-argstype(::Model{A,B}) where {A,B} = A
-bodytype(::Model{A,B}) where {A,B} = B
+argstype(::Model{A,B,M}) where {A,B,M} = A
+bodytype(::Model{A,B,M}) where {A,B,M} = B
 
-argstype(::Type{Model{A,B}}) where {A,B} = A
-bodytype(::Type{Model{A,B}}) where {A,B} = B
+argstype(::Type{Model{A,B,M}}) where {A,B,M} = A
+bodytype(::Type{Model{A,B,M}}) where {A,B,M} = B
 
 getmodule(::Type{Model{A,B,M}}) where {A,B,M} = from_type(M)
 getmodule(::Model{A,B,M}) where {A,B,M} = from_type(M)
diff --git a/src/core/utils.jl b/src/core/utils.jl
index 6b0481b1..fde5721f 100644
--- a/src/core/utils.jl
+++ b/src/core/utils.jl
@@ -17,8 +17,23 @@ argtuple(m) = arguments(m) |> astuple
 astuple(x) = Expr(:tuple,x...)
 astuple(x::Symbol) = Expr(:tuple,x)
 
+
+
+export arguments
+arguments(m::Model) = m.args
+arguments(d::JointDistribution) = d.args
+
+export sampled
+sampled(m::Model) = keys(m.dists) |> collect
+
+export assigned
+assigned(m::Model) = keys(m.vals) |> collect
+
+export parameters
+parameters(m::Model) = union(assigned(m), sampled(m))
+
 export variables
-variables(m::Model) = m.args ∪ keys(m.vals) ∪ keys(m.dists)
+variables(m::Model) = union(arguments(m), parameters(m))
 
 function variables(expr :: Expr) 
     leaf(x::Symbol) = begin
@@ -34,29 +49,13 @@ end
 variables(s::Symbol) = [s]
 variables(x) = []
 
-
-export arguments
-arguments(m) = m.args
-
-export stochastic
-stochastic(m::Model) = keys(m.dists)
-
-export bound
-bound(m::Model) = keys(m.vals)
-
-export bodyVariables
-bodyVariables(m::Model) = setdiff(variables(m), arguments(m))
-
-# TODO: Fix these broken methods
-# export observed
-# observed(m::Model) = keys(m.data)
-
-# export parameters
-# parameters(m::Model) = setdiff(
-#     stochastic(m), 
-#     observed(m)
-# )
-
+for f in [:arguments, :assigned, :sampled, :parameters, :variables]
+    @eval function $f(m::Model, nt::NamedTuple) 
+        vs = $f(m)
+        isempty(vs) && return NamedTuple()
+        return select(nt, $f(m))
+    end
+end
 
 
 export foldall
diff --git a/src/examples/repeated-measurements.jl b/src/examples/repeated-measurements.jl
index 7ea1ee2d..75dea82d 100644
--- a/src/examples/repeated-measurements.jl
+++ b/src/examples/repeated-measurements.jl
@@ -2,8 +2,8 @@ using Revise
 using Soss
 
 m = @model begin
-    n = 10 # number of samples
-    m = 10 # reviewers for each example
+    n = 100 # number of samples
+    m = 100 # reviewers for each example
 
     p_bad ~ Beta(1,3) |> iid(n)
     
@@ -20,16 +20,21 @@ using Random
 truth = rand(m());
 symlogpdf(m())
 
-@time logpdf(m(),truth, codegen)
-@time logpdf(m(),truth, codegen)
+using BenchmarkTools
+@btime logpdf(m(),truth)
+@btime logpdf(m(),truth, codegen)
 
+f1 = Soss._codegen(m, true);
+f2 = Soss._codegen(m,false);
 
+@btime f1((),truth)
+@btime f2((),truth)
 
 codegen(m(),truth)
 
 
 
-logpdf(m(), merge(truth, (p_bad=shuffle(truth.p_bad),)))
+logpdf(m(), merge(truth, (p_bad=shuffle(truth.p_bad),)), codegen)
 
 
 @time result = dynamicHMC(m(), (y=truth.y,), codegen) ;
diff --git a/src/examples/strata.jl b/src/examples/strata.jl
index f5e2676c..247639ea 100644
--- a/src/examples/strata.jl
+++ b/src/examples/strata.jl
@@ -109,7 +109,7 @@ scatter!(truth.x,truth.y, legend=false, c=1)
 plot!(xx, truth.α .+ truth.β .* xx, legend=false, lw=3, c=2)
 savefig("fit2.png")
 
-pred2 = predictive(m2, setdiff(stochastic(m2), [:y])...)
+pred2 = predictive(m2, setdiff(sampled(m2), [:y])...)
 
 post2pred= [rand(pred2(θ)((x=x,))) for θ ∈ post2]  |> particles
 
diff --git a/src/fromcube.jl b/src/fromcube.jl
index f594bcb2..35952cc3 100644
--- a/src/fromcube.jl
+++ b/src/fromcube.jl
@@ -44,7 +44,7 @@ function sourceFromcube(m::Model)
     argsExpr = Expr(:tuple,freeVariables(m)...)
 
     stochExpr = begin
-        vals = map(stochastic(m)) do x Expr(:(=), x,x) end
+        vals = map(sampled(m)) do x Expr(:(=), x,x) end
         Expr(:tuple, vals...)
     end
     
diff --git a/src/importance.jl b/src/importance.jl
index f9935392..87d1556d 100644
--- a/src/importance.jl
+++ b/src/importance.jl
@@ -47,7 +47,7 @@ function sourceImportanceSample()
         kwargsExpr = Expr(:tuple,kwargs...)
 
         stochExpr = begin
-            vals = map(stochastic(m)) do x Expr(:(=), x,x) end
+            vals = map(sampled(m)) do x Expr(:(=), x,x) end
             Expr(:tuple, vals...)
         end
 
diff --git a/src/inference/dynamicHMC.jl b/src/inference/dynamicHMC.jl
index 9d6116e7..b1ba9ac4 100644
--- a/src/inference/dynamicHMC.jl
+++ b/src/inference/dynamicHMC.jl
@@ -1,44 +1,88 @@
-using TransformVariables, LogDensityProblems, DynamicHMC, 
-    Distributions, Statistics, StatsFuns, ForwardDiff
+using TransformVariables,
+      LogDensityProblems,
+      DynamicHMC,
+      Distributions,
+      Statistics,
+      StatsFuns,
+      ForwardDiff
+import LogDensityProblems: ADgradient
 
 using Random
 
 export dynamicHMC
 
-# import Flux
-function dynamicHMC(m :: JointDistribution, _data, method=logpdf, N=1000::Int) 
+function dynamicHMC(
+    rng::AbstractRNG,
+    m::JointDistribution,
+    _data,
+    N::Int = 1000;
+    method = logpdf,
+    ad_backend = Val(:ForwardDiff),
+    reporter = DynamicHMC.NoProgressReport(),
+    kwargs...,
+)
     ℓ(pars) = logpdf(m, merge(pars, _data), method)
+    t = xform(m, _data)
+    P = LogDensityProblems.TransformedLogDensity(t, ℓ)
+    ∇P = LogDensityProblems.ADgradient(ad_backend, P)
 
-    t = xform(m,_data)
-    P = TransformedLogDensity(t, ℓ)
-    ∇P = ADgradient(Val(:ForwardDiff), P)
-    results = mcmc_with_warmup(MersenneTwister(), ∇P, N; reporter=DynamicHMC.NoProgressReport());
-    samples = TransformVariables.transform.(parent(∇P).transformation, results.chain)
+    results = DynamicHMC.mcmc_with_warmup(
+        rng,
+        ∇P,
+        N;
+        reporter = reporter,
+        kwargs...,
+    )
+    samples = TransformVariables.transform.(t, results.chain)
+    return samples
 end
 
+function dynamicHMC(
+    rng::AbstractRNG,
+    m::JointDistribution,
+    _data,
+    ::Val{Inf};
+    method = logpdf,
+    ad_backend = Val(:ForwardDiff),
+    reporter = DynamicHMC.NoProgressReport(),
+    kwargs...,
+)
+    ℓ(pars) = logpdf(m, merge(pars, _data), method)
+    t = xform(m, _data)
+    P = LogDensityProblems.TransformedLogDensity(t, ℓ)
+    ∇P = LogDensityProblems.ADgradient(ad_backend, P)
 
-function dynamicHMC(m :: JointDistribution, _data, ::Val{Inf}) 
-    ℓ(pars) = logpdf(m, merge(pars, _data))
-
-    t = xform(m,_data)
-    P = TransformedLogDensity(t, ℓ)
-    ∇P = ADgradient(Val(:ForwardDiff), P)
-    
-    # initialization
-    rng = MersenneTwister()
-    results = DynamicHMC.mcmc_keep_warmup(rng, ∇P, 0; reporter = NoProgressReport())
-    steps = DynamicHMC.mcmc_steps(results.sampling_logdensity, results.final_warmup_state)
+    results = DynamicHMC.mcmc_keep_warmup(
+        rng,
+        ∇P,
+        0;
+        reporter = reporter,
+        kwargs...,
+    )
+    steps = DynamicHMC.mcmc_steps(
+        results.sampling_logdensity,
+        results.final_warmup_state,
+    )
+    return results, steps
+end
 
-    (results, steps)
+function dynamicHMC(m::JointDistribution, args...; kwargs...)
+    return dynamicHMC(Random.GLOBAL_RNG, m, args...; kwargs...)
 end
 
 
 using ResumableFunctions
 
 export stream
-@resumable function stream(f::typeof(dynamicHMC), m :: JointDistribution, _data::NamedTuple) 
+
+@resumable function stream(
+    rng::AbstractRNG,
+    f::typeof(dynamicHMC),
+    m::JointDistribution,
+    _data::NamedTuple,
+)
     t = xform(m, _data)
-    (results, steps) = dynamicHMC(m, _data, Val(Inf))
+    (results, steps) = dynamicHMC(rng, m, _data, Val(Inf))
     Q = results.final_warmup_state.Q
     while true
         Q, tree_stats = DynamicHMC.mcmc_next_step(steps, Q)
@@ -46,3 +90,11 @@ export stream
     end
 end
 
+function stream(
+    f::typeof(dynamicHMC),
+    m::JointDistribution,
+    _data::NamedTuple;
+    kwargs...,
+)
+    return stream(Random.GLOBAL_RNG, f, m, _data; kwargs...)
+end
diff --git a/src/particles.jl b/src/particles.jl
index da067d57..624704de 100644
--- a/src/particles.jl
+++ b/src/particles.jl
@@ -1,4 +1,4 @@
-
+const DEFAULT_SAMPLE_SIZE = 1000
 
 export sourceParticles
 
@@ -9,9 +9,7 @@ end
 export particles
 
 particles(v::Vector{<:Real}) = Particles(v)
-particles(d) = begin
-    Particles(1000, d)
-end
+particles(d, N::Int=DEFAULT_SAMPLE_SIZE) = Particles(N, d)
 
 # particles(n :: NUTS_result) = particles(n.samples)
 using IterTools
@@ -33,18 +31,19 @@ export parts
 
 # Just a little helper function for particles
 # https://github.com/baggepinnen/MonteCarloMeasurements.jl/issues/22
-parts(d; N=1000) = particles(m)
-parts(x::Normal{P} where {P <: AbstractParticles}; N=1000) = Particles(length(x.μ), Normal()) * x.σ + x.μ
-parts(x::Sampleable{F,S}; N=1000) where {F,S} = Particles(N,x)
-parts(x::Integer; N=1000) = parts(float(x))
-parts(x::Real; N=1000) = parts(repeat([x],N))
-parts(x::AbstractArray; N=1000) = Particles(x)
-parts(p::Particles; N=1000) = p 
-parts(d::For; N=1000) = parts.(d.f.(d.θ...))
-
-
+parts(d, N::Int=DEFAULT_SAMPLE_SIZE) = particles(d, N)
+parts(x::Normal{P} where {P <: AbstractParticles}, N::Int=DEFAULT_SAMPLE_SIZE) = Particles(length(x.μ), Normal()) * x.σ + x.μ
+parts(x::Sampleable{F,S}, N::Int=DEFAULT_SAMPLE_SIZE) where {F,S} = Particles(N,x)
+parts(x::Integer, N::Int=DEFAULT_SAMPLE_SIZE) = parts(float(x))
+parts(x::Real, N::Int=DEFAULT_SAMPLE_SIZE) = parts(repeat([x],N))
+parts(x::AbstractArray, N::Int=DEFAULT_SAMPLE_SIZE) = Particles(x)
+parts(p::Particles, N::Int=DEFAULT_SAMPLE_SIZE) = p
+parts(d::For, N::Int=DEFAULT_SAMPLE_SIZE) = parts.(d.f.(d.θ...), N)
+function parts(d::For{F,Tuple{I}}, N::Int=DEFAULT_SAMPLE_SIZE) where {F <: Function, I <: Integer}
+    parts.(d.f.(Base.OneTo.(d.θ)...), N)
+end
 
-parts(d::iid; N=1000) = map(1:d.size) do j parts(d.dist) end
+parts(d::iid, N::Int=DEFAULT_SAMPLE_SIZE) = map(1:d.size) do j parts(d.dist, N) end
 # size
 # dist 
 
@@ -58,29 +57,29 @@ parts(d::iid; N=1000) = map(1:d.size) do j parts(d.dist) end
 # promote_rule(::Type{B}, ::Type{A}) where {A <: Real, B <: AbstractParticles{T,N}} where {T} = AbsractParticles{promote_type(A,T),N} where {N}
 
 
-@inline function particles(m::JointDistribution)
-    return _particles(getmoduletypencoding(m.model), m.model, m.args)
+@inline function particles(m::JointDistribution, N::Int=DEFAULT_SAMPLE_SIZE)
+    return _particles(getmoduletypencoding(m.model), m.model, m.args, Val(N))
 end
 
-@gg M function _particles(_::Type{M}, _m::Model, _args) where M <: TypeLevel{Module}
+@gg M function _particles(_::Type{M}, _m::Model, _args, _n::Val{_N}) where {M <: TypeLevel{Module},_N}
     Expr(:let,
         Expr(:(=), :M, from_type(M)),
-        type2model(_m) |> sourceParticles() |> loadvals(_args, NamedTuple()))
+        sourceParticles()(type2model(_m), _n) |> loadvals(_args, NamedTuple()))
 end
 
-@gg M function _particles(_::Type{M}, _m::Model, _args::NamedTuple{()}) where M <: TypeLevel{Module}
+@gg M function _particles(_::Type{M}, _m::Model, _args::NamedTuple{()}, _n::Val{_N}) where {M <: TypeLevel{Module},_N}
     Expr(:let,
         Expr(:(=), :M, from_type(M)),
-        type2model(_m) |> sourceParticles())
+        sourceParticles()(type2model(_m), _n))
 end
 
 export sourceParticles
 function sourceParticles() 
-    function(m::Model)
         
+    function(m::Model, ::Type{Val{_N}}) where {_N}
         _m = canonical(m)
         proc(_m, st::Assign)  = :($(st.x) = $(st.rhs))
-        proc(_m, st::Sample)  = :($(st.x) = parts($(st.rhs)))
+        proc(_m, st::Sample)  = :($(st.x) = parts($(st.rhs), $_N))
         proc(_m, st::Return)  = :(return $(st.rhs))
         proc(_m, st::LineNumber) = nothing
 
diff --git a/src/symbolic/codegen.jl b/src/symbolic/codegen.jl
index f292b236..ea6c4012 100644
--- a/src/symbolic/codegen.jl
+++ b/src/symbolic/codegen.jl
@@ -28,7 +28,7 @@ function _codegen(m :: Model, expand_sums=true)
         pushfirst!(code.args, :($v = getproperty(_args, $vname)))
     end
 
-    for v in stochastic(m)
+    for v in sampled(m)
         vname = QuoteNode(v)
         pushfirst!(code.args, :($v = getproperty(_data, $vname)))
     end
diff --git a/src/symbolic/symbolic.jl b/src/symbolic/symbolic.jl
index 04531bbb..b9063939 100644
--- a/src/symbolic/symbolic.jl
+++ b/src/symbolic/symbolic.jl
@@ -80,7 +80,7 @@ end
 
 
 
-"Half" distributions
+# "Half" distributions
 for dist in [:Normal, :Cauchy]
     let half = Symbol(:Half, dist)
         @eval begin
@@ -101,7 +101,7 @@ end
 
 
 export sym
-sym(s::Symbol) = sympy.IndexedBase(s, real=true)
+sym(s::Symbol) = sympy.symbols(s, real=true)
 sym(s) = Base.convert(Sym, s)
 function sym(expr::Expr)
     @match expr begin
@@ -120,86 +120,6 @@ end
 
 
 
-# export symlogpdf
-# # function symlogpdf(m::Model)
-# #     m = canonical(m)
-
-# #     result = @q begin
-# #         ctx = Dict()
-
-# #         ℓ = zero(Soss.Sym)
-# #     end
-
-# #     exprs = []
-
-# #     for st in map(v -> findStatement(m,v), toposortvars(m))
-# #         push!(exprs, symlogpdf(st))
-# #     end
-
-
-# #     append!(result.args, exprs)
-
-
-# #     # result
-
-# #     push!(result.args, :(ctx,ℓ))
-
-# #     result
-
-# #     # expandSums(ℓ)
-# # end
-
-# # function symlogpdf(st::Soss.Sample)
-# #     d = st.rhs
-# #     x = st.x
-# #     :(ℓ += $(symlogpdf(d,x)))
-# # end
-
-
-# # function symlogpdf(st::Soss.Assign)
-# #     val = st.rhs
-# #     x = st.x
-# #     :(ctx[$(QuoteNode(x))] = Expr(:$, $val))
-# #     # :($x = $val)
-# # end
-
-
-
-
-# # function symlogpdf(d::Expr, x::Symbol)
-# #     @match d begin
-# #         :(iid($n,$dist)) => begin
-# #             j = symbols(:j, cls=sympy.Idx)
-# #             dist = sym(dist)
-# #             x = sympy.IndexedBase(x)
-# #             n = sym(n)
-# #             :(Soss.sympy.Sum(logpdf($dist,$x[$j]), ($j,1,$n)))
-# #         end
-
-# #         :(For($f, (1:$n,))) => begin
-# #             n = sym(n)
-# #             @match f begin
-# #                 :(($j,) -> begin $lineno; $dist end) => begin
-# #                             j = symbols(j) # , cls=sympy.Idx)
-# #                             # @show j
-# #                             dist = sym(dist)
-# #                             # @show dist
-# #                             x = sympy.IndexedBase(x)
-# #                             return :(Soss.sympy.Sum(logpdf($dist,$x[$j]), ($j,1,$n)))
-# #                 end
-
-
-# #                 f => begin
-# #                     @show f
-# #                     error("symlogpdf: bad argument")
-# #                 end
-# #             end
-
-# #         end
-
-# #         _ => :(logpdf($(sym(d)), $(sym(x))))
-# #     end
-# # end
 
 export expandSums
 function expandSums(s::Sym)
@@ -221,19 +141,15 @@ end
 
 function expandSum(s::Sym, limits::Sym...)
 
-    hasIdx(s) || return maybeSum(s, limits...)
-
     func = s.func
     args = s.args
     args == () && return maybeSum(s,limits...)
 
-    if func == sympy.Add
-        return func([expandSum(t, limits...) for t in args]...)
-    elseif func == sympy.Mul
-        return expandMulSum(args, limits...)
-    else
-        return sympy.Sum(s, limits...)
-    end
+
+    func == sympy.Add && return func([expandSum(t, limits...) for t in args]...)
+    func == sympy.Mul && return expandMulSum(args, limits...)
+
+    return maybeSum(s, limits...)
 
 end
 
@@ -242,10 +158,13 @@ end
 function expandMulSum(factors::NTuple{N,Sym}, limits::Sym...) where {N}
     limits == () && return prod(factors)
 
+    p = sympy.expand_mul(prod(factors))
+    p.func == sympy.Mul || return expandSum(p, limits...)
+    factors = p.args
     for fac in factors
         for lim in limits
             (ix, ixlo, ixhi) = lim.args
-            if !insym(ix, fac)
+            if ix ∉ atoms(fac)
                 inSummand = prod(allbut(factors, fac))
                 inSum = expandSum(inSummand, lim)
                 outLims = allbut(limits, lim)
@@ -257,21 +176,20 @@ function expandMulSum(factors::NTuple{N,Sym}, limits::Sym...) where {N}
 end
 
 function atoms(s::Sym)
-    result = free_symbols(s)
-    union(result, map(x -> x.args, result)...)
-end
-
-function insym(j::Sym, s::Sym)
-    j ∈ atoms(s)
-    # for t in s.args
-    #     if j==t || in(j,t)
-    #         return true
-    #     end
-    # end
-    # return false
+    s.func == sympy.Symbol && return [s]
+    s.func == sympy.Idx && return [s]
+    isempty(free_symbols(s)) && return []
+
+    s.func == sympy.Sum && begin
+        summand = s.args[1]
+        ixs = [j.args[1] for j in s.args[2:end]]
+        bounds = union([j.args[2:3] for j in s.args[2:end]]...)
+        return setdiff(union(atoms(summand), atoms.(bounds)...), ixs)
+    end 
+    result = union(map(atoms, s.args)...)
+    return result
 end
 
-hasIdx(s::Sym) = any(startswith.(getproperty.(Soss.atoms(s), :name), "_j"))
 
 function allbut(tup, x)
     result = filter(collect(tup)) do v
@@ -286,7 +204,9 @@ function maybeSum(t::Sym, limits::Sym...)
 
     for lim in limits
         (ix, ixlo, ixhi) = lim.args
-        insym(ix, t) || return maybeSum(t * (ixhi - ixlo + 1), allbut(limits, lim)...)
+        ix ∈ atoms(t) || begin
+            return maybeSum(t * (ixhi - ixlo + 1), allbut(limits, lim)...)
+        end
     end
 
     return sympy.Sum(t, limits...)
@@ -297,60 +217,37 @@ export marginal
 function marginal(ℓ,v)
     f = ℓ.func
     f == sympy.Add || return ℓ
-    newargs = filter(t -> sym(v) in t, collect(ℓ.args))
+    newargs = filter(t -> sym(v) in atoms(t), collect(ℓ.args))
     foldl(+,newargs)
 end
 
-# marginal(m::Model, v) = marginal(m |> symlogpdf, v)
-
-# # We should be able to reason about a marginal from its derivative
-# export dmarginal
-# function dmarginal(ℓ, v)
-#     @as x ℓ begin
-#         marginal(x,sym(v))
-#         diff(x, sym(v))
-#         expand(x)
-#         sympy.collect(x, sym(v))
-#     end
-# end
-
-# dmarginal(m::Model, v) = dmarginal(m |> symlogpdf, v)
-
-
-
+export score
+score(ℓ::Sym, v) = diff(ℓ, sym(v))
 
-# logpdf(Normal(sym(:μ),sym(:σ)), :x) |> SymPy.cse
+score(m::Model, v) = score(symlogpdf(m), v)
 
+marginal(m::Model, v) = marginal(m |> symlogpdf, v)
 
-
-# # macro symdist(n, dist)
-# #     p = Expr(:tuple,gensym.(Symbol.(:p,1:n)))
-# #     @q begin
-# #         $dist($p)
-# #     end
-# # end
-# # @macroexpand @symdist(4,Normal)
+symvar(st) = :($sympy.IndexedBase($(st.x)))
 
 function symvar(st::Sample)
-    st.rhs.args[1] ∈ [:For, :iid] && return IndexedBase(st.x)
-    return sym(st.x)
+    st.rhs.args[1] ∈ [:For, :iid] && return :($sympy.IndexedBase($(st.x)))
+    return :($sym($(st.x)))
 end
+    
 
 export sourceSymlogpdf
 function sourceSymlogpdf()
     function(_m::Model)
         function proc(_m, st :: Assign)
-            # :($(st.x) = $(st.rhs))
             x = st.x
             xname = QuoteNode(x)
-            :($x = $sympy.IndexedBase($xname))
+            return :($x = $sympy.IndexedBase($xname))
         end
 
         function proc(_m, st :: Sample)
-            @q begin
-                _ℓ += symlogpdf($(st.rhs), $(st.x))
+            s = :(_ℓ += symlogpdf($(st.rhs), $(symvar(st))))
             end
-        end
         proc(_m, st :: Return)     = nothing
         proc(_m, st :: LineNumber) = nothing
 
@@ -359,9 +256,10 @@ function sourceSymlogpdf()
                 _ℓ = 0.0
             end
 
-            for x in variables(_m)
+            for x in arguments(_m)
                 xname = QuoteNode(x)
-                push!(q.args, :($x = $sympy.IndexedBase($xname)))
+                xsym = :($sympy.IndexedBase($xname))
+                push!(q.args, :($x = $xsym))
             end
 
             for st in map(v -> findStatement(_m,v), toposortvars(_m))
@@ -374,7 +272,7 @@ function sourceSymlogpdf()
                     , :($sympy.IndexedBase($xname))
                     , :($sym($xname))
                 )
-                # push!(q.args, :($x = $xsym))
+                push!(q.args, :($x = $xsym))
             end
 
             @q begin
@@ -384,8 +282,6 @@ function sourceSymlogpdf()
             end
         end
 
-
-
         buildSource(_m, proc, wrap) |> flatten
     end
 end
@@ -424,12 +320,6 @@ symlogpdf(d::Cauchy, x::Sym) = symlogpdf(Cauchy(sym(d.μ),sym(d.σ)), x)
 
 symlogpdf(d::Beta, x::Sym) = symlogpdf(Beta(sym(d.α),sym(d.β)), x)
 
-# @generated function symlogpdf(d,x::Sym)
-#     quote
-#         args = propertynames(d)
-
-#     end
-# end
 
 logpdf(d::Sym, x::Sym) = symlogpdf(d,x)
 
@@ -441,12 +331,14 @@ end
 
 symlogpdf(d,x::Sym) = logpdf(d,x)
 
-function symlogpdf(m::JointDistribution)
-    return _symlogpdf(getmoduletypencoding(m.model), m.model)
+function symlogpdf(d::JointDistribution, simplify=true)
+    symlogpdf(d.model)
 end
 
-function symlogpdf(m::Model)
-    return _symlogpdf(getmoduletypencoding(m), m)
+function symlogpdf(m::Model, simplify=true)
+    s = _symlogpdf(getmoduletypencoding(m), m)
+    simplify && return foldConstants(expandSums(s))
+    return s
 end
 
 @gg M function _symlogpdf(_::Type{M}, _m::Model) where M <: TypeLevel{Module}
@@ -456,41 +348,13 @@ end
 end
 
 
-# # s = symlogpdf(normalModel).args[7].args[3]
-
-# # export fexpr
-# # fexpr = quote
-# #     f = function(μ,σ,x)
-# #         a = $(codegen(s))
-# #         return a
-# #     end
-# # end
-
-# # julia> i,j = sympy.symbols("i j", integer=True)
-# # (i, j)
-
-# # julia> x = sympy.IndexedBase("x")
-# # x
-
-# # julia> a = sympy.Sum(x[i], (i, 1, j))
-# #   j
-# #  ___
-# #  ╲
-# #   ╲   x[i]
-# #   ╱
-# #  ╱
-# #  ‾‾‾
-# # i = 1
-
-# # julia> SymPy.walk_expression(a)
-# # :(Sum(Indexed(IndexedBase(x), i), (:i, 1, :j)))
-
 
 export tolatex
 function tolatex(ℓ::SymPy.Sym)
     r = r"_j(?<num>\d+)"
     s = s"j_{\g<num>}"
-    Base.replace(sympy.latex(ℓ), r => s)
+    result = Base.replace(sympy.latex(ℓ), r => s)
+    print()
 end
 
 
diff --git a/src/transforms/predictive.jl b/src/transforms/predictive.jl
index 94634d6e..6acca0aa 100644
--- a/src/transforms/predictive.jl
+++ b/src/transforms/predictive.jl
@@ -34,4 +34,34 @@ function predictive(m::Model, x :: Symbol)
     end
 end
 
-predictive(m::Model, x::Symbol, xs...) = predictive(predictive(m,x), xs...)
\ No newline at end of file
+predictive(m::Model, x::Symbol, xs...) = predictive(predictive(m,x), xs...)
+
+
+export predict
+
+function predict(d::JointDistribution, post::Vector{NamedTuple{N,T}}) where {N,T}
+    args = d.args
+    m = d.model
+    pred = predictive(m, keys(post[1])...)
+    map(nt -> rand(pred(merge(args,nt))), post)
+end
+
+function predict(m::Model, post::Vector{NamedTuple{N,T}}) where {N,T}
+    pred = predictive(m, keys(post[1])...)
+    map(nt -> rand(pred(nt)), post)
+end
+
+
+# TODO: These don't yet work properly t on particles
+
+function predict(d::JointDistribution, post::NamedTuple{N,T}) where {N,T}
+    args = d.args
+    m = d.model
+    pred = predictive(m, keys(post)...)
+    rand(pred(merge(args,post)))
+end
+
+function predict(m::Model, post::NamedTuple{N,T}) where {N,T}
+    pred = predictive(m, keys(post)...)
+    rand(pred(post))
+end
diff --git a/writing/2019-JOSS/paper.bib b/writing/2019-JOSS/paper.bib
new file mode 100644
index 00000000..d8585cf1
--- /dev/null
+++ b/writing/2019-JOSS/paper.bib
@@ -0,0 +1,48 @@
+@inproceedings{Cusumano-Towner:2019,
+ author = {Cusumano-Towner, Marco F. and Saad, Feras A. and Lew, Alexander K. and Mansinghka, Vikash K.},
+ title = {Gen: A General-purpose Probabilistic Programming System with Programmable Inference},
+ booktitle = {Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation},
+ series = {PLDI 2019},
+ year = {2019},
+ isbn = {978-1-4503-6712-7},
+ location = {Phoenix, AZ, USA},
+ pages = {221--236},
+ numpages = {16},
+ url = {http://doi.acm.org/10.1145/3314221.3314642},
+ doi = {10.1145/3314221.3314642},
+ acmid = {3314642},
+ publisher = {ACM},
+ address = {New York, NY, USA},
+ keywords = {Markov chain Monte Carlo, Probabilistic programming, sequential Monte Carlo, variational inference},
+} 
+
+@software{Blaom:2019,
+  author       = {Anthony Blaom and
+                  Franz Kiraly and
+                  Thibaut Lienart and
+                  Sebastian Vollmer},
+  title        = {alan-turing-institute/MLJ.jl: v0.5.3},
+  month        = nov,
+  year         = 2019,
+  publisher    = {Zenodo},
+  version      = {v0.5.3},
+  doi          = {10.5281/zenodo.3541506},
+  url          = {https://doi.org/10.5281/zenodo.3541506}
+}
+
+@article{Meurer:2017,
+ title = {SymPy: symbolic computing in Python},
+ author = {Meurer, Aaron and Smith, Christopher P. and Paprocki, Mateusz and \v{C}ert\'{i}k, Ond\v{r}ej and Kirpichev, Sergey B. and Rocklin, Matthew and Kumar, AMiT and Ivanov, Sergiu and Moore, Jason K. and Singh, Sartaj and Rathnayake, Thilina and Vig, Sean and Granger, Brian E. and Muller, Richard P. and Bonazzi, Francesco and Gupta, Harsh and Vats, Shivam and Johansson, Fredrik and Pedregosa, Fabian and Curry, Matthew J. and Terrel, Andy R. and Rou\v{c}ka, \v{S}t\v{e}p\'{a}n and Saboo, Ashutosh and Fernando, Isuru and Kulal, Sumith and Cimrman, Robert and Scopatz, Anthony},
+ year = 2017,
+ month = jan,
+ keywords = {Python, Computer algebra system, Symbolics},
+ abstract = {
+            SymPy is an open source computer algebra system written in pure Python. It is built with a focus on extensibility and ease of use, through both interactive and programmatic applications. These characteristics have led SymPy to become a popular symbolic library for the scientific Python ecosystem. This paper presents the architecture of SymPy, a description of its features, and a discussion of select submodules. The supplementary material provide additional examples and further outline details of the architecture and features of SymPy.
+         },
+ volume = 3,
+ pages = {e103},
+ journal = {PeerJ Computer Science},
+ issn = {2376-5992},
+ url = {https://doi.org/10.7717/peerj-cs.103},
+ doi = {10.7717/peerj-cs.103}
+}
\ No newline at end of file
diff --git a/writing/2019-JOSS/paper.md b/writing/2019-JOSS/paper.md
new file mode 100644
index 00000000..0c17e885
--- /dev/null
+++ b/writing/2019-JOSS/paper.md
@@ -0,0 +1,67 @@
+---
+title: 'Soss: Code Generation for Probabilistic Programming in Julia'
+tags:
+  - Julia language
+  - probabilistic programming
+  - Bayesian statistics
+  - code generation
+  - metaprogramming
+authors:
+  - name: Chad Scherrer
+    orcid: 0000-0002-1490-0304
+    affiliation: 1 # "1, 2" # (Multiple affiliations must be quoted)
+  - name: Taine Zhao
+    affiliation: 2
+affiliations:
+ - name: Metis
+   index: 1
+ - name: Department of Computer Science, University of Tsukuba
+   index: 2
+date: 7 Dec 2019
+bibliography: paper.bib
+---
+
+# Summary
+
+Probabilistic programming is a rapidly growing field, but is still far from mainstream use, due at least in part to a common diconnect between performance and ease of use. Soss aims to achieve the best of both worlds, by offering a simple mathematical syntax and specialized code generation behind the scenes.
+
+For example, here's a simple Gaussian model:
+
+```julia
+m = @model sigma,N begin
+    mu ~ Cauchy(0,1)
+    y ~ For(N) do j
+            Normal(mu,sigma)
+        end
+    end
+```
+
+Given this, a user can do things like
+
+- Specify the `sigma` and `N` arguments, and "forward sample" from the model (`rand`)
+- Compute the log-density (`logpdf`)
+- Call to external inference libraries that use these or other included methods
+- Build new models from `m`, for example using a known value for `mu` or computing the Markov blanket at a node
+- Find the symbolic log-density, using `SymPy.jl`
+- Use the result of symbolic simplification to generated optimized code, often with significant performance benefits
+
+At the time of this writing, Soss can connect (through the main library or optional add-ons) with Gen `[Cusumano-Towner:2019]`,  SymPy `[Meurer:2017]`, and MLJ `[Blaom:2019]`.
+
+<!-- 
+Citations to entries in paper.bib should be in
+[rMarkdown](http://rmarkdown.rstudio.com/authoring_bibliographies_and_citations.html)
+format.
+
+For a quick reference, the following citation commands can be used:
+- `@author:2001`  ->  "Author et al. (2001)"
+- `[@author:2001]` -> "(Author et al., 2001)"
+- `[@author1:2001; @author2:2001]` -> "(Author1 et al., 2001; Author2 et al., 2002)" -->
+
+# Acknowledgements
+
+The authors are grateful for 
+
+We acknowledge contributions from Brigitta Sipocz, Syrtis Major, and Semyeong
+Oh, and support from Kathryn Johnston during the genesis of this project.
+
+# References
\ No newline at end of file