Skip to content

SPDX-License-Identifier: AGPL-3.0-or-later

Physics — Stochastic PDE Solvers + Wolfram Hypergraph

Physics simulation using stochastic methods: heat equation via reflected Brownian Feynman-Kac paths and Wolfram Physics Project hypergraph evolution.

FeynmanKacHeatSolver — 1D Heat Equation via Reflected Brownian Motion

Solves the 1D heat equation using the Feynman-Kac connection between diffusion PDEs and Brownian motion. Walkers follow Euler-Maruyama Brownian increments with variance 2 * diffusivity * dt and exact reflective Neumann boundaries on [0, length].

This is not a clipped lattice random walk. Boundary reflection uses triangle-wave folding with period 2 * length, so arbitrarily large stochastic increments remain inside the physical domain without changing the reflected transition kernel.

Parameter Meaning
length Positive domain length L
diffusivity Non-negative heat-equation coefficient alpha
num_walkers Positive number of Monte Carlo walkers
dt Positive Euler-Maruyama timestep
seed Integer seed for reproducible trajectories

Methods:

  • set_initial_delta(x_0) — initialize all walkers at a point in [0, L]
  • set_initial_distribution(f, n_grid) — sample a non-negative finite density
  • step(n_substeps) — advance reflected Brownian paths
  • evolve_to(T) — advance monotonically to a finite target time
  • get_density(n_bins) — probability density histogram integrating to one
  • expectation(observable) — Monte Carlo Feynman-Kac expectation

StochasticHeatSolver is retained as a backwards-compatible alias for FeynmanKacHeatSolver.

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; each edge must be a non-empty tuple of unique non-negative integers
max_node_id Non-negative integer at least as large as every node appearing in edges

Methods:

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

The implementation validates the hypergraph before rewrite and before dimension estimation. Malformed edges, repeated nodes inside one hyperedge, negative node identifiers, or stale max_node_id values are rejected rather than silently corrupting the emergent topology.

Usage

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

# Heat equation
solver = FeynmanKacHeatSolver(
    length=1.0,
    diffusivity=0.1,
    num_walkers=10000,
    dt=1e-3,
    seed=42,
)
solver.set_initial_delta(0.5)
solver.evolve_to(0.2)
profile = solver.get_density(n_bins=64)
u_cos = solver.expectation(lambda x: np.cos(np.pi * x))

# 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
 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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
@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 __post_init__(self) -> None:
        self.edges = self._validated_edges(self.edges)
        self.max_node_id = self._validated_max_node_id(self.max_node_id, self.edges)

    @staticmethod
    def _validated_edges(edges: List[Tuple[int, ...]]) -> List[Tuple[int, ...]]:
        if not isinstance(edges, list):
            raise ValueError("edges must be a list of integer tuples")
        validated: list[tuple[int, ...]] = []
        for edge in edges:
            if not isinstance(edge, tuple) or len(edge) == 0:
                raise ValueError("edges must contain non-empty integer tuples")
            if any(isinstance(node, bool) or not isinstance(node, int) for node in edge):
                raise ValueError("edge nodes must be integers")
            if any(node < 0 for node in edge):
                raise ValueError("edge nodes must be non-negative")
            if len(set(edge)) != len(edge):
                raise ValueError("hyperedges must not repeat nodes")
            validated.append(edge)
        return validated

    @staticmethod
    def _validated_max_node_id(max_node_id: int, edges: List[Tuple[int, ...]]) -> int:
        if isinstance(max_node_id, bool) or not isinstance(max_node_id, int):
            raise ValueError("max_node_id must be a non-negative integer")
        if max_node_id < 0:
            raise ValueError("max_node_id must be a non-negative integer")
        observed = max((node for edge in edges for node in edge), default=-1)
        if observed > max_node_id:
            raise ValueError("max_node_id must be at least the largest node in edges")
        return max_node_id

    @staticmethod
    def _validated_steps(steps: int) -> int:
        if isinstance(steps, bool) or not isinstance(steps, int):
            raise ValueError("steps must be a non-negative integer")
        if steps < 0:
            raise ValueError("steps must be a non-negative integer")
        return steps

    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)
        """
        steps = self._validated_steps(steps)
        self.edges = self._validated_edges(self.edges)
        self.max_node_id = self._validated_max_node_id(self.max_node_id, self.edges)
        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 = self._validated_edges(new_edges)  # type: ignore[arg-type]
            self.max_node_id = self._validated_max_node_id(self.max_node_id, self.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.
        """
        self.edges = self._validated_edges(self.edges)
        self.max_node_id = self._validated_max_node_id(self.max_node_id, self.edges)
        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
 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
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)
    """
    steps = self._validated_steps(steps)
    self.edges = self._validated_edges(self.edges)
    self.max_node_id = self._validated_max_node_id(self.max_node_id, self.edges)
    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 = self._validated_edges(new_edges)  # type: ignore[arg-type]
        self.max_node_id = self._validated_max_node_id(self.max_node_id, self.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
Python
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
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.
    """
    self.edges = self._validated_edges(self.edges)
    self.max_node_id = self._validated_max_node_id(self.max_node_id, self.edges)
    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))