Implementing an OMMX Adapter

Implementing an OMMX Adapter#

As mentioned in Solve with multiple adapters and compare the results, OMMX Adapters have a common API. This common API is realized by inheriting the abstract base classes provided by the OMMX Python SDK. OMMX provides two abstract base classes depending on the type of adapter:

  • ommx.adapter.SolverAdapter: An abstract base class for optimization solvers that return one solution

  • ommx.adapter.SamplerAdapter: An abstract base class for sampling-based optimization solvers

Solvers that produce multiple solutions can be automatically treated as solvers returning a single solution by selecting the best sample. Therefore, SamplerAdapter inherits SolverAdapter. If you are unsure which one to implement, consider the number of solutions: if the solver returns one solution, use SolverAdapter; if it returns multiple solutions, use SamplerAdapter. For example, exact solvers like PySCIPOpt should use SolverAdapter, while samplers like OpenJij should use SamplerAdapter.

In OMMX, a class inheriting ommx.adapter.SolverAdapter is called a Solver Adapter and one inheriting ommx.adapter.SamplerAdapter is called a Sampler Adapter. For clear explaination in this chapter, the software that the adapter wraps (such as PySCIPOpt or OpenJij) is referred as “backend solver”.

Adapter Workflow#

The adapter process can be roughly divided into these 3 steps:

  1. Convert ommx.v1.Instance into a format the backend solver can understand

  2. Run the backend solver to obtain a solution

  3. Convert the backend solver’s output into ommx.v1.Solution or ommx.v1.SampleSet

Because the step 2 is nothing but the usage of the backend solver, we assume you to known it well. This tutorial explains steps 1 and 3.

Many backend solvers are designed to receive only the minimum necessary information to represent an optimization problem in a form suitable for their algorithms, whereas ommx.v1.Instance contains more information, assuming optimization as part of data analysis. Therefore, step 1 involves discarding much of this information. Additionally, OMMX manages decision variables and constraints with IDs that are not necessarily sequential, while some backend solvers manage them by names or sequential numbers. This correspondence is needed in step 3, so the adapter must manage it.

Conversely, in step 3, ommx.v1.Solution or ommx.v1.SampleSet, because these stores information same as ommx.v1.Instance, cannot be constructed solely from the backend solver’s output. Instead, the adapter will construct ommx.v1.State or ommx.v1.Samples from the backend solver’s output and the information from step 1, then convert it to ommx.v1.Solution or ommx.v1.SampleSet using ommx.v1.Instance.

Implementing a Solver Adapter#

Here, we will implement a Solver Adapter using PySCIPOpt as an example. For a complete example, refer to ommx-pyscipopt-adapter.

For this tutorial, we will proceed in the following order to make it easier to execute step by step:

  • Implement functions to construct a PySCIPOpt model from ommx.v1.Instance one by one.

  • Finally, combine these functions into the OMMXPySCIPOptAdapter class.

Custom Exception#

First, it is good to define custom exceptions. This makes it easier for users to understand which part is causing the problem when an exception occurs.

class OMMXPySCIPOptAdapterError(Exception):
    pass

OMMX can store a wide range of optimization problems, so there may be cases where the backend solver does not support the problem. In such cases, throw an error.

Setting Decision Variables#

PySCIPOpt manages decision variables by name, so register the OMMX decision variable IDs as strings. This allows you to reconstruct ommx.v1.State from PySCIPOpt decision variables in the decode_to_state function mentioned later. Note that the appropriate method depends on the backend solver’s implementation. The important thing is to retain the information needed to convert to ommx.v1.State after obtaining the solution.

import pyscipopt
from ommx.v1 import Instance, Solution, DecisionVariable, Constraint, State, Optimality, Function

