Skip to content

test: fix DelayParentScope test #1214

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 14 commits into from
Apr 10, 2025
Merged

test: fix DelayParentScope test #1214

Show file tree
Hide file tree
Changes from 9 commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
690ce43
test: fix `DelayParentScope` test
AayushSabharwal Apr 7, 2025
8e0bf3b
test: import `getname`
AayushSabharwal Apr 7, 2025
ddff0a3
test: remove reliance of tests on `Initial` parameters
AayushSabharwal Apr 9, 2025
bb374d7
fixup! test: remove reliance of tests on `Initial` parameters
AayushSabharwal Apr 9, 2025
ea38a83
spatial reactions fix
TorkelE Apr 9, 2025
8c30383
mark broken test as broken
TorkelE Apr 10, 2025
a52dd97
up
TorkelE Apr 10, 2025
4fb8a7c
disable jacobian tests
TorkelE Apr 10, 2025
aa24938
mark another broken test as broken
TorkelE Apr 10, 2025
de4506e
add jac eval
TorkelE Apr 10, 2025
cb5ee42
approx in jac evaluations
TorkelE Apr 10, 2025
21298b9
Update test/reactionsystem_core/reactionsystem.jl
TorkelE Apr 10, 2025
0e74125
Update test/reactionsystem_core/reactionsystem.jl
TorkelE Apr 10, 2025
86bd7c7
Update test_functions.jl
TorkelE Apr 10, 2025
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
18 changes: 11 additions & 7 deletions src/spatial_reaction_systems/spatial_reactions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,28 +103,32 @@ function check_spatial_reaction_validity(rs::ReactionSystem, tr::TransportReacti
if any(isequal(tr.species, s) && !isequivalent(tr.species, s) for s in species(rs))
error("A transport reaction used a species, $(tr.species), with metadata not matching its lattice reaction system. Please fetch this species from the reaction system and use it during transport reaction creation.")
end
# No `for` loop, just weird formatting by the formatter.
if any(isequal(rs_p, tr_p) && !isequivalent(rs_p, tr_p)
for rs_p in parameters(rs), tr_p in Symbolics.get_variables(tr.rate))
for rs_p in parameters(rs), tr_p in Symbolics.get_variables(tr.rate))
error("A transport reaction used a parameter with metadata not matching its lattice reaction system. Please fetch this parameter from the reaction system and use it during transport reaction creation.")
end

# Checks that no edge parameter occurs among rates of non-spatial reactions.
# No `for` loop, just weird formatting by the formatter.
if any(!isempty(intersect(Symbolics.get_variables(r.rate), edge_parameters))
for r in reactions(rs))
for r in reactions(rs))
error("Edge parameter(s) were found as a rate of a non-spatial reaction.")
end
end

# Since MTK's "isequal" ignores metadata, we have to use a special function that accounts for this.
# This is important because whether something is an edge parameter is defined in metadata.
# MTK potentially start adding there own metadata to system symbolic variables, to prevent breakage
# for this we now also have al list of `ignored_metadata` (which MTK will add to `ReactionSystem`)
# metadata. Long-term the solution is to permit the creation of spatial reactions within
# a `ReactionSystem` when it is created.
const ep_metadata = Catalyst.EdgeParameter => true
function isequivalent(sym1, sym2)
function isequivalent(sym1, sym2; ignored_metadata = [MT.SymScope])
isequal(sym1, sym2) || (return false)
if any((md1 != ep_metadata) && !(md1 in sym2.metadata) for md1 in sym1.metadata)
if any((md1 != ep_metadata) && (md1[1] ∉ ignored_metadata) && (md1 ∉ sym2.metadata)
for md1 in sym1.metadata)
return false
elseif any((md2 != ep_metadata) && !(md2 in sym1.metadata) for md2 in sym2.metadata)
elseif any((md2 != ep_metadata) && (md2[1] ∉ ignored_metadata) && (md2 ∉ sym1.metadata)
for md2 in sym2.metadata)
return false
elseif typeof(sym1) != typeof(sym2)
return false
Expand Down
10 changes: 5 additions & 5 deletions test/compositional_modelling/component_based_model_creation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@

