Identifiability from `ODEProblem` defined via function

[iliailmer]

The following code

using StructuralIdentifiability, OrdinaryDiffEq, ModelingToolkit
function f(du,u,p,t)
  b,c,a,beta,g,delta,sigma = p
  du[1] = -b * u[1] + 1/(c + u[4])
  du[2] = a * u[1] - beta * u[2]
  du[3] = g * u[2] - delta * u[3]
  du[4] = sigma * u[4] * (g * u[2] - delta * u[3])/u[3]
end
prob = ODEProblem(f,ones(4),(0.0,1.0),ones(7))
sys = modelingtoolkitize(prob)
st = states(sys)
params = parameters(sys)
@variables t y(t) x[1:length(st)](t) 
@parameters mu[1:length(params)]

states_map = Dict(s=>x[i] for (i, s) ∈ enumerate(st))
param_map = Dict(p=>mu[i] for (i, p) ∈ enumerate(params))

eqs = substitute(substitute(equations(sys), states_map), param_map)

sys = ODESystem(eqs, t, name=:Name)
measured_quantities = [y~states_map[us[1]]]
@time global_id = assess_identifiability(sys, measured_quantities=measured_quantities)

Produces this error at the very end

ERROR: ArgumentError: The polynomial contains a variable mu[2] not present in the new ring 9733449230697594235070667*mu[2]^2*mu[1]^5*mu[4]*mu[7]^2 + 9733449230697594235070667*mu[2]^2*mu[1]^5*mu[6]*mu[7]^2 - 48667246153487971175353335*mu[2]^2*mu[1]^4*mu[4]*mu[6]*mu[7]^3 - 12977932307596792313427556*mu[2]^2*mu[1]^4*mu[4]*mu[6]*mu[7]^2 + 520606133531898376008304868289618073659180126368695604752803844382800307368205982265758604268378386364287597134108550309252380302258675*mu[7]^2

This is due to a Singular issue, as in #43. I could not find a workaround so far. If it's possible to automate this part

@variables t y(t) x[1:length(st)](t) 
@parameters mu[1:length(params)]

avoiding indexed notation, i.e. x[1].

[pogudingleb]

@iliailmer Thank you for reporting this! Just to clarify: is the problem with the symbols [ and ] in the variable name? Or I missing something?

[pogudingleb]

So the problem still persists although the dependence on Singular has been eliminated. Apparently, the variables of the form x[i] are not converted properly. This would be great to fix.

[iliailmer]

@pogudingleb I just checked this code example and it works 😄

using StructuralIdentifiability, ModelingToolkit
function f(du, u, p, t)
  b, c, a, beta, g, delta, sigma = p
  du[1] = -b * u[1] + 1 / (c + u[4])
  du[2] = a * u[1] - beta * u[2]
  du[3] = g * u[2] - delta * u[3]
  du[4] = sigma * u[4] * (g * u[2] - delta * u[3]) / u[3]
end
prob = ODEProblem(f, ones(4), (0.0, 1.0), ones(7))
sys = modelingtoolkitize(prob)
st = states(sys)
params = parameters(sys)

measured_quantities = [y ~ st[1]]
@time global_id = assess_identifiability(sys, measured_quantities=measured_quantities)

The output:

Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, Symbol} with 7 entries:
  α₃ => :nonidentifiable
  α₇ => :globally
  α₆ => :locally
  α₄ => :locally
  α₁ => :globally
  α₂ => :globally
  α₅ => :nonidentifiable

[pogudingleb]

Thanks, @iliailmer ! So we are done here