Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add functions to find optimal partitions using ILP solvers. #3

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions demos/approximate_algorithms.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
#!python3

"""
Demonstrates how to use approximation algorithms
such as Karmarkar-Karp and Greedy.
"""

from numberpartitioning import karmarkar_karp, greedy, complete_greedy
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could be useful if the examples here were included as docstrings on the relevant functions or modules. This way, they would be easy to include in auto-generated documentation, and it would be possible to run them as tests locally and in CI.

E.g. something like the example sections in the methods in https://github.com/scipy/scipy/blob/master/scipy/sparse/csgraph/_flow.pyx works well.

numbers = [4, 6, 7, 5, 8]
print("Karmarkar-Karp with k=2", karmarkar_karp(numbers, num_parts=2))
print("Karmarkar-Karp with k=3", karmarkar_karp(numbers, num_parts=3))
print("Greedy with k=2", greedy(numbers, num_parts=2))
print("Greedy with k=3", greedy(numbers, num_parts=3))

print("More tests")
print(greedy([1,2,3,4,5,9], num_parts=2))
print(greedy([1,2,3,4,5,9], num_parts=3))
33 changes: 33 additions & 0 deletions demos/complete_algorithms.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
#!python3

"""
Demonstrates how to use complete algorithms
such as Complete Greedy.
"""

from numberpartitioning import complete_greedy

print("Complete Greedy with k=2, default objective:")
for p in complete_greedy([4, 6, 7, 5, 8], num_parts=2):
print (p)

# Define different objectives (all of them for minimization)
def largest_sum(partition):
return max([sum(part) for part in partition])
def minus_smallest_sum(partition):
return -min([sum(part) for part in partition])
def largest_diff(partition):
sums = [sum(part) for part in partition]
return max(sums)-min(sums)

print("An example from Walter (2013), 'Comparing the minimum completion times of two longest-first scheduling-heuristics'.")
walter_numbers = [46, 39, 27, 26, 16, 13, 10]
print("Complete Greedy with k=3, min-diff objective (default):")
for p in complete_greedy(walter_numbers, num_parts=3):
print (p)
print("Complete Greedy with k=3, min-max objective:")
for p in complete_greedy(walter_numbers, num_parts=3, objective=largest_sum):
print (p)
print("Complete Greedy with k=3, max-min objective:")
for p in complete_greedy(walter_numbers, num_parts=3, objective=minus_smallest_sum):
print (p)
13 changes: 13 additions & 0 deletions demos/ilp.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
#!python3

"""
A demo program for finding an optimal partition using ILP.
"""

from numberpartitioning.ilp import *

print("An example from Walter (2013), 'Comparing the minimum completion times of two longest-first scheduling-heuristics'.")
walter_numbers = [46, 39, 27, 26, 16, 13, 10]
print("Min-diff objective:", min_diff_partition(walter_numbers, num_parts=3))
print("Min-max objective:", min_max_partition(walter_numbers, num_parts=3))
print("Max-min objective:", max_min_partition(walter_numbers, num_parts=3))
2 changes: 2 additions & 0 deletions requirements.txt
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
cvxpy

# Dev
build

Expand Down
267 changes: 267 additions & 0 deletions src/numberpartitioning/ilp.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,267 @@
#!python3

from typing import List, Optional, Callable
from numbers import Number

from numberpartitioning.common import Partition, PartitioningResult

import cvxpy
from numberpartitioning.solve import minimize


