Skip to content

Stochastic Doctor

Bitstream-level stochastic correlation analysis and diagnostics engine. Detects correlation anomalies, precision drift, and produces full JSON-serializable audit reports per network layer.

Rust Acceleration

This module includes optional Rust PyO3 acceleration via the stochastic_doctor_core crate. When the compiled extension is present, all hot-path operations dispatch to Rust automatically.

Set SC_NEUROCORE_NO_RUST=1 to force the pure Python fallback.

Performance (Python vs Rust PyO3)

Benchmarked on x86_64, Python 3.12, Rust 1.86 (release):

SCC — Single Pair

Stream Length Python Rust Speedup
100 12 µs 0.3 µs 35×
1,000 29 µs 0.9 µs 32×
10,000 42 µs 4.6 µs 9×
100,000 128 µs 40 µs 3.2×
1,000,000 1,297 µs 369 µs 3.5×

Batch SCC — N×N Pairwise Matrix (stream_len=2048)

Neurons Pairs Python Rust Speedup
4 6 0.10 ms 0.007 ms 15×
8 28 0.40 ms 0.025 ms 16×
16 120 1.6 ms 0.10 ms 16×
32 496 7.0 ms 0.46 ms 15×
64 2,016 26.8 ms 1.5 ms 18×

Precision Estimation

Stream Length Python Rust Speedup
100 3 µs 0.2 µs 14×
1,000 4 µs 0.3 µs 12×
10,000 7 µs 0.8 µs 9×
100,000 32 µs 7 µs 5×
1,000,000 304 µs 61 µs 5×

Criterion Benchmarks (Pure Rust)

Function Input Time
scc_bytes 100 50 ns
scc_bytes 1,000 397 ns
scc_bytes 10,000 3.9 µs
scc_bytes 100,000 48 µs
scc_packed 256 bits 15 ns
scc_packed 1,024 bits 36 ns
scc_packed 8,192 bits 219 ns
scc_packed 65,536 bits 1.7 µs
scc_batch 4 streams 4.6 µs
scc_batch 8 streams 22 µs
scc_batch 16 streams 97 µs
scc_batch 32 streams 398 µs
precision_packed 256 bits 5.3 ns
precision_packed 65,536 bits 519 ns
histogram_u64 64 words 110 ns
histogram_u64 4,096 words 4.6 µs
drift_detector 1,000 obs 4.8 µs

Quick Start

Python
from sc_neurocore.stochastic_doctor.diagnostics import (
    StochasticDoctor, DriftDetector, BitstreamAuditReport,
)

import numpy as np

doc = StochasticDoctor(correlation_threshold=0.3, critical_threshold=0.7)

# Audit a layer of bitstreams
bitstreams = np.random.randint(0, 2, size=(8, 2048), dtype=np.uint8)
report = doc.audit_layer("V1_Cortex", bitstreams)

print(report.status)         # OK / WARNING / CRITICAL
print(report.max_correlation)
print(report.to_json())

# Drift monitoring
dd = DriftDetector(alpha=0.1, threshold=0.3)
for scc_value in correlation_stream:
    if dd.observe(scc_value):
        print(f"Drift detected! EMA={dd.ema:.4f}")

Architecture

Text Only
diagnostics.py ──┬── PyO3 (stochastic_doctor_core.so)  ← primary
                 ├── Python/NumPy fallback              ← secondary
                 └── SC_NEUROCORE_NO_RUST=1             ← force Python

sc_neurocore.stochastic_doctor.diagnostics

Bitstream-level stochastic correlation analysis and diagnostics.

Extends sc_neurocore.doctor with bitstream-specific metrics:

  • SCC: Stochastic Cross-Correlation (Alaghi & Hayes, 2013)
  • Precision estimation: empirical P̂ with variance bound σ² = p(1-p)/N
  • Drift detection: EMA-based correlation drift monitoring
  • Activity histograms: per-word popcount distribution

Includes optional Rust FFI acceleration via stochastic_doctor_core.

