Keep developing benchmarks

Author: staadecker

Snakefile

@@ -22,12 +22,20 @@ def generate_all_runs():
     return runs
 
 rule all:
+    input:
+        "benchmarks/results.csv"
+
+rule collect_benchmarks:
     input:
         [f"benchmarks/{problem}/results/{library}_{solver}_{size}.tsv"
          for problem, size, library, solver in generate_all_runs()]
+    output:
+        "benchmarks/results.csv"
+    script:
+        "benchmarks/collect_benchmarks.py"
 
 rule run_benchmark:
     benchmark:
         repeat("benchmarks/{problem}/results/{library}_{solver}_{size}.tsv", config["repeat"])
     shell:
-        "python benchmarks/run_benchmarks {wildcards.problem} --library {wildcards.library} --solver {wildcards.solver} --size {wildcards.size}"
+        "python benchmarks/run_benchmark.py {wildcards.problem} --library {wildcards.library} --solver {wildcards.solver} --size {wildcards.size}"

benchmarks/.gitignore

@@ -0,0 +1,3 @@
+*/data/**
+*/results/*.tsv
+*.csv

benchmarks/collect_benchmarks.py

@@ -0,0 +1,34 @@
+from pathlib import Path
+
+import polars as pl
+from energy_model.scripts.util import mock_snakemake
+
+
+def collect_benchmarks(input_files):
+    dfs = []
+    for input_file in input_files:
+        input_file = Path(input_file)
+        problem = input_file.parent.parent.name
+        name = input_file.name.partition(".")[0].split("_")
+        library = name[0]
+        solver = name[1]
+        size = name[2]
+        dfs.append(
+            pl.read_csv(input_file, separator="\t")
+            .mean()
+            .with_columns(
+                problem=pl.lit(problem),
+                library=pl.lit(library),
+                solver=pl.lit(solver),
+                size=pl.lit(size),
+            )
+        )
+    return pl.concat(dfs)
+
+
+if __name__ == "__main__":
+    if "snakemake" not in globals() or hasattr(snakemake, "mock"):  # noqa: F821
+        snakemake = mock_snakemake("collect_benchmarks")
+
+    result = collect_benchmarks(snakemake.input)
+    result.write_csv(snakemake.output[0])

benchmarks/config.yaml

@@ -1,4 +1,4 @@
-repeat: 3
+repeat: 1
 problems:
   facility_problem:
     size:
@@ -20,4 +20,4 @@ libraries:
   # - gurobipy
   # - jump
   # - linopy
-  # - pyoptinterface
+  - pyoptinterface

benchmarks/energy_model/.gitignore

@@ -0,0 +1,44 @@
+# Facility Location Problem
+
+This is the exact same problem as described in the original [JuMP paper](https://epubs.siam.org/doi/10.1137/15M1020575) (see Section 8.2). It is re-explained here for clarity.
+
+## Problem description
+
+The purpose of this problem is to choose the location of $F$ facilities such that the maximum distance between any customer and its nearest facility is minimized. The $C^2$ Customers are evenly spread out over a square grid (see image below). Facilities can be placed anywhere in the grid.
+
+TODO add image
+
+
+## Problem formulation
+
+To simplify the formulation, let us scale our square grid such that its axes span go from $0$ to $1$.
+
+We define a variable $y_{f,ax}$ which is the location of every facility ($f$) for each axis ($ax$) (either x-axis or y-axis).
+
+$$ 0 \leq y_{f,ax} \leq 1$$
+
+Next, we define the distance between every facility-customer pair, $s_{i,j,f}$.
+ $$r^x_{i,j,f} = (i / N - y_{f,1}) \qquad \forall i, j, f$$
+$$r^y_{i,j,f} = (j / N - y_{f,2}) \qquad \forall i,j,f$$
+
+$$s_{i,j,f} ^2 = (r_{i,j,f}^x) ^ 2 + (r^y_{i,j,f}) ^2 \qquad \forall i,j,f$$
+$$ s_{i,j,f} \geq 0 \qquad \forall i,j,f$$
+
+We'd like to minimize the maximum distance between any facility and its nearest neighbor $d$.
+
+$$\text{min} \quad d$$
+
+At first, we might try to simply define $d$ as follows:
+
+$$ d^{max} \geq s_{i,j,f} \qquad \forall i,j,f$$
+
+However this is not quite right. We only wish to ensure the maximum distance is bounded by the distances between every customer and its _nearest_ facility. To do this, suppose we had a binary variable, $z_{i,j,f}$ that equals $1$ for customer-facility pairs that are relevant (when the facility is the nearest to the customer) and $0$ otherwise. We can now rewrite the above constraint as follows to fix our issue.
+
+$$ d^{max} \geq s_{i,j,f} - \sqrt{2} (1 - z_{i,j,f}) \qquad \forall i,j,f $$
+
+Why does this work? When $z_{i,j,f} = 0$ the right hand side is necessarily zero or less (since $\sqrt{2}$ is the largest possible distance in a 1-by-1 square) and the constraint is thus non-binding as desired.
+
+Finally, how do we define this magical $z_{i,j,f}$ variable? We simply ensure that exactly one customer-facility pair exists for each customer. The optimization will automatically choose the nearest facility if it matters since it wishes to minimize the objective.
+
+$$\sum_{f} z_{i,j,f} = 1 \qquad \forall i,j$$
+

benchmarks/facility_problem/README.md

@@ -34,7 +34,6 @@ def solve(solver, size=None, G=None, F=None, is_benchmarking=True, use_var_names
         model.x_axis, model.y_axis, model.facilities, vtype="binary"
     )
     model.con_only_one_closest = model.is_closest.sum("f") == 1
-
     model.dist_x = Variable(model.x_axis, model.facilities)
     model.dist_y = Variable(model.y_axis, model.facilities)
     model.con_dist_x = model.dist_x == customer_position_x.to_expr().add_dim(

benchmarks/facility_problem/bm_pyoframe.py

@@ -0,0 +1,78 @@
+# Copyright (c) 2022: Miles Lubin and contributors
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+# See https://github.com/jump-dev/JuMPPaperBenchmarks
+
+import pyomo.environ as pyo
+from pyomo.opt import SolverFactory
+
+
+def solve(solver, size, is_benchmarking=True):
+    if type(size) == int:
+        size = (size, size)
+    G, F = size
+    model = pyo.ConcreteModel()
+    model.G = G
+    model.F = F
+    model.Grid = pyo.RangeSet(0, model.G)
+    model.Facs = pyo.RangeSet(1, model.F)
+    model.Dims = pyo.RangeSet(1, 2)
+    model.y = pyo.Var(model.Facs, model.Dims, bounds=(0.0, 1.0))
+    model.s = pyo.Var(model.Grid, model.Grid, model.Facs, bounds=(0.0, None))
+    model.z = pyo.Var(model.Grid, model.Grid, model.Facs, within=pyo.Binary)
+    model.r = pyo.Var(model.Grid, model.Grid, model.Facs, model.Dims)
+    model.d = pyo.Var()
+    model.obj = pyo.Objective(expr=1.0 * model.d)
+
+    def assmt_rule(mod, i, j):
+        return sum([mod.z[i, j, f] for f in mod.Facs]) == 1
+
+    model.assmt = pyo.Constraint(model.Grid, model.Grid, rule=assmt_rule)
+    M = 2 * 1.414
+
+    def quadrhs_rule(mod, i, j, f):
+        return mod.s[i, j, f] == mod.d + M * (1 - mod.z[i, j, f])
+
+    model.quadrhs = pyo.Constraint(
+        model.Grid, model.Grid, model.Facs, rule=quadrhs_rule
+    )
+
+    def quaddistk1_rule(mod, i, j, f):
+        return mod.r[i, j, f, 1] == (1.0 * i) / mod.G - mod.y[f, 1]
+
+    model.quaddistk1 = pyo.Constraint(
+        model.Grid, model.Grid, model.Facs, rule=quaddistk1_rule
+    )
+
+    def quaddistk2_rule(mod, i, j, f):
+        return mod.r[i, j, f, 2] == (1.0 * j) / mod.G - mod.y[f, 2]
+
+    model.quaddistk2 = pyo.Constraint(
+        model.Grid, model.Grid, model.Facs, rule=quaddistk2_rule
+    )
+
+    def quaddist_rule(mod, i, j, f):
+        return mod.r[i, j, f, 1] ** 2 + mod.r[i, j, f, 2] ** 2 <= mod.s[i, j, f] ** 2
+
+    model.quaddist = pyo.Constraint(
+        model.Grid, model.Grid, model.Facs, rule=quaddist_rule
+    )
+    opt = SolverFactory(solver)
+    if is_benchmarking:
+        opt.options["timelimit"] = 0.0
+        opt.options["presolve"] = False
+    try:
+        if solver == "gurobi_persistent":
+            opt.set_instance(model)
+            opt.solve(tee=True)
+        else:
+            opt.solve(model, tee=True)
+    except ValueError as e:
+        if "bad status: aborted" not in str(e):
+            raise e
+    return model
+
+
+if __name__ == "__main__":
+    solve("gurobi", 5)

benchmarks/facility_problem/bm_pyomo.py

@@ -0,0 +1,75 @@
+# Copyright (c) 2023: Yue Yang
+# https://github.com/metab0t/PyOptInterface_benchmark
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+
+import pyoptinterface as poi
+from pyoptinterface import gurobi, copt
+import numpy as np
+
+
+def solve(solver, size):
+    model_constructor = {
+        "gurobi": gurobi.Model,
+        "copt": copt.Model,
+    }.get(solver, None)
+    if model_constructor is None:
+        raise ValueError(f"Unknown solver {solver}")
+
+    m = model_constructor()
+    if type(size) == int:
+        size = (size, size)
+    G, F = size
+    # Create variables
+    y = add_ndarray_variable(m, (F, 2), lb=0.0, ub=1.0)
+    s = add_ndarray_variable(m, (G + 1, G + 1, F), lb=0.0)
+    z = add_ndarray_variable(m, (G + 1, G + 1, F), domain=poi.VariableDomain.Binary)
+    r = add_ndarray_variable(m, (G + 1, G + 1, F, 2))
+    d = m.add_variable()
+
+    # Set objective
+    m.set_objective(d * 1.0)
+
+    # Add constraints
+    for i in range(G + 1):
+        for j in range(G + 1):
+            expr = poi.quicksum(z[i, j, :])
+            m.add_linear_constraint(expr, poi.Eq, 1.0)
+
+    M = 2 * 1.414
+    for i in range(G + 1):
+        for j in range(G + 1):
+            for f in range(F):
+                expr = s[i, j, f] - d - M * (1 - z[i, j, f])
+                m.add_linear_constraint(expr, poi.Eq, 0.0)
+                expr = r[i, j, f, 0] - i / G + y[f, 0]
+                m.add_linear_constraint(expr, poi.Eq, 0.0)
+                expr = r[i, j, f, 1] - j / G + y[f, 1]
+                m.add_linear_constraint(expr, poi.Eq, 0.0)
+                m.add_second_order_cone_constraint(
+                    [s[i, j, f], r[i, j, f, 0], r[i, j, f, 1]]
+                )
+
+    # Optimize model
+    m.set_model_attribute(poi.ModelAttribute.Silent, True)
+    m.set_model_attribute(poi.ModelAttribute.TimeLimitSec, 0.0)
+    # close presolve for gurobi and open for copt
+    solver_name = m.get_model_attribute(poi.ModelAttribute.SolverName)
+    if solver_name.lower() == "gurobi":
+        m.set_raw_parameter("Presolve", 0)
+    elif solver_name.lower() == "copt":
+        m.set_raw_parameter("Presolve", 1)
+    m.optimize()
+
+
+def add_ndarray_variable(m, shape, **kwargs):
+    array = np.empty(shape, dtype=object)
+    array_flat = array.flat
+    for i in range(array.size):
+        array_flat[i] = m.add_variable(**kwargs)
+    return array
+
+
+if __name__ == "__main__":
+    solve("gurobi", 5)

benchmarks/facility_problem/bm_pyoptinterface.py

@@ -20,12 +20,14 @@ def get_problems() -> Collection[str]:
         prog="run_benchmark", description="Run a benchmark with a given solver."
     )
     parser.add_argument("problem", choices=get_problems())
-    parser.add_argument("--solver", choices=SUPPORTED_SOLVERS)
+    parser.add_argument("--solver", choices=[s.name for s in SUPPORTED_SOLVERS])
     parser.add_argument("--size", type=int)
     parser.add_argument("--library")
     args = parser.parse_args()
 
     # import solve function dynamically
-    solve_func = importlib.import_module(f"benchmarks.bm_{args.library.lower()}").solve
+    solve_func = importlib.import_module(
+        f"{args.problem}.bm_{args.library.lower()}"
+    ).solve
 
     solve_func(solver=args.solver, size=args.size)