Harden simulation and reference comparison validation (#29)

- Fail run_simulate when sol.t is empty or the system has no states
  or observed variables, preventing empty CSVs and downstream panics
- Include observed variables in the simulation CSV output so models
  without state variables (e.g. BusUsage) still produce usable results
- Fail compare_with_reference (instead of skip) when reference signals
  are absent from the simulation; write NaN sim/relerr columns for
  signals with no simulation counterpart
- Validate upfront that all comparisonSignals.txt entries exist in the
  reference CSV, and error if the simulation time interval does not
  cover the reference interval
- Extend the per-model status line to show a CMP OK/FAIL phase and
  signal counts including skipped signals

Co-authored-by: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: AnHeuermann

src/compare.jl

@@ -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).
@@ -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.
@@ -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
 
@@ -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))
@@ -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

src/pipeline.jl

@@ -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"

src/simulate.jl

@@ -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
@@ -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