diff --git a/.github/workflows/test_symmip.yaml b/.github/workflows/test_symmip.yaml new file mode 100644 index 00000000..972165e0 --- /dev/null +++ b/.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 diff --git a/README.rst b/README.rst index 4f54295d..5c0074bc 100644 --- a/README.rst +++ b/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 diff --git a/docs/installation.rst b/docs/installation.rst index 41248260..a0c915d6 100644 --- a/docs/installation.rst +++ b/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 diff --git a/examples/mip/bilinear.py b/examples/mip/bilinear.py new file mode 100644 index 00000000..d782ded1 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/examples/mip/mip1.py b/examples/mip/mip1.py new file mode 100644 index 00000000..fa7bf086 --- /dev/null +++ b/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') diff --git a/examples/mip/multiscenario.py b/examples/mip/multiscenario.py new file mode 100644 index 00000000..1c4404fb --- /dev/null +++ b/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) diff --git a/examples/mip/sudoku.py b/examples/mip/sudoku.py new file mode 100644 index 00000000..c9ed24d4 --- /dev/null +++ b/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) diff --git a/examples/mip/sudoku1 b/examples/mip/sudoku1 new file mode 100644 index 00000000..dcc2098f --- /dev/null +++ b/examples/mip/sudoku1 @@ -0,0 +1,9 @@ +.284763.. +...839.2. +7..512.8. +..179..4. +3........ +..9...1.. +.5..8.... +..692...5 +..2645..8 diff --git a/examples/mip/sudoku_alt.py b/examples/mip/sudoku_alt.py new file mode 100644 index 00000000..e4c553bd --- /dev/null +++ b/examples/mip/sudoku_alt.py @@ -0,0 +1,64 @@ +# Inspired by https://www.gurobi.com/documentation/9.5/examples/sudoku_py.html +# Almost identical to sudoku.py, except this example uses constraints +# to fix the known datapoints rather than bounds. This demonstrates +# how to provide constraints to only a subset of variables. + +import math + +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) + +# 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)) +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. +# 5. Fix known values. +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)], + *[Eq(x[val_i, val_j, int(char) - 1], 1) + for val_i, line in enumerate(grid) + for val_j, char in enumerate(line) + if char != "."] +] + +# 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) diff --git a/setup.cfg b/setup.cfg index bdf35ed6..aebc6c19 100644 --- a/setup.cfg +++ b/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 \ No newline at end of file diff --git a/symfit/api.py b/symfit/api.py index 7c10f641..12a2a714 100644 --- a/symfit/api.py +++ b/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 * \ No newline at end of file diff --git a/symfit/core/argument.py b/symfit/core/argument.py index ef0bba06..208e4efc 100644 --- a/symfit/core/argument.py +++ b/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,13 @@ 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: + if np.any(min > max) 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 +140,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__ diff --git a/symfit/core/models.py b/symfit/core/models.py index dbbe855e..b2c6c052 100644 --- a/symfit/core/models.py +++ b/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 diff --git a/symfit/core/support.py b/symfit/core/support.py index 58a4376f..1b1dd97d 100644 --- a/symfit/core/support.py +++ b/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): diff --git a/symfit/symmip/__init__.py b/symfit/symmip/__init__.py new file mode 100644 index 00000000..2871c077 --- /dev/null +++ b/symfit/symmip/__init__.py @@ -0,0 +1 @@ +from .mip import MIP \ No newline at end of file diff --git a/symfit/symmip/backend/__init__.py b/symfit/symmip/backend/__init__.py new file mode 100644 index 00000000..18f5225f --- /dev/null +++ b/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 diff --git a/symfit/symmip/backend/gurobi.py b/symfit/symmip/backend/gurobi.py new file mode 100644 index 00000000..22afd4a0 --- /dev/null +++ b/symfit/symmip/backend/gurobi.py @@ -0,0 +1,60 @@ +from dataclasses import dataclass, field + +import gurobipy as gp + +from sympy.printing.pycode import PythonCodePrinter + + +VTYPES = {'binary': 'B', 'integer': 'I', 'real': 'C'} + + +class GurobiPrinter(PythonCodePrinter): + modules = gp + + def _print_Sum(self, expr): + loops = ( + 'for {i} in range({a}, {b}+1)'.format( + i=self._print(i), + a=self._print(a), + b=self._print(b)) + for i, a, b in expr.limits) + return '(gurobipy.quicksum({function} {loops}))'.format( + function=self._print(expr.function), + loops=' '.join(loops)) + + +@dataclass +class GurobiBackend: + model: gp.Model = field(default_factory=gp.Model) + printer: PythonCodePrinter = field(default=GurobiPrinter) + + def add_var(self, *, vtype, **kwargs): + return self.model.addVar(vtype=VTYPES[vtype], **kwargs) + + def add_constr(self, constrexpr): + return self.model.addConstr(constrexpr) + + @property + def objective(self): + return self.model.getObjective() + + @objective.setter + def objective(self, obj): + self.model.setObjective(obj) + + def update(self): + self.model.update() + + def optimize(self): + self.model.optimize() + + @property + def objective_value(self): + return self.model.objVal + + @property + def solution(self): + return {v.VarName: v.X for v in self.model.getVars()} + + def get_value(self, v): + return v.X \ No newline at end of file diff --git a/symfit/symmip/backend/scipopt.py b/symfit/symmip/backend/scipopt.py new file mode 100644 index 00000000..78678ca8 --- /dev/null +++ b/symfit/symmip/backend/scipopt.py @@ -0,0 +1,49 @@ +from dataclasses import dataclass, field + +import pyscipopt as scip + +from sympy.printing.pycode import PythonCodePrinter + + +VTYPES = {'binary': 'B', 'integer': 'I', 'real': 'C'} + + +class SCIPOptPrinter(PythonCodePrinter): + modules = scip + + +@dataclass +class SCIPOptBackend: + model: scip.Model = field(default_factory=scip.Model) + printer: PythonCodePrinter = field(default=SCIPOptPrinter) + + def add_var(self, *, vtype, **kwargs): + return self.model.addVar(vtype=VTYPES[vtype], **kwargs) + + def add_constr(self, constrexpr): + return self.model.addCons(constrexpr) + + @property + def objective(self): + return self.model.getObjective() + + @objective.setter + def objective(self, obj): + self.model.setObjective(obj) + + def update(self): + return + + def optimize(self): + self.model.optimize() + + @property + def objective_value(self): + return self.model.getObjVal() + + @property + def solution(self): + return {v.name: self.model.getVal(v) for v in self.model.getVars()} + + def get_value(self, v): + return self.model.getVal(v) \ No newline at end of file diff --git a/symfit/symmip/mip.py b/symfit/symmip/mip.py new file mode 100644 index 00000000..778e3fb9 --- /dev/null +++ b/symfit/symmip/mip.py @@ -0,0 +1,214 @@ +from dataclasses import dataclass, field +from typing import List, Dict, Protocol +from functools import cached_property, reduce +from collections import defaultdict, namedtuple +import itertools + +from sympy import Expr, Rel, Eq, Idx, Indexed, Symbol, lambdify +from sympy.printing.pycode import PythonCodePrinter +import numpy as np + +from symfit.symmip.backend import SCIPOptBackend +from symfit.core.support import key2str + + +VTYPES = ['binary', 'integer', 'real'] # ToDo: Enum? + + +class MIPBackend(Protocol): + @property + def model(self): + pass + + @property + def printer(self) -> PythonCodePrinter: + pass + + def add_var(self, *args, **kwargs): + pass + + def add_constr(self, *args, **kwargs): + pass + + @property + def objective(self): + pass + + @objective.setter + def objective(self, obj): + pass + + def update(self, *args, **kwargs): + pass + + def optimize(self, *args, **kwargs): + pass + + @property + def objective_value(self): + pass + + @property + def solution(self) -> dict: + pass + + def get_value(self, v): + pass + + +@dataclass +class MIPResult: + best_vals: Dict + objective_value: float + + @property + def base2indexed(self): + res = {v.base.label: v for v in self.best_vals if isinstance(v, Indexed)} + res.update({v.base: v for v in self.best_vals if isinstance(v, Indexed)}) + return res + + def __getitem__(self, v): + if not isinstance(v, Indexed) and v in self.base2indexed: + v = self.base2indexed[v] + return self.best_vals[v] + + + +@dataclass +class MIP: + objective: Expr = field(default=None) + constraints: List[Rel] = field(default_factory=list) + backend: MIPBackend = field(default=SCIPOptBackend) + data: Dict = field(default_factory=dict) + + indices: List[Idx] = field(init=False) + variables: List[Symbol] = field(init=False) + indexed_variables: List[Indexed] = field(init=False) + + def __post_init__(self): + # Prepare constraints + self.constraints = [c if isinstance(c, Rel) else Eq(c, 0) + for c in self.constraints] + free_symbols = reduce(set.union, + (c.free_symbols for c in self.constraints), + self.objective.free_symbols if self.objective else set()) + + self.indices = {s for s in free_symbols if isinstance(s, Idx)} + self.indexed_variables = {s for s in free_symbols + if isinstance(s, Indexed) + if all(isinstance(idx, Idx) for idx in s.indices)} + labels = {s.base.label for s in self.indexed_variables} + self.variables = {s for s in free_symbols if not isinstance(s, (Idx, Indexed)) if s not in labels} + self.backend = self.backend() + self._make_mip_model() + + @property + def dependent_vars(self): + return (self.variables | self.indexed_variables) - self.independent_vars + + @property + def independent_vars(self): + return set(self.data) + + @staticmethod + def param_vtype(p: "Parameter"): + """ Return the gurobi vtype for the parameter `p`. """ + for vtype in VTYPES: + if p.assumptions0.get(vtype, False): + return vtype + return 'real' + + def _make_mip_model(self): + # Dictionary mapping symbolic entities to the corresponding MIP variable + self.mip_vars = defaultdict(dict) + + # If data was provided, then use its value instead of making it a Variable. + for p, data in self.data.items(): + if isinstance(p, Indexed): + self.mip_vars[p.base.label] = data + else: + self.mip_vars[p] = data + + for p in self.dependent_vars: + if not isinstance(p, Indexed): + kwargs = {'name': p.name} + if p.max is not None: + kwargs['ub'] = p.max + if p.min is not None: + kwargs['lb'] = p.min + kwargs['vtype'] = self.param_vtype(p) + self.mip_vars[p] = self.backend.add_var(**kwargs) + continue + + param = p.base.label + for indices in itertools.product(*(range(idx.lower, idx.upper + 1) for idx in p.indices)): + kwargs = {'name': f"{param.name}_{'_'.join(str(i) for i in indices)}"} + if param.max is not None: + kwargs['ub'] = param.max[indices] + if param.min is not None: + kwargs['lb'] = param.min[indices] + if param.value is not None: + try: + kwargs['obj'] = param.value[indices] + except TypeError: + pass + kwargs['vtype'] = self.param_vtype(param) + self.mip_vars[param][indices if len(indices) > 1 else indices[0]] = self.backend.add_var(**kwargs) + + all_vars = key2str(self.mip_vars) + # Translate the objective and constraints from sympy to their equivalent for the MIP solver. + if self.objective: + obj_func = lambdify(self.mip_vars.keys(), self.objective, printer=self.backend.printer, modules=self.backend.printer.modules) + obj = obj_func(**all_vars) + self.backend.objective = obj + + # For constraints the idea is the same as above, but a bit more involved since + # constraints can have free indices. + if self.constraints: + free_indices_for_constraints = defaultdict(list) + for constraint in self.constraints: + # Find the free indices in this constraint, and store them in alphabetical order. + free_indices = tuple(sorted((constraint.free_symbols & self.indices), key=lambda s: s.name)) + free_indices_for_constraints[free_indices].append(constraint) + + # TODO: Group the free indices in such a way that we can iterate over them even more efficiently. + for free_indices, constraints in free_indices_for_constraints.items(): + # Create the callable for these constraints to translate to MIP solver. + # Important: the args to these callables are free indices first, then the rest! + args = free_indices + tuple(self.mip_vars.keys()) + args = tuple(s.name for s in args) # To string, because sympy is a Dummy. + constrs_func = lambdify( + args, constraints, printer=self.backend.printer, modules=self.backend.printer.modules + ) + + # For all free indices, call the created function to make the MIP equivalent statement. + for idxs in itertools.product(*[range(i.lower, i.upper + 1) for i in free_indices]): + constrs = constrs_func(*idxs, **all_vars) + for constr in constrs: + self.backend.add_constr(constr) + + self.backend.update() + + def execute(self, *args, **kwargs) -> MIPResult: + self.backend.optimize(*args, **kwargs) + + # Extract the optimized values corresponding to the dependent variables. + dependent_vars = sorted(self.dependent_vars, key=lambda x: x.name) + best_vals = {} + for v in dependent_vars: + if isinstance(v, Indexed): + vtype = self.param_vtype(v.base.label) + if vtype == 'binary': + vals = np.zeros(v.shape, dtype=bool) + elif vtype == 'real': + vals = np.zeros(v.shape, dtype=float) + elif vtype == 'integer': + vals = np.zeros(v.shape, dtype=int) + for idxs in itertools.product(*[range(i.lower, i.upper + 1) for i in v.indices]): + vals[idxs] = self.backend.get_value(self.mip_vars[v.base.label][idxs if len(idxs) > 1 else idxs[0]]) + + best_vals[v] = vals + else: + best_vals[v] = self.backend.get_value(self.mip_vars[v]) + + return MIPResult(objective_value=self.backend.objective_value, best_vals=best_vals) diff --git a/tests/test_mip.py b/tests/test_mip.py new file mode 100644 index 00000000..4d871a97 --- /dev/null +++ b/tests/test_mip.py @@ -0,0 +1,34 @@ +import pytest +import numpy as np + +from symfit import Parameter, parameters, Eq, MIP, Fit + + +def test_parameter_array_bounds(): + lb = np.zeros(10) + ub = np.ones(10) + with pytest.raises(ValueError): + x = Parameter('x', min=ub, max=lb) + x = Parameter('x', min=lb, max=ub) + + +def test_bilinear(): + # Solve the bilinear example with MIP and compare to Fit. + 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() + fit = Fit(- objective, constraints=constraints) + fit_result = fit.execute() + + assert mip_result[x] == pytest.approx(fit_result.value(x), abs=1e-6) + assert mip_result[y] == pytest.approx(fit_result.value(y), abs=1e-6) + assert mip_result[z] == pytest.approx(fit_result.value(z), abs=1e-6) + assert mip_result.objective_value == pytest.approx(fit_result.objective_value, abs=1e-6) \ No newline at end of file