Problem when trying to asses identifiability of ModelingToolkit and setting observables using `measured_quantities` kwarg.

[TorkelE]

I am trying to write an extension for Catalyst/StructuralIdentifiability that would allow user to input by Catalyst defined ReactionSystems to assess_identifiability and assess_local_identifiability.

Basically, ReactionSystems can be converted to ODESystems, so I would simply use this conversion. However, generally, Catalyst does not have observables directly defined, so I wanted the user to input these using measured_quantities. Typically, they would just be species of the system, so the user would do something like

rs = @reaction_network begin
    (p,d), 0 <--> X
    (k1,k2), X <--> Y
end
assess_identifiability(rs; measured_quantities=[:Y])

and the appropriate conversion would be done under the hood. Also, we could do some automatic removal of conserved quantities to e.g. deal with cases like: https://github.com/SciML/StructuralIdentifiability.jl/issues/189

However, I was trying to prototype things following something like: https://github.com/SciML/StructuralIdentifiability.jl/issues/133 but fails. Here is my minimal example:

# Fectch packages.
using Catalyst, StructuralIdentifiability

# Create model and convert to ODESystem.
rn = @reaction_network begin
    (k1,k2), X <--> Y
end
osys = convert(ODESystem,rn)

# Create equation for observables (X only).
@variables t y(t)
eq = y~states(osys)[1]

# Do identifiability analysis.
assess_identifiability(osys; measured_quantities=[eq])

which gives

[ Info: Preproccessing `ModelingToolkit.AbstractTimeDependentSystem` object
[ Info: Assessing local identifiability
[ Info: Local identifiability assessed in 0.002107838 seconds
[ Info: Assessing global identifiability
[ Info: Functions to check involve states
[ Info: Computing IO-equations
┌ Info: Computed in 0.001850315 seconds
│   :ioeq_time = :ioeq_time
└   ioeq_time = 0.001850315
[ Info: Computing Wronskians
┌ Info: Computed in 0.001188107 seconds
│   :wrnsk_time = :wrnsk_time
└   wrnsk_time = 0.001188107
[ Info: Dimensions of the Wronskians [2]
┌ Info: Ranks of the Wronskians computed in 1.259e-5 seconds
│   :rank_time = :rank_time
└   rank_times = 1.259e-5
[ Info: Simplifying identifiable functions
ERROR: type VSCodeLogger has no field min_level
Stacktrace:
  [1] getproperty
    @ Base ./Base.jl:37 [inlined]
  [2] groebner_loglevel()
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/utils.jl:28
  [3] discover_shape!(state::ParamPunPam.GroebnerState{…}, modular::ParamPunPam.ModularTracker{…}; η::Int64)
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:181
  [4] discover_shape!
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:163 [inlined]
  [5] _paramgb(blackbox::StructuralIdentifiability.IdealMQS{…}, ordering::Groebner.InputOrdering, up_to_degree::Tuple{…}, estimate_degrees::Bool, assess_correctness::Bool, rational_interpolator::Symbol, polynomial_interpolator::Symbol)
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:123
  [6] paramgb(blackbox::StructuralIdentifiability.IdealMQS{…}; kwargs::@Kwargs{…})
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:97
  [7] paramgb
    @ ParamPunPam ~/.julia/packages/ParamPunPam/nTcWq/src/groebner/paramgb.jl:58 [inlined]
  [8] macro expansion
    @ StructuralIdentifiability ./timing.jl:395 [inlined]
  [9] groebner_basis_coeffs(rff::StructuralIdentifiability.RationalFunctionField{…}; seed::Int64, ordering::Groebner.InputOrdering, up_to_degree::Tuple{…}, rational_interpolator::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/RationalFunctionFields/RationalFunctionField.jl:328
 [10] simplified_generating_set(rff::StructuralIdentifiability.RationalFunctionField{…}; p::Float64, seed::Int64, simplify::Symbol, check_variables::Bool, rational_interpolator::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/RationalFunctionFields/RationalFunctionField.jl:529
 [11] macro expansion
    @ StructuralIdentifiability ./timing.jl:395 [inlined]
 [12] initial_identifiable_functions(ode::ODE{…}; p::Float64, known::Vector{…}, with_states::Bool, var_change_policy::Symbol, rational_interpolator::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:139
 [13] initial_identifiable_functions
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:83 [inlined]
 [14] check_identifiability(ode::ODE{…}, funcs_to_check::Vector{…}; known::Vector{…}, p::Float64, var_change_policy::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:181
 [15] check_identifiability
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:163 [inlined]
 [16] assess_global_identifiability(ode::ODE{…}, funcs_to_check::Vector{…}, known::Vector{…}, p::Float64; var_change::Symbol)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:280
 [17] assess_global_identifiability
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/global_identifiability.jl:268 [inlined]
 [18] macro expansion
    @ StructuralIdentifiability ./timing.jl:395 [inlined]
 [19] assess_identifiability(ode::ODE{Nemo.QQMPolyRingElem}; funcs_to_check::Vector{Nemo.QQMPolyRingElem}, p::Float64)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/StructuralIdentifiability.jl:132
 [20] assess_identifiability(ode::ODESystem; measured_quantities::Vector{Equation}, funcs_to_check::Vector{Any}, p::Float64)
    @ StructuralIdentifiability ~/.julia/packages/StructuralIdentifiability/ewu7r/src/StructuralIdentifiability.jl:184
 [21] top-level scope
    @ ~/Desktop/Julia Playground/Catalyst playground/playground.jl:115
Some type information was truncated. Use `show(err)` to see complete types.

[sumiya11]

Thanks for reporting this (and thanks for your work, supporting reaction networks would be very nice!).

I think it's a bug in ParamPunPam.jl. I do current_logger().min_level, which does not work for VSCodeLogger.

Do you know what are more generic ways of getting the logging level? Probably, Logging.min_enabled_level ?

[sumiya11]

@TorkelE I just registered ParamPunPam.jl version 0.2.3, it should fix the problem

[TorkelE]

Thanks :)

Will have a look tomorrow morning.

Also, and sorry for hijacking this issue for an unrelated issue. If I know the value of some parameters, e.g. k1 = 1.2, is there a way to input this information? Had a look but didn't see. (If not, I can make a routine in Catalyst that briefly substitutes this value in before passing the problem on to StructuralIdentifiability)

[sumiya11]

Thank you!

You can try set_parameter_values https://docs.sciml.ai/StructuralIdentifiability/stable/utils/ode/#StructuralIdentifiability.set_parameter_values-Union{Tuple{P},%20Tuple{T},%20Tuple{ODE{P},%20Dict{P,%20T}}}%20where%20{T%3C:AbstractAlgebra.FieldElem,%20P%3C:AbstractAlgebra.MPolyElem{T}}

[TorkelE]

That was what I was looking for, thanks :)

[TorkelE]

It works now, thanks a lot for the quick update :)

(there is some versioning problems preventing me to get the correct latest version on the 1.10 beta, but I think that should be fixed when it is actually released, if not I will update you)