First version of symmip, which adds tentative MIP support

.github/workflows/test_symmip.yaml

@@ -0,0 +1,34 @@
+name: pytest_and_symmip
+env:
+  version: 8.0.3
+on: [push, pull_request]
+
+jobs:
+  build:
+
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest]
+        python-version: [3.7, 3.8, 3.9, "3.10", 3.11]
+    steps:
+    - uses: actions/checkout@v3
+    - name: Install dependencies (SCIPOptSuite)
+      run: |
+        wget --quiet --no-check-certificate https://github.com/scipopt/scip/releases/download/$(echo "v${{env.version}}" | tr -d '.')/SCIPOptSuite-${{ env.version }}-Linux-ubuntu.deb
+        sudo apt-get update && sudo apt install -y ./SCIPOptSuite-${{ env.version }}-Linux-ubuntu.deb
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v4
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install wheel
+        pip install pytest
+        pip install -r requirements.txt
+        pip install matplotlib
+        pip install -e .[symmip]
+    - name: Test MIP
+      run:
+        pytest tests/test_mip.py

README.rst

@@ -2,7 +2,15 @@
    :target: https://zenodo.org/badge/latestdoi/24005390
 .. image:: https://coveralls.io/repos/github/tBuLi/symfit/badge.svg?branch=master
    :target: https://coveralls.io/github/tBuLi/symfit?branch=master
-
+.. image:: https://img.shields.io/pypi/v/symfit?label=pypi%20package
+   :alt: PyPI
+   :target:`https://pypi.org/project/symfit/`
+.. image:: https://img.shields.io/pypi/dm/symfit
+   :alt: PyPI - Downloads
+   :target:`https://pypi.org/project/symfit/`
+.. image:: https://img.shields.io/conda/dn/conda-forge/symfit?color=brightgreen&label=downloads&logo=conda-forge
+   :alt: Conda
+   :target: https://anaconda.org/conda-forge/symfit
 
 Please cite this DOI if ``symfit`` benefited your publication. Building this has been a lot of work, and as young researchers your citation means a lot to us.
 Martin Roelfs & Peter C Kroon, symfit. doi:10.5281/zenodo.1133336

docs/installation.rst

@@ -11,6 +11,15 @@ from your terminal. If you prefer to use `conda`, run ::
 instead. Lastly, if you prefer to install manually you can download
 the source from https://github.com/tBuLi/symfit.
 
+symmip module
+--------------
+To use `symfit`'s :class:`~symfit.symmip.mip.MIP` object for mixed integer programming (MIP) and
+mixed integer nonlinear programming (MINLP), you need to have a suitable backend installed.
+Because this is an optional feature, no such solver is installed by default.
+In order to install the non-commercial SCIPOpt package, install `symfit` by running
+
+  pip install symfit[symmip]
+
 Contrib module
 --------------
 To also install the dependencies of 3rd party contrib modules such as

examples/mip/bilinear.py

@@ -0,0 +1,26 @@
+# Inspired by https://www.gurobi.com/documentation/9.5/examples/bilinear_py.html#subsubsection:bilinear.py
+#
+# This example formulates and solves the following simple bilinear model:
+#  maximize    x
+#  subject to  x + y + z <= 10
+#              x * y <= 2         (bilinear inequality)
+#              x * z + y * z = 1  (bilinear equality)
+#              x, y, z non-negative (x integral in second version)
+
+from symfit import parameters, Eq
+from symfit import MIP
+
+# Create variables
+x, y, z = parameters('x, y, z', min=0)
+
+objective = 1.0 * x
+constraints = [
+    x + y + z <= 10,
+    x*y <= 2,
+    Eq(x*z + y*z, 1),
+]
+
+mip = MIP(- objective, constraints=constraints)
+mip_result = mip.execute()
+
+print(mip_result)

examples/mip/mip1.py

@@ -0,0 +1,35 @@
+# Inspired by https://www.gurobi.com/documentation/9.5/examples/mip1_py.html
+#
+# Solve the following MIP:
+#  maximize
+#        x +   y + 2 z
+#  subject to
+#        x + 2 y + 3 z <= 4
+#        x +   y       >= 1
+#        x, y, z binary
+
+from symfit import parameters, MIP
+from symfit.symmip.backend import SCIPOptBackend, GurobiBackend
+
+x, y, z = parameters('x, y, z', binary=True, min=0, max=1)
+
+objective = x + y + 2 * z
+constraints = [
+    x + 2 * y + 3 * z <= 4,
+    x + y >= 1
+]
+
+# We know solve this problem with different backends.
+for backend in [SCIPOptBackend, GurobiBackend]:
+    print(f'Run with {backend=}:')
+    fit = MIP(objective, constraints=constraints, backend=backend)
+    fit_result = fit.execute()
+
+    print(f"Optimal objective value: {fit_result.objective_value}")
+    print(
+        f"Solution values: "
+        f"x={fit_result[x]}, "
+        f"y={fit_result[y]}, "
+        f"z={fit_result[z]}"
+    )
+    print(fit_result, end='\n\n')

