Using DynamicHMC for a model defined in Python

Specific Domains
1

Hi I’ve been trying to use DynamicHMC to sample my model written in Python. The problem I’m encountering is that I don’t know how to structure the output of my log posterior function in julia correctly, so the MCMC sampler knows what to do with it.

In the examples provided it is suggested to use the following code
t = problem_transformation(p)
P = TransformedLogDensity(t, p)
∇P = ADgradient(:ForwardDiff, P);
https://www.tamaspapp.eu/DynamicHMCExamples.jl/latest/example_linear_regression/

However, because I use python output, the gradient of the likelihood is not calculated automatically anymore. So currently I only compute the log posterior without the gradients. And even if I would calculate these gradients myself. I am not sure how to give this all the proper structure so the sampler mcmc_with_warmup() can use it.

For context the function that computes the logposterior in the example linked above returns output such as (I printed it):
Dual{ForwardDiff.Tag{Base.Fix1{typeof(LogDensityProblems.logdensity), TransformedLogDensity{TransformVariables.TransformTuple{@NamedTuple{β::TransformVariables.ArrayTransformation{TransformVariables.Identity, 1}, σ::TransformVariables.ShiftedExp{true, Int64}}}, LinearRegressionProblem{Vector{Float64}, Matrix{Float64}, Float64}}}, Float64}}(-82.04794174652842,27.668542512145063,-19.965222890254054,-5.894398700546053,10.918972471706919)

This is my julia code:

ENV["JULIA_CONDAPKG_BACKEND"] = "Null" 
import PythonCall
py_module = PythonCall.Core.pyimport("models_douwe") 
using TransformVariables, TransformedLogDensities, LogDensityProblems, LogDensityProblemsAD,
    DynamicHMC, DynamicHMC.Diagnostics, SimpleUnPack, Statistics, Random, Distributions, 
    ThreadTools, MCMCDiagnosticTools, BenchmarkTools

struct Model_Structure
    y::Vector # hydraulic head observations
    model_nr::Int # identifier of the model defined in FloPy
    locations::Tuple # measurement locations
end

function create_model_structure(model_nr::Int, locations::Tuple)
    y_py = py_module.observations(model_nr, locations)
    y_jl = PythonCall.Core.pyconvert(Tuple, y_py)
    y = collect(y_jl)
    p = Model_Structure(y, model_nr, locations)
    return p
end

function (problem::Model_Structure)(θ) 
    @unpack hk, vk = θ 
    @unpack y, model_nr, locations = problem
    #create input suitable for python 
    if typeof(hk) == Float64
        all_k = (hk, vk) 
    else
        all_k = (hk.value, vk.value) 
    end
    #calculate likelihood
    model_output_py = py_module.model_output(model_nr, locations, K = all_k) 
    model_output_jl = PythonCall.pyconvert(Tuple, model_output_py) # julia Tuple
    model_output = collect(model_output_jl) # hydrualic head measurements, as a vector
    ϵ_distribution = Normal(0, 0.1)
    ℓ_error = loglikelihood(ϵ_distribution, y - model_output)
    # calculate prior density
    ℓ_hk = loglikelihood(Normal(10, 10), hk) 
    ℓ_vk = loglikelihood(Normal(10, 10), vk)
    return ℓ_error + ℓ_hk + ℓ_vk
end

function problem_transformation(p::Model_Structure)
    as((hk = asℝ₊, vk = asℝ₊)) # hydraulic conductivity cannot be negative
end

p = create_model_structure(1, LOCATIONS) # redefine for clarity
trans = problem_transformation(p) #transforms all k's to real positive numbers
P = TransformedLogDensity(trans, p) 
∇P = ADgradient(:ForwardDiff, P)

results = mcmc_with_warmup(Random.default_rng(), ∇P, 1000)

Any suggestions on correctly incorporating the gradients in my case?