Compute cut differential on CPU

Author: WrathfulSpatula

pyproject.toml

@@ -10,7 +10,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pyqrackising"
-version = "9.7.8"
+version = "9.8.0"
 requires-python = ">=3.8"
 description = "Fast MAXCUT, TSP, and sampling heuristics from near-ideal transverse field Ising model (TFIM)"
 readme = {file = "README.txt", content-type = "text/markdown"}

pyqrackising/maxcut_tfim_util.py

@@ -278,7 +278,7 @@ def convert_bool_to_uint(samples):
 
 @njit
 def compute_energy(sample, G_m, n_qubits):
-    energy = 0
+    energy = 0.0
     for u in range(n_qubits):
         u_bit = sample[u]
         for v in range(u + 1, n_qubits):
@@ -299,9 +299,93 @@ def compute_cut(sample, G_m, n_qubits):
     return cut
 
 
+@njit
+def compute_energy_diff(u, sample, G_m, n_qubits):
+    energy = 0.0
+    u_bit = sample[u]
+    for v in range(u):
+        val = 2 * G_m[u, v]
+        energy += val if u_bit == sample[v] else -val
+    for v in range(u + 1, n_qubits):
+        val = 2 * G_m[u, v]
+        energy += val if u_bit == sample[v] else -val
+
+    return -energy
+
+
+@njit
+def compute_cut_diff(u, sample, G_m, n_qubits):
+    energy = 0.0
+    u_bit = sample[u]
+    for v in range(u):
+        val = G_m[u, v]
+        energy += -val if u_bit == sample[v] else val
+    for v in range(u + 1, n_qubits):
+        val = G_m[u, v]
+        energy += -val if u_bit == sample[v] else val
+
+    return energy
+
+
+@njit
+def compute_energy_diff_2(k, l, sample, G_m, n_qubits):
+    if l < k:
+        t = k
+        k = l
+        l = t
+    energy = 0.0
+    k_bit = sample[k]
+    l_bit = sample[l]
+    for v in range(k):
+        val = 2 * G_m[k, v]
+        energy += val if k_bit == sample[v] else -val
+        val = 2 * G_m[l, v]
+        energy += val if l_bit == sample[v] else -val
+    for v in range(k + 1, l):
+        val = 2 * G_m[k, v]
+        energy += val if k_bit == sample[v] else -val
+        val = 2 * G_m[l, v]
+        energy += val if l_bit == sample[v] else -val
+    for v in range(l + 1, n_qubits):
+        val = 2 * G_m[k, v]
+        energy += val if k_bit == sample[v] else -val
+        val = 2 * G_m[l, v]
+        energy += val if l_bit == sample[v] else -val
+
+    return -energy
+
+
+@njit
+def compute_cut_diff_2(k, l, sample, G_m, n_qubits):
+    if l < k:
+        t = k
+        k = l
+        l = t
+    energy = 0.0
+    k_bit = sample[k]
+    l_bit = sample[l]
+    for v in range(k):
+        val = G_m[k, v]
+        energy += -val if k_bit == sample[v] else val
+        val = G_m[l, v]
+        energy += -val if l_bit == sample[v] else val
+    for v in range(k + 1, l):
+        val = G_m[k, v]
+        energy += -val if k_bit == sample[v] else val
+        val = G_m[l, v]
+        energy += -val if l_bit == sample[v] else val
+    for v in range(l + 1, n_qubits):
+        val = G_m[k, v]
+        energy += -val if k_bit == sample[v] else val
+        val = G_m[l, v]
+        energy += -val if l_bit == sample[v] else val
+
+    return energy
+
+
 @njit
 def compute_energy_sparse(sample, G_data, G_rows, G_cols, n_qubits):
-    energy = 0
+    energy = 0.0
     for u in range(n_qubits):
         u_bit = sample[u]
         for col in range(G_rows[u], G_rows[u + 1]):
@@ -327,7 +411,7 @@ def compute_cut_sparse(sample, G_data, G_rows, G_cols, n_qubits):
 
 @njit
 def compute_energy_streaming(sample, G_func, nodes, n_qubits):
-    energy = 0
+    energy = 0.0
     for u in range(n_qubits):
         u_bit = sample[u]
         for v in range(u + 1, n_qubits):

pyqrackising/spin_glass_solver.py

@@ -1,5 +1,5 @@
 from .maxcut_tfim import maxcut_tfim
-from .maxcut_tfim_util import compute_cut, compute_energy, get_cut, gray_code_next, gray_mutation, heuristic_threshold, int_to_bitstring, make_G_m_buf, make_best_theta_buf, make_best_theta_buf_64, opencl_context, setup_opencl
+from .maxcut_tfim_util import compute_cut, compute_cut_diff, compute_cut_diff_2, compute_energy, compute_energy_diff, compute_energy_diff_2, get_cut, gray_code_next, gray_mutation, heuristic_threshold, int_to_bitstring, make_G_m_buf, make_best_theta_buf, make_best_theta_buf_64, opencl_context, setup_opencl
 import networkx as nx
 import numpy as np
 from numba import njit, prange
@@ -28,12 +28,12 @@ def run_single_bit_flips(best_theta, is_spin_glass, G_m):
         for i in prange(n):
             state = best_theta.copy()
             state[i] = not state[i]
