SciML / StructuralIdentifiability.jl

Fast and automatic structural identifiability software for ODE systems
https://docs.sciml.ai/StructuralIdentifiability/stable/
MIT License
109 stars 17 forks source link

Error when attempting identifiability analysis on system with parameter in exponent. #190

Closed TorkelE closed 1 year ago

TorkelE commented 1 year ago

I have a case where I, instead of writing p1 writes 10^p1. The expressions are equivalent, but the rescaling makes sense in some cases (e.g. in the second case the value cannot be negative). However, the second case causes and error, see the following example:

using StructuralIdentifiability, ModelingToolkit

@parameters k1
@variables t LPd(t) ArBr(t) ArR(t)
D = Differential(t)
measured_quantities = [ArR]

@named ode1 = ODESystem([
    D(LPd) ~ -k1*ArBr*LPd,
    D(ArBr) ~ -k1*ArBr*LPd,
    D(ArR) ~ k1*ArBr*LPd
], t)
global_id = assess_identifiability(ode1, measured_quantities=measured_quantities)   # Works fine.
assess_local_identifiability(ode1, measured_quantities=measured_quantities)   # Works fine.

@named ode2 = ODESystem([
    D(LPd) ~ -10^k1*ArBr*LPd,
    D(ArBr) ~ -10^k1*ArBr*LPd,
    D(ArR) ~ 10^k1*ArBr*LPd
], t)
global_id = assess_identifiability(ode2, measured_quantities=measured_quantities)   # Throws an error.
assess_local_identifiability(ode2, measured_quantities=measured_quantities)   # Throws an error.

the error message is:

ERROR: MethodError: no method matching isless(::Int64, ::Nemo.QQMPolyRingElem)

Closest candidates are:
  isless(::Number, ::Unitful.AbstractQuantity)
   @ Unitful ~/.julia/packages/Unitful/orvol/src/quantities.jl:268
  isless(::Integer, ::Nemo.QQFieldElem)
   @ Nemo ~/.julia/packages/Nemo/EuCgH/src/flint/fmpq.jl:367
  isless(::Integer, ::Nemo.ZZRingElem)
   @ Nemo ~/.julia/packages/Nemo/EuCgH/src/flint/fmpz.jl:753
  ...

Stacktrace:
  [1] <(x::Int64, y::Nemo.QQMPolyRingElem)
    @ Base ./operators.jl:343
  [2] <=(x::Int64, y::Nemo.QQMPolyRingElem)
    @ Base ./operators.jl:392
  [3] >=(x::Nemo.QQMPolyRingElem, y::Int64)
    @ Base ./operators.jl:416
  [4] eval_at_nemo(e::SymbolicUtils.BasicSymbolic{Real}, vals::Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/util.jl:348
  [5] (::StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}})(a::SymbolicUtils.BasicSymbolic{Real})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/util.jl:342
  [6] iterate
    @ ./generator.jl:47 [inlined]
  [7] collect_to!(dest::Vector{Int64}, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}}, offs::Int64, st::Int64)
    @ Base ./array.jl:840
  [8] collect_to_with_first!(dest::Vector{Int64}, v1::Int64, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}}, st::Int64)
    @ Base ./array.jl:818
  [9] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:812
 [10] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, StructuralIdentifiability.var"#41#42"{Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem}}})
    @ Base ./array.jl:711
 [11] map(f::Function, A::Vector{Any})
    @ Base ./abstractarray.jl:3261
 [12] eval_at_nemo(e::SymbolicUtils.BasicSymbolic{Real}, vals::Dict{SymbolicUtils.BasicSymbolic{Real}, Nemo.QQMPolyRingElem})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/util.jl:342
 [13] __preprocess_ode(de::ODESystem, measured_quantities::Vector{Tuple{String, SymbolicUtils.BasicSymbolic{Real}}})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/ODE.jl:549
 [14] preprocess_ode(de::ODESystem, measured_quantities::Vector{Num})
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/ODE.jl:465
 [15] assess_local_identifiability(ode::ODESystem; measured_quantities::Vector{Num}, funcs_to_check::Vector{Array}, p::Float64, type::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/rYPpk/src/local_identifiability.jl:194
 [16] top-level scope
    @ ~/Desktop/Research Projects/New-Active-CRN-Learning/playground.jl:135
pogudingleb commented 1 year ago

@TorkelE Yes, right now we cannot handle expressions which are not rational functions (that is fractions of polynomials), see also #144 This indeed could be one of the simple cases for which we could have a preprocessing step, for now the solution would be to write p1 instead of 10^p1

TorkelE commented 1 year ago

Got it, thanks