Refactor continous discretization outside of specific comparison methods

Author: BenChung

ext/ModelTestingCalibration.jl

@@ -15,7 +15,6 @@ function ModelTesting.validate(model::AbstractTimeDependentSystem, data; pem_coe
     else
         model_transformations = []
     end
-    println(experiment_kwargs)
     experiment = Experiment(data, model; 
         model_transformations = model_transformations, 
         filter(arg->first(arg) != :model_transformations, experiment_kwargs)...)

src/test/continuous/delta_sol.jl

@@ -16,101 +16,79 @@ function make_cols(names, rows)
     return cols
 end
 
-function compare_solutions(
-    (ref_name, reference)::Pair{Symbol, <:SciMLBase.AbstractTimeseriesSolution},
-    sols::Vector{<:Pair{Symbol, <:SciMLBase.AbstractTimeseriesSolution}};
-    knots=nothing)
-    if isnothing(knots)
-        knots = reference.t
+function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution, time_ref::SciMLBase.AbstractTimeseriesSolution)
+    container = 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"
     end
-    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
-
-    if knots != reference.t && !reference.dense
-        @assert reference.dense "Interpolated evaluation points require a dense reference solution"
+    return discretize_solution(solution, time_ref[first(ref_t_vars)])
+end
+function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution, time_ref::DataFrame)
+    @assert "timestamp" ∈ names(time_ref) "The dataset B must contain a column named `timestamp`"
+    return discretize_solution(solution, collect(time_ref[!, "timestamp"]))
+end
+function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution)
+    container = 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
-
-    ref_t_vars = independent_variable_symbols(reference_container)
+    return discretize_solution(solution, solution[ref_t_var] )
+end
+function discretize_solution(solution::SciMLBase.AbstractTimeseriesSolution, time_ref::AbstractArray)
+    container = symbolic_container(solution)
+    measured = measured_values(container)
+    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
     ref_t_var = first(ref_t_vars)
-    if knots == reference.t
-        ref_soln = reference[measured]
-        cols = make_cols([Symbol(ref_t_var); namespace_symbol.((ref_name,), measured)], reference[[ref_t_var; measured]])
-    else 
-        ref_soln = reference(knots, idxs=measured)
-        cols = make_cols([Symbol(ref_t_var); namespace_symbol.((ref_name,), measured)], reference(knots, idxs=[ref_t_var; measured]))
-    end
-    testsols = Dict{Symbol, Any}()
-    results[:metrics] = testsols
-    for (name, test) in sols 
-        test_results = Dict{Symbol, Any}()
-        if test.dense
-            test_soln = test(knots, idxs=measured)
-        elseif knots == sol.t
-            test_soln = test[measured]
-        else 
-            @warn "Comparison of $ref_name to $name is invalid (cannot match timebases); only computing final state error"
-            continue
-        end
-        test_results[:final] = recursive_mean(abs.(ref_soln[end] - test_soln[end]))
-        compute_error_metrics(test_results, ref_soln, test_soln)
-        append!(cols, make_cols(namespace_symbol.((name,), measured), test_soln))
-        testsols[name] = test_results
+    
+    matching_timebase = solution[ref_t_var] == time_ref
+
+    if matching_timebase # if the time_ref match the timebase of the problem then use the value at the nodes regardless of if it's dense or sparse
+        cols = make_cols(String["timestamp"; measured_names(measured)], solution[[ref_t_var; measured]])
+    elseif solution.dense # continious-time solution, use the interpolant    
+        cols = make_cols(String["timestamp"; measured_names(measured)], solution(time_ref, idxs=[ref_t_var; measured]))
+    else
+        throw("Cannot discretize_solution sparse solution about a different timebase.")
     end
-    results[:data] = DataFrame(cols)
-    return results
+    return DataFrame(cols)
 end
-export compare_solutions
+export discretize_solution
 
 function compare_dense_solutions(
-    (ref_name, reference)::Pair{Symbol, <:SciMLBase.AbstractTimeseriesSolution},
-    sols::Vector{<:Pair{Symbol, <:SciMLBase.AbstractTimeseriesSolution}};
-    integrator=Tsit5(),
-    metric=abs
+    reference::SciMLBase.AbstractTimeseriesSolution,
+    sol::SciMLBase.AbstractTimeseriesSolution;
+    integrator=Tsit5()
 )
     results = Dict{Symbol, Any}()
     reference_container = symbolic_container(reference)
-    containers = symbolic_container.(last.(sols))
+    containers = symbolic_container(sol)
 
     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"
+    sol_measured = measured_values(containers)
+    @assert _symbolic_subset(measured_reference, sol_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
-
+    @assert sol.dense "Test solution must be dense in order to use continous-time comparison"
+    @assert timebounds(sol) == timebounds(reference) "Test solution has time range $(timebounds(sol)) which differs from the reference $(timebounds(reference))"
 
     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
+        du .= abs.(reference(t, idxs=measured) .- sol(t, idxs=measured))
     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))
+    func = ODEFunction(compare!; sys = SymbolCache(collect(Symbol.(measured)), [], ref_t_var))
+    prob = ODEProblem(func, zeros(length(measured)), timebounds(reference))
     soln = solve(prob, integrator)
     return soln
 end

src/test/discrete/single_shooting.jl

@@ -5,36 +5,19 @@ function validate(model, data)
     throw("Discrete comparison currently requires JSMO!")
 end
 
