From 32ad4f5e0c7d0f4cdb51d39248cb5d3d51a1174a Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Fri, 10 Nov 2023 00:49:50 +0100 Subject: [PATCH 01/24] getting started --- src/ODE.jl | 20 ++++++---- src/StructuralIdentifiability.jl | 28 ++++++------- src/identifiable_functions.jl | 6 +-- src/logging.jl | 68 ++++++++++++++++++++++++++++++++ src/precompile.jl | 23 ++++++----- 5 files changed, 105 insertions(+), 40 deletions(-) create mode 100644 src/logging.jl diff --git a/src/ODE.jl b/src/ODE.jl index fc05b91d9..945ae53b4 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -25,7 +25,7 @@ struct ODE{P} # P is the type of polynomials in the rhs of the ODE system # x_eqs is a dictionary x_i => f_i(x, u, params) # y_eqs is a dictionary y_i => g_i(x, u, params) if isempty(y_eqs) - @info "Could not find output variables in the model." + info_si("Could not find output variables in the model.") end poly_ring = parent(first(vcat(y_vars, x_vars))) u_vars = inputs @@ -458,14 +458,19 @@ macro ODEmodel(ex::Expr...) ) end end + + info_si("Summary of the model:") + info_si("State variables: " * join(map(string, collect(x_vars)), ", ")) + info_si("Parameters: " * join(map(string, collect(params)), ", ")) + info_si("Inputs: " * join(map(string, collect(u_vars)), ", ")) + info_si("Outputs: " * join(map(string, collect(y_vars)), ", ")) logging_exprs = [ - :(@info "Summary of the model:"), - :(@info "State variables: " * $(join(map(string, collect(x_vars)), ", "))), - :(@info "Parameters: " * $(join(map(string, collect(params)), ", "))), - :(@info "Inputs: " * $(join(map(string, collect(u_vars)), ", "))), - :(@info "Outputs: " * $(join(map(string, collect(y_vars)), ", "))), + :(StructuralIdentifiability.info_si("Summary of the model:")), + :(StructuralIdentifiability.info_si("State variables: " * $(join(map(string, collect(x_vars)), ", ")))), + :(StructuralIdentifiability.info_si("Parameters: " * $(join(map(string, collect(params)), ", ")))), + :(StructuralIdentifiability.info_si("Inputs: " * $(join(map(string, collect(u_vars)), ", ")))), + :(StructuralIdentifiability.info_si("Outputs: " * $(join(map(string, collect(y_vars)), ", ")))), ] - # creating the ode object ode_expr = :(StructuralIdentifiability.ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}( $vx, @@ -570,7 +575,6 @@ function __preprocess_ode( de::ModelingToolkit.AbstractTimeDependentSystem, measured_quantities::Array{<:Tuple{String, <:SymbolicUtils.BasicSymbolic}}, ) - @info "Preproccessing `ModelingToolkit.AbstractTimeDependentSystem` object" polytype = StructuralIdentifiability.Nemo.fmpq_mpoly fractype = StructuralIdentifiability.Nemo.Generic.Frac{polytype} diff_eqs = diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 2312e88ca..75a9417db 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -25,7 +25,7 @@ using ModelingToolkit export ODE, @ODEmodel, preprocess_ode # assessing identifiability -export assess_local_identifiability, assess_global_identifiability, assess_identifiability +export assess_local_identifiability, assess_identifiability # auxuliary function export set_parameter_values @@ -48,16 +48,7 @@ export find_submodels # finding identifiabile reparametrizations export reparametrize_global -# would be great to merge with the Julia logger -const _runtime_logger = Dict( - :id_calls_to_gb => 0, - :id_groebner_time => 0.0, - :id_inclusion_check_mod_p => 0, - :id_npoints_degree => 0, - :id_npoints_interpolation => 0, - :id_beautifulization => 0, -) - +include("logging.jl") include("util.jl") include("power_series_utils.jl") include("ODE.jl") @@ -99,7 +90,9 @@ function assess_identifiability( ode::ODE{P}; funcs_to_check = Vector(), p::Float64 = 0.99, + loglevel = Logging.Info ) where {P <: MPolyElem{fmpq}} + restart_logging(loglevel = loglevel) p_glob = 1 - (1 - p) * 0.9 p_loc = 1 - (1 - p) * 0.1 @@ -107,7 +100,7 @@ function assess_identifiability( funcs_to_check = vcat(ode.parameters, ode.x_vars) end - @info "Assessing local identifiability" + info_si("Assessing local identifiability") trbasis = Array{fmpq_mpoly, 1}() runtime = @elapsed local_result = assess_local_identifiability( ode, @@ -116,8 +109,8 @@ function assess_identifiability( type = :SE, trbasis = trbasis, ) - @info "Local identifiability assessed in $runtime seconds" - @debug "Trasncendence basis to be specialized is $trbasis" + info_si("Local identifiability assessed in $runtime seconds") + debug_si("Trasncendence basis to be specialized is $trbasis") _runtime_logger[:loc_time] = runtime loc_id = [local_result[each] for each in funcs_to_check] @@ -128,10 +121,10 @@ function assess_identifiability( end end - @info "Assessing global identifiability" + info_si("Assessing global identifiability") runtime = @elapsed global_result = assess_global_identifiability(ode, locally_identifiable, trbasis, p_glob) - @info "Global identifiability assessed in $runtime seconds" + info_si("Global identifiability assessed in $runtime seconds") _runtime_logger[:glob_time] = runtime result = Dict{Any, Symbol}() @@ -169,6 +162,7 @@ function assess_identifiability( measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = [], p = 0.99, + loglevel = Logging.Info ) if isempty(measured_quantities) measured_quantities = get_measured_quantities(ode) @@ -181,7 +175,7 @@ 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_, p = p, loglevel = loglevel) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) return out_dict diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 8045fccf3..feb404e42 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -49,12 +49,10 @@ function find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + loglevel = Logging.Info, ) where {T <: MPolyElem{fmpq}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) - _runtime_logger[:id_npoints_degree] = 0 - _runtime_logger[:id_npoints_interpolation] = 0 - _runtime_logger[:id_beautifulization] = 0.0 runtime_start = time_ns() if isempty(ode.parameters) && !with_states @warn """ @@ -139,6 +137,7 @@ function find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, + loglevel = Logging.Info, ) Random.seed!(seed) if isempty(measured_quantities) @@ -152,6 +151,7 @@ function find_identifiable_functions( seed = seed, with_states = with_states, rational_interpolator = rational_interpolator, + loglevel = loglevel, ) result = [parent_ring_change(f, ode.poly_ring) for f in result] nemo2mtk = Dict(v => Num(k) for (k, v) in conversion) diff --git a/src/logging.jl b/src/logging.jl new file mode 100644 index 000000000..946163c21 --- /dev/null +++ b/src/logging.jl @@ -0,0 +1,68 @@ +# Logging used by StructuralIdentifiability.jl +# Consists of two parts: +# - si_logger (using Logging package) for logs +# - runtime_logger dictionary for recording runtime info +# +# Defines functions to refresh both + specific logging functions +# info_si, debug_si, warn_si, and error_si calling @info, @debug, @warn, and @error, respectively + +# not the exhaustive list, just the ones to be initialized +const _runtime_rubrics = ( + :id_calls_to_gb, + :id_groebner_time, + :id_inclusion_check_mod_p, + :id_npoints_degree, + :id_npoints_interpolation, + :id_beautifulization, +) + +const _runtime_logger = Dict( + :id_calls_to_gb => 0, + :id_groebner_time => 0.0, + :id_inclusion_check_mod_p => 0, + :id_npoints_degree => 0, + :id_npoints_interpolation => 0, + :id_beautifulization => 0, +) + +const _si_logger = Ref{Logging.ConsoleLogger}( + Logging.ConsoleLogger( + Logging.Info, + show_limited = false, + ) +) + +function restart_logging(; loglevel = Logging.Info) + _si_logger[] = Logging.ConsoleLogger(loglevel, show_limited = false) + for r in _runtime_rubrics + _runtime_logger[r] = 0 + end +end + +function info_si(s) + with_logger(_si_logger[]) do + @info s + flush(stdout) + end +end + +function debug_si(s) + with_logger(_si_logger[]) do + @debug s + flush(stdout) + end +end + +function warn_si(s) + with_logger(_si_logger[]) do + @warn s + flush(stdout) + end +end + +function error_si(s) + with_logger(_si_logger[]) do + @error s + flush(stdout) + end +end diff --git a/src/precompile.jl b/src/precompile.jl index ba52f3346..e26244c5f 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -12,17 +12,16 @@ @compile_workload begin # all calls in this block will be precompiled, regardless of whether # they belong to your package or not (on Julia 1.8 and higher) - with_logger(ConsoleLogger(stdout, Logging.Warn)) do - ode = @ODEmodel( - x1'(t) = -(a01 + a21) * x1(t) + a12 * x2(t) + u(t), - x2'(t) = a21 * x1(t) - a12 * x2(t) - x3(t) / b, - x3'(t) = x3(t), - y(t) = x2(t) - ) - assess_identifiability(ode) - assess_identifiability(de; measured_quantities = [x0]) - assess_identifiability(de; measured_quantities = [y ~ x0]) - find_identifiable_functions(ode, with_states = true) - end + restart_logging(loglevel = Logging.Warn) + ode = @ODEmodel( + x1'(t) = -(a01 + a21) * x1(t) + a12 * x2(t) + u(t), + x2'(t) = a21 * x1(t) - a12 * x2(t) - x3(t) / b, + x3'(t) = x3(t), + y(t) = x2(t) + ) + assess_identifiability(ode) + assess_identifiability(de; measured_quantities = [x0]) + assess_identifiability(de; measured_quantities = [y ~ x0]) + find_identifiable_functions(ode, with_states = true) end end From ad5d6865467b2f00b89ced9ec4346d525e6fe9f8 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Fri, 10 Nov 2023 00:54:43 +0100 Subject: [PATCH 02/24] adding scope --- src/precompile.jl | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/src/precompile.jl b/src/precompile.jl index e26244c5f..af84e317b 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -10,18 +10,20 @@ eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1] de = ODESystem(eqs, t, name = :Test) @compile_workload begin - # all calls in this block will be precompiled, regardless of whether - # they belong to your package or not (on Julia 1.8 and higher) - restart_logging(loglevel = Logging.Warn) - ode = @ODEmodel( - x1'(t) = -(a01 + a21) * x1(t) + a12 * x2(t) + u(t), - x2'(t) = a21 * x1(t) - a12 * x2(t) - x3(t) / b, - x3'(t) = x3(t), - y(t) = x2(t) - ) - assess_identifiability(ode) - assess_identifiability(de; measured_quantities = [x0]) - assess_identifiability(de; measured_quantities = [y ~ x0]) - find_identifiable_functions(ode, with_states = true) + with_logger(Logging.ConsoleLogger(Logging.Warn)) do + # all calls in this block will be precompiled, regardless of whether + # they belong to your package or not (on Julia 1.8 and higher) + restart_logging(loglevel = Logging.Warn) + ode = @ODEmodel( + x1'(t) = -(a01 + a21) * x1(t) + a12 * x2(t) + u(t), + x2'(t) = a21 * x1(t) - a12 * x2(t) - x3(t) / b, + x3'(t) = x3(t), + y(t) = x2(t) + ) + assess_identifiability(ode) + assess_identifiability(de; measured_quantities = [x0]) + assess_identifiability(de; measured_quantities = [y ~ x0]) + find_identifiable_functions(ode, with_states = true) + end end end From 8ed0c9e470ca96690dbda624295c63e5f58d2603 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Fri, 10 Nov 2023 12:20:09 +0100 Subject: [PATCH 03/24] precompilation fixed --- src/ODE.jl | 5 ----- src/precompile.jl | 1 + 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 945ae53b4..887cdf95c 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -459,11 +459,6 @@ macro ODEmodel(ex::Expr...) end end - info_si("Summary of the model:") - info_si("State variables: " * join(map(string, collect(x_vars)), ", ")) - info_si("Parameters: " * join(map(string, collect(params)), ", ")) - info_si("Inputs: " * join(map(string, collect(u_vars)), ", ")) - info_si("Outputs: " * join(map(string, collect(y_vars)), ", ")) logging_exprs = [ :(StructuralIdentifiability.info_si("Summary of the model:")), :(StructuralIdentifiability.info_si("State variables: " * $(join(map(string, collect(x_vars)), ", ")))), diff --git a/src/precompile.jl b/src/precompile.jl index af84e317b..415dc0577 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -24,6 +24,7 @@ assess_identifiability(de; measured_quantities = [x0]) assess_identifiability(de; measured_quantities = [y ~ x0]) find_identifiable_functions(ode, with_states = true) + restart_logging(loglevel = Logging.Info) end end end From 8825005e7b84720192a786c557bbc6c3fad4379a Mon Sep 17 00:00:00 2001 From: gpogudin Date: Fri, 10 Nov 2023 15:11:28 +0100 Subject: [PATCH 04/24] applying logging --- src/RationalFunctionFields/IdealMQS.jl | 32 ++++--- .../RationalFunctionField.jl | 76 ++++++++-------- src/RationalFunctionFields/normalforms.jl | 42 ++++----- src/StructuralIdentifiability.jl | 5 +- src/discrete.jl | 24 ++--- src/elimination.jl | 49 ++++------- src/global_identifiability.jl | 35 ++++---- src/identifiable_functions.jl | 9 +- src/io_equation.jl | 88 ++++++++----------- src/lincomp.jl | 4 +- src/local_identifiability.jl | 40 +++++---- src/parametrizations.jl | 67 +++++++------- src/pb_representation.jl | 2 +- src/power_series_utils.jl | 2 +- src/precompile.jl | 8 +- src/primality_check.jl | 10 +-- src/util.jl | 17 ++-- src/wronskian.jl | 10 +-- 18 files changed, 244 insertions(+), 276 deletions(-) diff --git a/src/RationalFunctionFields/IdealMQS.jl b/src/RationalFunctionFields/IdealMQS.jl index bd4d7aabd..8a2e760bc 100644 --- a/src/RationalFunctionFields/IdealMQS.jl +++ b/src/RationalFunctionFields/IdealMQS.jl @@ -63,11 +63,11 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal # We prepare and store them to construct ideal specializations later @assert !isempty(funcs_den_nums) @assert sat_var_position in (:first, :last) - ordering !== :degrevlex && (@warn "Ordering is not degrevlex but $ordering") + ordering !== :degrevlex && (warn_si("Ordering is not degrevlex but $ordering")) ring = parent(first(first(funcs_den_nums))) - @debug "Constructing the MQS ideal in $ring" + debug_si("Constructing the MQS ideal in $ring") K, n = base_ring(ring), nvars(ring) - @debug "Finding pivot polynomials" + debug_si("Finding pivot polynomials") # In the component f1,f2,... find the polynomial with the minimal total # degree and length. Such element will serve as a normalizing term for # the component @@ -78,17 +78,17 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal map(plist -> findmin(p -> (total_degree(p), length(p)), plist), funcs_den_nums) pivots_indices = map(last, pivots) for plist in funcs_den_nums - @debug "\tDegrees in this list are $(map(total_degree, plist))" + debug_si("\tDegrees in this list are $(map(total_degree, plist))") end - @debug "\tDegrees and lengths are $(map(first, pivots))" + debug_si("\tDegrees and lengths are $(map(first, pivots))") den_lcm = mapreduce( i -> funcs_den_nums[i][pivots_indices[i]], lcm, 1:length(funcs_den_nums), ) - @debug "Rational functions common denominator is of degree $(total_degree(den_lcm)) and of length $(length(den_lcm))" + debug_si("Rational functions common denominator is of degree $(total_degree(den_lcm)) and of length $(length(den_lcm))") is_constant(den_lcm) && - (@debug "Common denominator of the field generators is constant") + (debug_si("Common denominator of the field generators is constant")) existing_varnames = map(String, symbols(ring)) ystrs = ["y$i" for i in 1:length(existing_varnames)] @assert !(sat_varname in ystrs) "The name of the saturation variable collided with a primary variable" @@ -99,8 +99,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal length(ystrs) + 1 end varnames = append_at_index(ystrs, sat_var_index, sat_varname) - @debug "Saturating variable is $sat_varname, index is $sat_var_index" - flush(stdout) + debug_si("Saturating variable is $sat_varname, index is $sat_var_index") R_sat, v_sat = Nemo.PolynomialRing(K, varnames, ordering = ordering) # Saturation t_sat = v_sat[sat_var_index] @@ -146,8 +145,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal end end parent_ring_param, _ = PolynomialRing(ring, varnames, ordering = ordering) - @debug "Constructed MQS ideal in $R_sat with $(length(nums_qq) + 1) elements" - flush(stdout) + debug_si("Constructed MQS ideal in $R_sat with $(length(nums_qq) + 1) elements") @assert length(pivots_indices) == length(dens_indices) == length(dens_qq) @assert length(pivots_indices) == length(funcs_den_nums) @@ -200,12 +198,12 @@ function fractionfree_generators_raw(mqs::IdealMQS) old_varnames = map(i -> "y$i", 1:length(varnames)) new_varnames = map(i -> "라$i", 1:(length(varnames) + 1)) if !isempty(intersect(old_varnames, new_varnames)) - @warn "Intersection in two sets of variables!" varnames new_varnames + warn_si("Intersection in two sets of variables! $varnames $new_varnames") end # NOTE: new variables go first! big_ring, big_vars = PolynomialRing(K, vcat(new_varnames, old_varnames), ordering = :lex) - @info "" mqs.sat_var_index varnames ring_params parent(mqs.sat_qq) + info_si("$(mqs.sat_var_index) $(varnames) $ring_params $(parent(mqs.sat_qq))") nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq nums_y = map(num -> parent_ring_change(num, big_ring, matching = :byindex), nums_qq) dens_y = map(den -> parent_ring_change(den, big_ring, matching = :byindex), dens_qq) @@ -232,10 +230,10 @@ function ParamPunPam.reduce_mod_p!( mqs::IdealMQS, ff::Field, ) where {Field <: Union{Nemo.GaloisField, Nemo.GaloisFmpzField}} - @debug "Reducing MQS ideal modulo $(ff)" + debug_si("Reducing MQS ideal modulo $(ff)") # If there is a reduction modulo this field already, if haskey(mqs.cached_nums_gf, ff) - @debug "Cache hit with $(ff)!" + debug_si("Cache hit with $(ff)!") return nothing end nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq @@ -255,7 +253,7 @@ function ParamPunPam.specialize_mod_p( saturated = true, ) where {T <: Union{gfp_elem, gfp_fmpz_elem}} K_1 = parent(first(point)) - @debug "Evaluating MQS ideal over $K_1 at $point" + debug_si("Evaluating MQS ideal over $K_1 at $point") @assert haskey(mqs.cached_nums_gf, K_1) nums_gf, dens_gf, sat_gf = mqs.cached_nums_gf[K_1], mqs.cached_dens_gf[K_1], mqs.cached_sat_gf[K_1] @@ -286,7 +284,7 @@ function ParamPunPam.specialize_mod_p( end function specialize(mqs::IdealMQS, point::Vector{Nemo.fmpq}; saturated = true) - @debug "Evaluating MQS ideal over QQ at $point" + debug_si("Evaluating MQS ideal over QQ at $point") nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq dens_indices = mqs.dens_indices K = base_ring(mqs) diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index f1316ead7..8835e19d9 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -100,13 +100,11 @@ function field_contains( if isempty(ratfuncs) return Bool[] end - @debug "Finding pivot polynomials" - flush(stdout) + debug_si("Finding pivot polynomials") pivots = map(plist -> plist[findmin(map(total_degree, plist))[2]], ratfuncs) - @debug "\tDegrees are $(map(total_degree, pivots))" - flush(stdout) + debug_si("\tDegrees are $(map(total_degree, pivots))") - @debug "Estimating the sampling bound" + debug_si("Estimating the sampling bound") # uses Theorem 3.3 from https://arxiv.org/pdf/2111.00991.pdf # the comments below use the notation from the theorem ring = parent(first(first(ratfuncs))) @@ -123,14 +121,13 @@ function field_contains( extra_degree = total_degree(den_lcm) - total_degree(field.mqs.dens_qq[i]) degree = max(degree, extra_degree + maximum(total_degree, plist)) end - @debug "\tBound for the degrees is $degree" - flush(stdout) + debug_si("\tBound for the degrees is $degree") total_vars = foldl( union, map(plist -> foldl(union, map(poly -> Set(vars(poly)), plist)), field.dennums), ) - @debug "\tThe total number of variables in $(length(total_vars))" + debug_si("\tThe total number of variables in $(length(total_vars))") sampling_bound = BigInt( 3 * @@ -138,19 +135,16 @@ function field_contains( (length(ratfuncs)) * ceil(1 / (1 - p)), ) - @debug "\tSampling from $(-sampling_bound) to $(sampling_bound)" - flush(stdout) + debug_si("\tSampling from $(-sampling_bound) to $(sampling_bound)") mqs = field.mqs param_ring = ParamPunPam.parent_params(mqs) point = map(v -> Nemo.QQ(rand((-sampling_bound):sampling_bound)), gens(param_ring)) mqs_specialized = specialize(mqs, point) - @debug "Computing Groebner basis ($(length(mqs_specialized)) equations)" - flush(stdout) + debug_si("Computing Groebner basis ($(length(mqs_specialized)) equations)") mqs_ratfuncs = specialize(IdealMQS(ratfuncs), point; saturated = false) @assert parent(first(mqs_specialized)) == parent(first(mqs_ratfuncs)) - @debug "Starting the groebner basis computation" - flush(stdout) + debug_si("Starting the groebner basis computation") gb = groebner(mqs_specialized) result = map(iszero, normalform(gb, mqs_ratfuncs)) return result @@ -226,13 +220,13 @@ function beautifuly_generators( fracs = filter(!is_rational_func_const, fracs) fracs = unique(fracs) if isempty(fracs) - @debug "The set of generators is empty" + debug_si("The set of generators is empty") return fracs end # Remove redundant pass if discard_redundant sort!(fracs, lt = rational_function_cmp) - @debug "The pool of fractions:\n$(join(map(repr, fracs), ",\n"))" + debug_si("The pool of fractions:\n$(join(map(repr, fracs), ",\n"))") if reversed_order non_redundant = collect(1:length(fracs)) for i in length(fracs):-1:1 @@ -242,9 +236,9 @@ function beautifuly_generators( end result = check_field_membership_mod_p(fracs[setdiff(non_redundant, i)], [func]) - @debug "Simplification: inclusion check" func result + debug_si("Simplification: inclusion check $func $result") if result[1] - @debug "The function $func is discarded" + debug_si("The function $func is discarded") setdiff!(non_redundant, i) end end @@ -254,14 +248,14 @@ function beautifuly_generators( for i in 2:length(fracs) func = fracs[i] result = check_field_membership_mod_p(fracs[non_redundant], [func]) - @debug "Simplification: inclusion check" func result + debug_si("Simplification: inclusion check $func $result") if !result[1] - @debug "The function $func is included in the set of generators" + debug_si("The function $func is included in the set of generators") push!(non_redundant, i) end end end - @debug "Out of $(length(fracs)) simplified generators there are $(length(non_redundant)) non redundant" + debug_si("Out of $(length(fracs)) simplified generators there are $(length(non_redundant)) non redundant") fracs = fracs[non_redundant] end sort!(fracs, lt = rational_function_cmp) @@ -321,7 +315,7 @@ function groebner_basis_coeffs( gb, fracs, new_rff = nothing, nothing, nothing # Check if the basis is in cache if haskey(mqs.cached_groebner_bases, (ordering, up_to_degree)) - @debug "Cache hit with ($ordering, $up_to_degree)!" + debug_si("Cache hit with ($ordering, $up_to_degree)!") gb = mqs.cached_groebner_bases[ordering, up_to_degree] basis_coeffs = map(collect ∘ coefficients, gb) fracs = collect(mapreduce(Set, union!, basis_coeffs)) @@ -331,7 +325,7 @@ function groebner_basis_coeffs( current_degrees = (2, 2) two_sided_inclusion = false while !two_sided_inclusion && all(current_degrees .<= up_to_degree) - @debug "Computing GB with parameters up to degrees $(current_degrees)" + debug_si("Computing GB with parameters up to degrees $(current_degrees)") runtime = @elapsed gb = ParamPunPam.paramgb( mqs, up_to_degree = current_degrees, @@ -343,12 +337,12 @@ function groebner_basis_coeffs( _runtime_logger[:id_npoints_interpolation] += ParamPunPam._runtime_data[:npoints_interpolation] _runtime_logger[:id_groebner_time] += runtime - @debug "Groebner basis computed in $runtime seconds" + debug_si("Groebner basis computed in $runtime seconds") basis_coeffs = map(collect ∘ coefficients, gb) basis_coeffs_set = mapreduce(Set, union!, basis_coeffs) fracs = collect(basis_coeffs_set) - @debug "Generators up to degrees $(current_degrees) are" fracs - @debug "Checking two-sided inclusion modulo a prime" + debug_si("Generators up to degrees $(current_degrees) are $fracs") + debug_si("Checking two-sided inclusion modulo a prime") time_start = time_ns() # Check inclusion: in new_rff = RationalFunctionField(fracs) @@ -359,10 +353,10 @@ function groebner_basis_coeffs( runtime = (time_ns() - time_start) / 1e9 _runtime_logger[:id_inclusion_check_mod_p] += runtime two_sided_inclusion = two_sided_inclusion && all(inclusion) - @debug "Inclusion checked in $(runtime) seconds. Result: $two_sided_inclusion" + debug_si("Inclusion checked in $(runtime) seconds. Result: $two_sided_inclusion") current_degrees = current_degrees .* 2 end - @debug "The coefficients of the Groebner basis are presented by $(length(fracs)) rational functions" + debug_si("The coefficients of the Groebner basis are presented by $(length(fracs)) rational functions") new_rff.mqs.cached_groebner_bases[ordering, up_to_degree] = gb rff.mqs.cached_groebner_bases[ordering, up_to_degree] = gb return new_rff @@ -388,7 +382,7 @@ function generating_sets_fan( time_start = time_ns() vars = gens(parent(rff.mqs)) nbases = length(vars) - @info "Computing $nbases Groebner bases for degrees $up_to_degree for block orderings" + info_si("Computing $nbases Groebner bases for degrees $up_to_degree for block orderings") ordering_to_generators = Dict() if code == 0 return ordering_to_generators @@ -414,7 +408,7 @@ function generating_sets_fan( # n1, n2 = div(n, 2), n - div(n, 2) n1, n2 = n - 1, 1 ord = DegRevLex(vars_shuffled[1:n1]) * DegRevLex(vars_shuffled[(n1 + 1):end]) - @debug "Computing GB for ordering" ord + debug_si("Computing GB for ordering $ord") new_rff = groebner_basis_coeffs( gb_rff, seed = seed, @@ -431,7 +425,7 @@ function generating_sets_fan( n = length(vars_shuffled) n1, n2 = max(n - 2, 1), min(2, length(vars) - 1) ord = DegRevLex(vars_shuffled[1:n1]) * DegRevLex(vars_shuffled[(n1 + 1):end]) - @debug "Computing GB for ordering" ord + debug_si("Computing GB for ordering $ord") new_rff = groebner_basis_coeffs( gb_rff, seed = seed, @@ -449,7 +443,7 @@ function generating_sets_fan( n1 = div(n, 2) n2 = n - n1 ord = DegRevLex(vars_shuffled[1:n1]) * DegRevLex(vars_shuffled[(n1 + 1):end]) - @debug "Computing GB for ordering" ord + debug_si("Computing GB for ordering $ord") new_rff = groebner_basis_coeffs( gb_rff, seed = seed, @@ -461,7 +455,7 @@ function generating_sets_fan( end end _runtime_logger[:id_gbfan_time] = (time_ns() - time_start) / 1e9 - @info "Computed Groebner bases in $((time_ns() - time_start) / 1e9) seconds" + info_si("Computed Groebner bases in $((time_ns() - time_start) / 1e9) seconds") return ordering_to_generators end @@ -494,7 +488,7 @@ function simplified_generating_set( check_variables = false, # almost always slows down and thus turned off rational_interpolator = :VanDerHoevenLecerf, ) - @info "Simplifying generating set. Simplification level: $simplify" + info_si("Simplifying generating set. Simplification level: $simplify") _runtime_logger[:id_groebner_time] = 0.0 _runtime_logger[:id_calls_to_gb] = 0 _runtime_logger[:id_inclusion_check_mod_p] = 0.0 @@ -558,23 +552,23 @@ function simplified_generating_set( append!(new_fracs, generators) end new_fracs_unique = unique(new_fracs) - @debug """ + debug_si(""" Final cleaning and simplification of generators. - Out of $(length(new_fracs)) fractions $(length(new_fracs_unique)) are syntactically unique.""" + 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" + debug_si("Checking inclusion with probability $p") runtime = @elapsed result = issubfield(rff, RationalFunctionField(new_fracs), p) _runtime_logger[:id_inclusion_check] = runtime if !result - @warn "Field membership check failed. Error will follow." + warn_si("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" - @debug "Out of $(length(rff.mqs.nums_qq)) initial generators there are $(length(new_fracs)) indepdendent" + info_si("Inclusion checked with probability $p in $(_runtime_logger[:id_inclusion_check]) seconds") + debug_si("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 - @debug "The ranking of the new set of generators is $ranking" + debug_si("The ranking of the new set of generators is $ranking") return new_fracs end diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index 02aafcbbd..c7f2c7c9a 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -63,27 +63,27 @@ function local_normal_forms( monoms_ff = Vector{elem_type(ring_ff)}(undef, 0) xs_ff = cut_at_index(xs_ff, mqs.sat_var_index) pivot_vectors = map(f -> exponent_vector(f, 1), xs_ff) - @debug """ + debug_si(""" variables in the finite field: $(xs_ff) gb parent: $(ring_ff) specialized gb: $(gb_ff_spec) - Evaluation point: $point_ff_ext""" + Evaluation point: $point_ff_ext""") # Compute the normal forms of all monomials of degrees up to `up_to_degree` for deg in 1:up_to_degree for combination in Combinatorics.with_replacement_combinations(pivot_vectors, deg) exp_vect = sum(combination) if in(stop_vectors, exp_vect) - @debug "Skipping exponent vector $exp_vect" + debug_si("Skipping exponent vector $exp_vect") continue end monom_ff = ring_ff([one(finite_field)], [exp_vect]) monom_ff_spec = evaluate(monom_ff, point_ff_ext) monom_mqs_ff_spec = monom_ff - monom_ff_spec divisors_ff, nf_ff = divrem(monom_mqs_ff_spec, gb_ff_spec) - @debug """ + debug_si(""" The normal form of $monom_mqs_ff_spec is: normalform = $nf_ff - divisors = $divisors_ff""" + divisors = $divisors_ff""") push!(monoms_ff, monom_ff) push!(normal_forms_ff, nf_ff) end @@ -161,7 +161,7 @@ function linear_relations_over_a_field(polys, preimages) for ind in zero_inds push!(relations, preimages[ind]) end - @debug "Zeroed polynomials are" preimages[zero_inds] + debug_si("Zeroed polynomials are $(preimages[zero_inds])") permutation = setdiff(permutation, zero_inds) # Sort, the first monom is the smallest lead_monoms = map(f -> iszero(f) ? one(f) : leading_monomial(f), polys) @@ -303,9 +303,9 @@ function linear_relations_between_normal_forms( nparams = nvars(ring_param) finite_field = Nemo.GF(2^30 + 3) ParamPunPam.reduce_mod_p!(mqs, finite_field) - @info "Computing normal forms of degree $up_to_degree in $nparams variables" - @debug """Variables ($nparams in total): $xs_param - Modulo: $finite_field""" + info_si("Computing normal forms of degree $up_to_degree in $nparams variables") + debug_si("""Variables ($nparams in total): $xs_param + Modulo: $finite_field""") # We first compute relations between the normal forms of linear monomials. # Then, we use this knowledge to drop out some monomials of higher degrees. tref = TinyRowEchelonForm{Int}() @@ -315,7 +315,7 @@ function linear_relations_between_normal_forms( for i in 1:length(normal_forms_ff_1) !iszero(normal_forms_ff_1[i]) && continue !(length(monoms_ff_1[i]) == 1) && continue - @debug "Registering existing monomial $(monoms_ff_1[i])" + debug_si("Registering existing monomial $(monoms_ff_1[i])") push!(relations_ff_1, monoms_ff_1[i]) push!(tref, exponent_vector(monoms_ff_1[i], 1)) end @@ -325,17 +325,17 @@ function linear_relations_between_normal_forms( while true iters += 1 point = ParamPunPam.distinct_nonzero_points(finite_field, nvars(ring_param)) - @debug "Used specialization points: $iters" - @debug "Computing normal forms to to degree $up_to_degree" + debug_si("Used specialization points: $iters") + debug_si("Computing normal forms to to degree $up_to_degree") gb_ff, normal_forms_ff, monoms_ff = local_normal_forms(mqs, finite_field, up_to_degree, point, stop_vectors = tref) if isempty(normal_forms_ff) break end - @debug "Computing relations of $(length(normal_forms_ff)) normal forms" + debug_si("Computing relations of $(length(normal_forms_ff)) normal forms") relations_ff, normal_forms_ff, monoms_ff = linear_relations_over_a_field(normal_forms_ff, monoms_ff) - @debug "Obtained $(length(relations_ff)) local relations" + debug_si("Obtained $(length(relations_ff)) local relations") if iters == 1 # first point in the sequence, take all relations complete_intersection_relations_ff = relations_ff @@ -365,9 +365,9 @@ function linear_relations_between_normal_forms( push!(zeroed_relations_inds, i) end end - @debug """ + debug_si(""" Relations in the previous intersection: $(length(complete_intersection_relations_ff)) - Vanished at the current point: $(length(zeroed_relations_inds))""" + Vanished at the current point: $(length(zeroed_relations_inds))""") non_zeroed_relations_inds = setdiff(collect(1:n_relations_ff), zeroed_relations_inds) zeroed_relations_from_complete_intersection = @@ -383,14 +383,14 @@ function linear_relations_between_normal_forms( non_zeroed_relations_from_complete_intersection, zeroed_relations_from_complete_intersection, ) - @debug "There are $(length(complete_intersection_relations_ff)) relations in the intersection" + debug_si("There are $(length(complete_intersection_relations_ff)) relations in the intersection") m_relations_ff = length(complete_intersection_relations_ff) if n_relations_ff == m_relations_ff || isempty(complete_intersection_relations_ff) break end end union!(complete_intersection_relations_ff, relations_ff_1) - @debug "Reconstructing relations to rationals" + debug_si("Reconstructing relations to rationals") relations_qq = Vector{Generic.Frac{elem_type(ring_param)}}( undef, length(complete_intersection_relations_ff), @@ -400,10 +400,10 @@ function linear_relations_between_normal_forms( success, relation_qq = ParamPunPam.rational_reconstruct_polynomial(ring, relation_ff) if !success - @warn """ + warn_si(""" Failed to reconstruct the $i-th relation. Error will follow. relation: $relation_ff - modulo: $finite_field""" + modulo: $finite_field""") throw(ErrorException("Rational reconstruction failed.")) end relation_qq_param = evaluate( @@ -412,6 +412,6 @@ function linear_relations_between_normal_forms( ) relations_qq[i] = relation_qq_param // one(relation_qq_param) end - @info "Used $iters specializations in $((time_ns() - time_start) / 1e9) seconds, found $(length(complete_intersection_relations_ff)) relations" + info_si("Used $iters specializations in $((time_ns() - time_start) / 1e9) seconds, found $(length(complete_intersection_relations_ff)) relations") relations_qq end diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 75a9417db..a36e0361f 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -90,7 +90,7 @@ function assess_identifiability( ode::ODE{P}; funcs_to_check = Vector(), p::Float64 = 0.99, - loglevel = Logging.Info + loglevel = Logging.Info, ) where {P <: MPolyElem{fmpq}} restart_logging(loglevel = loglevel) p_glob = 1 - (1 - p) * 0.9 @@ -108,6 +108,7 @@ function assess_identifiability( p = p_loc, type = :SE, trbasis = trbasis, + loglevel = loglevel, ) info_si("Local identifiability assessed in $runtime seconds") debug_si("Trasncendence basis to be specialized is $trbasis") @@ -162,7 +163,7 @@ function assess_identifiability( measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = [], p = 0.99, - loglevel = Logging.Info + loglevel = Logging.Info, ) if isempty(measured_quantities) measured_quantities = get_measured_quantities(ode) diff --git a/src/discrete.jl b/src/discrete.jl index 628797e98..1c29e43b0 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -72,12 +72,12 @@ function differentiate_sequence_solution( inputs::Dict{P, Array{T, 1}}, num_terms::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - @debug "Computing the power series solution of the system" + debug_si("Computing the power series solution of the system") seq_sol = sequence_solution(dds, params, ic, inputs, num_terms) generalized_params = vcat(dds.x_vars, dds.parameters) bring = base_ring(dds.poly_ring) - @debug "Solving the variational system at the solution" + debug_si("Solving the variational system at the solution") part_diffs = Dict( (x, p) => derivative(dds.x_equations[x], p) for x in dds.x_vars for p in generalized_params @@ -126,10 +126,10 @@ function differentiate_sequence_output( inputs::Dict{P, Array{T, 1}}, num_terms::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - @debug "Computing partial derivatives of the solution" + debug_si("Computing partial derivatives of the solution") seq_sol, sol_diff = differentiate_sequence_solution(dds, params, ic, inputs, num_terms) - @debug "Evaluating the partial derivatives of the outputs" + debug_si("Evaluating the partial derivatives of the outputs") generalized_params = vcat(dds.x_vars, dds.parameters) part_diffs = Dict( (y, p) => derivative(dds.y_equations[y], p) for y in dds.y_vars for @@ -195,7 +195,7 @@ function _assess_local_identifiability_discrete( ) where {P <: MPolyElem{Nemo.fmpq}} bring = base_ring(dds.poly_ring) - @debug "Extending the model" + debug_si("Extending the model") dds_ext = add_outputs(dds, Dict("loc_aux_$i" => f for (i, f) in enumerate(funcs_to_check))) @@ -206,7 +206,7 @@ function _assess_local_identifiability_discrete( known_ic = dds_ext.x_vars end - @debug "Computing the observability matrix" + debug_si("Computing the observability matrix") prec = length(dds.x_vars) + length(dds.parameters) # Computing the bound from the Schwartz-Zippel-DeMilo-Lipton lemma @@ -221,7 +221,7 @@ function _assess_local_identifiability_discrete( Jac_degree += 2 * deg_y * prec end D = Int(ceil(Jac_degree * length(funcs_to_check) / (1 - p))) - @debug "Sampling range $D" + debug_si("Sampling range $D") # Parameter values are the same across all the replicas params_vals = Dict(p => bring(rand(1:D)) for p in dds_ext.parameters) @@ -231,11 +231,11 @@ function _assess_local_identifiability_discrete( u => [bring(rand(1:D)) for i in 1:prec] for u in dds_ext.u_vars ) - @debug "Computing the output derivatives" + debug_si("Computing the output derivatives") output_derivatives = differentiate_sequence_output(dds_ext, params_vals, ic, inputs, prec) - @debug "Building the matrices" + debug_si("Building the matrices") Jac = zero( Nemo.MatrixSpace( bring, @@ -265,7 +265,7 @@ function _assess_local_identifiability_discrete( end end - @debug "Computing the result" + debug_si("Computing the result") base_rank = LinearAlgebra.rank(Jac) result = Dict{Any, Bool}() for i in 1:length(funcs_to_check) @@ -314,7 +314,7 @@ function assess_local_identifiability( ) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(dds)) - @info "Measured quantities are not provided, trying to find the outputs in input dynamical system." + info_si("Measured quantities are not provided, trying to find the outputs in input dynamical system.") measured_quantities = filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(dds), @@ -342,7 +342,7 @@ function assess_local_identifiability( nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) if length(known_ic) > 0 - @warn "Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !" + warn_si("Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !") end return out_dict end diff --git a/src/elimination.jl b/src/elimination.jl index c36b2db63..3c68b735a 100644 --- a/src/elimination.jl +++ b/src/elimination.jl @@ -32,8 +32,7 @@ function det_minor_expansion_inner( cache[discarded] = result end if length(discarded[1]) < 3 - @debug "Discarded: $discarded" - flush(stdout) + debug_si("Discarded: $discarded") end return result end @@ -202,12 +201,11 @@ function Base.iterate( i::Int = 1, ) where {P <: MPolyElem{<:FieldElem}} if i > length(gpg.cached_points) - @debug "Generating new point on the variety" + debug_si("Generating new point on the variety") sample_max = i * 50 result = undef while true - @debug "Preparing initial condition" - flush(stdout) + debug_si("Preparing initial condition") base_field = base_ring(gpg.big_ring) param_values = Dict{P, Int}(p => rand(1:sample_max) for p in gpg.ode.parameters) initial_conditions = @@ -216,8 +214,7 @@ function Base.iterate( u => [rand(1:sample_max) for _ in 1:(gpg.precision)] for u in gpg.ode.u_vars ) - @debug "Computing a power series solution" - flush(stdout) + debug_si("Computing a power series solution") ps_solution = undef try ps_solution = power_series_solution( @@ -228,13 +225,11 @@ function Base.iterate( gpg.precision, ) catch e - @debug "$e" - flush(stdout) + debug_si("$e") continue end - @debug "Constructing the point" - flush(stdout) + debug_si("Constructing the point") result = Dict{P, gpg.number_type}( switch_ring(p, gpg.big_ring) => base_field(c) for (p, c) in param_values ) @@ -283,7 +278,6 @@ function choose( # get accounts for the fact that the big ring may contain some auxiliary variables, e.g. rand_proj_var point = [get(p, v, zero(base_ring(parent(polys[1])))) for v in vars] polys = filter(e -> (evaluate(e, point) == 0), polys) - flush(stdout) if length(polys) <= 1 break end @@ -322,14 +316,12 @@ function eliminate_var( lg = coeff(g, [var_elim], [Nemo.degree(g, var_elim)]) (flag, q) = divides(lg, lf) if flag - @debug "\t Decreasing degree with linear combination" - flush(stdout) + debug_si("\t Decreasing degree with linear combination") g = g - q * f * var_elim^(Nemo.degree(g, var_elim) - Nemo.degree(f, var_elim)) elseif (Nemo.degree(g, var_elim) == Nemo.degree(f, var_elim)) (flag, q) = divides(lf, lg) if flag - @debug "\t Decreasing degree with linear combination" - flush(stdout) + debug_si("\t Decreasing degree with linear combination") f = f - q * g else break @@ -376,32 +368,28 @@ function eliminate_var( resultant = f else if Nemo.degree(f, var_elim) > 1 - @debug "Calculating Bezout Matrix" - flush(stdout) + debug_si("Calculating the Bezout Matrix") M = Bezout_matrix(f, g, var_elim) else - @debug "Calculating Sylvester matrix" - flush(stdout) + debug_si("Calculating the Sylvester matrix") M = Sylvester_matrix(f, g, var_elim) end - @debug "Simplifying the matrix" - flush(stdout) + debug_si("Simplifying the matrix") M_simp, matrix_factors = simplify_matrix(M) - @debug "Removed factors $(map(length, matrix_factors))" + debug_si("Removed factors $(map(length, matrix_factors))") M_size = zero(Nemo.MatrixSpace(Nemo.ZZ, ncols(M_simp), ncols(M_simp))) for i in 1:ncols(M_simp) for j in 1:ncols(M_simp) M_size[i, j] = length(M_simp[i, j]) end end - @debug "\t Matrix size: \n $M_size" - @debug "\t Computing determinant" - flush(stdout) + debug_si("\t Matrix size: \n $M_size") + debug_si("\t Computing determinant") resultant = det_minor_expansion(M_simp) end - @debug "Degrees are", [(v, Nemo.degree(resultant, v)) for v in vars(resultant)] + debug_si("Degrees are $([(v, Nemo.degree(resultant, v)) for v in vars(resultant)])") # Step 4: Eliminate extra factors - @debug "Eliminating extra factors" + debug_si("Eliminating extra factors") factors = fast_factor(resultant) for mfac in matrix_factors for fac in fast_factor(mfac) @@ -413,9 +401,8 @@ function eliminate_var( res = choose(factors, generic_point_generator) for f in factors if f != res - @debug "\t \t Size of extra factor: $(length(f))" - @debug "\t \t It is $f" - flush(stdout) + debug_si("\t \t Size of extra factor: $(length(f))") + debug_si("\t \t It is $f") end end return res diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index f33a57f20..d18fc2983 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -27,15 +27,14 @@ function extract_identifiable_functions_raw( bring, _ = Nemo.PolynomialRing(base_ring(ode.poly_ring), varnames) if with_states - @debug "Computing Lie derivatives" + debug_si("Computing Lie derivatives") for f in states_generators(ode, io_equations) num, den = unpack_fraction(parent_ring_change(f, bring)) push!(coeff_lists[:with_states], [den, num]) end end - @debug "Extracting coefficients" - flush(stdout) + debug_si("Extracting coefficients") if !isempty(ode.parameters) nonparameters = filter( v -> !(var_to_str(v) in map(var_to_str, ode.parameters)), @@ -57,7 +56,7 @@ function extract_identifiable_functions_raw( push!(coeff_lists[:no_states], as_list) end else - @debug "Known quantity $p cannot be casted and is thus dropped" + debug_si("Known quantity $p cannot be casted and is thus dropped") end end @@ -106,32 +105,32 @@ function initial_identifiable_functions( var_change_policy = :default, rational_interpolator = :VanDerHoevenLecerf, ) where {T} - @info "Computing IO-equations" + info_si("Computing IO-equations") ioeq_time = @elapsed io_equations = find_ioequations(ode; var_change_policy = var_change_policy) - @debug "Sizes: $(map(length, values(io_equations)))" - @info "Computed IO-equations in $ioeq_time seconds" + debug_si("Sizes: $(map(length, values(io_equations)))") + info_si("Computed IO-equations in $ioeq_time seconds") _runtime_logger[:ioeq_time] = ioeq_time if isempty(ode.parameters) - @info "No parameters, so Wronskian computation is not needed" + info_si("No parameters, so Wronskian computation is not needed") else - @info "Computing Wronskians" + info_si("Computing Wronskians") wrnsk_time = @elapsed wrnsk = wronskian(io_equations, ode) - @info "Computed Wronskians in $wrnsk_time seconds" + info_si("Computed Wronskians in $wrnsk_time seconds") _runtime_logger[:wrnsk_time] = wrnsk_time dims = map(ncols, wrnsk) - @info "Dimensions of the Wronskians $dims" + info_si("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" + debug_si("Dimensions of the Wronskians $dims") + debug_si("Ranks of the Wronskians $wranks") + info_si("Ranks of the Wronskians computed in $rank_times seconds") _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)" + warn_si("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) or transforming all the parameters to states with zero derivative") end end @@ -143,7 +142,7 @@ function initial_identifiable_functions( ) if with_states && !isempty(ode.parameters) - @debug "Generators of identifiable functions involve states, the parameter-only part is getting simplified" + debug_si("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( base_ring(bring), @@ -200,7 +199,7 @@ function check_identifiability( for f in funcs_to_check num, den = unpack_fraction(f) if !all(v -> v in ode.parameters, union(vars(num), vars(den))) - @info "Functions to check involve states" + info_si("Functions to check involve states") states_needed = true break end @@ -306,7 +305,7 @@ function assess_global_identifiability( ) where {P <: MPolyElem{fmpq}} submodels = find_submodels(ode) if length(submodels) > 0 - @info "Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`" + info_si("Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`") end result = check_identifiability( diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index feb404e42..43a66090b 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -51,15 +51,16 @@ function find_identifiable_functions( rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, ) where {T <: MPolyElem{fmpq}} + restart_logging(loglevel=loglevel) Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) runtime_start = time_ns() if isempty(ode.parameters) && !with_states - @warn """ + warn_si(""" There are no parameters in the given ODE, thus no identifiabile functions. Use `find_identifiable_functions` with keyword `with_states=true` to - compute the functions with the ODE states included.""" + compute the functions with the ODE states included.""") bring = parent(ode) id_funcs = [one(bring)] return id_funcs @@ -90,7 +91,7 @@ function find_identifiable_functions( id_funcs_fracs = [parent_ring_change(f, parent(ode)) for f in id_funcs_fracs] _runtime_logger[:id_total] = (time_ns() - runtime_start) / 1e9 _runtime_logger[:are_id_funcs_polynomial] = all(isone ∘ denominator, id_funcs_fracs) - @info "The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds" + info_si("The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds") return id_funcs_fracs end @@ -104,6 +105,8 @@ system. This functions takes the following optional arguments: - `measured_quantities` - the output functions of the model. +- `loglevel` - the verbosity of the logging + (can be Logging.Error, Logging.Warn, Logging.Info, Logging.Debug) ## Example diff --git a/src/io_equation.jl b/src/io_equation.jl index 19f8cd87e..92a77a9dc 100644 --- a/src/io_equation.jl +++ b/src/io_equation.jl @@ -146,7 +146,7 @@ function find_ioprojections( coef * evaluate(y_eq, [y], [zero(ring)]) // derivative(y_eq, y) end projection_equation, _ = unpack_fraction(projection_equation) - @debug "Extra projection equation $projection_equation" + debug_si("Extra projection equation $projection_equation") end while true @@ -157,9 +157,9 @@ function find_ioprojections( if isempty(var_degs) break end - @debug "Current degrees of io-equations $var_degs" - @debug "Orders: $y_orders" - @debug "Sizes: $(Dict(y => length(eq) for (y, eq) in y_equations))" + debug_si("Current degrees of io-equations $var_degs") + debug_si("Orders: $y_orders") + debug_si("Sizes: $(Dict(y => length(eq) for (y, eq) in y_equations))") # choosing the output to prolong outputs_with_scores = [ @@ -171,26 +171,22 @@ function find_ioprojections( d[1], ) for d in var_degs ] - @debug "Scores: $outputs_with_scores" + debug_si("Scores: $outputs_with_scores") y_prolong = sort(outputs_with_scores)[1][end] y_orders[y_prolong] += 1 - @debug "Prolonging output $y_prolong" - flush(stdout) + debug_si("Prolonging output $y_prolong") # Calculate the Lie derivative of the io_relation - @debug "Prolonging" - flush(stdout) + debug_si("Prolonging") next_y_equation = diff_poly(y_equations[y_prolong], derivation) for (x, eq) in x_equations - @debug "Eliminating the derivative of $x" - flush(stdout) + debug_si("Eliminating the derivative of $x") next_y_equation = eliminate_var(eq, next_y_equation, derivation[x], point_generator) end for (y, eq) in y_equations if y != y_prolong - @debug "Eliminating the leader of the equation for $y" - flush(stdout) + debug_si("Eliminating the leader of the equation for $y") # an ugly way of gettin the leader, to replace next_y_equation = eliminate_var( next_y_equation, @@ -208,8 +204,7 @@ function find_ioprojections( our_choice = sort(var_degs_next)[1] var_elim_deg, var_elim = our_choice[1], our_choice[3] - @debug "Elimination of $var_elim, $(length(x_equations)) left" - flush(stdout) + debug_si("Elimination of $var_elim, $(length(x_equations)) left") # Possible variable change for Axy + Bx + p(y) (x = var_elim) if auto_var_change && (var_elim_deg == 1) @@ -220,8 +215,7 @@ function find_ioprojections( A, B = simplify_frac(A, B) if isempty(filter!(v -> (v in keys(x_equations)), vars(A))) && (B != 0) # && (length(coeffs(A))==1) # variable change x_i' --> x_i' - (B/A)', x_i --> x_i - B/A - @debug "\t Applying variable change: $(x) --> $(x) - ( $B )/( $A )" - flush(stdout) + debug_si("\t Applying variable change: $(x) --> $(x) - ( $B )/( $A )") dB = diff_poly(B, derivation) dA = diff_poly(A, derivation) numer_d, denom_d = simplify_frac(A * dB - dA * B, A * A) @@ -236,8 +230,7 @@ function find_ioprojections( generator_var_change(point_generator, x, A * x + B, A) # Change current system - @debug "Change in the system" - flush(stdout) + debug_si("Change in the system") x_equations[x] = make_substitution( x_equations[x], derivation[x], @@ -245,29 +238,23 @@ function find_ioprojections( denom_d, ) for xx in keys(x_equations) - @debug "\t Change in the equation for $xx" - flush(stdout) + debug_si("\t Change in the equation for $xx") x_equations[xx] = make_substitution(x_equations[xx], x, A * x - B, A) end - @debug "Change in the outputs" - flush(stdout) + debug_si("Change in the outputs") for y in keys(y_equations) - @debug "\t Change in the output $y" - flush(stdout) + debug_si("\t Change in the output $y") y_equations[y] = make_substitution(y_equations[y], x, A * x - B, A) end - @debug "\t Change in the prolonged equation" - flush(stdout) + debug_si("\t Change in the prolonged equation") next_y_equation = make_substitution(next_y_equation, x, A * x - B, A) # recalibrate system - @debug "Unmixing the derivatives" - flush(stdout) + debug_si("Unmixing the derivatives") for xx in setdiff(keys(x_equations), [x]) - @debug "\t Unmixing $xx" - flush(stdout) + debug_si("\t Unmixing $xx") x_equations[x] = eliminate_var( x_equations[x], x_equations[xx], @@ -276,11 +263,10 @@ function find_ioprojections( ) end # change the projection - @debug "Change of variables in the extra projection" + debug_si("Change of variables in the extra projection") projection_equation = make_substitution(projection_equation, x, A * x - B, A) - @debug "Change of variables performed" - flush(stdout) + debug_si("Change of variables performed") break end end @@ -289,20 +275,16 @@ function find_ioprojections( # Eliminate var_elim from the system delete!(x_equations, var_elim) - @debug "Elimination in states" - flush(stdout) + debug_si("Elimination in states") for (x, eq) in x_equations - @debug "\t Elimination in the equation for $x" - flush(stdout) + debug_si("\t Elimination in the equation for $x") x_equations[x] = eliminate_var(eq, y_equations[y_prolong], var_elim, point_generator) end - @debug "Elimination in y_equations" - flush(stdout) + debug_si("Elimination in y_equations") for y in keys(y_equations) if y != y_prolong - @debug "Elimination in the output $y" - flush(stdout) + debug_si("Elimination in the output $y") y_equations[y] = eliminate_var( y_equations[y], y_equations[y_prolong], @@ -311,22 +293,20 @@ function find_ioprojections( ) end end - @debug "\t Elimination in the extra projection" + debug_si("\t Elimination in the extra projection") projection_equation = eliminate_var( projection_equation, y_equations[y_prolong], var_elim, point_generator, ) - @debug "\t Elimination in the prolonged equation" - flush(stdout) + debug_si("\t Elimination in the prolonged equation") y_equations[y_prolong] = eliminate_var( y_equations[y_prolong], next_y_equation, var_elim, point_generator, ) - flush(stdout) end io_projections = Dict( @@ -350,6 +330,7 @@ Finds the input-output equations of an ODE system Input: - `ode` - the ODE system - `var_change_policy` - whether to perform automatic variable change, can be one of `:default`, `:yes`, `:no` +- `loglevel` - logging level (default: Logging.Info) Output: - a dictionary from “leaders” to the corresponding input-output equations; if an extra projection is needed, @@ -358,6 +339,7 @@ Output: function find_ioequations( ode::ODE{P}; var_change_policy = :default, + loglevel = Logging.Info ) where {P <: MPolyElem{<:FieldElem}} # Setting the var_change policy if (var_change_policy == :yes) || @@ -367,29 +349,29 @@ function find_ioequations( (var_change_policy == :default && length(ode.y_vars) >= 3) auto_var_change = false else - @error "Unknown var_change policy $var_change_policy" + error_si("Unknown var_change policy $var_change_policy") return end io_projections, _, _ = find_ioprojections(ode, auto_var_change, nothing) ring = parent(first(values(io_projections))) - @debug "Check whether the original projections are enough" + debug_si("Check whether the original projections are enough") if length(io_projections) == 1 || check_primality(io_projections) - @debug "The projections generate an ideal with a single components of highest dimension, returning" + debug_si("The projections generate an ideal with a single components of highest dimension, returning") return io_projections end sampling_range = 5 while true - @debug "There are several components of the highest dimension, trying to isolate one" + debug_si("There are several components of the highest dimension, trying to isolate one") extra_projection = sum(rand(1:sampling_range) * v for v in keys(io_projections)) - @debug "Extra projections: $extra_projection" + debug_si("Extra projections: $extra_projection") new_projections, _, projection_equation = find_ioprojections(ode, auto_var_change, extra_projection) - @debug "Check primality" + debug_si("Check primality") if check_primality(io_projections, [projection_equation]) - @debug "Single component of highest dimension isolated, returning" + debug_si("Single component of highest dimension isolated, returning") io_projections[str_to_var(PROJECTION_VARNAME, parent(projection_equation))] = projection_equation break diff --git a/src/lincomp.jl b/src/lincomp.jl index 72dd6a327..3e394b46e 100644 --- a/src/lincomp.jl +++ b/src/lincomp.jl @@ -29,8 +29,10 @@ function linear_compartment_model( graph::Vector{Vector{Int}}, inputs::Vector{Int}, outputs::Vector{Int}, - leaks::Vector{Int}, + leaks::Vector{Int}; + loglevel = Logging.Info, ) + restart_logging(loglevel = loglevel) n = length(graph) x_vars_names = ["x$i" for i in 1:n] y_vars_names = ["y$i" for i in outputs] diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index c99baf6ba..e32e0abfb 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -31,14 +31,14 @@ function differentiate_solution( inputs::Dict{P, Array{T, 1}}, prec::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - @debug "Computing the power series solution of the system" + debug_si("Computing the power series solution of the system") ps_sol = power_series_solution(ode, params, ic, inputs, prec) ps_ring = parent(first(values(ps_sol))) for p in ode.parameters ps_sol[p] = ps_ring(params[p]) end - @debug "Building the variational system at the solution" + debug_si("Building the variational system at the solution") # Y' = AY + B vars = vcat(ode.x_vars, ode.parameters) SA = AbstractAlgebra.MatrixSpace(ps_ring, length(ode.x_vars), length(ode.x_vars)) @@ -60,7 +60,7 @@ function differentiate_solution( initial_condition[i, i] = 1 end - @debug "Solving the variational system and forming the output" + debug_si("Solving the variational system and forming the output") sol_var_system = ps_matrix_linear_de(A, B, initial_condition, prec) return ( ps_sol, @@ -87,14 +87,14 @@ function differentiate_output( inputs::Dict{P, Array{T, 1}}, prec::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - @debug "Computing partial derivatives of the solution" + debug_si("Computing partial derivatives of the solution") ps_sol, sol_diff = differentiate_solution(ode, params, ic, inputs, prec) ps_ring = parent(first(values(ps_sol))) for p in ode.parameters ps_sol[p] = ps_ring(params[p]) end - @debug "Evaluating the partial derivatives of the outputs" + debug_si("Evaluating the partial derivatives of the outputs") result = Dict() for (y, g) in ode.y_equations result[y] = Dict() @@ -169,10 +169,12 @@ function assess_local_identifiability( funcs_to_check = Array{}[], p::Float64 = 0.99, type = :SE, + loglevel = Logging.Info, ) + restart_logging(logelevel = loglevel) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) - @info "Measured quantities are not provided, trying to find the outputs in input ODE." + info_si("Measured quantities are not provided, trying to find the outputs in input ODE.") measured_quantities = filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(ode), @@ -200,6 +202,7 @@ function assess_local_identifiability( funcs_to_check = funcs_to_check_, p = p, type = type, + loglevel = loglevel ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) @@ -210,6 +213,7 @@ function assess_local_identifiability( funcs_to_check = funcs_to_check_, p = p, type = type, + loglevel = loglevel, ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) @@ -234,7 +238,9 @@ function assess_local_identifiability( p::Float64 = 0.99, type = :SE, trbasis = nothing, + loglevel = Logging.Info ) where {P <: MPolyElem{Nemo.fmpq}} + restart_logging(loglevel = loglevel) if isempty(funcs_to_check) funcs_to_check = ode.parameters if type == :SE @@ -248,7 +254,7 @@ function assess_local_identifiability( num, den = unpack_fraction(f) for v in vcat(vars(num), vars(den)) if !(v in ode.parameters) - @error "Multi-experiment identifiability is not properly defined for the states" + error_si("Multi-experiment identifiability is not properly defined for the states") throw(ArgumentError("State variable $v appears in $f")) end end @@ -256,7 +262,7 @@ function assess_local_identifiability( end # Computing the prime using Proposition 3.3 from https://doi.org/10.1006/jsco.2002.0532 - @debug "Computing the prime number" + debug_si("Computing the prime number") d, h = 1, 1 for f in vcat(collect(values(ode.x_equations)), collect(values(ode.y_equations))) df, hf = get_degree_and_coeffsize(f) @@ -281,17 +287,17 @@ function assess_local_identifiability( 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" + debug_si("The prime is $prime") F = Nemo.GF(prime) - @debug "Extending the model" + debug_si("Extending the model") ode_ext = add_outputs(ode, Dict("loc_aux_$i" => f for (i, f) in enumerate(funcs_to_check))) - @debug "Reducing the system modulo prime" + debug_si("Reducing the system modulo prime") ode_red = reduce_ode_mod_p(ode_ext, prime) - @debug "Computing the observability matrix (and, if ME, the bound)" + debug_si("Computing the observability matrix (and, if ME, the bound)") prec = length(ode.x_vars) + length(ode.parameters) # Parameter values are the same across all the replicas @@ -311,10 +317,10 @@ function assess_local_identifiability( u => [F(rand(1:prime)) for i in 1:prec] for u in ode_red.u_vars ) - @debug "Computing the output derivatives" + debug_si("Computing the output derivatives") output_derivatives = differentiate_output(ode_red, params_vals, ic, inputs, prec) - @debug "Building the matrices" + debug_si("Building the matrices") newJac = vcat(Jac, zero(Nemo.MatrixSpace(F, length(ode.x_vars), ncols(Jac)))) newJac = hcat( newJac, @@ -353,7 +359,7 @@ function assess_local_identifiability( end if !isnothing(trbasis) - @debug "Transcendence basis computation requested" + debug_si("Transcendence basis computation requested") reverted_Jac = zero(Nemo.MatrixSpace(F, size(Jac)[2], size(Jac)[1])) for i in 1:size(Jac)[1] for j in 1:size(Jac)[2] @@ -393,10 +399,10 @@ function assess_local_identifiability( for i in trbasis_indices_states push!(trbasis, ode.x_vars[i]) end - @debug "Transcendence basis $trbasis" + debug_si("Transcendence basis $trbasis") end - @debug "Computing the result" + debug_si("Computing the result") base_rank = LinearAlgebra.rank(Jac) result = Dict{Any, Bool}() for i in 1:length(funcs_to_check) diff --git a/src/parametrizations.jl b/src/parametrizations.jl index 8d6f78cae..f9019765d 100644 --- a/src/parametrizations.jl +++ b/src/parametrizations.jl @@ -53,35 +53,35 @@ function check_constructive_field_membership( end # Compute the remainders module the MQS. # These belong to K(T)[x]. - @debug """ - Reducing $tagged_num, $tagged_den""" + debug_si(""" + Reducing $tagged_num, $tagged_den""") _, num_rem = divrem(tagged_num, tagged_mqs) _, den_rem = divrem(tagged_den, tagged_mqs) - @debug """ + debug_si(""" Normal forms modulo MQS: Num: $(num_rem) - Den: $(den_rem)""" + Den: $(den_rem)""") common_factor = gcd(num_rem, den_rem) num_rem = divexact(num_rem, common_factor) den_rem = divexact(den_rem, common_factor) # If norm_form(Num) // norm_form(Den) does not belongs to K(T), then # the fraction does not belong to the field if iszero(den_rem) - @warn """ + warn_si(""" The element $tagged_num // $tagged_den is not in the sub-field Normal form, numerator: $num_rem Normal form, denominator: $den_rem Common factor: $(common_factor) - """ + """) return false, zero(ring_of_tags) // one(ring_of_tags) end if total_degree(num_rem) > 0 || total_degree(den_rem) > 0 - @warn """ + warn_si(""" The element $tagged_num // $tagged_den is not in the sub-field Normal form, numerator: $num_rem Normal form, denominator: $den_rem Common factor: $(common_factor) - """ + """) return false, zero(ring_of_tags) // one(ring_of_tags) end # Now we know that the remainders belong to K(T). @@ -95,11 +95,11 @@ function check_constructive_field_membership( _, num_den_factored = divrem(num_den, tag_relations) num_factored = num_num_factored // denominator(num) den_factored = num_den_factored // denominator(den) - @debug """ + debug_si(""" After factoring out relations: Num: $(num_factored) Den: $(den_factored) - """ + """) if iszero(den_factored) return false, zero(ring_of_tags) // one(ring_of_tags) end @@ -168,12 +168,12 @@ function check_constructive_field_membership( gen_tag_names(length(fracs_gen), "Tag") end sat_string = gen_tag_name("Sat") - @debug """ + debug_si(""" Tags: $(join(map(x -> string(x[1]) * " -> " * string(x[2]), zip(fracs_gen, tag_strings)), "\t\n")) Saturation tag: $sat_string - """ + """) poly_ring_tag, vars_tag = PolynomialRing(K, vcat(sat_string, orig_strings, tag_strings)) sat_var = vars_tag[1] orig_vars = vars_tag[2:(nvars(poly_ring) + 1)] @@ -198,11 +198,11 @@ function check_constructive_field_membership( # - the GB of the MQS in K(T)[x][t]. # ord = Lex() ord = DegRevLex([sat_var]) * DegRevLex(orig_vars) * DegRevLex(tag_vars) - @debug """ + debug_si(""" Tagged MQS ideal: $tagged_mqs Monom ordering: - $(ord)""" + $(ord)""") tagged_mqs_gb = groebner(tagged_mqs, ordering = ord, homogenize = :no) # Relations between tags in K[T] relations_between_tags = filter( @@ -212,12 +212,12 @@ function check_constructive_field_membership( # The basis in K[T][x] tagged_mqs_gb = setdiff(tagged_mqs_gb, relations_between_tags) tagged_mqs_gb = filter(poly -> isempty(intersect(vars(poly), [sat_var])), tagged_mqs_gb) - @debug """ + debug_si(""" Tagged MQS GB: $tagged_mqs_gb Relations between tags: $relations_between_tags - """ + """) # Reduce the fractions with respect to the MQS ideal. # # NOTE: reduction actually happens in K(T)[x]. So we map polynomials to the @@ -225,10 +225,10 @@ function check_constructive_field_membership( ring_of_tags, tags = PolynomialRing(K, tag_strings) tag_to_gen = Dict(tags[i] => fracs_gen[i] for i in 1:length(fracs_gen)) if !isempty(intersect(tag_strings, orig_strings)) - @warn """ + warn_si(""" There is an intersection between the names of the tag variables and the original variables. Tags: $tag_strings - Original vars: $orig_strings""" + Original vars: $orig_strings""") end parametric_ring, _ = PolynomialRing(FractionField(ring_of_tags), orig_strings, ordering = :degrevlex) @@ -238,18 +238,18 @@ function check_constructive_field_membership( Dict(gens(poly_ring_tag)[2:(nvars(poly_ring) + 1)] .=> gens(parametric_ring)), Dict(gens(poly_ring_tag)[(nvars(poly_ring) + 2):end] .=> gens(ring_of_tags)), ) - @debug """ + debug_si(""" Variable mapping: $param_var_mapping Parametric ring: $parametric_ring - """ + """) tagged_mqs_gb_param = map( poly -> crude_parent_ring_change(poly, parametric_ring, param_var_mapping), tagged_mqs_gb, ) tagged_mqs_gb_param = map(f -> divexact(f, leading_coefficient(f)), tagged_mqs_gb_param) - @debug "Tagged parametric mqs: $tagged_mqs_gb_param" + debug_si("Tagged parametric mqs: $tagged_mqs_gb_param") # Reduce each fraction var_mapping = Dict(gens(poly_ring) .=> gens(parametric_ring)) memberships = Vector{Bool}(undef, length(to_be_reduced)) @@ -303,7 +303,7 @@ function reparametrize_with_respect_to(ode, new_states, new_params) # Compute the new dynamics in terms of the original variables. # Paying attenton to the order.. new_vector_field = vector_field_along(ode.x_equations, new_states) - @debug "New vector field:\n$new_vector_field" + debug_si("New vector field:\n$new_vector_field") new_states = collect(keys(new_vector_field)) new_dynamics = [new_vector_field[new_state] for new_state in new_states] # Express the new dynamics in terms of new states and new parameters. @@ -322,14 +322,14 @@ function reparametrize_with_respect_to(ode, new_states, new_params) gen_tag_names(length(ode.u_vars), "Input"), gen_tag_names(length(ode.y_vars), "Output"), ) - @info """ + info_si(""" Tag names: $tag_names Generating functions: $generating_funcs Functions to be reduced: $to_be_reduced_funcs - """ + """) membership, new_dynamics_all, implicit_relations, new_vars = check_constructive_field_membership( RationalFunctionField(generating_funcs), @@ -351,20 +351,20 @@ function reparametrize_with_respect_to(ode, new_states, new_params) if !isempty(ode.u_vars) new_inputs = tag_inputs end - @info """ + info_si(""" New state dynamics: $new_dynamics_states New output dynamics: $new_dynamics_outputs New inputs: - $new_inputs""" + $new_inputs""") # Construct the new vector field. new_vars_vector_field = empty(ode.x_equations) for i in 1:length(new_states) state = tags[i] new_vars_vector_field[state] = new_dynamics_states[i] end - @info "Converting variable names to human-readable ones" + info_si("Converting variable names to human-readable ones") internal_variable_names = map(i -> "X$i", 1:length(new_states)) parameter_variable_names = map(i -> "a$i", 1:length(new_params)) input_variable_names = map(i -> "u$i", 1:length(tag_inputs)) @@ -464,27 +464,28 @@ Dict{Nemo.QQMPolyRingElem, AbstractAlgebra.Generic.Frac{Nemo.QQMPolyRingElem}} w Notice that the `new_ode` is fully identifiabile, and has `1` less parameter compared to the original one. """ -function reparametrize_global(ode::ODE{P}; p = 0.99, seed = 42) where {P} +function reparametrize_global(ode::ODE{P}; p = 0.99, seed = 42, loglevel = Logging.Info) where {P} + restart_logging(loglevel = loglevel) Random.seed!(seed) id_funcs = find_identifiable_functions(ode, with_states = true, simplify = :strong, p = p) ode_ring = parent(ode) @assert base_ring(parent(first(id_funcs))) == ode_ring - @info "Constructing a new parametrization" + info_si("Constructing a new parametrization") contains_states(poly::MPolyElem) = any(x -> degree(poly, x) > 0, ode.x_vars) contains_states(func) = contains_states(numerator(func)) || contains_states(denominator(func)) id_funcs_contains_states = filter(contains_states, id_funcs) - @info """ + info_si(""" Original states: $(ode.x_vars) Original params: $(ode.parameters) - Identifiable and contain states: $(id_funcs_contains_states)""" + Identifiable and contain states: $(id_funcs_contains_states)""") new_states = id_funcs_contains_states new_params = setdiff(id_funcs, id_funcs_contains_states) - @info """ + info_si(""" Reparametrizing with respect to: New states: $new_states - New params: $new_params""" + New params: $new_params""") new_vector_field, new_inputs, new_outputs, new_vars, implicit_relations = reparametrize_with_respect_to(ode, new_states, new_params) new_ring = parent(first(keys(new_vector_field))) diff --git a/src/pb_representation.jl b/src/pb_representation.jl index 3594b143a..d4d3cc0fa 100644 --- a/src/pb_representation.jl +++ b/src/pb_representation.jl @@ -91,7 +91,7 @@ function common_ring(poly::MPolyElem, pbr::PBRepresentation) if d != nothing max_ords[d[1]] = max(d[2], max_ords[d[1]]) elseif !(var_to_str(v) in pbr.param_names) - @warn "New variable $(var_to_str(v)), treating as a scalar parameter" + warn_si("New variable $(var_to_str(v)), treating as a scalar parameter") push!(new_params, var_to_str(v)) end end diff --git a/src/power_series_utils.jl b/src/power_series_utils.jl index 3124d47c0..276bf3e78 100644 --- a/src/power_series_utils.jl +++ b/src/power_series_utils.jl @@ -261,7 +261,7 @@ function ps_ode_solution( cur_prec = 1 while cur_prec < prec - @debug "\t Computing power series solution, currently at precision $cur_prec" + debug_si("\t Computing power series solution, currently at precision $cur_prec") new_prec = min(prec, 2 * cur_prec) for i in 1:length(x_vars) set_precision!(solution[x_vars[i]], new_prec) diff --git a/src/precompile.jl b/src/precompile.jl index 415dc0577..3136b1026 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -20,10 +20,10 @@ x3'(t) = x3(t), y(t) = x2(t) ) - assess_identifiability(ode) - assess_identifiability(de; measured_quantities = [x0]) - assess_identifiability(de; measured_quantities = [y ~ x0]) - find_identifiable_functions(ode, with_states = true) + assess_identifiability(ode, loglevel = Logging.Warn) + assess_identifiability(de; measured_quantities = [x0], loglevel = Logging.Warn) + assess_identifiability(de; measured_quantities = [y ~ x0], loglevel = Logging.Warn) + find_identifiable_functions(ode, with_states = true, loglevel = Logging.Warn) restart_logging(loglevel = Logging.Info) end end diff --git a/src/primality_check.jl b/src/primality_check.jl index 220188a9f..253e579cc 100644 --- a/src/primality_check.jl +++ b/src/primality_check.jl @@ -6,8 +6,8 @@ function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) dim = length(basis) S = Nemo.MatrixSpace(Nemo.QQ, dim, dim) matrices = [] - @debug "" J basis - @debug "Dim is $dim" + debug_si("$J $basis") + debug_si("Dim is $dim") for v in gens(parent(first(J))) M = zero(S) for (i, vec) in enumerate(basis) @@ -17,13 +17,13 @@ function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) end end push!(matrices, M) - @debug "Multiplication by $v" M + debug_si("Multiplication by $v: $M") end generic_multiplication = sum(Nemo.QQ(rand(1:100)) * M for M in matrices) - @debug "" generic_multiplication + debug_si("$generic_multiplication") R, t = Nemo.PolynomialRing(Nemo.QQ, "t") - @debug "" Nemo.charpoly(R, generic_multiplication) + debug_si("$(Nemo.charpoly(R, generic_multiplication))") return Nemo.isirreducible(Nemo.charpoly(R, generic_multiplication)) end diff --git a/src/util.jl b/src/util.jl index af9d95c2d..c7b279a49 100644 --- a/src/util.jl +++ b/src/util.jl @@ -167,13 +167,11 @@ function make_substitution( d = Nemo.degree(f, var_sub) result = 0 - @debug "Substitution in a polynomial of degree $d" - flush(stdout) + debug_si("Substitution in a polynomial of degree $d") for i in 0:d - @debug "\t Degree $i" - flush(stdout) + debug_si("\t Degree $i") result += coeff(f, [var_sub], [i]) * (val_numer^i) * (val_denom^(d - i)) - @debug "\t Intermediate result of size $(length(result))" + debug_si("\t Intermediate result of size $(length(result))") end return result end @@ -313,9 +311,7 @@ function uncertain_factorization(f::MPolyElem{fmpq}) mainvar_coeffs = [coeff(f, [main_var], [i]) for i in 0:d] gcd_coef = mainvar_coeffs[end] for i in d:-1:1 - #@info "Degrees $(total_degree(gcd_coef)) and $(total_degree(mainvar_coeffs[i]))" gcd_coef = gcd(gcd_coef, mainvar_coeffs[i]) - #@info "Time $tm and new degree $(total_degree(gcd_coef))" end f = divexact(f, gcd_coef) @@ -330,11 +326,10 @@ function uncertain_factorization(f::MPolyElem{fmpq}) if degree(factor_out, main_var) != degree(gcd(f_uni, derivative(f_uni))) continue end - @debug "Nonsquarefree poly, dividing by $factor_out" + debug_si("Nonsquarefree poly, dividing by $factor_out") f = divexact(f, factor_out) f_uni = divexact(f_uni, gcd(f_uni, derivative(f_uni))) end - # end is_irr = Nemo.isirreducible(f_uni) break end @@ -526,7 +521,7 @@ function eval_at_nemo(e::Union{Float16, Float32, Float64}, vals::Dict) else out = rationalize(e) end - @warn "Floating point value $e will be converted to $(out)." + warn_si("Floating point value $e will be converted to $(out).") return out end @@ -571,7 +566,7 @@ end function get_measured_quantities(ode::ModelingToolkit.ODESystem) if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) - @info "Measured quantities are not provided, trying to find the outputs in input ODE." + info_si("Measured quantities are not provided, trying to find the outputs in input ODE.") return filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(ode), diff --git a/src/wronskian.jl b/src/wronskian.jl index 9885f0f88..644eaa5af 100644 --- a/src/wronskian.jl +++ b/src/wronskian.jl @@ -198,9 +198,9 @@ Output: Computes the Wronskians of io_equations """ function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} - @debug "Compressing monomials" + debug_si("Compressing monomials") termlists = [monomial_compress(ioeq, ode)[2] for ioeq in values(io_equations)] - @debug "Matrix sizes $(map(length, termlists))" + debug_si("Matrix sizes $(map(length, termlists))") # estimating the required order of truncation ord = max(map(length, termlists)...) + length(ode.x_vars) @@ -214,7 +214,7 @@ function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} [map(p -> parent_ring_change(p, polyring_red), tlist) for tlist in termlists] ode_red = reduce_ode_mod_p(ode, PRIME) - @debug "Computing power series solution up to order $ord" + debug_si("Computing power series solution up to order $ord") ps = power_series_solution( ode_red, Dict(p => rand(1:100) for p in ode_red.parameters), @@ -222,7 +222,7 @@ function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} Dict(u => [rand(1:100) for i in 0:ord] for u in ode_red.u_vars), ord, ) - @debug "Computing the derivatives of the solution" + debug_si("Computing the derivatives of the solution") ps_ext = Dict{MPolyElem, Nemo.gfp_abs_series}()# Generic.AbsSeries}() for v in vcat(ode_red.y_vars, ode_red.u_vars) cur_ps = ps[v] @@ -232,7 +232,7 @@ function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} end end - @debug "Constructing Wronskians" + debug_si("Constructing Wronskians") result = [] for (i, tlist) in enumerate(termlists) n = length(tlist) From 0a4e694769861b390470284d84e77764add8caa6 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Fri, 10 Nov 2023 21:55:33 +0100 Subject: [PATCH 05/24] fixing typo, tests now pass --- src/local_identifiability.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index e32e0abfb..392e1421d 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -171,7 +171,7 @@ function assess_local_identifiability( type = :SE, loglevel = Logging.Info, ) - restart_logging(logelevel = loglevel) + restart_logging(loglevel = loglevel) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) info_si("Measured quantities are not provided, trying to find the outputs in input ODE.") From c2be531572867d6f8f549591d17dd82fff5e6386 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Fri, 10 Nov 2023 22:19:19 +0100 Subject: [PATCH 06/24] formatting --- src/ODE.jl | 16 +++++++--- src/RationalFunctionFields/IdealMQS.jl | 4 ++- .../RationalFunctionField.jl | 28 ++++++++++++----- src/RationalFunctionFields/normalforms.jl | 16 +++++++--- src/StructuralIdentifiability.jl | 7 ++++- src/discrete.jl | 8 +++-- src/global_identifiability.jl | 12 +++++-- src/identifiable_functions.jl | 6 ++-- src/io_equation.jl | 14 ++++++--- src/local_identifiability.jl | 12 ++++--- src/logging.jl | 12 +++---- src/parametrizations.jl | 31 ++++++++++++------- src/precompile.jl | 6 +++- src/util.jl | 4 ++- 14 files changed, 121 insertions(+), 55 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 887cdf95c..98835b10d 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -461,10 +461,18 @@ macro ODEmodel(ex::Expr...) logging_exprs = [ :(StructuralIdentifiability.info_si("Summary of the model:")), - :(StructuralIdentifiability.info_si("State variables: " * $(join(map(string, collect(x_vars)), ", ")))), - :(StructuralIdentifiability.info_si("Parameters: " * $(join(map(string, collect(params)), ", ")))), - :(StructuralIdentifiability.info_si("Inputs: " * $(join(map(string, collect(u_vars)), ", ")))), - :(StructuralIdentifiability.info_si("Outputs: " * $(join(map(string, collect(y_vars)), ", ")))), + :(StructuralIdentifiability.info_si( + "State variables: " * $(join(map(string, collect(x_vars)), ", ")), + )), + :(StructuralIdentifiability.info_si( + "Parameters: " * $(join(map(string, collect(params)), ", ")), + )), + :(StructuralIdentifiability.info_si( + "Inputs: " * $(join(map(string, collect(u_vars)), ", ")), + )), + :(StructuralIdentifiability.info_si( + "Outputs: " * $(join(map(string, collect(y_vars)), ", ")), + )), ] # creating the ode object ode_expr = :(StructuralIdentifiability.ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}( diff --git a/src/RationalFunctionFields/IdealMQS.jl b/src/RationalFunctionFields/IdealMQS.jl index 8a2e760bc..e1d9057dd 100644 --- a/src/RationalFunctionFields/IdealMQS.jl +++ b/src/RationalFunctionFields/IdealMQS.jl @@ -86,7 +86,9 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal lcm, 1:length(funcs_den_nums), ) - debug_si("Rational functions common denominator is of degree $(total_degree(den_lcm)) and of length $(length(den_lcm))") + debug_si( + "Rational functions common denominator is of degree $(total_degree(den_lcm)) and of length $(length(den_lcm))", + ) is_constant(den_lcm) && (debug_si("Common denominator of the field generators is constant")) existing_varnames = map(String, symbols(ring)) diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index 8835e19d9..113d59ee6 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -255,7 +255,9 @@ function beautifuly_generators( end end end - debug_si("Out of $(length(fracs)) simplified generators there are $(length(non_redundant)) non redundant") + debug_si( + "Out of $(length(fracs)) simplified generators there are $(length(non_redundant)) non redundant", + ) fracs = fracs[non_redundant] end sort!(fracs, lt = rational_function_cmp) @@ -356,7 +358,9 @@ function groebner_basis_coeffs( debug_si("Inclusion checked in $(runtime) seconds. Result: $two_sided_inclusion") current_degrees = current_degrees .* 2 end - debug_si("The coefficients of the Groebner basis are presented by $(length(fracs)) rational functions") + debug_si( + "The coefficients of the Groebner basis are presented by $(length(fracs)) rational functions", + ) new_rff.mqs.cached_groebner_bases[ordering, up_to_degree] = gb rff.mqs.cached_groebner_bases[ordering, up_to_degree] = gb return new_rff @@ -382,7 +386,9 @@ function generating_sets_fan( time_start = time_ns() vars = gens(parent(rff.mqs)) nbases = length(vars) - info_si("Computing $nbases Groebner bases for degrees $up_to_degree for block orderings") + info_si( + "Computing $nbases Groebner bases for degrees $up_to_degree for block orderings", + ) ordering_to_generators = Dict() if code == 0 return ordering_to_generators @@ -552,9 +558,11 @@ function simplified_generating_set( append!(new_fracs, generators) end new_fracs_unique = unique(new_fracs) - debug_si(""" - Final cleaning and simplification of generators. - Out of $(length(new_fracs)) fractions $(length(new_fracs_unique)) are syntactically unique.""") + debug_si( + """ +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_si("Checking inclusion with probability $p") @@ -564,8 +572,12 @@ function simplified_generating_set( warn_si("Field membership check failed. Error will follow.") throw("The new subfield generators are not correct.") end - info_si("Inclusion checked with probability $p in $(_runtime_logger[:id_inclusion_check]) seconds") - debug_si("Out of $(length(rff.mqs.nums_qq)) initial generators there are $(length(new_fracs)) indepdendent") + info_si( + "Inclusion checked with probability $p in $(_runtime_logger[:id_inclusion_check]) seconds", + ) + debug_si( + "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 debug_si("The ranking of the new set of generators is $ranking") diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index c7f2c7c9a..23efd94eb 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -365,9 +365,11 @@ function linear_relations_between_normal_forms( push!(zeroed_relations_inds, i) end end - debug_si(""" - Relations in the previous intersection: $(length(complete_intersection_relations_ff)) - Vanished at the current point: $(length(zeroed_relations_inds))""") + debug_si( + """ + Relations in the previous intersection: $(length(complete_intersection_relations_ff)) + Vanished at the current point: $(length(zeroed_relations_inds))""", + ) non_zeroed_relations_inds = setdiff(collect(1:n_relations_ff), zeroed_relations_inds) zeroed_relations_from_complete_intersection = @@ -383,7 +385,9 @@ function linear_relations_between_normal_forms( non_zeroed_relations_from_complete_intersection, zeroed_relations_from_complete_intersection, ) - debug_si("There are $(length(complete_intersection_relations_ff)) relations in the intersection") + debug_si( + "There are $(length(complete_intersection_relations_ff)) relations in the intersection", + ) m_relations_ff = length(complete_intersection_relations_ff) if n_relations_ff == m_relations_ff || isempty(complete_intersection_relations_ff) break @@ -412,6 +416,8 @@ function linear_relations_between_normal_forms( ) relations_qq[i] = relation_qq_param // one(relation_qq_param) end - info_si("Used $iters specializations in $((time_ns() - time_start) / 1e9) seconds, found $(length(complete_intersection_relations_ff)) relations") + info_si( + "Used $iters specializations in $((time_ns() - time_start) / 1e9) seconds, found $(length(complete_intersection_relations_ff)) relations", + ) relations_qq end diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index a36e0361f..aacae3d56 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -176,7 +176,12 @@ 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, loglevel = loglevel) + result = assess_identifiability( + ode, + funcs_to_check = funcs_to_check_, + p = p, + loglevel = loglevel, + ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) return out_dict diff --git a/src/discrete.jl b/src/discrete.jl index 1c29e43b0..b964570ac 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -314,7 +314,9 @@ function assess_local_identifiability( ) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(dds)) - info_si("Measured quantities are not provided, trying to find the outputs in input dynamical system.") + info_si( + "Measured quantities are not provided, trying to find the outputs in input dynamical system.", + ) measured_quantities = filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(dds), @@ -342,7 +344,9 @@ function assess_local_identifiability( nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) if length(known_ic) > 0 - warn_si("Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !") + warn_si( + "Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !", + ) end return out_dict end diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index d18fc2983..557f7abec 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -130,7 +130,9 @@ function initial_identifiable_functions( _runtime_logger[:rank_time] = rank_times if any([dim != rk + 1 for (dim, rk) in zip(dims, wranks)]) - warn_si("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) or transforming all the parameters to states with zero derivative") + warn_si( + "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) or transforming all the parameters to states with zero derivative", + ) end end @@ -142,7 +144,9 @@ function initial_identifiable_functions( ) if with_states && !isempty(ode.parameters) - debug_si("Generators of identifiable functions involve states, the parameter-only part is getting simplified") + debug_si( + "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( base_ring(bring), @@ -305,7 +309,9 @@ function assess_global_identifiability( ) where {P <: MPolyElem{fmpq}} submodels = find_submodels(ode) if length(submodels) > 0 - info_si("Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`") + info_si( + "Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`", + ) end result = check_identifiability( diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 43a66090b..2a49ce7e2 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -51,7 +51,7 @@ function find_identifiable_functions( rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, ) where {T <: MPolyElem{fmpq}} - restart_logging(loglevel=loglevel) + restart_logging(loglevel = loglevel) Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) runtime_start = time_ns() @@ -91,7 +91,9 @@ function find_identifiable_functions( id_funcs_fracs = [parent_ring_change(f, parent(ode)) for f in id_funcs_fracs] _runtime_logger[:id_total] = (time_ns() - runtime_start) / 1e9 _runtime_logger[:are_id_funcs_polynomial] = all(isone ∘ denominator, id_funcs_fracs) - info_si("The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds") + info_si( + "The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds", + ) return id_funcs_fracs end diff --git a/src/io_equation.jl b/src/io_equation.jl index 92a77a9dc..82048f8e2 100644 --- a/src/io_equation.jl +++ b/src/io_equation.jl @@ -215,7 +215,9 @@ function find_ioprojections( A, B = simplify_frac(A, B) if isempty(filter!(v -> (v in keys(x_equations)), vars(A))) && (B != 0) # && (length(coeffs(A))==1) # variable change x_i' --> x_i' - (B/A)', x_i --> x_i - B/A - debug_si("\t Applying variable change: $(x) --> $(x) - ( $B )/( $A )") + debug_si( + "\t Applying variable change: $(x) --> $(x) - ( $B )/( $A )", + ) dB = diff_poly(B, derivation) dA = diff_poly(A, derivation) numer_d, denom_d = simplify_frac(A * dB - dA * B, A * A) @@ -339,7 +341,7 @@ Output: function find_ioequations( ode::ODE{P}; var_change_policy = :default, - loglevel = Logging.Info + loglevel = Logging.Info, ) where {P <: MPolyElem{<:FieldElem}} # Setting the var_change policy if (var_change_policy == :yes) || @@ -358,13 +360,17 @@ function find_ioequations( debug_si("Check whether the original projections are enough") if length(io_projections) == 1 || check_primality(io_projections) - debug_si("The projections generate an ideal with a single components of highest dimension, returning") + debug_si( + "The projections generate an ideal with a single components of highest dimension, returning", + ) return io_projections end sampling_range = 5 while true - debug_si("There are several components of the highest dimension, trying to isolate one") + debug_si( + "There are several components of the highest dimension, trying to isolate one", + ) extra_projection = sum(rand(1:sampling_range) * v for v in keys(io_projections)) debug_si("Extra projections: $extra_projection") new_projections, _, projection_equation = diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 392e1421d..93fe08952 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -174,7 +174,9 @@ function assess_local_identifiability( restart_logging(loglevel = loglevel) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) - info_si("Measured quantities are not provided, trying to find the outputs in input ODE.") + info_si( + "Measured quantities are not provided, trying to find the outputs in input ODE.", + ) measured_quantities = filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(ode), @@ -202,7 +204,7 @@ function assess_local_identifiability( funcs_to_check = funcs_to_check_, p = p, type = type, - loglevel = loglevel + loglevel = loglevel, ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) @@ -238,7 +240,7 @@ function assess_local_identifiability( p::Float64 = 0.99, type = :SE, trbasis = nothing, - loglevel = Logging.Info + loglevel = Logging.Info, ) where {P <: MPolyElem{Nemo.fmpq}} restart_logging(loglevel = loglevel) if isempty(funcs_to_check) @@ -254,7 +256,9 @@ function assess_local_identifiability( num, den = unpack_fraction(f) for v in vcat(vars(num), vars(den)) if !(v in ode.parameters) - error_si("Multi-experiment identifiability is not properly defined for the states") + error_si( + "Multi-experiment identifiability is not properly defined for the states", + ) throw(ArgumentError("State variable $v appears in $f")) end end diff --git a/src/logging.jl b/src/logging.jl index 946163c21..a51940ff6 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -8,8 +8,8 @@ # not the exhaustive list, just the ones to be initialized const _runtime_rubrics = ( - :id_calls_to_gb, - :id_groebner_time, + :id_calls_to_gb, + :id_groebner_time, :id_inclusion_check_mod_p, :id_npoints_degree, :id_npoints_interpolation, @@ -25,12 +25,8 @@ const _runtime_logger = Dict( :id_beautifulization => 0, ) -const _si_logger = Ref{Logging.ConsoleLogger}( - Logging.ConsoleLogger( - Logging.Info, - show_limited = false, - ) -) +const _si_logger = + Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info, show_limited = false)) function restart_logging(; loglevel = Logging.Info) _si_logger[] = Logging.ConsoleLogger(loglevel, show_limited = false) diff --git a/src/parametrizations.jl b/src/parametrizations.jl index f9019765d..b1aa527ba 100644 --- a/src/parametrizations.jl +++ b/src/parametrizations.jl @@ -168,12 +168,14 @@ function check_constructive_field_membership( gen_tag_names(length(fracs_gen), "Tag") end sat_string = gen_tag_name("Sat") - debug_si(""" - Tags: - $(join(map(x -> string(x[1]) * " -> " * string(x[2]), zip(fracs_gen, tag_strings)), "\t\n")) - Saturation tag: - $sat_string - """) + debug_si( + """ +Tags: +$(join(map(x -> string(x[1]) * " -> " * string(x[2]), zip(fracs_gen, tag_strings)), "\t\n")) +Saturation tag: +$sat_string +""", + ) poly_ring_tag, vars_tag = PolynomialRing(K, vcat(sat_string, orig_strings, tag_strings)) sat_var = vars_tag[1] orig_vars = vars_tag[2:(nvars(poly_ring) + 1)] @@ -225,10 +227,12 @@ function check_constructive_field_membership( ring_of_tags, tags = PolynomialRing(K, tag_strings) tag_to_gen = Dict(tags[i] => fracs_gen[i] for i in 1:length(fracs_gen)) if !isempty(intersect(tag_strings, orig_strings)) - warn_si(""" - There is an intersection between the names of the tag variables and the original variables. - Tags: $tag_strings - Original vars: $orig_strings""") + warn_si( + """ + There is an intersection between the names of the tag variables and the original variables. + Tags: $tag_strings + Original vars: $orig_strings""", + ) end parametric_ring, _ = PolynomialRing(FractionField(ring_of_tags), orig_strings, ordering = :degrevlex) @@ -464,7 +468,12 @@ Dict{Nemo.QQMPolyRingElem, AbstractAlgebra.Generic.Frac{Nemo.QQMPolyRingElem}} w Notice that the `new_ode` is fully identifiabile, and has `1` less parameter compared to the original one. """ -function reparametrize_global(ode::ODE{P}; p = 0.99, seed = 42, loglevel = Logging.Info) where {P} +function reparametrize_global( + ode::ODE{P}; + p = 0.99, + seed = 42, + loglevel = Logging.Info, +) where {P} restart_logging(loglevel = loglevel) Random.seed!(seed) id_funcs = diff --git a/src/precompile.jl b/src/precompile.jl index 3136b1026..9b20445ab 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -22,7 +22,11 @@ ) assess_identifiability(ode, loglevel = Logging.Warn) assess_identifiability(de; measured_quantities = [x0], loglevel = Logging.Warn) - assess_identifiability(de; measured_quantities = [y ~ x0], loglevel = Logging.Warn) + assess_identifiability( + de; + measured_quantities = [y ~ x0], + loglevel = Logging.Warn, + ) find_identifiable_functions(ode, with_states = true, loglevel = Logging.Warn) restart_logging(loglevel = Logging.Info) end diff --git a/src/util.jl b/src/util.jl index c7b279a49..e51b3e313 100644 --- a/src/util.jl +++ b/src/util.jl @@ -566,7 +566,9 @@ end function get_measured_quantities(ode::ModelingToolkit.ODESystem) if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) - info_si("Measured quantities are not provided, trying to find the outputs in input ODE.") + info_si( + "Measured quantities are not provided, trying to find the outputs in input ODE.", + ) return filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(ode), From 265adc654cca718b3269368af79509b13ba10bba Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Sat, 11 Nov 2023 23:20:07 +0300 Subject: [PATCH 07/24] add TimerOutputs --- Project.toml | 1 + .../IdentifiableFunctions/experiments.jl | 9 ++- .../RationalFunctionField.jl | 12 ++-- src/RationalFunctionFields/normalforms.jl | 2 +- src/StructuralIdentifiability.jl | 7 ++- src/elimination.jl | 2 +- src/global_identifiability.jl | 8 +-- src/identifiable_functions.jl | 3 + src/io_equation.jl | 4 +- src/local_identifiability.jl | 4 ++ src/logging.jl | 55 +++++++++++++++++++ src/states.jl | 2 +- src/wronskian.jl | 2 +- 13 files changed, 90 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 2d698faa2..d3dbfe86b 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" [compat] AbstractAlgebra = "0.13, 0.18, 0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33" diff --git a/benchmarking/IdentifiableFunctions/experiments.jl b/benchmarking/IdentifiableFunctions/experiments.jl index 72dada132..e63eb2124 100644 --- a/benchmarking/IdentifiableFunctions/experiments.jl +++ b/benchmarking/IdentifiableFunctions/experiments.jl @@ -988,14 +988,17 @@ begin ) using Nemo, Logging - using JuliaInterpreter + # using JuliaInterpreter Groebner = StructuralIdentifiability.Groebner # ParamPunPam = StructuralIdentifiability.ParamPunPam Base.global_logger(ConsoleLogger(Logging.Info)) end -dennums, ring = - StructuralIdentifiability.initial_identifiable_functions(Pivastatin, p = 0.99); +begin + StructuralIdentifiability.find_identifiable_functions(Pivastatin); + StructuralIdentifiability.print_timings_table(); +end + fracs = StructuralIdentifiability.dennums_to_fractions(dennums); rff = StructuralIdentifiability.RationalFunctionField(dennums); diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index 113d59ee6..2f8378bfa 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -92,7 +92,7 @@ Output: - a list `L[i]` of bools of length `length(rat_funcs)` such that `L[i]` is true iff the i-th function belongs to `field` """ -function field_contains( +@timeit _to function field_contains( field::RationalFunctionField{T}, ratfuncs::Vector{Vector{T}}, p, @@ -176,7 +176,7 @@ end # ------------------------------------------------------------------------------ -function check_field_membership_mod_p!( +@timeit _to function check_field_membership_mod_p!( generators::RationalFunctionField{T}, tobereduced::RationalFunctionField{T}, ) where {T} @@ -209,7 +209,7 @@ Applies the following passes: 1. Filter constants, 2. Remove redundant generators. """ -function beautifuly_generators( +@timeit _to function beautifuly_generators( rff::RationalFunctionField; discard_redundant = true, reversed_order = false, @@ -303,7 +303,7 @@ end - `up_to_degree`: a tuple of integers, the degrees of numerator and denominator. The result is correct up to the requested degrees. """ -function groebner_basis_coeffs( +@timeit _to function groebner_basis_coeffs( rff::RationalFunctionField; seed = 42, ordering = Groebner.InputOrdering(), @@ -377,7 +377,7 @@ Returns a set of Groebner bases for multiple different rankings of variables. - Keyword `up_to_degree`: a tuple of integers, max. degrees of numerators and denominators. Result is correct up to the requested degrees. """ -function generating_sets_fan( +@timeit _to function generating_sets_fan( rff::RationalFunctionField{T}, code::Integer; seed = 42, @@ -486,7 +486,7 @@ end Returns a simplified set of generators for `rff`. Result is correct (in Monte-Carlo sense) with probability at least `p`. """ -function simplified_generating_set( +@timeit _to function simplified_generating_set( rff::RationalFunctionField; p = 0.99, seed = 42, diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index 23efd94eb..0bf45773f 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -290,7 +290,7 @@ Relations may include monomials up to the total degree `up_to_degree`. Note: uses Monte-Carlo probabilistic algorithm. The probability of correctness is not specified but is assumed to be close to 1. """ -function linear_relations_between_normal_forms( +@timeit _to function linear_relations_between_normal_forms( fracs::Vector{Generic.Frac{T}}, up_to_degree::Integer; seed = 42, diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index aacae3d56..fd4db817f 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -10,6 +10,7 @@ using Logging using MacroTools using Primes using Random +using TimerOutputs # Algebra packages using AbstractAlgebra @@ -93,6 +94,8 @@ function assess_identifiability( loglevel = Logging.Info, ) where {P <: MPolyElem{fmpq}} restart_logging(loglevel = loglevel) + reset_timings() + p_glob = 1 - (1 - p) * 0.9 p_loc = 1 - (1 - p) * 0.1 @@ -110,9 +113,9 @@ function assess_identifiability( trbasis = trbasis, loglevel = loglevel, ) - info_si("Local identifiability assessed in $runtime seconds") + # info_si("Local identifiability assessed in $runtime seconds") debug_si("Trasncendence basis to be specialized is $trbasis") - _runtime_logger[:loc_time] = runtime + # _runtime_logger[:loc_time] = runtime loc_id = [local_result[each] for each in funcs_to_check] locally_identifiable = Array{Any, 1}() diff --git a/src/elimination.jl b/src/elimination.jl index 3c68b735a..bab6aa9cc 100644 --- a/src/elimination.jl +++ b/src/elimination.jl @@ -301,7 +301,7 @@ Input: Output: - `polynomial` - the desired factor of the resultant of `f` and `g` """ -function eliminate_var( +@timeit _to function eliminate_var( f::P, g::P, var_elim::P, diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 557f7abec..88f78e7a1 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -12,7 +12,7 @@ are identifiable functions containing or not the state variables - a polynomial ring containing all returned identifiable functions (parameters or parameters + states) """ -function extract_identifiable_functions_raw( +@timeit _to function extract_identifiable_functions_raw( io_equations::Dict{P, P}, ode::ODE{P}, known::Array{P, 1}, @@ -97,7 +97,7 @@ The function returns a tuple containing the following: - a list of identifiable functions (as pairs [num, denum]) - the ring containing all these functuons (either parameters only of with states) """ -function initial_identifiable_functions( +@timeit _to function initial_identifiable_functions( ode::ODE{T}; p::Float64, known::Array{T, 1} = Array{T, 1}(), @@ -192,7 +192,7 @@ Input: Output: a list L of booleans with L[i] being the identifiability status of the i-th function to check """ -function check_identifiability( +@timeit _to function check_identifiability( ode::ODE{P}, funcs_to_check::Array{<:Any, 1}; known::Array{P, 1} = Array{P, 1}(), @@ -300,7 +300,7 @@ Output: Checks global identifiability of functions of parameters specified in `funcs_to_check`. """ -function assess_global_identifiability( +@timeit _to function assess_global_identifiability( ode::ODE{P}, funcs_to_check::Array{<:Any, 1}, known::Array{P, 1} = Array{P, 1}(), diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 2a49ce7e2..c3dcdf3a4 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -52,6 +52,8 @@ function find_identifiable_functions( loglevel = Logging.Info, ) where {T <: MPolyElem{fmpq}} restart_logging(loglevel = loglevel) + reset_timings() + Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) runtime_start = time_ns() @@ -94,6 +96,7 @@ function find_identifiable_functions( info_si( "The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds", ) + return id_funcs_fracs end diff --git a/src/io_equation.jl b/src/io_equation.jl index 82048f8e2..06e5ae0a4 100644 --- a/src/io_equation.jl +++ b/src/io_equation.jl @@ -101,7 +101,7 @@ Output: - generic point generator for the model (including the derivatives; mostly for testing) - an extra projection (if `extra_projection` was provided) """ -function find_ioprojections( +@timeit _to function find_ioprojections( ode::ODE{P}, auto_var_change::Bool, extra_projection = nothing, @@ -338,7 +338,7 @@ Output: - a dictionary from “leaders” to the corresponding input-output equations; if an extra projection is needed, it will be the value corresponding to `rand_proj_var` """ -function find_ioequations( +@timeit _to function find_ioequations( ode::ODE{P}; var_change_policy = :default, loglevel = Logging.Info, diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 93fe08952..0d4b83c8a 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -235,6 +235,8 @@ Call this function if you have a specific collection of parameters of which you If the type is `:ME`, states are not allowed to appear in the `funcs_to_check`. """ function assess_local_identifiability( + # TODO: report this bug (?) to TimerOutputs.jl + # @timeit _to function assess_local_identifiability( ode::ODE{P}; funcs_to_check::Array{<:Any, 1} = Array{Any, 1}(), p::Float64 = 0.99, @@ -243,6 +245,8 @@ function assess_local_identifiability( loglevel = Logging.Info, ) where {P <: MPolyElem{Nemo.fmpq}} restart_logging(loglevel = loglevel) + reset_timings() + if isempty(funcs_to_check) funcs_to_check = ode.parameters if type == :SE diff --git a/src/logging.jl b/src/logging.jl index a51940ff6..4418808c2 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -62,3 +62,58 @@ function error_si(s) flush(stdout) end end + +### +# Timings for StructuralIdentifiability.jl. +# Uses the the @timeit macro from TimerOutputs.jl. +# +# Provides a couple of useful functions: +# - enable_timer, for enabling or disabling the timer in SI +# - reset_timings, for clearing the collected timings +# - print_timings_table, for printing the table with timings + +timeit_debug_enabled() = false + +const _to = TimerOutputs.TimerOutput() + +""" + enable_timer(flag) + +If `flag` is `true`, enables the global timer. Otherwise, disables it. +""" +function enable_timer(flag::Bool) + if flag + enable_timer!(_to) + else + disable_timer!(_to) + end + nothing +end + +""" + reset_timings() + +Resets the global timer. +""" +function reset_timings() + TimerOutputs.reset_timer!(_to) + nothing +end + +""" + print_timings_table() + +Prints the table with collected timings data to `stdout`. +""" +function print_timings_table() + iostream = stdout + TimerOutputs.print_timer( + iostream, + _to, + allocations = true, + sortby = :time, + linechars = :ascii, + compact = false, + title = "SI.jl", + ) +end diff --git a/src/states.jl b/src/states.jl index 56fbd2c14..0936770ec 100644 --- a/src/states.jl +++ b/src/states.jl @@ -76,7 +76,7 @@ Output: all identifiabile functions over the field generated by the inputs and the identifiable functions of parameters only """ -function states_generators( +@timeit _to function states_generators( ode::ODE{P}, io_equations::Dict{P, P}, ) where {P <: MPolyElem{<:FieldElem}} diff --git a/src/wronskian.jl b/src/wronskian.jl index 644eaa5af..c36205449 100644 --- a/src/wronskian.jl +++ b/src/wronskian.jl @@ -197,7 +197,7 @@ Output: Computes the Wronskians of io_equations """ -function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} +@timeit _to function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} debug_si("Compressing monomials") termlists = [monomial_compress(ioeq, ode)[2] for ioeq in values(io_equations)] debug_si("Matrix sizes $(map(length, termlists))") From bf7a857da07dc06ddc5916ea30cb37282b812f60 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Sat, 11 Nov 2023 23:22:27 +0300 Subject: [PATCH 08/24] disable the global timer by default --- src/logging.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/logging.jl b/src/logging.jl index 4418808c2..988144e1d 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -90,6 +90,9 @@ function enable_timer(flag::Bool) nothing end +# By default, the timer is disabled +enable_timer(false) + """ reset_timings() From 51cf5907388f4f2c79ec06dc3f1c6fe50074032c Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Sun, 12 Nov 2023 01:24:41 +0100 Subject: [PATCH 09/24] wrapping to with_logger --- src/ODE.jl | 22 ++--- src/RationalFunctionFields/IdealMQS.jl | 32 ++++---- .../RationalFunctionField.jl | 82 ++++++++----------- src/RationalFunctionFields/normalforms.jl | 48 +++++------ src/StructuralIdentifiability.jl | 45 ++++++++-- src/discrete.jl | 28 +++---- src/elimination.jl | 42 ++++++---- src/global_identifiability.jl | 41 ++++------ src/identifiable_functions.jl | 52 ++++++++++-- src/io_equation.jl | 73 ++++++++--------- src/lincomp.jl | 4 +- src/local_identifiability.jl | 78 ++++++++++++------ src/logging.jl | 30 ------- src/parametrizations.jl | 81 ++++++++++-------- src/pb_representation.jl | 2 +- src/power_series_utils.jl | 2 +- src/primality_check.jl | 10 +-- src/util.jl | 14 ++-- src/wronskian.jl | 10 +-- 19 files changed, 365 insertions(+), 331 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index 98835b10d..52cad64d4 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -25,7 +25,7 @@ struct ODE{P} # P is the type of polynomials in the rhs of the ODE system # x_eqs is a dictionary x_i => f_i(x, u, params) # y_eqs is a dictionary y_i => g_i(x, u, params) if isempty(y_eqs) - info_si("Could not find output variables in the model.") + @info "Could not find output variables in the model." end poly_ring = parent(first(vcat(y_vars, x_vars))) u_vars = inputs @@ -460,19 +460,13 @@ macro ODEmodel(ex::Expr...) end logging_exprs = [ - :(StructuralIdentifiability.info_si("Summary of the model:")), - :(StructuralIdentifiability.info_si( - "State variables: " * $(join(map(string, collect(x_vars)), ", ")), - )), - :(StructuralIdentifiability.info_si( - "Parameters: " * $(join(map(string, collect(params)), ", ")), - )), - :(StructuralIdentifiability.info_si( - "Inputs: " * $(join(map(string, collect(u_vars)), ", ")), - )), - :(StructuralIdentifiability.info_si( - "Outputs: " * $(join(map(string, collect(y_vars)), ", ")), - )), + # :(with_logger(StructuralIdentifiability._si_logger[]) do), + :(@info "Summary of the model:"), + :(@info "State variables: " * $(join(map(string, collect(x_vars)), ", "))), + :(@info "Parameters: " * $(join(map(string, collect(params)), ", "))), + :(@info "Inputs: " * $(join(map(string, collect(u_vars)), ", "))), + :(@info "Outputs: " * $(join(map(string, collect(y_vars)), ", "))), + # :(end), ] # creating the ode object ode_expr = :(StructuralIdentifiability.ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}( diff --git a/src/RationalFunctionFields/IdealMQS.jl b/src/RationalFunctionFields/IdealMQS.jl index e1d9057dd..6d8ddb81c 100644 --- a/src/RationalFunctionFields/IdealMQS.jl +++ b/src/RationalFunctionFields/IdealMQS.jl @@ -63,11 +63,11 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal # We prepare and store them to construct ideal specializations later @assert !isempty(funcs_den_nums) @assert sat_var_position in (:first, :last) - ordering !== :degrevlex && (warn_si("Ordering is not degrevlex but $ordering")) + ordering !== :degrevlex && (@warn "Ordering is not degrevlex but $ordering") ring = parent(first(first(funcs_den_nums))) - debug_si("Constructing the MQS ideal in $ring") + @debug "Constructing the MQS ideal in $ring" K, n = base_ring(ring), nvars(ring) - debug_si("Finding pivot polynomials") + @debug "Finding pivot polynomials" # In the component f1,f2,... find the polynomial with the minimal total # degree and length. Such element will serve as a normalizing term for # the component @@ -78,19 +78,17 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal map(plist -> findmin(p -> (total_degree(p), length(p)), plist), funcs_den_nums) pivots_indices = map(last, pivots) for plist in funcs_den_nums - debug_si("\tDegrees in this list are $(map(total_degree, plist))") + @debug "\tDegrees in this list are $(map(total_degree, plist))" end - debug_si("\tDegrees and lengths are $(map(first, pivots))") + @debug "\tDegrees and lengths are $(map(first, pivots))" den_lcm = mapreduce( i -> funcs_den_nums[i][pivots_indices[i]], lcm, 1:length(funcs_den_nums), ) - debug_si( - "Rational functions common denominator is of degree $(total_degree(den_lcm)) and of length $(length(den_lcm))", - ) + @debug "Rational functions common denominator is of degree $(total_degree(den_lcm)) and of length $(length(den_lcm))" is_constant(den_lcm) && - (debug_si("Common denominator of the field generators is constant")) + (@debug "Common denominator of the field generators is constant") existing_varnames = map(String, symbols(ring)) ystrs = ["y$i" for i in 1:length(existing_varnames)] @assert !(sat_varname in ystrs) "The name of the saturation variable collided with a primary variable" @@ -101,7 +99,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal length(ystrs) + 1 end varnames = append_at_index(ystrs, sat_var_index, sat_varname) - debug_si("Saturating variable is $sat_varname, index is $sat_var_index") + @debug "Saturating variable is $sat_varname, index is $sat_var_index" R_sat, v_sat = Nemo.PolynomialRing(K, varnames, ordering = ordering) # Saturation t_sat = v_sat[sat_var_index] @@ -147,7 +145,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal end end parent_ring_param, _ = PolynomialRing(ring, varnames, ordering = ordering) - debug_si("Constructed MQS ideal in $R_sat with $(length(nums_qq) + 1) elements") + @debug "Constructed MQS ideal in $R_sat with $(length(nums_qq) + 1) elements" @assert length(pivots_indices) == length(dens_indices) == length(dens_qq) @assert length(pivots_indices) == length(funcs_den_nums) @@ -200,12 +198,12 @@ function fractionfree_generators_raw(mqs::IdealMQS) old_varnames = map(i -> "y$i", 1:length(varnames)) new_varnames = map(i -> "라$i", 1:(length(varnames) + 1)) if !isempty(intersect(old_varnames, new_varnames)) - warn_si("Intersection in two sets of variables! $varnames $new_varnames") + @warn "Intersection in two sets of variables! $varnames $new_varnames" end # NOTE: new variables go first! big_ring, big_vars = PolynomialRing(K, vcat(new_varnames, old_varnames), ordering = :lex) - info_si("$(mqs.sat_var_index) $(varnames) $ring_params $(parent(mqs.sat_qq))") + @info "$(mqs.sat_var_index) $(varnames) $ring_params $(parent(mqs.sat_qq))" nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq nums_y = map(num -> parent_ring_change(num, big_ring, matching = :byindex), nums_qq) dens_y = map(den -> parent_ring_change(den, big_ring, matching = :byindex), dens_qq) @@ -232,10 +230,10 @@ function ParamPunPam.reduce_mod_p!( mqs::IdealMQS, ff::Field, ) where {Field <: Union{Nemo.GaloisField, Nemo.GaloisFmpzField}} - debug_si("Reducing MQS ideal modulo $(ff)") + @debug "Reducing MQS ideal modulo $(ff)" # If there is a reduction modulo this field already, if haskey(mqs.cached_nums_gf, ff) - debug_si("Cache hit with $(ff)!") + @debug "Cache hit with $(ff)!" return nothing end nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq @@ -255,7 +253,7 @@ function ParamPunPam.specialize_mod_p( saturated = true, ) where {T <: Union{gfp_elem, gfp_fmpz_elem}} K_1 = parent(first(point)) - debug_si("Evaluating MQS ideal over $K_1 at $point") + @debug "Evaluating MQS ideal over $K_1 at $point" @assert haskey(mqs.cached_nums_gf, K_1) nums_gf, dens_gf, sat_gf = mqs.cached_nums_gf[K_1], mqs.cached_dens_gf[K_1], mqs.cached_sat_gf[K_1] @@ -286,7 +284,7 @@ function ParamPunPam.specialize_mod_p( end function specialize(mqs::IdealMQS, point::Vector{Nemo.fmpq}; saturated = true) - debug_si("Evaluating MQS ideal over QQ at $point") + @debug "Evaluating MQS ideal over QQ at $point" nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq dens_indices = mqs.dens_indices K = base_ring(mqs) diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index 2f8378bfa..e68f7f5f2 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -100,11 +100,11 @@ Output: if isempty(ratfuncs) return Bool[] end - debug_si("Finding pivot polynomials") + @debug "Finding pivot polynomials" pivots = map(plist -> plist[findmin(map(total_degree, plist))[2]], ratfuncs) - debug_si("\tDegrees are $(map(total_degree, pivots))") + @debug "\tDegrees are $(map(total_degree, pivots))" - debug_si("Estimating the sampling bound") + @debug "Estimating the sampling bound" # uses Theorem 3.3 from https://arxiv.org/pdf/2111.00991.pdf # the comments below use the notation from the theorem ring = parent(first(first(ratfuncs))) @@ -121,13 +121,13 @@ Output: extra_degree = total_degree(den_lcm) - total_degree(field.mqs.dens_qq[i]) degree = max(degree, extra_degree + maximum(total_degree, plist)) end - debug_si("\tBound for the degrees is $degree") + @debug "\tBound for the degrees is $degree" total_vars = foldl( union, map(plist -> foldl(union, map(poly -> Set(vars(poly)), plist)), field.dennums), ) - debug_si("\tThe total number of variables in $(length(total_vars))") + @debug "\tThe total number of variables in $(length(total_vars))" sampling_bound = BigInt( 3 * @@ -135,16 +135,16 @@ Output: (length(ratfuncs)) * ceil(1 / (1 - p)), ) - debug_si("\tSampling from $(-sampling_bound) to $(sampling_bound)") + @debug "\tSampling from $(-sampling_bound) to $(sampling_bound)" mqs = field.mqs param_ring = ParamPunPam.parent_params(mqs) point = map(v -> Nemo.QQ(rand((-sampling_bound):sampling_bound)), gens(param_ring)) mqs_specialized = specialize(mqs, point) - debug_si("Computing Groebner basis ($(length(mqs_specialized)) equations)") + @debug "Computing Groebner basis ($(length(mqs_specialized)) equations)" mqs_ratfuncs = specialize(IdealMQS(ratfuncs), point; saturated = false) @assert parent(first(mqs_specialized)) == parent(first(mqs_ratfuncs)) - debug_si("Starting the groebner basis computation") + @debug "Starting the groebner basis computation" gb = groebner(mqs_specialized) result = map(iszero, normalform(gb, mqs_ratfuncs)) return result @@ -220,13 +220,13 @@ Applies the following passes: fracs = filter(!is_rational_func_const, fracs) fracs = unique(fracs) if isempty(fracs) - debug_si("The set of generators is empty") + @debug "The set of generators is empty" return fracs end # Remove redundant pass if discard_redundant sort!(fracs, lt = rational_function_cmp) - debug_si("The pool of fractions:\n$(join(map(repr, fracs), ",\n"))") + @debug "The pool of fractions:\n$(join(map(repr, fracs), ",\n"))" if reversed_order non_redundant = collect(1:length(fracs)) for i in length(fracs):-1:1 @@ -236,9 +236,9 @@ Applies the following passes: end result = check_field_membership_mod_p(fracs[setdiff(non_redundant, i)], [func]) - debug_si("Simplification: inclusion check $func $result") + @debug "Simplification: inclusion check $func $result" if result[1] - debug_si("The function $func is discarded") + @debug "The function $func is discarded" setdiff!(non_redundant, i) end end @@ -248,16 +248,14 @@ Applies the following passes: for i in 2:length(fracs) func = fracs[i] result = check_field_membership_mod_p(fracs[non_redundant], [func]) - debug_si("Simplification: inclusion check $func $result") + @debug "Simplification: inclusion check $func $result" if !result[1] - debug_si("The function $func is included in the set of generators") + @debug "The function $func is included in the set of generators" push!(non_redundant, i) end end end - debug_si( - "Out of $(length(fracs)) simplified generators there are $(length(non_redundant)) non redundant", - ) + @debug "Out of $(length(fracs)) simplified generators there are $(length(non_redundant)) non redundant" fracs = fracs[non_redundant] end sort!(fracs, lt = rational_function_cmp) @@ -317,7 +315,7 @@ end gb, fracs, new_rff = nothing, nothing, nothing # Check if the basis is in cache if haskey(mqs.cached_groebner_bases, (ordering, up_to_degree)) - debug_si("Cache hit with ($ordering, $up_to_degree)!") + @debug "Cache hit with ($ordering, $up_to_degree)!" gb = mqs.cached_groebner_bases[ordering, up_to_degree] basis_coeffs = map(collect ∘ coefficients, gb) fracs = collect(mapreduce(Set, union!, basis_coeffs)) @@ -327,7 +325,7 @@ end current_degrees = (2, 2) two_sided_inclusion = false while !two_sided_inclusion && all(current_degrees .<= up_to_degree) - debug_si("Computing GB with parameters up to degrees $(current_degrees)") + @debug "Computing GB with parameters up to degrees $(current_degrees)" runtime = @elapsed gb = ParamPunPam.paramgb( mqs, up_to_degree = current_degrees, @@ -339,12 +337,12 @@ end _runtime_logger[:id_npoints_interpolation] += ParamPunPam._runtime_data[:npoints_interpolation] _runtime_logger[:id_groebner_time] += runtime - debug_si("Groebner basis computed in $runtime seconds") + @debug "Groebner basis computed in $runtime seconds" basis_coeffs = map(collect ∘ coefficients, gb) basis_coeffs_set = mapreduce(Set, union!, basis_coeffs) fracs = collect(basis_coeffs_set) - debug_si("Generators up to degrees $(current_degrees) are $fracs") - debug_si("Checking two-sided inclusion modulo a prime") + @debug "Generators up to degrees $(current_degrees) are $fracs" + @debug "Checking two-sided inclusion modulo a prime" time_start = time_ns() # Check inclusion: in new_rff = RationalFunctionField(fracs) @@ -355,12 +353,10 @@ end runtime = (time_ns() - time_start) / 1e9 _runtime_logger[:id_inclusion_check_mod_p] += runtime two_sided_inclusion = two_sided_inclusion && all(inclusion) - debug_si("Inclusion checked in $(runtime) seconds. Result: $two_sided_inclusion") + @debug "Inclusion checked in $(runtime) seconds. Result: $two_sided_inclusion" current_degrees = current_degrees .* 2 end - debug_si( - "The coefficients of the Groebner basis are presented by $(length(fracs)) rational functions", - ) + @debug "The coefficients of the Groebner basis are presented by $(length(fracs)) rational functions" new_rff.mqs.cached_groebner_bases[ordering, up_to_degree] = gb rff.mqs.cached_groebner_bases[ordering, up_to_degree] = gb return new_rff @@ -386,9 +382,7 @@ Returns a set of Groebner bases for multiple different rankings of variables. time_start = time_ns() vars = gens(parent(rff.mqs)) nbases = length(vars) - info_si( - "Computing $nbases Groebner bases for degrees $up_to_degree for block orderings", - ) + @info "Computing $nbases Groebner bases for degrees $up_to_degree for block orderings" ordering_to_generators = Dict() if code == 0 return ordering_to_generators @@ -414,7 +408,7 @@ Returns a set of Groebner bases for multiple different rankings of variables. # n1, n2 = div(n, 2), n - div(n, 2) n1, n2 = n - 1, 1 ord = DegRevLex(vars_shuffled[1:n1]) * DegRevLex(vars_shuffled[(n1 + 1):end]) - debug_si("Computing GB for ordering $ord") + @debug "Computing GB for ordering $ord" new_rff = groebner_basis_coeffs( gb_rff, seed = seed, @@ -431,7 +425,7 @@ Returns a set of Groebner bases for multiple different rankings of variables. n = length(vars_shuffled) n1, n2 = max(n - 2, 1), min(2, length(vars) - 1) ord = DegRevLex(vars_shuffled[1:n1]) * DegRevLex(vars_shuffled[(n1 + 1):end]) - debug_si("Computing GB for ordering $ord") + @debug "Computing GB for ordering $ord" new_rff = groebner_basis_coeffs( gb_rff, seed = seed, @@ -449,7 +443,7 @@ Returns a set of Groebner bases for multiple different rankings of variables. n1 = div(n, 2) n2 = n - n1 ord = DegRevLex(vars_shuffled[1:n1]) * DegRevLex(vars_shuffled[(n1 + 1):end]) - debug_si("Computing GB for ordering $ord") + @debug "Computing GB for ordering $ord" new_rff = groebner_basis_coeffs( gb_rff, seed = seed, @@ -461,7 +455,7 @@ Returns a set of Groebner bases for multiple different rankings of variables. end end _runtime_logger[:id_gbfan_time] = (time_ns() - time_start) / 1e9 - info_si("Computed Groebner bases in $((time_ns() - time_start) / 1e9) seconds") + @info "Computed Groebner bases in $((time_ns() - time_start) / 1e9) seconds" return ordering_to_generators end @@ -494,7 +488,7 @@ Result is correct (in Monte-Carlo sense) with probability at least `p`. check_variables = false, # almost always slows down and thus turned off rational_interpolator = :VanDerHoevenLecerf, ) - info_si("Simplifying generating set. Simplification level: $simplify") + @info "Simplifying generating set. Simplification level: $simplify" _runtime_logger[:id_groebner_time] = 0.0 _runtime_logger[:id_calls_to_gb] = 0 _runtime_logger[:id_inclusion_check_mod_p] = 0.0 @@ -558,29 +552,23 @@ Result is correct (in Monte-Carlo sense) with probability at least `p`. append!(new_fracs, generators) end new_fracs_unique = unique(new_fracs) - debug_si( - """ + @debug """ Final cleaning and simplification of generators. -Out of $(length(new_fracs)) fractions $(length(new_fracs_unique)) are syntactically unique.""", - ) +Out of $(length(new_fracs)) fractions $(length(new_fracs_unique)) are syntactically unique.""" runtime = @elapsed new_fracs = beautifuly_generators(RationalFunctionField(new_fracs_unique)) - debug_si("Checking inclusion with probability $p") + @debug "Checking inclusion with probability $p" runtime = @elapsed result = issubfield(rff, RationalFunctionField(new_fracs), p) _runtime_logger[:id_inclusion_check] = runtime if !result - warn_si("Field membership check failed. Error will follow.") + @warn "Field membership check failed. Error will follow." throw("The new subfield generators are not correct.") end - info_si( - "Inclusion checked with probability $p in $(_runtime_logger[:id_inclusion_check]) seconds", - ) - debug_si( - "Out of $(length(rff.mqs.nums_qq)) initial generators there are $(length(new_fracs)) indepdendent", - ) + @info "Inclusion checked with probability $p 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 - debug_si("The ranking of the new set of generators is $ranking") + @debug "The ranking of the new set of generators is $ranking" return new_fracs end diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index 0bf45773f..c095a00be 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -63,27 +63,27 @@ function local_normal_forms( monoms_ff = Vector{elem_type(ring_ff)}(undef, 0) xs_ff = cut_at_index(xs_ff, mqs.sat_var_index) pivot_vectors = map(f -> exponent_vector(f, 1), xs_ff) - debug_si(""" + @debug """ variables in the finite field: $(xs_ff) gb parent: $(ring_ff) specialized gb: $(gb_ff_spec) - Evaluation point: $point_ff_ext""") + Evaluation point: $point_ff_ext""" # Compute the normal forms of all monomials of degrees up to `up_to_degree` for deg in 1:up_to_degree for combination in Combinatorics.with_replacement_combinations(pivot_vectors, deg) exp_vect = sum(combination) if in(stop_vectors, exp_vect) - debug_si("Skipping exponent vector $exp_vect") + @debug "Skipping exponent vector $exp_vect" continue end monom_ff = ring_ff([one(finite_field)], [exp_vect]) monom_ff_spec = evaluate(monom_ff, point_ff_ext) monom_mqs_ff_spec = monom_ff - monom_ff_spec divisors_ff, nf_ff = divrem(monom_mqs_ff_spec, gb_ff_spec) - debug_si(""" + @debug """ The normal form of $monom_mqs_ff_spec is: normalform = $nf_ff - divisors = $divisors_ff""") + divisors = $divisors_ff""" push!(monoms_ff, monom_ff) push!(normal_forms_ff, nf_ff) end @@ -161,7 +161,7 @@ function linear_relations_over_a_field(polys, preimages) for ind in zero_inds push!(relations, preimages[ind]) end - debug_si("Zeroed polynomials are $(preimages[zero_inds])") + @debug "Zeroed polynomials are $(preimages[zero_inds])" permutation = setdiff(permutation, zero_inds) # Sort, the first monom is the smallest lead_monoms = map(f -> iszero(f) ? one(f) : leading_monomial(f), polys) @@ -303,9 +303,9 @@ is not specified but is assumed to be close to 1. nparams = nvars(ring_param) finite_field = Nemo.GF(2^30 + 3) ParamPunPam.reduce_mod_p!(mqs, finite_field) - info_si("Computing normal forms of degree $up_to_degree in $nparams variables") - debug_si("""Variables ($nparams in total): $xs_param - Modulo: $finite_field""") + @info "Computing normal forms of degree $up_to_degree in $nparams variables" + @debug """Variables ($nparams in total): $xs_param + Modulo: $finite_field""" # We first compute relations between the normal forms of linear monomials. # Then, we use this knowledge to drop out some monomials of higher degrees. tref = TinyRowEchelonForm{Int}() @@ -315,7 +315,7 @@ is not specified but is assumed to be close to 1. for i in 1:length(normal_forms_ff_1) !iszero(normal_forms_ff_1[i]) && continue !(length(monoms_ff_1[i]) == 1) && continue - debug_si("Registering existing monomial $(monoms_ff_1[i])") + @debug "Registering existing monomial $(monoms_ff_1[i])" push!(relations_ff_1, monoms_ff_1[i]) push!(tref, exponent_vector(monoms_ff_1[i], 1)) end @@ -325,17 +325,17 @@ is not specified but is assumed to be close to 1. while true iters += 1 point = ParamPunPam.distinct_nonzero_points(finite_field, nvars(ring_param)) - debug_si("Used specialization points: $iters") - debug_si("Computing normal forms to to degree $up_to_degree") + @debug "Used specialization points: $iters" + @debug "Computing normal forms to to degree $up_to_degree" gb_ff, normal_forms_ff, monoms_ff = local_normal_forms(mqs, finite_field, up_to_degree, point, stop_vectors = tref) if isempty(normal_forms_ff) break end - debug_si("Computing relations of $(length(normal_forms_ff)) normal forms") + @debug "Computing relations of $(length(normal_forms_ff)) normal forms" relations_ff, normal_forms_ff, monoms_ff = linear_relations_over_a_field(normal_forms_ff, monoms_ff) - debug_si("Obtained $(length(relations_ff)) local relations") + @debug "Obtained $(length(relations_ff)) local relations" if iters == 1 # first point in the sequence, take all relations complete_intersection_relations_ff = relations_ff @@ -365,11 +365,9 @@ is not specified but is assumed to be close to 1. push!(zeroed_relations_inds, i) end end - debug_si( - """ + @debug """ Relations in the previous intersection: $(length(complete_intersection_relations_ff)) - Vanished at the current point: $(length(zeroed_relations_inds))""", - ) + Vanished at the current point: $(length(zeroed_relations_inds))""" non_zeroed_relations_inds = setdiff(collect(1:n_relations_ff), zeroed_relations_inds) zeroed_relations_from_complete_intersection = @@ -385,16 +383,14 @@ is not specified but is assumed to be close to 1. non_zeroed_relations_from_complete_intersection, zeroed_relations_from_complete_intersection, ) - debug_si( - "There are $(length(complete_intersection_relations_ff)) relations in the intersection", - ) + @debug "There are $(length(complete_intersection_relations_ff)) relations in the intersection" m_relations_ff = length(complete_intersection_relations_ff) if n_relations_ff == m_relations_ff || isempty(complete_intersection_relations_ff) break end end union!(complete_intersection_relations_ff, relations_ff_1) - debug_si("Reconstructing relations to rationals") + @debug "Reconstructing relations to rationals" relations_qq = Vector{Generic.Frac{elem_type(ring_param)}}( undef, length(complete_intersection_relations_ff), @@ -404,10 +400,10 @@ is not specified but is assumed to be close to 1. success, relation_qq = ParamPunPam.rational_reconstruct_polynomial(ring, relation_ff) if !success - warn_si(""" + @warn """ Failed to reconstruct the $i-th relation. Error will follow. relation: $relation_ff - modulo: $finite_field""") + modulo: $finite_field""" throw(ErrorException("Rational reconstruction failed.")) end relation_qq_param = evaluate( @@ -416,8 +412,6 @@ is not specified but is assumed to be close to 1. ) relations_qq[i] = relation_qq_param // one(relation_qq_param) end - info_si( - "Used $iters specializations in $((time_ns() - time_start) / 1e9) seconds, found $(length(complete_intersection_relations_ff)) relations", - ) + @info "Used $iters specializations in $((time_ns() - time_start) / 1e9) seconds, found $(length(complete_intersection_relations_ff)) relations" relations_qq end diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index fd4db817f..ba146daf6 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -95,7 +95,20 @@ function assess_identifiability( ) 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 + ) + end +end +function _assess_identifiability( + ode::ODE{P}; + funcs_to_check = Vector(), + p::Float64 = 0.99, +) where {P <: MPolyElem{fmpq}} + p_glob = 1 - (1 - p) * 0.9 p_loc = 1 - (1 - p) * 0.1 @@ -103,18 +116,17 @@ function assess_identifiability( funcs_to_check = vcat(ode.parameters, ode.x_vars) end - info_si("Assessing local identifiability") + @info "Assessing local identifiability" trbasis = Array{fmpq_mpoly, 1}() - runtime = @elapsed local_result = assess_local_identifiability( + runtime = @elapsed local_result = _assess_local_identifiability( ode, funcs_to_check = funcs_to_check, p = p_loc, type = :SE, trbasis = trbasis, - loglevel = loglevel, ) - # info_si("Local identifiability assessed in $runtime seconds") - debug_si("Trasncendence basis to be specialized is $trbasis") + @debug "Local identifiability assessed in $runtime seconds" + @debug "Trasncendence basis to be specialized is $trbasis" # _runtime_logger[:loc_time] = runtime loc_id = [local_result[each] for each in funcs_to_check] @@ -125,10 +137,10 @@ function assess_identifiability( end end - info_si("Assessing global identifiability") + @info "Assessing global identifiability" runtime = @elapsed global_result = assess_global_identifiability(ode, locally_identifiable, trbasis, p_glob) - info_si("Global identifiability assessed in $runtime seconds") + @info "Global identifiability assessed in $runtime seconds" _runtime_logger[:glob_time] = runtime result = Dict{Any, Symbol}() @@ -167,6 +179,22 @@ function assess_identifiability( funcs_to_check = [], p = 0.99, loglevel = Logging.Info, +) + restart_logging(loglevel = loglevel) + with_logger(_si_logger[]) do + return _assess_identifiability(ode, + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + p = p, + ) + end +end + +function _assess_identifiability( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = [], + p = 0.99, ) if isempty(measured_quantities) measured_quantities = get_measured_quantities(ode) @@ -179,11 +207,10 @@ function assess_identifiability( end funcs_to_check_ = [eval_at_nemo(each, conversion) for each in funcs_to_check] - result = assess_identifiability( + result = _assess_identifiability( ode, funcs_to_check = funcs_to_check_, p = p, - loglevel = loglevel, ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) diff --git a/src/discrete.jl b/src/discrete.jl index b964570ac..628797e98 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -72,12 +72,12 @@ function differentiate_sequence_solution( inputs::Dict{P, Array{T, 1}}, num_terms::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - debug_si("Computing the power series solution of the system") + @debug "Computing the power series solution of the system" seq_sol = sequence_solution(dds, params, ic, inputs, num_terms) generalized_params = vcat(dds.x_vars, dds.parameters) bring = base_ring(dds.poly_ring) - debug_si("Solving the variational system at the solution") + @debug "Solving the variational system at the solution" part_diffs = Dict( (x, p) => derivative(dds.x_equations[x], p) for x in dds.x_vars for p in generalized_params @@ -126,10 +126,10 @@ function differentiate_sequence_output( inputs::Dict{P, Array{T, 1}}, num_terms::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - debug_si("Computing partial derivatives of the solution") + @debug "Computing partial derivatives of the solution" seq_sol, sol_diff = differentiate_sequence_solution(dds, params, ic, inputs, num_terms) - debug_si("Evaluating the partial derivatives of the outputs") + @debug "Evaluating the partial derivatives of the outputs" generalized_params = vcat(dds.x_vars, dds.parameters) part_diffs = Dict( (y, p) => derivative(dds.y_equations[y], p) for y in dds.y_vars for @@ -195,7 +195,7 @@ function _assess_local_identifiability_discrete( ) where {P <: MPolyElem{Nemo.fmpq}} bring = base_ring(dds.poly_ring) - debug_si("Extending the model") + @debug "Extending the model" dds_ext = add_outputs(dds, Dict("loc_aux_$i" => f for (i, f) in enumerate(funcs_to_check))) @@ -206,7 +206,7 @@ function _assess_local_identifiability_discrete( known_ic = dds_ext.x_vars end - debug_si("Computing the observability matrix") + @debug "Computing the observability matrix" prec = length(dds.x_vars) + length(dds.parameters) # Computing the bound from the Schwartz-Zippel-DeMilo-Lipton lemma @@ -221,7 +221,7 @@ function _assess_local_identifiability_discrete( Jac_degree += 2 * deg_y * prec end D = Int(ceil(Jac_degree * length(funcs_to_check) / (1 - p))) - debug_si("Sampling range $D") + @debug "Sampling range $D" # Parameter values are the same across all the replicas params_vals = Dict(p => bring(rand(1:D)) for p in dds_ext.parameters) @@ -231,11 +231,11 @@ function _assess_local_identifiability_discrete( u => [bring(rand(1:D)) for i in 1:prec] for u in dds_ext.u_vars ) - debug_si("Computing the output derivatives") + @debug "Computing the output derivatives" output_derivatives = differentiate_sequence_output(dds_ext, params_vals, ic, inputs, prec) - debug_si("Building the matrices") + @debug "Building the matrices" Jac = zero( Nemo.MatrixSpace( bring, @@ -265,7 +265,7 @@ function _assess_local_identifiability_discrete( end end - debug_si("Computing the result") + @debug "Computing the result" base_rank = LinearAlgebra.rank(Jac) result = Dict{Any, Bool}() for i in 1:length(funcs_to_check) @@ -314,9 +314,7 @@ function assess_local_identifiability( ) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(dds)) - info_si( - "Measured quantities are not provided, trying to find the outputs in input dynamical system.", - ) + @info "Measured quantities are not provided, trying to find the outputs in input dynamical system." measured_quantities = filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(dds), @@ -344,9 +342,7 @@ function assess_local_identifiability( nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) if length(known_ic) > 0 - warn_si( - "Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !", - ) + @warn "Since known initial conditions were provided, identifiability of states (e.g., `x(t)`) is at t = 0 only !" end return out_dict end diff --git a/src/elimination.jl b/src/elimination.jl index bab6aa9cc..6798889e9 100644 --- a/src/elimination.jl +++ b/src/elimination.jl @@ -32,7 +32,8 @@ function det_minor_expansion_inner( cache[discarded] = result end if length(discarded[1]) < 3 - debug_si("Discarded: $discarded") + @debug "Discarded: $discarded" + flush(_si_logger[].stream) end return result end @@ -201,11 +202,11 @@ function Base.iterate( i::Int = 1, ) where {P <: MPolyElem{<:FieldElem}} if i > length(gpg.cached_points) - debug_si("Generating new point on the variety") + @debug "Generating new point on the variety" sample_max = i * 50 result = undef while true - debug_si("Preparing initial condition") + @debug "Preparing initial condition" base_field = base_ring(gpg.big_ring) param_values = Dict{P, Int}(p => rand(1:sample_max) for p in gpg.ode.parameters) initial_conditions = @@ -214,7 +215,7 @@ function Base.iterate( u => [rand(1:sample_max) for _ in 1:(gpg.precision)] for u in gpg.ode.u_vars ) - debug_si("Computing a power series solution") + @debug "Computing a power series solution" ps_solution = undef try ps_solution = power_series_solution( @@ -225,11 +226,11 @@ function Base.iterate( gpg.precision, ) catch e - debug_si("$e") + @debug "$e" continue end - debug_si("Constructing the point") + @debug "Constructing the point" result = Dict{P, gpg.number_type}( switch_ring(p, gpg.big_ring) => base_field(c) for (p, c) in param_values ) @@ -316,12 +317,12 @@ Output: lg = coeff(g, [var_elim], [Nemo.degree(g, var_elim)]) (flag, q) = divides(lg, lf) if flag - debug_si("\t Decreasing degree with linear combination") + @debug "\t Decreasing degree with linear combination" g = g - q * f * var_elim^(Nemo.degree(g, var_elim) - Nemo.degree(f, var_elim)) elseif (Nemo.degree(g, var_elim) == Nemo.degree(f, var_elim)) (flag, q) = divides(lf, lg) if flag - debug_si("\t Decreasing degree with linear combination") + @debug "\t Decreasing degree with linear combination" f = f - q * g else break @@ -368,28 +369,33 @@ Output: resultant = f else if Nemo.degree(f, var_elim) > 1 - debug_si("Calculating the Bezout Matrix") + @debug "Calculating the Bezout Matrix" + flush(_si_logger[].stream) M = Bezout_matrix(f, g, var_elim) else - debug_si("Calculating the Sylvester matrix") + @debug "Calculating the Sylvester matrix" + flush(_si_logger[].stream) M = Sylvester_matrix(f, g, var_elim) end - debug_si("Simplifying the matrix") + @debug "Simplifying the matrix" + flush(_si_logger[].stream) M_simp, matrix_factors = simplify_matrix(M) - debug_si("Removed factors $(map(length, matrix_factors))") + @debug "Removed factors $(map(length, matrix_factors))" M_size = zero(Nemo.MatrixSpace(Nemo.ZZ, ncols(M_simp), ncols(M_simp))) for i in 1:ncols(M_simp) for j in 1:ncols(M_simp) M_size[i, j] = length(M_simp[i, j]) end end - debug_si("\t Matrix size: \n $M_size") - debug_si("\t Computing determinant") + @debug "\t Matrix size: \n $M_size" + @debug "\t Computing determinant" + flush(_si_logger[].stream) resultant = det_minor_expansion(M_simp) end - debug_si("Degrees are $([(v, Nemo.degree(resultant, v)) for v in vars(resultant)])") + @debug "Degrees are $([(v, Nemo.degree(resultant, v)) for v in vars(resultant)])" # Step 4: Eliminate extra factors - debug_si("Eliminating extra factors") + @debug "Eliminating extra factors" + flush(_si_logger[].stream) factors = fast_factor(resultant) for mfac in matrix_factors for fac in fast_factor(mfac) @@ -401,8 +407,8 @@ Output: res = choose(factors, generic_point_generator) for f in factors if f != res - debug_si("\t \t Size of extra factor: $(length(f))") - debug_si("\t \t It is $f") + @debug "\t \t Size of extra factor: $(length(f))" + @debug "\t \t It is $f" end end return res diff --git a/src/global_identifiability.jl b/src/global_identifiability.jl index 88f78e7a1..5f61644e0 100644 --- a/src/global_identifiability.jl +++ b/src/global_identifiability.jl @@ -27,14 +27,14 @@ are identifiable functions containing or not the state variables bring, _ = Nemo.PolynomialRing(base_ring(ode.poly_ring), varnames) if with_states - debug_si("Computing Lie derivatives") + @debug "Computing Lie derivatives" for f in states_generators(ode, io_equations) num, den = unpack_fraction(parent_ring_change(f, bring)) push!(coeff_lists[:with_states], [den, num]) end end - debug_si("Extracting coefficients") + @debug "Extracting coefficients" if !isempty(ode.parameters) nonparameters = filter( v -> !(var_to_str(v) in map(var_to_str, ode.parameters)), @@ -56,7 +56,7 @@ are identifiable functions containing or not the state variables push!(coeff_lists[:no_states], as_list) end else - debug_si("Known quantity $p cannot be casted and is thus dropped") + @debug "Known quantity $p cannot be casted and is thus dropped" end end @@ -105,34 +105,33 @@ The function returns a tuple containing the following: var_change_policy = :default, rational_interpolator = :VanDerHoevenLecerf, ) where {T} - info_si("Computing IO-equations") + @info "Computing IO-equations" ioeq_time = @elapsed io_equations = find_ioequations(ode; var_change_policy = var_change_policy) - debug_si("Sizes: $(map(length, values(io_equations)))") - info_si("Computed IO-equations in $ioeq_time seconds") + @debug "Sizes: $(map(length, values(io_equations)))" + @info "Computed IO-equations in $ioeq_time seconds" _runtime_logger[:ioeq_time] = ioeq_time if isempty(ode.parameters) - info_si("No parameters, so Wronskian computation is not needed") + @info "No parameters, so Wronskian computation is not needed" else - info_si("Computing Wronskians") + @info "Computing Wronskians" + flush(_si_logger[].stream) wrnsk_time = @elapsed wrnsk = wronskian(io_equations, ode) - info_si("Computed Wronskians in $wrnsk_time seconds") + @info "Computed Wronskians in $wrnsk_time seconds" _runtime_logger[:wrnsk_time] = wrnsk_time dims = map(ncols, wrnsk) - info_si("Dimensions of the Wronskians $dims") + @info "Dimensions of the Wronskians $dims" rank_times = @elapsed wranks = map(rank, wrnsk) - debug_si("Dimensions of the Wronskians $dims") - debug_si("Ranks of the Wronskians $wranks") - info_si("Ranks of the Wronskians computed in $rank_times seconds") + @debug "Dimensions of the Wronskians $dims" + @debug "Ranks of the Wronskians $wranks" + @info "Ranks of the Wronskians computed in $rank_times seconds" _runtime_logger[:rank_time] = rank_times if any([dim != rk + 1 for (dim, rk) in zip(dims, wranks)]) - warn_si( - "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) or transforming all the parameters to states with zero derivative", - ) + @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) or transforming all the parameters to states with zero derivative" end end @@ -144,9 +143,7 @@ The function returns a tuple containing the following: ) if with_states && !isempty(ode.parameters) - debug_si( - "Generators of identifiable functions involve states, the parameter-only part is getting simplified", - ) + @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( base_ring(bring), @@ -203,7 +200,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i for f in funcs_to_check num, den = unpack_fraction(f) if !all(v -> v in ode.parameters, union(vars(num), vars(den))) - info_si("Functions to check involve states") + @info "Functions to check involve states" states_needed = true break end @@ -309,9 +306,7 @@ Checks global identifiability of functions of parameters specified in `funcs_to_ ) where {P <: MPolyElem{fmpq}} submodels = find_submodels(ode) if length(submodels) > 0 - info_si( - "Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`", - ) + @info "Note: the input model has nontrivial submodels. If the computation for the full model will be too heavy, you may want to try to first analyze one of the submodels. They can be produced using function `find_submodels`" end result = check_identifiability( diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index c3dcdf3a4..378876880 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -53,16 +53,34 @@ function find_identifiable_functions( ) where {T <: MPolyElem{fmpq}} restart_logging(loglevel = loglevel) reset_timings() + with_logger(_si_logger[]) do + return _find_identifiable_functions(ode, + p = p, + seed = seed, + with_states = with_states, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + end +end +function _find_identifiable_functions( + ode::ODE{T}; + p::Float64 = 0.99, + seed = 42, + with_states = false, + simplify = :standard, + rational_interpolator = :VanDerHoevenLecerf, +) where {T <: MPolyElem{fmpq}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) runtime_start = time_ns() if isempty(ode.parameters) && !with_states - warn_si(""" + @warn """ There are no parameters in the given ODE, thus no identifiabile functions. Use `find_identifiable_functions` with keyword `with_states=true` to - compute the functions with the ODE states included.""") + compute the functions with the ODE states included.""" bring = parent(ode) id_funcs = [one(bring)] return id_funcs @@ -93,9 +111,7 @@ function find_identifiable_functions( id_funcs_fracs = [parent_ring_change(f, parent(ode)) for f in id_funcs_fracs] _runtime_logger[:id_total] = (time_ns() - runtime_start) / 1e9 _runtime_logger[:are_id_funcs_polynomial] = all(isone ∘ denominator, id_funcs_fracs) - info_si( - "The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds", - ) + @info "The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds" return id_funcs_fracs end @@ -146,20 +162,42 @@ function find_identifiable_functions( simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, loglevel = Logging.Info, +) + restart_logging(loglevel = loglevel) + reset_timings() + with_logger(_si_logger[]) do + return _find_identifiable_functions(ode, + measured_quantities = measured_quantities, + p = p, + seed = seed, + with_states = with_states, + simplify = simplify, + rational_interpolator = rational_interpolator, + ) + end +end + +function _find_identifiable_functions( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + p::Float64 = 0.99, + seed = 42, + with_states = false, + simplify = :standard, + rational_interpolator = :VanDerHoevenLecerf, ) Random.seed!(seed) if isempty(measured_quantities) measured_quantities = get_measured_quantities(ode) end ode, conversion = preprocess_ode(ode, measured_quantities) - result = find_identifiable_functions( + result = _find_identifiable_functions( ode, simplify = simplify, p = p, seed = seed, with_states = with_states, rational_interpolator = rational_interpolator, - loglevel = loglevel, ) result = [parent_ring_change(f, ode.poly_ring) for f in result] nemo2mtk = Dict(v => Num(k) for (k, v) in conversion) diff --git a/src/io_equation.jl b/src/io_equation.jl index 06e5ae0a4..ed1b1b1c4 100644 --- a/src/io_equation.jl +++ b/src/io_equation.jl @@ -146,7 +146,7 @@ Output: coef * evaluate(y_eq, [y], [zero(ring)]) // derivative(y_eq, y) end projection_equation, _ = unpack_fraction(projection_equation) - debug_si("Extra projection equation $projection_equation") + @debug "Extra projection equation $projection_equation" end while true @@ -157,9 +157,10 @@ Output: if isempty(var_degs) break end - debug_si("Current degrees of io-equations $var_degs") - debug_si("Orders: $y_orders") - debug_si("Sizes: $(Dict(y => length(eq) for (y, eq) in y_equations))") + @debug "Current degrees of io-equations $var_degs" + @debug "Orders: $y_orders" + @debug "Sizes: $(Dict(y => length(eq) for (y, eq) in y_equations))" + flush(_si_logger[].stream) # choosing the output to prolong outputs_with_scores = [ @@ -171,22 +172,22 @@ Output: d[1], ) for d in var_degs ] - debug_si("Scores: $outputs_with_scores") + @debug "Scores: $outputs_with_scores" y_prolong = sort(outputs_with_scores)[1][end] y_orders[y_prolong] += 1 - debug_si("Prolonging output $y_prolong") + @debug "Prolonging output $y_prolong" # Calculate the Lie derivative of the io_relation - debug_si("Prolonging") + @debug "Prolonging" next_y_equation = diff_poly(y_equations[y_prolong], derivation) for (x, eq) in x_equations - debug_si("Eliminating the derivative of $x") + @debug "Eliminating the derivative of $x" next_y_equation = eliminate_var(eq, next_y_equation, derivation[x], point_generator) end for (y, eq) in y_equations if y != y_prolong - debug_si("Eliminating the leader of the equation for $y") + @debug "Eliminating the leader of the equation for $y" # an ugly way of gettin the leader, to replace next_y_equation = eliminate_var( next_y_equation, @@ -204,7 +205,7 @@ Output: our_choice = sort(var_degs_next)[1] var_elim_deg, var_elim = our_choice[1], our_choice[3] - debug_si("Elimination of $var_elim, $(length(x_equations)) left") + @debug "Elimination of $var_elim, $(length(x_equations)) left" # Possible variable change for Axy + Bx + p(y) (x = var_elim) if auto_var_change && (var_elim_deg == 1) @@ -215,9 +216,7 @@ Output: A, B = simplify_frac(A, B) if isempty(filter!(v -> (v in keys(x_equations)), vars(A))) && (B != 0) # && (length(coeffs(A))==1) # variable change x_i' --> x_i' - (B/A)', x_i --> x_i - B/A - debug_si( - "\t Applying variable change: $(x) --> $(x) - ( $B )/( $A )", - ) + @debug "\t Applying variable change: $(x) --> $(x) - ( $B )/( $A )" dB = diff_poly(B, derivation) dA = diff_poly(A, derivation) numer_d, denom_d = simplify_frac(A * dB - dA * B, A * A) @@ -232,7 +231,7 @@ Output: generator_var_change(point_generator, x, A * x + B, A) # Change current system - debug_si("Change in the system") + @debug "Change in the system" x_equations[x] = make_substitution( x_equations[x], derivation[x], @@ -240,23 +239,23 @@ Output: denom_d, ) for xx in keys(x_equations) - debug_si("\t Change in the equation for $xx") + @debug "\t Change in the equation for $xx" x_equations[xx] = make_substitution(x_equations[xx], x, A * x - B, A) end - debug_si("Change in the outputs") + @debug "Change in the outputs" for y in keys(y_equations) - debug_si("\t Change in the output $y") + @debug "\t Change in the output $y" y_equations[y] = make_substitution(y_equations[y], x, A * x - B, A) end - debug_si("\t Change in the prolonged equation") + @debug "\t Change in the prolonged equation" next_y_equation = make_substitution(next_y_equation, x, A * x - B, A) # recalibrate system - debug_si("Unmixing the derivatives") + @debug "Unmixing the derivatives" for xx in setdiff(keys(x_equations), [x]) - debug_si("\t Unmixing $xx") + @debug "\t Unmixing $xx" x_equations[x] = eliminate_var( x_equations[x], x_equations[xx], @@ -265,10 +264,10 @@ Output: ) end # change the projection - debug_si("Change of variables in the extra projection") + @debug "Change of variables in the extra projection" projection_equation = make_substitution(projection_equation, x, A * x - B, A) - debug_si("Change of variables performed") + @debug "Change of variables performed" break end end @@ -277,16 +276,16 @@ Output: # Eliminate var_elim from the system delete!(x_equations, var_elim) - debug_si("Elimination in states") + @debug "Elimination in states" for (x, eq) in x_equations - debug_si("\t Elimination in the equation for $x") + @debug "\t Elimination in the equation for $x" x_equations[x] = eliminate_var(eq, y_equations[y_prolong], var_elim, point_generator) end - debug_si("Elimination in y_equations") + @debug "Elimination in y_equations" for y in keys(y_equations) if y != y_prolong - debug_si("Elimination in the output $y") + @debug "Elimination in the output $y" y_equations[y] = eliminate_var( y_equations[y], y_equations[y_prolong], @@ -295,14 +294,14 @@ Output: ) end end - debug_si("\t Elimination in the extra projection") + @debug "\t Elimination in the extra projection" projection_equation = eliminate_var( projection_equation, y_equations[y_prolong], var_elim, point_generator, ) - debug_si("\t Elimination in the prolonged equation") + @debug "\t Elimination in the prolonged equation" y_equations[y_prolong] = eliminate_var( y_equations[y_prolong], next_y_equation, @@ -351,33 +350,29 @@ Output: (var_change_policy == :default && length(ode.y_vars) >= 3) auto_var_change = false else - error_si("Unknown var_change policy $var_change_policy") + @error "Unknown var_change policy $var_change_policy" return end io_projections, _, _ = find_ioprojections(ode, auto_var_change, nothing) ring = parent(first(values(io_projections))) - debug_si("Check whether the original projections are enough") + @debug "Check whether the original projections are enough" if length(io_projections) == 1 || check_primality(io_projections) - debug_si( - "The projections generate an ideal with a single components of highest dimension, returning", - ) + @debug "The projections generate an ideal with a single components of highest dimension, returning" return io_projections end sampling_range = 5 while true - debug_si( - "There are several components of the highest dimension, trying to isolate one", - ) + @debug "There are several components of the highest dimension, trying to isolate one" extra_projection = sum(rand(1:sampling_range) * v for v in keys(io_projections)) - debug_si("Extra projections: $extra_projection") + @debug "Extra projections: $extra_projection" new_projections, _, projection_equation = find_ioprojections(ode, auto_var_change, extra_projection) - debug_si("Check primality") + @debug "Check primality" if check_primality(io_projections, [projection_equation]) - debug_si("Single component of highest dimension isolated, returning") + @debug "Single component of highest dimension isolated, returning" io_projections[str_to_var(PROJECTION_VARNAME, parent(projection_equation))] = projection_equation break diff --git a/src/lincomp.jl b/src/lincomp.jl index 3e394b46e..aef320884 100644 --- a/src/lincomp.jl +++ b/src/lincomp.jl @@ -29,10 +29,8 @@ function linear_compartment_model( graph::Vector{Vector{Int}}, inputs::Vector{Int}, outputs::Vector{Int}, - leaks::Vector{Int}; - loglevel = Logging.Info, + leaks::Vector{Int} ) - restart_logging(loglevel = loglevel) n = length(graph) x_vars_names = ["x$i" for i in 1:n] y_vars_names = ["y$i" for i in outputs] diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 0d4b83c8a..994a89938 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -31,14 +31,14 @@ function differentiate_solution( inputs::Dict{P, Array{T, 1}}, prec::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - debug_si("Computing the power series solution of the system") + @debug "Computing the power series solution of the system" ps_sol = power_series_solution(ode, params, ic, inputs, prec) ps_ring = parent(first(values(ps_sol))) for p in ode.parameters ps_sol[p] = ps_ring(params[p]) end - debug_si("Building the variational system at the solution") + @debug "Building the variational system at the solution" # Y' = AY + B vars = vcat(ode.x_vars, ode.parameters) SA = AbstractAlgebra.MatrixSpace(ps_ring, length(ode.x_vars), length(ode.x_vars)) @@ -60,7 +60,7 @@ function differentiate_solution( initial_condition[i, i] = 1 end - debug_si("Solving the variational system and forming the output") + @debug "Solving the variational system and forming the output" sol_var_system = ps_matrix_linear_de(A, B, initial_condition, prec) return ( ps_sol, @@ -87,14 +87,14 @@ function differentiate_output( inputs::Dict{P, Array{T, 1}}, prec::Int, ) where {T <: Generic.FieldElem, P <: MPolyElem{T}} - debug_si("Computing partial derivatives of the solution") + @debug "Computing partial derivatives of the solution" ps_sol, sol_diff = differentiate_solution(ode, params, ic, inputs, prec) ps_ring = parent(first(values(ps_sol))) for p in ode.parameters ps_sol[p] = ps_ring(params[p]) end - debug_si("Evaluating the partial derivatives of the outputs") + @debug "Evaluating the partial derivatives of the outputs" result = Dict() for (y, g) in ode.y_equations result[y] = Dict() @@ -172,11 +172,26 @@ function assess_local_identifiability( loglevel = Logging.Info, ) restart_logging(loglevel = loglevel) + with_logger(_si_logger[]) do + return _assess_local_identifiability(ode, + measured_quantities = measured_quantities, + funcs_to_check = funcs_to_check, + p = p, + type = type, + ) + end +end + +function _assess_local_identifiability( + ode::ModelingToolkit.ODESystem; + measured_quantities = Array{ModelingToolkit.Equation}[], + funcs_to_check = Array{}[], + p::Float64 = 0.99, + type = :SE, +) if length(measured_quantities) == 0 if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) - info_si( - "Measured quantities are not provided, trying to find the outputs in input ODE.", - ) + @info "Measured quantities are not provided, trying to find the outputs in input ODE." measured_quantities = filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(ode), @@ -199,23 +214,21 @@ function assess_local_identifiability( funcs_to_check_ = [eval_at_nemo(x, conversion) for x in funcs_to_check] if isequal(type, :SE) - result = assess_local_identifiability( + result = _assess_local_identifiability( ode, funcs_to_check = funcs_to_check_, p = p, type = type, - loglevel = loglevel, ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) return out_dict elseif isequal(type, :ME) - result, bd = assess_local_identifiability( + result, bd = _assess_local_identifiability( ode, funcs_to_check = funcs_to_check_, p = p, type = type, - loglevel = loglevel, ) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) @@ -246,7 +259,24 @@ function assess_local_identifiability( ) where {P <: MPolyElem{Nemo.fmpq}} restart_logging(loglevel = loglevel) reset_timings() + with_logger(_si_logger[]) do + return _assess_local_identifiability(ode, + funcs_to_check = funcs_to_check, + p = p, + type = type, + trbasis = trbasis, + ) + end +end +function _assess_local_identifiability( + ode::ODE{P}; + funcs_to_check::Array{<:Any, 1} = Array{Any, 1}(), + p::Float64 = 0.99, + type = :SE, + trbasis = nothing, +) where {P <: MPolyElem{Nemo.fmpq}} + if isempty(funcs_to_check) funcs_to_check = ode.parameters if type == :SE @@ -260,9 +290,7 @@ function assess_local_identifiability( num, den = unpack_fraction(f) for v in vcat(vars(num), vars(den)) if !(v in ode.parameters) - error_si( - "Multi-experiment identifiability is not properly defined for the states", - ) + @error "Multi-experiment identifiability is not properly defined for the states" throw(ArgumentError("State variable $v appears in $f")) end end @@ -270,7 +298,7 @@ function assess_local_identifiability( end # Computing the prime using Proposition 3.3 from https://doi.org/10.1006/jsco.2002.0532 - debug_si("Computing the prime number") + @debug "Computing the prime number" d, h = 1, 1 for f in vcat(collect(values(ode.x_equations)), collect(values(ode.y_equations))) df, hf = get_degree_and_coeffsize(f) @@ -295,17 +323,17 @@ function assess_local_identifiability( 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_si("The prime is $prime") + @debug "The prime is $prime" F = Nemo.GF(prime) - debug_si("Extending the model") + @debug "Extending the model" ode_ext = add_outputs(ode, Dict("loc_aux_$i" => f for (i, f) in enumerate(funcs_to_check))) - debug_si("Reducing the system modulo prime") + @debug "Reducing the system modulo prime" ode_red = reduce_ode_mod_p(ode_ext, prime) - debug_si("Computing the observability matrix (and, if ME, the bound)") + @debug "Computing the observability matrix (and, if ME, the bound)" prec = length(ode.x_vars) + length(ode.parameters) # Parameter values are the same across all the replicas @@ -325,10 +353,10 @@ function assess_local_identifiability( u => [F(rand(1:prime)) for i in 1:prec] for u in ode_red.u_vars ) - debug_si("Computing the output derivatives") + @debug "Computing the output derivatives" output_derivatives = differentiate_output(ode_red, params_vals, ic, inputs, prec) - debug_si("Building the matrices") + @debug "Building the matrices" newJac = vcat(Jac, zero(Nemo.MatrixSpace(F, length(ode.x_vars), ncols(Jac)))) newJac = hcat( newJac, @@ -367,7 +395,7 @@ function assess_local_identifiability( end if !isnothing(trbasis) - debug_si("Transcendence basis computation requested") + @debug "Transcendence basis computation requested" reverted_Jac = zero(Nemo.MatrixSpace(F, size(Jac)[2], size(Jac)[1])) for i in 1:size(Jac)[1] for j in 1:size(Jac)[2] @@ -407,10 +435,10 @@ function assess_local_identifiability( for i in trbasis_indices_states push!(trbasis, ode.x_vars[i]) end - debug_si("Transcendence basis $trbasis") + @debug "Transcendence basis $trbasis" end - debug_si("Computing the result") + @debug "Computing the result" base_rank = LinearAlgebra.rank(Jac) result = Dict{Any, Bool}() for i in 1:length(funcs_to_check) diff --git a/src/logging.jl b/src/logging.jl index 988144e1d..5d58bcf56 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -3,8 +3,6 @@ # - si_logger (using Logging package) for logs # - runtime_logger dictionary for recording runtime info # -# Defines functions to refresh both + specific logging functions -# info_si, debug_si, warn_si, and error_si calling @info, @debug, @warn, and @error, respectively # not the exhaustive list, just the ones to be initialized const _runtime_rubrics = ( @@ -35,34 +33,6 @@ function restart_logging(; loglevel = Logging.Info) end end -function info_si(s) - with_logger(_si_logger[]) do - @info s - flush(stdout) - end -end - -function debug_si(s) - with_logger(_si_logger[]) do - @debug s - flush(stdout) - end -end - -function warn_si(s) - with_logger(_si_logger[]) do - @warn s - flush(stdout) - end -end - -function error_si(s) - with_logger(_si_logger[]) do - @error s - flush(stdout) - end -end - ### # Timings for StructuralIdentifiability.jl. # Uses the the @timeit macro from TimerOutputs.jl. diff --git a/src/parametrizations.jl b/src/parametrizations.jl index b1aa527ba..56311fd60 100644 --- a/src/parametrizations.jl +++ b/src/parametrizations.jl @@ -53,35 +53,35 @@ function check_constructive_field_membership( end # Compute the remainders module the MQS. # These belong to K(T)[x]. - debug_si(""" - Reducing $tagged_num, $tagged_den""") + @debug """ + Reducing $tagged_num, $tagged_den""" _, num_rem = divrem(tagged_num, tagged_mqs) _, den_rem = divrem(tagged_den, tagged_mqs) - debug_si(""" + @debug """ Normal forms modulo MQS: Num: $(num_rem) - Den: $(den_rem)""") + Den: $(den_rem)""" common_factor = gcd(num_rem, den_rem) num_rem = divexact(num_rem, common_factor) den_rem = divexact(den_rem, common_factor) # If norm_form(Num) // norm_form(Den) does not belongs to K(T), then # the fraction does not belong to the field if iszero(den_rem) - warn_si(""" + @warn """ The element $tagged_num // $tagged_den is not in the sub-field Normal form, numerator: $num_rem Normal form, denominator: $den_rem Common factor: $(common_factor) - """) + """ return false, zero(ring_of_tags) // one(ring_of_tags) end if total_degree(num_rem) > 0 || total_degree(den_rem) > 0 - warn_si(""" + @warn """ The element $tagged_num // $tagged_den is not in the sub-field Normal form, numerator: $num_rem Normal form, denominator: $den_rem Common factor: $(common_factor) - """) + """ return false, zero(ring_of_tags) // one(ring_of_tags) end # Now we know that the remainders belong to K(T). @@ -95,11 +95,11 @@ function check_constructive_field_membership( _, num_den_factored = divrem(num_den, tag_relations) num_factored = num_num_factored // denominator(num) den_factored = num_den_factored // denominator(den) - debug_si(""" + @debug """ After factoring out relations: Num: $(num_factored) Den: $(den_factored) - """) + """ if iszero(den_factored) return false, zero(ring_of_tags) // one(ring_of_tags) end @@ -168,14 +168,12 @@ function check_constructive_field_membership( gen_tag_names(length(fracs_gen), "Tag") end sat_string = gen_tag_name("Sat") - debug_si( - """ + @debug """ Tags: $(join(map(x -> string(x[1]) * " -> " * string(x[2]), zip(fracs_gen, tag_strings)), "\t\n")) Saturation tag: $sat_string -""", - ) +""" poly_ring_tag, vars_tag = PolynomialRing(K, vcat(sat_string, orig_strings, tag_strings)) sat_var = vars_tag[1] orig_vars = vars_tag[2:(nvars(poly_ring) + 1)] @@ -200,11 +198,11 @@ $sat_string # - the GB of the MQS in K(T)[x][t]. # ord = Lex() ord = DegRevLex([sat_var]) * DegRevLex(orig_vars) * DegRevLex(tag_vars) - debug_si(""" + @debug """ Tagged MQS ideal: $tagged_mqs Monom ordering: - $(ord)""") + $(ord)""" tagged_mqs_gb = groebner(tagged_mqs, ordering = ord, homogenize = :no) # Relations between tags in K[T] relations_between_tags = filter( @@ -214,12 +212,12 @@ $sat_string # The basis in K[T][x] tagged_mqs_gb = setdiff(tagged_mqs_gb, relations_between_tags) tagged_mqs_gb = filter(poly -> isempty(intersect(vars(poly), [sat_var])), tagged_mqs_gb) - debug_si(""" + @debug """ Tagged MQS GB: $tagged_mqs_gb Relations between tags: $relations_between_tags - """) + """ # Reduce the fractions with respect to the MQS ideal. # # NOTE: reduction actually happens in K(T)[x]. So we map polynomials to the @@ -227,12 +225,10 @@ $sat_string ring_of_tags, tags = PolynomialRing(K, tag_strings) tag_to_gen = Dict(tags[i] => fracs_gen[i] for i in 1:length(fracs_gen)) if !isempty(intersect(tag_strings, orig_strings)) - warn_si( - """ + @warn """ There is an intersection between the names of the tag variables and the original variables. Tags: $tag_strings - Original vars: $orig_strings""", - ) + Original vars: $orig_strings""" end parametric_ring, _ = PolynomialRing(FractionField(ring_of_tags), orig_strings, ordering = :degrevlex) @@ -242,18 +238,18 @@ $sat_string Dict(gens(poly_ring_tag)[2:(nvars(poly_ring) + 1)] .=> gens(parametric_ring)), Dict(gens(poly_ring_tag)[(nvars(poly_ring) + 2):end] .=> gens(ring_of_tags)), ) - debug_si(""" + @debug """ Variable mapping: $param_var_mapping Parametric ring: $parametric_ring - """) + """ tagged_mqs_gb_param = map( poly -> crude_parent_ring_change(poly, parametric_ring, param_var_mapping), tagged_mqs_gb, ) tagged_mqs_gb_param = map(f -> divexact(f, leading_coefficient(f)), tagged_mqs_gb_param) - debug_si("Tagged parametric mqs: $tagged_mqs_gb_param") + @debug "Tagged parametric mqs: $tagged_mqs_gb_param" # Reduce each fraction var_mapping = Dict(gens(poly_ring) .=> gens(parametric_ring)) memberships = Vector{Bool}(undef, length(to_be_reduced)) @@ -307,7 +303,7 @@ function reparametrize_with_respect_to(ode, new_states, new_params) # Compute the new dynamics in terms of the original variables. # Paying attenton to the order.. new_vector_field = vector_field_along(ode.x_equations, new_states) - debug_si("New vector field:\n$new_vector_field") + @debug "New vector field:\n$new_vector_field" new_states = collect(keys(new_vector_field)) new_dynamics = [new_vector_field[new_state] for new_state in new_states] # Express the new dynamics in terms of new states and new parameters. @@ -326,14 +322,14 @@ function reparametrize_with_respect_to(ode, new_states, new_params) gen_tag_names(length(ode.u_vars), "Input"), gen_tag_names(length(ode.y_vars), "Output"), ) - info_si(""" + @info """ Tag names: $tag_names Generating functions: $generating_funcs Functions to be reduced: $to_be_reduced_funcs - """) + """ membership, new_dynamics_all, implicit_relations, new_vars = check_constructive_field_membership( RationalFunctionField(generating_funcs), @@ -355,20 +351,20 @@ function reparametrize_with_respect_to(ode, new_states, new_params) if !isempty(ode.u_vars) new_inputs = tag_inputs end - info_si(""" + @info """ New state dynamics: $new_dynamics_states New output dynamics: $new_dynamics_outputs New inputs: - $new_inputs""") + $new_inputs""" # Construct the new vector field. new_vars_vector_field = empty(ode.x_equations) for i in 1:length(new_states) state = tags[i] new_vars_vector_field[state] = new_dynamics_states[i] end - info_si("Converting variable names to human-readable ones") + @info "Converting variable names to human-readable ones" internal_variable_names = map(i -> "X$i", 1:length(new_states)) parameter_variable_names = map(i -> "a$i", 1:length(new_params)) input_variable_names = map(i -> "u$i", 1:length(tag_inputs)) @@ -475,26 +471,39 @@ function reparametrize_global( loglevel = Logging.Info, ) where {P} restart_logging(loglevel = loglevel) + with_logger(_si_logger[]) do + return _reparametrize_global(ode, + p = p, seed = seed, + ) + end +end + + +function _reparametrize_global( + ode::ODE{P}; + p = 0.99, + seed = 42, +) where {P} Random.seed!(seed) id_funcs = find_identifiable_functions(ode, with_states = true, simplify = :strong, p = p) ode_ring = parent(ode) @assert base_ring(parent(first(id_funcs))) == ode_ring - info_si("Constructing a new parametrization") + @info "Constructing a new parametrization" contains_states(poly::MPolyElem) = any(x -> degree(poly, x) > 0, ode.x_vars) contains_states(func) = contains_states(numerator(func)) || contains_states(denominator(func)) id_funcs_contains_states = filter(contains_states, id_funcs) - info_si(""" + @info """ Original states: $(ode.x_vars) Original params: $(ode.parameters) - Identifiable and contain states: $(id_funcs_contains_states)""") + Identifiable and contain states: $(id_funcs_contains_states)""" new_states = id_funcs_contains_states new_params = setdiff(id_funcs, id_funcs_contains_states) - info_si(""" + @info """ Reparametrizing with respect to: New states: $new_states - New params: $new_params""") + New params: $new_params""" new_vector_field, new_inputs, new_outputs, new_vars, implicit_relations = reparametrize_with_respect_to(ode, new_states, new_params) new_ring = parent(first(keys(new_vector_field))) diff --git a/src/pb_representation.jl b/src/pb_representation.jl index d4d3cc0fa..3594b143a 100644 --- a/src/pb_representation.jl +++ b/src/pb_representation.jl @@ -91,7 +91,7 @@ function common_ring(poly::MPolyElem, pbr::PBRepresentation) if d != nothing max_ords[d[1]] = max(d[2], max_ords[d[1]]) elseif !(var_to_str(v) in pbr.param_names) - warn_si("New variable $(var_to_str(v)), treating as a scalar parameter") + @warn "New variable $(var_to_str(v)), treating as a scalar parameter" push!(new_params, var_to_str(v)) end end diff --git a/src/power_series_utils.jl b/src/power_series_utils.jl index 276bf3e78..3124d47c0 100644 --- a/src/power_series_utils.jl +++ b/src/power_series_utils.jl @@ -261,7 +261,7 @@ function ps_ode_solution( cur_prec = 1 while cur_prec < prec - debug_si("\t Computing power series solution, currently at precision $cur_prec") + @debug "\t Computing power series solution, currently at precision $cur_prec" new_prec = min(prec, 2 * cur_prec) for i in 1:length(x_vars) set_precision!(solution[x_vars[i]], new_prec) diff --git a/src/primality_check.jl b/src/primality_check.jl index 253e579cc..e5655fd19 100644 --- a/src/primality_check.jl +++ b/src/primality_check.jl @@ -6,8 +6,8 @@ function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) dim = length(basis) S = Nemo.MatrixSpace(Nemo.QQ, dim, dim) matrices = [] - debug_si("$J $basis") - debug_si("Dim is $dim") + @debug "$J $basis" + @debug "Dim is $dim" for v in gens(parent(first(J))) M = zero(S) for (i, vec) in enumerate(basis) @@ -17,13 +17,13 @@ function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) end end push!(matrices, M) - debug_si("Multiplication by $v: $M") + @debug "Multiplication by $v: $M" end generic_multiplication = sum(Nemo.QQ(rand(1:100)) * M for M in matrices) - debug_si("$generic_multiplication") + @debug generic_multiplication R, t = Nemo.PolynomialRing(Nemo.QQ, "t") - debug_si("$(Nemo.charpoly(R, generic_multiplication))") + @debug "$(Nemo.charpoly(R, generic_multiplication))" return Nemo.isirreducible(Nemo.charpoly(R, generic_multiplication)) end diff --git a/src/util.jl b/src/util.jl index e51b3e313..a55d5571c 100644 --- a/src/util.jl +++ b/src/util.jl @@ -167,11 +167,11 @@ function make_substitution( d = Nemo.degree(f, var_sub) result = 0 - debug_si("Substitution in a polynomial of degree $d") + @debug "Substitution in a polynomial of degree $d" for i in 0:d - debug_si("\t Degree $i") + @debug "\t Degree $i" result += coeff(f, [var_sub], [i]) * (val_numer^i) * (val_denom^(d - i)) - debug_si("\t Intermediate result of size $(length(result))") + @debug "\t Intermediate result of size $(length(result))" end return result end @@ -326,7 +326,7 @@ function uncertain_factorization(f::MPolyElem{fmpq}) if degree(factor_out, main_var) != degree(gcd(f_uni, derivative(f_uni))) continue end - debug_si("Nonsquarefree poly, dividing by $factor_out") + @debug "Nonsquarefree poly, dividing by $factor_out" f = divexact(f, factor_out) f_uni = divexact(f_uni, gcd(f_uni, derivative(f_uni))) end @@ -521,7 +521,7 @@ function eval_at_nemo(e::Union{Float16, Float32, Float64}, vals::Dict) else out = rationalize(e) end - warn_si("Floating point value $e will be converted to $(out).") + @warn "Floating point value $e will be converted to $(out)." return out end @@ -566,9 +566,7 @@ end function get_measured_quantities(ode::ModelingToolkit.ODESystem) if any(ModelingToolkit.isoutput(eq.lhs) for eq in ModelingToolkit.equations(ode)) - info_si( - "Measured quantities are not provided, trying to find the outputs in input ODE.", - ) + @info "Measured quantities are not provided, trying to find the outputs in input ODE." return filter( eq -> (ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(ode), diff --git a/src/wronskian.jl b/src/wronskian.jl index c36205449..ac5750a12 100644 --- a/src/wronskian.jl +++ b/src/wronskian.jl @@ -198,9 +198,9 @@ Output: Computes the Wronskians of io_equations """ @timeit _to function wronskian(io_equations::Dict{P, P}, ode::ODE{P}) where {P <: MPolyElem} - debug_si("Compressing monomials") + @debug "Compressing monomials" termlists = [monomial_compress(ioeq, ode)[2] for ioeq in values(io_equations)] - debug_si("Matrix sizes $(map(length, termlists))") + @debug "Matrix sizes $(map(length, termlists))" # estimating the required order of truncation ord = max(map(length, termlists)...) + length(ode.x_vars) @@ -214,7 +214,7 @@ Computes the Wronskians of io_equations [map(p -> parent_ring_change(p, polyring_red), tlist) for tlist in termlists] ode_red = reduce_ode_mod_p(ode, PRIME) - debug_si("Computing power series solution up to order $ord") + @debug "Computing power series solution up to order $ord" ps = power_series_solution( ode_red, Dict(p => rand(1:100) for p in ode_red.parameters), @@ -222,7 +222,7 @@ Computes the Wronskians of io_equations Dict(u => [rand(1:100) for i in 0:ord] for u in ode_red.u_vars), ord, ) - debug_si("Computing the derivatives of the solution") + @debug "Computing the derivatives of the solution" ps_ext = Dict{MPolyElem, Nemo.gfp_abs_series}()# Generic.AbsSeries}() for v in vcat(ode_red.y_vars, ode_red.u_vars) cur_ps = ps[v] @@ -232,7 +232,7 @@ Computes the Wronskians of io_equations end end - debug_si("Constructing Wronskians") + @debug "Constructing Wronskians" result = [] for (i, tlist) in enumerate(termlists) n = length(tlist) From 93dbef9ee9c56329b1db26159bce8205cfdb7ff5 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Sun, 12 Nov 2023 14:34:48 +0300 Subject: [PATCH 10/24] adjust the groebner.jl logging level --- .../RationalFunctionField.jl | 8 ++++---- src/RationalFunctionFields/normalforms.jl | 2 +- src/logging.jl | 9 +++++++++ src/parametrizations.jl | 18 ++++++++---------- src/primality_check.jl | 6 +++--- 5 files changed, 25 insertions(+), 18 deletions(-) diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index e68f7f5f2..9c8e68fbf 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -145,8 +145,8 @@ Output: mqs_ratfuncs = specialize(IdealMQS(ratfuncs), point; saturated = false) @assert parent(first(mqs_specialized)) == parent(first(mqs_ratfuncs)) @debug "Starting the groebner basis computation" - gb = groebner(mqs_specialized) - result = map(iszero, normalform(gb, mqs_ratfuncs)) + gb = groebner(mqs_specialized, loglevel = _groebner_loglevel[]) + result = map(iszero, normalform(gb, mqs_ratfuncs, loglevel = _groebner_loglevel[])) return result end @@ -191,8 +191,8 @@ end polys_specialized = ParamPunPam.specialize_mod_p(mqs_tobereduced, point, saturated = false) @assert parent(first(gens_specialized)) == parent(first(polys_specialized)) - gb = groebner(gens_specialized) - nf = normalform(gb, polys_specialized) + gb = groebner(gens_specialized, loglevel = _groebner_loglevel[]) + nf = normalform(gb, polys_specialized, loglevel = _groebner_loglevel[]) result = map(iszero, nf) return result end diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index c095a00be..80c6329c0 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -56,7 +56,7 @@ function local_normal_forms( @assert parent(first(point)) == finite_field point_ff_ext = append_at_index(point, mqs.sat_var_index, one(finite_field)) gens_ff_spec = specialize_mod_p(mqs, point) - gb_ff_spec = Groebner.groebner(gens_ff_spec) + gb_ff_spec = Groebner.groebner(gens_ff_spec, loglevel = _groebner_loglevel[]) ring_ff = parent(gb_ff_spec[1]) xs_ff = gens(ring_ff) normal_forms_ff = Vector{elem_type(ring_ff)}(undef, 0) diff --git a/src/logging.jl b/src/logging.jl index 5d58bcf56..7e8bdbe41 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -25,12 +25,21 @@ const _runtime_logger = Dict( const _si_logger = Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info, show_limited = false)) +const _groebner_loglevel = Ref{Int}(0) function restart_logging(; loglevel = Logging.Info) _si_logger[] = Logging.ConsoleLogger(loglevel, show_limited = false) for r in _runtime_rubrics _runtime_logger[r] = 0 end + if loglevel < Logging.Info + _groebner_loglevel[] = -2 + elseif loglevel < Logging.Warn + _groebner_loglevel[] = 0 + else + _groebner_loglevel[] = 10 + end + return nothing end ### diff --git a/src/parametrizations.jl b/src/parametrizations.jl index 56311fd60..c712b3328 100644 --- a/src/parametrizations.jl +++ b/src/parametrizations.jl @@ -203,7 +203,12 @@ $sat_string $tagged_mqs Monom ordering: $(ord)""" - tagged_mqs_gb = groebner(tagged_mqs, ordering = ord, homogenize = :no) + tagged_mqs_gb = groebner( + tagged_mqs, + ordering = ord, + homogenize = :no, + loglevel = _groebner_loglevel[], + ) # Relations between tags in K[T] relations_between_tags = filter( poly -> isempty(intersect(vars(poly), vcat(sat_var, orig_vars))), @@ -472,18 +477,11 @@ function reparametrize_global( ) where {P} restart_logging(loglevel = loglevel) with_logger(_si_logger[]) do - return _reparametrize_global(ode, - p = p, seed = seed, - ) + return _reparametrize_global(ode, p = p, seed = seed) end end - -function _reparametrize_global( - ode::ODE{P}; - p = 0.99, - seed = 42, -) where {P} +function _reparametrize_global(ode::ODE{P}; p = 0.99, seed = 42) where {P} Random.seed!(seed) id_funcs = find_identifiable_functions(ode, with_states = true, simplify = :strong, p = p) diff --git a/src/primality_check.jl b/src/primality_check.jl index e5655fd19..ae98f006b 100644 --- a/src/primality_check.jl +++ b/src/primality_check.jl @@ -1,8 +1,8 @@ # ------------------------------------------------------------------------------ function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) - J = Groebner.groebner(J) - basis = Groebner.kbase(J) + J = Groebner.groebner(J, loglevel = _groebner_loglevel[]) + basis = Groebner.kbase(J, loglevel = _groebner_loglevel[]) dim = length(basis) S = Nemo.MatrixSpace(Nemo.QQ, dim, dim) matrices = [] @@ -11,7 +11,7 @@ function check_primality_zerodim(J::Array{fmpq_mpoly, 1}) for v in gens(parent(first(J))) M = zero(S) for (i, vec) in enumerate(basis) - image = Groebner.normalform(J, v * vec) + image = Groebner.normalform(J, v * vec, loglevel = _groebner_loglevel[]) for (j, base_vec) in enumerate(basis) M[i, j] = Nemo.QQ(coeff(image, base_vec)) end From cd54b1dd7726dd60c1c32bf55edeca66578f064d Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Sun, 12 Nov 2023 14:38:15 +0300 Subject: [PATCH 11/24] format --- benchmarking/IdentifiableFunctions/experiments.jl | 4 ++-- src/ODE.jl | 4 ++-- src/StructuralIdentifiability.jl | 15 ++++----------- src/identifiable_functions.jl | 10 ++++++---- src/lincomp.jl | 2 +- src/local_identifiability.jl | 7 ++++--- 6 files changed, 19 insertions(+), 23 deletions(-) diff --git a/benchmarking/IdentifiableFunctions/experiments.jl b/benchmarking/IdentifiableFunctions/experiments.jl index e63eb2124..f5b3077a9 100644 --- a/benchmarking/IdentifiableFunctions/experiments.jl +++ b/benchmarking/IdentifiableFunctions/experiments.jl @@ -995,8 +995,8 @@ begin end begin - StructuralIdentifiability.find_identifiable_functions(Pivastatin); - StructuralIdentifiability.print_timings_table(); + StructuralIdentifiability.find_identifiable_functions(Pivastatin) + StructuralIdentifiability.print_timings_table() end fracs = StructuralIdentifiability.dennums_to_fractions(dennums); diff --git a/src/ODE.jl b/src/ODE.jl index 52cad64d4..a69b59308 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -460,13 +460,13 @@ macro ODEmodel(ex::Expr...) end logging_exprs = [ - # :(with_logger(StructuralIdentifiability._si_logger[]) do), + # :(with_logger(StructuralIdentifiability._si_logger[]) do), :(@info "Summary of the model:"), :(@info "State variables: " * $(join(map(string, collect(x_vars)), ", "))), :(@info "Parameters: " * $(join(map(string, collect(params)), ", "))), :(@info "Inputs: " * $(join(map(string, collect(u_vars)), ", "))), :(@info "Outputs: " * $(join(map(string, collect(y_vars)), ", "))), - # :(end), + # :(end), ] # creating the ode object ode_expr = :(StructuralIdentifiability.ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}( diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index ba146daf6..f2eb868a6 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -96,10 +96,7 @@ function assess_identifiability( 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, p = p) end end @@ -108,7 +105,6 @@ function _assess_identifiability( funcs_to_check = Vector(), p::Float64 = 0.99, ) where {P <: MPolyElem{fmpq}} - p_glob = 1 - (1 - p) * 0.9 p_loc = 1 - (1 - p) * 0.1 @@ -182,7 +178,8 @@ function assess_identifiability( ) restart_logging(loglevel = loglevel) with_logger(_si_logger[]) do - return _assess_identifiability(ode, + return _assess_identifiability( + ode, measured_quantities = measured_quantities, funcs_to_check = funcs_to_check, p = p, @@ -207,11 +204,7 @@ 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_, p = p) nemo2mtk = Dict(funcs_to_check_ .=> funcs_to_check) out_dict = Dict(nemo2mtk[param] => result[param] for param in funcs_to_check_) return out_dict diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 378876880..52a6782c7 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -54,7 +54,8 @@ function find_identifiable_functions( restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do - return _find_identifiable_functions(ode, + return _find_identifiable_functions( + ode, p = p, seed = seed, with_states = with_states, @@ -71,7 +72,7 @@ function _find_identifiable_functions( with_states = false, simplify = :standard, rational_interpolator = :VanDerHoevenLecerf, -) where {T <: MPolyElem{fmpq}} +) where {T <: MPolyElem{fmpq}} Random.seed!(seed) @assert simplify in (:standard, :weak, :strong, :absent) runtime_start = time_ns() @@ -112,7 +113,7 @@ function _find_identifiable_functions( _runtime_logger[:id_total] = (time_ns() - runtime_start) / 1e9 _runtime_logger[:are_id_funcs_polynomial] = all(isone ∘ denominator, id_funcs_fracs) @info "The search for identifiable functions concluded in $(_runtime_logger[:id_total]) seconds" - + return id_funcs_fracs end @@ -166,7 +167,8 @@ function find_identifiable_functions( restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do - return _find_identifiable_functions(ode, + return _find_identifiable_functions( + ode, measured_quantities = measured_quantities, p = p, seed = seed, diff --git a/src/lincomp.jl b/src/lincomp.jl index aef320884..72dd6a327 100644 --- a/src/lincomp.jl +++ b/src/lincomp.jl @@ -29,7 +29,7 @@ function linear_compartment_model( graph::Vector{Vector{Int}}, inputs::Vector{Int}, outputs::Vector{Int}, - leaks::Vector{Int} + leaks::Vector{Int}, ) n = length(graph) x_vars_names = ["x$i" for i in 1:n] diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 994a89938..578240d96 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -173,7 +173,8 @@ function assess_local_identifiability( ) restart_logging(loglevel = loglevel) with_logger(_si_logger[]) do - return _assess_local_identifiability(ode, + return _assess_local_identifiability( + ode, measured_quantities = measured_quantities, funcs_to_check = funcs_to_check, p = p, @@ -260,7 +261,8 @@ function assess_local_identifiability( restart_logging(loglevel = loglevel) reset_timings() with_logger(_si_logger[]) do - return _assess_local_identifiability(ode, + return _assess_local_identifiability( + ode, funcs_to_check = funcs_to_check, p = p, type = type, @@ -276,7 +278,6 @@ function _assess_local_identifiability( type = :SE, trbasis = nothing, ) where {P <: MPolyElem{Nemo.fmpq}} - if isempty(funcs_to_check) funcs_to_check = ode.parameters if type == :SE From aded2d66fa34168d0beb2c46eeaa0e70d3e27c6d Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Sun, 12 Nov 2023 14:34:57 +0100 Subject: [PATCH 12/24] fixing old julia versions --- src/logging.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/logging.jl b/src/logging.jl index 7e8bdbe41..e6c2500ef 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -23,8 +23,12 @@ const _runtime_logger = Dict( :id_beautifulization => 0, ) -const _si_logger = +const _si_logger = @static if VERSION >= v"1.7.0" Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info, show_limited = false)) +else + Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info)) +end + const _groebner_loglevel = Ref{Int}(0) function restart_logging(; loglevel = Logging.Info) From c54df7f3d3f8d32600e79857ca26fa34f3fc9afc Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Sun, 12 Nov 2023 14:54:59 +0100 Subject: [PATCH 13/24] trying to fix for old Julia --- src/logging.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/logging.jl b/src/logging.jl index e6c2500ef..2879341f6 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -26,7 +26,7 @@ const _runtime_logger = Dict( const _si_logger = @static if VERSION >= v"1.7.0" Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info, show_limited = false)) else - Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info)) + Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger()) end const _groebner_loglevel = Ref{Int}(0) From 2c11a6e5f2a9baa9802c6aed81ff4873f20daa70 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Sun, 12 Nov 2023 15:23:01 +0100 Subject: [PATCH 14/24] v tretiy raz zabrosil starik nevod --- src/logging.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/logging.jl b/src/logging.jl index 2879341f6..79ecb2ff2 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -32,7 +32,11 @@ end const _groebner_loglevel = Ref{Int}(0) function restart_logging(; loglevel = Logging.Info) - _si_logger[] = Logging.ConsoleLogger(loglevel, show_limited = false) + _si_logger[] = @static if VERSION >= v"1.7.0" + Logging.ConsoleLogger(loglevel, show_limited = false) + else + Logging.ConsoleLogger() + end for r in _runtime_rubrics _runtime_logger[r] = 0 end From 769d8ee2b8bffc0862e6cc9613be5445395c1994 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Sun, 12 Nov 2023 18:00:06 +0300 Subject: [PATCH 15/24] fix timeit issue --- src/local_identifiability.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 578240d96..83a1396fc 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -183,7 +183,7 @@ function assess_local_identifiability( end end -function _assess_local_identifiability( +@timeit _to function _assess_local_identifiability( ode::ModelingToolkit.ODESystem; measured_quantities = Array{ModelingToolkit.Equation}[], funcs_to_check = Array{}[], @@ -249,8 +249,6 @@ Call this function if you have a specific collection of parameters of which you If the type is `:ME`, states are not allowed to appear in the `funcs_to_check`. """ function assess_local_identifiability( - # TODO: report this bug (?) to TimerOutputs.jl - # @timeit _to function assess_local_identifiability( ode::ODE{P}; funcs_to_check::Array{<:Any, 1} = Array{Any, 1}(), p::Float64 = 0.99, From 9b4aaf39d67d84fa92c2ece418ce42ab6da4d611 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Sun, 12 Nov 2023 21:19:12 +0100 Subject: [PATCH 16/24] yet another attempt on 1.6 --- src/logging.jl | 4 ++-- src/precompile.jl | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/logging.jl b/src/logging.jl index 79ecb2ff2..0056a1ecd 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -26,7 +26,7 @@ const _runtime_logger = Dict( const _si_logger = @static if VERSION >= v"1.7.0" Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info, show_limited = false)) else - Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger()) + Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(stderr, Logging.Info)) end const _groebner_loglevel = Ref{Int}(0) @@ -35,7 +35,7 @@ function restart_logging(; loglevel = Logging.Info) _si_logger[] = @static if VERSION >= v"1.7.0" Logging.ConsoleLogger(loglevel, show_limited = false) else - Logging.ConsoleLogger() + Logging.ConsoleLogger(stderr, loglevel) end for r in _runtime_rubrics _runtime_logger[r] = 0 diff --git a/src/precompile.jl b/src/precompile.jl index 9b20445ab..fb2c2094f 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -10,10 +10,10 @@ eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1] de = ODESystem(eqs, t, name = :Test) @compile_workload begin - with_logger(Logging.ConsoleLogger(Logging.Warn)) do + restart_logging(loglevel = Logging.Warn) + with_logger(_si_logger[]) do # all calls in this block will be precompiled, regardless of whether # they belong to your package or not (on Julia 1.8 and higher) - restart_logging(loglevel = Logging.Warn) ode = @ODEmodel( x1'(t) = -(a01 + a21) * x1(t) + a12 * x2(t) + u(t), x2'(t) = a21 * x1(t) - a12 * x2(t) - x3(t) / b, @@ -28,7 +28,7 @@ loglevel = Logging.Warn, ) find_identifiable_functions(ode, with_states = true, loglevel = Logging.Warn) - restart_logging(loglevel = Logging.Info) end + restart_logging(loglevel = Logging.Info) end end From 72633290af4be6a87461fbf0f75810d7a5a0c8f0 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Mon, 13 Nov 2023 00:40:07 +0300 Subject: [PATCH 17/24] run CI --- src/logging.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/logging.jl b/src/logging.jl index 0056a1ecd..7f993ba1e 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -26,12 +26,13 @@ const _runtime_logger = Dict( const _si_logger = @static if VERSION >= v"1.7.0" Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info, show_limited = false)) else - Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(stderr, Logging.Info)) + Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info)) end const _groebner_loglevel = Ref{Int}(0) function restart_logging(; loglevel = Logging.Info) + @assert loglevel isa Base.CoreLogging.LogLevel _si_logger[] = @static if VERSION >= v"1.7.0" Logging.ConsoleLogger(loglevel, show_limited = false) else From 0bd7609a142348fae845cb06754c3ef9635eaca2 Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Mon, 13 Nov 2023 00:49:34 +0300 Subject: [PATCH 18/24] revert previous commit; try to fix @ODEmodel logging --- src/ODE.jl | 18 +++++++++--------- src/logging.jl | 2 +- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/ODE.jl b/src/ODE.jl index a69b59308..c7d15afd7 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -459,15 +459,15 @@ macro ODEmodel(ex::Expr...) end end - logging_exprs = [ - # :(with_logger(StructuralIdentifiability._si_logger[]) do), - :(@info "Summary of the model:"), - :(@info "State variables: " * $(join(map(string, collect(x_vars)), ", "))), - :(@info "Parameters: " * $(join(map(string, collect(params)), ", "))), - :(@info "Inputs: " * $(join(map(string, collect(u_vars)), ", "))), - :(@info "Outputs: " * $(join(map(string, collect(y_vars)), ", "))), - # :(end), - ] + logging_exprs = [:( + with_logger(StructuralIdentifiability._si_logger[]) do + @info "Summary of the model:" + @info "State variables: " * $(join(map(string, collect(x_vars)), ", ")) + @info "Parameters: " * $(join(map(string, collect(params)), ", ")) + @info "Inputs: " * $(join(map(string, collect(u_vars)), ", ")) + @info "Outputs: " * $(join(map(string, collect(y_vars)), ", ")) + end + )] # creating the ode object ode_expr = :(StructuralIdentifiability.ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}( $vx, diff --git a/src/logging.jl b/src/logging.jl index 7f993ba1e..f00a2c763 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -26,7 +26,7 @@ const _runtime_logger = Dict( const _si_logger = @static if VERSION >= v"1.7.0" Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info, show_limited = false)) else - Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(Logging.Info)) + Ref{Logging.ConsoleLogger}(Logging.ConsoleLogger(stderr, Logging.Info)) end const _groebner_loglevel = Ref{Int}(0) From 18af6da62e2ea097463179f20731b81032c9694d Mon Sep 17 00:00:00 2001 From: Alexander Demin Date: Mon, 13 Nov 2023 01:09:17 +0300 Subject: [PATCH 19/24] add tests for logging; try to fix @ODEmodel logging --- src/ODE.jl | 22 +++++++++++++--------- test/logging.jl | 19 +++++++++++++++++++ 2 files changed, 32 insertions(+), 9 deletions(-) create mode 100644 test/logging.jl diff --git a/src/ODE.jl b/src/ODE.jl index c7d15afd7..d677479b0 100644 --- a/src/ODE.jl +++ b/src/ODE.jl @@ -459,15 +459,19 @@ macro ODEmodel(ex::Expr...) end end - logging_exprs = [:( - with_logger(StructuralIdentifiability._si_logger[]) do - @info "Summary of the model:" - @info "State variables: " * $(join(map(string, collect(x_vars)), ", ")) - @info "Parameters: " * $(join(map(string, collect(params)), ", ")) - @info "Inputs: " * $(join(map(string, collect(u_vars)), ", ")) - @info "Outputs: " * $(join(map(string, collect(y_vars)), ", ")) - end - )] + logging_exprs = [ + :( + StructuralIdentifiability.Logging.with_logger( + StructuralIdentifiability._si_logger[], + ) do + @info "Summary of the model:" + @info "State variables: " * $(join(map(string, collect(x_vars)), ", ")) + @info "Parameters: " * $(join(map(string, collect(params)), ", ")) + @info "Inputs: " * $(join(map(string, collect(u_vars)), ", ")) + @info "Outputs: " * $(join(map(string, collect(y_vars)), ", ")) + end + ), + ] # creating the ode object ode_expr = :(StructuralIdentifiability.ODE{StructuralIdentifiability.Nemo.fmpq_mpoly}( $vx, diff --git a/test/logging.jl b/test/logging.jl new file mode 100644 index 000000000..b4c5e9fb0 --- /dev/null +++ b/test/logging.jl @@ -0,0 +1,19 @@ +# An attempt to test logging +using Logging + +@testset "Logging" begin + ode = StructuralIdentifiability.@ODEmodel(x'(t) = x^42, y(t) = x) + + # Some logs + @test_logs StructuralIdentifiability.assess_identifiability(ode) + # No logs + @test_logs min_level = Logging.Warn StructuralIdentifiability.assess_identifiability( + ode, + loglevel = Logging.Warn, + ) + # Many logs + @test_logs StructuralIdentifiability.assess_identifiability( + ode, + loglevel = Logging.Debug, + ) +end From e8375acb793f395afbfe4d103c73816a6c07d523 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Mon, 13 Nov 2023 14:31:01 +0100 Subject: [PATCH 20/24] fixed for 1.6 (but nonnegative Groebner) --- src/StructuralIdentifiability.jl | 8 ++++++++ src/logging.jl | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index f2eb868a6..382b4fbc9 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -73,6 +73,14 @@ include("pb_representation.jl") include("submodels.jl") include("discrete.jl") +function __init__() + _si_logger[] = @static if VERSION >= v"1.7.0" + Logging.ConsoleLogger(Logging.Info, show_limited = false) + else + Logging.ConsoleLogger(stderr, Logging.Info) + end +end + """ assess_identifiability(ode; funcs_to_check = [], p=0.99) diff --git a/src/logging.jl b/src/logging.jl index f00a2c763..915f3c7af 100644 --- a/src/logging.jl +++ b/src/logging.jl @@ -42,7 +42,7 @@ function restart_logging(; loglevel = Logging.Info) _runtime_logger[r] = 0 end if loglevel < Logging.Info - _groebner_loglevel[] = -2 + _groebner_loglevel[] = 0 elseif loglevel < Logging.Warn _groebner_loglevel[] = 0 else From 075ffcf9ae00a3dbfd217e869da0556b8399b6af Mon Sep 17 00:00:00 2001 From: gpogudin Date: Mon, 13 Nov 2023 15:37:08 +0100 Subject: [PATCH 21/24] updating docs --- docs/src/identifiability/identifiability.md | 6 ------ docs/src/tutorials/discrete_time.md | 5 +++++ docs/src/tutorials/identifiability.md | 9 ++++++++- docs/src/tutorials/identifiable_functions.md | 3 +++ src/StructuralIdentifiability.jl | 8 +++++--- src/identifiable_functions.jl | 1 + src/local_identifiability.jl | 5 +++-- 7 files changed, 25 insertions(+), 12 deletions(-) diff --git a/docs/src/identifiability/identifiability.md b/docs/src/identifiability/identifiability.md index 3aa66d8a9..cb172fd15 100644 --- a/docs/src/identifiability/identifiability.md +++ b/docs/src/identifiability/identifiability.md @@ -12,12 +12,6 @@ assess_identifiability assess_local_identifiability ``` -## Assessing Global Identifiability - -```@docs -assess_global_identifiability -``` - ## Finding Identifiable Functions ```@docs diff --git a/docs/src/tutorials/discrete_time.md b/docs/src/tutorials/discrete_time.md index 4c160a079..22485d502 100644 --- a/docs/src/tutorials/discrete_time.md +++ b/docs/src/tutorials/discrete_time.md @@ -71,6 +71,11 @@ assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [ principle it may produce incorrect result but the probability of correctness of the returned result is guaranteed to be at least `p` (in fact, the employed bounds are quite conservative, so in practice incorrect result is almost never produced). +As other main functions in the package, `assess_local_identifiability` accepts an optional parameter `loglevel` (default: `Logging.Info`) +to adjust the verbosity of logging. + + + The implementation is based on a version of the observability rank criterion and will be described in a forthcoming paper. [^1]: > S. Nõmm, C. Moog, [*Identifiability of discrete-time nonlinear systems*](https://doi.org/10.1016/S1474-6670(17)31245-4), IFAC Proceedings Volumes, 2004. diff --git a/docs/src/tutorials/identifiability.md b/docs/src/tutorials/identifiability.md index 242b377c6..1beb767b4 100644 --- a/docs/src/tutorials/identifiability.md +++ b/docs/src/tutorials/identifiability.md @@ -62,9 +62,13 @@ Function `assess_local_identifiability` has several optional parameters - `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 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 + + - `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). + - `loglevel` (default `Logging.Info`). The minimal level of logging messages to be displayed. Available options: `Logging.Debug`, + `Logging.Info`, `Logging.Warn`, and `Logging.Error`. + Note that the scaling symmetry given above suggests that $b x_2(t)$ may in fact be identifiable. This can be checked using `funcs_to_check` parameter: ```@example local @@ -108,6 +112,9 @@ Similarly to `assess_local_identifiability`, this function has optional paramete 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. + - `loglevel` (default `Logging.Info`). The minimal level of logging messages to be displayed. Available options: `Logging.Debug`, + `Logging.Info`, `Logging.Warn`, and `Logging.Error`. + Using `funcs_to_check` parameter, one can further inverstigate the nature of the lack of identifiability in the model at hand. For example, for the Goodwin oscillator, we can check if `beta + delta` and `beta * delta` are identifiable: diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md index d9f6ef8d2..c7e830dca 100644 --- a/docs/src/tutorials/identifiable_functions.md +++ b/docs/src/tutorials/identifiable_functions.md @@ -40,4 +40,7 @@ By default, `find_identifiable_functions` tries to simplify the output functions the degree of simplification. The default value is `:standard` but one could use `:strong` to try to simplify further (at the expense of heavier computation) or use `:weak` to simplify less (but compute faster). +As `assess_identifiability` and `assess_local_identifiability`, `find_identifiable_functions` accepts an optional parameter `loglevel` (default: `Logging.Info`) +to adjust the verbosity of logging. + [^1]: > Y. Lecourtier, F. Lamnabhi-Lagarrigue, and E. Walter, [*A method to prove that nonlinear models can be unidentifiable*](https://doi.org/10.1109/CDC.1987.272467), Proceedings of the 26th Conference on Decision and Control, December 1987, 2144-2145; diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl index 382b4fbc9..227d12164 100644 --- a/src/StructuralIdentifiability.jl +++ b/src/StructuralIdentifiability.jl @@ -82,14 +82,15 @@ function __init__() end """ - assess_identifiability(ode; funcs_to_check = [], p=0.99) + assess_identifiability(ode; funcs_to_check = [], p=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. - +- `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`. The function returns a dictionary from the functions to check to their identifiability properties @@ -166,13 +167,14 @@ function _assess_identifiability( end """ - assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], p = 0.99) + assess_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=[], p = 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. +- `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`. diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl index 52a6782c7..18cde16c6 100644 --- a/src/identifiable_functions.jl +++ b/src/identifiable_functions.jl @@ -21,6 +21,7 @@ This functions takes the following optional arguments: - `p`: 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) ## Example diff --git a/src/local_identifiability.jl b/src/local_identifiability.jl index 83a1396fc..df3ef4e04 100644 --- a/src/local_identifiability.jl +++ b/src/local_identifiability.jl @@ -143,7 +143,7 @@ end # ------------------------------------------------------------------------------ """ - function assess_local_identifiability(ode::ModelingToolkit.ODESystem; measured_quantities=Array{ModelingToolkit.Equation}[], funcs_to_check=Array{}[], p::Float64=0.99, type=:SE) + 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) Input: - `ode` - the ODESystem object from ModelingToolkit @@ -151,6 +151,7 @@ Input: - `funcs_to_check` - functions of parameters for which to check identifiability - `p` - 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) Output: - for `type=:SE`, the result is a dictionary from each parameter to boolean; @@ -239,7 +240,7 @@ end # ------------------------------------------------------------------------------ """ - assess_local_identifiability(ode::ODE{P}; funcs_to_check::Array{<: Any, 1}, p::Float64=0.99, type=:SE) where P <: MPolyElem{Nemo.fmpq} + 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} Checks the local identifiability/observability of the functions in `funcs_to_check`. The result is correct with probability at least `p`. From d426726421f3982e49dd6b78fe5eacb88688ae50 Mon Sep 17 00:00:00 2001 From: gpogudin Date: Mon, 13 Nov 2023 16:41:41 +0100 Subject: [PATCH 22/24] Cleaning up the examples: stage 2 --- docs/src/tutorials/discrete_time.md | 2 -- docs/src/tutorials/identifiability.md | 11 +++----- examples/CRN.jl | 7 +----- examples/Crauste.jl | 11 ++++---- examples/Fujita.jl | 11 ++++---- examples/Goodwin.jl | 11 ++++---- examples/HIV.jl | 11 ++++---- examples/LV.jl | 11 ++++---- examples/MAPK.jl | 7 +----- examples/NFkB.jl | 5 ---- examples/Pharm.jl | 7 +----- examples/QWWC.jl | 6 ----- examples/SEAIJRC.jl | 11 ++++---- examples/SIAR.jl | 10 +++----- examples/SIRS_forced.jl | 9 +++---- examples/SIWR.jl | 7 +----- examples/SIWR_multioutput.jl | 7 +----- examples/arabidopsis.jl | 36 --------------------------- examples/multiout_counterexample.jl | 7 +----- 19 files changed, 48 insertions(+), 139 deletions(-) delete mode 100644 examples/arabidopsis.jl diff --git a/docs/src/tutorials/discrete_time.md b/docs/src/tutorials/discrete_time.md index 22485d502..4890643b5 100644 --- a/docs/src/tutorials/discrete_time.md +++ b/docs/src/tutorials/discrete_time.md @@ -74,8 +74,6 @@ assess_local_identifiability(sir; measured_quantities = [I], funcs_to_check = [ As other main functions in the package, `assess_local_identifiability` accepts an optional parameter `loglevel` (default: `Logging.Info`) to adjust the verbosity of logging. - - The implementation is based on a version of the observability rank criterion and will be described in a forthcoming paper. [^1]: > S. Nõmm, C. Moog, [*Identifiability of discrete-time nonlinear systems*](https://doi.org/10.1016/S1474-6670(17)31245-4), IFAC Proceedings Volumes, 2004. diff --git a/docs/src/tutorials/identifiability.md b/docs/src/tutorials/identifiability.md index 1beb767b4..d89ceb5f4 100644 --- a/docs/src/tutorials/identifiability.md +++ b/docs/src/tutorials/identifiability.md @@ -62,12 +62,10 @@ Function `assess_local_identifiability` has several optional parameters - `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 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 + - `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). - - - `loglevel` (default `Logging.Info`). The minimal level of logging messages to be displayed. Available options: `Logging.Debug`, - `Logging.Info`, `Logging.Warn`, and `Logging.Error`. + - `loglevel` (default `Logging.Info`). The minimal level of logging messages to be displayed. Available options: `Logging.Debug`, + `Logging.Info`, `Logging.Warn`, and `Logging.Error`. Note that the scaling symmetry given above suggests that $b x_2(t)$ may in fact be identifiable. This can be checked using `funcs_to_check` parameter: @@ -111,8 +109,7 @@ Similarly to `assess_local_identifiability`, this function has optional paramete 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. - - - `loglevel` (default `Logging.Info`). The minimal level of logging messages to be displayed. Available options: `Logging.Debug`, + - `loglevel` (default `Logging.Info`). The minimal level of logging messages to be displayed. Available options: `Logging.Debug`, `Logging.Info`, `Logging.Warn`, and `Logging.Error`. Using `funcs_to_check` parameter, one can further inverstigate the nature of the lack of identifiability in the model at hand. diff --git a/examples/CRN.jl b/examples/CRN.jl index 568b14df6..e6db44e48 100644 --- a/examples/CRN.jl +++ b/examples/CRN.jl @@ -1,13 +1,8 @@ # Conradi, C., Shiu, A., # Dynamics of post-translational modification systems: recent progress and future directions # Eq. 3.4 -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t), x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t), @@ -19,4 +14,4 @@ ode = @ODEmodel( y2(t) = x2(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) diff --git a/examples/Crauste.jl b/examples/Crauste.jl index 56c5051cb..f592e6bf3 100644 --- a/examples/Crauste.jl +++ b/examples/Crauste.jl @@ -1,10 +1,5 @@ -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - ode = @ODEmodel( N'(t) = -N(t) * mu_N - N(t) * P(t) * delta_NE, E'(t) = @@ -17,4 +12,8 @@ ode = @ODEmodel( y3(t) = M(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) + +# Not everything is identifiable, so we may wonder what are the identifiable functions + +println(find_identifiable_functions(ode, with_states = true)) diff --git a/examples/Fujita.jl b/examples/Fujita.jl index 9901b1b80..df0588ebb 100644 --- a/examples/Fujita.jl +++ b/examples/Fujita.jl @@ -2,13 +2,8 @@ # Decoupling of receptor and downstream signals in the Akt pathway by its low-pass filter characteristics # https://doi.org/10.1126/scisignal.2000810 -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - ode = @ODEmodel( EGFR'(t) = EGFR_turnover * pro_EGFR(t) + EGF_EGFR(t) * reaction_1_k2 - @@ -42,4 +37,8 @@ ode = @ODEmodel( y3(t) = a3 * pS6(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) + +# Not everything is identifiable, so we may wonder what are the identifiable functions + +println(find_identifiable_functions(ode, with_states = true)) diff --git a/examples/Goodwin.jl b/examples/Goodwin.jl index a6716bc0d..edf767137 100644 --- a/examples/Goodwin.jl +++ b/examples/Goodwin.jl @@ -2,13 +2,8 @@ # Oscillatory behavior in enzymatic control processes # https://doi.org/10.1016/0065-2571(65)90067-1 -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( x1'(t) = -b * x1(t) + 1 / (c + x4(t)), x2'(t) = alpha * x1(t) - beta * x2(t), @@ -17,4 +12,8 @@ ode = @ODEmodel( y(t) = x1(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) + +# Not everything is identifiable, so we may wonder what are the identifiable functions + +println(find_identifiable_functions(ode, with_states = true)) diff --git a/examples/HIV.jl b/examples/HIV.jl index 5386098af..a94fa23da 100644 --- a/examples/HIV.jl +++ b/examples/HIV.jl @@ -2,13 +2,8 @@ # Specific therapy regimes could lead to long-term immunological control of HIV # https://doi.org/10.1073/pnas.96.25.14464 # Page 1 -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( x'(t) = lm - d * x(t) - beta * x(t) * v(t), y'(t) = beta * x(t) * v(t) - a * y(t), @@ -19,4 +14,8 @@ ode = @ODEmodel( y2(t) = z(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) + +# Not everything is identifiable, so we may wonder what are the identifiable functions + +println(find_identifiable_functions(ode, with_states = true)) diff --git a/examples/LV.jl b/examples/LV.jl index 192756114..8767202b4 100644 --- a/examples/LV.jl +++ b/examples/LV.jl @@ -1,15 +1,14 @@ # Lotka-Volterra model -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - ode = @ODEmodel( x1'(t) = a * x1(t) - b * x1(t) * x2(t) + u(t), x2'(t) = -c * x2(t) + d * x1(t) * x2(t), y(t) = x1(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) + +# Not everything is identifiable, so we may wonder what are the identifiable functions + +println(find_identifiable_functions(ode, with_states = true)) diff --git a/examples/MAPK.jl b/examples/MAPK.jl index 8eeb77ee4..c533d79ed 100644 --- a/examples/MAPK.jl +++ b/examples/MAPK.jl @@ -1,13 +1,8 @@ # A. K. Manrai and J. Gunawardena. # The geometry of multisite phosphorylation # https://dx.doi.org/10.1529%2Fbiophysj.108.140632 -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( KS00'(t) = -a00 * K(t) * S00(t) + @@ -52,4 +47,4 @@ ode = @ODEmodel( y5(t) = S11(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) diff --git a/examples/NFkB.jl b/examples/NFkB.jl index 548883659..263ee7236 100644 --- a/examples/NFkB.jl +++ b/examples/NFkB.jl @@ -2,13 +2,8 @@ # An iterative identification procedure for dynamic modeling of biochemical networks # https://doi.org/10.1186/1752-0509-4-11 # The model is out of reach at the moment but can be analyzed by SIAN (https://github.com/pogudingleb/SIAN) -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( x1'(t) = k_prod - k_deg * x1(t) - k1 * x1(t) * u(t), x2'(t) = diff --git a/examples/Pharm.jl b/examples/Pharm.jl index ad621fe15..794232383 100644 --- a/examples/Pharm.jl +++ b/examples/Pharm.jl @@ -1,13 +1,8 @@ # Demignot, S., D., D., # Effect of prosthetic sugar groups on the pharmacokinetics of glucose-oxidase # https://pubmed.ncbi.nlm.nih.gov/2855567/ -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - ode = @ODEmodel( x0'(t) = a1 * (x1(t) - x0(t)) - (ka * n * x0(t)) / (kc * ka + kc * x2(t) + ka * x0(t)), @@ -18,4 +13,4 @@ ode = @ODEmodel( y1(t) = x0(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) diff --git a/examples/QWWC.jl b/examples/QWWC.jl index f298842b4..88aa1cebb 100644 --- a/examples/QWWC.jl +++ b/examples/QWWC.jl @@ -1,11 +1,5 @@ -using Logging - using StructuralIdentifiability -#QWWC -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - ode = @ODEmodel( x'(t) = a * (y(t) - x(t)) + y(t) * z(t), y'(t) = b * (x(t) + y(t)) - x(t) * z(t), diff --git a/examples/SEAIJRC.jl b/examples/SEAIJRC.jl index 35032df85..42a987f8e 100644 --- a/examples/SEAIJRC.jl +++ b/examples/SEAIJRC.jl @@ -2,13 +2,8 @@ # Assessing parameter identifiability in compartmental dynamic models using a computational approach: application to infectious disease transmission models # https://doi.org/10.1186/s12976-018-0097-6 # Model 2 -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - ode = @ODEmodel( S'(t) = -b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv, E'(t) = b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv - k * E(t), @@ -20,4 +15,8 @@ ode = @ODEmodel( y2(t) = Ninv ) -@time println(assess_global_identifiability(ode)) +println(assess_identifiability(ode)) + +# Not everything is identifiable, so we may wonder what are the identifiable functions + +println(find_identifiable_functions(ode, with_states = true)) diff --git a/examples/SIAR.jl b/examples/SIAR.jl index c8300a470..1b9493dd4 100644 --- a/examples/SIAR.jl +++ b/examples/SIAR.jl @@ -1,12 +1,10 @@ # SIAR model -# https://arxiv.org/pdf/2012.00443.pdf -using Logging - +# L. Gallo, M.Frasca, V. Latora, G. Russo +# Lack of practical identifiability may hamper reliable predictions in COVID-19 epidemic models +# Science Advances, DOI: 10.1126/sciadv.abg5234 +# Equation (21) using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( S'(t) = -S(t) * (betaIa * Ia(t) + betaIs * Is(t) + betaH * H(t) + betaT * T(t)) * Ninv, diff --git a/examples/SIRS_forced.jl b/examples/SIRS_forced.jl index 6cddef64f..a9814809e 100644 --- a/examples/SIRS_forced.jl +++ b/examples/SIRS_forced.jl @@ -2,13 +2,8 @@ # Parameter Estimation of Some Epidemic Models. The Case of Recurrent Epidemics Caused by Respiratory Syncytial Virus # doi.org/10.1007/s11538-009-9429-3 # Equations (7)-(11) -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( s'(t) = mu - mu * s(t) - b0 * (1 + b1 * x1(t)) * i(t) * s(t) + g * r(t), i'(t) = b0 * (1 + b1 * x1(t)) * i(t) * s(t) - (nu + mu) * i(t), @@ -20,3 +15,7 @@ ode = @ODEmodel( ) @time println(assess_identifiability(ode)) + +# Not everything is identifiable, so we may wonder what are the identifiable functions + +println(find_identifiable_functions(ode, with_states = true)) diff --git a/examples/SIWR.jl b/examples/SIWR.jl index 8b141176b..7e4c25da4 100644 --- a/examples/SIWR.jl +++ b/examples/SIWR.jl @@ -2,13 +2,8 @@ # Model distinguishability and inference robustness in mechanisms of cholera transmission and loss of immunity # https://doi.org/10.1016/j.jtbi.2017.01.032 # Eq. (3) -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t), I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t), @@ -17,4 +12,4 @@ ode = @ODEmodel( y(t) = k * I(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) diff --git a/examples/SIWR_multioutput.jl b/examples/SIWR_multioutput.jl index 1cc1b2d80..fbce7cc46 100644 --- a/examples/SIWR_multioutput.jl +++ b/examples/SIWR_multioutput.jl @@ -2,13 +2,8 @@ # Model distinguishability and inference robustness in mechanisms of cholera transmission and loss of immunity # https://doi.org/10.1016/j.jtbi.2017.01.032 # Eq. (3) + extra output -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - ode = @ODEmodel( x_0'(t) = mu - bi * x_0(t) * x_1(t) - bw * x_0(t) * x_2(t) - mu * x_0(t) + a * x_3(t), @@ -19,4 +14,4 @@ ode = @ODEmodel( y2(t) = x_0(t) + x_1(t) + x_3(t) ) -@time println(assess_global_identifiability(ode)) +println(assess_identifiability(ode)) diff --git a/examples/arabidopsis.jl b/examples/arabidopsis.jl deleted file mode 100644 index b627f7568..000000000 --- a/examples/arabidopsis.jl +++ /dev/null @@ -1,36 +0,0 @@ -include("io_equation.jl") - -logger = Logging.SimpleLogger(stdout, Logging.Debug) -global_logger(logger) - -varnames = vcat(["x$i" for i in 1:7], ["p$i" for i in 1:27], ["u"]) - -R, Rvars = PolynomialRing(QQ, varnames) -xs = Rvars[1:7] -ps = Rvars[8:(end - 1)] -u = Rvars[end] - -ode = ODE( - Dict( - xs[1] => - ps[1] * xs[6] // (ps[3] + xs[6]) - ps[5] * xs[1] // (ps[12] + xs[1]) + - ps[26] * xs[7] * u, - xs[2] => - ps[19] * xs[1] - ps[22] * xs[2] + ps[23] * xs[3] - - ps[6] * xs[2] // (ps[13] + xs[2]), - xs[3] => ps[22] * xs[2] - ps[23] * xs[3] - ps[7] * xs[3] // (ps[14] + xs[3]), - xs[4] => - ps[2] * ps[4]^2 // (ps[4]^2 + xs[3]^2) - ps[8] * xs[4] // (ps[15] + xs[4]), - xs[5] => - ps[20] * xs[4] - ps[24] * xs[5] + ps[25] * xs[6] - - ps[9] * xs[5] // (ps[16] + xs[5]), - xs[6] => ps[24] * xs[5] - ps[25] * xs[6] - ps[10] * xs[6] // (ps[17] + xs[6]), - xs[7] => - ps[21] - ps[11] * xs[7] // (ps[18] + xs[7]) - (ps[21] + ps[27] * xs[7]) * u, - ), - [u], -) - -@time io_equations = find_ioequations(ode, [xs[1], xs[4]]) - -print(map(length, io_equations), "\n") diff --git a/examples/multiout_counterexample.jl b/examples/multiout_counterexample.jl index a70a0ed72..5bad882c5 100644 --- a/examples/multiout_counterexample.jl +++ b/examples/multiout_counterexample.jl @@ -1,11 +1,6 @@ # an artificial example illustrating the necessity of an extra projection in the multi-output case -using Logging - using StructuralIdentifiability -logger = Logging.SimpleLogger(stdout, Logging.Info) -global_logger(logger) - ode = @ODEmodel( x1'(t) = (1 + x1(t)^2) // 2, x2'(t) = (1 - x1(t)^2) // (1 + x1(t)^2), @@ -13,4 +8,4 @@ ode = @ODEmodel( y2(t) = x2(t) ) -@time println(assess_global_identifiability(ode)) +@time println(assess_identifiability(ode)) From 6bd56f262d2740e312c416f6ae853ee0e20166e6 Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Mon, 13 Nov 2023 21:00:21 +0100 Subject: [PATCH 23/24] fixing a bug with states + random projection --- examples/multiout_counterexample.jl | 2 +- src/states.jl | 8 ++++++-- test/identifiability.jl | 4 ++-- 3 files changed, 9 insertions(+), 5 deletions(-) diff --git a/examples/multiout_counterexample.jl b/examples/multiout_counterexample.jl index 5bad882c5..a57ac9cfc 100644 --- a/examples/multiout_counterexample.jl +++ b/examples/multiout_counterexample.jl @@ -8,4 +8,4 @@ ode = @ODEmodel( y2(t) = x2(t) ) -@time println(assess_identifiability(ode)) +println(assess_identifiability(ode)) diff --git a/src/states.jl b/src/states.jl index 0936770ec..c706a9ab5 100644 --- a/src/states.jl +++ b/src/states.jl @@ -83,8 +83,12 @@ identifiable functions of parameters only y_to_ord = Dict{P, Int}() ynames = [var_to_str(y) for y in ode.y_vars] for (leader, ioeq) in io_equations - (y_str, ord) = decompose_derivative(var_to_str(leader), ynames) - y_to_ord[str_to_var(y_str, parent(ode))] = ord + decomposition = decompose_derivative(var_to_str(leader), ynames) + # decomposition will be `nothing` for extra random projection + if !isnothing(decomposition) + (y_str, ord) = decomposition + y_to_ord[str_to_var(y_str, parent(ode))] = ord + end end result = Array{Generic.Frac{P}, 1}() diff --git a/test/identifiability.jl b/test/identifiability.jl index 5b3380d16..361cb5cd5 100644 --- a/test/identifiability.jl +++ b/test/identifiability.jl @@ -130,8 +130,8 @@ y1(t) = 2 * x1(t) // (b * (1 + x1(t)^2)), y2(t) = x2(t) ) - funcs_to_test = [b] - correct = [:globally] + funcs_to_test = [b, x1, x2] + correct = [:globally, :globally, :globally] push!( test_cases, Dict( From 6c7171263a0fe341d02d28c8efa7538464fc085c Mon Sep 17 00:00:00 2001 From: pogudingleb Date: Mon, 13 Nov 2023 22:58:51 +0100 Subject: [PATCH 24/24] updating the examples --- examples/SIAR.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/SIAR.jl b/examples/SIAR.jl index 1b9493dd4..1e5f93b67 100644 --- a/examples/SIAR.jl +++ b/examples/SIAR.jl @@ -26,4 +26,4 @@ ode = @ODEmodel( y5(t) = Ninv ) -@time println(assess_identifiability(ode)) +println(find_identifiable_functions(ode, funcs_to_check = ode.parameters))