diff --git a/experiments/DAMPED_COUPLED_TRIAL_2026-09-18.md b/experiments/DAMPED_COUPLED_TRIAL_2026-09-18.md new file mode 100644 index 0000000..e66d84f --- /dev/null +++ b/experiments/DAMPED_COUPLED_TRIAL_2026-09-18.md @@ -0,0 +1,67 @@ +# Damped fully coupled trial correction + +This investigation branches from `46249d0` onto +`research/damped-coupled-trial-20260918`. It tests two separate mechanisms: +bounded damping of one fully coupled correction, and reserving that correction +solve for a trial whose merit is reasonably close to its affine prediction. +The previous full-step policy remains available and the default solver is unchanged. + +## Algorithm and architecture + +At accepted state x, capture F0=F(x), the actual Newton direction p and J0 p. +For outer step alpha and prepared candidate y (raw or gravity-corrected), set +L=F0+alpha J0 p and solve J(y)c=L-F(y). The unscaled c must satisfy both the +existing independent full-J defect check and the existing norm cap. + +The optional search prepares y+beta c with beta=1,1/2,1/4,1/8, always starting +from y. One auxiliary solve supplies all four candidates. The first improving +candidate passing ordinary outer acceptance ends the search. Otherwise the +lowest-merit admissible improvement is retained; an earlier selected candidate +is freshly prepared, and complete failure restores y with its current inverse. +No accepted-state buffer is written by this transaction. + +For the optional attempt selector, D=M(F0)-M(L) must be finite and positive. +The solve is eligible when (M(F(y))-M(F0))/D <= k, with default k=1. A deferral +spends no auxiliary budget. This is a bounded-work heuristic, not a convergence +guarantee. It requires a deterministic metric of the supplied residual with +fixed parameters; the default normalized residual metric satisfies this contract. + +Algorithms are selected by types: `FullCoupledStep`/`BacktrackedCoupledStep` +and `EveryRejectedCoupledAttempt`/`NearAcceptanceCoupledAttempt`. Runtime +iteration counts, reduction factors and tolerances are validated collectively. +`BacktrackedCoupledCurvatureCorrection` composes damping and near-acceptance +selection in the existing typed `NewtonNumerics` recipe. The original +`CoupledCurvatureCorrection` remains the full-step/every-rejected control. +Neither changes the physics-facing model specification or central-density policy. + +An auxiliary report describes c, not beta c. For a damped candidate the exact +remainder identity is + + F(y+beta c)-L = (1-beta)(F(y)-L) + beta (J(y)c-(L-F(y))) + E_beta, + E_beta = F(y+beta c)-F(y)-beta J(y)c. + +Scaling c therefore does not imply solving the original auxiliary RHS to its +original tolerance. The actual freshly prepared residual decides acceptance. + +## Reproduction and evidence + +The versioned `damped_coupled_screen_2026-09-18.json` compares none/full/damped/near +for the late n=1.5,J=0 snapshot on raw and gravity-corrected paths, plus a seed +control. No central constraint is imposed. Each sequence captures once, shares +one auxiliary-solve budget, tests every halving from alpha=1 through 1/4096 and +stops at its first accepted candidate. The prior independent-alpha screen did +not represent this shared budget and omitted some intermediate step lengths. + +The campaign runner builds a clean committed revision, copies an immutable +executable and hashes mesh, inputs and outputs. The experiment calls the same +production transaction and production gravity correction as Newton. It records +every candidate evaluation, fresh selected-state replay and unchanged accepted +buffers/normalization. Only the accepted trial is replayed for quadrature-14 and +quadrature-18 physical observations, outside the timed selection phase. + +Operation timings exclude artifact output and physical replay; search wall time +also includes artifact output. Callback vector copies remain in the measured +transaction. The no-correction control captures the same tangent for a comparable +screen, so capture cost is reported separately and is not ordinary Newton work +for that control. All builds and numerical cases run serially; artifact audits +run after timed FEM work. diff --git a/experiments/check_damped_coupled_trial.py b/experiments/check_damped_coupled_trial.py new file mode 100644 index 0000000..b5eb903 --- /dev/null +++ b/experiments/check_damped_coupled_trial.py @@ -0,0 +1,407 @@ +#!/usr/bin/env python3 +"""Audit a sequential frozen line search, separately from numerical success. + +Only saved-vector arithmetic is independent here: no FEM action or physical +quadrature is repeated. A shared attempt budget and first-accept stopping rule +are checked explicitly; these frozen searches are not Newton trajectories. +""" +import argparse +import json +import math +from pathlib import Path + +from check_nonlinear_step import EPS, close, digest, dot, layout_for, minus, norm, rows, sanitize, vector +from check_coupled_trial import inspect_physical, near_vector + +P = 'damped_coupled_' +STATUS = ('not_attempted', 'original_accepted', 'budget_exhausted', 'kept', 'linear_failure', + 'norm_limit', 'non_finite', 'inadmissible', 'not_improved', 'prediction_deferred') +BASE = {'metadata.txt', 'coupled_trial_physical.csv'} | { + P + name + '.csv' for name in ('summary', 'trials', 'candidates', 'evaluations', 'layout')} | { + P + name + '.txt' for name in ('base', 'base_physical', 'base_residual', 'direction', 'base_action', + 'state_factors', 'residual_factors', 'restored_residual', 'restored_action')} + + +def same_number(a, b, scale=0.): + return (math.isnan(a) and math.isnan(b)) or (not math.isfinite(a) and a == b) or close(a, b, scale) + + +def merit(values): + return .5 * math.fsum(value * value for value in values) + + +def eligibility(base_merit, prediction_merit, original_merit, maximum_excess): + decrease = base_merit - prediction_merit + ratio = (original_merit - base_merit) / decrease if math.isfinite(decrease) and decrease > 0. else math.nan + return math.isfinite(ratio) and ratio <= maximum_excess, decrease, ratio + + +def correction_mechanism(prediction, before, corrected, action, rhs, beta): + remainder = minus(before, prediction) + defect = minus(action, rhs) + own = [z - y - beta * a for z, y, a in zip(corrected, before, action)] + after = minus(corrected, prediction) + reconstructed = [(1. - beta) * r + beta * d + e for r, d, e in zip(remainder, defect, own)] + return {'remainder_before_norm': norm(remainder), 'remainder_after_norm': norm(after), + 'scaled_linear_defect_norm': abs(beta) * norm(defect), 'own_taylor_remainder_norm': norm(own), + 'decomposition_error_norm': norm(minus(after, reconstructed))} + + +def inspect_case(folder): + try: + return inspect_case_artifacts(Path(folder)) + except (OSError, ValueError, KeyError, TypeError, IndexError, OverflowError, ZeroDivisionError) as error: + return {'passed': False, 'observations_valid': False, 'required_artifacts': sorted(BASE), + 'checks': [{'check': 'complete_readable_artifacts', 'passed': False, 'error': str(error)}], + 'trials': [], 'physical_observations': []} + + +def inspect_case_artifacts(folder): + checks, required = [], set(BASE) + def check(name, condition): + checks.append({'check': name, 'passed': bool(condition)}) + def flag(row, key, label): + value = float(row[key]); check(label + key + ':boolean', value in (0., 1.)); return value == 1. + def get(name, available=True, prefix=P): + path = folder / (prefix + name + '.txt') + if available: + required.add(path.name); return vector(path) + check(name + ':unavailable_file_absent', not path.exists()); return None + summary_rows = rows(folder / (P + 'summary.csv')) + summary = {row['metric']: row['value'] if row['metric'] == 'mode' else float(row['value']) for row in summary_rows} + check('unique_summary', len(summary) == len(summary_rows)) + metadata = dict(line.split('=', 1) for line in (folder / 'metadata.txt').read_text().splitlines() if '=' in line) + mode = summary['mode']; check('mode', mode in ('none', 'full', 'damped', 'near')) + check('metadata_scope', metadata['damped_coupled_probe'] == '1' and metadata['coupled_trial_alphas_are_independent'] == '0' and + metadata['coupled_trial_admissibility'] == 'sampled_enthalpy_without_repair' and metadata['coupled_search'] == mode) + for key, source in {'outer_linear_tolerance': 'linear_tol', 'relative_tolerance': 'coupled_tolerance', + 'maximum_iterations': 'coupled_iterations', 'restart': 'coupled_restart', + 'maximum_correction_ratio': 'coupled_maximum_ratio', 'attempt_budget': 'coupled_attempts_per_iteration', + 'beta_reduction': 'coupled_beta_reduction', 'maximum_excess_ratio': 'coupled_excess_ratio'}.items(): + check('metadata:' + key, summary[key] == float(metadata[source])) + for key, expected in {'schema_version': 1, 'attempt_budget': 1, 'shared_iteration_budget': 1, + 'accepted_commits': 0, 'probe_exception': 0, 'accepted_buffers_bitwise_identical': 1, + 'direction_buffers_bitwise_identical': 1, 'normalization_factors_bitwise_identical': 1, + 'fresh_operator_and_inverse_current': 1}.items(): + check(key, summary[key] == expected) + for name in ('original_accepted', 'budget_exhausted', 'kept', 'prediction_deferred'): + check('status_mapping:' + name, summary['status_' + name] == STATUS.index(name)) + for value, name in enumerate(('evaluated', 'inadmissible', 'non_finite')): + check('candidate_status:' + name, summary['candidate_status_' + name] == value) + gravity_path = flag(summary, 'gravity_path', 'summary:') + check('gravity_policy', gravity_path == metadata['strategy'].startswith('curvature-')) + maximum_trials = int(summary['maximum_trials']) + maximum_steps = int(summary['maximum_beta_steps']) + reduction = summary['beta_reduction'] + check('outer_grid_control', 1 <= maximum_trials <= 30 and maximum_trials == summary['maximum_trials']) + check('beta_grid_control', maximum_steps >= 1 and maximum_steps == summary['maximum_beta_steps'] and 0 < reduction < 1) + check('full_step_policy', mode not in ('none', 'full') or maximum_steps == 1) + check('beta_steps_metadata', mode in ('none', 'full') or maximum_steps == int(metadata['coupled_beta_steps'])) + check('finite_controls', all(math.isfinite(summary[key]) and summary[key] > 0. for key in + ('relative_tolerance', 'maximum_correction_ratio', 'maximum_excess_ratio', 'outer_linear_tolerance'))) + check('nonnegative_absolute_tolerance', math.isfinite(summary['absolute_tolerance']) and summary['absolute_tolerance'] >= 0.) + check('sufficient_decrease', 0 < summary['sufficient_decrease'] < 1) + x, physical, f, p, j, sf, rf = [get(name) for name in + ('base', 'base_physical', 'base_residual', 'direction', 'base_action', 'state_factors', 'residual_factors')] + size = len(x) + check('dimensions', all(len(value) == size for value in (physical, f, p, j, sf, rf))) + check('positive_normalization', all(value > 0. for value in sf + rf)) + check('base_normalization', near_vector(x, [value * factor for value, factor in zip(physical, sf)], norm(x))) + layout = rows(folder / (P + 'layout.csv')) + value_layout = layout_for(layout, 'value', size); residual_layout = layout_for(layout, 'residual', size) + gravity_values = {i for name, (offset, count) in value_layout.items() if name in ('gravity_gradient', 'gravity_potential') + for i in range(offset, offset + count)} + gravity_rows = [i for name, (offset, count) in residual_layout.items() if name in ('gravity_gradient_relation', 'poisson_balance') + for i in range(offset, offset + count)] + fn, pn, jn, m0 = norm(f), norm(p), norm(j), merit(f) + check('finite_nonzero_base', math.isfinite(fn) and fn > 0. and math.isfinite(pn) and pn > 0.) + defect = [a + b for a, b in zip(f, j)] + check('outer_true_gate', norm(defect) <= summary['outer_linear_tolerance'] * fn and summary['outer_verified_gate'] == 1) + for key, value in {'base_norm': fn, 'base_merit': m0, 'outer_relative_defect': norm(defect) / fn}.items(): + check(key, close(summary[key], value, value)) + for name, expected in (('residual', f), ('action', j)): + restored = get('restored_' + name) + error = norm(minus(restored, expected)) + check('restored_' + name, near_vector(restored, expected, norm(expected)) and + close(summary['restored_' + name + '_error'], error, norm(expected))) + + trial_rows = rows(folder / (P + 'trials.csv')) + candidate_rows = rows(folder / (P + 'candidates.csv')) + evaluation_rows = rows(folder / (P + 'evaluations.csv')) + check('sequential_trials', [int(row['trial']) for row in trial_rows] == list(range(len(trial_rows))) and + 1 <= len(trial_rows) <= maximum_trials) + candidate_map = {(int(row['trial']), int(row['candidate'])): row for row in candidate_rows} + evaluation_map = {(int(row['trial']), int(row['evaluation'])): row for row in evaluation_rows} + check('unique_candidates', len(candidate_map) == len(candidate_rows)) + check('unique_evaluations', len(evaluation_map) == len(evaluation_rows)) + seen_candidates, seen_evaluations, outcomes = set(), set(), [] + attempts, deferred, iterations, operation_seconds = 0, 0, 0, 0. + accepted_trial, physical_trials = -1, [] + for row in trial_rows: + i, alpha = int(row['trial']), float(row['alpha']) + label, suffix = str(i) + ':', '_a' + str(i) + check(label + 'alpha', alpha == math.ldexp(1., -i)) + flags = {key: flag(row, key, label) for key in ('raw_prepared', 'gravity_gate', 'usable', 'before_accepted', + 'attempted', 'aux_gate', 'aux_reported_gate', 'aux_verification_applied', 'rhs_saved', 'correction_saved', + 'verified_vectors_saved', 'selected_prepared', 'selected_accepted', 'selected_original')} + status = int(row['status']); check(label + 'status', 0 <= status < len(STATUS)) + name = STATUS[status] + raw = get('raw' + suffix); prediction = get('prediction' + suffix) + raw_f = get('raw_residual' + suffix, flags['raw_prepared']) + check(label + 'raw_state', near_vector(raw, [a + alpha * b for a, b in zip(x, p)], norm(x) + alpha * pn)) + check(label + 'prediction', near_vector(prediction, [a + alpha * b for a, b in zip(f, j)], fn + alpha * jn)) + check(label + 'preparation_status', (int(row['raw_reason']) == -1 and int(row['raw_stage']) == -1) + if flags['raw_prepared'] else (int(row['raw_reason']) >= 0 and int(row['raw_stage']) >= 0)) + check(label + 'usable_requires_preparation', not flags['usable'] or (flags['raw_prepared'] and flags['gravity_gate'])) + if not gravity_path: + check(label + 'raw_path', flags['gravity_gate'] and flags['usable'] == flags['raw_prepared']) + before = get('before' + suffix, flags['usable']) + before_f = get('before_residual' + suffix, flags['usable']) + gravity_correction = get('gravity_correction' + suffix, gravity_path and flags['gravity_gate']) + if flags['usable'] and not gravity_path: + check(label + 'raw_baseline', before == raw and before_f == raw_f) + if gravity_correction is not None: + check(label + 'gravity_support', len(gravity_correction) == size and bool(gravity_values) and bool(gravity_rows) and + all(value == 0. for index, value in enumerate(gravity_correction) if index not in gravity_values)) + check(label + 'gravity_displacement', near_vector(before, [a + b for a, b in zip(raw, gravity_correction)], norm(raw) + norm(gravity_correction))) + target = norm([before_f[index] - prediction[index] for index in gravity_rows]) + threshold = float(row['gravity_target_threshold']); disagreement = float(row['gravity_affine_disagreement']) + check(label + 'gravity_target', close(float(row['gravity_target_defect']), target, fn + alpha * jn)) + check(label + 'gravity_gate', math.isfinite(threshold) and threshold >= 0. and target <= threshold + 512 * EPS * fn and + math.isfinite(disagreement) and 0. <= disagreement <= threshold) + selected = get('selected' + suffix, flags['selected_prepared']) + selected_f = get('selected_residual' + suffix, flags['selected_prepared']) + check(label + 'selected_prepared', flags['selected_prepared'] == flags['usable']) + cutoff = (1. - summary['sufficient_decrease'] * alpha) * m0 + for stage, value in (('before', before_f), ('selected', selected_f)): + if value is not None: + check(label + stage + '_norm', close(float(row[stage + '_norm']), norm(value), norm(value))) + check(label + stage + '_merit', close(float(row[stage + '_merit']), merit(value), merit(value))) + else: + check(label + stage + '_unavailable', math.isnan(float(row[stage + '_norm'])) and math.isnan(float(row[stage + '_merit']))) + before_accepted = flags['usable'] and merit(before_f) <= cutoff + check(label + 'before_armijo', flags['before_accepted'] == before_accepted) + eligible, predicted_decrease, excess = eligibility(m0, merit(prediction), merit(before_f), summary['maximum_excess_ratio']) \ + if flags['usable'] else (False, math.nan, math.nan) + for key, expected in (('predicted_decrease', predicted_decrease), ('excess_ratio', excess)): + if math.isfinite(float(row[key])): + check(label + key, close(float(row[key]), expected, abs(expected))) + should_attempt = flags['usable'] and mode != 'none' and not before_accepted and attempts < summary['attempt_budget'] and \ + (mode != 'near' or eligible) + check(label + 'attempt_decision', flags['attempted'] == should_attempt) + if mode == 'near' and flags['usable'] and not before_accepted and attempts < summary['attempt_budget']: + check(label + 'gate_recorded', same_number(float(row['predicted_decrease']), predicted_decrease, m0) and + same_number(float(row['excess_ratio']), excess, abs(excess))) + if not should_attempt: + expected = 'not_attempted' if not flags['usable'] or mode == 'none' else ( + 'original_accepted' if before_accepted else ('budget_exhausted' if attempts >= summary['attempt_budget'] else 'prediction_deferred')) + check(label + 'skip_status', name == expected) + deferred += expected == 'prediction_deferred' + if expected == 'prediction_deferred': + check(label + 'deferred_gate_recorded', same_number(float(row['predicted_decrease']), predicted_decrease, m0) and + same_number(float(row['excess_ratio']), excess, abs(excess))) + attempts += should_attempt + check(label + 'shared_budget', int(row['attempts_total']) == attempts <= summary['attempt_budget']) + rhs = get('rhs' + suffix, flags['rhs_saved']); correction = get('correction' + suffix, flags['correction_saved']) + action = get('action' + suffix, flags['verified_vectors_saved']); saved_defect = get('defect' + suffix, flags['verified_vectors_saved']) + check(label + 'save_flags', should_attempt or not any(flags[key] for key in ('rhs_saved', 'correction_saved', 'verified_vectors_saved'))) + check(label + 'finite_verification_saved', not (flags['aux_verification_applied'] and math.isfinite(float(row['aux_defect']))) or flags['verified_vectors_saved']) + core = int(row['aux_core_status']) + if should_attempt: + iterations += int(row['aux_iterations']) + check(label + 'iteration_limit', 0 <= int(row['aux_iterations']) <= summary['maximum_iterations']) + check(label + 'reported_status', flags['aux_reported_gate'] == (int(row['aux_status']) == 0)) + cap = summary['maximum_correction_ratio'] * alpha * pn + check(label + 'unscaled_cap', same_number(float(row['correction_cap']), cap, cap)) + if rhs is not None: + check(label + 'rhs', near_vector(rhs, minus(prediction, before_f), norm(prediction) + norm(before_f))) + if math.isfinite(float(row['aux_threshold'])): + check(label + 'rhs_norm', close(float(row['aux_rhs_norm']), norm(rhs), norm(rhs))) + check(label + 'threshold', close(float(row['aux_threshold']), max(summary['absolute_tolerance'], + summary['relative_tolerance'] * norm(rhs)), float(row['aux_threshold']))) + if flags['verified_vectors_saved']: + check(label + 'verified_availability', flags['aux_verification_applied'] and rhs is not None and correction is not None) + delta = minus(action, rhs); dn, cn = norm(delta), norm(correction) + check(label + 'saved_defect', near_vector(saved_defect, delta, norm(action) + norm(rhs))) + check(label + 'true_defect', same_number(float(row['aux_defect']), dn, norm(rhs))) + check(label + 'true_gate', flags['aux_gate'] == (math.isfinite(dn) and dn <= float(row['aux_threshold']))) + check(label + 'correction_norm', same_number(float(row['correction_norm']), cn, cn)) + expected_core = 3 if not math.isfinite(cn) else (1 if not flags['aux_reported_gate'] or not flags['aux_gate'] else (2 if cn > cap else 0)) + check(label + 'core_decision', core == expected_core) + if core == 0: + check(label + 'admissible_requires_gates', flags['aux_reported_gate'] and flags['aux_gate'] and + flags['verified_vectors_saved'] and float(row['correction_norm']) <= cap) + else: + check(label + 'failed_core_status', name == {1: 'linear_failure', 2: 'norm_limit', 3: 'non_finite'}.get(core)) + check(label + 'failed_core_no_preparation', int(row['preparations']) == 0 and int(row['restorations']) == 0) + else: + check(label + 'no_aux_certificate', core == -1 and int(row['aux_status']) == -1 and not flags['aux_gate'] and + not flags['aux_reported_gate'] and not flags['aux_verification_applied'] and int(row['aux_iterations']) == 0) + check(label + 'no_aux_preparations', int(row['preparations']) == 0 and int(row['restorations']) == 0) + + count = int(row['candidate_count']) + check(label + 'candidate_budget', 0 <= count <= maximum_steps and (should_attempt and core == 0 or count == 0)) + crows = [candidate_map[(i, k)] for k in range(count)] + seen_candidates.update((i, k) for k in range(count)) + erows = [value for key, value in evaluation_map.items() if key[0] == i] + check(label + 'evaluation_order', [int(value['evaluation']) for value in erows] == list(range(len(erows)))) + observed, mechanisms = [], [] + for erow in erows: + e = int(erow['evaluation']); seen_evaluations.add((i, e)); tag = suffix + '_e' + str(e) + ex, ef = get('evaluated_state' + tag), get('evaluated_residual' + tag) + beta = float(erow['beta']) + check(label + f'e{e}:requires_correction', should_attempt and core == 0 and correction is not None) + scale = norm(before) + norm(correction) + check(label + f'e{e}:displacement', near_vector(ex, [b + beta * c for b, c in zip(before, correction)], scale)) + allowed_beta = [0.] + [reduction ** k for k in range(count)] + check(label + f'e{e}:known_beta_or_restore', any(near_vector(ex, + [b + value * c for b, c in zip(before, correction)], scale) for value in allowed_beta)) + check(label + f'e{e}:norm', same_number(float(erow['residual_norm']), norm(ef), norm(ef))) + check(label + f'e{e}:merit', same_number(float(erow['merit']), merit(ef), merit(ef))) + check(label + f'e{e}:armijo', flag(erow, 'accepted', label) == (merit(ef) <= cutoff)) + observed.append((beta, ex, ef)) + mechanisms.append({'evaluation': e, 'beta': beta, **correction_mechanism(prediction, before_f, ef, action, rhs, beta)}) + evaluated = [] + first_accepted = None + for k, crow in enumerate(crows): + beta = reduction ** k + check(label + f'c{k}:beta', float(crow['beta']) == beta) + cstatus = int(crow['status']); check(label + f'c{k}:status', cstatus in (0, 1, 2)) + if cstatus == 0: + matches = [(ex, ef) for observed_beta, ex, ef in observed if near_vector(ex, + [b + beta * c for b, c in zip(before, correction)], norm(before) + norm(correction))] + check(label + f'c{k}:fresh_evaluation', bool(matches)) + ex, ef = matches[0] + cm = merit(ef); evaluated.append((beta, ex, ef, cm)) + check(label + f'c{k}:merit', close(float(crow['merit']), cm, cm)) + accepted = cm <= cutoff + check(label + f'c{k}:armijo', flag(crow, 'accepted', label) == accepted) + if accepted and first_accepted is None: first_accepted = k + else: + check(label + f'c{k}:rejected', not flag(crow, 'accepted', label)) + check(label + f'c{k}:time', math.isfinite(float(crow['preparation_seconds'])) and float(crow['preparation_seconds']) >= 0.) + if should_attempt and core == 0: + check(label + 'beta_stop_rule', count == (first_accepted + 1 if first_accepted is not None else maximum_steps)) + selected_beta = float(row['selected_beta']) + if flags['usable']: + better_candidates = [item for item in evaluated if item[3] < merit(before_f)] + best = min(better_candidates, key=lambda item: item[3]) if better_candidates else None + if name == 'kept': + fresh = [ef for beta, ex, ef in observed if ex == selected] + check(label + 'kept_best', best is not None and selected_beta == best[0] and selected == best[1]) + check(label + 'selected_fresh_residual', bool(fresh) and near_vector(selected_f, fresh[-1], norm(fresh[-1])) and merit(selected_f) < merit(before_f)) + else: + check(label + 'fallback_original', math.isnan(selected_beta) and selected == before and near_vector(selected_f, before_f, norm(before_f))) + check(label + 'discarded_best_requires_failed_replay', best is None or + (name in ('inadmissible', 'non_finite', 'not_improved') and int(row['restorations']) > 0)) + check(label + 'selected_original', flags['selected_original'] == (selected == before)) + else: + check(label + 'unavailable_selection', math.isnan(selected_beta) and not flags['selected_original']) + accepted_now = flags['usable'] and merit(selected_f) <= cutoff + check(label + 'selected_armijo', flags['selected_accepted'] == accepted_now) + if accepted_now: + check(label + 'first_accept_stops', i == len(trial_rows) - 1 and accepted_trial == -1) + accepted_trial = i + pv = get('selected_physical' + suffix, prefix='coupled_trial_') + check(label + 'physical_normalization', near_vector(selected, [v * a for v, a in zip(pv, sf)], norm(selected))) + physical_trials.append({'sample': str(i), 'selected_prepared': '1', 'selected_accepted': '1', 'selected_norm': row['selected_norm']}) + for key in ('raw_seconds', 'gravity_seconds', 'correction_seconds', 'solve_seconds', 'refresh_seconds', 'verification_seconds', 'preparation_seconds'): + check(label + key, math.isfinite(float(row[key])) and float(row[key]) >= 0.) + operation_seconds += sum(float(row[key]) for key in ('raw_seconds', 'gravity_seconds', 'correction_seconds')) + check(label + 'correction_time_scope', float(row['correction_seconds']) + 1.e-6 >= + sum(float(row[key]) for key in ('solve_seconds', 'refresh_seconds', 'verification_seconds', 'preparation_seconds'))) + outcomes.append({'trial': i, 'alpha': alpha, 'status': name, 'attempted': should_attempt, 'attempts_total': attempts, + 'before_accepted': before_accepted, 'selected_accepted': accepted_now, 'selected_beta': selected_beta, + 'before_norm': float(row['before_norm']), 'selected_norm': float(row['selected_norm']), + 'predicted_decrease': predicted_decrease, 'excess_ratio': excess, 'candidate_mechanisms': mechanisms}) + check('candidate_coverage', seen_candidates == set(candidate_map)) + check('evaluation_coverage', seen_evaluations == set(evaluation_map)) + check('accepted_trial', summary['accepted_trial'] == accepted_trial) + if accepted_trial < 0: + check('exhausted_trial_grid', len(trial_rows) == maximum_trials) + else: + last = trial_rows[-1] + for key, value in {'accepted_alpha': float(last['alpha']), 'accepted_norm': float(last['selected_norm']), + 'accepted_beta': float(last['selected_beta'])}.items(): check(key, same_number(summary[key], value, value)) + for key, value in {'attempts': attempts, 'aux_iterations': iterations, 'deferred_attempts': deferred}.items(): + check(key, summary[key] == value) + check('operation_time_total', close(summary['search_operation_seconds'], operation_seconds, operation_seconds)) + check('wall_time_scope', math.isfinite(summary['search_wall_seconds']) and summary['search_wall_seconds'] + 1.e-6 >= operation_seconds) + physical_observations = inspect_physical(folder, physical_trials, metadata, check) + check('no_unexplained_vectors', {path.name for pattern in (P + '*.txt', 'coupled_trial_*.txt') for path in folder.glob(pattern)} == + {name for name in required if name.endswith('.txt') and name != 'metadata.txt'}) + passed = all(item['passed'] for item in checks) + return sanitize({'passed': passed, 'observations_valid': passed, 'required_artifacts': sorted(required), 'checks': checks, + 'summary': summary, 'metadata': metadata, 'trials': outcomes if passed else [], + 'physical_observations': physical_observations if passed else [], + 'scope': 'Sequential frozen search with one shared auxiliary-solve budget. Saved-vector arithmetic, not independent FEM or Newton convergence.'}) + + +def inspect_campaign(campaign): + campaign = Path(campaign).resolve(); checks, results = [], {} + def check(name, value): checks.append({'check': name, 'passed': bool(value)}) + try: + provenance = json.loads((campaign / 'provenance.json').read_text()) + cases = provenance['cases'] + check('executable_hash', digest(campaign / 'plateau_diagnostic') == provenance['binary_sha256']) + check('mesh_hash', digest(Path(provenance['mesh'])) == provenance['mesh_sha256']) + check('case_manifest', json.loads((campaign / 'cases.json').read_text()) == cases) + check('unique_cases', bool(cases) and len(cases) == len({case['id'] for case in cases})) + for case in cases: + name, arguments = case['id'], case['arguments']; folder = campaign / name + record = json.loads((campaign / (name + '_run.json')).read_text()) + check(name + ':completed', record.get('exit_code') == 0 and bool(record.get('finished_utc')) and not record.get('timed_out', False)) + check(name + ':command', record.get('id') == name and record.get('command') == + [str(campaign / 'plateau_diagnostic'), *arguments, '--mesh', provenance['mesh'], '--output', str(folder)]) + expected_inputs = {} + for option, saved in (('--probe-state', 'base_physical'), ('--probe-direction', 'direction')): + path = Path(arguments[arguments.index(option) + 1]) + path = path if path.is_absolute() else Path(provenance['source']) / path + expected_inputs[str(path.resolve())] = saved + check(name + ':required_input_hashes', expected_inputs.keys() <= record.get('input_state_hashes', {}).keys()) + for path, value in record.get('input_state_hashes', {}).items(): check(name + ':input:' + path, digest(Path(path)) == value) + for filename, value in record.get('output_hashes', {}).items(): check(name + ':output:' + filename, digest(folder / filename) == value) + for path, saved in expected_inputs.items(): + source, target = vector(Path(path)), vector(folder / (P + saved + '.txt')) + check(name + ':input_binding:' + saved, source == target if saved == 'direction' else near_vector(source, target, norm(source))) + result = inspect_case(folder); results[name] = result + check(name + ':required_output_hashes', set(result['required_artifacts']) <= record.get('output_hashes', {}).keys()) + check(name + ':probe_scope', '--damped-coupled-probe' in arguments) + metadata = result.get('metadata', {}) + for option, key in {'--index': 'n', '--angular-momentum': 'J', '--linear-tol': 'linear_tol', + '--coupled-tol': 'coupled_tolerance', '--coupled-iterations': 'coupled_iterations', + '--coupled-restart': 'coupled_restart', '--coupled-max-ratio': 'coupled_maximum_ratio', + '--coupled-beta-steps': 'coupled_beta_steps', '--coupled-beta-reduction': 'coupled_beta_reduction', + '--coupled-excess-ratio': 'coupled_excess_ratio', '--coupled-attempts': 'coupled_attempts_per_iteration', + '--gravity-tol': 'gravity_tolerance', '--gravity-iterations': 'gravity_iterations', '--gravity-restart': 'gravity_restart'}.items(): + if option in arguments: + check(name + ':requested:' + key, key in metadata and float(metadata[key]) == float(arguments[arguments.index(option) + 1])) + for option, key in (('--coupled-search', 'coupled_search'), ('--strategy', 'strategy')): + if option in arguments: check(name + ':requested:' + key, metadata.get(key) == arguments[arguments.index(option) + 1]) + if '--no-central' in arguments: check(name + ':central_constraint', metadata.get('central_phase') == '0') + passed = all(item['passed'] for item in checks) and all(result['passed'] for result in results.values()) + if not all(item['passed'] for item in checks): + for result in results.values(): result.update(passed=False, observations_valid=False, trials=[], physical_observations=[]) + return {'passed': passed, 'source_commit': provenance['commit'], 'checker_sha256': digest(Path(__file__)), + 'checker_helpers_sha256': {name: digest(Path(__file__).with_name(name)) for name in + ('check_coupled_trial.py', 'check_nonlinear_step.py')}, + 'provenance_sha256': digest(campaign / 'provenance.json'), 'provenance_checks': checks, 'cases': results} + except (OSError, ValueError, KeyError, TypeError, IndexError) as error: + for result in results.values(): result.update(passed=False, observations_valid=False, trials=[], physical_observations=[]) + return {'passed': False, 'provenance_checks': checks + [{'check': 'complete_readable_campaign', 'passed': False, 'error': str(error)}], 'cases': results} + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('campaign', type=Path); parser.add_argument('--output', required=True, type=Path) + args = parser.parse_args() + if args.output.resolve() == args.campaign.resolve() or args.campaign.resolve() in args.output.resolve().parents: + parser.error('Write the audit outside the immutable campaign.') + result = inspect_campaign(args.campaign) + args.output.write_text(json.dumps(result, indent=2, allow_nan=False) + '\n') + print('PASS' if result['passed'] else 'FAIL', args.output) + return 0 if result['passed'] else 1 + + +if __name__ == '__main__': + raise SystemExit(main()) diff --git a/experiments/damped_coupled_screen_2026-09-18.json b/experiments/damped_coupled_screen_2026-09-18.json new file mode 100644 index 0000000..ce4d2d4 --- /dev/null +++ b/experiments/damped_coupled_screen_2026-09-18.json @@ -0,0 +1,425 @@ +[ + { + "id": "late_raw_none", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "none" + ] + }, + { + "id": "late_raw_full", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "full" + ] + }, + { + "id": "late_raw_damped", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "damped" + ] + }, + { + "id": "late_raw_near", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "near" + ] + }, + { + "id": "late_curvature_none", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "curvature-coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "none" + ] + }, + { + "id": "late_curvature_full", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "curvature-coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "full" + ] + }, + { + "id": "late_curvature_damped", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "curvature-coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "damped" + ] + }, + { + "id": "late_curvature_near", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "curvature-coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/late_v1/n1p5_late_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "near" + ] + }, + { + "id": "seed_curvature_near", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--max-newton", + "0", + "--damped-coupled-probe", + "--strategy", + "curvature-coupled", + "--probe-state", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/controls_v1/n1p5_seed_eta003/nonlinear_step_physical_base.txt", + "--probe-direction", + "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/nonlinear_step_2026-09-18/controls_v1/n1p5_seed_eta003/nonlinear_step_direction_full.txt", + "--linear-tol", + "0.03", + "--restart", + "160", + "--preconditioner-variant", + "gravity-schur", + "--schur-sweeps", + "1", + "--schur-relaxation", + "0.5", + "--gravity-tol", + "1e-10", + "--gravity-iterations", + "300", + "--gravity-restart", + "80", + "--coupled-tol", + "0.03", + "--coupled-iterations", + "160", + "--coupled-restart", + "160", + "--coupled-max-ratio", + "1", + "--coupled-attempts", + "1", + "--coupled-search", + "near" + ] + } +] diff --git a/experiments/plateau_damped_coupled_diagnostics.hpp b/experiments/plateau_damped_coupled_diagnostics.hpp new file mode 100644 index 0000000..7126544 --- /dev/null +++ b/experiments/plateau_damped_coupled_diagnostics.hpp @@ -0,0 +1,333 @@ +#pragma once + +#include "plateau_coupled_trial_diagnostics.hpp" + +namespace experiment::plateau_diagnostics { +struct DampedCoupledProbeOptions final { + CoupledTrialProbeOptions common{}; + // The uncorrected control still captures the same tangent so its measured + // search cost is comparable without hiding the optional feature's setup. + std::string mode{"near"}; + mean_field::solver::BacktrackedCoupledStep search{}; + mean_field::solver::NearAcceptanceCoupledAttempt attempt{}; + int maximumTrials{13}; +}; + +template +void ProbeDampedCoupledPath(State &state, const mfem::Vector &base, const mfem::Vector &direction, + const std::filesystem::path &output, DampedCoupledProbeOptions options, + Search search, Attempt attempt) { + using namespace mean_field; + namespace d = nonlinear_step_detail; + using Gravity = coupled_trial_detail::GravityWorkspace; + using GravityCorrection = solver::PreparedGravityCurvatureCorrection; + using Transaction = solver::detail::PreparedCoupledTrialTransaction; + const auto comm = state.Problem().GetCommunicator(); + const int size = state.Problem().StateSize(); + int ranks = 0; + MPI_Comm_size(comm, &ranks); + if (ranks != 1 || size != state.Problem().EquationSize() || base.Size() != size || direction.Size() != size || + !d::Finite(base) || !d::Finite(direction) || options.maximumTrials < 1 || options.maximumTrials > 30 || + !std::isfinite(options.common.outerLinearTolerance) || options.common.outerLinearTolerance <= 0. || + options.common.outerLinearTolerance >= 1. || !std::isfinite(options.common.sufficientDecrease) || + options.common.sufficientDecrease <= 0. || options.common.sufficientDecrease >= 1.) + throw std::invalid_argument("Sequential coupled probes require one rank and compatible finite vectors."); + const auto file = [&](const std::string &name) { return output / ("damped_coupled_" + name); }; + auto summary = d::Open(file("summary.csv")); + auto trials = d::Open(file("trials.csv")); + auto candidates = d::Open(file("candidates.csv")); + auto evaluations = d::Open(file("evaluations.csv")); + auto layout = d::Open(file("layout.csv")); + layout << "kind,block,offset,size\n"; + for (const auto &v : state.Problem().GetManifest().valueBlocks()) + layout << "value," << v.stableId << ',' << v.offset << ',' << v.size << '\n'; + for (const auto &v : state.Problem().GetManifest().residualBlocks()) + layout << "residual," << v.stableId << ',' << v.offset << ',' << v.size << '\n'; + layout.flush(); + summary << "metric,value\n"; + trials << "trial,alpha,raw_prepared,raw_reason,raw_stage,gravity_gate,usable,before_norm,before_merit," + "before_accepted,status,attempted,attempts_total,aux_iterations,aux_gate,aux_threshold,aux_defect," + "correction_norm,correction_cap,selected_beta,candidate_count,predicted_decrease,excess_ratio," + "selected_prepared,selected_norm,selected_merit,selected_accepted,selected_original," + "raw_seconds,gravity_seconds,correction_seconds,preparations,restorations,solve_seconds," + "refresh_seconds,verification_seconds,preparation_seconds,aux_core_status,aux_status,aux_reported_gate," + "aux_rhs_norm,aux_reported_true_norm,aux_verification_applied,rhs_saved,correction_saved,verified_vectors_saved," + "gravity_target_defect,gravity_target_threshold,gravity_affine_disagreement\n"; + candidates << "trial,candidate,beta,status,merit,accepted,preparation_seconds\n"; + evaluations << "trial,evaluation,beta,residual_norm,merit,accepted\n"; + const auto metric = [&](const char *key, auto value) { summary << key << ',' << value << '\n'; }; + const double nan = std::numeric_limits::quiet_NaN(); + const auto &normalization = state.normalizedOperator->GetNormalization(); + const mfem::Vector accepted(state.acceptedNormalizedState), acceptedF(state.acceptedNormalizedResidual), + acceptedPhysical(state.AcceptedPhysicalState()), savedDirection(state.normalizedCorrection), + savedPhysicalDirection(state.physicalCorrection), savedVolumeDirection(state.volumeDisplacementDirection), + stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors()); + const double acceptedJacobian = state.acceptedMinimumJacobianDeterminant; + const auto evaluate = [&](const mfem::Vector &v) { + return solver::nonlinear::getMetric(solver::nonlinear::NormalizedResidualMetric{}, v, comm); + }; + const auto prepare = [&](const mfem::Vector &x) { + state.trialNormalizedState = x; + state.normalizedOperator->DenormalizeState(x, state.candidatePhysicalState); + return state.PrepareTrial(); + }; + mfem::Vector baseF(size), action(size); + bool captured = false, finished = false; + const auto finish = [&]() { + finished = true; + if (!prepare(base)) throw std::runtime_error("Sequential coupled probe could not restore its base."); + state.physicalInverse->Refresh(); + mfem::Vector restoredAction(size); + state.normalizedOperator->Mult(direction, restoredAction); + const bool acceptedSame = d::Identical(accepted, state.acceptedNormalizedState) && + d::Identical(acceptedF, state.acceptedNormalizedResidual) && + d::Identical(acceptedPhysical, state.AcceptedPhysicalState()) && + acceptedJacobian == state.acceptedMinimumJacobianDeterminant; + const bool directionsSame = d::Identical(savedDirection, state.normalizedCorrection) && + d::Identical(savedPhysicalDirection, state.physicalCorrection) && + d::Identical(savedVolumeDirection, state.volumeDisplacementDirection); + const bool factorsSame = d::Identical(stateFactors, normalization.StateFactors()) && + d::Identical(residualFactors, normalization.ResidualFactors()); + metric("accepted_buffers_bitwise_identical", acceptedSame); + metric("direction_buffers_bitwise_identical", directionsSame); + metric("normalization_factors_bitwise_identical", factorsSame); + metric("restored_residual_error", captured ? d::Difference(baseF, state.trialNormalizedResidual) : nan); + metric("restored_action_error", captured ? d::Difference(action, restoredAction) : nan); + metric("fresh_operator_and_inverse_current", state.normalizedOperator->IsPrepared() && + state.normalizedInverse->IsCurrent() && state.physicalInverse->IsCurrent()); + d::Save(file("restored_residual.txt"), state.trialNormalizedResidual); + d::Save(file("restored_action.txt"), restoredAction); + summary.flush(); + if (!acceptedSame || !directionsSame || !factorsSame) + throw std::runtime_error("Sequential coupled probe mutated an accepted buffer or normalization."); + }; + metric("candidate_status_evaluated", int(solver::CoupledTrialCandidateStatus::evaluated)); + metric("candidate_status_inadmissible", int(solver::CoupledTrialCandidateStatus::inadmissible)); + metric("candidate_status_non_finite", int(solver::CoupledTrialCandidateStatus::non_finite)); + metric("schema_version", 1); + metric("mode", options.mode); + metric("gravity_path", options.common.gravityPreconditionPath); + metric("attempt_budget", 1); + metric("maximum_trials", options.maximumTrials); + metric("maximum_iterations", options.common.coupled.control.maximumIterations); + metric("restart", options.common.coupled.backend.restartLength); + metric("maximum_beta_steps", std::same_as ? 1 : options.search.maximumSteps); + metric("beta_reduction", options.search.reduction); + metric("maximum_excess_ratio", options.attempt.maximumExcessRatio); + metric("sufficient_decrease", options.common.sufficientDecrease); + metric("relative_tolerance", options.common.coupled.control.relativeTolerance); + metric("absolute_tolerance", options.common.coupled.control.absoluteTolerance); + metric("maximum_correction_ratio", options.common.coupled.maximumCorrectionRatio); + metric("accepted_commits", 0); + metric("shared_iteration_budget", 1); + metric("status_prediction_deferred", int(solver::CoupledTrialAttemptStatus::prediction_deferred)); + metric("status_original_accepted", int(solver::CoupledTrialAttemptStatus::original_accepted)); + metric("status_budget_exhausted", int(solver::CoupledTrialAttemptStatus::budget_exhausted)); + metric("status_kept", int(solver::CoupledTrialAttemptStatus::kept)); + try { + if (!prepare(base)) throw std::runtime_error("Sequential coupled probe base does not prepare."); + state.physicalInverse->Refresh(); + baseF = state.trialNormalizedResidual; + state.normalizedOperator->Mult(direction, action); + captured = true; + const auto baseMetric = evaluate(baseF); + mfem::Vector defect(action); + defect += baseF; + const bool outerGate = d::Norm(defect) <= options.common.outerLinearTolerance * baseMetric.residualNorm; + metric("base_norm", baseMetric.residualNorm); + metric("base_merit", baseMetric.merit); + metric("outer_relative_defect", d::Norm(defect) / baseMetric.residualNorm); + metric("outer_linear_tolerance", options.common.outerLinearTolerance); + metric("outer_verified_gate", outerGate); + d::Save(file("base.txt"), base); + d::Save(file("base_physical.txt"), state.candidatePhysicalState); + d::Save(file("base_residual.txt"), baseF); + d::Save(file("direction.txt"), direction); + d::Save(file("base_action.txt"), action); + d::Save(file("state_factors.txt"), stateFactors); + d::Save(file("residual_factors.txt"), residualFactors); + if (!outerGate) throw std::runtime_error("Sequential probe supplied direction fails the full-J gate."); + Transaction transaction(state, options.common.coupled, 1, search, attempt); + const auto captureStart = d::Clock::now(); + transaction.CaptureBase(baseF, direction); + std::optional gravity; + std::optional gravityCorrection; + if (options.common.gravityPreconditionPath) { + gravity.emplace(state, options.common.gravity); + gravityCorrection.emplace(gravity->subsystem); + gravity->subsystem.Refresh(); + gravityCorrection->CaptureBase(baseF, direction); + } + metric("capture_seconds", d::Seconds(captureStart)); + const auto searchStart = d::Clock::now(); + double operationSeconds = 0.; + int acceptedTrial = -1; + for (int trial = 0; trial < options.maximumTrials; ++trial) { + const double alpha = std::ldexp(1., -trial); + const std::string suffix = "_a" + std::to_string(trial) + ".txt"; + transaction.BeginTrial(); + if (gravity) gravity->BeginGravityAttempt(); + mfem::Vector raw(base), prediction(baseF); + raw.Add(alpha, direction); + prediction.Add(alpha, action); + d::Save(file("raw" + suffix), raw); + d::Save(file("prediction" + suffix), prediction); + const auto rawStart = d::Clock::now(); + auto preparation = prepare(raw); + const double rawSeconds = d::Seconds(rawStart); + const bool rawPrepared = bool(preparation); + if (rawPrepared) d::Save(file("raw_residual" + suffix), state.trialNormalizedResidual); + const int rawReason = rawPrepared ? -1 : int(preparation.error().reason); + const int rawStage = rawPrepared ? -1 : int(preparation.error().stage); + bool gravityGate = !gravity; + double gravitySeconds = 0.; + if (rawPrepared && gravity) { + const auto started = d::Clock::now(); + gravity->Refresh(state); + gravityCorrection->BuildLinearPrediction(alpha, gravity->desiredResidual); + gravity->lastSolve = gravityCorrection->SolveCorrection(alpha, state.trialNormalizedResidual, + gravity->gravityCorrection, options.common.gravity.control); + gravity->trialSolveSucceeded = gravity->lastSolve.Converged(); + if (gravity->trialSolveSucceeded) { + preparation = gravity->ApplyCorrection(state); + gravityGate = bool(preparation) && gravity->TrialSolveSucceeded(); + + } + gravitySeconds = d::Seconds(started); + if (gravityGate) d::Save(file("gravity_correction" + suffix), gravity->fullCorrection); + } + const bool usable = preparation.has_value() && gravityGate; + const mfem::Vector before(state.trialNormalizedState); + auto beforeMetric = solver::nonlinear::MetricEvaluation{nan, nan}; + if (usable) { + beforeMetric = evaluate(state.trialNormalizedResidual); + d::Save(file("before" + suffix), before); + d::Save(file("before_residual" + suffix), state.trialNormalizedResidual); + } + const auto accepts = [&](const auto &m) { + return std::isfinite(m.merit) && m.merit >= 0. && + m.merit <= (1. - options.common.sufficientDecrease * alpha) * baseMetric.merit; + }; + const auto better = [](const auto &a, const auto &b) { + return std::isfinite(a.merit) && a.merit >= 0. && std::isfinite(a.residualNorm) && + a.residualNorm >= 0. && a.merit < b.merit; + }; + auto selected = beforeMetric; + int observed = 0; + struct Observation { mfem::Vector x, f; double beta; solver::nonlinear::MetricEvaluation metric; }; + std::vector observations; + const auto observe = [&](const mfem::Vector &v) { + const auto m = evaluate(v); + // Prediction metrics use other vectors. Copy observations in + // memory here; writing large artifacts must not inflate the + // timed correction's preparation/selection cost. + if (&v == &state.trialNormalizedResidual) { + mfem::Vector delta(state.trialNormalizedState); + delta -= before; + const auto &c = transaction.GetCorrection(); + const double c2 = c * c; + observations.push_back({state.trialNormalizedState, v, c2 > 0. ? (delta * c) / c2 : 0., m}); + } + return m; + }; + const auto previous = transaction.GetDiagnostics(); + double correctionSeconds = 0.; + if (usable && options.mode != "none") { + const auto started = d::Clock::now(); + transaction.TryImprove(state, alpha, preparation, selected, observe, accepts, better, + [](const mfem::Vector &) { return true; }); + correctionSeconds = d::Seconds(started); + } + operationSeconds += rawSeconds + gravitySeconds + correctionSeconds; + const auto &diag = transaction.GetDiagnostics(); + const auto &report = diag.currentSolve; + const bool rhsSaved = report && d::Finite(transaction.GetCore().GetRightHandSide()); + const bool correctionSaved = report && d::Finite(transaction.GetCorrection()); + const bool verifiedSaved = report && report->verificationApplied && + d::Finite(transaction.GetCore().GetVerifiedAction()) && d::Finite(transaction.GetCore().GetVerifiedDefect()); + if (report) { + if (rhsSaved) d::Save(file("rhs" + suffix), transaction.GetCore().GetRightHandSide()); + if (correctionSaved) d::Save(file("correction" + suffix), transaction.GetCorrection()); + if (verifiedSaved) { + d::Save(file("action" + suffix), transaction.GetCore().GetVerifiedAction()); + d::Save(file("defect" + suffix), transaction.GetCore().GetVerifiedDefect()); + } + } + for (const auto &o : observations) { + const auto tag = "_a" + std::to_string(trial) + "_e" + std::to_string(observed) + ".txt"; + d::Save(file("evaluated_state" + tag), o.x); + d::Save(file("evaluated_residual" + tag), o.f); + evaluations << trial << ',' << observed++ << ',' << o.beta << ',' << o.metric.residualNorm + << ',' << o.metric.merit << ',' << accepts(o.metric) << '\n'; + } + int candidate = 0; + for (const auto &c : diag.candidates) + candidates << trial << ',' << candidate++ << ',' << c.beta << ',' << int(c.status) << ',' + << c.merit << ',' << c.accepted << ',' << c.preparationSeconds << '\n'; + if (usable && preparation) { + selected = evaluate(state.trialNormalizedResidual); + d::Save(file("selected" + suffix), state.trialNormalizedState); + d::Save(file("selected_residual" + suffix), state.trialNormalizedResidual); + } + const bool acceptedNow = usable && preparation && accepts(selected); + trials << trial << ',' << alpha << ',' << rawPrepared << ',' << rawReason << ',' << rawStage << ',' + << gravityGate << ',' << usable << ',' << beforeMetric.residualNorm << ',' << beforeMetric.merit << ',' + << (usable && accepts(beforeMetric)) << ',' << int(diag.status) << ',' << bool(report) << ',' + << diag.attempts << ',' << (report ? report->linearSolve.iterations : 0) << ',' + << (report && report->verifiedGate) << ',' << (report ? report->convergenceThreshold : nan) << ',' + << (report ? report->verifiedTrueResidualNorm : nan) << ',' << (report ? report->correctionNorm : nan) << ',' + << (report ? report->maximumCorrectionNorm : nan) << ',' << diag.selectedBeta << ',' << diag.candidateCount << ',' + << diag.predictedDecrease << ',' << diag.predictionExcessRatio << ',' << (usable && bool(preparation)) << ',' + << selected.residualNorm << ',' << selected.merit << ',' << acceptedNow << ',' + << (usable && d::Identical(before, state.trialNormalizedState)) << ',' << rawSeconds << ',' + << gravitySeconds << ',' << correctionSeconds << ',' << diag.preparations - previous.preparations << ',' + << diag.restorations - previous.restorations << ',' << diag.solveSeconds - previous.solveSeconds << ',' + << diag.refreshSeconds - previous.refreshSeconds << ',' << diag.verificationSeconds - previous.verificationSeconds << ',' + << diag.preparationSeconds - previous.preparationSeconds << ',' + << (report ? int(report->status) : -1) << ',' << (report ? int(report->linearSolve.status) : -1) << ',' + << (report && report->linearSolve.Converged()) << ',' << (report ? report->linearSolve.rightHandSideNorm : nan) << ',' + << (report ? report->linearSolve.trueResidualNorm : nan) << ',' << (report && report->verificationApplied) << ',' + << rhsSaved << ',' << correctionSaved << ',' << verifiedSaved << ',' + << (gravity && rawPrepared ? gravity->lastPreparedDefect : nan) << ',' + << (gravity && rawPrepared ? gravity->lastPreparedTargetThreshold : nan) << ',' + << (gravity && rawPrepared ? gravity->lastAffineDisagreement : nan) << '\n'; + trials.flush(); candidates.flush(); evaluations.flush(); + if (acceptedNow) { + acceptedTrial = trial; + metric("accepted_alpha", alpha); + metric("accepted_norm", selected.residualNorm); + metric("accepted_beta", diag.selectedBeta); + d::Save(output / ("coupled_trial_selected_physical" + suffix), state.candidatePhysicalState); + break; + } + } + metric("search_wall_seconds", d::Seconds(searchStart)); + metric("search_operation_seconds", operationSeconds); + metric("accepted_trial", acceptedTrial); + const auto &diag = transaction.GetDiagnostics(); + metric("attempts", diag.attempts); + metric("aux_iterations", diag.iterations); + metric("deferred_attempts", diag.deferredAttempts); + metric("probe_exception", 0); + } catch (...) { + const auto error = std::current_exception(); + metric("probe_exception", 1); + if (!finished) finish(); + std::rethrow_exception(error); + } + finish(); +} + +template +void ProbeDampedCoupled(State &state, const mfem::Vector &base, const mfem::Vector &direction, + const std::filesystem::path &output, DampedCoupledProbeOptions options) { + using namespace mean_field::solver; + if (options.mode == "none" || options.mode == "full") + ProbeDampedCoupledPath(state, base, direction, output, options, FullCoupledStep{}, EveryRejectedCoupledAttempt{}); + else if (options.mode == "damped") + ProbeDampedCoupledPath(state, base, direction, output, options, options.search, EveryRejectedCoupledAttempt{}); + else if (options.mode == "near") + ProbeDampedCoupledPath(state, base, direction, output, options, options.search, options.attempt); + else throw std::invalid_argument("Unknown sequential coupled probe mode."); +} +} // namespace experiment::plateau_diagnostics diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 941558c..2f8e3ec 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -26,6 +26,7 @@ import mean_field; #include "plateau_nonlinear_step_diagnostics.hpp" #include "plateau_coupled_trial_diagnostics.hpp" #include "plateau_coupled_trial_physical.hpp" +#include "plateau_damped_coupled_diagnostics.hpp" #include "plateau_gravity_application_diagnostics.hpp" #include "plateau_positivity_strategy_diagnostics.hpp" #include "plateau_polynomial_certificate.hpp" @@ -40,6 +41,10 @@ struct Options { bool strategyProbe = false, positivityProbe = false, positivityCertificate = false; bool preconditionerProbe = false, nonlinearStepProbe = false; bool coupledTrialProbe = false, coupledSeedControl = false; + bool dampedCoupledProbe = false; + std::string coupledSearch = "near"; + int coupledBetaSteps = 4; + double coupledBetaReduction = .5, coupledExcessRatio = 1.; double coupledTolerance = .03, coupledMaximumRatio = 1.; int coupledIterations = 160, coupledRestart = 160, coupledAttempts = 1; std::string preconditionerVariant = "independent"; @@ -108,6 +113,16 @@ Options Parse(int argc, char **argv) { o.preconditionerProbe = true; else if (a == "--nonlinear-step-probe") o.nonlinearStepProbe = true; + else if (a == "--damped-coupled-probe") + o.dampedCoupledProbe = true; + else if (a == "--coupled-search") + o.coupledSearch = v(); + else if (a == "--coupled-beta-steps") + o.coupledBetaSteps = std::stoi(v()); + else if (a == "--coupled-beta-reduction") + o.coupledBetaReduction = std::stod(v()); + else if (a == "--coupled-excess-ratio") + o.coupledExcessRatio = std::stod(v()); else if (a == "--coupled-trial-probe") o.coupledTrialProbe = true; else if (a == "--coupled-seed-control") @@ -158,17 +173,17 @@ Options Parse(int argc, char **argv) { throw std::invalid_argument("Invalid options"); if (o.probeRestart != 40 && o.probeRestart != 160) throw std::invalid_argument("--probe-restart must be 40 or 160"); - if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.coupledTrialProbe) && (o.probeState.empty() || o.probeDirection.empty())) + if ((o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.coupledTrialProbe || o.dampedCoupledProbe) && (o.probeState.empty() || o.probeDirection.empty())) throw std::invalid_argument("Focused probes require --probe-state and --probe-direction"); if (o.positivityCertificate && o.probeState.empty()) throw std::invalid_argument("--positivity-certificate requires --probe-state"); if (o.positivityCertificate && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) + o.strategyProbe || o.positivityProbe || o.preconditionerProbe || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options"); if (o.preconditionerProbe) { experiment::plateau_diagnostics::ParsePreconditionerBenchmarkVariant(o.preconditionerVariant); if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || o.coupledTrialProbe || !o.probeDirection.empty() || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || o.coupledTrialProbe || o.dampedCoupledProbe || !o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline") throw std::invalid_argument("--preconditioner-probe requires --max-newton 0 and no nonlinear/probe strategies"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) @@ -176,22 +191,25 @@ Options Parse(int argc, char **argv) { } if (o.nonlinearStepProbe) { if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.coupledTrialProbe || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.coupledTrialProbe || o.dampedCoupledProbe || !o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline") throw std::invalid_argument("--nonlinear-step-probe requires --max-newton 0 and no other probes or strategy; it compares raw and curvature paths itself"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) throw std::invalid_argument("A nonlinear-step probe requires a finite positive --linear-tol"); } - if (o.coupledTrialProbe) { + if (o.coupledTrialProbe || o.dampedCoupledProbe) { if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || - o.nonlinearStepProbe || !o.replayCandidate.empty() || !o.captureIterations.empty() || - o.coupledAttempts != 1 || (o.strategy != "coupled" && o.strategy != "curvature-coupled")) + o.nonlinearStepProbe || (o.coupledTrialProbe && o.dampedCoupledProbe) || !o.replayCandidate.empty() || !o.captureIterations.empty() || + o.coupledAttempts != 1 || (o.strategy != "coupled" && o.strategy != "curvature-coupled" && o.strategy != "backtracked-coupled" && o.strategy != "curvature-backtracked-coupled")) throw std::invalid_argument("--coupled-trial-probe requires --max-newton 0 and strategy coupled or curvature-coupled, without other probes"); + if (o.coupledTrialProbe && o.strategy != "coupled" && o.strategy != "curvature-coupled") + throw std::invalid_argument("The independent coupled probe requires its original full-step strategy"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.) throw std::invalid_argument("A coupled trial probe requires a finite positive outer --linear-tol"); - } else if (o.coupledSeedControl) { + } + if (o.coupledSeedControl && !o.coupledTrialProbe) { throw std::invalid_argument("--coupled-seed-control requires --coupled-trial-probe"); } if (!o.preconditionerProbe) { @@ -200,12 +218,12 @@ Options Parse(int argc, char **argv) { throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, stellar-ldu, or stellar-ldu-gravity-schur"); if (o.preconditionerVariant != "independent") { if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || - o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe) || - (!o.probeDirection.empty() && !o.coupledTrialProbe) || + o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe && !o.coupledTrialProbe && !o.dampedCoupledProbe) || + (!o.probeDirection.empty() && !o.coupledTrialProbe && !o.dampedCoupledProbe) || !o.replayCandidate.empty()) throw std::invalid_argument("A custom Newton preconditioner cannot be combined with frozen or certificate probes"); if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" && - o.strategy != "positivity-curvature" && o.strategy != "coupled" && o.strategy != "curvature-coupled") + o.strategy != "positivity-curvature" && o.strategy != "coupled" && o.strategy != "curvature-coupled" && o.strategy != "backtracked-coupled" && o.strategy != "curvature-backtracked-coupled") throw std::invalid_argument("Unsupported custom Newton preconditioner strategy"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0 || o.schurCorrectionSweeps < 0 || !std::isfinite(o.schurRelaxation) || o.schurRelaxation <= 0.0) @@ -216,13 +234,19 @@ Options Parse(int argc, char **argv) { throw std::invalid_argument("A supplied direction or candidate requires --probe-state"); if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" && o.strategy != "reduced" && o.strategy != "positivity-curvature" && o.strategy != "positivity-reduced" && - o.strategy != "coupled" && o.strategy != "curvature-coupled") + o.strategy != "coupled" && o.strategy != "curvature-coupled" && o.strategy != "backtracked-coupled" && o.strategy != "curvature-backtracked-coupled") throw std::invalid_argument("Unknown numerical strategy"); if (!std::isfinite(o.gravityTolerance) || o.gravityTolerance <= 0.0 || o.gravityIterations < 1 || o.gravityRestart < 1) throw std::invalid_argument("Invalid gravity solve controls"); if (!std::isfinite(o.coupledTolerance) || o.coupledTolerance <= 0. || o.coupledIterations < 1 || o.coupledRestart < 1 || !std::isfinite(o.coupledMaximumRatio) || o.coupledMaximumRatio <= 0. || o.coupledAttempts < 1) throw std::invalid_argument("Invalid coupled correction controls"); + if (o.coupledSearch != "none" && o.coupledSearch != "full" && o.coupledSearch != "damped" && o.coupledSearch != "near") + throw std::invalid_argument("Unknown coupled search mode"); + if (!o.dampedCoupledProbe && o.coupledSearch != "near") + throw std::invalid_argument("--coupled-search selects controls only in --damped-coupled-probe"); + solver::BacktrackedCoupledStep{.maximumSteps=o.coupledBetaSteps, .reduction=o.coupledBetaReduction}.Validate(); + solver::NearAcceptanceCoupledAttempt{.maximumExcessRatio=o.coupledExcessRatio}.Validate(); if (std::filesystem::exists(o.output)) throw std::runtime_error("Output directory already exists"); std::filesystem::create_directories(o.output); @@ -488,7 +512,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nM=1\nR_reference=1\nxi1=" << scales.xi << "\nmu1=" << scales.mu << "\ncentral_phase=" << !o.noCentral << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart << "\nfrozen_probe_backend=" << (o.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" : - o.coupledTrialProbe ? "supplied_outer_direction_verified_with_full_J_and_auxiliary_fgmres" : + (o.coupledTrialProbe || o.dampedCoupledProbe) ? "supplied_outer_direction_verified_with_full_J_and_auxiliary_fgmres" : o.nonlinearStepProbe ? "prepared_fgmres_cold_with_independent_true_gate" : o.positivityCertificate ? "polynomial_certificate_no_linear_solve" : o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve") @@ -498,6 +522,11 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart << "\npreconditioner_probe=" << o.preconditionerProbe << "\nnonlinear_step_probe=" << o.nonlinearStepProbe + << "\ndamped_coupled_probe=" << o.dampedCoupledProbe + << "\ncoupled_search=" << o.coupledSearch + << "\ncoupled_beta_steps=" << o.coupledBetaSteps + << "\ncoupled_beta_reduction=" << o.coupledBetaReduction + << "\ncoupled_excess_ratio=" << o.coupledExcessRatio << "\ncoupled_trial_probe=" << o.coupledTrialProbe << "\ncoupled_seed_control=" << o.coupledSeedControl << "\ncoupled_tolerance=" << o.coupledTolerance @@ -513,6 +542,34 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; metadata.flush(); + if (o.dampedCoupledProbe) { + metadata << "coupled_trial_input=physical_snapshot_and_supplied_normalized_direction\n" + << "coupled_trial_admissibility=sampled_enthalpy_without_repair\n" + << "coupled_trial_alphas_are_independent=0\n" + << "probe_state=" << o.probeState << "\nprobe_direction=" << o.probeDirection << '\n'; + metadata.flush(); + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) { + const auto physical = LoadVector(o.probeState, s.AcceptedPhysicalState().Size()); + mfem::Vector base; + s.normalizedOperator->NormalizeState(physical, base); + const auto direction = LoadVector(o.probeDirection, base.Size()); + experiment::plateau_diagnostics::ProbeDampedCoupled(s, base, direction, o.output, { + .common = {.outerLinearTolerance = o.linearTolerance, + .gravity = {.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14, + .maximumIterations = o.gravityIterations}, + .backend = {.restartLength = o.gravityRestart, .printLevel = -1}}, + .coupled = {.control = {.relativeTolerance = o.coupledTolerance, .absoluteTolerance = 1e-14, + .maximumIterations = o.coupledIterations}, + .backend = {.restartLength = o.coupledRestart, .printLevel = -1}, + .maximumCorrectionRatio = o.coupledMaximumRatio}, + .gravityPreconditionPath = o.strategy.starts_with("curvature-")}, + .mode = o.coupledSearch, + .search = {.maximumSteps = o.coupledBetaSteps, .reduction = o.coupledBetaReduction}, + .attempt = {.maximumExcessRatio = o.coupledExcessRatio}}); + experiment::plateau_diagnostics::MeasureCoupledTrialSnapshots(s, f, o.output, o.n, scales.K); + }); + return 0; + } if (o.coupledTrialProbe) { metadata << "coupled_trial_input=physical_snapshot_and_supplied_normalized_direction\n" << "coupled_trial_admissibility=sampled_enthalpy_without_repair\n" @@ -675,6 +732,10 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { "preparation_seconds_total,refresh_seconds_total,original_merit,corrected_merit,has_solve," "linear_status,linear_iterations,rhs_norm,reported_true_norm,verified_true_norm,threshold,verified_gate," "correction_norm,maximum_correction_norm\n"; + auto coupledSearchLog = File(o.output / "coupled_search.csv"); + coupledSearchLog << "iteration,trial,alpha,status,deferred_total,selected_beta,candidate_count,base_merit,predicted_merit,predicted_decrease,excess_ratio\n"; + auto coupledCandidatesLog = File(o.output / "coupled_candidates.csv"); + coupledCandidatesLog << "iteration,trial,candidate,beta,status,merit,accepted,preparation_seconds\n"; auto stateConsistency = File(o.output / "trial_state_consistency.csv"); stateConsistency << "iteration,trial,normalized_physical_roundtrip_error,relative_error\n"; mfem::Vector captureBase; @@ -721,6 +782,16 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { coupledLog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << (e.numericalStrategies != nullptr); if (e.numericalStrategies) { const auto &c = e.numericalStrategies->coupledTrial; + coupledSearchLog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << int(c.status) << ',' + << c.deferredAttempts << ',' << c.selectedBeta << ',' << c.candidateCount << ',' << c.baseMerit << ',' + << c.linearPredictionMerit << ',' << c.predictedDecrease << ',' << c.predictionExcessRatio << '\n'; + int candidateIndex = 0; + for (const auto &candidate : c.candidates) + coupledCandidatesLog << e.iteration << ',' << e.trial << ',' << candidateIndex++ << ',' + << candidate.beta << ',' << int(candidate.status) << ',' << candidate.merit << ',' + << candidate.accepted << ',' << candidate.preparationSeconds << '\n'; + coupledSearchLog.flush(); + coupledCandidatesLog.flush(); coupledLog << ',' << static_cast(c.status) << ',' << c.attempts << ',' << c.kept << ',' << c.fallbacks << ',' << c.iterations << ',' << c.preparations << ',' << c.restorations << ',' << c.solveSeconds << ',' << c.verificationSeconds << ',' << c.preparationSeconds << ',' @@ -836,6 +907,12 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { auto report = [&] { if (o.strategy == "positivity") return run(nl::NewtonNumerics{nl::PolynomialEnthalpy{}, nl::MonolithicGravity{}}); if (o.strategy == "curvature") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}}); + const nl::BacktrackedCoupledCurvatureCorrection backtracked{ + .correction = coupled.correction, .maximumAttemptsPerIteration = o.coupledAttempts, + .search = {.maximumSteps = o.coupledBetaSteps, .reduction = o.coupledBetaReduction}, + .attempt = {.maximumExcessRatio = o.coupledExcessRatio}}; + if (o.strategy == "backtracked-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::MonolithicGravity{}, backtracked}); + if (o.strategy == "curvature-backtracked-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}, backtracked}); if (o.strategy == "coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::MonolithicGravity{}, coupled}); if (o.strategy == "curvature-coupled") return run(nl::NewtonNumerics{nl::SampledEnthalpy{}, nl::GravityCurvatureCorrection{gravity}, coupled}); if constexpr (FullDefaultStack) @@ -872,6 +949,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { metric("last_lifted_relative_true_residual", strategyDiagnostics.lastLiftedRelativeTrueResidual); const auto &coupledDiagnostics = strategyDiagnostics.coupledTrial; metric("coupled_attempts", coupledDiagnostics.attempts); + metric("coupled_deferred_attempts", coupledDiagnostics.deferredAttempts); metric("coupled_kept", coupledDiagnostics.kept); metric("coupled_fallbacks", coupledDiagnostics.fallbacks); metric("coupled_iterations", coupledDiagnostics.iterations); diff --git a/experiments/test_damped_coupled_trial_checker.py b/experiments/test_damped_coupled_trial_checker.py new file mode 100644 index 0000000..cc39eb1 --- /dev/null +++ b/experiments/test_damped_coupled_trial_checker.py @@ -0,0 +1,248 @@ +#!/usr/bin/env python3 +"""Small synthetic searches exercise certificates, selection and provenance.""" +import csv +import json +import math +from pathlib import Path +import tempfile +import unittest + +import check_damped_coupled_trial as c +from test_coupled_trial_checker import csv_file, vec_file + + +class DampedCheckerTests(unittest.TestCase): + def setUp(self): + self.temp = tempfile.TemporaryDirectory(); self.addCleanup(self.temp.cleanup) + self.root = Path(self.temp.name).resolve(); self.folder = self.root / 'synthetic'; self.folder.mkdir() + self.make_case() + + def make_case(self, kind='damped'): + for path in self.folder.iterdir(): path.unlink() + pfx, folder, nan = c.P, self.folder, math.nan + mode = 'near' if kind in ('near', 'raw_failure') else ('none' if kind == 'none' else ('full' if kind in ('full', 'budget', 'linear_failure', 'cap') else 'damped')) + x, p, f = [1., 2., 3.], [-1., -1., -1.], [1., 1., 1.] + sf, rf = [2.] * 3, [3.] * 3 + save = lambda name, value: vec_file(folder / (pfx + name + '.txt'), value) + for name, values in {'base': x, 'base_physical': [v / 2 for v in x], 'base_residual': f, 'direction': p, + 'base_action': p, 'state_factors': sf, 'residual_factors': rf, 'restored_residual': f, 'restored_action': p}.items(): save(name, values) + layout = [{'kind': typ, 'block': name, 'offset': i, 'size': 1} for typ, names in ( + ('value', ('density', 'shape', 'enthalpy')), ('residual', ('closure', 'shape', 'hydrostatic'))) for i, name in enumerate(names)] + csv_file(folder / (pfx + 'layout.csv'), layout) + max_steps = 1 if mode in ('none', 'full') else 4 + cap_ratio = .01 if kind == 'cap' else 1. + records, candidates, evaluations = [], [], [] + attempts = iterations = deferred = 0 + plans = ['attempt'] + if kind == 'none': plans = ['accepted'] + elif kind in ('full', 'linear_failure', 'cap', 'best'): plans += ['accepted'] + elif kind == 'budget': plans += ['budget', 'accepted'] + elif kind == 'near': plans = ['deferred', 'attempt'] + elif kind == 'raw_failure': plans = ['raw_failure', 'accepted'] + for i, plan in enumerate(plans): + alpha, suffix = 2. ** -i, '_a' + str(i) + raw = [v + alpha * w for v, w in zip(x, p)]; prediction = [1. - alpha] * 3 + save('raw' + suffix, raw); save('prediction' + suffix, prediction) + usable = plan != 'raw_failure'; attempted = plan == 'attempt'; before_accepted = plan == 'accepted' + before_f = [.5] * 3 if before_accepted else ([3., 0., 0.] if plan == 'deferred' else [2., 0., 0.]) + if usable: + save('raw_residual' + suffix, before_f); save('before' + suffix, raw); save('before_residual' + suffix, before_f) + status = {'accepted': 'not_attempted' if mode == 'none' else 'original_accepted', 'budget': 'budget_exhausted', + 'deferred': 'prediction_deferred', 'raw_failure': 'not_attempted'}.get(plan, 'kept') + selected, selected_f, selected_beta = raw, before_f, nan + core, count, preps, restores = -1, 0, 0, 0 + correction = [-.5, 0., 0.]; rhs = [a - b for a, b in zip(prediction, before_f)] + action = [0., 0., 0.] if kind == 'linear_failure' else rhs + defect = c.minus(action, rhs); gate = c.norm(defect) <= .03 * c.norm(rhs) + if attempted: + attempts += 1; iterations += 2 + core = 1 if kind == 'linear_failure' else (2 if kind == 'cap' else 0) + for name, values in {'rhs': rhs, 'correction': correction, 'action': action, 'defect': defect}.items(): save(name + suffix, values) + if core: + status = 'linear_failure' if core == 1 else 'norm_limit' + else: + proposal_norms = [2.1] if mode == 'full' else ([1.8, 1.9, 2., 2.1] if kind == 'best' else [2.1, 1.]) + observed = [] + for k, value in enumerate(proposal_norms): + beta = .5 ** k; state = [v + beta * w for v, w in zip(raw, correction)]; residual = [value, 0., 0.] + failed = kind == 'candidate_failure' and k == 0 + cm = .5 * value * value; accepted = cm <= (1 - 1.e-4 * alpha) * 1.5 + candidates.append({'trial': i, 'candidate': k, 'beta': beta, 'status': 1 if failed else 0, + 'merit': nan if failed else cm, 'accepted': int(accepted and not failed), 'preparation_seconds': .05}) + count += 1; preps += 1 + if not failed: observed.append((beta, state, residual)) + if kind == 'best': + selected_beta = 1.; selected = [v + w for v, w in zip(raw, correction)]; selected_f = [1.8, 0., 0.] + observed.append((1., selected, selected_f)); restores = 1; preps += 1 + elif mode == 'full': + status = 'not_improved'; observed.append((0., raw, before_f)); restores = 1; preps += 1 + else: + selected_beta, selected, selected_f = observed[-1] + for e, (beta, state, residual) in enumerate(observed): + tag = suffix + '_e' + str(e) + save('evaluated_state' + tag, state); save('evaluated_residual' + tag, residual) + cm = c.merit(residual) + evaluations.append({'trial': i, 'evaluation': e, 'beta': beta, 'residual_norm': c.norm(residual), + 'merit': cm, 'accepted': int(cm <= (1 - 1.e-4 * alpha) * 1.5)}) + deferred += plan == 'deferred' + accepted_now = usable and c.merit(selected_f) <= (1 - 1.e-4 * alpha) * 1.5 + if usable: save('selected' + suffix, selected); save('selected_residual' + suffix, selected_f) + if accepted_now: vec_file(folder / ('coupled_trial_selected_physical' + suffix + '.txt'), [v / 2 for v in selected]) + gate_used = mode == 'near' and usable and not before_accepted and plan != 'budget' + _, decrease, excess = c.eligibility(1.5, c.merit(prediction), c.merit(before_f), 1.) + solve, verification = (.1, .01) if attempted else (0., 0.) + refresh, prep_seconds = (.02 * (1 + restores) if attempted else 0.), .05 * preps + correction_seconds = solve + verification + refresh + prep_seconds + (.01 if attempted else 0.) + records.append({'trial': i, 'alpha': alpha, 'raw_prepared': int(usable), 'raw_reason': -1 if usable else 2, + 'raw_stage': -1 if usable else 6, 'gravity_gate': 1, 'usable': int(usable), + 'before_norm': c.norm(before_f) if usable else nan, 'before_merit': c.merit(before_f) if usable else nan, + 'before_accepted': int(before_accepted), 'status': c.STATUS.index(status), 'attempted': int(attempted), + 'attempts_total': attempts, 'aux_iterations': 2 if attempted else 0, 'aux_gate': int(attempted and gate), + 'aux_threshold': .03 * c.norm(rhs) if attempted else nan, 'aux_defect': c.norm(defect) if attempted else nan, + 'correction_norm': .5 if attempted else nan, 'correction_cap': cap_ratio * alpha * math.sqrt(3) if attempted else nan, + 'selected_beta': selected_beta, 'candidate_count': count, 'predicted_decrease': decrease if gate_used else nan, + 'excess_ratio': excess if gate_used else nan, 'selected_prepared': int(usable), + 'selected_norm': c.norm(selected_f) if usable else nan, 'selected_merit': c.merit(selected_f) if usable else nan, + 'selected_accepted': int(accepted_now), 'selected_original': int(usable and selected == raw), + 'raw_seconds': .1, 'gravity_seconds': 0., 'correction_seconds': correction_seconds, 'preparations': preps, + 'restorations': restores, 'solve_seconds': solve, 'refresh_seconds': refresh, 'verification_seconds': verification, + 'preparation_seconds': prep_seconds, 'aux_core_status': core, 'aux_status': 0 if attempted else -1, + 'aux_reported_gate': int(attempted), 'aux_rhs_norm': c.norm(rhs) if attempted else nan, + 'aux_reported_true_norm': 0. if attempted else nan, 'aux_verification_applied': int(attempted), + 'rhs_saved': int(attempted), 'correction_saved': int(attempted), 'verified_vectors_saved': int(attempted), + 'gravity_target_defect': nan, 'gravity_target_threshold': nan, 'gravity_affine_disagreement': nan}) + csv_file(folder / (pfx + 'trials.csv'), records) + csv_file(folder / (pfx + 'candidates.csv'), candidates, ['trial', 'candidate', 'beta', 'status', 'merit', 'accepted', 'preparation_seconds']) + csv_file(folder / (pfx + 'evaluations.csv'), evaluations, ['trial', 'evaluation', 'beta', 'residual_norm', 'merit', 'accepted']) + last = records[-1] + summary = {'schema_version': 1, 'mode': mode, 'gravity_path': 0, 'attempt_budget': 1, 'maximum_trials': 13, + 'maximum_iterations': 160, 'restart': 160, 'maximum_beta_steps': max_steps, 'beta_reduction': .5, + 'maximum_excess_ratio': 1., 'sufficient_decrease': 1.e-4, 'relative_tolerance': .03, 'absolute_tolerance': 1.e-14, + 'maximum_correction_ratio': cap_ratio, 'accepted_commits': 0, 'shared_iteration_budget': 1, + 'base_norm': math.sqrt(3), 'base_merit': 1.5, 'outer_relative_defect': 0., 'outer_linear_tolerance': .03, + 'outer_verified_gate': 1, 'capture_seconds': .1, 'search_wall_seconds': 10., + 'search_operation_seconds': sum(row['raw_seconds'] + row['correction_seconds'] for row in records), + 'accepted_trial': len(records) - 1, 'accepted_alpha': last['alpha'], 'accepted_norm': last['selected_norm'], + 'accepted_beta': last['selected_beta'], 'attempts': attempts, 'aux_iterations': iterations, 'deferred_attempts': deferred, + 'probe_exception': 0, 'accepted_buffers_bitwise_identical': 1, 'direction_buffers_bitwise_identical': 1, + 'normalization_factors_bitwise_identical': 1, 'restored_residual_error': 0., 'restored_action_error': 0., + 'fresh_operator_and_inverse_current': 1} + summary.update({'status_' + name: c.STATUS.index(name) for name in ('original_accepted', 'budget_exhausted', 'kept', 'prediction_deferred')}) + summary.update({'candidate_status_' + name: i for i, name in enumerate(('evaluated', 'inadmissible', 'non_finite'))}) + csv_file(folder / (pfx + 'summary.csv'), [{'metric': key, 'value': value} for key, value in summary.items()]) + metadata = {'n': 1.5, 'K': 1., 'M': 1., 'J': 0., 'central_phase': 0, 'damped_coupled_probe': 1, + 'coupled_trial_alphas_are_independent': 0, 'coupled_trial_admissibility': 'sampled_enthalpy_without_repair', + 'coupled_search': mode, 'coupled_beta_steps': 4, 'coupled_beta_reduction': .5, 'coupled_excess_ratio': 1., + 'linear_tol': .03, 'coupled_tolerance': .03, 'coupled_iterations': 160, 'coupled_restart': 160, + 'coupled_maximum_ratio': cap_ratio, 'coupled_attempts_per_iteration': 1, 'strategy': 'coupled'} + (folder / 'metadata.txt').write_text(''.join(f'{key}={value}\n' for key, value in metadata.items())) + self.physical(last) + self.make_campaign(mode) + + def physical(self, selected): + records = [] + def add(sample, source, q, key, value): records.append({'sample': sample, 'source': source, 'quadrature': q, 'metric': key, 'value': value}) + for key in ('accepted_buffers_unchanged_before_restoration', 'accepted_buffers_unchanged_after_restoration', + 'normalization_factors_unchanged', 'direction_buffers_unchanged', 'accepted_operator_and_inverse_current'): + add(-1, 'restoration', 0, key, 1) + add(-1, 'replay', 0, 'selected_snapshots_measured', 1) + i = selected['trial']; source = f'coupled_trial_selected_physical_a{i}.txt' + add(i, source, 0, 'fresh_preparation_passed', 1) + for q in (14, 18): + values = {'quadrature_order': q, 'polytropic_index': 1.5, 'polytropic_constant': 1., 'target_mass': 1., + 'stellar_samples': 10, 'valid_stellar_samples': 10, 'invalid_stellar_samples': 0, 'negative_density_samples': 0, + 'negative_enthalpy_samples': 0, 'invalid_density_eos_samples': 0, 'invalid_enthalpy_eos_samples': 0, + 'volume': 1., 'complete_volume_measurement': 1, 'pressure_from_enthalpy_defined_everywhere': 1, + 'pressure_from_density_defined_everywhere': 1, 'saved_prepared_l2_difference': 0., 'saved_prepared_maximum_difference': 0., + 'saved_prepared_maximum_scaled_difference': 0., 'normalized_bordered_residual': selected['selected_norm'], + 'normalized_unbordered_residual': selected['selected_norm'], 'normalized_border_action': 0., 'central_border': 0., + 'mass': 1., 'mass_absolute_error': 0., 'minimum_density': 0., 'minimum_enthalpy': 0., + 'binding_energy': -3., 'force_binding_energy': -3., 'kinetic_energy': 0., 'pressure_integral': 1., 'enthalpy_pressure_integral': 1.} + values.update({prefix + force + 'virial_' + name: 0. for prefix in ('', 'enthalpy_') for force in ('', 'force_') for name in ('signed', 'error')}) + for key, value in values.items(): add(i, source, q, key, value) + csv_file(self.folder / 'coupled_trial_physical.csv', records) + + def make_campaign(self, mode): + root = self.root + state, direction, mesh = root / 'state.txt', root / 'direction.txt', root / 'mesh' + vec_file(state, [.5, 1., 1.5]); vec_file(direction, [-1.] * 3); mesh.write_text('synthetic mesh') + (root / 'plateau_diagnostic').write_text('synthetic executable') + arguments = ['--damped-coupled-probe', '--coupled-search', mode, '--strategy', 'coupled', '--no-central', + '--probe-state', str(state), '--probe-direction', str(direction)] + cases = [{'id': 'synthetic', 'arguments': arguments}] + (root / 'cases.json').write_text(json.dumps(cases)) + (root / 'provenance.json').write_text(json.dumps({'cases': cases, 'commit': 'synthetic', 'source': str(root), + 'mesh': str(mesh), 'mesh_sha256': c.digest(mesh), 'binary_sha256': c.digest(root / 'plateau_diagnostic')})) + (root / 'synthetic_run.json').write_text(json.dumps({'id': 'synthetic', 'finished_utc': 'synthetic', 'exit_code': 0, + 'command': [str(root / 'plateau_diagnostic'), *arguments, '--mesh', str(mesh), '--output', str(self.folder)], + 'input_state_hashes': {str(path): c.digest(path) for path in (state, direction)}, + 'output_hashes': {path.name: c.digest(path) for path in self.folder.iterdir()}})) + + def audit(self): return c.inspect_case(self.folder) + def assert_passes(self): + result = self.audit(); self.assertTrue(result['passed'], [row for row in result['checks'] if not row['passed']]) + def mutate_row(self, file, key, value, row=0): + path = self.folder / (c.P + file + '.csv'); records = c.rows(path); records[row][key] = value; csv_file(path, records) + + def add_gravity_path(self): + folder = self.folder; delta = [0., .01, .02] + for name in ('before_a0', 'selected_a0', 'evaluated_state_a0_e0', 'evaluated_state_a0_e1'): + path = folder / (c.P + name + '.txt'); values = c.vector(path) + vec_file(path, [a + b for a, b in zip(values, delta)]) + vec_file(folder / (c.P + 'gravity_correction_a0.txt'), delta) + vec_file(folder / (c.P + 'raw_residual_a0.txt'), [2., .2, -.1]) + selected = c.vector(folder / (c.P + 'selected_a0.txt')) + vec_file(folder / 'coupled_trial_selected_physical_a0.txt', [value / 2 for value in selected]) + layout = [{'kind': kind, 'block': name, 'offset': i, 'size': 1} for kind, names in ( + ('value', ('density', 'gravity_gradient', 'gravity_potential')), + ('residual', ('closure', 'gravity_gradient_relation', 'poisson_balance'))) for i, name in enumerate(names)] + csv_file(folder / (c.P + 'layout.csv'), layout) + for key, value in (('gravity_target_defect', 0.), ('gravity_target_threshold', 1.e-12), ('gravity_affine_disagreement', 0.)): + self.mutate_row('trials', key, value) + summary = c.rows(folder / (c.P + 'summary.csv')) + next(row for row in summary if row['metric'] == 'gravity_path')['value'] = 1 + csv_file(folder / (c.P + 'summary.csv'), summary) + path = folder / 'metadata.txt'; path.write_text(path.read_text().replace('strategy=coupled', 'strategy=curvature-coupled')) + + def test_damped_half_step(self): self.assert_passes() + def test_near_gate_preserves_budget(self): self.make_case('near'); self.assert_passes() + def test_full_fallback(self): self.make_case('full'); self.assert_passes() + def test_shared_budget_exhaustion(self): self.make_case('budget'); self.assert_passes() + def test_best_beta_replay(self): self.make_case('best'); self.assert_passes() + def test_unprepared_beta_skipped(self): self.make_case('candidate_failure'); self.assert_passes() + def test_no_correction_control(self): self.make_case('none'); self.assert_passes() + def test_pressure_failure_preserves_budget(self): self.make_case('raw_failure'); self.assert_passes() + def test_failed_full_true_gate(self): self.make_case('linear_failure'); self.assert_passes() + def test_full_cap_rejects_before_damping(self): self.make_case('cap'); self.assert_passes() + def test_gravity_path_uses_full_coupled_target(self): self.add_gravity_path(); self.assert_passes() + def test_corrupt_gravity_target_fails(self): + self.add_gravity_path(); self.mutate_row('trials', 'gravity_target_defect', 1.); self.assertFalse(self.audit()['passed']) + def test_provenance_passes(self): self.assertTrue(c.inspect_campaign(self.root)['passed']) + def test_missing_true_action_fails(self): (self.folder / (c.P + 'action_a0.txt')).unlink(); self.assertFalse(self.audit()['passed']) + def test_scaled_cannot_replace_verified_c(self): + vec_file(self.folder / (c.P + 'correction_a0.txt'), [-.25, 0., 0.]); self.assertFalse(self.audit()['passed']) + def test_corrupt_rhs_fails(self): + vec_file(self.folder / (c.P + 'rhs_a0.txt'), [-1., 0., 0.]); self.assertFalse(self.audit()['passed']) + def test_reset_budget_fails(self): + self.make_case('budget'); self.mutate_row('trials', 'attempts_total', 0, 1); self.assertFalse(self.audit()['passed']) + def test_deferred_gate_corruption_fails(self): + self.make_case('near'); self.mutate_row('trials', 'excess_ratio', .5); self.assertFalse(self.audit()['passed']) + def test_wrong_beta_fails(self): + self.mutate_row('candidates', 'beta', .25, 1); self.assertFalse(self.audit()['passed']) + def test_later_beta_after_acceptance_fails(self): + self.mutate_row('candidates', 'accepted', 1, 0); self.assertFalse(self.audit()['passed']) + def test_selected_stale_residual_fails(self): + vec_file(self.folder / (c.P + 'selected_residual_a0.txt'), [2.1, 0., 0.]); self.assertFalse(self.audit()['passed']) + def test_missing_physical_replay_fails(self): + (self.folder / 'coupled_trial_physical.csv').unlink(); self.assertFalse(self.audit()['passed']) + def test_campaign_hash_corruption_suppresses_observations(self): + self.mutate_row('candidates', 'merit', 99.); result = c.inspect_campaign(self.root) + self.assertFalse(result['passed']); self.assertFalse(result['cases']['synthetic']['observations_valid']) + def test_positive_predicted_decrease_required(self): self.assertFalse(c.eligibility(1., 1., 1., 1.)[0]) + def test_beta_remainder_identity(self): + result = c.correction_mechanism([1.], [2.], [1.7], [-.9], [-1.], .5) + self.assertLess(result['decomposition_error_norm'], 1.e-15) + self.assertAlmostEqual(result['own_taylor_remainder_norm'], .15) + + +if __name__ == '__main__': unittest.main() diff --git a/libmeanfield/interface/solver/coupled_trial_correction.cppm b/libmeanfield/interface/solver/coupled_trial_correction.cppm index cd85809..4a8da99 100644 --- a/libmeanfield/interface/solver/coupled_trial_correction.cppm +++ b/libmeanfield/interface/solver/coupled_trial_correction.cppm @@ -12,6 +12,7 @@ module; #include #include #include +#include #include #include @@ -21,6 +22,33 @@ export module mean_field:solver.coupled_trial_correction; export import :solver.linear_backend; export namespace mean_field::solver { + struct FullCoupledStep final { }; + struct BacktrackedCoupledStep final { + int maximumSteps{4}; + double reduction{.5}; + + void Validate() const { + if (maximumSteps <= 0 || !std::isfinite(reduction) || reduction <= 0.0 || reduction >= 1.0 || + std::pow(reduction, maximumSteps - 1) == 0.0) + throw std::invalid_argument("Coupled damping requires positive resolved steps and a reduction strictly between zero and one."); + } + }; + struct EveryRejectedCoupledAttempt final { }; + struct NearAcceptanceCoupledAttempt final { + double maximumExcessRatio{1.0}; + + void Validate() const { + if (!std::isfinite(maximumExcessRatio) || maximumExcessRatio <= 0.0) + throw std::invalid_argument("A coupled prediction gate requires a finite positive excess ratio."); + } + }; + template + concept CoupledCorrectionSearch = std::same_as || + std::same_as; + template + concept CoupledAttemptCriterion = std::same_as || + std::same_as; + struct CoupledTrialCorrectionOptions final { LinearSolveControl control{.relativeTolerance = .03, .absoluteTolerance = 1.e-14, .maximumIterations = 160}; linear::FGMRESOptions backend{.restartLength = 160, .printLevel = -1}; @@ -52,7 +80,16 @@ export namespace mean_field::solver { enum class CoupledTrialAttemptStatus : std::uint8_t { not_attempted, original_accepted, budget_exhausted, kept, linear_failure, norm_limit, - non_finite, inadmissible, not_improved + non_finite, inadmissible, not_improved, prediction_deferred + }; + + enum class CoupledTrialCandidateStatus : std::uint8_t { evaluated, inadmissible, non_finite }; + struct CoupledTrialCandidateReport final { + double beta{1.0}; + CoupledTrialCandidateStatus status{CoupledTrialCandidateStatus::inadmissible}; + double merit{std::numeric_limits::quiet_NaN()}; + bool accepted{false}; + double preparationSeconds{0.0}; }; struct CoupledTrialDiagnostics final { @@ -64,6 +101,16 @@ export namespace mean_field::solver { std::optional currentSolve{}; double originalMerit{std::numeric_limits::quiet_NaN()}; double correctedMerit{std::numeric_limits::quiet_NaN()}; + std::uint64_t deferredAttempts{0}; + int candidateCount{0}; + double selectedBeta{std::numeric_limits::quiet_NaN()}; + double baseMerit{std::numeric_limits::quiet_NaN()}; + double linearPredictionMerit{std::numeric_limits::quiet_NaN()}; + double predictedDecrease{std::numeric_limits::quiet_NaN()}; + double predictionExcessRatio{std::numeric_limits::quiet_NaN()}; + // Search candidates exclude the possible final replay of an earlier + // best beta. Preparations/restorations still account for that work. + std::vector candidates{}; }; } @@ -194,6 +241,7 @@ export namespace mean_field::solver { [[nodiscard]] double GetBaseDirectionNorm() const noexcept { return m_directionNorm; } [[nodiscard]] MPI_Comm GetCommunicator() const noexcept { return m_communicator; } [[nodiscard]] const mfem::Vector &GetRightHandSide() const noexcept { return m_rightHandSide; } + [[nodiscard]] const mfem::Vector &GetBaseResidual() const noexcept { return m_baseResidual; } [[nodiscard]] const mfem::Vector &GetVerifiedAction() const noexcept { return m_verifiedAction; } [[nodiscard]] const mfem::Vector &GetVerifiedDefect() const noexcept { return m_defect; } [[nodiscard]] const auto &GetLastReport() const noexcept { return m_lastReport; } @@ -302,23 +350,24 @@ export namespace mean_field::solver { // execute the same transaction while their typed callbacks choose the metric // and enthalpy policy. This helper never reads or writes accepted-state buffers. export namespace mean_field::solver::detail { - template class PreparedCoupledTrialTransaction final { + template + class PreparedCoupledTrialTransaction final { using Operation = std::remove_cvref_t().normalizedOperator)>; using Inverse = std::remove_cvref_t().normalizedInverse)>; using Clock = std::chrono::steady_clock; public: PreparedCoupledTrialTransaction(State &state, const CoupledTrialCorrectionOptions &options, - int maximumAttemptsPerIteration = 1) + int maximumAttemptsPerIteration = 1, Search search = {}, Attempt attempt = {}) : m_core(*state.normalizedOperator, *state.normalizedInverse, state.Problem().GetCommunicator(), options), - m_maximumAttempts(maximumAttemptsPerIteration), m_originalState(state.normalizedOperator->Width()), - m_originalResidual(state.normalizedOperator->Height()), m_originalPhysical(state.normalizedOperator->Width()), - m_correction(state.normalizedOperator->Width()) { - CoupledRequire(maximumAttemptsPerIteration > 0, m_core.GetCommunicator(), "A coupled trial attempt budget must be positive."); - int minimum = 0, maximum = 0; - if (MPI_Allreduce(&maximumAttemptsPerIteration, &minimum, 1, MPI_INT, MPI_MIN, m_core.GetCommunicator()) != MPI_SUCCESS || - MPI_Allreduce(&maximumAttemptsPerIteration, &maximum, 1, MPI_INT, MPI_MAX, m_core.GetCommunicator()) != MPI_SUCCESS) - throw std::runtime_error("Coupled trial attempt budgets could not be compared."); - if (minimum != maximum) throw std::invalid_argument("Coupled trial attempt budgets must agree across ranks."); + m_maximumAttempts(maximumAttemptsPerIteration), m_search(search), m_attempt(attempt), + m_originalState(state.normalizedOperator->Width()), m_originalResidual(state.normalizedOperator->Height()), + m_originalPhysical(state.normalizedOperator->Width()), m_correction(state.normalizedOperator->Width()), + m_selectedCorrection(state.normalizedOperator->Width()), + m_prediction(std::same_as ? state.normalizedOperator->Height() : 0) { + ValidatePoliciesCollectively(); + m_selectedCorrection = 0.0; + m_diagnostics.candidates.reserve(MaximumSteps()); } void CaptureBase(const mfem::Vector &residual, const mfem::Vector &direction) { m_core.CaptureBase(residual, direction); @@ -327,11 +376,19 @@ export namespace mean_field::solver::detail { void BeginTrial() noexcept { m_diagnostics.status = CoupledTrialAttemptStatus::not_attempted; m_diagnostics.currentSolve.reset(); - m_diagnostics.originalMerit = m_diagnostics.correctedMerit = std::numeric_limits::quiet_NaN(); + m_diagnostics.originalMerit = m_diagnostics.correctedMerit = m_diagnostics.selectedBeta = + m_diagnostics.baseMerit = m_diagnostics.linearPredictionMerit = m_diagnostics.predictedDecrease = + m_diagnostics.predictionExcessRatio = std::numeric_limits::quiet_NaN(); + m_diagnostics.candidateCount = 0; + m_diagnostics.candidates.clear(); + m_selectedCorrection = 0.0; } [[nodiscard]] const CoupledTrialDiagnostics &GetDiagnostics() const noexcept { return m_diagnostics; } [[nodiscard]] const auto &GetCore() const noexcept { return m_core; } + // The report verifies this unscaled vector. Damping changes the trial + // map, not the vector to which its original linear certificate applies. [[nodiscard]] const mfem::Vector &GetCorrection() const noexcept { return m_correction; } + [[nodiscard]] const mfem::Vector &GetSelectedCorrection() const noexcept { return m_selectedCorrection; } template void TryImprove(State &state, double alpha, Preparation &preparation, Metric &metric, @@ -339,8 +396,7 @@ export namespace mean_field::solver::detail { const auto communicator = m_core.GetCommunicator(); CoupledRequire(preparation.has_value() && state.normalizedOperator->IsPrepared() && CoupledFinite(state.trialNormalizedState) && CoupledFinite(state.trialNormalizedResidual) && - std::isfinite(metric.merit), communicator, - "Coupled correction requires a successfully prepared original trial."); + ValidMetric(metric), communicator, "Coupled correction requires a successfully prepared original trial."); m_diagnostics.originalMerit = metric.merit; if (CoupledAll(accepts(metric), communicator)) { m_diagnostics.status = CoupledTrialAttemptStatus::original_accepted; @@ -350,6 +406,27 @@ export namespace mean_field::solver::detail { m_diagnostics.status = CoupledTrialAttemptStatus::budget_exhausted; return; } + if constexpr (std::same_as) { + // These are evaluations of saved residual vectors, not new + // physical candidates. No preparation, inverse refresh or + // attempt budget is spent on the prediction gate. + m_core.BuildLinearPrediction(alpha, m_prediction); + const auto baseMetric = evaluate(m_core.GetBaseResidual()); + const auto predictionMetric = evaluate(m_prediction); + m_diagnostics.baseMerit = baseMetric.merit; + m_diagnostics.linearPredictionMerit = predictionMetric.merit; + m_diagnostics.predictedDecrease = baseMetric.merit - predictionMetric.merit; + const double decrease = m_diagnostics.predictedDecrease; + const double ratio = decrease > 0.0 ? (metric.merit - baseMetric.merit) / decrease : + std::numeric_limits::quiet_NaN(); + m_diagnostics.predictionExcessRatio = ratio; + if (!CoupledAll(ValidMetric(baseMetric) && ValidMetric(predictionMetric) && std::isfinite(decrease) && + decrease > 0.0 && std::isfinite(ratio) && ratio <= m_attempt.maximumExcessRatio, communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::prediction_deferred; + ++m_diagnostics.deferredAttempts; + return; + } + } m_originalState = state.trialNormalizedState; m_originalResidual = state.trialNormalizedResidual; m_originalPhysical = state.candidatePhysicalState; @@ -368,14 +445,14 @@ export namespace mean_field::solver::detail { case CoupledTrialCorrectionStatus::norm_limit: m_diagnostics.status = CoupledTrialAttemptStatus::norm_limit; break; default: m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; break; } - // No candidate mutation has occurred; its prepared operator and - // freshly refreshed inverse still describe the original trial. + // No candidate mutation has occurred. In particular, a small + // beta cannot rescue an unverified or over-cap full correction. return; } const Metric originalMetric = metric; - bool restorationStarted = false; - auto restore = [&]() { - restorationStarted = true; + bool originalRestorationStarted = false; + auto restoreOriginal = [&]() { + originalRestorationStarted = true; state.trialNormalizedState = m_originalState; state.candidatePhysicalState = m_originalPhysical; auto restored = Prepare(state); @@ -387,56 +464,155 @@ export namespace mean_field::solver::detail { metric = evaluate(state.trialNormalizedResidual); ++m_diagnostics.fallbacks; }; - try { - state.trialNormalizedState += m_correction; + const auto formCandidate = [&](double beta) { + state.trialNormalizedState = m_originalState; + state.trialNormalizedState.Add(beta, m_correction); state.normalizedOperator->DenormalizeState(state.trialNormalizedState, state.candidatePhysicalState); - if (!CoupledAll(CoupledFinite(state.trialNormalizedState) && CoupledFinite(state.candidatePhysicalState), communicator)) { - m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; - restore(); - return; - } - if (!CoupledAll(checkState(state.trialNormalizedState), communicator)) { - m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; - restore(); - return; - } - auto corrected = Prepare(state); - if (!CoupledAll(corrected.has_value(), communicator)) { - m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; - restore(); - return; - } - if (!CoupledAll(CoupledFinite(state.trialNormalizedResidual), communicator)) { - m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; - restore(); - return; - } - const auto correctedMetric = evaluate(state.trialNormalizedResidual); - m_diagnostics.correctedMerit = correctedMetric.merit; - if (!CoupledAll(better(correctedMetric, originalMetric), communicator)) { + }; + try { + std::optional bestMetric; + double bestBeta = 0.0; + int bestStep = -1, lastStep = -1; + double beta = 1.0; + for (int step = 0; step < MaximumSteps(); ++step, beta *= Reduction()) { + lastStep = step; + formCandidate(beta); + m_diagnostics.candidates.push_back({.beta = beta}); + ++m_diagnostics.candidateCount; + auto &candidate = m_diagnostics.candidates.back(); + if (!CoupledAll(CoupledFinite(state.trialNormalizedState) && CoupledFinite(state.candidatePhysicalState), communicator)) { + candidate.status = CoupledTrialCandidateStatus::non_finite; + m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; + continue; + } + if (!CoupledAll(checkState(state.trialNormalizedState), communicator)) { + candidate.status = CoupledTrialCandidateStatus::inadmissible; + m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; + continue; + } + const double preparationStart = m_diagnostics.preparationSeconds; + auto corrected = Prepare(state); + candidate.preparationSeconds = m_diagnostics.preparationSeconds - preparationStart; + if (!CoupledAll(corrected.has_value(), communicator)) { + candidate.status = CoupledTrialCandidateStatus::inadmissible; + m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; + continue; + } + if (!CoupledAll(CoupledFinite(state.trialNormalizedResidual), communicator)) { + candidate.status = CoupledTrialCandidateStatus::non_finite; + m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; + continue; + } + const auto correctedMetric = evaluate(state.trialNormalizedResidual); + m_diagnostics.correctedMerit = correctedMetric.merit; + if (!CoupledAll(ValidMetric(correctedMetric), communicator)) { + candidate.status = CoupledTrialCandidateStatus::non_finite; + m_diagnostics.status = CoupledTrialAttemptStatus::non_finite; + continue; + } + preparation = std::move(corrected); + candidate.status = CoupledTrialCandidateStatus::evaluated; + candidate.merit = correctedMetric.merit; + const bool improves = CoupledAll(better(correctedMetric, originalMetric), communicator); + candidate.accepted = improves && CoupledAll(accepts(correctedMetric), communicator); + if (improves && (!bestMetric || candidate.accepted || CoupledAll(better(correctedMetric, *bestMetric), communicator))) { + bestMetric = correctedMetric; + bestBeta = beta; + bestStep = step; + } m_diagnostics.status = CoupledTrialAttemptStatus::not_improved; - restore(); + // An improved but unacceptable full correction must not + // suppress the remaining, cheaper physical preparations. + if (candidate.accepted) break; + } + if (!bestMetric) { + restoreOriginal(); return; } - preparation = std::move(corrected); - metric = correctedMetric; + if (bestStep != lastStep) { + // The current prepared operator belongs to the last beta, + // which can even be inadmissible. Replaying the selected + // state restores its residual, geometry and generation. + formCandidate(bestBeta); + if (!CoupledAll(checkState(state.trialNormalizedState), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; + restoreOriginal(); + return; + } + auto restored = Prepare(state); + ++m_diagnostics.restorations; + if (!CoupledAll(restored.has_value() && CoupledFinite(state.trialNormalizedResidual), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::inadmissible; + restoreOriginal(); + return; + } + const auto restoredMetric = evaluate(state.trialNormalizedResidual); + if (!CoupledAll(ValidMetric(restoredMetric), communicator) || + !CoupledAll(better(restoredMetric, originalMetric), communicator)) { + m_diagnostics.status = CoupledTrialAttemptStatus::not_improved; + restoreOriginal(); + return; + } + preparation = std::move(restored); + bestMetric = restoredMetric; + Refresh(state); + } + metric = *bestMetric; + m_selectedCorrection = m_correction; + m_selectedCorrection *= bestBeta; + m_diagnostics.selectedBeta = bestBeta; + m_diagnostics.correctedMerit = metric.merit; m_diagnostics.status = CoupledTrialAttemptStatus::kept; ++m_diagnostics.kept; - // This is a nonlinear candidate, not an affine gravity completion. - // Its fresh full residual is judged by the caller's normal metric. - // As for an ordinary prepared trial, the caller refreshes the - // inverse only if it commits or needs another auxiliary solve. + // A current selected trial retains the ordinary lifecycle: + // its inverse is refreshed only on commit or another solve. } catch (...) { - // Numerical rejection is handled above. For collectively - // propagated exceptions, restore the original prepared trial - // before rethrowing. A failed restoration itself propagates to - // the owning Newton context, which restores the accepted base. - if (!restorationStarted) restore(); + // Only collectively propagated callback exceptions can be + // rolled back collectively. Failure to restore propagates to + // the owning Newton context's accepted-base restoration. + if (!originalRestorationStarted) restoreOriginal(); throw; } } private: + template static bool ValidMetric(const Metric &metric) { + bool valid = std::isfinite(metric.merit) && metric.merit >= 0.0; + if constexpr (requires { metric.residualNorm; }) + valid = valid && std::isfinite(metric.residualNorm) && metric.residualNorm >= 0.0; + return valid; + } + [[nodiscard]] int MaximumSteps() const noexcept { + if constexpr (std::same_as) return m_search.maximumSteps; + else return 1; + } + [[nodiscard]] double Reduction() const noexcept { + if constexpr (std::same_as) return m_search.reduction; + else return 1.0; + } + void ValidatePoliciesCollectively() const { + const auto communicator = m_core.GetCommunicator(); + bool valid = m_maximumAttempts > 0; + try { + if constexpr (std::same_as) m_search.Validate(); + if constexpr (std::same_as) m_attempt.Validate(); + } catch (const std::invalid_argument &) { valid = false; } + CoupledRequire(valid, communicator, "Coupled trial search or attempt controls are invalid on at least one rank."); + double excessRatio = 0.0; + if constexpr (std::same_as) excessRatio = m_attempt.maximumExcessRatio; + const std::array integers{m_maximumAttempts, std::same_as, + MaximumSteps(), std::same_as}; + const std::array reals{Reduction(), excessRatio}; + std::array minimumIntegers{}, maximumIntegers{}; + std::array minimumReals{}, maximumReals{}; + if (MPI_Allreduce(integers.data(), minimumIntegers.data(), 4, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(integers.data(), maximumIntegers.data(), 4, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), minimumReals.data(), 2, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), maximumReals.data(), 2, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) + throw std::runtime_error("Coupled trial search and attempt controls could not be compared."); + if (minimumIntegers != maximumIntegers || minimumReals != maximumReals) + throw std::invalid_argument("Coupled trial search and attempt controls must agree across ranks."); + } void Refresh(State &state) { const auto start = Clock::now(); state.physicalInverse->Refresh(); @@ -451,7 +627,9 @@ export namespace mean_field::solver::detail { } PreparedCoupledTrialCorrection m_core; int m_maximumAttempts, m_attemptsThisIteration{0}; - mfem::Vector m_originalState, m_originalResidual, m_originalPhysical, m_correction; + [[no_unique_address]] Search m_search; + [[no_unique_address]] Attempt m_attempt; + mfem::Vector m_originalState, m_originalResidual, m_originalPhysical, m_correction, m_selectedCorrection, m_prediction; CoupledTrialDiagnostics m_diagnostics{}; }; } diff --git a/libmeanfield/interface/solver/newton_numerics.cppm b/libmeanfield/interface/solver/newton_numerics.cppm index 618ccc3..0f44ad1 100644 --- a/libmeanfield/interface/solver/newton_numerics.cppm +++ b/libmeanfield/interface/solver/newton_numerics.cppm @@ -53,6 +53,25 @@ export namespace mean_field::solver::nonlinear { } }; + struct BacktrackedCoupledCurvatureCorrection final { + CoupledTrialCorrectionOptions correction{}; + int maximumAttemptsPerIteration{1}; + BacktrackedCoupledStep search{}; + NearAcceptanceCoupledAttempt attempt{}; + + void Validate() const { + correction.Validate(); + search.Validate(); + attempt.Validate(); + if (maximumAttemptsPerIteration <= 0) + throw std::invalid_argument("A coupled trial policy requires a positive attempt budget."); + } + }; + + template + concept CoupledTrialNumerics = std::same_as || + std::same_as; + template concept EnthalpyNumerics = std::same_as || std::same_as; @@ -64,7 +83,7 @@ export namespace mean_field::solver::nonlinear { template concept TrialPathNumerics = std::same_as || - std::same_as; + CoupledTrialNumerics; // Types select algorithms; tolerances and the mesh-dependent work sizes // remain runtime data. Adding a strategy cannot silently change a model's diff --git a/libmeanfield/interface/solver/stellar_newton_strategy.cppm b/libmeanfield/interface/solver/stellar_newton_strategy.cppm index 07170b5..a33e8f2 100644 --- a/libmeanfield/interface/solver/stellar_newton_strategy.cppm +++ b/libmeanfield/interface/solver/stellar_newton_strategy.cppm @@ -44,8 +44,8 @@ namespace mean_field::solver::detail { using Enthalpy = typename Recipe::EnthalpyType; using Gravity = typename Recipe::GravityType; using Trial = typename Recipe::TrialType; - std::array integers{}; - std::array reals{}; + std::array integers{}; + std::array reals{}; integers[0] = std::same_as; integers[1] = std::same_as ? 1 : std::same_as ? 2 : 0; @@ -61,8 +61,9 @@ namespace mean_field::solver::detail { reals[0] = options.control.relativeTolerance; reals[1] = options.control.absoluteTolerance; } - integers[8] = std::same_as; - if constexpr (std::same_as) { + integers[8] = std::same_as ? 1 : + std::same_as ? 2 : 0; + if constexpr (nonlinear::CoupledTrialNumerics) { const auto &options = recipe.trial.correction; integers[9] = options.control.maximumIterations; integers[10] = options.backend.restartLength; @@ -72,13 +73,18 @@ namespace mean_field::solver::detail { reals[3] = options.control.absoluteTolerance; reals[4] = options.maximumCorrectionRatio; } + if constexpr (std::same_as) { + integers[13] = recipe.trial.search.maximumSteps; + reals[5] = recipe.trial.search.reduction; + reals[6] = recipe.trial.attempt.maximumExcessRatio; + } integers[7] = std::all_of(reals.begin(), reals.end(), [](double value) { return std::isfinite(value); }); - std::array minimumIntegers{}, maximumIntegers{}; - std::array minimumReals{}, maximumReals{}; - if (MPI_Allreduce(integers.data(), minimumIntegers.data(), 13, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS || - MPI_Allreduce(integers.data(), maximumIntegers.data(), 13, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS || - MPI_Allreduce(reals.data(), minimumReals.data(), 5, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS || - MPI_Allreduce(reals.data(), maximumReals.data(), 5, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) { + std::array minimumIntegers{}, maximumIntegers{}; + std::array minimumReals{}, maximumReals{}; + if (MPI_Allreduce(integers.data(), minimumIntegers.data(), 14, MPI_INT, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(integers.data(), maximumIntegers.data(), 14, MPI_INT, MPI_MAX, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), minimumReals.data(), 7, MPI_DOUBLE, MPI_MIN, communicator) != MPI_SUCCESS || + MPI_Allreduce(reals.data(), maximumReals.data(), 7, MPI_DOUBLE, MPI_MAX, communicator) != MPI_SUCCESS) { throw std::runtime_error("Newton could not compare numerical strategies across ranks."); } if (minimumIntegers != maximumIntegers || minimumReals != maximumReals || @@ -86,7 +92,7 @@ namespace mean_field::solver::detail { throw std::invalid_argument("Newton numerical strategies and controls must agree on every rank."); } if constexpr (!std::same_as) { recipe.gravity.gravity.Validate(); } - if constexpr (std::same_as) { recipe.trial.Validate(); } + if constexpr (nonlinear::CoupledTrialNumerics) { recipe.trial.Validate(); } } template class PreparedGravityStrategy; @@ -431,6 +437,11 @@ namespace mean_field::solver::detail { : transaction(state, policy.correction, policy.maximumAttemptsPerIteration) { } PreparedCoupledTrialTransaction transaction; }; + template struct PreparedTrialPath { + PreparedTrialPath(State &state, const nonlinear::BacktrackedCoupledCurvatureCorrection &policy) + : transaction(state, policy.correction, policy.maximumAttemptsPerIteration, policy.search, policy.attempt) { } + PreparedCoupledTrialTransaction transaction; + }; template class PreparedNewtonStrategy final : public NewtonStrategyStorage { @@ -439,7 +450,7 @@ namespace mean_field::solver::detail { using Gravity = typename Recipe::GravityType; using Trial = typename Recipe::TrialType; static constexpr bool altersDirection = std::same_as; - static constexpr bool correctsTrial = std::same_as; + static constexpr bool correctsTrial = nonlinear::CoupledTrialNumerics; PreparedNewtonStrategy(State &state, const Recipe &recipe) : m_enthalpy(state, recipe.enthalpy), m_gravity(state, recipe.gravity), m_trial(state, recipe.trial) { } diff --git a/tests/solver/coupled_trial_correction.cpp b/tests/solver/coupled_trial_correction.cpp index 4e0f411..1ad4676 100644 --- a/tests/solver/coupled_trial_correction.cpp +++ b/tests/solver/coupled_trial_correction.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -162,6 +163,12 @@ namespace coupled_trial_correction_test { static_assert(CoupledTrialInverse); static_assert(!CoupledTrialOperation); static_assert(!CoupledTrialInverse); + static_assert(CoupledCorrectionSearch); + static_assert(CoupledCorrectionSearch); + static_assert(!CoupledCorrectionSearch); + static_assert(CoupledAttemptCriterion); + static_assert(CoupledAttemptCriterion); + static_assert(!CoupledAttemptCriterion); static_assert(LinearBackendRuntimeAvailableFor); using ScriptedCorrection = PreparedCoupledTrialCorrection; } // namespace coupled_trial_correction_test @@ -435,6 +442,11 @@ namespace coupled_trial_correction_test { int preparations{0}; bool quarticBarrier{false}, rejectCorrectedPreparation{false}, nonfiniteCorrectedResidual{false}; bool throwCorrectedPreparation{false}; + bool cubicSearchModel{false}; + bool rejectPositivePreparation{false}; + double barrierMultiplier{1.0}; + int throwAtPreparation{-1}; + std::vector preparedCoordinates; ToyState() { trialNormalizedState = 0.0; @@ -449,23 +461,37 @@ namespace coupled_trial_correction_test { [[nodiscard]] std::expected PrepareTrial() { ++preparations; candidatePhysicalState = trialNormalizedState; + preparedCoordinates.push_back(trialNormalizedState(0)); + if (preparations == throwAtPreparation) { + operation.prepared = false; + throw std::runtime_error("Synthetic delayed preparation failure."); + } if (throwCorrectedPreparation && trialNormalizedState(0) < -2.0) { operation.prepared = false; throw std::runtime_error("Synthetic corrected preparation failure."); } - if (rejectCorrectedPreparation && trialNormalizedState(0) < -2.0) { + if ((rejectCorrectedPreparation && trialNormalizedState(0) < -2.0) || + (rejectPositivePreparation && trialNormalizedState(0) > 0.0)) { operation.prepared = false; return std::unexpected(7); } trialNormalizedResidual = coupledResidual(trialNormalizedState); auto matrix = coupledJacobian(trialNormalizedState); + if (cubicSearchModel) { + // At alpha=.5, c=(1.5,0,0) overshoots badly at full scale. + // Its quarter step improves the ordinary residual merit. + const double x = trialNormalizedState(0); + trialNormalizedResidual(0) = 1.0 + x + 10.0*x*x + 11.0*x*x*x; + matrix = identity(3); + matrix(0, 0) = 1.0 + 20.0*x + 33.0*x*x; + } if (quarticBarrier) { // The extra term and its derivative vanish at the base and at // alpha=.5. It makes the solved correction harmful without // supplying the transaction an inconsistent trial Jacobian. const double x = trialNormalizedState(0); - trialNormalizedResidual(0) += x*x*(x+.5)*(x+.5); - matrix(0, 0) += 2.0*x*(x+.5)*(2.0*x+.5); + trialNormalizedResidual(0) += barrierMultiplier*x*x*(x+.5)*(x+.5); + matrix(0, 0) += barrierMultiplier*2.0*x*(x+.5)*(2.0*x+.5); } operation.Prepare(matrix); if (nonfiniteCorrectedResidual && trialNormalizedState(0) < -2.0) @@ -789,3 +815,365 @@ TEST_CASE("Exceptions after coupled mutation restore the original prepared candi REQUIRE(preparation.has_value()); checkRestoredTrial(state, originalState, originalResidual, originalAction); } + +TEST_CASE("Damping one verified coupled correction accepts a nonlinear step that its full scale misses", + "[solver][coupled_trial_correction][backtracked]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + ToyState state; + state.cubicSearchModel = true; + SECTION("Finite full and half corrections overshoot the nonlinear residual") { } + SECTION("Rejected full and half preparations are followed by a valid smaller correction") { + state.rejectPositivePreparation = true; + } + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + auto metric = toyMetric(state.trialNormalizedResidual); + REQUIRE(metric.merit == Approx(1.5703125)); + const auto accepts = [](const auto &candidate) { return candidate.merit <= 1.5*(1.0-1.e-4*.5); }; + REQUIRE_FALSE(accepts(metric)); + const int initialActions = state.operation.applications, initialInverses = state.inverse.applications; + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, accepts, toyImproves, admissibleToy); + const auto &diagnostics = transaction.GetDiagnostics(); + REQUIRE(diagnostics.currentSolve.has_value()); + REQUIRE(diagnostics.currentSolve->Admissible()); + REQUIRE(diagnostics.currentSolve->verifiedGate); + CHECK(diagnostics.status == CoupledTrialAttemptStatus::kept); + CHECK(diagnostics.selectedBeta == .25); + CHECK(diagnostics.attempts == 1); + CHECK(diagnostics.kept == 1); + CHECK(diagnostics.fallbacks == 0); + CHECK(diagnostics.candidateCount == 3); + REQUIRE(diagnostics.candidates.size() == 3); + CHECK(diagnostics.candidates[0].beta == 1.0); + CHECK_FALSE(diagnostics.candidates[0].accepted); + if (state.rejectPositivePreparation) { + CHECK(diagnostics.candidates[0].status == CoupledTrialCandidateStatus::inadmissible); + CHECK(diagnostics.candidates[1].status == CoupledTrialCandidateStatus::inadmissible); + } else CHECK(diagnostics.candidates[0].merit > diagnostics.originalMerit); + CHECK(diagnostics.candidates[1].beta == .5); + CHECK_FALSE(diagnostics.candidates[1].accepted); + CHECK(diagnostics.candidates[2].beta == .25); + CHECK(diagnostics.candidates[2].accepted); + CHECK(diagnostics.preparations == 3); + CHECK(diagnostics.restorations == 0); + CHECK(transaction.GetCorrection()(0) == Approx(1.5)); + CHECK(transaction.GetSelectedCorrection()(0) == Approx(.375)); + CHECK(state.trialNormalizedState(0) == Approx(-.125)); + CHECK(state.trialNormalizedResidual(0) == Approx(1.009765625)); + CHECK(accepts(metric)); + // Candidate preparations change J but must not launch another solve or + // verify a scaled correction against a different candidate Jacobian. + CHECK(state.operation.applications-initialActions == diagnostics.currentSolve->linearSolve.operatorApplications+1); + CHECK(state.inverse.applications-initialInverses == diagnostics.currentSolve->linearSolve.inversePreconditionerApplications); + CHECK(diagnostics.iterations == static_cast(diagnostics.currentSolve->linearSolve.iterations)); + CHECK(state.operation.IsPrepared()); + CHECK_FALSE(state.inverse.IsCurrent()); + checkOuterData(state); + transaction.BeginTrial(); + CHECK(std::isnan(transaction.GetDiagnostics().selectedBeta)); + CHECK(transaction.GetDiagnostics().candidateCount == 0); + CHECK(transaction.GetDiagnostics().candidates.empty()); + CHECK_FALSE(transaction.GetDiagnostics().currentSolve.has_value()); + CHECK(transaction.GetSelectedCorrection().Norml2() == 0.0); +} + +TEST_CASE("A bounded coupled search freshly restores its best earlier beta when no beta passes Armijo", + "[solver][coupled_trial_correction][backtracked][rollback]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + ToyState state; + state.cubicSearchModel = true; + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + auto metric = toyMetric(state.trialNormalizedResidual); + // This stringent but valid Armijo parameter excludes the quarter step; + // it remains the best full-merit improvement among all four proposals. + const auto accepts = [](const auto &candidate) { return candidate.merit <= 1.5*(1.0-.99*.5); }; + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, accepts, toyImproves, admissibleToy); + const auto &diagnostics = transaction.GetDiagnostics(); + CHECK(diagnostics.status == CoupledTrialAttemptStatus::kept); + CHECK(diagnostics.selectedBeta == .25); + CHECK(diagnostics.candidateCount == 4); + REQUIRE(diagnostics.candidates.size() == 4); + for (const auto &candidate : diagnostics.candidates) CHECK_FALSE(candidate.accepted); + CHECK(diagnostics.candidates[2].merit < diagnostics.candidates[3].merit); + CHECK_FALSE(accepts(metric)); + CHECK(metric.merit < diagnostics.originalMerit); + CHECK(diagnostics.preparations == 5); + CHECK(diagnostics.restorations == 1); + CHECK(diagnostics.fallbacks == 0); + REQUIRE(state.preparedCoordinates.size() == 6); + CHECK(state.preparedCoordinates[4] == Approx(-.3125)); + CHECK(state.preparedCoordinates[5] == Approx(-.125)); + CHECK(state.trialNormalizedState(0) == Approx(-.125)); + CHECK(state.candidatePhysicalState(0) == Approx(-.125)); + CHECK(state.trialNormalizedResidual(0) == Approx(1.009765625)); + CHECK(state.operation.Matrix()(0, 0) == Approx(1.0+20.0*(-.125)+33.0*.125*.125)); + CHECK(state.operation.IsPrepared()); + checkOuterData(state); +} + +TEST_CASE("Every rejected damping proposal falls back to the original prepared trial", + "[solver][coupled_trial_correction][backtracked][rollback]") { + using namespace coupled_trial_correction_test; + ToyState state; + bool rejectAdmissibility = false, invalidMetric = false; + auto expected = CoupledTrialAttemptStatus::not_improved; + SECTION("Every nonlinear correction has worse full merit") { + state.quarticBarrier = true; + state.barrierMultiplier = 100.0; + } + SECTION("Every beta fails admissibility before preparation") { + rejectAdmissibility = true; + expected = CoupledTrialAttemptStatus::inadmissible; + } + SECTION("A finite merit cannot hide a nonfinite metric norm") { + invalidMetric = true; + expected = CoupledTrialAttemptStatus::non_finite; + } + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + const mfem::Vector originalState(state.trialNormalizedState), originalResidual(state.trialNormalizedResidual); + mfem::Vector originalAction(3); + state.operation.Mult(state.normalizedCorrection, originalAction); + auto metric = toyMetric(originalResidual); + const double originalMerit = metric.merit; + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, + [&](const mfem::Vector &residual) { + if (invalidMetric && state.trialNormalizedState(0) != originalState(0)) + return nonlinear::MetricEvaluation{.residualNorm = std::numeric_limits::quiet_NaN(), .merit = .1}; + return toyMetric(residual); + }, toyAccepts, toyImproves, [rejectAdmissibility](const mfem::Vector &) { return !rejectAdmissibility; }); + const auto &diagnostics = transaction.GetDiagnostics(); + CHECK(diagnostics.status == expected); + REQUIRE(diagnostics.currentSolve.has_value()); + CHECK(diagnostics.currentSolve->Admissible()); + CHECK(diagnostics.attempts == 1); + CHECK(diagnostics.candidateCount == 4); + CHECK(diagnostics.kept == 0); + CHECK(diagnostics.fallbacks == 1); + CHECK(diagnostics.restorations == 1); + CHECK(diagnostics.preparations == (rejectAdmissibility ? 1 : 5)); + CHECK(std::isnan(diagnostics.selectedBeta)); + CHECK(transaction.GetSelectedCorrection().Norml2() == 0.0); + CHECK(metric.merit == originalMerit); + REQUIRE(preparation.has_value()); + checkRestoredTrial(state, originalState, originalResidual, originalAction); +} + +TEST_CASE("Exceptions during later damping or best-candidate replay restore the original trial", + "[solver][coupled_trial_correction][backtracked][rollback][exception]") { + using namespace coupled_trial_correction_test; + ToyState state; + state.cubicSearchModel = true; + enum class FailureSite { later_admissibility, later_metric, replay_preparation, replay_metric }; + FailureSite site = FailureSite::later_admissibility; + SECTION("The second admissibility callback throws") { site = FailureSite::later_admissibility; } + SECTION("The second candidate metric throws") { site = FailureSite::later_metric; } + SECTION("Repreparing the earlier best invalidates J and throws") { site = FailureSite::replay_preparation; } + SECTION("The freshly replayed best metric throws") { site = FailureSite::replay_metric; } + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + const mfem::Vector originalState(state.trialNormalizedState), originalResidual(state.trialNormalizedResidual); + mfem::Vector originalAction(3); + state.operation.Mult(state.normalizedCorrection, originalAction); + auto metric = toyMetric(originalResidual); + const double originalMerit = metric.merit; + if (site == FailureSite::replay_preparation) state.throwAtPreparation = state.preparations+5; + int metricCalls = 0, admissibilityCalls = 0; + transaction.BeginTrial(); + const auto attempt = [&] { + transaction.TryImprove(state, .5, preparation, metric, + [&](const mfem::Vector &residual) { + ++metricCalls; + if ((site == FailureSite::later_metric && metricCalls == 2) || + (site == FailureSite::replay_metric && metricCalls == 5)) + throw std::runtime_error("Synthetic damping metric failure."); + return toyMetric(residual); + }, [](const auto &candidate) { return candidate.merit <= 1.5*(1.0-.99*.5); }, toyImproves, + [&](const mfem::Vector &) { + if (++admissibilityCalls == 2 && site == FailureSite::later_admissibility) + throw std::runtime_error("Synthetic later admissibility failure."); + return true; + }); + }; + CHECK_THROWS_AS(attempt(), std::runtime_error); + CHECK(transaction.GetDiagnostics().kept == 0); + CHECK(std::isnan(transaction.GetDiagnostics().selectedBeta)); + CHECK(transaction.GetSelectedCorrection().Norml2() == 0.0); + CHECK(metric.merit == originalMerit); + REQUIRE(preparation.has_value()); + checkRestoredTrial(state, originalState, originalResidual, originalAction); +} + +TEST_CASE("Near-acceptance selection defers work without spending the captured direction budget", + "[solver][coupled_trial_correction][backtracked][budget]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + ToyState state; + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + auto metric = toyMetric(state.trialNormalizedResidual); + const int actions = state.operation.applications, inverses = state.inverse.applications, preparations = state.preparations; + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, + [](const auto &candidate) { return candidate.merit <= 1.5*(1.0-1.e-4*.5); }, toyImproves, admissibleToy); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::prediction_deferred); + CHECK(transaction.GetDiagnostics().baseMerit == Approx(1.5)); + CHECK(transaction.GetDiagnostics().linearPredictionMerit == Approx(.375)); + CHECK(transaction.GetDiagnostics().predictedDecrease == Approx(1.125)); + CHECK(transaction.GetDiagnostics().predictionExcessRatio == Approx(26.0/9.0)); + CHECK(transaction.GetDiagnostics().deferredAttempts == 1); + CHECK(transaction.GetDiagnostics().attempts == 0); + CHECK_FALSE(transaction.GetDiagnostics().currentSolve.has_value()); + CHECK(state.operation.applications == actions); + CHECK(state.inverse.applications == inverses); + CHECK(state.preparations == preparations); + + preparation = state.PrepareAlpha(.25); + metric = toyMetric(state.trialNormalizedResidual); + REQUIRE_FALSE(toyAccepts(metric)); + transaction.BeginTrial(); + CHECK(std::isnan(transaction.GetDiagnostics().predictionExcessRatio)); + transaction.TryImprove(state, .25, preparation, metric, toyMetric, toyAccepts, toyImproves, admissibleToy); + CHECK(transaction.GetDiagnostics().predictionExcessRatio == Approx(1.0/84.0)); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::kept); + CHECK(transaction.GetDiagnostics().selectedBeta == 1.0); + CHECK(transaction.GetDiagnostics().attempts == 1); + CHECK(transaction.GetDiagnostics().deferredAttempts == 1); + CHECK(metric.merit == Approx(.84375)); + CHECK(toyAccepts(metric)); + + preparation = state.PrepareAlpha(.25); + metric = toyMetric(state.trialNormalizedResidual); + transaction.BeginTrial(); + transaction.TryImprove(state, .25, preparation, metric, toyMetric, toyAccepts, toyImproves, admissibleToy); + CHECK(transaction.GetDiagnostics().status == CoupledTrialAttemptStatus::budget_exhausted); + CHECK(transaction.GetDiagnostics().attempts == 1); + CHECK_FALSE(transaction.GetDiagnostics().currentSolve.has_value()); + CHECK(transaction.GetDiagnostics().candidateCount == 0); + CHECK(transaction.GetDiagnostics().candidates.empty()); + CHECK(std::isnan(transaction.GetDiagnostics().selectedBeta)); + CHECK(transaction.GetSelectedCorrection().Norml2() == 0.0); + checkOuterData(state); +} + +TEST_CASE("Near-acceptance selection skips acceptable trials and requires a predicted decrease", + "[solver][coupled_trial_correction][backtracked][budget]") { + using namespace coupled_trial_correction_test; + ToyState state; + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + mfem::Vector direction(state.normalizedCorrection); + bool accepted = false; + SECTION("An already acceptable original never starts auxiliary work") { accepted = true; } + SECTION("A zero predicted decrease cannot define a near-acceptance allowance") { direction = 0.0; } + transaction.CaptureBase(state.acceptedResidual, direction); + auto preparation = state.PrepareAlpha(accepted ? .05 : 0.0); + auto metric = toyMetric(state.trialNormalizedResidual); + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, toyMetric, + [](const auto &candidate) { return candidate.merit <= 1.5*(1.0-1.e-4*.5); }, toyImproves, admissibleToy); + CHECK(transaction.GetDiagnostics().status == (accepted ? CoupledTrialAttemptStatus::original_accepted : + CoupledTrialAttemptStatus::prediction_deferred)); + CHECK(transaction.GetDiagnostics().attempts == 0); + CHECK_FALSE(transaction.GetDiagnostics().currentSolve.has_value()); + CHECK(transaction.GetDiagnostics().candidateCount == 0); + CHECK(transaction.GetSelectedCorrection().Norml2() == 0.0); + if (!accepted) CHECK(transaction.GetDiagnostics().predictedDecrease == 0.0); + checkOuterData(state); +} + +TEST_CASE("Damped coupled candidates make collective admissibility and metric decisions", + "[solver][coupled_trial_correction][backtracked][mpi]") { + using namespace coupled_trial_correction_test; + using Catch::Approx; + int rank = 0; + REQUIRE(MPI_Comm_rank(MPI_COMM_WORLD, &rank) == MPI_SUCCESS); + ToyState state; + state.problem.communicator = MPI_COMM_WORLD; + state.cubicSearchModel = true; + bool invalidMetric = false, strictAcceptance = false; + SECTION("One rank rejects the first two admissibility checks") { } + SECTION("One rank reports a nonfinite quarter-step metric") { invalidMetric = true; } + SECTION("One rank requires the full bounded search before keeping its best") { strictAcceptance = true; } + solver::detail::PreparedCoupledTrialTransaction transaction(state, exactOptions()); + transaction.CaptureBase(state.acceptedResidual, state.normalizedCorrection); + auto preparation = state.PrepareAlpha(.5); + auto metric = toyMetric(state.trialNormalizedResidual); + transaction.BeginTrial(); + transaction.TryImprove(state, .5, preparation, metric, + [&](const mfem::Vector &residual) { + if (invalidMetric && rank == 0 && std::abs(state.trialNormalizedState(0)+.125) < 1.e-12) + return nonlinear::MetricEvaluation{.residualNorm = std::numeric_limits::quiet_NaN(), .merit = .1}; + return toyMetric(residual); + }, [&](const auto &candidate) { + return candidate.merit <= 1.5*(1.0-(strictAcceptance && rank == 0 ? .99 : 1.e-4)*.5); + }, toyImproves, [&](const mfem::Vector &candidate) { + return invalidMetric || strictAcceptance || rank != 0 || candidate(0) <= 0.0; + }); + const auto &diagnostics = transaction.GetDiagnostics(); + REQUIRE(diagnostics.currentSolve.has_value()); + REQUIRE(diagnostics.currentSolve->Admissible()); + CHECK(diagnostics.status == CoupledTrialAttemptStatus::kept); + CHECK(diagnostics.attempts == 1); + CHECK(diagnostics.selectedBeta == (invalidMetric ? .125 : .25)); + CHECK(diagnostics.candidateCount == (invalidMetric || strictAcceptance ? 4 : 3)); + REQUIRE(diagnostics.candidates.size() == diagnostics.candidateCount); + if (invalidMetric) CHECK(diagnostics.candidates[2].status == CoupledTrialCandidateStatus::non_finite); + else if (strictAcceptance) CHECK_FALSE(diagnostics.candidates[2].accepted); + else { + CHECK(diagnostics.candidates[0].status == CoupledTrialCandidateStatus::inadmissible); + CHECK(diagnostics.candidates[1].status == CoupledTrialCandidateStatus::inadmissible); + CHECK(diagnostics.preparations == 1); + } + CHECK(state.trialNormalizedState(0) == Approx(invalidMetric ? -.3125 : -.125)); + checkOuterData(state); +} + +TEST_CASE("Damped coupled search controls reject invalid or inconsistent work budgets collectively", + "[solver][coupled_trial_correction][backtracked][configuration][mpi]") { + using namespace coupled_trial_correction_test; + int rank = 0, ranks = 0; + REQUIRE(MPI_Comm_rank(MPI_COMM_WORLD, &rank) == MPI_SUCCESS); + REQUIRE(MPI_Comm_size(MPI_COMM_WORLD, &ranks) == MPI_SUCCESS); + ToyState state; + state.problem.communicator = MPI_COMM_WORLD; + BacktrackedCoupledStep search; + NearAcceptanceCoupledAttempt selection; + SECTION("A rank-local zero candidate budget cannot bound a useful search") { if (rank == 0) search.maximumSteps = 0; } + SECTION("A rank-local nondecreasing scale is invalid") { if (rank == 0) search.reduction = 1.0; } + SECTION("A rank-local zero scale is invalid") { if (rank == 0) search.reduction = 0.0; } + SECTION("A rank-local nonfinite scale is invalid") { if (rank == 0) search.reduction = std::numeric_limits::quiet_NaN(); } + SECTION("A nominally finite budget whose last scale underflows is invalid") { + if (rank == 0) search.maximumSteps = 4096; + } + SECTION("A rank-local zero acceptance allowance is invalid") { if (rank == 0) selection.maximumExcessRatio = 0.0; } + SECTION("A rank-local nonfinite allowance is invalid") { if (rank == 0) selection.maximumExcessRatio = std::numeric_limits::infinity(); } + SECTION("Different valid candidate budgets cannot share a collective search") { + if (ranks < 2) return; + search.maximumSteps = rank == 0 ? 4 : 3; + } + SECTION("Different valid reductions cannot share a collective search") { + if (ranks < 2) return; + search.reduction = rank == 0 ? .5 : .25; + } + SECTION("Different valid allowances cannot select one global attempt") { + if (ranks < 2) return; + selection.maximumExcessRatio = rank == 0 ? 1.0 : 2.0; + } + const auto create = [&] { + return solver::detail::PreparedCoupledTrialTransaction(state, exactOptions(), 1, search, selection); + }; + CHECK_THROWS_AS(create(), std::invalid_argument); +} diff --git a/tests/solver/stellar_numerical_strategies.cpp b/tests/solver/stellar_numerical_strategies.cpp index ca187ec..680fa5a 100644 --- a/tests/solver/stellar_numerical_strategies.cpp +++ b/tests/solver/stellar_numerical_strategies.cpp @@ -47,6 +47,21 @@ namespace stellar_numerical_strategy_test { static_assert(!MakesTrialNumerics); static_assert(!MakesTrialNumerics); + static_assert(nonlinear::CoupledTrialNumerics); + static_assert(nonlinear::CoupledTrialNumerics); + static_assert(!nonlinear::CoupledTrialNumerics); + static_assert(!nonlinear::TrialPathNumerics); + static_assert(MakesTrialNumerics); + static_assert(MakesTrialNumerics); + static_assert(!MakesTrialNumerics); + static_assert(!MakesTrialNumerics); + static_assert(std::same_as); [[nodiscard]] double matchingCentralDensity() { @@ -186,6 +201,11 @@ TEST_CASE("Typed numerical policies preserve the existing Newton configuration", const auto selected = nonlinear::Newton(nonlinear::NewtonOptions{}, nonlinear::NormalizedResidualMetric{}, nonlinear::NewtonNumerics{nonlinear::SampledEnthalpy{}, nonlinear::GravityCurvatureCorrection{}}); CHECK(selected.numerics().gravity.gravity.control.relativeTolerance == 1.e-10); + const nonlinear::BacktrackedCoupledCurvatureCorrection damped; + CHECK(damped.search.maximumSteps == 4); + CHECK(damped.search.reduction == .5); + CHECK(damped.attempt.maximumExcessRatio == 1.0); + CHECK(damped.maximumAttemptsPerIteration == 1); } TEST_CASE("Production gravity strategies retain optional model borders and restore rejected trials", @@ -300,25 +320,28 @@ TEST_CASE("Polynomial admissibility is invoked by each production gravity strate check(nonlinear::EliminatedGravity{permissiveGravity()}); REQUIRE(context.isReady()); checkPreserved(original, snapshot(context)); - { + const auto checkTrial = [&](auto correction) { const auto report = evaluate(context, nonlinear::NewtonNumerics{polynomial, - nonlinear::MonolithicGravity{}, nonlinear::CoupledCurvatureCorrection{}}); + nonlinear::MonolithicGravity{}, correction}); CHECK_FALSE(report.converged()); CHECK(report.diagnostics().numericalStrategies.polynomialChecks > 0); CHECK(report.diagnostics().acceptedNonlinearIterations == 0); CHECK(report.lastAcceptedCheckpointView().valid()); - } + }; + checkTrial(nonlinear::CoupledCurvatureCorrection{}); + REQUIRE(context.isReady()); + checkPreserved(original, snapshot(context)); + checkTrial(nonlinear::BacktrackedCoupledCurvatureCorrection{}); REQUIRE(context.isReady()); checkPreserved(original, snapshot(context)); } TEST_CASE("Coupled trial selection keeps the outer linear report and rejected stellar checkpoint", - "[solver][numerical_strategies][coupled_trial_correction][rollback]") { + "[solver][numerical_strategies][coupled_trial_correction][backtracked][rollback]") { using namespace stellar_numerical_strategy_test; auto context = makeContext(); const auto original = snapshot(context); - const auto check = [&](auto gravity) { - const nonlinear::CoupledCurvatureCorrection correction{}; + const auto check = [&](auto gravity, auto correction) { const auto report = evaluate(context, nonlinear::NewtonNumerics{ nonlinear::SampledEnthalpy{}, gravity, correction}); CHECK_FALSE(report.converged()); @@ -329,18 +352,36 @@ TEST_CASE("Coupled trial selection keeps the outer linear report and rejected st CHECK(outer.control.absoluteTolerance == 1.e100); CHECK(report.diagnostics().numericalStrategies.lastLinearReportSpace == nonlinear::LinearReportSpace::full_system); const auto &auxiliary = report.diagnostics().numericalStrategies.coupledTrial; - REQUIRE(auxiliary.currentSolve.has_value()); - CHECK(auxiliary.attempts == 1); CHECK(auxiliary.kept == 0); - CHECK(auxiliary.fallbacks == 1); - CHECK(auxiliary.currentSolve->linearSolve.control.relativeTolerance == correction.correction.control.relativeTolerance); - CHECK(auxiliary.currentSolve->linearSolve.control.absoluteTolerance == correction.correction.control.absoluteTolerance); + if constexpr (std::same_as) { + REQUIRE(auxiliary.currentSolve.has_value()); + CHECK(auxiliary.attempts == 1); + CHECK(auxiliary.fallbacks == 1); + CHECK(auxiliary.currentSolve->linearSolve.control.relativeTolerance == correction.correction.control.relativeTolerance); + CHECK(auxiliary.currentSolve->linearSolve.control.absoluteTolerance == correction.correction.control.absoluteTolerance); + } else { + // The captured zero outer correction predicts no decrease. The + // optional selector must leave the auxiliary budget untouched. + CHECK(auxiliary.status == solver::CoupledTrialAttemptStatus::prediction_deferred); + CHECK(auxiliary.deferredAttempts == 1); + CHECK(auxiliary.attempts == 0); + CHECK_FALSE(auxiliary.currentSolve.has_value()); + CHECK(auxiliary.candidateCount == 0); + CHECK(auxiliary.candidates.empty()); + CHECK(std::isnan(auxiliary.selectedBeta)); + } CHECK(report.lastAcceptedCheckpointView().valid()); }; - check(nonlinear::MonolithicGravity{}); + check(nonlinear::MonolithicGravity{}, nonlinear::CoupledCurvatureCorrection{}); REQUIRE(context.isReady()); checkPreserved(original, snapshot(context)); - check(nonlinear::GravityCurvatureCorrection{permissiveGravity()}); + check(nonlinear::GravityCurvatureCorrection{permissiveGravity()}, nonlinear::CoupledCurvatureCorrection{}); + REQUIRE(context.isReady()); + checkPreserved(original, snapshot(context)); + check(nonlinear::MonolithicGravity{}, nonlinear::BacktrackedCoupledCurvatureCorrection{}); + REQUIRE(context.isReady()); + checkPreserved(original, snapshot(context)); + check(nonlinear::GravityCurvatureCorrection{permissiveGravity()}, nonlinear::BacktrackedCoupledCurvatureCorrection{}); REQUIRE(context.isReady()); checkPreserved(original, snapshot(context)); }