cscherrer / soss.jl

Probabilistic programming via source rewriting

Soss

Soss is a library for probabilistic programming.

Let's look at an example. First we'll load things:

using MeasureTheory
using Soss

MeasureTheory.jl is designed specifically with PPLs like Soss in mind, though you can also use Distributions.jl.

Now for a model. Here's a linear regression:

m = @model x begin
    α ~ Lebesgue(ℝ)
    β ~ Normal()
    σ ~ Exponential()
    y ~ For(x) do xj
        Normal(α + β * xj, σ)
    end
    return y
end

Next we'll generate some fake data to work with. For x-values, let's use

x = randn(20)

Now loosely speaking, Lebesgue(ℝ) is uniform over the real numbers, so we can't really sample from it. Instead, let's transform the model and make α an argument:

julia> predα = predictive(m, :α)
@model (x, α) begin
        σ ~ Exponential()
        β ~ Normal()
        y ~ For(x) do xj
                Normal(α + β * xj, σ)
            end
        return y
    end

Now we can do

julia> y = rand(predα(x=x,α=10.0))
20-element Vector{Float64}:
 10.554133456468438
  9.378065258831002
 12.873667041657287
  8.940799408080496
 10.737189595204965
  9.500536439014208
 11.327606120726893
 10.899892855024445
 10.18488773139243
 10.386969795947177
 10.382195272387214
  8.358407507910297
 10.727173015711768
 10.452311211064654
 11.076232496702387
 11.362009520020141
  9.539433052406448
 10.61851691333643
 11.586170856832645
  9.197496058151618

Now for inference! Let's use DynamicHMC, which we have wrapped in SampleChainsDynamicHMC.

julia> using SampleChainsDynamicHMC
[ Info: Precompiling SampleChainsDynamicHMC [6d9fd711-e8b2-4778-9c70-c1dfb499d4c4]

julia> post = sample(m(x=x) | (y=y,), dynamichmc())
4000-element MultiChain with 4 chains and schema (σ = Float64, β = Float64, α = Float64)
(σ = 1.0±0.15, β = 0.503±0.26, α = 10.2±0.25)

How is Soss different from Turing?

First, a fine point: When people say "the Turing PPL" they usually mean what's technically called "DynamicPPL".

• In Soss, models are first class, and can be composed or nested. For example, you can define a model and later nest it inside another model, and inference will handle both together. DynamicPPL can also handle nested models (see this PR) though I'm not aware of a way to combine independently-defined DynamicPPL models for a single inference pass.
• Soss has been updated to use MeasureTheory.jl, though everything from Distributions.jl is still available.
• Soss allows model transformations. This can be used, for example, to easily express predictive distributions or Markov blanket as a new model.
• Most of the focus of Soss is at the syntactic level; inference works in terms of "primitives" that transform the model's abstract syntax tree (AST) to new code. This adds the same benefits as using Julia's macros and generated functions, as opposed to higher-order functions alone.
• Soss can evaluate log-densities symbolically, which can then be used to produce optimized evaluations for much faster inference. This capability is in relatively early stages, and will be made more robust in our ongoing development.
• The Soss team is much smaller than that of DynamicPPL. But I hope that will change (contributors welcome!)

Soss and DynamicPPL are both maturing and becoming more complete, so the above will change over time. It's also worth noting that we (the Turing team and I) hope to move toward a natural way of using these systems together to arrive at the best of both.

How can I get involved?

I'm glad you asked! Lots of things:

• Contribute documentation or tests
• Ask questions on Discourse or Zulip
• File issues for bugs (or other problems) or feature requests
• Use Soss in your applications, teaching, or blogging
• Get involved in other libraries in the Soss ecosystem:
  • SossMLJ
  • SossGen (needs updating)
  • SampleChains
  • SampleChainsDynamicHMC
  • TupleVectors
  • NestedTuples
  • MeasureTheory
  • SymbolicCodegen

For more details, please see the documentation.

Stargazers over time

soss.jl's Issues

Variational inference

End-to-end

Given

• A model m
• Values x0 for some subset x of the parameters. Call the rest par
• A way of comparing x0 to another instance x1 of the same shape

Compute

• mForward by observeing par and returning x
• mBack by observein x and returning par

Then...

1. Use some inference procedure to sample par1 = mBack(x0)
2. Run mForward(par1) to find x1
3. Run the provided function comparing x0 and x1

sourceXform leaves some symbols

In the cars model,

julia> cars_formula
@model speed begin
        σ ~ Uniform(0, 50)
        β ~ Normal(0, 10)
        α ~ Normal(0, 100)
        μ = α .+ β .* speed
        N = length(speed)
        dist ~ For(1:N) do n
                Normal(μ[n], σ)
            end
    end

If we try to observe β, we get an error:

julia> xform(cars_formula(speed=cars.Speed),(β=0.0,))
ERROR: MethodError: no method matching *(::Symbol, ::Int64)

The root of this is sourceXform:

julia> sourceXform((β=0.0,))(cars_formula)
quote
    _result = NamedTuple()
    σ = rand(Uniform(0, 50))
    _t = xform(Uniform(0, 50))
    _result = merge(_result, (σ = _t,))
    β = (:β,)
    α = rand(Normal(0, 100))
    _t = xform(Normal(0, 100))
    _result = merge(_result, (α = _t,))
    μ = α .+ β .* speed
    N = length(speed)
    dist = rand(For(((n,)->begin
                        #= REPL[8]:8 =#
                        Normal(μ[n], σ)
                    end), 1:N))
    _t = xform(For(((n,)->begin
                        #= REPL[8]:8 =#
                        Normal(μ[n], σ)
                    end), 1:N))
    _result = merge(_result, (dist = _t,))
    as(_result)
end

Note the β = (:β,) line, which is causing the trouble

Use ModelingToolkit.jl?

It looks right up Soss.jl's alley! As far as representing and manipulating static representations of models.

Dynamic dispatch in `logpdf`

As of faa3c24,

julia> @code_warntype logpdf(m(x=x),truth)
Variables
  #self#::Core.Compiler.Const(Distributions.logpdf, false)
  m::Soss.JointDistribution{NamedTuple{(:x,),Tuple{Array{Float64,1}}},NamedTuple{(:x,),T} where T<:Tuple,TypeEncoding(begin
    σ ~ HalfNormal()
    β ~ Normal()
    α ~ Normal()
    yhat = α .+ β .* x
    n = length(x)
    y ~ For(n) do j
            Normal(yhat[j], σ)
        end
end),TypeEncoding(Main)}
  x::NamedTuple{(:x, :σ, :β, :α, :yhat, :n, :y),Tuple{Array{Float64,1},Float64,Float64,Float64,Array{Float64,1},Int64,Array{Float64,1}}}

Body::Any
1 ─ %1 = Soss.from_type($(Expr(:static_parameter, 4)))::Module
│   %2 = Base.getproperty(m, :model)::Model{NamedTuple{(:x,),T} where T<:Tuple,TypeEncoding(begin
    σ ~ HalfNormal()
    β ~ Normal()
    α ~ Normal()
    yhat = α .+ β .* x
    n = length(x)
    y ~ For(n) do j
            Normal(yhat[j], σ)
        end
end),TypeEncoding(Main)}
│   %3 = Base.getproperty(m, :args)::NamedTuple{(:x,),Tuple{Array{Float64,1}}}
│   %4 = Soss._logpdf(%1, %2, %3, x)::Any
└──      return %4

Note the Any inferences; those are bad ;)

