This notebook explores a question posed on OR Stack Exchange about grouping machines so as to make the group productivity values approximately equal.

As always, the first step is to load libraries.

library(GA, quietly = TRUE)                # for running GAs
  ____    _    
 / ___|  / \     Genetic 
| |  _  / _ \    Algorithms
| |_| |/ ___ \   
 \____/_/   \_\  version 3.2.2
Type 'citation("GA")' for citing this R package in publications.

Attaching package: ‘GA’

The following object is masked from ‘package:utils’:

    de
library(ROI, quietly = TRUE)               # general solver interface
ROI: R Optimization Infrastructure
Registered solver plugins: nlminb, cplex.
Default solver: auto.
library(ROI.plugin.cplex, quietly = TRUE)  # connects the ROI package to CPLEX
library(ompr, quietly = TRUE)              # MIP modeling package
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
library(ompr.roi, quietly = TRUE)          # connects OMPR to ROI
library(dplyr, quietly = TRUE)             # for data frame operations

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union
library(ggplot2, quietly = TRUE)           # to plot results

Data

We start by generating some random data, using dimensions from the original post.

M <- 21        # number of machines
G <- 7         # number of groups
S <- M / G     # size of a group
set.seed(246)  # seed the random number generator
productivity <- runif(M, min = 0, max = 10)  # productivity values
pAvg <- S * mean(productivity)               # mean group productivity

Solutions

A solution will be a permutation of the sequence \(1,\dots,M\), where the first \(S\) elements will form the first machine group, the next \(S\) elements will form the second group, etc.

For convenience, we will create an index vector indicating the group for each element in a permutation.

gIndex <- rep(1:G, each = S)

Criteria

We will seek to minimize either the range or the mean absolute deviation (MAD) of the group productivity levels. Other possible criteria include minimizing the mean squared deviation (MSD, which is likely to produce results similar but not necessarily identical to those using MAD), minimizing the highest group productivity and maximizing the lowest group productivity.

We first define a helper function that converts a solution into a vector of group productivity values.

# Input: a permutation of the machine indices 1,...,M
# Output: a vector of dimension G, where each entry is the productivity of a group
groupVals <- function(x) {
  productivity[x] |> rowsum(group = gIndex, reorder = F)
}

Next, we define functions that calculate the range and mean absolute deviation of a solution.

# Input: a permutation of the machine indices 1,...,M
# Output: the range of the group scores
pRange <- function(x) x |> groupVals() |> range() |> diff()

# Input: a permutation of the machine indices 1,...,M
# Output: the mean absolute deviation of group scores from the overall mean
pMAD <- function(x) {
  # Calculate the group values.
  y <- x |> groupVals()
  # Subtract the overall average, take absolute values and sum.
  (y - pAvg) |> abs() |> sum() / length(x)
}

We will use a data frame to store the group productivity values obtained in each solution, and another data frame to store the range and MAD for each solution.

results <- data.frame("Method" = NA, "Criterion" = NA, "Group" = NA, "Productivity" = NA)
metrics <- data.frame("Method" = NA, "Criterion" = NA, "Range" = NA, "MAD" = NA)

Author’s heuristic

The author suggests a greedy heuristic in which the \(G\) most productive machines are assigned to groups of size 1 and thereafter machines are assigned to the least productive group with space available.

author <- function() {
  p <- rep(0, G)   # group productivity values
  s <- rep(0, G)   # group sizes
  y <- rep(NA, M)  # y[i] is group index for machine i
  # Process machines in decreasing productivity order.
  for (i in sort(productivity, decreasing = T, index.return = T)$ix) {
    # Find the least productive group.
    g <- which.min(p)
    # Assign machine i to group g.
    y[i] <- g
    # Update the group size and productivity.
    s[g] <- s[g] + 1
    p[g] <- p[g] + productivity[i]
    # If the group is full, set its productivity to NA, preventing it from further selection.
    if (s[g] == S) {
      p[g] <- NA
    }
  }
  # Translate the group index vector into a permutation.
  (1:M)[sort(y, index.return = T)$ix]
}

We test the heuristic on the sample data.

sol <- author()
r <- pRange(sol)
m <- pMAD(sol)
results <- results |> add_row(Method = "Author", Criterion = "None", Group = 1:7, Productivity = groupVals(sol))
metrics <- metrics |> add_row(Method = "Author", Criterion = "None", Range = r, MAD = m)
# Summarize it.
cat("Author's heuristic solution has range = ", r, " and  MAD = ", m, ".")
Author's heuristic solution has range =  2.034865  and  MAD =  0.256263 .

Genetic algorithms

As one comparison, we try a genetic algorithms. Since the GA is designed to maximize rather than minimize, we use as fitness value the difference between an upper bound on either range or MAD and the actual range or MAD. The upper bound for range is simply range in machine productivity values multiplied by the group size. The upper bound for MAD is the product of group size with the larger of the absolute difference between the overall average and either the highest or lowest productivity of any machine.

rangeUB <-  productivity |> range() |> diff() * S
madUB <- S * max(max(productivity) - pAvg, pAvg - min(productivity))

The first GA minimizes range. We will use default values for some parameters, but increase the iteration limit and population size.

ga <- ga(type = "permutation",
         fitness = function(x) rangeUB - pRange(x),
         lower = rep.int(1, M),
         upper = rep.int(M, M),
         popSize = 200,
         maxiter = 500,
         monitor = FALSE)
# Extract the solution.
sol <- ga@solution[1,]
r <- pRange(sol)
m <- pMAD(sol)
# Store the results.
results <- results |> add_row(Method = "GA", Criterion = "Range", Group = 1:7, Productivity = groupVals(sol))
metrics <- metrics |> add_row(Method = "GA", Criterion = "Range", Range = r, MAD = m)
# Summarize it.
cat("The GA solution for minimizing range has range = ", r, " and MAD = ", m, ".")
The GA solution for minimizing range has range =  0.6431715  and MAD =  0.05772201 .

The second GA minimizes MAD.

ga <- ga(type = "permutation",
         fitness = function(x) madUB - pMAD(x),
         lower = rep.int(1, M),
         upper = rep.int(M, M),
         popSize = 200,
         maxiter = 500,
         monitor = FALSE)