AuditSeverity

Bases: Enum

Audit finding severity levels.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
63
64
65
66
67
68
class AuditSeverity(Enum):
    """Audit finding severity levels."""

    OK = "OK"
    WARNING = "WARNING"
    CRITICAL = "CRITICAL"

BitstreamAuditFinding dataclass

Single finding from a bitstream audit.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
71
72
73
74
75
76
77
78
79
@dataclass
class BitstreamAuditFinding:
    """Single finding from a bitstream audit."""

    category: str
    severity: AuditSeverity
    message: str
    metric: float = 0.0
    neuron_pair: Optional[tuple] = None

BitstreamAuditReport dataclass

Full bitstream-level audit report for a network layer.

JSON-serializable via to_json() for pipeline integration.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
 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
@dataclass
class BitstreamAuditReport:
    """Full bitstream-level audit report for a network layer.

    JSON-serializable via ``to_json()`` for pipeline integration.
    """

    layer: str
    stream_length: int
    num_neurons: int
    max_correlation: float = 0.0
    mean_precision: float = 0.0
    precision_variance: float = 0.0
    hot_neurons: List[tuple] = field(default_factory=list)
    findings: List[BitstreamAuditFinding] = field(default_factory=list)
    status: AuditSeverity = AuditSeverity.OK

    def to_dict(self) -> Dict:
        """Serialize to a plain dict (JSON-compatible)."""
        d = asdict(self)
        d["status"] = self.status.value
        d["findings"] = [{**asdict(f), "severity": f.severity.value} for f in self.findings]
        return d

    def to_json(self, indent: int = 2) -> str:
        """Serialize to JSON string."""
        return json.dumps(self.to_dict(), indent=indent)

to_dict()

Serialize to a plain dict (JSON-compatible).

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
 99
100
101
102
103
104
def to_dict(self) -> Dict:
    """Serialize to a plain dict (JSON-compatible)."""
    d = asdict(self)
    d["status"] = self.status.value
    d["findings"] = [{**asdict(f), "severity": f.severity.value} for f in self.findings]
    return d

to_json(indent=2)

Serialize to JSON string.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
106
107
108
def to_json(self, indent: int = 2) -> str:
    """Serialize to JSON string."""
    return json.dumps(self.to_dict(), indent=indent)

DriftDetector

Exponential moving average drift detector for SCC monitoring.

Tracks the running EMA of SCC values. Flags a drift event when |EMA| exceeds the threshold.

Parameters

alpha : float EMA smoothing factor (0.0–1.0; lower = smoother). threshold : float Absolute SCC value above which to flag drift.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
class DriftDetector:
    """Exponential moving average drift detector for SCC monitoring.

    Tracks the running EMA of SCC values. Flags a drift event when
    |EMA| exceeds the threshold.

    Parameters
    ----------
    alpha : float
        EMA smoothing factor (0.0–1.0; lower = smoother).
    threshold : float
        Absolute SCC value above which to flag drift.
    """

    def __init__(self, alpha: float = 0.1, threshold: float = 0.3):
        self.alpha = max(0.0, min(1.0, alpha))
        self.threshold = max(0.0, min(1.0, threshold))
        self.ema: float = 0.0
        self.active: bool = False
        self._history: List[float] = []

    def observe(self, scc_value: float) -> bool:
        """Feed a new SCC observation. Returns True if drift detected."""
        self.ema = self.alpha * scc_value + (1.0 - self.alpha) * self.ema
        self.active = abs(self.ema) > self.threshold
        self._history.append(self.ema)
        return self.active

    def reset(self) -> None:
        """Reset detector state."""
        self.ema = 0.0
        self.active = False
        self._history.clear()

    @property
    def history(self) -> List[float]:
        """EMA history for plotting/logging."""
        return self._history

history property

EMA history for plotting/logging.

observe(scc_value)

