"""Hopf oscillator for amplitude dynamics and state transitions (Section 22)."""
from typing import List, Optional, Tuple
import numpy as np
from mercurial.params.empirical import HopfParams
[docs]
class HopfOscillator:
"""
Single Hopf oscillator: dz/dt = (α + iω)z - β|z|²z + noise.
"""
def __init__(
self, params: Optional[HopfParams] = None, frequency_band: str = "alpha", **kwargs
):
if params is None:
params = HopfParams()
self.beta = kwargs.get("beta", params.beta)
self.sigma = kwargs.get("sigma", params.sigma)
# Set omega based on frequency band
if frequency_band == "alpha":
omega = params.omega_alpha
elif frequency_band == "beta":
omega = params.omega_beta
elif frequency_band == "gamma":
omega = params.omega_gamma
else:
omega = float(frequency_band) * 2 * np.pi
self.omega = kwargs.get("omega", omega)
self.alpha = kwargs.get("alpha", params.alpha_rest)
self.z = 0.05 + 0.05j
self.rng = np.random.default_rng()
[docs]
def derivative(self, z: complex) -> complex:
return (self.alpha + 1j * self.omega) * z - self.beta * (z * np.conj(z)) * z
[docs]
def step(self, dt: float, external_input: complex = 0.0) -> None:
k1 = self.derivative(self.z) + external_input
z_pred = self.z + k1 * dt
k2 = self.derivative(z_pred) + external_input
dz = (k1 + k2) / 2.0
noise = self.sigma * np.sqrt(dt) * (self.rng.normal() + 1j * self.rng.normal())
self.z += dz * dt + noise
[docs]
def get_amplitude(self) -> float:
return np.abs(self.z)
[docs]
def get_phase(self) -> float:
return np.angle(self.z)
[docs]
def set_bifurcation(self, alpha: float) -> None:
self.alpha = alpha
[docs]
class HopfNetwork:
"""
Network of Hopf oscillators with Kuramoto phase coupling and diffusive amplitude coupling.
"""
def __init__(
self,
n_oscillators: int,
alphas: Optional[List[float]] = None,
omegas: Optional[List[float]] = None,
beta: float = 1.0,
sigma: float = 0.05,
kuramoto_coupling: float = 0.0,
diffusive_coupling: float = 0.0,
):
self.n = n_oscillators
self.beta = beta
self.sigma = sigma
self.K = kuramoto_coupling
self.eps = diffusive_coupling
# Default frequencies (gamma band)
if alphas is None:
alphas = [-0.05] * n_oscillators
if omegas is None:
omegas = [2 * np.pi * 40.0] * n_oscillators
self.oscillators = [
HopfOscillator(alpha=alphas[i], omega=omegas[i], beta=beta, sigma=sigma)
for i in range(n_oscillators)
]
[docs]
def step(self, dt: float) -> None:
phases = np.array([osc.get_phase() for osc in self.oscillators])
z_vals = np.array([osc.z for osc in self.oscillators])
# Kuramoto coupling
if self.K > 0:
mean_field = np.mean(np.exp(1j * phases))
for i, osc in enumerate(self.oscillators):
coupling = 1j * self.K * np.imag(mean_field * np.exp(-1j * phases[i])) * osc.z
osc.step(dt, external_input=coupling)
else:
for osc in self.oscillators:
osc.step(dt)
# Diffusive amplitude coupling
if self.eps > 0:
z_vals = np.array([osc.z for osc in self.oscillators])
for i, osc in enumerate(self.oscillators):
diff = np.sum(z_vals - z_vals[i]) / self.n
osc.step(dt, external_input=self.eps * diff)
[docs]
def get_order_parameter(self) -> Tuple[float, float]:
phases = np.array([osc.get_phase() for osc in self.oscillators])
z = np.mean(np.exp(1j * phases))
return np.abs(z), np.angle(z)
[docs]
def get_mean_amplitude(self) -> float:
return np.mean([osc.get_amplitude() for osc in self.oscillators])