Na wypadek awarii

W ramach wakacyjnego relaksu proponuję zadanie z zawodów Central European Programming Contest organizowanych osiem lat temu przez Instytut Informatyki Uniwersytetu Wrocławskiego. Zadanie jest oczywiście z geometrii, ale tym razem w tylko dwóch wymiarach: dostajemy n \leq 10^5 punktów na płaszczyźnie, dla każdego z nich chcemy znaleźć najbliższego sąsiada. Odległość jest standardowa, to znaczy euklidesowa.

Oryginalną treść można zobaczyć na przykład tutaj, można tam także sprawdzić swoje rozwiązanie. Jedynym szczegółem, na który być może warto zwrócić uwagę, jest obietnica, że podane punkty nie powtarzają się (oraz n\geq 2, więc zadanie jest dobrze zdefiniowane).

Każdy choćby średnio doświadczony zawodnik po przeczytaniu treści pewnie od razu pomyślał o klasycznej wersji, w której szukamy pary najbliższych punktów. Jeśli potrafilibyśmy szybko znaleźć najbliższego sąsiada dla każdego punktu to oczywiście łatwo moglibyśmy wyznaczyć taką parę, warto więc przypomnieć sobie jak rozwiązuję się tę prostszą (i bardziej znaną) wersję w czasie \mathcal{O}(n\log n).

Są przynajmniej dwa sposoby. Obydwa korzystają z podobnej obserwacji, ale jeden stosuje strategię dziel-i-zwyciężaj, a drugi zamiatanie. Ten drugi jest pewnie nieco łatwiejszy w implementacji, jednak dla nas bardziej przydatny będzie ten pierwszy.

Podzielmy płaszczyznę na dwie części prostą pionową x=t tak, aby w każdej z otrzymanych dwóch części znalazła się mniej więcej połowa wszystkich punktów. (Załóżmy na razie, że wszystkie punkty mają różne współrzędne x i taki podział jest możliwy.) Rekurencyjnie znajdźmy parę najbliższych punktów na lewo od prostej oraz parę najbliższych punktów na prawo, niech \delta oznacza minimum ze znalezionych w ten sposób odległości. Być może \delta jest już poprawną odpowiedzią dla całego problemu, jednak może się także zdarzyć, że ta odpowiedź jest wyznaczona przez punkty leżące po różnych stronach prostej. Zauważmy jednak, że w takim przypadku jeden z tych punktów jest postaci (x_1,y_1) gdzie x_1 \in (t-\delta,t), a drugi (x_2,y_2) dla x_2 \in (t,t+\delta). Innymi słowy: możemy spokojnie usunąć wszystkie punkty leżące poza pasem szerokości 2\delta. Co dalej?

Ustalmy (x_1,y_1). Zauważmy, że wszystkie punkty (x_2,y_2), których wybranie mogłoby prowadzić do polepszenia odpowiedzi, spełniają nie tylko x_2 \in (t,t+\delta) lecz również y_2 \in (y_1-\delta,y_1+\delta). Ile może być takich punktów?

wwypadku

Nie za dużo: co najwyżej jeden w każdym z kwadratów \frac{\delta}{2} \times \frac{\delta}{2} na powyższym rysunku. Gdyby bowiem dwa z nich leżały w tym samym kwadracie to byłyby w odległości nie większej niż \frac{\delta}{2}\sqrt{2} <  \delta. Fajnie, pokazaliśmy więc, że dla każdego punktu po lewej mamy tylko kilku sensownych kandydatów na punkt po prawej, który należy rozważyć jako być może poprawiający odpowiedź. Ale jak szybko dobrać się do tych kilku kandydatów?

