Source code for eqc_models.base.polynomial

# (C) Quantum Computing Inc., 2024.
from typing import Tuple, Union, List
import numpy as np
from eqc_models.base.base import EqcModel
from eqc_models.base.operators import Polynomial, QUBO, OperatorNotAvailableError
from eqc_models.base.constraints import ConstraintsMixIn
class PolynomialMixin:
    """This class provides an instance method and property that
    manage polynomial models.
    """
    @property
    def H(self) -> Tuple[np.ndarray, np.ndarray]:
        """ 
        Hamiltonian specified as a polynomial : coefficients, indices 
        
        indices are of the format [0, idx-1, ..., idx-d] which must be non-decreasing
        and each idx-j is a 1-based index of the variable which is a power in the 
        term. For a polynomial where the highest degree is 3 and specifying a term 
        such as x_1x_2, the index array is [0, 1, 2]. Another example, x_1^2x_2 is
        [1, 1, 2].
        """
        return self.coefficients, self.indices
    
    @H.setter
    def H(self, value : Tuple[np.ndarray, np.ndarray]):
        """ Set H directly as coefficients, indices """
        coefficients, indices = value
        self.coefficients = coefficients
        self.indices = indices
    @property
    def sparse(self) -> Tuple[np.ndarray, np.ndarray]:
        return self.H
    
    def evaluate(self, solution : np.ndarray) -> float:
        """ 
        Evaluate polynomial at solution 
        
        :solution: 1-d numpy array with the same length as the number of variables
        returns a floating point value
        """
        value = self.polynomial.evaluate(np.array(solution))
        return value
    
    @property
    def dynamic_range(self) -> float:
        """
        Dynamic range is a measure in decibels of the ratio of the largest
        magnitude coefficient in a problem to the smallest non-zero magnitude
        coefficient. 
        The possible range of values are all greater than or equal to 0. The
        calculation is performed by finding the lowest non-zero of the 
        absolute value of all the coefficients, which could be empty. In that
        case, the dynamic range is undefined, so an exception is raised. If
        it is positive, then the maximum of the absolute values is divided
        by the lowest. The base-10 logarithm of that value is taken and mul-
        tiplied by 10. This is the dynamic range.
        Returns
        ----------
        
        float
        
        """
        H = self.H
        coefficients = np.array(H[0])
        try:
            lowest = np.min(np.abs(coefficients[coefficients!=0]))
        except IndexError:
            raise ValueError("Dynamic range of a Hamiltonian of all 0 is undefined")
        highest  = np.max(np.abs(coefficients))
        return 10*np.log10(highest / lowest)
