Dear all, I am to fit a model that uses an interaction surface in Turing.jl, the interaction surface is a function that is outside the model but gets called inside the model to estimate the competitive equivalence of individuals based on their size and density. It has a series of parameters that define that interaction, but for some reason, I am not able to make it work in the model. I was able to do it in stan, but I would really like to be able to do it in Turing. Below there is a minimal reproducible example:

```
# for the sake of time and learning I am performing the IPM analyses for the mesocosms in Julia
# load the modules
using Chain
using DataFrames
using Plots
using Random
using StatsPlots
using Statistics
using LinearAlgebra
using Turing
using DynamicHMC
using Distributions
using MCMCChains
function α_surface!(target_tank::Int64, tank_var::Vector{Int64},
spp::Vector{Int64}, z::Vector{Float64},
ϕ_i, ϕ_j, η_ij, ρ_i, ρ_j, ι_ij, H, ω )
sz = z[findall(tank_var .== target_tank)]# get the body size values per tank or channel
sp = spp[findall(tank_var .== target_tank)] # get the species for each tank
complength = zeros(length(sz),length(sz))
for j in eachindex(sz)
for k in eachindex(sz)
if(sp[j] == 1 && sp[k] == 1)
complength[j,k] = exp( ϕ_i * (sz[k] - ω) - ϕ_i * (sz[j] - ω)) * exp(-((ρ_i * (sz[k] - ω) - ρ_i * (sz[j] - ω))^2)/(2 * H^2))
end
if(sp[j] == 1 && sp[k] == 0)
complength[j,k] = exp(η_ij + ϕ_j * (sz[k] - ω) - ϕ_i * (sz[j] - ω)) * exp(-((ι_ij + ρ_j * (sz[k] - ω) - ρ_i * (sz[j] - ω))^2/(2 * H^2)))
end
if(sp[j] == 0 && sp[k] == 0)
complength[j,k] = (exp( ϕ_j * (sz[k] - ω) - ϕ_j * (sz[j] - ω)) )* exp(-((ρ_j * (sz[k] - ω) - ρ_j * (sz[j] - ω))^2)/(2 * H^2))
end
if(sp[j] == 0 && sp[k] == 1)
complength[j,k] = exp(-η_ij + ϕ_j * (sz[k] - ω) - ϕ_i .* (sz[j] - ω)) * exp(-((-1*ι_ij + ρ_j * (sz[k] - ω) - ρ_i * (sz[j] - ω))^2)/(2 * H^2))
end
end
end
return TArray(sum(complength, dims= 2)[:,1])
end
# # function to get the density for all tanks or channels
function α_density(tank_var::Vector{Int64}, spp::Vector{Int64}, z::Vector{Float64},
ϕ_i, ϕ_j, η_ij, ρ_i, ρ_j, ι_ij, H, ω )
d1 = []
for c in 1:maximum(tank_var)
d2 = α_surface!(c, tank_var, spp, z, ϕ_i, ϕ_j, η_ij, ρ_i, ρ_j, ι_ij, H, ω )
d1 = append!(d1, d2)
end
return d1
end
df = DataFrame()
for i in 1:50
append!(df, DataFrame(size = rand(Normal(32, 10), rand(5:10)), spp= 1, tank = i ))
append!(df, DataFrame(size = rand(Normal(40, 10), rand(5:10)), spp= 0, tank = i ))
end
# simulate data
α_density(df.tank, df.spp, df.size, 0.05,0,0,0,0,0,1,10)
N = α_density(df.tank, df.spp, df.size, 0.02,0.02,0,0,0,0,1,10)
df.growth .= 0.397 .+ 0.2*df.spp .- 0.762 .* log.(df.size ./ 10) .+
0.355.* log.(df.size.^2 ./ 10^2) .- (0.073) .* N .+ 0.084 .*N
df.growth .= exp.(rand.(Normal.(df.growth, 0.04)))
scatter(df.size, df.growth, group = df.spp)
α_density(df.tank, df.spp, df.size, 0.05,0,0,0,0,0,1,10)
@model function Model_Inter(growth, spp::Vector{Int64}, int_size::Vector{Float64}, ChanID, center)
size = int_size
z = (log.(int_size ./ center))
z2 = (log.(int_size.^2 ./ center^2))
## Total num of y
N = length(growth)
## Separate σ priors for each actor
σ_ch ~ truncated(Cauchy(0, 1), 0, 10)
## Number of unique mesocosms in the data set
N_ch = length(unique(ChanID)) #7
# --------------Survival model
## Vector of mesocosms (1,..,7) which we'll set priors on
α_ch ~ filldist(Normal(0, σ_ch), N_ch)
#λ ~ TArray( )#filldist(Normal(), N) # a latent variable
## Prior for intercept, prosoc_left, and the interaction
α ~ Normal(0, 10)
βz ~ Normal(0, 10)
βz2 ~ Normal(0, 10)
βspp ~ Normal(0, 10)
βN ~ Normal(0, 10)
βzN ~ Normal(0, 10)
ϕ_i ~ Normal(0.495, 0.109)
# ϕ_j ~ Normal(0.867, 0.102)
# η_ij ~ Normal(-0.020, 0.010)
# ρ_i ~ Normal(0.10, 0.021)
# ρ_j ~ Normal( 0.03, 0.004)
# ι_ij ~ Normal((3.71-3.12), 0.165)
σ ~ truncated(Cauchy(0,1),0,10)
H= 1.
ω= center
Den = tzeros(N)
Den .= α_density(ChanID, spp, size, ϕ_i,0,0,0,0,0,H,ω)
μ = α .+ α_ch[ChanID] .+ βspp .* spp .+ βz .* z .+ βz2 .* z2 .+
βN .* Den .+
βzN .* Den
growth .~ MvNormal(μ, σ)
end;
chns = sample(
Model_Inter(df.growth,df.spp, df.size, df.tank, 10.0 ),
NUTS(),
100
)
```

Thanks for the help