Add polynomial expression optimizer

.github/workflows/run_examples.yml

@@ -22,6 +22,7 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         pip install .
+        pip install .[docplex]
     - name: Run examples
       run: |
         cd examples

examples/docplex_poly.py

@@ -0,0 +1,18 @@
+from docplex.mp.model import Model
+from p_kit.library.docplex_optimizer import DocplexOptimizer
+from p_kit.solver.csd_solver import CaSuDaSolver
+from p_kit.solver.annealing import constant, execute
+from p_kit.visualization.histplot import histplot
+
+# Build the docplex model: minimize 3*x^2 + 2*x + 5
+mdl = Model(name="poly")
+x = mdl.integer_var(lb=0, ub=15, name="x")
+mdl.minimize(3 * x * x + 2 * x + 5)
+
+# Convert to p-bit circuit (n_bits auto-inferred from ub=15 → 4 bits)
+circuit = DocplexOptimizer(mdl)
+
+solver = CaSuDaSolver(Nt=int(1e4), dt=0.1, i0=1.0, expected_mean=0, seed=42)
+samples = execute(solver, circuit, constant, n_shots=100, n_last_samples=10, n_jobs=-1)
+
+histplot(samples, decode=lambda m: circuit.decode(m)['x'])

examples/poly.py

@@ -0,0 +1,26 @@
+import numpy as np
+from p_kit.library.poly import PolyOptimizer
+from p_kit.solver.csd_solver import CaSuDaSolver
+from p_kit.solver.annealing import constant, execute
+from p_kit.visualization.histplot import histplot
+
+# Minimize 3*x^2 + 2*x + 5  over integers [0, 15] (n_bits=4)
+# Global minimum on non-negative integers: x=0  (f=5)
+coeffs = {
+    ('x', 'x'): 3,
+    ('x',): 2,
+    (): 5,
+}
+circuit = PolyOptimizer(coeffs, variables=['x'], n_bits=4, minimize=True)
+
+solver = CaSuDaSolver(Nt=int(1e4), dt=0.1, i0=1.0, expected_mean=0, seed=42)
+samples = execute(solver, circuit, constant, n_shots=100, n_last_samples=10, n_jobs=-1)
+
+def f(x):
+    return 3 * x ** 2 + 2 * x + 5
+
+results = [circuit.decode(row) for row in samples]
+best = min(results, key=lambda r: f(r['x']))
+print(f"Best x = {best['x']},  f(x) = {f(best['x'])}")
+
+histplot(samples, decode=lambda m: circuit.decode(m)['x'])

examples/tsp.py

@@ -37,7 +37,7 @@
 circuit = TSP(city_graph=city_graph, tsp_modifier=0.33)
 
 # Simulated Annealing Schedule
-solver = CaSuDaSolver(Nt=int(1e5), dt=1 / (2 * len(city_graph)), i0=0, expected_mean=0, seed=None)
+solver = CaSuDaSolver(Nt=int(1e5), dt=1 / (2 * len(city_graph)), i0=0, expected_mean=0, seed=None, device="cpu")
 
 # Increase n_shots to have more reliable results
 n_shots = 200

p_kit/library/docplex_optimizer.py