Oddly enough, there's no such issue if you use codegen:

julia> @code_warntype logpdf(m(x=x),truth, codegen)
Variables
  #self#::Core.Compiler.Const(Distributions.logpdf, false)
  m::Soss.JointDistribution{NamedTuple{(:x,),Tuple{Array{Float64,1}}},NamedTuple{(:x,),T} where T<:Tuple,TypeEncoding(begin
    σ ~ HalfNormal()
    β ~ Normal()
    α ~ Normal()
    yhat = α .+ β .* x
    n = length(x)
    y ~ For(n) do j
            Normal(yhat[j], σ)
        end
end),TypeEncoding(Main)}
  x::NamedTuple{(:x, :σ, :β, :α, :yhat, :n, :y),Tuple{Array{Float64,1},Float64,Float64,Float64,Array{Float64,1},Int64,Array{Float64,1}}}
  #unused#::Core.Compiler.Const(Soss.codegen, false)

Body::Float64
1 ─ %1 = Soss.from_type($(Expr(:static_parameter, 4)))::Module
│   %2 = Base.getproperty(m, :model)::Model{NamedTuple{(:x,),T} where T<:Tuple,TypeEncoding(begin
    σ ~ HalfNormal()
    β ~ Normal()
    α ~ Normal()
    yhat = α .+ β .* x
    n = length(x)
    y ~ For(n) do j
            Normal(yhat[j], σ)
        end
end),TypeEncoding(Main)}
│   %3 = Base.getproperty(m, :args)::NamedTuple{(:x,),Tuple{Array{Float64,1}}}
│   %4 = Soss.codegen(%1, %2, %3, x)::Float64
└──      return %4

Both of these use GeneralizedGenerated in very similar ways

Implement dummy types for distributions

