[WIP] Develop concept

Example/Using_Custom_VI/demonstarte.jl

@@ -0,0 +1,57 @@
+
+## It is not yet a package...
+using Distributions
+using AdvancedPS
+using Libtask
+using BenchmarkTools
+using NamedTupleTools
+include("AdvancedPS/Example/Using_Custom_VI/test_interface.jl")
+# Define a short model.
+# The syntax is rather simple. Observations need to be reported with report_observation.
+# Transitions must be reported using report_transition.
+# The trace contains the variables which we want to infer using particle gibbs.
+# Thats all!
+n = 3000
+
+y = Vector{Float64}(undef,n-1)
+for i =1:n-1
+    y[i] = 0
+end
+
+function task_f(y)
+    var = initialize()
+    x = zeros(n,1)
+    r = vectorize(Normal(), rand(Normal()))
+    x[1,:] = update_var(var, 1,  r)
+    report_transition!(var,0.0,0.0)
+    for i = 2:n
+        # Sampling
+        r =  vectorize(Normal(),rand(Normal()))
+        x[i,:] = update_var(var, i, r)
+        logγ = logpdf(Normal(),x[i,1]) #γ(x_t|x_t-1)
+        logp = logpdf(Normal(),x[i,1])                # p(x_t|x_t-1)
+        report_transition!(var,logp,logγ)
+        #Proposal and Resampling
+        logpy = logpdf(Normal(x[i,1], 1.0), y[i-1])
+        var = report_observation!(var,logpy)
+    end
+end
+
+
+
+
+
+task = create_task(task_f, y)
+model = PFModel(task)
+tcontainer =  Container(zeros(n,1),Vector{Bool}(undef,n),zeros(n),0)
+
+
+
+alg = AdvancedPS.SMCAlgorithm()
+uf = AdvancedPS.SMCUtilityFunctions(deepcopy, set_retained_vns_del_by_spl!, empty!, tonamedtuple)
+@btime sample(model, alg, uf, tcontainer, 100)
+
+
+alg = AdvancedPS.PGAlgorithm(AdvancedPS.resample_systematic, 1.0, 10)
+uf = AdvancedPS.PGUtilityFunctions(deepcopy, set_retained_vns_del_by_spl!, empty!, tonamedtuple)
+@btime chn2 =sample(model, alg, uf, tcontainer, 5)

Example/Using_Custom_VI/test_interface.jl

@@ -0,0 +1,83 @@
+## Some function which make the model easier to define.
+
+
+# A very light weight container for a state space model
+# The state space is one dimensional
+# Note that this is only for a very simple demonstration.
+
+
+
+# This is a very shallow container solely for testing puropose
+mutable struct Container
+    x::Array{Float64,2}
+    marked::Vector{Bool}
+    produced_at::Vector{Int64}
+    num_produce::Float64
+end
+
+function Base.deepcopy(vi::Container)
+    return Container(deepcopy(vi.x),deepcopy(vi.marked),deepcopy(vi.produced_at),deepcopy(vi.num_produce))
+end
+
+
+function set_retained_vns_del_by_spl!(container::Container)
+    for i in 1:length(container.marked)
+        if container.marked[i]
+            if container.produced_at[i] >= container.num_produce
+                container.marked[i] = false
+            end
+        end
+    end
+    container
+end
+
+
+
+# This is important for initalizaiton
+const initialize = current_trace
+
+function report_observation!(trace, logp::Float64)
+    produce(logp)
+    trace = current_trace()
+end
+
+# logγ corresponds to the proposal distributoin we are sampling from.
+function report_transition!(trace,logp::Float64,logγ::Float64)
+    trace.taskinfo.logp += logp - logγ
+    trace.taskinfo.logpseq += logp
+end
+
+function update_var(trace, vn::Int64, r::Vector{Float64})
+    if !trace.vi.marked[vn]
+        trace.vi.x[vn,:] = r
+        trace.vi.marked[vn] = true
+        trace.vi.produced_at[vn] = trace.vi.num_produce
+        return r
+    end
+    return trace.vi.x[vn,:]
+end
+
+
+# The reason for this is that we need to pass it!
+function copy_container(vi::Container)
+    Container(deepcopy(vi.x),deepcopy(vi.marked),deepcopy(vi.produced_at),copy(vi.num_produce))
+end
+
+function create_task(f::Function, args...)
+    return CTask(() ->  begin new_vi=f(args...); produce(Val{:done}); new_vi; end )
+end
+
+function tonamedtuple(vi::Container)
+    tnames = Tuple([Symbol("x$i") for i in 1:size(vi.x)[1]])
+    tvalues = Tuple([vi.x[i,:] for i in 1:size(vi.x)[1]])
+    return namedtuple(tnames, tvalues)
+end
+function Base.empty!(vi::Container)
+    for i in 1:length(vi.marked)
+        vi.marked[i] = false
+    end
+    vi
+end
+vectorize(d::UnivariateDistribution, r::Real) = [r]
+vectorize(d::MultivariateDistribution, r::AbstractVector{<:Real}) = copy(r)
+vectorize(d::MatrixDistribution, r::AbstractMatrix{<:Real}) = copy(vec(r))

