Excel ANOVA Format in R

A while back, I wrote a Shiny application for a colleague who teaches basic linear regression in a few courses. She recently asked me to add an ANOVA table. I thought the change would be trivial -- R has at least two functions in the stats package, anova() and aov(), for generating ANOVA output -- but there was a catch. Because the students start in Excel before switching to the Shiny app, she wanted me to mimic Excel's output.

The following is Excel's version of an ANOVA table for a small model with independent variables X1 and X2 and dependent variable Y. 

	df	Sum of Squares	Mean Square	F	Significance F
Regression	2	193.4887725	96.74438625	334.348849898829	6.64178639346989E-12
Residual	13	3.76157124999999	0.289351634615384
Total	15	197.25034375

I'm not sure that the column headings are gospel -- my colleague might have tweaked them a bit -- but the important thing is that there are three rows: "Regression" (variation attributable to the model as a whole); "Residual" (variation attributable to the residuals or "errors"); and "Total", variation of the dependent variable around its mean. Contrast that with the output from R's anova() function for the same model.

	Df	Sum Sq	Mean Sq	F value	Pr(>F)
X1	1	183.527	183.527	634.270	2.037e-12	***
X2	1	9.962	9.962	34.428	5.524e-05	***
Residuals	13	3.762	0.289

Note that R splits the explained variation among the individual predictor variables and omits the total variation around the dependent variable mean. (It also omits a whole lot of decimal places in which I would place little faith.)

I search around for an R function that would generate an ANOVA table in something like Excel's format but came up empty. So I wrote my own, which I will share here in case anyone else has the same need. To combine the individual predictor rows into the first row of the Excel version, I just add the degrees of freedom and sum of squares and then recompute the mean square value, F statistic and probability. To create the bottom row of the Excel version, I regress the dependent variable around just a constant term (lm(Y ~ 1)) and then apply the anova() function to that. It might seem a bit roundabout, but it generates an output that is easily grafted onto the rest of the ANOVA table using rbind(). Here's the R code. 

#
# Generate an ANOVA table for a model in the format used by Excel.
#
# @param model the fitted model
# @param data the date set used for fitting
#
# @return the formatted ANOVA table
#
excelANOVA <- function(model, data) {
  # Get the dependent variable.
  dv <- (model |> formula() |> all.vars())[1]
  # Regress the d.v. against just a constant term.
  model2 <- lm(paste0(dv, " ~ 1"), data = data)
  # Run ANOVAs for both models.
  a <- anova(model)
  b <- anova(model2)
  # Get the number of predictors and the row count of the original ANOVA.
  r <- nrow(a) # total rows
  p <- r - 1   # number of predictors
  # Put the degrees of freedom, sum of squares and mean squares for the full model in row p.
  a[p, 1] <- p                   # df
  a[p, 2] <- sum(a[1:p, 2])      # ss
  a[p, 3] <- a[p, 2]/p           # ms
  # Drop the rows for individual predictors, leaving just model and residual rows.
  a <- a[p:r, ]
  # Fix the F statistic and p-value.
  a[1, 4] <- a[1, 3] / a[2, 3]                  # F statistic
  a[1, 5] <- 1 - pf(a[1, 4], a[1, 1], a[2, 1])  # p-value
  # Add the row for the constant model.
  a <- rbind(a, b)
  # Remove the MSE from the total (constant model), since Excel does not show that.
  a[3, 3] <- NA
  # Tweak the row and column labels to suit the boss.
  rownames(a) <- c("Regression", "Residual", "Total")
  colnames(a) <- c("df", "Sum of Squares", "Mean Square", "F", "Significance F")
  # Return the ANOVA table.
  return(a)
}

Here is the output it produces on the test instance.

	df	Sum of Squares	Mean Square	F	Significance F
Regression	2	193.489	96.744	334.35	6.6418e-12
Residual	13	3.762	0.289
Total	15	197.250

The default behavior of lm() is to encode the data in the lm object, so it might seem that the second argument to the function is redundant. I included it in case someone is working with a large dataset and suppresses the embedding of the data in the model object. One could tweak this to make the second argument optional, with the default being to grab the data from the first argument.

CRAN Repository Selection

RStudio gives you some control over the installation of packages via the menu path Tools > Global Options... > Packages. On the rare occasions I went there, I got a warning message at the top of the panel: "CRAN repositories modified outside package preferences". I'm not sure if this means that changes made in the settings dialog will not have any effect at all, or whether their effect will be limited to that RStudio session. Either way, until today I paid it no mind.

For reasons best captured by the saying "idle hands are the devil's workshop", I decided this morning to try to fix the problem. A Google search convinced me that other people are encountering this but was a bit shy of specifics on how to fix it. I did learn that running getOption("repos") in the RStudio console (or in an R session outside RStudio) would tell me where R would look for packages for installation or update. On my system (Linux Mint), it was defaulting to https://cloud.r-project.org. So the next step was to figure out where that was being set.

After considerable poking around, I found that code at the bottom of /etc/R/Rprofile.site was setting the option. After editing the file (as root) and commenting out that code, I restarted RStudio and found that the error message was gone. RStudio now points to  https://cran.rstudio.com/ (which may or may not be a mirror of the first site).

So what happens now if I run R from a terminal (outside RStudio) and try to install/update packages? With the original option no longer set, R pops up a list of mirrors and has me pick one. It's an extra step I could live without, but since I normally install packages from RStudio and not a plain R session, I'm content to leave the change in place. I suspect, though, that the next time I upgrade R itself I will have to remember to edit the global Rprofile.site file again. 

Identifying Extreme Points

I found a recent blog post by Erwin Kalvelagen titled "Convex hull models" rather interesting. Although the title mentions convex hulls, what Erwin actually discusses is finding the extreme points of the convex hull of a finite point set. (The full convex hull would also include bounding hyperplanes, which is a separate, arguably thornier, matter.) Erwin presents both a mixed integer nonlinear programming (MINLP) model and two mixed integer linear programming (MILP) models. In computational tests, they all solved at the root node, leading Erwin to make the following observation: "[a]s we consistently get integer models that need 0 nodes it is likely this model is really an LP." He tested the LP relaxation of one of his models, and it seemed to work.

My interest stems in part from the fact that after contemplating the problem a bit I assumed it could be solved as an LP. Assume we have a finite set $P=\lbrace p_1,\dots,p_m\rbrace\subset \mathbb{R}^n$ of points. Our mission is to determine which are the extreme points of $\mathrm{conv}(P),$ the convex hull of $P.$ A starting point is the recognition that the extreme points of $\mathrm{conv}(P)$ must belong to $P,$ and that an extreme point cannot be a convex combination of other points in $P.$

Now consider this system of linear equations and inequalities: 

$$p_i  = \sum^{m}_{k=1}\lambda_{k}p_k\quad (1)$$$$\sum^{m}_{k=1}\lambda_{k} = 1\quad (2)$$$$\lambda_{k} \ge0\ \ \forall k \quad (3).$$

