tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add bounded coupled damping and prediction-based attempt selection

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-18 10:45:45 -04:00
parent 46249d0758
commit 422c46bcbf
11 changed files with 2296 additions and 101 deletions
Download Patch File Download Diff File Expand all files Collapse all files

67
experiments/DAMPED_COUPLED_TRIAL_2026-09-18.md Normal file
View File

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

407
experiments/check_damped_coupled_trial.py Normal file
View File

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

425
experiments/damped_coupled_screen_2026-09-18.json Normal file
View File

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

333
experiments/plateau_damped_coupled_diagnostics.hpp Normal file
View File

@@ -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 <class State, class Search, class Attempt>
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<State>;
using GravityCorrection = solver::PreparedGravityCurvatureCorrection<typename Gravity::Subsystem>;
using Transaction = solver::detail::PreparedCoupledTrialTransaction<State, Search, Attempt>;
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<double>::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<Search, solver::FullCoupledStep> ? 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> gravity;
std::optional<GravityCorrection> 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<Observation> 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 <class State>
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

104
experiments/plateau_diagnostic.cpp
View File

@@ -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<int>(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);

248
experiments/test_damped_coupled_trial_checker.py Normal file
View File

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

296
libmeanfield/interface/solver/coupled_trial_correction.cppm
View File

