Clean up and test compare interface

Author: BenChung

src/test/compare.jl

@@ -6,23 +6,28 @@ struct DefaultComparison
                 :L∞ => (delta, t) -> norm.(delta, Inf), 
                 :L1 => (delta, t) -> norm.(delta, 1), 
                 :L2 => (delta, t) -> norm.(delta, 2), 
-                :rms => (delta, t) -> sqrt.(1/length(t) .* sum.(map(d-> d .^ 2, delta)))]))
+                :rms => (delta, t) -> sqrt.(1/length(t) .* sum.(map(d-> d .^ 2, delta))),
+                :final => (delta, t) -> last.(delta)]))
         end
         return new(field_cmp)
     end
 end
-function (d::DefaultComparison)(c, b, t, n, r) 
+function (d::DefaultComparison)(c, names, b, t, n, r) 
     delta = map((o, re) -> o .- re, n, r)
     cmps = [name => cmper(delta, t) for (name, cmper) in d.field_cmp]
-    return DataFrame([:var => b, cmps...,
+    return DataFrame([:var => names, cmps...,
         :observed => SymbolicIndexingInterface.is_observed.(c, b)])
 end
 
 
 function compare(
-    new_sol::SciMLBase.AbstractTimeseriesSolution, basis::Vector, 
-    reference::DataFrame, reference_basis::Vector,
-    cmp; warn_observed=true)
+    new_sol::SciMLBase.AbstractTimeseriesSolution, 
+    reference::DataFrame,
+    over::Vector{<:Union{Pair{<:Any, String}, Pair{<:Any, Pair{String, String}}}},
+    cmp=DefaultComparison(); warn_observed=true)
+    basis = first.(over)
+    reference_basis = map(e -> e isa Pair ? first(e) : e, last.(over))
+    output_names = map((n1, n2) -> n2 isa Pair ? last(n2) : string(n1), basis, last.(over))
     @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"
@@ -31,7 +36,7 @@ function compare(
         dat = new_sol(reference[:, :timestamp], idxs=basis)
         obs = [[dat[j][i] for j=1:nrow(reference)] for i in eachindex(basis)]
         ref = collect.(eachcol(reference[:, Not(:timestamp)]))
-        cmp(new_container, basis, reference[:, :timestamp], obs, ref)
+        cmp(new_container, output_names, basis, reference[:, :timestamp], obs, ref)
     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
@@ -54,6 +59,6 @@ function compare(
         end
         matching_sym 
     end for ref_name in setdiff(names(reference), ["timestamp"])] 
-    return compare(new_sol, eval_syms, reference, setdiff(names(reference), ["timestamp"]), cmp)
+    return compare(new_sol, reference, eval_syms .=> setdiff(names(reference), ["timestamp"]), cmp)
 end
 export compare

test/block_modeling.jl

@@ -28,23 +28,26 @@ using SymbolicIndexingInterface
         sys = structural_simplify(rc_model)
         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)
-        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))
+        prob2 = ODEProblem(sys, Pair[capacitor.C => 0.9], (0, 10.0)) # this converges to nearly the same solution but is a little too fast
         sol2 = solve(prob2, Tsit5())
-        prob3 = ODEProblem(sys, Pair[capacitor.C => 5.0], (0, 10.0))
+        prob3 = ODEProblem(sys, Pair[capacitor.C => 5.0], (0, 10.0)) # this doesn't stabilize in the allotted time
         sol3 = solve(prob3, Tsit5())
 
         d1 = discretize_solution(sol1, sol1)
-
         ds1 = discretize_solution(sol1, sol1; measured=SymbolicIndexingInterface.all_variable_symbols(sys))
-        println(compare(sol3, ds1))
-        println(compare(sol2, ds1))
-        println(compare(sol1, ds1))
+
+        @test sum(compare(sol1, ds1; warn_observed=false)[:, :L∞]) < 0.01
+        @test sum(compare(sol2, ds1; warn_observed=false)[:, :L∞]) > sum(compare(sol1, ds1; warn_observed=false)[:, :L∞])
+        @test sum(compare(sol3, ds1; warn_observed=false)[:, :L∞]) > sum(compare(sol2, ds1; warn_observed=false)[:, :L∞])
+        @test sum(compare(sol1, ds1; warn_observed=false)[:, :L∞]) > sum(compare(sol1, d1; warn_observed=false)[:, :L∞])
+        @test sum(compare(sol2, ds1; warn_observed=false)[:, :L∞]) > sum(compare(sol2, d1; warn_observed=false)[:, :L∞])
+        @test sum(compare(sol3, ds1; warn_observed=false)[:, :L∞]) > sum(compare(sol3, d1; warn_observed=false)[:, :L∞])
+
+        # construct a fictional power measurement
+        power_synth = select(ds1, :timestamp, ["capacitor₊v(t)", "capacitor₊i(t)"] => ((v, i) -> v .* i) => "power")
+        @test compare(sol1, power_synth, [capacitor.i * capacitor.v => "power" => "power"])[1, "L∞"] < 0.01
+        @test compare(sol2, power_synth, [capacitor.i * capacitor.v => "power" => "power"])[1, "L∞"] < 0.05
+        @test compare(sol3, power_synth, [capacitor.i * capacitor.v => "power" => "power"])[1, "L∞"] < 0.3
+        
     end
 end