use prob_thres

Author: TorkelE

docs/src/inverse_problems/structural_identifiability.md

@@ -82,12 +82,12 @@ 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:
+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 `prob_threshold` (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
-assess_identifiability(goodwind_oscillator; measured_quantities=[:M], p=0.999, loglevel=Logging.Error)
+assess_identifiability(goodwind_oscillator; measured_quantities=[:M], prob_threshold=0.999, loglevel=Logging.Error)
 nothing # hide
 ```
-giving a minimum bound of $99.9\%$ chance of correctness. In practise, the bounds used by StructuralIdentifiability are very conservative, which means that while the minimum guaranteed probability of correctness in the default case is $99\%$, in practise it is much higher. While increasing the value of `p` increases the certainty of correctness, it will also increase the time required to assess identifiability.
+giving a minimum bound of $99.9\%$ chance of correctness. In practise, the bounds used by StructuralIdentifiability are very conservative, which means that while the minimum guaranteed probability of correctness in the default case is $99\%$, in practise it is much higher. While increasing the value of `prob_threshold` increases the certainty of correctness, it will also increase the time required to assess identifiability.
 
 ## Local identifiability analysis
 Local identifiability can be assessed through the `assess_local_identifiability` function. While this is already determined by `assess_identifiability`, assessing local identifiability only has the advantage that it is easier to compute. Hence, there might be models where global identifiability analysis fails (or takes a prohibitively long time), where instead `assess_local_identifiability` can be used. This function takes the same inputs as `assess_identifiability` and returns, for each quantity, `true` if it is locally identifiable (or `false` if it is not). Here, for the Goodwind oscillator, we assesses it for local identifiability only:

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

@@ -64,14 +64,14 @@ Notes:
 """
 function SI.assess_local_identifiability(rs::ReactionSystem, args...; measured_quantities = [], 
                                          known_p = [], funcs_to_check = Vector(), remove_conserved = true, 
-                                         ignore_no_measured_warn=false, kwargs...)
+                                         ignore_no_measured_warn=false, prob_threshold = 0.99, kwargs...)
     # Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form.
     osys, conseqs, vars = make_osys(rs; remove_conserved)
     measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn)
     funcs_to_check = make_ftc(funcs_to_check, conseqs, vars)
 
     # Computes identifiability and converts it to a easy to read form.
-    out = SI.assess_local_identifiability(osys, args...; measured_quantities, funcs_to_check, kwargs...)
+    out = SI.assess_local_identifiability(osys, args...; measured_quantities, funcs_to_check, p=prob_threshold, kwargs...)
     return make_output(out, funcs_to_check, reverse.(conseqs))
 end
 
@@ -102,14 +102,14 @@ Notes:
 """
 function SI.assess_identifiability(rs::ReactionSystem, args...; measured_quantities = [], known_p = [], 
                                    funcs_to_check = Vector(), remove_conserved = true, ignore_no_measured_warn=false, 
-                                   kwargs...)
+                                   prob_threshold = 0.99, kwargs...)
     # Creates a ODESystem, list of measured quantities, and functions to check, of SI's preferred form.
     osys, conseqs, vars = make_osys(rs; remove_conserved)
     measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn)
     funcs_to_check = make_ftc(funcs_to_check, conseqs, vars)
 
     # Computes identifiability and converts it to a easy to read form.
-    out = SI.assess_identifiability(osys, args...; measured_quantities, funcs_to_check, kwargs...)
+    out = SI.assess_identifiability(osys, args...; measured_quantities, funcs_to_check, p=prob_threshold, kwargs...)
     return make_output(out, funcs_to_check, reverse.(conseqs))
 end
 
@@ -140,13 +140,13 @@ Notes:
 """
 function SI.find_identifiable_functions(rs::ReactionSystem, args...; measured_quantities = [], 
                                         known_p = [], remove_conserved = true, ignore_no_measured_warn=false, 
-                                        kwargs...)
+                                        prob_threshold = 0.99, kwargs...)
     # Creates a ODESystem, and list of measured quantities, of SI's preferred form.
     osys, conseqs = make_osys(rs; remove_conserved)
     measured_quantities = make_measured_quantities(rs, measured_quantities, known_p, conseqs; ignore_no_measured_warn)
 
     # Computes identifiable functions and converts it to a easy to read form.
-    out = SI.find_identifiable_functions(osys, args...; measured_quantities, kwargs...)
+    out = SI.find_identifiable_functions(osys, args...; measured_quantities, p=prob_threshold, kwargs...)
     return vector_subs(out, reverse.(conseqs))
 end
 

test/extensions/structural_identifiability.jl

@@ -203,6 +203,16 @@ let
     @test make_si_ode(gw_osc_complt; measured_quantities=[gw_osc_complt.M*gw_osc_complt.E]) isa ODE
 end
 
+# Check that `prob_threshold` alternative kwarg works.
+let 
+    rs = @reaction_network begin
+        p, X --> 0
+    end
+
+    assess_identifiability(rs_catalyst; measured_quantities=[rs.X], prob_thres=0.9)
+    assess_identifiability(rs_catalyst; measured_quantities=[rs.X], prob_thres=0.999)
+end
+
 # Tests for hierarchical model with conservation laws at both top and internal levels.
 let
     # Identifiability analysis for Catalyst model.