BENCH: optimize: milp: add MILP benchmarks

benchmarks/benchmarks/optimize_linprog.py

@@ -78,7 +78,7 @@ def setup(self, meth, prob):
                 raise NotImplementedError("skipped")
 
         dims, obj = prob
-        self.A_eq, self.b_eq, self.c, numbers = magic_square(dims)
+        self.A_eq, self.b_eq, self.c, numbers, _ = magic_square(dims)
         self.fun = None
 
     def time_magic_square(self, meth, prob):

benchmarks/benchmarks/optimize_milp.py

@@ -0,0 +1,71 @@
+import os
+
+import numpy as np
+from numpy.testing import assert_allclose
+
+from .common import Benchmark, safe_import
+
+with safe_import():
+    from scipy.optimize import milp
+
+with safe_import():
+    from scipy.optimize.tests.test_linprog import magic_square
+
+
+# MIPLIB 2017 benchmarks included with permission of the authors
+milp_problems = ["piperout-27"]
+
+
+class MilpMiplibBenchmarks(Benchmark):
+    params = [milp_problems]
+    param_names = ['problem']
+
+    def setup(self, prob):
+        if not hasattr(self, 'data'):
+            dir_path = os.path.dirname(os.path.realpath(__file__))
+            datafile = os.path.join(dir_path, "linprog_benchmark_files",
+                                    "milp_benchmarks.npz")
+            self.data = np.load(datafile, allow_pickle=True)
+
+        c, A_ub, b_ub, A_eq, b_eq, bounds, integrality = self.data[prob]
+
+        lb = [l for l, u in bounds]
+        ub = [u for l, u in bounds]
+
+        cons = []
+        if A_ub is not None:
+            cons.append((A_ub, -np.inf, b_ub))
+        if A_eq is not None:
+            cons.append((A_eq, b_eq, b_eq))
+
+        self.c = c
+        self.constraints = cons
+        self.bounds = (lb, ub)
+        self.integrality = integrality
+
+    def time_milp(self, prob):
+        res = milp(c=self.c, constraints=self.constraints, bounds=self.bounds,
+                   integrality=self.integrality)
+        assert res.success
+
+
+class MilpMagicSquare(Benchmark):
+
+    params = [[3, 4, 5, 6]]
+    param_names = ['size']
+
+    def setup(self, n):
+        A_eq, b_eq, self.c, self.numbers, self.M = magic_square(n)
+        self.constraints = (A_eq, b_eq, b_eq)
+
+    def time_magic_square(self, n):
+        res = milp(c=self.c*0, constraints=self.constraints,
+                      bounds=(0, 1), integrality=True)
+        assert res.status==0
+        x = np.round(res.x)
+        s = (self.numbers.flatten() * x).reshape(n**2, n, n)
+        square = np.sum(s, axis=0)
+        assert_allclose(square.sum(axis=0), self.M)
+        assert_allclose(square.sum(axis=1), self.M)
+        assert_allclose(np.diag(square).sum(), self.M)
+        assert_allclose(np.diag(square[:, ::-1]).sum(), self.M)