Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
45 changes: 45 additions & 0 deletions compare.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,45 @@
using Revise, StructuralIdentifiability, ParamPunPam, RationalFunctionFields
import StructuralIdentifiability: parent_ring_change, lie_derivative
import RationalFunctionFields: fields_equal

cmp = (f, g) -> begin
f = parent_ring_change(f, parent(ode))
g = parent_ring_change(g, parent(ode))
df = lie_derivative(f, ode)
dg = lie_derivative(g, ode)
default_cmp = RationalFunctionFields.rational_function_cmp
if iszero(df) && iszero(dg)
default_cmp(f, g)
else
default_cmp(df, dg)
end

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[Runic] reported by reviewdog 🐶

Suggested change
end
end

end

ode = @ODEmodel(
x1'(t) =
(p1 * x4(t)) - (p3 * x1(t)) -
p4 * ((x1(t)^2 / (p5 + x1(t))) * (1 + (p6 * u1(t) / (p7 + u1(t))))),
x2'(t) =
p8 - (p9 * x2(t)) -
p10 * (
(x1(t) * x2(t) / (p11 + x2(t))) * (1 + (p12 * u1(t) / (p13 + u1(t))))
),
x3'(t) =
p14 - (p15 * x3(t)) -
p16 * x1(t) * x3(t) * (1 - p18 * u1(t)) / (p17 + x3(t)),
x4'(t) =
p20 - p21 * (1 - p24) * (1 - p25) / ((p22^4) + 1) - (p20 * x4(t)) +
(p21 * (x3(t)^4)) +
(1 + p23 * u1(t)) * (1 - p24 * x1(t)) * (1 - p25 * x2(t)) /
(p22^4 + x3(t)^4),
y1(t) = x1(t),
y2(t) = x2(t),
y3(t) = x3(t),
y4(t) = x4(t)
)
Comment on lines +19 to +39

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[Runic] reported by reviewdog 🐶

Suggested change
x1'(t) =
(p1 * x4(t)) - (p3 * x1(t)) -
p4 * ((x1(t)^2 / (p5 + x1(t))) * (1 + (p6 * u1(t) / (p7 + u1(t))))),
x2'(t) =
p8 - (p9 * x2(t)) -
p10 * (
(x1(t) * x2(t) / (p11 + x2(t))) * (1 + (p12 * u1(t) / (p13 + u1(t))))
),
x3'(t) =
p14 - (p15 * x3(t)) -
p16 * x1(t) * x3(t) * (1 - p18 * u1(t)) / (p17 + x3(t)),
x4'(t) =
p20 - p21 * (1 - p24) * (1 - p25) / ((p22^4) + 1) - (p20 * x4(t)) +
(p21 * (x3(t)^4)) +
(1 + p23 * u1(t)) * (1 - p24 * x1(t)) * (1 - p25 * x2(t)) /
(p22^4 + x3(t)^4),
y1(t) = x1(t),
y2(t) = x2(t),
y3(t) = x3(t),
y4(t) = x4(t)
)
x1'(t) =
(p1 * x4(t)) - (p3 * x1(t)) -
p4 * ((x1(t)^2 / (p5 + x1(t))) * (1 + (p6 * u1(t) / (p7 + u1(t))))),
x2'(t) =
p8 - (p9 * x2(t)) -
p10 * (
(x1(t) * x2(t) / (p11 + x2(t))) * (1 + (p12 * u1(t) / (p13 + u1(t))))
),
x3'(t) =
p14 - (p15 * x3(t)) -
p16 * x1(t) * x3(t) * (1 - p18 * u1(t)) / (p17 + x3(t)),
x4'(t) =
p20 - p21 * (1 - p24) * (1 - p25) / ((p22^4) + 1) - (p20 * x4(t)) +
(p21 * (x3(t)^4)) +
(1 + p23 * u1(t)) * (1 - p24 * x1(t)) * (1 - p25 * x2(t)) /
(p22^4 + x3(t)^4),
y1(t) = x1(t),
y2(t) = x2(t),
y3(t) = x3(t),
y4(t) = x4(t)
)


id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[Runic] reported by reviewdog 🐶

Suggested change
id1 = find_identifiable_functions(ode, with_states=true, cmp = cmp)
id1 = find_identifiable_functions(ode, with_states = true, cmp = cmp)


id2 = find_identifiable_functions(ode, with_states=true, cmp = RationalFunctionFields.rational_function_cmp)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[Runic] reported by reviewdog 🐶