This expresses point $p_i$ as a convex combination of the points in $P$ including itself. So a trivial solution is $\lambda_i = 1$ and $\lambda_j = 0$ for all $j\neq i.$ If a solution exists with $\lambda_i=0$ then $p_i$ is a convex combination of the other points and therefore not extreme. If there is no solution with $\lambda_i=0,$ then $p_i$ is not a convex combination of other points and therefore must be extreme. This logic, by the way, is also embedded in Erwin's models. 

Suppose that we find a solution with $0 < \lambda_i < 1.$ Combining the $p_i$ terms on the left side of (1), we have $$(1-\lambda_i) p_i = \sum_{k \neq i} \lambda_k p_k$$ and therefore $$p_i = \sum_{k \neq i} \frac{\lambda_k}{1-\lambda_i} p_k.$$ Since $$\sum_{k=1}^m \lambda_k = 1 \implies \sum_{k\neq i}  \frac{\lambda_k}{1-\lambda_i} = 1,$$ we know that if a solution to (1)-(3) exists with $\lambda_i < 1$ then a solution exists with $\lambda_i = 0.$

So we can solve the problem of minimizing $\lambda_i$ subject to the constraints above. If the minimum value is 0, $p_i$ is not an extreme point. If the minimum value is strictly greater than 0, then it must be 1, which means that $p_i$ cannot be written as a convex combination of the other points and therefore must be extreme.

Along the lines of what Erwin did, we can write a single LP that identifies all extreme points. We add a second index to $\lambda,$ so that $\lambda_{i,k}$ is the weight assigned to $p_k$ when writing $p_i$ as a convex combination of the points. We then stack the constraints for each $p_i$ into a single model and minimize $\sum_i \lambda_{i,i}.$ Positive values of the objective terms identify extreme points. That said, I find it perhaps more convenient to solve a sequence of smaller LPs, one for each point in $P.$

As a proof of concept, I wrote an R Markdown document that generates a random set of points in two dimensions and uses the single-point model to find extreme points. I stuck to two dimensions to allow for easy plotting of the results. It uses the lpSolve package, a wrapper for the open-source lp_solve LP/MIP solver. (It also needs the ggplot2 package for plotting purposes.) You can download the .Rmd file and try it for yourself (assuming that you have R and the two aforementioned packages installed). On my fairly humble PC, tow-dimensional test cases with 1,000 points solved in the blink of an eye. Of course, solution times increase as dimension and sample size increase. It is an empirical question (for which I have no answer) whether solving for all points at once in a single model is faster or slower than solving for each point individually. To belabor the obvious, doing all in a single model will require much more memory.

R Upgrade Glitch

Today I was able to download and install the next version (4.6) of the R programming language. Almost exactly a year ago I posted about retaining libraries during an R upgrade. So I followed my own advice (not always a good idea but appropriate in this case), copied the version 4.5 user packages over to 4.6, and then checked for any needing an update. The only one flagged in my case was the devtools package ... but when I tried to update it, the update failed.

The initial error message pointed to a missing symbol definition in the rlang package. There was no update available for rlang, so I uninstalled it, reinstalled it (no glitches) and then tried devtools again. This time, I got a missing symbol error relating to a different package (vctrs). Encouraged by the fact that it at least was not complaining about rlang anymore, I uninstalled and reinstalled vctrs, then ran the devtools update yet again. This time the update succeeded.

I'm not sure I took the most efficient path here, but at least I made it to the goal line.

 

Shiny: Highlighting DT Rows

I am working on a Shiny application involving data tables (displayed using the DT library). There is one particular table in which I would like my code to highlight rows meeting certain criteria, which the application will evaluate. Not knowing how to set a background color for an entire row in a table, I went down a rabbit hole with Shiny Assistant, the relatively knew AI bot trained to help with Shiny coding. It took way too many attempts, but eventually I stumbled on prompts that got me more or less what I need.

It turns out there are some tricky bits, so I put together a small demonstration application. The demo app displays a portion of the mtcars dataset (which installs with R) and lets the user select, via check boxes, which rows to highlight. You can download the source code here. I will just mention a few things about the code.

• The check boxes to select which rows are highlighted belong to a single checkboxGroupInput control named "selected". Initially I used observeEvent() to watch for changes to input$selected, which almost worked. The problem was that when the user removed the last check mark, signalling that no rows should be highlighted, the observer did not see an event and the table continued to have one highlighted row. Switching to observe()fixed the part about not seeing an event. I have no explanation as to why. See the next paragraph for the rest of the fix.
• Early on, the table itself would not display until there was at least one row selected for highlighting. After switching to observe(), removing the last checkmark still failed to clear all formatting. I was forced to use an if-else approach (the highlightRows() function), in which the code applies the row formatting only if at least one row is selected. I am still in the dark as to why this is necessary.
• The formatting code ends up being translated to JavaScript to communicate with the DataTables library. DataTables uses zero-based indexing, which explains why the specification of which columns to highlight (in this case, all of them) is 0:(ncol(df) - 1) rather than 1:ncol(df).
• Using reactiveVal() may be overkill in the demo application, since input$selected is itself reactive (I think) (maybe). I used it in part because I intend to use it in the larger application I am building, and I wanted to be sure it worked.

 

Finding a Row that Moved

This post deals with programmatically scrolling a table in a Shiny app (using the DT library) to a specific row. 

I've been working on a Shiny app involving a number of tables, which are displayed using the R DT package. The tables include controls to let the user sort them by any column, which plays a role in today's post. Controls external to the table let the user insert, delete and edit records. The reason I use external controls and not the editing capabilities built into DT is that there are a number of limitations on edits that the application code enforces.

Before getting into the weeds, I need to establish a bit of terminology. The underlying data in any table has associated row numbers, just as in a spreadsheet. Those do not change when the user sorts the displayed table. I'll refer to that as the "row index", and to the position (1 if first, 2 if second, etc.) of a row in the displayed table as the "display index". If the seventeenth row of a table winds up at the top after sorting, it has row index 17 and display index 1. (Sorry if I'm belaboring the obvious.)

Here is the problem I ran into. Suppose that a user selects a row and edits it. When the user commits the change, DT updates the display, and the edited row may wander quite far in either direction and in particular leave the portion of the table being displayed. So I want the application to automatically scroll the table to the new location of the edited record. I initially assumed this would be easy. It was not. In fact, 20+ attempts with the Shiny Assistant coding AI and 10+ attempts with Claude AI failed to solve the problem (although they did produce some insights, and a bit of code described below).

We can split the problem into two parts: finding the new location of the record, and then scrolling to it. It turns out that the latter is the easier part. At least it's the easier part if you know JavaScript, which I do not. One of the bots (I forget which) provided me with the following code to embed in an R function. The arguments to the R function are outputID (the name of the table in the UI) and row (the row number to which to scroll).

