OSCR

Neural sampling from cognitive maps enables goal-directed imagination and planning.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § Methods › Details of our model for spatial goal-directed imagination › Using the cognitive map for generating imagined paths ↔ gcml_grid_cell.ipynb, lines 1–61 · score 0.75 · constant speed, random walk, turn noise, forbidden, resampling, Invalid
  2. [2] § Methods › Details of our model for spatial goal-directed imagination › Generating a cognitive map based on grid cells ↔ gcml_grid_cell.ipynb, lines 142–285 · score 0.53 · event, plateau, arena, field, Pi, firing

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 874 lines · 26 KB · MIT · 2 matches

  1. # %%
  2. import numpy as np
  3. import matplotlib.pyplot as plt
  4. # Environment settings
  5. env_min = 0
  6. env_max = 4
  8. # Random walk settings
  9. n_steps = 50000
  10. dt = 0.05
  11. speed = 0.1 # constant speed
  12. turn_noise = 0.1 # radians per step
  13. n_traj = 1 # number of trajectories to plot
  14. # Helper: is a point valid?
  15. def is_valid(pos):
  16. x, y = pos
  17. if not (env_min <= x <= env_max and env_min <= y <= env_max):
  18. return False
  19. return True
  20. # Generate trajectories
  21. trajectories = []
  22. for _ in range(n_traj):
  23. pos = np.array([np.random.uniform(0, 4), np.random.uniform(0, 4)])
  24. angle = np.random.uniform(0, 2*np.pi)
  25. path = [pos.copy()]
  26. for _ in range(n_steps):
  27. # Add small random turn
  28. angle += np.random.randn() * turn_noise
  29. velocity = speed * np.array([np.cos(angle), np.sin(angle)])
  30. new_pos = pos + velocity * dt
  31. if is_valid(new_pos):
  32. pos = new_pos
  33. else:
  34. # Resample a new random direction if invalid move
  35. angle = np.random.uniform(0, 2*np.pi)
  36. path.append(pos.copy())
  37. trajectories.append(np.array(path))
  38. # Plot
  39. fig, ax = plt.subplots(figsize=(6, 6),dpi=300)
  40. for traj in trajectories:
  41. # Plot trajectory
  42. ax.plot(traj[:,0], traj[:,1], color='gray', alpha=0.5, linewidth=1)
  43. # for i in range(1, len(traj)):
  44. # ax.plot(traj[i-1:i+1,0], traj[i-1:i+1,1], color=‘gray’, alpha=0.5, linewidth=1)
  45. # Mark start
  46. start = traj[0]
  47. ax.plot(start[0], start[1], 'o', markersize=20, c='black')
  48. # Mark end
  49. end = traj[-1]
  50. ax.plot(end[0], end[1], 's', markersize=20, c='black')
  51. # Forbidden region
  52. # ax.add_patch(plt.Rectangle((0, 0), forbidden_x, forbidden_y,
  53. # color='red', alpha=0.2, zorder=0))
  54. # Formatting
  55. ax.set_xlim(env_min, env_max)
  56. ax.set_ylim(env_min, env_max)
  57. ax.set_aspect('equal')
  58. ax.set_xticks([])
  59. ax.set_yticks([])
  60. # ax.set_title(‘Random Walk with Forbidden Region’)
  61. plt.show()
  63. # %%
  65. import numpy as np
  66. import matplotlib.pyplot as plt
  68. # Environment settings
  69. env_min = 0
  70. env_max = 4
  71. forbidden_x = 2
  72. forbidden_y = 2
  74. # Random walk settings
  75. n_steps = 50000
  76. dt = 0.05
  77. speed = 0.1 # constant speed
  78. turn_noise = 0.1 # radians per step
  79. n_traj = 1 # number of trajectories to plot
  81. # Helper: is a point valid?
  82. def is_valid(pos):
  83. x, y = pos
  84. if not (env_min <= x <= env_max and env_min <= y <= env_max):
  85. return False
  86. if x < forbidden_x and y < forbidden_y:
  87. return False
  88. return True
  90. # Generate trajectories
  91. trajectories = []
  92. for _ in range(n_traj):
  93. pos = np.array([np.random.uniform(2, 4), np.random.uniform(2, 4)]) # start outside forbidden
  94. angle = np.random.uniform(0, 2*np.pi)
  95. path = [pos.copy()]
  97. for _ in range(n_steps):
  98. # Add small random turn
  99. angle += np.random.randn() * turn_noise
  100. velocity = speed * np.array([np.cos(angle), np.sin(angle)])
  101. new_pos = pos + velocity * dt
  103. if is_valid(new_pos):
  104. pos = new_pos
  105. else:
  106. # Resample a new random direction if invalid move
  107. angle = np.random.uniform(0, 2*np.pi)
  109. path.append(pos.copy())
  111. trajectories.append(np.array(path))
  113. # Plot
  114. fig, ax = plt.subplots(figsize=(6, 6),dpi=300)
  115. for traj in trajectories:
  116. # Plot trajectory
  117. ax.plot(traj[:,0], traj[:,1], color='gray', alpha=0.99, linewidth=1)
  119. # for i in range(1, len(traj)):
  120. # ax.plot(traj[i-1:i+1,0], traj[i-1:i+1,1], color='gray', alpha=0.5, linewidth=1)
  121. # Mark start
  122. start = traj[0]
  123. # ax.plot(start[0], start[1], 'o', markersize=20, c='black')
  124. # Mark end
  125. end = traj[-1]
  126. # ax.plot(end[0], end[1], 's', markersize=20, c='black')
  128. # Forbidden region
  129. # ax.add_patch(plt.Rectangle((0, 0), forbidden_x, forbidden_y,
  130. # color='red', alpha=0.2, zorder=0))
  132. # Formatting
  133. ax.set_xlim(env_min, env_max)
  134. ax.set_ylim(env_min, env_max)
  135. ax.set_aspect('equal')
  136. ax.set_xticks([])
  137. ax.set_yticks([])
  138. # ax.set_title('Random Walk with Forbidden Region')
  140. plt.show()
  142. # %%
  143. import numpy as np
  144. import matplotlib.pyplot as plt
  146. # -------------------------
  147. # Function Definitions
  148. # -------------------------
  150. def is_valid(pos, env_min=0, env_max=4):
  151. x, y = pos
  152. return env_min <= x <= env_max and env_min <= y <= env_max
  154. def generate_grid_cells(n_cells=1000, arena_size=4, seed=None):
  155. if seed is not None:
  156. np.random.seed(seed)
  157. scales = np.random.uniform(0.02, 8, size=n_cells)
  158. phases = np.random.uniform(0, arena_size, size=(n_cells, 2))
  159. orientations = np.random.uniform(0, np.pi/3, size=n_cells)
  160. return scales, phases, orientations
  162. def grid_activation(x, y, scales, phases, orientations):
  163. positions = np.array([x, y])
  164. activations = []
  165. angles = np.array([0, np.pi / 3, 2 * np.pi / 3])
  166. for scale, phase, orientation in zip(scales, phases, orientations):
  167. # Rotate and compute grating
  168. R = np.array([[np.cos(orientation), -np.sin(orientation)],
  169. [np.sin(orientation), np.cos(orientation)]])
  170. pos_rot = R @ (positions - phase)
  171. proj = pos_rot[0] * np.cos(angles) + pos_rot[1] * np.sin(angles)
  172. grating = np.sum(np.cos((4 * np.pi / (scale * np.sqrt(3))) * proj))
  173. activations.append((2/3) * grating)
  174. return np.array(activations)
  176. def grid_activation_vectorized(X, Y, scales, phases, orientations):
  177. H, W = X.shape
  178. n_cells = len(scales)
  179. pos = np.stack([X.ravel(), Y.ravel()], axis=1) # (H*W, 2)
  180. activations = np.zeros((H * W, n_cells))
  181. angles = np.array([0, np.pi / 3, 2 * np.pi / 3])
  183. for i in range(n_cells):
  184. scale = scales[i]
  185. phase = phases[i]
  186. orientation = orientations[i]
  187. c, s = np.cos(orientation), np.sin(orientation)
  188. R = np.array([[c, -s], [s, c]])
  189. rel_pos = pos - phase
  190. pos_rot = rel_pos @ R.T
  191. proj = np.outer(pos_rot[:, 0], np.cos(angles)) + np.outer(pos_rot[:, 1], np.sin(angles))
  192. grating = np.sum(np.cos((4 * np.pi / (scale * np.sqrt(3))) * proj), axis=1)
  193. activations[:, i] = (2/3) * grating
  195. return activations.reshape(H, W, n_cells)
  197. # -------------------------
  198. # Random Walk Generation
  199. # -------------------------
  201. env_min, env_max = 0, 4
  202. n_steps = 50000
  203. dt = 0.05
  204. speed = 0.2
  205. turn_noise = 0.1
  207. pos = np.random.uniform(env_min, env_max, size=2)
  208. angle = np.random.uniform(0, 2*np.pi)
  209. trajectory = [pos.copy()]
  211. for _ in range(n_steps):
  212. angle += np.random.randn() * turn_noise
  213. vel = speed * np.array([np.cos(angle), np.sin(angle)])
  214. new_pos = pos + vel * dt
  215. if is_valid(new_pos, env_min, env_max):
  216. pos = new_pos
  217. else:
  218. angle = np.random.uniform(0, 2*np.pi)
  219. trajectory.append(pos.copy())
  221. trajectory = np.array(trajectory)
  223. # -------------------------
  224. # Grid Cell Setup
  225. # -------------------------
  227. n_cells = 1000
  228. scales, phases, orientations = generate_grid_cells(n_cells, env_max)
  230. # -------------------------
  231. # Plateau Events & Weights
  232. # -------------------------
  234. n_plateau = 1000
  235. plateau_times = np.sort(np.random.choice(np.arange(len(trajectory)), size=n_plateau, replace=False))
  237. weights = np.zeros((n_plateau, n_cells))
  238. for i, t in enumerate(plateau_times):
  239. x, y = trajectory[t]
  240. weights[i] = grid_activation(x, y, scales, phases, orientations)
  242. # -------------------------
  243. # Plot 1: Random Walk + Plateau Fires
  244. # -------------------------
  246. fig, ax = plt.subplots(figsize=(6, 6), dpi=300)
  247. ax.plot(trajectory[:, 0], trajectory[:, 1], color='gray', alpha=0.5, linewidth=1)
  248. ax.scatter(trajectory[plateau_times, 0], trajectory[plateau_times, 1],
  249. c='black', s=10, label='Plateau Fires')
  250. ax.plot(trajectory[0, 0], trajectory[0, 1], 'o', markersize=12, c='black')
  251. ax.plot(trajectory[-1, 0], trajectory[-1, 1], 's', markersize=12, c='black')
  252. ax.set_xlim(env_min, env_max)
  253. ax.set_ylim(env_min, env_max)
  254. ax.set_aspect('equal')
  255. ax.axis('off')
  256. plt.show()
  258. # -------------------------
  259. # Compute Reactive Fields
  260. # -------------------------
  262. grid_res = 50
  263. x_vals = np.linspace(env_min, env_max, grid_res)
  264. y_vals = np.linspace(env_min, env_max, grid_res)
  265. X, Y = np.meshgrid(x_vals, y_vals)
  267. grid_embeddings = grid_activation_vectorized(X, Y, scales, phases, orientations)
  268. # resulting shape: (grid_res, grid_res, n_cells)
  270. # reactive_fields(h, w, neuron) = sum_i grid_embeddings(h,w,i) * weights(neuron,i)
  271. reactive_fields = np.tensordot(grid_embeddings, weights, axes=([2], [1])) # shape: (H, W, n_plateau)
  273. # -------------------------
  274. # Plot 2: First 100 Reactive Fields
  275. # -------------------------
  277. fig2, axes = plt.subplots(10, 10, figsize=(20, 20))
  278. for idx, ax in enumerate(axes.flatten()):
  279. im = ax.imshow(reactive_fields[:, :, idx],
  280. extent=[env_min, env_max, env_min, env_max],
  281. origin='lower')
  282. ax.axis('off')
  284. plt.tight_layout()
  285. plt.show()
  288. # %%
  289. from sklearn.decomposition import PCA
  290. import numpy as np
  291. import matplotlib.pyplot as plt
  293. # Identify top 1% largest-scale grid cells
  294. threshold = np.percentile(scales,85)
  295. large_indices = np.where(scales >= threshold)[0]
  297. # Extract embeddings for only those cells and flatten
  298. H, W, _ = grid_embeddings.shape
  299. emb_large = grid_embeddings[:, :, large_indices].reshape(H * W, len(large_indices))
  301. # PCA on filtered embeddings
  302. pca_large = PCA(n_components=2)
  303. pc_large = pca_large.fit_transform(emb_large)
  305. # Plot
  306. plt.figure(figsize=(6, 6), dpi=300)
  307. plt.scatter(pc_large[:, 0], pc_large[:, 1], s=5)
  308. # plt.title('PCA of Top 1%‑Scale Grid Cell Embeddings')
  309. # plt.xlabel('PC1')
  310. plt.axis('off')
  311. # plt.ylabel('PC2')
  312. plt.tight_layout()
  313. plt.show()
  316. # %%
  317. from sklearn.decomposition import PCA
  318. import numpy as np
  319. import matplotlib.pyplot as plt
  320. from matplotlib.colors import Normalize, ListedColormap
  322. # Identify top 15% largest-scale grid cells (you had 85th percentile)
  323. threshold = np.percentile(scales, 85)
  324. large_indices = np.where(scales >= threshold)[0]
  326. # Extract embeddings for only those cells and flatten
  327. H, W, _ = grid_embeddings.shape
  328. emb_large = grid_embeddings[:, :, large_indices].reshape(H * W, len(large_indices))
  330. # PCA on filtered embeddings
  331. pca_large = PCA(n_components=2)
  332. pc_large = pca_large.fit_transform(emb_large)
  334. # --- Colormap: Blues truncated to [0.1, 0.3] ---
  335. blues = plt.cm.Blues
  336. trunc_blues = ListedColormap(blues(np.linspace(0.7, 0.8, 256)))
  338. # Normalize by x-values (PC1)
  339. norm = Normalize(vmin=pc_large[:, 0].min(), vmax=pc_large[:, 0].max())
  341. # Plot
  342. plt.figure(figsize=(6, 6), dpi=300)
  343. plt.scatter(pc_large[:, 0], pc_large[:, 1],
  344. c=pc_large[:, 0], cmap=trunc_blues, norm=norm,
  345. s=40, linewidths=0)
  346. plt.axis('off')
  347. plt.tight_layout()
  348. plt.show()
  351. # %%
  352. import numpy as np
  353. import matplotlib.pyplot as plt
  354. from sklearn.decomposition import PCA
  356. grid_res = 50
  357. x_vals = np.linspace(env_min, env_max, grid_res)
  358. y_vals = np.linspace(env_min, env_max, grid_res)
  359. X, Y = np.meshgrid(x_vals, y_vals)
  361. # Grid embedding and coords from workspace
  362. H, W, n_cells = grid_embeddings.shape
  363. dx = X[0, 1] - X[0, 0]
  364. dy = Y[1, 0] - Y[0, 0]
  366. # Top 10% largest-scale cells
  367. threshold_10 = np.percentile(scales, 80)
  368. top10_idx = np.where(scales >= threshold_10)[0]
  369. n_top = len(top10_idx)
  371. # Extract embeddings for top cells
  372. emb_top10 = grid_embeddings.take(top10_idx, axis=2)
  374. # Compute finite differences
  375. partial_x = (emb_top10[:, 1:, :] - emb_top10[:, :-1, :]) / dx
  376. partial_y = (emb_top10[1:, :, :] - emb_top10[:-1, :, :]) / dy
  378. # Average over space
  379. avg_px = partial_x.reshape(-1, n_top).mean(axis=0)
  380. avg_py = partial_y.reshape(-1, n_top).mean(axis=0)
  382. # W matrix (#top,4): [Up, Down, Left, Right]
  383. W = np.stack([ avg_py, -avg_py, -avg_px, avg_px ], axis=1)
  385. # Build flattened embedding matrix manually (H*W, n_top)
  386. emb_flat_top = np.stack([emb_top10[:, :, i].ravel() for i in range(n_top)], axis=1)
  388. # PCA on embeddings of top cells
  389. pca_top10 = PCA(n_components=2)
  390. pc_top10 = pca_top10.fit_transform(emb_flat_top)
  392. # Project action vectors
  393. vectors = pca_top10.components_.dot(W) # (2,4)
  395. # Plot PCA scatter + arrows
  396. plt.figure(figsize=(6, 6),dpi=300)
  397. plt.scatter(pc_top10[:, 0], pc_top10[:, 1], s=0.2,c='black')
  399. actions = ['Up', 'Down', 'Left', 'Right']
  400. for i, action in enumerate(actions):
  401. vx, vy = vectors[0, i], vectors[1, i]
  402. plt.arrow(20, 0, vx, vy, head_width=0.4,fc='black')
  403. plt.text(vx * 1.3, vy * 1.3, action)
  405. # plt.title('PCA of Top-10%-Scale Grid Embeddings with Action Vectors')
  406. plt.xlabel('PC1')
  407. plt.ylabel('PC2')
  408. plt.axis('equal')
  409. plt.tight_layout()
  410. plt.axis('off')
  411. plt.show()
  413. print("W matrix shape:", W.shape)
  416. # %%
  417. import numpy as np
  418. import matplotlib.pyplot as plt
  419. import matplotlib.patches as patches
  421. # -------------------------
  422. # Helpers from earlier
  423. # -------------------------
  424. def grid_activation(x, y, scales, phases, orientations):
  425. pos = np.array([x, y])
  426. angles = np.array([0, np.pi/3, 2*np.pi/3])
  427. acts = []
  428. for scale, phase, orient in zip(scales, phases, orientations):
  429. R = np.array([[np.cos(orient), -np.sin(orient)],
  430. [np.sin(orient), np.cos(orient)]])
  431. pr = R @ (pos - phase)
  432. proj = pr[0]*np.cos(angles) + pr[1]*np.sin(angles)
  433. grating = np.sum(np.cos((4*np.pi/(scale*np.sqrt(3))) * proj))
  434. acts.append((2/3)*grating)
  435. return np.array(acts)
  437. def is_inside_rect(pos, rect):
  438. x, y, l, orient = rect
  439. if orient == 'h':
  440. x0, y0 = x - l/2, y - 0.05
  441. x1, y1 = x + l/2, y + 0.05
  442. else:
  443. x0, y0 = x - 0.05, y - l/2
  444. x1, y1 = x + 0.05, y + l/2
  445. return (x0 <= pos[0] <= x1) and (y0 <= pos[1] <= y1)
  447. def seg_intersect(A, B, C, D):
  448. def ccw(P, Q, R):
  449. return (R[1]-P[1])*(Q[0]-P[0]) > (Q[1]-P[1])*(R[0]-P[0])
  450. return ccw(A, C, D) != ccw(B, C, D) and ccw(A, B, C) != ccw(A, B, D)
  452. def segment_intersects_rect(p1, p2, rect):
  453. x0, x1, y0, y1 = (rect[0] - rect[2]/2, rect[0] + rect[2]/2,
  454. rect[1] - (0.01 if rect[3]=='h' else rect[2]/2),
  455. rect[1] + (0.01 if rect[3]=='h' else rect[2]/2))
  456. if x0 <= p2[0] <= x1 and y0 <= p2[1] <= y1:
  457. return True
  458. edges = [
  459. (np.array([x0, y0]), np.array([x1, y0])),
  460. (np.array([x1, y0]), np.array([x1, y1])),
  461. (np.array([x1, y1]), np.array([x0, y1])),
  462. (np.array([x0, y1]), np.array([x0, y0])),
  463. ]
  464. for C, D in edges:
  465. if seg_intersect(p1, p2, C, D):
  466. return True
  467. return False
  469. def is_valid_position(pos, walls):
  470. if not (0 <= pos[0] <= 4 and 0 <= pos[1] <= 4):
  471. return False
  472. for w in walls:
  473. if is_inside_rect(pos, w):
  474. return False
  475. return True
  477. def compute_valid_step(pos, vel, dt, walls):
  478. full = vel * dt
  479. cand = pos + full
  480. if is_valid_position(cand, walls) and not any(segment_intersects_rect(pos, cand, w) for w in walls):
  481. return cand
  482. # slide
  483. for dx, dy in [(full[0], 0), (0, full[1])]:
  484. cand2 = pos + np.array([dx, dy])
  485. if is_valid_position(cand2, walls) and not any(segment_intersects_rect(pos, cand2, w) for w in walls):
  486. return cand2
  487. # fallback smaller
  488. for f in [0.5, 0.25, 0.1]:
  489. cand3 = pos + full * f
  490. if is_valid_position(cand3, walls) and not any(segment_intersects_rect(pos, cand3, w) for w in walls):
  491. return cand3
  492. return pos
  494. # -------------------------
  495. # Setup grid code & W
  496. # -------------------------
  497. # Use pre-generated scales, phases, orientations, X, Y from workspace
  499. # Select top 10%
  500. threshold = np.percentile(scales, 90)
  501. top_idx = np.where(scales >= threshold)[0]
  502. sc_top = scales[top_idx]
  503. ph_top = phases[top_idx]
  504. or_top = orientations[top_idx]
  505. n_top = len(top_idx)
  507. # Compute W by finite differences averaged
  508. # reuse grid_embeddings, X, Y
  509. H, Wg, _ = grid_embeddings.shape
  510. dx = X[0,1] - X[0,0]
  511. dy = Y[1,0] - Y[0,0]
  512. emb_top = grid_embeddings[:,:,top_idx]
  513. pd_x = (emb_top[:,1:,:] - emb_top[:,:-1,:]) / dx
  514. pd_y = (emb_top[1:,:,:] - emb_top[:-1,:,:]) / dy
  515. avg_px = pd_x.reshape(-1, n_top).mean(axis=0)
  516. avg_py = pd_y.reshape(-1, n_top).mean(axis=0)
  517. Wmat = np.stack([avg_py, -avg_py, -avg_px, avg_px], axis=1)
  519. # -------------------------
  520. # Navigation via grid code
  521. # -------------------------
  522. walls = [] # no walls
  523. start = np.array([0.5, 2.5])
  524. goal = np.array([2.5, 0.5])
  525. noise_level = 80
  528. dt = 0.05
  529. threshold = 0.1
  530. max_steps = 2000
  532. # Precompute goal embedding
  533. goal_emb = grid_activation(goal[0], goal[1], sc_top, ph_top, or_top)
  535. paths = []
  536. for _ in range(5):
  537. pos = start.copy()
  538. path = [pos.copy()]
  539. for _ in range(max_steps):
  540. # current embedding
  541. cur_emb = grid_activation(pos[0], pos[1], sc_top, ph_top, or_top)
  542. delta = goal_emb - cur_emb # high-D target diff
  543. scores = delta @ Wmat + noise_level * np.random.randn(4) # (4,) for [Up,Down,Left,Right]
  544. # compute movement
  545. step = np.array([scores[3] - scores[2], scores[0] - scores[1]])
  546. # normalize to unit-length
  547. if np.linalg.norm(step) > 1e-6:
  548. step = step / np.linalg.norm(step)
  549. pos = compute_valid_step(pos, step, dt, walls)
  550. path.append(pos.copy())
  551. if np.linalg.norm(pos - goal) < threshold:
  552. break
  553. paths.append(np.array(path))
  555. # -------------------------
  556. # Plotting
  557. # -------------------------
  558. fig, ax = plt.subplots(figsize=(6,6))
  560. ax.plot(start[0], start[1], 'o', ms=20, c='black')
  561. ax.plot(goal[0], goal[1], '*', ms=30, c='black')
  562. colors = plt.cm.rainbow(np.linspace(0,1,len(paths)))
  563. for p, c in zip(paths, colors):
  564. ax.plot(p[:,0], p[:,1], color=c, lw=2, alpha=0.8)
  565. ax.set_xlim(0,4); ax.set_ylim(0,4); ax.set_aspect('equal')
  567. ax.set_xticks([])
  568. ax.set_yticks([])
  569. # ax.axis('off')
  570. plt.show()
  573. # %%
  574. probs.max()
  576. # %%
  577. import numpy as np
  578. import matplotlib.pyplot as plt
  580. # — your existing grid_activation, scales, phases, orientations, and BTSP weights —
  581. # grid_activation(x,y,scales,phases,orientations) → (n_grid_cells,)
  582. # weights: shape (n_place_cells, n_grid_cells)
  583. # paths: list of trajectories, each shape [T,2]
  585. # 1) Grab first trajectory
  586. traj = paths[0]
  587. T = traj.shape[0]
  588. n_pc = weights.shape[0]
  590. # 2) Sample spikes from learned place-cell activations
  591. spikes = np.zeros((T, n_pc), dtype=int)
  592. for t, (x, y) in enumerate(traj):
  593. gact = grid_activation(x, y, scales, phases, orientations) # pre-syn grid code
  594. p_act = weights @ gact # learned place-cell drive
  595. exps = np.exp(p_act - p_act.max())
  596. probs = exps / exps.sum() * 1 + 0.0005*np.random.randn(exps.shape[0]) # softmax → firing probs
  597. spikes[t] = (np.random.rand(n_pc) < probs).astype(int)
  599. # 3) Pick top-50 most active, rank by mean spike time
  600. counts = spikes.sum(axis=0)
  601. avg_time = (np.arange(T)[:,None] * spikes).sum(axis=0) / (counts + 1e-9)
  603. top50 = np.argsort(counts)[-50:] # highest-firing 50
  604. order = np.argsort(avg_time[top50]) # earliest mean time first
  605. selected = top50[order]
  607. # 4) Raster plot (spikes only)
  608. fig, ax = plt.subplots(figsize=(3,5), dpi=300)
  609. for r, neuron in enumerate(selected):
  610. times = np.where(spikes[:, neuron])[0]
  611. ax.vlines(times, r + 0.2, r + 0.8, color='black', linewidth=0.7)
  613. ax.set_xlim(0, T)
  614. ax.set_ylim(0, 50)
  615. # ax.set_xlabel('Time step')
  616. ticks = ax.get_xticks() # e.g. [0, 100, 200, …]
  617. ax.set_xticks(ticks)
  618. ax.set_xticklabels([(t/100) for t in ticks])
  620. ax.set_yticks([0,10,20,30,40,50]) # hide neuron labels
  621. plt.tight_layout()
  622. plt.show()
  625. # %%
  626. ticks
  628. # %%
  629. fig, ax = plt.subplots(figsize=(1.5,2.5), dpi=300)
  630. for r, neuron in enumerate(selected):
  631. times = np.where(spikes[:, neuron])[0]
  632. ax.vlines(times, r + 0.2, r + 0.8, color='black', linewidth=0.4)
  634. ax.set_xlim(0, T)
  635. ax.set_ylim(0, 50)
  636. # ax.set_xlabel('Time step')
  637. ticks = ax.get_xticks() # e.g. [0, 100, 200, …]
  638. # ticks = [0,100,200,300,400,500]
  639. ax.set_xticks(ticks)
  640. ax.set_xticklabels([(t/100) for t in ticks])
  642. ax.set_yticks([0,10,20,30,40,50]) # hide neuron labels
  643. plt.tight_layout()
  644. plt.show()
  647. # %%
  648. import numpy as np
  649. import matplotlib.pyplot as plt
  650. import matplotlib.patches as patches
  652. # -------------------------
  653. # Helpers from earlier
  654. # -------------------------
  655. def grid_activation(x, y, scales, phases, orientations):
  656. pos = np.array([x, y])
  657. angles = np.array([0, np.pi/3, 2*np.pi/3])
  658. acts = []
  659. for scale, phase, orient in zip(scales, phases, orientations):
  660. R = np.array([[np.cos(orient), -np.sin(orient)],
  661. [np.sin(orient), np.cos(orient)]])
  662. pr = R @ (pos - phase)
  663. proj = pr[0]*np.cos(angles) + pr[1]*np.sin(angles)
  664. grating = np.sum(np.cos((4*np.pi/(scale*np.sqrt(3))) * proj))
  665. acts.append((2/3)*grating)
  666. return np.array(acts)
  669. def rect_bounds(rect, eps=0.01):
  670. x, y, l, orient = rect
  671. if orient == 'h':
  672. # horizontal: thin in y
  673. return (x - l/2, x + l/2,
  674. y - eps, y + eps)
  675. else:
  676. # vertical: thin in x
  677. return (x - eps, x + eps,
  678. y - l/2, y + l/2)
  680. def is_inside_rect(pos, rect):
  681. x0, x1, y0, y1 = rect_bounds(rect)
  682. return (x0 <= pos[0] <= x1) and (y0 <= pos[1] <= y1)
  684. def segment_intersects_rect(p1, p2, rect):
  685. x0, x1, y0, y1 = rect_bounds(rect)
  686. # If endpoint lies within rect, it intersects
  687. if x0 <= p2[0] <= x1 and y0 <= p2[1] <= y1:
  688. return True
  689. # Otherwise check each edge
  690. edges = [
  691. (np.array([x0, y0]), np.array([x1, y0])),
  692. (np.array([x1, y0]), np.array([x1, y1])),
  693. (np.array([x1, y1]), np.array([x0, y1])),
  694. (np.array([x0, y1]), np.array([x0, y0]))
  695. ]
  696. for C, D in edges:
  697. if seg_intersect(p1, p2, C, D):
  698. return True
  699. return False
  701. # def is_inside_rect(pos, rect):
  702. # x, y, l, orient = rect
  703. # if orient == 'h':
  704. # x0, y0 = x - l/2, y - 0.05
  705. # x1, y1 = x + l/2, y + 0.05
  706. # else:
  707. # x0, y0 = x - 0.05, y - l/2
  708. # x1, y1 = x + 0.05, y + l/2
  709. # return (x0 <= pos[0] <= x1) and (y0 <= pos[1] <= y1)
  711. def ccw(A, B, C):
  712. return (C[1]-A[1])*(B[0]-A[0]) > (B[1]-A[1])*(C[0]-A[0])
  714. def seg_intersect(A, B, C, D):
  715. return ccw(A, C, D) != ccw(B, C, D) and ccw(A, B, C) != ccw(A, B, D)
  717. # def segment_intersects_rect(p1, p2, rect):
  718. # x0, x1, y0, y1 = (
  719. # rect[0] - rect[2]/2, rect[0] + rect[2]/2,
  720. # rect[1] - (0.01 if rect[3]=='h' else rect[2]/2),
  721. # rect[1] + (0.01 if rect[3]=='h' else rect[2]/2)
  722. # )
  723. # if x0 <= p2[0] <= x1 and y0 <= p2[1] <= y1:
  724. # return True
  725. # edges = [
  726. # (np.array([x0, y0]), np.array([x1, y0])),
  727. # (np.array([x1, y0]), np.array([x1, y1])),
  728. # (np.array([x1, y1]), np.array([x0, y1])),
  729. # (np.array([x0, y1]), np.array([x0, y0])),
  730. # ]
  731. # for C, D in edges:
  732. # if seg_intersect(p1, p2, C, D):
  733. # return True
  734. return False
  736. def is_valid_position(pos, walls):
  737. if not (0 <= pos[0] <= 4 and 0 <= pos[1] <= 4):
  738. return False
  739. for w in walls:
  740. if is_inside_rect(pos, w):
  741. return False
  742. return True
  744. def compute_valid_step(pos, vel, dt, walls):
  745. full = vel * dt
  746. cand = pos + full
  747. if is_valid_position(cand, walls) and not any(segment_intersects_rect(pos, cand, w) for w in walls):
  748. return cand
  749. # slide
  750. for dx, dy in [(full[0], 0), (0, full[1])]:
  751. cand2 = pos + np.array([dx, dy])
  752. if is_valid_position(cand2, walls) and not any(segment_intersects_rect(pos, cand2, w) for w in walls):
  753. return cand2
  754. # fallback smaller
  755. for f in [0.5, 0.25, 0.1]:
  756. cand3 = pos + full * f
  757. if is_valid_position(cand3, walls) and not any(segment_intersects_rect(pos, cand3, w) for w in walls):
  758. return cand3
  759. return pos
  761. # -------------------------
  762. # Setup grid code & W
  763. # -------------------------
  764. # Top 10% scales
  765. threshold = np.percentile(scales, 90)
  766. top_idx = np.where(scales >= threshold)[0]
  767. sc_top = scales[top_idx]
  768. ph_top = phases[top_idx]
  769. or_top = orientations[top_idx]
  770. n_top = len(top_idx)
  772. # Finite diff to build Wmat
  773. H, Wg, _ = grid_embeddings.shape
  774. dx = X[0,1] - X[0,0]
  775. dy = Y[1,0] - Y[0,0]
  776. emb_top = grid_embeddings[:,:,top_idx]
  777. pd_x = (emb_top[:,1:,:] - emb_top[:,:-1,:]) / dx
  778. pd_y = (emb_top[1:,:,:] - emb_top[:-1,:,:]) / dy
  779. avg_px = pd_x.reshape(-1, n_top).mean(axis=0)
  780. avg_py = pd_y.reshape(-1, n_top).mean(axis=0)
  781. Wmat = np.stack([avg_py, -avg_py, -avg_px, avg_px], axis=1)
  783. # -------------------------
  784. # Navigation with repulsion
  785. # -------------------------
  786. walls = [
  787. (1,1.5,1,'h'),
  788. (1.5,1,1,'v'),
  789. (1.5,2,1,'v'),
  790. (2,2.5,1,'h'),
  791. (3,2.5,1,'h'),
  792. (2.5,3,1,'v')
  793. ]
  794. # Precompute repulsion points
  795. rep_points = []
  796. for x, y, l, orient in walls:
  797. rep_points.append(np.array([x, y]))
  798. if orient == 'h':
  799. rep_points += [np.array([x-l/2,y]), np.array([x+l/2,y])]
  800. else:
  801. rep_points += [np.array([x,y-l/2]), np.array([x,y+l/2])]
  803. start = np.array([0.2, 0.2])
  804. goal = np.array([2.0, 2.0])
  805. k_wall = 0.08
  806. noise_level = 30.0
  807. dt = 0.05
  808. threshold_goal = 0.1
  809. max_steps = 1000
  810. n_paths = 5
  811. max_dist = np.linalg.norm(start - goal)
  812. goal_emb = grid_activation(goal[0], goal[1], sc_top, ph_top, or_top)
  814. def compute_velocity(pos):
  815. # Repulsive 2D
  816. f_rep = np.zeros(2)
  817. for wc in rep_points:
  818. delta = pos - wc
  819. dsq = np.sum(delta**2)
  820. if dsq > 1e-4:
  821. f_rep += k_wall * delta / dsq
  822. rep_scale = np.clip(np.linalg.norm(pos-goal) / max_dist, 0, 1)
  823. f_rep *= rep_scale
  825. # Grid-code guidance
  826. cur_emb = grid_activation(pos[0], pos[1], sc_top, ph_top, or_top)
  827. delta_emb = goal_emb - cur_emb
  828. scores = delta_emb @ Wmat + noise_level * np.random.randn(4) # [Up,Down,Left,Right]
  829. step = np.array([scores[3]-scores[2], scores[0]-scores[1]])
  830. if np.linalg.norm(step) > 1e-6:
  831. step = step / np.linalg.norm(step)
  833. return step + f_rep # combine
  835. # Generate paths
  836. paths = []
  837. for _ in range(n_paths):
  838. pos = start.copy()
  839. path = [pos.copy()]
  840. for _ in range(max_steps):
  841. vel = compute_velocity(pos)
  842. pos = compute_valid_step(pos, vel, dt, walls)
  843. path.append(pos.copy())
  844. if np.linalg.norm(pos - goal) < threshold_goal:
  845. break
  846. paths.append(np.array(path))
  848. # -------------------------
  849. # Plotting
  850. # -------------------------
  851. fig, ax = plt.subplots(figsize=(6,6))
  852. # draw walls
  853. for x, y, l, orient in walls:
  854. if orient == 'h':
  855. ax.add_patch(patches.Rectangle((x-l/2, y-0.02), l, 0.04, color='black'))
  856. else:
  857. ax.add_patch(patches.Rectangle((x-0.02, y-l/2), 0.04, l, color='black'))
  858. # start & goal
  859. ax.plot(start[0], start[1], 'o', ms=12, c='black')
  860. ax.plot(goal[0], goal[1], '*', ms=15, c='black')
  861. # paths
  862. colors = plt.cm.rainbow(np.linspace(0,1,n_paths))
  863. for p, c in zip(paths, colors):
  864. ax.plot(p[:,0], p[:,1], color=c, lw=2, alpha=0.8)
  865. ax.set_xlim(0,4); ax.set_ylim(0,4); ax.set_aspect('equal')
  866. ax.set_xticks([]); ax.set_yticks([])
  867. plt.show()
  870. # %%
  871. vel
  873. # %%