Suggested change
id2 = find_identifiable_functions(ode, with_states=true, cmp = RationalFunctionFields.rational_function_cmp)
id2 = find_identifiable_functions(ode, with_states = true, cmp = RationalFunctionFields.rational_function_cmp)


@assert fields_equal(RationalFunctionField(id1), RationalFunctionField(id2), 0.99)
4 changes: 3 additions & 1 deletion src/global_identifiability.jl
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ The function returns a tuple containing the following:
with_states::Bool = false,
var_change_policy = :default,
rational_interpolator = :VanDerHoevenLecerf,
cmp = RationalFunctionFields.rational_function_cmp
) where {T}
@info "Computing IO-equations"
ioeq_time = @elapsed io_equations =
Expand Down Expand Up @@ -148,6 +149,7 @@ The function returns a tuple containing the following:
seed = 42,
simplify = :standard,
rational_interpolator = rational_interpolator,
cmp = cmp
)
dennums_simplified = fractions_to_dennums(no_states_simplified)
# switch back the ring
Expand Down Expand Up @@ -202,7 +204,7 @@ Output: a list L of booleans with L[i] being the identifiability status of the i
known = known,
prob_threshold = half_p,
var_change_policy = var_change_policy,
with_states = states_needed,
with_states = states_needed
)

funcs_to_check = Vector{Generic.FracFieldElem{P}}(
Expand Down
8 changes: 7 additions & 1 deletion src/identifiable_functions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ function find_identifiable_functions(
simplify = :standard,
rational_interpolator = :VanDerHoevenLecerf,
loglevel = Logging.Info,
cmp = RationalFunctionFields.rational_function_cmp,
) where {T <: MPolyRingElem{QQFieldElem}}
restart_logging(loglevel = loglevel)
reset_timings()
Expand All @@ -66,6 +67,7 @@ function find_identifiable_functions(
with_states = with_states,
simplify = simplify,
rational_interpolator = rational_interpolator,
cmp = cmp
)
else
id_funcs = _find_identifiable_functions_kic(
Expand All @@ -75,6 +77,7 @@ function find_identifiable_functions(
seed = seed,
simplify = simplify,
rational_interpolator = rational_interpolator,
cmp = cmp
)
# renaming variables from `x(t)` to `x(0)`
return replace_with_ic(ode, id_funcs)
Expand All @@ -89,6 +92,7 @@ function _find_identifiable_functions(
with_states = false,
simplify = :standard,
rational_interpolator = :VanDerHoevenLecerf,
cmp = RationalFunctionFields.rational_function_cmp
) where {T <: MPolyRingElem{QQFieldElem}}
Random.seed!(seed)
@assert simplify in (:standard, :weak, :strong, :absent)
Expand All @@ -109,6 +113,7 @@ function _find_identifiable_functions(
prob_threshold = half_p,
with_states = with_states,
rational_interpolator = rational_interpolator,
cmp = cmp
)
# If simplification is needed
if simplify !== :absent
Expand All @@ -122,7 +127,8 @@ function _find_identifiable_functions(
seed = seed,
simplify = simplify,
rational_interpolator = rational_interpolator,
priority_variables = [parent_ring_change(p, bring) for p in ode.parameters],
# priority_variables = [parent_ring_change(p, bring) for p in ode.parameters],
cmp = cmp
)
else
id_funcs_fracs = RationalFunctionFields.dennums_to_fractions(id_funcs)
Expand Down
3 changes: 3 additions & 0 deletions src/known_ic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ function _find_identifiable_functions_kic(
seed = 42,
simplify = :standard,
rational_interpolator = :VanDerHoevenLecerf,
cmp = RationalFunctionFields.rational_function_cmp
) where {T <: MPolyRingElem{Nemo.QQFieldElem}}
Random.seed!(seed)
@assert simplify in (:standard, :weak, :strong, :absent)
Expand All @@ -45,6 +46,7 @@ function _find_identifiable_functions_kic(
simplify = :absent,
rational_interpolator = rational_interpolator,
seed = seed,
cmp = cmp
)

id_funcs = vcat(id_funcs_general, [f // one(parent(ode)) for f in known_ic])
Expand All @@ -56,6 +58,7 @@ function _find_identifiable_functions_kic(
seed = seed,
simplify = simplify,
rational_interpolator = rational_interpolator,
cmp = cmp
)
end

Expand Down
Loading