Skip to content

Physics — Stochastic PDE Solvers + Wolfram Hypergraph

Physics simulation using stochastic methods: heat equation via random walks (Feynman-Kac) and Wolfram Physics Project hypergraph evolution.

StochasticHeatSolver — 1D Heat Equation via Random Walks

Solves the 1D heat equation using the Feynman-Kac connection between diffusion PDEs and Brownian motion. N random walkers perform discrete random walks on a 1D lattice; their density at time t approximates the temperature profile u(x,t).

Walker dynamics: at each step, move -1, 0, or +1 with probabilities [0.25, 0.5, 0.25]. Reflective boundary conditions (walkers clipped to [0, length-1]).

Parameter Meaning
length Lattice size (spatial resolution)
num_walkers Number of random walkers (more = smoother profile)
alpha Diffusion coefficient (reserved for future dt scaling)

Methods: step() (advance one timestep), get_temperature_profile() → normalized density histogram.

WolframHypergraph — Discrete Space-Time Evolution

Simulates the Wolfram Physics Project hypergraph rewrite system. The universe is a set of hyperedges (relations between nodes). Evolution applies a rewrite rule:

{{x, y}, {y, z}} → {{x, z}, {x, w}, {y, w}}

This is triangle completion with a new node w. The graph grows at each step, and the effective spatial dimension emerges from the topology.

Parameter Meaning
edges Initial hyperedge list (list of int tuples)
max_node_id Highest existing node ID

Methods:

  • evolve(steps) — Apply rewrite rule for N steps
  • dimension_estimate() — Estimate effective dimension via BFS neighborhood growth: measures how |B(r)| scales with r, fits d from V(r) ~ r^d

Usage

Python
from sc_neurocore.physics.heat import StochasticHeatSolver
from sc_neurocore.physics.wolfram_hypergraph import WolframHypergraph
import numpy as np

# Heat equation
solver = StochasticHeatSolver(length=100, num_walkers=10000, alpha=0.1)
solver.walkers[:] = 50  # point source at center
for _ in range(200):
    solver.step()
profile = solver.get_temperature_profile()  # diffused Gaussian-like

# Wolfram hypergraph
wh = WolframHypergraph(
    edges=[(0, 1), (1, 2), (2, 3), (3, 4)],
    max_node_id=4,
)
wh.evolve(steps=10)
print(f"Edges: {len(wh.edges)}, Nodes: {wh.max_node_id}")
print(f"Effective dimension: {wh.dimension_estimate():.2f}")

sc_neurocore.physics

sc_neurocore.physics -- Tier: research (experimental / research).

WolframHypergraph dataclass

Simulates the Wolfram Physics Project Hypergraph. Universe is a set of relations (Hyperedges).