def set_decision_variables(model: pyscipopt.Model, instance: Instance) -> dict[str, pyscipopt.Variable]:
    """Add decision variables to the model and create a mapping from variable names to variables"""
    # Create PySCIPOpt variables from OMMX decision variable information
    for var in instance.raw.decision_variables:
        if var.kind == DecisionVariable.BINARY:
            model.addVar(name=str(var.id), vtype="B")
        elif var.kind == DecisionVariable.INTEGER:
            model.addVar(
                name=str(var.id), vtype="I", lb=var.bound.lower, ub=var.bound.upper
            )
        elif var.kind == DecisionVariable.CONTINUOUS:
            model.addVar(
                name=str(var.id), vtype="C", lb=var.bound.lower, ub=var.bound.upper
            )
        else:
            # Throw an error if an unsupported decision variable type is encountered
            raise OMMXPySCIPOptAdapterError(
                f"Unsupported decision variable kind: id: {var.id}, kind: {var.kind}"
            )

    # If the objective is quadratic, add an auxiliary variable for linearization
    if instance.raw.objective.HasField("quadratic"):
        model.addVar(
            name="auxiliary_for_linearized_objective", vtype="C", lb=None, ub=None
        )

    # Create a dictionary to access the variables added to the model
    return {var.name: var for var in model.getVars()}

Converting ommx.v1.Function to pyscipopt.Expr#

Implement a function to convert ommx.v1.Function to pyscipopt.Expr. Since ommx.v1.Function only has the OMMX decision variable IDs, you need to obtain the PySCIPOpt variables from the IDs using the variable name and variable mapping created in set_decision_variables.

def make_linear_expr(function: Function, varname_map: dict) -> pyscipopt.Expr:
    """Helper function to generate a linear expression"""
    linear = function.linear
    return (
        pyscipopt.quicksum(
            term.coefficient * varname_map[str(term.id)]
            for term in linear.terms
        )
        + linear.constant
    )

def make_quadratic_expr(function: Function, varname_map: dict) -> pyscipopt.Expr:
    """Helper function to generate a quadratic expression"""
    quad = function.quadratic
    quad_terms = pyscipopt.quicksum(
        varname_map[str(row)] * varname_map[str(column)] * value
        for row, column, value in zip(quad.rows, quad.columns, quad.values)
    )

    linear_terms = pyscipopt.quicksum(
        term.coefficient * varname_map[str(term.id)]
        for term in quad.linear.terms
    )

    constant = quad.linear.constant

    return quad_terms + linear_terms + constant

Setting Objective Function and Constraints#

Add the objective function and constraints to the pyscipopt.Model. This part requires knowledge of what and how the backend solver supports. For example, in the following code, since PySCIPOpt cannot directly handle quadratic objective functions, an auxiliary variable is introduced according to the PySCIPOpt documentation.

import math

def set_objective(model: pyscipopt.Model, instance: Instance, varname_map: dict):
    """Set the objective function for the model"""
    objective = instance.raw.objective

    if instance.sense == Instance.MAXIMIZE:
        sense = "maximize"
    elif instance.sense == Instance.MINIMIZE:
        sense = "minimize"
    else:
        raise OMMXPySCIPOptAdapterError(
            f"Sense not supported: {instance.sense}"
        )

    if objective.HasField("constant"):
        model.setObjective(objective.constant, sense=sense)
    elif objective.HasField("linear"):
        expr = make_linear_expr(objective, varname_map)
        model.setObjective(expr, sense=sense)
    elif objective.HasField("quadratic"):
        # Since PySCIPOpt doesn't support quadratic objectives directly, linearize using an auxiliary variable
        auxilary_var = varname_map["auxiliary_for_linearized_objective"]

        # Set the auxiliary variable as the objective
        model.setObjective(auxilary_var, sense=sense)

        # Add a constraint for the auxiliary variable
        expr = make_quadratic_expr(objective, varname_map)
        if sense == "minimize":
            constr_expr = auxilary_var >= expr
        else:  # sense == "maximize"
            constr_expr = auxilary_var <= expr

        model.addCons(constr_expr, name="constraint_for_linearized_objective")
    else:
        raise OMMXPySCIPOptAdapterError(
            "The objective function must be `constant`, `linear`, or `quadratic`."
        )