-            energies[i] = compute_energy(state, G_m, n)
+            energies[i] = compute_energy_diff(i, state, G_m, n)
     else:
         for i in prange(n):
             state = best_theta.copy()
             state[i] = not state[i]
-            energies[i] = compute_cut(state, G_m, n)
+            energies[i] = compute_cut_diff(i, state, G_m, n)
 
     best_index = np.argmax(energies)
     best_energy = energies[best_index]
@@ -71,7 +71,7 @@ def run_double_bit_flips(best_theta, is_spin_glass, G_m, thread_count):
                 state[i] = not state[i]
                 state[j] = not state[j]
 
-                states[t], energies[t] = state, compute_energy(state, G_m, n)
+                states[t], energies[t] = state, compute_energy_diff_2(i, j, state, G_m, n)
 
                 s += thread_batch
     else:
@@ -93,7 +93,7 @@ def run_double_bit_flips(best_theta, is_spin_glass, G_m, thread_count):
                 state[i] = not state[i]
                 state[j] = not state[j]
 
-                states[t], energies[t] = state, compute_cut(state, G_m, n)
+                states[t], energies[t] = state, compute_cut_diff_2(i, j, state, G_m, n)
 
                 s += thread_batch
 
@@ -184,7 +184,7 @@ def run_gray_optimization(best_theta, iterators, energies, gray_iterations, thre
     return best_energy, best_state
 
 
-def run_bit_flips_opencl(is_double, n, kernel, best_energy, theta, theta_buf, G_m_buf, is_segmented, local_size, global_size, args_buf, local_energy_buf, local_index_buf, max_energy_host, max_index_host, max_energy_buf, max_index_buf):
+def run_bit_flips_opencl(is_double, n, kernel, theta, theta_buf, G_m_buf, is_segmented, local_size, global_size, args_buf, local_energy_buf, local_index_buf, max_energy_host, max_index_host, max_energy_buf, max_index_buf):
     queue = opencl_context.queue
 
     # Set kernel args
@@ -227,9 +227,7 @@ def run_bit_flips_opencl(is_double, n, kernel, best_energy, theta, theta_buf, G_
 
     if energy <= 0.0:
         # No improvement: we can exit early
-        return best_energy, theta
-
-    energy += best_energy
+        return 0.0, theta
 
     # We need the best index
     queue.finish()
@@ -425,23 +423,23 @@ def spin_glass_solver(
         # Single bit flips with O(n^2)
         if is_opencl:
             theta_buf = make_best_theta_buf(best_theta)
-            energy, state = run_bit_flips_opencl(False, n_qubits, single_bit_flips_kernel, max_energy, best_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
+            energy, state = run_bit_flips_opencl(False, n_qubits, single_bit_flips_kernel, best_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
         else:
             energy, state = run_single_bit_flips(best_theta, is_spin_glass, G_m)
-        if energy > max_energy:
-            max_energy = energy
+        if energy > 0.0:
+            max_energy += energy
             best_theta = state
             improved = True
             continue
 
         # Double bit flips with O(n^3)
         if is_opencl:
             # theta_buf has not changed
-            energy, state = run_bit_flips_opencl(True, n_qubits, double_bit_flips_kernel, max_energy, best_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
+            energy, state = run_bit_flips_opencl(True, n_qubits, double_bit_flips_kernel, best_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
         else:
             energy, state = run_double_bit_flips(best_theta, is_spin_glass, G_m, thread_count)
-        if energy > max_energy:
-            max_energy = energy
+        if energy > 0.0:
+            max_energy += energy
             best_theta = state
             improved = True
             continue
@@ -472,23 +470,23 @@ def spin_glass_solver(
         # Single bit flips with O(n^2)
         if is_opencl:
             theta_buf = make_best_theta_buf(reheat_theta)
-            energy, state = run_bit_flips_opencl(False, n_qubits, single_bit_flips_kernel, max_energy, reheat_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
+            energy, state = run_bit_flips_opencl(False, n_qubits, single_bit_flips_kernel, reheat_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
         else:
             energy, state = run_single_bit_flips(reheat_theta, is_spin_glass, G_m)
-        if energy > max_energy:
-            max_energy = energy
+        if energy > 0.0:
+            max_energy += energy
             best_theta = state
             improved = True
             continue
 
         # Double bit flips with O(n^3)
         if is_opencl:
             # theta_buf has not changed
-            energy, state = run_bit_flips_opencl(True, n_qubits, double_bit_flips_kernel, max_energy, reheat_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
+            energy, state = run_bit_flips_opencl(True, n_qubits, double_bit_flips_kernel, reheat_theta, theta_buf, G_m_buf, is_segmented, *opencl_args)
         else:
             energy, state = run_double_bit_flips(reheat_theta, is_spin_glass, G_m, thread_count)
-        if energy > max_energy:
-            max_energy = energy
+        if energy > 0.0:
+            max_energy += energy
             best_theta = state
             improved = True
 

setup.py

@@ -7,7 +7,7 @@
 
 setup(
     name='pyqrackising',
-    version='9.7.8',
+    version='9.8.0',
     author='Dan Strano',
     author_email='[email protected]',
     description='Fast MAXCUT, TSP, and sampling heuristics from near-ideal transverse field Ising model (TFIM)',