SciML / DiffEqBayes.jl

Extension functionality which uses Stan.jl, DynamicHMC.jl, and Turing.jl to estimate the parameters to differential equations and perform Bayesian probabilistic scientific machine learning
https://docs.sciml.ai/DiffEqBayes/stable/
Other
121 stars 29 forks source link

Parameter estimation with Turing and HMC #76

Closed denisshepelin closed 5 years ago

denisshepelin commented 5 years ago

Dear all, I'm trying to build a Turing-first simple example with estimation of ODE parameters using HMC sampler.

My environment with Julia 1.1.1 :

julia> Pkg.status()
    Status `~/.julia/environments/v1.1/Project.toml`
  [ebbdde9d] DiffEqBayes v1.1.0
  [aae7a2af] DiffEqFlux v0.5.0
  [41bf760c] DiffEqSensitivity v3.2.3
  [0c46a032] DifferentialEquations v6.4.0
  [31c24e10] Distributions v0.20.0
  [bbc10e6e] DynamicHMC v1.0.5
  [587475ba] Flux v0.8.3
  [f6369f11] ForwardDiff v0.10.3
  [7073ff75] IJulia v1.18.1
  [6fdf6af0] LogDensityProblems v0.8.2
  [65888b18] ParameterizedFunctions v4.1.1
  [d96e819e] Parameters v0.10.3
  [91a5bcdd] Plots v0.25.1
  [731186ca] RecursiveArrayTools v0.20.0
  [f3b207a7] StatsPlots v0.11.0
  [84d833dd] TransformVariables v0.3.3
  [fce5fe82] Turing v0.6.17

So far the code looks like this:

using DifferentialEquations, Plots, ParameterizedFunctions, Random, Distributions, Turing, StatsPlots

predator_prey = @ode_def begin
    du = (alpha - beta*v) * u
    dv = (-gamma + delta * u) * v
    end alpha beta gamma delta

#  "Experimental" data generation
u0 = [10.0, 10.0]
tspan = (0.0, 40.0)
p = [1.1,0.4,0.4, 0.1]
tp= collect(range(0,stop=40))
prob = ODEProblem(predator_prey, u0, tspan, p)
sol = solve(prob, saveat=tp)

# optional plot
# plot(sol, linestyle=:dot)
# modify data with adding gaussian noise

d = Normal(0, 0.5)
newsol = sol + rand(d, size(sol))

function test_f(p, tp_save)
  _prob = remake(prob;u0=convert.(eltype(p),prob.u0),p=p)
  solve(_prob,Tsit5(),saveat=tp_save;abstol=1e-6, reltol= 1e-6)
end

plot(sol)
scatter!(tp,newsol[1,:])
scatter!(tp,newsol[2,:])

# define statistical model
@model bayes_predator_prey(x, tp_x) = begin
  alpha ~ Normal(1.0, 0.2)
  beta ~ Normal(0.5, 0.2)
  gamma ~ Normal(0.5, 0.2)
  delta ~ Normal(0.1, 0.1)
  # gather parameters and solve equation
  p = [alpha, beta, gamma ,delta]
  sol_tmp = test_f(p, tp_x)
  # hope that it's likelihood
  N = length(tp_x)

  fill_length = length(tp_x) - length(sol_tmp.u)

  for i in 1:fill_length
    if eltype(sol_tmp.u) <: Number
      push!(sol_tmp.u, Inf)
    else
      push!(sol_tmp.u, fill(Inf, size(sol_tmp[1])))
    end
  end

  for i in 1:N
    x[:,i] ~ MvNormal(sol_tmp.u[i], [0.5,0.5])
  end
end

# run the sampling
chain = sample(bayes_predator_prey(newsol, tp), HMC(1000, 0.1, 5))

This code is a very close adaptation of how Turing based inference works in DiffEqBayes but with different sampler (HMC instead of IS). Unfortunately this code is not working at the stage of sampling. The error is the following:

       chain = sample(bayes_predator_prey(newsol, tp), HMC(1000, 0.1, 5))
ERROR: MethodError: no method matching OrdinaryDiffEq.Tsit5ConstantCache(::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Float64, ::Rational{Int64}, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::Rational{Int64}, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::Rational{Int64}, ::BigFloat, ::BigFloat, ::BigFloat, ::Float64, ::BigFloat, ::BigFloat, ::Float64, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::BigFloat, ::Float64, ::Int64, ::Float64)
Closest candidates are:
  OrdinaryDiffEq.Tsit5ConstantCache(::T2, ::T2, ::T2, ::T2, ::T2, ::T2, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T, ::T) where {T, T2} at /Users/denshe/.julia/packages/OrdinaryDiffEq/tat1c/src/tableaus/low_order_rk_tableaus.jl:498
