DataDrivenDiffEq.jl for estimating DAE system of equations from data

I have been experimenting DataDrivenDiffEq.jl for retrieval of system of differential equations from a dataset( where data is generated by solving that differential equation system). But unlike ODESystem, I find it cumbersome to estimate DAESystem of equations from underlying data.

I have a simple classical case of index 1 DAE - Roberts Equation , whose solution data is used to retrieve corresponding DAESystem, which gives wrong estimates . Here is the code for the same.

using DataDrivenDiffEq
using ModelingToolkit
using OrdinaryDiffEq
using LinearAlgebra , DifferentialEquations

function f(du, u, p, t)
    du[1] = -p[1]*u[1] + p[2]*u[2]*u[3]
    du[2] = p[1]*u[1] - p[2]*u[2]*u[3] - p[3]*u[2]*u[2]
    du[3] = u[1] + u[2] + u[3] - 1.
M = [1 0 0; 0 1 0; 0 0 0.]
p = [0.04, 10^4, 3e7]
u0 = [1.,0.,0.]
tspan = (0., 1e6)
prob = ODEProblem(ODEFunction(f, mass_matrix = M), u0, tspan, p)
sol1 = solve(prob,saveat=10.0, Rodas5())

ddprob = ContinuousDataDrivenProblem(sol1)
@variables t x(t) y(t) 
u = [x;y]

basis = Basis(polynomial_basis(u,3), u, iv = t)
opt = STLSQ()
ddsol = solve(ddprob,DMDSVD(), normalize = true)
print(ddsol, Val{true})

Is there a way to retreive original DAESystem providing appropriate basis in this case?
Thanking in advanceā€¦

Did you try the implicit SINDy? The mass matrix form will be much more stable.