Implement compare

Author: BenChung

Project.toml

@@ -22,6 +22,7 @@ ModelTestingCalibration = "JuliaSimModelOptimizer"
 
 [compat]
 julia = "1.10"
+SymbolicIndexingInterface = "^0.3.15"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/ModelTesting.jl

@@ -8,6 +8,8 @@ include("test/continuous/delta_sol.jl")
 include("test/discrete/single_shooting.jl")
 include("test/discrete/merge.jl")
 include("test/metric/metric.jl")
+include("test/discrete/compare.jl")
+include("test/metric/spec.jl")
 include("test/instantaneous/instant.jl")
 include("problem_config/problem.jl")
 end

src/test/continuous/delta_sol.jl

@@ -17,7 +17,7 @@ function make_cols(names, rows)
 end
 
 function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution, time_ref::SciMLBase.AbstractTimeseriesSolution; all_observed=false)
-    container = symbolic_container(time_ref)
+    container = SymbolicIndexingInterface.symbolic_container(time_ref)
     ref_t_vars = independent_variable_symbols(container)
     if length(ref_t_vars) > 1
         @error "PDE solutions not currently supported; only one iv is allowed"
@@ -29,15 +29,15 @@ function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution, tim
     return discretize_solution(solution, collect(time_ref[!, "timestamp"]); all_observed=all_observed)
 end
 function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution; all_observed=false)
-    container = symbolic_container(solution)
+    container = SymbolicIndexingInterface.symbolic_container(solution)
     ref_t_vars = independent_variable_symbols(container)
     if length(ref_t_vars) > 1
         @error "PDE solutions not currently supported; only one iv is allowed"
     end
-    return discretize_solution(solution, solution[ref_t_var]; all_observed=all_observed )
+    return discretize_solution(solution, solution[first(ref_t_vars)]; all_observed=all_observed )
 end
 function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution, time_ref::AbstractArray; measured=nothing, all_observed=false)
-    container = symbolic_container(solution)
+    container = SymbolicIndexingInterface.symbolic_container(solution)
     if isnothing(measured)
         if all_observed
             measured = all_variable_symbols(container)
@@ -72,8 +72,8 @@ function compare_dense_solutions(
     integrator=Tsit5()
 )
     results = Dict{Symbol, Any}()
-    reference_container = symbolic_container(reference)
-    containers = symbolic_container(sol)
+    reference_container = SymbolicIndexingInterface.symbolic_container(reference)
+    containers = SymbolicIndexingInterface.symbolic_container(sol)
 
     measured_reference = isnothing(reference_measured) ? measured_values(reference_container) : reference_measured
     sol_measured = isnothing(solution_measured) ? measured_values(containers) : sol_measured

src/test/discrete/compare.jl

@@ -0,0 +1,46 @@
+function compare(
+    new_sol::SciMLBase.AbstractTimeseriesSolution, basis::Vector, 
+    reference::DataFrame, reference_basis::Vector,
+    cmp::Function, red::Function; init=nothing)
+    @assert "timestamp" ∈ names(reference) "The dataset must contain a column named `timestamp`"
+    new_container = SymbolicIndexingInterface.symbolic_container(new_sol)
+    @assert all(SymbolicIndexingInterface.is_observed.((new_container, ), basis) .| SymbolicIndexingInterface.is_variable.((new_container, ), basis)) "All basis symbols must be observed in the new system"
+    @assert all(b ∈ names(reference) for b in reference_basis) "The reference basis must be a subset of the columns in the reference data"
+    if new_sol.dense
+        foldl(red, cmp(row[:timestamp], new_sol(row[:timestamp], idxs=basis), row[reference_basis]) for row in eachrow(reference); init=isnothing(init) ? nothing : init(basis))
+    else 
+        foldl(red, cmp(row[:timestamp], new_sol[row[:timestamp], basis], row[reference_basis]) for row in eachrow(reference); init=isnothing(init) ? nothing : init(basis))
+    end
+end
+
+function compare(
+    new_sol::SciMLBase.AbstractTimeseriesSolution,
+    reference::DataFrame,
+    cmp::Function, red::Function; init=nothing, to_name=string, warn_observed=true)
+    @assert "timestamp" ∈ names(reference) "The dataset must contain a column named `timestamp`"
+
+    new_container = SymbolicIndexingInterface.symbolic_container(new_sol)
+    all_available = SymbolicIndexingInterface.all_variable_symbols(new_container)
+    all_available_syms = Dict(to_name.(all_available) .=> all_available)
+    eval_syms = [ begin
+        @assert ref_name ∈ keys(all_available_syms) "Reference value $ref_name not found in the new solution (either as state or observed)"
+        matching_sym = all_available_syms[ref_name]
+        if warn_observed && SymbolicIndexingInterface.is_observed(new_container, matching_sym)
+            @warn "The variable $matching_sym is observed in the new solution while the past $ref_name is provided explicitly; problem structure may have changed"
+        end
+        matching_sym 
+    end for ref_name in setdiff(names(reference), ["timestamp"])] 
+    return compare(new_sol, eval_syms, reference, setdiff(names(reference), ["timestamp"]), cmp, red; init=init)
+end
+
+function compare(
+    new_sol::SciMLBase.AbstractTimeseriesSolution,
+    reference::DataFrame; warn_observed=true)
+    @assert "timestamp" ∈ names(reference) "The dataset must contain a column named `timestamp`"
+    return compare(new_sol, reference, 
+        (t, n, r) -> abs.(collect(n) .- collect(r)), 
+        (acc, nv) -> begin acc[:, :l∞] = max.(collect(acc[:, :l∞]), collect(nv)); acc end; 
+        init=(v) -> DataFrame(:name => v, :l∞ => zeros.(length(v))), 
+        warn_observed=warn_observed)
+end
+export compare

test/block_modeling.jl

@@ -65,17 +65,24 @@ using ModelingToolkitStandardLibrary.Blocks: Constant
         prob1 = ODEProblem(sys, Pair[], (0, 10.0))
         sol1 = solve(prob1, Tsit5())
 
+        #=
         ref = DataFrame(:timestamp => sol1.t, Symbol(resistor.v) => sol1[resistor.v])
-        measure, measured_sys = compare_data(resistor.v, ref)(sys)
+        measure, measured_sys = compare_data([resistor.v], ref)(sys)
         prob_measured = ODEProblem(measured_sys, Pair[], (0, 10.0))
         sol_measured = solve(prob_measured, Tsit5())
-
+=#
         prob2 = ODEProblem(sys, Pair[capacitor.C => 0.9], (0, 10.0))
         sol2 = solve(prob2, Tsit5())
         prob3 = ODEProblem(sys, Pair[capacitor.C => 5.0], (0, 10.0))
         sol3 = solve(prob3, Tsit5())
 
         d1 = discretize_solution(sol1, sol1)
+
+        ds1 = discretize_solution(sol1, sol1; measured=states(sys))
+        println(compare(sol3, ds1))
+        println(compare(sol2, ds1))
+        println(compare(sol1, ds1))
+
         d2 = discretize_solution(sol2, sol1)
         d3 = discretize_solution(sol3, sol1)
         results_good = compare_discrete(sys, d1, d2)