Vec assume

appveyor.yml

@@ -1,11 +1,11 @@
 environment:
   matrix:
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5.0-win64.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x64/0.5/julia-0.5.2-win64.exe"
     MINGW_DIR: mingw64
     # MINGW_URL: https://sourceforge.net/projects/mingw-w64/files/Toolchains%20targetting%20Win64/Personal%20Builds/mingw-builds/5.3.0/threads-win32/seh/x86_64-5.3.0-release-win32-seh-rt_v4-rev0.7z/download
     MINGW_URL: http://mlg.eng.cam.ac.uk/hong/x86_64-5.3.0-release-win32-seh-rt_v4-rev0.7z
     MINGW_ARCHIVE: x86_64-5.3.0-release-win32-seh-rt_v4-rev0.7z
-  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5.0-win32.exe"
+  - JULIA_URL: "https://julialang-s3.julialang.org/bin/winnt/x86/0.5/julia-0.5.2-win32.exe"
     MINGW_DIR: mingw32
     # MINGW_URL: https://sourceforge.net/projects/mingw-w64/files/Toolchains%20targetting%20Win32/Personal%20Builds/mingw-builds/5.3.0/threads-win32/dwarf/i686-5.3.0-release-win32-dwarf-rt_v4-rev0.7z/download
     MINGW_URL: http://mlg.eng.cam.ac.uk/hong/i686-5.3.0-release-win32-dwarf-rt_v4-rev0.7z
@@ -20,8 +20,8 @@ notifications:
 cache:
 # Cache large downloads to avoid network unreliability
   - "%MINGW_ARCHIVE%"
-  - C:\Users\appveyor\.julia\v0.5
-  
+  #- C:\Users\appveyor\.julia\v0.5
+
 install:
   - ps: "[System.Net.ServicePointManager]::SecurityProtocol = [System.Net.SecurityProtocolType]::Tls12"
 # Download and install MinGW
