Finding largest empty shapes

Data

I generated $n=100$ data points $({\color{darkblue}x}_i,{\color{darkblue}y}_i), i=1,\dots,n$ drawn from a uniform distribution \[\begin{align}&{\color{darkblue}x}_i \sim U(0,10) \\ & {\color{darkblue}y}_i \sim U(0,10)\end{align} \]

It is always a good idea to have a careful look at the data. I.e., print and plot it (for very large datasets, print or plot a subset), and collect some elementary statistics (min, max, average, etc.). In practical modeling, it is surprising how many problems arise from simple data errors.

In the rest of this post, the dataset is formed by \[\begin{align}& {\color{darkblue}p}_{i,c} \in [0,{\color{darkblue}{\mathit{maxSize}}}] \\ & i \in \{1,\dots,{\color{darkblue}n}\} \\ & c \in \{x,y\} \\&{\color{darkblue}{\mathit{maxSize}}}=10 \end{align}\]

Data

----     26 PARAMETER p  points

                   x           y

point1         4.973       9.639
point2         8.739       8.443
point3         1.912       4.945
point4         2.346       3.810
point5         6.009       1.515
point6         5.077       4.633
point7         3.591       1.526
point8         3.664       8.260
point9         0.239       0.593
point10        2.848       1.140
point11        8.905       3.743
point12        2.341       1.477
point13        0.191       0.784
point14        9.280       5.550
point15        6.164       3.605
point16        5.214       4.320
point17        7.499       8.174
point18        0.160       3.729
point19        6.971       7.669
point20        8.875       8.117
point21        7.897       8.036
point22        0.166       3.932
point23        9.601       7.215
point24        8.147       8.375
point25        7.745       5.283
point26        6.557       3.924
point27        4.851       4.026
point28        4.346       1.935
point29        1.290       4.876
point30        7.912       6.856
point31        3.026       0.001
point32        4.463       8.148
point33        2.080       3.360
point34        2.816       1.565
point35        6.203       8.599
point36        6.153       3.857
point37        7.633       7.315
point38        0.798       0.829
point39        0.192       2.454
point40        3.341       3.038
point41        8.349       3.863
point42        0.531       5.056
point43        1.919       2.345
point44        9.347       7.861
point45        9.435       1.718
point46        1.817       2.026
point47        5.323       5.065
point48        2.108       6.382
point49        2.959       5.533
point50        0.342       8.075
point51        3.551       9.887
point52        8.136       7.153
point53        2.435       5.417
point54        5.019       1.108
point55        1.686       2.869
point56        8.416       8.060
point57        5.631       1.570
point58        2.715       0.339
point59        7.184       7.050
point60        1.577       7.522
point61        9.554       7.855
point62        2.566       0.880
point63        7.997       2.193
point64        9.556       8.034
point65        4.599       2.022
point66        8.957       5.241
point67        8.720       9.062
point68        3.519       3.229
point69        9.245       6.341
point70        4.141       2.687
point71        2.817       7.412
point72        6.571       4.559
point73        3.806       4.961
point74        8.999       5.323
point75        0.602       3.960
point76        4.302       9.600
point77        2.717       8.777
point78        1.872       4.952
point79        3.734       4.279
point80        9.519       5.120
point81        4.618       4.190
point82        0.915       9.232
point83        9.528       6.382
point84        5.353       8.983
point85        0.462       4.249
point86        8.164       4.834
point87        8.054       3.067
point88        8.594       7.318
point89        5.609       3.905
point90        2.741       4.446
point91        9.810       2.776
point92        7.988       9.737
point93        5.810       7.965
point94        8.287       7.178
point95        4.339       2.033
point96        2.708       6.648
point97        1.108       8.467
point98        2.175       2.350
point99        8.997       8.511
point100       8.928       2.236

1. Largest empty square problem

A nice little problem is: find the largest square that covers no points. We want this:

Modeling this will help us with some more complicated problems later on. A high-level model can look like:

High-level model
\[\begin{align}\max_{{\color{darkred}x},{\color{darkred}y},{\color{darkred}s}}\>&{\color{darkred}{\mathit{Area}}}={\color{darkred}s}^2 \\ & {\color{darkblue}p}_{i,x} \le {\color{darkred}x} \>\mathbf{or}\> {\color{darkblue}p}_{i,x} \ge {\color{darkred}x}+{\color{darkred}s} \> \mathbf{or}\> {\color{darkblue}p}_{i,y} \le {\color{darkred}y} \>\mathbf{or}\> {\color{darkblue}p}_{i,y} \ge {\color{darkred}y}+{\color{darkred}s} && \forall i \\ & {\color{darkred}x},{\color{darkred}y} \in [0,{\color{darkblue}{\mathit{maxSize}}} - {\color{darkred}s}] \\ & {\color{darkred}s} \ge 0 \\ & i \in \{1,\dots,{\color{darkblue}n}\} \end{align}\]

High-level model

\[\begin{align}\max_{{\color{darkred}x},{\color{darkred}y},{\color{darkred}s}}\>&{\color{darkred}{\mathit{Area}}}={\color{darkred}s}^2 \\ & {\color{darkblue}p}_{i,x} \le {\color{darkred}x} \>\mathbf{or}\> {\color{darkblue}p}_{i,x} \ge {\color{darkred}x}+{\color{darkred}s} \> \mathbf{or}\> {\color{darkblue}p}_{i,y} \le {\color{darkred}y} \>\mathbf{or}\> {\color{darkblue}p}_{i,y} \ge {\color{darkred}y}+{\color{darkred}s} && \forall i \\ & {\color{darkred}x},{\color{darkred}y} \in [0,{\color{darkblue}{\mathit{maxSize}}} - {\color{darkred}s}] \\ & {\color{darkred}s} \ge 0 \\ & i \in \{1,\dots,{\color{darkblue}n}\} \end{align}\]

