From 30517c0243144c2bbd08deaf3178265b77e42e3b Mon Sep 17 00:00:00 2001
From: "google-labs-jules[bot]"
 <161369871+google-labs-jules[bot]@users.noreply.github.com>
Date: Wed, 8 Jul 2026 02:16:48 +0000
Subject: [PATCH] =?UTF-8?q?=E2=9A=A1=20Bolt:=20optimize=20topoplot=20grid?=
 =?UTF-8?q?=20interpolation?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

- Optimized griddata_v4 in topoplot.py by replacing nested loops with NumPy vectorization.
- Achieved ~2.6x-4.8x speedup in biharmonic spline interpolation.
- Verified numerical parity with existing tests.
- Updated Bolt journal in .jules/bolt.md.

Co-authored-by: suraj-ranganath <14310165+suraj-ranganath@users.noreply.github.com>
---
 .jules/bolt.md                                |  3 ++
 src/eegprep/functions/sigprocfunc/topoplot.py | 31 +++++++++++--------
 2 files changed, 21 insertions(+), 13 deletions(-)
 create mode 100644 .jules/bolt.md

diff --git a/.jules/bolt.md b/.jules/bolt.md
new file mode 100644
index 00000000..2338c1d1
--- /dev/null
+++ b/.jules/bolt.md
@@ -0,0 +1,3 @@
+## 2026-07-08 - [Optimization of topoplot grid interpolation]
+**Learning:** Replaced a double-nested loop in `griddata_v4` with vectorized NumPy broadcasting and matrix multiplication (@). This is particularly effective for 2D grid evaluations in interpolation functions where query points are numerous.
+**Action:** Always look for nested loops over query grids in signal processing and visualization functions and replace them with 3D broadcasting + @ to shift computation to BLAS.
diff --git a/src/eegprep/functions/sigprocfunc/topoplot.py b/src/eegprep/functions/sigprocfunc/topoplot.py
index a63aff4c..32a25b5f 100644
--- a/src/eegprep/functions/sigprocfunc/topoplot.py
+++ b/src/eegprep/functions/sigprocfunc/topoplot.py
@@ -45,19 +45,24 @@ def griddata_v4(x, y, v, xq, yq):
         # If still singular, use pseudoinverse as last resort
         weights = np.linalg.pinv(g_reg) @ v
 
-    # Initialize output array
-    m, n = xq.shape
-    vq = np.zeros_like(xq)
-
-    # Evaluate at requested points
-    xy = xy[:, None]  # Make it column vector for broadcasting
-    for i in range(m):
-        for j in range(n):
-            d = np.abs(xq[i, j] + 1j * yq[i, j] - xy.ravel())
-            with np.errstate(divide='ignore', invalid='ignore'):
-                g = (d**2) * (np.log(d) - 1)  # Green's function
-            g[d == 0] = 0  # Handle Green's function at zero
-            vq[i, j] = np.dot(g, weights)
+    # Evaluate at requested points (vectorized)
+    # Combine xq and yq into complex numbers
+    q = xq + 1j * yq
+
+    # Calculate distances from all query points to all electrode points
+    # q has shape (m, n), xy has shape (k,)
+    # Resulting d will have shape (m, n, k)
+    d = np.abs(q[:, :, np.newaxis] - xy[np.newaxis, np.newaxis, :])
+
+    with np.errstate(divide='ignore', invalid='ignore'):
+        g = (d**2) * (np.log(d) - 1)  # Green's function
+
+    g[d == 0] = 0  # Handle Green's function at zero
+
+    # Weights has shape (k,)
+    # vq[i, j] = sum(g[i, j, k] * weights[k])
+    # This is equivalent to matrix multiplication g @ weights
+    vq = g @ weights
 
     return vq