Smooth transition from exact MAXCUT (orig. from Elara, custom GPT) to heuristic

WrathfulSpatula · WrathfulSpatula · commit d775876ffe3a · 2025-11-10T16:49:45.000-05:00
diff --git a/pyproject.toml b/pyproject.toml
@@ -10,7 +10,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pyqrackising"
-version = "9.7.0"
+version = "9.7.1"
 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"}
diff --git a/pyqrackising/maxcut_tfim.py b/pyqrackising/maxcut_tfim.py
@@ -3,7 +3,7 @@
 from numba import njit, prange
 import os
 
-from .maxcut_tfim_util import compute_cut, compute_energy, convert_bool_to_uint, get_cut, get_cut_base, init_thresholds, make_G_m_buf, make_theta_buf, maxcut_hamming_cdf, opencl_context, sample_mag, setup_opencl, bit_pick
+from .maxcut_tfim_util import compute_cut, compute_energy, convert_bool_to_uint, get_cut, get_cut_base, heuristic_threshold, init_thresholds, make_G_m_buf, make_theta_buf, maxcut_hamming_cdf, opencl_context, sample_mag, setup_opencl, bit_pick
 
 IS_OPENCL_AVAILABLE = True
 try:
@@ -242,6 +242,60 @@ def run_cut_opencl(best_energy, samples, G_m_buf, is_segmented, local_size, glob
     return samples[max_index_host[best_x]], energy
 
 
+@njit
+def exact_maxcut(G):
+    """Brute-force exact MAXCUT solver using Numba JIT."""
+    n = G.shape[0]
+    max_cut = -1.0
+    best_mask = 0
+
+    # Enumerate all 2^n possible bitstrings
+    for mask in range(1 << n):
+        cut = 0.0
+        for i in range(n):
+            bi = (mask >> i) & 1
+            for j in range(i + 1, n):
+                if bi != ((mask >> j) & 1):
+                    cut += G[i, j]
+        if cut > max_cut:
+            max_cut = cut
+            best_mask = mask
+
+    # Reconstruct best bitstring
+    best_bits = np.zeros(n, dtype=np.bool_)
+    for i in range(n):
+        best_bits[i] = (best_mask >> i) & 1
+
+    return best_bits, max_cut
+
+
+@njit
+def exact_spin_glass(G):
+    """Brute-force exact spin-glass solver using Numba JIT."""
+    n = G.shape[0]
+    max_cut = -1.0
+    best_mask = 0
+
+    # Enumerate all 2^n possible bitstrings
+    for mask in range(1 << n):
+        cut = 0.0
+        for i in range(n):
+            bi = (mask >> i) & 1
+            for j in range(i + 1, n):
+                val = G[i, j]
+                cut += val if bi == ((mask >> j) & 1) else -val
+        if cut > max_cut:
+            max_cut = cut
+            best_mask = mask
+
+    # Reconstruct best bitstring
+    best_bits = np.zeros(n, dtype=np.bool_)
+    for i in range(n):
+        best_bits[i] = (best_mask >> i) & 1
+
+    return best_bits, max_cut
+
+
 def maxcut_tfim(
     G,
     quality=None,
@@ -264,22 +318,15 @@ def maxcut_tfim(
 
     n_qubits = len(G_m)
 
-    if n_qubits < 3:
-        empty = [nodes[0]]
-        empty.clear()
-
-        if n_qubits == 0:
-            return "", 0, (empty, empty.copy())
-
-        if n_qubits == 1:
-            return "0", 0, (nodes, empty)
+    if n_qubits < heuristic_threshold:
+        best_solution, best_value = exact_spin_glass(G_m) if is_spin_glass else exact_maxcut(G_m)
+        bit_string, l, r = get_cut(best_solution, nodes, n_qubits)
 
-        if n_qubits == 2:
-            weight = G_m[0, 1]
-            if weight < 0.0:
-                return "00", 0, (nodes, empty)
+        if best_value < 0.0:
+            # Best cut is trivial partition, all/empty
+            return '0' * n_qubits, 0.0, (nodes, [])
 
-            return "01", weight, ([nodes[0]], [nodes[1]])
+        return bit_string, best_value, (l, r)
 
     if quality is None:
         quality = 6
diff --git a/pyqrackising/maxcut_tfim_sparse.py b/pyqrackising/maxcut_tfim_sparse.py
@@ -3,7 +3,7 @@
 from numba import njit, prange
 import os
 
-from .maxcut_tfim_util import binary_search, compute_cut_sparse, compute_energy_sparse, convert_bool_to_uint, get_cut, get_cut_base, init_thresholds, make_G_m_csr_buf, make_theta_buf, maxcut_hamming_cdf, opencl_context, sample_mag, setup_opencl, bit_pick, to_scipy_sparse_upper_triangular
+from .maxcut_tfim_util import binary_search, compute_cut_sparse, compute_energy_sparse, convert_bool_to_uint, get_cut, get_cut_base, heuristic_threshold_sparse, init_thresholds, make_G_m_csr_buf, make_theta_buf, maxcut_hamming_cdf, opencl_context, sample_mag, setup_opencl, bit_pick, to_scipy_sparse_upper_triangular
 
 IS_OPENCL_AVAILABLE = True
 try:
@@ -256,6 +256,69 @@ def run_cut_opencl(best_energy, samples, G_data_buf, G_rows_buf, G_cols_buf, is_
     return samples[max_index_host[best_x]], energy
 
 
+@njit
+def exact_maxcut(G_data, G_rows, G_cols):
+    """Brute-force exact MAXCUT solver using Numba JIT."""
+    n = G_rows.shape[0] - 1
+    max_cut = -1.0
+    best_mask = 0
+
+    # Enumerate all 2^n possible bitstrings
+    for mask in range(1 << n):
+        cut = 0.0
+        for i in range(n):
+            bi = (mask >> i) & 1
+            for j in range(i + 1, n):
+                if bi != ((mask >> j) & 1):
+                    u, v = (i, j) if i < j else (j, i)
+                    start = G_rows[u]
+                    end = G_rows[u + 1]
+                    k = binary_search(G_cols[start:end], v) + start
+                    if k < end:
+                        cut += G_data[k]
+        if cut > max_cut:
+            max_cut = cut
+            best_mask = mask
+
+    # Reconstruct best bitstring
+    best_bits = np.zeros(n, dtype=np.bool_)
+    for i in range(n):
+        best_bits[i] = (best_mask >> i) & 1
+
+    return best_bits, max_cut
+
+
+@njit
+def exact_spin_glass(G_data, G_rows, G_cols):
+    """Brute-force exact spin-glass solver using Numba JIT."""
+    n = G_rows.shape[0] - 1
+    max_cut = -1.0
+    best_mask = 0
+
+    # Enumerate all 2^n possible bitstrings
+    for mask in range(1 << n):
+        cut = 0.0
+        for i in range(n):
+            bi = (mask >> i) & 1
+            for j in range(i + 1, n):
+                u, v = (i, j) if i < j else (j, i)
+                start = G_rows[u]
+                end = G_rows[u + 1]
+                k = binary_search(G_cols[start:end], v) + start
+                if k < end:
+                    val = G_data[k]
+                    cut += val if bi == ((mask >> j) & 1) else -val
+        if cut > max_cut:
+            max_cut = cut
+            best_mask = mask
+
+    # Reconstruct best bitstring
+    best_bits = np.zeros(n, dtype=np.bool_)
+    for i in range(n):
+        best_bits[i] = (best_mask >> i) & 1
+
+    return best_bits, max_cut
+
 
 def maxcut_tfim_sparse(
     G,
@@ -281,19 +344,15 @@ def maxcut_tfim_sparse(
         nodes = list(range(n_qubits))
         G_m = G
 
-    if n_qubits < 3:
-        if n_qubits == 0:
-            return "", 0, ([], [])
-
-        if n_qubits == 1:
-            return "0", 0, (nodes, [])
+    if n_qubits < heuristic_threshold_sparse:
+        best_solution, best_value = exact_spin_glass(G_m.data, G_m.indptr, G_m.indices) if is_spin_glass else exact_maxcut(G_m.data, G_m.indptr, G_m.indices)
+        bit_string, l, r = get_cut(best_solution, nodes, n_qubits)
 
-        if n_qubits == 2:
-            weight = G_m[0, 1]
-            if weight < 0.0:
-                return "00", 0, (nodes, [])
+        if best_value < 0.0:
+            # Best cut is trivial partition, all/empty
+            return '0' * n_qubits, 0.0, (nodes, [])
 
-            return "01", weight, ([nodes[0]], [nodes[1]])
+        return bit_string, best_value, (l, r)
 
     if quality is None:
         quality = 6
diff --git a/pyqrackising/maxcut_tfim_streaming.py b/pyqrackising/maxcut_tfim_streaming.py
@@ -3,7 +3,7 @@
 from numba import njit, prange
 import os
 
-from .maxcut_tfim_util import compute_cut_streaming, compute_energy_streaming, get_cut, get_cut_base, init_thresholds, maxcut_hamming_cdf, opencl_context, sample_mag, bit_pick
+from .maxcut_tfim_util import compute_cut_streaming, compute_energy_streaming, get_cut, get_cut_base, heuristic_threshold, init_thresholds, maxcut_hamming_cdf, opencl_context, sample_mag, bit_pick
 
 
 epsilon = opencl_context.epsilon
@@ -148,6 +148,58 @@ def find_G_min(G_func, nodes, n_nodes):
     return G_min
 
 
+@njit
+def exact_maxcut(G_func, n):
+    """Brute-force exact MAXCUT solver using Numba JIT."""
+    max_cut = -1.0
+    best_mask = 0
+
+    # Enumerate all 2^n possible bitstrings
+    for mask in range(1 << n):
+        cut = 0.0
+        for i in range(n):
+            bi = (mask >> i) & 1
+            for j in range(i + 1, n):
+                if bi != ((mask >> j) & 1):
+                    cut += G_func(i, j)
+        if cut > max_cut:
+            max_cut = cut
+            best_mask = mask
+
+    # Reconstruct best bitstring
+    best_bits = np.zeros(n, dtype=np.bool_)
+    for i in range(n):
+        best_bits[i] = (best_mask >> i) & 1
+
+    return best_bits, max_cut
+
+
+@njit
+def exact_spin_glass(G_func, n):
+    """Brute-force exact spin-glass solver using Numba JIT."""
+    max_cut = -1.0
+    best_mask = 0
+
+    # Enumerate all 2^n possible bitstrings
+    for mask in range(1 << n):
+        cut = 0.0
+        for i in range(n):
+            bi = (mask >> i) & 1
+            for j in range(i + 1, n):
+                val = G_func(i, j)
+                cut += val if bi == ((mask >> j) & 1) else -val
+        if cut > max_cut:
+            max_cut = cut
+            best_mask = mask
+
+    # Reconstruct best bitstring
+    best_bits = np.zeros(n, dtype=np.bool_)
+    for i in range(n):
+        best_bits[i] = (best_mask >> i) & 1
+
+    return best_bits, max_cut
+
+
 def maxcut_tfim_streaming(
     G_func,
     nodes,
@@ -161,19 +213,15 @@ def maxcut_tfim_streaming(
     wgs = opencl_context.work_group_size
     n_qubits = len(nodes)
 
-    if n_qubits < 3:
-        if n_qubits == 0:
-            return "", 0, ([], [])
-
-        if n_qubits == 1:
-            return "0", 0, (nodes, [])
+    if n_qubits < heuristic_threshold:
+        best_solution, best_value = exact_spin_glass(G_func, n_qubits) if is_spin_glass else exact_maxcut(G_func, n_qubits)
+        bit_string, l, r = get_cut(best_solution, nodes, n_qubits)
 
-        if n_qubits == 2:
-            weight = G_func(nodes[0], nodes[1])
-            if weight < 0.0:
-                return "00", 0, (nodes, [])
+        if best_value < 0.0:
+            # Best cut is trivial partition, all/empty
+            return '0' * n_qubits, 0.0, (nodes, [])
 
-            return "01", weight, ([nodes[0]], [nodes[1]])
+        return bit_string, best_value, (l, r)
 
     if quality is None:
         quality = 6
diff --git a/pyqrackising/maxcut_tfim_util.py b/pyqrackising/maxcut_tfim_util.py
@@ -121,6 +121,8 @@ def __init__(self, a, b, g, d, e, f, c, q, i, j, k, l, m, n, o, p, x, y, z, w):
     print("PyOpenCL not installed. (If you have any OpenCL accelerator devices with available ICDs, you might want to optionally install pyopencl.)")
 
 opencl_context = OpenCLContext(compute_units, IS_OPENCL_AVAILABLE, work_group_size, dtype, epsilon, max_alloc, ctx, queue, calculate_cut_kernel, calculate_cut_sparse_kernel, calculate_cut_segmented_kernel, calculate_cut_sparse_segmented_kernel, single_bit_flips_kernel, single_bit_flips_sparse_kernel, single_bit_flips_segmented_kernel, single_bit_flips_sparse_segmented_kernel, double_bit_flips_kernel, double_bit_flips_sparse_kernel, double_bit_flips_segmented_kernel, double_bit_flips_sparse_segmented_kernel)
+heuristic_threshold = 24
+heuristic_threshold_sparse = 23
 
 
 def setup_opencl(l, g, args_np):
diff --git a/pyqrackising/spin_glass_solver.py b/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, int_to_bitstring, make_G_m_buf, make_best_theta_buf, opencl_context, setup_opencl
+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, opencl_context, setup_opencl
 import networkx as nx
 import numpy as np
 from numba import njit, prange
@@ -296,6 +296,9 @@ def spin_glass_solver(
 
             return "01", weight, ([nodes[0]], [nodes[1]]), -weight
 
+    if n_qubits < heuristic_threshold:
+        best_guess = None
+
     bitstring = ""
     if isinstance(best_guess, str):
         bitstring = best_guess
@@ -305,16 +308,27 @@ def spin_glass_solver(
         bitstring = "".join(["1" if b else "0" for b in best_guess])
     else:
         bitstring, cut_value, _ = maxcut_tfim(G_m, quality=quality, shots=shots, is_spin_glass=is_spin_glass, anneal_t=anneal_t, anneal_h=anneal_h, repulsion_base=repulsion_base, is_maxcut_gpu=is_maxcut_gpu, is_nested=True)
+
     best_theta = np.array([b == "1" for b in list(bitstring)], dtype=np.bool_)
+    max_energy = compute_energy(best_theta, G_m, n_qubits) if is_spin_glass else cut_value
+
+    if n_qubits < heuristic_threshold:
+        bitstring, l, r = get_cut(best_theta, nodes, n_qubits)
+        if is_spin_glass:
+            cut_value = compute_cut(best_theta, G_m, n_qubits)
+            min_energy = -max_energy
+        else:
+            cut_value = max_energy
+            min_energy = compute_energy(best_theta, G_m, n_qubits)
+
+        return bitstring, float(cut_value), (l, r), float(min_energy)
 
     if gray_iterations is None:
         gray_iterations = n_qubits * n_qubits
 
     if gray_seed_multiple is None:
         gray_seed_multiple = os.cpu_count()
 
-    max_energy = compute_energy(best_theta, G_m, n_qubits) if is_spin_glass else cut_value
-
     is_opencl = is_maxcut_gpu and IS_OPENCL_AVAILABLE
 
     if is_opencl:
diff --git a/pyqrackising/spin_glass_solver_sparse.py b/pyqrackising/spin_glass_solver_sparse.py
diff --git a/pyqrackising/spin_glass_solver_streaming.py b/pyqrackising/spin_glass_solver_streaming.py
diff --git a/setup.py b/setup.py