Stacktrace:
 [1] OrdinaryDiffEq.Tsit5ConstantCache(::Type, ::Type) at /Users/denshe/.julia/packages/OrdinaryDiffEq/tat1c/src/tableaus/low_order_rk_tableaus.jl:647
 [2] alg_cache(::Tsit5, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}, ::Type, ::Type, ::Type, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}, ::getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr}, ::Float64, ::Float64, ::Float64, ::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}, ::Bool, ::Type{Val{true}}) at /Users/denshe/.julia/packages/OrdinaryDiffEq/tat1c/src/caches/low_order_rk_caches.jl:348
 [3] #__init#297(::Array{Int64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Bool, ::Nothing, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Float64, ::Float64, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Int64, ::Float64, ::Float64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::typeof(LinearAlgebra.opnorm), ::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), ::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/denshe/.julia/packages/OrdinaryDiffEq/tat1c/src/solve.jl:242
 [4] (::getfield(DiffEqBase, Symbol("#kw##__init")))(::NamedTuple{(:saveat, :abstol, :reltol),Tuple{Array{Int64,1},Float64,Float64}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0
 [5] #__solve#296(::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:saveat, :abstol, :reltol),Tuple{Array{Int64,1},Float64,Float64}}}, ::Function, ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/denshe/.julia/packages/OrdinaryDiffEq/tat1c/src/solve.jl:6
 [6] (::getfield(DiffEqBase, Symbol("#kw##__solve")))(::NamedTuple{(:saveat, :abstol, :reltol),Tuple{Array{Int64,1},Float64,Float64}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0 (repeats 5 times)
 [7] #solve#456(::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:saveat, :abstol, :reltol),Tuple{Array{Int64,1},Float64,Float64}}}, ::Function, ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5) at /Users/denshe/.julia/packages/DiffEqBase/6ewdP/src/solve.jl:39
 [8] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:saveat, :abstol, :reltol),Tuple{Array{Int64,1},Float64,Float64}}, ::typeof(solve), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1},getfield(Main, Symbol("##363")){getfield(Main, Symbol("##3#7")),getfield(Main, Symbol("##4#8")),getfield(Main, Symbol("##5#9")),Nothing,Nothing,getfield(Main, Symbol("##6#10")),Expr,Expr},Nothing,DiffEqBase.StandardODEProblem}, ::Tsit5) at ./none:0
 [9] test_f(::Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}, ::Array{Int64,1}) at ./REPL[11]:5
 [10] macro expansion at ./REPL[15]:9 [inlined]
 [11] (::getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}})(::Turing.Core.RandomVariables.UntypedVarInfo, ::Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}, ::Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/core/compiler.jl:388
 [12] #call#3 at /Users/denshe/.julia/packages/Turing/RZOZ8/src/Turing.jl:62 [inlined]
 [13] Model at /Users/denshe/.julia/packages/Turing/RZOZ8/src/Turing.jl:62 [inlined]
 [14] runmodel! at /Users/denshe/.julia/packages/Turing/RZOZ8/src/core/RandomVariables.jl:82 [inlined]
 [15] f at /Users/denshe/.julia/packages/Turing/RZOZ8/src/core/ad.jl:109 [inlined]
 [16] vector_mode_dual_eval(::getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}, ::Array{Real,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}}) at /Users/denshe/.julia/packages/ForwardDiff/N0wMF/src/apiutils.jl:37
 [17] vector_mode_gradient!(::Array{Real,1}, ::getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}, ::Array{Real,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}}) at /Users/denshe/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:103
 [18] gradient! at /Users/denshe/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:35 [inlined]
 [19] gradient!(::Array{Real,1}, ::getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}, ::Array{Real,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4,Array{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#26")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}},Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},Real},Real,4},1}}) at /Users/denshe/.julia/packages/ForwardDiff/N0wMF/src/gradient.jl:33
 [20] gradient_logp_forward(::Array{Real,1}, ::Turing.Core.RandomVariables.UntypedVarInfo, ::Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}}, ::Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/core/ad.jl:116
 [21] gradient_logp at /Users/denshe/.julia/packages/Turing/RZOZ8/src/core/ad.jl:80 [inlined]
 [22] #2 at /Users/denshe/.julia/packages/Turing/RZOZ8/src/inference/support/hmc_core.jl:12 [inlined]
 [23] #_leapfrog#26(::getfield(Turing.Inference, Symbol("##6#7")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}, ::getfield(Turing.Inference, Symbol("##8#9")){Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}, ::Function, ::Array{Real,1}, ::Array{Float64,1}, ::Int64, ::Float64, ::getfield(Turing.Inference, Symbol("##2#3")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}},Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}}}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/inference/support/hmc_core.jl:147
 [24] (::getfield(Turing.Inference, Symbol("#kw##_leapfrog")))(::NamedTuple{(:rev_func, :log_func),Tuple{getfield(Turing.Inference, Symbol("##6#7")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},getfield(Turing.Inference, Symbol("##8#9")){Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}}}, ::typeof(Turing.Inference._leapfrog), ::Array{Real,1}, ::Array{Float64,1}, ::Int64, ::Float64, ::Function) at ./none:0
 [25] #_hmc_step#27(::Function, ::Function, ::Function, ::Array{Real,1}, ::Float64, ::getfield(Turing.Inference, Symbol("##4#5")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}},Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}}}, ::Function, ::getfield(Turing.Inference, Symbol("##16#17")), ::Int64, ::Float64, ::getfield(Turing.Inference, Symbol("##14#15")){Int64}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/inference/support/hmc_core.jl:195
 [26] (::getfield(Turing.Inference, Symbol("#kw##_hmc_step")))(::NamedTuple{(:rev_func, :log_func),Tuple{getfield(Turing.Inference, Symbol("##6#7")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},getfield(Turing.Inference, Symbol("##8#9")){Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}}}, ::typeof(Turing.Inference._hmc_step), ::Array{Real,1}, ::Float64, ::Function, ::Function, ::Function, ::Int64, ::Float64, ::getfield(Turing.Inference, Symbol("##14#15")){Int64}) at ./none:0
 [27] (::getfield(Turing.Inference, Symbol("#kw##hmc_step")))(::NamedTuple{(:rev_func, :log_func),Tuple{getfield(Turing.Inference, Symbol("##6#7")){Turing.Core.RandomVariables.UntypedVarInfo,Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}},getfield(Turing.Inference, Symbol("##8#9")){Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}}}}, ::typeof(Turing.Inference.hmc_step), ::Array{Real,1}, ::Float64, ::Function, ::Function, ::Function, ::Float64, ::HMC{Turing.Core.ForwardDiffAD{40},Any}, ::getfield(Turing.Inference, Symbol("##14#15")){Int64}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/inference/hmc.jl:67
 [28] step(::Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}}, ::Turing.Sampler{HMC{Turing.Core.ForwardDiffAD{40},Any}}, ::Turing.Core.RandomVariables.UntypedVarInfo, ::Val{false}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/inference/hmc.jl:233
 [29] macro expansion at ./util.jl:213 [inlined]
 [30] #sample#37(::Bool, ::Nothing, ::Int64, ::Nothing, ::Function, ::Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}}, ::HMC{Turing.Core.ForwardDiffAD{40},Any}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/inference/hmc.jl:150
 [31] sample(::Turing.Model{Tuple{:alpha,:beta,:gamma,:delta},Tuple{:x},getfield(Main, Symbol("###inner_function#371#18")){Array{Int64,1}},NamedTuple{(:x,),Tuple{Array{Float64,2}}},NamedTuple{(:x,),Tuple{Symbol}}}, ::HMC{Turing.Core.ForwardDiffAD{40},Any}) at /Users/denshe/.julia/packages/Turing/RZOZ8/src/inference/hmc.jl:102
 [32] top-level scope at none:0

