BoundsError when inverting a model with Turing

I am trying to infer the parameters (a11, a12, …, c21, c22) of the following model:

@model function invert(u, y, pE, pC, options)
    
    # Specify priors
    h = [Normal(pE.H[i], pC.H[i]) for i in eachindex(pE.H)]
    hc = [Normal(pE.HC[i], pC.HC[i]) for i in eachindex(pE.HC)]
    H = reshape(rand.(h), 2, 7)
    HC = reshape(rand.(hc), 2, 2)
    
    a11 ~ Normal(pE.A[1,1], pC.A[1,1])
    a12 ~ Normal(pE.A[1,2], pC.A[1,2])
    a21 ~ Normal(pE.A[2,1], pC.A[2,1])
    a22 ~ Normal(pE.A[2,2], pC.A[2,2])
    
    c11 ~ Normal(pE.C[1,1], pC.C[1,1])
    c12 ~ Normal(pE.C[1,2], pC.C[1,2])
    c21 ~ Normal(pE.C[2,1], pC.C[2,1])
    c22 ~ Normal(pE.C[2,2], pC.C[2,2])
    
    x0 = zeros(options.ℓ, 1) # initial condition
    z0 = zeros(options.ℓ, 6) # initial condition
    states = [x0, z0]
    P = [a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC]
    tspan = (0.0, options.duration)
    predicted = forward_model(P, states, options)  
 
    y[:,1] ~ Normal(predicted[:,1]) 
    y[:,2] ~ Normal(predicted[:,2])
 
end

m = invert(U, y, pE, pC, options)
chain = mapreduce(c -> sample(m, NUTS(.65),1000), chainscat, 1:3) 

but get the error message:

BoundsError: attempt to access 1-element Vector{Matrix{Float64}} at index [2]

The error occurs in a function called g, which is called by the function forward_model:

function forward_model(params, states, options)
    a11, a12, a21, a22, c11, c12, c21, c22, u, H, HC = params
    x0, z0 = states
    
    ts = range(0.05, step=options.dt, stop=options.duration) # sample times
    tspan = (0.0, options.duration) # time span
    
    p = [a11, a12, a21, a22, c11, c12, c21, c22, u, options.dt] # parameters
    ode = ODEProblem(f, x0, tspan, p)
    sol = solve(ode, Tsit5(), saveat=ts, tstops=ts, dt=1e-2)
    
    ph = [H, HC, sol.u, options] # parameters
    odeh = ODEProblem(h, z0, tspan, ph)
    solh = solve(odeh, Tsit5(), saveat=ts, tstops=ts, dt = 1e-2)
    
    # Model signal
    y = g(solh.u, H, options)
end

function g(z, H, opt)
    # Define constants
    TE = 0.04
    V₀ = 100*0.04
    r₀  = 110
    ϑ₀ = 84.8
    E₀ = opt.Pf_H[:,5] .* exp.(H[:,5]) 
    εₕ = opt.Pf_H[2,6] .* exp.(H[2,6])
    
    # Compute coefficients
    k₁ = 4.3.*ϑ₀.*E₀.*TE
    k₂ = εₕ.*r₀.*E₀.*TE
    k₃ = 1 - εₕ
    
    # Output equation
    v = map(i -> (z[i][:,3]), 1:12000)
    v = transpose(reduce(hcat, v))
    q = map(i -> (z[i][:,4]), 1:12000)
    q = transpose(reduce(hcat, q))
    y = V₀ * (k₁'.*(1 .- q) + k₂'.*(1 .- q./v) + k₃'.*(1 .- v))
end

For the sake of simplicity I have left out all other functions (f, h, etc.); I think that these should not be necessary for answering my question. The input z into function g is the solution of an ODE problem and has the following form:

12000-element Vector{Matrix{Float64}}:
 [0.0 1.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0]
 [0.0 1.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0]
 ⋮

Generating data with the function forward_model works. The BoundsError during the inference procedure occurs in the line:

v = map(i -> (z[i][:,3]), 1:12000)

(in function g). I can see that the error occurs because the input argument y (in function invert) is a Matrix{Float64} of size (12000, 2), while z (in function g) is a Vector{Matrix{Float64}} of size (12000,), but my attempts to resolve the problem remain unsuccessful.

I have tried, for example, to convert the data in the function invert:

d = [c[:] for c in eachrow(y)]
d[:,1] ~ Normal(predicted[:,1]) 
d[:,2] ~ Normal(predicted[:,2])

but this converts d to a Vector{Vector{Float64}} and leads to the same BoundsError. I also tried

d = [c for c in eachrow(y)]

and

d = [(y)]

but these are also not correct. How can I convert y to a Vector{Matrix{Float64}}? Or is there a better alternative solution?

Ok, the following converts the data to a Vector{Matrix{Float64}} of size (12000,):

d = [hcat(c[:]) for c in eachrow(data)]
d[:,1] ~ Normal(predicted[:,1]) 
d[:,2] ~ Normal(predicted[:,2])

but this was not the solution after all - I still get the same BoundsError in the same place.

Any suggestions of what else may be wrong with my code/approach would be very welcome.