def max_min_partition(
numbers: List[int], num_parts: int = 2, return_indices: bool = False,
copies = 1, num_smallest_parts:int = 1,
) -> PartitioningResult:
"""
Produce a partition that maximizes the smallest sum,
by solving an integer linear program (ILP).
Credit: Rob Pratt, https://or.stackexchange.com/a/6115/2576
Uses CVXPY as an ILP solver.

Parameters
----------
numbers
The list of numbers to be partitioned.
num_parts
The desired number of parts in the partition. Default: 2.
return_indices
If True, the elements of the parts are the indices of the corresponding entries
of numbers; if False (default), the elements are the numbers themselves.
copies
how many copies are there from each number.
Can be either a single integer (the same #copies for all numbers),
or a list of integers (a different #copies for each number).
Default: 1
num_smallest_parts
number of smallest parts whose sum should be maximized.
Default is 1, which means to just maximize the smallest part-sum.
A value of 3, for example, means to maximize the sum of the three smallest part-sums.

Returns
-------
A partition represented by a ``PartitioningResult``.

>>> max_min_partition([10,20,40,0], num_parts=1, return_indices=True)
PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])

>>> max_min_partition([10,20,40,0], num_parts=2, return_indices=False)
PartitioningResult(partition=[[10, 20, 0], [40]], sizes=[30.0, 40.0])

>>> result = max_min_partition([10,20,40,0], num_parts=3)
>>> int(min(result.sizes))
10

>>> result = max_min_partition([10,20,40,0], num_parts=4)
>>> int(min(result.sizes))
0

>>> result = max_min_partition([10,20,40,0], num_parts=5)
>>> int(min(result.sizes))
0

>>> max_min_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])

>>> max_min_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])

>>> max_min_partition([10,20,40,1], num_parts=3, num_smallest_parts=2)
PartitioningResult(partition=[[], [10, 20, 1], [40]], sizes=[0.0, 31.0, 40.0])
"""
def minimization_objective(parts_sums):
return - sum(parts_sums[0:num_smallest_parts])
return find_optimal_partition(numbers, minimization_objective, num_parts, return_indices, copies)



def min_max_partition(
numbers: List[int], num_parts: int = 2, return_indices: bool = False,
copies = 1, num_largest_parts:int = 1,
) -> PartitioningResult:
"""
Produce a partition that minimizes the largest sum,
by solving an integer linear program (ILP).

Parameters
----------
numbers
The list of numbers to be partitioned.
num_parts
The desired number of parts in the partition. Default: 2.
return_indices
If True, the elements of the parts are the indices of the corresponding entries
of numbers; if False (default), the elements are the numbers themselves.
copies
how many copies are there from each number.
Can be either a single integer (the same #copies for all numbers),
or a list of integers (a different #copies for each number).
Default: 1
num_largest_parts
number of largest parts whose sum should be minimized.
Default is 1, which means to just minimize the largest part-sum.
A value of 3, for example, means to minimize the sum of the three largest part-sums.

Returns
-------
A partition represented by a ``PartitioningResult``.

>>> min_max_partition([10,20,40,0], num_parts=1, return_indices=True)
PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])

>>> min_max_partition([10,20,40,0], num_parts=2, return_indices=False)
PartitioningResult(partition=[[10, 20], [40, 0]], sizes=[30.0, 40.0])

>>> result = min_max_partition([10,20,40,0], num_parts=3)
>>> int(max(result.sizes))
40

>>> result = min_max_partition([10,20,40,0], num_parts=5)
>>> int(max(result.sizes))
40

>>> min_max_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])

>>> min_max_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])

>>> min_max_partition([10,20,40,1], num_parts=3, num_largest_parts=2)
PartitioningResult(partition=[[10, 1], [20], [40]], sizes=[11.0, 20.0, 40.0])
"""
def minimization_objective(parts_sums):
return sum(parts_sums[-num_largest_parts:])
return find_optimal_partition(numbers, minimization_objective, num_parts, return_indices, copies)