# Fetch packages.
using Catalyst, LinearAlgebra, OrdinaryDiffEqTsit5, SciMLNLSolve, Test
using ModelingToolkit: nameof
using ModelingToolkit: nameof, getname

# Sets the default `t` to use.
t = default_t()
Expand Down Expand Up @@ -507,12 +507,12 @@ let
@variables x3(t) x4(t) x5(t)
x2 = ParentScope(x2)
x3 = ParentScope(ParentScope(x3))
x4 = DelayParentScope(x4, 2)
x4 = DelayParentScope(x4)
x5 = GlobalScope(x5)
@parameters p1 p2 p3 p4 p5
p2 = ParentScope(p2)
p3 = ParentScope(ParentScope(p3))
p4 = DelayParentScope(p4, 2)
p4 = DelayParentScope(p4)
p5 = GlobalScope(p5)
rxs = [Reaction(p1, nothing, [x1]), Reaction(p2, [x2], nothing),
D(x3) ~ p3, D(x4) ~ p4, D(x5) ~ p5]
Expand All @@ -530,11 +530,11 @@ let
sys3 = sys3 ∘ sys2
@test length(unknowns(sys3)) == 4
@test any(isequal(x3), unknowns(sys3))
@test any(isequal(x4), unknowns(sys3))
@test any(endswith("x4") ∘ string ∘ getname, unknowns(sys3))
@test length(species(sys3)) == 2
@test length(parameters(sys3)) == 4
@test any(isequal(p3), parameters(sys3))
@test any(isequal(p4), parameters(sys3))
@test any(endswith("p4") ∘ string ∘ getname, parameters(sys3))
sys4 = complete(sys3)
@test length(unknowns(sys3)) == 4
@test length(parameters(sys4)) == 5
Expand Down
13 changes: 6 additions & 7 deletions test/reactionsystem_core/reactionsystem.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,9 +118,7 @@ let
kvals = Float64.(1:length(k))
def_p = [k => kvals]
def_u0 = [A => 0.5, B => 1.0, C => 1.5, D => 2.0]
inits = (MT.Initial(A), MT.Initial(B), MT.Initial(C), MT.Initial(D)) .=> 0
defs = merge(Dict(def_p), Dict(def_u0))
fulldefs = merge(copy(defs), Dict(inits))

@named rs = ReactionSystem(rxs, t, [A, B, C, D], [k]; defaults = defs)
rs = complete(rs)
Expand All @@ -133,7 +131,8 @@ let

# these systems add initial conditions to the defaults
@test ModelingToolkit.get_defaults(odesys) ==
ModelingToolkit.get_defaults(sdesys) == fulldefs
ModelingToolkit.get_defaults(sdesys)
@test issubset(defs, ModelingToolkit.get_defaults(odesys))

u0map = [A => 5.0]
kvals[1] = 5.0
Expand Down Expand Up @@ -465,9 +464,9 @@ let
rs = complete(rs)
@test all(eq -> eq isa Reaction, ModelingToolkit.get_eqs(rs)[1:4])
osys = complete(convert(ODESystem, rs))
inits = Initial.((B, C, D, E))
@test issetequal(MT.get_unknowns(osys), [B, C, D, E])
@test issetequal(MT.get_ps(osys), [k1, k2, A, inits...])
_ps = filter(x -> !iscall(x) || !(operation(x) isa Initial), MT.get_ps(osys))
@test issetequal(_ps, [k1, k2, A])

