physarum

Physarum polycephalum (slime mould) on the CPU, two ways. The agent model grows the organic transport networks the organism forms across a space. The flow solver finds shortest paths and builds networks by the Tero adaptive-conductivity model. Both vectorise with NEON on aarch64, with a portable scalar fallback; the agent update is also threaded and tile-sorted.

1200x1200 field, 680k agents, on a Raspberry Pi 4.

Network solver

Tero et al. (2010, Rules for biologically inspired adaptive network design): current flows between terminals on a masked grid; each step reinforces high-flux edges and lets the rest decay, so the conductivities converge onto the shortest path (a single source and sink) or an efficient network (several terminals). The inner step is a weighted-graph-Laplacian solve, exposed as the flow_cg kernel op (Jacobi-preconditioned conjugate gradients over a NEON 5-point mat-vec).

import torch
from kernels import get_kernel

physarum = get_kernel("phanerozoic/physarum", version=1, trust_remote_code=True)

open_mask = physarum.maze(127, 127, seed=5)      # 1 = open, 0 = wall
solver = physarum.PhysarumFlow(open_mask)
solver.solve([(2, 2), (124, 124)], [1, -1])      # a source and a sink
path = solver.path((2, 2), (124, 124))           # [(x, y), ...], the shortest route

solve(terminals, inject) takes any list of terminals and per-terminal net current; .network(thresh) returns the boolean mask of the converged network and .path(a, b) the route between two terminals through it. On recursive-division mazes the recovered path length equals the breadth-first-search optimum. A multi-terminal run connects every terminal with a network of spanning-tree length.

Single-core solve time (200 adaptation steps, flow_cg):

CG is sequential and the grids are modest, so the solve is vectorised but not threaded; the mat-vec is a 5-point weighted stencil in gather form.

vs classical algorithms

benchmark.py compares the solver to breadth-first search, Dijkstra, and the terminals' minimum spanning tree on a 4-connected grid, averaged over random configurations:

The flow model converges to the shortest path for a single source and sink, so it reproduces Dijkstra's length exactly (Dijkstra is faster in time). solve connects several terminals with a network of spanning-tree length; solve_robust reinforces the flux over every terminal pair to build a network that survives any single edge cut. benchmark.py reports the measured lengths against BFS, Dijkstra, and the MST across random configurations.

Agent model

Agents deposit a chemoattractant onto a scalar field and steer up its gradient; the field is diffused and decayed each step. The emergent network is the structure the organism forms spanning a space. Per step (Jones 2010, Characteristics of pattern formation and evolution in approximations of Physarum transport networks): each agent samples the field at three sensors (ahead, and +/- sense_ang, sense_dist px out), rotates by turn toward the strongest, moves speed px, and deposits deposit; the field is then convolved with a separable 3-tap [1 2 1]/4 kernel and multiplied by decay.

sim = physarum.Physarum(width=1024, height=1024, agents=200_000, seed=0)
for _ in range(400):
    sim.step()
img = sim.image()            # [H, W, 3] uint8, inferno colormap

sim.trail is the [H, W] float32 field. The agent update runs across cores via at::get_num_threads() (set with torch.set_num_threads); agents are counting-sorted into 32px tiles every eight steps so consecutive and same-core agents sense a hot window rather than scattering three gathers across the field. The step is deterministic for any thread count.

1024x1024 field, 200k agents, single core and four cores (NEON, tile-sorted):

API

Scope

• float32, CPU. NEON on aarch64; portable scalar elsewhere.
• Flow solver: 4-connected grid, one node grounded for the gauge; deterministic.
• Agent model: threaded update (torch.set_num_threads), serial deposit and diffusion, deterministic for any thread count.
• Behavioural and phenomenological models of Physarum, not biophysical ones.

Build from source

get_kernel loads a prebuilt aarch64 or x86_64 variant for torch 2.11-2.13. For a torch version outside that set, clone the repo and call load_local.load(), which JIT-compiles the source with torch.utils.cpp_extension and returns the same module.

License

Apache-2.0.