From 25106135b32571c41e3b303aa31dbe8de36f5995 Mon Sep 17 00:00:00 2001 From: WojtAcht Date: Mon, 8 Jun 2026 00:52:46 +0200 Subject: [PATCH 1/2] feat: add NK Landscape --- .gitignore | 1 + README.md | 2 +- docs/api/discrete.md | 5 + docs/api/index.md | 1 + docs/index.md | 2 +- docs/user-guide/discrete.md | 23 ++++ src/lonkit/__init__.py | 9 +- src/lonkit/discrete/__init__.py | 9 +- src/lonkit/discrete/problems/__init__.py | 2 + src/lonkit/discrete/problems/bitstring.py | 125 ++++++++++++++++++++++ tests/discrete/test_problems.py | 98 +++++++++++++++++ 11 files changed, 273 insertions(+), 4 deletions(-) diff --git a/.gitignore b/.gitignore index 74bbabe..d5ff491 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,4 @@ site/ CLAUDE.md __pycache__/ .DS_Store +.gstack/ diff --git a/README.md b/README.md index 9ebdb79..d7b7a05 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ lonkit is a Python library for constructing, analyzing, and visualizing Local Op ## Features - **Basin-Hopping Sampling**: Efficient exploration of continuous fitness landscapes using configurable Basin-Hopping -- **Iterated Local Search**: Discrete LON construction via ILS with built-in problems (Number Partitioning, OneMax) +- **Iterated Local Search**: Discrete LON construction via ILS with built-in problems (Number Partitioning, NK Landscape, OneMax) - **LON Construction**: Automatic construction of Local Optima Networks from sampling data - **CMLON Support**: Compressed Monotonic LONs for cleaner landscape analysis - **Rich Metrics**: Compute landscape metrics including funnel analysis and neutrality diff --git a/docs/api/discrete.md b/docs/api/discrete.md index 3a22412..367e918 100644 --- a/docs/api/discrete.md +++ b/docs/api/discrete.md @@ -17,6 +17,11 @@ show_root_heading: true show_source: true +::: lonkit.discrete.problems.bitstring.NKLandscape + options: + show_root_heading: true + show_source: true + ::: lonkit.discrete.problems.bitstring.OneMax options: show_root_heading: true diff --git a/docs/api/index.md b/docs/api/index.md index 2693b7e..6e04f04 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -36,6 +36,7 @@ Iterated Local Search sampling and built-in discrete problems. - [`DiscreteProblem`](discrete.md#lonkit.discrete.problems.problem.DiscreteProblem) - Abstract base for discrete problems - [`BitstringProblem`](discrete.md#lonkit.discrete.problems.bitstring.BitstringProblem) - Base class for bitstring problems - [`NumberPartitioning`](discrete.md#lonkit.discrete.problems.bitstring.NumberPartitioning) - Number Partitioning Problem +- [`NKLandscape`](discrete.md#lonkit.discrete.problems.bitstring.NKLandscape) - Kauffman's NK Landscape benchmark - [`OneMax`](discrete.md#lonkit.discrete.problems.bitstring.OneMax) - OneMax benchmark - [`ILSSampler`](discrete.md#lonkit.discrete.sampling.ILSSampler) - ILS sampling class - [`ILSSamplerConfig`](discrete.md#lonkit.discrete.sampling.ILSSamplerConfig) - ILS configuration diff --git a/docs/index.md b/docs/index.md index 7ffaaf6..0226300 100644 --- a/docs/index.md +++ b/docs/index.md @@ -32,7 +32,7 @@ Local Optima Networks (LONs) are graph-based models that capture the global stru --- - Discrete LON construction via ILS with built-in problems (Number Partitioning, OneMax) and custom problem support + Discrete LON construction via ILS with built-in problems (Number Partitioning, NK Landscape, OneMax) and custom problem support - **LON Construction** diff --git a/docs/user-guide/discrete.md b/docs/user-guide/discrete.md index cd212d1..289b862 100644 --- a/docs/user-guide/discrete.md +++ b/docs/user-guide/discrete.md @@ -48,6 +48,29 @@ problem = NumberPartitioning(n=4, weights=[7, 5, 6, 4]) **Fitness:** `|sum(A) - sum(B)|` (minimization, optimal = 0) +### NKLandscape + +**NKLandscape** is Kauffman's NK benchmark: maximize the average contribution +of `N` bit-dependent fitness components, where each component depends on its own +bit and `K` interacting bits. + +```python +from lonkit import NKLandscape + +problem = NKLandscape(n=20, k=4, instance_seed=1) +``` + +**Parameters:** + +- `n`: Length of the bitstring +- `k`: Number of epistatic interactions per position. Must be in `[0, n - 1]`. +- `instance_seed`: Seed for generating the fixed landscape instance +- `neighbor_model`: `"adjacent"` for cyclic adjacent interactions or `"random"` for random distinct neighbors (default: `"adjacent"`) +- `n_perturbation_flips`: Number of random bit flips per perturbation (default: 2) +- `first_improvement`: Use first-improvement hill climbing (default: True) + +**Fitness:** Average contribution in `[0, 1)` (maximization) + ### OneMax **OneMax** is a simple benchmark: maximize the number of 1-bits in a bitstring. diff --git a/src/lonkit/__init__.py b/src/lonkit/__init__.py index d71eca9..4ad4dac 100644 --- a/src/lonkit/__init__.py +++ b/src/lonkit/__init__.py @@ -5,7 +5,13 @@ compute_lon, ) from lonkit.continuous.step_size import StepSizeEstimator, StepSizeEstimatorConfig, StepSizeResult -from lonkit.discrete.problems import BitstringProblem, DiscreteProblem, NumberPartitioning, OneMax +from lonkit.discrete.problems import ( + BitstringProblem, + DiscreteProblem, + NKLandscape, + NumberPartitioning, + OneMax, +) from lonkit.discrete.sampling import ILSResult, ILSSampler, ILSSamplerConfig from lonkit.lon import CMLON, LON, LONConfig from lonkit.visualization import LONVisualizer @@ -24,6 +30,7 @@ "ILSSamplerConfig", "LONConfig", "LONVisualizer", + "NKLandscape", "NumberPartitioning", "OneMax", "StepSizeEstimator", diff --git a/src/lonkit/discrete/__init__.py b/src/lonkit/discrete/__init__.py index 1e3ac5c..cb0738a 100644 --- a/src/lonkit/discrete/__init__.py +++ b/src/lonkit/discrete/__init__.py @@ -1,4 +1,10 @@ -from lonkit.discrete.problems import BitstringProblem, DiscreteProblem, NumberPartitioning, OneMax +from lonkit.discrete.problems import ( + BitstringProblem, + DiscreteProblem, + NKLandscape, + NumberPartitioning, + OneMax, +) from lonkit.discrete.sampling import ILSResult, ILSSampler, ILSSamplerConfig __all__ = [ @@ -7,6 +13,7 @@ "ILSResult", "ILSSampler", "ILSSamplerConfig", + "NKLandscape", "NumberPartitioning", "OneMax", ] diff --git a/src/lonkit/discrete/problems/__init__.py b/src/lonkit/discrete/problems/__init__.py index 881a894..8d7a2bc 100644 --- a/src/lonkit/discrete/problems/__init__.py +++ b/src/lonkit/discrete/problems/__init__.py @@ -1,5 +1,6 @@ from lonkit.discrete.problems.bitstring import ( BitstringProblem, + NKLandscape, NumberPartitioning, OneMax, ) @@ -8,6 +9,7 @@ __all__ = [ "BitstringProblem", "DiscreteProblem", + "NKLandscape", "NumberPartitioning", "OneMax", ] diff --git a/src/lonkit/discrete/problems/bitstring.py b/src/lonkit/discrete/problems/bitstring.py index db7b8a3..fda9759 100644 --- a/src/lonkit/discrete/problems/bitstring.py +++ b/src/lonkit/discrete/problems/bitstring.py @@ -1,5 +1,6 @@ import math import warnings +from typing import Literal import numpy as np @@ -225,6 +226,130 @@ def evaluate(self, solution: list[int]) -> float: return float(abs(cost_a - cost_b)) +class NKLandscape(BitstringProblem): + """ + Kauffman's NK Landscape. + + A tunable family of rugged fitness landscapes over bitstrings of length + ``N``. Each of the ``N`` positions contributes a fitness component that + depends on its own bit plus ``K`` other ("epistatic") bits. The overall + fitness is the average of the ``N`` contributions:: + + fitness(x) = (1 / N) * sum_i f_i(x_i, x_{neighbors(i)}) + + Each component ``f_i`` is defined by a lookup table of ``2^(K+1)`` values + drawn uniformly from ``[0, 1)`` and indexed by the ``(K + 1)``-bit pattern + formed by bit ``i`` followed by its ``K`` neighbors. + + ``K`` tunes the ruggedness: + - ``K = 0`` gives a smooth, single-optimum landscape (each bit independent). + - ``K = N - 1`` gives a maximally rugged, random landscape. + + This is a **maximization** problem (optimal fitness close to 1.0). + + The instance (neighbor structure and contribution tables) is fixed at + construction time from ``instance_seed``, so two problems built with the + same parameters are identical. + + Args: + n: Length of the bitstring. Must be > 0. + k: Number of epistatic interactions per position. Must be in [0, n-1]. + instance_seed: Seed for generating the neighbor structure (random model) + and the random contribution tables. Required for reproducibility. + neighbor_model: ``"adjacent"`` (each position interacts with its ``K`` + cyclically-following neighbors) or ``"random"`` (``K`` distinct + positions drawn uniformly at random, excluding the position itself). + Default: ``"adjacent"``. + n_perturbation_flips: Number of random flips per perturbation (default: 2). + first_improvement: If True, local search uses first-improvement + hill climbing (stochastic -- scan order randomized each pass). + If False, uses best-improvement (deterministic). Default: True. + """ + + @property + def minimize(self) -> bool: + return False + + def __init__( + self, + n: int, + k: int, + instance_seed: int, + neighbor_model: Literal["adjacent", "random"] = "adjacent", + n_perturbation_flips: int = 2, + first_improvement: bool = True, + ): + super().__init__(n, n_perturbation_flips, first_improvement) + if instance_seed is None: + raise ValueError("instance_seed is required for NKLandscape") + if k < 0 or k > n - 1: + raise ValueError(f"k must be in [0, {n - 1}], got {k}") + if neighbor_model not in ("adjacent", "random"): + raise ValueError( + f"neighbor_model must be 'adjacent' or 'random', got {neighbor_model!r}" + ) + self.k = k + self.instance_seed = instance_seed + self.neighbor_model = neighbor_model + + rng = np.random.default_rng(instance_seed) + self.neighbors = self._build_neighbors(rng) + # Contribution tables: one row of 2^(k+1) random values per position. + self.tables = rng.random(size=(n, 1 << (k + 1))) + # For each bit, which positions' contributions depend on it (for delta). + self._affected: list[list[int]] = [[] for _ in range(n)] + for pos in range(n): + self._affected[pos].append(pos) + for j in self.neighbors[pos]: + self._affected[j].append(pos) + + def _build_neighbors(self, rng: np.random.Generator) -> list[list[int]]: + """Return the list of K epistatic neighbors for each position.""" + if self.neighbor_model == "adjacent": + return [ + [(i + offset) % self.n for offset in range(1, self.k + 1)] for i in range(self.n) + ] + # random model: K distinct positions, excluding i itself + neighbors: list[list[int]] = [] + for i in range(self.n): + candidates = [j for j in range(self.n) if j != i] + chosen = rng.choice(candidates, size=self.k, replace=False) + neighbors.append(sorted(int(j) for j in chosen)) + return neighbors + + def _contribution_index(self, solution: list[int], pos: int) -> int: + """Build the table index from bit `pos` followed by its neighbors.""" + idx = solution[pos] + for j in self.neighbors[pos]: + idx = (idx << 1) | solution[j] + return idx + + def evaluate(self, solution: list[int]) -> float: + total = 0.0 + for pos in range(self.n): + total += self.tables[pos][self._contribution_index(solution, pos)] + return total / self.n + + def delta_evaluate(self, solution: list[int], index: int) -> float | None: + """ + Delta evaluation: flipping bit `index` only changes the + contributions of positions that depend on it. + """ + old_total = sum( + float(self.tables[pos][self._contribution_index(solution, pos)]) + for pos in self._affected[index] + ) + solution[index] = 1 - solution[index] + try: + new_total = sum( + float(self.tables[pos][self._contribution_index(solution, pos)]) + for pos in self._affected[index] + ) + finally: + solution[index] = 1 - solution[index] + return (new_total - old_total) / self.n + + class OneMax(BitstringProblem): """ OneMax problem: maximize the number of 1-bits. diff --git a/tests/discrete/test_problems.py b/tests/discrete/test_problems.py index d528ca1..feafc5e 100644 --- a/tests/discrete/test_problems.py +++ b/tests/discrete/test_problems.py @@ -1,6 +1,8 @@ import numpy as np +import pytest from lonkit.discrete.problems.bitstring import ( + NKLandscape, NumberPartitioning, OneMax, ) @@ -71,3 +73,99 @@ def test_compare_minimization(self): assert p.compare(0, 5) == 1 assert p.compare(5, 0) == -1 assert p.compare(3, 3) == 0 + + +class TestNKLandscape: + def test_minimize_is_false(self): + p = NKLandscape(n=8, k=2, instance_seed=0) + assert p.minimize is False + + def test_fitness_in_unit_interval(self): + p = NKLandscape(n=10, k=3, instance_seed=1) + rng = np.random.default_rng(7) + for _ in range(50): + sol = p.random_solution(rng) + fit = p.evaluate(sol) + assert 0.0 <= fit < 1.0 + + def test_evaluate_is_deterministic(self): + p = NKLandscape(n=8, k=2, instance_seed=3) + sol = [0, 1, 1, 0, 1, 0, 0, 1] + assert p.evaluate(sol) == p.evaluate(sol) + + def test_same_seed_same_instance(self): + p1 = NKLandscape(n=8, k=2, instance_seed=42) + p2 = NKLandscape(n=8, k=2, instance_seed=42) + rng = np.random.default_rng(0) + for _ in range(20): + sol = p1.random_solution(rng) + assert p1.evaluate(sol) == p2.evaluate(sol) + + def test_different_seed_different_instance(self): + p1 = NKLandscape(n=8, k=2, instance_seed=1) + p2 = NKLandscape(n=8, k=2, instance_seed=2) + sol = [0, 1, 0, 1, 0, 1, 0, 1] + assert p1.evaluate(sol) != p2.evaluate(sol) + + def test_adjacent_neighbors_structure(self): + p = NKLandscape(n=5, k=2, instance_seed=0, neighbor_model="adjacent") + assert p.neighbors[0] == [1, 2] + assert p.neighbors[4] == [0, 1] # cyclic wrap-around + + def test_random_neighbors_excludes_self_and_distinct(self): + p = NKLandscape(n=8, k=3, instance_seed=5, neighbor_model="random") + for i in range(8): + assert i not in p.neighbors[i] + assert len(set(p.neighbors[i])) == 3 + + @pytest.mark.parametrize("neighbor_model", ["adjacent", "random"]) + def test_delta_evaluate_matches_full_evaluation(self, neighbor_model): + p = NKLandscape(n=10, k=3, instance_seed=11, neighbor_model=neighbor_model) + rng = np.random.default_rng(2) + for _ in range(20): + sol = p.random_solution(rng) + base = p.evaluate(sol) + for i in range(p.n): + delta = p.delta_evaluate(sol, i) + flipped = list(sol) + flipped[i] = 1 - flipped[i] + expected = p.evaluate(flipped) - base + assert delta == pytest.approx(expected, abs=1e-12) + + def test_delta_evaluate_does_not_modify_solution(self): + p = NKLandscape(n=6, k=2, instance_seed=4) + sol = [1, 0, 1, 0, 1, 0] + original = list(sol) + p.delta_evaluate(sol, 3) + assert sol == original + + def test_local_search_returns_local_optimum(self): + p = NKLandscape(n=10, k=2, instance_seed=9) + rng = np.random.default_rng(3) + sol, fit = p.local_search(p.random_solution(rng), rng) + # No single bit flip should improve a local optimum. + for i in range(p.n): + assert not p.is_better(fit + p.delta_evaluate(sol, i), fit) + + def test_k_zero_is_smooth(self): + # With k=0 each bit is independent; greedy local search reaches the + # global optimum regardless of starting point. + p = NKLandscape(n=8, k=0, instance_seed=6) + rng = np.random.default_rng(1) + _, fit_a = p.local_search([0] * 8, rng) + _, fit_b = p.local_search([1] * 8, rng) + assert fit_a == pytest.approx(fit_b) + + def test_invalid_k_raises(self): + with pytest.raises(ValueError, match="k must be in"): + NKLandscape(n=4, k=4, instance_seed=0) + with pytest.raises(ValueError, match="k must be in"): + NKLandscape(n=4, k=-1, instance_seed=0) + + def test_instance_seed_is_required(self): + with pytest.raises(ValueError, match="instance_seed is required"): + NKLandscape(n=4, k=1, instance_seed=None) # type: ignore[arg-type] + + def test_invalid_neighbor_model_raises(self): + with pytest.raises(ValueError, match="neighbor_model"): + NKLandscape(n=4, k=1, instance_seed=0, neighbor_model="circular") From 5dad3af79da2475f4cc106c6350e02b5d2aa3b51 Mon Sep 17 00:00:00 2001 From: WojtAcht Date: Mon, 8 Jun 2026 01:07:21 +0200 Subject: [PATCH 2/2] feat: vectorize `evaluate` method --- src/lonkit/discrete/problems/bitstring.py | 27 +++++++++++++---------- tests/discrete/test_problems.py | 1 + 2 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/lonkit/discrete/problems/bitstring.py b/src/lonkit/discrete/problems/bitstring.py index fda9759..b6911ec 100644 --- a/src/lonkit/discrete/problems/bitstring.py +++ b/src/lonkit/discrete/problems/bitstring.py @@ -294,6 +294,12 @@ def __init__( rng = np.random.default_rng(instance_seed) self.neighbors = self._build_neighbors(rng) + self._dependencies = np.asarray( + [[pos, *neighbors] for pos, neighbors in enumerate(self.neighbors)], + dtype=np.intp, + ) + self._index_weights = 1 << np.arange(k, -1, -1, dtype=np.int64) + self._table_rows = np.arange(n, dtype=np.intp) # Contribution tables: one row of 2^(k+1) random values per position. self.tables = rng.random(size=(n, 1 << (k + 1))) # For each bit, which positions' contributions depend on it (for delta). @@ -324,27 +330,24 @@ def _contribution_index(self, solution: list[int], pos: int) -> int: idx = (idx << 1) | solution[j] return idx + def _contribution(self, solution: list[int], pos: int) -> float: + """Return the Python float contribution for position `pos`.""" + return float(self.tables[pos][self._contribution_index(solution, pos)]) + def evaluate(self, solution: list[int]) -> float: - total = 0.0 - for pos in range(self.n): - total += self.tables[pos][self._contribution_index(solution, pos)] - return total / self.n + bits = np.asarray(solution, dtype=np.int64)[self._dependencies] + indices = bits @ self._index_weights + return float(self.tables[self._table_rows, indices].mean()) def delta_evaluate(self, solution: list[int], index: int) -> float | None: """ Delta evaluation: flipping bit `index` only changes the contributions of positions that depend on it. """ - old_total = sum( - float(self.tables[pos][self._contribution_index(solution, pos)]) - for pos in self._affected[index] - ) + old_total = sum(self._contribution(solution, pos) for pos in self._affected[index]) solution[index] = 1 - solution[index] try: - new_total = sum( - float(self.tables[pos][self._contribution_index(solution, pos)]) - for pos in self._affected[index] - ) + new_total = sum(self._contribution(solution, pos) for pos in self._affected[index]) finally: solution[index] = 1 - solution[index] return (new_total - old_total) / self.n diff --git a/tests/discrete/test_problems.py b/tests/discrete/test_problems.py index feafc5e..ce9bf90 100644 --- a/tests/discrete/test_problems.py +++ b/tests/discrete/test_problems.py @@ -86,6 +86,7 @@ def test_fitness_in_unit_interval(self): for _ in range(50): sol = p.random_solution(rng) fit = p.evaluate(sol) + assert type(fit) is float assert 0.0 <= fit < 1.0 def test_evaluate_is_deterministic(self):