The problem treated here involves assigning disjoint sets of items to a smaller number of collections such that (a) each set goes into one collection in its entirety (i.e., all elements of a set end up in the same collection) and (b) the collection sizes are as close to each other as is possible. The source of the problem is a question on Mathematics Stack Exchange.
Before beginning, we load a few packages. The OMPR package provides a domain specific language (DSL) for writing linear and integer programming models in R. It uses the ROI package to interface to a solver, which here will be CPLEX. (This in turn requires that the Rcplex package be installed.) The dplyr package is used for data manipulation, and the tictoc package for calculating execution times.
library(ROI)ROI: R Optimization Infrastructure
Registered solver plugins: nlminb, cplex.
Default solver: auto.library(ROI.plugin.cplex)
library(ompr)Registered S3 method overwritten by 'data.table':
method from
print.data.table library(ompr.roi)
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, unionlibrary(tictoc)We’ll pick some arbitrary problem dimensions and randomly assign items to sets. Note: The distribution of items across sets is fairly uniform due to the use of the sample() function with default probabilities. Creating a distribution of set sizes with heavier tales might change the performance of any or all of the solution methods shown here.
nItems <- 5723 # total number of items
nSets <- 137 # the number of original sets
nCollections <- 10 # the number of collections to create
# Set a seed for reproducibility.
set.seed(81620)
# Define set memberships for the items.
# To ensure no empty sets, we start with an item in each set and randomly assign the remaining items.
setID <- c(1:nSets, sample(1:nSets, nItems - nSets, replace = TRUE))
# Convert the set IDs into cardinalities for each set.
cardinality <- as.vector(table(setID))To catch possible errors, we will use a function that takes as input a vector of collection assignments (where the i-th vector entry is the index of the collection into which set i is assigned) and outputs the sizes of the collections.
sizeOf <- function(x) {
suppressMessages({ # suppresses an annoying warning
data.frame(c = cardinality, p = x) %>%
group_by(p) %>%
summarize(psize = sum(c)) %>%
pull(psize)
})
}For convenience, we define another function that gives the difference between the largest and smallest sizes of the collections.
spreadOf <- function(x) {
temp <- range(sizeOf(x))
temp[2] - temp[1]
}The greedy heuristic assigns sets in descending size order, putting each set in the currently smallest collection.
tic("greedy")
# Store the collection index of each set in a vector (initially zeroes).
collection.index <- rep(0, nSets)
# Keep a running list of collection sizes.
collection.size <- rep(0, nCollections)
# Get the descending cardinality sort order for the sets.
sort.index <- order(cardinality, decreasing = TRUE)
# Assign sets to collections.
for (i in 1:nSets) {
j <- sort.index[i] # index of the set to be assigned
k <- which.min(collection.size) # index of the currently smallest collection
collection.index[j] <- k
collection.size[k] <- collection.size[k] + cardinality[j] # update the size
}
toc()greedy: 0.038 sec elapsedShow the range of collection sizes.
cat("Collection size range: ", range(collection.size), "\n")Collection size range: 552 582 cat("Spread = ", spreadOf(collection.index), "\n")Spread = 30 Show the individual collection sizes and the spread.
sizeOf(collection.index) [1] 582 582 582 581 580 578 552 552 552 582The pairwise interchange heuristic attempts to shrink the range in collection sizes by swapping a member of the smallest collection (A) with a member of a larger collection (B) so as to increase the size of A without shrinking B so much that it is as small as, or smaller than, the starting size of A.
tic("Improvement heuristic")
retry <- TRUE
# Make a working copy of the first heuristic's solution.
collection.2.index <- collection.index
collection.2.size <- sizeOf(collection.2.index)
spread <- spreadOf(collection.2.index) # to be reduced/minimized
while (retry) {
retry <- FALSE
# Find the smallest collection (i) and its size (s).
i <- which.min(collection.2.size)
i.size <- collection.2.size[i]
# Identify the smallest set in collection i. This is the target for swapping.
in.i <- which(collection.2.index == i) # sets in collection i
g <- in.i[which.min(cardinality[in.i])] # the set to swap
g.size <- cardinality[g] # the size of the set to be moved
# Identify the indices of larger collections and sort them into descending size order.
bigger <- which(collection.2.size > i.size)
temp <- collection.2.size[bigger]
temp <- order(temp, decreasing = TRUE)
bigger <- bigger[temp]
# Process the larger collections in descending size order until a valid swap is found.
for (j in bigger) {
j.size <- collection.2.size[j]
# Get the difference in collection sizes, which is a bound on the difference in size between the two sets being swapped.
bd <- collection.2.size[j] - i.size # bound on swap size
# Get the indices of the sets in collection j.
available <- which(collection.2.index == j)
# Select the largest set that is bigger than g and but not bigger than g plus the bound on the swap difference.
temp <- cardinality[available]
available <- available[which(temp > g.size & temp < g.size + bd)]
# Check whether any set is eligible to be swapped.
if (length(available) >= 1) {
# Find the largest eligible set.
g2 <- available[which.max(cardinality[available])]
g2.size <- cardinality[g2]
# Swap g and g2.
collection.2.index[g] <- j
collection.2.index[g2] <- i
collection.2.size[i] <- i.size - g.size + g2.size
collection.2.size[j] <- j.size + g.size - g2.size
retry <- TRUE
break # Break out of the for loop and start over, looking for another swap.
}
}
}
toc()Improvement heuristic: 0.111 sec elapsedcat("Improvement heuristic spread = ", spreadOf(collection.2.index), "\n")Improvement heuristic spread = 3 sizeOf(collection.2.index) [1] 572 572 574 573 572 571 572 572 572 573Assign sets to collections to minimize the range of collection sizes, using a mixed integer linear program.
tic("IP")
model <-
MIPModel() %>%
# x[i, j] = 1 iff set i is assigned to collection j.
add_variable(x[i, j], i = 1:nSets, j = 1:nCollections, type = "binary") %>%
# size[j] is the size of collection j.
add_variable(size[j], j = 1:nCollections, lb = 0, ub = nItems) %>%
# least and most are the sizes of the smallest and largest collections.
add_variable(least, lb = 0, ub = nItems) %>%
add_variable(most, lb = 0, ub = nItems) %>%
# The objective is to minimize the spread between smallest and largest collections.
set_objective(most - least, "min") %>%
# Every set must be assigned to exactly one collection.
add_constraint(sum_expr(x[i, j], j = 1:nCollections) == 1, i = 1:nSets) %>%
# The size of a collection is the sum of the cardinalities of the member sets.
add_constraint(sum_expr(cardinality[i] * x[i, j], i = 1:nSets) == size[j], j = 1:nCollections) %>%
# Every collection size is between least and most.
add_constraint(size[j] <= most, j = 1:nCollections) %>%
add_constraint(size[j] >= least, j = 1:nCollections) %>%
# To mitigate symmetry with respect to the indexing of the collections, we require the collection sizes to be nondecreasing.
add_constraint(size[j] >= size[j - 1], j = 2:nCollections)
adding constraints [====>-------------------------------] 14% eta 1s
adding constraints [====>-------------------------------] 15% eta 1s
adding constraints [=====>------------------------------] 15% eta 1s
adding constraints [=====>------------------------------] 16% eta 1s
adding constraints [=====>------------------------------] 17% eta 1s
adding constraints [=====>------------------------------] 18% eta 1s
adding constraints [======>-----------------------------] 18% eta 1s
adding constraints [======>-----------------------------] 19% eta 1s
adding constraints [======>-----------------------------] 20% eta 1s
adding constraints [=======>----------------------------] 21% eta 1s
adding constraints [=======>----------------------------] 22% eta 1s
adding constraints [=======>----------------------------] 23% eta 1s
adding constraints [========>---------------------------] 24% eta 1s
adding constraints [========>---------------------------] 25% eta 1s
adding constraints [========>---------------------------] 26% eta 1s
adding constraints [=========>--------------------------] 27% eta 1s
adding constraints [=========>--------------------------] 28% eta 1s
adding constraints [==========>-------------------------] 29% eta 1s
adding constraints [==========>-------------------------] 30% eta 1s
adding constraints [==========>-------------------------] 31% eta 1s
adding constraints [===========>------------------------] 32% eta 1s
adding constraints [===========>------------------------] 33% eta 1s
adding constraints [===========>------------------------] 34% eta 1s
adding constraints [============>-----------------------] 35% eta 1s
adding constraints [============>-----------------------] 36% eta 1s
adding constraints [============>-----------------------] 37% eta 1s
adding constraints [=============>----------------------] 38% eta 1s
adding constraints [=============>----------------------] 39% eta 1s
adding constraints [=============>----------------------] 40% eta 1s
adding constraints [==============>---------------------] 41% eta 1s
adding constraints [==============>---------------------] 42% eta 1s
adding constraints [===============>--------------------] 43% eta 1s
adding constraints [===============>--------------------] 44% eta 1s
adding constraints [===============>--------------------] 45% eta 1s
adding constraints [================>-------------------] 46% eta 1s
adding constraints [================>-------------------] 47% eta 1s
adding constraints [================>-------------------] 48% eta 1s
adding constraints [=================>------------------] 49% eta 1s
adding constraints [=================>------------------] 50% eta 1s
adding constraints [=================>------------------] 51% eta 1s
adding constraints [==================>-----------------] 52% eta 0s
adding constraints [==================>-----------------] 53% eta 0s
adding constraints [==================>-----------------] 54% eta 0s
adding constraints [===================>----------------] 55% eta 0s
adding constraints [===================>----------------] 56% eta 0s
adding constraints [===================>----------------] 57% eta 0s
adding constraints [====================>---------------] 58% eta 0s
adding constraints [====================>---------------] 59% eta 0s
adding constraints [=====================>--------------] 60% eta 0s
adding constraints [=====================>--------------] 61% eta 0s
adding constraints [=====================>--------------] 62% eta 0s
adding constraints [======================>-------------] 63% eta 0s
adding constraints [======================>-------------] 64% eta 0s
adding constraints [======================>-------------] 65% eta 0s
adding constraints [=======================>------------] 66% eta 0s
adding constraints [=======================>------------] 67% eta 0s
adding constraints [=======================>------------] 68% eta 0s
adding constraints [========================>-----------] 69% eta 0s
adding constraints [========================>-----------] 70% eta 0s
adding constraints [========================>-----------] 71% eta 0s
adding constraints [=========================>----------] 72% eta 0s
adding constraints [=========================>----------] 73% eta 0s
adding constraints [==========================>---------] 74% eta 0s
adding constraints [==========================>---------] 75% eta 0s
adding constraints [==========================>---------] 76% eta 0s
adding constraints [===========================>--------] 77% eta 0s
adding constraints [===========================>--------] 78% eta 0s
adding constraints [===========================>--------] 79% eta 0s
adding constraints [============================>-------] 80% eta 0s
adding constraints [============================>-------] 81% eta 0s
adding constraints [============================>-------] 82% eta 0s
adding constraints [=============================>------] 82% eta 0s
adding constraints [=============================>------] 83% eta 0s
adding constraints [=============================>------] 84% eta 0s
adding constraints [=============================>------] 85% eta 0s
adding constraints [==============================>-----] 85% eta 0s
adding constraints [==============================>-----] 86% eta 0s
adding constraints [==============================>-----] 87% eta 0s
adding constraints [===============================>----] 88% eta 0s
adding constraints [===============================>----] 89% eta 0s
adding constraints [===============================>----] 90% eta 0s
adding constraints [================================>---] 91% eta 0s
adding constraints [================================>---] 92% eta 0s
adding constraints [================================>---] 93% eta 0s
adding constraints [=================================>--] 93% eta 0s
adding constraints [=================================>--] 94% eta 0s
adding constraints [=================================>--] 95% eta 0s
adding constraints [=================================>--] 96% eta 0s
adding constraints [==================================>-] 96% eta 0s
adding constraints [==================================>-] 97% eta 0s
adding constraints [==================================>-] 98% eta 0s
adding constraints [==================================>-] 99% eta 0s
adding constraints [===================================>] 99% eta 0s
adding constraints [====================================] 100% eta 0s
adding constraints [==========>-------------------------] 30% eta 1s
adding constraints [=============>----------------------] 40% eta 1s
adding constraints [=================>------------------] 50% eta 0s
adding constraints [=====================>--------------] 60% eta 0s
adding constraints [========================>-----------] 70% eta 0s
adding constraints [============================>-------] 80% eta 0s
adding constraints [===============================>----] 90% eta 0s
adding constraints [====================================] 100% eta 0s
Solve the model.
modelObject <<- as_ROI_model(model)
solution <- model %>% solve_model(with_ROI(solver = "cplex"))CPLEX environment opened
Closed CPLEX environmenttoc()IP: 20.943 sec elapsedDisplay the sizes of the smallest and largest collections.
cat("Smallest collection size = ", get_solution(solution, least), "\n")Smallest collection size = 572 cat("Largest collection size = ", get_solution(solution, most), "\n")Largest collection size = 573 Get the collection assignments as a vector.
temp <-
get_solution(solution, x[i, j]) %>%
filter(value == 1) %>%
arrange(i) %>%
select(j)
collection.3.index <- temp$jShow the individual collection sizes and their spread.
sizeOf(collection.3.index) [1] 572 572 572 572 572 572 572 573 573 573spreadOf(collection.3.index)[1] 1