-function compare_discrete_to_continous(
-    sol::SciMLBase.AbstractTimeseriesSolution,
-    data::DataFrame)
-    @assert "timestamp" ∈ names(data) "The data must contain a column named `timestamp`"
-    sort!(data, [:timestamp])
-    @assert data[!, :timestamp] == sol.t "The data's discretization points and the solution time discretization must match"
-    
-    results = Dict{Symbol, Any}()
-    reference_container = symbolic_container(sol)
-    measured_reference = measured_values(reference_container)
-    @assert length(measured_reference) > 0 "At least one variable must be marked as measured"
-    measured = measured_reference
-    
-    measured_names = string.(measured)
-    if !(all(name->name ∈ names(data), measured_names))
-        error("Measured data points must exist in both model solution & test data; measured parameters in solution: $(measured_names) vs. in data $(names(data))")
-    end
-
-    data_matrix = collect.(eachrow(data[!, string.(measured)]))  # lame and slow
-    solution_data = sol[measured]
-    metrics = Dict{Symbol, Any}()
-    compute_error_metrics(metrics, solution_data, data_matrix)
-    metrics[:final] = recursive_mean(abs.(solution_data[end] - collect(data[end, string.(measured)])))
-    results[:metrics] = metrics
-    cols = []
-    push!(cols, :timestamp => data[:, :timestamp])
-    append!(cols, make_cols(namespace_symbol.((:simulated,), measured), solution_data))
-    append!(cols, make_cols(namespace_symbol.((:data,), measured), data_matrix))
-    results[:data] = DataFrame(cols)
-    return results
+function compare_discrete(sys, data_a::DataFrame, data_b::DataFrame)
+    @assert "timestamp" ∈ names(data_a) "The dataset A must contain a column named `timestamp`"
+    @assert "timestamp" ∈ names(data_b) "The dataset B must contain a column named `timestamp`"
 
+    measured = measured_values(sys)
+    measured_cols = measured_names(measured)
+    @assert all(c->c ∈ names(data_a), measured_cols) "All measured values must exist in both datasets (missing value in A)"
+    @assert all(c->c ∈ names(data_b), measured_cols) "All measured values must exist in both datasets (missing value in B)"
+    
+    test_results = Dict{Symbol, Any}()
+    delta_sol = data_a[:, measured_cols] .- data_b[:, measured_cols]
+    test_results[:final] = recursive_mean(abs.(collect(delta_sol[end, :])))
+    compute_error_metrics(test_results, collect.(eachrow(data_a[:, measured_cols])), collect.(eachrow(data_b[:, measured_cols])))
+    return test_results
 end
-export validate, compare_discrete_to_continous
+export validate, compare_discrete_to_continous, compare_discrete

src/test/measured.jl

@@ -11,3 +11,12 @@ end
 function measured_values(sys, v=all_variable_symbols(sys))
     filter(x -> ismeasured(x, false), v)
 end
+
+function measured_system_values(sys)
+    reference_container = symbolic_container(sys)
+    return measured_values(reference_container)
+end
+
+function measured_names(measured)
+    return string.(measured)
+end

test/timeseries.jl

@@ -17,25 +17,30 @@ tol = 1e-6
         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_solutions(:reference=>sol1, [:good=>sol2, :bad=>sol3])
-        
-        @test results[:metrics][:good][:l∞] < tol
-        @test results[:metrics][:good][:l2] < tol
-        @test results[:metrics][:good][:final] < tol
-        @test results[:metrics][:bad][:l∞] > tol
-        @test results[:metrics][:bad][:l2] > tol
-        @test results[:metrics][:bad][:final] > tol
+        d1 = discretize_solution(sol1, sol1)
+        d2 = discretize_solution(sol2, sol1)
+        d3 = discretize_solution(sol3, sol1)
+        results_good = compare_discrete(fol, d1, d2)
+        results_bad = compare_discrete(fol, d1, d3)
+        @test results_good[:l∞] < tol
+        @test results_good[:l2] < tol
+        @test results_good[:final] < tol
+        @test results_bad[:l∞] > tol
+        @test results_bad[:l2] > tol
+        @test results_bad[:final] > tol
 
         knots = collect(1:0.1:10)
-        results = compare_solutions(:reference=>sol1, [:good=>sol2, :bad=>sol3]; knots=knots)
-        
-        @test nrow(results[:data]) == length(knots)
-        @test results[:metrics][:good][:l∞] < tol
-        @test results[:metrics][:good][:l2] < tol
-        @test results[:metrics][:good][:final] < tol
-        @test results[:metrics][:bad][:l∞] > tol
-        @test results[:metrics][:bad][:l2] > tol
-        @test results[:metrics][:bad][:final] > tol
+        d1 = discretize_solution(sol1, knots)
+        d2 = discretize_solution(sol2, knots)
+        d3 = discretize_solution(sol3, knots)
+        results_good = compare_discrete(fol, d1, d2)
+        results_bad = compare_discrete(fol, d1, d3)
+        @test results_good[:l∞] < tol
+        @test results_good[:l2] < tol
+        @test results_good[:final] < tol
+        @test results_bad[:l∞] > tol
+        @test results_bad[:l2] > tol
+        @test results_bad[:final] > tol
     end
     @testset "Model-Model Continous Comparison" begin 
         @variables t 
@@ -52,11 +57,12 @@ tol = 1e-6
         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
+        results_good = compare_dense_solutions(sol1, sol2)
+        results_bad = compare_dense_solutions(sol1, sol3)
+        @test results_good[Symbol("x(t)")][end] < tol
+        @test results_good[Symbol("y(t)")][end] < tol
+        @test results_bad[Symbol("x(t)")][end] > tol
+        @test results_bad[Symbol("y(t)")][end] > tol
     end
 end