[docs]
class PolynomialModel(PolynomialMixin, EqcModel):
    """
    Polynomial model base class.
    Parameters
    ------------
    coefficients: An array of polynomial coeffients.
    indices: An array of polynomial indices.
    Examples
    ------------
    >>> coeffs = np.array([1, 2, 3])
    >>> indices = np.array([[0, 0, 1], [0, 1, 1], [1, 1, 1]])
    >>> from eqc_models.base.polynomial import PolynomialModel        
    >>> polynomial = PolynomialModel(coeffs, indices)
    >>> solution = np.array([1, 1, 1])
    >>> value = polynomial.evaluate(solution)
    >>> int(value)
    6
    >>> polynomial.H # doctest: +ELLIPSIS, +NORMALIZE_WHITESPACE
    (array([1, 2, 3]), array([[0, 0, 1],
       [0, 1, 1],
       [1, 1, 1]]))    
    """
    def __init__(self, coefficients : Union[List, np.ndarray], indices : Union[List, np.ndarray]) -> None:
        self.coefficients = coefficients
        self.indices = indices
    @property
    def polynomial(self) -> Polynomial:
        coefficients, indices = self.H
        return Polynomial(coefficients=coefficients, indices=indices)
    @property
    def qubo(self) -> QUBO:
        try:
            if np.all([len(self.polynomial.indices[i]) == 2 for i in range(len(self.polynomial.indices))]):
                bin_n = 0
                bits = []
                C, J = self._quadratic_polynomial_to_qubo_coefficients(self.polynomial.coefficients,
                                                                       self.polynomial.indices, self.n)
                # upper_bound is an array of the maximum values each variable can take
                upper_bound = self.upper_bound
                if np.sum(upper_bound) != upper_bound.shape[0]:
                    for i in range(upper_bound.shape[0]):
                        bits.append(1 + np.floor(np.log2(upper_bound[i])))
                        bin_n += bits[-1]
                    bin_n = int(bin_n)
                    Q = np.zeros((bin_n, bin_n), dtype=np.float32)
                    powers = [2 ** np.arange(bit_count) for bit_count in bits]
                    blocks = []
                    linear_blocks = []
                    for i in range(len(powers)):
                        # add the linear terms to the diagonal
                        linear_blocks.append(C[i] * powers[i])
                        row = []
                        for j in range(len(powers)):
                            mult = J[i, j]
                            block = np.outer(powers[i], powers[j])
                            block *= mult
                            row.append(block)
                        blocks.append(row)
                    Q[:, :] = np.block(blocks)
                    linear_operator = np.hstack(linear_blocks)
                    Q += np.diag(linear_operator)
                else:
                    # in this case, the fomulation already has only binary variables
                    Q = np.zeros_like(J)
                    Q[:, :] = J
                    Q += np.diag(np.squeeze(C))
                return QUBO(Q)
            else:
                raise OperatorNotAvailableError("QUBO operator not available")
        except OperatorNotAvailableError as e:
            print(e)
    def _quadratic_polynomial_to_qubo_coefficients(self, coefficients, indices, num_variables):
        """
        Transform polynomial into linear and quadratic qubo coefficient arrays.
        Parameters
        ----------
        coefficients : List[float]
            Coefficients of the polynomial terms, sorted according to the assumed format.
        indices : List[Tuple[int, int]]
            Sorted list of variable indices for the polynomial terms.
        num_variables : int
            Total number of variables in the polynomial.
        Returns
        -------
        linear_coefficients : np.ndarray
            1D array of linear coefficients.
        quadratic_coefficients : np.ndarray
            2D array of quadratic coefficients.
        """
        # Initialize arrays
        linear_coefficients = np.zeros(num_variables, dtype=np.float32)
        quadratic_coefficients = np.zeros((num_variables, num_variables), dtype=np.float32)
        # Populate arrays
        for coeff, index in zip(coefficients, indices):
            if index[0] == 0:  # Linear term
                linear_coefficients[index[1] - 1] += coeff
            else:  # Quadratic term
                i, j = index
                if i == j:
                    quadratic_coefficients[i - 1, j - 1] += coeff
                elif i != j:  # Symmetric terms
                    quadratic_coefficients[i - 1, j - 1] += coeff / 2
                    quadratic_coefficients[j - 1, i - 1] += coeff / 2
        return linear_coefficients, quadratic_coefficients