@@ -0,0 +1,76 @@
+import math
+from p_kit.library.poly import PolyOptimizer
+
+
+class DocplexOptimizer(PolyOptimizer):
+    """Build a PolyOptimizer from a DOcplex linear/quadratic model.
+
+    Only the objective function is encoded into J/h — constraints are ignored.
+    Objective must be linear or quadratic (degree ≤ 2).
+    The optimization sense (minimize/maximize) is read from the model.
+
+    Parameters
+    ----------
+    model : docplex.mp.model.Model
+    n_bits : int, optional
+        Bit-width per variable. If None, auto-inferred from variable bounds:
+        binary variables → 1, integer [0, ub] → ceil(log2(ub + 1)).
+        All variables share the same n_bits; for mixed models, pass an explicit value.
+
+    Notes
+    -----
+    .. versionadded:: 0.0.1
+    """
+
+    def __init__(self, model, n_bits=None):
+        coeffs, variables = DocplexOptimizer._parse(model)
+        minimize = model.is_minimized()
+        if n_bits is None:
+            n_bits = DocplexOptimizer._infer_n_bits(model, variables)
+        super().__init__(coeffs, variables, n_bits=n_bits, minimize=minimize)
+
+    @staticmethod
+    def _parse(model):
+        obj = model.objective_expr
+        coeffs = {}
+        variables = [v.name for v in model.iter_variables()]
+
+        # linear_part exists on QuadExpr; fall back to obj itself for LinearExpr
+        linear_part = getattr(obj, 'linear_part', obj)
+        if hasattr(linear_part, 'iter_terms'):
+            for var, coeff in linear_part.iter_terms():
+                coeffs[(var.name,)] = coeffs.get((var.name,), 0) + coeff
+        elif hasattr(linear_part, 'name'):
+            coeffs[(linear_part.name,)] = coeffs.get((linear_part.name,), 0) + 1.0
+        const = getattr(linear_part, 'constant', 0)
+        if const:
+            coeffs[()] = const
+
+        if hasattr(obj, 'iter_quad_triplets'):
+            for v1, v2, coeff in obj.iter_quad_triplets():
+                key = tuple(sorted([v1.name, v2.name]))
+                coeffs[key] = coeffs.get(key, 0) + coeff
+
+        return coeffs, variables
+
+    @staticmethod
+    def _infer_n_bits(model, variables):
+        var_map = {v.name: v for v in model.iter_variables()}
+        max_bits = 1
+        for name in variables:
+            v = var_map[name]
+            if v.is_binary():
+                bits = 1
+            elif v.is_integer():
+                ub = v.ub
+                if math.isinf(ub):
+                    raise ValueError(
+                        f"Variable '{name}' has no finite upper bound; specify n_bits explicitly."
+                    )
+                bits = max(1, math.ceil(math.log2(ub + 1)))
+            else:
+                raise ValueError(
+                    f"Variable '{name}' is continuous; only binary and integer variables are supported."
+                )
+            max_bits = max(max_bits, bits)
+        return max_bits

p_kit/library/poly.py

