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) ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