# Extract the solution.
sol <- ga@solution[1,]
r <- pRange(sol)
m <- pMAD(sol)
# Store the results.
results <- results |> add_row(Method = "GA", Criterion = "MAD", Group = 1:7, Productivity = groupVals(sol))
metrics <- metrics |> add_row(Method = "GA", Criterion = "MAD", Range = r, MAD = m)
# Summarize it.
cat("The GA solution for minimizing MAD has range = ", r, " and MAD = ", m, ".")
The GA solution for minimizing MAD has range =  0.6639519  and MAD =  0.06434813 .

MIP models

We can find a provably optimal solution using MIP models, which we will code using the ompr package and solve with CPLEX. We start with a model for minimizing range.

mip <- MIPModel() |>
          # Use an indicator variable for each combination of machine and group.
          add_variable(x[m, g], m = 1:M, g = 1:G, type = "binary") |>
          # Use a nonnegative variable for the production level of each group.
          add_variable(p[g], g = 1:G, type = "continuous", lb = 0) |>
          # Use nonnegative variables for the smallest and largest group productivity levels.
          add_variable(zlo, type = "continuous", lb = 0) |>
          add_variable(zhi, type = "continuous", lb = 0) |>
          # Limit every machine to a single group.
          add_constraint(sum_expr(x[m, g], g = 1:G) == 1, m = 1:M) |>
          # Limit every group to S machines.
          add_constraint(sum_expr(x[m, g], m = 1:M) == S, g = 1:G) |>
          # Define the group production levels.
          add_constraint(sum_expr(productivity[m] * x[m, g], m = 1:M) - p[g] == 0
                         , g = 1:G) |>
          # Define the range limit variables.
          add_constraint(p[g] - zhi <= 0, g = 1:G) |>
          add_constraint(p[g] - zlo >= 0, g = 1:G) |>
          # Set the objective.
          set_objective(zhi - zlo, sense = "min")
# Solve the MIP.
mipsol <- mip |> solve_model(with_ROI(solver = "cplex"))
CPLEX environment opened
Closed CPLEX environment
# Decode the solution.
sol <- mipsol |> get_solution(x[m, g]) |> filter(value > 0.9) |> arrange(g) |>
                 select(m) |> as.vector()
# Change from tibble to vector.
sol <- sol$m
r <- pRange(sol)
m <- pMAD(sol)
# Store the results.
results <- results |> add_row(Method = "MIP", Criterion = "Range", Group = 1:7, Productivity = groupVals(sol))
metrics <- metrics |> add_row(Method = "MIP", Criterion = "Range", Range = r, MAD = m)
# Summarize it.
cat("The MILP solution for minimizing range has range = ", r, " and MAD = ", m, ".")
The MILP solution for minimizing range has range =  0.4230259  and MAD =  0.03940527 .

Now we try a MIP model to minimize the sum of absolute deviations (which will minimize MAD, since they differ by a constant factor.)

mip <- MIPModel() |>
          # Use an indicator variable for each combination of machine and group.
          add_variable(x[m, g], m = 1:M, g = 1:G, type = "binary") |>
          # Use a nonnegative variable for the production level of each group.
          add_variable(p[g], g = 1:G, type = "continuous", lb = 0) |>
          # Use a nonnegative variable for the absolute deviation of each group's productivity from the overall mean.
          add_variable(z[g], g = 1:G, type = "continuous", lb = 0) |>
          # Limit every machine to a single group.
          add_constraint(sum_expr(x[m, g], g = 1:G) == 1, m = 1:M) |>
          # Limit every group to S machines.
          add_constraint(sum_expr(x[m, g], m = 1:M) == S, g = 1:G) |>
          # Define the group production levels.
          add_constraint(sum_expr(productivity[m] * x[m, g], m = 1:M) - p[g] == 0
                         , g = 1:G) |>
          # Define the absolute value variables.
          add_constraint(p[g] - pAvg - z[g] <= 0, g = 1:G) |>
          add_constraint(p[g] - pAvg + z[g] >= 0, g = 1:G) |>
          # Set the objective.
          set_objective(sum_expr(z[g], g = 1:G), sense = "min")
# Solve the MIP.
mipsol <- mip |> solve_model(with_ROI(solver = "cplex"))
CPLEX environment opened
Closed CPLEX environment
# Decode the solution.
sol <- mipsol |> get_solution(x[m, g]) |> filter(value > 0.9) |> arrange(g) |>
                 select(m) |> as.vector()
# Change from tibble to vector.
sol <- sol$m
r <- pRange(sol)
m <- pMAD(sol)
# Store the results.
results <- results |> add_row(Method = "MIP", Criterion = "MAD", Group = 1:7, Productivity = groupVals(sol))
metrics <- metrics |> add_row(Method = "MIP", Criterion = "MAD", Range = r, MAD = m)
# Summarize it.
cat("The MILP solution for minimizing MAD has range = ", r, " and MAD = ", m, ".")
The MILP solution for minimizing MAD has range =  0.4950586  and MAD =  0.0385512 .

We can drop the dummy rows used to initialize the data frames.

results <- results |> filter(!is.na(Method))
metrics <- metrics |> filter(!is.na(Method))

Results

The overall performance results are summarized below.

metrics

Percentage optimality gaps are as follows.

metrics |>
  mutate(Range = 100 * (Range / min(Range) - 1)) |>
  rename(`Range % Gap` = Range) |>
  mutate(MAD = 100 * (MAD / min(MAD) - 1)) |>
  rename(`MAD % Gap` = MAD)

The author’s heuristic, while fast, is far from optimal in multiple test runs. The GAs also tend to have substantial gaps, although increasing population size and iteration limits might help them.

Finally, we can visualize the group productivity levels of the various solutions.

ggplot(results, aes(x = Group, y = Criterion)) +
  facet_grid(Method ~ ., scales = "free") +
  geom_point(size = 100 * abs(results$Productivity - pAvg) / pAvg) +
  scale_x_continuous(breaks = 1:G) +
  ggtitle("% Deviation from Group Mean") +
  theme(plot.title = element_text(hjust = 0.5))