Posortujmy wszystkie punkty z lewej (x_1,y_1), dla których x_1 > t-\delta, względem współrzędnej y. Podobnie, posortujmy wszystkie punkty z prawej (x_2,y_2), dla których x_2 < t+\delta, względem współrzędnej y. Teraz zauważmy, że dla każdego punktu z lewej wystarczy znaleźć pierwszy punkt z prawej na naszej posortowanej liście, dla którego y_2 \geq y_1: wszyscy kandydaci mogą być wtedy wygenerowani biorąc jego kilku poprzedników i kilku następników. Jak szybko znaleźć ten pierwszy punkt z prawej? Wystarczy przesuwać dwa palce po obydwu posortowanych listach, co wymaga w sumie tylko czasu liniowego. Ignorując więc konieczność posortowania punktów, czas działania całego rozwiązania jest opisany przez rekurencję T(n)=2T(n/2)+\mathcal{O}(n), czyli \mathcal{O}(n\log n). Niestety, nie do końca możemy zignorować tę konieczność. Na szczęście dość łatwo jednak uniknąć zwiększenie złożoności do \mathcal{O}(n\log^2 n) stosując prosty trik: sortujemy wszystkie punkty raz a dobrze względem współrzędnej y. Następnie w każdym wywołaniu rekurencyjnym zakładamy, że dostaliśmy tak posortowane punkty. Łatwo widać, że możemy utrzymywać taki porządek dzieląc punkty na lewe i prawe. Jest jeszcze jeden szczegół: chcemy wybrać prostą pionową dzielącą punkty na dwie mniej więcej równe grupy. W tym celu można pewnie wybrać medianę wszystkich współrzędnych x, co jest trywialne jeśli są one posortowane względem współrzędnej x. Tak naprawdę przekazujemy więc dwie listy, obydwie zawierające ten sam zbiór punktów: jedna jest posortowana względem współrzędnej x, a druga względem współrzędnej y.

Skoro przypomnieliśmy już sobie jak rozwiązać prostszą wersję zadania, wróćmy do tej trudniejszej. Naturalne wydaje się zastosowanie podobnej metody, to znaczy wybranie prostej pionowej x=t dzielącej punkty na dwie mniej więcej równe grupy i uruchomienie się rekurencyjnie dla wszystkich punktów po lewej i wszystkich punktów po prawej. Teraz jednak wynikiem takiego wywołania nie jest jedna liczba: dla każdego punktu p=(x,y) dostajemy jego odległość do najbliższego sąsiada \delta_p. Być może wywołanie rekurencyjne wystarczy, żeby poprawnie wyznaczyć \delta_p, może być jednak i tak, że najbliższy sąsiad punktu po lewej jest punktem po prawej (lub na odwrót). Sytuacja jest symetryczna, skupmy się więc na „naprawieniu” sytuacji dla każdego punktu p=(x_1,y_1) po lewej. Chcemy więc sprawdzić, czy przypadkiem nie ma punktu po prawej p'=(x_2,y_2), dla którego \delta_p > d(p,p'). Ale jak?

Kuszące wydaje się ograniczenie liczby kandydatów, których należy rozważyć dla każdego takiego punktu z lewej. Zauważmy więc, że (x_2,y_2) jest sensownym kandydatem tylko jeśli przecięcie okręgu o środku w (x_1,y_1) i promieniu \delta_p z prostą x=t zawiera punkt (t,y_2) tak jak na poniższej ilustracji.

wwypadku1

Czy potrafilibyśmy ograniczyć liczbę takich kandydatów? Wydaje się, że nie bardzo: może się zdarzyć, że jeden punkt z lewej będzie miał ich wielu. Ale może dałoby się ograniczyć ich sumaryczną liczbę? Lub, mówiąc inaczej, może dałoby się ograniczyć liczbę punktów z lewej, dla których dany punkt z prawej jest sensownym kandydatem? Spróbujmy!

Załóżmy, że punkt z prawej strony (x_2,y_2) jest sensownym kandydatem dla 4 punktów z lewej strony. Podzielmy lewą część płaszczyzny na 3 kawałki tak jak poniżej.

wwypadku21

