Skip to content

Commit

Permalink
updating to the newest AA and Nemo
Browse files Browse the repository at this point in the history
  • Loading branch information
pogudingleb committed Feb 23, 2024
1 parent 9453d5f commit 79e433b
Show file tree
Hide file tree
Showing 26 changed files with 91 additions and 82 deletions.
2 changes: 1 addition & 1 deletion ext/ModelingToolkitSIExt.jl
Original file line number Diff line number Diff line change
Expand Up @@ -165,7 +165,7 @@ function __mtk_to_si(
measured_quantities::Array{<:Tuple{String, <:SymbolicUtils.BasicSymbolic}},
)
polytype = StructuralIdentifiability.Nemo.QQMPolyRingElem
fractype = StructuralIdentifiability.Nemo.Generic.Frac{polytype}
fractype = StructuralIdentifiability.Nemo.Generic.FracFieldElem{polytype}
diff_eqs =
filter(eq -> !(ModelingToolkit.isoutput(eq.lhs)), ModelingToolkit.equations(de))
# performing full structural simplification
Expand Down
28 changes: 14 additions & 14 deletions src/ODE.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,14 +11,14 @@ struct ODE{P} # P is the type of polynomials in the rhs of the ODE system
y_vars::Array{P, 1}
u_vars::Array{P, 1}
parameters::Array{P, 1}
x_equations::Dict{P, <:Union{P, Generic.Frac{P}}}
y_equations::Dict{P, <:Union{P, Generic.Frac{P}}}
x_equations::Dict{P, <:ExtendedFraction{P}}
y_equations::Dict{P, <:ExtendedFraction{P}}

function ODE{P}(
x_vars::Array{P, 1},
y_vars::Array{P, 1},
x_eqs::Dict{P, <:Union{P, Generic.Frac{P}}},
y_eqs::Dict{P, <:Union{P, Generic.Frac{P}}},
x_eqs::Dict{P, <:ExtendedFraction{P}},
y_eqs::Dict{P, <:ExtendedFraction{P}},
inputs::Array{P, 1},
) where {P <: MPolyRingElem{<:FieldElem}}
# Initialize ODE
Expand All @@ -37,8 +37,8 @@ struct ODE{P} # P is the type of polynomials in the rhs of the ODE system
end

