diff --git a/docs/src/tutorials/discrete_time.md b/docs/src/tutorials/discrete_time.md index 6978fb4a7..212f59db8 100644 --- a/docs/src/tutorials/discrete_time.md +++ b/docs/src/tutorials/discrete_time.md @@ -67,8 +67,8 @@ The `assess_local_identifiability` function has three important keyword argument assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [β * S]) ``` - - `p` is the probability of correctness (default value `0.99`, i.e., 99%). The underlying algorithm is a Monte-Carlo algorithm, so in - principle it may produce incorrect result but the probability of correctness of the returned result is guaranteed to be at least `p` + - `prob_threshold` is the probability of correctness (default value `0.99`, i.e., 99%). The underlying algorithm is a Monte-Carlo algorithm, so in + principle it may produce incorrect result but the probability of correctness of the returned result is guaranteed to be at least `prob_threshold` (in fact, the employed bounds are quite conservative, so in practice incorrect result is almost never produced). - `known_ic` is a list of the states for which initial conditions are known. In this case, the identifiability results will be valid not diff --git a/docs/src/tutorials/identifiability.md b/docs/src/tutorials/identifiability.md index d89ceb5f4..c88ba3d94 100644 --- a/docs/src/tutorials/identifiability.md +++ b/docs/src/tutorials/identifiability.md @@ -59,8 +59,8 @@ Function `assess_local_identifiability` has several optional parameters - `funcs_to_check` a list of specific functions of parameters and states to check identifiability for (see an example below). If not provided, the identifiability is assessed for all parameters and states. - - `p` (default $0.99$) is the probability of correctness. The algorithm can, in theory, produce wrong result, but the probability that it is correct - is guaranteed to be at least `p`. However, the probability bounds we use are quite conservative, so the actual probability of correctness is + - `prob_threshold` (default $0.99$, i.e. 99%) is the probability of correctness. The algorithm can, in theory, produce wrong result, but the probability that it is correct + is guaranteed to be at least `prob_threshold`. However, the probability bounds we use are quite conservative, so the actual probability of correctness is likely to be much higher. - `type` (default `:SE`). By default, the algorithm checks the standard single-experiment identifiability. If one sets `type = :ME`, then the algorithm checks multi-experiment identifiability, that is, identifiability from several experiments with independent initial conditions (the algorithm from [^2] is used). @@ -105,7 +105,7 @@ Similarly to `assess_local_identifiability`, this function has optional paramete more involved than for the parameters, so one may want to call the function with `funcs_to_check = ode.parameters` if the call `assess_identifiability(ode)` takes too long. - - `p` (default $0.99$) is the probability of correctness. Same story as above: the probability estimates are very conservative, so the actual + - `prob_threshold` (default $0.99$, i.e. 99%) is the probability of correctness. Same story as above: the probability estimates are very conservative, so the actual error probability is much lower than 1%. Also, currently, the probability of correctness does not include the probability of correctness of the modular reconstruction for Groebner bases. This probability is ensured by an additional check modulo a large prime, and can be neglected for practical purposes. diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index 9c8e68fbf..cf38a9da4 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -76,17 +76,17 @@ end # ------------------------------------------------------------------------------ """ - field_contains(field, ratfuncs, p) + field_contains(field, ratfuncs, prob_threshold) Checks whether given rational function field `field` contains given rational functions `ratfuncs` (represented as a list of lists). The result is correct with -probability at least `p` +probability at least `prob_threshold` Inputs: - `field` - a rational function field - `ratfuncs` - a list of lists of polynomials. Each of the lists, say, `[f1, ..., fn]`, defines generators `f2/f1, ..., fn/f1`. -- `p` real number from (0, 1) +- `prob_threshold` real number from (0, 1) Output: - a list `L[i]` of bools of length `length(rat_funcs)` such that `L[i]` is true iff @@ -95,7 +95,7 @@ Output: @timeit _to function field_contains( field::RationalFunctionField{T}, ratfuncs::Vector{Vector{T}}, - p, + prob_threshold, ) where {T} if isempty(ratfuncs) return Bool[] @@ -133,7 +133,7 @@ Output: 3 * BigInt(degree)^(length(total_vars) + 3) * (length(ratfuncs)) * - ceil(1 / (1 - p)), + ceil(1 / (1 - prob_threshold)), ) @debug "\tSampling from $(-sampling_bound) to $(sampling_bound)" @@ -153,24 +153,36 @@ end function field_contains( field::RationalFunctionField{T}, ratfuncs::Vector{Generic.Frac{T}}, - p, + prob_threshold, ) where {T} - return field_contains(field, fractions_to_dennums(ratfuncs), p) + return field_contains(field, fractions_to_dennums(ratfuncs), prob_threshold) end -function field_contains(field::RationalFunctionField{T}, polys::Vector{T}, p) where {T} +function field_contains( + field::RationalFunctionField{T}, + polys::Vector{T}, + prob_threshold, +) where {T} id = one(parent(first(polys))) - return field_contains(field, [[id, p] for p in polys], p) + return field_contains(field, [[id, p] for p in polys], prob_threshold) end # ------------------------------------------------------------------------------ -function issubfield(F::RationalFunctionField{T}, E::RationalFunctionField{T}, p) where {T} - return all(field_contains(E, F.dennums, p)) +function issubfield( + F::RationalFunctionField{T}, + E::RationalFunctionField{T}, + prob_threshold, +) where {T} + return all(field_contains(E, F.dennums, prob_threshold)) end -function fields_equal(F::RationalFunctionField{T}, E::RationalFunctionField{T}, p) where {T} - new_p = 1 - (1 - p) / 2 +function fields_equal( + F::RationalFunctionField{T}, + E::RationalFunctionField{T}, + prob_threshold, +) where {T} + new_p = 1 - (1 - prob_threshold) / 2 return issubfield(F, E, new_p) && issubfield(E, F, new_p) end @@ -475,14 +487,14 @@ function monomial_generators_up_to_degree( end """ - simplified_generating_set(rff; p = 0.99, seed = 42) + simplified_generating_set(rff; prob_threshold = 0.99, seed = 42) Returns a simplified set of generators for `rff`. -Result is correct (in Monte-Carlo sense) with probability at least `p`. +Result is correct (in the Monte-Carlo sense) with probability at least `prob_threshold`. """ @timeit _to function simplified_generating_set( rff::RationalFunctionField; - p = 0.99, + prob_threshold = 0.99, seed = 42, simplify = :standard, check_variables = false, # almost always slows down and thus turned off @@ -500,8 +512,8 @@ Result is correct (in Monte-Carlo sense) with probability at least `p`. # Checking membership of particular variables and adding them to the field if check_variables vars = gens(poly_ring(rff)) - containment = field_contains(rff, vars, (1.0 + p) / 2) - p = (1.0 + p) / 2 + containment = field_contains(rff, vars, (1.0 + prob_threshold) / 2) + prob_threshold = (1.0 + prob_threshold) / 2 if all(containment) return [v // one(poly_ring(rff)) for v in vars] end @@ -557,14 +569,15 @@ Final cleaning and simplification of generators. Out of $(length(new_fracs)) fractions $(length(new_fracs_unique)) are syntactically unique.""" runtime = @elapsed new_fracs = beautifuly_generators(RationalFunctionField(new_fracs_unique)) - @debug "Checking inclusion with probability $p" - runtime = @elapsed result = issubfield(rff, RationalFunctionField(new_fracs), p) + @debug "Checking inclusion with probability $prob_threshold" + runtime = + @elapsed result = issubfield(rff, RationalFunctionField(new_fracs), prob_threshold) _runtime_logger[:id_inclusion_check] = runtime if !result @warn "Field membership check failed. Error will follow." throw("The new subfield generators are not correct.") end - @info "Inclusion checked with probability $p in $(_runtime_logger[:id_inclusion_check]) seconds" + @info "Inclusion checked with probability $prob_threshold in $(_runtime_logger[:id_inclusion_check]) seconds" @debug "Out of $(length(rff.mqs.nums_qq)) initial generators there are $(length(new_fracs)) indepdendent" ranking = generating_set_rank(new_fracs) _runtime_logger[:id_ranking] = ranking diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 165570f3a..86c83993e 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -82,40 +82,44 @@ function __init__() end """ - assess_identifiability(ode; funcs_to_check = [], p=0.99, loglevel=Logging.Info) + assess_identifiability(ode; funcs_to_check = [], prob_threshold=0.99, loglevel=Logging.Info) Input: - `ode` - the ODE model - `funcs_to_check` - list of functions to check identifiability for; if empty, all parameters and states are taken -- `p` - probability of correctness. +- `prob_threshold` - probability of correctness. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) Assesses identifiability of a given ODE model. The result is guaranteed to be correct with the probability -at least `p`. +at least `prob_threshold`. The function returns an (ordered) dictionary from the functions to check to their identifiability properties (one of `:nonidentifiable`, `:locally`, `:globally`). """ function assess_identifiability( ode::ODE{P}; funcs_to_check = Vector(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, loglevel = Logging.Info, ) where {P <: MPolyElem{fmpq}} restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do - return _assess_identifiability(ode, funcs_to_check = funcs_to_check, p = p) + return _assess_identifiability( + ode, + funcs_to_check = funcs_to_check, + prob_threshold = prob_threshold, + ) end end function _assess_identifiability( ode::ODE{P}; funcs_to_check = Vector(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, ) where {P <: MPolyElem{fmpq}} - p_glob = 1 - (1 - p) * 0.9 - p_loc = 1 - (1 - p) * 0.1 + p_glob = 1 - (1 - prob_threshold) * 0.9 + p_loc = 1 - (1 - prob_threshold) * 0.1 if isempty(funcs_to_check) funcs_to_check = vcat(ode.x_vars, ode.parameters) @@ -126,7 +130,7 @@ function _assess_identifiability( runtime = @elapsed local_result = _assess_local_identifiability( ode, funcs_to_check = funcs_to_check, - p = p_loc, + prob_threshold = p_loc, type = :SE, trbasis = trbasis, ) @@ -167,23 +171,23 @@ function _assess_identifiability( end """ - assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], p = 0.99, loglevel=Logging.Info) + assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], prob_threshold = 0.99, loglevel=Logging.Info) Input: - `ode` - the ModelingToolkit.ODESystem object that defines the model - `measured_quantities` - the output functions of the model - `funcs_to_check` - functions of parameters for which to check the identifiability -- `p` - probability of correctness. +- `prob_threshold` - probability of correctness. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) Assesses identifiability (both local and global) of a given ODE model (parameters detected automatically). The result is guaranteed to be correct with the probability -at least `p`. +at least `prob_threshold`. """ function assess_identifiability( ode::ModelingToolkit.ODESystem; measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = [], - p = 0.99, + prob_threshold = 0.99, loglevel = Logging.Info, ) restart_logging(loglevel = loglevel) @@ -192,7 +196,7 @@ function assess_identifiability( ode, measured_quantities = measured_quantities, funcs_to_check = funcs_to_check, - p = p, + prob_threshold = prob_threshold, ) end end @@ -201,7 +205,7 @@ function _assess_identifiability( ode::ModelingToolkit.ODESystem; measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = [], - p = 0.99, + prob_threshold = 0.99, ) if isempty(measured_quantities) measured_quantities = get_measured_quantities(ode) @@ -214,7 +218,11 @@ function _assess_identifiability( end funcs_to_check_ = [eval_at_nemo(each, conversion) for each in funcs_to_check] - result = _assess_identifiability(ode, funcs_to_check = funcs_to_check_, p = p) + result = _assess_identifiability( + ode, + funcs_to_check = funcs_to_check_, + prob_threshold = prob_threshold, + ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) return out_dict diff --git a/src/discrete.jl b/src/discrete.jl index a8bf8da70..193d4dcb3 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -178,11 +178,11 @@ function _degree_with_common_denom(polys) end """ - _assess_local_identifiability_discrete_aux(dds::ODE{P}, funcs_to_check::Array{<: Any, 1}, known_ic, p::Float64=0.99) where P <: MPolyElem{Nemo.fmpq} + _assess_local_identifiability_discrete_aux(dds::ODE{P}, funcs_to_check::Array{<: Any, 1}, known_ic, prob_threshold::Float64=0.99) where P <: MPolyElem{Nemo.fmpq} Checks the local identifiability/observability of the functions in `funcs_to_check` treating `dds` as a discrete-time system with **shift** instead of derivative in the right-hand side. -The result is correct with probability at least `p`. +The result is correct with probability at least `prob_threshold`. `known_ic` can take one of the following * `:none` - no initial conditions are assumed to be known * `:all` - all initial conditions are assumed to be known @@ -192,7 +192,7 @@ function _assess_local_identifiability_discrete_aux( dds::ODE{P}, funcs_to_check::Array{<:Any, 1}, known_ic = :none, - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, ) where {P <: MPolyElem{Nemo.fmpq}} bring = base_ring(dds.poly_ring) @@ -222,7 +222,7 @@ function _assess_local_identifiability_discrete_aux( else Jac_degree += 2 * deg_y * prec end - D = Int(ceil(Jac_degree * length(funcs_to_check) / (1 - p))) + D = Int(ceil(Jac_degree * length(funcs_to_check) / (1 - prob_threshold))) @debug "Sampling range $D" # Parameter values are the same across all the replicas @@ -293,26 +293,26 @@ end measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=Array{}[], known_ic=Array{}[], - p::Float64=0.99) + prob_threshold::Float64=0.99) Input: - `dds` - the DiscreteSystem object from ModelingToolkit (with **difference** operator in the right-hand side) - `measured_quantities` - the measurable outputs of the model - `funcs_to_check` - functions of parameters for which to check identifiability (all parameters and states if not specified) - `known_ic` - functions (of states and parameter) whose initial conditions are assumed to be known -- `p` - probability of correctness +- `prob_threshold` - probability of correctness Output: - the result is an (ordered) dictionary from each function to to boolean; -The result is correct with probability at least `p`. +The result is correct with probability at least `prob_threshold`. """ function assess_local_identifiability( dds::ModelingToolkit.DiscreteSystem; measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = Array{}[], known_ic = Array{}[], - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, loglevel = Logging.Info, ) restart_logging(loglevel = loglevel) @@ -322,7 +322,7 @@ function assess_local_identifiability( measured_quantities = measured_quantities, funcs_to_check = funcs_to_check, known_ic = known_ic, - p = p, + prob_threshold = prob_threshold, ) end end @@ -332,7 +332,7 @@ function _assess_local_identifiability( measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = Array{}[], known_ic = Array{}[], - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, ) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(dds)) @@ -373,8 +373,12 @@ function _assess_local_identifiability( funcs_to_check_ = [eval_at_nemo(x, conversion) for x in funcs_to_check] known_ic_ = [eval_at_nemo(x, conversion) for x in known_ic] - result = - _assess_local_identifiability_discrete_aux(dds_aux, funcs_to_check_, known_ic_, p) + result = _assess_local_identifiability_discrete_aux( + dds_aux, + funcs_to_check_, + known_ic_, + prob_threshold, + ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = OrderedDict(nemo2mtk[param] => result[param] for param in funcs_to_check_) if length(known_ic) > 0 diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 5f61644e0..984747e3b 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -88,7 +88,7 @@ These are presented by the coefficients of the IO-equations. ## Options This function takes the following optional arguments: -- `p`: probability of correctness +- `prob_threshold`: probability of correctness - `with_states`: Also report the identifiabile functions in states. Default is `false`. If this is `true`, the identifiable functions involving parameters only will be simplified @@ -99,7 +99,7 @@ The function returns a tuple containing the following: """ @timeit _to function initial_identifiable_functions( ode::ODE{T}; - p::Float64, + prob_threshold::Float64, known::Array{T, 1} = Array{T, 1}(), with_states::Bool = false, var_change_policy = :default, @@ -157,7 +157,7 @@ The function returns a tuple containing the following: _runtime_logger[:check_time] = @elapsed no_states_simplified = simplified_generating_set( RationalFunctionField(id_funcs_no_states_param), - p = p, + prob_threshold = prob_threshold, seed = 42, simplify = :standard, rational_interpolator = rational_interpolator, @@ -178,13 +178,13 @@ end # ------------------------------------------------------------------------------ """ - check_identifiability(ode, funcs_to_check; known, p, var_change_policy) + check_identifiability(ode, funcs_to_check; known, prob_threshold, var_change_policy) Input: - `ode` - the ODE model - `funcs_to_check` - the functions to check identifiability for - `known` - a list of functions in states which are assumed to be known and generic -- `p` - probability of correctness +- `prob_threshold` - probability of correctness - `var_change` - a policy for variable change (`:default`, `:yes`, `:no`), affects only the runtime Output: a list L of booleans with L[i] being the identifiability status of the i-th function to check @@ -193,7 +193,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i ode::ODE{P}, funcs_to_check::Array{<:Any, 1}; known::Array{P, 1} = Array{P, 1}(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, var_change_policy = :default, ) where {P <: MPolyElem{fmpq}} states_needed = false @@ -209,11 +209,11 @@ Output: a list L of booleans with L[i] being the identifiability status of the i return [true for _ in funcs_to_check] end - half_p = 0.5 + p / 2 + half_p = 0.5 + prob_threshold / 2 identifiable_functions_raw, bring = initial_identifiable_functions( ode, known = known, - p = p, + prob_threshold = half_p, var_change_policy = var_change_policy, with_states = states_needed, ) @@ -234,26 +234,26 @@ end function check_identifiability( ode::ODE{P}; known::Array{P, 1} = Array{P, 1}(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, var_change_policy = :default, ) where {P <: MPolyElem{fmpq}} return check_identifiability( ode, ode.parameters, known = known, - p = p, + prob_threshold = prob_threshold, var_change_policy = var_change_policy, ) end #------------------------------------------------------------------------------ """ - assess_global_identifiability(ode::ODE{P}, p::Float64=0.99; var_change=:default) where P <: MPolyElem{fmpq} + assess_global_identifiability(ode::ODE{P}, prob_threshold::Float64=0.99; var_change=:default) where P <: MPolyElem{fmpq} Input: - `ode` - the ODE model - `known` - a list of functions in states which are assumed to be known and generic -- `p` - probability of correctness +- `prob_threshold` - probability of correctness - `var_change` - a policy for variable change (`:default`, `:yes`, `:no`), affects only the runtime Output: @@ -264,14 +264,14 @@ Checks global identifiability for parameters of the model provided in `ode`. Cal function assess_global_identifiability( ode::ODE{P}, known::Array{P, 1} = Array{P, 1}(), - p::Float64 = 0.99; + prob_threshold::Float64 = 0.99; var_change = :default, ) where {P <: MPolyElem{fmpq}} result_list = assess_global_identifiability( ode, ode.parameters, known, - p; + prob_threshold; var_change = var_change, ) @@ -281,13 +281,13 @@ end #------------------------------------------------------------------------------ """ - assess_global_identifiability(ode, [funcs_to_check, p=0.99, var_change=:default]) + assess_global_identifiability(ode, [funcs_to_check, prob_threshold=0.99, var_change=:default]) Input: - `ode` - the ODE model - `funcs_to_check` - rational functions in parameters - `known` - function in parameters that are assumed to be known and generic -- `p` - probability of correctness +- `prob_threshold` - probability of correctness - `var_change` - a policy for variable change (`:default`, `:yes`, `:no`), affects only the runtime @@ -301,7 +301,7 @@ Checks global identifiability of functions of parameters specified in `funcs_to_ ode::ODE{P}, funcs_to_check::Array{<:Any, 1}, known::Array{P, 1} = Array{P, 1}(), - p::Float64 = 0.99; + prob_threshold::Float64 = 0.99; var_change = :default, ) where {P <: MPolyElem{fmpq}} submodels = find_submodels(ode) @@ -313,7 +313,7 @@ Checks global identifiability of functions of parameters specified in `funcs_to_ ode, funcs_to_check, known = known, - p = p, + prob_threshold = prob_threshold, var_change_policy = var_change, ) diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 4edb29b83..79db0ffbd 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -18,7 +18,7 @@ This functions takes the following optional arguments: - `:strong`: Strong simplification. This option is the slowest, but the output functions are nice and simple. - `:absent`: No simplification. -- `p`: A float in the range from 0 to 1, the probability of correctness. Default +- `prob_threshold`: A float in the range from 0 to 1, the probability of correctness. Default is `0.99`. - `seed`: The rng seed. Default value is `42`. - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) @@ -45,7 +45,7 @@ find_identifiable_functions(ode) """ function find_identifiable_functions( ode::ODE{T}; - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, seed = 42, with_states = false, simplify = :standard, @@ -57,7 +57,7 @@ function find_identifiable_functions( with_logger(_si_logger[]) do return _find_identifiable_functions( ode, - p = p, + prob_threshold = prob_threshold, seed = seed, with_states = with_states, simplify = simplify, @@ -68,7 +68,7 @@ end function _find_identifiable_functions( ode::ODE{T}; - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, seed = 42, with_states = false, simplify = :standard, @@ -87,10 +87,10 @@ function _find_identifiable_functions( id_funcs = [one(bring)] return id_funcs end - half_p = 0.5 + p / 2 + half_p = 0.5 + prob_threshold / 2 id_funcs, bring = initial_identifiable_functions( ode, - p = half_p, + prob_threshold = half_p, with_states = with_states, rational_interpolator = rational_interpolator, ) @@ -102,7 +102,7 @@ function _find_identifiable_functions( end id_funcs_fracs = simplified_generating_set( RationalFunctionField(id_funcs), - p = half_p, + prob_threshold = half_p, seed = seed, simplify = simplify, rational_interpolator = rational_interpolator, @@ -158,7 +158,7 @@ find_identifiable_functions(de, measured_quantities = [y1 ~ x0]) function find_identifiable_functions( ode::ModelingToolkit.ODESystem; measured_quantities = Array{ModelingToolkit.Equation}[], - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, seed = 42, with_states = false, simplify = :standard, @@ -171,7 +171,7 @@ function find_identifiable_functions( return _find_identifiable_functions( ode, measured_quantities = measured_quantities, - p = p, + prob_threshold = prob_threshold, seed = seed, with_states = with_states, simplify = simplify, @@ -183,7 +183,7 @@ end function _find_identifiable_functions( ode::ModelingToolkit.ODESystem; measured_quantities = Array{ModelingToolkit.Equation}[], - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, seed = 42, with_states = false, simplify = :standard, @@ -197,7 +197,7 @@ function _find_identifiable_functions( result = _find_identifiable_functions( ode, simplify = simplify, - p = p, + prob_threshold = prob_threshold, seed = seed, with_states = with_states, rational_interpolator = rational_interpolator, diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 3d680c336..ae2a3cbba 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -143,13 +143,13 @@ end # ------------------------------------------------------------------------------ """ - function assess_local_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=Array{}[], p::Float64=0.99, type=:SE, loglevel=Logging.Info) + function assess_local_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=Array{}[], prob_threshold::Float64=0.99, type=:SE, loglevel=Logging.Info) Input: - `ode` - the ODESystem object from ModelingToolkit - `measured_quantities` - the measurable outputs of the model - `funcs_to_check` - functions of parameters for which to check identifiability -- `p` - probability of correctness +- `prob_threshold` - probability of correctness - `type` - identifiability type (`:SE` for single-experiment, `:ME` for multi-experiment) - `loglevel` - the minimal level of log messages to display (`Logging.Info` by default) @@ -159,7 +159,7 @@ Output: The function determines local identifiability of parameters in `funcs_to_check` or all possible parameters if `funcs_to_check` is empty -The result is correct with probability at least `p`. +The result is correct with probability at least `prob_threshold`. `type` can be either `:SE` (single-experiment identifiability) or `:ME` (multi-experiment identifiability). The return value is a tuple consisting of the array of bools and the number of experiments to be performed. @@ -168,7 +168,7 @@ function assess_local_identifiability( ode::ModelingToolkit.ODESystem; measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = Array{}[], - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, type = :SE, loglevel = Logging.Info, ) @@ -178,7 +178,7 @@ function assess_local_identifiability( ode, measured_quantities = measured_quantities, funcs_to_check = funcs_to_check, - p = p, + prob_threshold = prob_threshold, type = type, ) end @@ -188,7 +188,7 @@ end ode::ModelingToolkit.ODESystem; measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = Array{}[], - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, type = :SE, ) if length(measured_quantities) == 0 @@ -219,7 +219,7 @@ end result = _assess_local_identifiability( ode, funcs_to_check = funcs_to_check_, - p = p, + prob_threshold = prob_threshold, type = type, ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) @@ -230,7 +230,7 @@ end result, bd = _assess_local_identifiability( ode, funcs_to_check = funcs_to_check_, - p = p, + prob_threshold = prob_threshold, type = type, ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) @@ -242,9 +242,9 @@ end # ------------------------------------------------------------------------------ """ - assess_local_identifiability(ode::ODE{P}; funcs_to_check::Array{<: Any, 1}, p::Float64=0.99, type=:SE, loglevel=Logging.Info) where P <: MPolyElem{Nemo.fmpq} + assess_local_identifiability(ode::ODE{P}; funcs_to_check::Array{<: Any, 1}, prob_threshold::Float64=0.99, type=:SE, loglevel=Logging.Info) where P <: MPolyElem{Nemo.fmpq} -Checks the local identifiability/observability of the functions in `funcs_to_check`. The result is correct with probability at least `p`. +Checks the local identifiability/observability of the functions in `funcs_to_check`. The result is correct with probability at least `prob_threshold`. Call this function if you have a specific collection of parameters of which you would like to check local identifiability. @@ -254,7 +254,7 @@ If the type is `:ME`, states are not allowed to appear in the `funcs_to_check`. function assess_local_identifiability( ode::ODE{P}; funcs_to_check::Array{<:Any, 1} = Array{Any, 1}(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, type = :SE, trbasis = nothing, loglevel = Logging.Info, @@ -265,7 +265,7 @@ function assess_local_identifiability( return _assess_local_identifiability( ode, funcs_to_check = funcs_to_check, - p = p, + prob_threshold = prob_threshold, type = type, trbasis = trbasis, ) @@ -275,7 +275,7 @@ end function _assess_local_identifiability( ode::ODE{P}; funcs_to_check::Array{<:Any, 1} = Array{Any, 1}(), - p::Float64 = 0.99, + prob_threshold::Float64 = 0.99, type = :SE, trbasis = nothing, ) where {P <: MPolyElem{Nemo.fmpq}} @@ -307,7 +307,7 @@ function _assess_local_identifiability( d = max(d, df) h = max(h, hf) end - p_per_func = 1 - (1 - p) / length(funcs_to_check) + p_per_func = 1 - (1 - prob_threshold) / length(funcs_to_check) mu = ceil(1 / (1 - sqrt(p_per_func))) n = length(ode.x_vars) diff --git a/src/parametrizations.jl b/src/parametrizations.jl index aef46cc12..a8b6e6ec4 100644 --- a/src/parametrizations.jl +++ b/src/parametrizations.jl @@ -432,7 +432,7 @@ Returns a tuple (`new_ode`, `new_vars`, `implicit_relations`), such that: The function accepts the following optional arguments. - `seed`: A float in the range from 0 to 1, random seed (default is `seed = 42`). -- `p`: The probability of correctness (default is `p = 0.99`). +- `prob_threshold`: The probability of correctness (default is `prob_threshold = 0.99`). ## Example @@ -471,20 +471,24 @@ compared to the original one. """ function reparametrize_global( ode::ODE{P}; - p = 0.99, + prob_threshold = 0.99, seed = 42, loglevel = Logging.Info, ) where {P} restart_logging(loglevel = loglevel) with_logger(_si_logger[]) do - return _reparametrize_global(ode, p = p, seed = seed) + return _reparametrize_global(ode, prob_threshold = prob_threshold, seed = seed) end end -function _reparametrize_global(ode::ODE{P}; p = 0.99, seed = 42) where {P} +function _reparametrize_global(ode::ODE{P}; prob_threshold = 0.99, seed = 42) where {P} Random.seed!(seed) - id_funcs = - find_identifiable_functions(ode, with_states = true, simplify = :strong, p = p) + id_funcs = find_identifiable_functions( + ode, + with_states = true, + simplify = :strong, + prob_threshold = prob_threshold, + ) ode_ring = parent(ode) @assert base_ring(parent(first(id_funcs))) == ode_ring @info "Constructing a new parametrization" diff --git a/test/local_identifiability_me.jl b/test/local_identifiability_me.jl index b64c0bb3a..01365120c 100644 --- a/test/local_identifiability_me.jl +++ b/test/local_identifiability_me.jl @@ -213,7 +213,7 @@ end result = assess_local_identifiability( case[:ode], funcs_to_check = case[:funcs], - p = 0.932, + prob_threshold = 0.932, type = :ME, ) @test result == case[:correct] diff --git a/test/mtk_compat.jl b/test/mtk_compat.jl index 0927500b9..ccbdaf287 100644 --- a/test/mtk_compat.jl +++ b/test/mtk_compat.jl @@ -149,7 +149,7 @@ de; measured_quantities = measured_quantities, funcs_to_check = funcs_to_check, - p = 0.99, + prob_threshold = 0.99, type = :ME, ), ) @@ -195,7 +195,7 @@ assess_local_identifiability( de; funcs_to_check = funcs_to_check, - p = 0.99, + prob_threshold = 0.99, type = :ME, ), ) @@ -232,7 +232,7 @@ de; measured_quantities = measured_quantities, funcs_to_check = funcs_to_check, - p = 0.99, + prob_threshold = 0.99, type = :ME, ), )