cheap way to compute a pretty good traversal to visit 2d points

traversal.js

// public domain or Unlicense or MIT licensed

// operates in place
function sortPoints(points) {
    let bb = bbox(points);
    let bands = numberOfBands(points, bb);
    sortYbandsX(points, bb, bands);
}

function sortYbandsX(points, bbox, bands) {
    let [_1, _2, ymin, ymax] = bbox;
    let yspan = ymax - ymin + 1;
    let band = (p) => Math.floor((p.y - ymin) / yspan * bands);
    let cmp = (a, b) => {
        let band_a = band(a);
        let band_b = band(b);
        if (band_a < band_b) return -1;
        if (band_a > band_b) return 1;
        // bands are equal, sort ltr on even bands and rtl on odd bands
        let sign = band_a % 2 === 0 ? 1 : -1;
        if (a.x < b.x) return sign * -1;
        if (a.x > b.x) return sign * 1;
        // if x's are equal, sort by y
        return a.y - b.y;

    };
    points.sort(cmp);
    return points;
}

function numberOfBands(points, bbox) {
    const maxBands = 20;
    let [xmin, xmax, ymin, ymax] = bbox;
    let bestScore = Infinity;
    let P = points.length;
    let xspan = xmax - xmin;
    let yspan = ymax - ymin;
    for (let b = 1; b < maxBands; b++) {
        let p = P / b;  // avg number of points per band
        let h = yspan / b;  // avg height of band
        let w = xspan / p;  // avg x spacing between points
        // each move will avg h/2 in vertical and w in width
        // can skip the sqrt
        let c = Math.pow(h/2, 2) + Math.pow(w, 2);
        if (c > bestScore) return b - 1;
        bestScore = c;
    }
    return maxBands;
}

function bbox(points) {
    return points.reduce(
        ([xmin, xmax, ymin, ymax], p) => [
            Math.min(xmin, p.x),
            Math.max(xmax, p.x),
            Math.min(ymin, p.y),
            Math.max(ymax, p.y),
        ],
        [Infinity, -Infinity, Infinity, -Infinity]
    );
}

traversal.py

# public domain or Unlicense or MIT licensed

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import time
import json
import subprocess

def plot(pointss, strategies, outname='/tmp/tplot.png'):
    fig, axs = plt.subplots(len(strategies), len(pointss), figsize=(6 * len(pointss), 6 * len(strategies)), squeeze=False)
    for i, strategy in enumerate(strategies):
        for j, points in enumerate(pointss):
            t0 = time.time()
            xy = strategy(points)
            t1 = time.time()
            axs[i][j].plot(xy[:, 0], xy[:, 1], '.')
            axs[i][j].plot(xy[:, 0], xy[:, 1], alpha=0.5)
            dist = np.linalg.norm(xy[:-1] - xy[1:], axis=1).sum()
            axs[i][j].set_title('{}, dist: {:.2f}, elapsed: {:.2f}ms'.format(strategy.__name__, dist, (t1 - t0) * 1000))
            axs[i][j].set_aspect('equal')

    plt.tight_layout()
    plt.savefig(outname)
    plt.close()

def strat_noop(points):
    return points

def strat_random(points):
    I = np.arange(len(points))
    np.random.shuffle(I)
    return points[I]

def strat_sort_x(points):
    I = np.argsort(points[:, 0])
    return points[I]

def strat_sort_y_band_x(points, bands=10):
    ymin = points[:, 1].min()
    ymax = points[:, 1].max()
    yspan = ymax - ymin
    yspan += 1
    def key(xy):
        x, y = xy
        band = int(np.floor((y - ymin) / yspan * bands))
        if band % 2 == 1:
            x *= -1
        return band, x, y

    return np.array(sorted(points, key=key))

def strat_sort_y_band_x_bands_auto(points):
    xmin = points[:, 0].min()
    xmax = points[:, 0].max()
    ymin = points[:, 1].min()
    ymax = points[:, 1].max()
    yspan = ymax - ymin
    xspan = xmax - xmin
    P = len(points)
    bs = np.arange(1, 20)
    p = P / bs
    l = yspan / bs
    move_per_p = (l/2) ** 2 + (xspan/p) ** 2
    c = P * move_per_p

    bands = bs[np.argmin(c)]
    if bands > 1:
        xy = strat_sort_y_band_x(points, bands=bands-1)
        dist = np.linalg.norm(xy[:-1] - xy[1:], axis=1).sum()
        print('bands-1', dist, 'pred', c[bands-1-1])

    if bands < bs[-1]:
        xy = strat_sort_y_band_x(points, bands=bands+1)
        dist = np.linalg.norm(xy[:-1] - xy[1:], axis=1).sum()
        print('bands+1', dist, 'pred', c[bands+1-1])

    xy = strat_sort_y_band_x(points, bands=bands)
    dist = np.linalg.norm(xy[:-1] - xy[1:], axis=1).sum()
    print('bands', bands, dist, 'pred', c[bands-1])

    return xy

def strat_tsp_christofides(points):
    G = nx.Graph()
    for i in range(len(points)):
        for j in range(i+1, len(points)):
            G.add_edge(i, j, weight=np.linalg.norm(points[i] - points[j]))
    from networkx.algorithms.approximation import christofides
    nodes = christofides(G)
    nodes.pop()
    return points[nodes]

def gen_pad_points(centers=50, npins=(1, 4), pin_spacing=1, board=(50, 50)):
    ret = []
    centers = np.random.random(size=(centers, 2))
    centers[:, 0] *= board[0]
    centers[:, 1] *= board[1]
    for center in centers:
        pins = np.random.randint(*npins)
        offsets = np.zeros((pins, 2))
        axis = int(np.random.random() < 0.5)
        offsets[:, axis] = (np.arange(pins) - pins//2) * pin_spacing
        points = center + offsets
        np.random.shuffle(points)
        ret.extend(points)

    return np.array(ret)

def js(points):
    with open('traversal.js') as fh:
        code = fh.read()
    code += '''
let points = JSON.parse(require('fs').readFileSync('/dev/stdin').toString());
sortPoints(points);
console.log(JSON.stringify(points));
    '''
    stdin = json.dumps([{'x': x, 'y': y} for x, y in points])
    ret = subprocess.run(['node', '-e', code], input=stdin, text=True, capture_output=True)
    if ret.returncode != 0:
        print(ret.stdout)
        print(ret.stderr)
        raise ValueError('fail')
    if ret.stderr:
        print(ret.stderr)
    j = json.loads(ret.stdout)
    return np.array([[d['x'], d['y']] for d in j])

np.random.seed(42)

points = [
    # np.array([[0, 0], [0, 10], [10, 0], [-10, 0], [0, -10]]),
    # np.random.random(size=(30, 2)),
    # np.random.random(size=(40, 2)) * 10 - 5,
    # np.round(np.random.random(size=(40, 2)) * 10, 0),
    gen_pad_points(),
    gen_pad_points(board=(100, 50)),
    gen_pad_points(board=(50, 100)),
    ]

strats = [
    strat_noop,
    strat_random,
    strat_sort_x,
    strat_sort_y_band_x_bands_auto,
    strat_tsp_christofides,
    js,
    ]

plot(points, strats, outname='/tmp/tplot.svg')