@@ -12,6 +12,7 @@ module;
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <vector>
#include <mfem.hpp>
#include <mpi.h>
@@ -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 <typename Candidate>
concept CoupledCorrectionSearch = std::same_as<Candidate, FullCoupledStep> ||
std::same_as<Candidate, BacktrackedCoupledStep>;
template <typename Candidate>
concept CoupledAttemptCriterion = std::same_as<Candidate, EveryRejectedCoupledAttempt> ||
std::same_as<Candidate, NearAcceptanceCoupledAttempt>;
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<double>::quiet_NaN()};
bool accepted{false};
double preparationSeconds{0.0};
};
struct CoupledTrialDiagnostics final {
@@ -64,6 +101,16 @@ export namespace mean_field::solver {
std::optional<CoupledTrialCorrectionReport> currentSolve{};
double originalMerit{std::numeric_limits<double>::quiet_NaN()};
double correctedMerit{std::numeric_limits<double>::quiet_NaN()};
std::uint64_t deferredAttempts{0};
int candidateCount{0};
double selectedBeta{std::numeric_limits<double>::quiet_NaN()};
double baseMerit{std::numeric_limits<double>::quiet_NaN()};
double linearPredictionMerit{std::numeric_limits<double>::quiet_NaN()};
double predictedDecrease{std::numeric_limits<double>::quiet_NaN()};
double predictionExcessRatio{std::numeric_limits<double>::quiet_NaN()};
// Search candidates exclude the possible final replay of an earlier
// best beta. Preparations/restorations still account for that work.
std::vector<CoupledTrialCandidateReport> 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 <typename State> class PreparedCoupledTrialTransaction final {
template <typename State, CoupledCorrectionSearch Search = FullCoupledStep,
CoupledAttemptCriterion Attempt = EveryRejectedCoupledAttempt>
class PreparedCoupledTrialTransaction final {
using Operation = std::remove_cvref_t<decltype(*std::declval<State>().normalizedOperator)>;
using Inverse = std::remove_cvref_t<decltype(*std::declval<State>().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<Attempt, NearAcceptanceCoupledAttempt> ? 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<double>::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<double>::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 <typename Preparation, typename Metric, typename Evaluate, typename Accepts, typename Better, typename CheckState>
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<Attempt, NearAcceptanceCoupledAttempt>) {
// 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<double>::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<Metric> 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 <typename Metric> 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<Search, BacktrackedCoupledStep>) return m_search.maximumSteps;
else return 1;
}
[[nodiscard]] double Reduction() const noexcept {
if constexpr (std::same_as<Search, BacktrackedCoupledStep>) 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<Search, BacktrackedCoupledStep>) m_search.Validate();
if constexpr (std::same_as<Attempt, NearAcceptanceCoupledAttempt>) 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<Attempt, NearAcceptanceCoupledAttempt>) excessRatio = m_attempt.maximumExcessRatio;
const std::array<int, 4> integers{m_maximumAttempts, std::same_as<Search, BacktrackedCoupledStep>,
MaximumSteps(), std::same_as<Attempt, NearAcceptanceCoupledAttempt>};
const std::array<double, 2> reals{Reduction(), excessRatio};
std::array<int, 4> minimumIntegers{}, maximumIntegers{};
std::array<double, 2> 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<Operation, Inverse> 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{};
};
}

21
libmeanfield/interface/solver/newton_numerics.cppm
View File

@@ -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 <typename Candidate>
concept CoupledTrialNumerics = std::same_as<Candidate, CoupledCurvatureCorrection> ||
std::same_as<Candidate, BacktrackedCoupledCurvatureCorrection>;
template <typename Candidate>
concept EnthalpyNumerics = std::same_as<Candidate, SampledEnthalpy> ||
std::same_as<Candidate, PolynomialEnthalpy>;
@@ -64,7 +83,7 @@ export namespace mean_field::solver::nonlinear {
template <typename Candidate>
concept TrialPathNumerics = std::same_as<Candidate, UncorrectedTrial> ||
std::same_as<Candidate, CoupledCurvatureCorrection>;
CoupledTrialNumerics<Candidate>;
// Types select algorithms; tolerances and the mesh-dependent work sizes
// remain runtime data. Adding a strategy cannot silently change a model's

35
libmeanfield/interface/solver/stellar_newton_strategy.cppm
View File

@@ -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<int, 13> integers{};
std::array<double, 5> reals{};
std::array<int, 14> integers{};
std::array<double, 7> reals{};
integers[0] = std::same_as<Enthalpy, nonlinear::PolynomialEnthalpy>;
integers[1] = std::same_as<Gravity, nonlinear::GravityCurvatureCorrection> ? 1 :
std::same_as<Gravity, nonlinear::EliminatedGravity> ? 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<Trial, nonlinear::CoupledCurvatureCorrection>;
if constexpr (std::same_as<Trial, nonlinear::CoupledCurvatureCorrection>) {
integers[8] = std::same_as<Trial, nonlinear::CoupledCurvatureCorrection> ? 1 :
std::same_as<Trial, nonlinear::BacktrackedCoupledCurvatureCorrection> ? 2 : 0;
if constexpr (nonlinear::CoupledTrialNumerics<Trial>) {
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<Trial, nonlinear::BacktrackedCoupledCurvatureCorrection>) {
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<int, 13> minimumIntegers{}, maximumIntegers{};
std::array<double, 5> 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<int, 14> minimumIntegers{}, maximumIntegers{};
std::array<double, 7> 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<Gravity, nonlinear::MonolithicGravity>) { recipe.gravity.gravity.Validate(); }
if constexpr (std::same_as<Trial, nonlinear::CoupledCurvatureCorrection>) { recipe.trial.Validate(); }
if constexpr (nonlinear::CoupledTrialNumerics<Trial>) { recipe.trial.Validate(); }
}
template <typename State, typename Policy> class PreparedGravityStrategy;
@@ -431,6 +437,11 @@ namespace mean_field::solver::detail {
: transaction(state, policy.correction, policy.maximumAttemptsPerIteration) { }
PreparedCoupledTrialTransaction<State> transaction;
};
template <typename State> struct PreparedTrialPath<State, nonlinear::BacktrackedCoupledCurvatureCorrection> {
PreparedTrialPath(State &state, const nonlinear::BacktrackedCoupledCurvatureCorrection &policy)
: transaction(state, policy.correction, policy.maximumAttemptsPerIteration, policy.search, policy.attempt) { }
PreparedCoupledTrialTransaction<State, BacktrackedCoupledStep, NearAcceptanceCoupledAttempt> transaction;
};
template <typename State, nonlinear::NewtonNumericalConfiguration Recipe>
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<Enthalpy, nonlinear::PolynomialEnthalpy>;
static constexpr bool correctsTrial = std::same_as<Trial, nonlinear::CoupledCurvatureCorrection>;
static constexpr bool correctsTrial = nonlinear::CoupledTrialNumerics<Trial>;
PreparedNewtonStrategy(State &state, const Recipe &recipe)
: m_enthalpy(state, recipe.enthalpy), m_gravity(state, recipe.gravity), m_trial(state, recipe.trial) { }

394
tests/solver/coupled_trial_correction.cpp
View File

@@ -7,6 +7,7 @@
#include <memory>
#include <stdexcept>
#include <utility>
#include <vector>
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
@@ -162,6 +163,12 @@ namespace coupled_trial_correction_test {
static_assert(CoupledTrialInverse<ExactInverse>);
static_assert(!CoupledTrialOperation<mfem::Operator>);
static_assert(!CoupledTrialInverse<mfem::Solver>);
static_assert(CoupledCorrectionSearch<FullCoupledStep>);
static_assert(CoupledCorrectionSearch<BacktrackedCoupledStep>);
static_assert(!CoupledCorrectionSearch<int>);
static_assert(CoupledAttemptCriterion<EveryRejectedCoupledAttempt>);
static_assert(CoupledAttemptCriterion<NearAcceptanceCoupledAttempt>);
static_assert(!CoupledAttemptCriterion<int>);
static_assert(LinearBackendRuntimeAvailableFor<TestBackend, DenseOperation, ExactInverse>);
using ScriptedCorrection = PreparedCoupledTrialCorrection<DenseOperation, ExactInverse, TestBackend>;
} // 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<double> preparedCoordinates;
ToyState() {
trialNormalizedState = 0.0;
@@ -449,23 +461,37 @@ namespace coupled_trial_correction_test {
[[nodiscard]] std::expected<void, int> 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<ToyState, BacktrackedCoupledStep> 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<std::uint64_t>(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<ToyState, BacktrackedCoupledStep> 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<ToyState, BacktrackedCoupledStep> 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<double>::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<ToyState, BacktrackedCoupledStep> 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<ToyState, BacktrackedCoupledStep,
NearAcceptanceCoupledAttempt> 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<ToyState, BacktrackedCoupledStep,
NearAcceptanceCoupledAttempt> 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<ToyState, BacktrackedCoupledStep> 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<double>::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<double>::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<double>::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);
}

67
tests/solver/stellar_numerical_strategies.cpp
View File

@@ -47,6 +47,21 @@ namespace stellar_numerical_strategy_test {
static_assert(!MakesTrialNumerics<nonlinear::PolynomialEnthalpy, nonlinear::EliminatedGravity,
nonlinear::CoupledCurvatureCorrection>);
static_assert(!MakesTrialNumerics<nonlinear::SampledEnthalpy, nonlinear::MonolithicGravity, UnregisteredTrial>);
static_assert(nonlinear::CoupledTrialNumerics<nonlinear::CoupledCurvatureCorrection>);
static_assert(nonlinear::CoupledTrialNumerics<nonlinear::BacktrackedCoupledCurvatureCorrection>);
static_assert(!nonlinear::CoupledTrialNumerics<nonlinear::UncorrectedTrial>);
static_assert(!nonlinear::TrialPathNumerics<solver::BacktrackedCoupledStep>);
static_assert(MakesTrialNumerics<nonlinear::SampledEnthalpy, nonlinear::MonolithicGravity,
nonlinear::BacktrackedCoupledCurvatureCorrection>);
static_assert(MakesTrialNumerics<nonlinear::PolynomialEnthalpy, nonlinear::GravityCurvatureCorrection,
nonlinear::BacktrackedCoupledCurvatureCorrection>);
static_assert(!MakesTrialNumerics<nonlinear::SampledEnthalpy, nonlinear::EliminatedGravity,
nonlinear::BacktrackedCoupledCurvatureCorrection>);
static_assert(!MakesTrialNumerics<nonlinear::PolynomialEnthalpy, nonlinear::EliminatedGravity,
nonlinear::BacktrackedCoupledCurvatureCorrection>);
static_assert(std::same_as<typename decltype(nonlinear::NewtonNumerics{nonlinear::SampledEnthalpy{},
nonlinear::MonolithicGravity{}, nonlinear::BacktrackedCoupledCurvatureCorrection{}})::TrialType,
nonlinear::BacktrackedCoupledCurvatureCorrection>);
[[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<false>();
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<decltype(correction), nonlinear::CoupledCurvatureCorrection>) {
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));
}