Distancia de Manhattan

Definición

Para puntos $p$ y $q$ en un plano, podemos definir la distancia entre ellos como la suma de las diferencias entre sus coordenadas $x$ e $y$ :

$d(p,q) = |x_p - x_q| + |y_p - y_q|$

Definida de esta manera, la distancia corresponde a la llamada geometría de Manhattan (taxicab) , en la que los puntos se consideran intersecciones en una ciudad bien diseñada, como Manhattan, donde solo se puede mover por las calles en horizontal o en vertical, como se muestra en la imagen de abajo:

Estas imágenes muestran algunos de los caminos más cortos de un punto negro al otro, todos ellos de longitud $12$ .

Hay algunos trucos y algoritmos interesantes que se pueden hacer con esta distancia, y mostraremos algunos de ellos aquí.

Par de puntos más lejanos en distancia de Manhattan

Dados $n$ puntos $P$ , queremos encontrar el par de puntos $p,q$ que están más alejados, es decir, maximizar $|x_p - x_q| + |y_p - y_q|$ .

Pensemos primero en una dimensión, de modo que $y=0$ . La observación principal es que podemos hacer fuerza bruta sobre si $|x_p - x_q|$ es igual a $x_p - x_q$ o a $-x_p + x_q$ , porque si “erramos el signo” del valor absoluto, solo obtendremos un valor más pequeño, de modo que no puede afectar la respuesta. Más formalmente, se cumple que:

$|x_p - x_q| = \max(x_p - x_q, -x_p + x_q)$

Así, por ejemplo, podemos intentar que $p$ sea tal que $x_p$ tenga el signo más, y entonces $q$ debe tener el signo menos. De este modo queremos encontrar:

$\max\limits_{p, q \in P}(x_p + (-x_q)) = \max\limits_{p \in P}(x_p) + \max\limits_{q \in P}( - x_q ).$

Nótese que podemos extender esta idea más allá a 2 (¡o más!) dimensiones. Para $d$ dimensiones, debemos hacer fuerza bruta sobre $2^d$ valores posibles de los signos. Por ejemplo, si estamos en $2$ dimensiones y forzamos que $p$ tenga ambos signos más, queremos encontrar:

$\max\limits_{p, q \in P} [(x_p + (-x_q)) + (y_p + (-y_q))] = \max\limits_{p \in P}(x_p + y_p) + \max\limits_{q \in P}(-x_q - y_q).$

Como hicimos independientes a $p$ y $q$ , ahora es fácil encontrar el $p$ y el $q$ que maximizan la expresión.

El código de abajo generaliza esto a $d$ dimensiones y corre en $O(n \cdot 2^d \cdot d)$ .


long long ans = 0;
for (int msk = 0; msk < (1 << d); msk++) {
    long long mx = LLONG_MIN, mn = LLONG_MAX;
    for (int i = 0; i < n; i++) {
        long long cur = 0;
        for (int j = 0; j < d; j++) {
            if (msk & (1 << j)) cur += p[i][j];
            else cur -= p[i][j];
        }
        mx = max(mx, cur);
        mn = min(mn, cur);
    }
    ans = max(ans, mx - mn);
}

Rotación de los puntos y distancia de Chebyshev

Es bien sabido que, para todo $m, n \in \mathbb{R}$ ,

$|m| + |n| = \text{max}(|m + n|, |m - n|).$

Para demostrarlo, solo hay que analizar los signos de $m$ y $n$ . Y se deja como ejercicio.

Podemos aplicar esta ecuación a la fórmula de la distancia de Manhattan para descubrir que

$d((x_1, y_1), (x_2, y_2)) = |x_1 - x_2| + |y_1 - y_2| = \text{max}(|(x_1 + y_1) - (x_2 + y_2)|, |(y_1 - x_1) - (y_2 - x_2)|).$

La última expresión de la ecuación anterior es la distancia de Chebyshev de los puntos $(x_1 + y_1, y_1 - x_1)$ y $(x_2 + y_2, y_2 - x_2)$ . Esto significa que, después de aplicar la transformación

$\alpha : (x, y) \to (x + y, y - x),$