Basically: each point should be below, above, on the left, or right of the square. A practical implementation, using binary variables and big-M constraints to implement the $\mathbf{or}$ disjunction, can look like:

MIP Implementation
\[\begin{align}\max_{{\color{darkred}x},{\color{darkred}y},{\color{darkred}s}}\>&{\color{darkred}s} \\ & {\color{darkblue}p}_{i,x} \le {\color{darkred}x} + {\color{darkblue}M} {\color{darkred}\delta}_{i,1} && \forall i \\ & p_{i,x} \ge {\color{darkred}x}+{\color{darkred}s} - {\color{darkblue}M}{\color{darkred}\delta}_{i,2} && \forall i \\ & p_{i,y} \le {\color{darkred}y} +{\color{darkblue}M} {\color{darkred}\delta}_{i,3}&& \forall i \\ & p_{i,y} \ge {\color{darkred}y}+{\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,4} && \forall i \\ & \sum_k {\color{darkred}\delta}_{i,k} = 3 && \forall i\\ & {\color{darkred}x}+ {\color{darkred}s}\le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}y} + {\color{darkred}s} \le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}x},{\color{darkred}y},{\color{darkred}s} \ge 0 \\ & {\color{darkred}\delta}_{i,k} \in \{0,1\} \\ & {\color{darkblue}M} = {\color{darkblue}{\mathit{maxSize}}} \\ & i \in \{1,\dots,{\color{darkblue}n}\} \end{align}\]

MIP Implementation

\[\begin{align}\max_{{\color{darkred}x},{\color{darkred}y},{\color{darkred}s}}\>&{\color{darkred}s} \\ & {\color{darkblue}p}_{i,x} \le {\color{darkred}x} + {\color{darkblue}M} {\color{darkred}\delta}_{i,1} && \forall i \\ & p_{i,x} \ge {\color{darkred}x}+{\color{darkred}s} - {\color{darkblue}M}{\color{darkred}\delta}_{i,2} && \forall i \\ & p_{i,y} \le {\color{darkred}y} +{\color{darkblue}M} {\color{darkred}\delta}_{i,3}&& \forall i \\ & p_{i,y} \ge {\color{darkred}y}+{\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,4} && \forall i \\ & \sum_k {\color{darkred}\delta}_{i,k} = 3 && \forall i\\ & {\color{darkred}x}+ {\color{darkred}s}\le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}y} + {\color{darkred}s} \le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}x},{\color{darkred}y},{\color{darkred}s} \ge 0 \\ & {\color{darkred}\delta}_{i,k} \in \{0,1\} \\ & {\color{darkblue}M} = {\color{darkblue}{\mathit{maxSize}}} \\ & i \in \{1,\dots,{\color{darkblue}n}\} \end{align}\]

Notes:

Maximizing ${\color{darkred}s}^2$ would make the problem nonconvex nonlinear. It is much easier to maximize ${\color{darkred}s}$. We can do this because this is a monotonic transformation. The model is now linear, so we can use any MIP solver.
$({\color{darkred}x},{\color{darkred}y})$ indicates the left-lower corner of the square in this model. Obviously, we could also have used the center of the square as the anchor point for the location.
The ${\color{darkred}\delta}$ variables relax the corresponding constraint when ${\color{darkred}\delta}_{i,k}=1$. We want at least one constraint to hold, so one ${\color{darkred}\delta}_{i,k}$ should be zero. Note that we could have used \[\sum_k {\color{darkred}\delta}_{i,k} \le 3\] instead of the equality. Is that better (i.e., faster)? I don't know. I would guess it does not make much difference.
We want the size of the big-${\color{darkblue}M}$ constants to be relatively small. If your data has very large coordinates, simply scale the points.
This solves quite easily. I see:

--- Cplex Time: 0.06sec (det. 56.38 ticks)

Proven optimal solution
MIP Solution: 2.900147 (484 iterations, 0 nodes)

so no branch & bound nodes were needed. This is a very good performance.
The solution is not unique: we can shift the square around a bit.

2. Largest diamond problem

Here we rotate the square by 45°. The result is a diamond shape.

This problem has a special interpretation. If we take $(x,y)$ as the center, then the rim is at the L1 or Manhattan distance from the center to the closest data point. The Manhattan distance between points $(x_1,y_1)$ and $(x_2,y_2)$ is defined by: \[d = |x_1-x_2| + |y_1-y_2|\] We will try to maximize the smallest L1 distance between the center of the diamond and any data point.

The optimization model can look like:

High-Level Model
\[\begin{align}\max_{{\color{darkred}x},{\color{darkred}y},{\color{darkred}d}}\>&{\color{darkred}d} \\ & {\color{darkred}d} \le \| {\color{darkblue}p}_{i,x} - {\color{darkred}x} \| +\| {\color{darkblue}p}_{i,y} - {\color{darkred}y} \| && \forall i\\ &{\color{darkblue}p}_{i,c} -{\color{darkred}d} \ge 0 && \forall i,c\\&{\color{darkblue}p}_{i,c} +{\color{darkred}d} \le{\color{darkblue}{\mathit{maxSize}}}&& \forall i,c\\ & c \in \{x,y\} \\ & i \in \{1,\dots,{\color{darkblue}n}\} \end{align}\]