Feed a new SCC observation. Returns True if drift detected.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
174
175
176
177
178
179
def observe(self, scc_value: float) -> bool:
    """Feed a new SCC observation. Returns True if drift detected."""
    self.ema = self.alpha * scc_value + (1.0 - self.alpha) * self.ema
    self.active = abs(self.ema) > self.threshold
    self._history.append(self.ema)
    return self.active

reset()

Reset detector state.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
181
182
183
184
185
def reset(self) -> None:
    """Reset detector state."""
    self.ema = 0.0
    self.active = False
    self._history.clear()

StochasticDoctor

Bitstream-level stochastic diagnostics engine.

Parameters

correlation_threshold : float |SCC| above this triggers a WARNING (default 0.3). critical_threshold : float |SCC| above this triggers a CRITICAL (default 0.7).

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
class StochasticDoctor:
    """Bitstream-level stochastic diagnostics engine.

    Parameters
    ----------
    correlation_threshold : float
        |SCC| above this triggers a WARNING (default 0.3).
    critical_threshold : float
        |SCC| above this triggers a CRITICAL (default 0.7).
    """

    def __init__(
        self,
        correlation_threshold: float = 0.3,
        critical_threshold: float = 0.7,
    ):
        self.correlation_threshold = correlation_threshold
        self.critical_threshold = critical_threshold

    def compute_correlation(self, a: np.ndarray, b: np.ndarray) -> float:
        """Compute SCC between two bitstreams."""
        return compute_scc(a, b)

    def estimate_precision(self, bitstream: np.ndarray) -> tuple:
        """Estimate probability and variance bound for a bitstream.

        Uses Rust PyO3 acceleration when available.

        Returns
        -------
        (probability, variance_bound)
        """
        if _HAS_PYO3 and _sdc_rust is not None:
            bs = require_c_contiguous(bitstream, "bitstream", np.uint8)
            return _sdc_rust.py_precision_bytes(bs)
        bs = np.ascontiguousarray(bitstream, dtype=np.uint8)
        n = len(bs)
        if n == 0:
            return (0.0, 0.0)
        p = float(np.mean(bs))
        variance = p * (1.0 - p) / n
        return (p, variance)

    def compute_histogram(self, bitstream: np.ndarray, word_size: int = 64) -> np.ndarray:
        """Compute per-word popcount histogram.

        Splits the bitstream into chunks of ``word_size`` and counts
        the popcount of each chunk. Returns a histogram with bins 0..word_size.
        Uses Rust PyO3 acceleration when available.

        Parameters
        ----------
        bitstream : ndarray
            1D array of 0/1 values.
        word_size : int
            Number of bits per word (default 64).
        """
        if _HAS_PYO3 and _sdc_rust is not None:
            bs = require_c_contiguous(bitstream, "bitstream", np.uint8)
            return np.asarray(_sdc_rust.py_histogram(bs, word_size))
        bs = np.ascontiguousarray(bitstream, dtype=np.uint8)
        n = len(bs)
        hist = np.zeros(word_size + 1, dtype=np.int64)
        for start in range(0, n, word_size):
            chunk = bs[start : start + word_size]
            pc = int(np.sum(chunk))
            hist[pc] += 1
        return hist

    def audit_layer(self, layer_id: str, bitstreams: np.ndarray) -> BitstreamAuditReport:
        """Audit a full layer of bitstreams.

        Parameters
        ----------
        layer_id : str
            Human-readable layer identifier.
        bitstreams : ndarray
            Shape (num_neurons, stream_length), each element 0 or 1.

        Returns
        -------
        BitstreamAuditReport
        """
        num_neurons, stream_len = bitstreams.shape
        report = BitstreamAuditReport(
            layer=layer_id,
            stream_length=stream_len,
            num_neurons=num_neurons,
        )

        # Precision analysis
        precisions = []
        for i in range(num_neurons):
            p, var = self.estimate_precision(bitstreams[i])
            precisions.append(p)

        report.mean_precision = float(np.mean(precisions))
        report.precision_variance = float(np.var(precisions))

        # Pairwise SCC analysis
        max_corr = 0.0
        hot_pairs: List[tuple] = []

        for i in range(num_neurons):
            for j in range(i + 1, num_neurons):
                scc_val = self.compute_correlation(bitstreams[i], bitstreams[j])
                abs_scc = abs(scc_val)

                if abs_scc > abs(max_corr):
                    max_corr = scc_val

                if abs_scc > self.critical_threshold:
                    hot_pairs.append((i, j, scc_val))
                    report.findings.append(
                        BitstreamAuditFinding(
                            category="critical_correlation",
                            severity=AuditSeverity.CRITICAL,
                            message=f"Neurons ({i},{j}): SCC={scc_val:.4f} exceeds critical threshold",
                            metric=scc_val,
                            neuron_pair=(i, j),
                        )
                    )
                elif abs_scc > self.correlation_threshold:
                    hot_pairs.append((i, j, scc_val))
                    report.findings.append(
                        BitstreamAuditFinding(
                            category="high_correlation",
                            severity=AuditSeverity.WARNING,
                            message=f"Neurons ({i},{j}): SCC={scc_val:.4f} exceeds warning threshold",
                            metric=scc_val,
                            neuron_pair=(i, j),
                        )
                    )

        report.max_correlation = max_corr
        report.hot_neurons = hot_pairs

        # Overall status
        if any(f.severity == AuditSeverity.CRITICAL for f in report.findings):
            report.status = AuditSeverity.CRITICAL
        elif any(f.severity == AuditSeverity.WARNING for f in report.findings):
            report.status = AuditSeverity.WARNING
        else:
            report.status = AuditSeverity.OK

        return report

