dynamic_inference update

REQUIRE

@@ -1,5 +1,5 @@
 julia 0.7
-DiffEqBase 1.7.0
+DiffEqBase 
 #Mamba
 Stan
 Distributions
@@ -10,8 +10,9 @@ RecursiveArrayTools
 ParameterizedFunctions
 OrdinaryDiffEq
 Parameters
-#DiffWrappers
-ContinuousTransformations 1.0.0
-DynamicHMC 0.1.1
+DynamicHMC 
 Distances
 ApproxBayes
+TransformVariables
+LogDensityProblems
+Flux

src/DiffEqBayes.jl

@@ -1,9 +1,9 @@
 module DiffEqBayes
 using DiffEqBase, Stan, Distributions, Turing, MacroTools
 using OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
-using DynamicHMC, ContinuousTransformations
+using DynamicHMC, TransformVariables, LogDensityProblems
 using Parameters, Distributions, Optim
-using Distances, ApproxBayes
+using Distances, ApproxBayes, Flux
 
 include("stan_inference.jl")
 include("turing_inference.jl")

src/dynamichmc_inference.jl

@@ -8,9 +8,10 @@ end
 
 function (P::DynamicHMCPosterior)(a)
     @unpack alg, problem, likelihood, priors, kwargs = P
-
-    prob = problem_new_parameters(problem, a)
-    sol = solve(prob, alg; kwargs...)
+    u0_new = Flux.Tracker.TrackedReal.(problem.u0)
+    tspan_new = Flux.Tracker.TrackedReal.(problem.tspan)
+    prob = remake(problem,u0 = u0_new,tspan = tspan_new,p=a.a)
+    sol = solve(prob, alg,dtmin = 1e-4; kwargs...)
     if any((s.retcode != :Success for s in sol))
         ℓ = -Inf
     else
@@ -33,7 +34,6 @@ function dynamichmc_inference(prob::DiffEqBase.DEProblem, alg, t, data, priors,
                               kwargs...)
     likelihood = sol -> sum( sum(logpdf.(Normal(0.0, σ), sol(t) .- data[:, i]))
                              for (i, t) in enumerate(t) )
-
     dynamichmc_inference(prob, alg, likelihood, priors, transformations;
                          ϵ=ϵ, initial=initial, num_samples=num_samples,
                          kwargs...)
@@ -43,39 +43,12 @@ function dynamichmc_inference(prob::DiffEqBase.DEProblem, alg, likelihood, prior
                               ϵ=0.001, initial=Float64[], num_samples=1000,
                               kwargs...)
     P = DynamicHMCPosterior(alg, prob, likelihood, priors, kwargs)
+    PT = TransformedLogDensity(transformations, P)
+    println(typeof(PT))
+    PTG = FluxGradientLogDensity(PT);
 
-    transformations_tuple = Tuple(transformations)
-    parameter_transformation = TransformationTuple(transformations_tuple) # assuming a > 0
-    PT = TransformLogLikelihood(P, parameter_transformation)
-    PTG = ForwardGradientWrapper(PT, zeros(length(priors)));
-
-    lower_bound = Float64[]
-    upper_bound = Float64[]
-
-    for i in priors
-        push!(lower_bound, minimum(i))
-        push!(upper_bound, maximum(i))
-    end
-
-    # If no initial position is given use local minimum near expectation of priors.
-    if length(initial) == 0
-        for i in priors
-            push!(initial, mean(i))
-        end
-        initial_opt = Optim.minimizer(optimize(a -> -P(a),lower_bound,upper_bound,initial,Fminbox(GradientDescent())))
-    end
-
-    initial_inverse_transformed = Float64[]
-    for i in 1:length(initial_opt)
-       para = TransformationTuple(transformations[i])
-       push!(initial_inverse_transformed,inverse(para, (initial_opt[i], ))[1])
-    end
-    #println(initial_inverse_transformed)
-    sample, NUTS_tuned = NUTS_init_tune_mcmc(PTG,
-                                             initial_inverse_transformed,
-                                             num_samples, ϵ=ϵ)
-
-    posterior = ungrouping_map(Vector, get_transformation(PT) ∘ get_position, sample)
+    chain, NUTS_tuned = NUTS_init_tune_mcmc(PTG,num_samples, ϵ=ϵ)
+    posterior = transform.(Ref(PTG.transformation), get_position.(chain));
 