High-Level Model

\[\begin{align}\max_{{\color{darkred}x},{\color{darkred}y},{\color{darkred}d}}\>&{\color{darkred}d} \\ & {\color{darkred}d} \le | {\color{darkblue}p}_{i,x} - {\color{darkred}x} | +| {\color{darkblue}p}_{i,y} - {\color{darkred}y} | && \forall i\\ &{\color{darkblue}p}_{i,c} -{\color{darkred}d} \ge 0 && \forall i,c\\&{\color{darkblue}p}_{i,c} +{\color{darkred}d} \le{\color{darkblue}{\mathit{maxSize}}}&& \forall i,c\\ & c \in \{x,y\} \\ & i \in \{1,\dots,{\color{darkblue}n}\} \end{align}\]

This model looks very different from the previous one. We can implement absolute values using variable splitting, in which a variable is split into its positive and negative parts. In convex cases (like minimizing the absolute value), we don't need binary variables to ensure that one of the positive and negative variables is zero. In this case, however, as we are essentially maximizing the absolute value, we need to address this extra complication.

MIP Implementation
\[\begin{align}\max\>&{\color{darkred}d} \\ & {\color{darkred}\Delta}_{i,c}^+ -{\color{darkred}\Delta}_{i,c}^- = {\color{darkblue}p}_{i,c} - {\color{darkred}{\mathit{center}}}_c && \forall i,c \\ &{\color{darkred}\Delta}_{i,c}^+ \le {\color{darkblue}M} {\color{darkred}\delta}_{i,c}&& \forall i,c \\ & {\color{darkred}\Delta}_{i,c}^- \le {\color{darkblue}M} (1-{\color{darkred}\delta}_{i,c})&& \forall i,c \\ & {\color{darkred}d} \le \sum_c \left( {\color{darkred}\Delta}_{i,c}^+ +{\color{darkred}\Delta}_{i,c}^- \right) && \forall i\\ &{\color{darkblue}p}_{i,c} -{\color{darkred}d} \ge 0 && \forall i,c\\&{\color{darkblue}p}_{i,c} +{\color{darkred}d} \le{\color{darkblue}{\mathit{maxSize}}}&& \forall i,c \\ & {\color{darkred}\Delta}_{i,c}^+, {\color{darkred}\Delta}_{i,c}^- \ge 0 \\ & {\color{darkred}\delta}_{i,c} \in \{0,1\} \\ & c \in \{x,y\}\\ & i \in \{1,\dots,{\color{darkblue}n}\}\\ & {\color{darkblue}M} = {\color{darkblue}{\mathit{maxSize}}}\end{align}\]

MIP Implementation

\[\begin{align}\max\>&{\color{darkred}d} \\ & {\color{darkred}\Delta}_{i,c}^+ -{\color{darkred}\Delta}_{i,c}^- = {\color{darkblue}p}_{i,c} - {\color{darkred}{\mathit{center}}}_c && \forall i,c \\ &{\color{darkred}\Delta}_{i,c}^+ \le {\color{darkblue}M} {\color{darkred}\delta}_{i,c}&& \forall i,c \\ & {\color{darkred}\Delta}_{i,c}^- \le {\color{darkblue}M} (1-{\color{darkred}\delta}_{i,c})&& \forall i,c \\ & {\color{darkred}d} \le \sum_c \left( {\color{darkred}\Delta}_{i,c}^+ +{\color{darkred}\Delta}_{i,c}^- \right) && \forall i\\ &{\color{darkblue}p}_{i,c} -{\color{darkred}d} \ge 0 && \forall i,c\\&{\color{darkblue}p}_{i,c} +{\color{darkred}d} \le{\color{darkblue}{\mathit{maxSize}}}&& \forall i,c \\ & {\color{darkred}\Delta}_{i,c}^+, {\color{darkred}\Delta}_{i,c}^- \ge 0 \\ & {\color{darkred}\delta}_{i,c} \in \{0,1\} \\ & c \in \{x,y\}\\ & i \in \{1,\dots,{\color{darkblue}n}\}\\ & {\color{darkblue}M} = {\color{darkblue}{\mathit{maxSize}}}\end{align}\]

Notes:

It is important that the big-Ms are relatively small.
This is also an easy MIP model. It was solved during preprocessing:

--- Cplex Time: 0.09sec (det. 119.55 ticks)

Proven optimal solution
MIP Solution: 1.972049 (438 iterations, 0 nodes)
Here $({\color{darkred}x},{\color{darkred}y})$ is the center of the diamond shape and not the lower-left corner as in the previous model.
This model finds the point that is farthest from the closest data point, measured as an L1 distance. We only need to allow the point to be close to the border of the $[0,{\color{darkblue}{\mathit{maxSize}}}] \times [0,{\color{darkblue}{\mathit{maxSize}}}]$ box. In the above model we want to have the whole diamond inside the outer box. So we are a bit too strict for the "point" interpretation of the model.

3. Largest empty rectangle problem

