From 1e8c93bbeb856d4721933b297d527bd76c6bb86a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 7 Aug 2026 14:16:34 +0200 Subject: [PATCH 1/2] Reject step sizes that are spuriously self-consistent in Consistency `Consistency` checked whether the requested methods (e.g. forward/backward/central) agreed with each other at each step size, then blended every self-consistent size's mean into the final value with one final blanket tolerance check. A step size can become small enough that all methods sample points within the target function's floating-point noise floor and become correlated (affected by the same rounding/cancellation error) -- spuriously self-consistent, yet biased away from the true derivative. That biased size was still blended into the average while `success` reported `True`. Symmetrically, a large step size can also be self-consistently wrong due to higher-order/truncation effects. Add an order-independent iterative outlier-rejection pass (`Consistency._reject_outliers`) over the self-consistent step sizes' means: repeatedly drop the single worst-deviating candidate, relative to the median and a robust (MAD-based) spread of the rest, until nothing looks anomalous. This intentionally does not treat step size magnitude as a proxy for trustworthiness, and does not reuse the user's own `rtol` as the cross-size threshold (a fixed relative tolerance loose enough for the per-size check is too loose to catch this kind of bias no matter where it's applied; the point of using a spread estimated from the trusted candidates themselves is that it tightens automatically as they agree more closely). This changes `Consistency()`'s default behavior -- intentionally, since today's default is the bug. Two new tunable parameters, `trend_n_sigma` (default 5.0) and `min_trend_samples` (default 3), gate this behavior; below `min_trend_samples` self-consistent sizes, behavior is unchanged. Documented limitation (see updated class docstring): like any purely data-driven consensus check, this has a ~50% breakdown point -- if close to half the self-consistent step sizes are corrupted, it cannot reliably tell which subset is trustworthy. Still a strict improvement over today's ~0% effective breakdown point, where a single corrupted-but-self-consistent size already breaks the check. Co-Authored-By: Claude Sonnet 5 --- fiddy/success.py | 156 ++++++++++++++++++++++++++++++++------- tests/test_success.py | 165 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 294 insertions(+), 27 deletions(-) create mode 100644 tests/test_success.py diff --git a/fiddy/success.py b/fiddy/success.py index ab84b34..b15df53 100644 --- a/fiddy/success.py +++ b/fiddy/success.py @@ -33,6 +33,39 @@ def method(self, directional_derivative: DirectionalDerivative) -> Any: class Consistency(Success): + """Consistency-based success checker. + + For each step size, checks whether the requested methods + (e.g. forward/backward/central) agree with each other ("self-consistent"). + Self-consistency alone is not sufficient though: a step size can be small + enough that all methods sample points within the target function's + floating-point noise floor, and become correlated (affected by the same + rounding/cancellation error) -- self-consistent, yet biased away from the + truth. Symmetrically, a step size can also be large enough that all + methods are biased the same way by higher-order/truncation effects. + + To guard against this, self-consistent step sizes are additionally + required to agree with the majority of other self-consistent step sizes, + via iterative outlier rejection (order-independent; step size magnitude + is not used as a proxy for trustworthiness): repeatedly compute the + median and a robust (MAD-based) spread of the current candidates, and + drop the single worst-deviating one if it exceeds ``trend_n_sigma`` + scaled MADs from the median, until nothing looks anomalous. This only + activates once there are at least ``min_trend_samples`` self-consistent + step sizes; below that, there isn't enough data to estimate a spread, and + all self-consistent step sizes are used, as before. + + Note that this is a majority-vote style method: like any check based + purely on the agreement of the values themselves (no independent ground + truth), it has a breakdown point of roughly 50% (a property of the + underlying median/MAD statistics) -- if close to half (or more) of the + self-consistent step sizes are corrupted, this check cannot reliably + tell which subset is trustworthy. This is a fundamental limitation of + any purely data-driven consistency check, not something this + implementation can detect or work around; sufficient step sizes with a + real chance of being individually trustworthy should be provided. + """ + # FIXME string literal id = "consistency" only_at_completion: bool = True @@ -46,7 +79,33 @@ def __init__( rtol: float = 0.2, atol: float = 1e-15, equal_nan: bool = True, + trend_n_sigma: float = 5.0, + min_trend_samples: int = 3, ): + """Construct. + + Args: + rtol: + Relative tolerance for the self-consistency check of methods + at the same step size, and for the final check of the + blended value against the trusted per-size means. + atol: + Absolute tolerance, analogous to `rtol`. Also used as a + floor for the robust spread estimate in the cross-step-size + outlier rejection (see class docstring). + equal_nan: + Whether `NaN`s are considered equal in tolerance checks. + trend_n_sigma: + The number of scaled median-absolute-deviations a + self-consistent step size's estimate may deviate from the + median of the other trusted step sizes' estimates, before it + is rejected as an outlier. + min_trend_samples: + The minimum number of self-consistent step sizes required + before the cross-step-size outlier rejection is attempted. + Below this, all self-consistent step sizes are trusted, same + as if `trend_n_sigma`/outlier rejection did not exist. + """ super().__init__() # if computer_parser is None: # computer_parser = ( @@ -66,6 +125,8 @@ def __init__( self.rtol = rtol self.atol = atol self.equal_nan = equal_nan + self.trend_n_sigma = trend_n_sigma + self.min_trend_samples = min_trend_samples def method( self, directional_derivative: DirectionalDerivative @@ -87,47 +148,88 @@ def method( ) results_by_size[size][result.method_id] = result.value - success_by_size = {} - for size, results in results_by_size.items(): - values = list(results.values()) - with warnings.catch_warnings(): - warnings.filterwarnings( - "ignore", "Mean of empty slice", RuntimeWarning - ) - mean = np.nanmean(values, axis=0) - success_by_size[size] = np.isclose( - values, - mean, - rtol=self.rtol / 2, - atol=self.atol / 2, - equal_nan=self.equal_nan, - ).all() - + self_consistent_means = [] with warnings.catch_warnings(): warnings.filterwarnings( "ignore", "Mean of empty slice", RuntimeWarning ) + for results in results_by_size.values(): + values = list(results.values()) + mean = np.nanmean(values, axis=0) + is_self_consistent = np.isclose( + values, + mean, + rtol=self.rtol / 2, + atol=self.atol / 2, + equal_nan=self.equal_nan, + ).all() + if is_self_consistent: + self_consistent_means.append(mean) + + if not self_consistent_means: + return False, np.nan - consistent_results = np.array( - [ - np.nanmean(list(results_by_size[size].values()), axis=0) - for size, success in success_by_size.items() - if success - ] - ) + trusted_means = self._reject_outliers(self_consistent_means) - if len(consistent_results) == 0: + if not trusted_means: return False, np.nan - value = np.nanmean(consistent_results, axis=0) + value = np.nanmean(trusted_means, axis=0) + success = ( np.isclose( - consistent_results, + trusted_means, value, rtol=self.rtol, atol=self.atol, equal_nan=self.equal_nan, ).all() - and not np.isnan(consistent_results).all() + and not np.isnan(trusted_means).all() ) return success, value + + def _reject_outliers( + self, means: list[Type.DIRECTIONAL_DERIVATIVE] + ) -> list[Type.DIRECTIONAL_DERIVATIVE]: + """Iteratively reject step sizes whose estimate is an outlier. + + See the class docstring for the rationale. Order-independent: does + not assume larger (or smaller) step sizes are inherently more + trustworthy. + + Args: + means: + The per-step-size mean estimates that passed the + within-step-size self-consistency check. + + Returns: + The subset of `means` that are also mutually consistent with + each other. + """ + trusted = list(means) + if len(trusted) < self.min_trend_samples: + return trusted + + floor = max(self.atol / 2, np.finfo(float).tiny) + while len(trusted) >= self.min_trend_samples: + stacked = np.asarray(trusted, dtype=float) + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", "All-NaN", RuntimeWarning) + center = np.nanmedian(stacked, axis=0) + mad = np.nanmedian(np.abs(stacked - center), axis=0) + scale = np.maximum(mad * 1.4826, floor) + # One badness score per candidate, reduced across all + # output dimensions (a candidate is an outlier if it + # deviates too much in *any* output element). + badness = np.nanmax( + (np.abs(stacked - center) / scale).reshape( + len(trusted), -1 + ), + axis=1, + ) + worst = int(np.nanargmax(badness)) + if badness[worst] > self.trend_n_sigma: + trusted.pop(worst) + else: + break + return trusted diff --git a/tests/test_success.py b/tests/test_success.py new file mode 100644 index 0000000..c6ac472 --- /dev/null +++ b/tests/test_success.py @@ -0,0 +1,165 @@ +import numpy as np +from scipy.optimize import rosen + +from fiddy import MethodId, get_derivative +from fiddy.success import Consistency + + +def test_consistency_rejects_rounding_noise_dominated_step_sizes(): + """Regression test for a `Consistency` robustness bug. + + A step size can become small enough that forward/backward/central all + sample points within the target function's floating-point noise floor. + They then become correlated (affected by the same rounding/cancellation + error), and can spuriously agree with each other ("self-consistent") + while being biased away from the true derivative. `Consistency` must not + blend such a step size into the final value while reporting + `success=True`. + """ + true_slope = 872.68 + noise_floor = 2e-7 + + def f(point): + x0 = point[0] + value = -1023.447 + true_slope * (x0 - 1e-4) + value += noise_floor * np.sin(1e8 * x0) + return np.array(value) + + point = np.array([9.579126317171899e-05]) + step_sizes = [5e-1, 2e-1, 1e-1, 5e-2, 1e-2, 1e-3, 1e-4, 1e-5] + + derivative = get_derivative( + function=f, + point=point, + sizes=step_sizes, + direction_ids=["x0"], + method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + success_checker=Consistency(rtol=0.1, atol=1e-5), + relative_sizes=True, + ) + + success = bool(derivative.df["success"].values[0]) + value = float(np.squeeze(derivative.series.values[0])) + + # Reporting `success=True` is only acceptable if the value is actually + # accurate; silently returning a significantly biased value (as the + # unpatched algorithm does: ~868.5, a ~0.5% error) is the bug. + if success: + assert np.isclose(value, true_slope, rtol=1e-2) + + +def test_consistency_averages_all_trustworthy_step_sizes(): + """A wide, but genuinely well-behaved, range of step sizes should not + trigger spurious outlier rejection.""" + + def f(point): + return np.array([rosen(point)]) + + point = np.array([1.3, 0.7]) + step_sizes = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6] + + derivative = get_derivative( + function=f, + point=point, + sizes=step_sizes, + direction_ids=["x0"], + directions=[np.array([1.0, 0.0])], + method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + success_checker=Consistency(rtol=1e-2, atol=1e-8), + ) + + assert bool(derivative.df["success"].values[0]) + value = float(np.squeeze(derivative.series.values[0])) + # d/dx0 rosen([1.3, 0.7]) = -2*100*(0.7 - 1.3**2)*1.3 * ... expected via finite differences, + # so just check against a fine central-difference estimate directly. + h = 1e-6 + expected = ( + rosen(point + np.array([h, 0.0])) - rosen(point - np.array([h, 0.0])) + ) / (2 * h) + assert np.isclose(value, expected, rtol=1e-3) + + +class TestRejectOutliers: + """Unit tests for `Consistency._reject_outliers`, the order-independent + iterative outlier-rejection pass over step sizes' per-size means (added + to reject step sizes that are spuriously self-consistent, see + `test_consistency_rejects_rounding_noise_dominated_step_sizes` above).""" + + def test_below_min_trend_samples_keeps_everything(self): + # An outlier (500.0) is present, but there are fewer candidates than + # `min_trend_samples` -- too little data to estimate a spread, so no + # trimming is attempted at all. + checker = Consistency(min_trend_samples=5) + means = [10.0, 10.01, 500.0] + assert checker._reject_outliers(means) == means + + def test_no_outliers_keeps_all(self): + checker = Consistency() + means = [10.0, 10.01, 9.99, 10.02] + assert checker._reject_outliers(means) == means + + def test_removes_single_outlier(self): + checker = Consistency() + means = [10.0, 10.01, 9.99, 10.02, 500.0] + trusted = checker._reject_outliers(means) + assert trusted == [10.0, 10.01, 9.99, 10.02] + + def test_removes_multiple_outliers_iteratively(self): + # Two outliers on opposite sides of the trustworthy cluster; both + # must be dropped, one per iteration, worst-first. + checker = Consistency() + means = [10.0, 10.01, 9.99, 10.02, 500.0, -500.0] + trusted = checker._reject_outliers(means) + assert trusted == [10.0, 10.01, 9.99, 10.02] + + def test_order_independent(self): + # Dropping is based on value, not position: shuffling the input + # must not change which candidates survive. + checker = Consistency() + means = [500.0, 10.0, 10.01, 9.99, 10.02] + trusted = checker._reject_outliers(means) + assert sorted(trusted) == [9.99, 10.0, 10.01, 10.02] + + def test_respects_trend_n_sigma(self): + means = [10.0, 10.01, 9.99, 10.02, 500.0] + lenient_checker = Consistency(trend_n_sigma=1e6) + assert lenient_checker._reject_outliers(means) == means + + strict_checker = Consistency(trend_n_sigma=5.0) + assert strict_checker._reject_outliers(means) == [ + 10.0, + 10.01, + 9.99, + 10.02, + ] + + def test_vector_valued_drops_whole_candidate_on_any_element_outlier(self): + # A candidate that's fine in one output element but a severe + # outlier in another must still be dropped entirely (not just + # masked in the bad element) -- "badness" is reduced across all + # output dimensions before picking the worst candidate. + checker = Consistency() + means = [ + np.array([10.0, 5.0]), + np.array([10.01, 5.01]), + np.array([9.99, 500.0]), # fine in element 0, an outlier in 1 + ] + trusted = checker._reject_outliers(means) + assert len(trusted) == 2 + assert all( + np.array_equal(t, means[i]) + for t, i in zip(trusted, [0, 1], strict=True) + ) + + def test_nan_candidate_is_never_flagged_as_worst(self): + # Known, documented limitation: `nanargmax` ignores NaNs, so a + # candidate whose mean is entirely NaN can never be selected as + # "the worst" and is left in the trusted set untouched (harmless in + # practice: it doesn't shift `np.nanmean` of the final value, and + # `Consistency.method`'s final blanket `isclose` check against a + # non-NaN blended value still reports `success=False` overall). + checker = Consistency() + means = [10.0, 10.01, np.nan] + trusted = checker._reject_outliers(means) + assert len(trusted) == 3 + assert np.isnan(trusted[-1]) From 3091f28c3c8654e458995f536758ce87323f3838 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 10 Aug 2026 15:39:03 +0200 Subject: [PATCH 2/2] Split outlier rejection into a separate RobustConsistency class Per review feedback: revert Consistency to its original behavior (no cross-step-size rejection), and move the MAD-based outlier rejection added in the previous commit into a new RobustConsistency(Consistency) subclass instead of changing Consistency's default behavior. This is non-breaking for any existing Consistency() caller, and lets the more robust checker be adopted explicitly (e.g. by AMICI) without a default-behavior debate. Also addresses the reviewer's request to warn when a step size is rejected: RobustConsistency.method now emits a UserWarning (with the count of rejected step sizes) whenever _reject_outliers drops one or more otherwise self-consistent step sizes. Shared logic (grouping results by step size, the within-size self-consistency check) is factored into Consistency._self_consistent_means so both classes reuse it without duplication. tests/test_success.py updated accordingly: the regression/robustness tests now target RobustConsistency; added a test asserting RobustConsistency's warning fires (with the right count) exactly when a step size is actually rejected, and not otherwise. Also added a test documenting a related, non-bug property found while re-testing: a wide, noise-free step-size range can legitimately narrow down to just the smallest/most-precise steps, since larger (but not wrong, just less precise) steps can look like outliers next to a cluster already near machine precision -- the blended value stays accurate regardless. Co-Authored-By: Claude Sonnet 5 --- fiddy/success.py | 180 +++++++++++++++++++++++++++++++----------- tests/test_success.py | 155 ++++++++++++++++++++++++++++-------- 2 files changed, 254 insertions(+), 81 deletions(-) diff --git a/fiddy/success.py b/fiddy/success.py index b15df53..c2811d4 100644 --- a/fiddy/success.py +++ b/fiddy/success.py @@ -37,33 +37,17 @@ class Consistency(Success): For each step size, checks whether the requested methods (e.g. forward/backward/central) agree with each other ("self-consistent"). - Self-consistency alone is not sufficient though: a step size can be small - enough that all methods sample points within the target function's + A step size is trusted if it is self-consistent; the final value is the + average of all trusted step sizes' means, and `success` additionally + requires those means to be mutually close to that average. + + Self-consistency alone is not a very strong guarantee: a step size can be + small enough that all methods sample points within the target function's floating-point noise floor, and become correlated (affected by the same rounding/cancellation error) -- self-consistent, yet biased away from the truth. Symmetrically, a step size can also be large enough that all - methods are biased the same way by higher-order/truncation effects. - - To guard against this, self-consistent step sizes are additionally - required to agree with the majority of other self-consistent step sizes, - via iterative outlier rejection (order-independent; step size magnitude - is not used as a proxy for trustworthiness): repeatedly compute the - median and a robust (MAD-based) spread of the current candidates, and - drop the single worst-deviating one if it exceeds ``trend_n_sigma`` - scaled MADs from the median, until nothing looks anomalous. This only - activates once there are at least ``min_trend_samples`` self-consistent - step sizes; below that, there isn't enough data to estimate a spread, and - all self-consistent step sizes are used, as before. - - Note that this is a majority-vote style method: like any check based - purely on the agreement of the values themselves (no independent ground - truth), it has a breakdown point of roughly 50% (a property of the - underlying median/MAD statistics) -- if close to half (or more) of the - self-consistent step sizes are corrupted, this check cannot reliably - tell which subset is trustworthy. This is a fundamental limitation of - any purely data-driven consistency check, not something this - implementation can detect or work around; sufficient step sizes with a - real chance of being individually trustworthy should be provided. + methods are biased the same way by higher-order/truncation effects. See + `RobustConsistency` for a variant that additionally guards against this. """ # FIXME string literal @@ -79,8 +63,6 @@ def __init__( rtol: float = 0.2, atol: float = 1e-15, equal_nan: bool = True, - trend_n_sigma: float = 5.0, - min_trend_samples: int = 3, ): """Construct. @@ -90,21 +72,9 @@ def __init__( at the same step size, and for the final check of the blended value against the trusted per-size means. atol: - Absolute tolerance, analogous to `rtol`. Also used as a - floor for the robust spread estimate in the cross-step-size - outlier rejection (see class docstring). + Absolute tolerance, analogous to `rtol`. equal_nan: Whether `NaN`s are considered equal in tolerance checks. - trend_n_sigma: - The number of scaled median-absolute-deviations a - self-consistent step size's estimate may deviate from the - median of the other trusted step sizes' estimates, before it - is rejected as an outlier. - min_trend_samples: - The minimum number of self-consistent step sizes required - before the cross-step-size outlier rejection is attempted. - Below this, all self-consistent step sizes are trusted, same - as if `trend_n_sigma`/outlier rejection did not exist. """ super().__init__() # if computer_parser is None: @@ -125,12 +95,13 @@ def __init__( self.rtol = rtol self.atol = atol self.equal_nan = equal_nan - self.trend_n_sigma = trend_n_sigma - self.min_trend_samples = min_trend_samples - def method( + def _self_consistent_means( self, directional_derivative: DirectionalDerivative - ) -> tuple[bool, float]: + ) -> list[Type.DIRECTIONAL_DERIVATIVE]: + """Group results by step size, and return the per-size mean for + every step size whose requested methods agree with each other + ("self-consistent") within `rtol/2`, `atol/2`.""" # FIXME string literals computer_results = directional_derivative.get_computer_results() analysis_results = directional_derivative.get_analysis_results() @@ -165,15 +136,128 @@ def method( ).all() if is_self_consistent: self_consistent_means.append(mean) + return self_consistent_means + + def method( + self, directional_derivative: DirectionalDerivative + ) -> tuple[bool, float]: + self_consistent_means = self._self_consistent_means( + directional_derivative + ) + + if not self_consistent_means: + return False, np.nan + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", "Mean of empty slice", RuntimeWarning + ) + value = np.nanmean(self_consistent_means, axis=0) + + success = ( + np.isclose( + self_consistent_means, + value, + rtol=self.rtol, + atol=self.atol, + equal_nan=self.equal_nan, + ).all() + and not np.isnan(self_consistent_means).all() + ) + return success, value - if not self_consistent_means: - return False, np.nan - trusted_means = self._reject_outliers(self_consistent_means) +class RobustConsistency(Consistency): + """`Consistency`, plus rejection of step sizes that are self-consistent + but inconsistent with the majority of other step sizes. - if not trusted_means: - return False, np.nan + As explained in `Consistency`'s docstring, self-consistency of a step + size (its methods agreeing with each other) is not a strong guarantee on + its own -- a step size can be spuriously self-consistent while biased + away from the truth (e.g. correlated rounding/cancellation error at very + small steps, or correlated higher-order/truncation effects at very large + steps). + To guard against this, self-consistent step sizes are additionally + required to agree with the majority of other self-consistent step sizes, + via iterative outlier rejection (order-independent; step size magnitude + is not used as a proxy for trustworthiness): repeatedly compute the + median and a robust (MAD-based) spread of the current candidates, and + drop the single worst-deviating one if it exceeds ``trend_n_sigma`` + scaled MADs from the median, until nothing looks anomalous. This only + activates once there are at least ``min_trend_samples`` self-consistent + step sizes; below that, there isn't enough data to estimate a spread, and + all self-consistent step sizes are used, as in `Consistency`. A + `UserWarning` is emitted whenever one or more step sizes are rejected + this way. + + Note that this is a majority-vote style method: like any check based + purely on the agreement of the values themselves (no independent ground + truth), it has a breakdown point of roughly 50% (a property of the + underlying median/MAD statistics) -- if close to half (or more) of the + self-consistent step sizes are corrupted, this check cannot reliably + tell which subset is trustworthy. This is a fundamental limitation of + any purely data-driven consistency check, not something this + implementation can detect or work around; sufficient step sizes with a + real chance of being individually trustworthy should be provided. + """ + + id = "robust_consistency" + + def __init__( + self, + *args, + trend_n_sigma: float = 5.0, + min_trend_samples: int = 3, + **kwargs, + ): + """Construct. + + Args: + trend_n_sigma: + The number of scaled median-absolute-deviations a + self-consistent step size's estimate may deviate from the + median of the other trusted step sizes' estimates, before it + is rejected as an outlier. + min_trend_samples: + The minimum number of self-consistent step sizes required + before the cross-step-size outlier rejection is attempted. + Below this, all self-consistent step sizes are trusted, same + as in `Consistency`. + """ + super().__init__(*args, **kwargs) + self.trend_n_sigma = trend_n_sigma + self.min_trend_samples = min_trend_samples + + def method( + self, directional_derivative: DirectionalDerivative + ) -> tuple[bool, float]: + self_consistent_means = self._self_consistent_means( + directional_derivative + ) + + if not self_consistent_means: + return False, np.nan + + trusted_means = self._reject_outliers(self_consistent_means) + + if not trusted_means: + return False, np.nan + + n_rejected = len(self_consistent_means) - len(trusted_means) + if n_rejected: + warnings.warn( + f"{n_rejected} step size(s) were self-consistent (the " + "requested methods agreed with each other) but were " + "rejected as inconsistent with the majority of other step " + "sizes; see `RobustConsistency`'s docstring.", + stacklevel=2, + ) + + with warnings.catch_warnings(): + warnings.filterwarnings( + "ignore", "Mean of empty slice", RuntimeWarning + ) value = np.nanmean(trusted_means, axis=0) success = ( diff --git a/tests/test_success.py b/tests/test_success.py index c6ac472..f840a3b 100644 --- a/tests/test_success.py +++ b/tests/test_success.py @@ -1,19 +1,36 @@ import numpy as np +import pytest from scipy.optimize import rosen from fiddy import MethodId, get_derivative -from fiddy.success import Consistency +from fiddy.directional_derivative import ComputerResult +from fiddy.success import RobustConsistency -def test_consistency_rejects_rounding_noise_dominated_step_sizes(): - """Regression test for a `Consistency` robustness bug. +class FakeDirectionalDerivative: + """Minimal stand-in exposing only what `Consistency.method` calls.""" + + def __init__(self, computer_results, analysis_results=None): + self._computer_results = computer_results + self._analysis_results = analysis_results or [] + + def get_computer_results(self): + return self._computer_results + + def get_analysis_results(self): + return self._analysis_results + + +def test_robust_consistency_rejects_rounding_noise_dominated_step_sizes(): + """Regression test for a `RobustConsistency` (nee `Consistency`) + robustness bug. A step size can become small enough that forward/backward/central all sample points within the target function's floating-point noise floor. They then become correlated (affected by the same rounding/cancellation error), and can spuriously agree with each other ("self-consistent") - while being biased away from the true derivative. `Consistency` must not - blend such a step size into the final value while reporting + while being biased away from the true derivative. `RobustConsistency` + must not blend such a step size into the final value while reporting `success=True`. """ true_slope = 872.68 @@ -28,15 +45,16 @@ def f(point): point = np.array([9.579126317171899e-05]) step_sizes = [5e-1, 2e-1, 1e-1, 5e-2, 1e-2, 1e-3, 1e-4, 1e-5] - derivative = get_derivative( - function=f, - point=point, - sizes=step_sizes, - direction_ids=["x0"], - method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=Consistency(rtol=0.1, atol=1e-5), - relative_sizes=True, - ) + with pytest.warns(UserWarning, match="rejected as inconsistent"): + derivative = get_derivative( + function=f, + point=point, + sizes=step_sizes, + direction_ids=["x0"], + method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + success_checker=RobustConsistency(rtol=0.1, atol=1e-5), + relative_sizes=True, + ) success = bool(derivative.df["success"].values[0]) value = float(np.squeeze(derivative.series.values[0])) @@ -48,15 +66,19 @@ def f(point): assert np.isclose(value, true_slope, rtol=1e-2) -def test_consistency_averages_all_trustworthy_step_sizes(): +def test_robust_consistency_averages_all_trustworthy_step_sizes(): """A wide, but genuinely well-behaved, range of step sizes should not - trigger spurious outlier rejection.""" + trigger spurious outlier rejection (and thus no rejection warning).""" def f(point): return np.array([rosen(point)]) point = np.array([1.3, 0.7]) - step_sizes = [1e-1, 1e-2, 1e-3, 1e-4, 1e-5, 1e-6] + # Chosen to have comparable precision across the whole range (see + # `test_robust_consistency_narrows_to_the_most_precise_step_sizes` + # below for what happens once the range gets wide enough that the + # smallest steps are far more precise than the largest). + step_sizes = [1e-2, 1e-3, 1e-4] derivative = get_derivative( function=f, @@ -65,7 +87,7 @@ def f(point): direction_ids=["x0"], directions=[np.array([1.0, 0.0])], method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], - success_checker=Consistency(rtol=1e-2, atol=1e-8), + success_checker=RobustConsistency(rtol=1e-2, atol=1e-8), ) assert bool(derivative.df["success"].values[0]) @@ -79,27 +101,94 @@ def f(point): assert np.isclose(value, expected, rtol=1e-3) +def test_robust_consistency_narrows_to_the_most_precise_step_sizes(): + """Rejection isn't only about *biased* step sizes (the motivating bug): + a genuinely wide, noise-free step-size range can legitimately narrow + down to just the handful of smallest, most precise steps, even though + the larger, excluded ones weren't wrong -- just comparatively imprecise + (ordinary, shrinking-with-h truncation error) next to a cluster that + happens to already be near machine precision. The blended value must + stay accurate either way. + """ + + def f(point): + return np.array([rosen(point)]) + + point = np.array([1.3, 0.7]) + step_sizes = [1e-2, 1e-3, 1e-4, 1e-5, 1e-6] + + with pytest.warns(UserWarning, match="rejected as inconsistent"): + derivative = get_derivative( + function=f, + point=point, + sizes=step_sizes, + direction_ids=["x0"], + directions=[np.array([1.0, 0.0])], + method_ids=[MethodId.CENTRAL, MethodId.FORWARD, MethodId.BACKWARD], + success_checker=RobustConsistency(rtol=1e-2, atol=1e-8), + ) + + assert bool(derivative.df["success"].values[0]) + value = float(np.squeeze(derivative.series.values[0])) + h = 1e-6 + expected = ( + rosen(point + np.array([h, 0.0])) - rosen(point - np.array([h, 0.0])) + ) / (2 * h) + assert np.isclose(value, expected, rtol=1e-6) + + +def test_robust_consistency_warns_when_rejecting_step_sizes(): + """`RobustConsistency` should tell the user when it rejects a step size + that looked self-consistent on its own -- otherwise a legitimate-seeming + result could silently vanish from the blend without a trace.""" + results = [ + ComputerResult( + method_id="central", value=10.0, metadata={"size_absolute": 1.0} + ), + ComputerResult( + method_id="central", value=10.01, metadata={"size_absolute": 0.5} + ), + ComputerResult( + method_id="central", value=9.99, metadata={"size_absolute": 0.1} + ), + ComputerResult( + method_id="central", + value=500.0, + metadata={"size_absolute": 0.01}, + ), + ] + + checker = RobustConsistency() + with pytest.warns(UserWarning, match="1 step size"): + success, value = checker.method(FakeDirectionalDerivative(results)) + + assert success + assert np.isclose(value, np.mean([10.0, 10.01, 9.99])) + + class TestRejectOutliers: - """Unit tests for `Consistency._reject_outliers`, the order-independent - iterative outlier-rejection pass over step sizes' per-size means (added - to reject step sizes that are spuriously self-consistent, see - `test_consistency_rejects_rounding_noise_dominated_step_sizes` above).""" + """Unit tests for `RobustConsistency._reject_outliers`, the + order-independent iterative outlier-rejection pass over step sizes' + per-size means (added to reject step sizes that are spuriously + self-consistent, see + `test_robust_consistency_rejects_rounding_noise_dominated_step_sizes` + above).""" def test_below_min_trend_samples_keeps_everything(self): # An outlier (500.0) is present, but there are fewer candidates than # `min_trend_samples` -- too little data to estimate a spread, so no # trimming is attempted at all. - checker = Consistency(min_trend_samples=5) + checker = RobustConsistency(min_trend_samples=5) means = [10.0, 10.01, 500.0] assert checker._reject_outliers(means) == means def test_no_outliers_keeps_all(self): - checker = Consistency() + checker = RobustConsistency() means = [10.0, 10.01, 9.99, 10.02] assert checker._reject_outliers(means) == means def test_removes_single_outlier(self): - checker = Consistency() + checker = RobustConsistency() means = [10.0, 10.01, 9.99, 10.02, 500.0] trusted = checker._reject_outliers(means) assert trusted == [10.0, 10.01, 9.99, 10.02] @@ -107,7 +196,7 @@ def test_removes_single_outlier(self): def test_removes_multiple_outliers_iteratively(self): # Two outliers on opposite sides of the trustworthy cluster; both # must be dropped, one per iteration, worst-first. - checker = Consistency() + checker = RobustConsistency() means = [10.0, 10.01, 9.99, 10.02, 500.0, -500.0] trusted = checker._reject_outliers(means) assert trusted == [10.0, 10.01, 9.99, 10.02] @@ -115,17 +204,17 @@ def test_removes_multiple_outliers_iteratively(self): def test_order_independent(self): # Dropping is based on value, not position: shuffling the input # must not change which candidates survive. - checker = Consistency() + checker = RobustConsistency() means = [500.0, 10.0, 10.01, 9.99, 10.02] trusted = checker._reject_outliers(means) assert sorted(trusted) == [9.99, 10.0, 10.01, 10.02] def test_respects_trend_n_sigma(self): means = [10.0, 10.01, 9.99, 10.02, 500.0] - lenient_checker = Consistency(trend_n_sigma=1e6) + lenient_checker = RobustConsistency(trend_n_sigma=1e6) assert lenient_checker._reject_outliers(means) == means - strict_checker = Consistency(trend_n_sigma=5.0) + strict_checker = RobustConsistency(trend_n_sigma=5.0) assert strict_checker._reject_outliers(means) == [ 10.0, 10.01, @@ -138,7 +227,7 @@ def test_vector_valued_drops_whole_candidate_on_any_element_outlier(self): # outlier in another must still be dropped entirely (not just # masked in the bad element) -- "badness" is reduced across all # output dimensions before picking the worst candidate. - checker = Consistency() + checker = RobustConsistency() means = [ np.array([10.0, 5.0]), np.array([10.01, 5.01]), @@ -156,9 +245,9 @@ def test_nan_candidate_is_never_flagged_as_worst(self): # candidate whose mean is entirely NaN can never be selected as # "the worst" and is left in the trusted set untouched (harmless in # practice: it doesn't shift `np.nanmean` of the final value, and - # `Consistency.method`'s final blanket `isclose` check against a - # non-NaN blended value still reports `success=False` overall). - checker = Consistency() + # `RobustConsistency.method`'s final blanket `isclose` check against + # a non-NaN blended value still reports `success=False` overall). + checker = RobustConsistency() means = [10.0, 10.01, np.nan] trusted = checker._reject_outliers(means) assert len(trusted) == 3