Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
TorkelE opened this issue Jul 3, 2023 · 2 comments

Comments

@TorkelE
Copy link
Member

TorkelE commented Jul 3, 2023

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
Copy link
Collaborator

@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
Copy link
Member Author

TorkelE commented Jul 4, 2023

Got it, thanks

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants