OMMX Adapterを実装する#
複数のAdapterで最適化問題を解いて結果を比較するで触れた通り、OMMX Adapterは共通化されたAPIを持っています。この共通化されたAPIは、OMMX Python SDKが用意している抽象基底クラスを継承することで実現されています。OMMXはAdapterの性質に応じて二つの抽象基底クラスを用意しています。
ommx.adapter.SolverAdapter: 一つの解を返す最適化ソルバーのための抽象基底クラスommx.adapter.SamplerAdapter: サンプリングベースの最適化ソルバーのための抽象基底クラス
複数の解が得られるソルバーは、特に複数得られたサンプルのうち最善のものを選択することによって、自動的に単一の解を返すソルバーと見なす事ができるため、SamplerAdapter は SolverAdapter を継承しています。Adapterを作るときにどちらを実装するか悩んだら、出力される解の数を見て、一つの解を返すなら SolverAdapter、複数の解を返すなら SamplerAdapter を継承すると良いでしょう。たとえば PySCIPOpt などの厳密解を1つ返す最適化ソルバーはSolverAdapter を使い、OpenJij などの複数の解を返すサンプラーはSamplerAdapter を使います。
OMMXでは ommx.adapter.SolverAdapter を継承したクラスを Solver Adapter、ommx.adapter.SamplerAdapter を継承したクラスを Sampler Adapter と呼びます。
またここでの説明のため、PySCIPOptやOpenJijのようにAdapterがラップしようとしているソフトウェアのことをバックエンドソルバーと呼びます。
Adapterの処理の流れ#
Adapterの処理は大雑把にいうと次の3ステップからなります：
ommx.v1.Instanceをバックエンドソルバーが読める形式に変換するバックエンドソルバーを実行して解を取得する
バックエンドソルバーの解を
ommx.v1.Solutionやommx.v1.SampleSetに変換して返す
2.はバックエンドソルバーの使い方そのものなので、これは既知としてこのチュートリアルでは扱いません。ここでは 1. と 3. をどのように実装するかを説明します。
多くのバックエンドソルバーが最適化問題を表すための必要な最小限の情報だけを、アルゴリズムに応じた形で受け取るように作られているのに比べて、ommx.v1.Instance はデータ分析の一部として最適化を行うことを想定しているためより多くの情報を持っています。なのでステップ 1. では多くの情報を削ぎ落とすことになります。またOMMXでは決定変数や制約条件は連番とは限らないIDで管理されていますが、バックエンドソルバーによっては名前で管理されいたり、連番で管理されていることがあります。この対応関係は 3. の処理で必要になるのでAdapterが管理しておく必要があります。
逆にステップ 3. では ommx.v1.Solution や ommx.v1.SampleSet はバックエンドソルバーの出力だけからは構築できないので、まずバックエンドソルバーの返した解と 1. の時の情報から ommx.v1.State あるいは ommx.v1.Samples を構築し、それを ommx.v1.Instance を使って ommx.v1.Solution や ommx.v1.SampleSet に変換します。
Solver Adapterを実装する#
ここでは PySCIPOpt を例としてSolver Adapterを実装してみましょう。なお完全な例は ommx-pyscipopt-adapter を確認してください。
ここではチュートリアルということで、順番に実行しやすいように以下のように作業します。
ommx.v1.Instanceから PySCIPOpt のモデルを構築するための関数を順番に実装していきます。最後にこれらの関数を
OMMXPySCIPOptAdapterクラスとしてまとめます
カスタム例外#
まずカスタム例外を定義しておくといいでしょう。これによりユーザーは例外が発生したときに、どの部分が問題を引き起こしているのかを理解しやすくなります。
class OMMXPySCIPOptAdapterError(Exception):
pass
OMMXは広いクラスの最適化問題を保存できるようになっているので、バックエンドソルバーが対応していない問題が入力されるケースがあります。その場合はエラーを投げるようにしてください。
決定変数を設定する#
PySCIPOptは決定変数を名前で管理するので、OMMXの決定変数のIDを文字列にして名前として登録します。これにより後述する decode_to_state においてPySCIPOptの決定変数から ommx.v1.State を復元することができます。これはバックエンドソルバーの実装に応じて適切な方法が変わることに注意してください。重要なのは解を得た後に ommx.v1.State に変換するための情報を保持することです。
import pyscipopt
from ommx.v1 import Instance, Solution, DecisionVariable, Constraint, State, Optimality, Function
def set_decision_variables(
model: pyscipopt.Model, # チュートリアルのために状態を引数で受け取っているがclassで管理するのが一般的
instance: Instance
) -> dict[str, pyscipopt.Variable]:
"""
モデルに決定変数を追加し、変数名のマッピングを作成して返す
"""
# OMMXの決定変数の情報からPySCIPOptの変数を作成
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:
# 未対応の決定変数の種類がある場合はエラー
raise OMMXPySCIPOptAdapterError(
f"Unsupported decision variable kind: "
f"id: {var.id}, kind: {var.kind}"
)
# 目的関数が2次の場合、線形化のために補助変数を追加
if instance.raw.objective.HasField("quadratic"):
model.addVar(
name="auxiliary_for_linearized_objective", vtype="C", lb=None, ub=None
)
# モデルに追加された変数へアクセスするための辞書を作成
return {var.name: var for var in model.getVars()}
ommx.v1.Function を pyscipopt.Expr に変換する#
ommx.v1.Function を pyscipopt.Expr に変換するための関数を実装します。ommx.v1.Function はOMMXの決定変数のIDしか持っていないので、IDからPySCIPOpt側の変数を取得する必要があり、そのために set_decision_variables で作成した変数名と変数のマッピングを使います。
def make_linear_expr(function: Function, varname_map: dict) -> pyscipopt.Expr:
"""線形式を生成するヘルパー関数"""
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:
"""2次式を生成するヘルパー関数"""
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
目的関数と制約条件を設定する#
pyscipopt.Model に目的関数と制約条件を追加します。この部分はバックエンドソルバーが何をどのようにサポートしているのかの知識が必要になります。例えば以下のコードでは、PySCIPOptが目的関数として2次式を直接扱うことができないため、PySCIPOptのドキュメントに従って補助変数を導入しています。
import math
def set_objective(model: pyscipopt.Model, instance: Instance, varname_map: dict):
"""モデルに目的関数を設定"""
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"):
# PySCIPOptでは2次の目的関数を直接サポートしていないため、補助変数を使って線形化
auxilary_var = varname_map["auxiliary_for_linearized_objective"]
# 補助変数を目的関数として設定
model.setObjective(auxilary_var, sense=sense)
# 補助変数に対する制約を追加
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`, `quadratic`."
)
def set_constraints(model: pyscipopt.Model, instance: Instance, varname_map: dict):
"""モデルに制約条件を設定"""
# 通常の制約条件を処理
for constraint in instance.raw.constraints:
# 制約関数の種類に基づいて式を生成
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"):
# 定数制約の場合、実行可能かどうかをチェック
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 was found: id {constraint.id}"
)
else:
raise OMMXPySCIPOptAdapterError(
f"Constraints must be either `constant`, `linear` or `quadratic`."
f"id: {constraint.id}, "
f"type: {constraint.function.WhichOneof('function')}"
)
# 制約種別（等式/不等式）に基づいて制約を追加
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: "
f"id: {constraint.id}, equality: {constraint.equality}"
)
# 制約をモデルに追加
model.addCons(constr_expr, name=str(constraint.id))
また、バックエンドソルバーが特殊な制約条件（例: SOS制約 など）をサポートしている場合は、それに対応するための関数を追加する必要があります。
以上で ommx.v1.Instance から pycscipopt.Model が構築できるようになりました。
得られた解を ommx.v1.State に変換する#
次に、PySCIPOptのモデルを解いて得られた解を ommx.v1.State に変換する関数を実装します。まず解けているかを確認します。SCIPには最適性を保証する機能や解が非有界であることを検知する機能があるので、それらを検知していたら対応した例外を投げます。これもバックエンドソルバーに依存します。
Warning
特に ommx.adapter.InfeasibleDetected は解がInfeasibleではなくて最適化問題自体がInfeasible、つまり 一つも解を持ち得ないことが保証できた という意味であることに注意してください。ヒューリスティックソルバーが一つも実行可能解を見つけられなかった場合にこれを使ってはいけません。
from ommx.adapter import InfeasibleDetected, UnboundedDetected
def decode_to_state(model: pyscipopt.Model, instance: Instance) -> State:
"""最適化済みのPySCIPOpt Modelからommx.v1.Stateを作成する"""
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:
# 最適解を取得
sol = model.getBestSol()
# 変数名と変数のマッピングを作成
varname_map = {var.name: var for var in model.getVars()}
# 変数IDと値のマッピングを持つStateを作成
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()}]"
)
ommx.adapter.SolverAdapter を継承した class を作る#
最後に、Adapter毎のAPIを揃えるために ommx.adapter.SolverAdapter を継承したクラスを作成します。これは @abstractmethod を含む次のような抽象基底クラスです：
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
この抽象基底クラスは以下の2通りのユースケースを想定しています:
バックエンドソルバーのパラメータなどを調整しない場合は、
solveクラスメソッドを使う。バックエンドソルバーのパラメータなどを調整する場合は、
solver_inputを使ってバックエンドソルバーの入力用のデータ構造（今回はpyscipopt.Model）を取得し、調整した後にバックエンドソルバーへ入力し、最後にバックエンドソルバーの出力をdecodeで変換する。
ここまでで用意した関数を使って次のように実装することができます：
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()
# 関数を使用してモデルを構築
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:
"""
PySCIPoptを使ってommx.v1.Instanceを解き、ommx.v1.Solutionを返す
"""
adapter = cls(ommx_instance)
model = adapter.solver_input
model.optimize()
return adapter.decode(model)
@property
def solver_input(self) -> pyscipopt.Model:
"""生成されたPySCIPOptモデルを返す"""
return self.model
def decode(self, data: pyscipopt.Model) -> Solution:
"""
最適化後のpyscipopt.ModelとOMMX Instanceからommx.v1.Solutionを生成する
"""
# 解の状態をチェック
if data.getStatus() == "infeasible":
raise InfeasibleDetected("Model was infeasible")
if data.getStatus() == "unbounded":
raise UnboundedDetected("Model was unbounded")
# 解を状態に変換
state = decode_to_state(data, self.instance)
# インスタンスを使用して解を評価
solution = self.instance.evaluate(state)
# 最適性ステータスを設定
if data.getStatus() == "optimal":
solution.raw.optimality = Optimality.OPTIMALITY_OPTIMAL
return solution
これでSolver Adapter完成です 🎉
Note
Pythonは継承したクラスでパラメータ引数を追加してもいいので、次のように追加のパラメータを定義することもできます。ただし、これによってバックエンドソルバーの様々な機能が使えるようになる一方、他のAdapterとの互換性が損なわれるので、Adapterを作る際には慎重に検討してください。
@classmethod
def solve(
cls,
ommx_instance: Instance,
*,
timeout: Optional[int] = None,
) -> Solution:
Solver Adapterを使ってナップザック問題を解く#
動作確認のため、これを使ってナップザック問題を解いてみましょう
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)
Sampler Adapterを実装する#
次にOpenJijを使ったSampler Adapterを作ってみましょう。OpenJijには Simulated Annealing (SA) による openjij.SASamplerと Simulated Quantum Annealing (SQA) による openjij.SQASampler が含まれています。このチュートリアルでは、 SASampler を例に説明します。
このチュートリアルでは簡単のためにOpenJijに渡すパラメータは省略しています。詳しくは ommx-openjij-adapter の実装を参照してください。また OpenJij Adapterの使い方については OMMX AdapterでQUBOからサンプリングする を参照してください。
openjij.Response から ommx.v1.Samples への変換#
OpenJijは決定変数をOMMXと同様に連番とは限らないIDで管理しているので、PySCIPOptの時のようにIDの対応表を作る必要はありません。
OpenJijのサンプル結果は openjij.Response として得られるので、これを ommx.v1.Samples に変換する関数を実装します。OpenJijは同じサンプルが得られた時、それが発生した回数を num_occurrence として返します。一方 ommx.v1.Samples は個々のサンプルが固有のサンプルIDをもち、同じ値を持つサンプルは SamplesEntry として圧縮されます。この差異を埋めるための変換が必要なことに注意します。
import openjij as oj
from ommx.v1 import Instance, SampleSet, Solution, Samples, State
def decode_to_samples(response: oj.Response) -> Samples:
# サンプルIDを生成
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))
# `num_occurrences`をサンプルIDリストにエンコード
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)
IDの対応を考えなくて良いため、この段階では ommx.v1.Instance やその情報を抽出した対応表などが必要ないことに注意してください。
ommx.adapter.SamplerAdapter を継承したクラスの実装#
PySCIPOptの時は SolverAdapter を継承しましたが、今回は SamplerAdapter を継承します。これは次のように3つの @abstractmethod を持っています。
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 は SolverAdapter を継承しているので solve などの @abstractmethod も実装する必要と思うかもしれません。しかし、これらについては sample を使って最善のサンプルを返すという機能が SamplerAdapter に実装されているため、sample だけを実装すれば十分です。自分でより効率の良い実装を行いたい場合は solve をオーバーライドしてください。
from ommx.adapter import SamplerAdapter
class OMMXOpenJijSAAdapter(SamplerAdapter):
"""
Sampling QUBO with Simulated Annealing (SA) by `openjij.SASampler`
"""
# SampleSetに変換する必要があるので、Instanceを保持
ommx_instance: Instance
def __init__(self, ommx_instance: Instance):
self.ommx_instance = ommx_instance
# サンプリングを行う
def _sample(self) -> oj.Response:
sampler = oj.SASampler()
# QUBOの辞書形式に変換
# InstanceがQUBO形式でなければここでエラーになる
qubo, _offset = self.ommx_instance.as_qubo_format()
return sampler.sample_qubo(qubo)
# サンプリングを行う共通のメソッド
@classmethod
def sample(cls, ommx_instance: Instance) -> SampleSet:
adapter = cls(ommx_instance)
response = adapter._sample()
return adapter.decode_to_sampleset(response)
# このAdapterでは `SamplerInput` は QUBO形式の辞書を使うことにする
@property
def sampler_input(self) -> dict[tuple[int, int], float]:
qubo, _offset = self.ommx_instance.as_qubo_format()
return qubo
# OpenJijのResponseをSampleSetに変換
def decode_to_sampleset(self, data: oj.Response) -> SampleSet:
samples = decode_to_samples(data)
# ここで `ommx.v1.Instance` が保持している情報が必要になる
return self.ommx_instance.evaluate_samples(samples)
まとめ#
このチュートリアルでは、PySCIPOptと接続するSolver Adapterの実装とOpenJijと接続するSampler Adapterの実装を通して、OMMX Adapterの実装方法について学びました。以下がOMMX Adapterを実装する際の重要なポイントです：
OMMX Adapterは
SolverAdapterまたはSamplerAdapterの抽象基底クラスを継承することで実装します実装の主なステップは以下の通りです：
ommx.v1.Instanceをバックエンドソルバーが理解できる形式に変換するバックエンドソルバーを実行して解を取得する
バックエンドソルバーの出力を
ommx.v1.Solutionやommx.v1.SampleSetに変換する
各バックエンドソルバーの特性や制限を理解し、適切に処理する必要があります
IDの管理や変数の対応付けなど、バックエンドソルバーとOMMXの橋渡しに注意を払う必要があります
独自のバックエンドソルバーをOMMXと接続したい場合は、このチュートリアルを参考に実装すると良いでしょう。このチュートリアルに従ってOMMX Adapterを実装することで、様々なバックエンドソルバーでの最適化を共通化されたAPIで利用できるようになります。
より詳しい実装例については、ommx-pyscipopt-adapterやommx-openjij-adapterなどのリポジトリを参照してください。