"""Test 2D neural field dynamics with spatial interactions."""
import matplotlib.pyplot as plt
import numpy as np
from mercurial.simulation.engine import SimulationEngine
# Create engine
engine = SimulationEngine(dim=10)
# Define a simple external stimulus: Gaussian blob at centre
[docs]
def gaussian_stimulus(t, nx=64, ny=64, sigma=5.0):
x = np.linspace(-nx / 2, nx / 2, nx)
y = np.linspace(-ny / 2, ny / 2, ny)
X, Y = np.meshgrid(x, y)
# Stimulus amplitude modulated by time (pulse)
amp = 0.2 * np.exp(-((t - 0.5) ** 2) / 0.1)
return amp * np.exp(-(X**2 + Y**2) / (2 * sigma**2))
# Run simulation
t_span = (0.0, 2.0)
results = engine.run_neural_field(
t_span,
dt=0.005,
nx=64,
ny=64,
dx=0.1,
external_stimulus=lambda t: gaussian_stimulus(t, nx=64, ny=64),
save_history_every=20,
)
# Plot final field
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.imshow(results["final_E"], cmap="hot", origin="lower")
plt.colorbar(label="Excitatory activity")
plt.title("Final Excitatory Field")
plt.subplot(1, 2, 2)
plt.imshow(results["final_I"], cmap="cool", origin="lower")
plt.colorbar(label="Inhibitory activity")
plt.title("Final Inhibitory Field")
plt.tight_layout()
plt.show()
# Plot a slice through the centre over time
E_slice = [E[32, :] for E in results["E_history"]]
plt.figure()
plt.imshow(
np.array(E_slice).T,
aspect="auto",
cmap="hot",
extent=[0, results["t"][-1], 0, 64],
origin="lower",
)
plt.xlabel("Time (s)")
plt.ylabel("Spatial index (y)")
plt.title("Excitatory Activity Slice (x = centre)")
plt.colorbar()
plt.show()
print(f"Final order parameter R: {results['order_parameter']:.4f}")