W jednej z tych części mamy 2 punkty z lewej strony, oznaczmy je p=(x_1,y_1) oraz p'=(x_1',y_1'). Zauważmy, że kąt wyznaczony przez p,(x_2,y_2),p' jest nie większy niż 60 stopni, oznaczmy go przez \alpha. Tak naprawdę chcielibyśmy założyć, że \alpha < 60^{\circ}: wymaga to zdefiniowania kawałków tak, aby były jednostronnie domknięte i założenia, żew wszystkie punkty po lewej są ściśle na lewo od prostej x=t. Mamy więc trójkąt o dwóch bokach długości \delta_p,\delta_{p'} i kącie \alpha < 60^{\circ} między nimi. Długość trzeciego boku to po prostu odległość między p oraz p', która nie może być mniejsza niż \max(\delta_p,\delta_{p'}) (bo przecież w przeciwnym wypadku wyznaczone odległości \delta_p,\delta_{p'} nie byłyby poprawne!). Korzystając z prawa kosinusów i zakładając bez straty ogólności, że \delta_p \leq \delta_{p'} dostajemy więc, że \delta^2_{p'} \leq \delta^2_p+\delta^2_{p'}-2\delta_p\delta_{p'}\cos \alpha, więc po uproszczeniu \cos \alpha \leq \frac{1}{2}. To jest jednak sprzeczność z założeniem, że \alpha \leq 60^{\circ}.

No dobrze, czyli pokazaliśmy, że sumaryczna liczba kandydatów jest liniowa. Ale, znów, jak ich wygenerować? Prawie tak samo jak poprzednio: sortujemy punkty po lewej i punkty po prawej względem ich współrzędnych y. Dla każdego punktu po lewej p=(x_1,y_1) przeglądamy wszystkie punkty po prawej, których współrzędne y mieszczą się w odpowiednim zakresie. Zauważmy, że y_1 należy do tego zakresu, więc znając punkt po prawej z najmniejszą współrzędną y_2 \geq y_1 jesteśmy w stanie wygenerować wszystkie takie punkty w czasie proporcjonalnym do ich licby. Tak jak w prostszej wersji zadania, możemy uniknąć sortowania utrzymując dwie kopie listy punktów i otrzymać rozwiązanie, które działa w czasie \mathcal{O}(n\log n).

Wypada jeszcze wrócić do założenia, że wszystkie punkty mają różne współrzędne x. Jeśli tak nie jest, nie zawsze istnieje prosta pionowa dzieląca je na dwie z grubsza równe grupy. Prosta nie musi być jednak pionowa: można ją lekko "kopnąć". Tak naprawdę oznacza to, że można posortować punkty leksykograficznie i wywołać się rekurencyjnie dla pierwszej i drugiej połówki tak posortowanej listy.

Warto także zastanowić się, czy nasze rozwiązanie jest optymalne. Nietrudno pokazać, że i owszem, jeśli dopuszczamy tylko podstawowe operacje arytmetyczne dla współrzędnych (a konkretniej: testy liniowe). Jeśli jednak dopuścimy także używanie podłogi (zaokrąglania) to można skonstruować algorytm, który działa w czasie \mathcal{O}(n\log\log n). Jeśli zaś dopuścimy także randomizację: istnieje algorytm liniowy! Czy którekolwiek z tych szybszych rozwiązań uogólnia się do naszej wersji zadania?

4 myśli nt. „Na wypadek awarii

  1. „Chcemy więc sprawdzić, czy przypadkiem nie ma punktu po prawej p’=(x_2,y_2), dla którego \delta_p < d(p,p')"

    Chyba nierówność w drugą stronę?

  2. Hi, thanks for the post.
    Could you please clarify how to generate all the candidates in linear time?
    For each point p on the left, I cannot filter all on the right using it delta[p]. This would spam a quadratic algorithm, right?
    I am probably missing something and the translation isn’t helping.
    :/

    • For every point (x_p,y_p) on the left, we want to extract all points (x,y) on the right, such that y belongs to the interval [y_p – delta_p, y_p + delta_p]. (This is really the definition of a candidate, and we have proved that there wouldn’t be too many of them.) So now it is enough to sort all points on the right with respect to their y coordinates. Also sort points on the left with respect to their value of y_p – delta_p. Then scan both sorted lists at the same time, so that while processing the point (x_p,y_p) we know the point on the right with the smallest y coordinate, such that y >= y_p – delta_p. If y > y_p + delta_p, (x_p,y_p) does not generate any candidates. Otherwise, we start extracting the candidates from the sorted list of point on the right, starting from (x,y) (who is a candidate).

      I hope this helps!

Dodaj komentarz