Source code for eqc_models.base.operators

"""
QUBO and Polynomial operators are used in EQC Models to pass the appropriate 
format to the solver. All models must output one or both of QUBO or Polynomial
types. Each Solver checks for the appropriate type. If a model does not provide
the type, then it must raise OperatorNotAvailableError.

>>> Q = np.array([[1, 2], [0, 1]])
>>> qubo = QUBO(Q)
>>> (qubo.Q == np.array([[1, 1], [1, 1]])).all()
True
>>> coefficients = [1, 1, 2]
>>> indices = [(1, 1), (2, 2), (1, 2)]
>>> poly = Polynomial(coefficients, indices)
>>> indices = [(1, 1), (2, 2), (2, 1)]
>>> poly = Polynomial(coefficients, indices)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: Input data to polynomial is not correct: index order must be monotonic
>>> s = np.array([1, 1])
>>> (poly.evaluate(s)==qubo.evaluate(s)).all()
True
"""
from typing import List
from dataclasses import dataclass
import numpy as np
# polyeval module is a Cython module useful for speeding up computation 
# of an operator's value at a solution. If the Cython module is not 
# available, the poly_eval variable will be set to None, triggering the
# use of a pure Python method for evaluation
try:
    from .polyeval import poly_eval
except ImportError:
    poly_eval = None

class OperatorNotAvailableError(Exception):
    """ General error to raise when an operator type is not implemented """

[docs] @dataclass class QUBO: """ Contains a QUBO in a symmetric matrix If the matrix Q is not symmetric already, it will be after __post_init__ Parameters ---------- qubo : np.array 2d symmetric matrix of values which describe a quadratic unconstrained binary optimization problem. """ Q: np.ndarray def __post_init__(self): Q2 = (self.Q + self.Q.T) / 2 self.Q = Q2
[docs] def evaluate(self, solution : np.ndarray): return solution.T@self.Q@solution
[docs] @dataclass class Polynomial: """ Represents an operator and evalute the operator at a point or set of points. The operator must be a polynomial, possibly multivariate, with powers of up to 5, at the time of this version. The representation of a polynomial uses a sparse method with two components per term. A term is described by the coefficient and a tuple of integers which indicate the variable indexes of the term. The coefficients can be any value which fits in 32-bit floating point representation, but the dynamic range of the coefficients should be within the limit of the hardware's sensitivity for best results. The term index tuple length must be consistent across all terms. If a term does not have a variable index for all positions, such as with a term which is the square of a variable when other terms have third-order powers, then there must be a placeholder of 0 for the unused power. The variable indexes must be in the tuple in ascending order. Here are some examples (suppose the max degree is 4): - :math:`x_1^2`: :code:`(0, 0, 1, 1)` - :math:`x_1 x_2 x_3`: :code:`(0, 1, 2, 3)` - :math:`x_2^2 x_3^2`: :code:`(2, 2, 3, 3)` while it does not affect the optimiztion, a constant term can be applied to the polynomial by using an index of all zeros :code:`(0, 0, 0, 0)`. When listing the coefficients, the position in the array must correspond to the position in the array of indexes. Also, the indices must be ordered with linear terms first, quadratic terms next and so forth. A polynomial operator does not have an explicit domain. It could be evaluated on an array of any real numbers. Parameters ---------- coefficients : List, np.array Floating point values for the coefficients of a polynomial. Must correspond to the entries in the indices array. indices : List, np.array List of tuples or 2d np.array with integer values describing the term which the corresponding coefficient value is used for. Examples -------- >>> coefficients = [-1, -1, 2] >>> indices = [(0, 1), (0, 2), (1, 2)] >>> polynomial = Polynomial(coefficients, indices) >>> test_solution = np.array([1, 1]) >>> polynomial.evaluate(test_solution) [0] >>> test_solution = np.array([1, 0]) >>> polynomial.evaluate(test_solution) [-1] >>> test_solution = np.array([5, -1]) >>> polynomial.evaluate(test_solution) [-14] >>> test_solution = np.array([2.5, -2.5]) >>> polynomial.evaluate(test_solution) [-12.5] """ coefficients : List indices : List def __post_init__(self): issues = set() degree_count = None if len(self.coefficients)!=len(self.indices): issues.add("coefficients and indices must be the same length") # ensure indices are not numpy self.indices = [[int(val) for val in index] for index in self.indices] for i in range(len(self.coefficients)): if degree_count is None: degree_count = len(self.indices[i]) elif len(self.indices[i])!=degree_count: issues.add("term rank is not consistent") for j in range(1, degree_count): if self.indices[i][j] < self.indices[i][j-1]: issues.add("index order must be monotonic") try: coeff = float(self.coefficients[i]) except TypeError: issues.add("coefficient data types must be coercible to float type") except ValueError: issues.add("coefficient values must be coercible to float values") if len(issues) > 0: msg = "Input data to polynomial is not correct: " msg += "; ".join([issue for issue in issues]) raise ValueError(msg)
[docs] def pure_evaluate(self, solution : np.ndarray) -> np.ndarray: """ Evaluation in pure python Parameters ---------- solution : np.array Solution to evaluate, is optionally 2-d, which results in multiple evaluations Returns ------- np.array """ if len(solution.shape) == 1: solution = solution.reshape((1, solution.shape[0])) objective = [0 for k in range(solution.shape[0])] for k in range(solution.shape[0]): for i in range(len(self.coefficients)): term = self.coefficients[i] for j in self.indices[i]: if j > 0: term *= solution[k, j-1] objective[k] += term return objective
[docs] def evaluate(self, solution: np.ndarray): """ Evaluate the polynomial at the solution point. If the Cython module is available, use that for speedup, otherwise evaluate with Python loops. Parameters ---------- solution : np.array Solution to evaluate. Is optinoally 2-d, which results in multiple exaluations. Returns ------- 1-d array of values which match the coerced dtype of the inputs. """ if poly_eval is None: return self.pure_evaluate(solution) else: return poly_eval(np.array(self.coefficients, dtype=np.float64), np.array(self.indices, dtype=np.int64), solution)