js_code <- sprintf(
  "setTimeout(function() {
     var table = $('#%s table').DataTable();
     // Scroll to the row.
     setTimeout(function() {
       var rowNode = table.row(%d).node();
       if (rowNode) {
         $(rowNode).addClass('highlight-row');
         rowNode.scrollIntoView({behavior: 'smooth', block: 'center'});
       }
     }, 300);
   }, 200);",
  outputId,
  row - 1
)
runjs(js_code)

I'm not positive, but I believe you need to load the shinyjs R library for this to work. (My app was already using it.) The values 200 and 300 in the code are apparently delays (in milliseconds) intended to allow time for DT to finish rendering the table before scrolling occurs.

On to what I originally thought would be the easier part: finding the target record. What the JS code is looking for in the second argument (row number) is actually what I'm calling the display index. My code knows the row index of the target record. Assume the output ID of the table is "A". After some digging, I discovered that input$A_rows_all contains the vector of row indices in display order (meaning the first entry is the row index of the record with display index 1 etc.). So if variable r contains the row index of the record we want, which(input$A_rows_all == r) will give us the display index of that row. Simple enough ... except that it did not work.

Poking around my code, I discovered that after the user committed a change the table display would update "immediately" but that input$A_rows_all retained its pre-change value. This turns out to be a synchronization issue -- DT updates the table and then updates the rows_all vector at it's own pace, while my function attempts to move immediately from the line committing the change to the line looking for the new display index.

The answer was to wrap the code doing the scrolling in an event observer: observeEvent(input$A_rows_all, { ... }) where "..." is the code that finds the target row and scrolls to it. This effectively pauses my code until DT has gotten around to updating the rows_all vector. 

Queueing Cuts in XpressMP

In my previous post, I described doing Benders decomposition with the XpressMP optimizer (Java API), including some sample code. As with previous code solving the same sample problem using CPLEX, a callback was required. The callback takes candidate solutions from the master problem, solves the corresponding LP subproblem, and either accepts the candidate or rejects it by adding a Benders feasibility cut or a Benders optimality cut. Well, that's how it works with CPLEX. The Xpress callback allowed me to add Benders cuts when the candidate solution was the optimal solution to a node LP in the master search tree, but not when the candidate was found via heuristics. (With CPLEX, the source of the solution did not matter.)

Being obstinate, I tried to work around that limitation. I updated my Java source code (still in the same repository) with a second Benders approach using a cut queue. In this approach, the PreIntsol callback solves the LP subproblem and, if warranted, cranks out a Benders cut. If the source was a node optimum, the PreIntsol callback adds the cut as before. If the source was a heuristic, the PreIntsol callback adds the cut to a queue within the Java code. A second callback, which implements the XpressProblem.CallbackAPI.CutRoundCallback interface, is called when Xpress is generating cuts. If any cuts are in the queue, this callback removes them from the queue and adds them via calls to XpressProblem.addManagedCut(). Thus, in theory, the opportunity to generate Benders cuts based on heuristic solutions does not go to waste.

So did this help? Somewhat, though not as much as I had expected. The first thing I discovered is that running either of the Benders models with a single thread was faster than running them with five threads. I suspect this is partly due to the test problem being very easy and partly due to multiple threads doing redundant work (for instance, churning out the same heuristic solution multiple times). 

 With a single thread, Benders with the queue was a wee bit faster than Benders without the queue. In one run, the original Benders model needed 597 ms. versus 578 ms. for the Benders run with the cut queue. (For comparison purposes, the basic MIP model with no decomposition needed 162 ms.) This might understate the advantage of the model with cut queue very slightly, since it does a bit more printing, and when you are dealing with differences in the millisecond range maybe time printing matters (?). The original Benders model had its callback invoked 28 times with integer-feasible node solutions, resulting in 25 optimality cuts and one feasibility cut. It was invoked 74 times with heuristic solutions. The version with the cut queue was invoked only six times with integer-feasible solutions but 74 times with heuristic solutions. It generated 56 optimality cuts and 21 feasibility cuts. The higher cut counts make sense, since it can generate cuts from heuristic solutions. What is a bit perplexing to me is that it encountered the same number of heuristic solutions, suggesting that the extra cuts blocked some node optima from occurring but had no effect on Xpress's heuristics.

With five threads, the story changes.  The run times were 168 ms. for the MIP model (no change after accounting for randomness), 817 ms. using the queue and 2557 ms. without the queue. So the queue cut Benders run time by about 2/3 or so. The original Benders callback was called 84 times with integer-feasible node solutions and 142 times with heuristic solutions, versus 17 times with integer-feasible node solutions and 87 times with heuristic solutions for the version with the queue. Note that, in this case, the number of heuristic solutions also went down.

So is the cut queue worth the effort of programming it (and the extra time devoted to invoke the cut callback repeatedly)? Maybe. It appears to depend at least in part on the number of threads being used, although one or two runs on one test case is far from conclusive. 

Benders Decomposition in XpressMP

Thirteen and a half years ago (yikes!) I posted about doing Benders decomposition in CPLEX, including source code (Java) for solving a small fixed-cost transportation problem using Benders. The source code changed a couple of years later and I moved it to a new repository. I probably should state up front that when I allude to Benders decomposition I mean what is known in some quarters as "one tree Benders" -- solving the problem in a single pass, using callbacks, rather than Jacques Benders's original approach, which solved the master problem to "optimality", added cuts and then started the master problem solution again from scratch.

I decided to try applying Benders to the same problem using the Java API to FICO's XpressMP solver. It turns out that implementing callbacks in XpressMP is a bit of an adventure, in part due to its design and in part due to rather sparse documentation for the Java API. (Someone at FICO advised me to consult the documentation for the C API whenever I could not find some detail in the Java documentation.) Eventually, I managed to get a working implementation of Benders for the same example used in those old posts. I will not repeat the formulation here. You can find it in the first of the posts above. It involves deciding which of a given set of warehouses to open and how much to ship from each opened warehouse to each client. Clients have demands that must be met, warehouses have capacities, there are fixed charges to open warehouses and per-unit flow costs between warehouses and customers. The warehouse-customer links do not have any capacity limits. In the Benders decomposition, the master problem (MIP) contains binary variables indicating which warehouses are opened plus a surrogate variable for the total flow cost. It initially has no constraints. The subproblem (LP) contains continuous flow variables for each combination of warehouse and customer, demand constraints for the customers and capacity constraints for the warehouses.

The Java code for my implementation can be found in this repository. The "Code" drop-down menu provides options for downloading the code. Running it requires that you have a recent version of XpressMP installed. Outside that (and having Java, of course) there are no software requirements. The code first runs a single MIP model as a benchmark, then runs the Benders model.

Callbacks

The Benders implementation requires only one callback in the master problem. It implements the XpressProblem.CallbackAPI.PreIntsolCallback interface. According to the C documentation, Xpress calls it "when an integer solution is found by heuristics or during the branch and bound search, but before it is accepted by the Optimizer". One of the arguments (soltype) tells you whether the solution is an integer-feasible node optimum (soltype = 0) or a heuristic solution (soltype = 1). The distinction is important because Xpress will allow you to add cuts when it has a node solution but will not allow you to add cuts when it has a heuristic solution.

