diff --git a/CHANGELOG.md b/CHANGELOG.md index 130b46c..67af694 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,8 @@ ### Bug Fixes +- Continuous sampling now handles failed local searches explicitly. A local search fails when `minimize` raises `ValueError` or returns a non-finite fitness. A failed initial local search fails its run, and a failed search after a perturbation is rejected and counts as a non-improving iteration. Previously, runs starting at a `nan` minimum silently vanished, and an `inf` minimum made `LON.from_trace_data()` fail after all sampling had finished. Results with `success=False` and a finite fitness are still used as local optima. Each failed local search emits a warning with its context and reason (#44). +- `StepSizeEstimator` divides each sample's escape count by its number of evaluated perturbations instead of `n_perturbations`, so failed local searches no longer bias the escape rate down. Samples whose initial local search fails are skipped, and `RuntimeError` is raised if no sample can be evaluated, instead of returning `nan` (#44). - Runs without accepted transitions are no longer dropped from the LON and its performance metrics. `BasinHoppingSampler` and `ILSSampler` now record such a run as a single trace row from its initial local optimum to itself. `LON.from_trace_data()` drops self-loops, so the optimum becomes a vertex without edges and the run counts towards `success` and `deviation`. Basin-Hopping runs with a non-finite initial minimum are still skipped. Seeded traces and LONs change only when some run had no accepted transitions. Set `preserve_runs_without_transitions=False` in either sampler config to restore the earlier behavior, excluding these runs from both the graph and performance metrics (#39). - Raised the minimum joblib version to 1.4.0, which introduced the `generator_unordered` mode required by parallel Basin-Hopping sampling. diff --git a/docs/user-guide/sampling.md b/docs/user-guide/sampling.md index b254f96..c1652b8 100644 --- a/docs/user-guide/sampling.md +++ b/docs/user-guide/sampling.md @@ -168,6 +168,23 @@ config = BasinHoppingSamplerConfig( ) ``` +### Failed Local Searches + +A local search fails when `minimize` raises `ValueError` or returns a non-finite +fitness (`nan` or `inf`). Such a result never enters the LON: + +- if the initial local search of a run fails, the whole run fails and is left out + of the trace and the performance metrics; +- if a local search after a perturbation fails, it is rejected and counts as a + non-improving iteration. + +Results with `success=False` but a finite fitness are used as local optima. On +non-smooth functions such as Ackley, L-BFGS-B often reports `ABNORMAL` at genuine +optima. + +Each failed local search emits a warning with the run and iteration (for +perturbations), the starting point, and the reason for failure. + ## Discrete Sampling with ILS Independent ILS runs can use multiple processes via `n_jobs`: diff --git a/src/lonkit/continuous/sampling.py b/src/lonkit/continuous/sampling.py index ffa53b1..52f6fb9 100644 --- a/src/lonkit/continuous/sampling.py +++ b/src/lonkit/continuous/sampling.py @@ -6,7 +6,7 @@ import joblib import numpy as np import pandas as pd -from scipy.optimize import minimize +from scipy.optimize import OptimizeResult, minimize from tqdm import tqdm from lonkit.lon import LON, LONConfig @@ -53,6 +53,11 @@ class BasinHoppingSamplerConfig: preserve_runs_without_transitions: Keep runs with no accepted moves as self-loop trace rows. If False, omit them from the trace, LON, and performance metrics, matching the earlier behavior. Default: True. + + A local search fails if ``minimize`` raises ``ValueError`` or returns a non-finite + ``fun``. A failed initial search skips the run; a failed perturbation is rejected + and counts as a non-improving iteration. Failures emit warnings. Results with + ``success=False`` and finite fitness are still used as local optima. """ n_runs: int = 100 @@ -182,6 +187,34 @@ def _hash_solution(self, x: np.ndarray) -> str: return hash_str + def _local_search( + self, + func: Callable[[np.ndarray], float], + x0: np.ndarray, + bounds_array: np.ndarray | None, + ) -> tuple[OptimizeResult | None, str | None]: + """ + Run ``minimize`` from ``x0``. + + Returns: + Tuple of (result, failure). ``result`` is ``None`` if ``minimize`` raised + ``ValueError``. ``failure`` is ``None`` if the result can be used as a local + optimum, otherwise a description of why it cannot. + """ + try: + res = minimize( + func, + x0, + method=self.config.minimizer_method, + options=self.config.minimizer_options, + bounds=bounds_array if self.config.bounded else None, + ) + except ValueError as e: + return None, f"ValueError: {e}" + if not np.isfinite(res.fun): + return res, f"non-finite fitness {res.fun}" + return res, None + def _single_bh_run( self, run: int, @@ -206,26 +239,20 @@ def _single_bh_run( Tuple of (records, nfev) for this run. """ nfev = 0 - - try: - res = minimize( - func, - initial_point, - method=self.config.minimizer_method, - options=self.config.minimizer_options, - bounds=bounds_array if self.config.bounded else None, - ) - except ValueError as e: + res, failure = self._local_search(func, initial_point, bounds_array) + if res is not None: + nfev += res.nfev + if failure is not None: warnings.warn( - f"Run {run}: initial minimize failed with ValueError: {e}. " + f"Run {run}: initial local search failed ({failure}). " f"Starting point: {initial_point}. Skipping run.", stacklevel=3, ) - return [], 0 + return [], nfev + assert res is not None current_x = res.x current_f = res.fun - nfev += res.nfev records: list[dict] = [] iters_without_improvement = 0 @@ -241,24 +268,15 @@ def _single_bh_run( break x_perturbed = self._perturbation(current_x, p, rng, bounds_array) - try: - res = minimize( - func, - x_perturbed, - method=self.config.minimizer_method, - options=self.config.minimizer_options, - bounds=bounds_array if self.config.bounded else None, - ) - except ValueError as e: - # L-BFGS-B can produce internal iterates that slightly - # violate bounds, causing approx_derivative to fail. - # Skip this perturbation and try the next one. + res, failure = self._local_search(func, x_perturbed, bounds_array) + if failure is not None: warnings.warn( f"Run {run}, iteration {iter_index}: minimize after perturbation " - f"failed with ValueError: {e}. " + f"failed ({failure}). " f"Perturbed point: {x_perturbed}. Skipping perturbation.", stacklevel=3, ) + if res is None: iters_without_improvement += 1 iter_index += 1 continue @@ -266,6 +284,7 @@ def _single_bh_run( new_x = res.x new_f = res.fun nfev += res.nfev + accepted = failure is None and new_f <= current_f records.append( { @@ -275,26 +294,27 @@ def _single_bh_run( "current_f": current_f, "new_x": new_x.copy(), "new_f": new_f, - "accepted": new_f <= current_f, + "accepted": accepted, } ) + # A failed local search counts as a non-improving iteration if self.config.n_iter_no_change is not None: - if new_f < current_f: + if failure is None and new_f < current_f: iters_without_improvement = 0 else: iters_without_improvement += 1 # Acceptance criterion (minimization: accept if better or equal) - if new_f <= current_f: + if accepted: current_x = new_x.copy() current_f = new_f iter_index += 1 - # No perturbation produced a record: preserve the finite initial result + # No perturbation produced a record: preserve the initial result # for _construct_trace_data without inventing an accepted move. - if self.config.preserve_runs_without_transitions and not records and np.isfinite(current_f): + if self.config.preserve_runs_without_transitions and not records: records.append( { "run": run, @@ -363,6 +383,9 @@ def _parallel_bh( The ``joblib.Parallel`` runner needs to use the ``_run_single_bh_in_worker`` function, because of pickling - the pickled objects cannot share the same memory to achieve true parallelism between the processes, so we need to create separate ``BasinHoppingSampler`` object for each run. + + Warnings raised in worker processes do not reach the parent process, so workers + return them and they are re-emitted here. """ parallel_runner = joblib.Parallel( n_jobs=effective_n_jobs, @@ -388,7 +411,9 @@ def _parallel_bh( result_iter = ( tqdm(parallel_results, total=self.config.n_runs) if verbose else parallel_results ) - for completed_runs, (_, records, nfev) in enumerate(result_iter, 1): + for completed_runs, (_, records, nfev, caught) in enumerate(result_iter, 1): + for message in caught: + warnings.warn(message, stacklevel=2) if progress_callback: progress_callback(completed_runs, self.config.n_runs) raw_records.extend(records) @@ -512,13 +537,13 @@ def _construct_trace_data(self, raw_records: list[dict]) -> pd.DataFrame: if self.config.preserve_runs_without_transitions: # A run with no accepted moves keeps its initial minimum as a self-loop row # (dropped from the graph by LON.from_trace_data()), so it still counts - # towards success/deviation. Non-finite minima are skipped. + # towards success/deviation. runs_in_trace = {row["run"] for row in trace_records} first_records: dict[int, dict] = {} for rec in raw_records: first_records.setdefault(rec["run"], rec) for run, rec in first_records.items(): - if run in runs_in_trace or not np.isfinite(rec["current_f"]): + if run in runs_in_trace: continue node = self._hash_solution( self._round_value(rec["current_x"], self.config.coordinate_precision) @@ -650,7 +675,7 @@ def _run_single_bh_in_worker( bounds_array: np.ndarray | None, config: BasinHoppingSamplerConfig, seed: np.random.SeedSequence, -) -> tuple[int, list[dict], int]: +) -> tuple[int, list[dict], int, list[Warning | str]]: """ Module-level worker for a single Basin-Hopping run. @@ -669,12 +694,15 @@ def _run_single_bh_in_worker( seed: ``SeedSequence`` child spawned for this run. Returns: - Tuple of (run, records, nfev) for this run. + Tuple of (run, records, nfev, warnings) for this run. Warnings are returned + instead of emitted, because they would not reach the parent process. """ sampler = BasinHoppingSampler(config) rng = np.random.default_rng(seed) - records, nfev = sampler._single_bh_run(run, func, initial_point, p, bounds_array, rng) - return run, records, nfev + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + records, nfev = sampler._single_bh_run(run, func, initial_point, p, bounds_array, rng) + return run, records, nfev, [w.message for w in caught] def compute_lon( diff --git a/src/lonkit/continuous/step_size.py b/src/lonkit/continuous/step_size.py index 08794ab..24cdc61 100644 --- a/src/lonkit/continuous/step_size.py +++ b/src/lonkit/continuous/step_size.py @@ -3,7 +3,6 @@ from dataclasses import dataclass import numpy as np -from scipy.optimize import minimize from lonkit.continuous.sampling import BasinHoppingSampler, BasinHoppingSamplerConfig @@ -33,6 +32,10 @@ class StepSizeEstimatorConfig: Default: `None`. bounded: Whether to enforce domain bounds during perturbation. Default: `True`. seed: Random seed for reproducibility. Default: `None`. + + Failed local searches (see ``BasinHoppingSamplerConfig``) are skipped: a sample whose + initial local search fails is dropped, and failed perturbations are excluded from + that sample's escape rate. """ n_samples: int = 100 @@ -125,6 +128,13 @@ def _compute_escape_rate( Compute the average escape rate for a given step size. Cost per call: `n_samples` + `n_samples * n_perturbations` minimizations. + + Each sample's escape rate is the fraction of escapes among its evaluated + perturbations, i.e. those whose local search did not fail. + + Raises: + RuntimeError: If no sample has a valid initial local search and at least + one evaluated perturbation. """ bounds_array = domain_array if self.config.bounded else None p = step_size * np.abs(domain_array[:, 1] - domain_array[:, 0]) @@ -134,53 +144,46 @@ def _compute_escape_rate( for _ in range(self.config.n_samples): x0 = rng.uniform(domain_array[:, 0], domain_array[:, 1]) - try: - res = minimize( - func, - x0, - method=self.config.minimizer_method, - options=self.config.minimizer_options, - bounds=bounds_array if self.config.bounded else None, - ) - except ValueError as e: + res, failure = sampler._local_search(func, x0, bounds_array) + if failure is not None: warnings.warn( - f"Initial minimize failed with ValueError: {e}. " - f"Starting point: {x0}. Skipping run.", + f"Initial local search failed ({failure}). " + f"Starting point: {x0}. Skipping sample.", stacklevel=3, ) continue + assert res is not None optimum = res.x optimum_rounded = sampler._round_value(optimum, self.config.coordinate_precision) optimum_hash = sampler._hash_solution(optimum_rounded) escapes = 0 + evaluated = 0 for _ in range(self.config.n_perturbations): x_perturbed = sampler._perturbation(optimum, p, rng, bounds_array) - try: - res_perturbed = minimize( - func, - x_perturbed, - method=self.config.minimizer_method, - options=self.config.minimizer_options, - bounds=bounds_array, - ) - except ValueError as e: - # L-BFGS-B can produce internal iterates that slightly - # violate bounds, causing approx_derivative to fail. - # Skip this perturbation and try the next one. + res_perturbed, failure = sampler._local_search(func, x_perturbed, bounds_array) + if failure is not None: + # Excluded from this sample's escape rate warnings.warn( - f"Minimize after perturbation " - f"failed with ValueError: {e}. " + f"Minimize after perturbation failed ({failure}). " f"Perturbed point: {x_perturbed}. Skipping perturbation.", stacklevel=3, ) continue + assert res_perturbed is not None + evaluated += 1 new_hash = sampler._hash_solution(res_perturbed.x) if new_hash != optimum_hash: escapes += 1 - escape_rates.append(escapes / self.config.n_perturbations) + if evaluated: + escape_rates.append(escapes / evaluated) + if not escape_rates: + raise RuntimeError( + f"Could not estimate the escape rate for step size {step_size}: " + "every sample failed its initial local search or all of its perturbations." + ) return float(np.mean(escape_rates)) def estimate( diff --git a/tests/continuous/test_sampling.py b/tests/continuous/test_sampling.py index c8c2360..528f853 100644 --- a/tests/continuous/test_sampling.py +++ b/tests/continuous/test_sampling.py @@ -1,6 +1,7 @@ from unittest.mock import patch import numpy as np +import pandas as pd import pytest from scipy.optimize import OptimizeResult @@ -71,27 +72,27 @@ def test_mixed_runs_keep_accepted_transitions(self, preserve: bool) -> None: } @pytest.mark.parametrize("preserve", [True, False]) - @pytest.mark.parametrize("initial_fitness", [0.0, float("nan")]) - def test_all_perturbations_fail(self, initial_fitness: float, preserve: bool) -> None: + def test_all_perturbations_fail(self, preserve: bool) -> None: sampler = BasinHoppingSampler( BasinHoppingSamplerConfig( n_runs=1, max_iter=2, seed=SEED, preserve_runs_without_transitions=preserve ) ) - initial_result = OptimizeResult(x=np.zeros(2), fun=initial_fitness, nfev=1) + initial_result = OptimizeResult(x=np.zeros(2), fun=0.0, nfev=1) with ( patch( "lonkit.continuous.sampling.minimize", side_effect=[initial_result, ValueError("failed"), ValueError("failed")], ), - pytest.warns(UserWarning, match="minimize after perturbation failed") as warnings, + pytest.warns(UserWarning) as warnings, ): result = sampler.sample(sphere, DOMAIN_2D, initial_points=np.zeros((1, 2))) - assert len(warnings) == 2 + messages = [str(w.message) for w in warnings] + assert sum("minimize after perturbation failed" in m for m in messages) == 2 assert result.nfev == 1 lon = sampler.sample_to_lon(result) - if not preserve or not np.isfinite(initial_fitness): + if not preserve: assert result.raw_records == [] assert result.trace_df.empty assert lon.n_vertices == 0 @@ -144,3 +145,107 @@ def nan_for_negative_x0(x: np.ndarray) -> float: nan_for_negative_x0, DOMAIN_2D, initial_points=np.array([[-2.0, 0.0], [2.0, 0.0]]) ) assert set(result.trace_df["run"]) == {2} + + +def _res(x: float, fun: float, success: bool = True) -> OptimizeResult: + return OptimizeResult( + x=np.full(2, x), + fun=fun, + nfev=1, + success=success, + status=0 if success else 2, + message="ok" if success else "ABNORMAL", + ) + + +class TestMinimizerFailures: + def test_unconverged_finite_result_is_used(self) -> None: + sampler = BasinHoppingSampler(BasinHoppingSamplerConfig(n_runs=1, max_iter=1, seed=SEED)) + results = [_res(2.0, 8.0), _res(1.0, 2.0, success=False)] + with patch("lonkit.continuous.sampling.minimize", side_effect=results): + result = sampler.sample(sphere, DOMAIN_2D, initial_points=np.zeros((1, 2))) + + (record,) = result.raw_records + assert record["accepted"] + assert result.nfev == 2 + assert sampler.sample_to_lon(result).n_vertices == 2 + + @pytest.mark.parametrize("bad", [float("nan"), float("inf")]) + def test_non_finite_initial_minimum_fails_run(self, bad: float) -> None: + def f(x: np.ndarray) -> float: + return bad if x[0] > 3 else float(np.sum(x**2)) + + sampler = BasinHoppingSampler( + BasinHoppingSamplerConfig(n_runs=4, n_iter_no_change=5, seed=1) + ) + points = np.array([[4.0, 1.0], [5.0, 2.0], [-1.0, -1.0], [-2.0, 1.0]]) + with pytest.warns(UserWarning, match="initial local search failed"): + result = sampler.sample(f, [(-5, 5)] * 2, initial_points=points) + + assert set(result.trace_df["run"]) == {3, 4} + assert np.isfinite(result.trace_df[["fit1", "fit2"]].to_numpy()).all() + assert all(r["run"] in (3, 4) for r in result.raw_records) + sampler.sample_to_lon(result) # trace is valid + + def test_non_finite_mid_run_is_never_accepted(self) -> None: + sampler = BasinHoppingSampler(BasinHoppingSamplerConfig(n_runs=1, max_iter=2, seed=SEED)) + results = [_res(2.0, 8.0), _res(1.0, float("-inf")), _res(1.0, 2.0)] + with ( + patch("lonkit.continuous.sampling.minimize", side_effect=results), + pytest.warns(UserWarning, match="minimize after perturbation failed"), + ): + result = sampler.sample(sphere, DOMAIN_2D, initial_points=np.zeros((1, 2))) + + assert [r["accepted"] for r in result.raw_records] == [False, True] + assert result.trace_df[["fit1", "fit2"]].to_numpy().tolist() == [[8.0, 2.0]] + + def test_initial_value_error_reports_failed_run(self) -> None: + sampler = BasinHoppingSampler(BasinHoppingSamplerConfig(n_runs=2, max_iter=1, seed=SEED)) + results = [ValueError("bad"), _res(2.0, 8.0), _res(1.0, 2.0)] + with ( + patch("lonkit.continuous.sampling.minimize", side_effect=results), + pytest.warns(UserWarning, match="Run 1: initial local search failed"), + ): + result = sampler.sample(sphere, DOMAIN_2D, initial_points=np.zeros((2, 2))) + assert result.trace_df["run"].tolist() == [2] + + @pytest.mark.parametrize("bad", [float("nan"), float("inf"), float("-inf")]) + def test_failed_perturbations_stop_after_no_change_limit(self, bad: float) -> None: + sampler = BasinHoppingSampler( + BasinHoppingSamplerConfig(n_runs=1, n_iter_no_change=2, max_iter=None, seed=SEED) + ) + results = [_res(2.0, 8.0), _res(1.0, bad), _res(1.0, bad)] + with ( + patch("lonkit.continuous.sampling.minimize", side_effect=results) as minimize, + pytest.warns(UserWarning, match="minimize after perturbation failed") as warnings, + ): + result = sampler.sample(sphere, DOMAIN_2D, initial_points=np.zeros((1, 2))) + + assert minimize.call_count == 3 + assert len(warnings) == 2 + assert not any(r["accepted"] for r in result.raw_records) + assert result.nfev == 3 + assert result.trace_df["fit1"].tolist() == [8.0] + assert sampler.sample_to_lon(result).n_vertices == 1 + + def test_failed_runs_match_across_n_jobs(self) -> None: + def f(x: np.ndarray) -> float: + return float("nan") if x[0] > 3 else float(np.sum(x**2)) + + results = [] + messages = [] + points = np.array([[4.0, 1.0], [-1.0, -1.0]]) + for n_jobs in (1, 2): + config = BasinHoppingSamplerConfig( + n_runs=2, n_iter_no_change=3, seed=SEED, n_jobs=n_jobs + ) + with pytest.warns(UserWarning) as caught: + results.append( + BasinHoppingSampler(config).sample(f, DOMAIN_2D, initial_points=points) + ) + messages.append(sorted(str(w.message) for w in caught)) + assert any("Run 1: initial local search failed" in m for m in messages[1]) + assert messages[0] == messages[1] + pd.testing.assert_frame_equal(results[0].trace_df, results[1].trace_df) + assert results[0].nfev == results[1].nfev + assert set(results[0].trace_df["run"]) == {2} diff --git a/tests/continuous/test_step_size.py b/tests/continuous/test_step_size.py new file mode 100644 index 0000000..d215474 --- /dev/null +++ b/tests/continuous/test_step_size.py @@ -0,0 +1,56 @@ +from unittest.mock import patch + +import numpy as np +import pytest +from scipy.optimize import OptimizeResult + +from lonkit import StepSizeEstimator, StepSizeEstimatorConfig +from tests.conftest import DOMAIN_2D, SEED, sphere + + +def _res(x: float, fun: float = 0.0, success: bool = True) -> OptimizeResult: + return OptimizeResult(x=np.full(2, x), fun=fun, nfev=1, success=success, message="") + + +class TestEscapeRateWithFailures: + def test_failed_perturbations_are_excluded_from_denominator(self) -> None: + estimator = StepSizeEstimator( + StepSizeEstimatorConfig(n_samples=1, n_perturbations=4, seed=SEED) + ) + # Two of four perturbations fail; one of the two evaluated escapes. + results = [_res(0.0), ValueError("bad"), _res(1.0, float("nan")), _res(0.0), _res(1.0)] + with ( + patch("lonkit.continuous.sampling.minimize", side_effect=results), + pytest.warns(UserWarning, match="Skipping perturbation"), + ): + rate = estimator._compute_escape_rate( + sphere, np.array(DOMAIN_2D), 0.1, estimator._make_sampler() + ) + assert rate == 0.5 + + def test_sample_without_evaluated_perturbations_is_skipped(self) -> None: + estimator = StepSizeEstimator( + StepSizeEstimatorConfig(n_samples=2, n_perturbations=1, seed=SEED) + ) + results = [_res(0.0), ValueError("bad"), _res(0.0), _res(1.0)] + with ( + patch("lonkit.continuous.sampling.minimize", side_effect=results), + pytest.warns(UserWarning), + ): + rate = estimator._compute_escape_rate( + sphere, np.array(DOMAIN_2D), 0.1, estimator._make_sampler() + ) + assert rate == 1.0 + + def test_all_samples_failing_raises(self) -> None: + def nan_everywhere(x: np.ndarray) -> float: + return float("nan") + + estimator = StepSizeEstimator( + StepSizeEstimatorConfig(n_samples=3, n_perturbations=2, seed=SEED) + ) + with ( + pytest.warns(UserWarning, match="Initial local search failed"), + pytest.raises(RuntimeError, match="Could not estimate the escape rate"), + ): + estimator.estimate(nan_everywhere, DOMAIN_2D)