examples/mip/multiscenario.py

@@ -0,0 +1,66 @@
+"""
+Inspired by https://www.gurobi.com/documentation/9.5/examples/multiscenario_py.html.
+
+For now this symfit equivalent only solves a single scenario, as we do not (currently) support Gurobi's
+multiscenario feature. However, this is still a nice example to demonstrate some of symfit's features.
+"""
+
+from symfit import MIP, IndexedBase, Eq, Idx, Parameter, symbols, Sum, pprint
+from symfit.symmip.backend import SCIPOptBackend, GurobiBackend
+import numpy as np
+
+# Warehouse demand in thousands of units
+data_demand = np.array([15, 18, 14, 20])
+
+# Plant capacity in thousands of units
+data_capacity = np.array([20, 22, 17, 19, 18])
+
+# Fixed costs for each plant
+data_fixed_costs = np.array([12000, 15000, 17000, 13000, 16000])
+
+# Transportation costs per thousand units
+data_trans_costs = np.array(
+    [[4000, 2000, 3000, 2500, 4500],
+     [2500, 2600, 3400, 3000, 4000],
+     [1200, 1800, 2600, 4100, 3000],
+     [2200, 2600, 3100, 3700, 3200]]
+)
+
+# Indices over the plants and warehouses
+plant = Idx('plant', range=len(data_capacity))
+warehouse = Idx('warehouse', range=len(data_demand))
+
+# Indexed variables. Initial values become coefficients in the objective function.
+open = IndexedBase(Parameter('Open', binary=True))
+transport = IndexedBase(Parameter('Transport'))
+fixed_costs = IndexedBase(Parameter('fixed_costs'))
+trans_cost = IndexedBase(Parameter('trans_cost'))
+capacity = IndexedBase(Parameter('capacity'))
+demand = IndexedBase(Parameter('demand'))
+
+objective = Sum(fixed_costs[plant] * open[plant], plant) + Sum(trans_cost[warehouse, plant] * transport[warehouse, plant], warehouse, plant)
+constraints = [
+    Sum(transport[warehouse, plant], warehouse) <= capacity[plant] * open[plant],
+    Eq(Sum(transport[warehouse, plant], plant), demand[warehouse])
+]
+
+print('Objective:')
+pprint(objective, wrap_line=False)
+print('\nSubject to:')
+for constraint in constraints:
+    pprint(constraint, wrap_line=False)
+print('\n\n')
+
+data = {
+    fixed_costs[plant]: data_fixed_costs,
+    trans_cost[warehouse, plant]: data_trans_costs,
+    capacity[plant]: data_capacity,
+    demand[warehouse]: data_demand,
+}
+
+# We know solve this problem with different backends.
+for backend in [SCIPOptBackend, GurobiBackend]:
+    print(f'Run with {backend=}:')
+    mip = MIP(objective, constraints=constraints, data=data, backend=backend)
+    results = mip.execute()
+    print(results)

examples/mip/sudoku.py

@@ -0,0 +1,68 @@
+# Inspired by https://www.gurobi.com/documentation/9.5/examples/sudoku_py.html
+
+import math
+
+import numpy as np
+
+from symfit import Parameter, symbols, IndexedBase, Idx, Sum, Eq
+from symfit.symmip.backend import SCIPOptBackend, GurobiBackend
+from symfit import MIP
+
+with open('sudoku1') as f:
+    grid = f.read().split()
+
+for line in grid:
+    print(line)
+
+n = len(grid[0])
+s = math.isqrt(n)
+
+# Fix variables associated with cells whose values are pre-specified
+# lb = np.array([[0 if char == "." else 1 for j, char in enumerate(line)] for i, line in enumerate(grid)])
+lb = np.zeros((n, n, n), dtype=int)
+for i in range(n):
+    for j in range(n):
+        if grid[i][j] != '.':
+            v = int(grid[i][j]) - 1
+            lb[i, j, v] = 1
+ub = np.ones_like(lb)
+
+# Prepare the boolean parameters for our sudoku board.
+# Because every position on the board can have only one value,
+# we make a binary Indexed symbol x[i,j,v], where i is the column,
+# j is the row, and v is the value in the (i, j) position.
+x = IndexedBase(Parameter('x', binary=True, min=lb, max=ub))
+i, j, v = symbols('i, j, v', cls=Idx, range=n)
+x_ijv = x[i, j, v]
+
+# Add the sudoku constraints:
+#   1. Each cell must take exactly one value: Sum(x[i,j,v], v) == 1
+#   2. Each value is used exactly once per row: Sum(x[i,j,v], i) == 1
+#   3. Each value is used exactly once per column: Sum(x[i,j,v], j) == 1
+#   4. Each value is used exactly once per 3x3 subgrid.
+constraints = [
+    Eq(Sum(x[i, j, v], v), 1),
+    Eq(Sum(x[i, j, v], i), 1),
+    Eq(Sum(x[i, j, v], j), 1),
+    *[Eq(Sum(x[i, j, v], (i, i_lb, i_lb + s - 1), (j, j_lb, j_lb + s - 1)), 1)
+      for i_lb in range(0, n, s)
+      for j_lb in range(0, n, s)]
+]
+
+# We know solve this problem with different backends.
+for backend in [SCIPOptBackend, GurobiBackend]:
+    print(f'Run with {backend=}:')
+    fit = MIP(constraints=constraints, backend=backend)
+    result = fit.execute()
+
+    print('')
+    print('Solution:')
+    print('')
+    solution = result[x]
+    for i in range(n):
+        sol = ''
+        for j in range(n):
+            for v in range(n):
+                if solution[i, j, v] > 0.5:
+                    sol += str(v+1)
+        print(sol)