A relatively small extension is to find the largest empty rectangle. For this, we have two different sides ${\color{darkred}s}_x$ and ${\color{darkred}s}_y$. If we try to maximize the area ${\color{darkred}s}_x \cdot {\color{darkred}s}_y$ directly, we end up with a nonconvex quadratic model. However, there is a trick that makes this a convex problem. A constraint like $\sum_i x_i^2 \le y\cdot z$ with $y,z\ge 0$ is a rotated second order cone. Solvers like Cplex, Gurobi and others can handle these directly. So we use this special quadratic constraint and formulate:

MISOCP Implementation
\[\begin{align}\max \>&{\color{darkred}t} \\ & {\color{darkred}t}^2 \le {\color{darkred}s}_x \cdot {\color{darkred}s}_y \\ & {\color{darkblue}p}_{i,x} \le {\color{darkred}x} + {\color{darkblue}M} {\color{darkred}\delta}_{i,1} && \forall i \\ & p_{i,x} \ge {\color{darkred}x}+{\color{darkred}s}_x - {\color{darkblue}M}{\color{darkred}\delta}_{i,2} && \forall i \\ & p_{i,y} \le {\color{darkred}y} +{\color{darkblue}M} {\color{darkred}\delta}_{i,3}&& \forall i \\ & p_{i,y} \ge {\color{darkred}y}+{\color{darkred}s}_y - {\color{darkblue}M} {\color{darkred}\delta}_{i,4} && \forall i \\ & \sum_k {\color{darkred}\delta}_{i,k} = 3 && \forall i\\ & {\color{darkred}x} + {\color{darkred}s}_x \le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}y} + {\color{darkred}s}_y\le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}x},{\color{darkred}y},{\color{darkred}s}_x,{\color{darkred}s}_y\ge 0 \\&{\color{darkred}t}\> \mathbf{free}\\ & {\color{darkred}\delta}_{i,k} \in \{0,1\} \\ & {\color{darkblue}M} = {\color{darkblue}{\mathit{maxSize}}} \end{align}\]

MISOCP Implementation

\[\begin{align}\max \>&{\color{darkred}t} \\ & {\color{darkred}t}^2 \le {\color{darkred}s}_x \cdot {\color{darkred}s}_y \\ & {\color{darkblue}p}_{i,x} \le {\color{darkred}x} + {\color{darkblue}M} {\color{darkred}\delta}_{i,1} && \forall i \\ & p_{i,x} \ge {\color{darkred}x}+{\color{darkred}s}_x - {\color{darkblue}M}{\color{darkred}\delta}_{i,2} && \forall i \\ & p_{i,y} \le {\color{darkred}y} +{\color{darkblue}M} {\color{darkred}\delta}_{i,3}&& \forall i \\ & p_{i,y} \ge {\color{darkred}y}+{\color{darkred}s}_y - {\color{darkblue}M} {\color{darkred}\delta}_{i,4} && \forall i \\ & \sum_k {\color{darkred}\delta}_{i,k} = 3 && \forall i\\ & {\color{darkred}x} + {\color{darkred}s}_x \le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}y} + {\color{darkred}s}_y\le {\color{darkblue}{\mathit{maxSize}}} \\& {\color{darkred}x},{\color{darkred}y},{\color{darkred}s}_x,{\color{darkred}s}_y\ge 0 \\&{\color{darkred}t}\> \mathbf{free}\\ & {\color{darkred}\delta}_{i,k} \in \{0,1\} \\ & {\color{darkblue}M} = {\color{darkblue}{\mathit{maxSize}}} \end{align}\]

Notes:

Much is borrowed from the first linear model.
Quadratic models, including SOCP models, typically require a bit more time to solve than linear models. Here we see:

--- Cplex Time: 1.09sec (det. 1570.77 ticks)

Solution satisfies tolerances
MIP Solution: 3.410958 (49129 iterations, 3404 nodes)
Of course, we could have solved this with the simpler objective: $\max {\color{darkred}s}_x \cdot {\color{darkred}s}_y$ using a global solver.
This is called a convex problem. The binary variables obviously introduce a non-convexity. But, in mathematical programming, we say an integer problem is convex if relaxations are convex.
Mosek refuses this model. I have not investigated this further. Other conic solvers accept this.

The solution looks like:

4. Largest empty circle problem

Here we want to find the largest circle with center ${\color{darkred}x},{\color{darkred}y}$ and radius ${\color{darkred}r}$ that does not cover any of the points ${\color{darkblue}p}_i$.

This is the same as finding the point that is farthest away from the closest point. Sometimes this is called the "toxic waste dump problem." It is often solved using Voronoi diagrams. Formulated as an optimization problem, this is, unfortunately, nonconvex:

Nonconvex NLP Implementation
\[\begin{align}\max \>&{\color{darkred}r} \\ & {\color{darkred}\Delta}_{i,c} ={\color{darkblue}p}_{i,c} -{\color{darkred}{\mathit{center}}}_{c}\\ & \sum_c {\color{darkred}\Delta}_{i,c}^2 \ge {\color{darkred}r}^2 && \forall i \\ & {\color{darkred}{\mathit{center}}}_{c} -{\color{darkred}r} \ge 0 \\ & {\color{darkred}{\mathit{center}}}_{c} +{\color{darkred}r} \le {\color{darkblue}{\mathit{maxSize}}} \\ & c \in \{x,y\} \end{align}\]

Nonconvex NLP Implementation