def set_constraints(model: pyscipopt.Model, instance: Instance, varname_map: dict):
    """Set the constraints for the model"""
    # Process regular constraints
    for constraint in instance.raw.constraints:
        # Generate an expression based on the type of constraint function
        if constraint.function.HasField("linear"):
            expr = make_linear_expr(constraint.function, varname_map)
        elif constraint.function.HasField("quadratic"):
            expr = make_quadratic_expr(constraint.function, varname_map)
        elif constraint.function.HasField("constant"):
            # For constant constraints, check feasibility
            if constraint.equality == Constraint.EQUAL_TO_ZERO and math.isclose(
                constraint.function.constant, 0, abs_tol=1e-6
            ):
                continue
            elif (
                constraint.equality == Constraint.LESS_THAN_OR_EQUAL_TO_ZERO
                and constraint.function.constant <= 1e-6
            ):
                continue
            else:
                raise OMMXPySCIPOptAdapterError(
                    f"Infeasible constant constraint found: id {constraint.id}"
                )
        else:
            raise OMMXPySCIPOptAdapterError(
                f"Constraints must be either `constant`, `linear` or `quadratic`. id: {constraint.id}, type: {constraint.function.WhichOneof('function')}"
            )

        # Add constraints based on the type (equality/inequality)
        if constraint.equality == Constraint.EQUAL_TO_ZERO:
            constr_expr = expr == 0
        elif constraint.equality == Constraint.LESS_THAN_OR_EQUAL_TO_ZERO:
            constr_expr = expr <= 0
        else:
            raise OMMXPySCIPOptAdapterError(
                f"Not supported constraint equality: id: {constraint.id}, equality: {constraint.equality}"
            )

        # Add the constraint to the model
        model.addCons(constr_expr, name=str(constraint.id))

Also, if the backend solver supports special constraints (e.g., SOS constraints), you need to add functions to handle them.

Now, we can construct a pycscipopt.Model from ommx.v1.Instance.

Converting Obtained Solutions to ommx.v1.State#

Next, implement a function to convert the solution obtained by solving the PySCIPOpt model to ommx.v1.State. First, check if it is solved. SCIP has functions to guarantee optimality and detect unbounded solutions, so throw corresponding exceptions if detected. This also depends on the backend solver.

Warning

Note that ommx.adapter.InfeasibleDetected means that the optimization problem itself is infeasible, i.e., it is guaranteed to have no solutions. Do not use this when a heuristic solver fails to find any feasible solutions.

from ommx.adapter import InfeasibleDetected, UnboundedDetected

def decode_to_state(model: pyscipopt.Model, instance: Instance) -> State:
    """Create an ommx.v1.State from an optimized PySCIPOpt Model"""
    if model.getStatus() == "unknown":
        raise OMMXPySCIPOptAdapterError(
            "The model may not be optimized. [status: unknown]"
        )

    if model.getStatus() == "infeasible":
        raise InfeasibleDetected("Model was infeasible")

    if model.getStatus() == "unbounded":
        raise UnboundedDetected("Model was unbounded")

    try:
        # Get the best solution
        sol = model.getBestSol()
        # Create a mapping from variable names to variables
        varname_map = {var.name: var for var in model.getVars()}
        # Create a State with a mapping from variable IDs to their values
        return State(
            entries={
                var.id: sol[varname_map[str(var.id)]]
                for var in instance.raw.decision_variables
            }
        )
    except Exception:
        raise OMMXPySCIPOptAdapterError(
            f"There is no feasible solution. [status: {model.getStatus()}]"
        )

Creating a Class that Inherits ommx.adapter.SolverAdapter#

Finally, create a class that inherits ommx.adapter.SolverAdapter to standardize the API for each adapter. This is an abstract base class with @abstractmethod as follows:

class SolverAdapter(ABC):
    @abstractmethod
    def __init__(self, ommx_instance: Instance):
        pass

    @classmethod
    @abstractmethod
    def solve(cls, ommx_instance: Instance) -> Solution:
        pass

    @property
    @abstractmethod
    def solver_input(self) -> SolverInput:
        pass

    @abstractmethod
    def decode(self, data: SolverOutput) -> Solution:
        pass

