From a67b9550612f544751048f5e55d2e42daf529915 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Fri, 29 Sep 2023 07:55:35 +0200 Subject: [PATCH] systems without parameters fix --- src/StructuralIdentifiability.jl | 11 +++---- src/global_identifiability.jl | 56 ++++++++++++++++++-------------- src/local_identifiability.jl | 1 + test/identifiable_functions.jl | 12 +++---- 4 files changed, 42 insertions(+), 38 deletions(-) diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index d88e302b3..bedb3a4b8 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -81,8 +81,6 @@ include("pb_representation.jl") include("submodels.jl") include("discrete.jl") -# TODO: handle finding identifiabile functions in case there are no parameters - """ assess_identifiability(ode::ODE{P}, p::Float64=0.99) where P <: MPolyElem{fmpq} @@ -108,11 +106,10 @@ Input: Assesses identifiability of a given ODE model. The result is guaranteed to be correct with the probability at least `p`. -If `funcs_to_check` are given, then the function will assess the identifiability of the provided functions -and return a list of the same length with each element being one of `:nonidentifiable`, `:locally`, `:globally`. - -If `funcs_to_check` are not given, the function will assess identifiability of the parameters, and the result will -be a dictionary from the parameters to their identifiability properties (again, one of `:nonidentifiable`, `:locally`, `:globally`). +If `funcs_to_check` are given, then the function will assess the identifiability of the provided functions. +Otherwise, it will assess identifiability of all the parameters and states. +The function returns a dictionary from the functions to check to their identifiability properties +(one of `:nonidentifiable`, `:locally`, `:globally`). """ function assess_identifiability( ode::ODE{P}, diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index d44bfa53a..01f49c689 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -22,14 +22,16 @@ function extract_identifiable_functions_raw( @debug "Extracting coefficients" flush(stdout) - nonparameters = filter( - v -> !(var_to_str(v) in map(var_to_str, ode.parameters)), - gens(parent(first(values(io_equations)))), - ) - for eq in values(io_equations) - eq_coefs = collect(values(extract_coefficients(eq, nonparameters))) - eq_coefs = [parent_ring_change(c, bring) for c in eq_coefs] - push!(coeff_lists[:no_states], eq_coefs) + if !isempty(ode.parameters) + nonparameters = filter( + v -> !(var_to_str(v) in map(var_to_str, ode.parameters)), + gens(parent(first(values(io_equations)))), + ) + for eq in values(io_equations) + eq_coefs = collect(values(extract_coefficients(eq, nonparameters))) + eq_coefs = [parent_ring_change(c, bring) for c in eq_coefs] + push!(coeff_lists[:no_states], eq_coefs) + end end for p in known @@ -93,22 +95,26 @@ function initial_identifiable_functions( @info "Computed in $ioeq_time seconds" :ioeq_time ioeq_time _runtime_logger[:ioeq_time] = ioeq_time - @info "Computing Wronskians" - wrnsk_time = @elapsed wrnsk = wronskian(io_equations, ode) - @info "Computed in $wrnsk_time seconds" :wrnsk_time wrnsk_time - _runtime_logger[:wrnsk_time] = wrnsk_time - - dims = map(ncols, wrnsk) - @info "Dimensions of the Wronskians $dims" - - rank_times = @elapsed wranks = map(rank, wrnsk) - @debug "Dimensions of the Wronskians $dims" - @debug "Ranks of the Wronskians $wranks" - @info "Ranks of the Wronskians computed in $rank_times seconds" :rank_time rank_times - _runtime_logger[:rank_time] = rank_times - - if any([dim != rk + 1 for (dim, rk) in zip(dims, wranks)]) - @warn "One of the Wronskians has corank greater than one, so the results of the algorithm will be valid only for multiexperiment identifiability. If you still would like to assess single-experiment identifiability, we recommend using SIAN (https://github.com/alexeyovchinnikov/SIAN-Julia)" + if isempty(ode.parameters) + @info "No parameters, so Wronskian computation is not needed" + else + @info "Computing Wronskians" + wrnsk_time = @elapsed wrnsk = wronskian(io_equations, ode) + @info "Computed in $wrnsk_time seconds" :wrnsk_time wrnsk_time + _runtime_logger[:wrnsk_time] = wrnsk_time + + dims = map(ncols, wrnsk) + @info "Dimensions of the Wronskians $dims" + + rank_times = @elapsed wranks = map(rank, wrnsk) + @debug "Dimensions of the Wronskians $dims" + @debug "Ranks of the Wronskians $wranks" + @info "Ranks of the Wronskians computed in $rank_times seconds" :rank_time rank_times + _runtime_logger[:rank_time] = rank_times + + if any([dim != rk + 1 for (dim, rk) in zip(dims, wranks)]) + @warn "One of the Wronskians has corank greater than one, so the results of the algorithm will be valid only for multiexperiment identifiability. If you still would like to assess single-experiment identifiability, we recommend using SIAN (https://github.com/alexeyovchinnikov/SIAN-Julia)" + end end id_funcs, bring = extract_identifiable_functions_raw( @@ -118,7 +124,7 @@ function initial_identifiable_functions( with_states, ) - if with_states + if with_states && !isempty(ode.parameters) @debug "Generators of identifiable functions involve states, the parameter-only part is getting simplified" # NOTE: switching to a ring without states for a moment param_ring, _ = PolynomialRing( diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 790e7d208..91f4a60c2 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -297,6 +297,7 @@ function assess_local_identifiability( Dprime = D * (2 * log(n + ell + r + 1) + log(mu * D)) + 4 * (n + ell)^2 * ((n + m) * h + log(2 * n * D)) + Dprime = max(Dprime, 1.0) prime = Primes.nextprime(Int(ceil(2 * mu * Dprime))) @debug "The prime is $prime" F = Nemo.GF(prime) diff --git a/test/identifiable_functions.jl b/test/identifiable_functions.jl index 66bce75f0..854a55b18 100644 --- a/test/identifiable_functions.jl +++ b/test/identifiable_functions.jl @@ -759,16 +759,16 @@ ident_funcs = [ phi * Mar * siga2 * beta_SA - phi * siga2^2 * beta_SA - Mar * siga2 * beta_SA ) // (phi * M * siga2), ] -push!(test_cases, (ode = ode, ident_funcs = ident_funcs)) +# Really large and takes a lot of time, so commented +# push!(test_cases, (ode = ode, ident_funcs = ident_funcs)) ### # Cases with states -# TODO: uncomment when this is handled! -# ode = StructuralIdentifiability.@ODEmodel(x'(t) = x(t), y(t) = x(t)) -# T = typeof(x) -# ident_funcs = [x] -# push!(test_cases, (ode = ode, ident_funcs = ident_funcs)) +ode = StructuralIdentifiability.@ODEmodel(x'(t) = x(t), y(t) = x(t)) +T = typeof(x) +ident_funcs = [x] +push!(test_cases, (ode = ode, ident_funcs = ident_funcs, with_states = true)) ode = StructuralIdentifiability.@ODEmodel(x'(t) = a * x(t) + u(t), y(t) = x(t)) ident_funcs = [a, x]