\[\begin{align}\max \>&{\color{darkred}r} \\ & {\color{darkred}\Delta}_{i,c} ={\color{darkblue}p}_{i,c} -{\color{darkred}{\mathit{center}}}_{c}\\ & \sum_c {\color{darkred}\Delta}_{i,c}^2 \ge {\color{darkred}r}^2 && \forall i \\ & {\color{darkred}{\mathit{center}}}_{c} -{\color{darkred}r} \ge 0 \\ & {\color{darkred}{\mathit{center}}}_{c} +{\color{darkred}r} \le {\color{darkblue}{\mathit{maxSize}}} \\ & c \in \{x,y\} \end{align}\]

Notes:

We could add the bound ${\color{darkred}r}\ge 0$.
We could get rid of one nonlinearity by saying \[\begin{align}\max\>&{\color{darkred}R} \\ &\sum_c {\color{darkred}\Delta}_{i,c}^2 \ge {\color{darkred}R} \end{align}\] where ${\color{darkred}R}={\color{darkred}r}^2$.
I'd like the circle to be completely inside the box. If we are looking for just the farthest point, we can drop that.
I added bounds of $0$ and ${\color{darkblue}{\mathit{maxSize}}}$ to the variables ${\color{darkred}r}$ and ${\color{darkred}{\mathit{center}}}_{c}$. When using global solvers, it is a good idea to have reasonable bounds on the variables.
This problem was solved very fast using the global solver Baron:

Wall clock time: 0.66
Total CPU time used: 0.53

Total no. of BaR iterations: 33
Best solution found at node: 1
Max. no. of nodes in memory: 5

5. The five largest empty rectangles problem

Here we try to find the five largest different rectangles. The first question to ask: What is different? If we don't pay attention to this, we may find the original largest rectangle (or a very close alternative). We really want the next rectangle to be significantly different from the previous ones. I used the position of the left lower corner $(x,y)$ as a measure to differentiate rectangles: the position must be 3 units away from the positions of other rectangles, using the Manhattan distance. That minimum distance measure is, of course, subjective and dependent on the problem at hand. It required a few trials to get an appealing picture. This is an interesting modeling exercise. I get the following results:

Note that, as expected, each next rectangle we find is smaller than the one before. I used a loop for this, but we can also create a single, larger problem that finds all these rectangles in one solve.

In practical modeling, instead of just showing the optimal solution, it is often useful to present a few good solutions that are sufficiently different. In this case, we see some "holes" in the point coverage that I did not notice before.

Conclusions

The problem of shapes not covering points is basically inherently nonconvex. Finding shapes containing points is often convex, so it makes sense that the opposite problem is not. In some cases, we can handle the nonconvexity by using binary variables. In others, we really end up with problems that require global nonlinear solvers. All models solve very quickly.

I find it interesting that basically the same problem, finding the largest empty shape, is modeled differently depending very much on the shape and solved in very different ways: linear MIPs, an MISOCP model, and a nonconvex quadratic problem. In most cases, there were not that many concepts we could reuse.

Q: how would a model for a rotated square (with endogenous angle) look like? This would be a generalization of models 1 and 2.

GAMS Code

$onText

Find largest shapes that don't cover any points.

Shapes:

1. Square (MIP)

2. Diamond (MIP)

3. Rectangle (MISOCP)

4. Circle (nonconvex NLP)

5. 5 different "next best" rectangles

$offtext

*---------------------------------------------------------------

* options

*---------------------------------------------------------------

* 0: no plot

* 1: produce HTML plot)

$set htmlplot 1

option qcp=cplex,mip=cplex,miqcp=cplex,nlp=baron;

option seed=12345;

*---------------------------------------------------------------

* data

*---------------------------------------------------------------

scalar maxsize 'size of the box' / 10 /;

sets

i 'points' /point1*point100/

cr 'circle attributes' /x,y,r/

c(cr) 'coordinates' /x,y/

;

parameter p(i,c) 'points';

p(i,c) = uniform(0,maxsize);

display p;

*---------------------------------------------------------------

* 1. largest empty square model (MIP)

*---------------------------------------------------------------

set case 'for comparison' /less,greater/;

positive variable sq(*) 'square location and size';

binary variable d(i,c,case) 'd=1 means: relax no-overlap constraint';

variable z 'objective';

Equations

emptysq1(i,c) 'no-overlap constraint: less-than version'

emptysq2(i,c) 'no-overlap constraint: greater-than version'

sumd(i) '(at least) one no-overlap constraint must hold'

maxsq 'stay inside our box'

objsq 'objective'

;

objsq.. z =e= sq('side');

maxsq(c).. sq(c) + sq('s') =l= maxsize;

emptysq1(i,c).. p(i,c) =l= sq(c) + d(i,c,'less')*maxsize;

emptysq2(i,c).. p(i,c) =g= sq(c) + sq('side') - d(i,c,'greater')*maxsize;

sumd(i).. sum((c,case),d(i,c,case)) =e= 3;

model emptysquare /all/;

solve emptysquare maximizing z using mip;

parameter square(*) 'optimal values';

square(c) = sq.l(c);

square('side') = sq.l('side');

square('area') = sqr(sq.l('side'));

display square;

*---------------------------------------------------------------

* 2. largest empty diamond model (MIP)

*---------------------------------------------------------------

set pm 'var splitting for abs()' /'+','-'/;

positive variables

diamnd(*) 'center + dist of diamond'

