[R] Find a rectangle of maximal area

[Ray Brownrigg]

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)