diff --git a/experiments/audit_density_order.py b/experiments/audit_density_order.py index db1151c..ca49fec 100644 --- a/experiments/audit_density_order.py +++ b/experiments/audit_density_order.py @@ -527,6 +527,89 @@ def audit_closure_rules(table, support): for (order, points), count in sorted(counts.items())]} +def audit_line_search(table, history, initial_residual, strategy): + require(strategy == "baseline", "Trial audit currently covers the recorded baseline policy only") + # Reviewed against LineSearchTrialDisposition in solver/newton.cppm and + # PreparationStageName/line-search dispatch in solver/stellar_context.cppm. + dispositions = {0: "accepted", 1: "inadmissible_state", 2: "non_finite_state", + 3: "non_finite_residual", 4: "insufficient_decrease", 5: "auxiliary_solve_failure"} + stages = {"generated geometry", "gravity", "barotropic closure", "hydrostatic equilibrium", + "displacement residual", "pressure force", "gravity displacement force", "rotational displacement force", + "displacement composition", "mass normalization", "model specification", "unspecified preparation", "unknown preparation"} + grouped, source_counts, disposition_counts, joint_counts = defaultdict(list), defaultdict(int), defaultdict(int), defaultdict(int) + saved = [] + for raw in table: + r = {**numeric(raw, ("source",)), "source": raw["source"]} + iteration, trial, code = int(raw["iteration"]), int(raw["trial"]), int(raw["disposition"]) + require(iteration >= 0 and trial >= 0 and code in dispositions, "Invalid trial index/disposition") + require(math.isfinite(r["alpha"]) and 0 < r["alpha"] <= 1, "Invalid trial step length") + require(not saved or (iteration, trial) > (int(saved[-1]["iteration"]), int(saved[-1]["trial"])), "Trial log is not sequential") + require(trial == len(grouped[iteration]), "Missing/repeated trial within an iteration") + if grouped[iteration]: + require(r["alpha"] < grouped[iteration][-1]["alpha"], "Backtracking trial did not shorten the step") + require(grouped[iteration][-1]["disposition"] != 0, "Trial recorded after acceptance") + source = r["source"] + if code == 0: + require(source == "", "Accepted trial has a rejection source") + elif code == 4: + require(source == "residual metric", "Insufficient decrease source mismatch") + elif code == 1: + require(source in stages or source == "enthalpy polynomial", "Inadmissible-state source mismatch") + elif code == 2: + require(source in stages or source == "candidate state", "Nonfinite-state source mismatch") + elif code == 3: + require(source in stages or source in ("normalized residual", "residual metric"), "Nonfinite-residual source mismatch") + else: + require(False, "Baseline monolithic trial unexpectedly reports an auxiliary solve failure") + if code in (0, 4): + require(math.isfinite(r["F"]) and r["F"] >= 0 and math.isfinite(r["min_jacobian"]) and r["min_jacobian"] > 0, + "Metric-evaluated trial lacks finite residual/valid geometry") + elif source != "residual metric": + require(math.isnan(r["F"]), "Rejected preparation has a fabricated residual metric") + require(math.isfinite(r["min_jacobian"]) or math.isnan(r["min_jacobian"]), "Invalid optional trial geometry") + grouped[iteration].append(r) + saved.append(r) + source_counts[source] += 1 + disposition_counts[dispositions[code]] += 1 + joint_counts[(source, code)] += 1 + history_by_iteration = keyed(history, lambda row: int(row["iteration"])) + require(set(grouped) <= set(history_by_iteration), "Trials lack an after-iteration record") + previous = initial_residual + accepted_alphas = [] + for iteration, row in sorted(history_by_iteration.items()): + trials = grouped[iteration] + require(len(trials) == row["trials"], "Trial/iteration count mismatch") + accepted = [r for r in trials if r["disposition"] == 0] + require(len(accepted) == int(row["accepted"]), "Trial/iteration acceptance mismatch") + for trial in trials: + if trial["disposition"] in (0, 4): + limit = previous * math.sqrt(1 - 1e-4 * trial["alpha"]) + allowance = 128 * EPS * max(previous, trial["F"]) + if trial["disposition"] == 0: + require(trial["F"] <= limit + allowance, "Accepted baseline trial fails recorded merit inequality") + else: + require(trial["F"] >= limit - allowance, "Rejected baseline metric satisfies decrease away from rounding") + if accepted: + require(accepted[0] is trials[-1], "Accepted trial was not the last attempt") + close(accepted[0]["alpha"], row["alpha"], "Accepted alpha binding") + close(accepted[0]["F"], row["F"], "Accepted residual binding") + accepted_alphas.append(row["alpha"]) + else: + require(row["alpha"] == 0, "Failed iteration has a nonzero accepted step") + close(row["F"], previous, "Failed iteration residual restoration") + previous = row["F"] + return {"scope": "Recorded baseline trial stages and metric inequalities; rejected preparation has no evaluated residual metric. Geometry may shorten the backtracking ladder.", + "disposition_enum": {str(code): name for code, name in dispositions.items()}, + "total_trials": len(saved), "counts_by_source": dict(sorted(source_counts.items())), + "counts_by_disposition": dict(sorted(disposition_counts.items())), + "source_disposition_counts": [{"source": source, "disposition": dispositions[code], "disposition_code": code, "count": count} + for (source, code), count in sorted(joint_counts.items())], + "pressure_force_rejections": sum(r["source"] == "pressure force" and r["disposition"] != 0 for r in saved), + "residual_metric_rejections": sum(r["source"] == "residual metric" and r["disposition"] != 0 for r in saved), + "residual_metric_insufficient_decrease_rejections": sum(r["source"] == "residual metric" and r["disposition"] == 4 for r in saved), + "minimum_accepted_alpha": min(accepted_alphas) if accepted_alphas else None, "trials": saved} + + def audit_case(campaign, name): get = lambda filename: campaign.artifact(name, filename) record = campaign.run(name) @@ -621,6 +704,7 @@ def audit_case(campaign, name): close(row["linear_initial_over_rhs"], row["linear_initial_residual"] / row["linear_rhs_norm"], "Initial linear defect ratio") final_f = result["physical"]["18"]["metrics"]["normalized_bordered_residual"] initial_f = result["initial_physical"]["18"]["metrics"]["normalized_bordered_residual"] + result["line_search"] = audit_line_search(rows(get("trials.csv")), history, initial_f, meta["strategy"]) if history: close(history[-1]["F"], final_f, "Final residual/log binding") threshold = max(float(meta["atol"]), float(meta["rtol"]) * initial_f) @@ -655,6 +739,24 @@ def audit_case(campaign, name): class SelfTests(unittest.TestCase): + def test_trial_sources_and_sequential_binding(self): + trial_rows = [dict(iteration="0", trial="0", alpha="1", disposition="1", source="pressure force", F="nan", min_jacobian="nan"), + dict(iteration="0", trial="1", alpha="0.5", disposition="4", source="residual metric", F="1.1", min_jacobian=".8"), + dict(iteration="0", trial="2", alpha="0.25", disposition="0", source="", F=".8", min_jacobian=".9")] + history = [dict(iteration=0., trials=3., accepted=1., alpha=.25, F=.8)] + result = audit_line_search(trial_rows, history, 1., "baseline") + self.assertEqual(result["pressure_force_rejections"], 1) + self.assertEqual(result["residual_metric_rejections"], 1) + self.assertEqual(result["minimum_accepted_alpha"], .25) + for index, key, value in ((0, "F", "1"), (1, "F", ".7"), (1, "trial", "2"), + (2, "alpha", ".5"), (2, "source", "pressure force"), (2, "F", ".7")): + altered = [dict(row) for row in trial_rows] + altered[index][key] = value + with self.subTest(key=key), self.assertRaises(AssertionError): + audit_line_search(altered, history, 1., "baseline") + failed = [dict(iteration=0., trials=0., accepted=0., alpha=0., F=1.)] + self.assertIsNone(audit_line_search([], failed, 1., "baseline")["minimum_accepted_alpha"]) + def test_actual_closure_rule_support_and_tensor_counts(self): rows = [{"element": "2", "quadrature_order": "19", "quadrature_points": "1000"}, {"element": "5", "quadrature_order": "21", "quadrature_points": "1331"}] diff --git a/experiments/summarize_density_order.py b/experiments/summarize_density_order.py index 06c7e94..5f3145e 100644 --- a/experiments/summarize_density_order.py +++ b/experiments/summarize_density_order.py @@ -133,6 +133,12 @@ def build_rows(audits): # unbordered equations still carry the auxiliary border action. for key in ("normalized_unbordered_residual", "central_border", "normalized_border_action"): row["final_" + key] = finite(physical18[key], "final " + key, nonnegative=key != "central_border") + line_search = case["line_search"] + for key in ("pressure_force_rejections", "residual_metric_rejections"): + row[key] = finite(line_search[key], key, True) + row["minimum_accepted_alpha"] = line_search["minimum_accepted_alpha"] + if row["minimum_accepted_alpha"] is not None: + require(0 < finite(row["minimum_accepted_alpha"], "minimum accepted alpha") <= 1, "Invalid accepted step") for q in (14, 18): physical = case["physical"][str(q)]["metrics"] for source, label in (("enthalpy_virial_error", "h_energy_virial"), @@ -320,6 +326,7 @@ class Tests(unittest.TestCase): "final_residual": .01, "threshold": 1e-8, "accepted_steps": 1, "linear_iterations": 7, "iteration_seconds": 10., "iterations": [{"linear_seconds": 6.}]}, "observer_timings": {"observer_total": 200.}, "run": {"seconds": 215., "utc_minus_monotonic_seconds": .01}, + "line_search": {"pressure_force_rejections": 4, "residual_metric_rejections": 1, "minimum_accepted_alpha": 1/32}, "physical": {"14": {"metrics": metrics}, "18": {"metrics": metrics}}} audit = {"artifact_passed": True, "campaigns": {"/recorded/campaign": { "producer_commit": "test", "mesh_sha256": "test", "cases": {"sample": case}}}} @@ -336,6 +343,9 @@ class Tests(unittest.TestCase): self.assertEqual(row["final_normalized_unbordered_residual"], .04) self.assertEqual(row["final_central_border"], -.125) self.assertEqual(row["final_normalized_border_action"], .03) + self.assertEqual(row["pressure_force_rejections"], 4) + self.assertEqual(row["residual_metric_rejections"], 1) + self.assertEqual(row["minimum_accepted_alpha"], 1/32) def main(): @@ -364,6 +374,7 @@ def main(): "surface_minimum_enthalpy", "surface_negative_enthalpy_samples"]) write_csv(args.output / "coupled_outcomes.csv", coupled, ["case", "n", "J", "density_degree", "source_state_dofs", "density_dofs", "converged", "failure", "accepted_steps", "krylov_iterations", "newton_iteration_seconds", "linear_solve_seconds", + "pressure_force_rejections", "residual_metric_rejections", "minimum_accepted_alpha", "initial_residual", "final_residual", "final_over_initial_residual", "h_energy_virial_q18", "h_force_virial_q18", "final_normalized_unbordered_residual", "final_central_border", "final_normalized_border_action", "h_energy_virial_q14", "h_force_virial_q14", "negative_mass_over_abs_mass", "negative_volume_fraction", @@ -380,6 +391,7 @@ def main(): "Uncertainty scales are reconstruction indicators, not rigorous bounds; zero observed negativity is not a global positivity certificate.", "Newton seconds sum reported nonlinear-iteration times, excluding observer and setup costs; total elapsed time is kept separately.", "Final unbordered residual, signed central border and normalized border action come from q18 physical diagnostics; bordered convergence does not imply that these vanish.", + "Pressure-force and residual-metric rejection counts retain the recorded stage distinction; residual-metric counts include any nonfinite metric failures as well as insufficient decrease.", "Ratios to Q2 compare recorded endpoints and do not imply equal convergence or error tolerances were attained.", "Undefined values are null/blank; zero is kept exactly. No logarithmic plot floor is inserted.", "P(h) virial milestones do not establish pointwise density validity; density-based virials are retained separately."],