la distancia de Manhattan entre los puntos $p$ y $q$ se convierte en la distancia de Chebyshev entre $\alpha(p)$ y $\alpha(q)$ .

Además, podemos darnos cuenta de que $\alpha$ es una semejanza espiral (rotación del plano seguida de una dilatación respecto de un centro $O$ ) con centro $(0, 0)$ , ángulo de rotación de $45^{\circ}$ en sentido horario y dilatación por $\sqrt{2}$ .

Aquí hay una imagen para ayudar a visualizar la transformación:

Árbol de expansión mínima de Manhattan

El problema del MST de Manhattan consiste en, dados algunos puntos en el plano, encontrar las aristas que conectan todos los puntos y tienen una suma total mínima de pesos. El peso de una arista que conecta dos puntos es su distancia de Manhattan. Por simplicidad, asumimos que todos los puntos tienen ubicaciones distintas. Aquí mostramos una forma de encontrar el MST en $O(n \log{n})$ hallando para cada punto su vecino más cercano en cada octante, como se representa en la imagen de abajo. Esto nos dará $O(n)$ aristas candidatas, que, como mostramos más abajo, garantizarán que contienen el MST. El paso final es entonces usar algún MST estándar, por ejemplo, el algoritmo de Kruskal usando conjuntos disjuntos .

*Los 8 octantes relativos a un punto S*

El algoritmo mostrado aquí se presentó por primera vez en un artículo de H. Zhou, N. Shenoy y W. Nichollos (2002) . También hay otro algoritmo conocido que usa un enfoque de divide y vencerás de J. Stolfi , que también es muy interesante y solo difiere en la forma en que encuentran el vecino más cercano en cada octante. Ambos tienen la misma complejidad, pero el que se presenta aquí es más fácil de implementar y tiene un factor constante menor.

Primero, entendamos por qué basta considerar solo el vecino más cercano en cada octante. La idea es mostrar que para un punto $s$ y cualesquiera otros dos puntos $p$ y $q$ en el mismo octante, $d(p, q) < \max(d(s, p), d(s, q))$ . Esto es importante, porque muestra que si hubiera un MST donde $s$ está conectado tanto a $p$ como a $q$ , podríamos borrar una de estas aristas y añadir la arista $(p,q)$ , lo que disminuiría el costo total. Para demostrarlo, asumimos sin pérdida de generalidad que $p$ y $q$ están en el octante $R_1$ , que se define por: $x_s \leq x$ y $x_s - y_s > x - y$ , y luego hacemos un análisis por casos. La imagen de abajo da algo de intuición sobre por qué esto es cierto.

*Intuitivamente, la limitación del octante hace imposible que

p

q

estén ambos más cerca de

s

que el uno del otro*

Por tanto, la pregunta principal es cómo encontrar el vecino más cercano en cada octante para cada uno de los $n$ puntos.

Vecino más cercano en cada octante en O(n log n)

Por simplicidad nos concentramos en el octante NNE ( $R_1$ en la imagen de arriba). Todas las demás direcciones se pueden encontrar con el mismo algoritmo rotando la entrada.

Usaremos un enfoque de línea de barrido. Procesamos los puntos de suroeste a noreste, es decir, por $x + y$ no decreciente. También mantenemos un conjunto de puntos que aún no tienen su vecino más cercano, al que llamamos “conjunto activo”. Añadimos las imágenes de abajo para ayudar a visualizar el algoritmo.

*En negro con una flecha se puede ver la dirección de la línea de barrido. Todos los puntos por debajo de esta línea están en el conjunto activo, y los puntos de arriba aún no están procesados. En verde vemos los puntos que están en el octante del punto procesado. En rojo los puntos que no están en el octante buscado.*

*En esta imagen vemos el conjunto activo después de procesar el punto

p

. Nótese que los

2

puntos verdes de la imagen anterior tenían a

p

en su octante norte-noreste y ya no están en el conjunto activo, porque ya encontraron su vecino más cercano.*

