Skip to main content
NK landscapes, introduced by Stuart Kauffman, are a standard tool in complex systems research for studying how epistatic interactions between components affect the shape of a fitness landscape. A landscape with higher K has more interdependencies, creating a more rugged surface with many local optima that agents can get trapped in. The OptimizedNKLandscape class in examples/LandscapeConstruction_NK_Modular.py extends the classical NK model with a two-module structure and a parameter R that controls how insular each module’s dependencies are.

​
Parameters

ParameterTypeDescription
NintNumber of bits in each agent’s state. Must be even.
KintNumber of epistatic dependencies per bit. Must satisfy 0 <= K < N.
RfloatFraction of intra-module dependencies for module 1 (0–1). R=0 means all cross-module; R=1 means fully self-contained.
seedintOptional random seed for reproducibility.
use_cacheboolEnable LRU fitness caching (default True).
max_cache_sizeintMaximum number of cached fitness values (default 100000).
N must be even because the landscape is divided into exactly 2 modules of N/2 bits each. Passing an odd N raises a ValueError.

​
Module structure

The landscape is split into two modules of N/2 bits:
  • Module 0 (bits 0 to N/2 - 1): dependencies are chosen completely at random from all N bits.
  • Module 1 (bits N/2 to N - 1): R controls the fraction of dependencies that stay within module 1. A fraction (1 - R) are drawn from module 0 (cross-module).
This lets you study how modularity interacts with ruggedness: a high R landscape has two semi-independent search problems, while R=0 creates tight cross-module coupling.
The maximum valid R is constrained by K and module size. If R is too high to satisfy K intra-module dependencies, the constructor raises a ValueError with the maximum permissible value.

​
Fitness function

Fitness is calculated in two stages:
  1. Raw fitness — for each bit i, look up fitness_contributions[i] using the binary states of i and its K dependencies as an index. Average across all N bits.
  2. Normalization and transformation — the raw value is linearly rescaled to [0, 1] using the true min and max across all 2^N states, then raised to the power of 4:
normalized = (raw - true_min) / (true_max - true_min)
transformed = normalized ** 4
The ^4 transform concentrates the distribution near zero, making high-fitness states rare and meaningful.

​
Memory-efficient bounds calculation

Finding the true min and max requires evaluating all 2^N states. OptimizedNKLandscape does this with a streaming pass at construction time — it never stores all fitness values at once, processing states in batches of 10,000: 
for batch_start in range(0, 2**self.N, batch_size):
    for i in range(batch_start, batch_end):
        state = np.array([(i >> bit) & 1 for bit in range(self.N)], dtype=np.uint8)
        raw_fitness = self._get_raw_fitness(state)
        min_fitness = min(min_fitness, raw_fitness)
        max_fitness = max(max_fitness, raw_fitness)
This means construction is O(2^N) in time but O(N) in memory.

​
Key methods

Compute the normalized, transformed fitness for a binary state array.
landscape = OptimizedNKLandscape(N=20, K=4, R=0.5, seed=42)
state = np.array([1, 0, 1, 1, 0, ...], dtype=np.uint8)  # length N
fitness = landscape.get_fitness(state)  # float in [0, 1]
Results are cached with LRU eviction. When the cache reaches max_cache_size, the least-recently-used entry is evicted.
Generate a random binary state vector of length N.
state = landscape.random_state()  # np.ndarray of dtype uint8
Uses the landscape’s internal RandomState, so behavior is reproducible when seed is set.
Return all single-bit-flip neighbors of a state — a list of N arrays.
neighbors = landscape.get_neighbors(state)
# neighbors[i] is `state` with bit i flipped
Useful for greedy hill-climbing: evaluate all neighbors and move to the best.
Return a dict summarizing the dependency structure:
analysis = landscape.analyze_modularity()
# {
#   'total_dependencies': int,
#   'module0_intra_dependencies': int,
#   'module1_intra_dependencies': int,
#   'inter_module_dependencies': int,
#   'target_R': float,
#   'actual_R_module1': float,
#   'module1_intra_fraction': float,
#   'avg_dependencies_per_bit': float
# }
actual_R_module1 may differ slightly from target_R due to integer rounding when computing int(round(K * R)).
Plot the N × N binary interaction matrix as a heatmap with module boundaries highlighted. Requires matplotlib.
fig = landscape.visualize_interactions(save_path="interactions.png")
Module 0 is outlined in green, module 1 in orange, and the boundary between modules is marked with red lines.

​
Using with AgentModel

To simulate agents searching the NK landscape with AgentModel, store the landscape as a parameter and access it inside your timestep function.
1

Create the landscape and model

from emergent.main import AgentModel
from LandscapeConstruction_NK_Modular import OptimizedNKLandscape
import numpy as np

landscape = OptimizedNKLandscape(N=20, K=4, R=0.5, seed=0)

model = AgentModel()
model.update_parameters({
    "num_nodes": 20,
    "graph_type": "complete",
    "convergence_data_key": "fitness",
    "convergence_std_dev": 0.01,
})
model["landscape"] = landscape
2

Define the initial data function

Each node gets a random binary state and its fitness on the landscape.
def initial_data(m):
    ls = m["landscape"]
    state = ls.random_state()
    return {
        "state": state,
        "fitness": ls.get_fitness(state),
    }

model.set_initial_data_function(initial_data)
3

Define the timestep function

Each agent either explores (random single-bit flip) or copies the fittest neighbor.
import random

def timestep(m):
    graph = m.get_graph()
    ls = m["landscape"]
    velocity = m["velocity"]

    for node in graph.nodes():
        node_data = graph.nodes[node]
        if random.random() <= velocity:
            # Social learning: copy fittest neighbor
            neighbors = list(graph.neighbors(node))
            if neighbors:
                best = max(neighbors, key=lambda n: graph.nodes[n]["fitness"])
                if graph.nodes[best]["fitness"] > node_data["fitness"]:
                    node_data["state"] = graph.nodes[best]["state"].copy()
                    node_data["fitness"] = graph.nodes[best]["fitness"]
        else:
            # Individual exploration: flip a random bit
            new_state = node_data["state"].copy()
            bit = random.randint(0, ls.N - 1)
            new_state[bit] = 1 - new_state[bit]
            new_fitness = ls.get_fitness(new_state)
            if new_fitness > node_data["fitness"]:
                node_data["state"] = new_state
                node_data["fitness"] = new_fitness

model.set_timestep_function(timestep)
4

Run the simulation

model["velocity"] = 0.5
model.initialize_graph()
steps = model.run_to_convergence()
print(f"Converged in {steps} steps")

# Inspect final fitness values
graph = model.get_graph()
fitnesses = [graph.nodes[n]["fitness"] for n in graph.nodes()]
print(f"Mean fitness: {sum(fitnesses) / len(fitnesses):.4f}")

​
Cache management

For large N or long simulations, the cache can grow large. Call clear_cache() periodically to release memory: 
# Inside a long simulation loop
if t % 100 == 0:
    landscape.clear_cache()
Check current cache usage with get_cache_stats(): 
stats = landscape.get_cache_stats()
print(stats)
# {'cache_size': 43201, 'max_cache_size': 100000, 'cache_enabled': True}
For problems where N > 16, consider reducing max_cache_size and using more processes carefully. The simulation in AgentDecisionProcedures.py uses max_cache_size=50000 and caps parallel processes at 8 for N > 16.

Build docs developers (and LLMs) love

Get started for freeTalk to us