Implement continous-time solution comparison by integration measure(ref - sol)

Author: BenChung

src/test/continous/delta_sol.jl

@@ -71,6 +71,51 @@ function compare_solutions(
 end
 export compare_solutions
 
+function compare_dense_solutions(
+    (ref_name, reference)::Pair{Symbol, <:SciMLBase.AbstractTimeseriesSolution},
+    sols::Vector{<:Pair{Symbol, <:SciMLBase.AbstractTimeseriesSolution}};
+    integrator=Tsit5(),
+    metric=abs
+)
+    results = Dict{Symbol, Any}()
+    reference_container = symbolic_container(reference)
+    containers = symbolic_container.(last.(sols))
+
+    measured_reference = measured_values(reference_container)
+    sols_measured = measured_values.(containers)
+    @assert all(_symbolic_subset.((measured_reference,), sols_measured)) "Test solutions must expose a superset of the reference's variables for comparison"
+    @assert length(measured_reference) > 0 "Compared solutions must share at least one measured variable"
+    measured = measured_reference
+
+    timebounds(sol) = (sol.t[1], sol.t[end])
+    @assert reference.dense "Dense (integrated) comparision requires a dense reference solution"
+    for (test_name, test_sol) in sols
+        @assert test_sol.dense "Test solution $(test_name) must be dense in order to use continous-time comparison"
+        @assert timebounds(test_sol) == timebounds(reference) "Test solution $(test_name) has time range $(timebounds(test_sol)) which differs from the reference $(timebounds(reference))"
+    end
+
+
+    ref_t_vars = independent_variable_symbols(reference_container)
+    if length(ref_t_vars) > 1
+        @error "PDE solutions not currently supported; only one iv is allowed"
+    end
+    ref_t_var = first(ref_t_vars) 
+
+    state_size = length(measured)
+    function compare!(du, u, p, t) 
+        offs = 1
+        for (test_name, test_sol) in sols
+            du[offs:offs+state_size-1] .= abs.(reference(t, idxs=measured) .- test_sol(t, idxs=measured))
+            offs += state_size
+        end
+    end
+    func = ODEFunction(compare!; sys = SymbolCache(collect(Iterators.flatten([namespace_symbol.((test_name, ), measured) for (test_name, _) in sols])), [], ref_t_var))
+    prob = ODEProblem(func, zeros(length(sols) * length(measured)), timebounds(reference))
+    soln = solve(prob, integrator)
+    return soln
+end
+export compare_dense_solutions
+
 function compute_error_metrics(output_results, ref_soln, test_soln)
     delta = ref_soln - test_soln
     output_results[:l∞] = maximum(vecvecapply((x) -> abs.(x), delta))

test/timeseries.jl

@@ -2,7 +2,7 @@ Base.isapprox((v1, t1)::Tuple{Vector{Float64}, Float64}, (v2, t2)::Tuple{Vector{
 tol = 1e-8
 
 @testset "Continous" begin 
-    @testset "Model-Model Comparison" begin 
+    @testset "Model-Model Mixed (self-discretized) Comparison" begin 
         using ModelingToolkit: t_nounits as t, D_nounits as D
         @variables x(t) [measured=true] y(t) [measured=true]
         @parameters τ
@@ -35,4 +35,23 @@ tol = 1e-8
         @test results[:metrics][:bad][:l2] > tol
         @test results[:metrics][:bad][:final] > tol
     end
+    @testset "Model-Model Continous Comparison" begin 
+        using ModelingToolkit: t_nounits as t, D_nounits as D
+        @variables x(t) [measured=true] y(t) [measured=true]
+        @parameters τ
+        @named fol_sys = ODESystem([D(x) ~ (1 - y) / τ; y ~ x], t)
+        fol = structural_simplify(fol_sys)
+        prob1 = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
+        sol1 = solve(prob1, Tsit5(), reltol = 1e-8, abstol = 1e-8)
+        prob2 = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 3.0])
+        sol2 = solve(prob2, Tsit5(), reltol = 1e-8, abstol = 1e-8)
+        prob3 = ODEProblem(fol, [fol.x => 0.0], (0.0, 10.0), [fol.τ => 1.0])
+        sol3 = solve(prob3, Tsit5(), reltol = 1e-8, abstol = 1e-8)
+
+        results = compare_dense_solutions(:reference=>sol1, [:good=>sol2, :bad=>sol3])
+        @test results[Symbol("good/x(t)"), end] < tol
+        @test results[Symbol("good/y(t)"), end] < tol
+        @test results[Symbol("bad/x(t)"), end] > tol
+        @test results[Symbol("bad/y(t)"), end] > tol
+    end
 end