The callback includes an argument (pReject in my code) that you set to 0 to accept the solution or 1 to reject it. Somewhat curiously, if a node solution is infeasible or underestimates flow costs and you are able to add a Benders cut, you are supposed to "accept" the solution (set pReject to 0) in the callback (else the cuts are lost, if I understand this correctly). If a heuristic solution is infeasible or underestimates the flow cost, you reject it (set pReject to 1). The argument has class IntHolder, so you do not directly set it to 0 or 1; you instead set its value field to the desired value (for instance, pReject.value = 0 to accept the solution).

According to a source at FICO, there is currently no callback that allows you to add Benders cuts when dealing with a heuristic solution. All you can do is reject it if it is infeasible or underestimates the flow costs. This results in the master problem solver repeatedly generating (via heuristic) the same or similar infeasible solutions, which could be avoided if Benders cuts could be added.

While researching this, I came across a Benders example (in Python)  that also used a callback which in Java would implement the OptNodeCallback interface. This is apparently no longer necessary, due to changes to the solver.

Variable names

As with all other solvers, Xpress lets you assign names to variables.  In my code, I print out the Benders cuts being added. This is mainly for debugging purposes. I certainly would not do it in a production problem with hundreds or thousands of binary variables. It turns out this is not as easy as one might think. For rather arcane reasons that I will not get into (since I would get the explanation wrong), inside the callback the variables are named X[0], X[1], X[2] etc., regardless of whatever names you originally assigned them. This is a bit confusing (and makes debugging tricky, since it's unclear which of the original variables is now known as X[12]). Daniel Junglas at FICO provided me Java code for a workaround, which is included in my code. You will find it in the BendersModel.printCut method.

Thread safety

When using multiple threads, there is good news and bad news. The good news is that Xpress goes to great lengths to be thread-safe.  Part of their implementation of thread safety is the first argument of the callback function, which is a presolved (I believe) copy of the original (master) problem. This means that if the callback is running simultaneously in multiple threads, each thread is working with a different copy of the problem, avoiding collisions. 

The bad news is that this introduces some rather curious problems. In my case, problems arose when printing the cuts being added. The cuts were expressed in terms of the original variables (defined in the master problem and saved in an array for later use generating cuts). For whatever reason, the cuts themselves came out fine, and adding them worked, but when I went to print them (which implicitly meant invoking toString() on the cut expressions) Xpress threw exceptions and eventually crashed. This only happened when using multiple threads. If I throttled Xpress to a single thread, printing the cuts caused no problems. Fortunately, the workaround of the variable naming issue (implemented in the printCut method) also eliminated the exceptions occurring when I printed the cuts.

So, bottom line, I now have working code that solves a problem via Benders decomposition with the Java API to Xpress. You are free to download it if you wish, provided you promise not to laugh at it. Okay, not to laugh at it too much. 

XpressMP Java Examples

I've been looking at the examples that ship with the FICO XpressMP optimizer, specifically those using the Java API, and I ran into a bit of a conundrum. The examples live in the .../examples/solver/optimizer/java folder, which also contains a subfolder named objects. The file .../examples/solver/optimizer/index.html has a section title "Calling the library from Java" which lists the various example problems, gives a short description of what they exemplify (for instance, "Using Xpress callbacks"), and points to the corresponding source code files.

The first point of confusion for me was why the objects subfolder contained examples that in some cases seemed to duplicate what was in the parent folder.  The second point of confusion was why some example entries pointed to multiple Java files. For instance, the entry "Goal programming: Lexicographic goal programming using the Xpress multi-objective API"  points to a single file GoalProg.java, whereas the entry "The travelling salesman problem: Using Xpress callbacks" points to two files, TSP.java and TravelingSalesPerson.java. The latter is in the objects folder. So why two versions?

A kind soul at FICO cleared things up for me. The Java API is being transitioned from a "thin wrapper" around the C API (based on the XPRSprob class) to a more object oriented version (based on the XpressProblem class). The old API is not going away, but the new API is now the recommended one. Examples in the objects folder use the new API, while those in the parent folder use the old API.

One distinction between the two APIs becomes apparent when you look at the TSP examples, which employ callbacks. The distinction has to do with the standard practice of taking the user's original model and modifying it in the presolve stage. With the original API, you have access to the presolved model, and my understanding is that you can add a cut either to the original model or to the presolved model. Since those models differ substantially (including in regard to which variables they contain), if you add a cut to the presolved model you have to "crush" it (convert from original space to presolved space) yourself. In the newer API, I believe you only have access to the original model, and so you express cuts in terms of the original variables and Xpress handles the crushing ... which is fine with me. Barring some weird exigency, I'll be sticking to the new API.

 

A Shinyapps Misfire

Posit, the company behind the open source RStudio IDE that I use for R coding and many other R-related products and services, operates an online server for Shiny web applications at shinyapps.io. As with their other services (and as is common with many web services these days), they offer various pricing tiers, including a free tier that I use. The free tier, as one would expect, comes with a variety of limitations, including lack of direct customer support. (Free tier users are encouraged to use their community forums for help.) I use my shinyapps account primarily for developing applications for nonprofit clients of INFORMS Pro Bono Analytics (PBA). Those clients host the completed applications on their own accounts (whether free or paid).

Recently, my account was suspended after I did a test run of a program I was writing for a  PBA client. It took some digging to find out why, but once I did the support staff at Posit very quickly restored my account, for which I am grateful. It turns out that, in addition to limitations on number of applications hosted and usage amounts (time, memory, ...), there are certain things an application in a free account cannot do that an application in a paid account can do. As best I can tell, the following things can get your account suspended and/or cause your application to crash:

•  launching a headless browser (such as Google Chrome or Chromium); or
•  trying to launch background processes (parallel workers, system scripts, ...).

I'm not sure that launching background scripts/processes will result in account suspension, but I have it on good authority that they can cause applications to bomb or act funny (for instance, because the server silently blocked the background process). Also, the above list is by no means guaranteed to be comprehensive. There may be other ways to get your account suspended.

In my case, the Shiny application I was testing used the leaflet library to draw maps and the mapview library (which in turn uses the webshot2 library) to create PNG images of those maps. I did not realize that one of those (I think webshot2 but don't count on my being correct) uses a headless browser (Chromium?) to do the HTML to PNG conversion. OOPS!

Apparently there are legitimate reasons for Posit to be concerned about headless browsers -- ways they can be abused (?). I just wish that the response to their use in the free tier were an aborted application launch with explanation and not account suspension. Possibly my account was suspended because I tried multiple times to get the application running ... because I had no idea it was doing something not permitted.

Anyway, I'm back up and running, and now I know better than to try that approach to HTML -> PNG conversion. 

 

Timing is Everything?

I've been running the latest versions of Linux Mint (MATE desktop) on a somewhat long in the tooth HP desktop for years now. Part of my setup used the "Startup Applications" feature of the system menu to load the following applications:

• Firefox (web browser) with a half-dozen tabs preloaded;
• Thunderbird (email client) with connections to two email accounts (using IMAP);
• LibreOffice Calc (spreadsheet) with a particular file opened;
• Workrave (a timer that signals me when I need to take a break from banging on the keyboard); and recently
• Diodon (clipboard manager).

These all ran like a charm ... until I replaced the HP tower with a much faster and more powerful "mini" computer.

After setting up the new computer, all the apps loaded at boot as before (only much faster than with the HP). Firefox and Calc continued to run just fine, as did other applications that I launched manually as needed. That leaves three of the original five programs, which had issues.

• Thunderbird would intermittently (but frequently) fail to connect to the mail servers. It would grab the contents of mail folders when it started, but subsequent updates might or might not work, and sometimes I could not send a message (or even save it as a draft). Shutting down and restarting the application seemed to consistently fix the problem.
• Workrave mostly worked, but sometimes the break timer would count down to zero and then freeze rather than resetting, and sometimes I would walk away from the computer long enough to qualify as a break but Workrave would continue counting down as if I had been using the machine. Again, exiting and restarting the application fixed the problem.
• Diodon would initially fail to recognize my copying text. The copied text would be in the clipboard (i.e., I could paste it) but not in Diodon's popup clipboard history. Once again, restarting the application solved the problem.

Each application had a different problem, but the common thread was that the problems only occurred on the initial (automatic) load. I suspect (although I cannot be sure) that there is some timing issue during startup affecting things -- perhaps because the new computer is faster, perhaps because it has more cores and threads than its predecessor, perhaps due to something beyond my limited grasp of computer engineering. So I tried playing with the feature in  "Startup Applications" that lets you program in a delay in starting an application, in the hope that it was indeed a timing issue and that I could disrupt the issue by delaying starts. Here are my results.

• Thunderbird seems to be working properly with a 50 second delay. I tried 10 seconds initially, which worked sometimes but not always.
• Workrave works correctly for the most part with a 12 second delay. It has glitched once or twice since I made the change, but it is much better with the delay than without.
• Diodon did not work correctly with 8 or 15 second delays, but works fine with a 60 second delay.

The moral of the story here is that the startup delay setting apparently has a good reason for existing (and that weird things can happen due to "random" timing issues).

 

Building a Partial Tour

A user on Mathematics Stack Exchange posted a question that, while superficially similar to the prize collecting traveling salesman problem (TSP), also has substantial differences to it. You are given a rectangular grid where each cell contains a reward. The objective is to build a tour of some (but not all) cells, with movement restricted to adjacent cells (those sharing an edge with the current cell), subject to the restriction that you cannot visit two cells with identical rewards. The objective is to maximize the value of the cells visited. It differs from the prize collecting TSP in several ways:

• there is no designated start/end cell for the tour;
• there are no edge costs (movement incurs no objective penalty);
• there is no penalty for unvisited cells; and of course
• duplicate rewards are prohibited.

A key part of the user's question was how to handle subtour elimination. The user was looking at the classic Miller-Tucker-Zemlin (MTZ) approach for TSPs, which involves an auxiliary variable at each cell recording in essence the position of the cell (first, second, ...) in the tour. The implementation of MTZ involves constraints at each cell except the depot (tour origin/destination). In the current problem, there is not designated depot, and you cannot just pick a cell to serve as the depot because you do not know which cells will be in the optimal tour.

The workaround I suggested in a reply to the post is to add an artificial cell (the "depot"), connected to every cell in the original grid, and use that to implement the MTZ constraints. The one mildly tricky part is that the cell immediately preceding the depot and the cell immediately following the depot in the tour need to be adjacent (or the same, if the tour only visits one cell). That is easily handled at the cost of some extra constraints.

Assume the grid contains $m$ rows and $n$ columns. Let us start by numbering the cells $1,\dots,mn$ in raster scan order from upper left to lower right, and call the depot cell 0. Let $r_i$ denote the reward for visiting cell $i$ and $N_i$ denote the set of cells (including the depot) adjacent to cell $i$. Finally, let $R$ denote the set of possible rewards and $C_r$ the "cluster" of cells having the same reward $r\in R.$ We build the model as follows.

Variables

• $x_{ij}\in \{0,1\}$ will be 1 if the tour moves from cell $i$ to cell $j,$ 0 otherwise. 

• $y_{i} \in \{0,1\}$ will be 1 if the tour visits cell $i > 0,$ 0 if not.

• $z_{i} \ge 0$ is the auxiliary variable (position) for cell $i$ in the MTZ constraints (indexing where cell $i$ falls in the tour). 

If cell $i$ does not appear in the tour, the value of $z_{i}$ will be irrelevant (and the solver will likely set it to 0).

 

Objective

The objective is to maximize $$\sum_{i} r_i y_i.$$ 

 

Constraints

• Movements may only occur between adjacent cells (with the depot cell adjacent to every other cell).  $$x_{ij} = 0 \quad \forall i, \forall j\notin N_i.$$

 

• The depot must be exited and entered exactly once. $$\sum_{i > 0} x_{0i} = 1 = \sum_{i > 0} x_{i0}.$$

 

• Each cell other than the depot is entered and exited once if visited and zero times if not visited. $$\sum_{j} x_{ji} = y_i = \sum_{j} x_{ij} \quad \forall i>0.$$

 

• At most one cell with a given reward value can be visited. $$\sum_{i\in C_r} y_i \le 1 \quad \forall r\in R.$$

 

• MTZ: If the tour moves from cell $i$ to cell $j,$ the position value for cell $j$ must be at least one greater than the position value for cell $i$ (excluding $j=0,$ i.e., return to depot). $$x_{ij} = 1 \implies z_j \ge z_i + 1 \quad \forall i, \forall j > 0.$$

 

• If the tour goes $0 \rightarrow i \rightarrow \dots \rightarrow j \rightarrow 0,$ then cells $i$ and $j$ must be adjacent (or the same). $$x_{0i} + x_{j0} \le 1 \quad \forall i>0, j>0, i\neq j, j \notin N_i.$$

 

I tested the code on the example in the Math SE post, using Java and the FICO XpressMP MIP solver, and it produced what appears to be an optimal tour (confirmed by the original poster) in about two minutes on my somewhat archaic desktop. I posted the tour in a comment to my answer to the original question. The Java code is available (Creative Commons 4 open source license) from my Git repository.

Mint MATE Monitor Brightness

I ran into a rare but not unheard of problem this morning while working on my desktop computer (running Linux Mint with the MATE desktop, and using a Dell monitor). Due to a somewhat funky arrangement of solar position, clouds, and possibly other celestial entities, the natural light behind my monitor caused a glare that made my eyes hurt a bit and caused me to have trouble seeing the screen. (Some what I guess you might call after-images were messing with my vision.) To reduce the glare, I wanted to dial down the monitor brightness ... but could not figure out how. The buttons on the Dell monitor are unlabeled, and I was not overly excited about a bunch of trial-and-error button pushing. So I naturally went to the system's display controls ... which did not have a brightness setting.

When I had some free time, I started searching the web for how to adjust display brightness in MATE, where to find GUI tools or applets for display brightness, and so on. I found lots of seeming helpful posts that either (a) referred to controls that do not exist in my operating system, despite allegedly being specific to MATE and possibly to MATE on Linux Mint or (b) pointed me to an applet that did not work (in some cases possibly because it only works on laptops for some reason?).

So, bottom line, I eventually wrote a bash script with a primitive GUI to let me futz with the display brightness. It uses the xrandr  and zenity commands, either or both of which may need to be installed. (You can run "which xrandr" and "which zenity" in a terminal to see if you have them.) Once you have them, run "xrandr" in a terminal with no options to find out the name of your display (DP-2 in my case). Substitute that for DP-2 in the following script (and remember to make it executable with "chmod +x <your script>.sh", and you should be in business.

#!/bin/bash

# Adjust display brightness (1.0 being the default brightness).
# You can get the monitor name (DP-2 here) by running 'xrandr' in a terminal.

target=$(zenity --entry \
                --title="Set Display Brightness" \
                --text="Set a brightness level (base value = 1.0):" \
                --entry-text="1.0")
xrandr --output DP-2 --brightness $target
zenity --info --text="Click OK to close the terminal."

A Shinyvalidate Hack

Lately I've been coding a Shiny web application (in R) in which users are confronted by several forms. Since their inputs go into a database, some amount of input validation seems prudent (to put it mildly). I came across the shinyvalidate package, which is proving very useful. You can specify a variety of rules for input fields (starting with whether they are required and moving on to various filters on what is allowed). If the input fails to meet the validation rules, messages are printed in red under the input controls (something you've probably seen before if you ever tried to submit a web form with something missing or glaringly incorrect).

The remainder of this post probably won't make much sense unless you have some experience with shinyvalidate. (I'm not promising it will make sense if you do.) I ran into a bit of a snag with a few input fields where the validation rules involved a disjunction. In one case the fields was the URL of an organization's web site, and in the other it was their contact email address. Both fields are "required" in my application, with required in quotes because it is possible for an organization not to have a web site and/or not to have a contact email address. (For instance, some organizations have web sites but no email address, expecting people to use a contact form on their home page.) The approach I'm taking is to allow "none" as a valid response for both web and email addresses if there is none. So I want to enforce the rule that the input value is either "none" or a valid address.

The shinyvalidate package provides mechanisms for enforcing conjunctions of rules (satisfy all of these requirements) but apparently does not have an explicit mechanism for enforcing disjunctions (either this or that). It does let you create custom functions to validate inputs, but I wanted to use their built-in web (sv_url()) and email (sv_email()) functions so that I would not have to go down the regular expression rabbit hole. That meant writing a function that combined a simple test (x != "none", where x is the value submitted in a text field of the form) with their functions.

It took me a while to figure out that while sv_url and sv_email take a few arguments (including an optional customized error message), they do not actually take the form input (what I'm calling x) as an argument. Instead,  the output value of sv_url() or sv_email() is itself a function with the form input as its sole argument. Once that sank in, making the custom validators was trivial. Here is my URL validator.

urlOrNone <- function(x) {
  if (x != "none") sv_url("Please enter a valid URL, including the http/https prefix.")(x)
} 

If x either equals "none" or satisfies the sv_url validator, the return value is NULL, which allows the input to pass validation. Otherwise, the user gets the custom error message and the input is disallowed. 

Android Silliness Part II

Last year I documented at great length a mess I went through when a very straightforward morning radio alarm I had set on my bedroom smart speaker had to be replaced by a Google "automation". That post ended with my having successfully implemented a solution.

Well, that was then and this is now. Google recently introduced Android 16 and I recently upgraded my phone from a Pixel 6a to a Pixel 9a. One or both of those events killed my morning alarm. Instead of hearing Google's voice telling me "streaming NPR on WKAR" (my local NPR radio station), I heard "OK, here's a playlist from YouTube Music" followed by some noises that I suspect are used to soften up presumed terrorist before interrogating them.

This was baffling in several ways, starting with the fact that I have never, ever instructed Google to play random music from YouTube Music early in the morning. More baffling is that when I went into the Google automation settings, my morning routine was still there, still marked enabled, and still listed as doing the same things: playing the radio (station = WKAR FM) on the bedroom speaker. Manually triggering it, however, did not work.

I have no idea what the problem was (is), but I eventually found a work around. I changed the "action" setting from playing the radio (one of a menu of possibilities controlled by a multiple select input) to a "custom" action (selected via "Try adding your own (Experiment with custom actions)". That opens a text input where you can type pretty much anything you might verbally ask Google Assistant to do. In my case that meant "Stream WKAR FM radio". To date, that has worked, which I assume means Google has engineers working on a fix for it.

If any Google engineers are reading this, may I recommend consulting the "Redneck Repair Manual" (page 1): "If it ain't broke, don't fix it."

 

RStudio Republish Menu

I'm a big fan of the RStudio IDE for coding in R, but like all software it has one or two quirks. Someone recently asked about one on the Posit Community forum. It's one that has annoyed me a bit.

When editing certain kinds of files (Shiny applications, RMarkdown or Quarto documents, ...), the IDE gives you a drop-down menu to publish (or republish) the document. The menu lets you publish to an existing account or to define a new account and publish there. It lists accounts to which you previously published this app or document. There is also an option to clear the list of previous locations. Here comes the gripe: there is no option to delete selectively one or more previous locations from the list. All you can do is clear the entire list.

It turns out that you can partially but not completely winnow the list, but it takes a file browser or terminal/command prompt and a bit of exertion. Since I've only published Shiny apps, I'll document the steps for that. The process for a Quarto or RMarkdown document is presumably similar but perhaps not identical.

The first step is to navigate to the folder containing the document or application. It should contain a folder named "rsconnect", which is where you want to go next. In there, you should find a folder for each site (server) on which the app has been published. In each of those folders you should find a folder for each name/account under which the app was published on that server. Delete the folders corresponding to the locations you want to delete from the "republish" menu. If you are deleting locations on more than one server, repeat for each pertinent server folder.

For example, I wrote a Shiny application for a colleague to use in a course. The app is installed on shinyapps.io in two places, a paid account used by my colleague and my free developer account. The app files (ui.R and server.R -- it predates the option to use a single consolidated app.R file) live in file X on my PC. So I go there and drill down to X/rsconnect. Since both installations are on the same service, there is only one file there, X/rsconnect/shinyapps.io. Inside that folder are two folders, one bearing the name of the course account (call that A) and the other bearing the name of my developer account (call that B). If I want to remove just my developer account from the republish list, I delete folder X/rsconnect/shinyapps.io/B but leave A in place.

Note that this does not remove the app from any server on which it is current published. For that you need to log into the server and do something. On shinyapps.io, you find the app in your administrative panel, sleep it, archive it and then delete it.

 

Routing With Sequencing

The motivation for this post comes from a sequence of questions posted on OR Stack Exchange (including this one), having to do with a mixed integer programming model for routing an electronic vehicle (EV) serving various customers. One difference from the basic single vehicle routing models with which I'm familiar is that the EV has to visit a charging station periodically during the route. That is easy to accommodate. Where it gets funky is that the modeler needs to know within the model which customer was last on the route, because the EV is required to go to the nearest charging station after its last stop. I'll take that a step further and require the model to provide (via variables) the position of each node (first, second, third) in the route sequence. This might be useful if, for example, the model had to enforce a rule that customer X must be among the first three customers served. Identifying just the last customer node is easier, as I'll describe at the end.

Attempts by the author of the original question followed the usual pattern for vehicle routing. Assume that there is a single vehicle, each customer must be visited once, and there are no time windows complicating things. You have a digraph containing nodes for each customer and each charging station. You typically start by assigning a binary variable $x_{ij}$ to each arc $(i, j),$ taking value 1 if and only if the vehicle crosses that arc, and proceed from there.

To collect sequencing information, I would normally employ the Miller-Tucker-Zemlin formulation of subtour elimination constraints. The MTZ approach adds a nonnegative auxiliary variable $u_i$ for each node $i$ together with the constraints $$u_j \ge  u_i + x_{ij} - M(1 - x_{ij})$$ for each pair of distinct nodes $i\neq j.$ This says that if we cross arc $(i, j),$ the value of $u_j$ must be at least one higher than the value of $u_i,$ preventing any loops. If $n$ is the number of nodes and we are willing to start numbering with $u_s=0$ ($s$ being the starting node for the tour), we can choose $M=n-1.$ The MTZ constraints are intended to prevent subtours, but as a side effect the $u$ variables number the stops in the order they occur.

This would work for the EV problem if there were a rule that the EV cannot use the same charging station twice during a tour. If the vehicle can stop more than once at the same charging station (which I assume would normally be the case), we cannot use the MTZ constraints because a repeat visit to a charging station would create a subtour.  This also complicates (I think) the use of subtour elimination constraints to prevent disjoint subtours. Fortunately, there are at least two "reasonable" (in my opinion) workarounds, both using the MTZ constraints. Unfortunately both are clunky.

The first workaround is to create multiple clones of each charging station node. So if node $s$ represents a charging station, we introduce additional nodes $s', s'', s''' \dots$ that are all charging stations, all in the same location (meaning time / distance / charge consumption between node $i$ and any of the clones is the same). For any arcs $(i, s)$ and $(s, j)$ we add arcs $(i, s'), (s', j), (i, s''), (s'', j)$ etc. We do not require that every charging node be entered (unlike customer nodes, which must all be visited), but we do limit each charging node to at most one entry. That removes the threat of loops and let us use the MTZ approach. Besides making the digraph larger, this forces the modeler to guess how many clones of each charging node will be needed.

The other approach is change to a multigraph, with fewer nodes but more arcs. We include only customer nodes, plus dummy start and end nodes. Arcs from the start node to each customer (within EV range) are the same as before. The end node is a stand-in for the closest charging station to the last customer visited. The arc from any customer node to the end node has the time / distance / charge consumption required to reach the closest charging station. (The closest charging station to each customer node is computed before building the model.)

For each pair of customer nodes $i \neq j,$ the arc $(i, j)$ (if it exists) represents moving directly from $i$ to $j.$ For select charging nodes $s$ we add another arc from $i$ to $j$, which I will denote $<i, s, j>,$ that represents going from $i$ to $s,$ recharging, and then proceeding from $s$ to $j.$ Each of those arcs produces another MTZ constraint. As usual, every customer node should be entered/exited exactly once.

If the number of charging stations is small, we can create extra arcs  $<i, s, j>$ for every combination of two customers and a charging station, weeding out those that are infeasible (meaning an EV with a full charge could not get from $i$ to $s$ or from $s$ to $j$). To get a smaller model, we can throw out arcs that are dominated by other arcs. For $<i, s, j>$ to dominate $<i, s', j>,$ you would need power consumption from $i$ to $s$ to be no greater than power consumption from $i$ to $s'$ and power consumption from $s$ to $j$ to be no greater than power consumption from $s'$ to $j.$ (If other criteria, such as mileage or transit time, appear in the objective function then they would also factor into the determination of dominance.)

One advantage of this approach is that it would let us dispense with the $u$ variables and the MTZ constraints if the only reason for them was to enable constraints forcing the EV to end at the charging station closest to the last customer, since that is baked into the arcs leading to the dummy end node. If we want to use the cloned charging station approach, we can also use a dummy end node linked to each customer by an arc representing the link to that customer's nearest charging station to enforce the desired ending rule for the tour.

Retaining Libraries During R Upgrades

Today I was able to upgrade R to version 4.5, and was reminded in the process of a tedious "feature" of the upgrade.

R libraries are organized into two distinct groups. If you use RStudio, look at the "Packages" tab. You will see the heading "User Library" followed by whatever packages you have installed manually. Scroll down and you will get to another heading, "System Library", followed by another group of packages. These are the packages that were automatically installed when you installed R itself.

The upgrade from R 4.4 to R 4.5 was very easy, since it comes as a system package, at least on Ubuntu and Linux Mint (and presumably other Linux distributions). I'm not sure about Windows, macOS etc. The Mint update manager offered me updates to several system packages (r-base, r-base-core, r-base-html and r-recommended) among the morning's gaggle of updates. I just installed those and R 4.4.3 was replaced by R 4.5.0. That part could not be easier.

After the updates were done, I opened RStudio and looked to see if any packages there needed updates. The "System Library" group was there, and none of them needed updates (no shock since they had just been installed during the upgrade). The "User Library" did not exist. I should have known this was coming based on previous R upgrades, but I forgot.
You can of course reinstall all your previously installed libraries manually (if you can remember which ones they were), or you can just wait until something doesn't work due to a missing library and install it then. I prefer to reinstall them all at once, and I most definitely do not have the list memorized. The fix is easy if you know how to do it (and remember that you have to do it). 

The first step is to open a file manager and navigate to the system directory where your libraries for the previous R version are stored. They will still be there. If you do not know where they are hiding, you can run the command .libPaths() and get a list of the directories in which R will look for libraries. One of them will contain the R version number. (It is consistently the first entry in the list when I do this, but I do not know if that will always be true.) In my case, the entry is "/home/paul/R/x86_64-pc-linux-gnu-library/4.5", which means I want to open "/home/paul/R/x86_64-pc-linux-gnu-library" in the file manager. There I find two directories, one for the previous version ( "/home/paul/R/x86_64-pc-linux-gnu-library/4.4") with lots of subdirectories and one for the new version ( "/home/paul/R/x86_64-pc-linux-gnu-library/4.5") that is empty. All it takes is copying or moving the contents of the older directory to the newer one. Once you have confirmed that the new R version can see the libraries (for instance, by observing that the "User Library" section has returned to the "Packages" tab in RStudio), you can delete the folder for the older version.

With that done, you will want to check for updates to the "User Library" packages. Several that I had installed needed updates today after moving to R 4.5. Updating them is done in the usual way.

I wonder if either R or RStudio has a "inherit libraries from previous version" function stashed away that would automate this? If so, I haven't found it.

Boolean Grid II

As the title implies, this is a sequel to my previous post, to deal with a variety of odds and ends.

There's a saying that you cannot teach an old dog new tricks. That's untrue; it's just that dogs as old as I am are slow learners. When I first started tangling with integer programs, there was a rule that you never made the model dimensions (number of rows or columns) any larger than you absolutely had to, partly to conserve memory and partly so as not to slow down pivoting. Once solvers switched from Gaussian pivoting (the way I learned to do it by hand) to matrix factoring, keeping dimensions down took a back seat to reducing matrix density. 

Similarly, once upon a time I learned (the hard way) that symmetry in my model would slow down pruning of the search tree. In the previous post, I alluded to research on exploiting symmetry and said something about solvers having symmetry detection. Imre Polik mentioned in a comment that Xpress already detects the symmetry by default. CPLEX might need some encouragement to do so. Near the start of the Xpress output, it describes the model as symmetric and lists statistics on "orbits" (groups of variables whose values can be permuted). It is unclear to me whether Xpress exploits that information by using "orbital branching" (a relatively recent development) or in some other way. Early in the CPLEX output I see a message that it is "detecting symmetries", but the model dimensions do not change and there are no further mentions of symmetry.

Moving on, Imre suggested in his comment to the previous post that another possible antisymmetry constraint is to assert a dominance inequality between the number of  true values in the top row and the number in the left column. This is compatible with my original two constraints, and so I added a version of it (combining Imre's constraint with mine) to my code. Rob Pratt suggested yet another possibility, an inequality between the subdiagonal and superdiagonal. I'm not convinced that one plays nice with my original constraints, meaning that if you added Rob's constraint to my original two you might legislate the optimal solution out of existence, so I added it to my code by itself. Also, it only works when the grid is square. Finally, since both solvers have parameters to control how hard they work to detect symmetry, I added an option to my code to skip adding any constraints and just crank up the solver's response.

The results are summarized in the following graph, which shows the optimality gap (best bound at the left end, best incumbent at the right end) for each solver and modeling option. The dashed vertical line is the optimal solution (from Erwin's post).

 

The "Bilateral" and "Trilateral" models are my original two constraints and my two constraints plus Imre's, respectively. The "Diagonal" model uses Rob's constraint. "None" is just the model by itself and "Solver" is the unmodified model plus a solver parameter setting to get it to work harder detecting symmetries. Note that the results for Xpress with "None" and Xpress with "Solver" are pretty much identical, confirming Imre's assertion that Xpress would detect the symmetry on its own. CPLEX saw a bit of improvement in both incumbent and bound going from "None" to "Solver", so apparently the nudge helped there. Only two runs found the optimal solution, both times CPLEX with some help from antisymmetry constraints. None of the runs got the best bound anywhere near tight.

Returning to my "old dog, new tricks" theme, the takeaway for me is that before I go nuts try to constrain away symmetry in a model, I need to investigate whether the solver can recognize it and, if so, whether it can eliminate or even exploit the symmetry.

Lastly, I belatedly realized that Erwin got his proven optimum quickly because the modified the model to use equality constraints in the interior of the grid and inequalities only in the two outermost rows/columns on each edge. I added that to my code as well, and yes, it gets a proven optimum incredibly fast. I was a bit leery about assuming that redundant coverage would only be required near the boundary, but per some comments by Rob on Erwin's post, 227 is indeed the (known) optimal value for a 32x32 grid.

A Boolean Grid

A recent blog post by Erwin Kalvelagen discusses a very straightforward integer programming problem. You have a rectangular grid of boolean variables, where a variables neighbors are the variable immediately above, below, to the left or to the right of it. The sole constraint is that, for any cell, at least one of that cell's variable or its neighbor variables must be true. The objective is to minimize the number of true cells in the grid. Erwin coded the model in GAMS and ran it for a 32x32 grid. He reported that he got an incumbent value of 227 in about 65 seconds but had trouble getting to optimality. (This might be a good time to point out that Erwin's computer is probably better than mine, since he is a consultant.)

I was curious whether a couple of redundant constraints would help. The problem suffers (if that is the correct term) from symmetry. Draw a grid and color in the cells of an optimal solution. Now flip the grid, switching either top with bottom, left with right, or both. The colored cells still form an optimal solution. What is the harm of symmetry? Think about the branch-and-bound (or, if you prefer, branch-and-cut) algorithm and specifically how it prunes nodes based on bound. When you find a new incumbent solution, you prune any node whose bound is no better than that solution. Typically the node bound will be at least slightly loose (meaning, in a minimization problem, that the bound will be strictly less than the objective value of the best feasible solution lurking in that node). In the context of the current problem, when a feasible solution is found, there will be at least three other feasible solutions with the same objective value, obtained by reversing the indexing of rows, columns or both. Each of them will likely be in a node of the search tree with an objective value at least slightly better than their true value, meaning that none of those nodes can be pruned right away, even if they do not contain an even better solution.

So symmetry can slow pruning and also slow improvement of the best bound. There has been research on how to exploit symmetry in IP models, but as far as I know that work has to be baked into a solver to be used. At least some solvers have some built-in capability to recognize and deal with symmetry, but I'm not sure how well that works. I usually keep an eye out for symmetry and, if I think it might be slowing improvement of the bound, see if I can constrain away some of it.

In this case, the symmetry I identified can be removed by adding just two constraints. One is that the number of true cells in the top row should not exceed the number in the bottom row (so that the vertical flip is ruled out unless the top and bottom rows are tied). Similarly, the other is that the number of true cells in the left column should not exceed the number in the right column (ruling out the horizontal flip). Since these constraints shrink the feasible region, it would not be surprising if they slow down identification of improved incumbents. By enlarging the model, they also slow down (very slightly?) the rate at which nodes are solved. The hope is that faster bound improvements compensate for that.

I ran the 32x32 case (coded in Java) using two solvers, FICO Xpress MP (version 44.01.01) and IBM CPLEX (version 22.1.2), both with and without the antisymmetry constraints. Each run was limited to 15 minutes (wall clock time) and used default settings for all parameters. Here are the results, in the format "best solution/best bound".

 

	Xpress MP	CPLEX
With antisymmetry	231 / 215.73	227 / 214.40
Without antisymmetry	228 / 215.95	229 / 214.38

 

There is enough randomness in IP solvers that I would not read much into a single iteration of each model. CPLEX actually did a bit better on the incumbent with the antisymmetry constraints included, which surprises me. Xpress had a very slightly worse bound with them included, which also surprises me. The bottom line seems to be that the antisymmetry constraints do not help much (which I find disappointing) and that, as Erwin noted, the problem is a bit stubborn.

As always, my code is available for download.