Distributions.jl is not very flexible with respect to argument types. Implementing custom types, which are basically just wrappers for the arguments would allow more flexibility. For example, the following abuse of notation is not currently possible, since Poisson only allows scalar arguments (note that #14 is probably a better option for this)

simple_poisson = @model (α, β, traps) begin
    λ = α .+ β .* traps
    complaints ~ Poisson(λ)
end

These distribution types would also make it possible to separate the code generation for different back ends (Distributions.jl, Stan math, etc.) into different packages.

use comprehensions instead of `For(1:N)`?

Would it perhaps be more natural to use comprehensions to express multidimensional sampling? In other words, instead of writing

x ~ For(1:N) do n 
    Cauchy(0,100)
end

you would write

x ~ [Cauchy(0,100) for n=1:N]

Is there something about this construct that I'm missing?

Gradient computation

For Stan-like models, XGrad should be a good fit.

Translate Stan DiffEq code

Code is here

Clean up warnings

The warnings situation has improved thanks to baggepinnen/MonteCarloMeasurements.jl#29, but there are still a few lingering. Need to clean these up:

julia> using Soss
[ Info: Recompiling stale cache file /home/chad/.julia/compiled/v1.3/Soss/vXX0U.ji for Soss [8ce77f84-9b61-11e8-39ff-d17a774bf41c]
WARNING: Method definition invokefrozen(Any, Any, Any...) in module Soss at /home/chad/git/jl/Soss/src/utils.jl:438 overwritten at /home/chad/git/jl/Soss/src/utils.jl:444.
  ** incremental compilation may be fatally broken for this module **

WARNING: using LogDensityProblems.logdensity in module Soss conflicts with an existing identifier.
WARNING: using Lazy.flatten in module Soss conflicts with an existing identifier.
WARNING: Method definition in(SymPy.Sym, SymPy.Sym) in module SymPy at /home/chad/.julia/packages/SymPy/wd4d9/src/sets.jl:22 overwritten in module Soss at /home/chad/git/jl/Soss/src/symbolic.jl:151.
  ** incremental compilation may be fatally broken for this module **

Make `iid` take a shape instead of length

@vasishth asked on Twitter about implementation of models like these. The "varying intercepts, varying slopes" model in Listing 6 would be something like this:

# varying intercepts, varying slopes
vivs = @model N,so,subj,item  begin
    σu ~ HalfCauchy() |> iid(2)
    σw ~ HalfCauchy() |> iid(2)
    σe ~ HalfCauchy()
    u ~ Normal(0, σu) |> iid(2,J)
    w ~ Normal(0, σw) |> iid(2,J)
    rt ~ For(1:N) do i 
        μ = (β[1] 
            + u[1,subj[i]] 
            + w[1,item[i]]
            + so[i] * (β[2] 
                    + u[1,subj[2]] 
                    + w[1,item[2]]
                    )
            )
        LogNormal(μ,σe)
    end
end

Currently this doesn't work, because iid only takes a length. Let's change this by removing the Int shape methods and making this always a Tuple. For example, x ~ dist |> iid(3) should produce iid((3,),dist).

Down the road, we'll have iid methods to produce an iterator. We can either make the first argument a Union type or leave it unspecified.

Sampling variables used on rhs inside `For` seems to break GeneralizedGenerated

Starting with

using Soss
using RDatasets

cars = RDatasets.dataset("datasets", "cars")

cars_formula = @model speed begin
    σ ~ Uniform(0, 50)
    β ~ Normal(0, 10)
    α ~ Normal(0, 100)
    μ = α .+ β .* speed
    N = length(speed)
    dist ~ For(1:N) do n
            Normal(μ[n], σ)
        end
end

We get this

julia> @time nuts(cars_formula(speed=cars.Speed),(dist=cars.Dist,)) ;
 11.767596 seconds (205.76 M allocations: 7.263 GiB, 21.66% gc time)

julia> @time nuts(cars_formula(speed=cars.Speed),(β = 3.91,));
409.249730 seconds (7.01 G allocations: 285.622 GiB, 23.51% gc time)

julia> @time nuts(cars_formula(speed=cars.Speed),(α = -17.2,));
422.182354 seconds (7.46 G allocations: 303.788 GiB, 22.52% gc time)

julia> @time nuts(cars_formula(speed=cars.Speed),(σ = 15.9,));
ERROR: type Float64 has no field contents

Those times for α and β are ridiculous, but (I think) that's another question.

Here's what sourceLogpdf builds:

julia> Soss.sourceLogpdf()(cars_formula)
quote
    _ℓ = 0.0
    _ℓ += logpdf(Uniform(0, 50), σ)
    _ℓ += logpdf(Normal(0, 10), β)
    _ℓ += logpdf(Normal(0, 100), α)
    μ = α .+ β .* speed
    N = length(speed)
    _ℓ += logpdf(For(1:N) do n
                #= REPL[15]:8 =#
                Normal(μ[n], σ)
            end, dist)
    return _ℓ
end

Ducktype sampleable objects

The is a draft for adding arbitrary simulators here. Would it be possible to just treat everything on the right hand side of a tilde as implementing a rand method and just add custom types? Probably related to #42 .

HMC sampling

Actually probably the HMC-like NUTS, as implemented in DymanicHMC.jl

Leaky `M`

This works fine:

julia> m = @model J,μ begin
           x ~ For(J) do j
               Normal(μ[j],1)
           end
       end;

julia> rand(m(J=3,μ=randn(3)))
(J = 3, μ = [-1.0834777431935945, -0.25066269165681376, 1.2285812
342919462], x = [-1.2785760840079383, -0.10252931007812652, -0.83
70227365078675])

But when we try the same thing with M,

julia> m = @model M,μ begin
           x ~ For(M) do j
               Normal(μ[j],1)
           end
       end;

julia> rand(m(M=3,μ=randn(3)))
ERROR: type Int64 has no field rand
Stacktrace:
 [1] getproperty(::Int64, ::Symbol) at ./Base.jl:20
 [2] macro expansion at /home/chad/.julia/packages/GeneralizedGenerated/NDqgV/src/closure_conv.jl:121 [inlined]
 [3] _rand(::Type{TypeEncoding(Main)}, ::Model{NamedTuple{(:M, :μ),T} where T<:Tuple,TypeEncoding(begin
    x ~ For(M) do j
            Normal(μ[j], 1)
        end
end),TypeEncoding(Main)}, ::NamedTuple{(:M, :μ),Tuple{Int64,Array{Float64,1}}}) at /home/chad/.julia/packages/GeneralizedGenerated/NDqgV/src/closure_conv.jl:121
 [4] rand(::Soss.JointDistribution{NamedTuple{(:M, :μ),Tuple{Int64,Array{Float64,1}}},NamedTuple{(:M, :μ),T} where T<:Tuple,TypeEncoding(begin
    x ~ For(M) do j
            Normal(μ[j], 1)
        end
end),TypeEncoding(Main)}) at /home/chad/git/Soss.jl/src/primitives/rand.jl:7
 [5] top-level scope at REPL[23]:1

Implement log likelihood

Condition on all parameters (~ values)
Sum logpdf for all <~ values

Example in README doesn't work

I am just trying the example in the readme

hello = @model μ,x begin
              σ ~ HalfCauchy()
              x ~ Normal(μ,σ) |> iid
              end
data = (μ=1, x=[2,4,5])
nuts(hello, data=data).samples

But I get the error ERROR: type NamedTuple has no field x. Any idea?

Find a way around `invokelatest`

I've been using invokelatest for quite a while. This is definitely a hack. It evaluates in global scope, has a been of additional overhead, and is generally discouraged.

@thautwarm offered to have a look at this (THANK YOU!!!) to see if there might be a workaround.

Simplest place to start is here:
https://github.com/cscherrer/Soss.jl/blob/master/src/rand.jl

The current approach is

1. sourceRand generates the source code
2. makeRand uses invokelatest to create the function
3. rand is a little wrapper that's intended to be similar to Base.rand

For example, here's a little coin flip model:

julia> coin
@model flips begin
        pHeads ~ Beta(1, 1)
        flips ~ Bernoulli(pHeads) |> iid(20)
    end

Here's the source code:

julia> sourceRand(coin)
:(function ##rand#397(args...; kwargs...)
      @unpack () = kwargs
      pHeads = rand(Beta(1, 1))
      flips = rand(iid(20, Bernoulli(pHeads)))
      (flips = flips, pHeads = pHeads)
  end)

And sampling from it:

julia> rand(coin)
(flips = [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1], pHeads = 0.65847447186616)

Dependencies not propagated correctly in `predictive`

Consider this model:

julia> m = @model begin
       N ~ Poisson(10)
       x ~ Normal() |> iid(N)
       y ~ For(N) do j
               Normal(x,1)
           end
       end;

julia> predictive(m, :x)
@model x begin
        N ~ Poisson(10)
        y ~ For(N) do j
                Normal(x, 1)
            end
    end

Clearly in this case x determines N, though the original dependency is the reverse. Ideally we'd replace N ~ Poisson(10) with N = length(x). How to get to this is not yet clear.

Support for non-unit RealIntervals

@OdE33 noticed a model using σ ~ Uniform(0,50) throwing this error:

ERROR: LoadError: asTransform(RealInterval(0.0, 50.0)) not yet supported
Stacktrace:
 [1] error(::String) at ./error.jl:33
 [2] asTransform(::RealInterval) at /home/oliver/.julia/packages/Soss/bC3jk/src/xform.jl:76

Add parameter for choice of PRNG

Based on Tamas Papp's point here.

Some possibilities:

• Have a global default, and allow users to select others by calling a different method
• Pass PRNG as an optional keyword argument
• Have a configuration file that is used to set this option (and maybe others)

Fix log density computation

Using values outside Soss

This doesn't work:

julia> f(x) = x+1
f (generic function with 1 method)

julia> m = @model begin
           y = f(x)
           x ~ Normal(y,1)
       end;

julia> rand(m())
ERROR: UndefVarError: y not defined

I'm guessing it could be a very hard problem to fix. Maybe we should implement $ interpolation inside a @model?

Storing current module when calling `@model`

julia> m1 = @model begin
                 x ~ N()
             end;

There're use cases that we need to define custom global variables outside module Soss, and currently all global variables are treated as in module Soss, which prevents the expected behaviours.

To end this, I propose to add an extra field mod to struct Model.

JuliaStaging/GeneralizedGenerated.jl#24

Computing supports and variable transformations

Some statistical modeling approaches rely on transforming variables to ℝⁿ. For example, if y ~ Gamma() we may instead work in terms of log(y). We can reason about this example statically from the AST (just have a function on symbols that takes :Gamma to log), and in general if we had a fixed set of distributions this would be no problem. It's more difficult when we have something like y ~ f(x); we need to know x, evaluate f(x) to determine a distribution, and then find the transformation for that.

One way to address this (and the best I have thought of so far) is to a sort of abstract interpretation based on random sampling from the model. This could, for example, turn y ~ f(x) into

dist = f(x)
y = rand(dist)
xforms[:y] = getTransform(dist)

where xforms is a dictionary of the transformations that need to be done.

Because this determines transformations based on a random sample, it requires an additional restriction that supports must be static and cannot change from one execution to the next.

I'm using @thautwarm's MLStyle.jl quite a lot, and had thought a nice approach would be to use this to build a little interpreter to hook into.

But as I type this, the problem seems less challenging than it originally had. [Nice side effect of writing it down. Note to self: write more things down]

So for now I'll stick with the approach that has been working so far: generate code, @eval, and then invokelatest. There are clear problems with this -- global evaluation, execution overhead -- but we can keep hunting for a better approach and migrate everything to a new approach if an alternative becomes clear.

Parameterize `For`

TODO: edit this with more details - just a quick note to myself for now...

For a while, the For type was parameterized. At a point, this broke things and I had to remove it to get things working in time for JuliaCon.

Need to have a closer look at this, in particular because something like For{Normal{Float64}} would allow specialization of the loop we build in logpdf with no need for any symbolic manipulation. Should be able to make things crazy fast. :)

Settle syntax

For a while we had been passing observations to a model as keyword arguments, for example if we observe a and b and want to reason about the other parameters we would write

nuts(m; a=2, b=[2,3,4])

This has the advantage of not requiring many parentheses or extra keywords, but gets in the way if we have configuration parameters to adjust. So we changed it to

nuts(m, (a=2, b=[2,3,4]); kwargs...)

kwargs isn't used yet, but it's available. We should expect different models to have different configuration kwargs.

This feels "right", but I'm not yet sure how it should generalize to other computations. For example, say we want to sample from the model. We do

rand(m; a=2, b=[2,3,4])

This is different than for nuts, but also feels right - it's the usual signature for rand with some keyword arguments.

But maybe we wanted more than one sample. Then we can do

rand(m, numSamples::Int; a=2, b=[2,3,4])

This again is in line with the original methods for rand.

The original way of working around all of this was to call everything with a data keyword. So for example,

nuts(m; data=(a=2, b=[2,3,4]))

I had moved away from this because it felt clunky. But now I wonder if it might be the best way to make everything consistent. Any thoughts?

Separate the Soss compiler into its own package

Currently Soss bundles many dependencies in to make the full analysis convenient. Maybe it is possible to separate the compiler into a separate lightweight package that other packages can use more easily.

GeneralizedGenerated function not finding defined variable

Trying to sample from this model:

m = @model begin
    n = 50
    p_bad_ = 0.1
    m = 5 # reviewers for each example
    
    rcl_ ~ Uniform(0.35, 0.45)
    prc_ ~ Uniform(0.65, 0.75)

    tpr = rcl_
    fpr = p_bad_ * tpr * (1/prc_ - 1.) / (1. - p_bad_)
    tnr = 1 - fpr

    prc ~ Beta(5,2) |> iid(m)
    rcl ~ Beta(2,2) |> iid(m)
    p_bad ~ Beta(1,10)

    Good = For(1:m) do mj
        Bernoulli(p_bad * rcl[mj] * (1/prc[mj] - 1)/(1 - p_bad))
    end

    Bad = For(1:m) do mj
        Bernoulli(rcl[mj])
    end

    y ~ For(1:n) do nj
            For(1:m) do mj
                MixtureModel([Good, Bad], [1 - p_bad, p_bad])
            end
        end
end;

I get

julia> rand(m())
ERROR: UndefVarError: Bad not defined
Stacktrace:
 [1] macro expansion at /home/chad/.julia/packages/GeneralizedGenerated/X0aLF/src/runtime_funcs.jl:54 [inlined]
 [2] (::GeneralizedGenerated.RuntimeFn{list(GeneralizedGenerated.Argument(:m, nothing, GeneralizedGenerated.Unset()), GeneralizedGenerated.Argument(:nj, nothing, GeneralizedGenerated.Unset())),nil(GeneralizedGenerated.Argument),let
    (Soss).For((Soss).:(:)(1, m)) do _free, (Bad, p_bad, Good)
        #= /home/chad/.julia/packages/GeneralizedGenerated/X0aLF/src/closure_conv.jl:102 =#
        (GeneralizedGenerated.Closure){GeneralizedGenerated.RuntimeFn{list(GeneralizedGenerated.Argument(:Bad, nothing, GeneralizedGenerated.Unset()), GeneralizedGenerated.Argument(:p_bad, nothing, GeneralizedGenerated.Unset()), GeneralizedGenerated.Argument(:Good, nothing, GeneralizedGenerated.Unset()), GeneralizedGenerated.Argument(:mj, nothing, GeneralizedGenerated.Unset())),nil(GeneralizedGenerated.Argument),let
    (Soss).MixtureModel([Good, Bad], [(Soss).:-(1, p_bad), p_bad])
end}(), (typeof)((Bad, p_bad, Good))}(_free)
    end
end})(::Int64, ::Int64) at /home/chad/.julia/packages/GeneralizedGenerated/X0aLF/src/runtime_funcs.jl:47
 [3] #call#29 at /home/chad/.julia/packages/GeneralizedGenerated/X0aLF/src/closure.jl:6 [inlined]
 [4] Closure at /home/chad/.julia/packages/GeneralizedGenerated/X0aLF/src/closure.jl:6 [inlined]
 [5] iterate at ./generator.jl:47 [inlined]
 [6] _collect at ./array.jl:619 [inlined]
 [7] collect_similar at ./array.jl:548 [inlined]
 [8] map at ./abstractarray.jl:2073 [inlined]
 [9] rand at /home/chad/git/jl/Soss/src/for.jl:24 [inlined]
 [10] macro expansion at /home/chad/.julia/packages/GeneralizedGenerated/X0aLF/src/closure_conv.jl:102 [inlined]
 [11] _rand(::Model{NamedTuple{(),T} where T<:Tuple,begin
    rcl_ ~ Uniform(0.35, 0.45)
    tpr = rcl_
    prc_ ~ Uniform(0.65, 0.75)
    p_bad_ = 0.1
    p_bad ~ Beta(1, 10)
    n = 50
    m = 5
    rcl ~ Beta(2, 2) |> iid(m)
    prc ~ Beta(5, 2) |> iid(m)
    fpr = (p_bad_ * tpr * (1 / prc_ - 1.0)) / (1.0 - p_bad_)
    tnr = 1 - fpr
    Good = For(1:m) do mj
            Bernoulli((p_bad * rcl[mj] * (1 / prc[mj] - 1)) / (1 - p_bad))
        end
    Bad = For(1:m) do mj
            Bernoulli(rcl[mj])
        end
    y ~ For(1:n) do nj
            For(1:m) do mj
                MixtureModel([Good, Bad], [1 - p_bad, p_bad])
            end
        end
end}, ::NamedTuple{(),Tuple{}}) at /home/chad/.julia/packages/GeneralizedGenerated/X0aLF/src/closure_conv.jl:143
 [12] rand(::Soss.JointDistribution{NamedTuple{(),Tuple{}},NamedTuple{(),T} where T<:Tuple,begin
    rcl_ ~ Uniform(0.35, 0.45)
    tpr = rcl_
    prc_ ~ Uniform(0.65, 0.75)
    p_bad_ = 0.1
    p_bad ~ Beta(1, 10)
    n = 50
    m = 5
    rcl ~ Beta(2, 2) |> iid(m)
    prc ~ Beta(5, 2) |> iid(m)
    fpr = (p_bad_ * tpr * (1 / prc_ - 1.0)) / (1.0 - p_bad_)
    tnr = 1 - fpr
    Good = For(1:m) do mj
            Bernoulli((p_bad * rcl[mj] * (1 / prc[mj] - 1)) / (1 - p_bad))
        end
    Bad = For(1:m) do mj
            Bernoulli(rcl[mj])
        end
    y ~ For(1:n) do nj
            For(1:m) do mj
                MixtureModel([Good, Bad], [1 - p_bad, p_bad])
            end
        end
end}) at /home/chad/git/jl/Soss/src/rand.jl:8
 [13] top-level scope at REPL[43]:1

But the generated code looks like this, and Bad seems to be defined:

julia> sourceRand()(m)
quote
    rcl_ = rand(Uniform(0.35, 0.45))
    tpr = rcl_
    prc_ = rand(Uniform(0.65, 0.75))
    p_bad_ = 0.1
    p_bad = rand(Beta(1, 10))
    n = 50
    m = 5
    rcl = rand(iid(m, Beta(2, 2)))
    prc = rand(iid(m, Beta(5, 2)))
    fpr = (p_bad_ * tpr * (1 / prc_ - 1.0)) / (1.0 - p_bad_)
    tnr = 1 - fpr
    Good = For(((mj,)->begin
                    #= REPL[42]:15 =#
                    Bernoulli((p_bad * rcl[mj] * (1 / prc[mj] - 1)) / (1 - p_bad))
                end), 1:m)
    Bad = For(((mj,)->begin
                    #= REPL[42]:18 =#
                    Bernoulli(rcl[mj])
                end), 1:m)
    y = rand(For(((nj,)->begin
                        #= REPL[42]:21 =#
                        For(1:m) do mj
                            #= REPL[42]:22 =#
                            MixtureModel([Good, Bad], [1 - p_bad, p_bad])
                        end
                    end), 1:n))
    (n = n, p_bad_ = p_bad_, m = m, tpr = tpr, fpr = fpr, tnr = tnr, Good = Good, Bad = Bad, rcl_ = rcl_, prc_ = prc_, prc = prc, rcl = rcl, p_bad = p_bad, y = y)
end

Model chaining

Feed the return values of one model into the parameters of another

Differential Equations

Implement SDE, etc using DifferentialEquations.jl and (most likely) Stan-like approach

Implement `particles` for NamedTuple with array-valued field

Given

dept_id = [1,1,2,2,3,3,4,4,5,5,6,6]
male = [1,0,1,0,1,0,1,0,1,0,1,0]
applications = [825,108,560,25,325,593,417,375,191,393,373,341]
admit = [512,89,353,17,120,202,138,131,53,94,22,24]

we would like a varying intercepts model for each unique department

using Soss

varying_intercepts = @model (dept_id, applications, male) begin
    α ~ For(dept_id) do id
        Normal(0,10)
    end
    β ~ Normal(0,1)
    p = logistic.(α .+ β .* male)
    N = length(male)
    admit ~ For(1:N) do n
        Binomial(applications[n], p[n])
    end
end

errors

post = nuts(varying_intercepts(dept_id = dept_id, applications = applications, male = male), (admit = admit,)) |> particles

with output

(β = 0.00901 ± 0.97, α = ERROR: LoadError: MethodError: no method matching eps(::Type{Array{Float64,1}})
Closest candidates are:
  eps(!Matched::Type{Float16}) at float.jl:739
  eps(!Matched::Type{Float32}) at float.jl:740
  eps(!Matched::Type{Float64}) at float.jl:741
  ...
Stacktrace:
 [1] eps(::Particles{Array{Float64,1},1000}) at /home/oliver/.julia/packages/MonteCarloMeasurements/WZtSS/src/particles.jl:392

We would like this model to output the intercepts for each department.

reprex.zip

`For` on a sequence of ranges

It can be convenient to write things like

For(1:M,1:S) do (m,s)
    Normal(μ[m],σ[s])
end

This would specify a distribution over matrices, where entries are independent.

So For should be able to take a sequence of ranges.

Citation?

Is there a published paper or JOSS entry that can be cited? If so, would be great to have a CITATION.bib file in the package.

Some improvements to sampling

particles

Currently particles is used in 2 ways. The first is convert a sample into a Particle-based format. The second is to draw a sample and then convert into a Particle-based format.

I propose to split this into two functions. This will enable more intuitive use of the N keyword for sample size. It could also let us return the samples in multiple formats than just Particle.

dynamicHMC

dynamicHMC 1) conditions the JointDistribution on the data 2) wraps it in an ADgradient, calls either mcmc_with_warmup or mcmc_keep_warmup, and 3) transforms the samples. Perhaps a cleaner way to implement this is to overload mcmc_keep_warmup for JointDistribution inputs. We could still provide convenient dynamicHMC functions, but this would let a user pass a Soss model to DynamicHMC and be able to tune number of warmup steps an such. A possibly better alternative is to have a wrapper struct that implements the LogDensityProblems api that could be passed directly to DynamicHMC.

