# eqc_models.base

This subpackage contains the building blocks for formulating and solving models with EQC devices. Many well known models in the space of quantum computing use quadratic operators. Two generally accepted formats for quadratic models are Ising Hamiltonian and QUBO operators. Ising models are not explicitly supported because of the simple equivalence to QUBO models. Converters from Ising to QUBO are readily available and may be added as a utility function in the future, but for now are left out of the package. Other quadratic models, such as integer or floating point models are supported using the domain and upper bound of the model variables and the solver choice.

Polynomial models are used to support higher order interactions between variables. These models are defined over non-negative variables with the highest power interaction being the only restriction on variable interactions. This is also known as all-to-all connectivity.

## Subpackages

constraints

  This module defines support for linear equality and inequality constraints. It can be paired with the penalty multiplier algorithm to choose multiplier values which satisfy the enforcement of constraints as penalties. The mixin pattern is used to allow more complex classes to be built from these bases.

  - ConstraintMixIn

  - ConstraintModel

  - InequalitiesMixIn

  - InequalityConstraintModel

quadratic

  This module defines the methods required to convert a linear vector and quadratic matrix into the required operator, either QUBO or Polynomial, for solving with EQC. The support for constraints is also incorporated.

  - ConstrainedQuadraticModel

  - QuadraticModel

polynomial

  This module defines the classes and methods for conversion from arrays of coefficients and indices to operators for solving with EQC. The support for constriants is also incorporated.

  - PolynomialModel

  - ConstrainedPolynomialModel

base

  This module defines abstract base classes for models and solvers.

  - EqcModel

  - ModelSolver

operators

  This module defines the operator types to use for passing problems to EQC devices.

  - Polynomial

  - QUBO

class eqc_models.base.ConstraintsMixIn

  Bases: `object`

  This mixin class contains methods and attributes which transform linear constraints into penalties.

  lhs = None
  rhs = None
  alpha = 1
  linear_objective = None
  quad_objective = None
  property penalties: Tuple[ndarray, ndarray]

    Returns two numpy arrays, one linear and one quadratic pieces of an operator

  property penalty_multiplier: float
  property constraints
  property offset: float

    Calculate the offset due to the conversion of constraints to penalties

  evaluate(solution: ndarray, alpha: float | None = None, includeoffset: bool = False)

    Compute the objective value plus penalties for the given solution. Including the offset will ensure the penalty contribution is non-negative.

    Parameters:

      - solution (np.array) – The solution vector for the problem.

      - alpha (float) – Penalty multiplier, optional. This can be used to test different multipliers for determination of sufficiently large values.

  evaluatePenalties(solution) → float

    Evaluate penalty function without alpha or offset

    Parameters:

      solution (np.array) – The solution vector for the problem.

  checkPenalty(solution: ndarray)

    Get the penalty of the solution.

    Parameters:

      solution (np.array) – The solution vector for the problem.

class eqc_models.base.ConstraintModel

  Bases: `ConstraintsMixIn`, `EqcModel`

  Abstract class for representing linear constrained optimization problems as EQC models.

class eqc_models.base.ConstrainedQuadraticModel(C_obj: array, J_obj: array, lhs: array, rhs: array)

  Bases: `ConstraintsMixIn`, `QuadraticModel`

  Provides a constrained quadratic operator and device sum constraint support.

  Parameters:

    - J (Quadratic hamiltonian array.)

    - C (Linear hamiltonian array.)

    - lhs (Left hand side of the linear constraints.)

    - rhs (Right hand side of the linear constraints.)

  Examples

      >>> C = np.array([1, 2])
      >>> J = np.array([[2, 1], [1, 2]])
      >>> lhs = np.array([[1, 1], [1, 1]])
      >>> rhs = np.array([1, 1])
      >>> from eqc_models.base.quadratic import ConstrainedQuadraticModel
      >>> model = ConstrainedQuadraticModel(C, J, lhs, rhs)
      >>> model.penalty_multiplier = 1
      >>> model.upper_bound = np.array([1, 1])
      >>> qubo = model.qubo
      >>> qubo.Q
      array([[1., 3.],
             [3., 2.]])

  property H: Tuple[ndarray, ndarray]

    Return a pair of arrays, the linear and quadratic portions of a quadratic operator that has the objective plus penalty functions multiplied by the penalty multiplier. The linear terms are the first array and the quadratic are the second.

    Return type:

      np.ndarray, np.nedarray

  property penalty_multiplier: float
  property constraints
  evaluateObjective(solution: ndarray) → float

class eqc_models.base.QuadraticModel(C: ndarray, J: ndarray)

  Bases: `QuadraticMixIn`, `EqcModel`

  Provides a quadratic operator and device sum constraint support.

  Parameters:

    - J (Quadratic hamiltonian array.)

    - C (Linear hamiltonian array.)

  Examples

      >>> C = np.array([1, 2])
      >>> J = np.array([[2, 1], [1, 2]])
      >>> from eqc_models.base.quadratic import QuadraticModel
      >>> model = QuadraticModel(C, J)
      >>> model.upper_bound = np.array([1, 1])
      >>> qubo = model.qubo
      >>> qubo.Q
      array([[3., 1.],
             [1., 4.]])

