tpapp / DynamicHMC.jl

Implementation of robust dynamic Hamiltonian Monte Carlo methods (NUTS) in Julia.
Other
245 stars 21 forks source link

Epsilon adaptation is machine-dependent #38

Closed ChrisRackauckas closed 5 years ago

ChrisRackauckas commented 5 years ago

I'm having a hard time getting an MWE because right now it's inside of a private repo and the result is different on every run, but the setup is like this. It's a Bayesian estimation of some ODE system, with Random.seed!(1) before the call. However, Travis, my laptop, and my workstation all give different results. Somehow, Travis passes tests (yay!), my laptop runs but is slightly off on one parameter, while my workstation fails with the error pasted below. It seems to be tied to different steps being taken, and I'm not sure entirely what would cause it.

Bayesian models: Error During Test at /home/travis/.julia/packages/SafeTestsets/A83XK/src/SafeTestsets.jl:25
  Got exception outside of a @test
  LoadError: ArgumentError: 0 ≤ a must hold. Got
  0 => 0.0
  a => NaN
  Stacktrace:
   [1] macro expansion at /home/travis/.julia/packages/ArgCheck/BUMkA/src/checks.jl:165 [inlined]
   [2] adapt_stepsize(::DynamicHMC.DualAveragingParameters{Float64}, ::DynamicHMC.DualAveragingAdaptation{Float64}, ::Float64) at /home/travis/.julia/packages/DynamicHMC/pkv4F/src/stepsize.jl:162
   [3] tune(::DynamicHMC.NUTS{Array{Float64,1},Float64,Random.MersenneTwister,DynamicHMC.Hamiltonian{PuMaS.BayesLogDensity{PKPDModel{ParamSet{NamedTuple{(:θ, :Ω, :σ),Tuple{Constrained{MvNormal{Float64,PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{Float64,1}}},InverseWishart{Float64,PDMat{Float64,Array{Float64,2}}},RealDomain{Float64}}}},getfield(Main.##408, Symbol("##27#40")),getfield(Main.##408, Symbol("##28#41")),getfield(Main.##408, Symbol("##29#42")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(ModelingToolkit, Symbol("##156#157")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main.##408, Symbol("##30#43")),getfield(Main.##408, Symbol("##35#48"))},Population{Subject{StructArrays.StructArray{NamedTuple{(:dv,),Tuple{Float64}},1,NamedTuple{(:dv,),Tuple{Array{Float64,1}}}},NamedTuple{(:time, :WT, :SEX),Tuple{Array{Float64,1},Float64,Int64}},Array{PuMaS.Event{Float64,Float64,Float64,Float64,Float64,Float64},1},Float64}},Array{Float64,1},ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8},1}},DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}},Tuple{},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},DynamicHMC.GaussianKE{Symmetric{Float64,Array{Float64,2}},LowerTriangular{Float64,Array{Float64,2}}}},DynamicHMC.ReportIO{Base.TTY}}, ::DynamicHMC.StepsizeTuner) at /home/travis/.julia/packages/DynamicHMC/pkv4F/src/sampler.jl:74
   [4] tune(::DynamicHMC.NUTS{Array{Float64,1},Float64,Random.MersenneTwister,DynamicHMC.Hamiltonian{PuMaS.BayesLogDensity{PKPDModel{ParamSet{NamedTuple{(:θ, :Ω, :σ),Tuple{Constrained{MvNormal{Float64,PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{Float64,1}}},InverseWishart{Float64,PDMat{Float64,Array{Float64,2}}},RealDomain{Float64}}}},getfield(Main.##408, Symbol("##27#40")),getfield(Main.##408, Symbol("##28#41")),getfield(Main.##408, Symbol("##29#42")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(ModelingToolkit, Symbol("##156#157")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main.##408, Symbol("##30#43")),getfield(Main.##408, Symbol("##35#48"))},Population{Subject{StructArrays.StructArray{NamedTuple{(:dv,),Tuple{Float64}},1,NamedTuple{(:dv,),Tuple{Array{Float64,1}}}},NamedTuple{(:time, :WT, :SEX),Tuple{Array{Float64,1},Float64,Int64}},Array{PuMaS.Event{Float64,Float64,Float64,Float64,Float64,Float64},1},Float64}},Array{Float64,1},ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8},1}},DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}},Tuple{},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}},DynamicHMC.GaussianKE{Diagonal{Float64,Array{Float64,1}},LowerTriangular{Float64,Diagonal{Float64,Array{Float64,1}}}}},DynamicHMC.ReportIO{Base.TTY}}, ::DynamicHMC.TunerSequence{Tuple{DynamicHMC.StepsizeTuner,DynamicHMC.StepsizeCovTuner{Float64},DynamicHMC.StepsizeCovTuner{Float64},DynamicHMC.StepsizeCovTuner{Float64},DynamicHMC.StepsizeCovTuner{Float64},DynamicHMC.StepsizeCovTuner{Float64},DynamicHMC.StepsizeTuner}}) at /home/travis/.julia/packages/DynamicHMC/pkv4F/src/sampler.jl:262
   [5] #NUTS_init_tune_mcmc#20(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Random.MersenneTwister, ::PuMaS.BayesLogDensity{PKPDModel{ParamSet{NamedTuple{(:θ, :Ω, :σ),Tuple{Constrained{MvNormal{Float64,PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{Float64,1}}},InverseWishart{Float64,PDMat{Float64,Array{Float64,2}}},RealDomain{Float64}}}},getfield(Main.##408, Symbol("##27#40")),getfield(Main.##408, Symbol("##28#41")),getfield(Main.##408, Symbol("##29#42")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(ModelingToolkit, Symbol("##156#157")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main.##408, Symbol("##30#43")),getfield(Main.##408, Symbol("##35#48"))},Population{Subject{StructArrays.StructArray{NamedTuple{(:dv,),Tuple{Float64}},1,NamedTuple{(:dv,),Tuple{Array{Float64,1}}}},NamedTuple{(:time, :WT, :SEX),Tuple{Array{Float64,1},Float64,Int64}},Array{PuMaS.Event{Float64,Float64,Float64,Float64,Float64,Float64},1},Float64}},Array{Float64,1},ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8},1}},DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}},Tuple{},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::Int64) at /home/travis/.julia/packages/DynamicHMC/pkv4F/src/sampler.jl:285
   [6] NUTS_init_tune_mcmc at /home/travis/.julia/packages/DynamicHMC/pkv4F/src/sampler.jl:284 [inlined]
   [7] #NUTS_init_tune_mcmc#21 at /home/travis/.julia/packages/DynamicHMC/pkv4F/src/sampler.jl:294 [inlined]
   [8] NUTS_init_tune_mcmc(::PuMaS.BayesLogDensity{PKPDModel{ParamSet{NamedTuple{(:θ, :Ω, :σ),Tuple{Constrained{MvNormal{Float64,PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{Float64,1}}},InverseWishart{Float64,PDMat{Float64,Array{Float64,2}}},RealDomain{Float64}}}},getfield(Main.##408, Symbol("##27#40")),getfield(Main.##408, Symbol("##28#41")),getfield(Main.##408, Symbol("##29#42")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(ModelingToolkit, Symbol("##156#157")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main.##408, Symbol("##30#43")),getfield(Main.##408, Symbol("##35#48"))},Population{Subject{StructArrays.StructArray{NamedTuple{(:dv,),Tuple{Float64}},1,NamedTuple{(:dv,),Tuple{Array{Float64,1}}}},NamedTuple{(:time, :WT, :SEX),Tuple{Array{Float64,1},Float64,Int64}},Array{PuMaS.Event{Float64,Float64,Float64,Float64,Float64,Float64},1},Float64}},Array{Float64,1},ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(LogDensityProblems.logdensity),Float64},Float64,8},1}},DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}},Tuple{},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}}, ::Int64) at /home/travis/.julia/packages/DynamicHMC/pkv4F/src/sampler.jl:294
   [9] #fit#207(::Int64, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::PKPDModel{ParamSet{NamedTuple{(:θ, :Ω, :σ),Tuple{Constrained{MvNormal{Float64,PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{Float64,1}}},InverseWishart{Float64,PDMat{Float64,Array{Float64,2}}},RealDomain{Float64}}}},getfield(Main.##408, Symbol("##27#40")),getfield(Main.##408, Symbol("##28#41")),getfield(Main.##408, Symbol("##29#42")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(ModelingToolkit, Symbol("##156#157")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main.##408, Symbol("##30#43")),getfield(Main.##408, Symbol("##35#48"))}, ::Population{Subject{StructArrays.StructArray{NamedTuple{(:dv,),Tuple{Float64}},1,NamedTuple{(:dv,),Tuple{Array{Float64,1}}}},NamedTuple{(:time, :WT, :SEX),Tuple{Array{Float64,1},Float64,Int64}},Array{PuMaS.Event{Float64,Float64,Float64,Float64,Float64,Float64},1},Float64}}, ::PuMaS.BayesMCMC) at /home/travis/build/UMCTM/PuMaS.jl/src/models/bayes.jl:133
   [10] (::getfield(StatsBase, Symbol("#kw##fit")))(::NamedTuple{(:nsamples,),Tuple{Int64}}, ::typeof(fit), ::PKPDModel{ParamSet{NamedTuple{(:θ, :Ω, :σ),Tuple{Constrained{MvNormal{Float64,PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{Float64,1}}},InverseWishart{Float64,PDMat{Float64,Array{Float64,2}}},RealDomain{Float64}}}},getfield(Main.##408, Symbol("##27#40")),getfield(Main.##408, Symbol("##28#41")),getfield(Main.##408, Symbol("##29#42")),ODEProblem{Nothing,Tuple{Nothing,Nothing},false,Nothing,ODEFunction{false,getfield(ModelingToolkit, Symbol("##156#157")),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem},getfield(Main.##408, Symbol("##30#43")),getfield(Main.##408, Symbol("##35#48"))}, ::Population{Subject{StructArrays.StructArray{NamedTuple{(:dv,),Tuple{Float64}},1,NamedTuple{(:dv,),Tuple{Array{Float64,1}}}},NamedTuple{(:time, :WT, :SEX),Tuple{Array{Float64,1},Float64,Int64}},Array{PuMaS.Event{Float64,Float64,Float64,Float64,Float64,Float64},1},Float64}}, ::PuMaS.BayesMCMC) at ./none:0
   [11] top-level scope at none:0
   [12] include at ./boot.jl:326 [inlined]
   [13] include_relative(::Module, ::String) at ./loading.jl:1038
   [14] include at ./sysimg.jl:29 [inlined]
   [15] include(::String) at /home/travis/.julia/packages/SafeTestsets/A83XK/src/SafeTestsets.jl:23
   [16] top-level scope at /home/travis/build/UMCTM/PuMaS.jl/test/nlme/runtests.jl:16
   [17] top-level scope at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.1/Test/src/Test.jl:1083
   [18] top-level scope at /home/travis/build/UMCTM/PuMaS.jl/test/nlme/runtests.jl:16
   [19] eval(::Module, ::Any) at ./boot.jl:328
   [20] top-level scope at /home/travis/.julia/packages/SafeTestsets/A83XK/src/SafeTestsets.jl:23
   [21] top-level scope at util.jl:156
   [22] top-level scope at /home/travis/build/UMCTM/PuMaS.jl/test/nlme/runtests.jl:16
   [23] include at ./boot.jl:326 [inlined]
   [24] include_relative(::Module, ::String) at ./loading.jl:1038
   [25] include(::Module, ::String) at ./sysimg.jl:29
   [26] include(::String) at ./client.jl:403
   [27] top-level scope at /home/travis/build/UMCTM/PuMaS.jl/test/runtests.jl:54
   [28] top-level scope at util.jl:156
   [29] include at ./boot.jl:326 [inlined]
   [30] include_relative(::Module, ::String) at ./loading.jl:1038
   [31] include(::Module, ::String) at ./sysimg.jl:29
   [32] include(::String) at ./client.jl:403
   [33] top-level scope at none:0
   [34] eval(::Module, ::Any) at ./boot.jl:328
   [35] exec_options(::Base.JLOptions) at ./client.jl:243
   [36] _start() at ./client.jl:436
  in expression starting at /home/travis/build/UMCTM/PuMaS.jl/test/nlme/bayes.jl:108
tpapp commented 5 years ago

This looks similar to #37, are you using the latest tagged version?

If yes, and the error still persists, I will try to help you narrow it down.

ChrisRackauckas commented 5 years ago

Yup, seeing it on v1.0.3 and it's the same project @simonbyrne was working on. It seems to have fixed it in some cases, but not others. Curiously, I haven't seen it on Linux but do see it on Windows (both 64-bit), but that correlation might be a fluke since it's different computers as well.

tpapp commented 5 years ago

I have a hunch, can you please try ~#38~ #39?

tpapp commented 5 years ago

I would love to follow up on this, but I need some code I can run.

ChrisRackauckas commented 5 years ago

https://github.com/JuliaDiffEq/DiffEqBayes.jl/blob/master/test/dynamicHMC.jl is pretty much an MWE.

tpapp commented 5 years ago

thanks, that's on my agenda but I thought that was a different issue

tpapp commented 5 years ago

I think this is fixed by fixing that MWE according to

https://github.com/JuliaDiffEq/DiffEqBayes.jl/issues/63#issuecomment-521211764

Will keep this open to remind me to make that PR.

tpapp commented 5 years ago

I think this is closed (commenting in the issue above).