examples/mip/sudoku1

@@ -0,0 +1,9 @@
+.284763..
+...839.2.
+7..512.8.
+..179..4.
+3........
+..9...1..
+.5..8....
+..692...5
+..2645..8

setup.cfg

@@ -27,7 +27,10 @@ universal=1
 [extras]
 contrib = 
     matplotlib >= 2.0
+symmip =
+    pyscipopt
 # all should be a complete list of all dependencies of all other extras. How to
 # automate this?
 all =
     matplotlib >= 2.0
+    pyscipopt

symfit/api.py

@@ -15,5 +15,11 @@
 from symfit.core.argument import Variable, Parameter
 from symfit.core.support import variables, parameters, D
 
+try:
+    # MIP is an optional feature. If no solvers are installed this will raise an import error.
+    from symfit.symmip import MIP
+except ImportError:
+    pass
+
 # Expose the sympy API
 from sympy import *

symfit/core/argument.py

@@ -6,6 +6,8 @@
 import numbers
 import warnings
 
+import numpy as np
+
 from sympy.core.symbol import Symbol
 
 
@@ -93,6 +95,9 @@ class Parameter(Argument):
     _argument_name = 'par'
 
     def __new__(cls, name=None, value=1.0, min=None, max=None, fixed=False, **kwargs):
+        if 'binary' in kwargs:
+            kwargs['integer'] = kwargs['real'] = kwargs['nonnegative'] = True
+
         try:
             return super(Parameter, cls).__new__(cls, name, **kwargs)
         except TypeError as err:
@@ -115,13 +120,17 @@ def __init__(self, name=None, value=1.0, min=None, max=None, fixed=False, **assu
         self.value = value
         self.fixed = fixed
 
-        if min is not None and max is not None and min > max:
-            if not self.fixed:
+        if min is not None and max is not None:
+            test = min > max
+            if isinstance(test, np.ndarray):
+                test = test.any()
+
+            if test and not self.fixed:
                 raise ValueError('The value of `min` should be less than or'
                                  ' equal to the value of `max`.')
-        else:
-            self.min = min
-            self.max = max
+
+        self.min = min
+        self.max = max
 
     def __eq__(self, other):
         """
@@ -135,10 +144,14 @@ def __eq__(self, other):
         if not equal:
             return False
         else:
-            return (self.min == other.min and
-                    self.max == other.max and
+            min_eq = self.min == other.min
+            max_eq = self.max == other.max
+            value_eq = self.value == other.value
+
+            return (min_eq.all() if isinstance(min_eq, np.ndarray) else min_eq and
+                    max_eq.all() if isinstance(max_eq, np.ndarray) else max_eq and
                     self.fixed == other.fixed and
-                    self.value == other.value)
+                    value_eq.all() if isinstance(value_eq, np.ndarray) else value_eq)
 
     __hash__ = Argument.__hash__
 

symfit/core/models.py

@@ -198,8 +198,7 @@ def as_constraint(cls, constraint, model, constraint_type=None, **init_kwargs):
         if set(instance.params) <= set(model.params):
             instance.params = model.params
         else:
-            raise ModelError('The parameters of ``constraint`` have to be a '
-                             'subset of those of ``model``.')
+            model.params = sorted(set(model.params) | set(instance.params), key=lambda x: x.name)
 
         return instance
 

symfit/core/support.py

@@ -302,7 +302,7 @@ def key2str(target):
     :param target: `Mapping` to be made save
     :return: `Mapping` of str(symbol): value pairs.
     """
-    return target.__class__((str(symbol), value) for symbol, value in target.items())
+    return {str(symbol): value for symbol, value in target.items()}
 
 
 def D(*args, **kwargs):

symfit/symmip/__init__.py

@@ -0,0 +1 @@
+from .mip import MIP

symfit/symmip/backend/__init__.py

@@ -0,0 +1,7 @@
+try:
+    # As a commercial solver, gurobi is optional.
+    from .gurobi import GurobiBackend
+except ImportError:
+    pass
+
+from .scipopt import SCIPOptBackend