absdiffd(i,c,*) 'Δ(i,c) = |p(i,c)-center(c)|'

;

binary variable delta 'for variable splitting';

variable dmin 'Manhattan distance to closest point';

equation

eabsdiff1(i,c) '|p(i,c)-center(c)|'

eabsdiff2(i,c) '|p(i,c)-center(c)|'

eabsdiff3(i,c) '|p(i,c)-center(c)|'

smallestd 'smallest distance'

elbnd(c) 'diamond must be inside box'

eubnd(c) 'diamond must be inside box'

;

eabsdiff1(i,c).. absdiffd(i,c,'+')-absdiffd(i,c,'-') =e= p(i,c)-diamnd(c);

eabsdiff2(i,c).. absdiffd(i,c,'+') =l= delta(i,c)*Maxsize;

eabsdiff3(i,c).. absdiffd(i,c,'-') =l= (1-delta(i,c))*Maxsize;

smallestd(i).. dmin =l= sum((c,pm),absdiffd(i,c,pm));

elbnd(c).. diamnd(c) - dmin =g= 0;

eubnd(c).. diamnd(c) + dmin =l= maxsize;

model emptydiamond /all-emptysquare/;

solve emptydiamond maximizing dmin using mip;

parameter diamond(*) 'optimal values';

diamond(c) = diamnd.l(c);

diamond('L1 dist') = dmin.l;

diamond('diagonal') = 2*dmin.l;

diamond('side') = sqrt(2*sqr(dmin.l));

diamond('area') = 2*sqr(dmin.l);

display diamond;

*---------------------------------------------------------------

* largest empty rectangle model (MISOCP)

*---------------------------------------------------------------

set a 'attribute' /pos,side/;

positive variable rect(a,c) 'rectangle location and size';

binary variable d(i,c,case) 'd=1 means: relax no-overlap constraint';

variables t 'to be maximized';

Equations

maxrect(c) 'stay inside box'

emptyrect1(i,c) 'no-overlap: less-than version'

emptyrect2(i,c) 'no-overlap: greater-than version'

socp '2nd order cone constraint'

;

maxrect(c).. rect('pos',c) + rect('side',c) =l= maxsize;

emptyrect1(i,c).. p(i,c) =l= rect('pos',c) + d(i,c,'less')*maxsize;

emptyrect2(i,c).. p(i,c) =g= rect('pos',c) + rect('side',c) - d(i,c,'greater')*maxsize;

socp.. sqr(t) =l= rect('side','x')*rect('side','y');

model maxrectangle /maxrect,emptyrect1,emptyrect2,socp,sumd/;

solve maxrectangle maximizing t using miqcp;

parameter rectangle(a,c) 'optimal rectangle';

saverect(a,c) = rect.l(a,c);

display rectangle;

*---------------------------------------------------------------

* largest empty circle model (nonconvex NLP)

*---------------------------------------------------------------

Variables

circle(cr) 'center+radius'

diff(i,c) 'intermediate variables'

;

equations

calcdiff(i,c) 'compute diff'

emptycircle(i) 'no points inside circle'

objcircle 'objective'

bnd1(c) 'stay inside box'

bnd2(c) 'stay inside box'

;

objcircle.. z =e= circle('r');

calcdiff(i,c).. diff(i,c) =e= p(i,c)-circle(c);

emptycircle(i).. sum(c, sqr(diff(i,c))) =g= sqr(circle('r'));

bnd1(c).. circle(c) =g= circle('r');

bnd2(c).. circle(c) =l= maxsize-circle('r');

circle.lo(cr) = 0;

circle.up(cr) = maxsize;

model maxcircle /objcircle,calcdiff,emptycircle,bnd1,bnd2/;

solve maxcircle maximizing z using nlp;

display circle.l;

*---------------------------------------------------------------

* five best rectangles

*---------------------------------------------------------------

scalar difftol 'min difference between rectangles |Δx|+|Δy|>difftol' /3/;

set

run 'solve runs' /run1*run5/

prevrun(run) 'previous runs'

;

parameter prevrect(run,*,*);

binary variable bindiff(run,c);

positive variable absdiff(run,c,pm);

equation

difference1(run,c)

difference2(run,c)

difference3(run,c)

difference4(run)

;

difference1(prevrun,c).. absdiff(prevrun,c,'+')-absdiff(prevrun,c,'-') =e= prevrect(prevrun,'pos',c) - rect('pos',c);

difference2(prevrun,c).. absdiff(prevrun,c,'+') =l= bindiff(prevrun,c)*maxsize;

difference3(prevrun,c).. absdiff(prevrun,c,'-') =l= (1-bindiff(prevrun,c))*maxsize;

difference4(prevrun).. sum((c,pm),absdiff(prevrun,c,pm)) =g= difftol;

model maxrectangle2 /maxrectangle,difference1,difference2,difference3,difference4/;

prevrect(run,'pos',c) = 0;

prevrun(run) = no;

loop(run,

solve maxrectangle2 maximizing t using miqcp;

display rect.l;

prevrect(run,'pos',c) = rect.l('pos',c);

prevrect(run,'side',c) = rect.l('side',c);

prevrun(run) = yes;

display prevrect;

);

*---------------------------------------------------------------------

* visualization

*---------------------------------------------------------------------

$set html plot.html

$set data data.js

$if %htmlplot%==0 $goto skipplot

file fdata /%data%/; put fdata;

* points

