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 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/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/docplex_optimizer.py b/p_kit/library/docplex_optimizer.py new file mode 100644 index 0000000..54b407d --- /dev/null +++ b/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 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