[R] Find a rectangle of maximal area

[Ray Brownrigg]

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

> Thanks again, Hans Werner
>
> > I had a look at some of the references from Wikipedia, but they all
> > follow a similar pattern, one I have noticed in many computer science
> > journal articles:
> >
> >  1. State a problem that looks tricky.
> >  2. Say "We have an efficient algorithm for the problem stated in #1"
> >  3. Proceed to derive, using much algebra and theory, the efficient
> > algorithm. 4. Stop.
> >
> > The idea of actually producing some dirty, filthy, actual code to
> > implement their shiny algorithms never seems to cross their minds.
> >
> >  I also found a similar question from 2008 asked on the R-sig-geo
> > mailing list. That didn't get much help either!
> >
> > Sorry.
> >
> > Barry

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:j], 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]
      }
    }
  }
lines(c(mx1, mx2, mx2, mx1, mx1), c(my1, my1, my2, my2, my1), col=2)