systems without parameters fix

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},

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(

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)

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]