gcml_grid_cell.ipynb at commit ff76859, under MIT · at the source

Overview

Authors: Hui Lin1,2, Yukun Yang2, Rong Zhao1, Giovanni Pezzulo3, Wolfgang Maass2
ORCID iDs: Yukun Yang, Rong Zhao, Giovanni Pezzulo, Wolfgang Maass
  1. Department of Precision Instruments, Center for Brain-Inspired Computing Research (CBICR), Tsinghua University,Beijing, China
  2. Institute of Machine Learning and Neural Computation, Graz University of Technology,Graz, Austria
  3. Institute of Cognitive Sciences and Technologies, National Research Council,Rome, Italy
Institutions: Tsinghua University (China); Graz University of Technology (Austria); Institute of Cognitive Sciences and Technologies (Italy)
Journal: Nature machine intelligence, volume 8, issue 7, pages 1045-1065
Dates: received 6 August 2025; accepted 11 May 2026; published online 21 July 2026; in print 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s42256-026-01254-4 · PMID 42499994 · PMCID PMC13395624 · OpenAlex W7169839672
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: computational modeling (no new data) (modality), none (in silico) (organism), cognitive (subfield)
Methods: Connectivity, Smoothing, state filtering, decompositions, Graphs, Statistics, Machine learning, Single-unit activity, calcium imaging
Keywords: Computational science, Learning algorithms, Network models
Topic: Ferroelectric and Negative Capacitance Devices (Electrical and Electronic Engineering, Engineering), according to OpenAlex
Funding: Austrian Science Fund (10.55776/COE12); NSF (EFRI BRAID project 2318152); EC | EU Framework Programme for Research and Innovation H2020 | H2020 Priority Excellent Science | H2020 European Research Council (H2020 Excellent Science - European Research Council) (No. 820213 (ThinkAhead))
Citations: not cited yet (Europe PMC); 104 references in the paper

