technology from back to front

The right placement

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+1xi + d, but |aixi| 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 ai = aiid, and similarly xi = xiid. Note that xi+1xi + d if and only if xi+1xi, in other words, xi has to be a monotonically nondecreasing sequence. Note also that |aixi| = |aixi|. 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 ai sequence is monotonically increasing, and so the other sequences we average cling to it; the result is that ai = xi. 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.

by
Paul Crowley
on
11/01/07
  1. Jan Van lent
    on 11/01/07 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 )
    c
    xup = ( 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()

  2. Jan Van lent
    on 11/01/07 at 8:41 pm

    Should have put code tags around formulas and code in the previous comments.

  3. Jan Van lent
    on 11/01/07 at 8:38 pm

    Probably overkill, but you can try to formulate this as a linear programming problem.

    l1 norm
    min \sum
    i ti
    s.t.
    x
    min <= xi <= xmax
    x{i+1} – xi >= d
    -ti <= ai – xi <= ti

    l\infty norm
    min t
    s.t.
    x
    min <= 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/).

  4. 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.         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
    on 14/01/07 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
    on 14/01/07 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. Now I have three solutions:

    • the fast but not-quite-right solution I started with
    • the quadratic programming solution discussed in the comments
    • and a new linear-time algorithm

    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)]
    
    
 
 


6 − four =

2000-14 LShift Ltd, 1st Floor, Hoxton Point, 6 Rufus Street, London, N1 6PE, UK+44 (0)20 7729 7060   Contact us