From c187a467fab907b36cbdacdc79b15a5011f29b0c Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Wed, 26 Jul 2017 17:39:39 +0100 Subject: [PATCH 01/13] Fix return --- src/core/compiler.jl | 5 +---- src/samplers/gibbs.jl | 5 ++++- src/samplers/hmc.jl | 6 ++++++ src/samplers/is.jl | 3 +++ src/samplers/pgibbs.jl | 4 ++++ src/samplers/smc.jl | 8 +++++++- src/samplers/support/helper.jl | 14 +++++++++++++- 7 files changed, 38 insertions(+), 7 deletions(-) diff --git a/src/core/compiler.jl b/src/core/compiler.jl index ec2580183..43438f57a 100644 --- a/src/core/compiler.jl +++ b/src/core/compiler.jl @@ -220,10 +220,7 @@ macro model(fexpr) pop!(fbody_inner.args) for v = return_ex.args @assert typeof(v) == Symbol "Returned variable ($v) name must be a symbol." - vstr = string(v) - push!(fbody_inner.args, :(vn = Turing.VarName(:ret, Symbol($vstr*"_ret"), "", 1))) - # push!(fbody_inner.args, :(haskey(vi, vn) ? Turing.setval!(vi, $v, vn) : - # push!(vi, vn, $v, Distributions.Uniform(-Inf,+Inf), -1))) + push!(fbody_inner.args, :(if sampler != nothing sampler.info[:pred][Symbol($(string(v)))] = Turing.realpart($v) end)) end end push!(fbody_inner.args, Expr(:return, :vi)) diff --git a/src/samplers/gibbs.jl b/src/samplers/gibbs.jl index 0d751e809..6b67e1b09 100644 --- a/src/samplers/gibbs.jl +++ b/src/samplers/gibbs.jl @@ -98,9 +98,10 @@ sample(model::Function, alg::Gibbs; dprintln(2, "Gibbs stepping...") time_elapsed = zero(Float64) - lp = nothing; epsilon = nothing; lf_num = nothing + lp = nothing; epsilon = nothing; lf_num = nothing; last_spl = nothing for local_spl in spl.info[:samplers] + last_spl = local_spl # if PROGRESS && haskey(spl.info, :progress) local_spl.info[:progress] = spl.info[:progress] end dprintln(2, "$(typeof(local_spl)) stepping...") @@ -123,6 +124,7 @@ sample(model::Function, alg::Gibbs; if epsilon != nothing samples[i_thin].value[:epsilon] = epsilon end if lf_num != nothing samples[i_thin].value[:lf_num] = lf_num end end + update_pred(samples[i], local_spl) i_thin += 1 end time_elapsed += time_elapsed_thin @@ -147,6 +149,7 @@ sample(model::Function, alg::Gibbs; if lp != nothing samples[i].value[:lp] = lp end if epsilon != nothing samples[i].value[:epsilon] = epsilon end if lf_num != nothing samples[i].value[:lf_num] = lf_num end + update_pred(samples[i], last_spl) end if PROGRESS diff --git a/src/samplers/hmc.jl b/src/samplers/hmc.jl index d8b465aef..d583b057e 100644 --- a/src/samplers/hmc.jl +++ b/src/samplers/hmc.jl @@ -62,6 +62,9 @@ Sampler(alg::Hamiltonian) = begin # For caching gradient info[:grad_cache] = Dict{UInt64,Vector}() + + # For explicit return + info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -120,6 +123,9 @@ function sample{T<:Hamiltonian}(model::Function, alg::T; end samples[i].value[:elapsed] = time_elapsed samples[i].value[:lf_eps] = spl.info[:wum][:ϵ][end] + + update_pred(samples[i], spl) + if PROGRESS ProgressMeter.next!(spl.info[:progress]) end end diff --git a/src/samplers/is.jl b/src/samplers/is.jl index fe2337582..f2f9b1013 100644 --- a/src/samplers/is.jl +++ b/src/samplers/is.jl @@ -32,6 +32,8 @@ end Sampler(alg::IS) = begin info = Dict{Symbol, Any}() + # For explicit return + info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -43,6 +45,7 @@ sample(model::Function, alg::IS) = begin for i = 1:n vi = model(vi=VarInfo(), sampler=spl) samples[i] = Sample(vi) + update_pred(samples[i], spl) end le = logsum(map(x->x[:lp], samples)) - log(n) diff --git a/src/samplers/pgibbs.jl b/src/samplers/pgibbs.jl index 3ce98fb17..70a7e4e52 100644 --- a/src/samplers/pgibbs.jl +++ b/src/samplers/pgibbs.jl @@ -42,6 +42,8 @@ end Sampler(alg::PG) = begin info = Dict{Symbol, Any}() info[:logevidence] = [] + # For explicit return + info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -114,6 +116,8 @@ sample(model::Function, alg::PG; time_elapsed = @elapsed vi = step(model, spl, vi) push!(samples, Sample(vi)) samples[i].value[:elapsed] = time_elapsed + update_pred(samples[i], spl) + time_total += time_elapsed if PROGRESS && spl.alg.gid == 0 diff --git a/src/samplers/smc.jl b/src/samplers/smc.jl index 142b4430c..460f06892 100644 --- a/src/samplers/smc.jl +++ b/src/samplers/smc.jl @@ -38,6 +38,8 @@ end Sampler(alg::SMC) = begin info = Dict{Symbol, Any}() + # For explicit return + info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -55,6 +57,10 @@ function sample(model::Function, alg::SMC) resample!(particles,use_replay=spl.alg.use_replay) end end - res = Chain(getsample(particles)...) + w, samples = getsample(particles) + for s = samples + update_pred(s, spl) + end + res = Chain(w, samples) end diff --git a/src/samplers/support/helper.jl b/src/samplers/support/helper.jl index b1f4be7fe..314806a6e 100644 --- a/src/samplers/support/helper.jl +++ b/src/samplers/support/helper.jl @@ -6,7 +6,8 @@ @inline realpart(d::ForwardDiff.Dual) = d.value @inline realpart(ds::Union{Vector,SubArray}) = Float64[realpart(d) for d in ds] # NOTE: the function below is assumed to return a Vector now @inline realpart!(arr::Union{Array,SubArray}, ds::Union{Array,SubArray}) = for i = 1:length(ds) arr[i] = realpart(ds[i]) end -@inline realpart(ds::Matrix) = Float64[realpart(col) for col in ds] +@inline realpart{T<:Real}(ds::Matrix{T}) = Float64[realpart(col) for col in ds] +@inline realpart(ds::Matrix{Any}) = [realpart(col) for col in ds] @inline realpart(ds::Array) = map(d -> realpart(d), ds) # NOTE: this function is not optimized @inline dualpart(d::ForwardDiff.Dual) = d.partials.values @@ -76,3 +77,14 @@ end s end + +update_pred(sample::Sample, spl::Sampler) = begin + if ~isempty(spl.info[:pred]) + for sym in keys(spl.info[:pred]) + if ~haskey(sample.value, sym) + sample.value[sym] = spl.info[:pred][sym] + end + end + end + +end From 731b47fead5626ba09ebd7a4d479a83eadcc05ec Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Wed, 26 Jul 2017 21:15:14 +0100 Subject: [PATCH 02/13] Fix explict return with corountine --- src/core/compiler.jl | 2 +- src/core/util.jl | 2 +- src/core/varinfo.jl | 5 ++++- src/samplers/gibbs.jl | 4 +--- src/samplers/hmc.jl | 4 ---- src/samplers/is.jl | 3 --- src/samplers/pgibbs.jl | 3 --- src/samplers/smc.jl | 5 ----- src/samplers/support/helper.jl | 29 ++++++++++++++++++----------- test/compiler.jl/explicit_ret.jl | 18 ++++++++++++++++++ test/runtests.jl | 2 +- 11 files changed, 44 insertions(+), 33 deletions(-) create mode 100644 test/compiler.jl/explicit_ret.jl diff --git a/src/core/compiler.jl b/src/core/compiler.jl index 43438f57a..4866a019c 100644 --- a/src/core/compiler.jl +++ b/src/core/compiler.jl @@ -220,7 +220,7 @@ macro model(fexpr) pop!(fbody_inner.args) for v = return_ex.args @assert typeof(v) == Symbol "Returned variable ($v) name must be a symbol." - push!(fbody_inner.args, :(if sampler != nothing sampler.info[:pred][Symbol($(string(v)))] = Turing.realpart($v) end)) + push!(fbody_inner.args, :(if sampler != nothing vi.pred[Symbol($(string(v)))] = Turing.realpart($v) end)) end end push!(fbody_inner.args, Expr(:return, :vi)) diff --git a/src/core/util.jl b/src/core/util.jl index 1e9d7b404..522146ed5 100644 --- a/src/core/util.jl +++ b/src/core/util.jl @@ -20,7 +20,7 @@ end type NotImplementedException <: Exception end # Numerically stable sum of values represented in log domain. -function logsum(xs::Vector{Float64}) +logsum{T<:Real}(xs::Vector{T}) = begin largest = maximum(xs) ys = map(x -> exp(x - largest), xs) diff --git a/src/core/varinfo.jl b/src/core/varinfo.jl index 8efe758a0..3932af88f 100644 --- a/src/core/varinfo.jl +++ b/src/core/varinfo.jl @@ -30,7 +30,7 @@ copybyindex(vn::VarName, indexing::String) = VarName(vn.csym, vn.sym, indexing, ########### type VarInfo - idcs :: Dict{VarName, Int} + idcs :: Dict{VarName,Int} vns :: Vector{VarName} ranges :: Vector{UnitRange{Int}} vals :: Vector{Vector{Real}} @@ -38,12 +38,14 @@ type VarInfo gids :: Vector{Int} trans :: Vector{Vector{Bool}} logp :: Vector{Real} + pred :: Dict{Symbol,Any} index :: Int # index of current randomness num_produce :: Int # num of produce calls from trace, each produce corresponds to an observe. VarInfo() = begin vals = Vector{Vector{Real}}(); push!(vals, Vector{Real}()) trans = Vector{Vector{Real}}(); push!(trans, Vector{Real}()) logp = Vector{Real}(); push!(logp, zero(Real)) + pred = Dict{Symbol,Any}() new( Dict{VarName, Int}(), @@ -53,6 +55,7 @@ type VarInfo Vector{Distributions.Distribution}(), Vector{Int}(), trans, logp, + pred, 0, 0 ) diff --git a/src/samplers/gibbs.jl b/src/samplers/gibbs.jl index 6b67e1b09..0a2912f8f 100644 --- a/src/samplers/gibbs.jl +++ b/src/samplers/gibbs.jl @@ -98,7 +98,7 @@ sample(model::Function, alg::Gibbs; dprintln(2, "Gibbs stepping...") time_elapsed = zero(Float64) - lp = nothing; epsilon = nothing; lf_num = nothing; last_spl = nothing + lp = nothing; epsilon = nothing; lf_num = nothing for local_spl in spl.info[:samplers] last_spl = local_spl @@ -124,7 +124,6 @@ sample(model::Function, alg::Gibbs; if epsilon != nothing samples[i_thin].value[:epsilon] = epsilon end if lf_num != nothing samples[i_thin].value[:lf_num] = lf_num end end - update_pred(samples[i], local_spl) i_thin += 1 end time_elapsed += time_elapsed_thin @@ -149,7 +148,6 @@ sample(model::Function, alg::Gibbs; if lp != nothing samples[i].value[:lp] = lp end if epsilon != nothing samples[i].value[:epsilon] = epsilon end if lf_num != nothing samples[i].value[:lf_num] = lf_num end - update_pred(samples[i], last_spl) end if PROGRESS diff --git a/src/samplers/hmc.jl b/src/samplers/hmc.jl index d583b057e..71a9d1702 100644 --- a/src/samplers/hmc.jl +++ b/src/samplers/hmc.jl @@ -63,8 +63,6 @@ Sampler(alg::Hamiltonian) = begin # For caching gradient info[:grad_cache] = Dict{UInt64,Vector}() - # For explicit return - info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -124,8 +122,6 @@ function sample{T<:Hamiltonian}(model::Function, alg::T; samples[i].value[:elapsed] = time_elapsed samples[i].value[:lf_eps] = spl.info[:wum][:ϵ][end] - update_pred(samples[i], spl) - if PROGRESS ProgressMeter.next!(spl.info[:progress]) end end diff --git a/src/samplers/is.jl b/src/samplers/is.jl index f2f9b1013..fe2337582 100644 --- a/src/samplers/is.jl +++ b/src/samplers/is.jl @@ -32,8 +32,6 @@ end Sampler(alg::IS) = begin info = Dict{Symbol, Any}() - # For explicit return - info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -45,7 +43,6 @@ sample(model::Function, alg::IS) = begin for i = 1:n vi = model(vi=VarInfo(), sampler=spl) samples[i] = Sample(vi) - update_pred(samples[i], spl) end le = logsum(map(x->x[:lp], samples)) - log(n) diff --git a/src/samplers/pgibbs.jl b/src/samplers/pgibbs.jl index 70a7e4e52..154cf2d5c 100644 --- a/src/samplers/pgibbs.jl +++ b/src/samplers/pgibbs.jl @@ -42,8 +42,6 @@ end Sampler(alg::PG) = begin info = Dict{Symbol, Any}() info[:logevidence] = [] - # For explicit return - info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -116,7 +114,6 @@ sample(model::Function, alg::PG; time_elapsed = @elapsed vi = step(model, spl, vi) push!(samples, Sample(vi)) samples[i].value[:elapsed] = time_elapsed - update_pred(samples[i], spl) time_total += time_elapsed diff --git a/src/samplers/smc.jl b/src/samplers/smc.jl index 460f06892..1a1ad940c 100644 --- a/src/samplers/smc.jl +++ b/src/samplers/smc.jl @@ -38,8 +38,6 @@ end Sampler(alg::SMC) = begin info = Dict{Symbol, Any}() - # For explicit return - info[:pred] = Dict{Symbol,Any}() Sampler(alg, info) end @@ -58,9 +56,6 @@ function sample(model::Function, alg::SMC) end end w, samples = getsample(particles) - for s = samples - update_pred(s, spl) - end res = Chain(w, samples) end diff --git a/src/samplers/support/helper.jl b/src/samplers/support/helper.jl index 314806a6e..b4814d42f 100644 --- a/src/samplers/support/helper.jl +++ b/src/samplers/support/helper.jl @@ -57,9 +57,27 @@ end for vn in keys(vi) value[sym(vn)] = realpart(vi[vn]) end + # NOTE: do we need to check if lp is 0? value[:lp] = realpart(getlogp(vi)) + + + if ~isempty(vi.pred) + for sym in keys(vi.pred) + # if ~haskey(sample.value, sym) + value[sym] = vi.pred[sym] + # end + end + # TODO: check why 1. 2. cause errors + # TODO: which one is faster? + # 1. Using empty! + # empty!(vi.pred) + # 2. Reassign an enmtpy dict + # vi.pred = Dict{Symbol,Any}() + # 3. Do nothing? + end + Sample(0.0, value) end @@ -77,14 +95,3 @@ end s end - -update_pred(sample::Sample, spl::Sampler) = begin - if ~isempty(spl.info[:pred]) - for sym in keys(spl.info[:pred]) - if ~haskey(sample.value, sym) - sample.value[sym] = spl.info[:pred][sym] - end - end - end - -end diff --git a/test/compiler.jl/explicit_ret.jl b/test/compiler.jl/explicit_ret.jl new file mode 100644 index 000000000..c3ea4f0d8 --- /dev/null +++ b/test/compiler.jl/explicit_ret.jl @@ -0,0 +1,18 @@ +using Turing +using Distributions +using Base.Test + +@model test_ex_rt() = begin + x ~ Normal(10, 1) + y ~ Normal(x / 2, 1) + z = y + 1 + x, y, z +end + +mf = test_ex_rt() + +for alg = [HMC(1000, 0.2, 3), PG(20, 1000), SMC(10000), IS(10000), Gibbs(1000, PG(20, 1, :x), HMC(1, 0.2, 3, :y))] + chn = sample(mf, alg) + @test_approx_eq_eps mean(chn[:y]) 5.0 0.1 + @test_approx_eq_eps mean(chn[:z]) 6.0 0.1 +end diff --git a/test/runtests.jl b/test/runtests.jl index cec635dd1..a5ff90d9c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,7 +17,7 @@ testcases = Dict( "compiler.jl" => ["assume", "observe", "predict", "sample", "beta_binomial", "noparam", #"opt_param_of_dist", - "new_grammar", "newinterface", "noreturn", "forbid_global",], + "new_grammar", "newinterface", "noreturn", "forbid_global", "explicit_ret",], "container.jl" => ["copy_particle_container",], "varinfo.jl" => ["replay", "test_varname", "varinfo", "is_inside",], "io.jl" => ["chain_utility", "save_resume_chain",], From eefa8ca5c4ed4c63d9707d59e3c0c5b791243d6d Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Wed, 26 Jul 2017 21:18:24 +0100 Subject: [PATCH 03/13] Add helper for TArray --- src/samplers/support/helper.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/samplers/support/helper.jl b/src/samplers/support/helper.jl index b4814d42f..54e719a76 100644 --- a/src/samplers/support/helper.jl +++ b/src/samplers/support/helper.jl @@ -9,6 +9,7 @@ @inline realpart{T<:Real}(ds::Matrix{T}) = Float64[realpart(col) for col in ds] @inline realpart(ds::Matrix{Any}) = [realpart(col) for col in ds] @inline realpart(ds::Array) = map(d -> realpart(d), ds) # NOTE: this function is not optimized +@inline realpart(ds::TArray) = realpart(Array(ds)) @inline dualpart(d::ForwardDiff.Dual) = d.partials.values @inline dualpart(ds::Union{Array,SubArray}) = map(d -> dualpart(d), ds) From 06afdbdbc4c24d76c076e6cd56d42e5280279be3 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Wed, 26 Jul 2017 21:44:50 +0100 Subject: [PATCH 04/13] Update explicit ret test --- test/compiler.jl/explicit_ret.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/test/compiler.jl/explicit_ret.jl b/test/compiler.jl/explicit_ret.jl index c3ea4f0d8..bb4bda513 100644 --- a/test/compiler.jl/explicit_ret.jl +++ b/test/compiler.jl/explicit_ret.jl @@ -6,13 +6,16 @@ using Base.Test x ~ Normal(10, 1) y ~ Normal(x / 2, 1) z = y + 1 + x = x - 1 x, y, z end mf = test_ex_rt() -for alg = [HMC(1000, 0.2, 3), PG(20, 1000), SMC(10000), IS(10000), Gibbs(1000, PG(20, 1, :x), HMC(1, 0.2, 3, :y))] +for alg = [HMC(2000, 0.2, 3), PG(20, 2000), SMC(10000), IS(10000), Gibbs(2000, PG(20, 1, :x), HMC(1, 0.2, 3, :y))] + println("[explicit_ret.jl] testing $alg") chn = sample(mf, alg) - @test_approx_eq_eps mean(chn[:y]) 5.0 0.1 - @test_approx_eq_eps mean(chn[:z]) 6.0 0.1 + @test_approx_eq_eps mean(chn[:x]) 9.0 0.2 + @test_approx_eq_eps mean(chn[:y]) 5.0 0.2 + @test_approx_eq_eps mean(chn[:z]) 6.0 0.2 end From f4ca7bfc8a63e5a6825ec272e7dffed7be623b31 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 27 Jul 2017 15:47:12 +0100 Subject: [PATCH 05/13] Use online database for logging --- benchmarks/benchmarkhelper.jl | 12 ++++++++++-- benchmarks/bernoulli.run.jl | 2 ++ benchmarks/binormal.run.jl | 2 ++ benchmarks/gauss.run.jl | 2 ++ benchmarks/gdemo.run.jl | 2 ++ benchmarks/kid.run.jl | 1 + benchmarks/normal-mixture.run.jl | 2 ++ benchmarks/school8.run.jl | 2 ++ 8 files changed, 23 insertions(+), 2 deletions(-) diff --git a/benchmarks/benchmarkhelper.jl b/benchmarks/benchmarkhelper.jl index 3b53c0273..c7cf802e5 100644 --- a/benchmarks/benchmarkhelper.jl +++ b/benchmarks/benchmarkhelper.jl @@ -83,8 +83,16 @@ end print_log(logd::Dict, monitor=[]) = print(log2str(logd, monitor)) send_log(logd::Dict, monitor=[]) = begin - log_str = log2str(logd, monitor) - send_str(log_str, logd["name"]) + # log_str = log2str(logd, monitor) + # send_str(log_str, logd["name"]) + dir_old = pwd() + cd(Pkg.dir("Turing")) + commit_str = replace(split(readstring(pipeline(`git show --summary `, `grep "commit"`)), " ")[2], "\n", "") + cd(dir_old) + time_str = "$(Dates.format(now(), "dd-u-yyyy-HH-MM-SS"))" + logd["created"] = time_str + logd["commit"] = commit_str + post("https://api.mlab.com/api/1/databases/benchmark/collections/log?apiKey=Hak1H9--KFJz7aAx2rAbNNgub1KEylgN"; json=logd) end send_str(str::String, fname::String) = begin diff --git a/benchmarks/bernoulli.run.jl b/benchmarks/bernoulli.run.jl index 99ab56fa1..e0b899d2d 100644 --- a/benchmarks/bernoulli.run.jl +++ b/benchmarks/bernoulli.run.jl @@ -6,6 +6,8 @@ include(Pkg.dir("Turing")*"/benchmarks/benchmarkhelper.jl") include(Pkg.dir("Turing")*"/example-models/stan-models/bernoulli-stan.data.jl") include(Pkg.dir("Turing")*"/example-models/stan-models/bernoulli.model.jl") +tbenchmark("HMC(10, 0.25, 5)", "bermodel", "data=berstandata[1]") + bench_res = tbenchmark("HMC(1000, 0.25, 5)", "bermodel", "data=berstandata[1]") logd = build_logd("Bernoulli Model", bench_res...) diff --git a/benchmarks/binormal.run.jl b/benchmarks/binormal.run.jl index c1e4de83d..a19d80f0b 100644 --- a/benchmarks/binormal.run.jl +++ b/benchmarks/binormal.run.jl @@ -5,6 +5,8 @@ using Turing include(Pkg.dir("Turing")*"/benchmarks/benchmarkhelper.jl") include(Pkg.dir("Turing")*"/example-models/stan-models/binormal.model.jl") +tbenchmark("HMC(20, 0.5, 5)", "binormal", "") + bench_res = tbenchmark("HMC(2000, 0.5, 5)", "binormal", "") # chn = sample(binormal(), HMC(2000,0.5,5)) diff --git a/benchmarks/gauss.run.jl b/benchmarks/gauss.run.jl index f89b327c4..cc6fcd249 100644 --- a/benchmarks/gauss.run.jl +++ b/benchmarks/gauss.run.jl @@ -6,6 +6,8 @@ include(Pkg.dir("Turing")*"/benchmarks/benchmarkhelper.jl") include(Pkg.dir("Turing")*"/example-models/benchmarks/gdemo-stan.data.jl") include(Pkg.dir("Turing")*"/example-models/benchmarks/gdemo.model.jl") +tbenchmark("HMC(20, 0.1, 3)", "simplegaussmodel", "data=simplegaussstandata[1]") + bench_res = tbenchmark("HMC(2000, 0.1, 3)", "simplegaussmodel", "data=simplegaussstandata[1]") logd = build_logd("Simple Gaussian Model", bench_res...) logd["analytic"] = Dict("s" => 49/24, "m" => 7/6) diff --git a/benchmarks/gdemo.run.jl b/benchmarks/gdemo.run.jl index 8f93de4d3..c2cdd9f89 100644 --- a/benchmarks/gdemo.run.jl +++ b/benchmarks/gdemo.run.jl @@ -6,6 +6,8 @@ include(Pkg.dir("Turing")*"/benchmarks/benchmarkhelper.jl") include(Pkg.dir("Turing")*"/example-models/benchmarks/gauss.data.jl") include(Pkg.dir("Turing")*"/example-models/benchmarks/gauss.model.jl") +tbenchmark("PG(20, 20)", "gaussmodel", "gaussdata") + bench_res = tbenchmark("PG(20, 2000)", "gaussmodel", "gaussdata") logd = build_logd("Gaussian Model", bench_res...) print_log(logd) diff --git a/benchmarks/kid.run.jl b/benchmarks/kid.run.jl index 2586f42be..16b7a4161 100644 --- a/benchmarks/kid.run.jl +++ b/benchmarks/kid.run.jl @@ -200,6 +200,7 @@ using Turing end # chn = sample(kid_turing(data=kiddata[1]), HMC(2000, 0.0025, 10)) +tbenchmark("HMC(20, 0.0025, 10)", "kid_turing", "data=kiddata[1]") bench_res = tbenchmark("HMC(2000, 0.0025, 10)", "kid_turing", "data=kiddata[1]") chn = bench_res[4] logd = build_logd("Kid", bench_res...) diff --git a/benchmarks/normal-mixture.run.jl b/benchmarks/normal-mixture.run.jl index 5b2b711e3..7533bdc22 100644 --- a/benchmarks/normal-mixture.run.jl +++ b/benchmarks/normal-mixture.run.jl @@ -7,6 +7,8 @@ include(Pkg.dir("Turing")*"/example-models/stan-models/normal-mixture-stan.data. include(Pkg.dir("Turing")*"/example-models/stan-models/normal-mixture.model.jl") # NOTE: I only run a sub-set of the data as running the whole is quite slow +tbenchmark("Gibbs(10, HMC(1, 0.05, 1, :theta), PG(50, 1, :k), HMC(1, 0.2, 3, :mu))", "nmmodel", "simplenormalmixturestandata[1][\"y\"][1:100]") + bench_res = tbenchmark("Gibbs(1000, HMC(1, 0.05, 1, :theta), PG(50, 1, :k), HMC(1, 0.2, 3, :mu))", "nmmodel", "simplenormalmixturestandata[1][\"y\"][1:100]") logd = build_logd("Simple Gaussian Mixture Model", bench_res...) diff --git a/benchmarks/school8.run.jl b/benchmarks/school8.run.jl index d2b8a79bb..923f3919e 100644 --- a/benchmarks/school8.run.jl +++ b/benchmarks/school8.run.jl @@ -10,6 +10,8 @@ delete!(data, "tau") # chn = sample(school8(data=data), HMC(2000, 0.75, 5)) +tbenchmark("HMC(20, 0.75, 5)", "school8", "data=data") + bench_res = tbenchmark("HMC(2000, 0.75, 5)", "school8", "data=data") # bench_res[4].names = ["phi[1]", "phi[2]", "phi[3]", "phi[4]"] logd = build_logd("School 8", bench_res...) From d03d99bb7610c2b9ff95bbf4c2874b85ad5e3610 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Thu, 27 Jul 2017 20:00:13 +0100 Subject: [PATCH 06/13] Update benchmark script --- benchmarks/benchmark.jl | 58 ++++++++--------------------------- benchmarks/benchmarkhelper.jl | 44 ++++++++++++++++++++++++++ 2 files changed, 56 insertions(+), 46 deletions(-) diff --git a/benchmarks/benchmark.jl b/benchmarks/benchmark.jl index 0e3769242..f051fe0bc 100644 --- a/benchmarks/benchmark.jl +++ b/benchmarks/benchmark.jl @@ -3,49 +3,15 @@ using Turing using Stan # NOTE: put Stan models before Turing ones if you want to compare them in print_log -CONFIG = Dict( - "model-list" => [ - "gdemo-geweke", - #"normal-loc", - "normal-mixture", - "gdemo", - "gauss", - "bernoulli", - #"negative-binomial", - "school8", - "binormal", - "kid" - ], - - "test-level" => 2 # 1 = model lang, 2 = whole interface -) - -if CONFIG["test-level"] == 1 - - println("Turing compiler test started.") - - for model in CONFIG["model-list"] - println("Tesing `$model` ... ") - include("$(model).run.jl") - println("✓") - end - - println("Turing compiler test passed.") - -elseif CONFIG["test-level"] == 2 - - println("Turing benchmarking started.") - - for model in CONFIG["model-list"] - println("Benchmarking `$model` ... ") - job = `julia -e " cd(\"$(pwd())\");include(dirname(\"$(@__FILE__)\")*\"/benchmarkhelper.jl\"); - CMDSTAN_HOME = \"$CMDSTAN_HOME\"; - using Turing, Distributions, Stan; - include(dirname(\"$(@__FILE__)\")*\"/$(model).run.jl\") "` - println(job); run(job) - println("`$model` ✓") - end - - println("Turing benchmarking completed.") - -end +model_list = ["gdemo-geweke", + #"normal-loc", + "normal-mixture", + "gdemo", + "gauss", + "bernoulli", + #"negative-binomial", + "school8", + "binormal", + "kid"] + +benchmakr_turing(model_list) diff --git a/benchmarks/benchmarkhelper.jl b/benchmarks/benchmarkhelper.jl index c7cf802e5..134b68feb 100644 --- a/benchmarks/benchmarkhelper.jl +++ b/benchmarks/benchmarkhelper.jl @@ -103,3 +103,47 @@ send_str(str::String, fname::String) = begin time_str = "$(Dates.format(now(), "dd-u-yyyy-HH-MM-SS"))" post("http://80.85.86.210:1110"; files = [FileParam(str, "text","upfile","benchmark-$time_str-$commit_str-$fname.txt")]) end + + + +# using Requests +# import Requests: get, post, put, delete, options, FileParam +# import JSON + +gen_mkd_table_for_commit(commit) = begin + # commit = "f4ca7bfc8a63e5a6825ec272e7dffed7be623b31" + api_url = "https://api.mlab.com/api/1/databases/benchmark/collections/log?q={%22commit%22:%22$commit%22}&apiKey=Hak1H9--KFJz7aAx2rAbNNgub1KEylgN" + res = get(api_url) + # print(res) + + json = JSON.parse(readstring(res)) + # json[1] + + mkd = "| Model | Turing | Stan | Ratio |\n" + mkd *= "| ----- | ------ | ---- | ----- |\n" + for log in json + modelName = log["name"] + tt, ts = log["time"], log["time_stan"] + rt = tt / ts + tt, ts, rt = round(tt, 2), round(ts, 2), round(rt, 2) + mkd *= "|$modelName|$tt|$ts|$rt|\n" + end + + mkd +end + +benchmakr_turing(model_list) = begin + println("Turing benchmarking started.") + + for model in model_list + println("Benchmarking `$model` ... ") + job = `julia -e " cd(\"$(pwd())\");include(dirname(\"$(@__FILE__)\")*\"/benchmarkhelper.jl\"); + CMDSTAN_HOME = \"$CMDSTAN_HOME\"; + using Turing, Distributions, Stan; + include(dirname(\"$(@__FILE__)\")*\"/$(model).run.jl\") "` + println(job); run(job) + println("`$model` ✓") + end + + println("Turing benchmarking completed.") +end From 71ea5374ac3027d6405550e1bb90a4a9060fec16 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Fri, 28 Jul 2017 13:44:35 +0100 Subject: [PATCH 07/13] Fix missing file --- benchmarks/benchmark.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/benchmarks/benchmark.jl b/benchmarks/benchmark.jl index f051fe0bc..613d515ac 100644 --- a/benchmarks/benchmark.jl +++ b/benchmarks/benchmark.jl @@ -1,6 +1,7 @@ using Distributions using Turing using Stan +include(Pkg.dir("Turing")*"/benchmarks/benchmarkhelper.jl") # NOTE: put Stan models before Turing ones if you want to compare them in print_log model_list = ["gdemo-geweke", From 52958ab44558fa28f73e73b3e9a1432bad66947e Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Tue, 22 Aug 2017 01:36:58 +0800 Subject: [PATCH 08/13] disable explicit return --- src/core/compiler.jl | 35 +++++++++++++++++++++++------------ test/runtests.jl | 5 +++-- 2 files changed, 26 insertions(+), 14 deletions(-) diff --git a/src/core/compiler.jl b/src/core/compiler.jl index 4866a019c..28a8eb699 100644 --- a/src/core/compiler.jl +++ b/src/core/compiler.jl @@ -206,24 +206,35 @@ macro model(fexpr) # Modify fbody, so that we always return VarInfo fbody_inner = deepcopy(fbody) - return_ex = fbody.args[end] # get last statement of defined model + + return_ex = fbody.args[end] # get last statement of defined model if typeof(return_ex) == Symbol pop!(fbody_inner.args) - vstr = string(return_ex) - push!(fbody_inner.args, :(vn = Turing.VarName(:ret, Symbol($vstr*"_ret"), "", 1))) + # NOTE: code below is commented out to disable explict return + # vstr = string(return_ex) + # push!(fbody_inner.args, :(vn = Turing.VarName(:ret, Symbol($vstr*"_ret"), "", 1))) + # NOTE: code above is commented out to disable explict return + # push!(fbody_inner.args, :(haskey(vi, vn) ? Turing.setval!(vi, $return_ex, vn) : - # push!(vi, vn, $return_ex, Distributions.Uniform(-Inf,+Inf), -1))) + # push!(vi, vn, $return_ex, Distributions.Uniform(-Inf,+Inf), -1))) elseif return_ex.head == :return || return_ex.head == :tuple - if return_ex.head == :return && typeof(return_ex.args[1])!=Symbol && return_ex.args[1].head == :tuple - return_ex = return_ex.args[1] - end pop!(fbody_inner.args) - for v = return_ex.args - @assert typeof(v) == Symbol "Returned variable ($v) name must be a symbol." - push!(fbody_inner.args, :(if sampler != nothing vi.pred[Symbol($(string(v)))] = Turing.realpart($v) end)) - end + # NOTE: code below is commented out to disable explict return + # # Turn statement from return to tuple + # if return_ex.head == :return && typeof(return_ex.args[1]) != Symbol && return_ex.args[1].head == :tuple + # return_ex = return_ex.args[1] + # end + # + # # Replace :return or :tuple statement with corresponding operations on vi + # for v = return_ex.args + # @assert typeof(v) == Symbol "Returned variable ($v) name must be a symbol." + # push!(fbody_inner.args, :(if sampler != nothing vi.pred[Symbol($(string(v)))] = Turing.realpart($v) end)) + # end + # NOTE: code above is commented out to disable explict return end - push!(fbody_inner.args, Expr(:return, :vi)) + + push!(fbody_inner.args, Expr(:return, :vi)) # always return vi in the end of function body + dprintln(1, fbody_inner) fname_inner = Symbol("$(fname)_model") diff --git a/test/runtests.jl b/test/runtests.jl index a5ff90d9c..175e90a45 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,8 +16,9 @@ testcases = Dict( "ad.jl" => ["ad1", "ad2", "ad3", "pass_dual_to_dists",], "compiler.jl" => ["assume", "observe", "predict", "sample", "beta_binomial", "noparam", - #"opt_param_of_dist", - "new_grammar", "newinterface", "noreturn", "forbid_global", "explicit_ret",], + # "opt_param_of_dist", + # "explicit_ret", + "new_grammar", "newinterface", "noreturn", "forbid_global",], "container.jl" => ["copy_particle_container",], "varinfo.jl" => ["replay", "test_varname", "varinfo", "is_inside",], "io.jl" => ["chain_utility", "save_resume_chain",], From 9da504c90487791ba375c3a82e5a2d52c3bfafd7 Mon Sep 17 00:00:00 2001 From: Kai Xu Date: Tue, 22 Aug 2017 01:37:23 +0800 Subject: [PATCH 09/13] remove old code --- src/core/compiler.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/core/compiler.jl b/src/core/compiler.jl index 28a8eb699..554eb6407 100644 --- a/src/core/compiler.jl +++ b/src/core/compiler.jl @@ -214,9 +214,6 @@ macro model(fexpr) # vstr = string(return_ex) # push!(fbody_inner.args, :(vn = Turing.VarName(:ret, Symbol($vstr*"_ret"), "", 1))) # NOTE: code above is commented out to disable explict return - - # push!(fbody_inner.args, :(haskey(vi, vn) ? Turing.setval!(vi, $return_ex, vn) : - # push!(vi, vn, $return_ex, Distributions.Uniform(-Inf,+Inf), -1))) elseif return_ex.head == :return || return_ex.head == :tuple pop!(fbody_inner.args) # NOTE: code below is commented out to disable explict return From fdb884652a6f005494ad2a3dbe9f7af5586f038f Mon Sep 17 00:00:00 2001 From: Emile Mathieu Date: Thu, 24 Aug 2017 13:17:14 +0100 Subject: [PATCH 10/13] use SMC replay function and not default one --- src/samplers/smc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/samplers/smc.jl b/src/samplers/smc.jl index 1a1ad940c..e937fa689 100644 --- a/src/samplers/smc.jl +++ b/src/samplers/smc.jl @@ -52,7 +52,7 @@ function sample(model::Function, alg::SMC) while consume(particles) != Val{:done} ess = effectiveSampleSize(particles) if ess <= spl.alg.resampler_threshold * length(particles) - resample!(particles,use_replay=spl.alg.use_replay) + resample!(particles,smc_spl.alg.resampler,use_replay=smc_spl.alg.use_replay) end end w, samples = getsample(particles) From 2e972a147a11d8bfc0b73d7924889c6f2f61b407 Mon Sep 17 00:00:00 2001 From: Emile Mathieu Date: Thu, 24 Aug 2017 17:36:41 +0100 Subject: [PATCH 11/13] typo fix --- src/samplers/smc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/samplers/smc.jl b/src/samplers/smc.jl index e937fa689..dbc23dbfa 100644 --- a/src/samplers/smc.jl +++ b/src/samplers/smc.jl @@ -52,7 +52,7 @@ function sample(model::Function, alg::SMC) while consume(particles) != Val{:done} ess = effectiveSampleSize(particles) if ess <= spl.alg.resampler_threshold * length(particles) - resample!(particles,smc_spl.alg.resampler,use_replay=smc_spl.alg.use_replay) + resample!(particles,spl.alg.resampler,use_replay=smc_spl.alg.use_replay) end end w, samples = getsample(particles) From fc3b13026049588ea980bb3b4e16980275b22929 Mon Sep 17 00:00:00 2001 From: Emile Mathieu Date: Fri, 25 Aug 2017 11:03:50 +0100 Subject: [PATCH 12/13] typo fix, again... --- src/samplers/smc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/samplers/smc.jl b/src/samplers/smc.jl index dbc23dbfa..da520bab5 100644 --- a/src/samplers/smc.jl +++ b/src/samplers/smc.jl @@ -52,7 +52,7 @@ function sample(model::Function, alg::SMC) while consume(particles) != Val{:done} ess = effectiveSampleSize(particles) if ess <= spl.alg.resampler_threshold * length(particles) - resample!(particles,spl.alg.resampler,use_replay=smc_spl.alg.use_replay) + resample!(particles,spl.alg.resampler,use_replay=spl.alg.use_replay) end end w, samples = getsample(particles) From 8d0516f0b836cc4a6de18359bf8abfaff45b34b7 Mon Sep 17 00:00:00 2001 From: Emile Mathieu Date: Fri, 25 Aug 2017 11:07:01 +0100 Subject: [PATCH 13/13] vectorize assume for PG and SMC fix --- src/samplers/pgibbs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/samplers/pgibbs.jl b/src/samplers/pgibbs.jl index 154cf2d5c..768f77c59 100644 --- a/src/samplers/pgibbs.jl +++ b/src/samplers/pgibbs.jl @@ -166,7 +166,7 @@ assume{T<:Union{PG,SMC}}(spl::Sampler{T}, dist::Distribution, vn::VarName, _::Va end end -assume{T<:Union{PG,SMC}}(spl::Void, dists::Vector{T}, vn::VarName, var::Any, vi::VarInfo) = +assume{A<:Union{PG,SMC},D<:Distribution}(spl::Sampler{A}, dists::Vector{D}, vn::VarName, var::Any, vi::VarInfo) = error("[Turing] PG and SMC doesn't support vectorizing assume statement") observe{T<:Union{PG,SMC}}(spl::Sampler{T}, dist::Distribution, value, vi) =