Skip to content

SCPN Layers

Stochastic Coupled Phase-oscillator Network layer adapters.

Maps the 16-layer SCPN holonomic model (L1 Quantum through L16 Meta) into SC-NeuroCore's simulation engine. JAX-accelerated Kuramoto coupling, UPDE solvers, and phase coherence metrics.

Python
from sc_neurocore.scpn import KuramotoCoupling, PhaseCoherenceMetric

sc_neurocore.scpn

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

L1_QuantumLayer

Stochastic implementation of the Quantum Cellular Field.

Source code in src/sc_neurocore/scpn/layers/l1_quantum.py
Python
 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
173
174
175
176
177
178
179
180
181
182
183
184
185
class L1_QuantumLayer:
    """
    Stochastic implementation of the Quantum Cellular Field.
    """

    def __init__(self, params: Optional[L1_StochasticParameters] = None) -> None:
        self.params = params or L1_StochasticParameters()
        self._validate_params(self.params)
        self._rng = np.random.default_rng(self.params.rng_seed)

        # The Core Engine
        if self.params.backend == "simulated":
            self.quantum_core: Any = QuantumStochasticLayer(
                n_qubits=self.params.n_qubits,
                length=self.params.bitstream_length,
                rng_seed=self.params.rng_seed,
            )
        else:
            self.quantum_core = QuantumHardwareLayer(
                n_qubits=self.params.n_qubits,
                length=self.params.bitstream_length,
                backend_type=self.params.backend,
            )

        # State: Coherence represented as probabilities (0.0 to 1.0)
        # In SC, this will be encoded into bitstreams.
        # Initialize with max coherence (pure state)
        self.coherence_probs = np.ones(self.params.n_qubits) * 0.95

        # History for non-Markovian effects
        self.history: list[float] = []

    def step(
        self, dt: float, external_field: Optional[np.ndarray[Any, Any]] = None
    ) -> np.ndarray[Any, Any]:
        """
        Advance the layer by one time step.

        Args:
            dt: Time step in seconds.
            external_field: Optional coupling input from other layers (normalized 0-1).

        Returns:
            output_bitstreams: The stochastic state of the field.
        """
        self._validate_step_inputs(dt, external_field, self.params.n_qubits)
        # 1. Apply Decoherence (Classical Decay)
        # Adjusted by Non-Markovian factor
        effective_decay = self.params.decoherence_rate * dt / np.log10(self.params.F_non_Markov)
        self.coherence_probs *= 1.0 - effective_decay

        # 2. Apply External Coupling (e.g. from L2 Neurochemical)
        if external_field is not None:
            field = np.asarray(external_field, dtype=np.float64).reshape(self.params.n_qubits)
            # Mix the field: coherence is modulated by external input
            self.coherence_probs = (
                1 - self.params.coupling_strength
            ) * self.coherence_probs + self.params.coupling_strength * field

        # 3. Quantum Rotation via Stochastic Core
        # The core takes the probabilities, rotates them (simulating evolution),
        # and returns collapsed bitstreams.
        # We assume the 'probability' maps to the quantum phase/amplitude.

        # Generate input bitstreams from current probabilities
        rands = self._rng.random((self.params.n_qubits, self.params.bitstream_length))
        input_bits = (rands < self.coherence_probs[:, None]).astype(np.uint8)

        # Pass through Quantum Hybrid Layer
        output_bits = self.quantum_core.forward(input_bits)

        # 4. Update State from Measurement (Collapse/Update)
        # The output bits represent the measured state.
        # We update our internal probabilities based on the measurement (Bayesian update or similar)
        # For this simulation, we'll take the mean as the new base, but add some "Quantum Zeno" recovery
        measured_probs = np.mean(output_bits, axis=1)

        # "Zeno" effect: frequent measurement can freeze evolution or reset it.
        # Here we just blend it back.
        self.coherence_probs = 0.9 * self.coherence_probs + 0.1 * measured_probs

        res: np.ndarray[Any, Any] = output_bits
        return res

    def get_global_metric(self) -> float:
        """Return the global coherence metric (Phi-like)."""
        return float(np.mean(self.coherence_probs))

    @staticmethod
    def _validate_params(params: L1_StochasticParameters) -> None:
        if (
            not isinstance(params.n_qubits, int)
            or isinstance(params.n_qubits, bool)
            or params.n_qubits <= 0
        ):
            raise ValueError("n_qubits must be a positive integer")
        if (
            not isinstance(params.bitstream_length, int)
            or isinstance(params.bitstream_length, bool)
            or params.bitstream_length <= 0
        ):
            raise ValueError("bitstream_length must be a positive integer")
        if not math.isfinite(float(params.F_non_Markov)) or params.F_non_Markov <= 1.0:
            raise ValueError("F_non_Markov must be finite and greater than 1")
        if not math.isfinite(float(params.temperature)) or params.temperature <= 0.0:
            raise ValueError("temperature must be finite and positive")
        if (
            not math.isfinite(float(params.coupling_strength))
            or params.coupling_strength < 0.0
            or params.coupling_strength > 1.0
        ):
            raise ValueError("coupling_strength must be finite and within [0, 1]")
        if not math.isfinite(float(params.decoherence_rate)) or params.decoherence_rate < 0.0:
            raise ValueError("decoherence_rate must be finite and non-negative")
        if not isinstance(params.backend, str) or not params.backend:
            raise ValueError("backend must be a non-empty string")
        if params.rng_seed is not None:
            if isinstance(params.rng_seed, bool) or not isinstance(params.rng_seed, int):
                raise ValueError("rng_seed must be a non-negative integer or None")
            if params.rng_seed < 0:
                raise ValueError("rng_seed must be a non-negative integer or None")

    @staticmethod
    def _validate_step_inputs(
        dt: float, external_field: Optional[np.ndarray[Any, Any]], n_qubits: int
    ) -> None:
        if not math.isfinite(float(dt)) or dt <= 0.0:
            raise ValueError("dt must be finite and positive")
        if external_field is not None:
            field = np.asarray(external_field, dtype=np.float64)
            if field.size != n_qubits or not np.all(np.isfinite(field)):
                raise ValueError("external_field must contain one finite value per qubit")
            if np.any(field < 0.0) or np.any(field > 1.0):
                raise ValueError("external_field must be within [0, 1]")

step(dt, external_field=None)

Advance the layer by one time step.

Parameters:

Name Type Description Default
dt float

Time step in seconds.

 required
external_field Optional[ndarray[Any, Any]]

Optional coupling input from other layers (normalized 0-1).

 None

Returns:

Name Type Description
output_bitstreams ndarray[Any, Any]

The stochastic state of the field.
Source code in src/sc_neurocore/scpn/layers/l1_quantum.py
Python
 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
def step(
    self, dt: float, external_field: Optional[np.ndarray[Any, Any]] = None
) -> np.ndarray[Any, Any]:
    """
    Advance the layer by one time step.

    Args:
        dt: Time step in seconds.
        external_field: Optional coupling input from other layers (normalized 0-1).

    Returns:
        output_bitstreams: The stochastic state of the field.
    """
    self._validate_step_inputs(dt, external_field, self.params.n_qubits)
    # 1. Apply Decoherence (Classical Decay)
    # Adjusted by Non-Markovian factor
    effective_decay = self.params.decoherence_rate * dt / np.log10(self.params.F_non_Markov)
    self.coherence_probs *= 1.0 - effective_decay

    # 2. Apply External Coupling (e.g. from L2 Neurochemical)
    if external_field is not None:
        field = np.asarray(external_field, dtype=np.float64).reshape(self.params.n_qubits)
        # Mix the field: coherence is modulated by external input
        self.coherence_probs = (
            1 - self.params.coupling_strength
        ) * self.coherence_probs + self.params.coupling_strength * field

    # 3. Quantum Rotation via Stochastic Core
    # The core takes the probabilities, rotates them (simulating evolution),
    # and returns collapsed bitstreams.
    # We assume the 'probability' maps to the quantum phase/amplitude.

    # Generate input bitstreams from current probabilities
    rands = self._rng.random((self.params.n_qubits, self.params.bitstream_length))
    input_bits = (rands < self.coherence_probs[:, None]).astype(np.uint8)

    # Pass through Quantum Hybrid Layer
    output_bits = self.quantum_core.forward(input_bits)

    # 4. Update State from Measurement (Collapse/Update)
    # The output bits represent the measured state.
    # We update our internal probabilities based on the measurement (Bayesian update or similar)
    # For this simulation, we'll take the mean as the new base, but add some "Quantum Zeno" recovery
    measured_probs = np.mean(output_bits, axis=1)

    # "Zeno" effect: frequent measurement can freeze evolution or reset it.
    # Here we just blend it back.
    self.coherence_probs = 0.9 * self.coherence_probs + 0.1 * measured_probs

    res: np.ndarray[Any, Any] = output_bits
    return res

get_global_metric()

Return the global coherence metric (Phi-like).

Source code in src/sc_neurocore/scpn/layers/l1_quantum.py
Python
136
137
138
def get_global_metric(self) -> float:
    """Return the global coherence metric (Phi-like)."""
    return float(np.mean(self.coherence_probs))

L2_NeurochemicalLayer

Stochastic implementation of the Neurochemical Signaling Layer.

Models receptor-ligand binding, neurotransmitter dynamics, and second messenger cascades using bitstream representations.

Source code in src/sc_neurocore/scpn/layers/l2_neurochemical.py
Python
 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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
class L2_NeurochemicalLayer:
    """
    Stochastic implementation of the Neurochemical Signaling Layer.

    Models receptor-ligand binding, neurotransmitter dynamics, and
    second messenger cascades using bitstream representations.
    """

    # Neurotransmitter indices
    DA = 0  # Dopamine
    SEROTONIN = 1  # 5-HT
    NE = 2  # Norepinephrine
    ACH = 3  # Acetylcholine

    def __init__(self, params: Optional[L2_StochasticParameters] = None):
        self.params = params or L2_StochasticParameters()
        self._validate_params(self.params)
        self._rng = np.random.default_rng(self.params.rng_seed)

        # Receptor states: 0 = unbound, 1 = bound
        self.receptor_states = np.zeros(
            (self.params.n_neurotransmitter_types, self.params.n_receptors), dtype=np.float32
        )

        # Neurotransmitter concentrations (as probabilities for bitstream encoding)
        self.nt_concentrations = np.ones(self.params.n_neurotransmitter_types) * 0.5

        # Second messenger cascade state (cAMP, Ca2+, etc.)
        self.second_messenger_levels = np.zeros(self.params.n_neurotransmitter_types)

        # History for temporal dynamics
        self.history: list[Any] = []

    def step(
        self,
        dt: float,
        nt_release: Optional[np.ndarray[Any, Any]] = None,
        l1_input: Optional[np.ndarray[Any, Any]] = None,
    ) -> Dict[str, Any]:
        """
        Advance the layer by one time step.

        Args:
            dt: Time step in seconds.
            nt_release: Neurotransmitter release rates [4] (0-1 normalized).
            l1_input: Quantum layer input (coherence modulation).

        Returns:
            Dict with receptor_activity, second_messengers, output_bitstreams
        """
        self._validate_step_inputs(dt, nt_release, l1_input, self.params.n_neurotransmitter_types)
        # 1. Update neurotransmitter concentrations from release
        if nt_release is not None:
            self.nt_concentrations = np.clip(
                self.nt_concentrations
                + np.asarray(nt_release, dtype=np.float64) * dt
                - self.params.reuptake_rate * dt,
                0.0,
                1.0,
            )

        # 2. Receptor binding dynamics (stochastic)
        for nt_idx in range(self.params.n_neurotransmitter_types):
            nt_conc = self.nt_concentrations[nt_idx]

            # Binding: P(bind) = affinity * concentration * (1 - current_state)
            binding_prob = self.params.binding_affinity * nt_conc
            bind_mask = self._rng.random(self.params.n_receptors) < binding_prob * dt

            # Unbinding: P(unbind) = unbinding_rate * current_state
            unbind_mask = (
                self._rng.random(self.params.n_receptors) < self.params.unbinding_rate * dt
            )

            # Update states
            self.receptor_states[nt_idx] = np.where(
                bind_mask & (self.receptor_states[nt_idx] < 0.5), 1.0, self.receptor_states[nt_idx]
            )
            self.receptor_states[nt_idx] = np.where(
                unbind_mask & (self.receptor_states[nt_idx] > 0.5),
                0.0,
                self.receptor_states[nt_idx],
            )

        # 3. Second messenger cascade
        receptor_activity = np.mean(self.receptor_states, axis=1)
        self.second_messenger_levels = 0.9 * self.second_messenger_levels + 0.1 * receptor_activity

        # 4. Quantum coupling (L1 modulates receptor sensitivity)
        if l1_input is not None:
            quantum_mod = self._finite_mean(l1_input, "l1_input") * self.params.quantum_coupling
            self.receptor_states *= 1.0 + quantum_mod
            self.receptor_states[...] = np.clip(self.receptor_states, 0.0, 1.0).astype(
                np.float32, copy=False
            )

        # 5. Generate output bitstreams
        receptor_activity = np.mean(self.receptor_states, axis=1)
        output_probs = np.clip(receptor_activity, 0.0, 1.0)
        rands = self._rng.random(
            (self.params.n_neurotransmitter_types, self.params.bitstream_length)
        )
        output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
        genomic_drive = self.params.genomic_coupling * self.second_messenger_levels

        # Store history
        self.history.append(
            {
                "nt_concentrations": self.nt_concentrations.copy(),
                "receptor_activity": receptor_activity.copy(),
                "second_messengers": self.second_messenger_levels.copy(),
            }
        )
        if len(self.history) > 100:
            self.history.pop(0)

        return {
            "receptor_activity": receptor_activity,
            "second_messengers": self.second_messenger_levels.copy(),
            "genomic_drive": genomic_drive.copy(),
            "output_bitstreams": output_bitstreams,
            "nt_concentrations": self.nt_concentrations.copy(),
        }

    def release_neurotransmitter(self, nt_type: int, amount: float) -> None:
        """Trigger neurotransmitter release."""
        if not isinstance(nt_type, int) or isinstance(nt_type, bool):
            raise ValueError("nt_type must be a valid neurotransmitter index")
        if nt_type < 0 or nt_type >= self.params.n_neurotransmitter_types:
            raise ValueError("nt_type must be a valid neurotransmitter index")
        if not math.isfinite(float(amount)) or amount < 0.0:
            raise ValueError("amount must be finite and non-negative")
        self.nt_concentrations[nt_type] = np.clip(
            self.nt_concentrations[nt_type] + amount, 0.0, 1.0
        )

    def get_global_metric(self) -> float:
        """Return the global neurochemical activity metric."""
        return float(np.mean(self.receptor_states))

    def get_neuromodulation_state(self) -> Dict[str, float]:
        """Return named neurotransmitter levels for external use."""
        return {
            "dopamine": float(self.nt_concentrations[self.DA]),
            "serotonin": float(self.nt_concentrations[self.SEROTONIN]),
            "norepinephrine": float(self.nt_concentrations[self.NE]),
            "acetylcholine": float(self.nt_concentrations[self.ACH]),
        }

    @staticmethod
    def _validate_params(params: L2_StochasticParameters) -> None:
        if (
            not isinstance(params.n_receptors, int)
            or isinstance(params.n_receptors, bool)
            or params.n_receptors <= 0
        ):
            raise ValueError("n_receptors must be a positive integer")
        if (
            not isinstance(params.n_neurotransmitter_types, int)
            or isinstance(params.n_neurotransmitter_types, bool)
            or params.n_neurotransmitter_types <= 0
        ):
            raise ValueError("n_neurotransmitter_types must be a positive integer")
        if (
            not isinstance(params.bitstream_length, int)
            or isinstance(params.bitstream_length, bool)
            or params.bitstream_length <= 0
        ):
            raise ValueError("bitstream_length must be a positive integer")
        if (
            not math.isfinite(float(params.binding_affinity))
            or params.binding_affinity < 0.0
            or params.binding_affinity > 1.0
        ):
            raise ValueError("binding_affinity must be finite and within [0, 1]")
        for field_name in (
            "unbinding_rate",
            "diffusion_rate",
            "reuptake_rate",
            "quantum_coupling",
            "genomic_coupling",
        ):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value < 0.0:
                raise ValueError(f"{field_name} must be finite and non-negative")
        if params.rng_seed is not None:
            if isinstance(params.rng_seed, bool) or not isinstance(params.rng_seed, int):
                raise ValueError("rng_seed must be a non-negative integer or None")
            if params.rng_seed < 0:
                raise ValueError("rng_seed must be a non-negative integer or None")

    @classmethod
    def _validate_step_inputs(
        cls,
        dt: float,
        nt_release: Optional[np.ndarray[Any, Any]],
        l1_input: Optional[np.ndarray[Any, Any]],
        n_neurotransmitter_types: int,
    ) -> None:
        if not math.isfinite(float(dt)) or dt <= 0.0:
            raise ValueError("dt must be finite and positive")
        if nt_release is not None:
            release = np.asarray(nt_release, dtype=np.float64)
            if release.size != n_neurotransmitter_types or not np.all(np.isfinite(release)):
                raise ValueError("nt_release must contain one finite value per neurotransmitter")
            if np.any(release < 0.0) or np.any(release > 1.0):
                raise ValueError("nt_release must be within [0, 1]")
        if l1_input is not None:
            cls._finite_mean(l1_input, "l1_input")

    @staticmethod
    def _finite_mean(values: Any, name: str) -> float:
        arr = np.asarray(values, dtype=np.float64)
        if arr.size == 0 or not np.all(np.isfinite(arr)):
            raise ValueError(f"{name} must contain finite values")
        return float(np.mean(arr))