This abstract base class assumes the following two use cases:

  • If you do not adjust the backend solver’s parameters, use the solve class method.

  • If you adjust the backend solver’s parameters, use solver_input to get the data structure for the backend solver (in this case, pyscipopt.Model), adjust it, then input it to the backend solver, and finally convert the backend solver’s output using decode.

Using the functions prepared so far, you can implement it as follows:

from ommx.adapter import SolverAdapter

class OMMXPySCIPOptAdapter(SolverAdapter):
    def __init__(self, ommx_instance: Instance):
        self.instance = ommx_instance
        self.model = pyscipopt.Model()
        self.model.hideOutput()
        
        # Build the model with helper functions
        self.varname_map = set_decision_variables(self.model, self.instance)
        set_objective(self.model, self.instance, self.varname_map)
        set_constraints(self.model, self.instance, self.varname_map)

    @classmethod
    def solve(cls, ommx_instance: Instance) -> Solution:
        """Solve an ommx.v1.Instance using PySCIPopt and return an ommx.v1.Solution"""
        adapter = cls(ommx_instance)
        model = adapter.solver_input
        model.optimize()
        return adapter.decode(model)

    @property
    def solver_input(self) -> pyscipopt.Model:
        """Return the generated PySCIPopt model"""
        return self.model

    def decode(self, data: pyscipopt.Model) -> Solution:
        """Generate an ommx.v1.Solution from an optimized pyscipopt.Model and the OMMX Instance"""
        if data.getStatus() == "infeasible":
            raise InfeasibleDetected("Model was infeasible")

        if data.getStatus() == "unbounded":
            raise UnboundedDetected("Model was unbounded")

        # Convert the solution to state
        state = decode_to_state(data, self.instance)
        # Evaluate the state using the instance
        solution = self.instance.evaluate(state)

        # Set the optimality status if the model is optimal
        if data.getStatus() == "optimal":
            solution.raw.optimality = Optimality.OPTIMALITY_OPTIMAL

        return solution

This completes the Solver Adapter 🎉

Note

You can add parameter arguments in the inherited class in Python, so you can define additional parameters as follows. However, while this allows you to use various features of the backend solver, it may compromise compatibility with other adapters, so carefully consider when creating an adapter.

    @classmethod
    def solve(
        cls,
        ommx_instance: Instance,
        *,
        timeout: Optional[int] = None,
    ) -> Solution:

Solving a Knapsack Problem Using the Solver Adapter#

For verification, let’s solve a knapsack problem using this.

v = [10, 13, 18, 31, 7, 15]
w = [11, 25, 20, 35, 10, 33]
W = 47
N = len(v)

x = [
    DecisionVariable.binary(
        id=i,
        name="x",
        subscripts=[i],
    )
    for i in range(N)
]
instance = Instance.from_components(
    decision_variables=x,
    objective=sum(v[i] * x[i] for i in range(N)),
    constraints=[sum(w[i] * x[i] for i in range(N)) - W <= 0],
    sense=Instance.MAXIMIZE,
)

solution = OMMXPySCIPOptAdapter.solve(instance)

Implementing a Sampler Adapter#

Next, let’s create a Sampler Adapter using OpenJij. OpenJij includes openjij.SASampler for Simulated Annealing (SA) and openjij.SQASampler for Simulated Quantum Annealing (SQA). In this tutorial, we will use SASampler as an example.

For simplicity, this tutorial omits the parameters passed to OpenJij. For more details, refer to the implementation of ommx-openjij-adapter. For how to use the OpenJij Adapter, refer to Sampling from QUBO with OMMX Adapter.

Converting openjij.Response to ommx.v1.Samples#

OpenJij manages decision variables with IDs that are not necessarily sequential, similar to OMMX, so there is no need to create an ID correspondence table as in the case of PySCIPOpt.

The sample results from OpenJij are obtained as openjij.Response, so implement a function to convert this to ommx.v1.Samples. OpenJij returns the number of occurrences of the same sample as num_occurrence. On the other hand, ommx.v1.Samples has unique sample IDs for each sample, and the same value samples are compressed as SamplesEntry. Note that a conversion is needed to bridge this difference.

import openjij as oj
from ommx.v1 import Instance, SampleSet, Solution, Samples, State