Source code in src/sc_neurocore/physics/wolfram_hypergraph.py
Python
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
@dataclass
class WolframHypergraph:
    """
    Simulates the Wolfram Physics Project Hypergraph.
    Universe is a set of relations (Hyperedges).
    """

    edges: List[Tuple[int, ...]]
    max_node_id: int

    def evolve(self, steps: int = 1) -> None:
        """
        Applies a rewrite rule.
        Rule: {{x, y}, {y, z}} -> {{x, z}, {x, w}, {y, w}}
        (Triangle completion with new node w)
        """
        for _ in range(steps):
            new_edges = []
            matched_indices = set()

            # Naive pattern matching O(E^2)
            # Find (x, y) and (y, z)
            for i, e1 in enumerate(self.edges):
                if i in matched_indices:
                    continue
                if len(e1) != 2:
                    continue

                x, y = e1

                for j, e2 in enumerate(self.edges):
                    if i == j or j in matched_indices:
                        continue
                    if len(e2) != 2:
                        continue

                    if e2[0] == y:  # Found chain x->y->z
                        z = e2[1]

                        # Apply Rule
                        w = self.max_node_id + 1
                        self.max_node_id += 1

                        # New edges: {x,z}, {x,w}, {y,w}
                        new_edges.append((x, z))
                        new_edges.append((x, w))
                        new_edges.append((y, w))

                        matched_indices.add(i)
                        matched_indices.add(j)
                        break

            # Keep unmatched edges
            for k, e in enumerate(self.edges):
                if k not in matched_indices:
                    new_edges.append(e)  # type: ignore[arg-type]

            self.edges = new_edges  # type: ignore[assignment]

    def dimension_estimate(self) -> float:
        """Estimate effective dimension via BFS neighborhood growth.

        Builds an adjacency graph from hyperedges, picks a random node,
        measures how |B(r)| grows with r. Fits d from V(r) ~ r^d.
        Returns 0.0 if the graph is too small for meaningful estimation.
        """
        if len(self.edges) < 3:
            return 0.0

        adj: dict[int, set[int]] = {}
        for edge in self.edges:
            for node in edge:
                adj.setdefault(node, set())
            for i in range(len(edge)):
                for j in range(i + 1, len(edge)):
                    adj[edge[i]].add(edge[j])
                    adj[edge[j]].add(edge[i])

        nodes = list(adj.keys())
        if len(nodes) < 4:
            return 0.0

        start = nodes[len(nodes) // 2]
        visited = {start}
        frontier = {start}
        volumes = []

        for _ in range(min(10, len(nodes))):
            next_frontier: set[int] = set()
            for n in frontier:
                for nb in adj.get(n, set()):
                    if nb not in visited:
                        visited.add(nb)
                        next_frontier.add(nb)
            if not next_frontier:
                break
            frontier = next_frontier
            volumes.append(len(visited))

        if len(volumes) < 2:
            return 0.0

        import numpy as np

        r_vals = np.arange(1, len(volumes) + 1, dtype=np.float64)
        v_vals = np.array(volumes, dtype=np.float64)
        log_r = np.log(r_vals)
        log_v = np.log(np.clip(v_vals, 1, None))
        if log_r[-1] - log_r[0] < 1e-10:  # pragma: no cover
            return 0.0
        slope = (log_v[-1] - log_v[0]) / (log_r[-1] - log_r[0])
        return float(max(slope, 0.0))

evolve(steps=1)

Applies a rewrite rule. Rule: {{x, y}, {y, z}} -> {{x, z}, {x, w}, {y, w}} (Triangle completion with new node w)

Source code in src/sc_neurocore/physics/wolfram_hypergraph.py
Python
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
def evolve(self, steps: int = 1) -> None:
    """
    Applies a rewrite rule.
    Rule: {{x, y}, {y, z}} -> {{x, z}, {x, w}, {y, w}}
    (Triangle completion with new node w)
    """
    for _ in range(steps):
        new_edges = []
        matched_indices = set()

        # Naive pattern matching O(E^2)
        # Find (x, y) and (y, z)
        for i, e1 in enumerate(self.edges):
            if i in matched_indices:
                continue
            if len(e1) != 2:
                continue

            x, y = e1

            for j, e2 in enumerate(self.edges):
                if i == j or j in matched_indices:
                    continue
                if len(e2) != 2:
                    continue

                if e2[0] == y:  # Found chain x->y->z
                    z = e2[1]

                    # Apply Rule
                    w = self.max_node_id + 1
                    self.max_node_id += 1

                    # New edges: {x,z}, {x,w}, {y,w}
                    new_edges.append((x, z))
                    new_edges.append((x, w))
                    new_edges.append((y, w))

                    matched_indices.add(i)
                    matched_indices.add(j)
                    break

        # Keep unmatched edges
        for k, e in enumerate(self.edges):
            if k not in matched_indices:
                new_edges.append(e)  # type: ignore[arg-type]

        self.edges = new_edges  # type: ignore[assignment]

dimension_estimate()

Estimate effective dimension via BFS neighborhood growth.

Builds an adjacency graph from hyperedges, picks a random node, measures how |B(r)| grows with r. Fits d from V(r) ~ r^d. Returns 0.0 if the graph is too small for meaningful estimation.

Source code in src/sc_neurocore/physics/wolfram_hypergraph.py
Python
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
def dimension_estimate(self) -> float:
    """Estimate effective dimension via BFS neighborhood growth.

    Builds an adjacency graph from hyperedges, picks a random node,
    measures how |B(r)| grows with r. Fits d from V(r) ~ r^d.
    Returns 0.0 if the graph is too small for meaningful estimation.
    """
    if len(self.edges) < 3:
        return 0.0

    adj: dict[int, set[int]] = {}
    for edge in self.edges:
        for node in edge:
            adj.setdefault(node, set())
        for i in range(len(edge)):
            for j in range(i + 1, len(edge)):
                adj[edge[i]].add(edge[j])
                adj[edge[j]].add(edge[i])

    nodes = list(adj.keys())
    if len(nodes) < 4:
        return 0.0

    start = nodes[len(nodes) // 2]
    visited = {start}
    frontier = {start}
    volumes = []

    for _ in range(min(10, len(nodes))):
        next_frontier: set[int] = set()
        for n in frontier:
            for nb in adj.get(n, set()):
                if nb not in visited:
                    visited.add(nb)
                    next_frontier.add(nb)
        if not next_frontier:
            break
        frontier = next_frontier
        volumes.append(len(visited))

    if len(volumes) < 2:
        return 0.0

    import numpy as np

    r_vals = np.arange(1, len(volumes) + 1, dtype=np.float64)
    v_vals = np.array(volumes, dtype=np.float64)
    log_r = np.log(r_vals)
    log_v = np.log(np.clip(v_vals, 1, None))
    if log_r[-1] - log_r[0] < 1e-10:  # pragma: no cover
        return 0.0
    slope = (log_v[-1] - log_v[0]) / (log_r[-1] - log_r[0])
    return float(max(slope, 0.0))