-    return posterior, sample, NUTS_tuned
+    return posterior, chain, NUTS_tuned
 end

src/stan_inference.jl

@@ -93,7 +93,7 @@ function stan_inference(prob::DiffEqBase.DEProblem,t,data,priors = nothing;alg=:
       setup_params = string(setup_params,"row_vector<lower=0>[$length_of_y] sigma$(i-1);")
     end
   end
-  tuple_hyper_params = tuple_hyper_params[1:endof(tuple_hyper_params)-1]
+  tuple_hyper_params = tuple_hyper_params[1:length(tuple_hyper_params)-1]
   differential_equation = generate_differential_equation(f)
   priors_string = string(generate_priors(f,priors))
   stan_likelihood = stan_string(likelihood)

test/dynamicHMC.jl

@@ -1,8 +1,8 @@
 using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions, RecursiveArrayTools
-using DynamicHMC, DiffWrappers, ContinuousTransformations
+using DynamicHMC, TransformVariables
 using Parameters, Distributions, Optim
 
-f1 = @ode_def_nohes LotkaVolterraTest1 begin
+f1 = @ode_def LotkaVolterraTest1 begin
   dx = a*x - x*y
   dy = -3*y + x*y
 end a
@@ -15,9 +15,8 @@ prob1 = ODEProblem(f1,u0,tspan,p)
 t = collect(range(1,stop=10,length=10))   # observation times
 sol = solve(prob1,Tsit5())
 randomized = VectorOfArray([(sol(t[i]) + σ * randn(2)) for i in 1:length(t)])
-data = convert(Array,randomized)
-
-bayesian_result = dynamichmc_inference(prob1, Tsit5(), t, data, [Normal(1.5, 1)], [bridge(ℝ, ℝ⁺, )])
+data = convert(Array,randomized) 
+bayesian_result = dynamichmc_inference(prob1, Tsit5(), t, data, [Normal(1.5, 1)], as((a = asℝ₊,)))
 @test mean(bayesian_result[1][1]) ≈ 1.5 atol=1e-1
 
 # With hand-code likelihood function
@@ -33,11 +32,11 @@ likelihood = function (sol)
     end
     return l
 end
-bayesian_result = dynamichmc_inference(prob1, Tsit5(), likelihood, [Truncated(Normal(1.5, 1), 0, 2)], [bridge(ℝ, ℝ⁺, )])
+bayesian_result = dynamichmc_inference(prob1, Tsit5(), likelihood, [Truncated(Normal(1.5, 1), 0, 2)])
 @test mean(bayesian_result[1][1]) ≈ 1.5 atol=1e-1
 
 
-f1 = @ode_def_nohes LotkaVolterraTest4 begin
+f1 = @ode_def LotkaVolterraTest4 begin
   dx = a*x - b*x*y
   dy = -c*y + d*x*y
 end a b c d

test/stan.jl

@@ -2,7 +2,7 @@ using DiffEqBayes, OrdinaryDiffEq, ParameterizedFunctions,
       RecursiveArrayTools, Distributions, Test
 
 println("One parameter case")
-f1 = @ode_def_nohes LotkaVolterraTest1 begin
+f1 = @ode_def LotkaVolterraTest1 begin
   dx = a*x - x*y
   dy = -3y + x*y
 end a
@@ -25,7 +25,7 @@ theta1 = bayesian_result.chain_results[:,["theta.1"],:]
 
 
 println("Four parameter case")
-f1 = @ode_def_nohes LotkaVolterraTest4 begin
+f1 = @ode_def LotkaVolterraTest4 begin
   dx = a*x - b*x*y
   dy = -c*y + d*x*y
 end a b c d