# test nonlinear systems
u0 = [1.0, 2.0, 3.0, 4.0]
Expand Down Expand Up @@ -497,9 +496,9 @@ let
@named rs = ReactionSystem(rxs, t) # add constraint csys when supported!
rs = complete(rs)
ssys = complete(convert(SDESystem, rs))
inits = Initial.((B, C, D, E))
@test issetequal(MT.get_unknowns(ssys), [B, C, D, E])
@test issetequal(MT.get_ps(ssys), [A, k1, k2, inits...])
_ps = filter(x -> !iscall(x) || !(operation(x) isa Initial), MT.get_ps(ssys))
@test issetequal(_ps, [A, k1, k2])
du1 = zeros(4)
du2 = zeros(4)
sprob = SDEProblem(ssys, u0map, tspan, pmap; check_length = false)
Expand Down
2 changes: 1 addition & 1 deletion test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,7 @@ end

# Tests ODE, SDE, jump simulations, nonlinear solving, and steady state simulations.
@time @safetestset "ODE System Simulations" begin include("simulation_and_solving/simulate_ODEs.jl") end
@time @safetestset "Automatic Jacobian Construction" begin include("simulation_and_solving/jacobian_construction.jl") end
#@time @safetestset "Automatic Jacobian Construction" begin include("simulation_and_solving/jacobian_construction.jl") end
Copy link
Member

Choose a reason for hiding this comment

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

Is the CI failure just from testing floating point results are equal instead of agreeing to some tolerance? If so just update the test to test agreement to some reasonable tolerance. We generally shouldn’t test floating point results are exactly equal anyways.

Copy link
Member

Choose a reason for hiding this comment

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

So I think doing approx is a sensible solution. There is a discussion over at SciML/ModelingToolkit.jl#3554 (comment) and there might be some reason why the MTK people would want to have a look. From our purposes, using approx should be more than fine.

Copy link
Member

Choose a reason for hiding this comment

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

Yeah with how SciML/ModelingToolkit.jl#3554 was narrowed down, I thought at first it could be something a bit more fundamental but it was narrowed down to CSE, which is an expression reordering to save computations thing. For reference, common subexpression elimination is

x = (x + 1) * y + (x + 1) * z

you instead do

tmp = x + 1
x = tmp*y + tmp*z

Symbolics.jl build function now supports doing this, and MTK turns it on by default. You can see it changes the BCR sparse Jacobian build from 200 seconds to 20 seconds:

───────────────────────────────────────────────────────────────────────────
───────────────────────────────────
                                                                     Time  
                  Allocations      
                                                            ───────────────
────────   ────────────────────────
                     Tot / % measured:                            500s /  9
3.2%           34.4GiB /  90.1%    

Section                                             ncalls     time    %tot
     avg     alloc    %tot      avg
───────────────────────────────────────────────────────────────────────────
───────────────────────────────────
Compile jacobian - no CSE                                1     213s   45.7%
    213s   3.05GiB    9.8%  3.05GiB
Compile jacobian - CSE                                   1    77.8s   16.7%
   77.8s   0.96GiB    3.1%  0.96GiB
Calculate jacobian - without hashconsing                 1    73.3s   15.7%
   73.3s   10.5GiB   34.0%  10.5GiB
Calculate jacobian - hashconsing, without caching        1    72.1s   15.5%
   72.1s   10.5GiB   33.8%  10.5GiB
Calculate jacobian - hashconsing and caching             1    23.7s    5.1%
   23.7s   5.19GiB   16.8%  5.19GiB
Build jacobian - no CSE                                  1    4.42s    0.9%
   4.42s    615MiB    1.9%   615MiB
Build jacobian - CSE                                     1    1.65s    0.4%
   1.65s    162MiB    0.5%   162MiB
Compute jacobian - no CSE                                1   91.4μs    0.0%
  91.4μs      272B    0.0%     272B
Compute jacobian - CSE                                   1   46.6μs    0.0%
  46.6μs      272B    0.0%     272B
───────────────────────────────────────────────────────────────────────────
───────────────────────────────────

in the unreleased benchmarks SciML/SciMLBenchmarks.jl#1178

But yes, it is a floating point level difference and I don't think Catalyst tests should care about a 1e-16 from a reordering of arithmetic, but at the same time I do want to follow up with this in MTK just to find out why the in-place and out of place are not giving exactly the same ordering post LLVM optimizations.

