How to use DiffEqFlux on DAEProblem

I’m trying to use diffeqflux on differential algebraic equations. I know there is an example in the diffeqflux docs, but there they use constant mass matrices. I want to have the option to use non-constant matrices. Ultimately, I’d like to be able to solve a problem where the mass matrix M contains an unknown function represented by a neural network. I’m unsure if this will be possible so would appreciate insight on the matter. For now I’m just trying to make it work with basic parameter estimation

For sake of simplicity, I started with the first example on parameter estimation and rewrote the LV function to work with DAEProblems,

using DifferentialEquations, DiffEqFlux, Plots, Sundials
function LV!(out,du, u, p, t)
  x, y = u
  α, β, δ, γ = p
  out[1] = α*x - β*x*y - du[1]
  out[2] = -δ*y + γ*x*y - du[2]
end

u0 = [1.0, 1.0]

tspan = (0.0, 10.0)
tsteps = 0.0:0.1:10.0

p = [1.5, 1.0, 3.0, 1.0]

differential_vars = [true,true]
sol = Sundials.solve(DAEProblem(LV!, u0, u0, tspan, p,differential_vars=differential_vars),IDA(),rtol=1e-8)

function loss(p)
    sol = Sundials.solve(DAEProblem(LV!, u0, u0, tspan, p,differential_vars=differential_vars,saveat = tsteps),IDA(),rtol=1e-8)
  loss = sum(abs2, sol.-1)
  return loss, sol
end

callback = function (p, l, pred)
  display(l)
  plt = plot(pred, ylim = (0, 6))
  display(plt)
  return false
end

result_ode = DiffEqFlux.sciml_train(loss, p, cb = callback, maxiters = 100)

The last line gives me the following error:
Continuous adjoint sensitivities are only supported for ODE/SDE/RODE problems.

So I’m unsure if this is just not feasible or if I should be approaching it differently.

I’ve also tried to define a mass matrix and invert explicitly within the function (accepting the fact that I will later have to regularize the mass matrix when it contains the functions represented by NNs). Plugging the following into the diffeqflux parameter estimation example doesn’t work either:

function lotka_volterra!(du, u, p, t)
    x, y = u
    α, β, δ, γ = p
    M = [1.0 0.0;0.0 1.0]
    
    G = zeros(2)
    G[1] = α*x - β*x*y
    G[2] = -δ*y + γ*x*y
    
    G = M\G
    du[1],du[2]=G[1],G[2]
    
end

I get the error
ArgumentError: Converting an instance of ReverseDiff.TrackedReal{Float64, Float64, Nothing} to Float64 is not defined. Please use ReverseDiff.value instead.
I’m not very good at interpreting Julia errors yet so am unsure what the problem is here. I suspect the error is with the \ operator, but unsure how. Either way, my preference would be for using DAEProblem directly, but perhaps that’s not possible.

Thanks!