function ODE{P}(
x_eqs::Dict{P, <:Union{P, Generic.Frac{P}}},
y_eqs::Dict{P, <:Union{P, Generic.Frac{P}}},
x_eqs::Dict{P, <:ExtendedFraction{P}},
y_eqs::Dict{P, <:ExtendedFraction{P}},
inputs::Array{P, 1},
) where {P <: MPolyRingElem{<:FieldElem}}
x_vars = collect(keys(x_eqs))
Expand All @@ -62,19 +62,19 @@ function add_outputs(
new_ring, new_vars = Nemo.polynomial_ring(base_ring(ode.poly_ring), new_var_names)

new_x = Array{P, 1}([parent_ring_change(x, new_ring) for x in ode.x_vars])
new_x_eqs = Dict{P, Union{P, Generic.Frac{P}}}(
new_x_eqs = Dict{P, ExtendedFraction{P}}(
parent_ring_change(x, new_ring) => parent_ring_change(f, new_ring) for
(x, f) in ode.x_equations
)
new_y = Array{P, 1}([parent_ring_change(y, new_ring) for y in ode.y_vars])
for y in keys(extra_y)
push!(new_y, str_to_var(y, new_ring))
end
new_y_eqs = Dict{P, Union{P, Generic.Frac{P}}}(
new_y_eqs = Dict{P, ExtendedFraction{P}}(
parent_ring_change(y, new_ring) => parent_ring_change(g, new_ring) for
(y, g) in ode.y_equations
)
extra_y_eqs = Dict{P, Union{P, Generic.Frac{P}}}(
extra_y_eqs = Dict{P, ExtendedFraction{P}}(
str_to_var(y, new_ring) => parent_ring_change(g, new_ring) for (y, g) in extra_y
)
merge!(new_y_eqs, extra_y_eqs)
Expand Down Expand Up @@ -106,11 +106,11 @@ function set_parameter_values(
merge!(eval_dict, Dict(p => small_ring(val) for (p, val) in param_values))

return ODE{P}(
Dict{P, Union{P, Generic.Frac{P}}}(
Dict{P, ExtendedFraction{P}}(
eval_at_dict(v, eval_dict) => eval_at_dict(f, eval_dict) for
(v, f) in ode.x_equations
),
Dict{P, Union{P, Generic.Frac{P}}}(
Dict{P, ExtendedFraction{P}}(
eval_at_dict(v, eval_dict) => eval_at_dict(f, eval_dict) for
(v, f) in ode.y_equations
),
Expand Down Expand Up @@ -225,7 +225,7 @@ function _reduce_mod_p(poly::QQMPolyRingElem, p::Int)
return change_base_ring(Nemo.Native.GF(p), num) * (1 // Nemo.Native.GF(p)(den))
end

function _reduce_mod_p(rat::Generic.Frac{QQMPolyRingElem}, p::Int)
function _reduce_mod_p(rat::Generic.FracFieldElem{QQMPolyRingElem}, p::Int)
num, den = map(poly -> _reduce_mod_p(poly, p), [numerator(rat), denominator(rat)])
if den == 0
throw(Base.ArgumentError("Prime $p divides the denominator of $rat"))
Expand All @@ -248,8 +248,8 @@ function reduce_ode_mod_p(ode::ODE{<:MPolyRingElem{Nemo.QQFieldElem}}, p::Int)
new_inputs = map(u -> switch_ring(u, new_ring), ode.u_vars)
new_x = map(x -> switch_ring(x, new_ring), ode.x_vars)
new_y = map(y -> switch_ring(y, new_ring), ode.y_vars)
new_x_eqs = Dict{new_type, Union{new_type, Generic.Frac{new_type}}}()
new_y_eqs = Dict{new_type, Union{new_type, Generic.Frac{new_type}}}()
new_x_eqs = Dict{new_type, ExtendedFraction{new_type}}()
new_y_eqs = Dict{new_type, ExtendedFraction{new_type}}()
for (old, new) in Dict(ode.x_equations => new_x_eqs, ode.y_equations => new_y_eqs)
for (v, f) in old
new_v = switch_ring(v, new_ring)
Expand Down
6 changes: 3 additions & 3 deletions src/RationalFunctionFields/IdealMQS.jl
Original file line number Diff line number Diff line change
Expand Up @@ -100,7 +100,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal
end
varnames = append_at_index(ystrs, sat_var_index, sat_varname)
@debug "Saturating variable is $sat_varname, index is $sat_var_index"
R_sat, v_sat = Nemo.polynomial_ring(K, varnames, ordering = ordering)
R_sat, v_sat = Nemo.polynomial_ring(K, varnames, internal_ordering = ordering)
# Saturation
t_sat = v_sat[sat_var_index]
den_lcm_orig = den_lcm
Expand Down Expand Up @@ -144,7 +144,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal
push!(nums_qq, num)
end
end
parent_ring_param, _ = polynomial_ring(ring, varnames, ordering = ordering)
parent_ring_param, _ = polynomial_ring(ring, varnames, internal_ordering = ordering)
@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 @@ -202,7 +202,7 @@ function fractionfree_generators_raw(mqs::IdealMQS)
end
# NOTE: new variables go first!
big_ring, big_vars =
polynomial_ring(K, vcat(new_varnames, old_varnames), ordering = :lex)
polynomial_ring(K, vcat(new_varnames, old_varnames), internal_ordering = :lex)
@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)
Expand Down
6 changes: 3 additions & 3 deletions src/RationalFunctionFields/RationalFunctionField.jl
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ mutable struct RationalFunctionField{T}
function RationalFunctionField(polys::Vector{T}) where {T}
RationalFunctionField(polys .// one(parent(first(polys))))
end
function RationalFunctionField(fractions::Vector{Generic.Frac{T}}) where {T}
function RationalFunctionField(fractions::Vector{Generic.FracFieldElem{T}}) where {T}
RationalFunctionField(fractions_to_dennums(fractions))
end
function RationalFunctionField(dennums::Vector{Vector{T}}) where {T}
Expand Down Expand Up @@ -156,7 +156,7 @@ end

function field_contains(
field::RationalFunctionField{T},
ratfuncs::Vector{Generic.Frac{T}},
ratfuncs::Vector{Generic.FracFieldElem{T}},
prob_threshold,
) where {T}
return field_contains(field, fractions_to_dennums(ratfuncs), prob_threshold)
Expand Down Expand Up @@ -200,7 +200,7 @@ function check_algebraicity(field, ratfuncs, p)
eval_point = [Nemo.QQ(rand(1:D)) for x in base_vars]

# Filling the jacobain for generators
S = MatrixSpace(Nemo.QQ, length(base_vars), length(fgens) + 1)
S = matrix_space(Nemo.QQ, length(base_vars), length(fgens) + 1)
J = zero(S)
for (i, f) in enumerate(fgens)
for (j, x) in enumerate(base_vars)
Expand Down
4 changes: 2 additions & 2 deletions src/RationalFunctionFields/normalforms.jl
Original file line number Diff line number Diff line change
Expand Up @@ -291,7 +291,7 @@ Note: uses Monte-Carlo probabilistic algorithm. The probability of correctness
is not specified but is assumed to be close to 1.
"""
@timeit _to function linear_relations_between_normal_forms(
fracs::Vector{Generic.Frac{T}},
fracs::Vector{Generic.FracFieldElem{T}},
up_to_degree::Integer;
seed = 42,
) where {T}
Expand Down Expand Up @@ -391,7 +391,7 @@ is not specified but is assumed to be close to 1.
end
union!(complete_intersection_relations_ff, relations_ff_1)
@debug "Reconstructing relations to rationals"
relations_qq = Vector{Generic.Frac{elem_type(ring_param)}}(
relations_qq = Vector{Generic.FracFieldElem{elem_type(ring_param)}}(
undef,
length(complete_intersection_relations_ff),
)
Expand Down
2 changes: 1 addition & 1 deletion src/RationalFunctionFields/util.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ Output: an array of fractions
`[f2/f1, f3/f1, ..., g2/g1, g3/g1, ...]`
"""
function dennums_to_fractions(dennums::Vector{Vector{T}}) where {T}
fractions = Vector{AbstractAlgebra.Generic.Frac{T}}()
fractions = Vector{AbstractAlgebra.Generic.FracFieldElem{T}}()
for dni in dennums
den, nums = dni[1], dni[2:end]
isempty(nums) && continue
Expand Down
4 changes: 3 additions & 1 deletion src/StructuralIdentifiability.jl
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,8 @@ export find_submodels
# finding identifiabile reparametrizations
export reparametrize_global

ExtendedFraction{P} = Union{P, Generic.FracFieldElem{P}}

include("logging.jl")
include("util.jl")
include("power_series_utils.jl")
Expand Down Expand Up @@ -102,7 +104,7 @@ The function returns an (ordered) dictionary from the functions to check to thei
function assess_identifiability(
ode::ODE{P};
funcs_to_check = Vector(),
known_ic::Vector{<:Union{P, Generic.Frac{P}}} = Vector{Union{P, Generic.Frac{P}}}(),
known_ic::Vector{<:ExtendedFraction{P}} = Vector{ExtendedFraction{P}}(),
prob_threshold::Float64 = 0.99,
loglevel = Logging.Info,
) where {P <: MPolyRingElem{QQFieldElem}}
Expand Down
4 changes: 2 additions & 2 deletions src/discrete.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,8 @@ struct DDS{P} # P is the type of polynomials in the rhs of the DDS system
function DDS{P}(
x_vars::Array{P, 1},
y_vars::Array{P, 1},
x_eqs::Dict{P, <:Union{P, Generic.Frac{P}}},
y_eqs::Dict{P, <:Union{P, Generic.Frac{P}}},
x_eqs::Dict{P, <:ExtendedFraction{P}},
y_eqs::Dict{P, <:ExtendedFraction{P}},
u_vars::Array{P, 1},
) where {P <: MPolyRingElem{<:FieldElem}}
new{P}(ODE{P}(x_vars, y_vars, x_eqs, y_eqs, u_vars))
Expand Down
4 changes: 2 additions & 2 deletions src/global_identifiability.jl
Original file line number Diff line number Diff line change
Expand Up @@ -148,7 +148,7 @@ The function returns a tuple containing the following:
param_ring, _ = polynomial_ring(
base_ring(bring),
map(string, ode.parameters),
ordering = Nemo.ordering(bring),
internal_ordering = Nemo.internal_ordering(bring),
)
id_funcs_no_states_param = map(
polys -> map(poly -> parent_ring_change(poly, param_ring), polys),
Expand Down Expand Up @@ -218,7 +218,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i
with_states = states_needed,
)

funcs_to_check = Vector{Generic.Frac{P}}(
funcs_to_check = Vector{Generic.FracFieldElem{P}}(
map(f -> parent_ring_change(f, bring) // one(bring), funcs_to_check),
)

Expand Down
4 changes: 2 additions & 2 deletions src/identifiable_functions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -40,15 +40,15 @@ ode = @ODEmodel(
find_identifiable_functions(ode)
# prints
3-element Vector{AbstractAlgebra.Generic.Frac{Nemo.QQMPolyRingElem}}:
3-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}:
a12 + a01 + a21
a12*a01
```
"""
function find_identifiable_functions(
ode::ODE{T};
known_ic::Vector{<:Union{T, Generic.Frac{T}}} = Vector{Union{T, Generic.Frac{T}}}(),
known_ic::Vector{<:ExtendedFraction{T}} = Vector{ExtendedFraction{T}}(),
prob_threshold::Float64 = 0.99,
seed = 42,
with_states = false,
Expand Down
6 changes: 3 additions & 3 deletions src/input_macro.jl
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ function generate_model_code(type, ex::Expr...)
StructuralIdentifiability.Nemo.QQMPolyRingElem,
Union{
StructuralIdentifiability.Nemo.QQMPolyRingElem,
StructuralIdentifiability.AbstractAlgebra.Generic.Frac{
StructuralIdentifiability.AbstractAlgebra.Generic.FracFieldElem{
StructuralIdentifiability.Nemo.QQMPolyRingElem,
},
},
Expand All @@ -147,7 +147,7 @@ function generate_model_code(type, ex::Expr...)
StructuralIdentifiability.Nemo.QQMPolyRingElem,
Union{
StructuralIdentifiability.Nemo.QQMPolyRingElem,
StructuralIdentifiability.AbstractAlgebra.Generic.Frac{
StructuralIdentifiability.AbstractAlgebra.Generic.FracFieldElem{
StructuralIdentifiability.Nemo.QQMPolyRingElem,
},
},
Expand Down Expand Up @@ -304,4 +304,4 @@ macro DDSmodel(ex::Expr...)
return esc(generate_model_code(:dds, ex...))
end

#------------------------------------------------------------------------------
#-----------------------------------------------------------------------------i
2 changes: 1 addition & 1 deletion src/io_equation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ function diff_poly(poly::P, derivation::Dict{P, T}) where {P <: MPolyRingElem, T
return sum(derivative(poly, x) * xd for (x, xd) in derivation)
end

function diff_frac(frac::Generic.Frac{F}, derivation::Dict{P, T}) where {F, P, T}
function diff_frac(frac::Generic.FracFieldElem{F}, derivation::Dict{P, T}) where {F, P, T}
num, den = unpack_fraction(frac)
numd = den * diff_poly(num, derivation) - num * diff_poly(den, derivation)
dend = den^2
Expand Down
8 changes: 4 additions & 4 deletions src/known_ic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,12 +28,12 @@ This functions takes the following optional arguments:
"""
function _find_identifiable_functions_kic(
ode::ODE{T},
known_ic::Vector{<:Union{T, Generic.Frac{T}}};
known_ic::Vector{<:ExtendedFraction{T}};
prob_threshold::Float64 = 0.99,
seed = 42,
simplify = :standard,
rational_interpolator = :VanDerHoevenLecerf,
) where {T <: MPolyElem{fmpq}}
) where {T <: MPolyRingElem{Nemo.QQFieldElem}}
Random.seed!(seed)
@assert simplify in (:standard, :weak, :strong, :absent)
half_p = 0.5 + prob_threshold / 2
Expand Down Expand Up @@ -80,10 +80,10 @@ The function returns an (ordered) dictionary from the functions to check to thei
"""
function _assess_identifiability_kic(
ode::ODE{P},
known_ic::Vector{<:Union{P, Generic.Frac{P}}};
known_ic::Vector{<:ExtendedFraction{P}};
funcs_to_check = Vector(),
prob_threshold::Float64 = 0.99,
) where {P <: MPolyElem{fmpq}}
) where {P <: MPolyRingElem{Nemo.QQFieldElem}}
runtime_start = time_ns()
if length(funcs_to_check) == 0
funcs_to_check = vcat(ode.x_vars, ode.parameters)
Expand Down
11 changes: 4 additions & 7 deletions src/lincomp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,7 @@ function linear_compartment_model(
)
x_vars = @view vars[1:n]
x_equations =
Dict{QQMPolyRingElem, Union{QQMPolyRingElem, Generic.Frac{QQMPolyRingElem}}}(
x => R(0) for x in x_vars
)
Dict{QQMPolyRingElem, ExtendedFraction{QQMPolyRingElem}}(x => R(0) for x in x_vars)
for i in 1:n
for j in graph[i]
rate = str_to_var("a_$(j)_$(i)", R)
Expand All @@ -69,10 +67,9 @@ function linear_compartment_model(
end
end

y_equations =
Dict{QQMPolyRingElem, Union{QQMPolyRingElem, Generic.Frac{QQMPolyRingElem}}}(
str_to_var("y$i", R) => str_to_var("x$i", R) for i in outputs
)
y_equations = Dict{QQMPolyRingElem, ExtendedFraction{QQMPolyRingElem}}(
str_to_var("y$i", R) => str_to_var("x$i", R) for i in outputs
)

return ODE{QQMPolyRingElem}(
x_equations,
Expand Down
4 changes: 3 additions & 1 deletion src/local_identifiability.jl
Original file line number Diff line number Diff line change
Expand Up @@ -135,7 +135,9 @@ function get_degree_and_coeffsize(f::MPolyRingElem{Nemo.QQFieldElem})
return (total_degree(f), max_coef)
end

function get_degree_and_coeffsize(f::Generic.Frac{<:MPolyRingElem{Nemo.QQFieldElem}})
function get_degree_and_coeffsize(
f::Generic.FracFieldElem{<:MPolyRingElem{Nemo.QQFieldElem}},
)
num_deg, num_coef = get_degree_and_coeffsize(numerator(f))
den_deg, den_coef = get_degree_and_coeffsize(denominator(f))
return (max(num_deg, den_deg), max(num_coef, den_coef))
Expand Down
Loading

0 comments on commit 79e433b

Please sign in to comment.