@@ -0,0 +1,130 @@
+import numpy as np
+from p_kit.psl.p_circuit import PCircuit
+
+
+class PolyOptimizer(PCircuit):
+    """Encode a multivariate polynomial (linear + quadratic terms only) into J/h.
+
+    Each integer variable is binary-expanded: x = sum_k(2^k * b_k), b_k in {0,1},
+    then mapped to p-bits via b_k = (s_k + 1) / 2, s_k in {-1, +1}.
+
+    Parameters
+    ----------
+    coeffs : dict
+        Maps monomial tuples of variable names to float coefficients.
+        Order within a tuple does not matter; {('x','y'): c} == {('y','x'): c}.
+        - Constant:         {(): 5}           (ignored)
+        - Linear:           {('x',): 2}
+        - Quadratic (x^2):  {('x', 'x'): 3}
+        - Cross-term (x*y): {('x', 'y'): 1}
+        Degree > 2 raises ValueError.
+    variables : list of str
+        Variable names. Defines ordering and total p-bit count.
+    n_bits : int, default 4
+        Bit-width per variable. Each variable ranges over [0, 2^n_bits - 1].
+    minimize : bool, default True
+        If True, negate J and h so the solver minimizes the polynomial.
+        Set to False to maximize instead.
+
+    Attributes
+    ----------
+    h : np.array((n_pbits,))
+        Bias vector.
+    J : np.array((n_pbits, n_pbits))
+        Coupling matrix.
+
+    Notes
+    -----
+    .. versionadded:: 0.0.1
+    """
+
+    def __init__(self, coeffs, variables, n_bits=4, minimize=True):
+        self.variables = list(variables)
+        self.n_bits = n_bits
+        self.minimize = minimize
+        self.coeffs = self._normalize(coeffs)
+        PCircuit.__init__(self, len(variables) * n_bits)
+        self._encode()
+
+    @staticmethod
+    def _normalize(coeffs):
+        result = {}
+        for mono, val in coeffs.items():
+            if len(mono) > 2:
+                raise ValueError(
+                    f"Monomial {mono} has degree {len(mono)}; only linear and quadratic terms are supported."
+                )
+            key = tuple(sorted(mono))
+            result[key] = result.get(key, 0) + val
+        return result
+
+    def _idx(self, var, bit):
+        return self.variables.index(var) * self.n_bits + bit
+
+    def _encode(self):
+        n = self.n_bits
+        # Offset introduced by the (s+1)/2 binary expansion: x = (1/2)*sum_k 2^k*s_k + C
+        C = (2 ** n - 1) / 2
+
+        for mono, coeff in self.coeffs.items():
+            if len(mono) == 0:
+                continue
+
+            elif len(mono) == 1:
+                # Linear: coeff * x_v = (coeff/2) * sum_k 2^k * s_{v,k} + const
+                v = mono[0]
+                for k in range(n):
+                    self.h[self._idx(v, k)] += coeff / 2 * (2 ** k)
+
+            elif mono[0] == mono[1]:
+                # Quadratic x_v^2:
+                #   h contribution:  coeff * C * 2^k  per bit k
+                #   J contribution (j<k):  coeff/4 * 2^j * 2^k  (factor-of-2 symmetry accounted for)
+                v = mono[0]
+                for k in range(n):
+                    self.h[self._idx(v, k)] += coeff * C * (2 ** k)
+                for j in range(n):
+                    for k in range(j + 1, n):
+                        w = coeff / 4 * (2 ** j) * (2 ** k)
+                        ij, ik = self._idx(v, j), self._idx(v, k)
+                        self.J[ij, ik] += w
+                        self.J[ik, ij] += w
+
+            else:
+                # Cross-term x_u * x_v (u != v):
+                #   h contribution:  coeff * C/2 * 2^k  per bit k, for both variables
+                #   J contribution:  coeff/8 * 2^j * 2^k  (factor-of-2 symmetry accounted for)
+                u, v = mono
+                for k in range(n):
+                    wk = 2 ** k
+                    self.h[self._idx(u, k)] += coeff * C / 2 * wk
+                    self.h[self._idx(v, k)] += coeff * C / 2 * wk
+                for j in range(n):
+                    for k in range(n):
+                        w = coeff / 8 * (2 ** j) * (2 ** k)
+                        ij, ik = self._idx(u, j), self._idx(v, k)
+                        self.J[ij, ik] += w
+                        self.J[ik, ij] += w
+
+        if self.minimize:
+            self.J = -self.J
+            self.h = -self.h
+
+        np.fill_diagonal(self.J, 0)
+
+    def decode(self, samples):
+        """Convert p-bit samples {-1, +1} to integer variable values.
+
+        Parameters
+        ----------
+        samples : array-like of shape (n_pbits,)
+
+        Returns
+        -------
+        dict mapping variable name to integer value in [0, 2^n_bits - 1]
+        """
+        bits = ((np.asarray(samples) + 1) / 2).astype(int)
+        return {
+            var: int(sum(int(bits[self._idx(var, k)]) * (2 ** k) for k in range(self.n_bits)))
+            for var in self.variables
+        }

p_kit/visualization/histplot.py

@@ -5,11 +5,11 @@
 import numpy as np
 
 
-def histplot(output):
+def histplot(output, decode=m_to_string):
     ret = {}
 
     for m in output:
-        s = m_to_string(m)
+        s = str(decode(m))
         if s in ret:
             ret[s] = ret[s] + 1
         else:

p_kit/visualization/utils.py

@@ -18,7 +18,7 @@ def extract_tsp_path(sample_matrix):
         path.append(city)
     return path
 
-def tsp_hist(samples, city_graph):  
+def tsp_hist(samples, city_graph):
     hist = {}
     for i in range(len(samples)):
         s = samples[i, :].reshape((len(city_graph), len(city_graph[0])))

setup.py

@@ -45,6 +45,7 @@
       extras_require={
                       'tests': ['pytest', 'seaborn', 'flake8'],
                       'gpu': ['cupy-cuda13x'],
+                      'docplex': ['docplex'],
                       },
       zip_safe=False,
 )