Example/Using_Turing_VI/Turing_baseline.jl

@@ -0,0 +1,25 @@
+using Turing
+using BenchmarkTools
+n = 500
+
+y = Vector{Float64}(undef,n-1)
+for i =1:n-1
+    y[i] = 0
+end
+
+@model demo(y) = begin
+    x = TArray{Float64}(undef,n)
+    x[1] ~ Normal()
+    for i = 2:n
+        x[i] ~ Normal()
+        y[i-1] ~ Normal(x[i],1.0)
+    end
+end
+
+
+
+
+#@elapsed sample(demo(),PG(10),5)
+chn = @btime sample(demo(y), SMC(), 100)
+
+chn.na.axes

Example/Using_Turing_VI/demonstrate_with_Turing_VI.jl

@@ -0,0 +1,75 @@
+
+## It is not yet a package...
+using Distributions
+using AdvancedPS
+using BenchmarkTools
+using Libtask
+include("AdvancedPS/Example/Using_Turing_VI/turing_interface.jl")
+# Define a short model.
+# The syntax is rather simple. Observations need to be reported with report_observation.
+# Transitions must be reported using report_transition.
+# The trace contains the variables which we want to infer using particle gibbs.
+# Thats all!
+n = 500
+
+y = Vector{Float64}(undef,n-1)
+for i =1:n-1
+    y[i] = 0
+end
+
+function task_f(y)
+    var = initialize()
+    x = TArray{Float64}(undef,n)
+    vn = @varname x[1]
+    x[1] = update_var(var, vn, rand(Normal()))
+    report_transition!(var,0.0,0.0)
+    for i = 2:n
+        # Sampling
+        r =  rand(Normal())
+        vn = @varname x[i]
+        x[i] = update_var(var, vn, r)
+        logγ = logpdf(Normal(),x[i]) #γ(x_t|x_t-1)
+        logp = logγ             # p(x_t|x_t-1)
+        report_transition!(var,logp,logγ)
+        #Proposal and Resampling
+        logpy = logpdf(Normal(x[i], 1.0), y[i-1])
+        var = report_observation!(var,logpy)
+    end
+end
+
+
+
+
+
+
+
+task = create_task(task_f, y)
+model = PFModel(task)
+
+
+#################################################################
+# Get type stability!!                                          #
+#################################################################
+alg = AdvancedPS.SMCAlgorithm()
+uf = AdvancedPS.SMCUtilityFunctions(deepcopy, set_retained_vns_del_by_spl!, empty!, tonamedtuple)
+untypedcontainer = VarInfo()
+T = Trace{typeof(untypedcontainer),SMCTaskInfo{Float64}}
+particles = T[ Trace(untypedcontainer, task, SMCTaskInfo(), uf.copy) for _ =1:1]
+pc = ParticleContainer{typeof(particles[1])}(particles,zeros(1),0.0,0)
+AdvancedPS.sample!(pc, alg, uf, nothing)
+container = uf.empty!(TypedVarInfo(pc[1].vi))
+tcontainer = container
+
+
+
+
+alg = AdvancedPS.SMCAlgorithm()
+uf = AdvancedPS.SMCUtilityFunctions(deepcopy, set_retained_vns_del_by_spl!, empty!, tonamedtuple)
+@btime sample(model, alg, uf, tcontainer, 100)
+
+
+
+
+alg = AdvancedPS.PGAlgorithm(AdvancedPS.resample_systematic, 1.0, 10)
+uf = AdvancedPS.PGUtilityFunctions(deepcopy, set_retained_vns_del_by_spl!, empty!, tonamedtuple)
+@elapsed chn2 =sample(model, alg, uf, tcontainer, 5)