advancedHMC

advancedHMC does the same 3 steps as dynamicHMC, although it looks like it might actually do the transformation online. A third option for simplifying dynamicHMC would be to implement a single wrapper type that is compatible with AdvancedHMC.sample and DynamicHMC.mcmc_keep_warmup.

Flexible return values

Just as we can turn ~ values into parameters, we should be able to easily add values to the return value of any model.

Are recusions properly handled?

I accidentally included a recursion in a model today:

m = @model begin
    μ ~ Uniform(-2,3)
    x ~ Normal(μ,x)
end

nuts(m(), (x=0.2,)) |> particles

And it seemed to work fine. Cycles between variables currently break things because the model depends on being able to topologically sort dependencies.

There is an interpretation of this model that's consistent, it's just not clear yet whether Soss is doing that or some weird random thing.

Warnings on startup, ok to ignore?

I'm a Julia newcomer, apologies if this is a dumb question. I got a lot of warnings on loading Soss. OK to ignore them? Using julia 1.2.

julia> using Soss
[ Info: Precompiling Soss [8ce77f84-9b61-11e8-39ff-d17a774bf41c]
┌ Warning: lgamma(x::Real) is deprecated, use (logabsgamma(x))[1] instead.
│ caller = lstirling_asym(::BigFloat) at misc.jl:56
└ @ StatsFuns ~/.julia/packages/StatsFuns/2QE7p/src/misc.jl:56

WARNING: using Lazy.flatten in module Soss conflicts with an existing identifier.

julia>

julia> WARNING: both AdvancedHMC and Distributions export "Multinomial"; uses of it in module Soss must be qualified
WARNING: both SimpleGraphs and Distributions export "components"; uses of it in module Soss must be qualified
WARNING: both IterTools and DataFrames export "groupby"; uses of it in module Soss must be qualified
WARNING: both DataFrames and DataStructures export "head"; uses of it in module Soss must be qualified
WARNING: both SimpleGraphs and Distributions export "scale"; uses of it in module Soss must be qualified
WARNING: both DataFrames and NamedTupleTools export "select"; uses of it in module Soss must be qualified
WARNING: both DataFrames and SimplePosets export "stack"; uses of it in module Soss must be qualified
WARNING: both DataFrames and Lazy export "tail"; uses of it in module Soss must be qualified
julia>

Wrap `model` Expression in constructor

Wrapping will make it easy to distinguish a model from a generic expression. We can also have a different way of displaying a model, for example keeping For expressions in do form.

Sample from prior

• ~ rewritten to call rand
• <~ ignored

Precompiler: ERROR: UndefVarError: depwarn not defined

I'm unable to compile Soss or debug this error. Just playing around with Soss so I hope it's something simple. The modules are loaded using Python3 but the error seems to stem from .jl packages.

[ Info: Precompiling Soss [8ce77f84-9b61-11e8-39ff-d17a774bf41c]
WARNING: using Lazy.flatten in module Soss conflicts with an existing identifier.
WARNING: using SymPy.flatten in module Soss conflicts with an existing identifier.
WARNING: Method definition in(SymPy.Sym, SymPy.Sym) in module SymPy at /home/oliver/.julia/packages/SymPy/wd4d9/src/sets.jl:22 overwritten in module Soss at /home/oliver/.julia/packages/Soss/ZXfSY/src/symbolic.jl:139.
  ** incremental compilation may be fatally broken for this module **

ERROR: LoadError: LoadError: LoadError: UndefVarError: depwarn not defined
Stacktrace:
 [1] func_for_method_checked(::Method, ::Any) at ./reflection.jl:1032
 [2] my_code_typed(::Any, ::Any) at /home/oliver/.julia/packages/ResumableFunctions/ATGuu/src/utils.jl:63
 [3] my_code_typed(::Any) at /home/oliver/.julia/packages/ResumableFunctions/ATGuu/src/utils.jl:58
 [4] top-level scope at /home/oliver/.julia/packages/ResumableFunctions/ATGuu/src/utils.jl:35

System information:

Julia Version 1.3.0-alpha.0
Commit 6c11e7c2c4 (2019-07-23 01:46 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-9600K CPU @ 3.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)

Many thanks!

Better named for subsets of variables

I propose

• arguments for values required as inputs
• assigned for LHS of Assign (=) statements
• sampled for LHS of Sample (~) statements
• parameters(m) = union(assigned(m), sampled(m))
• variables(m) = union(arguments(m), parameters(m))

@sethaxen any feedback?

Simplify log-density

Often (but not always) we're only interested in the log-density up to an additive constant. In this case normlogpdf can have unnecessary computation. For example, for a normal distribution the definition is equivalent to

normlogpdf(μ::Real, σ::Real, x::Number) = -(abs2((x-μ)/σ) + log2π)/2 - log(σ)

Overhead from the additive constant log2π is low. But if σ is also fixed, we would often prefer

normlogpdf(μ::Real, σ::Real, x::Number) = -abs2((x-μ)/σ) / 2

Implement the TableTraits.jl interface

Ingesting data from the IterableTables interface would open a whole host of data sources from various Julia packages,

World age weirdness in `codegen`

The way SymPy is set up requires loading some code within an __init__(). Here's a simplified version of what I have:

function __init__()
    stats = PyCall.pyimport_conda("sympy.stats", "sympy")
    SymPy.import_from(stats)
    global stats = stats
    
    @eval begin
        function Distributions.Normal(μ::Sym, σ::Sym)
            println("Evaluating Normal(",μ," :: Sym, ", σ, " :: Sym)")
            stats.$dist(:dist, μ,σ) |> SymPy.density
        end

        Distributions.$dist(μ,σ) = $dist(promote(μ,σ)...)
    end    


end

For a given model, symlogpdf is a GG function that uses SymPy to find its logpdf. Somewhat surprisingly, it works just fine:

julia> m = @model begin
           x ~ Normal()
       end;

julia> truth = rand(m());

julia> symlogpdf(m)
Evaluating Normal(0.0 :: Sym, 1.00000000000000 :: Sym)
       2                       log(π)   log(2)
- 0.5⋅x  - 0.693147180559945 - ────── + ──────
                                 2        2   

The next GG function, codegen, builds code that includes a call to symlogpdf. The code isn't anything too unusual:

function codegen(m::JointDistribution,x)
    return _codegen(m.model, m.args, x)    
end

@gg function _codegen(_m::Model, _args, _data)  
    type2model(_m) |> sourceCodegen() |> loadvals(_args, _data)
end

export sourceCodegen
function sourceCodegen()
    function(m::Model)
        body = @q begin end

        for (x, rhs) in pairs(m.vals)
            push!(body.args, :($x = $rhs))
        end

        push!(body.args, codegen(symlogpdf(m)))
        return body
    end
end

There's a codegen(::Sym) method that does the heavy lifting. This breaks:

julia> codegen(m(),truth)
ERROR: MethodError: no method matching Normal(::SymPy.Sym, ::SymPy.Sym)
The applicable method may be too new: running in world age 26973, while current world is 27213.

It's very strange that the first works but the second doesn't.

@thautwarm I'm pretty stuck on this one, do you have any idea what might be going on here? All relevant code is in master.

underspecified argument types

There's a problem with xform that comes up in the "8 schools" model. Here's the model:

school8 = @model σ begin
  μ ~ Flat()
  τ ~ HalfFlat()
  η ~ Normal() |> iid(8)
  y ~ For(1:8) do j
      Normal(μ + τ * η[j], σ[j])
  end
end

Trying to xform this throws an error:

julia> xform(school8( σ=[15., 10., 16., 11.,  9., 11., 10., 18.]),( y = [28.,  8., -3.,  7., -1., 
 1., 18., 12.],))
ERROR: MethodError: no method matching as(::NamedTuple{(:τ, :μ, :η),Tuple{Nothing,Nothing,TransformVariables.ArrayTransform{TransformVariables.Identity,1}}})
Closest candidates are:
  as(::Type{Real}, ::TransformVariables.Infinity{false}, ::TransformVariables.Infinity{true}) at /home/chad/.julia/packages/TransformVariables/49XcW/src/scalar.jl:166
  as(::Type{Real}, ::Real, ::TransformVariables.Infinity{true}) at /home/chad/.julia/packages/TransformVariables/49XcW/src/scalar.jl:168
  as(::Type{Real}, ::TransformVariables.Infinity{false}, ::Real) at /home/chad/.julia/packages/TransformVariables/49XcW/src/scalar.jl:170
  ...
Stacktrace:
 [1] _xform(::Model{NamedTuple{(:σ,),T} where T<:Tuple,begin
    τ ~ HalfFlat()
    μ ~ Flat()
    η ~ Normal() |> iid(8)
    y ~ For(1:8) do j
            Normal(μ + τ * η[j], σ[j])
        end
end}, ::NamedTuple{(:σ,),Tuple{Array{Float64,1}}}, ::NamedTuple{(:y,),Tuple{Array{Float64,1}}}) at /home/chad/.julia/packages/GG/Fa6NT/src/closure_conv.jl:143
 [2] xform(::Soss.BoundModel{NamedTuple{(:σ,),T} where T<:Tuple,begin
    τ ~ HalfFlat()
    μ ~ Flat()
    η ~ Normal() |> iid(8)
    y ~ For(1:8) do j
            Normal(μ + τ * η[j], σ[j])
        end
end}, ::NamedTuple{(:y,),Tuple{Array{Float64,1}}}) at /home/chad/git/Soss.jl/src/xform.jl:13
 [3] top-level scope at REPL[35]:1

I think the problem comes from the types not being known statically:

julia> school8( σ=[15., 10., 16., 11.,  9., 11., 10., 18.]) |> typeofSoss.BoundModel{NamedTuple{(:σ,),T} where T<:Tuple,begin    τ ~ HalfFlat()
    μ ~ Flat()
    η ~ Normal() |> iid(8)
    y ~ For(1:8) do j
            Normal(μ + τ * η[j], σ[j])
        end
end}

@thautwarm I don't want all the data to be type-encoded, but I think we can change a BoundModel so its args has a more specific type than that of the corresponding Model. Does that sound right to you?

Soss + Stheno

I'm trying to make the basic features of Stheno play nicely with Soss, and struggling to figure out what I'm doing wrong. In particular, here's some code that I'm trying

using Soss, Stheno
m = @model x begin
    log_l ~ Normal(0, 1)
    f = Stheno.GP(eq(exp(-l)), Stheno.GPC())
    y ~ f(x, 0.01)
end
rand(m(x=randn(10)))

which produces the following error message:

ERROR: UndefVarError: Stheno not defined
Stacktrace:
 [1] getproperty at ./Base.jl:13 [inlined]
 [2] macro expansion at /home/wct23/.julia/packages/GeneralizedGenerated/x3uMp/src/closure_conv.jl:155 [inlined]
 [3] _rand(::Model{NamedTuple{(:x,),T} where T<:Tuple,begin
    log_l ~ Normal(0, 1)
    f = Stheno.GP(eq(exp(-l)), Stheno.GPC())
    y ~ f(x, 0.01)
end}, ::NamedTuple{(),Tuple{}}) at /home/wct23/.julia/packages/GeneralizedGenerated/x3uMp/src/closure_conv.jl:155
 [4] rand(::Soss.JointDistribution{NamedTuple{(),Tuple{}},NamedTuple{(:x,),T} where T<:Tuple,begin
    log_l ~ Normal(0, 1)
    f = Stheno.GP(eq(exp(-l)), Stheno.GPC())
    y ~ f(x, 0.01)
end}) at /home/wct23/.julia/packages/Soss/qX8ds/src/rand.jl:8
 [5] top-level scope at REPL[7]:1

Versions:
Julia: 1.3.0-rc3.0
Soss: 0.6.0
Stheno: master -- this commit

Add type parameters to distribution calls

Say we have some probabilities from a posterior

julia> p = Particles(Beta(3,5))
Part500(0.375 ± 0.161)

In posterior predictions, we might try to do something like this:

julia> Bernoulli(p)
ERROR: TypeError: non-boolean (Particles{Float64,500}) used in boolean context

The problem comes from @check_args in Distributions.jl:

struct Bernoulli{T<:Real} <: DiscreteUnivariateDistribution
    p::T

    Bernoulli{T}(p::T) where {T <: Real} = new{T}(p)
end

function Bernoulli(p::T; check_args=true) where {T <: Real}
    check_args && @check_args(Bernoulli, zero(p) <= p <= one(p))
    return Bernoulli{T}(p)
end

We can't just "always pass check_args=false, because that sometimes throws an error:

julia> Normal(post.p; check_args=false)
ERROR: MethodError: no method matching Normal(::Particles{Float64,1000}; check_args=false)

But this works just fine:

julia> Bernoulli{typeof(p)}(p)
Bernoulli{Particles{Float64,500}}(
p: 0.375 ± 0.16
)

rand, particles, etc should probably do this to allow particles to be passed in.

Readme

Hi,

in the Readme, it says

using Soss 

m = @model X begin
   β ~ Normal() |> iid(size(X,2))
   y ~ For(eachrow(X)) do x
      Normal(x' * β, 1)
   end
end;

and then

truth = rand(m(x=randn(6)));

This gives; ERROR: UndefVarError: X not defined

Error when running the example in the readme

Running the following with the current registry version of Soss and on the master branch:

using Soss, Random

m = @model X begin
    β ~ Normal() |> iid(size(X,2))
    y ~ For(eachrow(X)) do x
        Normal(x' * β, 1)
    end
end;

Random.seed!(3)
X = randn(6,2)
rand(m(X=X))

gives me the following error:

ERROR: UndefVarError: interpret not defined
Stacktrace:
 [1] type2model(::Type{Model{NamedTuple{(:X,),T} where T<:Tuple,TypeEncoding(begin
    β ~ Normal() |> iid(size(X, 2))
    y ~ For(eachrow(X)) do x
            Normal(x' * β, 1)
        end
end),TypeEncoding(Main)}}) at ~/.local/share/julia/packages/Soss/g2PnG/src/core/model.jl:35
 [2] macro expansion at ~/.local/share/julia/packages/Soss/g2PnG/src/primitives/rand.jl:16 [inlined]
 [3] #s54#71(::Any, ::Any, ::Any) at ~/.local/share/julia/packages/GeneralizedGenerated/NDqgV/src/closure_conv.jl:121
 [4] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any,N} where N) at ./boot.jl:524
 [5] rand(::Soss.JointDistribution{NamedTuple{(:X,),Tuple{Array{Float64,2}}},NamedTuple{(:X,),T} where T<:Tuple,TypeEncoding(begin
    β ~ Normal() |> iid(size(X, 2))
    y ~ For(eachrow(X)) do x
            Normal(x' * β, 1)
        end
end),TypeEncoding(Main)}) at ~/.local/share/julia/packages/Soss/g2PnG/src/primitives/rand.jl:8
 [6] top-level scope at REPL[14]:1