# [R] Find a rectangle of maximal area

Ray Brownrigg Ray.Brownrigg at ecs.vuw.ac.nz
Tue Mar 23 04:07:08 CET 2010

```Sorry, minor tweak to the algorithm in line 16 (thanks Barry).

Looks better now (at end again).

Ray

On Tue, 23 Mar 2010, Ray Brownrigg wrote:
> On Tue, 23 Mar 2010, Hans W Borchers wrote:
> > Barry Rowlingson <b.rowlingson <at> lancaster.ac.uk> writes:
> > > On Mon, Mar 22, 2010 at 4:28 PM, Hans W Borchers
> > >
> > > <hwborchers <at> googlemail.com> wrote:
> > > > Still I believe that a clever approach might be possible avoiding the
> > > > need to call a commercial solver. I am getting this hope from one of
> > > > Jon Bentley's articles in the series Programming Pearls.
> > >
> > > Is this the 'Largest Empty Rectangle' problem?
> > >
> > > http://en.wikipedia.org/wiki/Largest_empty_rectangle
> >
> > Dear Barry,
> >
> > thanks for this pointer. I never suspected this problem could have a name
> > of its own. Rethinking the many possible applications makes it clear: I
> > should have searched for it before.
> >
> > I looked in some of the references of the late 80s and found two
> > algorithms that appear to be appropriate for implementation in R. The
> > goal is to solve the problem for n=200 points in less than 10-15 secs.
>
> How about less than 2 seconds? [And 500 points in less than 15 seconds - on
> a 2-year-old DELL Optiplex GX755.]
>
> The implementation below (at end) loops over all 'feasible' pairs of x
> values, then selects the largest rectangle for each pair, subject to your
> specified constraints.  I have no idea if it implements a previously
> published algorithm.
>
> Other constraints are reasonably easily accommodated.
>
> HTH,
> Ray Brownrigg

N = 200
x <- runif(N)
y <- runif(N)
ox <- order(x)
x <- x[ox]; y <- y[ox]
x <- c(0, x, 1)
y <- c(0, y, 1)
plot(x, y, xlim=c(0, 1), ylim=c(0, 1), pch="*", col=2)
omaxa <- 0
for(i in 1:(length(x) - 1))
for(j in (i+1):length(x)) {
x1 <- x[i]
x2 <- x[j]
XX <- x2 - x1
if (XX > 0.5) break
yy <- c(0, y[(i + 1):(j - 1)], 1)
oyy <- order(yy)
yy <- yy[oyy]
dyy <- diff(yy)
whichdyy <- (dyy <= 0.5)  & (dyy >= 0.5*XX) & (dyy <= 2*XX)
wy1 <- yy[whichdyy]
if (length(wy1) > 0) {
wy2 <- yy[(1:length(dyy))[whichdyy] + 1]
k <- which.max(dyy[whichdyy])
maxa <- (x2 - x1)*(wy2[k] - wy1[k])
if (maxa > omaxa) {
omaxa <- maxa
mx1 <- x1
mx2 <- x2
my1 <- wy1[k]
my2 <- wy2[k]
}
}
}
print(omaxa)
lines(c(mx1, mx2, mx2, mx1, mx1), c(my1, my1, my2, my2, my1), col=2)

```