step(dt, nt_release=None, l1_input=None)

Advance the layer by one time step.

Parameters:

Name Type Description Default
dt float

Time step in seconds.

 required
nt_release Optional[ndarray[Any, Any]]

Neurotransmitter release rates [4] (0-1 normalized).

 None
l1_input Optional[ndarray[Any, Any]]

Quantum layer input (coherence modulation).

 None

Returns:

Type Description
Dict[str, Any]

Dict with receptor_activity, second_messengers, output_bitstreams
Source code in src/sc_neurocore/scpn/layers/l2_neurochemical.py
Python
 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
173
174
175
176
177
178
179
def step(
    self,
    dt: float,
    nt_release: Optional[np.ndarray[Any, Any]] = None,
    l1_input: Optional[np.ndarray[Any, Any]] = None,
) -> Dict[str, Any]:
    """
    Advance the layer by one time step.

    Args:
        dt: Time step in seconds.
        nt_release: Neurotransmitter release rates [4] (0-1 normalized).
        l1_input: Quantum layer input (coherence modulation).

    Returns:
        Dict with receptor_activity, second_messengers, output_bitstreams
    """
    self._validate_step_inputs(dt, nt_release, l1_input, self.params.n_neurotransmitter_types)
    # 1. Update neurotransmitter concentrations from release
    if nt_release is not None:
        self.nt_concentrations = np.clip(
            self.nt_concentrations
            + np.asarray(nt_release, dtype=np.float64) * dt
            - self.params.reuptake_rate * dt,
            0.0,
            1.0,
        )

    # 2. Receptor binding dynamics (stochastic)
    for nt_idx in range(self.params.n_neurotransmitter_types):
        nt_conc = self.nt_concentrations[nt_idx]

        # Binding: P(bind) = affinity * concentration * (1 - current_state)
        binding_prob = self.params.binding_affinity * nt_conc
        bind_mask = self._rng.random(self.params.n_receptors) < binding_prob * dt

        # Unbinding: P(unbind) = unbinding_rate * current_state
        unbind_mask = (
            self._rng.random(self.params.n_receptors) < self.params.unbinding_rate * dt
        )

        # Update states
        self.receptor_states[nt_idx] = np.where(
            bind_mask & (self.receptor_states[nt_idx] < 0.5), 1.0, self.receptor_states[nt_idx]
        )
        self.receptor_states[nt_idx] = np.where(
            unbind_mask & (self.receptor_states[nt_idx] > 0.5),
            0.0,
            self.receptor_states[nt_idx],
        )

    # 3. Second messenger cascade
    receptor_activity = np.mean(self.receptor_states, axis=1)
    self.second_messenger_levels = 0.9 * self.second_messenger_levels + 0.1 * receptor_activity

    # 4. Quantum coupling (L1 modulates receptor sensitivity)
    if l1_input is not None:
        quantum_mod = self._finite_mean(l1_input, "l1_input") * self.params.quantum_coupling
        self.receptor_states *= 1.0 + quantum_mod
        self.receptor_states[...] = np.clip(self.receptor_states, 0.0, 1.0).astype(
            np.float32, copy=False
        )

    # 5. Generate output bitstreams
    receptor_activity = np.mean(self.receptor_states, axis=1)
    output_probs = np.clip(receptor_activity, 0.0, 1.0)
    rands = self._rng.random(
        (self.params.n_neurotransmitter_types, self.params.bitstream_length)
    )
    output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
    genomic_drive = self.params.genomic_coupling * self.second_messenger_levels

    # Store history
    self.history.append(
        {
            "nt_concentrations": self.nt_concentrations.copy(),
            "receptor_activity": receptor_activity.copy(),
            "second_messengers": self.second_messenger_levels.copy(),
        }
    )
    if len(self.history) > 100:
        self.history.pop(0)

    return {
        "receptor_activity": receptor_activity,
        "second_messengers": self.second_messenger_levels.copy(),
        "genomic_drive": genomic_drive.copy(),
        "output_bitstreams": output_bitstreams,
        "nt_concentrations": self.nt_concentrations.copy(),
    }

release_neurotransmitter(nt_type, amount)

Trigger neurotransmitter release.

Source code in src/sc_neurocore/scpn/layers/l2_neurochemical.py
Python
181
182
183
184
185
186
187
188
189
190
191
def release_neurotransmitter(self, nt_type: int, amount: float) -> None:
    """Trigger neurotransmitter release."""
    if not isinstance(nt_type, int) or isinstance(nt_type, bool):
        raise ValueError("nt_type must be a valid neurotransmitter index")
    if nt_type < 0 or nt_type >= self.params.n_neurotransmitter_types:
        raise ValueError("nt_type must be a valid neurotransmitter index")
    if not math.isfinite(float(amount)) or amount < 0.0:
        raise ValueError("amount must be finite and non-negative")
    self.nt_concentrations[nt_type] = np.clip(
        self.nt_concentrations[nt_type] + amount, 0.0, 1.0
    )

get_global_metric()

Return the global neurochemical activity metric.

Source code in src/sc_neurocore/scpn/layers/l2_neurochemical.py
Python
193
194
195
def get_global_metric(self) -> float:
    """Return the global neurochemical activity metric."""
    return float(np.mean(self.receptor_states))

get_neuromodulation_state()

Return named neurotransmitter levels for external use.

Source code in src/sc_neurocore/scpn/layers/l2_neurochemical.py
Python
197
198
199
200
201
202
203
204
def get_neuromodulation_state(self) -> Dict[str, float]:
    """Return named neurotransmitter levels for external use."""
    return {
        "dopamine": float(self.nt_concentrations[self.DA]),
        "serotonin": float(self.nt_concentrations[self.SEROTONIN]),
        "norepinephrine": float(self.nt_concentrations[self.NE]),
        "acetylcholine": float(self.nt_concentrations[self.ACH]),
    }

L3_GenomicLayer

Stochastic implementation of the Genomic-Epigenomic Layer.

Models gene expression, epigenetic modifications, and bioelectric pattern formation using bitstream representations.

Source code in src/sc_neurocore/scpn/layers/l3_genomic.py
Python
 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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