put "p=["/;

loop(i,

put " {i:'",i.tl:0,"',x:",p(i,'x'):0:4,",y:",p(i,'y'):0:4,"},"/;

);

put "]"/;

put "pnts={n:",(card(i)):0:0,

",map:'[0,",maxsize:0:0,"]×[0,",maxsize:0:0,"]'",

",minx:",(smin(i,p(i,'x'))):0:3,

",miny:",(smin(i,p(i,'y'))):0:3,

",maxx:",(smax(i,p(i,'x'))):0:3,

",maxy:",(smax(i,p(i,'y'))):0:3,

",avgx:",(sum(i,p(i,'x'))/card(i)):0:3,

",avgy:",(sum(i,p(i,'y'))/card(i)):0:3,

"}"/;

put "sq={x:",square('x'):0:3,

",y:",square('y'):0:3,

",side:",square('side'):0:3,

",area:",(sqr(square('side'))):0:3,

"}"/;

put "d={x:",diamond('x'):0:3,

",y:",diamond('y'):0:3,

",area:",diamond('area'):0:3,

",side:",diamond('side'):0:3,

",diagonal:",diamond('diagonal'):0:3,

"}"/;

put "r={x:",rectangle('pos','x'):0:3,

",y:",rectangle('pos','y'):0:3,

",sx:",rectangle('side','x'):0:3,

",sy:",rectangle('side','y'):0:3,

",a:",(rectangle('side','x')*rectangle('side','y')):0:3,"}"/;

put "circle={x:",circle.l('x'):0:3,",y:",circle.l('y'):0:3,",r:",circle.l('r'):0:3,

",a:",(pi*sqr(circle.l('r'))):0:3,"}"/;

put "rs=["/;

loop(run,

put " {x:",prevrect(run,'pos','x'):0:3,

",y:",prevrect(run,'pos','y'):0:3,

",sx:",prevrect(run,'side','x'):0:3,

",sy:",prevrect(run,'side','y'):0:3,

",a:",(prevrect(run,'side','x')*prevrect(run,'side','y')):0:3,

"},"/;

);

put "]"/;

putclose;

$onecho > %html%

<html>

<style>

.brdr {

border: 1px solid black;

border-collapse: collapse;

}

th, td {

padding-left: 4px;

padding-right: 4px;

}

</style>

<tr><th colspan=2 class="brdr">points</th></tr>

</table>

</td><td>

</td></tr></table>

<h1>Largest Empty Square</h1>

<tr><th colspan=2 class="brdr">square</th></tr>

</table>

</td><td>

</td></tr></table>

<h1>Largest Empty Diamond</h1>

<tr><th colspan=2 class="brdr">diamond</th></tr>

<tr><td class="brdr">diagonal</td><td class="brdr" id="d_diag"></td></tr>

</table>

</td><td>

</td></tr></table>

<h1>Largest Empty Rectangle</h1>

<tr><th colspan=2 class="brdr">rectangle</td></tr>

</table>

</td><td>

</td></tr></table>

<h1>Largest Empty Circle</h1>

<tr><th colspan=2 class="brdr">circle</td></tr>

<tr><td class="brdr">radius</td><td class="brdr" id="c_r">0</td></tr>

</table>

</td><td>

</td></tr></table>

<h1>Five largest rectangles</h1>

<tr><th class="brdr">rectangles</th><th class="brdr">rect1</th><th class="brdr">rect2</th><th class="brdr">rect3</th><th class="brdr">rect4</th><th class="brdr">rect5</th></tr>

</table>

</td><td>

</td></tr></table>

document.getElementById("d_n").innerHTML = pnts['n'];

document.getElementById("d_map").innerHTML = pnts['map'];

document.getElementById("d_minx").innerHTML = pnts['minx'];

document.getElementById("d_maxx").innerHTML = pnts['maxx'];

document.getElementById("d_miny").innerHTML = pnts['miny'];

document.getElementById("d_maxy").innerHTML = pnts['maxy'];

document.getElementById("d_avgx").innerHTML = pnts['avgx'];

document.getElementById("d_avgy").innerHTML = pnts['avgy'];

document.getElementById("sq_x").innerHTML = sq['x'];

document.getElementById("sq_y").innerHTML = sq['y'];

document.getElementById("sq_s").innerHTML = sq['side'];

document.getElementById("sq_area").innerHTML = sq['area'];

document.getElementById("d_x").innerHTML = d['x'];

document.getElementById("d_y").innerHTML = d['y'];

document.getElementById("d_diag").innerHTML = d['diagonal'];

document.getElementById("d_side").innerHTML = d['side'];

document.getElementById("d_area").innerHTML = d['area'];

document.getElementById("r_x").innerHTML = r['x'];

document.getElementById("r_y").innerHTML = r['y'];

document.getElementById("r_s1").innerHTML = r['sx'];

document.getElementById("r_s2").innerHTML = r['sy'];

document.getElementById("r_area").innerHTML = r['a'];

document.getElementById("c_x").innerHTML = circle['x'];

document.getElementById("c_y").innerHTML = circle['y'];

document.getElementById("c_r").innerHTML = circle['r'];

document.getElementById("c_area").innerHTML = circle['a'];

document.getElementById("r_x1").innerHTML = rs[0]['x'];

document.getElementById("r_y1").innerHTML = rs[0]['y'];

document.getElementById("r_sx1").innerHTML = rs[0]['sx'];

