Comparing Packings

The Problem

Imagine we have generated two sets of points, like the following:

Packing-Example

Now we want some pairing of points on the left with points on the right, so that they “match up”. More specifically, we want a pairing that minimizes the total distance between paired points.

Table Form

number in A number in B
1 5
2 2
3 4
4 1
5 3

Image Form

Packing-Example

The Question

How do we find that pairing such that we minimize the sum of distances?

Mathematics

Let us define an ordering tuple to be a tuple . So and are ordering tuples, but is not.

Now suppose we have two sets of points, and , with .

We now want to find the ordering tuple that minimizes

where for each , we have a location and a “pairing” numbered at a location .

The Solution

Step 1: What does a solution even look like?

A solution is an ordering tuple, which in code can just be a list. In my code, I just use a wrapper around vec<uint>, with a distance. Here is a simpler version of it:

class Solution {
    public:
        vector<uint> assigned;
        flt distsq;
    [ ... ]
}

The distance (or distance squared, distsq) is simply from above.

Below, I’m going to use the form ([2, 0, 1], 10.3), where that corresponds to , with a total distance squared .

Step 2: Djikstra’s Algorithm

If you look at the form of the solution, you’ll notice that it has a lot in common with a path-finding algorithm: we want a list of “nodes” to visit, and the more nodes we visit, the greater our distance — is monotonically increasing. This isn’t quite a path-finding algorithm, because we want to visit all nodes (i.e. assign every particle in to a particle in ), so its close to the Traveling Salesman problem, but we can still use a similar solution.

We start with a Queue with one element in it:

([], 0.0)

Then we take out the first element, and “expand” it to try assigning every possible particle in to the particle 0 in , and put those back in our Queue, keeping it sorted by distance squared:

([3], 1.0),
([1], 9.0),
([2], 16.0),
([0], 25.0)

Now we repeat the process with the first element in our Queue:

([3, 2], 2.0),
([1], 9.0),
([2], 16.0),
([3, 0], 17.0),
([0], 25.0),
([3, 1], 29.0)

If we keep doing this, we might end up something like this:

([3, 2, 0, 1], 3.0),
([1], 9.0),
([2], 16.0),
([3, 0], 17.0),
([0], 25.0),
([3, 1], 29.0),
([3, 2, 1], 38.0)

And now we know the best pairing is , with . Conveniently enough, the remainder of our Queue is basically a proof that any other pairing would have a greater — for example, any pairing including has , so we don’t have to investigate all 6 ways that could be the case.

Extra Difficulties

Translations

Now let’s imagine our particles might not only be “unpaired”, but possibly also shifted – so now we want to find not just the best pairing, but the best (pairing, translation) to minimize : We now want to find the ordering tuple and vector that minimizes

It turns out, this is simple: you subtract off the “center of mass” difference:

So we use the algorithm above, but add in that calculation first before calculating .

Periodic Boundary Conditions

In my case, we often want to deal with particles in a box with periodic boundary conditions, i.e. a box where a particle that disappears on one side appears on the other. If we still have to calculate the minimal (pairing, translation), we now have a problem: you can’t use the center of mass trick above, because there is no well-defined center-of-mass like that with periodic boundary conditions.

To make a long story short, you can instead calculate in the following way, that is independent of translations:

where means “shortest distance given periodic boundary conditions in a box of shape ”. In practice, this is defined like so:

for(i=0; i<NUMBER_OF_DIMENSIONS; ++i){
    float dxi = r[i] - s[i];
    ominused[i] = dxi - remainder(dxi, L[i]);
}

Proving that the defined here for the periodic boundary conditions is equivalent to the previous one is well outside the scope of this post, but for anyone interested, see here, or the original LyX file.

Rotations

Now imagine that however you are generating these sets, sometimes they appear as rotated versions of each other, and possibly even a mirror.

This is also fairly simple to handle:

I have code that handles this here.