Skip to content

Commit

Permalink
wrapping to with_logger
Browse files Browse the repository at this point in the history
  • Loading branch information
@pogudingleb
pogudingleb committed Nov 12, 2023
1 parent bf7a857 commit 51cf590
Show file tree
Hide file tree
Showing 19 changed files with 365 additions and 331 deletions.
22 changes: 8 additions & 14 deletions src/ODE.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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}(
Expand Down
32 changes: 15 additions & 17 deletions src/RationalFunctionFields/IdealMQS.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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"
Expand All @@ -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]
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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)
Expand Down
82 changes: 35 additions & 47 deletions src/RationalFunctionFields/RationalFunctionField.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)))
Expand All @@ -121,30 +121,30 @@ 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 *
BigInt(degree)^(length(total_vars) + 3) *
(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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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))
Expand All @@ -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,
Expand All @@ -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: <simplified generators> in <original generators>
new_rff = RationalFunctionField(fracs)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down
Loading

0 comments on commit 51cf590

Please sign in to comment.