document.getElementById("r_sy1").innerHTML = rs[0]['sy'];

document.getElementById("r_a1").innerHTML = rs[0]['a'];

document.getElementById("r_x2").innerHTML = rs[1]['x'];

document.getElementById("r_y2").innerHTML = rs[1]['y'];

document.getElementById("r_sx2").innerHTML = rs[1]['sx'];

document.getElementById("r_sy2").innerHTML = rs[1]['sy'];

document.getElementById("r_a2").innerHTML = rs[1]['a'];

document.getElementById("r_x3").innerHTML = rs[2]['x'];

document.getElementById("r_y3").innerHTML = rs[2]['y'];

document.getElementById("r_sx3").innerHTML = rs[2]['sx'];

document.getElementById("r_sy3").innerHTML = rs[2]['sy'];

document.getElementById("r_a3").innerHTML = rs[2]['a'];

document.getElementById("r_x4").innerHTML = rs[3]['x'];

document.getElementById("r_y4").innerHTML = rs[3]['y'];

document.getElementById("r_sx4").innerHTML = rs[3]['sx'];

document.getElementById("r_sy4").innerHTML = rs[3]['sy'];

document.getElementById("r_a4").innerHTML = rs[3]['a'];

document.getElementById("r_x5").innerHTML = rs[4]['x'];

document.getElementById("r_y5").innerHTML = rs[4]['y'];

document.getElementById("r_sx5").innerHTML = rs[4]['sx'];

document.getElementById("r_sy5").innerHTML = rs[4]['sy'];

document.getElementById("r_a5").innerHTML = rs[4]['a'];

px = p.map((x) => x['x'])

py = p.map((x) => x['y'])

var data1 = {

x: px,

y: py,

mode: 'markers',

type: 'scatter',

name: 'data points',

color: 'darkblue'

};

var layout0 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

}

var layoutSquare = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

shapes: [{type:'rect',

xref:'x',

yref:'y',

x0:sq['x'],

y0:sq['y'],

x1:sq['x']+sq['side'],

y1:sq['y']+sq['side'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

name:'largest square',

showlegend:true

}

]

}

var path = `M ${d['x']-0.5*d['diagonal']} ${d['y']} L ${d['x']} ${d['y']+0.5*d['diagonal']} L ${d['x']+0.5*d['diagonal']} ${d['y']} L ${d['x']} ${d['y']-0.5*d['diagonal']} Z`;

var layoutDiamond = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

shapes: [{type:'path',

path:path,

xref:'x',

yref:'y',

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

name:'largest diamond',

showlegend:true

}

]

}

var layout2 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

shapes: [{type:'rect',

xref:'x',

yref:'y',

x0:r['x'],

y0:r['y'],

x1:r['x']+r['sx'],

y1:r['y']+r['sy'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

name:'largest rectangle',

showlegend:true

}

]

}

var layout3 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

shapes: [{type:'circle',

xref:'x',

yref:'y',

x0:circle['x']-circle['r'],

y0:circle['y']-circle['r'],

x1:circle['x']+circle['r'],

y1:circle['y']+circle['r'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

name:'largest circle',

showlegend:true

}

]

}

var layout4 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

shapes: [{type:'rect',

xref:'x',

yref:'y',

x0:rs[0]['x'],

y0:rs[0]['y'],

x1:rs[0]['x']+rs[0]['sx'],

y1:rs[0]['y']+rs[0]['sy'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

name:'largest rectangle',

showlegend:true

{type:'rect',

xref:'x',

yref:'y',

x0:rs[1]['x'],

y0:rs[1]['y'],

x1:rs[1]['x']+rs[1]['sx'],

y1:rs[1]['y']+rs[1]['sy'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

{type:'rect',

xref:'x',

yref:'y',

x0:rs[2]['x'],

y0:rs[2]['y'],

x1:rs[2]['x']+rs[2]['sx'],

y1:rs[2]['y']+rs[2]['sy'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

{type:'rect',

xref:'x',

yref:'y',

x0:rs[3]['x'],

y0:rs[3]['y'],

x1:rs[3]['x']+rs[3]['sx'],

y1:rs[3]['y']+rs[3]['sy'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

{type:'rect',

xref:'x',

yref:'y',

x0:rs[4]['x'],

y0:rs[4]['y'],

x1:rs[4]['x']+rs[4]['sx'],

y1:rs[4]['y']+rs[4]['sy'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 255, 255, 0.5)',

]

}

var trc1 = [data1];

var options1 = {staticPlot: true, displayModeBar: false, responsive: false};

Plotly.newPlot('plotDiv0', trc1, layout0, options1);

Plotly.newPlot('plotDiv1', trc1, layoutSquare, options1);

Plotly.newPlot('plotDiv1a', trc1, layoutDiamond, options1);

Plotly.newPlot('plotDiv2', trc1, layout2, options1);

Plotly.newPlot('plotDiv3', trc1, layout3, options1);

Plotly.newPlot('plotDiv4', trc1, layout4, options1);

</script>

</html>

$offecho

executetool 'win32.ShellExecute "%html%"';

$label skipplot

Yet Another Math Programming Consultant

Monday, May 11, 2026

Largest empty shapes

Finding largest empty shapes

Data

1. Largest empty square problem

2. Largest diamond problem

3. Largest empty rectangle problem

4. Largest empty circle problem

5. The five largest empty rectangles problem

Conclusions