I've tried to fill issue on the Turing side (https://github.com/TuringLang/Turing.jl/issues/810) but the error seem to be DifferentialEquations.jl specific. I will appreciate any help. Thank you

Vaibhavdixit02 commented 5 years ago

Hi @denisshepelin you are trying to redo a lot of what this package is meant for actually. BTW we moved away from using IS as the default sampler to NUTS now https://github.com/JuliaDiffEq/DiffEqBayes.jl/blob/master/src/turing_inference.jl#L4, that said you can pass the sampler you want with the sampler arg in turing_inference. Using turing_inference would make your life a bit easier I guess. https://github.com/JuliaDiffEq/DiffEqTutorials.jl/blob/master/tutorials/models/06-pendulum_bayesian_inference.jmd is a good reference on how to use it.

turing_inference(prob, Tsit5(), tp, newsol, [Normal(1.0, 0.2), Normal(0.5, 0.2),Normal(0.5, 0.2), Normal(0.1, 0.1)];syms=[:alpha,:beta, :gamma, :delta],sampler=HMC(1000, 0.1, 5)) this works for me.

You should try ]up once as some of your packages are out of date, that will probably fix the problem.

denisshepelin commented 5 years ago

I've just noted recent update of turing_inference, tried it on my simple example and it worked very nicely!

Vaibhavdixit02 commented 5 years ago

Cool! I'll close this issue then. Do ping in Slack or open an issue if you face any more problems.