def min_diff_partition(
numbers: List[int], num_parts: int = 2, return_indices: bool = False,
copies = 1
) -> PartitioningResult:
"""
Produce a partition that minimizes the difference between the largest and smallest sum,
by solving an integer linear program (ILP).

Parameters
----------
numbers
The list of numbers to be partitioned.
num_parts
The desired number of parts in the partition. Default: 2.
return_indices
If True, the elements of the parts are the indices of the corresponding entries
of numbers; if False (default), the elements are the numbers themselves.
copies
how many copies are there from each number.
Can be either a single integer (the same #copies for all numbers),
or a list of integers (a different #copies for each number).
Default: 1

Returns
-------
A partition represented by a ``PartitioningResult``.

>>> min_diff_partition([10,20,40,0], num_parts=1, return_indices=True)
PartitioningResult(partition=[[0, 1, 2, 3]], sizes=[70.0])

>>> min_diff_partition([10,20,40,0], num_parts=2, return_indices=False)
PartitioningResult(partition=[[10, 20], [40, 0]], sizes=[30.0, 40.0])

>>> min_diff_partition([10,20,40,1], num_parts=2, copies=2, return_indices=True)
PartitioningResult(partition=[[0, 1, 2, 3], [0, 1, 2, 3]], sizes=[71.0, 71.0])

>>> min_diff_partition([10,20,40,0], num_parts=2, copies=[2,1,1,0], return_indices=True)
PartitioningResult(partition=[[0, 0, 1], [2]], sizes=[40.0, 40.0])
"""
def minimization_objective(parts_sums):
return parts_sums[-1] - parts_sums[0]
return find_optimal_partition(numbers, minimization_objective, num_parts, return_indices, copies)


def find_optimal_partition(
numbers: List[int],
minimization_objective: Optional[Callable[[list], float]],
num_parts: int = 2, return_indices: bool = False, copies = 1
) -> PartitioningResult:
"""
Produce a partition that minimizes the given objective,
by solving an integer linear program (ILP).
Credit: Rob Pratt, https://or.stackexchange.com/a/6115/2576
Uses CVXPY as an ILP solver.

Parameters
----------
numbers
The list of numbers to be partitioned.
minimization_objective
A callable for constructing the objective function to be minimized.
Gets as input a list of part-sums, sorted from small to large.
Each of the part-sums is a cvxpy expression.
Returns as output an expression that should be minimized.
See max_min_partition for usage examples.
num_parts
The desired number of parts in the partition. Default: 2.
return_indices
If True, the elements of the parts are the indices of the corresponding entries
of numbers; if False (default), the elements are the numbers themselves.
copies
how many copies are there from each number.
Can be either a single integer (the same #copies for all numbers),
or a list of integers (a different #copies for each number).
Default: 1

Returns
-------
A partition representing by a ``PartitioningResult``.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So this always runs to optimality right? Given the NP-hard nature of the problem, that might be a bit of an ask. Would it be a lot of work to use callbacks to generate solutions along the way, or maybe make it easy to allow it to return non-optimak results (through time constraints, allowed gap, etc.)?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The code uses an existing ILP solver, so - as far as I know - it cannot return solutions along the way. On the positive side, ILP solvers use many tricks that make the solution very efficient in practice.


"""
parts = range(num_parts)
items = range(len(numbers))
if isinstance(copies, Number):
copies = [copies]*len(numbers)

counts:dict = {
item:
[cvxpy.Variable(integer=True) for part in parts]
for item in items
} # counts[i][j] determines how many times item i appears in part j.
parts_sums = [
sum([counts[item][part]*numbers[item] for item in items])
for part in parts]

# Construct the list of constraints:
constraints = []

# The counts must be non-negative:
constraints += [counts[item][part] >= 0 for part in parts for item in items]

# Each item must be in exactly one part:
constraints += [sum([counts[item][part] for part in parts]) == copies[item] for item in items]

# Parts must be in ascending order of their sum (a symmetry-breaker):
constraints += [parts_sums[part+1] >= parts_sums[part] for part in range(num_parts-1)]

objective = minimization_objective(parts_sums)
minimize(objective, constraints)

partition = [
sum([int(counts[item][part].value)*[item]
for item in items if counts[item][part].value>=1], [])
for part in parts
]
sums = [parts_sums[part].value for part in parts]
partition:Partition = [[] for _ in parts]
for part in parts:
for item in items:
count = int(counts[item][part].value)
if count>=1:
partition[part] += count * [item if return_indices else numbers[item]]
return PartitioningResult(partition, sums)




if __name__ == "__main__":
import doctest
(failures,tests) = doctest.testmod(report=True)
print ("{} failures, {} tests".format(failures,tests))
Loading