Miscommunication: keep the reversed attraction/repulsion

Author: WrathfulSpatula

pyproject.toml

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pyqrackising"
-version = "9.15.1"
+version = "9.15.2"
 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/generate_tfim_samples.py

@@ -49,7 +49,8 @@
 epsilon = opencl_context.epsilon
 
 
-# Fixed-Hamming-weight Kawasaki swaps with acceptance probability, by Elara
+# ── Geometry ───────────────────────────────────────────────────────────────────
+
 def random_fixed_hamming_state(n_qubits, h):
     bits = np.zeros(n_qubits, dtype=np.bool_)
     bits[:h] = True
@@ -59,139 +60,122 @@ def random_fixed_hamming_state(n_qubits, h):
         state <<= 1
         if b:
             state |= 1
-
     return state
 
 
 def build_neighbors(n_rows, n_cols):
     L = n_rows * n_cols
     neighbors = [[] for _ in range(L)]
-
     for i in range(n_rows):
         for j in range(n_cols):
             idx = i * n_cols + j
-
             right = i * n_cols + ((j + 1) % n_cols)
-            left = i * n_cols + ((j - 1) % n_cols)
-            down = ((i + 1) % n_rows) * n_cols + j
-            up = ((i - 1) % n_rows) * n_cols + j
-
+            left  = i * n_cols + ((j - 1) % n_cols)
+            down  = ((i + 1) % n_rows) * n_cols + j
+            up    = ((i - 1) % n_rows) * n_cols + j
             neighbors[idx] = [right, left, down, up]
-
     return neighbors
 
 
-def count_like_edges(state, neighbors):
-    like = 0
+def count_unlike_edges(state, neighbors):
+    unlike = 0
     visited = set()
-
     for i, nbrs in enumerate(neighbors):
         si = (state >> i) & 1
         for j in nbrs:
             if (j, i) in visited:
                 continue
-            sj = (state >> j) & 1
-            if si == sj:
-                like += 1
+            if si != ((state >> j) & 1):
+                unlike += 1
             visited.add((i, j))
-
-    return like
+    return unlike
 
 
 def delta_like_edges(state, i, j, neighbors):
     si = (state >> i) & 1
     sj = (state >> j) & 1
-
     delta = 0
-
     for idx in [i, j]:
         s_old = (state >> idx) & 1
-
         for nbr in neighbors[idx]:
             if nbr == i or nbr == j:
                 continue
-
-            s_nbr = (state >> nbr) & 1
-
+            s_nbr    = (state >> nbr) & 1
             old_like = s_old == s_nbr
             new_spin = sj if idx == i else si
             new_like = new_spin == s_nbr
-
-            delta += int(new_like) - int(old_like)
-
-    # handle edge between i and j if neighbors
+            delta   += int(new_like) - int(old_like)
     if j in neighbors[i]:
-        old_like = si == sj
-        new_like = sj == si
-        delta += int(new_like) - int(old_like)
-
+        delta += int(sj == si) - int(si == sj)
     return delta
 
 
+# ── Spatial sampler: repulsive Kawasaki ───────────────────────────────────────
+# Electrons carry intrinsic magnetic dipole moments; antiparallel (unlike)
+# neighbours have lower magnetic potential energy.  The acceptance criterion
+# therefore unconditionally favours unlike neighbours, independent of J.
+
 def sample_fixed_hamming_weight(h_weight, count, n_rows, n_cols, burnin=10):
     if count == 0:
         return []
 
-    n_qubits = n_rows * n_cols
-    neighbors = build_neighbors(n_rows, n_cols)
-
-    like_edges = 0
-    while like_edges == 0:
-        state = random_fixed_hamming_state(n_qubits, h_weight)
-        like_edges = count_like_edges(state, neighbors)
-
-    samples = []
-
-    burn = burnin * n_qubits
+    n_qubits     = n_rows * n_cols
+    neighbors    = build_neighbors(n_rows, n_cols)
+    total_edges  = 2 * n_qubits   # 4 neighbours per site, each edge once → 2n
+
+    unlike_edges = 0
+    attempts     = 0
+    while unlike_edges == 0:
+        state        = random_fixed_hamming_state(n_qubits, h_weight)
+        unlike_edges = count_unlike_edges(state, neighbors)
+        attempts    += 1
+        if attempts > 1000:
+            return [state] * count
+
+    samples  = []
+    burn     = burnin * n_qubits
     thinning = n_qubits
 
-    ones = set([i for i in range(n_qubits) if (state >> i) & 1])
-    zeros = set([i for i in range(n_qubits) if not (state >> i) & 1])
+    ones  = set(i for i in range(n_qubits) if     (state >> i) & 1)
+    zeros = set(i for i in range(n_qubits) if not (state >> i) & 1)
 
     for step in range(burn + thinning * count):
-        # choose random 1-bit and 0-bit
-
         i = random.choice(tuple(ones))
         j = random.choice(tuple(zeros))
 
-        delta = delta_like_edges(state, i, j, neighbors)
-        new_like = like_edges + delta
+        delta_like   = delta_like_edges(state, i, j, neighbors)
+        new_unlike   = unlike_edges - delta_like   # unlike = total - like
 
-        if new_like > 0:
-            accept_prob = new_like / like_edges
-            if np.random.random() < accept_prob:
+        if new_unlike > 0:
+            # Repulsive: accept proportional to new_unlike / unlike_edges
+            if np.random.random() < new_unlike / unlike_edges:
                 state ^= 1 << i