compute_correlation(a, b)

Compute SCC between two bitstreams.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
217
218
219
def compute_correlation(self, a: np.ndarray, b: np.ndarray) -> float:
    """Compute SCC between two bitstreams."""
    return compute_scc(a, b)

estimate_precision(bitstream)

Estimate probability and variance bound for a bitstream.

Uses Rust PyO3 acceleration when available.

Returns

(probability, variance_bound)

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
def estimate_precision(self, bitstream: np.ndarray) -> tuple:
    """Estimate probability and variance bound for a bitstream.

    Uses Rust PyO3 acceleration when available.

    Returns
    -------
    (probability, variance_bound)
    """
    if _HAS_PYO3 and _sdc_rust is not None:
        bs = require_c_contiguous(bitstream, "bitstream", np.uint8)
        return _sdc_rust.py_precision_bytes(bs)
    bs = np.ascontiguousarray(bitstream, dtype=np.uint8)
    n = len(bs)
    if n == 0:
        return (0.0, 0.0)
    p = float(np.mean(bs))
    variance = p * (1.0 - p) / n
    return (p, variance)

compute_histogram(bitstream, word_size=64)

Compute per-word popcount histogram.

Splits the bitstream into chunks of word_size and counts the popcount of each chunk. Returns a histogram with bins 0..word_size. Uses Rust PyO3 acceleration when available.

Parameters

bitstream : ndarray 1D array of 0/1 values. word_size : int Number of bits per word (default 64).

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
def compute_histogram(self, bitstream: np.ndarray, word_size: int = 64) -> np.ndarray:
    """Compute per-word popcount histogram.

    Splits the bitstream into chunks of ``word_size`` and counts
    the popcount of each chunk. Returns a histogram with bins 0..word_size.
    Uses Rust PyO3 acceleration when available.

    Parameters
    ----------
    bitstream : ndarray
        1D array of 0/1 values.
    word_size : int
        Number of bits per word (default 64).
    """
    if _HAS_PYO3 and _sdc_rust is not None:
        bs = require_c_contiguous(bitstream, "bitstream", np.uint8)
        return np.asarray(_sdc_rust.py_histogram(bs, word_size))
    bs = np.ascontiguousarray(bitstream, dtype=np.uint8)
    n = len(bs)
    hist = np.zeros(word_size + 1, dtype=np.int64)
    for start in range(0, n, word_size):
        chunk = bs[start : start + word_size]
        pc = int(np.sum(chunk))
        hist[pc] += 1
    return hist

