Claude does not understand Gray code search at all

Author: WrathfulSpatula

pyqrackising/maxcut_tfim_util.py

@@ -416,6 +416,55 @@ def compute_cut_sparse(sample, G_data, G_rows, G_cols, n_qubits):
     return cut
 
 
+@njit(cache=True)
+def compute_cut_diff_sparse(u, sample, G_data, G_rows, G_cols, n_qubits):
+    energy = 0.0
+    u_bit = sample[u]
+    for col in range(G_rows[u], G_rows[u + 1]):
+        v = G_cols[col]
+        val = G_data[col]
+        energy += -val if u_bit == sample[v] else val
+
+    return energy
+
+
+@njit(cache=True)
+def compute_energy_diff_sparse(u, sample, G_data, G_rows, G_cols, n_qubits):
+    energy = 0.0
+    u_bit = sample[u]
+    for col in range(G_rows[u], G_rows[u + 1]):
+        v = G_cols[col]
+        val = G_data[col]
+        energy += -val if u_bit == sample[v] else val
+
+    return energy
+
+
+@njit(cache=True)
+def compute_cut_diff_2_sparse(k, l, sample, G_data, G_rows, G_cols, n_qubits):
+    if l < k:
+        t = k
+        k = l
+        l = t
+    energy = 0.0
+    k_bit = sample[k]
+    l_bit = sample[l]
+    for col in range(G_rows[k], G_rows[k + 1]):
+        v = G_cols[col]
+        if v == l:
+            continue
+        val = G_data[col]
+        energy += -val if k_bit == sample[v] else val
+    for col in range(G_rows[l], G_rows[l + 1]):
+        v = G_cols[col]
+        if v == k:
+            continue
+        val = G_data[col]
+        energy += -val if l_bit == sample[v] else val
+
+    return energy
+
+
 @njit(cache=True)
 def compute_energy_streaming(sample, G_func, nodes, n_qubits):
     energy = 0.0

pyqrackising/spin_glass_solver_sparse.py