[docs]
class ConstrainedPolynomialModel(ConstraintsMixIn, PolynomialModel):
    """
    Constrained Polynomial model base class.
    Parameters
    ------------
    coefficients: An array of polynomial coeffients.
    indices: An array of polynomial indices.
    lhs: Left hand side of the linear constraints.
    rhs: Right hand side of the linear constraints.
    """
    def __init__(self, coefficients : Union[List, np.ndarray], indices : Union[List, np.ndarray],
                 lhs : np.ndarray, rhs: np.ndarray):
        self.coefficients = np.array(coefficients)
        self.indices = np.array(indices).astype(np.int64)
        self.max_order = self.indices.shape[1]
        self.lhs = lhs
        self.rhs = rhs
    @property
    def penalties(self):
        """ 
        Penalty terms specified as a polynomial: coefficients, indices 
        
        indices are of the format [0, idx-1, ..., idx-d] which must be non-decreasing
        and each idx-j is a 1-based index of the variable which is a power in the 
        term. For a polynomial where the highest degree is 3 and specifying a term 
        such as x_1x_2, the index array is [0, 1, 2]. Another example, x_1^2x_2 is
        [1, 1, 2].
        Only linear equality constraints are supported. Translate Ax=b into 
        penalties using the superclass.
        """
        indices = []
        coefficients = []
        def lpad(index):
            missing = self.max_order - len(index)
            if missing > 0:
                index = (0,) * missing + index
            assert len(index) > 0
            return np.array(index)
        Pl, Pq = super(ConstrainedPolynomialModel, self).penalties
        for i in range(Pl.shape[0]):
            if Pl[i] != 0:
                indices.append(lpad((0, i+1)))
                coefficients.append(Pl[i])
            for j in range(i, Pq.shape[1]):
                if Pq[i, j] != 0:
                    indices.append(lpad((i+1, j+1)))
                    value = Pq[i, j]
                    if i!=j:
                        value += Pq[j, i]
                    coefficients.append(value)
        return coefficients, indices
[docs]
    def evaluatePenalties(self, solution : np.ndarray, include_offset=False) -> float:
        """
        Take the polynomial form of the penalties from the penalties property
        and evaluate the solution. The offset can be included by passing a True
        value to the `include_offset` keyword argument.
        Parameters
        -----------
        solution : np.ndarray
            Solution to evaluate for a penalty value
        include_offset : bool
            Optional argument indicating whether or not to include the offset value.
        Returns
        ---------
        Penalty value : float
        Examples 
        ---------
        >>> coeff = np.array([-1.0, -1.0])
        >>> indices = np.array([(0, 1), (0, 2)])
        >>> lhs = np.array([[1.0, 1.0]])
        >>> rhs = np.array([1.0])
        >>> model = ConstrainedPolynomialModel(coeff, indices, lhs, rhs)
        >>> sol = np.array([1.0, 1.0])
        >>> lhs@sol - rhs
        array([1.])
        >>> model.evaluatePenalties(sol)+model.offset
        1.0
        >>> model.evaluatePenalties(sol)
        0.0
        >>> model.evaluatePenalties(sol, include_offset=True)
        1.0
        
        """
# get the coefficients and indices for the penalty polynomial
# use the Polynomial operator to evaluate the solution
        coefficients, indices = self.penalties
        solution = np.array(solution, dtype=np.float64)
        polynomial = Polynomial(coefficients, indices)
        if include_offset:
            value = self.offset
        else:
            value = 0
        value += polynomial.evaluate(solution)
        return value
[docs]
    def evaluateObjective(self, solution : np.ndarray) -> float:
        """
        Take the polynomial coeff and indices from constructor and evalute the 
        solution with it.
        Parameters
        -----------
        solution : np.ndarray
            Soluttion to evaluate the objective value
        Returns
        --------
        objective value : float
        """
        coefficients = self.coefficients
        indices = self.indices
        solution = np.array(solution, dtype=np.float64)
        polynomial = Polynomial(coefficients, indices)
        return polynomial.evaluate(solution)
    @property
    def H(self) -> Tuple[np.ndarray, np.ndarray]:
        """ Provide the sparse format for the Hamiltonian """
        p_coeff, p_indices = self.penalties
        coefficients, indices = self.coefficients, self.indices
        terms = {}
        alpha = self.alpha
        for index, coeff in zip(p_indices, p_coeff):
            index = tuple(index)
            assert len(index) > 1
            if index not in terms:
                terms[index] = alpha * coeff
            else:
                terms[index] += alpha * coeff
        for index, coeff in zip(indices, coefficients):
            index = tuple(index)
            if index not in terms:
                terms[index] = coeff
            else:
                terms[index] += coeff
        indices = [index for index in terms.keys()]
        indices.sort()
        coefficients = [terms[tuple(index)] for index in indices]
        return coefficients, indices