audit_layer(layer_id, bitstreams)

Audit a full layer of bitstreams.

Parameters

layer_id : str Human-readable layer identifier. bitstreams : ndarray Shape (num_neurons, stream_length), each element 0 or 1.

Returns

BitstreamAuditReport

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
def audit_layer(self, layer_id: str, bitstreams: np.ndarray) -> BitstreamAuditReport:
    """Audit a full layer of bitstreams.

    Parameters
    ----------
    layer_id : str
        Human-readable layer identifier.
    bitstreams : ndarray
        Shape (num_neurons, stream_length), each element 0 or 1.

    Returns
    -------
    BitstreamAuditReport
    """
    num_neurons, stream_len = bitstreams.shape
    report = BitstreamAuditReport(
        layer=layer_id,
        stream_length=stream_len,
        num_neurons=num_neurons,
    )

    # Precision analysis
    precisions = []
    for i in range(num_neurons):
        p, var = self.estimate_precision(bitstreams[i])
        precisions.append(p)

    report.mean_precision = float(np.mean(precisions))
    report.precision_variance = float(np.var(precisions))

    # Pairwise SCC analysis
    max_corr = 0.0
    hot_pairs: List[tuple] = []

    for i in range(num_neurons):
        for j in range(i + 1, num_neurons):
            scc_val = self.compute_correlation(bitstreams[i], bitstreams[j])
            abs_scc = abs(scc_val)

            if abs_scc > abs(max_corr):
                max_corr = scc_val

            if abs_scc > self.critical_threshold:
                hot_pairs.append((i, j, scc_val))
                report.findings.append(
                    BitstreamAuditFinding(
                        category="critical_correlation",
                        severity=AuditSeverity.CRITICAL,
                        message=f"Neurons ({i},{j}): SCC={scc_val:.4f} exceeds critical threshold",
                        metric=scc_val,
                        neuron_pair=(i, j),
                    )
                )
            elif abs_scc > self.correlation_threshold:
                hot_pairs.append((i, j, scc_val))
                report.findings.append(
                    BitstreamAuditFinding(
                        category="high_correlation",
                        severity=AuditSeverity.WARNING,
                        message=f"Neurons ({i},{j}): SCC={scc_val:.4f} exceeds warning threshold",
                        metric=scc_val,
                        neuron_pair=(i, j),
                    )
                )

    report.max_correlation = max_corr
    report.hot_neurons = hot_pairs

    # Overall status
    if any(f.severity == AuditSeverity.CRITICAL for f in report.findings):
        report.status = AuditSeverity.CRITICAL
    elif any(f.severity == AuditSeverity.WARNING for f in report.findings):
        report.status = AuditSeverity.WARNING
    else:
        report.status = AuditSeverity.OK

    return report

compute_scc(a, b)

Compute SCC between two bitstreams.

Uses Rust PyO3 acceleration when available, falls back to pure Python. Set SC_NEUROCORE_NO_RUST=1 to force Python path.

Source code in src/sc_neurocore/stochastic_doctor/diagnostics.py
Python
133
134
135
136
137
138
139
140
141
142
143
144
145
def compute_scc(a: np.ndarray, b: np.ndarray) -> float:
    """Compute SCC between two bitstreams.

    Uses Rust PyO3 acceleration when available, falls back to pure Python.
    Set ``SC_NEUROCORE_NO_RUST=1`` to force Python path.
    """
    if _HAS_PYO3 and _sdc_rust is not None:
        a = require_c_contiguous(a, "a", np.uint8)
        b = require_c_contiguous(b, "b", np.uint8)
        return float(_sdc_rust.py_scc_bytes(a, b))
    a = np.ascontiguousarray(a, dtype=np.uint8)
    b = np.ascontiguousarray(b, dtype=np.uint8)
    return _scc_python(a, b)