def decode_to_samples(response: oj.Response) -> Samples:
    # Generate sample IDs
    sample_id = 0
    entries = []

    num_reads = len(response.record.num_occurrences)
    for i in range(num_reads):
        sample = response.record.sample[i]
        state = State(entries=zip(response.variables, sample))
        # Encode `num_occurrences` into a list of sample IDs
        ids = []
        for _ in range(response.record.num_occurrences[i]):
            ids.append(sample_id)
            sample_id += 1
        entries.append(Samples.SamplesEntry(state=state, ids=ids))
    return Samples(entries=entries)

Note that at this stage, ommx.v1.Instance or its extracted correspondence table is not needed because there is no need to consider ID correspondence.

Implementing a Class that Inherits ommx.adapter.SamplerAdapter#

In the case of PySCIPOpt, we inherited SolverAdapter, but this time we will inherit SamplerAdapter. This has three @abstractmethod as follows:

class SamplerAdapter(SolverAdapter):
    @classmethod
    @abstractmethod
    def sample(cls, ommx_instance: Instance) -> SampleSet:
        pass

    @property
    @abstractmethod
    def sampler_input(self) -> SamplerInput:
        pass

    @abstractmethod
    def decode_to_sampleset(self, data: SamplerOutput) -> SampleSet:
        pass

SamplerAdapter inherits from SolverAdapter, so you might think you need to implement solve and other @abstractmethod. However, since SamplerAdapter has a function to return the best sample using sample, it is sufficient to implement only sample. If you want to implement a more efficient implementation yourself, override solve.

from ommx.adapter import SamplerAdapter

class OMMXOpenJijSAAdapter(SamplerAdapter):
    """
    Sampling QUBO with Simulated Annealing (SA) by `openjij.SASampler`
    """

    # Retain the Instance because it is required to convert to SampleSet
    ommx_instance: Instance
    
    def __init__(self, ommx_instance: Instance):
        self.ommx_instance = ommx_instance

    # Perform sampling
    def _sample(self) -> oj.Response:
        sampler = oj.SASampler()
        # Convert to QUBO dictionary format
        # If the Instance is not in QUBO format, an error will be raised here
        qubo, _offset = self.ommx_instance.as_qubo_format()
        return sampler.sample_qubo(qubo)

    # Common method for performing sampling
    @classmethod
    def sample(cls, ommx_instance: Instance) -> SampleSet:
        adapter = cls(ommx_instance)
        response = adapter._sample()
        return adapter.decode_to_sampleset(response)
    
    # In this adapter, `SamplerInput` uses a QUBO dictionary
    @property
    def sampler_input(self) -> dict[tuple[int, int], float]:
        qubo, _offset = self.ommx_instance.as_qubo_format()
        return qubo
   
    # Convert OpenJij Response to a SampleSet
    def decode_to_sampleset(self, data: oj.Response) -> SampleSet:
        samples = decode_to_samples(data)
        # The information stored in `ommx.v1.Instance` is required here
        return self.ommx_instance.evaluate_samples(samples)

Summary#

In this tutorial, we learned how to implement an OMMX Adapter by connecting to PySCIPOpt as a Solver Adapter and OpenJij as a Sampler Adapter. Here are the key points when implementing an OMMX Adapter:

  1. Implement an OMMX Adapter by inheriting the abstract base class SolverAdapter or SamplerAdapter.

  2. The main steps of the implementation are as follows:

    • Convert ommx.v1.Instance into a format that the backend solver can understand.

    • Run the backend solver to obtain a solution.

    • Convert the backend solver’s output into ommx.v1.Solution or ommx.v1.SampleSet.

  3. Understand the characteristics and limitations of each backend solver and handle them appropriately.

  4. Pay attention to managing IDs and mapping variables to bridge the backend solver and OMMX.

If you want to connect your own backend solver to OMMX, refer to this tutorial for implementation. By implementing an OMMX Adapter following this tutorial, you can use optimization with various backend solvers through a common API.

For more detailed implementation examples, refer to the repositories such as ommx-pyscipopt-adapter and ommx-openjij-adapter.