Source code for eqc_models.algorithms.penaltymultiplier

# (C) Quantum Computing Inc., 2024.
from typing import Union
import logging
import numpy as np
from eqc_models.base.constraints import ConstraintModel
from eqc_models.base.base import ModelSolver
from eqc_models.algorithms.base import Algorithm
log = logging.getLogger(name=__name__)
[docs]
class PenaltyMultiplierAlgorithm(Algorithm):
    """
    Parameters
    ----------
    model : ConstraintModel
        Instance of a model to search out a penalty multiplier value, must be constrained model.
    solver : ModelSolver subclass
        Instance of a solver class to use to run the algorithm.
    Properties
    ----------
    upper_bound : float
        Upper bound value for the objective function, this need not be a least upper bound, 
        but the tighter the value, the more efficient the search
    solutions : List
        The solutions found during the algorithm run
    alphas : List
        The values of multiplier found at each algorithm iteration
    penalties : List
        The values for penalties found at each algorithm iteration. A penalty of 0
        indicates algorithm termination.
    penalty_threshold : float
        This value is the cutoff for penalty values that are threated as 0. Default
        is 1e-6. 
    progress_threshold : float
        This value is the cutoff for checking for progress on reducing the penalty
        value. The current penalty is compared to the average of the previous two
        and if the absolute difference is less than this value, then the algorithm
        stops, reporting progress has not been made. The default value is 1e-4.
    dynamic_range : List
        The values for the dynamic range of the unconstrained problem formulation,
        which is useful for identifying difficulty in representation of the problem
        on the analog device.
    The penalty multiplier search algorithm uses an infeasible solution to select the next 
    value for the penalty multiplier. The algorithm depends upon good solutions and only
    guarantees termination when the solution found for a given multiplier is optimal. For
    this reason, the implementation will terminate when no progress is made, thus making
    it a heuristic. Providing an exact solver for the solver instance will guarantee that 
    the algorithm is correct and the penalty mulktiplier found is the minimal multiplier
    capable of enforcing the condition that an unconstrained objective value for a feasible
    solution is less than an unconstrained objective value for an infeasible solution.
    This example uses the quadratic assignment problem and the known multiplier to test 
    the implementation of the algorithm.
    >>> from eqc_models.solvers.qciclient import Dirac3IntegerCloudSolver
    >>> from eqc_models.assignment.qap import QAPModel
    >>> A = np.array([[0, 5, 8, 0, 1],
    ...               [0, 0, 0, 10, 15],
    ...               [0, 0, 0, 13, 18],
    ...               [0, 0, 0, 0, 0.],
    ...               [0, 0, 0, 1, 0.]])
    >>> B = np.array([[0, 8.54, 6.4, 10, 8.94],
    ...               [8.54, 0, 4.47, 5.39, 6.49],
    ...               [6.4, 4.47, 0, 3.61, 3.0],
    ...               [10, 5.39, 3.61, 0, 2.0],
    ...               [8.94, 6.49, 3.0, 2.0, 0.]])
    >>> C = np.array([[2, 3, 6, 3, 7],
    ...               [3, 9, 2, 5, 9],
    ...               [2, 6, 4, 1, 2],
    ...               [7, 5, 8, 5, 7],
    ...               [1, 9, 2, 9, 2.]])
    >>> model = QAPModel(A, B, C)
    >>> solver = Dirac3IntegerCloudSolver() # must be configured with environment variables
    >>> algo = PenaltyMultiplierAlgorithm(model, solver)
    >>> algo.upper_bound = 330.64
    >>> algo.run(relaxation_schedule=2, num_samples=5) # doctest: +ELLIPSIS
    2... RUNNING... COMPLETED...
    >>> algo.alphas[-1] # doctest: +SKIP
    106.25
    >>> algo.penalties[-1] # doctest: +SKIP
    0.0
 
    """
    def __init__(self, model : ConstraintModel, solver : ModelSolver, penalty_threshold:float=1e-6,
                 progress_threshold : float = 1e-4):
        self.model = model
        self.solver = solver
        # ub = np.max(model.quad_objective)
        # if ub < np.max(model.linear_objective):
        #     ub = np.max(model.linear_objective)
        #     ub *= model.sum_constraint
        # else:
        #     ub *= model.sum_constraint ** 2
        self.ub = None # ub
        self.solutions = None
        self.penalties = None
        self.alphas = None
        self.dynamic_range = None
        self.responses = None
        self.penalty_threshold = penalty_threshold
        self.progress_threshold = progress_threshold
    @property
    def upper_bound(self) -> float:
        return self.ub
    @upper_bound.setter
    def upper_bound(self, value : float):
        self.ub = value
[docs]
    def run(self, initial_alpha : float=None, initial_solution : np.array = None, **kwargs):
        """ Start with a guess at alpha, iterate until alpha is sufficiently large """
        self.solutions = solutions = []
        self.penalties = penalties = []
        self.alphas = alphas = []
        self.dynamic_range = dynamic_range = []
        self.responses = responses = []
        self.energies = energies = []
        model = self.model
        solver = self.solver
        offset = model.offset
        ub = self.ub
        if initial_alpha is None and offset > 0:
            alpha = ub / offset
            log.info("UPPER BOUND %f OFFSET %f  -> ALPHA %f", 
                     ub, offset, alpha)
            if alpha < 1:
                alpha = 1
        elif initial_alpha is not None:
            alpha = initial_alpha
        else:
            log.info("No tricks for initial alpha, setting to 1")
            alpha = 1
        if initial_solution is not None:
            log.debug("INITIAL SOLUTION GIVEN")
            obj_val = model.evaluate(initial_solution, alpha, True)
            penalty = model.evaluatePenalties(initial_solution) + offset
            log.info("INITIAL SOLUTION OBJECTIVE %f PENALTY %f", obj_val, penalty)
            if obj_val < ub:
                alpha += (ub - obj_val) / penalty
            log.info("INITIAL SOLUTION DETERMINED ALPHA %f", alpha)
        else:
            penalty = None
            
        while penalty is None or penalty > self.penalty_threshold:
            log.info("NEW RUN")
            log.info("SETTING MULTIPLIER %f", alpha)
            model.penalty_multiplier = float(alpha)
            log.info("GOT MULTIPLIER %f NEW OFFSET %f", model.penalty_multiplier, 
                     model.penalty_multiplier * model.offset)
            dynamic_range.append(float(model.dynamic_range))
            log.info("CALLING SOLVE WITH ALPHA %f DYNAMIC RANGE %f", alpha, dynamic_range[-1])
            alphas.append(float(alpha))
            response = solver.solve(model, **kwargs)
            responses.append(response)
            results = response["results"]
            solution = np.array(results["solutions"][0])
            solutions.append(solution)
            penalty = model.evaluatePenalties(solution) + offset
            penalties.append(float(penalty))
            obj_val = model.evaluate(solution, alpha, True)
            less_offset = model.evaluate(solution, alpha, False)
            energies.append(results["energies"][0])
            log.info("NEW SOLUTION OBJECTIVE %f LESS OFFSET %f ENERGY %f PENALTY %f", 
                     obj_val, less_offset, energies[-1], penalty)
            if penalty < self.penalty_threshold:
                pass
            elif obj_val < ub:
                alpha += (ub - obj_val) / penalty
            if penalty > self.penalty_threshold and abs(sum(penalties[-2:])/2-penalty) < self.progress_threshold:
                log.warn("SUFFICIENT PROGRESS NOT MADE FOR THREE ITERATIONS, QUITTING")
                break