Copy link
Member

Choose a reason for hiding this comment

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

I had seen the recent CSE stuff and figured that was the cause, but yes we shouldn’t generally expect floating point equality except for results from the exact same function. I figured it could be that something was slightly different in the codegen for in and out of place.

@time @safetestset "SDE System Simulations" begin include("simulation_and_solving/simulate_SDEs.jl") end
@time @safetestset "Jump System Simulations" begin include("simulation_and_solving/simulate_jumps.jl") end
@time @safetestset "Nonlinear and SteadyState System Solving" begin include("simulation_and_solving/solve_nonlinear.jl") end
Expand Down
2 changes: 1 addition & 1 deletion test/simulation_and_solving/jacobian_construction.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,7 @@ let
return J
end
@test eval_jac(oprob_jac, false) == eval_jac(sprob_jac, false) == eval_jac(nlprob_jac, false)
@test_broken eval_jac(oprob_sjac, true) == eval_jac(sprob_sjac, true) == eval_jac(nlprob_sjac, true) # https://github.com/SciML/ModelingToolkit.jl/issues/3527
@test eval_jac(oprob_sjac, true) == eval_jac(sprob_sjac, true) == eval_jac(nlprob_sjac, true)
end
end

Expand Down
2 changes: 1 addition & 1 deletion test/simulation_and_solving/simulate_ODEs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -193,7 +193,7 @@ let
ps = [:k1 => 2.0f0, :k2 => 3.0f0]
oprob = ODEProblem(rn, u0, 1.0f0, ps)
osol = solve(oprob, Tsit5())
@test eltype(osol[:X1]) == eltype(osol[:X2]) == typeof(oprob[:X1]) == typeof(oprob[:X2]) == Float32
@test_broken eltype(osol[:X1]) == eltype(osol[:X2]) == typeof(oprob[:X1]) == typeof(oprob[:X2]) == Float32 # https://github.com/SciML/ModelingToolkit.jl/issues/3553
@test eltype(osol.t) == typeof(oprob.tspan[1]) == typeof(oprob.tspan[2]) == Float32
end

Expand Down
6 changes: 3 additions & 3 deletions test/simulation_and_solving/simulate_SDEs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -384,7 +384,7 @@ end
### Other Tests ###

# Checks that solution values have types consistent with their input types.
# Check that both float types are preserved in the solution (and problems), while integers are
# Check that both float types are preserved in the solution (and problems), while integers are
# promoted to floats.
# Checks that the time types are correct (`Float64` by default or possibly `Float32`), however,
# type conversion only occurs in the solution, and integer types are preserved in problems.
Expand All @@ -400,7 +400,7 @@ let
@test eltype(ssol[:X1]) == eltype(ssol[:X2]) == typeof(sprob[:X1]) == typeof(sprob[:X2]) == Float64
@test eltype(ssol.t) == typeof(sprob.tspan[1]) == typeof(sprob.tspan[2]) == Float64

# Checks that `Int64` values are promoted to `Float64`.
# Checks that `Int64` values are promoted to `Float64`.
u0 = [:X1 => 1, :X2 => 3]
ps = [:k1 => 2, :k2 => 3]
sprob = SDEProblem(rn, u0, 1, ps)
Expand All @@ -413,7 +413,7 @@ let
ps = [:k1 => 2.0f0, :k2 => 3.0f0]
sprob = SDEProblem(rn, u0, 1.0f0, ps)
ssol = solve(sprob, ISSEM())
@test eltype(ssol[:X1]) == eltype(ssol[:X2]) == typeof(sprob[:X1]) == typeof(sprob[:X2]) == Float32
@test_broken eltype(ssol[:X1]) == eltype(ssol[:X2]) == typeof(sprob[:X1]) == typeof(sprob[:X2]) == Float32
@test eltype(ssol.t) == typeof(sprob.tspan[1]) == typeof(sprob.tspan[2]) == Float32
end

Expand Down
Loading