Abstract

Artificial intelligence systems are becoming more intelligent, but at a very high cost in terms of energy consumption and training requirements. By contrast, our brains only require 20 W of energy, they learn online and they can instantly adjust to changing contingencies. This begs the question what data structures, algorithms and learning methods enable brains to achieve that, and whether these can be ported into artificial devices. We are addressing this question for a core feature of intelligence: the capacity to plan and solve problems, including new problems that involve states that were never encountered before. Here we examine three tools that brains are likely to use for achieving that: cognitive maps, stochastic computing and compositional coding. We integrate these tools into a transparent neural network model, and demonstrate its power for flexible planning and problem-solving. Importantly, this approach is suitable for implementation by in-memory computing and other energy-efficient neuromorphic hardware. In particular, it only requires self-supervised local synaptic plasticity that is suited for on-chip learning. Hence, a core feature of brain intelligence—the capacity to generate solutions to problems that were never encountered before—does not require deep neural networks or large language models, and can be implemented in energy-efficient edge devices.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

LH-cbicr/GCML

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: ff76859b71a2bc2056b50f5e052475351c007f76, 1 April 2026
Languages: Jupyter (3), Python (1)
Size: 10 files, 4 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, environment (environment.yml), 3 notebooks
Not found: CITATION.cff, tests, continuous integration, documentation
Tools: Matplotlib (4 files), NumPy (4 files), scikit-learn (4 files), NetworkX (3 files), PyTorch (3 files), h5py (2 files), seaborn (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
6 files
At the source: github.com/LH-cbicr/GCML

Zenodo 19370442

License: MIT
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Size: 1 file
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: Matplotlib (4 files), NumPy (4 files), scikit-learn (4 files), NetworkX (3 files), PyTorch (3 files), h5py (2 files), seaborn (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)
  • 27 September 2026: the link answers (HTTP 200)
6 files
At the source:

Code availability

The code used for training and evaluating the GCML is publicly available via GitHub at https://github.com/LH-cbicr/GCML and via Zenodo at 10.5281/zenodo.19370442 (ref. 104).

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 8 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

All data used in the experiments were generated synthetically and are publicly available via GitHub at https://github.com/LH-cbicr/GCML and via Zenodo at 10.5281/zenodo.19370442 (ref. 104).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Publisher: n/a → Nature Portfolio

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 3 keywords, 3 funders, 87 references.

Cite

This paper

Lin, H., Yang, Y., Zhao, R., Pezzulo, G., & Maass, W. (2026). Neural sampling from cognitive maps enables goal-directed imagination and planning. Nature machine intelligence, 8(7), 1045-1065. https://doi.org/10.1038/s42256-026-01254-4

BibTeX

@article{lin2026neural,
author = {Lin, Hui and Yang, Yukun and Zhao, Rong and Pezzulo, Giovanni and Maass, Wolfgang},
title = {{Neural sampling from cognitive maps enables goal-directed imagination and planning}},
journal = {Nature machine intelligence},
year = {2026},
month = jul,
volume = {8},
number = {7},
pages = {1045--1065},
publisher = {Nature Portfolio},
issn = {2522-5839},
doi = {10.1038/s42256-026-01254-4},
url = {https://doi.org/10.1038/s42256-026-01254-4},
pmid = {42499994},
pmcid = {PMC13395624}
}

RIS

TY - JOUR
AU - Lin, Hui
AU - Yang, Yukun
AU - Zhao, Rong
AU - Pezzulo, Giovanni
AU - Maass, Wolfgang
TI - Neural sampling from cognitive maps enables goal-directed imagination and planning
T2 - Nature machine intelligence
J2 - Nat Mach Intell
PY - 2026
DA - 2026/07/21
VL - 8
IS - 7
SP - 1045
EP - 1065
SN - 2522-5839
PB - Nature Portfolio
DO - 10.1038/s42256-026-01254-4
UR - https://doi.org/10.1038/s42256-026-01254-4
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s42256-026-01254-4",
"type": "article-journal",
"title": "Neural sampling from cognitive maps enables goal-directed imagination and planning",
"container-title": "Nature machine intelligence",
"author": [
{
"family": "Lin",
"given": "Hui"
},
{
"family": "Yang",
"given": "Yukun"
},
{
"family": "Zhao",
"given": "Rong"
},
{
"family": "Pezzulo",
"given": "Giovanni"
},
{
"family": "Maass",
"given": "Wolfgang"
}
],
"container-title-short": "Nat Mach Intell",
"volume": "8",
"issue": "7",
"page": "1045-1065",
"DOI": "10.1038/s42256-026-01254-4",
"PMID": "42499994",
"PMCID": "PMC13395624",
"ISSN": "2522-5839",
"publisher": "Nature Portfolio",
"URL": "https://doi.org/10.1038/s42256-026-01254-4",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
21
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1371/journal.pbio.3003755 [code]
Action information is integrated into entorhinal representations of conceptual space and is reflected in eye movements.
Journal: PLoS biology
In common: seaborn, Matplotlib, NumPy, cognitive, 16 references
[2] doi:10.1038/s41467-026-76102-5 [code]
The inherent capacity of neurons to learn order relations and support abstract reasoning.
Journal: Nature communications
In common: PyTorch, scikit-learn, Matplotlib, 1 other tool, computational modeling (no new data), cognitive, 6 references, 2 authors
[3] doi:10.1016/j.cub.2026.05.068 [code]
An abstract relational map emerges in the human medial prefrontal cortex with consolidation.
Journal: Current biology : CB
In common: seaborn, Matplotlib, NumPy, 9 references
[4] doi:10.1002/hipo.70131 [code]
Decoding Medial Entorhinal Cortical Dynamics Produces Planning-Like Alternations in Hippocampal theta Sequences.
Journal: Hippocampus
In common: h5py, scikit-learn, Matplotlib, 1 other tool, 7 references
[5] doi:10.1038/s41467-026-74357-6 [code]
Hippocampo-neocortical interaction as compressive retrieval-augmented generation.
Journal: Nature communications
In common: NetworkX, PyTorch, seaborn, 3 other tools, computational modeling (no new data), 4 references
[6] doi:10.1371/journal.pcbi.1013487 [code]
Flexible navigation with neuromodulated cognitive maps.
Journal: PLoS computational biology
In common: Matplotlib, NumPy, none (in silico), 6 references
[7] doi:10.1038/s41467-026-74358-5 [code]
Brain-inspired spatial intelligence for embodied agents.
Journal: Nature communications
In common: NetworkX, h5py, PyTorch, 4 other tools, cognitive, 2 references
[8] doi:10.1126/sciadv.aeg6797 [code]
Dorsoventral gradient of theta sweeps in the medial entorhinal cortex.
Journal: Science advances
In common: NetworkX, h5py, seaborn, 3 other tools, none (in silico), 2 references
[9] doi:10.1038/s41593-026-02232-0 [code]
Entorhinal cortex represents task-relevant remote locations independently of CA1.
Journal: Nature neuroscience
In common: NetworkX, h5py, PyTorch, 4 other tools, 2 references
[10] doi:10.1038/s41593-026-02365-2 [code]
Hippocampal theta sweeps indicate goal direction during navigation.
Journal: Nature neuroscience
In common: Matplotlib, NumPy, computational modeling (no new data), 5 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.