feat: add NK Landscape

.gitignore

@@ -7,3 +7,4 @@ site/
 CLAUDE.md
 __pycache__/
 .DS_Store
+.gstack/

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

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

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

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**
 

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.

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",

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",
 ]

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",
 ]

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,133 @@ 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)
+        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).
+        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 _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:
+        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(self._contribution(solution, pos) for pos in self._affected[index])
+        solution[index] = 1 - solution[index]
+        try:
+            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
+
+
 class OneMax(BitstringProblem):
     """
     OneMax problem: maximize the number of 1-bits.

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,100 @@ 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 type(fit) is float
+            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")