add paoa

README.md

@@ -123,4 +123,29 @@ class LogicCircuit:
         self.and_gate.output.connect(self.or_gate.input1)
 ```
 
+### Higher-level Library
+
+`p_kit.library` provides ready-made circuits for common problem types.
+
+| Class | Purpose | Extra p-bits mean… |
+|---|---|---|
+| `TSP` | Travelling Salesman Problem | one p-bit per (city, order) pair |
+| `PolyOptimizer` | Classical polynomial minimisation (QUBO) | binary bits for integer precision |
+| `TransverseFieldIsing` | Quantum-inspired annealing (TFIM emulation) | Trotter replicas in imaginary time |
+
+Both `PolyOptimizer` and `TransverseFieldIsing` can target the same QUBO/Ising problems, but they correspond to different parts of QAOA:
+
+- **`PolyOptimizer`** encodes only the QAOA *cost Hamiltonian* H_C = -Σ J_ij σ^z_i σ^z_j. It uses pure Boltzmann sampling to find the minimum — equivalent to QAOA with zero mixer strength (the classical baseline).
+
+- **`TransverseFieldIsing`** emulates the full QAOA *dynamics*: the cost Hamiltonian (intra-replica couplings) plus the mixer Hamiltonian H_B = -Γ Σ σ^x_i (inter-replica Trotter couplings). The `gamma` parameter is the QAOA mixer strength and `n_replicas` corresponds to the number of QAOA layers. This is what p-kit calls **PAOA** — a probabilistic emulation of QAOA for stoquastic Hamiltonians, running at room temperature on classical hardware (see [Camsari et al., 2019](https://arxiv.org/abs/1809.04028), Fig. 9).
+
+To use `TransverseFieldIsing` with a QUBO problem:
+
+```python
+from p_kit.library.quantum import TransverseFieldIsing
+
+# Q is your (n, n) QUBO matrix
+circuit = TransverseFieldIsing.from_qubo(Q, gamma=0.5, beta=2.0, n_replicas=20)
+```
+
 For more information, visit our [documentation](https://github.com/IBM/p-kit/wiki).

examples/docplex_poly.py

@@ -1,5 +1,8 @@
-from docplex.mp.model import Model
-from p_kit.library.docplex_optimizer import DocplexOptimizer
+try:
+    from docplex.mp.model import Model
+    from p_kit.library.docplex_optimizer import DocplexOptimizer
+except ImportError:
+    raise SystemExit(0)  # skip gracefully if docplex is not installed
 from p_kit.solver.csd_solver import CaSuDaSolver
 from p_kit.solver.annealing import constant, execute
 from p_kit.visualization.histplot import histplot

examples/quantum_tfim.py

@@ -0,0 +1,69 @@
+# Transverse-Field Ising Model (TFIM) emulated with p-bits
+# via the Suzuki-Trotter path integral decomposition.
+#
+# Ref: Camsari et al., arXiv:1809.04028, Fig. 9
+#
+# A stoquastic n-qubit TFIM is mapped to a p-circuit with
+# n * n_replicas p-bits. Intra-replica couplings reproduce
+# the Ising Hamiltonian; inter-replica couplings encode the
+# transverse field via the Trotter coupling K.
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from p_kit.library.quantum import TransverseFieldIsing
+from p_kit.solver.annealing import linear, execute
+from p_kit.solver.csd_solver import CaSuDaSolver
+
+# --- Problem definition --------------------------------------------------
+# 3-qubit ferromagnetic TFIM: J_ij = 1 (spins want to align)
+n = 3
+J_q = np.array([
+    [0., 1., 1.],
+    [1., 0., 1.],
+    [1., 1., 0.],
+])
+h_q = np.zeros(n)
+
+# Transverse field strength (quantum fluctuations).
+# gamma=0  → classical Ising (all replicas lock together)
+# gamma>>0 → strong quantum fluctuations
+gamma = 0.5
+beta = 2.0    # inverse temperature
+n_replicas = 20
+
+# --- Build p-circuit from Trotter decomposition --------------------------
+circuit = TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=n_replicas)
+print(f"Qubits: {n}  |  Replicas: {n_replicas}  |  Total p-bits: {circuit.n_pbits}")
+print(f"Trotter coupling K = {-0.5 * np.log(np.tanh(beta * gamma / n_replicas)):.4f}")
+
+# --- Solve with annealing ------------------------------------------------
+solver = CaSuDaSolver(Nt=int(5e4), dt=0.1, i0=0, seed=42)
+n_shots = 100
+samples = execute(solver, circuit, linear, n_shots=n_shots, n_last_samples=50, n_jobs=-1)
+
+# Recover qubit magnetizations by averaging over replicas
+# samples shape: (n_shots * n_last_samples, n * n_replicas)
+qubit_m = samples.reshape(-1, n_replicas, n)  # (samples, replicas, qubits)
+qubit_m_avg = qubit_m.mean(axis=1)            # (samples, qubits): per-qubit mean over replicas
+
+# --- Plot ----------------------------------------------------------------
+fig, axes = plt.subplots(1, 2, figsize=(12, 4))
+fig.suptitle(f"TFIM emulated via Suzuki-Trotter  (n={n}, R={n_replicas}, γ={gamma}, β={beta})")
+
+# Left: histogram of qubit-0 averaged magnetization
+axes[0].hist(qubit_m_avg[:, 0], bins=30, density=True, color="steelblue", edgecolor="white")
+axes[0].set_xlabel("⟨m₀⟩ (averaged over replicas)")
+axes[0].set_ylabel("Probability density")
+axes[0].set_title("Qubit-0 magnetization distribution")
+
+# Right: inter-qubit correlation <m_i * m_j>
+corr = np.mean(qubit_m_avg[:, :, None] * qubit_m_avg[:, None, :], axis=0)
+im = axes[1].imshow(corr, vmin=-1, vmax=1, cmap="RdBu_r")
+axes[1].set_title("Qubit-qubit correlations ⟨mᵢ mⱼ⟩")
+axes[1].set_xticks(range(n))
+axes[1].set_yticks(range(n))
+plt.colorbar(im, ax=axes[1])
+
+plt.tight_layout()
+plt.show()

p_kit/library/__init__.py

@@ -1,3 +1,6 @@
 from . import tsp
+from . import poly
+from . import docplex_optimizer
+from . import quantum
 
-__all__ = ["tsp"]
+__all__ = ["tsp", "poly", "docplex_optimizer", "quantum"]

p_kit/library/quantum.py

@@ -0,0 +1,96 @@
+import numpy as np
+from p_kit.psl.p_circuit import PCircuit
+
+
+class TransverseFieldIsing(PCircuit):
+    """Maps a stoquastic transverse-field Ising model to a p-bit circuit
+    via the Suzuki-Trotter (path integral) decomposition.
+
+    An n-qubit quantum system becomes an n * n_replicas p-bit circuit.
+    P-bit (i, tau) maps to index tau * n + i (replica-major ordering).
+
+    This is the probabilistic analog of QAOA (PAOA): intra-replica couplings
+    encode the QAOA cost Hamiltonian H_C = -sum_ij J_ij sz_i sz_j, while
+    inter-replica Trotter couplings encode the QAOA mixer H_B = -gamma sum_i sx_i.
+    ``gamma`` is the mixer strength and ``n_replicas`` the number of QAOA layers.
+    Use ``PolyOptimizer`` instead for the classical baseline (gamma=0 equivalent).
+
+    Parameters
+    ----------
+    j_q : np.ndarray, shape (n, n)
+        Symmetric Ising coupling matrix of the quantum system (QAOA cost H_C).
+    h_q : np.ndarray, shape (n,)
+        Longitudinal field (bias) of the quantum system.
+    gamma : float
+        Transverse field / QAOA mixer strength. Use 0 for a purely classical Ising model.
+    beta : float, default 1.0
+        Inverse temperature (1/T).
+    n_replicas : int, default 10
+        Number of Trotter replicas (QAOA layers). Higher values give a better approximation.
+
+    Notes
+    -----
+    Ref: Camsari et al., arXiv:1809.04028, Fig. 9
+
+    .. versionadded:: 0.0.1
+    """
+
+    def __init__(self, j_q, h_q, gamma, beta=1.0, n_replicas=10):
+        n = len(h_q)
+        R = n_replicas
+        super().__init__(n * R)
+
+        J = np.zeros((n * R, n * R))
+        h_vec = np.zeros(n * R)
+
+        # Intra-replica: spatial couplings and bias scaled by beta/R
+        scale = beta / R
+        for tau in range(R):
+            s = tau * n
+            J[s:s + n, s:s + n] += scale * j_q
+            h_vec[s:s + n] += scale * h_q
+
+        # Inter-replica: Trotter coupling K = -(1/2) ln(tanh(beta*gamma/R))
+        # K > 0 ferromagnetically couples same qubit across adjacent replicas.
+        # gamma=0 means classical limit (no inter-replica coupling).
+        if gamma > 0:
+            K = -0.5 * np.log(np.tanh(beta * gamma / R))
+            seen = set()
+            for tau in range(R):
+                tau_next = (tau + 1) % R
+                pair = (min(tau, tau_next), max(tau, tau_next))
+                if pair in seen:
+                    continue
+                seen.add(pair)
+                idx = np.arange(n)
+                J[tau * n + idx, tau_next * n + idx] += K
+                J[tau_next * n + idx, tau * n + idx] += K
+
+        np.fill_diagonal(J, 0)
+        self.J = J
+        self.h = h_vec
+
+    @classmethod
+    def from_qubo(cls, qubo, gamma, beta=1.0, n_replicas=10):
+        """Build a TransverseFieldIsing circuit from a QUBO matrix.
+
+        Minimises x^T Q x with x ∈ {0, 1} by mapping to the Ising
+        spin variables s ∈ {-1, +1} via x_i = (1 + s_i) / 2.
+
+        Parameters
+        ----------
+        qubo : np.ndarray, shape (n, n)
+            QUBO matrix (upper or full triangle).
+        gamma : float
+            Transverse field strength.
+        beta : float, default 1.0
+            Inverse temperature (1/T).
+        n_replicas : int, default 10
+            Number of Trotter replicas.
+        """
+        qubo = np.asarray(qubo, dtype=float)
+        qubo_sym = qubo + qubo.T - np.diag(np.diag(qubo))   # symmetrise, count diagonal once
+        j_q = -qubo_sym / 4
+        np.fill_diagonal(j_q, 0)
+        h_q = -np.sum(qubo_sym, axis=1) / 4
+        return cls(j_q, h_q, gamma, beta, n_replicas)

tests/test_quantum.py

@@ -0,0 +1,82 @@
+import numpy as np
+import pytest
+from p_kit.library.quantum import TransverseFieldIsing
+from p_kit.solver import CaSuDaSolver
+
+
+def _make(n=3, R=4, gamma=0.5, beta=1.0, seed=0):
+    rng = np.random.default_rng(seed)
+    J_q = rng.uniform(-1, 1, (n, n))
+    J_q = (J_q + J_q.T) / 2
+    np.fill_diagonal(J_q, 0)
+    h_q = rng.uniform(-1, 1, n)
+    return TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=R), n, R
+
+
+def test_shape():
+    c, n, R = _make()
+    assert c.n_pbits == n * R
+    assert c.J.shape == (n * R, n * R)
+    assert c.h.shape == (n * R,)
+
+
+def test_symmetry():
+    c, _, _ = _make()
+    np.testing.assert_array_almost_equal(c.J, c.J.T)
+
+
+def test_zero_diagonal():
+    c, n, R = _make()
+    np.testing.assert_array_equal(np.diag(c.J), np.zeros(n * R))
+
+
+def test_intra_replica_scaling():
+    n, R, beta = 2, 3, 2.0
+    J_q = np.array([[0.0, 1.0], [1.0, 0.0]])
+    h_q = np.array([0.5, -0.5])
+    c = TransverseFieldIsing(J_q, h_q, gamma=0.0, beta=beta, n_replicas=R)
+    scale = beta / R
+    for tau in range(R):
+        s = tau * n
+        np.testing.assert_array_almost_equal(c.J[s:s + n, s:s + n], scale * J_q)
+        np.testing.assert_array_almost_equal(c.h[s:s + n], scale * h_q)
+
+
+def test_no_inter_replica_when_gamma_zero():
+    n, R = 2, 3
+    J_q = np.array([[0.0, 1.0], [1.0, 0.0]])
+    h_q = np.zeros(n)
+    c = TransverseFieldIsing(J_q, h_q, gamma=0.0, n_replicas=R)
+    for tau in range(R):
+        tau_next = (tau + 1) % R
+        assert c.J[tau * n, tau_next * n] == 0.0
+
+
+def test_inter_replica_coupling_value():
+    n, R = 2, 4
+    beta, gamma = 1.0, 0.5
+    J_q = np.zeros((n, n))
+    h_q = np.zeros(n)
+    c = TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=R)
+    K = -0.5 * np.log(np.tanh(beta * gamma / R))
+    # qubit 0: replica 0 connected to replica 1 at index [0, n]
+    assert c.J[0, n] == pytest.approx(K)
+    assert c.J[n, 0] == pytest.approx(K)
+
+
+def test_r2_no_double_counting():
+    """For R=2 the ring has one unique temporal bond per qubit — not two."""
+    n, R = 2, 2
+    beta, gamma = 1.0, 0.5
+    J_q = np.zeros((n, n))
+    h_q = np.zeros(n)
+    c = TransverseFieldIsing(J_q, h_q, gamma=gamma, beta=beta, n_replicas=R)
+    K = -0.5 * np.log(np.tanh(beta * gamma / R))
+    assert c.J[0, n] == pytest.approx(K)   # not 2*K
+
+
+def test_solve():
+    c, _, _ = _make(n=2, R=4)
+    solver = CaSuDaSolver(Nt=100, dt=0.1667, i0=0.8, seed=42)
+    _, all_m, _ = solver.solve(c)
+    assert all_m.shape[-1] == c.n_pbits