Cuando añadimos un punto nuevo $p$ , para cada punto $s$ que lo tiene en su octante podemos asignar de forma segura $p$ como el vecino más cercano. Esto es cierto porque su distancia es $d(p,s) = |x_p - x_s| + |y_p - y_s| = (x_p + y_p) - (x_s + y_s)$ , porque $p$ está en el octante norte-noreste. Como todos los puntos siguientes no tendrán un valor menor de $x + y$ por el paso de ordenamiento, está garantizado que $p$ tiene la menor distancia. Luego podemos quitar todos esos puntos del conjunto activo, y finalmente añadir $p$ al conjunto activo.

La siguiente pregunta es cómo encontrar de forma eficiente qué puntos $s$ tienen a $p$ en el octante norte-noreste. Es decir, qué puntos $s$ satisfacen:

$x_s \leq x_p$
$x_p - y_p < x_s - y_s$

Como ningún punto del conjunto activo está en la región $R_1$ de otro, también tenemos que para dos puntos $q_1$ y $q_2$ del conjunto activo, $x_{q_1} \neq x_{q_2}$ y su orden implica $x_{q_1} < x_{q_2} \implies x_{q_1} - y_{q_1} \leq x_{q_2} - y_{q_2}$ .

Se puede intentar visualizar esto en las imágenes de arriba pensando en el orden de $x - y$ como una “línea de barrido” que va del noroeste al sureste, de modo que es perpendicular a la que está dibujada.

Esto significa que si mantenemos el conjunto activo ordenado por $x$ los candidatos $s$ están colocados de forma consecutiva. Entonces podemos encontrar el mayor $x_s \leq x_p$ y procesar los puntos en orden decreciente de $x$ hasta que se rompa la segunda condición $x_p - y_p < x_s - y_s$ (en realidad podemos permitir que $x_p - y_p = x_s - y_s$ y eso trata el caso de puntos con coordenadas iguales). Nótese que, como quitamos del conjunto justo después de procesar, esto tendrá una complejidad amortizada de $O(n \log(n))$ . Ahora que tenemos el punto más cercano en la dirección noreste, rotamos los puntos y repetimos. Es posible mostrar que de esta forma también encontramos el punto más cercano en la dirección suroeste, de modo que podemos repetir solo 4 veces, en lugar de 8.

En resumen:

Ordenamos los puntos por $x + y$ en orden no decreciente;
Para cada punto, iteramos sobre el conjunto activo empezando por el punto con el mayor $x$ tal que $x \leq x_p$ , y rompemos el bucle si $x_p - y_p \geq x_s - y_s$ . Para cada punto válido $s$ añadimos la arista $(s,p, d(s,p))$ a nuestra lista;
Añadimos el punto $p$ al conjunto activo;
Rotamos los puntos y repetimos hasta iterar sobre todos los octantes.
Aplicamos el algoritmo de Kruskal sobre la lista de aristas para obtener el MST.

Abajo se puede encontrar una implementación, basada en la de KACTL .


struct point {
    long long x, y;
};

// Devuelve una lista de aristas en el formato (weight, u, v). 
// Pasar esta lista al algoritmo de Kruskal da el MST de Manhattan.
vector<tuple<long long, int, int>> manhattan_mst_edges(vector<point> ps) {
    vector<int> ids(ps.size());
    iota(ids.begin(), ids.end(), 0);
    vector<tuple<long long, int, int>> edges;
    for (int rot = 0; rot < 4; rot++) { // para cada rotación
        sort(ids.begin(), ids.end(), [&](int i, int j){
            return (ps[i].x + ps[i].y) < (ps[j].x + ps[j].y);
        });
        map<int, int, greater<int>> active; // (xs, id)
        for (auto i : ids) {
            for (auto it = active.lower_bound(ps[i].x); it != active.end();
            active.erase(it++)) {
                int j = it->second;
                if (ps[i].x - ps[i].y > ps[j].x - ps[j].y) break;
                assert(ps[i].x >= ps[j].x && ps[i].y >= ps[j].y);
                edges.push_back({(ps[i].x - ps[j].x) + (ps[i].y - ps[j].y), i, j});
            }
            active[ps[i].x] = i;
        }
        for (auto &p : ps) { // rotar
            if (rot & 1) p.x *= -1;
            else swap(p.x, p.y);
        }
    }
    return edges;
}