The right placement
January 11th, 2007 Paul Crowley
I wrote a short script to display the results of the 2004 US presidential elections as a bar chart, shown below (click for full-size image, which will make more sense):
The width of each state indicates the number of electoral college votes it has; the height the extent to which votes for one party exceeded votes for the other. The script is in Python and uses Cairo for drawing. The most challenging part of writing the script was…
…placing the labels that indicate the states.
You can’t just place the labels above/below the states, because some of the states are too narrow; the labels would overlap each other. So some labels need to be moved a little to one side or the other for them all to fit. It took me a little while to find a neat, satisfying algorithm for solving the problem. I’m sure it’s well known - it’s too simple not to be, and I’m sure I’m not the first with this problem - but it was fun to find.
Formally, we have some minimum distance between labels d, as well as smallest and largest positions they can be in (xmin and xmax). Each label wants to be in a position ai, so we want to choose each xi such that xi+1 ≥ xi + d, but |ai - xi| is small. The first requirement is strict and formally defined; the second is a fuzzier goodness factor.
My solution looks something like this:
offsets = [s.desired_x - i * d for i, s in enumerate(states)]
max_offset = x_max - d * (len(states) -1)
for i, s in enumerate(states):
new_offset = (max(offsets[:i+1]) + min(offsets[i:])) / 2.0
s.label_x = i * d + max(x_min, min(max_offset, new_offset))
How does it work?
Define a‘i = ai - id, and similarly x‘i = xi - id. Note that xi+1 ≥ xi + d if and only if x‘i+1 ≥ x‘i, in other words, x‘i has to be a monotonically nondecreasing sequence. Note also that |a‘i - x‘i| = |ai - xi|. So we’ve reduced the problem to finding a nondecreasing sequence that’s “close” in some sense to an arbitrary sequence. We make two non-decreasing sequences - “smallest point ahead of me” and “largest point behind me” - and average the two, which is guaranteed to yield another non-decreasing sequence. We then clamp this sequence between suitable min/max bounds, which again is guaranteed to yield a non-decreasing sequence. So it’s not hard to prove that labels never overlap.
How does it do on the goodness factor? As can be seen from the example, not too badly. Where there is plenty of room, the a‘i sequence is monotonically increasing, and so the other sequences we average cling to it; the result is that a‘i = x‘i. Where they are bunched together, the labels are spread out above, centered on the center of the labels, except at the ends where they fan out into the space available.
It’s not perfect; look in particular at New Jersey, which is making more room for its neighbours than it needs to. I may yet be able to improve on this. But it’s pleasing when you replace a nastily hacked-up thirty-line succesive approximation solution with a five-line solution that runs faster and produces more aesthetically pleasing results.
Entry Filed under: Technology, Politics
8 Comments Add your own
1. Jan Van lent | January 11th, 2007 at 8:38 pm
Probably overkill, but you can try to formulate this as a linear programming problem.
l1 norm
min \sumi ti
s.t.
xmin <= xi <= xmax
x{i+1} - xi >= d
-ti <= ai - xi <= ti
l\infty norm
min t
s.t.
xmin <= xi <= xmax
x{i+1} - xi >= d
-t <= ai - xi <= t
The l_2 norm would give a quadratic programming problem.
You can solve linear programming (and other) problems in python using for example cvxopt (http://www.ee.ucla.edu/~vandenbe/cvxopt/).
2. Jan Van lent | January 11th, 2007 at 8:39 pm
Untested python code for cvxopt:
d = 0.1
a = arange(n) + uniform(-d, d, n)
x = variable(n, ‘x’)
cxlow = ( xmin <= x )
cxup = ( x <= xmax )
c_d = [ x[i+1] - x[i] >= d for i in range(n-1) ]
t = variable(n, ‘t’)
ctlow = ( -t <= a - x )
ctup = ( a - x <= t )
lp1 = op(t, [ cxlow, cxup ] + cd + [ ctlow, ct_up ])
tinf = variable(1, ‘tinf’)
ctinflow = ( -tinf <= a - x )
ctinfup = ( a - x <= tinf )
lpinf = op(t, [ cxlow, cxup ] + cd + [ ctinflow, ctinf_up ])
lp_1.solve()
print x.value(), t.value()
3. Jan Van lent | January 11th, 2007 at 8:41 pm
Should have put code tags around formulas and code in the previous comments.
4. Paul Crowley | January 12th, 2007 at 5:56 pm
Yes, that seems a plausible approach. It occurred to me to cast this as a quadratic programming problem, but I hadn’t realised that such a convenient library for it existed - thanks for the pointer.
5. Paul Crowley | January 14th, 2007 at 3:24 pm
res = cvxopt.solvers.qp(
cvxopt.base.matrix([[1.0*int(i == j) for j in range(len(xpos))] for i in range(len(xpos))]),
cvxopt.base.matrix([[-x for x in xpos]]),
cvxopt.base.matrix([[0.0 for j in range(i)] + [-1.0, 1.0] + [0.0 for j in range(len(xpos)-i-1)] for i in range(len(xpos))]),
cvxopt.base.matrix([0.0 - xmin] + [-lwidth for i in range(len(xpos)-1)] + [xmax]))
6. Joachim Dahl | January 14th, 2007 at 4:49 pm
Hi Paul,
Specifying this problem as sparse will be much faster for
anything but toy-problems:
P = spmatrix(1.0, range(n), range(n))
q = matrix(xpos)
G = spmatrix(-1,range(n),range(n),(n+1,n)) + \
spmatrix( 1,range(1,n+1),range(n),(n+1,n))
h = matrix([0.0 - xmin] + [-lwidth for i in range(len(q)-1)] + [xmax]);
res = solvers.qp(P, q, G, h)
Best
Joachim
7. Joachim Dahl | January 14th, 2007 at 7:59 pm
You could also write your own KKT solver (ยง8.7 in the CVXOPT documentation). Essentially you would solve
the positive definite system
(I + G’diag(d)G)dx = r
where ‘d’ is a given positive scaling vector. That’s a
tridiagonal system, that can be solved very efficiently;
the next release of CVXOPT will include the banded
solvers in LAPACK suited for something like this.
Joachim
8. Paul Crowley | January 16th, 2007 at 11:49 am
Now I have three solutions:
I think the new algorithm produces the same optimal solution defined by the quadratic programming algorithm, but my current proof of it is so ugly that it’s too much work to set it out here and I don’t trust it. I think I’m on the tail of a more elegant proof now - I’ll make a new blog entry if I find it. In any case, this algorithm is much faster than the CVXOPT-based solution and produces results that look indistinguishable, without the artifacts of the first algorithm.
Fast algorithm:
def place_labels(xmin, xmax, lwidth, xpos): offsets = [x - i * lwidth for i, x in enumerate(xpos)] max_offset = xmax - lwidth * (len(xpos) -1) return [i * lwidth + max(xmin, min(max_offset, (max(offsets[:i+1]) + min(offsets[i:])) / 2.0)) for i in range(len(xpos))]CVXOPT-based algorithm:
import sys sys.path.append('/home/paul/path/lib/python') import cvxopt.base import cvxopt.solvers cvxopt.solvers.options['show_progress'] = False cvxopt.solvers.options['maxiters'] = 5000 #cvxopt.solvers.options['feastol'] = 1E-15 def place_labels(xmin, xmax, lwidth, xpos): n = len(xpos) res = cvxopt.solvers.qp( cvxopt.base.spmatrix(1.0, range(n), range(n)), -cvxopt.base.matrix(xpos), cvxopt.base.spmatrix(-1, range(n), range(n), (n+1, n)) + cvxopt.base.spmatrix(1, range(1, n+1), range(n), (n+1, n)), cvxopt.base.matrix([0.0 - xmin] + [-lwidth for i in range(len(xpos)-1)] + [xmax])) if res['status'] != 'optimal': raise Exception("Label optimizer failed: " + repr(res)) return [x for x in res['x']]New fast and perhaps optimal algorithm:
class PointGroup(object): def __init__(self, kind, n, x): self.kind = kind self.n = n self.x = x def place_labels(xmin, xmax, lwidth, xpos): groups = [PointGroup('l', 0, xmin)] + [PointGroup('m', 1, x) for x in xpos] + [PointGroup('r', 0, xmax + lwidth)] res = [] for g in groups: res.append(g) while len(res) >= 2 and res[-2].x + lwidth * res[-2].n >= res[-1].x: left = res[-2]; right = res[-1] del res[-2:] if left.kind == 'l' and right.kind == 'm': res.append(PointGroup('l', left.n +right.n, left.x)) elif left.kind == 'm' and right.kind == 'r': res.append(PointGroup('r', left.n + right.n, right.x - lwidth * left.n)) elif left.kind == 'm' and right.kind == 'm': res.append(PointGroup('m', left.n + right.n, (left.x * left.n + (right.x - left.n * lwidth) * right.n) / float(left.n + right.n))) else: raise Exception("Bonk!") return [g.x + lwidth * i for g in res for i in range(g.n)]Leave a Comment
Some HTML allowed:
<a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <code> <em> <i> <strike> <strong>
Trackback this post | Subscribe to the comments via RSS Feed