From 236b310f69f721454088f07f119d981ad15a62ed Mon Sep 17 00:00:00 2001 From: Gregoire CATTAN Date: Fri, 22 May 2026 16:38:36 +0200 Subject: [PATCH 1/5] add polynomial optimizer --- examples/poly.py | 26 +++++++ examples/tsp.py | 2 +- p_kit/library/poly.py | 130 ++++++++++++++++++++++++++++++++ p_kit/visualization/histplot.py | 4 +- p_kit/visualization/utils.py | 2 +- 5 files changed, 160 insertions(+), 4 deletions(-) create mode 100644 examples/poly.py create mode 100644 p_kit/library/poly.py diff --git a/examples/poly.py b/examples/poly.py new file mode 100644 index 0000000..dec0d51 --- /dev/null +++ b/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']) diff --git a/examples/tsp.py b/examples/tsp.py index 4f05266..4d57265 100644 --- a/examples/tsp.py +++ b/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 diff --git a/p_kit/library/poly.py b/p_kit/library/poly.py new file mode 100644 index 0000000..25732ee --- /dev/null +++ b/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 Date: Fri, 22 May 2026 16:54:33 +0200 Subject: [PATCH 2/5] add docplex optimizer --- examples/docplex_poly.py | 18 ++++++++ p_kit/library/docplex_optimizer.py | 72 ++++++++++++++++++++++++++++++ setup.py | 1 + 3 files changed, 91 insertions(+) create mode 100644 examples/docplex_poly.py create mode 100644 p_kit/library/docplex_optimizer.py diff --git a/examples/docplex_poly.py b/examples/docplex_poly.py new file mode 100644 index 0000000..cf66738 --- /dev/null +++ b/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']) diff --git a/p_kit/library/docplex_optimizer.py b/p_kit/library/docplex_optimizer.py new file mode 100644 index 0000000..c951bda --- /dev/null +++ b/p_kit/library/docplex_optimizer.py @@ -0,0 +1,72 @@ +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) + for var, coeff in linear_part.iter_terms(): + coeffs[(var.name,)] = coeffs.get((var.name,), 0) + coeff + if linear_part.constant: + coeffs[()] = linear_part.constant + + 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 diff --git a/setup.py b/setup.py index 2d0cd22..273f39c 100644 --- a/setup.py +++ b/setup.py @@ -45,6 +45,7 @@ extras_require={ 'tests': ['pytest', 'seaborn', 'flake8'], 'gpu': ['cupy-cuda13x'], + 'docplex': ['docplex'], }, zip_safe=False, ) From 79c28bba6e481ad403f31d1e7fc06150c0684d63 Mon Sep 17 00:00:00 2001 From: Gregoire CATTAN Date: Fri, 22 May 2026 16:59:29 +0200 Subject: [PATCH 3/5] push test --- tests/test_poly.py | 165 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 165 insertions(+) create mode 100644 tests/test_poly.py diff --git a/tests/test_poly.py b/tests/test_poly.py new file mode 100644 index 0000000..d5f4942 --- /dev/null +++ b/tests/test_poly.py @@ -0,0 +1,165 @@ +import numpy as np +import pytest +from functools import partial + +from tests.conftest import requires_module +from p_kit.library.poly import PolyOptimizer + +requires_docplex = partial(requires_module, name="docplex") + + +# --- PolyOptimizer --- + +def test_n_pbits(): + circuit = PolyOptimizer({('x',): 1}, ['x'], n_bits=4) + assert circuit.n_pbits == 4 + + +def test_n_pbits_multivar(): + circuit = PolyOptimizer({('x', 'y'): 1}, ['x', 'y'], n_bits=3) + assert circuit.n_pbits == 6 + + +def test_j_symmetric(): + circuit = PolyOptimizer({('x', 'x'): 3}, ['x'], n_bits=4, minimize=False) + assert np.allclose(circuit.J, circuit.J.T) + + +def test_j_diagonal_zero(): + circuit = PolyOptimizer({('x', 'x'): 3}, ['x'], n_bits=4, minimize=False) + assert np.allclose(np.diag(circuit.J), 0) + + +def test_linear_h(): + # 2*x with n_bits=1: h[0] = 2/2 * 2^0 = 1 + circuit = PolyOptimizer({('x',): 2}, ['x'], n_bits=1, minimize=False) + assert np.isclose(circuit.h[0], 1.0) + + +def test_quadratic_h(): + # 3*x^2 with n_bits=2, C=1.5: h[k] = 3 * 1.5 * 2^k + circuit = PolyOptimizer({('x', 'x'): 3}, ['x'], n_bits=2, minimize=False) + assert np.isclose(circuit.h[0], 4.5) + assert np.isclose(circuit.h[1], 9.0) + + +def test_quadratic_J(): + # 3*x^2 with n_bits=2: J[0,1] = 3/4 * 2^0 * 2^1 = 1.5 + circuit = PolyOptimizer({('x', 'x'): 3}, ['x'], n_bits=2, minimize=False) + assert np.isclose(circuit.J[0, 1], 1.5) + + +def test_minimize_negates(): + c_min = PolyOptimizer({('x',): 2, ('x', 'x'): 3}, ['x'], n_bits=2, minimize=True) + c_max = PolyOptimizer({('x',): 2, ('x', 'x'): 3}, ['x'], n_bits=2, minimize=False) + assert np.allclose(c_min.h, -c_max.h) + assert np.allclose(c_min.J, -c_max.J) + + +def test_monomial_order_invariant(): + c1 = PolyOptimizer({('x', 'y'): 1}, ['x', 'y'], n_bits=2, minimize=False) + c2 = PolyOptimizer({('y', 'x'): 1}, ['x', 'y'], n_bits=2, minimize=False) + assert np.allclose(c1.J, c2.J) + assert np.allclose(c1.h, c2.h) + + +def test_degree_3_raises(): + with pytest.raises(ValueError): + PolyOptimizer({('x', 'x', 'x'): 1}, ['x'], n_bits=4) + + +def test_decode_zero(): + circuit = PolyOptimizer({('x',): 1}, ['x'], n_bits=4) + assert circuit.decode(np.array([-1, -1, -1, -1]))['x'] == 0 + + +def test_decode_one(): + circuit = PolyOptimizer({('x',): 1}, ['x'], n_bits=4) + assert circuit.decode(np.array([1, -1, -1, -1]))['x'] == 1 + + +def test_decode_eight(): + circuit = PolyOptimizer({('x',): 1}, ['x'], n_bits=4) + assert circuit.decode(np.array([-1, -1, -1, 1]))['x'] == 8 + + +def test_decode_multivar(): + circuit = PolyOptimizer({('x', 'y'): 1}, ['x', 'y'], n_bits=2) + # x bits: [1,-1] → x=1; y bits: [-1,1] → y=2 + result = circuit.decode(np.array([1, -1, -1, 1])) + assert result['x'] == 1 + assert result['y'] == 2 + + +# --- DocplexOptimizer --- + +@requires_docplex +def test_docplex_linear(): + from docplex.mp.model import Model + from p_kit.library.docplex_optimizer import DocplexOptimizer + mdl = Model() + x = mdl.integer_var(lb=0, ub=3, name='x') + mdl.minimize(2 * x) + circuit = DocplexOptimizer(mdl) + expected = PolyOptimizer({('x',): 2}, ['x'], n_bits=2, minimize=True) + assert np.allclose(circuit.h, expected.h) + assert np.allclose(circuit.J, expected.J) + + +@requires_docplex +def test_docplex_quadratic(): + from docplex.mp.model import Model + from p_kit.library.docplex_optimizer import DocplexOptimizer + mdl = Model() + x = mdl.integer_var(lb=0, ub=3, name='x') + mdl.minimize(3 * x * x + 2 * x) + circuit = DocplexOptimizer(mdl) + expected = PolyOptimizer({('x', 'x'): 3, ('x',): 2}, ['x'], n_bits=2, minimize=True) + assert np.allclose(circuit.h, expected.h) + assert np.allclose(circuit.J, expected.J) + + +@requires_docplex +def test_docplex_infer_n_bits_integer(): + from docplex.mp.model import Model + from p_kit.library.docplex_optimizer import DocplexOptimizer + mdl = Model() + x = mdl.integer_var(lb=0, ub=15, name='x') + mdl.minimize(x) + circuit = DocplexOptimizer(mdl) + assert circuit.n_bits == 4 + + +@requires_docplex +def test_docplex_infer_n_bits_binary(): + from docplex.mp.model import Model + from p_kit.library.docplex_optimizer import DocplexOptimizer + mdl = Model() + x = mdl.binary_var(name='x') + mdl.minimize(x) + circuit = DocplexOptimizer(mdl) + assert circuit.n_bits == 1 + assert circuit.n_pbits == 1 + + +@requires_docplex +def test_docplex_maximize(): + from docplex.mp.model import Model + from p_kit.library.docplex_optimizer import DocplexOptimizer + mdl = Model() + x = mdl.integer_var(lb=0, ub=3, name='x') + mdl.maximize(2 * x) + circuit = DocplexOptimizer(mdl) + expected = PolyOptimizer({('x',): 2}, ['x'], n_bits=2, minimize=False) + assert np.allclose(circuit.h, expected.h) + + +@requires_docplex +def test_docplex_continuous_raises(): + from docplex.mp.model import Model + from p_kit.library.docplex_optimizer import DocplexOptimizer + mdl = Model() + x = mdl.continuous_var(lb=0, ub=1, name='x') + mdl.minimize(x) + with pytest.raises(ValueError, match="continuous"): + DocplexOptimizer(mdl) From ba78c73bfe67f208ae170a2355a83b55c3c77918 Mon Sep 17 00:00:00 2001 From: Gregoire CATTAN Date: Fri, 22 May 2026 18:19:16 +0200 Subject: [PATCH 4/5] update tests --- p_kit/library/docplex_optimizer.py | 12 ++++++++---- tests/test_poly.py | 2 +- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/p_kit/library/docplex_optimizer.py b/p_kit/library/docplex_optimizer.py index c951bda..54b407d 100644 --- a/p_kit/library/docplex_optimizer.py +++ b/p_kit/library/docplex_optimizer.py @@ -37,10 +37,14 @@ def _parse(model): # linear_part exists on QuadExpr; fall back to obj itself for LinearExpr linear_part = getattr(obj, 'linear_part', obj) - for var, coeff in linear_part.iter_terms(): - coeffs[(var.name,)] = coeffs.get((var.name,), 0) + coeff - if linear_part.constant: - coeffs[()] = linear_part.constant + 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(): diff --git a/tests/test_poly.py b/tests/test_poly.py index d5f4942..55f3169 100644 --- a/tests/test_poly.py +++ b/tests/test_poly.py @@ -2,7 +2,7 @@ import pytest from functools import partial -from tests.conftest import requires_module +from conftest import requires_module from p_kit.library.poly import PolyOptimizer requires_docplex = partial(requires_module, name="docplex") From 65fc1300adb15563655164fc2b909ee80a1e39a7 Mon Sep 17 00:00:00 2001 From: gcattan Date: Fri, 22 May 2026 18:27:11 +0200 Subject: [PATCH 5/5] Update run_examples.yml --- .github/workflows/run_examples.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/run_examples.yml b/.github/workflows/run_examples.yml index 1766fcd..0ebaf01 100644 --- a/.github/workflows/run_examples.yml +++ b/.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