class eqc_models.base.PolynomialModel(coefficients: List | ndarray, indices: List | ndarray)

  Bases: `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
      (array([1, 2, 3]), array([[0, 0, 1],
         [0, 1, 1],
         [1, 1, 1]]))

  property polynomial: Polynomial
  property qubo: QUBO

class eqc_models.base.ConstrainedPolynomialModel(coefficients: List | ndarray, indices: List | ndarray, lhs: ndarray, rhs: ndarray)

  Bases: `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.)

  property penalties

    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.

    Type:

      Penalty terms specified as a polynomial

  evaluatePenalties(solution: 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

    Return type:

      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

  evaluateObjective(solution: 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

    Return type:

      float

  property H: Tuple[ndarray, ndarray]

    Provide the sparse format for the Hamiltonian

class eqc_models.base.InequalitiesMixin

  Bases: `object`

  This mixin enables inequality constraints by automatically generating slack variables for each inequality

  This mixin adds a senses attribute which has a value for each constraint. The values are one of ‘EQ’, ‘LE’ or ‘GE’ for equal to, less than or equal to or greater than or equal to. The effect of the value is to control whether a slack is added and what the sign of the slack variable in the constraint is. Negative is used for GE, positive is used for LE and all slack variables get a coefficient magnitude of 1.

  The constraints are modified on demand, so the class members, lhs and rhs remain unmodified.

  property senses: List[str]

    Comparison operator by constraint

  property num_slacks: int

    The number of slack variables. Will match the number of inequality constraints.

    Returns:

      number

    Return type:

      int

  property constraints: Tuple[ndarray, ndarray]

    Get the general form of the constraints, add slacks where needed and return a standard, equality constraint form.

class eqc_models.base.InequalityConstraintModel

  Bases: `InequalitiesMixin`, `ConstraintModel`

  Abstract class for a linear constrained optimization model with inequality constraints

class eqc_models.base.EqcModel

  Bases: `object`

  EqcModel subclasses must provide these properties/methods.

  Decode:

    takes a raw solution and translates it into the original problem formulation

  H:

    property which returns a Hamiltonian operator

  Upper_bound:

    Let D be an array of length n which contains the largest possible value allowed for x[i], which is the variable at index i, 0<=i<n. This means that a x[i] is in the domain [0,D[i]]. If the solution type of x[i] is integer, then x[i] is in the set of integers, Z, and also 0<=x[i]<=floor(D[i]).

  Qudit_limits:

    maximum value permitted for each qudit

      >>> model = EqcModel()
      >>> ub = np.array([1, 1.5, 2])
      >>> model.upper_bound = ub
      Traceback (most recent call last):
      ...
      ValueError: ...
      >>> model.upper_bound = np.ones((3,))
      >>> (model.upper_bound==np.ones((3,))).all()
      True
      >>> model.upper_bound = 2*np.ones((3,))
      >>> (model.upper_bound==2).all()
      True

  decode(solution: ndarray, from_encoding: str | None = None) → ndarray

    Manipulate the solution to match the variable count. Optionally, a log-encoded solution vector can be translated into a vector with integral values.

    Parameters:

      - solution (np.ndarray) – 1d solution vector

      - from_encoding (str) – string indicating if the vector is an encoding of particular type. The text ‘qubo’ is the only option that does anything for now.

      - EqcModel() (>>> model =)

      - np.array([1 (>>> solution =)

      - 2

      - 5])

      - ub (>>> model.upper_bound =)

      - 1 (>>> model.machine_slacks =)

      - np.array([1

      - 1

      - 1

      - 1

      - 0

      - 1

      - 3 (0]) # this solution has a)

      - bound (which violates the upper)

      - care (but decode doesn't)

      - model.decode(solution (>>>)

      - "qubo")

      - array([1

      - 3

      - 5])

  property upper_bound: array

    An array of upper bound values for every variable in the model. Must be integer.

  property domains: array
  property is_discrete: List

    An array of boolean values indicating if variable should be restricted to a discrete domain

  property n: int

    Return the number of variables

  property H

    Hamiltonian operator of unknown type

  property sparse: Tuple[ndarray, ndarray]
  property machine_slacks

    Number of slack qudits to add to the model

  evaluateObjective(solution: ndarray) → float
  createConfigElements() → Dict
  createBenchmarkConfig(fname: str) → None
  property dynamic_range: float
  property polynomial: Polynomial
  property qubo: QUBO

class eqc_models.base.ModelSolver

  Bases: `object`

  Provide a common interface for solver implementations. Store a model, implement a solve method.

  solve(model: EqcModel, *args, **kwargs) → Dict

class eqc_models.base.QUBO(Q: ndarray)

  Bases: `object`

  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: ndarray
  evaluate(solution: ndarray)

class eqc_models.base.Polynomial(coefficients: List, indices: List)

  Bases: `object`

  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):

  - \(x_1^2\): `(0, 0, 1, 1)`

  - \(x_1 x_2 x_3\): `(0, 1, 2, 3)`

  - \(x_2^2 x_3^2\): `(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 `(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
  pure_evaluate(solution: ndarray) → ndarray

    Evaluation in pure python

    Parameters:

      solution (np.array) – Solution to evaluate, is optionally 2-d, which results in multiple evaluations

    Return type:

      np.array

  evaluate(solution: 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.

    Return type:

      1-d array of values which match the coerced dtype of the inputs.

  to_json(save_file: str | TextIO, sum_constraint: float | None = None, levels: int | List | None = None) → None

    Stream this operator to a json file for running directly on Dirac-3.

    Parameters:

      - save_file:: – path (string) or IO object for streaming the json encoding.

      - sum_constraint:: – scalar value indicating the sum constraint for the problem. This indicates a continuous valued problem.

      - levels:: – Either an integer or iterable of integers indicating the number of levels of each variable.

class eqc_models.base.SolutionResults(solutions: ndarray, energies: ndarray, counts: ndarray, objectives: ndarray, run_time: ndarray, preprocessing_time: int, postprocessing_time: ndarray, penalties: ndarray | None = None, device: str | None = None, raw_solutions: ndarray | None = None, time_units: str = 'ns', response: Any | None = None)

  Bases: `object`

  The class is meant to provide a uniform interface to results, no matter the method of running the job. If available, the metrics are reported in nanoseconds.

  ## Properties

  solutionsnp.ndarray

    2-d array of solution vectors

  energiesnp.ndarray

    1-d array of energies computed from the device for each sample

  countsnp.ndarray

    1-d array of counts the particular sample occurred during sampling

  objectivesnp.ndarray

    1-d array of objective values. Is None if the model does not provide a separate objective function

  run_timenp.ndarray

    1-d array of runtimes reported by the device.

  preprocessing_timeint

    Single value for time spent preprocessing before sampling occurs.

  postprocessing_timenp.ndarray

    1-d array of time spent post-processing samples.

  penaltiesnp.ndarray

    1-d array of penalty values for each sample. Is None if the model does not have constraints.

  devicestr

    String that represents the device used to solve the model.

  raw_solutionsnp.ndarray

    Numpy array of the solutions as returned by the solver device.

  calibration_timefloat

    Total time spend during job exectution where the device performed calibration. The calibration is not directly affected by the job submission and the time is not included in run_time.

  time_unitsstr

    String indicator of the unit of time reported in the metrics. Only ns is supported at this time.

  solutions: ndarray
  energies: ndarray
  counts: ndarray
  objectives: ndarray
  run_time: ndarray
  preprocessing_time: int
  postprocessing_time: ndarray
  penalties: ndarray = None
  device: str = None
  raw_solutions: ndarray = None
  time_units: str = 'ns'
  response: Any = None
  to_dict() → Dict

    Return a dictionary of values with the dataclass attributes as keys

  property device_time: ndarray

    1-d array of device usage computed from preprocessing, runtime and postprocessing time.

  property total_samples
  property best_energy
  classmethod determine_device_type(device_config)

    Use the device config object from a cloud response to get the device info. It will have a device and job type identifiers in it.

  classmethod from_cloud_response(model, response, solver)

    Fill in the details from the cloud

    Parameters:

      - model (eqc_models.base.EqcModel) – EqcModel object describing the problem solved in response

      - response (Dict) – Dictionary of the repsonse from the solver device.

      - solver (eqc_models.base.ModelSolver) – ModelSolver object which is used to obtain job metrics.

  classmethod from_eqcdirect_response(model, response, solver)

    Fill in details from the response dictionary and possibly the solver device.

    Parameters:

      - model (eqc_models.base.EqcModel) – EqcModel object describing the problem solved in response

      - response (Dict) – Dictionary of the repsonse from the solver device.

      - solver (eqc_models.base.ModelSolver) – ModelSolver object which is used to obtain device information.

  get(key, default)

# eqc_models.solvers

class eqc_models.solvers.Dirac3DirectSolver

  Bases: `ModelSolver`

  ip_addr: str = 'localhost'
  port: str = '50051'
  cert_file: str = None
  solve(model, sum_constraint, relaxation_schedule=1, num_samples=1, mean_photon_number=None, quantum_fluctuation_coefficient=None) → SolutionResults
  connect(ip_addr: str, port: str, cert_file: str | None = None)
  property client: EqcClient

    Return a new client from eqc-direct based on class config.

  makeResults(model: EqcModel, response: Dict)

    Builds the results object

eqc_models.solvers.Dirac3ContinuousDirectSolver

  alias of `Dirac3DirectSolver`

class eqc_models.solvers.Dirac3IntegerDirectSolver

  Bases: `Dirac3DirectSolver`

  solve(model, relaxation_schedule=1, num_samples=1, mean_photon_number=None, quantum_fluctuation_coefficient=None) → SolutionResults

class eqc_models.solvers.Dirac1CloudSolver(url=None, api_token=None)

  Bases: `Dirac1Mixin`, `QuboSolverMixin`, `QciClientSolver`

  ## Overview

  Dirac1CloudSolver is a class that encapsulates the different calls to Qatalyst for Dirac-1 jobs, which are quadratic binary optimization problems.

  Examples

      >>> C = np.array([[-1], [-1]])
      >>> J = np.array([[0, 1.0], [1.0, 0]])
      >>> from eqc_models.base.quadratic import QuadraticModel
      >>> model = QuadraticModel(C, J)
      >>> model.upper_bound = np.array([1, 1]) # set the domain maximum per variable
      >>> solver = Dirac1CloudSolver()
      >>> results = solver.solve(model, num_samples=5)
      2... submitted... COMPLETED...
      >>> results["results"]["energies"][0] <= 1.0 # old style response
      True
      >>> results.energies[0] <= 1.0 # SolutionResults
      True

class eqc_models.solvers.Dirac3CloudSolver(url=None, api_token=None)

  Bases: `Dirac3Mixin`, `QciClientSolver`

  Dirac3CloudSolver is a class that encapsulates the different calls to Qatalyst for Dirac-3 jobs. Currently, there are two different jobs, one for integer and another for continuous solutions. Calling the solve method with different arguments controls which job is submitted. The continuous job requires `sum_constraint`. The integer job does not accept this parameter, so specifying a sum constraint forces the job type to be continuous, and not specifying it results in the integer job being called.

  ## Continuous Solver

  Utilizing Dirac-3 as a continuous solver involves encoding the variables in single time bins with the values of each determined by a normalized photon count value.

  ## Integer Solver

  Utilizing Dirac-3 as an integer solver involves encoding the variables in multiple time bins, each representing a certain value for that variable, or “qudit”.

  Examples

      >>> C = np.array([[1], [1]])
      >>> J = np.array([[-1.0, 0], [0, -1.0]])
      >>> from eqc_models.base.quadratic import QuadraticModel
      >>> model = QuadraticModel(C, J)
      >>> model.upper_bound = np.array([1, 1]) # set the domain maximum per variable
      >>> solver = Dirac3CloudSolver()
      >>> results = solver.solve(model, sum_constraint=1, relaxation_schedule=1)
      2... submitted... COMPLETED...
      >>> results["results"]["energies"][0] <= 1.0
      True
      >>> results.energies[0] <= 1.0
      True
      >>> C = np.array([-1, -1], dtype=np.float32)
      >>> J = np.array([[0, 1], [1, 0]], dtype=np.float32)
      >>> model = QuadraticModel(C, J)
      >>> model.upper_bound = np.array([1, 1]) # set the domain maximum per variable
      >>> results = solver.solve(model, relaxation_schedule=1)
      2... submitted... COMPLETED...
      >>> results["results"]["energies"][0] == -1.0
      True
      >>> results.energies[0] == -1.0
      True

  job_type = 'hamiltonian'
  job_params_names = ['num_samples', 'relaxation_schedule', 'mean_photon_number', 'quantum_fluctuation_coefficient', 'num_levels', 'sum_constraint']
  solve(model: EqcModel, name: str | None = None, tags: List | None = None, sum_constraint: float | None = None, relaxation_schedule: int | None = None, num_samples: int = 1, wait: bool = True, mean_photon_number: float | None = None, quantum_fluctuation_coefficient: int | None = None, **job_kwargs)
    Parameters:

      - model (EqcModel) – a model object which supplies a hamiltonian operator for the device to sample. Must support the polynomial operator property.

      - tags (List) – a list of strings to save with the job

      - sum_constraint (float) – a value which applies a constraint to the solution, forcing all variables to sum to this value, changes method to continuous solver

      - relaxation_schedule (int) – a predefined schedule indicator which sets parameters on the device to control the sampling through photon measurement

      - num_samples (int) – the number of samples to take, defaults to 1

      - wait (bool) – a flag for waiting for the response or letting it run asynchronously. Asynchronous runs must retrieve results directly using qci-client and the job_id.

      - mean_photon_number (float) – an optional decimal value which sets the average number of photons that are present in a given quantum state. Modify this value to control the relaxation schedule more precisely than the four presets given in schedules 1 through 4. Allowed values are decimals between 0.000133333333333 and 0.001.

      - quantum_fluctuation_coefficient (int) – an integer value which Sets the amount of loss introduced into the system for each loop during the measurement process. Modify this value to control the relaxation schedule more precisely than the four presets given in schedules 1 through 4. Allowed values range from 1 to 50.

class eqc_models.solvers.QciClientSolver(url=None, api_token=None)

  Bases: `QciClientMixin`, `ModelSolver`

  Parameters:

    - url (string) – optional value specifying the QCi API URL

    - api_token (string) – optional value specifying the authentication token for the QCi API

    - the (QCi API client wrapper for solving an EQC model. This class provides)

    - file (common method for uploading a file to the API for solving. Since the)

    - types (types change for the job)

    - are (the specific files required for the job)

    - method. (specified in subclasses within the uploadFiles)

  static uploadFile(file_data: ndarray, file_name: str | None = None, file_type: str | None = None, client: OptimizationClient | None = None) → str

    Upload the operator file, return the file ID.

    Parameters:

      - file_data (numpy array, dictionary or list) – contains file data to be uploaded

      - file_name (str) – Name of the file to be uploaded.

      - file_type (str) – Type of the file to be uploaded.

      - client (QciClient) – QciClient instance

  uploadJobFiles(client: OptimizationClient, model: EqcModel)
  checkModel(model)
    Parameters:

      - model (EqcModel) – Instance of a model to validate against the solver requirements

      - solver. (This method raises an exception if the model supplied does not meet the requirements for the)

      - the (One of the validations is that the model supplies the operator that the solver uses. This is)

      - supply (qubo operator for Dirac-1 and polynomial operator for Dirac-3. If the model does not)

      - operator (the)

      - explanation. (then the solver cannot accept it and the method will fail with an)

      - instance (Another validation is for the allowed upper bound of a model. For)

      - only (solvers which)

      - 1. (handle binary variables can only accept models with variabes having an upper bound of)

  solve(model: EqcModel, name: str | None = None, tags: List | None = None, num_samples: int = 1, wait=True, job_type=None, **job_kwargs) → SolutionResults
    Parameters:

      - model (EqcModel) – Instance of a model for solving.

      - name (str) – Name of the job; default is None.

      - tags (list) – A list of job tags; default is None.

      - num_samples (int) – Number of samples used; default is 1.

      - wait (bool) – The wait flag indicating whether to wait for the job to complete before returning the complete job data otherwise return a job ID as soon as a job is submitted; default is True.

      - job_type (str) – Type of the job; default is None. When None, it is constructed from the instance job_type property.

    Returns:

      - job response dictionary

      - This method takes the particulars of the instance model and handles

      - the QciClient.solve call.

  getResults(response: Dict) → Dict[str, List]

    Extract the results from response.

    Parameters:

      response (The responce from QciClient.)

    Return type:

      The results json object.

class eqc_models.solvers.Dirac3IntegerCloudSolver(url=None, api_token=None)

  Bases: `Dirac3Mixin`, `QciClientSolver`

      >>> C = np.array([-1, -1], dtype=np.float32)
      >>> J = np.array([[0, 1], [1, 0]], dtype=np.float32)
      >>> from eqc_models.base.quadratic import QuadraticModel
      >>> model = QuadraticModel(C, J)
      >>> model.upper_bound = np.array([1, 1]) # set the domain maximum per variable
      >>> solver = Dirac3IntegerCloudSolver()
      >>> model = QuadraticModel(C, J)
      >>> model.upper_bound = np.array([1, 1]) # set the domain maximum per variable
      >>> results = solver.solve(model, relaxation_schedule=1)
      2... submitted... COMPLETED...
      >>> results["results"]["energies"][0] == -1.0
      True
      >>> results.energies[0] == -1.0
      True

  job_type = 'hamiltonian-integer'
  job_params_names = ['num_samples', 'relaxation_schedule', 'mean_photon_number', 'quantum_fluctuation_coefficient', 'num_levels']
  solve(model: EqcModel, name: str | None = None, tags: List | None = None, relaxation_schedule: int | None = None, num_samples: int = 1, wait: bool = True, mean_photon_number: float | None = None, quantum_fluctuation_coefficient: int | None = None, **job_kwargs)
    Parameters:

      - model (EqcModel) – a model object which supplies a hamiltonian operator for the device to sample. Must support the polynomial operator property.

      - tags (List) – a list of strings to save with the job

      - relaxation_schedule (int) – a predefined schedule indicator which sets parameters on the device to control the sampling through photon measurement

      - num_samples (int) – the number of samples to take, defaults to 1

      - wait (bool) – a flag for waiting for the response or letting it run asynchronously. Asynchronous runs must retrieve results directly using qci-client and the job_id.

      - mean_photon_number (float) – an optional decimal value which sets the average number of photons that are present in a given quantum state. Modify this value to control the relaxation schedule more precisely than the four presets given in schedules 1 through 4. Allowed values are decimals between 0.000133333333333 and 0.001.

      - quantum_fluctuation_coefficient (int) – an integer value which Sets the amount of loss introduced into the system for each loop during the measurement process. Modify this value to control the relaxation schedule more precisely than the four presets given in schedules 1 through 4. Allowed values range from 1 to 50.

    Dirac3IntegerCloudSolver is a class that encapsulates the different calls to Qatalyst for Dirac-3 jobs. Utilizing Dirac-3 as an integer solver involves encoding the variables in multiple time bins, each representing a certain value for that variable, or “qudit”.

class eqc_models.solvers.Dirac3ContinuousCloudSolver(url=None, api_token=None)

  Bases: `Dirac3Mixin`, `QciClientSolver`

      >>> C = np.array([[1], [1]])
      >>> J = np.array([[-1.0, 0], [0, -1.0]])
      >>> from eqc_models.base.quadratic import QuadraticModel
      >>> model = QuadraticModel(C, J)
      >>> model.upper_bound = np.array([1, 1]) # set the domain maximum per variable
      >>> solver = Dirac3ContinuousCloudSolver()
      >>> results = solver.solve(model, sum_constraint=1, relaxation_schedule=1)
      2... submitted... COMPLETED...
      >>> results["results"]["energies"][0] <= 1.0
      True
      >>> results.energies[0] <= 1.0
      True

  job_type = 'hamiltonian'
  job_params_names = ['num_samples', 'relaxation_schedule', 'mean_photon_number', 'quantum_fluctuation_coefficient', 'sum_constraint']
  solve(model: EqcModel, name: str | None = None, tags: List | None = None, sum_constraint: float | None = None, relaxation_schedule: int | None = None, num_samples: int = 1, wait: bool = True, mean_photon_number: float | None = None, quantum_fluctuation_coefficient: int | None = None, **job_kwargs)
    Parameters:

      - model (EqcModel) – a model object which supplies a hamiltonian operator for the device to sample. Must support the polynomial operator property.

      - tags (List) – a list of strings to save with the job

      - sum_constraint (float) – a value which applies a constraint to the solution, forcing all variables to sum to this value

      - relaxation_schedule (int) – a predefined schedule indicator which sets parameters on the device to control the sampling through photon measurement

      - num_samples (int) – the number of samples to take, defaults to 1

      - wait (bool) – a flag for waiting for the response or letting it run asynchronously. Asynchronous runs must retrieve results directly using qci-client and the job_id.

      - mean_photon_number (float) – an optional decimal value which sets the average number of photons that are present in a given quantum state. Modify this value to control the relaxation schedule more precisely than the four presets given in schedules 1 through 4. Allowed values are decimals between 0.000133333333333 and 0.001.

      - quantum_fluctuation_coefficient (int) – an integer value which Sets the amount of loss introduced into the system for each loop during the measurement process. Modify this value to control the relaxation schedule more precisely than the four presets given in schedules 1 through 4. Allowed values range from 1 to 50.

class eqc_models.solvers.MIPMixin

  Bases: `object`

  Implements a solve method which intercepts the operator and builds a new model with added terms for restriction of discrete variables with an upper bound of 1 to take on only values 0 or 1 at minima.

  Following the submission of a new model with the added penalties, the solutions are updated to exclude added slack variables.

  This is only supported with continuous-capable devices.

  solve(model: EqcModel, *args, **kwargs)

class eqc_models.solvers.Dirac3MIPCloudSolver(url=None, api_token=None)

  Bases: `MIPMixin`, `Dirac3ContinuousCloudSolver`

class eqc_models.solvers.Dirac3MIPDirectSolver

  Bases: `MIPMixin`, `Dirac3DirectSolver`

# eqc_models.ml

class eqc_models.ml.QBoostClassifier(relaxation_schedule=2, num_samples=1, solver_access='cloud', api_url=None, api_token=None, ip_addr=None, port=None, lambda_coef=0, weak_cls_schedule=1, weak_cls_type='lg', weak_cls_params={}, weak_cls_strategy='multi_processing', weak_cls_num_jobs=None, weak_cls_pair_count=None)

  Bases: `ClassifierBase`

  An implementation of QBoost classifier that uses QCi’s Dirac-3.

  Parameters:

    - relaxation_schedule (Relaxation schedule used by Dirac-3;)

    - default (lg.)

    - num_samples (Number of samples used by Dirac-3; default: 1.)

    - solver_access (Solver access type: cloud or direct; default: cloud.)

    - api_url (API URL used when cloud access is used; default: None.)

    - api_token (API token used when cloud access is used; default: None.)

    - ip_addr (IP address of the device when direct access is used; default: None.)

    - port (Port number of the device when direct access is used; default: None.)

    - lambda_coef (A penalty multiplier; default: 0.)

    - weak_cls_schedule (Weak classifier schedule. Is either 1, 2,)

    - default

    - weak_cls_type (Type of weak classifier) –

      - dct: Decison tree classifier

      - nb: Naive Baysian classifier

      - lg: Logistic regression

      - gp: Gaussian process classifier

      - knn: K-nearest neighbors classifier

      - lda: Linear discriminant analysis classifier

      - qda: Quadratic discriminant analysis classifier

      - lgb: Light-GBM classifier

      - xgb: XGBoost classifier

    - default

    - weak_cls_params (Dict of weak classifier parameters. Default: {};)

    - parameters. (use default)

    - weak_cls_strategy (Computation strategy for weak classifier)

    - training

    - Default (weak_cls_strategy="multi_processing".)

    - "multi_processing".

    - weak_cls_num_jobs (Number of jobs when)

    - Default

    - cores. (available)

  Examples

      >>> from sklearn import datasets
      >>> from sklearn.preprocessing import MinMaxScaler
      >>> from sklearn.model_selection import train_test_split
      >>> iris = datasets.load_iris()
      >>> X = iris.data
      >>> y = iris.target
      >>> scaler = MinMaxScaler()
      >>> X = scaler.fit_transform(X)
      >>> for i in range(len(y)):
      ...     if y[i] == 0:
      ...         y[i] = -1
      ...     elif y[i] == 2:
      ...         y[i] = 1
      >>> X_train, X_test, y_train, y_test = train_test_split(
      ...     X,
      ...     y,
      ...     test_size=0.2,
      ...     random_state=42,
      ... )
      >>> from eqc_models.ml.classifierqboost import QBoostClassifier
      >>> obj = QBoostClassifier(
      ...     relaxation_schedule=2,
      ...     num_samples=1,
      ...     lambda_coef=0.0,
      ... )
      >>> from contextlib import redirect_stdout
      >>> import io
      >>> f = io.StringIO()
      >>> with redirect_stdout(f):
      ...    obj.fit(X_train, y_train)
      ...    y_train_prd = obj.predict(X_train)
      ...    y_test_prd = obj.predict(X_test)

  topNPairs(X, n)
  fit(X, y)

    Build a QBoost classifier from the training set (X, y).

    Parameters:

      - X ({array-like, sparse matrix} of shape (n_samples, n_features))

      - samples. (The training input)

      - y (array-like of shape (n_samples,))

      - values. (The target)

    Return type:

      Response of Dirac-3 in JSON format.

  predict_raw(X: array)

    Predict raw output of the classifier for input X.

    Parameters:

      X ({array-like, sparse matrix} of shape (n_samples, n_features))

    Returns:

      - y (ndarray of shape (n_samples,))

      - The predicted raw output of the classifier.

  predict(X: array, threshold=0)

    Predict classes for X.

    Parameters:

      - X ({array-like, sparse matrix} of shape (n_samples, n_features))

      - threshold (Prediction threshold.)

    Returns:

      - y (ndarray of shape (n_samples,))

      - The predicted classes.

  get_hamiltonian(X: array, y: array)
  convert_sol_to_params(sol)

class eqc_models.ml.QSVMClassifier(relaxation_schedule=1, num_samples=1, solver_access='cloud', api_url=None, api_token=None, ip_addr=None, port=None, lambda_coef=1.0)

  Bases: `ClassifierBase`

  An implementation of QSVM classifier that uses QCi’s Dirac-3.

  Parameters:

    - relaxation_schedule (Relaxation schedule used by Dirac-3; default:)

    - 2.

    - num_samples (Number of samples used by Dirac-3; default: 1.)

    - solver_access (Solver access type: cloud or direct; default: cloud.)

    - api_url (API URL used when cloud access is used; default: None.)

    - api_token (API token used when cloud access is used; default: None.)

    - ip_addr (IP address of the device when direct access is used; default: None.)

    - port (Port number of the device when direct access is used; default: None.)

    - lambda_coef (The penalty multipler)

  Examples

      >>> from sklearn import datasets
      >>> from sklearn.preprocessing import MinMaxScaler
      >>> from sklearn.model_selection import train_test_split
      >>> iris = datasets.load_iris()
      >>> X = iris.data
      >>> y = iris.target
      >>> scaler = MinMaxScaler()
      >>> X = scaler.fit_transform(X)
      >>> for i in range(len(y)):
      ...     if y[i] == 0:
      ...         y[i] = -1
      ...     elif y[i] == 2:
      ...         y[i] = 1
      >>> X_train, X_test, y_train, y_test = train_test_split(
      ...     X,
      ...     y,
      ...     test_size=0.2,
      ...     random_state=42,
      ... )
      >>> from eqc_models.ml.classifierqsvm import QSVMClassifier
      >>> obj = QSVMClassifier(
      ...     relaxation_schedule=2,
      ...     num_samples=1,
      ... )
      >>> from contextlib import redirect_stdout
      >>> import io
      >>> f = io.StringIO()
      >>> with redirect_stdout(f):
      ...    obj = obj.fit(X_train, y_train)
      ...    y_train_prd = obj.predict(X_train)
      ...    y_test_prd = obj.predict(X_test)

  fit(X, y)

    Build a QSVM classifier from the training set (X, y).

    Parameters:

      - X ({array-like, sparse matrix} of shape (n_samples, n_features))

      - samples. (The training input)

      - y (array-like of shape (n_samples,))

      - values. (The target)

    Return type:

      Response of Dirac-3 in JSON format.

  predict_raw(X: array)

    Predict classes for X.

    Parameters:

      X ({array-like, sparse matrix} of shape (n_samples, n_features))

    Returns:

      - y (ndarray of shape (n_samples,))

      - The predicted classes.

  predict(X: array)

    Predict classes for X.

    Parameters:

      X ({array-like, sparse matrix} of shape (n_samples, n_features))

    Returns:

      - y (ndarray of shape (n_samples,))

      - The predicted classes.

  get_hamiltonian(X: array, y: array)
  convert_sol_to_params(sol)

class eqc_models.ml.PCA(n_components=None, relaxation_schedule=2, num_samples=1, solver_access='cloud', api_url=None, api_token=None, ip_addr=None, port=None, mode='largest_components')

  Bases: `DecompBase`

  An implementation of Principal component analysis (PCA) that uses QCi’s Dirac-3.

  Linear dimensionality reduction using Singular Value Decomposition of the data to project it to a lower dimensional space.

  Parameters:

    - n_components (Number of components to keep; if n_components is not)

    - default (largest_components vs. smallest_components;)

    - relaxation_schedule (Relaxation schedule used by Dirac-3; default:)

    - 2.

    - num_samples (Number of samples used by Dirac-3; default: 1.)

    - solver_access (Solver access type: cloud or direct; default: cloud.)

    - api_url (API URL used when cloud access is used; default: None.)

    - api_token (API token used when cloud access is used; default: None.)

    - ip_addr (IP address of the device when direct access is used; default: None.)

    - port (Port number of the device when direct access is used; default: None.)

    - mode (Compute the largest or smallest principal components,)

    - default

    - largest_components.

  Examples

      >>> from sklearn import datasets
      >>> iris = datasets.load_iris()
      >>> X = iris.data
      >>> from sklearn.preprocessing import StandardScaler
      >>> scaler = StandardScaler()
      >>> X = scaler.fit_transform(X)
      >>> from eqc_models.ml.decomposition import PCA
      >>> from contextlib import redirect_stdout
      >>> import io
      >>> f = io.StringIO()
      >>> with redirect_stdout(f):
      ...    obj = PCA(
      ...        n_components=4,
      ...        relaxation_schedule=2,
      ...        num_samples=1,
      ...    )
      ...    X_pca = obj.fit_transform(X)

  fit(X)

    Build a PCA object from the training set X.

    Parameters:

      - X ({array-like, sparse matrix} of shape (n_samples, n_features))

      - samples. (The training input)

    Returns:

      - responses.

      - A dirct containing Dirac responses.

  transform(X: array)

    Apply dimensionality reduction to X.

    X is projected on the first principal components previously extracted from a training set.

    Parameters:

      - X (array-like of shape (n_samples, n_features))

      - data (New)

      - samples (where n_samples is the number of)

      - features. (and n_features is the number of)

    Returns:

      - X_new (array-like of shape (n_samples, n_components))

      - Projection of X in the first principal components, where n_samples

      - is the number of samples and n_components is the number of the components.

  fit_transform(X)

    Fit the model with X and apply the dimensionality reduction on X.

    Parameters:

      - X (array-like of shape (n_samples, n_features))

      - data (Training)

      - samples (where n_samples is the number of)

      - features. (and n_features is the number of)

    Returns:

      - X_new (ndarray of shape (n_samples, n_components))

      - Transformed values.

class eqc_models.ml.ReservoirForecastModel(ip_addr, num_nodes, feature_scaling, num_pads: int = 0, reg_coef: float = 0.0, device: str = 'EmuCore')

  Bases: `BaseForecastModel`, `QciReservoir`

  A reservoir based forecast model.

  Parameters:

    - ip_addr (The IP address of the device.)

    - num_nodes (Number of reservoir network nodes.)

    - feature_scaling (The factor used to scale the reservoir output.)

    - num_pads (Size of the pad used in the reservoir input;)

    - default (supported;)

    - reg_coef (L2 regularization coefficient for linear regression;)

    - default

    - device (The QCi reservoir device. Currently only 'EmuCore' is)

    - default

  Examples

      >>> MAX_TRAIN_DAY = 800
      >>> IP_ADDR = "172.22.19.49"
      >>> FEATURE_SCALING = 0.1
      >>> NUM_NODES = 1000
      >>> NUM_PADS = 100
      >>> LAGS = 2
      >>> from contextlib import redirect_stdout
      >>> import io
      >>> f = io.StringIO()
      >>> from eqc_models.ml import ReservoirForecastModel
      >>> with redirect_stdout(f):
      ...    model = ReservoirForecastModel(
      ...        ip_addr=IP_ADDR,
      ...        num_nodes=NUM_NODES,
      ...        feature_scaling=FEATURE_SCALING,
      ...        num_pads=NUM_PADS,
      ...        device="EmuCore",
      ...    )
      ...    model.fit(
      ...        data=train_df,
      ...        feature_fields=["norm_cell_prod"],
      ...        target_fields=["norm_cell_prod"],
      ...        lags=LAGS,
      ...        horizon_size=1,
      ...    )
      ...    y_train_pred = model.predict(train_df, mode="in_sample")
      ...    y_test_pred = model.predict(test_df, mode="in_sample")
      >>> model.close()

  close()
  fit(data: DataFrame, feature_fields: list, target_fields: list, lags: int = 0, horizon_size: int = 1)

    A function to train a forecast model.

    Parameters:

      - data (A pandas data frame that contain the time series.)

      - feature_fields (A list of fields in the data frame that are as)

      - reservoir. (inputs to the)

      - target_fields (A list of fields in teh data frame that are to be)

      - forecasted.

      - lags (Number of lags used; default = 0.)

      - horizon_size (Size of the horizon, e.g. number of forecast)

      - steps.

  predict(data: DataFrame, pad_mode: str = 'zero', mode: str = 'in_sample')

    A function to get predictions from forecast model.

    Parameters:

      - data (A pandas data frame that contain the time series.)

      - pad_mode (Mode of the reservoir input padding, either)

      - default (sample (used for testing);)

      - mode (A value of 'out_of_sample' predicts the horizon)

      - in (following the time series. A value of 'in_sample' predicts)

      - default

    Returns:

      The predictions

    Return type:

      numpy.array((horizon_size, num_dims)).

# eqc_models.allocation

class eqc_models.allocation.AllocationModel(resources: List, tasks: List, resource_usage: List, resource_limits: List, cost_per_task: List)

  Bases: `ConstrainedQuadraticModel`

  Parameters:

    - resources (List) – names of available resources.

    - tasks (List) – names of tasks.

    - resource_usage (List of Lists or 2D np.ndarray) – rows represent tasks and columns represent resources, specifying the amount of each resource required per task.

    - resource_limits (1D array or List) – specifying the limit on each resource.

    - cost_per_task (1D array List) – specifying the cost per task (or benefit with negative of value).

  penalty_multiplier

    value for weighting the penalties formed from the equality constraints

    Type:

      float

  qubo

    QUBO operator representation

    Type:

      eqc_models.base.operators.QUBO

  polynomial

    Polynomial operator representation

    Type:

      eqc_models.base.operators.Polynomial

  This class represents a resource allocation model for maximizing total benefit. In other words, Given a list of resources and a list of tasks, allocate the resources among the tasks so as to maximize the economic benefit.

  Here’s an example. Five tasks must share 4 resources. Each task can use a different amount of each resource.

  | Spam | Eggs | Coconut | Sparrow | Benefit
  Breakfast | 1 | 2 | 0 | 0 | 3
  Countryside Stroll | 0 | 0 | 1 | 0 | 1
  Storm Castle | 0 | 12 | 1 | 1 | 10
  Availability | 1 | 12 | 2 | 1 |

      >>> resources = ["Spam", "Eggs", "Coconut", "Sparrow"]
      >>> tasks = ["Breakfast", "Countryside Stroll", "Storm Castle"]
      >>> resource_usage = [[1, 2, 0, 0], [0, 0, 1, 0], [0, 12, 1, 1]]
      >>> resource_limits = [1, 12, 2, 1]
      >>> cost_per_task = [-3, -1, -10.]
      >>> allocation_model = AllocationModel(resources, tasks, resource_usage, resource_limits, cost_per_task)
      >>> allocation_model.penalty_multiplier = 1
      >>> C, J = allocation_model.H
      >>> C # -3 -2 * (12 * 2 + 1 * 1), -1 -2 * 2*1, -10 -2 * (12 * 12 + 1 * 2 + 1 * 1)
      ...
      array([ -53., -5.,     -304.])
      >>> J
      array([[ 5., 0., 24.],
             [ 0., 1.,  1.],
             [ 24., 1., 146.]])

  property upper_bound: array

    An array of upper bound values for every variable in the model. Must be integer.

  property variables
  add_task(task: str, task_resource_usage: List, task_cost: float)
    task: str

      Name of task to add

    task_resource_usage: List

      Quantity of resource used for task for all tasks

    task_cost: float

      Quantity indicating the cost or benefit (negative) of the task

    Add a task to the problem, modifying the resource usage and task cost arrays.

  property H: Tuple[ndarray, ndarray]

    Return linear, quadratic portions of the (quadratic) Hamiltonian

class eqc_models.allocation.AllocationModelX(resources: List, tasks: List, resource_usage: List, resource_limits: List, resource_rule: List[ResourceRuleEnum], cost_per_task: List)

  Bases: `InequalitiesMixin`, `AllocationModel`

  Parameters:

    - resources (List) – names of available resources.

    - tasks (List) – names of tasks.

    - resource_usage (List of Lists or 2D np.ndarray) – rows represent tasks and columns represent resources, specifying the amount of each resource required per task.

    - resource_limits (1D array or List) – specifying the limit on each resource.

    - resource_rule (List) – ResourceRuleEnum values for each resource

    - cost_per_task (1D array List) – specifying the cost per task (or benefit with negative of value).

  penalty_multiplier

    value for weighting the penalties formed from the equality constraints

    Type:

      float

  qubo

    QUBO oeprator representation

    Type:

      eqc_models.base.operators.QUBO

  polynomial

    Polynomial operator representation

    Type:

      eqc_models.base.operators.Polynomial

  variables

    names of variables formed from tasks and assignments

    Type:

      List

  This class represents a resource allocation model for maximizing total benefit. In other words, Given a list of resources and a list of tasks, allocate the resources among the tasks so as to maximize the economic benefit.

  Adds resource_rule as an argument. This must be a list of strings specifying constraints for each resource (LE, GE, or EQ).

  Here’s an example. Five tasks must share 4 resources. Each task can use a different amount of each resource.

  | Spam | Eggs | Coconut | Sparrow | Benefit
  Breakfast | 1 | 2 | 0 | 0 | 3
  Countryside Stroll | 0 | 0 | 1 | 0 | 1
  Storm Castle | 0 | 12 | 1 | 1 | 10
  Availability | 1 | 12 | 2 | 1 |

      >>> resources = ["Spam", "Eggs", "Coconut", "Sparrow"]
      >>> tasks = ["Breakfast", "Countryside Stroll", "Storm Castle"]
      >>> resource_usage = [[1, 2, 0, 0], [0, 0, 1, 0], [0, 12, 1, 1]]
      >>> resource_limits = [1, 12, 2, 1]
      >>> cost_per_task = [-3, -1, -10.]
      >>> resource_rules = [ResourceRuleEnum.MAXIMUM for i in range(len(resources))]
      >>> allocation_model = AllocationModelX(resources, tasks, resource_usage, resource_limits, resource_rules, cost_per_task)
      >>> allocation_model.penalty_multiplier = 1
      >>> C, J = allocation_model.H
      >>> C # -3 -2 * (12 * 2 + 1 * 1), -1 -2 * 2*1, -10 -2 * (12 * 12 + 1 * 2 + 1 * 1)
      ...
      array([ -53., -5., -304., -2., -24., -4., -2.])
      >>> J
      array([[  5.,  0.,  24., 1.,  2., 0., 0.],
             [  0.,  1.,   1., 0.,  0., 1., 0.],
             [ 24.,  1., 146., 0., 12., 1., 1.],
             [  1.,  0.,   0., 1.,  0., 0., 0.],
             [  2.,  0.,  12., 0.,  1., 0., 0.],
             [  0.,  1.,   1., 0.,  0., 1., 0.],
             [  0.,  0.,   1., 0.,  0., 0., 1.]])

  property upper_bound: array

    An array of upper bound values for every variable in the model. Must be integer.

  property linear_objective: ndarray

    Returns a 1D numpy array representing the linear part of the objective function (total profit).

  property quad_objective: ndarray

    Returns a 2D numpy array representing the quadratic part of the objective function (always zero in this case).

  property H

    Overrides the parent build method to incorporate slack variables based on resource_rule.

  checkResources(solution)
    Parameters:

      solution (List or np.ndarray) – solution vector to check for resource violations

    Returns:

      List of violations

    Return type:

      (name, rule, violation quantity, message)

  property n

    Return the number of variables

class eqc_models.allocation.PortMomentum(stocks: list, stock_data_dir: str, adj_date: str, lookback_days: int = 60, window_days: int = 30, window_overlap_days: int = 15, weight_upper_limit: float = 0.08, r_base: float = 0.00013698630136986303, alpha: float = 5.0, beta: float = 1.0, xi: float = 1.0)

  Bases: `PortBase`

  get_hamiltonian(return_df, min_date, max_date)

class eqc_models.allocation.ResourceRuleEnum(value)

  Bases: `Enum`

  Enumeration of the allowed resource rules, mapping to the mathematical expression:

  MAXIMUM -> LE (less than or equal to) MINIMUM -> GE (greater than or equal to) EXACT -> EQ (equal to)

  MAXIMUM = 'LE'
  MINIMUM = 'GE'
  EXACT = 'EQ'

# eqc_models.assignment

class eqc_models.assignment.QAPModel(A: ndarray, B: ndarray, C: ndarray)

  Bases: `ConstrainedQuadraticModel`

  Parameters:

    - A (np.array) – matrix of distances or costs between locations

    - B (np.array) – matrix of flows between facilities

    - C (np.array) – matrix pf fixed costs of assigning facilities to locations

  Formulates a quadratic programming model from three inputs. The inputs are composed of: a matrix which describes the unit cost of moving a single asset between locations. This matrix can describe asynchronous costs that differ by direction of flow; a matrix which describes the quantity of flow between facilities. This matrix describes the quantity in the direction of flow; a matrix which desribes the fixed cost of assigning a facility to a location.

  Using these flows and costs, an optimization problem is formulated to determine a solution which minimizes the total cost for positioning facilities at locations. The facilities do not have to be buildings or campuses and the locations do not have to be geographical locations. See Beckmann and and Koopmans (1957) for the original reference.

      >>> 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)
      >>> model.penalty_multiplier = 105.625
      >>> Q = model.qubo.Q
      >>> np.sum(Q)
      24318.03

class eqc_models.assignment.ResourceAssignmentModel(resources, tasks)

  Bases: `InequalitiesMixin`, `ConstrainedQuadraticModel`

  Resource assignment model

  Parameters:

    - resources (List)

    - tasks (List)

      >>> # name is not a required attribute of the resources or tasks
      >>> crews = [{"name": "Maintenance Crew 1", "skills": ["A", "F"], "capacity": 5, "cost": 4},
      ...          {"name": "Baggage Crew 1", "skills": ["B"], "capacity": 4, "cost": 1},
      ...          {"name": "Maintenance Crew 2", "skills": ["A", "F"], "capacity": 5, "cost": 2}]
      >>> tasks = [{"name": "Refuel", "skill_need": "F", "load": 3},
      ...          {"name": "Baggage", "skill_need": "B", "load": 1}]
      >>> model = ResourceAssignmentModel(crews, tasks)
      >>> assignments = model.createAssignmentVars()
      >>> assignments
      [{'resource': 0, 'task': 0}, {'resource': 1, 'task': 1}, {'resource': 2, 'task': 0}]
      >>> A, b, senses = model.constrainAssignments(assignments)
      >>> A
      array([[3., 0., 0.],
             [0., 1., 0.],
             [0., 0., 3.],
             [3., 0., 3.],
             [0., 3., 0.]], dtype=float32)
      >>> b
      array([5., 4., 5., 3., 3.], dtype=float32)
      >>> senses
      ['LE', 'LE', 'LE', 'EQ', 'EQ']
      >>> A, b = model.constraints
      >>> A
      array([[3., 0., 0., 1., 0., 0.],
             [0., 1., 0., 0., 1., 0.],
             [0., 0., 3., 0., 0., 1.],
             [3., 0., 3., 0., 0., 0.],
             [0., 3., 0., 0., 0., 0.]])

  classmethod checkTasks(tasks)
  createAssignmentVars()

    Examine all combinatins of possible crew-task assignments

  constrainAssignments(assignments: List) → List

    Examine the assignments to determine the necessary constraints to ensure feasibility of solution.

  property sum_constraint: int

    This value is a suggestion which should be used with a machine slack

  decode(solution: array) → List[Dict]

    Convert the binary solution into a list of tasks

# eqc_models.combinatorics

class eqc_models.combinatorics.SetCoverModel(subsets, weights)

  Bases: `ConstrainedQuadraticModel`

  Parameters:

    - subsets (List) – List of sets where the union of all sets is S

    - weights (List) – List of weights where each weight is the cost of choosing the subset corresponding to the index of the weight.

      >>> X = [set(['A', 'B']), set(['B', 'C']), set(['C'])]
      >>> weights = [2, 2, 1]
      >>> model = SetCoverModel(X, weights)
      >>> model.penalty_multiplier = 8
      >>> lhs, rhs = model.constraints
      >>> lhs
      array([[ 1,  0,  0,  0,  0],
             [ 1,  1,  0, -1,  0],
             [ 0,  1,  1,  0, -1]])
      >>> from eqc_models.solvers import Dirac3IntegerCloudSolver
      >>> solver = Dirac3IntegerCloudSolver()
      >>> response = solver.solve(model, relaxation_schedule=1, num_samples=5)
      20...
      >>> solutions = response["results"]["solutions"]
      >>> solutions[0]
      [1, 0, 1, 0, 0]
      >>> selections = model.decode(solutions[0])
      >>> assert {'B', 'A'} in selections
      >>> assert {'C'} in selections
      >>> assert len(selections) == 2

  decode(solution) → List

    Manipulate the solution to match the variable count. Optionally, a log-encoded solution vector can be translated into a vector with integral values.

    Parameters:

      - solution (np.ndarray) – 1d solution vector

      - from_encoding (str) – string indicating if the vector is an encoding of particular type. The text ‘qubo’ is the only option that does anything for now.

      - EqcModel() (>>> model =)

      - np.array([1 (>>> solution =)

      - 2

      - 5])

      - ub (>>> model.upper_bound =)

      - 1 (>>> model.machine_slacks =)

      - np.array([1

      - 1

      - 1

      - 1

      - 0

      - 1

      - 3 (0]) # this solution has a)

      - bound (which violates the upper)

      - care (but decode doesn't)

      - model.decode(solution (>>>)

      - "qubo")

      - array([1

      - 3

      - 5])

class eqc_models.combinatorics.SetPartitionModel(subsets: List[set], weights: List[float])

  Bases: `ConstrainedPolynomialModel`

  This class represents a set partitioning model for optimization problems that require selecting subsets to partition a universal set while minimizing an objective function defined by weights. Given a collection of subsets, a weight is assigned to each subset, and the goal is to determine an optimal selection of subsets to fully cover the universal set with minimized total weight.

  Parameters:

    - subsets (List of sets) – List of sets (subsets) defining the collection to partition. Each element in this list represents a subset containing elements from the universal set.

    - weights (List of floats) – List of weights corresponding to each subset. The length of this list should be equal to the number of subsets, with each weight indicating the cost or weight associated with selecting a particular subset.

  H

    Tuple containing the linear (h) and quadratic (J) coefficients for the Hamiltonian representation of the quadratic problem formulation.

    Type:

      tuple of np.ndarray

  penalty_multiplier

    Value for weighting the penalties formed from the equality constraints.

    Type:

      float

  polynomial

    Polynomial operator representation for the model, constructed from the penalty terms to encode the set partition constraints.

    Type:

      Polynomial

  linear_objective

    Array representing the linear objective function coefficients based on the weights of subsets.

    Type:

      np.ndarray

  quad_objective

    Quadratic objective function matrix initialized as zeros (no quadratic terms for objective).

    Type:

      np.ndarray

  constraints

    Matrix A and vector b representing constraints for the set partition problem: - A: Binary matrix indicating subset membership of elements in the universal set. - b: Vector of ones, enforcing full coverage of the universal set by the selected subsets.

    Type:

      tuple of np.ndarray

  universal_set

    Set containing all unique elements present across the input subsets, representing the elements that must be fully covered in the partition solution.

    Type:

      set

  Example

  Given a list of subsets representing a set partition problem, each with an associated weight:

      >>> import numpy as np
      >>> np.random.seed(21)
      >>> subsets = [{"A", "B", "C"}, {"D", "E", "F"}, {"A", "F"}, {"B", "E"}, {"C", "D"}, {"A"},
                     {"B"}, {"C", "D", "E"}, {"B", "C"}]
      >>> weights = [100 * np.random.random() for _ in subsets]

  We can construct and use the SetPartitionModel as follows:

      >>> model = SetPartitionModel(subsets=subsets, weights=weights)
      >>> solution = np.random.randint(0, 2, len(subsets))  # Random binary solution vector
      >>> print("Objective Value:", model.evaluateObjective(solution))  # Evaluate solution cost

  This approach builds the constraints matrix and penalties automatically, enabling efficient optimization using solvers like Dirac3CloudSolver.

  evaluateObjective(solution: ndarray) → float

    Evaluate the objective function by calculating the weighted sum of selected subsets.

    Parameters:

      solution (np.ndarray) – Binary array where each element indicates if a subset is selected (1) or not (0).

    Returns:

      The value of the objective function, representing the total weight of selected subsets.

    Return type:

      float

# eqc_models.graph

class eqc_models.graph.MaxCutModel(G: Graph)

  Bases: `TwoPartitionModel`

  decode(solution: ndarray) → ndarray

    Override the default decoding to use a the max cut metric to determine a solution

  property J: ndarray
  property C: ndarray
  property H

    Hamiltonian operator of unknown type

  partition(solution)

    Return a dictionary with the partition number of each node

  getCutSize(partition)
  costFunction()
    Parameters:

      None

    Returns:

      - C: linear operator (vector array of coefficients) for cost function

      - J: quadratic operator (N by N matrix array of coefficients ) for cost function

class eqc_models.graph.GraphPartitionModel(G: Graph, k: int = 2, weight: str = 'weight', alpha: float = 1.0, beta_obj: float = 1.0, gamma: float = 1.0)

  Bases: `GraphModel`

  A model for graph partitioning into k parts with objective and constraints derived from the Laplacian matrix and additional penalties for balance and constraints.

  set_and_validate_k()

    Sets k and encoding length for a graph problem

  initialize_model()

    Build the objective matrix and constraints for the k-partition problem.

  get_balanced_partition_term() → spmatrix

    Construct the quadratic penalty term for balanced partitions.

  get_constraints() → spmatrix

    Construct the quadratic penalty term for assignment constraints.

  get_two_partition_qubo() → spmatrix

    Construct the QUBO matrix for two partitions using adjacency and penalties.

  evaluate(solution: ndarray) → float

    Evaluate the objective function for a given solution.

  decode(solution: ndarray) → dict

    Decode the solution vector into a partition assignment.

  costFunction() → Tuple[ndarray, ndarray]

    Return the linear and quadratic components of the objective function.

class eqc_models.graph.EdgeMixin

  Bases: `object`

  property variables: List[str]

    Provide a variable name to index lookup; order enforced by sorting the list before returning

class eqc_models.graph.EdgeModel(G: Graph)

  Bases: `EdgeMixin`, `GraphModel`

  Create a model where the variables are edge-based

class eqc_models.graph.GraphModel(G: Graph)

  Bases: `QuadraticModel`

  property linear_objective

    Return linear terms as a vector.

  property quad_objective

    Return quadratic terms as a matrix.

  costFunction()
    Parameters:

      None

    Returns:

      - C: linear operator (vector array of coefficients) for cost function

      - J: quadratic operator (N by N matrix array of coefficients ) for cost function

class eqc_models.graph.NodeModel(G: Graph)

  Bases: `GraphModel`

  Base class for a model where the decision variables correspond to the graph nodes.

  property variables: List[str]

    Provide a variable name to index lookup; order enforced by sorting the list before returning

  modularity(partition: Set[Set]) → float

    Calculate modularity from a partition (set of communities)

# eqc_models.sequence

class eqc_models.sequence.MTZTSPModel(D: Dict[Tuple[int, int], float])

  Bases: `InequalitiesMixin`, `TSPModel`

  Using the Miller Tucker Zemlin (1960) formulation of TSP, create a model instance that contains variables for node order (to eliminate subtours) and edge choices.

      >>> D = {(1, 2): 1, (2, 1): 1, (1, 3): 2, (3, 1): 2, (2, 3): 3, (3, 2): 3}
      >>> model = MTZTSPModel(D)
      >>> model.penalty_multiplier = 10
      >>> model.distance(1, 2)
      1
      >>> model.distance(3, 2)
      3
      >>> solution = np.array([1, 0, 0, 1, 1, 0, 1, 2, 3, 2, 2, 2, 5])
      >>> model.cost(solution)
      6
      >>> lhs, rhs = model.constraints
      >>> lhs.shape
      (10, 13)
      >>> model.senses
      ['EQ', 'EQ', 'EQ', 'EQ', 'EQ', 'EQ', 'LE', 'LE', 'LE', 'LE']
      >>> rhs.shape
      (10,)
      >>> (lhs@solution - rhs == 0).all()
      True
      >>> model.upper_bound
      array([1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3])
      >>> poly = model.polynomial
      >>> poly.evaluate(solution) + model.penalty_multiplier * model.offset
      6.0
      >>> infeasible = np.array([0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 4, 2, 5])
      >>> Pl, Pq = -2 * rhs.T@lhs, lhs.T@lhs
      >>> Pl.T@solution + solution.T@Pq@solution + model.offset
      0.0
      >>> Pl.T@infeasible + infeasible.T@Pq@infeasible + model.offset > 0
      True
      >>> poly.evaluate(infeasible) + model.penalty_multiplier * model.offset > 0
      True
      >>> infeasible = np.array([1, 1, 0, 1, 0, 0, 0, 1, 2, 2, 2, 0, 3])
      >>> poly.evaluate(infeasible) + model.penalty_multiplier * model.offset > 6
      True
      >>> val1 = poly.evaluate(infeasible) + model.penalty_multiplier * model.offset
      >>> model.penalty_multiplier *= 2
      >>> poly2 = model.polynomial
      >>> val2 = poly2.evaluate(infeasible) + model.penalty_multiplier * model.offset
      >>> val1 < val2
      True

  cost(solution: ndarray) → float

    Solution cost is the sum of all D values where (i,j) is chosen

    :paramsolution: np.ndarray - An array of of 0,1 values which describe a route

      by setting variables representing active edges to 1.

    Returns: float

  property upper_bound: ndarray

    For all route variables, the domain is {0, 1}. The sequence variables can take on values in [0, 1, 2, …, N-1].

# eqc_models.algorithms

class eqc_models.algorithms.PenaltyMultiplierAlgorithm(model: ConstraintModel, solver: ModelSolver, penalty_threshold: float = 1e-06, progress_threshold: float = 0.0001)

  Bases: `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_boundfloat

    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

  solutionsList

    The solutions found during the algorithm run

  alphasList

    The values of multiplier found at each algorithm iteration

  penaltiesList

    The values for penalties found at each algorithm iteration. A penalty of 0 indicates algorithm termination.

  penalty_thresholdfloat

    This value is the cutoff for penalty values that are threated as 0. Default is 1e-6.

  progress_thresholdfloat

    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_rangeList

    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)
      2... RUNNING... COMPLETED...
      >>> algo.alphas[-1]
      106.25
      >>> algo.penalties[-1]
      0.0

  property upper_bound: float
  run(initial_alpha: float | None = None, initial_solution: array | None = None, **kwargs)

    Start with a guess at alpha, iterate until alpha is sufficiently large

class eqc_models.algorithms.ALMAlgorithm

  Bases: `object`

  Stateless ALM outer loop. Call run(model, registry, core, cfg, **core_kwargs).

  static lift_Qc_to_poly_terms(Q_native: ndarray, c_native: ndarray, blocks: List[ALMBlock]) → Tuple[List[Tuple[int, ...]], List[float]]

    Expand a quadratic/linear objective over m native discrete variables into the lifted (one-hot) s-space.

    Given Q_native and c_native over original m discrete vars (one block per var), and blocks (list[ALMBlock] in the same order as original vars), we enforce: - For variable i with level values b_i (length k_i), we create k_i lifted coords s_{i,a}. - Quadratic term: J_ij = Q_native[i,j] * (b_i b_j^T) contributes to pairs (s_{i,a}, s_{j,b}) - Linear term: C_i = c_native[i] * b_i contributes to s_{i,a}

    Returns polynomial (indices, coeffs) over concatenated s variables.

  static run(base_model: PolynomialModel, registry: ConstraintRegistry, solver, cfg: ALMConfig = ALMConfig(rho_h=1.0, rho_g=1.0, rho_min=0.001, rho_max=1000.0, rho_eq_family={}, rho_ineq_family={}, adapt_by_family=True, adapt=True, tau_up_h=1.05, tau_down_h=0.95, tau_up_g=0.9, tau_down_g=0.5, gamma_up=1.2, gamma_down=0.9, tol_h=1e-06, tol_g=1e-06, max_outer=100, use_stagnation_bump=False, patience_h=10, patience_g=10, stagnation_factor=0.001, ema_alpha=1.0, fd_eps=1e-06, act_tol=1e-10, bcd=BCDConfig(enabled=False, max_sweeps=25, batch_blocks=1, shuffle_blocks=False, repeats_per_block=1, tolerance=1e-08, patience=5, proximal_strength=0.0, warm_start=True, clip_min=None, clip_max=None, random_seed=None)), x0: ndarray | None = None, *, parse_output=None, verbose: bool = True, **solver_kwargs) → Dict[str, ndarray | Dict[int, float] | Dict]

    Solve with ALM. Keep all ALM state local to this call (no global side-effects). Handles three modes:

    > 1. No blocks -> continuous; use base_model as-is

    > 2. Blocks + native base_model -> lift to s-space

    > 3. Blocks + already-lifted base_model -> use as-is (compat)

    Returns:

      {

        “x”: final iterate, “decoded”: {start_idx_of_block: level_value, …} for lifted blocks, “decoded_debug”: {start_idx_of_block: native_device_value, …}, “hist”: { “eq_inf”: […], “ineq_inf”: […], “obj”: […], “x”: […] }

      }

class eqc_models.algorithms.ConstraintRegistry

  Bases: `object`

  Holds constraints and block metadata; keeps ALMAlgorithm stateless. Register constraints and (optional) lifted-discrete blocks here.

  add_equality(fun, jac=None, name='', family='')
  add_inequality(fun, jac=None, name='', family='')
  add_block(idx: Sequence[int], levels: ndarray, sum_to_one=True, one_hot=True)

class eqc_models.algorithms.ALMConfig(rho_h: float = 1.0, rho_g: float = 1.0, rho_min: float = 0.001, rho_max: float = 1000.0, rho_eq_family: Dict[str, float] = <factory>, rho_ineq_family: Dict[str, float] = <factory>, adapt_by_family: bool = True, adapt: bool = True, tau_up_h: float = 1.05, tau_down_h: float = 0.95, tau_up_g: float = 0.9, tau_down_g: float = 0.5, gamma_up: float = 1.2, gamma_down: float = 0.9, tol_h: float = 1e-06, tol_g: float = 1e-06, max_outer: int = 100, use_stagnation_bump: bool = False, patience_h: int = 10, patience_g: int = 10, stagnation_factor: float = 0.001, ema_alpha: float = 1.0, fd_eps: float = 1e-06, act_tol: float = 1e-10, bcd: eqc_models.algorithms.bcd.BCDConfig = <factory>)

  Bases: `object`

  rho_h: float = 1.0
  rho_g: float = 1.0
  rho_min: float = 0.001
  rho_max: float = 1000.0
  rho_eq_family: Dict[str, float]
  rho_ineq_family: Dict[str, float]
  adapt_by_family: bool = True
  adapt: bool = True
  tau_up_h: float = 1.05
  tau_down_h: float = 0.95
  tau_up_g: float = 0.9
  tau_down_g: float = 0.5
  gamma_up: float = 1.2
  gamma_down: float = 0.9
  tol_h: float = 1e-06
  tol_g: float = 1e-06
  max_outer: int = 100
  use_stagnation_bump: bool = False
  patience_h: int = 10
  patience_g: int = 10
  stagnation_factor: float = 0.001
  ema_alpha: float = 1.0
  fd_eps: float = 1e-06
  act_tol: float = 1e-10
  bcd: BCDConfig

class eqc_models.algorithms.ALMConstraint(kind: str, fun: Callable[[ndarray], ndarray], jac: Callable[[ndarray], ndarray] | None = None, name: str = '', family: str = '')

  Bases: `object`

  One constraint family; fun returns a vector; jac returns its Jacobian.

  kind: str
  fun: Callable[[ndarray], ndarray]
  jac: Callable[[ndarray], ndarray] | None = None
  name: str = ''
  family: str = ''

class eqc_models.algorithms.ALMBlock(idx: Sequence[int], levels: ndarray, enforce_sum_to_one: bool = True, enforce_one_hot: bool = True)

  Bases: `object`

  Lifted discrete variable block (optional).

  idx: Sequence[int]
  levels: ndarray
  enforce_sum_to_one: bool = True
  enforce_one_hot: bool = True

# eqc_models.algorithms.alm

## Augmented Lagrangian Method

### Overview

The augmented Lagrangian method algorithm in eqc-models is for automating parameter adjustments to achieve feasible solutions to constrained optimization models using entropy quantum computing devices. It supports fractional (continuous), binary and general integer variables. ALM also supports enforcement of one-hot constraints, inequality constraints and equality constraints without introduciton of slack variables.

This method requires a different approach to modeling than used in model implementations of constrained optimization models as direct subclasses of EqcModel. To use this approach, implement an objective function in one of the subclasses of EqcModel. The next step is to use a ConstraintRegistry object to describe the model constraints. There are many parameters available inside the ALMConfig class. Only a few of these parameters should be adjusted by novice users. Once the model, constraint registry and config objects are ready, pass them to an ALMAlgorithm instance with a solver and it will iterate over parameter adjustments until feasibility tolerance is found or maximum iteration count is exceeded.

The following example is a basic outline of how to use the ALM implementation to solve models.

    model = PolynomialModel()
    solver = Dirac3ContinuousCloudSolver()
    solver_kwargs = {} # add solver parameters as dict here
    config = ALMConfig() # add configuration parameters here
    registry = ConstraintRegistry()
    # add variable domains and constraints here
    algorithm = ALMAlgorithm()
    response = algorithm.run(model, registry, solver, config, **solver_kwargs) # keyword arguments are passed to the solver

### Model Creation

Definining the objective function for the optimization problem can be done in a number of ways. There are two base classes on which all other models are based and they define how function coefficients are specified. With a QuadraticModel, the coefficients are specified using a one-dimensional array for the linear terms and a two-dimensional, symmetric array or Hermitian matrix for quadratic terms. The PolynomialModel way uses two lists, one for indices of terms and another for the coefficients of terms. These terms can be up to degree five. The two classes just described are preferred for use with ALM because most other model classes incorporate constraints into the operator. In the future, ALM will provide a way to extract the objective function and constraints from the model separately, populating the ConstraintRegistry automatically. First, choose variable definitions by describing the variable names in a single list. Then organize the terms in ascending lexicographic order by index. Finally, create a list with coefficient entries corresponding to the index list. The example problem used here is a simple logical decision of one variable among 10. Those familiar with QUBO formulations will recognize the penalties in the quadratic coefficients. This first example sets the one-hot constraint as penalties in the objective function and adds block constraints to each of the variables to enforce 0 or 1. This has the effect of using two time bins per binary variable.

    >>> import warnings; warnings.filterwarnings("ignore") # silence warnings
    >>> n = 10
    >>> variables = [f"x_{i}" for i in range(n)]

Define the indices and coefficients. The solution is x_1 == 1 and all other x_i == 0. Variable indices sent to the device are 1-based.

    >>> indices = [(0, i+1) for i in range(n)] + [(i+1, j+1) for i in range(10) for j in range(i, n)]
    >>> coefficients = [-i if i==n-1 else 0 for i in range(n)] + [2 for i in range(n) for j in range(i, n)]

Define the model, setting the upper bound to 1 on all varaibles.

    >>> model = PolynomialModel(coefficients, indices)
    >>> model.variables = variables
    >>> model.variables[0] == "x_0"
    True
    >>> model.upper_bound = np.ones((n,))
    >>> model.is_discrete = [True for i in range(n)] # for future modeling support

Begin the ALM setup. One of the features of this algorithm in eqc-models is the ability to enforce discrete variables while using the continuous solver. This can allow for mixed integer formulations. For binary variables, add a block for each one so that it is forced to 0 or 1. The sum_to_one parameter is normally true and setting to false here is done because the simplicity of the problem allows a more efficient approach.

    >>> reg = ConstraintRegistry()
    >>> for i in range(n): reg.add_block([i], levels=[0, 1], sum_to_one=False)

Configure rho_h (equality constraints) to start at 1 and only allow 4 iterations. Also set tolerance high.

    >>> config = ALMConfig(tol_h=0.5, max_outer=4, rho_h=1, rho_g=1, gamma_up=2.0, gamma_down=1.0)
    >>> from eqc_models.solvers import Dirac3ContinuousCloudSolver
    >>> solver = Dirac3ContinuousCloudSolver()
    >>> algo = ALMAlgorithm()
    >>> results = algo.run(model, reg, solver, config, sum_constraint=n, relaxation_schedule=1)
     20...COMPLETE...
    [ALM 00] ...
    ...
    >>> results["decoded"] # had to hard-code to n==10
    {0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 0.0, 8: 0.0, 9: 1.0}

### Equality constraints

This ALM algorithm can seek feasible solutions with targeted multipliers rather than adding penalty multipliers into the objective coefficients and adjusting manually, which would be required for more complex problems. These constraints should be specified in the constraint registry. Specification of the binary variables using the add_block call sets up the domain. The objective function is formulated directly, purely as the cost of each variable, which is the negative value of the variable index.

    >>> indices = [(0, i+1) for i in range(n)]
    >>> coefficients = [-i for i in range(n)]
    >>> model = PolynomialModel(coefficients, indices)
    >>> model.upper_bound = np.ones((n, ))
    >>> model.is_discrete = [True for i in range(n)]
    >>> reg = ConstraintRegistry()
    >>> for i in range(n): reg.add_block([i], levels=[0, 1], sum_to_one=False)

Now, a new concept. Instead of placing the penalties explicitly in the objective, this statement enforces the one-hot constraint by defining a function which is 0 when true and either positive or negative when false. The lambda function is passed an encoded array for the variable values, so the array must be decoded before evaluating the function. To assist in this, a helper function called level_decoder is available for import. Use this decoder and the solution x to extract the values of the variables.

    >>> decoder = level_decoder(n, [0, 1])
    >>> test_x = np.zeros((20,))
    >>> test_x[::2] = 1
    >>> (decoder@test_x==0).all()
    True
    >>> config.max_outer = 15 # added complexity in the definition requires more steps
    >>> config.gamma_up = 1.2
    >>> config.gamma_down = 0.9
    >>> reg.add_equality(lambda x: np.sum(decoder@x) - 1)

Everything else is the same as in the previous example.

    >>> results = algo.run(model, reg, solver, config, sum_constraint=n, relaxation_schedule=1)
     20...COMPLETE...
    [ALM 00] ...
    ...
    >>> results["decoded"] # had to hard-code to n==10
    {0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 0.0, 8: 0.0, 9: 1.0}

### Inequality Constraints

Incorporating inequality constraints is not much different. Specify an inequality constraint by using the ConstraintRegistry.add_inequality method to pass a function which is zero or negative when the constraint is met and positive when the constraint is not met. This is straightforward for less than or equal to constraints and requires negating the function when the constraint is for greater than or equal. Here is an example showing a model where less than two variables are non zero. When adding the binary variable blocks, the sum_to_one parameter is left out, which results in the default of true. The use of the sum to one parameter gives the algorithm more control over the dual variable computations.

    >>> reg = ConstraintRegistry()
    >>> for i in range(n): reg.add_block([i], levels=[0, 1])
    >>> reg.add_inequality(lambda x: np.sum(decoder@x) - 2) # reuse the decoder from the previous example
    >>> config.gamma_up = 1.2 # slower adjustment
    >>> config.gamma_down = 0.8 # slower adjustment
    >>> results = algo.run(model, reg, solver, config, sum_constraint=n, relaxation_schedule=1)
     20...COMPLETED...
     [ALM 00] ...
     [ALM 05] ...
     ...
    >>> results["decoded"]
    {0: 0.0, 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0, 5: 0.0, 6: 0.0, 7: 0.0, 8: 1.0, 9: 1.0}

eqc_models.algorithms.alm.level_decoder(n, levels)

  Return an operator that transforms an encoded solution to a solution with a single value per variable. For instance, a sequence of two level variables with levels 0 and 1 has values [0, 1, 1, 0] is decoded into [1, 0] with the operator:

      [[ 0, 1],
       [ 0, 1]]

  and used with the expression op@variable.

      >>> decoder = level_decoder(3, [0, 1, 2])
      >>> values = np.array([1, 0, 0, 0, 1, 0, 0, 0, 1])
      >>> decoder@values
      array([0., 1., 2.])

class eqc_models.algorithms.alm.ALMConstraint(kind: str, fun: Callable[[ndarray], ndarray], jac: Callable[[ndarray], ndarray] | None = None, name: str = '', family: str = '')

  Bases: `object`

  One constraint family; fun returns a vector; jac returns its Jacobian.

  kind: str
  fun: Callable[[ndarray], ndarray]
  jac: Callable[[ndarray], ndarray] | None = None
  name: str = ''
  family: str = ''

class eqc_models.algorithms.alm.ALMBlock(idx: Sequence[int], levels: ndarray, enforce_sum_to_one: bool = True, enforce_one_hot: bool = True)

  Bases: `object`

  Lifted discrete variable block (optional).

  idx: Sequence[int]
  levels: ndarray
  enforce_sum_to_one: bool = True
  enforce_one_hot: bool = True

class eqc_models.algorithms.alm.ALMConfig(rho_h: float = 1.0, rho_g: float = 1.0, rho_min: float = 0.001, rho_max: float = 1000.0, rho_eq_family: Dict[str, float] = <factory>, rho_ineq_family: Dict[str, float] = <factory>, adapt_by_family: bool = True, adapt: bool = True, tau_up_h: float = 1.05, tau_down_h: float = 0.95, tau_up_g: float = 0.9, tau_down_g: float = 0.5, gamma_up: float = 1.2, gamma_down: float = 0.9, tol_h: float = 1e-06, tol_g: float = 1e-06, max_outer: int = 100, use_stagnation_bump: bool = False, patience_h: int = 10, patience_g: int = 10, stagnation_factor: float = 0.001, ema_alpha: float = 1.0, fd_eps: float = 1e-06, act_tol: float = 1e-10, bcd: eqc_models.algorithms.bcd.BCDConfig = <factory>)

  Bases: `object`

  rho_h: float = 1.0
  rho_g: float = 1.0
  rho_min: float = 0.001
  rho_max: float = 1000.0
  rho_eq_family: Dict[str, float]
  rho_ineq_family: Dict[str, float]
  adapt_by_family: bool = True
  adapt: bool = True
  tau_up_h: float = 1.05
  tau_down_h: float = 0.95
  tau_up_g: float = 0.9
  tau_down_g: float = 0.5
  gamma_up: float = 1.2
  gamma_down: float = 0.9
  tol_h: float = 1e-06
  tol_g: float = 1e-06
  max_outer: int = 100
  use_stagnation_bump: bool = False
  patience_h: int = 10
  patience_g: int = 10
  stagnation_factor: float = 0.001
  ema_alpha: float = 1.0
  fd_eps: float = 1e-06
  act_tol: float = 1e-10
  bcd: BCDConfig

class eqc_models.algorithms.alm.ConstraintRegistry

  Bases: `object`

  Holds constraints and block metadata; keeps ALMAlgorithm stateless. Register constraints and (optional) lifted-discrete blocks here.

  add_equality(fun, jac=None, name='', family='')
  add_inequality(fun, jac=None, name='', family='')
  add_block(idx: Sequence[int], levels: ndarray, sum_to_one=True, one_hot=True)

class eqc_models.algorithms.alm.ALMAlgorithm

  Bases: `object`

  Stateless ALM outer loop. Call run(model, registry, core, cfg, **core_kwargs).

  static lift_Qc_to_poly_terms(Q_native: ndarray, c_native: ndarray, blocks: List[ALMBlock]) → Tuple[List[Tuple[int, ...]], List[float]]

    Expand a quadratic/linear objective over m native discrete variables into the lifted (one-hot) s-space.

    Given Q_native and c_native over original m discrete vars (one block per var), and blocks (list[ALMBlock] in the same order as original vars), we enforce: - For variable i with level values b_i (length k_i), we create k_i lifted coords s_{i,a}. - Quadratic term: J_ij = Q_native[i,j] * (b_i b_j^T) contributes to pairs (s_{i,a}, s_{j,b}) - Linear term: C_i = c_native[i] * b_i contributes to s_{i,a}

    Returns polynomial (indices, coeffs) over concatenated s variables.

  static run(base_model: PolynomialModel, registry: ConstraintRegistry, solver, cfg: ALMConfig = ALMConfig(rho_h=1.0, rho_g=1.0, rho_min=0.001, rho_max=1000.0, rho_eq_family={}, rho_ineq_family={}, adapt_by_family=True, adapt=True, tau_up_h=1.05, tau_down_h=0.95, tau_up_g=0.9, tau_down_g=0.5, gamma_up=1.2, gamma_down=0.9, tol_h=1e-06, tol_g=1e-06, max_outer=100, use_stagnation_bump=False, patience_h=10, patience_g=10, stagnation_factor=0.001, ema_alpha=1.0, fd_eps=1e-06, act_tol=1e-10, bcd=BCDConfig(enabled=False, max_sweeps=25, batch_blocks=1, shuffle_blocks=False, repeats_per_block=1, tolerance=1e-08, patience=5, proximal_strength=0.0, warm_start=True, clip_min=None, clip_max=None, random_seed=None)), x0: ndarray | None = None, *, parse_output=None, verbose: bool = True, **solver_kwargs) → Dict[str, ndarray | Dict[int, float] | Dict]

    Solve with ALM. Keep all ALM state local to this call (no global side-effects). Handles three modes:

    > 1. No blocks -> continuous; use base_model as-is

    > 2. Blocks + native base_model -> lift to s-space

    > 3. Blocks + already-lifted base_model -> use as-is (compat)

    Returns:

      {

        “x”: final iterate, “decoded”: {start_idx_of_block: level_value, …} for lifted blocks, “decoded_debug”: {start_idx_of_block: native_device_value, …}, “hist”: { “eq_inf”: […], “ineq_inf”: […], “obj”: […], “x”: […] }

      }

# eqc_models.utilities

eqc_models.utilities.evaluate_polynomial(terms, solution)

eqc_models.utilities.read_coefficient_file(filename)

eqc_models.utilities.read_index_file(filename)

eqc_models.utilities.read_config_file(filename)

eqc_models.utilities.convert_hamiltonian_to_polynomial(A, B, C, D, num_vars)

  Converts a hamiltonian of up to fourth order to a polynomial.

  D_{ijkl} x_i x_j x_k x_l + C_{ijk} x_i x_j x_k + B_{ij} x_i x_j + A_i x_i

  Input:

  A: First order hamiltonian. B: Second order hamiltonian. C: Third order hamiltonian. D: Fourth order hamiltonian. num_vars: Number of variables.

  Output:

  polynomial_indices, polynomila_coefs: Indices and coefficients of polynomial that respresents the hamiltonian.

eqc_models.utilities.create_json_problem(A: array, B: array, C: array, D: array, num_vars: int, sum_constraint: float | None = None, num_levels: int | None = None)

  Converts a hamiltonian of up to fourth order to a polynomial.

  D_{ijkl} x_i x_j x_k x_l + C_{ijk} x_i x_j x_k + B_{ij} x_i x_j + A_i x_i

  Input:

  A: First order hamiltonian. B: Second order hamiltonian. C: Third order hamiltonian. D: Fourth order hamiltonian. num_vars: Number of variables.

  Output:

  Problem in json format.

# eqc_models.utilities.qplib

Read the file format from QPLIB (https://doi.org/10.1007/s12532-018-0147-4)

The “problem type” is a string of three characters.

The first character indicates the type of objective function used. It must be one of the following:

> L a linear objective function D a convex quadratic objective function whose Hessian is a diagonal matrix C a convex quadratic objective function Q a quadratic objective function whose Hessian may be indefinite

The second character indicates the types of variables that are present. It must be one of the following:

> C all the variables are continuous B all the variables are binary (0-1) M the variables are a mix of continuous and binary I all the variables are integer G the variables are a mix of continuous, binary and integer

The third character indicates the type of the (most extreme) constraint function used; other constraints may be of a lesser type. It must be one of the following:

> N there are no constraints B some of the variables lie between lower and upper bounds (box constraints) L the constraint functions are linear D the constraint functions are convex quadratics with diagonal Hessians C the constraint functions are convex quadratics Q the constraint functions are quadratics whose Hessians may be indefinite

Thus for continuous problems, we would have

> LCL a linear program LCC or LCQ a linear program with quadratic constraints CCB or QCB a bound-constrained quadratic program CCL or QCL a quadratic program CCC or CCQ or a quadratic program with quadratic constraints QCC or QCQ

class eqc_models.utilities.qplib.QBBModel(C: ndarray, J: ndarray)

  Bases: `QuadraticModel`

  QBB - Quadratic objective, binary variable type, unconstrained

  This model type is the same as a QUBO model, but the linear portion is specified separately, instead of in the diagonal.

  :param : :type : C: np.ndarray - linear portion of the objective :param : :type : J: np.ndarray - quadratic portion of the objective

class eqc_models.utilities.qplib.QGLModel(C_obj: array, J_obj: array, lhs: array, rhs: array)

  Bases: `ConstrainedQuadraticModel`

  QGL - Quadratic objective, General variable type, Linear constraints Handles QCL, QIL, QBL, QML problems

  :param : :type : C: np.ndarray - linear portion of the objective :param : :type : J: np.ndarray - quadratic portion of the objective :param : :type : A: np.ndarray - left hand side of the linear constraints :param : :type : b: np.ndarray - right hand side of the linear constraints

eqc_models.utilities.qplib.file_to_model(fp: ~typing.TextIO, DTYPE=<class 'numpy.float32'>) → QGLModel

eqc_models.utilities.qplib.file_to_polynomial(fp: ~typing.TextIO, DTYPE=<class 'numpy.float32'>, penalty_multiplier=1) → Tuple[List, List]

  Create a pair of lists describing a polynomial that represents the file contained in the qplib-formatted file descriptor fp. Sets a guess for the sum constraint when constructing the model, but it has no effect on the polynomial file output.

  Parameters :fp: File descriptor for qplib-formatted file :DTYPE: Default numeric datatype for non-integer numeric values

  Returns
    coefficients:, :indices:

eqc_models.utilities.qplib.process_file(fp, DTYPE=<class 'numpy.float32'>)

  Read the file line by line, saving the parts objective (C, Q) and constraints (A, b)

eqc_models.utilities.qplib.read_line(fp: TextIO) → Dict