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

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).

StochasticHeatSolver

Solves 1D Heat Equation using Stochastic Random Walks (Feynman-Kac).

Source code in src/sc_neurocore/physics/heat.py 
12
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
class StochasticHeatSolver:
    """
    Solves 1D Heat Equation using Stochastic Random Walks (Feynman-Kac).
    """

    def __init__(self, length: int, num_walkers: int, alpha: float):
        self.length = length
        self.walkers = np.random.randint(0, length, num_walkers)
        self.alpha = alpha

    def step(self) -> None:
        """
        Move walkers.
        """
        # Random step -1, 0, 1
        steps = np.random.choice([-1, 0, 1], size=len(self.walkers), p=[0.25, 0.5, 0.25])
        self.walkers += steps

        # Boundary conditions (Reflective)
        self.walkers = np.clip(self.walkers, 0, self.length - 1)

    def get_temperature_profile(self) -> np.ndarray[Any, Any]:
        """
        Convert walker density to temperature.
        """
        density, _ = np.histogram(self.walkers, bins=self.length, range=(0, self.length))
        return density / len(self.walkers)

step()

Move walkers.

Source code in src/sc_neurocore/physics/heat.py 
22
23
24
25
26
27
28
29
30
31
def step(self) -> None:
    """
    Move walkers.
    """
    # Random step -1, 0, 1
    steps = np.random.choice([-1, 0, 1], size=len(self.walkers), p=[0.25, 0.5, 0.25])
    self.walkers += steps

    # Boundary conditions (Reflective)
    self.walkers = np.clip(self.walkers, 0, self.length - 1)

get_temperature_profile()

Convert walker density to temperature.

Source code in src/sc_neurocore/physics/heat.py 
33
34
35
36
37
38
def get_temperature_profile(self) -> np.ndarray[Any, Any]:
    """
    Convert walker density to temperature.
    """
    density, _ = np.histogram(self.walkers, bins=self.length, range=(0, self.length))
    return density / len(self.walkers)

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 
 12
 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
@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)

            self.edges = new_edges

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

        self.edges = new_edges

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