@@ -43,9 +43,7 @@ build_script:
   - C:\projects\julia\bin\julia -e "versioninfo();
       Pkg.update();
       rm(Pkg.dir(\"Turing\"), force=true, recursive=true);
-      Pkg.clone(pwd(), \"Turing\"); 
-      Pkg.pin(\"Gadfly\", v\"0.6.0\");
-      using Gadfly;
+      Pkg.clone(pwd(), \"Turing\");
       Pkg.build(\"Turing\")"
 
 test_script:

benchmarks/MoC.run.jl

@@ -8,10 +8,14 @@ include(Pkg.dir("Turing")*"/example-models/stan-models/MoC.model.jl")
 
 bench_res = tbenchmark("HMCDA(1000, 0.65, 0.3)", "nbmodel", "data=nbstandata[1]")
 bench_res[4].names = ["phi[1]", "phi[2]", "phi[3]", "phi[4]"]
-logd = build_logd("Naive Bayes", bench_res...)
+logd = build_logd("Mixture-of-Categorical", bench_res...)
 
 include(Pkg.dir("Turing")*"/benchmarks/"*"MoC-stan.run.jl")
 logd["stan"] = stan_d
 logd["time_stan"] = nb_time
 
 print_log(logd)
+
+using Requests
+import Requests: get, post, put, delete, options, FileParam
+send_log(logd)

benchmarks/benchmarkhelper.jl

@@ -23,48 +23,65 @@ build_logd(name::String, engine::String, time, mem, tchain, _) = begin
 end
 
 # Log function
-print_log(logd::Dict, monitor=[]) = begin
-  println("/=======================================================================")
-  println("| Benchmark Result for >>> $(logd["name"]) <<<")
-  println("|-----------------------------------------------------------------------")
-  println("| Overview")
-  println("|-----------------------------------------------------------------------")
-  println("| Inference Engine  : $(logd["engine"])")
-  println("| Time Used (s)     : $(logd["time"])")
+log2str(logd::Dict, monitor=[]) = begin
+  str = ""
+  str *= ("/=======================================================================") * "\n"
+  str *= ("| Benchmark Result for >>> $(logd["name"]) <<<") * "\n"
+  str *= ("|-----------------------------------------------------------------------") * "\n"
+  str *= ("| Overview") * "\n"
+  str *= ("|-----------------------------------------------------------------------") * "\n"
+  str *= ("| Inference Engine  : $(logd["engine"])") * "\n"
+  str *= ("| Time Used (s)     : $(logd["time"])") * "\n"
   if haskey(logd, "time_stan")
-    println("|   -> time by Stan : $(logd["time_stan"])")
+    str *= ("|   -> time by Stan : $(logd["time_stan"])") * "\n"
   end
-  println("| Mem Alloc (bytes) : $(logd["mem"])")
+  str *= ("| Mem Alloc (bytes) : $(logd["mem"])") * "\n"
   if haskey(logd, "turing")
-    println("|-----------------------------------------------------------------------")
-    println("| Turing Inference Result")
-    println("|-----------------------------------------------------------------------")
+    str *= ("|-----------------------------------------------------------------------") * "\n"
+    str *= ("| Turing Inference Result") * "\n"
+    str *= ("|-----------------------------------------------------------------------") * "\n"
     for (v, m) = logd["turing"]
       if isempty(monitor) || v in monitor
-        println("| >> $v <<")
-        println("| mean = $(round(m, 3))")
+        str *= ("| >> $v <<") * "\n"
+        str *= ("| mean = $(round(m, 3))") * "\n"
         if haskey(logd, "analytic") && haskey(logd["analytic"], v)
-          print("|   -> analytic = $(round(logd["analytic"][v], 3)), ")
+          str *= ("|   -> analytic = $(round(logd["analytic"][v], 3)), ")
           diff = abs(m - logd["analytic"][v])
           diff_output = "diff = $(round(diff, 3))"
           if sum(diff) > 0.2
+            # TODO: try to fix this
             print_with_color(:red, diff_output*"\n")
+            str *= (diff_output) * "\n"
           else
-            println(diff_output)
+            str *= (diff_output) * "\n"
           end
         end
         if haskey(logd, "stan") && haskey(logd["stan"], v)
-          print("|   -> Stan     = $(round(logd["stan"][v], 3)), ")
+          str *= ("|   -> Stan     = $(round(logd["stan"][v], 3)), ")
           diff = abs(m - logd["stan"][v])
           diff_output = "diff = $(round(diff, 3))"
           if sum(diff) > 0.2
+            # TODO: try to fix this
             print_with_color(:red, diff_output*"\n")
+            str *= (diff_output) * "\n"
           else
-            println(diff_output)
+            str *= (diff_output) * "\n"
           end
         end
       end
     end
   end
-  println("\\=======================================================================")
+  str *= ("\\=======================================================================") * "\n"
+end
+
+print_log(logd::Dict, monitor=[]) = print(log2str(logd, monitor))
+
+send_log(logd::Dict, monitor=[]) = begin
+  log_str = log2str(logd, monitor)
+  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"))"
+  post("http://80.85.86.210:1110"; files = [FileParam(log_str, "text","upfile","benchmark-$time_str-$commit_str-$(logd["name"]).txt")])
 end

benchmarks/install_deps.jl

@@ -8,3 +8,5 @@ Pkg.add("UnicodePlots")
 
 Pkg.add("HDF5")
 Pkg.add("JLD")
+
+Pkg.add("Requests")

benchmarks/lda-stan.run.jl

@@ -7,7 +7,12 @@ include(Pkg.dir("Turing")*"/example-models/stan-models/lda-stan.data.jl")
 include(Pkg.dir("Turing")*"/example-models/stan-models/lda-stan.model.jl")
 
 stan_model_name = "LDA"
-ldastan = Stanmodel(Sample(save_warmup=true), name=stan_model_name, model=ldastanmodel, nchains=1);
+# ldastan = Stanmodel(Sample(save_warmup=true), name=stan_model_name, model=ldastanmodel, nchains=1);
+# To understand parameters, use: ?Stan.Static, ?Stan,Hmc
+ldastan = Stanmodel(Sample(algorithm=Stan.Hmc(Stan.Static(0.25),Stan.diag_e(),0.025,0.0),
+  save_warmup=true,adapt=Stan.Adapt(engaged=false)),
+  num_samples=2000, num_warmup=0, thin=1,
+  name=stan_model_name, model=ldastanmodel, nchains=1);
 
 rc, lda_stan_sim = stan(ldastan, ldastandata, CmdStanDir=CMDSTAN_HOME, summary=false);
 # lda_stan_sim.names

benchmarks/lda.run.jl

@@ -2,6 +2,9 @@ using Distributions
 using Turing
 using Stan
 
+using Requests
+import Requests: get, post, put, delete, options, FileParam
+
 include(Pkg.dir("Turing")*"/benchmarks/benchmarkhelper.jl")
 include(Pkg.dir("Turing")*"/example-models/stan-models/lda-stan.data.jl")
 include(Pkg.dir("Turing")*"/example-models/stan-models/lda.model.jl")
@@ -10,11 +13,14 @@ include(Pkg.dir("Turing")*"/benchmarks/"*"lda-stan.run.jl")
 
 setchunksize(60)
 
-for alg in ["HMC(1000, 0.25, 6)", "HMCDA(1000, 0.65, 1.5)", "NUTS(2000, 1000, 0.65)"]
-  bench_res = tbenchmark(alg, "ldamodel", "data=ldastandata[1]")
+#for alg in ["HMC(2000, 0.25, 10)", "HMCDA(1000, 0.65, 1.5)", "NUTS(2000, 1000, 0.65)"]
+
+for (modelc, modeln) in zip(["ldamodel_vec", "ldamodel"], ["LDA-vec", "LDA"])
+  bench_res = tbenchmark("HMC(2000, 0.025, 10)", modelc, "data=ldastandata[1]")
   bench_res[4].names = ["phi[1]", "phi[2]"]
-  logd = build_logd("LDA", bench_res...)
+  logd = build_logd(modeln, bench_res...)
   logd["stan"] = lda_stan_d
   logd["time_stan"] = lda_time
   print_log(logd)
+  send_log(logd)
 end

benchmarks/sv.run.jl

@@ -5,22 +5,16 @@ TPATH = Pkg.dir("Turing")
 
 include(TPATH*"/example-models/nips-2017/sv.model.jl")
 
-PR_for_distributions_is_merged = false
-
-if PR_for_distributions_is_merged
-
 using HDF5, JLD
-sv_data = load(TPATH*"/example-models/nips-2017/sv-data.jld")["data"]
+sv_data = load(TPATH*"/example-models/nips-2017/sv-data.jld.data")["data"]
 
 model_f = sv_model(data=sv_data[1])
-sample_n = 10000
+sample_n = 500
 
 setchunksize(550)
 chain_nuts = sample(model_f, NUTS(sample_n, 0.65))
-decribe(chain_nuts)
+describe(chain_nuts)
 
 setchunksize(5)
-chain_gibbs = sample(model_f, Gibbs(sample_n, PG(50,1,:h),NUTS(1000,0.65,:ϕ,:σ,:μ)))
-decribe(chain_gibbs)
-
-end
+chain_gibbs = sample(model_f, Gibbs(sample_n, PG(50,1,:h), NUTS(1000,0.65,:ϕ,:σ,:μ)))
+describe(chain_gibbs)

example-models/nips-2017/sv.model.jl

@@ -4,9 +4,6 @@
   σ ~ Truncated(Cauchy(0,5), 0, +Inf)
   μ ~ Cauchy(0, 10)
   h = tzeros(Real, T)
-  if σ / sqrt(1 - ϕ^2) <= 0
-    println(σ, ϕ)
-  end
   h[1] ~ Normal(μ, σ / sqrt(1 - ϕ^2))
   y[1] ~ Normal(0, exp(h[1] / 2))
   for t = 2:T

example-models/stan-models/lda.model.jl

@@ -44,10 +44,28 @@
   # end
 
   phi_dot_theta = [log([dot(map(p -> p[i], phi), theta[m]) for i = 1:V]) for m=1:M]
-  for n = 1:N
-    # phi_dot_theta = [dot(map(p -> p[i], phi), theta[doc[n]]) for i = 1:V]
-    # w[n] ~ Categorical(phi_dot_theta)
-    Turing.acclogp!(vi, phi_dot_theta[doc[n]][w[n]])
-  end
+  #for n = 1:N
+  #  # phi_dot_theta = [dot(map(p -> p[i], phi), theta[doc[n]]) for i = 1:V]
+  #  # w[n] ~ Categorical(phi_dot_theta)
+  #  Turing.acclogp!(vi, phi_dot_theta[doc[n]][w[n]])
+  #end
+  lp = mapreduce(n->phi_dot_theta[doc[n]][w[n]], +, 1:N)
+  Turing.acclogp!(vi, lp)
+
+end
+
+
+@model ldamodel_vec(K, V, M, N, w, doc, beta, alpha) = begin
+  theta = Matrix{Real}(K, M)
+  theta ~ [Dirichlet(alpha)]
+
+  phi = Matrix{Real}(V, K)
+  phi ~ [Dirichlet(beta)]
 
+  phi_dot_theta = log(phi * theta)
+  #for n = 1:N
+  #  Turing.acclogp!(vi, phi_dot_theta[w[n], doc[n]])
+  #end
+  lp = mapreduce(n->phi_dot_theta[w[n], doc[n]], +, 1:N)
+  Turing.acclogp!(vi, lp)
 end

src/core/io.jl

@@ -192,7 +192,7 @@ end
 save!(c::Chain, spl::Sampler, model::Function, vi::VarInfo) = begin
   c.info[:spl] = spl
   c.info[:model] = model
-  c.info[:vi] = vi
+  c.info[:vi] = deepcopy(vi)
 end
 
 resume(c::Chain, n_iter::Int) = begin

src/core/util.jl

@@ -25,9 +25,8 @@ macro VarName(ex::Union{Expr, Symbol})
   end
 end
 
-invlogit(x::Real) = one(x) / (one(x) + exp(-x))
-
-logit(x::Real) = log(x / (one(x) - x))
+invlogit{T<:Real}(x::Union{T,Vector{T},Matrix{T}}) = one(T) ./ (one(T) + exp(-x))
+logit{T<:Real}(x::Union{T,Vector{T},Matrix{T}}) = log(x ./ (one(T) - x))
 
 # More stable, faster version of rand(Categorical)
 function randcat(p::Vector{Float64})

src/core/varinfo.jl

@@ -143,7 +143,7 @@ end
 Base.getindex(vi::VarInfo, vns::Vector{VarName}) = begin
   @assert haskey(vi, vns[1]) "[Turing] attempted to replay unexisting variables in VarInfo"
   dist = getdist(vi, vns[1])
-  rs = reconstruct(dist, getval(vi, vns))
+  rs = reconstruct(dist, getval(vi, vns), length(vns))
   rs = istrans(vi, vns[1]) ? invlink(dist, rs) : rs
 end
 

src/samplers/gibbs.jl

@@ -76,7 +76,10 @@ sample(model::Function, alg::Gibbs;
   end
 
   # Init parameters
-  varInfo = model(); n = spl.alg.n_iters; i_thin = 1
+  varInfo = resume_from == nothing ?
+            model() :
+            resume_from.info[:vi]
+  n = spl.alg.n_iters; i_thin = 1
 
   # Gibbs steps
   spl.info[:progress] = ProgressMeter.Progress(n, 1, "[Gibbs] Sampling...", 0)

src/samplers/hmc.jl

@@ -95,7 +95,10 @@ function sample{T<:Hamiltonian}(model::Function, alg::T;
   for i = 1:n
     samples[i] = Sample(weight, Dict{Symbol, Any}())
   end
-  vi = model()
+
+  vi = resume_from == nothing ?
+       model() :
+       deepcopy(resume_from.info[:vi])
 
   if spl.alg.gid == 0
     link!(vi, spl)
@@ -107,7 +110,7 @@ function sample{T<:Hamiltonian}(model::Function, alg::T;
   for i = 1:n
     dprintln(2, "$alg_str stepping...")
 
-    time_elapsed = @elapsed vi = step(model, spl, vi, i==1)
+    time_elapsed = @elapsed vi = step(model, spl, vi, i == 1)
     time_total += time_elapsed
 
     if spl.info[:accept_his][end]     # accepted => store the new predcits
@@ -138,6 +141,8 @@ function sample{T<:Hamiltonian}(model::Function, alg::T;
   end
   c = Chain(0, samples)       # wrap the result by Chain
   if save_state               # save state
+    # Convert vi back to X if vi is required to be saved
+    if spl.alg.gid == 0 invlink!(vi, spl) end
     save!(c, spl, model, vi)
   end
 
@@ -152,18 +157,35 @@ assume{T<:Hamiltonian}(spl::Sampler{T}, dist::Distribution, vn::VarName, vi::Var
   r
 end
 
-assume{A<:Hamiltonian,D<:Distribution}(spl::Sampler{A}, dists::Vector{D}, vn::VarName, variable::Any, vi::VarInfo) = begin
+assume{A<:Hamiltonian,D<:Distribution}(spl::Sampler{A}, dists::Vector{D}, vn::VarName, var::Any, vi::VarInfo) = begin
   @assert length(dists) == 1 "[observe] Turing only support vectorizing i.i.d distribution"
   dist = dists[1]
-  n = size(variable)[end]
+  n = size(var)[end]
 
   vns = map(i -> copybyindex(vn, "[$i]"), 1:n)
 
-  rs = vi[vns]
+  rs = vi[vns]  # NOTE: inside Turing the Julia conversion should be sticked to
 
   acclogp!(vi, sum(logpdf(dist, rs, istrans(vi, vns[1]))))
 
-  rs
+  if isa(dist, UnivariateDistribution) || isa(dist, MatrixDistribution)
+    @assert size(var) == size(rs) "[assume] variable and random number dimension unmatched"
+    var = rs
+  elseif isa(dist, MultivariateDistribution)
+    if isa(var, Vector)
+      @assert length(var) == size(rs)[2] "[assume] variable and random number dimension unmatched"
+      for i = 1:n
+        var[i] = rs[:,i]
+      end
+    elseif isa(var, Matrix)
+      @assert size(var) == size(rs) "[assume] variable and random number dimension unmatched"
+      var = rs
+    else
+      error("[Turing] unsupported variable container")
+    end
+  end
+
+  var
 end
 
 observe{A<:Hamiltonian}(spl::Sampler{A}, d::Distribution, value::Any, vi::VarInfo) =