A fast Fortran solution may be found by:
require(Hmisc)
?largest.empty
Frank
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
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)
______________________________________________
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.
--
Frank E Harrell Jr Professor and Chairman School of Medicine
Department of Biostatistics Vanderbilt University
______________________________________________
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.