@@ -2,6 +2,8 @@
 from .maxcut_tfim_util import (
     compute_cut_sparse,
     compute_energy_sparse,
+    compute_cut_diff_sparse,
+    compute_cut_diff_2_sparse,
     get_cut,
     gray_code_next,
     gray_mutation,
@@ -35,31 +37,19 @@
 @njit(parallel=True, cache=True)
 def run_single_bit_flips(best_theta, is_spin_glass, G_data, G_rows, G_cols):
     n = len(best_theta)
-    deltas = np.empty(n, dtype=dtype)
-
+    energies = np.empty(n, dtype=dtype)
     for i in prange(n):
-        # Compute delta: walk only adjacency list of i (O(degree(i)))
-        old_bit = best_theta[i]
-        new_bit = not old_bit
-        delta = dtype(0)
-        for col in range(G_rows[i], G_rows[i + 1]):
-            v = G_cols[col]
-            val = G_data[col]
-            v_bit = best_theta[v]
-            old_cut = old_bit != v_bit
-            new_cut = new_bit != v_bit
-            if old_cut != new_cut:
-                delta += val if new_cut else -val
-        if is_spin_glass:
-            delta *= dtype(2)
-        deltas[i] = delta
-
-    best_index = np.argmax(deltas)
-    best_delta = deltas[best_index]
+        state = best_theta.copy()
+        state[i] = not state[i]
+        energies[i] = compute_cut_diff_sparse(i, state, G_data, G_rows, G_cols, n)
+    best_index = np.argmax(energies)
+    best_energy = energies[best_index]
+    if is_spin_glass:
+        best_energy *= 2.0
     best_state = best_theta.copy()
     best_state[best_index] = not best_state[best_index]
 
-    return best_delta, best_state
+    return best_energy, best_state
 
 
 @njit(parallel=True, cache=True)
@@ -69,7 +59,7 @@ def run_double_bit_flips(best_theta, is_spin_glass, G_data, G_rows, G_cols, thre
     thread_batch = (combo_count + thread_count - 1) // thread_count
 
     states = np.empty((thread_count, n), dtype=np.bool_)
-    deltas = np.empty(thread_count, dtype=dtype)
+    energies = np.empty(thread_count, dtype=dtype)
 
     for t in prange(thread_count):
         s = int(t)
@@ -85,49 +75,21 @@ def run_double_bit_flips(best_theta, is_spin_glass, G_data, G_rows, G_cols, thre
                     break
             j = int(c + i + 1)
 
-            # Delta: walk adjacency lists of i and j only
-            old_i = best_theta[i]
-            new_i = not old_i
-            old_j = best_theta[j]
-            new_j = not old_j
-            delta = dtype(0)
-
-            for col in range(G_rows[i], G_rows[i + 1]):
-                v = G_cols[col]
-                if v == j:
-                    continue   # k-l edge: both flip, XOR unchanged, zero delta
-                val = G_data[col]
-                v_bit = best_theta[v]
-                if (old_i != v_bit) != (new_i != v_bit):
-                    delta += val if (new_i != v_bit) else -val
-
-            for col in range(G_rows[j], G_rows[j + 1]):
-                v = G_cols[col]
-                if v == i:
-                    continue
-                val = G_data[col]
-                v_bit = best_theta[v]
-                if (old_j != v_bit) != (new_j != v_bit):
-                    delta += val if (new_j != v_bit) else -val
-
-            if is_spin_glass:
-                delta *= dtype(2)
-
-            states[t, i] = new_i
-            states[t, j] = new_j
-            # copy other bits (only need to track i,j changes)
-            for k in range(n):
-                if k != i and k != j:
-                    states[t, k] = best_theta[k]
-            deltas[t] = delta
+            state = best_theta.copy()
+            state[i] = not state[i]
+            state[j] = not state[j]
+
+            states[t], energies[t] = state, compute_cut_diff_2_sparse(i, j, state, G_data, G_rows, G_cols, n)
 
             s += thread_batch
 
-    best_index = np.argmax(deltas)
-    best_delta = deltas[best_index]
+    best_index = np.argmax(energies)
+    best_energy = energies[best_index]
+    if is_spin_glass:
+        best_energy *= 2.0
     best_state = states[best_index]
 
-    return best_delta, best_state
+    return best_energy, best_state
 
 
 @njit(parallel=True, cache=True)
@@ -167,50 +129,30 @@ def pick_gray_seeds(best_theta, thread_count, gray_seed_multiple, G_data, G_rows
 
 
 @njit(parallel=True, cache=True)
-def run_gray_optimization(
-    best_theta,
-    iterators,
-    energies,
-    gray_iterations,
-    thread_count,
-    is_spin_glass,
-    G_data,
-    G_rows,
-    G_cols,
-):
+def run_gray_optimization(best_theta, iterators, gray_iterations, thread_count, is_spin_glass, G_data, G_rows, G_cols):
     n = len(best_theta)
     thread_iterations = (gray_iterations + thread_count - 1) // thread_count
     blocks = (n + 63) >> 6
+    energies = np.empty(thread_count, dtype=dtype)
 
     for i in prange(thread_count):
         iterator = iterators[i]
-        best_energy = energies[i]
         for curr_idx in range(thread_iterations):
+            best_energy = 0.0
             for block in range(blocks):
                 flip_bit = gray_code_next(iterator, curr_idx, block << 6)
-                # Incremental delta: walk only adjacency list of flip_bit
-                new_bit = iterator[flip_bit]
-                old_bit = not new_bit
-                delta = dtype(0)
-                for col in range(G_rows[flip_bit], G_rows[flip_bit + 1]):
-                    v = G_cols[col]
-                    val = G_data[col]
-                    v_bit = iterator[v]
-                    old_cut = old_bit != v_bit
-                    new_cut = new_bit != v_bit
-                    if old_cut != new_cut:
-                        delta += val if new_cut else -val
-                if is_spin_glass:
-                    delta *= dtype(2)
-                if delta > dtype(0):
-                    best_energy += delta
+                energy = compute_cut_diff(flip_bit, iterator, G_data, G_rows, G_cols, n)
+                if energy > 0.0:
+                    best_energy += energy
                 else:
+                    # Revert iterator
                     iterator[flip_bit] = not iterator[flip_bit]
-        if best_energy > energies[i]:
-            energies[i] = best_energy
+            energies[i] += best_energy
 
     best_index = np.argmax(energies)
     best_energy = energies[best_index]
+    if is_spin_glass:
+        best_energy *= 2.0
     best_state = iterators[best_index]
 
     return best_energy, best_state
@@ -220,6 +162,7 @@ def run_bit_flips_opencl(
     is_double,
     n,
     kernel,
+    best_energy,
     theta,
     theta_buf,
     G_data_buf,
@@ -276,12 +219,12 @@ def run_bit_flips_opencl(
     # Queue read for results we might not need
     cl.enqueue_copy(queue, max_index_host, max_index_buf)
 
-    # Find global best delta
+    # Find global minimum
     best_x = np.argmax(max_energy_host)
     energy = max_energy_host[best_x]
 
     if energy <= 0.0:
-        # No improvement: best delta is non-positive, exit early.
+        # No improvement: we can exit early
         return 0.0, theta
 
     # We need the best index
@@ -360,9 +303,9 @@ def run_gray_search_opencl(
     best_x = np.argmax(max_energy_host)
     energy = max_energy_host[best_x]
 
-    if energy <= best_energy:
+    if energy <= 0.0:
         # No improvement: we can exit early
-        return best_energy, theta
+        return 0.0, theta
 
     # We need the best index
     queue.finish()
@@ -715,6 +658,7 @@ def spin_glass_solver_sparse(
                 False,
                 n_qubits,
                 single_bit_flips_kernel,
+                max_energy,
                 best_theta,
                 theta_buf,
                 G_data_buf,
@@ -725,8 +669,8 @@ def spin_glass_solver_sparse(
             )
         else:
             energy, state = run_single_bit_flips(best_theta, is_spin_glass, G_m.data, G_m.indptr, G_m.indices)
-        if energy > 0.0:
-            max_energy += energy
+        if energy > max_energy:
+            max_energy = energy
             best_theta = state
             improved = True
             continue
@@ -738,6 +682,7 @@ def spin_glass_solver_sparse(
                 True,
                 n_qubits,
                 double_bit_flips_kernel,
+                max_energy,
                 best_theta,
                 theta_buf,
                 G_data_buf,
@@ -748,8 +693,8 @@ def spin_glass_solver_sparse(
             )
         else:
             energy, state = run_double_bit_flips(best_theta, is_spin_glass, G_m.data, G_m.indptr, G_m.indices, thread_count)
-        if energy > 0.0:
-            max_energy += energy
+        if energy > max_energy:
+            max_energy = energy
             best_theta = state
             improved = True
             continue
@@ -815,6 +760,7 @@ def spin_glass_solver_sparse(
                 False,
                 n_qubits,
                 single_bit_flips_kernel,
+                max_energy,
                 reheat_theta,
                 theta_buf,
                 G_data_buf,
@@ -825,8 +771,8 @@ def spin_glass_solver_sparse(
             )
         else:
             energy, state = run_single_bit_flips(reheat_theta, is_spin_glass, G_m.data, G_m.indptr, G_m.indices)
-        if energy > 0.0:
-            max_energy += energy
+        if energy > max_energy:
+            max_energy = energy
             best_theta = state
             improved = True
             continue
@@ -841,6 +787,7 @@ def spin_glass_solver_sparse(
                 True,
                 n_qubits,
                 double_bit_flips_kernel,
+                max_energy,
                 reheat_theta,
                 theta_buf,
                 G_data_buf,
@@ -851,8 +798,8 @@ def spin_glass_solver_sparse(
             )
         else:
             energy, state = run_double_bit_flips(reheat_theta, is_spin_glass, G_m.data, G_m.indptr, G_m.indices, thread_count)
-        if energy > 0.0:
-            max_energy += energy
+        if energy > max_energy:
+            max_energy = energy
             best_theta = state
             improved = True