-                ones.remove(i)
-                zeros.add(i)
-
+                ones.remove(i);  zeros.add(i)
                 state ^= 1 << j
-                zeros.remove(j)
-                ones.add(j)
-
-                like_edges = new_like
+                zeros.remove(j); ones.add(j)
+                unlike_edges = new_unlike
 
         if step >= burn and (step - burn) % thinning == 0:
             samples.append(state)
 
     return samples
 
 
+# ── Numba helpers ──────────────────────────────────────────────────────────────
+
 @njit(cache=True)
 def factor_width(width, is_transpose=False):
     col_len = math.floor(math.sqrt(width))
     while ((width // col_len) * col_len) != width:
         col_len -= 1
     row_len = width // col_len
-
     return (col_len, row_len) if is_transpose else (row_len, col_len)
 
 
 @njit(cache=True)
 def comb(n, k):
-    if (k < 0) or (k > n):
-        return 0
-    if (k == 0) or (k == n):
-        return 1
+    if (k < 0) or (k > n): return 0
+    if (k == 0) or (k == n): return 1
     res = 1
     for i in range(1, k + 1):
         res = res * (n - k + i) // i
@@ -201,13 +185,12 @@ def comb(n, k):
 @njit(cache=True)
 def fix_cdf(hamming_prob):
     tot_prob = 0.0
-    n_bias = len(hamming_prob)
+    n_bias   = len(hamming_prob)
     cum_prob = np.empty(n_bias, dtype=np.float64)
     for i in range(n_bias):
-        tot_prob += hamming_prob[i]
-        cum_prob[i] = tot_prob
+        tot_prob    += hamming_prob[i]
+        cum_prob[i]  = tot_prob
     cum_prob[-1] = 1.0
-
     return cum_prob
 
 
@@ -232,6 +215,8 @@ def get_tfim_hamming_distribution(J=-1.0, h=2.0, z=4, theta=0.174532925199432957
     return probability_by_hamming_weight(J, h, z, theta, t, n_qubits + 1, normalized=True, omega=omega)
 
 
+# ── Public generators ──────────────────────────────────────────────────────────
+
 def generate_tfim_samples(J=-1.0, h=2.0, z=4, theta=0.174532925199432957, t=5, n_qubits=56, shots=100):
     samples = []
 
@@ -244,36 +229,34 @@ def generate_tfim_samples(J=-1.0, h=2.0, z=4, theta=0.174532925199432957, t=5, n
                 if np.random.random() < prob:
                     sample |= 1
             samples.append(sample)
-
         return samples
 
     n_rows, n_cols = factor_width(n_qubits)
 
-    # First dimension: Hamming weight (envelope applied inside)
-    bias = get_tfim_hamming_distribution(J=J, h=h, z=z, theta=theta, t=t, n_qubits=n_qubits)
+    # Dimension 1: Hamming weight with time-dependent envelope
+    bias   = get_tfim_hamming_distribution(J=J, h=h, z=z, theta=theta, t=t, n_qubits=n_qubits)
     counts = np.random.multinomial(shots, bias)
 
     samples += [0] * counts[0]
     samples += [(1 << n_qubits) - 1] * counts[-1]
     if n_qubits > 1:
         samples += [int(1 << np.random.randint(n_qubits)) for _ in range(counts[1])]
     if n_qubits > 2:
-        mask = (1 << n_qubits) - 1
+        mask     = (1 << n_qubits) - 1
         samples += [(mask ^ int(1 << np.random.randint(n_qubits))) for _ in range(counts[-2])]
 
-    # Second dimension: magnetic localization
+    # Dimension 2: spatial magnetic localization (repulsive dipole coupling)
     for h_weight in range(2, len(bias) - 2):
         samples += sample_fixed_hamming_weight(h_weight, counts[h_weight], n_rows, n_cols)
 
     np.random.shuffle(samples)
-
     return samples
 
 
 def generate_fermi_hubbard_samples(J=-1.0, h=2.0, z=4, theta=0.174532925199432957, t=5, n_qubits=56, shots=100, omega=1.5 * np.pi):
     shots_x, shots_y, shots_z = np.random.multinomial(shots, [1 / 3, 1 / 3, 1 / 3])
     return (
-        generate_tfim_samples(J=J, h=h, z=z, theta=theta, t=t, n_qubits=n_qubits, shots=shots_z)
+        generate_tfim_samples(J=J,  h=h,  z=z, theta=theta,              t=t, n_qubits=n_qubits, shots=shots_z)
         + generate_tfim_samples(J=-h, h=-J, z=z, theta=theta + np.pi / 2, t=t, n_qubits=n_qubits, shots=shots_x)
-        + generate_tfim_samples(J=J, h=h, z=z, theta=theta + np.pi / 2, t=t, n_qubits=n_qubits, shots=shots_y)
+        + generate_tfim_samples(J=J,  h=h,  z=z, theta=theta + np.pi / 2, t=t, n_qubits=n_qubits, shots=shots_y)
     )

setup.py

@@ -7,7 +7,7 @@
 
 setup(
     name="pyqrackising",
-    version="9.15.1",
+    version="9.15.2",
     author="Dan Strano",
     author_email="stranoj@gmail.com",
     description="Fast MAXCUT, TSP, and sampling heuristics from near-ideal transverse field Ising model (TFIM)",