Skip to content

Harden simulation and reference comparison validation #29

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
AnHeuermann merged 1 commit into from
Mar 18, 2026
Merged

Harden simulation and reference comparison validation #29

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
59 changes: 36 additions & 23 deletions src/compare.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -370,12 +370,17 @@ function compare_with_reference(
isempty(times) && return 0, 0, 0, ""

# Determine which signals to compare: prefer comparisonSignals.txt
sig_file = joinpath(dirname(ref_csv_path), "comparisonSignals.txt")
signals = if isfile(sig_file)
filter(s -> lowercase(s) != "time" && !isempty(s), strip.(readlines(sig_file)))
sig_file = joinpath(dirname(ref_csv_path), "comparisonSignals.txt")
using_sig_file = isfile(sig_file)
signals = if using_sig_file
sigs = filter(s -> lowercase(s) != "time" && !isempty(s), strip.(readlines(sig_file)))
sigs_missing = filter(s -> !haskey(ref_data, s), sigs)
isempty(sigs_missing) || error("Signal(s) listed in comparisonSignals.txt not present in reference CSV: $(join(sigs_missing, ", "))")
sigs
else
filter(k -> lowercase(k) != "time", collect(keys(ref_data)))
end
n_total = length(signals)

# ── Build variable accessor map ──────────────────────────────────────────────
# var_access: normalized name → Int (state index) or MTK symbolic (observed).
Expand All @@ -397,32 +402,40 @@ function compare_with_reference(
@warn "Could not enumerate observed variables: $(sprint(showerror, e))"
end

# Clip reference time to the simulation interval
# Verify the simulation covers the expected reference time interval.
# A large gap means the solver stopped early or started late.
isempty(sol.t) && return n_total, 0, 0, ""
t_start = sol.t[1]
t_end = sol.t[end]
ref_t_start = times[1]
ref_t_end = times[end]
if t_start > ref_t_start || t_end < ref_t_end
@error "Simulation interval [$(t_start), $(t_end)] does not cover " *
"reference interval [$(ref_t_start), $(ref_t_end)]"
return n_total, 0, 0, ""
end

# Clip reference time to the simulation interval
valid_mask = (times .>= t_start) .& (times .<= t_end)
t_ref = times[valid_mask]
isempty(t_ref) && return 0, 0, 0, ""
isempty(t_ref) && return n_total, 0, 0, ""

n_total = 0
n_pass = 0
pass_sigs = String[]
fail_sigs = String[]
skip_sigs = String[]
fail_scales = Dict{String,Float64}()

for sig in signals
haskey(ref_data, sig) || continue # signal absent from ref CSV entirely
signal_name = _normalize_var(sig)
ref_vals = ref_data[sig][valid_mask]

norm = _normalize_var(sig)
if !haskey(var_access, norm)
push!(skip_sigs, sig)
if !haskey(var_access, signal_name)
push!(fail_sigs, sig)
fail_scales[sig] = isempty(ref_vals) ? 0.0 : maximum(abs, ref_vals)
continue
end

accessor = var_access[norm]
ref_vals = ref_data[sig][valid_mask]
n_total += 1
accessor = var_access[signal_name]

# Peak |ref| — used as amplitude floor so relative error stays finite
# near zero crossings.
Expand All @@ -431,10 +444,10 @@ function compare_with_reference(
# Interpolate simulation at reference time points.
sim_vals = [_eval_sim(sol, accessor, t) for t in t_ref]

# If evaluation returned NaN (observed-var access failed), treat as skip.
# If evaluation returned NaN (observed-var access failed), treat as fail.
if any(isnan, sim_vals)
n_total -= 1
push!(skip_sigs, sig)
push!(fail_sigs, sig)
fail_scales[sig] = ref_scale
continue
end

Expand Down Expand Up @@ -468,9 +481,10 @@ function compare_with_reference(
for sig in fail_sigs
ref_vals = ref_data[sig][valid_mask]
r = ref_vals[ti]
s = _eval_sim(sol, var_access[_normalize_var(sig)], t)
acc = get(var_access, _normalize_var(sig), nothing)
s = acc === nothing ? NaN : _eval_sim(sol, acc, t)
ref_scale = get(fail_scales, sig, 0.0)
relerr = abs(s - r) / max(abs(r), ref_scale, settings.abs_tol)
relerr = isnan(s) ? NaN : abs(s - r) / max(abs(r), ref_scale, settings.abs_tol)
push!(row, @sprintf("%.10g", r),
@sprintf("%.10g", s),
@sprintf("%.6g", relerr))
Expand All @@ -481,12 +495,11 @@ function compare_with_reference(
end

# ── Write detail HTML whenever there is anything worth showing ───────────────
if !isempty(fail_sigs) || !isempty(skip_sigs)
if !isempty(fail_sigs)
write_diff_html(model_dir, model;
diff_csv_path = diff_csv,
pass_sigs = pass_sigs,
skip_sigs = skip_sigs)
pass_sigs = pass_sigs)
end

return n_total, n_pass, length(skip_sigs), diff_csv
return n_total, n_pass, 0, diff_csv
end
22 changes: 17 additions & 5 deletions src/pipeline.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,11 +189,23 @@ function main(;
result = test_model(omc, model, results_root, ref_root; csv_max_size_mb)
push!(results, result)

phase = result.sim_success ? "SIM OK" :
result.parse_success ? "SIM FAIL" :
result.export_success ? "PARSE FAIL" : "EXPORT FAIL"
cmp_info = result.cmp_total > 0 ?
" cmp=$(result.cmp_pass)/$(result.cmp_total)" : ""
phase = if result.sim_success && result.cmp_total > 0
result.cmp_pass == result.cmp_total ? "CMP OK" : "CMP FAIL"
elseif result.sim_success
"SIM OK"
elseif result.parse_success
"SIM FAIL"
elseif result.export_success
"PARSE FAIL"
else
"EXPORT FAIL"
end
cmp_info = if result.cmp_total > 0
skip_note = result.cmp_skip > 0 ? " skip=$(result.cmp_skip)" : ""
" cmp=$(result.cmp_pass)/$(result.cmp_total)$skip_note"
else
""
end
@info " → $phase export=$(round(result.export_time;digits=2))s" *
" parse=$(round(result.parse_time;digits=2))s" *
" sim=$(round(result.sim_time;digits=2))s$cmp_info"
Expand Down
28 changes: 23 additions & 5 deletions src/simulate.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,16 @@ function run_simulate(ode_prob, model_dir::String, model::String;
end
sim_time = time() - t0
if sol.retcode == ReturnCode.Success
sim_success = true
sys = sol.prob.f.sys
n_vars = length(ModelingToolkit.unknowns(sys))
n_obs = length(ModelingToolkit.observed(sys))
if isempty(sol.t)
sim_error = "Simulation produced no time points"
elseif n_vars == 0 && n_obs == 0
sim_error = "Simulation produced no output variables (no states or observed)"
else
sim_success = true
end
else
sim_error = "Solver returned: $(sol.retcode)"
end
Expand All @@ -50,21 +59,30 @@ function run_simulate(ode_prob, model_dir::String, model::String;
isempty(sim_error) || println(log_file, "\n--- Error ---\n$sim_error")
close(log_file)

# Write simulation results CSV (time + all state variables)
# Write simulation results CSV (time + state variables + observed variables)
if sim_success && sol !== nothing
short_name = split(model, ".")[end]
sim_csv = joinpath(model_dir, "$(short_name)_sim.csv")
try
sys = sol.prob.f.sys
vars = ModelingToolkit.unknowns(sys)
col_names = [_clean_var_name(string(v)) for v in vars]
sys = sol.prob.f.sys
vars = ModelingToolkit.unknowns(sys)
obs_eqs = ModelingToolkit.observed(sys)
obs_syms = [eq.lhs for eq in obs_eqs]
col_names = vcat(
[_clean_var_name(string(v)) for v in vars],
[_clean_var_name(string(s)) for s in obs_syms],
)
open(sim_csv, "w") do f
println(f, join(["time"; col_names], ","))
for (ti, t) in enumerate(sol.t)
row = [@sprintf("%.10g", t)]
for vi in eachindex(vars)
push!(row, @sprintf("%.10g", sol[vi, ti]))
end
for sym in obs_syms
val = try Float64(sol(t; idxs = sym)) catch; NaN end
push!(row, @sprintf("%.10g", val))
end
println(f, join(row, ","))
end
end
Expand Down
Loading