LS0tCnRpdGxlOiAiTWFjaGluZSBHcm91cHMiCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KClRoaXMgbm90ZWJvb2sgZXhwbG9yZXMgYSBbcXVlc3Rpb25dKGh0dHBzOi8vb3Iuc3RhY2tleGNoYW5nZS5jb20vcXVlc3Rpb25zLzc2MTEvZGl2aWRpbmctbWFjaGluZXMtaW50by1ncm91cHMtb2YtZXF1YWwtc2l6ZXMtc28tdGhhdC1lYWNoLWdyb3VwLWhhcy1hcHByb3hpbWF0ZWwpIHBvc2VkIG9uIE9SIFN0YWNrIEV4Y2hhbmdlIGFib3V0IGdyb3VwaW5nIG1hY2hpbmVzIHNvIGFzIHRvIG1ha2UgdGhlIGdyb3VwIHByb2R1Y3Rpdml0eSB2YWx1ZXMgYXBwcm94aW1hdGVseSBlcXVhbC4KCkFzIGFsd2F5cywgdGhlIGZpcnN0IHN0ZXAgaXMgdG8gbG9hZCBsaWJyYXJpZXMuCgpgYGB7cn0KbGlicmFyeShHQSwgcXVpZXRseSA9IFRSVUUpICAgICAgICAgICAgICAgICMgZm9yIHJ1bm5pbmcgR0FzCmxpYnJhcnkoUk9JLCBxdWlldGx5ID0gVFJVRSkgICAgICAgICAgICAgICAjIGdlbmVyYWwgc29sdmVyIGludGVyZmFjZQpsaWJyYXJ5KFJPSS5wbHVnaW4uY3BsZXgsIHF1aWV0bHkgPSBUUlVFKSAgIyBjb25uZWN0cyB0aGUgUk9JIHBhY2thZ2UgdG8gQ1BMRVgKbGlicmFyeShvbXByLCBxdWlldGx5ID0gVFJVRSkgICAgICAgICAgICAgICMgTUlQIG1vZGVsaW5nIHBhY2thZ2UKbGlicmFyeShvbXByLnJvaSwgcXVpZXRseSA9IFRSVUUpICAgICAgICAgICMgY29ubmVjdHMgT01QUiB0byBST0kKbGlicmFyeShkcGx5ciwgcXVpZXRseSA9IFRSVUUpICAgICAgICAgICAgICMgZm9yIGRhdGEgZnJhbWUgb3BlcmF0aW9ucwpsaWJyYXJ5KGdncGxvdDIsIHF1aWV0bHkgPSBUUlVFKSAgICAgICAgICAgIyB0byBwbG90IHJlc3VsdHMKYGBgCgojIERhdGEKCldlIHN0YXJ0IGJ5IGdlbmVyYXRpbmcgc29tZSByYW5kb20gZGF0YSwgdXNpbmcgZGltZW5zaW9ucyBmcm9tIHRoZSBvcmlnaW5hbCBwb3N0LgoKYGBge3J9Ck0gPC0gMjEgICAgICAgICMgbnVtYmVyIG9mIG1hY2hpbmVzCkcgPC0gNyAgICAgICAgICMgbnVtYmVyIG9mIGdyb3VwcwpTIDwtIE0gLyBHICAgICAjIHNpemUgb2YgYSBncm91cApzZXQuc2VlZCgyNDYpICAjIHNlZWQgdGhlIHJhbmRvbSBudW1iZXIgZ2VuZXJhdG9yCnByb2R1Y3Rpdml0eSA8LSBydW5pZihNLCBtaW4gPSAwLCBtYXggPSAxMCkgICMgcHJvZHVjdGl2aXR5IHZhbHVlcwpwQXZnIDwtIFMgKiBtZWFuKHByb2R1Y3Rpdml0eSkgICAgICAgICAgICAgICAjIG1lYW4gZ3JvdXAgcHJvZHVjdGl2aXR5CmBgYAoKIyBTb2x1dGlvbnMKCkEgc29sdXRpb24gd2lsbCBiZSBhIHBlcm11dGF0aW9uIG9mIHRoZSBzZXF1ZW5jZSAkMSxcZG90cyxNJCwgd2hlcmUgdGhlIGZpcnN0ICRTJCBlbGVtZW50cyB3aWxsIGZvcm0gdGhlIGZpcnN0IG1hY2hpbmUgZ3JvdXAsIHRoZSBuZXh0ICRTJCBlbGVtZW50cyB3aWxsIGZvcm0gdGhlIHNlY29uZCBncm91cCwgZXRjLgoKRm9yIGNvbnZlbmllbmNlLCB3ZSB3aWxsIGNyZWF0ZSBhbiBpbmRleCB2ZWN0b3IgaW5kaWNhdGluZyB0aGUgZ3JvdXAgZm9yIGVhY2ggZWxlbWVudCBpbiBhIHBlcm11dGF0aW9uLgoKYGBge3J9CmdJbmRleCA8LSByZXAoMTpHLCBlYWNoID0gUykKYGBgCgojIENyaXRlcmlhCgpXZSB3aWxsIHNlZWsgdG8gbWluaW1pemUgZWl0aGVyIHRoZSByYW5nZSBvciB0aGUgbWVhbiBhYnNvbHV0ZSBkZXZpYXRpb24gKE1BRCkgb2YgdGhlIGdyb3VwIHByb2R1Y3Rpdml0eSBsZXZlbHMuIE90aGVyIHBvc3NpYmxlIGNyaXRlcmlhIGluY2x1ZGUgbWluaW1pemluZyB0aGUgbWVhbiBzcXVhcmVkIGRldmlhdGlvbiAoTVNELCB3aGljaCBpcyBsaWtlbHkgdG8gcHJvZHVjZSByZXN1bHRzIHNpbWlsYXIgYnV0IG5vdCBuZWNlc3NhcmlseSBpZGVudGljYWwgdG8gdGhvc2UgdXNpbmcgTUFEKSwgbWluaW1pemluZyB0aGUgaGlnaGVzdCBncm91cCBwcm9kdWN0aXZpdHkgYW5kIG1heGltaXppbmcgdGhlIGxvd2VzdCBncm91cCBwcm9kdWN0aXZpdHkuCgpXZSBmaXJzdCBkZWZpbmUgYSBoZWxwZXIgZnVuY3Rpb24gdGhhdCBjb252ZXJ0cyBhIHNvbHV0aW9uIGludG8gYSB2ZWN0b3Igb2YgZ3JvdXAgcHJvZHVjdGl2aXR5IHZhbHVlcy4KCmBgYHtyfQojIElucHV0OiBhIHBlcm11dGF0aW9uIG9mIHRoZSBtYWNoaW5lIGluZGljZXMgMSwuLi4sTQojIE91dHB1dDogYSB2ZWN0b3Igb2YgZGltZW5zaW9uIEcsIHdoZXJlIGVhY2ggZW50cnkgaXMgdGhlIHByb2R1Y3Rpdml0eSBvZiBhIGdyb3VwCmdyb3VwVmFscyA8LSBmdW5jdGlvbih4KSB7CiAgcHJvZHVjdGl2aXR5W3hdIHw+IHJvd3N1bShncm91cCA9IGdJbmRleCwgcmVvcmRlciA9IEYpCn0KYGBgCgpOZXh0LCB3ZSBkZWZpbmUgZnVuY3Rpb25zIHRoYXQgY2FsY3VsYXRlIHRoZSByYW5nZSBhbmQgbWVhbiBhYnNvbHV0ZSBkZXZpYXRpb24gb2YgYSBzb2x1dGlvbi4KCmBgYHtyfQojIElucHV0OiBhIHBlcm11dGF0aW9uIG9mIHRoZSBtYWNoaW5lIGluZGljZXMgMSwuLi4sTQojIE91dHB1dDogdGhlIHJhbmdlIG9mIHRoZSBncm91cCBzY29yZXMKcFJhbmdlIDwtIGZ1bmN0aW9uKHgpIHggfD4gZ3JvdXBWYWxzKCkgfD4gcmFuZ2UoKSB8PiBkaWZmKCkKCiMgSW5wdXQ6IGEgcGVybXV0YXRpb24gb2YgdGhlIG1hY2hpbmUgaW5kaWNlcyAxLC4uLixNCiMgT3V0cHV0OiB0aGUgbWVhbiBhYnNvbHV0ZSBkZXZpYXRpb24gb2YgZ3JvdXAgc2NvcmVzIGZyb20gdGhlIG92ZXJhbGwgbWVhbgpwTUFEIDwtIGZ1bmN0aW9uKHgpIHsKICAjIENhbGN1bGF0ZSB0aGUgZ3JvdXAgdmFsdWVzLgogIHkgPC0geCB8PiBncm91cFZhbHMoKQogICMgU3VidHJhY3QgdGhlIG92ZXJhbGwgYXZlcmFnZSwgdGFrZSBhYnNvbHV0ZSB2YWx1ZXMgYW5kIHN1bS4KICAoeSAtIHBBdmcpIHw+IGFicygpIHw+IHN1bSgpIC8gbGVuZ3RoKHgpCn0KYGBgCgpXZSB3aWxsIHVzZSBhIGRhdGEgZnJhbWUgdG8gc3RvcmUgdGhlIGdyb3VwIHByb2R1Y3Rpdml0eSB2YWx1ZXMgb2J0YWluZWQgaW4gZWFjaCBzb2x1dGlvbiwgYW5kIGFub3RoZXIgZGF0YSBmcmFtZSB0byBzdG9yZSB0aGUgcmFuZ2UgYW5kIE1BRCBmb3IgZWFjaCBzb2x1dGlvbi4KCmBgYHtyfQpyZXN1bHRzIDwtIGRhdGEuZnJhbWUoIk1ldGhvZCIgPSBOQSwgIkNyaXRlcmlvbiIgPSBOQSwgIkdyb3VwIiA9IE5BLCAiUHJvZHVjdGl2aXR5IiA9IE5BKQptZXRyaWNzIDwtIGRhdGEuZnJhbWUoIk1ldGhvZCIgPSBOQSwgIkNyaXRlcmlvbiIgPSBOQSwgIlJhbmdlIiA9IE5BLCAiTUFEIiA9IE5BKQpgYGAKCiMgQXV0aG9yJ3MgaGV1cmlzdGljCgpUaGUgYXV0aG9yIHN1Z2dlc3RzIGEgZ3JlZWR5IGhldXJpc3RpYyBpbiB3aGljaCB0aGUgJEckIG1vc3QgcHJvZHVjdGl2ZSBtYWNoaW5lcyBhcmUgYXNzaWduZWQgdG8gZ3JvdXBzIG9mIHNpemUgMSBhbmQgdGhlcmVhZnRlciBtYWNoaW5lcyBhcmUgYXNzaWduZWQgdG8gdGhlIGxlYXN0IHByb2R1Y3RpdmUgZ3JvdXAgd2l0aCBzcGFjZSBhdmFpbGFibGUuCgpgYGB7cn0KYXV0aG9yIDwtIGZ1bmN0aW9uKCkgewogIHAgPC0gcmVwKDAsIEcpICAgIyBncm91cCBwcm9kdWN0aXZpdHkgdmFsdWVzCiAgcyA8LSByZXAoMCwgRykgICAjIGdyb3VwIHNpemVzCiAgeSA8LSByZXAoTkEsIE0pICAjIHlbaV0gaXMgZ3JvdXAgaW5kZXggZm9yIG1hY2hpbmUgaQogICMgUHJvY2VzcyBtYWNoaW5lcyBpbiBkZWNyZWFzaW5nIHByb2R1Y3Rpdml0eSBvcmRlci4KICBmb3IgKGkgaW4gc29ydChwcm9kdWN0aXZpdHksIGRlY3JlYXNpbmcgPSBULCBpbmRleC5yZXR1cm4gPSBUKSRpeCkgewogICAgIyBGaW5kIHRoZSBsZWFzdCBwcm9kdWN0aXZlIGdyb3VwLgogICAgZyA8LSB3aGljaC5taW4ocCkKICAgICMgQXNzaWduIG1hY2hpbmUgaSB0byBncm91cCBnLgogICAgeVtpXSA8LSBnCiAgICAjIFVwZGF0ZSB0aGUgZ3JvdXAgc2l6ZSBhbmQgcHJvZHVjdGl2aXR5LgogICAgc1tnXSA8LSBzW2ddICsgMQogICAgcFtnXSA8LSBwW2ddICsgcHJvZHVjdGl2aXR5W2ldCiAgICAjIElmIHRoZSBncm91cCBpcyBmdWxsLCBzZXQgaXRzIHByb2R1Y3Rpdml0eSB0byBOQSwgcHJldmVudGluZyBpdCBmcm9tIGZ1cnRoZXIgc2VsZWN0aW9uLgogICAgaWYgKHNbZ10gPT0gUykgewogICAgICBwW2ddIDwtIE5BCiAgICB9CiAgfQogICMgVHJhbnNsYXRlIHRoZSBncm91cCBpbmRleCB2ZWN0b3IgaW50byBhIHBlcm11dGF0aW9uLgogICgxOk0pW3NvcnQoeSwgaW5kZXgucmV0dXJuID0gVCkkaXhdCn0KYGBgCgpXZSB0ZXN0IHRoZSBoZXVyaXN0aWMgb24gdGhlIHNhbXBsZSBkYXRhLgoKYGBge3J9CnNvbCA8LSBhdXRob3IoKQpyIDwtIHBSYW5nZShzb2wpCm0gPC0gcE1BRChzb2wpCnJlc3VsdHMgPC0gcmVzdWx0cyB8PiBhZGRfcm93KE1ldGhvZCA9ICJBdXRob3IiLCBDcml0ZXJpb24gPSAiTm9uZSIsIEdyb3VwID0gMTo3LCBQcm9kdWN0aXZpdHkgPSBncm91cFZhbHMoc29sKSkKbWV0cmljcyA8LSBtZXRyaWNzIHw+IGFkZF9yb3coTWV0aG9kID0gIkF1dGhvciIsIENyaXRlcmlvbiA9ICJOb25lIiwgUmFuZ2UgPSByLCBNQUQgPSBtKQojIFN1bW1hcml6ZSBpdC4KY2F0KCJBdXRob3IncyBoZXVyaXN0aWMgc29sdXRpb24gaGFzIHJhbmdlID0gIiwgciwgIiBhbmQgIE1BRCA9ICIsIG0sICIuIikKYGBgCgojIEdlbmV0aWMgYWxnb3JpdGhtcwoKQXMgb25lIGNvbXBhcmlzb24sIHdlIHRyeSBhIGdlbmV0aWMgYWxnb3JpdGhtcy4gU2luY2UgdGhlIEdBIGlzIGRlc2lnbmVkIHRvIG1heGltaXplIHJhdGhlciB0aGFuIG1pbmltaXplLCB3ZSB1c2UgYXMgZml0bmVzcyB2YWx1ZSB0aGUgZGlmZmVyZW5jZSBiZXR3ZWVuIGFuIHVwcGVyIGJvdW5kIG9uIGVpdGhlciByYW5nZSBvciBNQUQgYW5kIHRoZSBhY3R1YWwgcmFuZ2Ugb3IgTUFELiBUaGUgdXBwZXIgYm91bmQgZm9yIHJhbmdlIGlzIHNpbXBseSByYW5nZSBpbiBtYWNoaW5lIHByb2R1Y3Rpdml0eSB2YWx1ZXMgbXVsdGlwbGllZCBieSB0aGUgZ3JvdXAgc2l6ZS4gVGhlIHVwcGVyIGJvdW5kIGZvciBNQUQgaXMgdGhlIHByb2R1Y3Qgb2YgZ3JvdXAgc2l6ZSB3aXRoIHRoZSBsYXJnZXIgb2YgdGhlIGFic29sdXRlIGRpZmZlcmVuY2UgYmV0d2VlbiB0aGUgb3ZlcmFsbCBhdmVyYWdlIGFuZCBlaXRoZXIgdGhlIGhpZ2hlc3Qgb3IgbG93ZXN0IHByb2R1Y3Rpdml0eSBvZiBhbnkgbWFjaGluZS4KCmBgYHtyfQpyYW5nZVVCIDwtICBwcm9kdWN0aXZpdHkgfD4gcmFuZ2UoKSB8PiBkaWZmKCkgKiBTCm1hZFVCIDwtIFMgKiBtYXgobWF4KHByb2R1Y3Rpdml0eSkgLSBwQXZnLCBwQXZnIC0gbWluKHByb2R1Y3Rpdml0eSkpCmBgYAoKVGhlIGZpcnN0IEdBIG1pbmltaXplcyByYW5nZS4gV2Ugd2lsbCB1c2UgZGVmYXVsdCB2YWx1ZXMgZm9yIHNvbWUgcGFyYW1ldGVycywgYnV0IGluY3JlYXNlIHRoZSBpdGVyYXRpb24gbGltaXQgYW5kIHBvcHVsYXRpb24gc2l6ZS4KCmBgYHtyfQpnYSA8LSBnYSh0eXBlID0gInBlcm11dGF0aW9uIiwKICAgICAgICAgZml0bmVzcyA9IGZ1bmN0aW9uKHgpIHJhbmdlVUIgLSBwUmFuZ2UoeCksCiAgICAgICAgIGxvd2VyID0gcmVwLmludCgxLCBNKSwKICAgICAgICAgdXBwZXIgPSByZXAuaW50KE0sIE0pLAogICAgICAgICBwb3BTaXplID0gMjAwLAogICAgICAgICBtYXhpdGVyID0gNTAwLAogICAgICAgICBtb25pdG9yID0gRkFMU0UpCiMgRXh0cmFjdCB0aGUgc29sdXRpb24uCnNvbCA8LSBnYUBzb2x1dGlvblsxLF0KciA8LSBwUmFuZ2Uoc29sKQptIDwtIHBNQUQoc29sKQojIFN0b3JlIHRoZSByZXN1bHRzLgpyZXN1bHRzIDwtIHJlc3VsdHMgfD4gYWRkX3JvdyhNZXRob2QgPSAiR0EiLCBDcml0ZXJpb24gPSAiUmFuZ2UiLCBHcm91cCA9IDE6NywgUHJvZHVjdGl2aXR5ID0gZ3JvdXBWYWxzKHNvbCkpCm1ldHJpY3MgPC0gbWV0cmljcyB8PiBhZGRfcm93KE1ldGhvZCA9ICJHQSIsIENyaXRlcmlvbiA9ICJSYW5nZSIsIFJhbmdlID0gciwgTUFEID0gbSkKIyBTdW1tYXJpemUgaXQuCmNhdCgiVGhlIEdBIHNvbHV0aW9uIGZvciBtaW5pbWl6aW5nIHJhbmdlIGhhcyByYW5nZSA9ICIsIHIsICIgYW5kIE1BRCA9ICIsIG0sICIuIikKYGBgClRoZSBzZWNvbmQgR0EgbWluaW1pemVzIE1BRC4KCmBgYHtyfQpnYSA8LSBnYSh0eXBlID0gInBlcm11dGF0aW9uIiwKICAgICAgICAgZml0bmVzcyA9IGZ1bmN0aW9uKHgpIG1hZFVCIC0gcE1BRCh4KSwKICAgICAgICAgbG93ZXIgPSByZXAuaW50KDEsIE0pLAogICAgICAgICB1cHBlciA9IHJlcC5pbnQoTSwgTSksCiAgICAgICAgIHBvcFNpemUgPSAyMDAsCiAgICAgICAgIG1heGl0ZXIgPSA1MDAsCiAgICAgICAgIG1vbml0b3IgPSBGQUxTRSkKIyBFeHRyYWN0IHRoZSBzb2x1dGlvbi4Kc29sIDwtIGdhQHNvbHV0aW9uWzEsXQpyIDwtIHBSYW5nZShzb2wpCm0gPC0gcE1BRChzb2wpCiMgU3RvcmUgdGhlIHJlc3VsdHMuCnJlc3VsdHMgPC0gcmVzdWx0cyB8PiBhZGRfcm93KE1ldGhvZCA9ICJHQSIsIENyaXRlcmlvbiA9ICJNQUQiLCBHcm91cCA9IDE6NywgUHJvZHVjdGl2aXR5ID0gZ3JvdXBWYWxzKHNvbCkpCm1ldHJpY3MgPC0gbWV0cmljcyB8PiBhZGRfcm93KE1ldGhvZCA9ICJHQSIsIENyaXRlcmlvbiA9ICJNQUQiLCBSYW5nZSA9IHIsIE1BRCA9IG0pCiMgU3VtbWFyaXplIGl0LgpjYXQoIlRoZSBHQSBzb2x1dGlvbiBmb3IgbWluaW1pemluZyBNQUQgaGFzIHJhbmdlID0gIiwgciwgIiBhbmQgTUFEID0gIiwgbSwgIi4iKQpgYGAKCiMgTUlQIG1vZGVscwoKV2UgY2FuIGZpbmQgYSBwcm92YWJseSBvcHRpbWFsIHNvbHV0aW9uIHVzaW5nIE1JUCBtb2RlbHMsIHdoaWNoIHdlIHdpbGwgY29kZSB1c2luZyB0aGUgYG9tcHJgIHBhY2thZ2UgYW5kIHNvbHZlIHdpdGggQ1BMRVguIFdlIHN0YXJ0IHdpdGggYSBtb2RlbCBmb3IgbWluaW1pemluZyByYW5nZS4KCmBgYHtyfQptaXAgPC0gTUlQTW9kZWwoKSB8PgogICAgICAgICAgIyBVc2UgYW4gaW5kaWNhdG9yIHZhcmlhYmxlIGZvciBlYWNoIGNvbWJpbmF0aW9uIG9mIG1hY2hpbmUgYW5kIGdyb3VwLgogICAgICAgICAgYWRkX3ZhcmlhYmxlKHhbbSwgZ10sIG0gPSAxOk0sIGcgPSAxOkcsIHR5cGUgPSAiYmluYXJ5IikgfD4KICAgICAgICAgICMgVXNlIGEgbm9ubmVnYXRpdmUgdmFyaWFibGUgZm9yIHRoZSBwcm9kdWN0aW9uIGxldmVsIG9mIGVhY2ggZ3JvdXAuCiAgICAgICAgICBhZGRfdmFyaWFibGUocFtnXSwgZyA9IDE6RywgdHlwZSA9ICJjb250aW51b3VzIiwgbGIgPSAwKSB8PgogICAgICAgICAgIyBVc2Ugbm9ubmVnYXRpdmUgdmFyaWFibGVzIGZvciB0aGUgc21hbGxlc3QgYW5kIGxhcmdlc3QgZ3JvdXAgcHJvZHVjdGl2aXR5IGxldmVscy4KICAgICAgICAgIGFkZF92YXJpYWJsZSh6bG8sIHR5cGUgPSAiY29udGludW91cyIsIGxiID0gMCkgfD4KICAgICAgICAgIGFkZF92YXJpYWJsZSh6aGksIHR5cGUgPSAiY29udGludW91cyIsIGxiID0gMCkgfD4KICAgICAgICAgICMgTGltaXQgZXZlcnkgbWFjaGluZSB0byBhIHNpbmdsZSBncm91cC4KICAgICAgICAgIGFkZF9jb25zdHJhaW50KHN1bV9leHByKHhbbSwgZ10sIGcgPSAxOkcpID09IDEsIG0gPSAxOk0pIHw+CiAgICAgICAgICAjIExpbWl0IGV2ZXJ5IGdyb3VwIHRvIFMgbWFjaGluZXMuCiAgICAgICAgICBhZGRfY29uc3RyYWludChzdW1fZXhwcih4W20sIGddLCBtID0gMTpNKSA9PSBTLCBnID0gMTpHKSB8PgogICAgICAgICAgIyBEZWZpbmUgdGhlIGdyb3VwIHByb2R1Y3Rpb24gbGV2ZWxzLgogICAgICAgICAgYWRkX2NvbnN0cmFpbnQoc3VtX2V4cHIocHJvZHVjdGl2aXR5W21dICogeFttLCBnXSwgbSA9IDE6TSkgLSBwW2ddID09IDAKICAgICAgICAgICAgICAgICAgICAgICAgICwgZyA9IDE6RykgfD4KICAgICAgICAgICMgRGVmaW5lIHRoZSByYW5nZSBsaW1pdCB2YXJpYWJsZXMuCiAgICAgICAgICBhZGRfY29uc3RyYWludChwW2ddIC0gemhpIDw9IDAsIGcgPSAxOkcpIHw+CiAgICAgICAgICBhZGRfY29uc3RyYWludChwW2ddIC0gemxvID49IDAsIGcgPSAxOkcpIHw+CiAgICAgICAgICAjIFNldCB0aGUgb2JqZWN0aXZlLgogICAgICAgICAgc2V0X29iamVjdGl2ZSh6aGkgLSB6bG8sIHNlbnNlID0gIm1pbiIpCiMgU29sdmUgdGhlIE1JUC4KbWlwc29sIDwtIG1pcCB8PiBzb2x2ZV9tb2RlbCh3aXRoX1JPSShzb2x2ZXIgPSAiY3BsZXgiKSkKIyBEZWNvZGUgdGhlIHNvbHV0aW9uLgpzb2wgPC0gbWlwc29sIHw+IGdldF9zb2x1dGlvbih4W20sIGddKSB8PiBmaWx0ZXIodmFsdWUgPiAwLjkpIHw+IGFycmFuZ2UoZykgfD4KICAgICAgICAgICAgICAgICBzZWxlY3QobSkgfD4gYXMudmVjdG9yKCkKIyBDaGFuZ2UgZnJvbSB0aWJibGUgdG8gdmVjdG9yLgpzb2wgPC0gc29sJG0KciA8LSBwUmFuZ2Uoc29sKQptIDwtIHBNQUQoc29sKQojIFN0b3JlIHRoZSByZXN1bHRzLgpyZXN1bHRzIDwtIHJlc3VsdHMgfD4gYWRkX3JvdyhNZXRob2QgPSAiTUlQIiwgQ3JpdGVyaW9uID0gIlJhbmdlIiwgR3JvdXAgPSAxOjcsIFByb2R1Y3Rpdml0eSA9IGdyb3VwVmFscyhzb2wpKQptZXRyaWNzIDwtIG1ldHJpY3MgfD4gYWRkX3JvdyhNZXRob2QgPSAiTUlQIiwgQ3JpdGVyaW9uID0gIlJhbmdlIiwgUmFuZ2UgPSByLCBNQUQgPSBtKQojIFN1bW1hcml6ZSBpdC4KY2F0KCJUaGUgTUlMUCBzb2x1dGlvbiBmb3IgbWluaW1pemluZyByYW5nZSBoYXMgcmFuZ2UgPSAiLCByLCAiIGFuZCBNQUQgPSAiLCBtLCAiLiIpCmBgYAoKTm93IHdlIHRyeSBhIE1JUCBtb2RlbCB0byBtaW5pbWl6ZSB0aGUgc3VtIG9mIGFic29sdXRlIGRldmlhdGlvbnMgKHdoaWNoIHdpbGwgbWluaW1pemUgTUFELCBzaW5jZSB0aGV5IGRpZmZlciBieSBhIGNvbnN0YW50IGZhY3Rvci4pCgpgYGB7cn0KbWlwIDwtIE1JUE1vZGVsKCkgfD4KICAgICAgICAgICMgVXNlIGFuIGluZGljYXRvciB2YXJpYWJsZSBmb3IgZWFjaCBjb21iaW5hdGlvbiBvZiBtYWNoaW5lIGFuZCBncm91cC4KICAgICAgICAgIGFkZF92YXJpYWJsZSh4W20sIGddLCBtID0gMTpNLCBnID0gMTpHLCB0eXBlID0gImJpbmFyeSIpIHw+CiAgICAgICAgICAjIFVzZSBhIG5vbm5lZ2F0aXZlIHZhcmlhYmxlIGZvciB0aGUgcHJvZHVjdGlvbiBsZXZlbCBvZiBlYWNoIGdyb3VwLgogICAgICAgICAgYWRkX3ZhcmlhYmxlKHBbZ10sIGcgPSAxOkcsIHR5cGUgPSAiY29udGludW91cyIsIGxiID0gMCkgfD4KICAgICAgICAgICMgVXNlIGEgbm9ubmVnYXRpdmUgdmFyaWFibGUgZm9yIHRoZSBhYnNvbHV0ZSBkZXZpYXRpb24gb2YgZWFjaCBncm91cCdzIHByb2R1Y3Rpdml0eSBmcm9tIHRoZSBvdmVyYWxsIG1lYW4uCiAgICAgICAgICBhZGRfdmFyaWFibGUoeltnXSwgZyA9IDE6RywgdHlwZSA9ICJjb250aW51b3VzIiwgbGIgPSAwKSB8PgogICAgICAgICAgIyBMaW1pdCBldmVyeSBtYWNoaW5lIHRvIGEgc2luZ2xlIGdyb3VwLgogICAgICAgICAgYWRkX2NvbnN0cmFpbnQoc3VtX2V4cHIoeFttLCBnXSwgZyA9IDE6RykgPT0gMSwgbSA9IDE6TSkgfD4KICAgICAgICAgICMgTGltaXQgZXZlcnkgZ3JvdXAgdG8gUyBtYWNoaW5lcy4KICAgICAgICAgIGFkZF9jb25zdHJhaW50KHN1bV9leHByKHhbbSwgZ10sIG0gPSAxOk0pID09IFMsIGcgPSAxOkcpIHw+CiAgICAgICAgICAjIERlZmluZSB0aGUgZ3JvdXAgcHJvZHVjdGlvbiBsZXZlbHMuCiAgICAgICAgICBhZGRfY29uc3RyYWludChzdW1fZXhwcihwcm9kdWN0aXZpdHlbbV0gKiB4W20sIGddLCBtID0gMTpNKSAtIHBbZ10gPT0gMAogICAgICAgICAgICAgICAgICAgICAgICAgLCBnID0gMTpHKSB8PgogICAgICAgICAgIyBEZWZpbmUgdGhlIGFic29sdXRlIHZhbHVlIHZhcmlhYmxlcy4KICAgICAgICAgIGFkZF9jb25zdHJhaW50KHBbZ10gLSBwQXZnIC0geltnXSA8PSAwLCBnID0gMTpHKSB8PgogICAgICAgICAgYWRkX2NvbnN0cmFpbnQocFtnXSAtIHBBdmcgKyB6W2ddID49IDAsIGcgPSAxOkcpIHw+CiAgICAgICAgICAjIFNldCB0aGUgb2JqZWN0aXZlLgogICAgICAgICAgc2V0X29iamVjdGl2ZShzdW1fZXhwcih6W2ddLCBnID0gMTpHKSwgc2Vuc2UgPSAibWluIikKIyBTb2x2ZSB0aGUgTUlQLgptaXBzb2wgPC0gbWlwIHw+IHNvbHZlX21vZGVsKHdpdGhfUk9JKHNvbHZlciA9ICJjcGxleCIpKQojIERlY29kZSB0aGUgc29sdXRpb24uCnNvbCA8LSBtaXBzb2wgfD4gZ2V0X3NvbHV0aW9uKHhbbSwgZ10pIHw+IGZpbHRlcih2YWx1ZSA+IDAuOSkgfD4gYXJyYW5nZShnKSB8PgogICAgICAgICAgICAgICAgIHNlbGVjdChtKSB8PiBhcy52ZWN0b3IoKQojIENoYW5nZSBmcm9tIHRpYmJsZSB0byB2ZWN0b3IuCnNvbCA8LSBzb2wkbQpyIDwtIHBSYW5nZShzb2wpCm0gPC0gcE1BRChzb2wpCiMgU3RvcmUgdGhlIHJlc3VsdHMuCnJlc3VsdHMgPC0gcmVzdWx0cyB8PiBhZGRfcm93KE1ldGhvZCA9ICJNSVAiLCBDcml0ZXJpb24gPSAiTUFEIiwgR3JvdXAgPSAxOjcsIFByb2R1Y3Rpdml0eSA9IGdyb3VwVmFscyhzb2wpKQptZXRyaWNzIDwtIG1ldHJpY3MgfD4gYWRkX3JvdyhNZXRob2QgPSAiTUlQIiwgQ3JpdGVyaW9uID0gIk1BRCIsIFJhbmdlID0gciwgTUFEID0gbSkKIyBTdW1tYXJpemUgaXQuCmNhdCgiVGhlIE1JTFAgc29sdXRpb24gZm9yIG1pbmltaXppbmcgTUFEIGhhcyByYW5nZSA9ICIsIHIsICIgYW5kIE1BRCA9ICIsIG0sICIuIikKYGBgCgpXZSBjYW4gZHJvcCB0aGUgZHVtbXkgcm93cyB1c2VkIHRvIGluaXRpYWxpemUgdGhlIGRhdGEgZnJhbWVzLgoKYGBge3J9CnJlc3VsdHMgPC0gcmVzdWx0cyB8PiBmaWx0ZXIoIWlzLm5hKE1ldGhvZCkpCm1ldHJpY3MgPC0gbWV0cmljcyB8PiBmaWx0ZXIoIWlzLm5hKE1ldGhvZCkpCmBgYAoKIyBSZXN1bHRzCgpUaGUgb3ZlcmFsbCBwZXJmb3JtYW5jZSByZXN1bHRzIGFyZSBzdW1tYXJpemVkIGJlbG93LgoKYGBge3J9Cm1ldHJpY3MKYGBgCgpQZXJjZW50YWdlIG9wdGltYWxpdHkgZ2FwcyBhcmUgYXMgZm9sbG93cy4KCmBgYHtyfQptZXRyaWNzIHw+CiAgbXV0YXRlKFJhbmdlID0gMTAwICogKFJhbmdlIC8gbWluKFJhbmdlKSAtIDEpKSB8PgogIHJlbmFtZShgUmFuZ2UgJSBHYXBgID0gUmFuZ2UpIHw+CiAgbXV0YXRlKE1BRCA9IDEwMCAqIChNQUQgLyBtaW4oTUFEKSAtIDEpKSB8PgogIHJlbmFtZShgTUFEICUgR2FwYCA9IE1BRCkKYGBgCgpUaGUgYXV0aG9yJ3MgaGV1cmlzdGljLCB3aGlsZSBmYXN0LCBpcyBmYXIgZnJvbSBvcHRpbWFsIGluIG11bHRpcGxlIHRlc3QgcnVucy4gVGhlIEdBcyBhbHNvIHRlbmQgdG8gaGF2ZSBzdWJzdGFudGlhbCBnYXBzLCBhbHRob3VnaCBpbmNyZWFzaW5nIHBvcHVsYXRpb24gc2l6ZSBhbmQgaXRlcmF0aW9uIGxpbWl0cyBtaWdodCBoZWxwIHRoZW0uCgpGaW5hbGx5LCB3ZSBjYW4gdmlzdWFsaXplIHRoZSBncm91cCBwcm9kdWN0aXZpdHkgbGV2ZWxzIG9mIHRoZSB2YXJpb3VzIHNvbHV0aW9ucy4KCmBgYHtyfQpnZ3Bsb3QocmVzdWx0cywgYWVzKHggPSBHcm91cCwgeSA9IENyaXRlcmlvbikpICsKICBmYWNldF9ncmlkKE1ldGhvZCB+IC4sIHNjYWxlcyA9ICJmcmVlIikgKwogIGdlb21fcG9pbnQoc2l6ZSA9IDEwMCAqIGFicyhyZXN1bHRzJFByb2R1Y3Rpdml0eSAtIHBBdmcpIC8gcEF2ZykgKwogIHNjYWxlX3hfY29udGludW91cyhicmVha3MgPSAxOkcpICsKICBnZ3RpdGxlKCIlIERldmlhdGlvbiBmcm9tIEdyb3VwIE1lYW4iKSArCiAgdGhlbWUocGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChoanVzdCA9IDAuNSkpCmBgYAoK