Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed Nov 14, 2023
1 parent b77a9eb commit ed8eb5c
Show file tree
Hide file tree
Showing 3 changed files with 12 additions and 5 deletions.
6 changes: 3 additions & 3 deletions docs/src/inverse_problems/structural_identifiability.md
Original file line number Diff line number Diff line change
Expand Up @@ -73,11 +73,11 @@ nothing # hide
### Assessing identifiability for specified quantities only
By default, StructuralIdentifiability assesses identifiability for all parameters and variables. It is, however, possible to designate precisely which quantities you want to check using the `funcs_to_check` option. This both includes selecting a smaller subset of parameters and variables to check, or defining customised expressions. Let us consider the Goodwind from previously, and say that we would like to check whether the production parameters ($pₘ$, $pₑ$, and $pₚ$) and the total amount of the three species ($P + M + E$) are identifiable quantities. Here, we would first unpack these (allowing us to form algebraic expressions) and then use the following code:
```example si1
@unpack pₘ, pₑ, and pₚ, M, E, P = goodwind_oscillator
assess_identifiability(goodwind_oscillator; funcs_to_check=[pₘ, pₑ, pₚ, M + E + P], loglevel=Logging.Error)
@unpack pₘ, pₑ, pₚ, M, E, P = goodwind_oscillator
assess_identifiability(goodwind_oscillator; measured_quantities=[:M], funcs_to_check=[pₘ, pₑ, pₚ, M + E + P], loglevel=Logging.Error)
nothing # hide
```


### Probability of correctness
The identifiability methods used can, in theory, produce erroneous results. However, it is possible to adjust the lower bound for the probability of correctness using the argument `p` (by default set to `0.99`, that is, at least a $99%$ chance of correctness). We can e.g. increase the bound through:
```example si2
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ si_ode(rs; measured_quantities = [:X], known_p = [:p])
```
"""
function Catalyst.make_si_ode(rs::ReactionSystem; measured_quantities = [], known_p = [], ignore_no_measured_warn=false)
rs = Catalyst.expand_registered_functions(rs)
known_quantities = make_measured_quantities(rs, measured_quantities, known_p; ignore_no_measured_warn)
return StructuralIdentifiability.preprocess_ode(convert(ODESystem, rs), known_quantities)[1]

Check warning on line 26 in ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

View check run for this annotation

Codecov / codecov/patch

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl#L23-L26

Added lines #L23 - L26 were not covered by tests
end
Expand All @@ -38,20 +39,23 @@ end

# Local identifiability.
function StructuralIdentifiability.assess_local_identifiability(rs::ReactionSystem, args...; measured_quantities = Num[], known_p = Num[], ignore_no_measured_warn=false, kwargs...)
rs = Catalyst.expand_registered_functions(rs)
known_quantities = make_measured_quantities(rs, measured_quantities, known_p; ignore_no_measured_warn)
osys = convert(ODESystem, rs)
return StructuralIdentifiability.assess_local_identifiability(osys, args...; measured_quantities=known_quantities, kwargs...)

Check warning on line 45 in ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

View check run for this annotation

Codecov / codecov/patch

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl#L41-L45

Added lines #L41 - L45 were not covered by tests
end

# Global identifiability.
function StructuralIdentifiability.assess_identifiability(rs::ReactionSystem, args...; measured_quantities = Num[], known_p = Num[], ignore_no_measured_warn=false, kwargs...)
rs = Catalyst.expand_registered_functions(rs)
known_quantities = make_measured_quantities(rs, measured_quantities, known_p; ignore_no_measured_warn)
osys = convert(ODESystem, rs)
return StructuralIdentifiability.assess_identifiability(osys, args...; measured_quantities=known_quantities, kwargs...)

Check warning on line 53 in ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

View check run for this annotation

Codecov / codecov/patch

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl#L49-L53

Added lines #L49 - L53 were not covered by tests
end

# Identifiable functions.
function StructuralIdentifiability.find_identifiable_functions(rs::ReactionSystem, args...; measured_quantities = Num[], known_p = Num[], ignore_no_measured_warn=false, kwargs...)
rs = Catalyst.expand_registered_functions(rs)
known_quantities = make_measured_quantities(rs, measured_quantities, known_p; ignore_no_measured_warn)
osys = convert(ODESystem, rs)
return StructuralIdentifiability.find_identifiable_functions(osys, args...; measured_quantities=known_quantities, kwargs...)

Check warning on line 61 in ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

View check run for this annotation

Codecov / codecov/patch

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl#L57-L61

Added lines #L57 - L61 were not covered by tests
Expand Down
7 changes: 5 additions & 2 deletions test/extensions/structural_identifiability.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,14 @@ using StructuralIdentifiability
### Helper Function ###

# Converts the output dicts from StructuralIdentifiability functions from "weird symbol => stuff" to "symbol => stuff" (the output have some strange meta data which prevents equality checks, this enables this).
# Structural identifiability also provides variables like x (rather than x(t)). This is a bug, but we have to convert to make it work (now just remove any (t) to make them all equal).
function sym_dict(dict_in)
dict_out = Dict{Symbol,Any}()
for key in keys(dict_in)
dict_out[Symbol(key)] = dict_in[key]
end
sym_key = Symbol(key)
sym_key = Symbol(replace(String(sym_key), "(t)" => ""))
dict_out[sym_key] = dict_in[key]
end
return dict_out
end

Expand Down

0 comments on commit ed8eb5c

Please sign in to comment.