class L3_GenomicLayer:
    """
    Stochastic implementation of the Genomic-Epigenomic Layer.

    Models gene expression, epigenetic modifications, and bioelectric
    pattern formation using bitstream representations.
    """

    def __init__(self, params: Optional[L3_StochasticParameters] = None):
        self.params = params or L3_StochasticParameters()
        self._validate_params(self.params)
        self._rng = np.random.default_rng(self.params.rng_seed)

        # Gene expression levels (mRNA proxy, 0-1 normalized)
        self.expression_levels = self._rng.random(self.params.n_genes) * 0.3

        # Protein concentrations
        self.protein_levels = self._rng.random(self.params.n_genes) * 0.2

        # Chromatin state: 0 = closed (silenced), 1 = open (active)
        self.chromatin_state = self._rng.random(self.params.n_genes) > 0.5
        self.chromatin_openness = self.chromatin_state.astype(np.float64)

        # Methylation pattern (0-1, higher = more methylated = silenced)
        self.methylation = self._rng.random(self.params.n_genes) * 0.3

        # Bioelectric membrane potential grid
        self.membrane_potential = np.ones(self.params.n_genes) * self.params.membrane_potential_rest

        # CISS spin state (for quantum-genomic coupling)
        self.spin_polarization = np.zeros(self.params.n_genes)

        # Sparse regulatory network adjacency.
        self.regulatory_matrix = self._init_regulatory_network()

    def _init_regulatory_network(self) -> np.ndarray[Any, Any]:
        """Initialize gene regulatory network with sparse connections."""
        # Sparse random regulatory matrix
        matrix = self._rng.random((self.params.n_genes, self.params.n_regulatory_elements))
        matrix = np.where(matrix > 0.9, matrix, 0)  # Sparse
        # Add some inhibitory connections
        matrix[:, : self.params.n_regulatory_elements // 3] *= -1
        return matrix

    def step(
        self,
        dt: float,
        l2_input: Optional[Dict[str, Any]] = None,
        bioelectric_signal: Optional[np.ndarray[Any, Any]] = None,
    ) -> Dict[str, Any]:
        """
        Advance the layer by one time step.

        Args:
            dt: Time step in seconds.
            l2_input: Neurochemical layer output (second messengers).
            bioelectric_signal: External bioelectric modulation.

        Returns:
            Dict with expression, protein_levels, chromatin_state, output_bitstreams
        """
        self._validate_step_inputs(dt, l2_input, bioelectric_signal, self.params.n_genes)
        # 1. Update chromatin state (epigenetic dynamics)
        # Methylation silences genes
        demeth_prob = self.params.demethylation_rate * dt
        meth_prob = self.params.methylation_rate * dt

        demeth_mask = self._rng.random(self.params.n_genes) < demeth_prob
        meth_mask = self._rng.random(self.params.n_genes) < meth_prob

        self.methylation = np.where(demeth_mask, self.methylation * 0.9, self.methylation)
        self.methylation = np.where(meth_mask, self.methylation + 0.1, self.methylation)
        self.methylation = np.clip(self.methylation, 0.0, 1.0)

        # Chromatin openness inversely related to methylation
        self.chromatin_openness = (
            1.0 - self.methylation + self._rng.normal(0, 0.05, self.params.n_genes)
        )
        self.chromatin_openness = np.clip(self.chromatin_openness, 0.0, 1.0)

        # 2. Gene expression (stochastic transcription)
        # Only open chromatin can be transcribed
        transcription_prob = self.params.transcription_rate * self.chromatin_openness * dt
        transcription = self._rng.random(self.params.n_genes) < transcription_prob

        self.expression_levels = np.where(
            transcription,
            self.expression_levels + 0.1,
            self.expression_levels - self.params.degradation_rate * dt,
        )
        self.expression_levels = np.clip(self.expression_levels, 0.0, 1.0)

        # 3. Translation to protein
        translation_prob = self.params.translation_rate * self.expression_levels * dt
        translation = self._rng.random(self.params.n_genes) < translation_prob

        self.protein_levels = np.where(
            translation,
            self.protein_levels + 0.05,
            self.protein_levels - self.params.degradation_rate * dt * 0.5,
        )
        self.protein_levels = np.clip(self.protein_levels, 0.0, 1.0)

        # 4. CISS effect (quantum spin filtering)
        # Spin polarization depends on DNA chirality and electron flow
        electron_flow = np.mean(self.expression_levels)  # Proxy for metabolic activity
        ciss_baseline = self.params.ciss_efficiency * self.params.dna_chirality * electron_flow
        spin_polarization = np.clip(
            np.full(self.params.n_genes, ciss_baseline, dtype=np.float64)
            + self._rng.normal(0, 0.1, self.params.n_genes),
            -1.0,
            1.0,
        )
        self.spin_polarization[...] = spin_polarization

        # 5. Neurochemical coupling (L2 input modulates expression)
        if l2_input is not None and "second_messengers" in l2_input:
            # cAMP from second messengers activates transcription factors
            camp_level = self._finite_mean(l2_input["second_messengers"], "second_messengers")
            activation_boost = camp_level * self.params.neurochemical_coupling
            self.expression_levels += activation_boost * dt
            self.expression_levels = np.clip(self.expression_levels, 0.0, 1.0)

        # 6. Bioelectric pattern formation
        if bioelectric_signal is not None:
            signal = self._bioelectric_signal(bioelectric_signal, self.params.n_genes)
            self.membrane_potential = 0.9 * self.membrane_potential + 0.1 * signal
        # Internal bioelectric dynamics (gap junction diffusion)
        diffusion = np.roll(self.membrane_potential, 1) - self.membrane_potential
        self.membrane_potential += diffusion * self.params.bioelectric_coupling * dt

        # 7. Generate output bitstreams
        output_probs = np.clip(self.protein_levels, 0.0, 1.0)
        rands = self._rng.random((self.params.n_genes, self.params.bitstream_length))
        output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
        cellular_drive = self.params.cellular_coupling * self.protein_levels

        return {
            "expression_levels": self.expression_levels.copy(),
            "protein_levels": self.protein_levels.copy(),
            "chromatin_openness": self.chromatin_openness.copy(),
            "methylation": self.methylation.copy(),
            "spin_polarization": self.spin_polarization.copy(),
            "membrane_potential": self.membrane_potential.copy(),
            "cellular_drive": cellular_drive.copy(),
            "output_bitstreams": output_bitstreams,
        }

    def get_global_metric(self) -> float:
        """Return the global genomic activity metric."""
        return float(np.mean(self.expression_levels))

    def get_ciss_coherence(self) -> float:
        """Return CISS spin coherence metric."""
        return float(np.abs(np.mean(self.spin_polarization)))

    @staticmethod
    def _validate_params(params: L3_StochasticParameters) -> None:
        if (
            not isinstance(params.n_genes, int)
            or isinstance(params.n_genes, bool)
            or params.n_genes <= 0
        ):
            raise ValueError("n_genes must be a positive integer")
        if (
            not isinstance(params.n_regulatory_elements, int)
            or isinstance(params.n_regulatory_elements, bool)
            or params.n_regulatory_elements <= 0
        ):
            raise ValueError("n_regulatory_elements must be a positive integer")
        if (
            not isinstance(params.bitstream_length, int)
            or isinstance(params.bitstream_length, bool)
            or params.bitstream_length <= 0
        ):
            raise ValueError("bitstream_length must be a positive integer")
        for field_name in (
            "transcription_rate",
            "translation_rate",
            "degradation_rate",
            "methylation_rate",
            "demethylation_rate",
            "histone_mod_rate",
            "bioelectric_coupling",
            "neurochemical_coupling",
            "cellular_coupling",
        ):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value < 0.0:
                raise ValueError(f"{field_name} must be finite and non-negative")
        if (
            not math.isfinite(float(params.ciss_efficiency))
            or params.ciss_efficiency < 0.0
            or params.ciss_efficiency > 1.0
        ):
            raise ValueError("ciss_efficiency must be finite and within [0, 1]")
        if not math.isfinite(float(params.dna_chirality)) or params.dna_chirality not in (
            -1.0,
            1.0,
        ):
            raise ValueError("dna_chirality must be finite and either -1.0 or 1.0")
        if not math.isfinite(float(params.membrane_potential_rest)):
            raise ValueError("membrane_potential_rest must be finite")
        if params.rng_seed is not None:
            if isinstance(params.rng_seed, bool) or not isinstance(params.rng_seed, int):
                raise ValueError("rng_seed must be a non-negative integer or None")
            if params.rng_seed < 0:
                raise ValueError("rng_seed must be a non-negative integer or None")

    @classmethod
    def _validate_step_inputs(
        cls,
        dt: float,
        l2_input: Optional[Dict[str, Any]],
        bioelectric_signal: Optional[np.ndarray[Any, Any]],
        n_genes: int,
    ) -> None:
        if not math.isfinite(float(dt)) or dt <= 0.0:
            raise ValueError("dt must be finite and positive")
        if l2_input is not None and "second_messengers" in l2_input:
            cls._finite_mean(l2_input["second_messengers"], "second_messengers")
        if bioelectric_signal is not None:
            cls._bioelectric_signal(bioelectric_signal, n_genes)

    @staticmethod
    def _finite_mean(values: Any, name: str) -> float:
        arr = np.asarray(values, dtype=np.float64)
        if arr.size == 0 or not np.all(np.isfinite(arr)):
            raise ValueError(f"{name} must contain finite values")
        return float(np.mean(arr))

    @staticmethod
    def _bioelectric_signal(values: Any, n_genes: int) -> np.ndarray[Any, Any]:
        signal = np.asarray(values, dtype=np.float64)
        if signal.size == 1 and np.all(np.isfinite(signal)):
            return np.full(n_genes, float(signal.reshape(-1)[0]), dtype=np.float64)
        if signal.size != n_genes or not np.all(np.isfinite(signal)):
            raise ValueError("bioelectric_signal must be finite and scalar or one value per gene")
        return signal.reshape(n_genes)

step(dt, l2_input=None, bioelectric_signal=None)

Advance the layer by one time step.

Parameters:

Name Type Description Default
dt float

Time step in seconds.

 required
l2_input Optional[Dict[str, Any]]

Neurochemical layer output (second messengers).

 None
bioelectric_signal Optional[ndarray[Any, Any]]

External bioelectric modulation.

 None

Returns:

Type Description
Dict[str, Any]

Dict with expression, protein_levels, chromatin_state, output_bitstreams
Source code in src/sc_neurocore/scpn/layers/l3_genomic.py
Python
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
def step(
    self,
    dt: float,
    l2_input: Optional[Dict[str, Any]] = None,
    bioelectric_signal: Optional[np.ndarray[Any, Any]] = None,
) -> Dict[str, Any]:
    """
    Advance the layer by one time step.

    Args:
        dt: Time step in seconds.
        l2_input: Neurochemical layer output (second messengers).
        bioelectric_signal: External bioelectric modulation.

    Returns:
        Dict with expression, protein_levels, chromatin_state, output_bitstreams
    """
    self._validate_step_inputs(dt, l2_input, bioelectric_signal, self.params.n_genes)
    # 1. Update chromatin state (epigenetic dynamics)
    # Methylation silences genes
    demeth_prob = self.params.demethylation_rate * dt
    meth_prob = self.params.methylation_rate * dt

    demeth_mask = self._rng.random(self.params.n_genes) < demeth_prob
    meth_mask = self._rng.random(self.params.n_genes) < meth_prob

    self.methylation = np.where(demeth_mask, self.methylation * 0.9, self.methylation)
    self.methylation = np.where(meth_mask, self.methylation + 0.1, self.methylation)
    self.methylation = np.clip(self.methylation, 0.0, 1.0)

    # Chromatin openness inversely related to methylation
    self.chromatin_openness = (
        1.0 - self.methylation + self._rng.normal(0, 0.05, self.params.n_genes)
    )
    self.chromatin_openness = np.clip(self.chromatin_openness, 0.0, 1.0)

    # 2. Gene expression (stochastic transcription)
    # Only open chromatin can be transcribed
    transcription_prob = self.params.transcription_rate * self.chromatin_openness * dt
    transcription = self._rng.random(self.params.n_genes) < transcription_prob

    self.expression_levels = np.where(
        transcription,
        self.expression_levels + 0.1,
        self.expression_levels - self.params.degradation_rate * dt,
    )
    self.expression_levels = np.clip(self.expression_levels, 0.0, 1.0)

    # 3. Translation to protein
    translation_prob = self.params.translation_rate * self.expression_levels * dt
    translation = self._rng.random(self.params.n_genes) < translation_prob

    self.protein_levels = np.where(
        translation,
        self.protein_levels + 0.05,
        self.protein_levels - self.params.degradation_rate * dt * 0.5,
    )
    self.protein_levels = np.clip(self.protein_levels, 0.0, 1.0)

    # 4. CISS effect (quantum spin filtering)
    # Spin polarization depends on DNA chirality and electron flow
    electron_flow = np.mean(self.expression_levels)  # Proxy for metabolic activity
    ciss_baseline = self.params.ciss_efficiency * self.params.dna_chirality * electron_flow
    spin_polarization = np.clip(
        np.full(self.params.n_genes, ciss_baseline, dtype=np.float64)
        + self._rng.normal(0, 0.1, self.params.n_genes),
        -1.0,
        1.0,
    )
    self.spin_polarization[...] = spin_polarization

    # 5. Neurochemical coupling (L2 input modulates expression)
    if l2_input is not None and "second_messengers" in l2_input:
        # cAMP from second messengers activates transcription factors
        camp_level = self._finite_mean(l2_input["second_messengers"], "second_messengers")
        activation_boost = camp_level * self.params.neurochemical_coupling
        self.expression_levels += activation_boost * dt
        self.expression_levels = np.clip(self.expression_levels, 0.0, 1.0)

    # 6. Bioelectric pattern formation
    if bioelectric_signal is not None:
        signal = self._bioelectric_signal(bioelectric_signal, self.params.n_genes)
        self.membrane_potential = 0.9 * self.membrane_potential + 0.1 * signal
    # Internal bioelectric dynamics (gap junction diffusion)
    diffusion = np.roll(self.membrane_potential, 1) - self.membrane_potential
    self.membrane_potential += diffusion * self.params.bioelectric_coupling * dt

    # 7. Generate output bitstreams
    output_probs = np.clip(self.protein_levels, 0.0, 1.0)
    rands = self._rng.random((self.params.n_genes, self.params.bitstream_length))
    output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
    cellular_drive = self.params.cellular_coupling * self.protein_levels

    return {
        "expression_levels": self.expression_levels.copy(),
        "protein_levels": self.protein_levels.copy(),
        "chromatin_openness": self.chromatin_openness.copy(),
        "methylation": self.methylation.copy(),
        "spin_polarization": self.spin_polarization.copy(),
        "membrane_potential": self.membrane_potential.copy(),
        "cellular_drive": cellular_drive.copy(),
        "output_bitstreams": output_bitstreams,
    }

get_global_metric()

Return the global genomic activity metric.

Source code in src/sc_neurocore/scpn/layers/l3_genomic.py
Python
217
218
219
def get_global_metric(self) -> float:
    """Return the global genomic activity metric."""
    return float(np.mean(self.expression_levels))

get_ciss_coherence()

Return CISS spin coherence metric.

Source code in src/sc_neurocore/scpn/layers/l3_genomic.py
Python
221
222
223
def get_ciss_coherence(self) -> float:
    """Return CISS spin coherence metric."""
    return float(np.abs(np.mean(self.spin_polarization)))

L4_CellularLayer

Stochastic implementation of the Cellular-Tissue Synchronization Layer.

Models collective cellular behavior, gap junction coupling, and tissue-level pattern formation using bitstream representations.

Source code in src/sc_neurocore/scpn/layers/l4_cellular.py
Python
 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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
class L4_CellularLayer:
    """
    Stochastic implementation of the Cellular-Tissue Synchronization Layer.

    Models collective cellular behavior, gap junction coupling, and
    tissue-level pattern formation using bitstream representations.
    """

    def __init__(self, params: Optional[L4_StochasticParameters] = None):
        self.params = params or L4_StochasticParameters()
        self._validate_params(self.params)
        self._rng = np.random.default_rng(self.params.rng_seed)
        self.n_cells = self.params.grid_size[0] * self.params.grid_size[1]

        # Oscillator phases (Kuramoto-like model)
        self.phases = self._rng.random(self.n_cells) * 2 * np.pi

        # Oscillator amplitudes
        self.amplitudes = np.ones(self.n_cells) * 0.5

        # Calcium concentrations
        self.calcium = self._rng.random(self.n_cells) * 0.3

        # Gap junction states (0 = closed, 1 = open)
        self.gap_junctions = self._init_gap_junctions()

        # Tissue activity pattern
        self.activity_pattern = np.zeros(self.n_cells)

        # Build neighbor connectivity matrix
        self.neighbors = self._build_neighbor_matrix()

    def _init_gap_junctions(self) -> np.ndarray[Any, Any]:
        """Initialize gap junction connectivity."""
        # Random initial state with bias toward open
        return (self._rng.random(self.n_cells) > 0.3).astype(np.float32)

    def _build_neighbor_matrix(self) -> np.ndarray[Any, Any]:
        """Build 2D grid neighbor connectivity matrix."""
        h, w = self.params.grid_size
        n = self.n_cells
        neighbors = np.zeros((n, n), dtype=np.float32)

        for i in range(n):
            row, col = i // w, i % w
            # 4-connectivity (von Neumann)
            for dr, dc in [(-1, 0), (1, 0), (0, -1), (0, 1)]:
                nr, nc = row + dr, col + dc
                if 0 <= nr < h and 0 <= nc < w:
                    j = nr * w + nc
                    neighbors[i, j] = 1.0

        return neighbors

    def step(
        self,
        dt: float,
        l3_input: Optional[Dict[str, Any]] = None,
        external_stimulus: Optional[np.ndarray[Any, Any]] = None,
    ) -> Dict[str, Any]:
        """
        Advance the layer by one time step.

        Args:
            dt: Time step in seconds.
            l3_input: Genomic layer output (protein levels).
            external_stimulus: External stimulation pattern.

        Returns:
            Dict with phases, calcium, synchronization, output_bitstreams
        """
        self._validate_step_inputs(dt, l3_input, external_stimulus, self.n_cells)
        # 1. Kuramoto oscillator dynamics
        # dθ/dt = ω + K/N * Σ sin(θ_j - θ_i)
        phase_diffs = np.sin(self.phases[None, :] - self.phases[:, None])
        coupling_term = (
            self.params.coupling_strength
            * np.sum(self.neighbors * phase_diffs, axis=1)
            / np.maximum(np.sum(self.neighbors, axis=1), 1)
        )

        noise = self.params.noise_amplitude * self._rng.normal(0, 1, self.n_cells)

        self.phases += (2 * np.pi * self.params.natural_frequency + coupling_term + noise) * dt
        self.phases = self.phases % (2 * np.pi)

        # 2. Calcium wave dynamics
        # Diffusion via gap junctions
        ca_diff = np.zeros(self.n_cells)
        for i in range(self.n_cells):
            neighbor_indices = np.where(self.neighbors[i] > 0)[0]
            if len(neighbor_indices) > 0:
                # Diffusion weighted by gap junction state
                for j in neighbor_indices:
                    gj_state = (self.gap_junctions[i] + self.gap_junctions[j]) / 2
                    ca_diff[i] += (
                        self.params.gap_junction_conductance
                        * gj_state
                        * (self.calcium[j] - self.calcium[i])
                    )

        self.calcium += (
            self.params.ca_diffusion_rate * ca_diff - self.params.ca_decay_rate * self.calcium
        ) * dt

        # Calcium-induced calcium release (CICR)
        cicr_mask = self.calcium > self.params.ca_release_threshold
        self.calcium = np.where(cicr_mask, self.calcium + 0.2, self.calcium)
        self.calcium = np.clip(self.calcium, 0.0, 1.0)

        # 3. Gap junction dynamics
        # Gap junctions modulated by calcium and coupling
        gj_noise = self.params.gap_junction_noise * self._rng.normal(0, 1, self.n_cells)
        self.gap_junctions = np.clip(
            self.gap_junctions + gj_noise * dt + 0.1 * (1 - self.calcium) * dt, 0.0, 1.0
        )

        # 4. Genomic input coupling (L3 proteins modulate oscillators)
        if l3_input is not None and "protein_levels" in l3_input:
            protein_mean = self._finite_mean(l3_input["protein_levels"], "protein_levels")
            self.amplitudes = np.clip(
                self.amplitudes + protein_mean * self.params.genomic_coupling * dt, 0.1, 1.0
            )

        # 5. External stimulus
        if external_stimulus is not None:
            self.calcium += external_stimulus[: self.n_cells] * dt
            self.calcium = np.clip(self.calcium, 0.0, 1.0)

        # 6. Compute activity pattern
        self.activity_pattern = self.amplitudes * (1 + np.cos(self.phases)) / 2

        # 7. Compute synchronization order parameter
        order_parameter = float(np.abs(np.mean(np.exp(1j * self.phases))))

        # 8. Generate output bitstreams
        output_probs = np.clip(self.activity_pattern, 0.0, 1.0)
        rands = self._rng.random((self.n_cells, self.params.bitstream_length))
        output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
        organismal_drive = self.params.organismal_coupling * order_parameter

        return {
            "phases": self.phases.copy(),
            "amplitudes": self.amplitudes.copy(),
            "calcium": self.calcium.copy(),
            "gap_junctions": self.gap_junctions.copy(),
            "activity_pattern": self.activity_pattern.copy(),
            "synchronization": order_parameter,
            "organismal_drive": organismal_drive,
            "output_bitstreams": output_bitstreams,
        }

    def get_global_metric(self) -> float:
        """Return the global synchronization metric (Kuramoto order parameter)."""
        return float(np.abs(np.mean(np.exp(1j * self.phases))))

    def get_tissue_pattern(self) -> np.ndarray[Any, Any]:
        """Return 2D tissue activity pattern."""
        return self.activity_pattern.reshape(self.params.grid_size)

    @staticmethod
    def _validate_params(params: L4_StochasticParameters) -> None:
        if (
            not isinstance(params.grid_size, tuple)
            or len(params.grid_size) != 2
            or any(
                not isinstance(dim, int) or isinstance(dim, bool) or dim <= 0
                for dim in params.grid_size
            )
        ):
            raise ValueError("grid_size must be a tuple of two positive integers")
        if (
            not isinstance(params.bitstream_length, int)
            or isinstance(params.bitstream_length, bool)
            or params.bitstream_length <= 0
        ):
            raise ValueError("bitstream_length must be a positive integer")
        if not math.isfinite(float(params.natural_frequency)) or params.natural_frequency <= 0.0:
            raise ValueError("natural_frequency must be finite and positive")
        for field_name in (
            "coupling_strength",
            "noise_amplitude",
            "gap_junction_noise",
            "ca_diffusion_rate",
            "ca_decay_rate",
            "genomic_coupling",
            "organismal_coupling",
        ):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value < 0.0:
                raise ValueError(f"{field_name} must be finite and non-negative")
        if (
            not math.isfinite(float(params.gap_junction_conductance))
            or params.gap_junction_conductance < 0.0
            or params.gap_junction_conductance > 1.0
        ):
            raise ValueError("gap_junction_conductance must be finite and within [0, 1]")
        if (
            not math.isfinite(float(params.ca_release_threshold))
            or params.ca_release_threshold < 0.0
            or params.ca_release_threshold > 1.0
        ):
            raise ValueError("ca_release_threshold must be finite and within [0, 1]")
        if params.rng_seed is not None:
            if isinstance(params.rng_seed, bool) or not isinstance(params.rng_seed, int):
                raise ValueError("rng_seed must be a non-negative integer or None")
            if params.rng_seed < 0:
                raise ValueError("rng_seed must be a non-negative integer or None")

    @classmethod
    def _validate_step_inputs(
        cls,
        dt: float,
        l3_input: Optional[Dict[str, Any]],
        external_stimulus: Optional[np.ndarray[Any, Any]],
        n_cells: int,
    ) -> None:
        if not math.isfinite(float(dt)) or dt <= 0.0:
            raise ValueError("dt must be finite and positive")
        if l3_input is not None and "protein_levels" in l3_input:
            cls._finite_mean(l3_input["protein_levels"], "protein_levels")
        if external_stimulus is not None:
            stimulus = np.asarray(external_stimulus, dtype=np.float64)
            if stimulus.size != n_cells or not np.all(np.isfinite(stimulus)):
                raise ValueError("external_stimulus must contain one finite value per cell")

    @staticmethod
    def _finite_mean(values: Any, name: str) -> float:
        arr = np.asarray(values, dtype=np.float64)
        if arr.size == 0 or not np.all(np.isfinite(arr)):
            raise ValueError(f"{name} must contain finite values")
        return float(np.mean(arr))

step(dt, l3_input=None, external_stimulus=None)

Advance the layer by one time step.

Parameters:

Name Type Description Default
dt float

Time step in seconds.

 required
l3_input Optional[Dict[str, Any]]

Genomic layer output (protein levels).

 None
external_stimulus Optional[ndarray[Any, Any]]

External stimulation pattern.

 None

Returns:

Type Description
Dict[str, Any]

Dict with phases, calcium, synchronization, output_bitstreams
Source code in src/sc_neurocore/scpn/layers/l4_cellular.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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
def step(
    self,
    dt: float,
    l3_input: Optional[Dict[str, Any]] = None,
    external_stimulus: Optional[np.ndarray[Any, Any]] = None,
) -> Dict[str, Any]:
    """
    Advance the layer by one time step.

    Args:
        dt: Time step in seconds.
        l3_input: Genomic layer output (protein levels).
        external_stimulus: External stimulation pattern.

    Returns:
        Dict with phases, calcium, synchronization, output_bitstreams
    """
    self._validate_step_inputs(dt, l3_input, external_stimulus, self.n_cells)
    # 1. Kuramoto oscillator dynamics
    # dθ/dt = ω + K/N * Σ sin(θ_j - θ_i)
    phase_diffs = np.sin(self.phases[None, :] - self.phases[:, None])
    coupling_term = (
        self.params.coupling_strength
        * np.sum(self.neighbors * phase_diffs, axis=1)
        / np.maximum(np.sum(self.neighbors, axis=1), 1)
    )

    noise = self.params.noise_amplitude * self._rng.normal(0, 1, self.n_cells)

    self.phases += (2 * np.pi * self.params.natural_frequency + coupling_term + noise) * dt
    self.phases = self.phases % (2 * np.pi)

    # 2. Calcium wave dynamics
    # Diffusion via gap junctions
    ca_diff = np.zeros(self.n_cells)
    for i in range(self.n_cells):
        neighbor_indices = np.where(self.neighbors[i] > 0)[0]
        if len(neighbor_indices) > 0:
            # Diffusion weighted by gap junction state
            for j in neighbor_indices:
                gj_state = (self.gap_junctions[i] + self.gap_junctions[j]) / 2
                ca_diff[i] += (
                    self.params.gap_junction_conductance
                    * gj_state
                    * (self.calcium[j] - self.calcium[i])
                )

    self.calcium += (
        self.params.ca_diffusion_rate * ca_diff - self.params.ca_decay_rate * self.calcium
    ) * dt

    # Calcium-induced calcium release (CICR)
    cicr_mask = self.calcium > self.params.ca_release_threshold
    self.calcium = np.where(cicr_mask, self.calcium + 0.2, self.calcium)
    self.calcium = np.clip(self.calcium, 0.0, 1.0)

    # 3. Gap junction dynamics
    # Gap junctions modulated by calcium and coupling
    gj_noise = self.params.gap_junction_noise * self._rng.normal(0, 1, self.n_cells)
    self.gap_junctions = np.clip(
        self.gap_junctions + gj_noise * dt + 0.1 * (1 - self.calcium) * dt, 0.0, 1.0
    )

    # 4. Genomic input coupling (L3 proteins modulate oscillators)
    if l3_input is not None and "protein_levels" in l3_input:
        protein_mean = self._finite_mean(l3_input["protein_levels"], "protein_levels")
        self.amplitudes = np.clip(
            self.amplitudes + protein_mean * self.params.genomic_coupling * dt, 0.1, 1.0
        )

    # 5. External stimulus
    if external_stimulus is not None:
        self.calcium += external_stimulus[: self.n_cells] * dt
        self.calcium = np.clip(self.calcium, 0.0, 1.0)

    # 6. Compute activity pattern
    self.activity_pattern = self.amplitudes * (1 + np.cos(self.phases)) / 2

    # 7. Compute synchronization order parameter
    order_parameter = float(np.abs(np.mean(np.exp(1j * self.phases))))

    # 8. Generate output bitstreams
    output_probs = np.clip(self.activity_pattern, 0.0, 1.0)
    rands = self._rng.random((self.n_cells, self.params.bitstream_length))
    output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
    organismal_drive = self.params.organismal_coupling * order_parameter

    return {
        "phases": self.phases.copy(),
        "amplitudes": self.amplitudes.copy(),
        "calcium": self.calcium.copy(),
        "gap_junctions": self.gap_junctions.copy(),
        "activity_pattern": self.activity_pattern.copy(),
        "synchronization": order_parameter,
        "organismal_drive": organismal_drive,
        "output_bitstreams": output_bitstreams,
    }

get_global_metric()

Return the global synchronization metric (Kuramoto order parameter).

Source code in src/sc_neurocore/scpn/layers/l4_cellular.py
Python
216
217
218
def get_global_metric(self) -> float:
    """Return the global synchronization metric (Kuramoto order parameter)."""
    return float(np.abs(np.mean(np.exp(1j * self.phases))))

get_tissue_pattern()

Return 2D tissue activity pattern.

Source code in src/sc_neurocore/scpn/layers/l4_cellular.py
Python
220
221
222
def get_tissue_pattern(self) -> np.ndarray[Any, Any]:
    """Return 2D tissue activity pattern."""
    return self.activity_pattern.reshape(self.params.grid_size)

L5_OrganismalLayer

Stochastic implementation of the Organismal-Psychoemotional Layer.

Models whole-organism integration, autonomic regulation, and emotional dynamics using bitstream representations.

Source code in src/sc_neurocore/scpn/layers/l5_organismal.py
Python
 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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
class L5_OrganismalLayer:
    """
    Stochastic implementation of the Organismal-Psychoemotional Layer.

    Models whole-organism integration, autonomic regulation, and
    emotional dynamics using bitstream representations.
    """

    # Emotional dimension indices
    VALENCE = 0  # Pleasant-Unpleasant
    AROUSAL = 1  # Activated-Deactivated
    DOMINANCE = 2  # Dominant-Submissive
    APPROACH = 3  # Approach-Avoid
    CERTAINTY = 4  # Certain-Uncertain
    ATTENTION = 5  # Focused-Diffuse
    FAIRNESS = 6  # Fair-Unfair
    SAFETY = 7  # Safe-Threatened

    def __init__(self, params: L5_StochasticParameters | None = None) -> None:
        self.params = params or L5_StochasticParameters()
        self._validate_params(self.params)
        self._rng = np.random.default_rng(self.params.rng_seed)

        # Emotional state vector
        self.emotional_state = np.array([0.5, 0.3, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6])

        # Autonomic nervous system state
        self.sympathetic = self.params.sympathetic_baseline
        self.parasympathetic = self.params.parasympathetic_baseline

        # Heart dynamics
        self.heart_rate = self.params.base_heart_rate
        self.hrv_phase = 0.0
        self.rr_intervals: list[float] = []

        # Interoceptive state (body sense)
        self.interoceptive_state = self._rng.random(self.params.n_autonomic_nodes) * 0.3

        # Attractor basins for emotional states
        self.attractors = self._init_emotional_attractors()

        # Time tracking
        self.time = 0.0

    def _init_emotional_attractors(self) -> np.ndarray:
        """Initialize emotional attractor states."""
        # Define stable emotional configurations
        attractors = np.array(
            [
                [0.8, 0.3, 0.6, 0.7, 0.7, 0.5, 0.6, 0.8],  # Joy/contentment
                [0.2, 0.8, 0.3, 0.2, 0.3, 0.8, 0.3, 0.2],  # Fear/anxiety
                [0.2, 0.7, 0.7, 0.8, 0.6, 0.7, 0.2, 0.4],  # Anger
                [0.3, 0.2, 0.2, 0.2, 0.4, 0.3, 0.5, 0.5],  # Sadness
                [0.5, 0.4, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6],  # Neutral
            ]
        )
        return attractors

    def step(
        self,
        dt: float,
        l4_input: dict[str, Any] | None = None,
        external_event: dict[str, Any] | None = None,
    ) -> dict[str, Any]:
        """
        Advance the layer by one time step.

        Args:
            dt: Time step in seconds.
            l4_input: Cellular layer output (synchronization).
            external_event: External emotional trigger {valence, arousal, ...}.

        Returns:
            Dict with emotional_state, autonomic, heart_rate, output_bitstreams
        """
        self._validate_step_inputs(dt, l4_input, external_event)
        self.time += dt

        # 1. Process external emotional events
        if external_event is not None:
            self._apply_external_event(external_event)

        # 2. Attractor dynamics (emotional states converge to stable patterns)
        # Find nearest attractor
        distances = np.linalg.norm(self.attractors - self.emotional_state, axis=1)
        nearest_attractor = self.attractors[np.argmin(distances)]

        # Pull toward attractor
        self.emotional_state += (
            self.params.attractor_strength * (nearest_attractor - self.emotional_state) * dt
        )

        # Add noise
        self.emotional_state += (
            self.params.emotional_noise * self._rng.normal(0, 1, self.params.n_emotional_dims) * dt
        )

        # Decay toward baseline
        baseline = np.array([0.5, 0.3, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6])
        self.emotional_state += self.params.emotional_decay * (baseline - self.emotional_state) * dt

        self.emotional_state = np.clip(self.emotional_state, 0.0, 1.0).astype(
            np.float64, copy=False
        )

        # 3. Autonomic nervous system dynamics
        # Sympathetic driven by arousal and threat
        target_symp = (
            self.emotional_state[self.AROUSAL] * 0.5 + (1 - self.emotional_state[self.SAFETY]) * 0.5
        )
        # Parasympathetic driven by valence and safety
        target_para = (
            self.emotional_state[self.VALENCE] * 0.3 + self.emotional_state[self.SAFETY] * 0.7
        )

        tau = self.params.autonomic_time_constant
        self.sympathetic += (target_symp - self.sympathetic) * dt / tau
        self.parasympathetic += (target_para - self.parasympathetic) * dt / tau

        self.sympathetic = np.clip(self.sympathetic, 0.0, 1.0)
        self.parasympathetic = np.clip(self.parasympathetic, 0.0, 1.0)

        # 4. Heart rate and HRV
        # RSA (Respiratory Sinus Arrhythmia)
        self.hrv_phase += 2 * np.pi * self.params.respiratory_frequency * dt
        rsa_component = self.params.hrv_amplitude * np.sin(self.hrv_phase) * self.parasympathetic

        # Sympathetic raises HR, parasympathetic lowers it
        target_hr = self.params.base_heart_rate + 20 * self.sympathetic - 15 * self.parasympathetic
        self.heart_rate += (target_hr - self.heart_rate) * dt * 0.5
        self.heart_rate += rsa_component * 10  # RSA effect

        # Track RR intervals
        rr = 60000.0 / self.heart_rate  # ms
        self.rr_intervals.append(rr)
        if len(self.rr_intervals) > 100:
            self.rr_intervals.pop(0)

        # 5. Cellular input coupling (L4 synchronization affects coherence)
        if l4_input is not None and "synchronization" in l4_input:
            sync = float(l4_input["synchronization"])
            # High cellular sync improves emotional stability
            self.emotional_state[self.CERTAINTY] += sync * self.params.cellular_coupling * dt
            self.emotional_state = np.clip(self.emotional_state, 0.0, 1.0).astype(
                np.float64, copy=False
            )

        # 6. Update interoceptive state
        self.interoceptive_state = (
            0.8 * self.interoceptive_state
            + 0.2
            * np.tile(
                [self.sympathetic, self.parasympathetic, self.heart_rate / 100],
                self.params.n_autonomic_nodes // 3 + 1,
            )[: self.params.n_autonomic_nodes]
        )

        # 7. Generate output bitstreams
        output_probs = np.concatenate(
            [self.emotional_state, [self.sympathetic, self.parasympathetic, self.heart_rate / 100]]
        )
        output_probs = np.tile(output_probs, self.params.n_autonomic_nodes // len(output_probs) + 1)
        output_probs = output_probs[: self.params.n_autonomic_nodes]
        output_probs = np.clip(output_probs, 0.0, 1.0)

        rands = self._rng.random((self.params.n_autonomic_nodes, self.params.bitstream_length))
        output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
        ecological_drive = self.params.ecological_coupling * self.emotional_state

        return {
            "emotional_state": self.emotional_state.copy(),
            "sympathetic": self.sympathetic,
            "parasympathetic": self.parasympathetic,
            "heart_rate": self.heart_rate,
            "hrv_rmssd": self._compute_rmssd(),
            "interoceptive_state": self.interoceptive_state.copy(),
            "ecological_drive": ecological_drive.copy(),
            "output_bitstreams": output_bitstreams,
        }

    def _compute_rmssd(self) -> float:
        """Compute RMSSD (root mean square of successive differences)."""
        if len(self.rr_intervals) < 2:
            return 0.0
        rr = np.array(self.rr_intervals)
        diff = np.diff(rr)
        return float(np.sqrt(np.mean(diff**2)))

    def get_global_metric(self) -> float:
        """Return the global organismal coherence metric."""
        # Combine HRV coherence with emotional stability
        hrv_coherence = self._compute_rmssd() / 100  # Normalize
        emotional_stability = 1.0 - np.std(self.emotional_state)
        return float(0.5 * hrv_coherence + 0.5 * emotional_stability)

    def get_emotional_valence(self) -> float:
        """Return current emotional valence."""
        return float(self.emotional_state[self.VALENCE])

    @staticmethod
    def _validate_params(params: L5_StochasticParameters) -> None:
        if params.n_emotional_dims != 8:
            raise ValueError("n_emotional_dims must be 8 for the fixed L5 emotional basis")
        if (
            not isinstance(params.n_autonomic_nodes, int)
            or isinstance(params.n_autonomic_nodes, bool)
            or params.n_autonomic_nodes <= 0
        ):
            raise ValueError("n_autonomic_nodes must be a positive integer")
        if (
            not isinstance(params.bitstream_length, int)
            or isinstance(params.bitstream_length, bool)
            or params.bitstream_length <= 0
        ):
            raise ValueError("bitstream_length must be a positive integer")
        for field_name in ("sympathetic_baseline", "parasympathetic_baseline"):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value < 0.0 or value > 1.0:
                raise ValueError(f"{field_name} must be finite and within [0, 1]")
        for field_name in ("autonomic_time_constant", "base_heart_rate", "respiratory_frequency"):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value <= 0.0:
                raise ValueError(f"{field_name} must be finite and positive")
        for field_name in (
            "hrv_amplitude",
            "emotional_decay",
            "emotional_noise",
            "attractor_strength",
            "cellular_coupling",
            "ecological_coupling",
        ):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value < 0.0:
                raise ValueError(f"{field_name} must be finite and non-negative")
        if params.rng_seed is not None:
            if isinstance(params.rng_seed, bool) or not isinstance(params.rng_seed, int):
                raise ValueError("rng_seed must be a non-negative integer or None")
            if params.rng_seed < 0:
                raise ValueError("rng_seed must be a non-negative integer or None")

    @classmethod
    def _validate_step_inputs(
        cls,
        dt: float,
        l4_input: dict[str, Any] | None,
        external_event: dict[str, Any] | None,
    ) -> None:
        if not math.isfinite(float(dt)) or dt <= 0.0:
            raise ValueError("dt must be finite and positive")
        if l4_input is not None and "synchronization" in l4_input:
            sync = float(l4_input["synchronization"])
            if not math.isfinite(sync) or sync < 0.0 or sync > 1.0:
                raise ValueError("synchronization must be finite and within [0, 1]")
        if external_event is not None:
            cls._validate_external_event(external_event)

    @classmethod
    def _validate_external_event(cls, external_event: dict[str, Any]) -> None:
        dimension_names = cls._dimension_names()
        for key, value in external_event.items():
            if key == "type":
                event_type = str(value)
                if event_type not in {"stress", "calm", "reward", "threat"}:
                    raise ValueError("external_event type must be stress, calm, reward, or threat")
                continue
            if key == "intensity":
                intensity = float(value)
                if not math.isfinite(intensity) or intensity < 0.0 or intensity > 1.0:
                    raise ValueError("intensity must be finite and within [0, 1]")
                continue
            if isinstance(key, int) and not isinstance(key, bool):
                if key < 0 or key >= 8 or not math.isfinite(float(value)):
                    raise ValueError("external_event dimension values must be finite")
                continue
            if isinstance(key, str) and key in dimension_names:
                if not math.isfinite(float(value)):
                    raise ValueError("external_event dimension values must be finite")
                continue
            raise ValueError(
                "external_event keys must be dimension indices, dimension names, or type/intensity"
            )

    def _apply_external_event(self, external_event: dict[str, Any]) -> None:
        event_type = str(external_event.get("type", ""))
        intensity = float(external_event.get("intensity", 1.0))
        if event_type in {"stress", "threat"}:
            self.emotional_state[self.VALENCE] -= 0.2 * intensity
            self.emotional_state[self.AROUSAL] += 0.4 * intensity
            self.emotional_state[self.SAFETY] -= 0.5 * intensity
            self.emotional_state[self.ATTENTION] += 0.2 * intensity
        elif event_type == "calm":
            self.emotional_state[self.VALENCE] += 0.2 * intensity
            self.emotional_state[self.AROUSAL] -= 0.25 * intensity
            self.emotional_state[self.SAFETY] += 0.4 * intensity
        elif event_type == "reward":
            self.emotional_state[self.VALENCE] += 0.3 * intensity
            self.emotional_state[self.APPROACH] += 0.25 * intensity
            self.emotional_state[self.AROUSAL] += 0.1 * intensity

        dimension_names = self._dimension_names()
        for key, value in external_event.items():
            if key in {"type", "intensity"}:
                continue
            if isinstance(key, int) and not isinstance(key, bool):
                self.emotional_state[key] += float(value) * 0.3
            elif isinstance(key, str) and key in dimension_names:
                self.emotional_state[dimension_names[key]] += float(value) * 0.3

    @classmethod
    def _dimension_names(cls) -> dict[str, int]:
        return {
            "valence": cls.VALENCE,
            "arousal": cls.AROUSAL,
            "dominance": cls.DOMINANCE,
            "approach": cls.APPROACH,
            "certainty": cls.CERTAINTY,
            "attention": cls.ATTENTION,
            "fairness": cls.FAIRNESS,
            "safety": cls.SAFETY,
        }

step(dt, l4_input=None, external_event=None)

Advance the layer by one time step.

Parameters:

Name Type Description Default
dt float

Time step in seconds.

 required
l4_input dict[str, Any] | None

Cellular layer output (synchronization).

 None
external_event dict[str, Any] | None

External emotional trigger {valence, arousal, ...}.

 None

Returns:

Type Description
dict[str, Any]

Dict with emotional_state, autonomic, heart_rate, output_bitstreams
Source code in src/sc_neurocore/scpn/layers/l5_organismal.py
Python
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
def step(
    self,
    dt: float,
    l4_input: dict[str, Any] | None = None,
    external_event: dict[str, Any] | None = None,
) -> dict[str, Any]:
    """
    Advance the layer by one time step.

    Args:
        dt: Time step in seconds.
        l4_input: Cellular layer output (synchronization).
        external_event: External emotional trigger {valence, arousal, ...}.

    Returns:
        Dict with emotional_state, autonomic, heart_rate, output_bitstreams
    """
    self._validate_step_inputs(dt, l4_input, external_event)
    self.time += dt

    # 1. Process external emotional events
    if external_event is not None:
        self._apply_external_event(external_event)

    # 2. Attractor dynamics (emotional states converge to stable patterns)
    # Find nearest attractor
    distances = np.linalg.norm(self.attractors - self.emotional_state, axis=1)
    nearest_attractor = self.attractors[np.argmin(distances)]

    # Pull toward attractor
    self.emotional_state += (
        self.params.attractor_strength * (nearest_attractor - self.emotional_state) * dt
    )

    # Add noise
    self.emotional_state += (
        self.params.emotional_noise * self._rng.normal(0, 1, self.params.n_emotional_dims) * dt
    )

    # Decay toward baseline
    baseline = np.array([0.5, 0.3, 0.5, 0.5, 0.5, 0.5, 0.5, 0.6])
    self.emotional_state += self.params.emotional_decay * (baseline - self.emotional_state) * dt

    self.emotional_state = np.clip(self.emotional_state, 0.0, 1.0).astype(
        np.float64, copy=False
    )

    # 3. Autonomic nervous system dynamics
    # Sympathetic driven by arousal and threat
    target_symp = (
        self.emotional_state[self.AROUSAL] * 0.5 + (1 - self.emotional_state[self.SAFETY]) * 0.5
    )
    # Parasympathetic driven by valence and safety
    target_para = (
        self.emotional_state[self.VALENCE] * 0.3 + self.emotional_state[self.SAFETY] * 0.7
    )

    tau = self.params.autonomic_time_constant
    self.sympathetic += (target_symp - self.sympathetic) * dt / tau
    self.parasympathetic += (target_para - self.parasympathetic) * dt / tau

    self.sympathetic = np.clip(self.sympathetic, 0.0, 1.0)
    self.parasympathetic = np.clip(self.parasympathetic, 0.0, 1.0)

    # 4. Heart rate and HRV
    # RSA (Respiratory Sinus Arrhythmia)
    self.hrv_phase += 2 * np.pi * self.params.respiratory_frequency * dt
    rsa_component = self.params.hrv_amplitude * np.sin(self.hrv_phase) * self.parasympathetic

    # Sympathetic raises HR, parasympathetic lowers it
    target_hr = self.params.base_heart_rate + 20 * self.sympathetic - 15 * self.parasympathetic
    self.heart_rate += (target_hr - self.heart_rate) * dt * 0.5
    self.heart_rate += rsa_component * 10  # RSA effect

    # Track RR intervals
    rr = 60000.0 / self.heart_rate  # ms
    self.rr_intervals.append(rr)
    if len(self.rr_intervals) > 100:
        self.rr_intervals.pop(0)

    # 5. Cellular input coupling (L4 synchronization affects coherence)
    if l4_input is not None and "synchronization" in l4_input:
        sync = float(l4_input["synchronization"])
        # High cellular sync improves emotional stability
        self.emotional_state[self.CERTAINTY] += sync * self.params.cellular_coupling * dt
        self.emotional_state = np.clip(self.emotional_state, 0.0, 1.0).astype(
            np.float64, copy=False
        )

    # 6. Update interoceptive state
    self.interoceptive_state = (
        0.8 * self.interoceptive_state
        + 0.2
        * np.tile(
            [self.sympathetic, self.parasympathetic, self.heart_rate / 100],
            self.params.n_autonomic_nodes // 3 + 1,
        )[: self.params.n_autonomic_nodes]
    )

    # 7. Generate output bitstreams
    output_probs = np.concatenate(
        [self.emotional_state, [self.sympathetic, self.parasympathetic, self.heart_rate / 100]]
    )
    output_probs = np.tile(output_probs, self.params.n_autonomic_nodes // len(output_probs) + 1)
    output_probs = output_probs[: self.params.n_autonomic_nodes]
    output_probs = np.clip(output_probs, 0.0, 1.0)

    rands = self._rng.random((self.params.n_autonomic_nodes, self.params.bitstream_length))
    output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
    ecological_drive = self.params.ecological_coupling * self.emotional_state

    return {
        "emotional_state": self.emotional_state.copy(),
        "sympathetic": self.sympathetic,
        "parasympathetic": self.parasympathetic,
        "heart_rate": self.heart_rate,
        "hrv_rmssd": self._compute_rmssd(),
        "interoceptive_state": self.interoceptive_state.copy(),
        "ecological_drive": ecological_drive.copy(),
        "output_bitstreams": output_bitstreams,
    }

get_global_metric()

Return the global organismal coherence metric.

Source code in src/sc_neurocore/scpn/layers/l5_organismal.py
Python
254
255
256
257
258
259
def get_global_metric(self) -> float:
    """Return the global organismal coherence metric."""
    # Combine HRV coherence with emotional stability
    hrv_coherence = self._compute_rmssd() / 100  # Normalize
    emotional_stability = 1.0 - np.std(self.emotional_state)
    return float(0.5 * hrv_coherence + 0.5 * emotional_stability)

get_emotional_valence()

Return current emotional valence.

Source code in src/sc_neurocore/scpn/layers/l5_organismal.py
Python
261
262
263
def get_emotional_valence(self) -> float:
    """Return current emotional valence."""
    return float(self.emotional_state[self.VALENCE])

L6_EcologicalLayer

Stochastic implementation of the Ecological-Planetary Layer.

Models planetary-scale electromagnetic fields, Schumann resonances, and biospheric network dynamics using bitstream representations.

Source code in src/sc_neurocore/scpn/layers/l6_ecological.py
Python
 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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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 L6_EcologicalLayer:
    """
    Stochastic implementation of the Ecological-Planetary Layer.

    Models planetary-scale electromagnetic fields, Schumann resonances,
    and biospheric network dynamics using bitstream representations.
    """

    def __init__(self, params: Optional[L6_StochasticParameters] = None):
        self.params = params or L6_StochasticParameters()
        self._validate_params(self.params)
        self._rng = np.random.default_rng(self.params.rng_seed)

        # Schumann resonance field (superposition of modes)
        self.schumann_phases = np.zeros(len(self.params.schumann_frequencies))
        self.schumann_amplitudes = np.ones(len(self.params.schumann_frequencies))

        # Geomagnetic field state
        self.geomag_field = np.ones(self.params.n_field_nodes) * self.params.geomag_baseline

        # Circadian phase
        self.circadian_phase = 0.0

        # Biospheric network state (collective field)
        self.biospheric_field = self._rng.random(self.params.n_field_nodes) * 0.3

        # Planetary consciousness coherence
        self.planetary_coherence = 0.5

        # History for temporal patterns
        self.history: List[Dict[str, Any]] = []

        # Time tracking
        self.time = 0.0

    def step(
        self,
        dt: float,
        l5_input: Optional[Dict[str, Any]] = None,
        solar_activity: float = 0.5,
        lunar_phase: float = 0.0,
    ) -> Dict[str, Any]:
        """
        Advance the layer by one time step.

        Args:
            dt: Time step in seconds.
            l5_input: Organismal layer output (emotional coherence).
            solar_activity: Solar activity index (0-1).
            lunar_phase: Lunar phase (0 to 2π).

        Returns:
            Dict with schumann_field, geomag, circadian, output_bitstreams
        """
        self._validate_step_inputs(dt, l5_input, solar_activity, lunar_phase)
        self.time += dt

        # 1. Schumann resonance dynamics
        for i, freq in enumerate(self.params.schumann_frequencies):
            self.schumann_phases[i] += 2 * np.pi * freq * dt
            self.schumann_phases[i] = self.schumann_phases[i] % (2 * np.pi)

        # Compute Schumann field as superposition
        schumann_signal = np.zeros(self.params.n_field_nodes)
        for i, freq in enumerate(self.params.schumann_frequencies):
            spatial_pattern = np.sin(np.linspace(0, 2 * np.pi * (i + 1), self.params.n_field_nodes))
            schumann_signal += (
                self.schumann_amplitudes[i]
                * self.params.schumann_amplitude
                * np.cos(self.schumann_phases[i])
                * spatial_pattern
            )

        # Add noise
        schumann_signal += self.params.schumann_noise * self._rng.normal(
            0, 1, self.params.n_field_nodes
        )

        # Normalize to [0, 1]
        schumann_field = (schumann_signal - schumann_signal.min()) / (
            schumann_signal.max() - schumann_signal.min() + 1e-8
        )

        # 2. Geomagnetic field dynamics
        # Solar activity modulates geomagnetic storms
        storm_factor = 1.0 + 0.5 * (solar_activity - 0.5)
        geomag_variation = (
            self.params.geomag_variation
            * storm_factor
            * self._rng.normal(0, 1, self.params.n_field_nodes)
        )
        self.geomag_field = np.clip(
            self.geomag_field + geomag_variation * dt,
            self.params.geomag_baseline * 0.5,
            self.params.geomag_baseline * 1.5,
        )

        # 3. Circadian rhythm
        self.circadian_phase += 2 * np.pi * dt / self.params.circadian_period
        self.circadian_phase = self.circadian_phase % (2 * np.pi)
        circadian_signal = 0.5 + self.params.circadian_amplitude * np.cos(self.circadian_phase)

        # 4. Biospheric network dynamics
        # Coupling between nodes
        network_coupling = np.zeros(self.params.n_field_nodes)
        for i in range(self.params.n_field_nodes):
            neighbors = [(i - 1) % self.params.n_field_nodes, (i + 1) % self.params.n_field_nodes]
            neighbor_mean = np.mean([self.biospheric_field[j] for j in neighbors])
            network_coupling[i] = neighbor_mean - self.biospheric_field[i]

        self.biospheric_field += (
            self.params.network_coupling * network_coupling
            + self.params.network_noise * self._rng.normal(0, 1, self.params.n_field_nodes)
        ) * dt

        # Modulate by Schumann and circadian
        self.biospheric_field *= (0.9 + 0.1 * schumann_field) * (0.8 + 0.2 * circadian_signal)
        self.biospheric_field = np.clip(self.biospheric_field, 0.0, 1.0)

        # 5. Organismal coupling (L5 collective emotional state affects field)
        if l5_input is not None:
            l5_effect = self._l5_organismal_effect(l5_input)
            if l5_effect != 0.0:
                self.biospheric_field += self.params.organismal_coupling * l5_effect * dt
                self.biospheric_field = np.clip(self.biospheric_field, 0.0, 1.0)

        # 6. Lunar phase modulation
        lunar_factor = 0.5 + 0.5 * np.cos(lunar_phase)
        self.schumann_amplitudes = np.ones(len(self.params.schumann_frequencies)) * (
            0.8 + 0.2 * lunar_factor
        )

        # 7. Compute planetary coherence
        self.planetary_coherence = float(
            np.abs(np.mean(np.exp(1j * 2 * np.pi * self.biospheric_field)))
        )

        # 8. Generate output bitstreams
        output_probs = self.biospheric_field * circadian_signal
        rands = self._rng.random((self.params.n_field_nodes, self.params.bitstream_length))
        output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
        symbolic_drive = self.params.symbolic_coupling * schumann_field

        # Store history
        result = {
            "schumann_field": schumann_field,
            "schumann_phases": self.schumann_phases.copy(),
            "geomag_field": self.geomag_field.copy(),
            "circadian_phase": self.circadian_phase,
            "circadian_signal": circadian_signal,
            "biospheric_field": self.biospheric_field.copy(),
            "planetary_coherence": self.planetary_coherence,
            "symbolic_drive": symbolic_drive.copy(),
            "output_bitstreams": output_bitstreams,
        }

        self.history.append(
            {
                "time": self.time,
                "coherence": self.planetary_coherence,
                "schumann_power": float(np.mean(schumann_field**2)),
            }
        )
        if len(self.history) > 100:
            self.history.pop(0)

        return result

    def get_global_metric(self) -> float:
        """Return the global planetary coherence metric."""
        return self.planetary_coherence

    def get_schumann_spectrum(self) -> Dict[float, float]:
        """Return current Schumann resonance spectrum."""
        return {
            freq: float(amp * np.cos(phase))
            for freq, amp, phase in zip(
                self.params.schumann_frequencies, self.schumann_amplitudes, self.schumann_phases
            )
        }

    def get_circadian_time(self) -> float:
        """Return current circadian time (0-24 hours)."""
        return (self.circadian_phase / (2 * np.pi)) * 24.0

    @staticmethod
    def _validate_params(params: L6_StochasticParameters) -> None:
        if (
            not isinstance(params.n_field_nodes, int)
            or isinstance(params.n_field_nodes, bool)
            or params.n_field_nodes <= 0
        ):
            raise ValueError("n_field_nodes must be a positive integer")
        if (
            not isinstance(params.bitstream_length, int)
            or isinstance(params.bitstream_length, bool)
            or params.bitstream_length <= 0
        ):
            raise ValueError("bitstream_length must be a positive integer")
        if (
            not isinstance(params.schumann_frequencies, tuple)
            or len(params.schumann_frequencies) == 0
            or any(
                not math.isfinite(float(freq)) or float(freq) <= 0.0
                for freq in params.schumann_frequencies
            )
        ):
            raise ValueError(
                "schumann_frequencies must be a non-empty tuple of positive finite values"
            )
        for field_name in (
            "schumann_amplitude",
            "schumann_noise",
            "geomag_variation",
            "circadian_amplitude",
            "network_coupling",
            "network_noise",
            "organismal_coupling",
            "symbolic_coupling",
        ):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value < 0.0:
                raise ValueError(f"{field_name} must be finite and non-negative")
        if not math.isfinite(float(params.geomag_baseline)) or params.geomag_baseline <= 0.0:
            raise ValueError("geomag_baseline must be finite and positive")
        if not math.isfinite(float(params.circadian_period)) or params.circadian_period <= 0.0:
            raise ValueError("circadian_period must be finite and positive")
        if params.circadian_amplitude > 0.5:
            raise ValueError("circadian_amplitude must keep circadian output within [0, 1]")
        if params.rng_seed is not None and (
            not isinstance(params.rng_seed, int)
            or isinstance(params.rng_seed, bool)
            or params.rng_seed < 0
        ):
            raise ValueError("rng_seed must be None or a non-negative integer")

    @classmethod
    def _validate_step_inputs(
        cls,
        dt: float,
        l5_input: Optional[Dict[str, Any]],
        solar_activity: float,
        lunar_phase: float,
    ) -> None:
        if not math.isfinite(float(dt)) or dt <= 0.0:
            raise ValueError("dt must be finite and positive")
        if not math.isfinite(float(solar_activity)) or solar_activity < 0.0 or solar_activity > 1.0:
            raise ValueError("solar_activity must be finite and within [0, 1]")
        if not math.isfinite(float(lunar_phase)):
            raise ValueError("lunar_phase must be finite")
        if l5_input is not None:
            cls._l5_organismal_effect(l5_input)

    @staticmethod
    def _finite_mean(values: Any, name: str) -> float:
        arr = np.asarray(values, dtype=np.float64)
        if arr.size == 0 or not np.all(np.isfinite(arr)):
            raise ValueError(f"{name} must contain finite values")
        return float(np.mean(arr))

    @classmethod
    def _l5_organismal_effect(cls, l5_input: Dict[str, Any]) -> float:
        if "ecological_drive" in l5_input:
            return cls._unit_mean(l5_input["ecological_drive"], "ecological_drive")
        if "emotional_state" in l5_input:
            emotional_coherence = cls._finite_mean(l5_input["emotional_state"], "emotional_state")
            return emotional_coherence - 0.5
        return 0.0

    @classmethod
    def _unit_mean(cls, values: Any, name: str) -> float:
        arr = np.asarray(values, dtype=np.float64)
        if arr.size == 0 or not np.all(np.isfinite(arr)):
            raise ValueError(f"{name} must contain finite values")
        if np.any(arr < 0.0) or np.any(arr > 1.0):
            raise ValueError(f"{name} values must be within [0, 1]")
        return float(np.mean(arr))

step(dt, l5_input=None, solar_activity=0.5, lunar_phase=0.0)

Advance the layer by one time step.

Parameters:

Name Type Description Default
dt float

Time step in seconds.

 required
l5_input Optional[Dict[str, Any]]

Organismal layer output (emotional coherence).

 None
solar_activity float

Solar activity index (0-1).

 0.5
lunar_phase float

Lunar phase (0 to 2π).

 0.0

Returns:

Type Description
Dict[str, Any]

Dict with schumann_field, geomag, circadian, output_bitstreams
Source code in src/sc_neurocore/scpn/layers/l6_ecological.py
Python
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
def step(
    self,
    dt: float,
    l5_input: Optional[Dict[str, Any]] = None,
    solar_activity: float = 0.5,
    lunar_phase: float = 0.0,
) -> Dict[str, Any]:
    """
    Advance the layer by one time step.

    Args:
        dt: Time step in seconds.
        l5_input: Organismal layer output (emotional coherence).
        solar_activity: Solar activity index (0-1).
        lunar_phase: Lunar phase (0 to 2π).

    Returns:
        Dict with schumann_field, geomag, circadian, output_bitstreams
    """
    self._validate_step_inputs(dt, l5_input, solar_activity, lunar_phase)
    self.time += dt

    # 1. Schumann resonance dynamics
    for i, freq in enumerate(self.params.schumann_frequencies):
        self.schumann_phases[i] += 2 * np.pi * freq * dt
        self.schumann_phases[i] = self.schumann_phases[i] % (2 * np.pi)

    # Compute Schumann field as superposition
    schumann_signal = np.zeros(self.params.n_field_nodes)
    for i, freq in enumerate(self.params.schumann_frequencies):
        spatial_pattern = np.sin(np.linspace(0, 2 * np.pi * (i + 1), self.params.n_field_nodes))
        schumann_signal += (
            self.schumann_amplitudes[i]
            * self.params.schumann_amplitude
            * np.cos(self.schumann_phases[i])
            * spatial_pattern
        )

    # Add noise
    schumann_signal += self.params.schumann_noise * self._rng.normal(
        0, 1, self.params.n_field_nodes
    )

    # Normalize to [0, 1]
    schumann_field = (schumann_signal - schumann_signal.min()) / (
        schumann_signal.max() - schumann_signal.min() + 1e-8
    )

    # 2. Geomagnetic field dynamics
    # Solar activity modulates geomagnetic storms
    storm_factor = 1.0 + 0.5 * (solar_activity - 0.5)
    geomag_variation = (
        self.params.geomag_variation
        * storm_factor
        * self._rng.normal(0, 1, self.params.n_field_nodes)
    )
    self.geomag_field = np.clip(
        self.geomag_field + geomag_variation * dt,
        self.params.geomag_baseline * 0.5,
        self.params.geomag_baseline * 1.5,
    )

    # 3. Circadian rhythm
    self.circadian_phase += 2 * np.pi * dt / self.params.circadian_period
    self.circadian_phase = self.circadian_phase % (2 * np.pi)
    circadian_signal = 0.5 + self.params.circadian_amplitude * np.cos(self.circadian_phase)

    # 4. Biospheric network dynamics
    # Coupling between nodes
    network_coupling = np.zeros(self.params.n_field_nodes)
    for i in range(self.params.n_field_nodes):
        neighbors = [(i - 1) % self.params.n_field_nodes, (i + 1) % self.params.n_field_nodes]
        neighbor_mean = np.mean([self.biospheric_field[j] for j in neighbors])
        network_coupling[i] = neighbor_mean - self.biospheric_field[i]

    self.biospheric_field += (
        self.params.network_coupling * network_coupling
        + self.params.network_noise * self._rng.normal(0, 1, self.params.n_field_nodes)
    ) * dt

    # Modulate by Schumann and circadian
    self.biospheric_field *= (0.9 + 0.1 * schumann_field) * (0.8 + 0.2 * circadian_signal)
    self.biospheric_field = np.clip(self.biospheric_field, 0.0, 1.0)

    # 5. Organismal coupling (L5 collective emotional state affects field)
    if l5_input is not None:
        l5_effect = self._l5_organismal_effect(l5_input)
        if l5_effect != 0.0:
            self.biospheric_field += self.params.organismal_coupling * l5_effect * dt
            self.biospheric_field = np.clip(self.biospheric_field, 0.0, 1.0)

    # 6. Lunar phase modulation
    lunar_factor = 0.5 + 0.5 * np.cos(lunar_phase)
    self.schumann_amplitudes = np.ones(len(self.params.schumann_frequencies)) * (
        0.8 + 0.2 * lunar_factor
    )

    # 7. Compute planetary coherence
    self.planetary_coherence = float(
        np.abs(np.mean(np.exp(1j * 2 * np.pi * self.biospheric_field)))
    )

    # 8. Generate output bitstreams
    output_probs = self.biospheric_field * circadian_signal
    rands = self._rng.random((self.params.n_field_nodes, self.params.bitstream_length))
    output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)
    symbolic_drive = self.params.symbolic_coupling * schumann_field

    # Store history
    result = {
        "schumann_field": schumann_field,
        "schumann_phases": self.schumann_phases.copy(),
        "geomag_field": self.geomag_field.copy(),
        "circadian_phase": self.circadian_phase,
        "circadian_signal": circadian_signal,
        "biospheric_field": self.biospheric_field.copy(),
        "planetary_coherence": self.planetary_coherence,
        "symbolic_drive": symbolic_drive.copy(),
        "output_bitstreams": output_bitstreams,
    }

    self.history.append(
        {
            "time": self.time,
            "coherence": self.planetary_coherence,
            "schumann_power": float(np.mean(schumann_field**2)),
        }
    )
    if len(self.history) > 100:
        self.history.pop(0)

    return result

get_global_metric()

Return the global planetary coherence metric.

Source code in src/sc_neurocore/scpn/layers/l6_ecological.py
Python
235
236
237
def get_global_metric(self) -> float:
    """Return the global planetary coherence metric."""
    return self.planetary_coherence

get_schumann_spectrum()

Return current Schumann resonance spectrum.

Source code in src/sc_neurocore/scpn/layers/l6_ecological.py
Python
239
240
241
242
243
244
245
246
def get_schumann_spectrum(self) -> Dict[float, float]:
    """Return current Schumann resonance spectrum."""
    return {
        freq: float(amp * np.cos(phase))
        for freq, amp, phase in zip(
            self.params.schumann_frequencies, self.schumann_amplitudes, self.schumann_phases
        )
    }

get_circadian_time()

Return current circadian time (0-24 hours).

Source code in src/sc_neurocore/scpn/layers/l6_ecological.py
Python
248
249
250
def get_circadian_time(self) -> float:
    """Return current circadian time (0-24 hours)."""
    return (self.circadian_phase / (2 * np.pi)) * 24.0

L7_SymbolicLayer

Stochastic implementation of the Geometric-Symbolic Layer.

Models sacred geometry patterns, symbolic resonances, and acupuncture point dynamics using bitstream representations.

Source code in src/sc_neurocore/scpn/layers/l7_symbolic.py
Python
 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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
class L7_SymbolicLayer:
    """
    Stochastic implementation of the Geometric-Symbolic Layer.

    Models sacred geometry patterns, symbolic resonances, and
    acupuncture point dynamics using bitstream representations.
    """

    # Platonic solid vertex counts
    PLATONIC_VERTICES = {
        "tetrahedron": 4,
        "cube": 8,
        "octahedron": 6,
        "dodecahedron": 20,
        "icosahedron": 12,
    }

    # Fibonacci sequence for alignment
    FIBONACCI = [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144]

    def __init__(self, params: Optional[L7_StochasticParameters] = None):
        self.params = params or L7_StochasticParameters()
        self._validate_params(self.params)
        self._rng = np.random.default_rng(self.params.rng_seed)

        # Symbol activation states
        self.symbol_activations = self._rng.random(self.params.n_symbols) * 0.3

        # Sacred geometry metrics
        self.phi_alignment = 0.5
        self.fibonacci_alignment = 0.5
        self.metatron_flow = 0.5
        self.platonic_coherence = 0.5
        self.e8_alignment = 0.5
        self.symbolic_health = 0.5

        # Meridian states (TCM)
        self.meridian_qi = np.ones(self.params.n_meridians) * 0.5

        # Acupuncture point activations
        self.acupoint_activations = np.zeros(self.params.n_acupoints)

        # Glyph vector (normalized output)
        self.glyph_vector = np.zeros(self.params.glyph_dimensions)

        # E8 lattice representation in the canonical 8D root-space basis.
        self.e8_state = self._rng.random(8) * 0.5

        # Time
        self.time = 0.0

    def step(
        self,
        dt: float,
        l6_input: Optional[Dict[str, Any]] = None,
        symbol_input: Optional[np.ndarray[Any, Any]] = None,
        acupoint_stimulus: Optional[Dict[int, float]] = None,
    ) -> Dict[str, Any]:
        """
        Advance the layer by one time step.

        Args:
            dt: Time step in seconds.
            l6_input: Ecological layer output (Schumann, circadian).
            symbol_input: External symbolic input vector.
            acupoint_stimulus: Dict of {point_id: intensity} for acupuncture.

        Returns:
            Dict with glyph_vector, meridian_qi, sacred_geometry, output_bitstreams
        """
        if not math.isfinite(float(dt)) or float(dt) <= 0.0:
            raise ValueError("dt must be finite and positive")
        self.time += dt

        # 1. Process symbol input
        if symbol_input is not None:
            symbol_values = self._symbol_input(symbol_input)
            self.symbol_activations = np.clip(
                self.symbol_activations + symbol_values * 0.2, 0.0, 1.0
            )

        # 2. Compute Phi (Golden Ratio) alignment
        # Check how close symbol ratios are to Phi
        sorted_activations = np.sort(self.symbol_activations)[::-1]
        if sorted_activations[1] > 0.01:
            ratios = sorted_activations[:-1] / (sorted_activations[1:] + 1e-8)
            phi_distances = np.abs(ratios - PHI)
            self.phi_alignment = float(np.exp(-np.mean(phi_distances)))
        else:
            self.phi_alignment = 0.5

        # 3. Compute Fibonacci alignment
        # Check if activation levels follow Fibonacci ratios
        fib_normalized = np.array(self.FIBONACCI[:8]) / self.FIBONACCI[7]
        top_8 = sorted_activations[:8]
        if np.max(top_8) > 0.01:
            top_8_norm = top_8 / (np.max(top_8) + 1e-8)
            fib_corr = np.corrcoef(top_8_norm, fib_normalized)[0, 1]
            self.fibonacci_alignment = float(max(0, (fib_corr + 1) / 2))
        else:
            self.fibonacci_alignment = 0.5

        # 4. Compute Metatron's Cube flow
        # Based on 13-sphere / 78-line connectivity pattern
        metatron_nodes = min(13, self.params.n_symbols)
        active_nodes = np.sum(self.symbol_activations[:metatron_nodes] > 0.5)
        self.metatron_flow = active_nodes / metatron_nodes
        # Add flow dynamics
        self.metatron_flow = 0.9 * self.metatron_flow + 0.1 * self._rng.random()

        # 5. Compute Platonic solid coherence
        platonic_metrics = []
        for solid, vertices in self.PLATONIC_VERTICES.items():
            solid_activations = self.symbol_activations[:vertices]
            coherence = np.std(solid_activations)  # Lower std = more coherent
            platonic_metrics.append(1.0 - coherence)
        self.platonic_coherence = float(np.mean(platonic_metrics))

        # 6. E8 lattice alignment against the full 240-root system.
        e8_norm = np.linalg.norm(self.e8_state)
        if e8_norm > 0:
            e8_unit = self.e8_state / e8_norm
            roots = self._e8_roots()
            root_units = roots / np.linalg.norm(roots, axis=1, keepdims=True)
            alignments = np.abs(root_units @ e8_unit)
            self.e8_alignment = float(np.max(alignments))
        else:
            self.e8_alignment = 0.5

        # Update E8 state with noise
        self.e8_state += 0.1 * self._rng.normal(0.0, 1.0, 8) * dt
        self.e8_state = np.clip(self.e8_state, -1, 1)

        # 7. Compute symbolic health
        self.symbolic_health = (
            self.params.phi_alignment_weight * self.phi_alignment
            + self.params.fibonacci_weight * self.fibonacci_alignment
            + self.params.metatron_weight * self.metatron_flow
            + self.params.platonic_weight * self.platonic_coherence
            + self.params.e8_weight * self.e8_alignment
        )

        # 8. Meridian Qi dynamics
        # Qi flows through meridians with circadian modulation
        qi_flow = np.roll(self.meridian_qi, 1) - self.meridian_qi
        self.meridian_qi += qi_flow * self.params.symbol_coupling * dt

        # Ecological coupling (structured L6 symbolic drive, with Schumann fallback)
        if l6_input is not None:
            l6_effect = self._l6_symbolic_effect(l6_input)
            if l6_effect != 0.0:
                self.meridian_qi *= 1.0 + self.params.ecological_coupling * l6_effect

        self.meridian_qi = np.clip(self.meridian_qi, 0.0, 1.0)

        # 9. Acupuncture point activation
        if acupoint_stimulus is not None:
            for point_id, intensity in self._acupoint_stimulus(acupoint_stimulus).items():
                self.acupoint_activations[point_id] = np.clip(
                    self.acupoint_activations[point_id] + intensity, 0.0, 1.0
                )

        # Decay acupoint activations
        self.acupoint_activations *= 1.0 - self.params.symbol_decay * dt

        # 10. Assemble glyph vector
        self.glyph_vector[...] = np.array(
            [
                self.phi_alignment,
                self.fibonacci_alignment,
                self.metatron_flow,
                self.platonic_coherence,
                self.e8_alignment,
                self.symbolic_health,
            ]
        )

        # 11. Symbol dynamics (decay and coupling)
        # Coupling: nearby symbols influence each other
        coupling = np.roll(self.symbol_activations, 1) + np.roll(self.symbol_activations, -1)
        self.symbol_activations += (
            self.params.symbol_coupling * (coupling / 2 - self.symbol_activations) * dt
        )
        # Decay
        self.symbol_activations *= 1.0 - self.params.symbol_decay * dt
        self.symbol_activations = np.clip(self.symbol_activations, 0.0, 1.0)

        # 12. Generate output bitstreams
        output_probs = np.concatenate(
            [self.symbol_activations, self.meridian_qi, self.glyph_vector]
        )
        output_probs = output_probs[: self.params.n_symbols]

        rands = self._rng.random((self.params.n_symbols, self.params.bitstream_length))
        output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)

        return {
            "glyph_vector": self.glyph_vector.copy(),
            "phi_alignment": self.phi_alignment,
            "fibonacci_alignment": self.fibonacci_alignment,
            "metatron_flow": self.metatron_flow,
            "platonic_coherence": self.platonic_coherence,
            "e8_alignment": self.e8_alignment,
            "symbolic_health": self.symbolic_health,
            "cosmic_phase_drive": self.params.cosmic_coupling * self.symbolic_health,
            "meridian_qi": self.meridian_qi.copy(),
            "acupoint_activations": self.acupoint_activations.copy(),
            "e8_state": self.e8_state.copy(),
            "output_bitstreams": output_bitstreams,
        }

    def get_global_metric(self) -> float:
        """Return the global symbolic coherence metric."""
        return self.symbolic_health

    def get_glyph_vector_normalized(self) -> np.ndarray[Any, Any]:
        """Return normalized glyph vector for external use."""
        return self.glyph_vector / (np.max(self.glyph_vector) + 1e-8)

    def stimulate_meridian(self, meridian_id: int, intensity: float) -> None:
        """Stimulate a specific meridian."""
        if not 0 <= meridian_id < self.params.n_meridians:
            raise ValueError("meridian_id must be in range")
        if not math.isfinite(float(intensity)):
            raise ValueError("intensity must be finite")
        self.meridian_qi[meridian_id] = np.clip(self.meridian_qi[meridian_id] + intensity, 0.0, 1.0)

    def get_acupoint_map(self) -> Dict[str, float]:
        """Return clinically common named acupoint activations."""
        named_points = {
            "LI4_Hegu": 4,
            "ST36_Zusanli": 36,
            "SP6_Sanyinjiao": 60,
            "PC6_Neiguan": 96,
            "LV3_Taichong": 120,
            "GV20_Baihui": 200,
            "CV4_Guanyuan": 250,
            "BL23_Shenshu": 300,
        }
        return {
            name: float(self.acupoint_activations[idx])
            for name, idx in named_points.items()
            if idx < self.params.n_acupoints
        }

    @staticmethod
    def _e8_roots() -> np.ndarray:
        roots: list[np.ndarray] = []
        for i in range(8):
            for j in range(i + 1, 8):
                for si in (-1.0, 1.0):
                    for sj in (-1.0, 1.0):
                        root = np.zeros(8, dtype=np.float64)
                        root[i] = si
                        root[j] = sj
                        roots.append(root)
        for mask in range(256):
            signs = np.array(
                [-1.0 if (mask >> bit) & 1 else 1.0 for bit in range(8)],
                dtype=np.float64,
            )
            if np.count_nonzero(signs < 0.0) % 2 == 0:
                roots.append(0.5 * signs)
        return np.asarray(roots, dtype=np.float64)

    @staticmethod
    def _validate_params(params: L7_StochasticParameters) -> None:
        if not isinstance(params.n_symbols, int) or isinstance(params.n_symbols, bool):
            raise ValueError("n_symbols must be a positive integer")
        if params.n_symbols < 2:
            raise ValueError("n_symbols must be at least two")
        if not isinstance(params.n_meridians, int) or isinstance(params.n_meridians, bool):
            raise ValueError("n_meridians must be a positive integer")
        if params.n_meridians <= 0:
            raise ValueError("n_meridians must be positive")
        if not isinstance(params.n_acupoints, int) or isinstance(params.n_acupoints, bool):
            raise ValueError("n_acupoints must be a positive integer")
        if params.n_acupoints <= 0:
            raise ValueError("n_acupoints must be positive")
        if not isinstance(params.bitstream_length, int) or isinstance(
            params.bitstream_length, bool
        ):
            raise ValueError("bitstream_length must be a positive integer")
        if params.bitstream_length <= 0:
            raise ValueError("bitstream_length must be positive")
        if params.glyph_dimensions != 6:
            raise ValueError("glyph_dimensions must be six for the current glyph contract")
        weights = np.array(
            [
                params.phi_alignment_weight,
                params.fibonacci_weight,
                params.metatron_weight,
                params.platonic_weight,
                params.e8_weight,
            ],
            dtype=np.float64,
        )
        if (
            not np.all(np.isfinite(weights))
            or np.any(weights < 0.0)
            or float(np.sum(weights)) <= 0.0
        ):
            raise ValueError("weights must be finite, non-negative, and sum positive")
        for field_name in (
            "symbol_decay",
            "symbol_coupling",
            "ecological_coupling",
            "cosmic_coupling",
        ):
            value = float(getattr(params, field_name))
            if not math.isfinite(value) or value < 0.0:
                raise ValueError(f"{field_name} must be finite and non-negative")
        if params.rng_seed is not None:
            if isinstance(params.rng_seed, bool) or not isinstance(params.rng_seed, int):
                raise ValueError("rng_seed must be a non-negative integer or None")
            if params.rng_seed < 0:
                raise ValueError("rng_seed must be a non-negative integer or None")

    def _symbol_input(self, symbol_input: np.ndarray[Any, Any]) -> np.ndarray:
        values = np.asarray(symbol_input, dtype=np.float64).reshape(-1)
        if values.size < self.params.n_symbols:
            raise ValueError("symbol_input must contain at least n_symbols values")
        trimmed = values[: self.params.n_symbols]
        if not np.all(np.isfinite(trimmed)):
            raise ValueError("symbol_input must contain only finite values")
        return trimmed

    @staticmethod
    def _finite_mean(value: Any, name: str) -> float:
        values = np.asarray(value, dtype=np.float64).reshape(-1)
        if values.size == 0:
            raise ValueError(f"{name} must contain at least one value")
        if not np.all(np.isfinite(values)):
            raise ValueError(f"{name} must contain only finite values")
        return float(np.mean(values))

    @classmethod
    def _l6_symbolic_effect(cls, l6_input: Dict[str, Any]) -> float:
        if "symbolic_drive" in l6_input:
            return cls._unit_mean(l6_input["symbolic_drive"], "symbolic_drive")
        if "schumann_field" in l6_input:
            return cls._finite_mean(l6_input["schumann_field"], "schumann_field") - 1.0
        return 0.0

    @staticmethod
    def _unit_mean(value: Any, name: str) -> float:
        values = np.asarray(value, dtype=np.float64).reshape(-1)
        if values.size == 0:
            raise ValueError(f"{name} must contain at least one value")
        if not np.all(np.isfinite(values)):
            raise ValueError(f"{name} must contain only finite values")
        if np.any(values < 0.0) or np.any(values > 1.0):
            raise ValueError(f"{name} values must be within [0, 1]")
        return float(np.mean(values))

    def _acupoint_stimulus(self, stimulus: Dict[int, float]) -> Dict[int, float]:
        validated: Dict[int, float] = {}
        for point_id, intensity in stimulus.items():
            if not isinstance(point_id, int) or isinstance(point_id, bool):
                raise ValueError("acupoint_stimulus keys must be integer point ids")
            if not 0 <= point_id < self.params.n_acupoints:
                raise ValueError("acupoint_stimulus point id out of range")
            intensity_value = float(intensity)
            if not math.isfinite(intensity_value):
                raise ValueError("acupoint_stimulus intensities must be finite")
            validated[point_id] = intensity_value
        return validated

step(dt, l6_input=None, symbol_input=None, acupoint_stimulus=None)

Advance the layer by one time step.

Parameters:

Name Type Description Default
dt float

Time step in seconds.

 required
l6_input Optional[Dict[str, Any]]

Ecological layer output (Schumann, circadian).

 None
symbol_input Optional[ndarray[Any, Any]]

External symbolic input vector.

 None
acupoint_stimulus Optional[Dict[int, float]]

Dict of {point_id: intensity} for acupuncture.

 None

Returns:

Type Description
Dict[str, Any]

Dict with glyph_vector, meridian_qi, sacred_geometry, output_bitstreams
Source code in src/sc_neurocore/scpn/layers/l7_symbolic.py
Python
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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
def step(
    self,
    dt: float,
    l6_input: Optional[Dict[str, Any]] = None,
    symbol_input: Optional[np.ndarray[Any, Any]] = None,
    acupoint_stimulus: Optional[Dict[int, float]] = None,
) -> Dict[str, Any]:
    """
    Advance the layer by one time step.

    Args:
        dt: Time step in seconds.
        l6_input: Ecological layer output (Schumann, circadian).
        symbol_input: External symbolic input vector.
        acupoint_stimulus: Dict of {point_id: intensity} for acupuncture.

    Returns:
        Dict with glyph_vector, meridian_qi, sacred_geometry, output_bitstreams
    """
    if not math.isfinite(float(dt)) or float(dt) <= 0.0:
        raise ValueError("dt must be finite and positive")
    self.time += dt

    # 1. Process symbol input
    if symbol_input is not None:
        symbol_values = self._symbol_input(symbol_input)
        self.symbol_activations = np.clip(
            self.symbol_activations + symbol_values * 0.2, 0.0, 1.0
        )

    # 2. Compute Phi (Golden Ratio) alignment
    # Check how close symbol ratios are to Phi
    sorted_activations = np.sort(self.symbol_activations)[::-1]
    if sorted_activations[1] > 0.01:
        ratios = sorted_activations[:-1] / (sorted_activations[1:] + 1e-8)
        phi_distances = np.abs(ratios - PHI)
        self.phi_alignment = float(np.exp(-np.mean(phi_distances)))
    else:
        self.phi_alignment = 0.5

    # 3. Compute Fibonacci alignment
    # Check if activation levels follow Fibonacci ratios
    fib_normalized = np.array(self.FIBONACCI[:8]) / self.FIBONACCI[7]
    top_8 = sorted_activations[:8]
    if np.max(top_8) > 0.01:
        top_8_norm = top_8 / (np.max(top_8) + 1e-8)
        fib_corr = np.corrcoef(top_8_norm, fib_normalized)[0, 1]
        self.fibonacci_alignment = float(max(0, (fib_corr + 1) / 2))
    else:
        self.fibonacci_alignment = 0.5

    # 4. Compute Metatron's Cube flow
    # Based on 13-sphere / 78-line connectivity pattern
    metatron_nodes = min(13, self.params.n_symbols)
    active_nodes = np.sum(self.symbol_activations[:metatron_nodes] > 0.5)
    self.metatron_flow = active_nodes / metatron_nodes
    # Add flow dynamics
    self.metatron_flow = 0.9 * self.metatron_flow + 0.1 * self._rng.random()

    # 5. Compute Platonic solid coherence
    platonic_metrics = []
    for solid, vertices in self.PLATONIC_VERTICES.items():
        solid_activations = self.symbol_activations[:vertices]
        coherence = np.std(solid_activations)  # Lower std = more coherent
        platonic_metrics.append(1.0 - coherence)
    self.platonic_coherence = float(np.mean(platonic_metrics))

    # 6. E8 lattice alignment against the full 240-root system.
    e8_norm = np.linalg.norm(self.e8_state)
    if e8_norm > 0:
        e8_unit = self.e8_state / e8_norm
        roots = self._e8_roots()
        root_units = roots / np.linalg.norm(roots, axis=1, keepdims=True)
        alignments = np.abs(root_units @ e8_unit)
        self.e8_alignment = float(np.max(alignments))
    else:
        self.e8_alignment = 0.5

    # Update E8 state with noise
    self.e8_state += 0.1 * self._rng.normal(0.0, 1.0, 8) * dt
    self.e8_state = np.clip(self.e8_state, -1, 1)

    # 7. Compute symbolic health
    self.symbolic_health = (
        self.params.phi_alignment_weight * self.phi_alignment
        + self.params.fibonacci_weight * self.fibonacci_alignment
        + self.params.metatron_weight * self.metatron_flow
        + self.params.platonic_weight * self.platonic_coherence
        + self.params.e8_weight * self.e8_alignment
    )

    # 8. Meridian Qi dynamics
    # Qi flows through meridians with circadian modulation
    qi_flow = np.roll(self.meridian_qi, 1) - self.meridian_qi
    self.meridian_qi += qi_flow * self.params.symbol_coupling * dt

    # Ecological coupling (structured L6 symbolic drive, with Schumann fallback)
    if l6_input is not None:
        l6_effect = self._l6_symbolic_effect(l6_input)
        if l6_effect != 0.0:
            self.meridian_qi *= 1.0 + self.params.ecological_coupling * l6_effect

    self.meridian_qi = np.clip(self.meridian_qi, 0.0, 1.0)

    # 9. Acupuncture point activation
    if acupoint_stimulus is not None:
        for point_id, intensity in self._acupoint_stimulus(acupoint_stimulus).items():
            self.acupoint_activations[point_id] = np.clip(
                self.acupoint_activations[point_id] + intensity, 0.0, 1.0
            )

    # Decay acupoint activations
    self.acupoint_activations *= 1.0 - self.params.symbol_decay * dt

    # 10. Assemble glyph vector
    self.glyph_vector[...] = np.array(
        [
            self.phi_alignment,
            self.fibonacci_alignment,
            self.metatron_flow,
            self.platonic_coherence,
            self.e8_alignment,
            self.symbolic_health,
        ]
    )

    # 11. Symbol dynamics (decay and coupling)
    # Coupling: nearby symbols influence each other
    coupling = np.roll(self.symbol_activations, 1) + np.roll(self.symbol_activations, -1)
    self.symbol_activations += (
        self.params.symbol_coupling * (coupling / 2 - self.symbol_activations) * dt
    )
    # Decay
    self.symbol_activations *= 1.0 - self.params.symbol_decay * dt
    self.symbol_activations = np.clip(self.symbol_activations, 0.0, 1.0)

    # 12. Generate output bitstreams
    output_probs = np.concatenate(
        [self.symbol_activations, self.meridian_qi, self.glyph_vector]
    )
    output_probs = output_probs[: self.params.n_symbols]

    rands = self._rng.random((self.params.n_symbols, self.params.bitstream_length))
    output_bitstreams = (rands < output_probs[:, None]).astype(np.uint8)

    return {
        "glyph_vector": self.glyph_vector.copy(),
        "phi_alignment": self.phi_alignment,
        "fibonacci_alignment": self.fibonacci_alignment,
        "metatron_flow": self.metatron_flow,
        "platonic_coherence": self.platonic_coherence,
        "e8_alignment": self.e8_alignment,
        "symbolic_health": self.symbolic_health,
        "cosmic_phase_drive": self.params.cosmic_coupling * self.symbolic_health,
        "meridian_qi": self.meridian_qi.copy(),
        "acupoint_activations": self.acupoint_activations.copy(),
        "e8_state": self.e8_state.copy(),
        "output_bitstreams": output_bitstreams,
    }

get_global_metric()

Return the global symbolic coherence metric.

Source code in src/sc_neurocore/scpn/layers/l7_symbolic.py
Python
281
282
283
def get_global_metric(self) -> float:
    """Return the global symbolic coherence metric."""
    return self.symbolic_health

get_glyph_vector_normalized()

Return normalized glyph vector for external use.

Source code in src/sc_neurocore/scpn/layers/l7_symbolic.py
Python
285
286
287
def get_glyph_vector_normalized(self) -> np.ndarray[Any, Any]:
    """Return normalized glyph vector for external use."""
    return self.glyph_vector / (np.max(self.glyph_vector) + 1e-8)

stimulate_meridian(meridian_id, intensity)

Stimulate a specific meridian.

Source code in src/sc_neurocore/scpn/layers/l7_symbolic.py
Python
289
290
291
292
293
294
295
def stimulate_meridian(self, meridian_id: int, intensity: float) -> None:
    """Stimulate a specific meridian."""
    if not 0 <= meridian_id < self.params.n_meridians:
        raise ValueError("meridian_id must be in range")
    if not math.isfinite(float(intensity)):
        raise ValueError("intensity must be finite")
    self.meridian_qi[meridian_id] = np.clip(self.meridian_qi[meridian_id] + intensity, 0.0, 1.0)

get_acupoint_map()

Return clinically common named acupoint activations.

Source code in src/sc_neurocore/scpn/layers/l7_symbolic.py
Python
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
def get_acupoint_map(self) -> Dict[str, float]:
    """Return clinically common named acupoint activations."""
    named_points = {
        "LI4_Hegu": 4,
        "ST36_Zusanli": 36,
        "SP6_Sanyinjiao": 60,
        "PC6_Neiguan": 96,
        "LV3_Taichong": 120,
        "GV20_Baihui": 200,
        "CV4_Guanyuan": 250,
        "BL23_Shenshu": 300,
    }
    return {
        name: float(self.acupoint_activations[idx])
        for name, idx in named_points.items()
        if idx < self.params.n_acupoints
    }

SCPNDatastream dataclass

In-memory representation of one deterministic SCPN stream.

Source code in src/sc_neurocore/scpn/datastream.py
Python
 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
@dataclass(frozen=True)
class SCPNDatastream:
    """In-memory representation of one deterministic SCPN stream."""

    dt_s: float
    seed: int
    probabilities: np.ndarray
    spike_train: np.ndarray
    omega_rad_s: np.ndarray
    knm: np.ndarray

    @property
    def n_steps(self) -> int:
        """Number of timesteps in the stream."""
        return int(self.spike_train.shape[0])

    @property
    def n_layers(self) -> int:
        """Number of SCPN layer channels in the stream."""
        return int(self.spike_train.shape[1])

    @property
    def firing_rates(self) -> np.ndarray:
        """Mean spike probability per layer over this stream window."""
        return self.spike_train.astype(np.float64).mean(axis=0)

    @property
    def rotation_angles_rad(self) -> np.ndarray:
        """``Ry`` angles for quantum-control bridges: firing rate times pi."""
        return self.firing_rates * np.pi

    @property
    def quantum_amplitudes(self) -> np.ndarray:
        """Real amplitude encoding of firing rates as ``[alpha, beta]`` pairs."""
        amplitudes = TensorStream.from_prob(self.firing_rates).to_quantum()
        return np.real(amplitudes).astype(np.float64)

    def to_json_dict(self) -> dict[str, Any]:
        """Serialise the datastream to a stable JSON-compatible mapping."""
        validate_scpn_datastream(self)
        return {
            "schema_version": SCHEMA_VERSION,
            "source_project": "sc-neurocore",
            "dt_s": self.dt_s,
            "seed": self.seed,
            "n_steps": self.n_steps,
            "n_layers": self.n_layers,
            "layer_ids": list(LAYER_IDS),
            "omega_rad_s": self.omega_rad_s.tolist(),
            "knm": self.knm.tolist(),
            "probabilities": self.probabilities.tolist(),
            "spike_train": self.spike_train.astype(int).tolist(),
            "firing_rates": self.firing_rates.tolist(),
            "rotation_angles_rad": self.rotation_angles_rad.tolist(),
            "quantum_amplitudes": self.quantum_amplitudes.tolist(),
        }

    @classmethod
    def from_json_dict(cls, payload: dict[str, Any]) -> SCPNDatastream:
        """Load and validate a datastream from a JSON-compatible mapping."""
        if not isinstance(payload, dict):
            raise ValueError("SCPN datastream JSON root must be an object")
        if payload.get("schema_version") != SCHEMA_VERSION:
            raise ValueError("unsupported SCPN datastream schema version")
        missing = [key for key in _REQUIRED_PAYLOAD_FIELDS if key not in payload]
        if missing:
            raise ValueError(f"SCPN datastream payload missing required fields: {missing}")

        spike_train = _numeric_array_from_payload(payload["spike_train"], key="spike_train")
        if not np.all((spike_train == 0.0) | (spike_train == 1.0)):
            raise ValueError("spike_train must be binary")

        stream = cls(
            dt_s=_float_scalar_from_payload(payload["dt_s"], key="dt_s"),
            seed=_int_scalar_from_payload(payload["seed"], key="seed"),
            probabilities=_numeric_array_from_payload(
                payload["probabilities"], key="probabilities"
            ),
            spike_train=spike_train.astype(np.uint8, copy=False),
            omega_rad_s=_numeric_array_from_payload(payload["omega_rad_s"], key="omega_rad_s"),
            knm=_numeric_array_from_payload(payload["knm"], key="knm"),
        )
        validate_scpn_datastream(stream)
        return stream

n_steps property

Number of timesteps in the stream.

n_layers property

Number of SCPN layer channels in the stream.

firing_rates property

Mean spike probability per layer over this stream window.

rotation_angles_rad property

Ry angles for quantum-control bridges: firing rate times pi.

quantum_amplitudes property

Real amplitude encoding of firing rates as [alpha, beta] pairs.

to_json_dict()

Serialise the datastream to a stable JSON-compatible mapping.

Source code in src/sc_neurocore/scpn/datastream.py
Python
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
def to_json_dict(self) -> dict[str, Any]:
    """Serialise the datastream to a stable JSON-compatible mapping."""
    validate_scpn_datastream(self)
    return {
        "schema_version": SCHEMA_VERSION,
        "source_project": "sc-neurocore",
        "dt_s": self.dt_s,
        "seed": self.seed,
        "n_steps": self.n_steps,
        "n_layers": self.n_layers,
        "layer_ids": list(LAYER_IDS),
        "omega_rad_s": self.omega_rad_s.tolist(),
        "knm": self.knm.tolist(),
        "probabilities": self.probabilities.tolist(),
        "spike_train": self.spike_train.astype(int).tolist(),
        "firing_rates": self.firing_rates.tolist(),
        "rotation_angles_rad": self.rotation_angles_rad.tolist(),
        "quantum_amplitudes": self.quantum_amplitudes.tolist(),
    }

from_json_dict(payload) classmethod

Load and validate a datastream from a JSON-compatible mapping.

Source code in src/sc_neurocore/scpn/datastream.py
Python
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
@classmethod
def from_json_dict(cls, payload: dict[str, Any]) -> SCPNDatastream:
    """Load and validate a datastream from a JSON-compatible mapping."""
    if not isinstance(payload, dict):
        raise ValueError("SCPN datastream JSON root must be an object")
    if payload.get("schema_version") != SCHEMA_VERSION:
        raise ValueError("unsupported SCPN datastream schema version")
    missing = [key for key in _REQUIRED_PAYLOAD_FIELDS if key not in payload]
    if missing:
        raise ValueError(f"SCPN datastream payload missing required fields: {missing}")

    spike_train = _numeric_array_from_payload(payload["spike_train"], key="spike_train")
    if not np.all((spike_train == 0.0) | (spike_train == 1.0)):
        raise ValueError("spike_train must be binary")

    stream = cls(
        dt_s=_float_scalar_from_payload(payload["dt_s"], key="dt_s"),
        seed=_int_scalar_from_payload(payload["seed"], key="seed"),
        probabilities=_numeric_array_from_payload(
            payload["probabilities"], key="probabilities"
        ),
        spike_train=spike_train.astype(np.uint8, copy=False),
        omega_rad_s=_numeric_array_from_payload(payload["omega_rad_s"], key="omega_rad_s"),
        knm=_numeric_array_from_payload(payload["knm"], key="knm"),
    )
    validate_scpn_datastream(stream)
    return stream

create_full_stack(params=None)

Create a complete 16-layer SCPN stack.

Source code in src/sc_neurocore/scpn/layers/__init__.py
Python
 99
100
101
102
def create_full_stack(params: Optional[dict[str, Any]] = None) -> dict[str, Any]:
    """Create a complete 16-layer SCPN stack."""
    params = params or {}
    return {key: cls(params.get(key)) for key, cls in LAYER_REGISTRY.items()}

run_integrated_step(layers, dt, inputs=None)

Run one integrated time step across all SCPN layers with inter-layer coupling.

Source code in src/sc_neurocore/scpn/layers/__init__.py
Python
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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
def run_integrated_step(
    layers: dict[str, Any], dt: float, inputs: Optional[dict[str, Any]] = None
) -> dict[str, Any]:
    """
    Run one integrated time step across all SCPN layers with inter-layer coupling.
    """
    inputs = inputs or {}
    outputs = {}

    # L1: Quantum (foundation)
    l1_bitstreams = layers["l1"].step(dt, external_field=inputs.get("l1_field"))
    outputs["l1"] = {
        "output_bitstreams": l1_bitstreams,
        "coherence": layers["l1"].get_global_metric(),
    }

    # L2: Neurochemical (receives L1 quantum modulation)
    l2_out = layers["l2"].step(dt, nt_release=inputs.get("nt_release"), l1_input=l1_bitstreams)
    outputs["l2"] = l2_out

    # L3: Genomic (receives L2 second messengers)
    l3_out = layers["l3"].step(dt, l2_input=l2_out, bioelectric_signal=inputs.get("bioelectric"))
    outputs["l3"] = l3_out

    # L4: Cellular (receives L3 protein modulation)
    l4_out = layers["l4"].step(
        dt, l3_input=l3_out, external_stimulus=inputs.get("cellular_stimulus")
    )
    outputs["l4"] = l4_out

    # L5: Organismal (receives L4 synchronization)
    l5_out = layers["l5"].step(dt, l4_input=l4_out, external_event=inputs.get("emotional_event"))
    outputs["l5"] = l5_out

    # L6: Ecological (receives L5 organismal state)
    l6_out = layers["l6"].step(
        dt,
        l5_input=l5_out,
        solar_activity=inputs.get("solar", 0.5),
        lunar_phase=inputs.get("lunar", 0.0),
    )
    outputs["l6"] = l6_out

    # L7: Symbolic (receives L6 Schumann/ecological)
    l7_out = layers["l7"].step(
        dt,
        l6_input=l6_out,
        symbol_input=inputs.get("symbols"),
        acupoint_stimulus=inputs.get("acupoints"),
    )
    outputs["l7"] = l7_out

    # L8: Phase Field / Cosmic (receives L7 symbolic)
    l8_out = layers["l8"].step(dt, l7_input=l7_out)
    outputs["l8"] = l8_out

    # L9: Memory (receives L8 cosmic alignment)
    l9_out = layers["l9"].step(dt, l8_input=l8_out)
    outputs["l9"] = l9_out

    # L10: Boundary (receives L9 retrieval quality)
    l10_out = layers["l10"].step(dt, l9_input=l9_out)
    outputs["l10"] = l10_out

    # L11: Morphic (receives L10 integrity)
    l11_out = layers["l11"].step(dt, l10_input=l10_out)
    outputs["l11"] = l11_out

    # L12: Quantum Info (receives L11 info saturation)
    l12_out = layers["l12"].step(dt, l11_input=l11_out)
    outputs["l12"] = l12_out

    # L13: Temporal (receives L12 coherence)
    l13_out = layers["l13"].step(dt, l12_input=l12_out)
    outputs["l13"] = l13_out

    # L14: Integration (receives metrics from all layers)
    all_metrics = get_global_metrics(layers)
    l14_out = layers["l14"].step(dt, layer_metrics=all_metrics, l13_input=l13_out)
    outputs["l14"] = l14_out

    # L15: Meta-cognitive (receives L14 integrated coherence)
    l15_out = layers["l15"].step(dt, l14_input=l14_out)
    outputs["l15"] = l15_out

    # L16: Director (receives L15 GCI)
    l16_out = layers["l16"].step(dt, l15_input=l15_out)
    outputs["l16"] = l16_out

    return outputs

get_global_metrics(layers)

Get global coherence metrics from all layers.

Source code in src/sc_neurocore/scpn/layers/__init__.py
Python
197
198
199
def get_global_metrics(layers: dict[str, Any]) -> dict[str, Any]:
    """Get global coherence metrics from all layers."""
    return {f"l{i}": layers[f"l{i}"].get_global_metric() for i in range(1, 17) if f"l{i}" in layers}

build_knm_matrix(n_layers=N_LAYERS)

Build the Knm inter-layer coupling matrix.

Construction: exponential decay baseline, calibration anchor overrides, cross-hierarchy boosts, symmetrisation, zero diagonal.

Source code in src/sc_neurocore/scpn/params.py
Python
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
def build_knm_matrix(n_layers: int = N_LAYERS) -> np.ndarray:
    """
    Build the Knm inter-layer coupling matrix.

    Construction: exponential decay baseline, calibration anchor overrides,
    cross-hierarchy boosts, symmetrisation, zero diagonal.
    """
    if not isinstance(n_layers, int) or isinstance(n_layers, bool) or n_layers <= 0:
        raise ValueError("n_layers must be a positive integer")

    K = np.zeros((n_layers, n_layers), dtype=np.float64)

    for n in range(n_layers):
        for m in range(n_layers):
            if n != m:
                K[n, m] = K_BASE * np.exp(-DECAY_ALPHA * abs(n - m))

    for (i, j), val in CALIBRATION_ANCHORS.items():
        if i <= n_layers and j <= n_layers:
            K[i - 1, j - 1] = val
            K[j - 1, i - 1] = val

    for (i, j), val in CROSS_BOOSTS.items():
        if i <= n_layers and j <= n_layers:
            K[i - 1, j - 1] = val
            K[j - 1, i - 1] = val

    K[...] = 0.5 * (K + K.T)
    np.fill_diagonal(K, 0.0)

    return K

generate_scpn_datastream(*, n_steps=32, dt_s=0.01, seed=1729, spike_floor=0.02, spike_ceiling=0.98)

Generate a deterministic 16-layer stream for inter-repository tests.

The probability envelope is a bounded phase oscillator driven by the canonical SCPN natural frequencies and a small normalised coupling bias derived from K_nm. The binary spike train is sampled from that envelope with a local RNG seeded by seed.

Source code in src/sc_neurocore/scpn/datastream.py
Python
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
173
174
175
176
177
178
def generate_scpn_datastream(
    *,
    n_steps: int = 32,
    dt_s: float = 0.01,
    seed: int = 1729,
    spike_floor: float = 0.02,
    spike_ceiling: float = 0.98,
) -> SCPNDatastream:
    """Generate a deterministic 16-layer stream for inter-repository tests.

    The probability envelope is a bounded phase oscillator driven by the
    canonical SCPN natural frequencies and a small normalised coupling
    bias derived from ``K_nm``. The binary spike train is sampled from
    that envelope with a local RNG seeded by ``seed``.
    """
    if n_steps < 1:
        raise ValueError("n_steps must be >= 1")
    if dt_s <= 0.0:
        raise ValueError("dt_s must be > 0")
    if not 0.0 <= spike_floor < spike_ceiling <= 1.0:
        raise ValueError("spike_floor and spike_ceiling must satisfy 0 <= floor < ceiling <= 1")

    omega = np.asarray(OMEGA_N, dtype=np.float64)
    knm = build_knm_matrix()
    coupling = knm.sum(axis=1)
    coupling_span = float(np.ptp(coupling))
    if coupling_span > 0.0:
        coupling_bias = (coupling - coupling.min()) / coupling_span
    else:
        coupling_bias = np.zeros_like(coupling)
    coupling_bias = coupling_bias - float(coupling_bias.mean())

    t = np.arange(n_steps, dtype=np.float64)[:, None] * dt_s
    phase = t * omega[None, :]
    probabilities = 0.5 + 0.25 * np.sin(phase) + 0.08 * coupling_bias[None, :]
    probabilities = np.clip(probabilities, spike_floor, spike_ceiling)

    rng = np.random.default_rng(seed)
    spike_train = (rng.random(probabilities.shape) < probabilities).astype(np.uint8)

    stream = SCPNDatastream(
        dt_s=dt_s,
        seed=seed,
        probabilities=probabilities,
        spike_train=spike_train,
        omega_rad_s=omega.copy(),
        knm=knm,
    )
    validate_scpn_datastream(stream)
    return stream

generate_scpn_datastream_payload(**kwargs)

Generate a JSON-compatible stream payload in one call.

Source code in src/sc_neurocore/scpn/datastream.py
Python
262
263
264
def generate_scpn_datastream_payload(**kwargs: Any) -> dict[str, Any]:
    """Generate a JSON-compatible stream payload in one call."""
    return generate_scpn_datastream(**kwargs).to_json_dict()

read_scpn_datastream(path)

Read a stream payload from JSON.

Source code in src/sc_neurocore/scpn/datastream.py
Python
217
218
219
220
221
222
def read_scpn_datastream(path: str | Path) -> SCPNDatastream:
    """Read a stream payload from JSON."""
    payload = json.loads(Path(path).read_text(encoding="utf-8"))
    if not isinstance(payload, dict):
        raise ValueError("SCPN datastream JSON root must be an object")
    return SCPNDatastream.from_json_dict(payload)

validate_scpn_datastream(stream)

Validate shape, bounds, and canonical matrix invariants.

Source code in src/sc_neurocore/scpn/datastream.py
Python
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
def validate_scpn_datastream(stream: SCPNDatastream) -> None:
    """Validate shape, bounds, and canonical matrix invariants."""
    if not np.isfinite(stream.dt_s) or stream.dt_s <= 0.0:
        raise ValueError("dt_s must be > 0")
    if stream.probabilities.shape != stream.spike_train.shape:
        raise ValueError("probabilities and spike_train must have matching shapes")
    if stream.probabilities.ndim != 2:
        raise ValueError("probabilities must be a 2-D array")
    if stream.spike_train.shape[1] != N_LAYERS:
        raise ValueError(f"spike_train must have {N_LAYERS} layer columns")
    if stream.omega_rad_s.shape != (N_LAYERS,):
        raise ValueError(f"omega_rad_s must have shape ({N_LAYERS},)")
    if stream.knm.shape != (N_LAYERS, N_LAYERS):
        raise ValueError(f"knm must have shape ({N_LAYERS}, {N_LAYERS})")
    if not np.isfinite(stream.probabilities).all():
        raise ValueError("probabilities must be finite")
    if not np.isfinite(stream.omega_rad_s).all():
        raise ValueError("omega_rad_s must be finite")
    if not np.isfinite(stream.knm).all():
        raise ValueError("knm must be finite")
    if not np.all((stream.probabilities >= 0.0) & (stream.probabilities <= 1.0)):
        raise ValueError("probabilities must be in [0, 1]")
    if not np.isin(stream.spike_train, [0, 1]).all():
        raise ValueError("spike_train must be binary")
    if not np.allclose(stream.knm, stream.knm.T, atol=1e-12):
        raise ValueError("knm must be symmetric")
    if not np.allclose(np.diag(stream.knm), 0.0, atol=1e-12):
        raise ValueError("knm diagonal must be zero")

write_scpn_datastream(path, stream)

Write a stream payload to JSON.

Source code in src/sc_neurocore/scpn/datastream.py
Python
211
212
213
214
def write_scpn_datastream(path: str | Path, stream: SCPNDatastream) -> None:
    """Write a stream payload to JSON."""
    path = Path(path)
    path.write_text(json.dumps(stream.to_json_dict(), indent=2, sort_keys=True) + "\n")