Поиск кратчайшего расстояния между точками в трёхмерном пространстве

Ставший уже традиционным ежемесячный конкурс по функциональному программированию в июне сего года собрал пока наибольшее количество участников с момента запуска оных конкурсов. Не знаю, что послужило причиной этого, но факт остаётся фактом — в конкурсе приняло участие 35 человек, представивших 36 решений, из которых большиинство было на языке Haskell (7 решений) и C++ (6 решений). Другими использованными языками программирования были: Assembler, C, C#, Clojure, D2, Erlang, F#, Java, Nemerle, OCaml, Perl, Python и Scala. И это очень приятно, что имеет место такая диверсификация языков и подходов.

Ну а в качестве задачи была предложена такая формулировка одной из традиционных задач для олимпиад по программированию:

Демиург задумал поместить новую разумную расу на планету во Вселенной так, чтобы она как можно быстрее освоила весь доступный космос. Для этих целей он собрал перечень всех планет, находящихся в поясе Златовласки, после чего решил найти среди них две, наиболее близко расположенные друг к другу. Разумная раса должна быть расположена на одной из таких планет. Это, по замыслу Демиурга, позволит максимально быстро обеспечить освоение космоса. Помогите Демиургу найти две планеты, находящиеся на наименьшем расстоянии друг от друга.

Как видно, это обычная задача на поиск минимального расстояния между точками в пространстве, однако сложность была увеличена тем, что пространство внезапно стало трёхмерным, в то время как в традиционном условии пространство двумерно. Но это нисколько не испугало конкурсантов, хотя и внесло в исходные коды для решения некоторые изменения известного алгоритма.

Далее в этой статье рассматривается несколько вопросов — генерация исходных данных, решение одним способом, решение другим способом, а также код на языке Haskell, завоевавший приз зрительских симпатий.

Генерация миллиона точек

Для генерации миллиона точек можно было бы воспользоваться простейшими функциями из стандартных модулей в поставке GHC, однако решено было написать один модуль, который использовался бы как при генерации, так и при поиске решения. Само собой разумеется, что речь идёт о модуле с описанием программных сущностей для представления точки в трёхмерном пространстве, а также всего остального, что надо будет для работы с точками. Определим этот модуль очень правильно, то есть со всеми необходимыми декларациями интерфейса.

Начнём:

module Point
(
  Point(..),
  readPoint,
  loadFromFile,
  distance,
  prettyPrintPoint
)
where

Здесь мы определяем интерфейс модуля, то есть все те программные сущности, которые будут видны в тех модулях, которые импортируют модуль Point. Для одноимённого алгебраического типа данных указано, что он экспортирует все имеющиеся в нём конструкторы. Ну и перечислены все утилитарные функции для работы с точками: чтение координат точки из строки, загрузка списка точек из файла, вычисление расстояния между двумя точками и вывод описания точки в «приятном виде» в строку. А где же простое преобразование описания точки в строку — может спросить кто-то? А оно будет определено в виде экземпляра класса Show, а, как известно, все экземпляры неявно экспортируются и импортируются.

import System.IO (Handle, hGetContents)

Для работы с файлами нам придётся импортировать пару программных сущностей из стандартного модуля для работы с системой ввода/вывода.

Теперь займёмся определением самого главного АТД, вокруг которого всё будет вертеться. Поскольку язык Haskell предлагает программистам истинный полиморфизм, нет никакого резона делать типы данных конкретными. Поэтому мы определим полиморфный тип для представления точки в трёхмерном пространстве. Параметризоваться, само собой разумеется, будет тип координаты. Делается это так:

data Point a = Point
               {
                 getX :: a,
                 getY :: a,
                 getZ :: a
               }

Заодно сразу определяются три функции для доступа (как по чтению, так и для записи) к компанентам типа — getX, getY и getZ. Конечно, префикс get здесь не очень годен, поскольку эти идентификаторы также являются «сеттерами», но использовать однобуквенные наименования полей в данном случае было бы слишком неосмотрительно — однобуквенные идентификаторы часто встречаются, а потому возможны коллизии имён, что, в свою очередь, может привести к нехорошим логическим ошибкам.

Теперь же определим экземпляр класса Show для этого типа. Этот экземпляр потребуется нам для записи списка координат в файл. Вот простейшее определение:

instance Show a => Show (Point a) where
  show (Point x y z) = show x ++ "t" ++ show y ++ "t" ++ show z

Как видно, мы просто выводим в строку три координаты, тип которых, кстати, также должны иметь экземпляр класса Show, разделяя их символом табуляции. Возможно, что лучшим определением функции show было бы show (Point x y z) = intercalate "t" $ map show [x, y, z] (или даже ещё мудрёнее — show p = intercalate "t" $ map (show . ($ p)) [getX, getY, getZ]), но насколько такое определение будет более понятным начинающим программистам на языке Haskell :)?

Ну и далее определения указанных ранее в интерфейсе модуля четырёх функций:

readPoint :: Read a => String -> Point a
readPoint = ((x:y:z:_) -> Point x y z) . map read . words

Данная функция преобразует заданную строку в описание точки типа Point. Строка, само собой разумеется, должна иметь формат «XXX YYY ZZZ», то есть три числа, разделённые пробельными символами (не обязательно табуляцией, можно и пробелами). Какие-либо знаки препинания не нужны, а числа XXX и т. д. должны иметь возможность преобразовываться из строки в значения при помощи соответствующих экземпляров класса Read. В противном случае произойдёт ошибка времени исполнения.

loadFromFile :: Read a => Handle -> IO [Point a]
loadFromFile h = do cnt <- hGetContents h
                    length cnt `seq` return (map readPoint $ lines cnt)

Как видно, эта функция получает на вход идентификатор открытого файла и возвращает список точек, при этом координаты точек должны иметь возможность быть преобразованными из строки. Функция предназначена для чтения списка точек из внешнего файла, причём она делает это энергичным образом (часто ленивость монады IO в языке Haskell приносит одни огорчения, поэтому приходится использовать функцию seq, которая энергично вычисляет значение своего первого аргумента, а возвращает значение второго). Поскольку для вычисления длины списка, его надо пробежаться от головы до самого последнего элемента, это даёт гарантию того, что всё содержимое файла быдет прочитано.

distance :: (Integral a, Floating d) => Point a -> Point a -> d
distance (Point x y z) (Point x' y' z')
  = sqrt $ fromIntegral $ (x - x')^2 + (y - y')^2 + (z - z')^2

Банальная функция, вычисляющая евклидово расстояние между двумя точками. Использование явного преобразования fromIntegral необходимо из-за типа функции sqrt, вычисляющей квадратный корень.

prettyPrintPoint :: Show a => Point a -> String
prettyPrintPoint (Point x y z) = "(" ++ show x ++ ", " ++ show y ++ ", " ++ show z ++ ")"

Наконец, функция для вывода описания точки в «приятном виде». Эта функция будет использоваться для вывода координат точек на экран для пользователя. Как видно, здесь используется математическая нотация — координаты заключаются в круглые скобки и разделяются запятыми.

В общем, модуль Point закончен. Сюда, конечно, можно было бы добавить кучу всего ещё, но для целей решения задач описанных программных сущностей вполне достаточно. Теперь можно перейти к модулю генерации случайного набора точек. Начнём, как обычно, с декларации модуля и секции импорта. Вот они:

module StarsGenerator where

import Control.Monad (replicateM)
import System.IO (Handle, IOMode(WriteMode), withFile, hPrint)
import System.Random (randomRIO)

import Point

Тут всё просто. Из стандартных модулей импортируем только то, что требуется в программе. Из написанного ранее модуля Point импортируем всё, что он экспортирует. Хотя можно было бы тоже импортировать только то, что требуется. Но тут это не так важно.

Теперь начнём разработку сверху вниз. Как обычно, с функции main. Вот её определение:

main :: FilePath -> Int -> Integer -> IO ()
main fn m n = replicateM m (generateStar n) >>= (withFile fn WriteMode . save)

Тут, конечно, печаль для неподготовленного новичка. Но если внимательно разобраться, то подобные определения будут расщёлкиваться, как орехи на завтрак. Оператор (>>=) разворачивает монаду IO с результата, возвращённого выражением replicateM m (generateStar n), и передаёт его выражению (withFile fn WriteMode . save). Последнее же и рассчитывает окончательный результат функции main.

Выражение replicateM m (generateStar n) выполняет функцию generateStar с параметром n ровно m раз и собирает результаты таких исполнений в список, оборачивая его в монаду (в данном случае, в монаду IO). Более сложным является выражение (withFile fn WriteMode . save). Чтобы его понять, легче всего обратиться к типам подвыражений. Тип слева от композиции (.) частично применённой стандартной функции withFile равен (Handle -> IO r) -> IO r. Тип функции save, которую мы ещё определим, равен [Point a] -> Handle -> IO (), но для того чтобы понять композицию, его лучше всего воспринимать как [Point a] -> (Handle -> IO ()). В итоге после композиции тип всего выражения становится равным [Point a] -> IO (). То есть всё это выражение справа от оператора связывания является функцией, которая принимает на вход список точек и ничего не возвращает.

Теперь определим недостающие функции generateStar и save. Вот определение первой:

generateStar :: Integer -> IO (Point Integer)
generateStar n = do x <- randomRIO (1, n)
                    y <- randomRIO (1, n)
                    z <- randomRIO (1, n)
                    return $ Point x y z

Ничего необычного. Параметр функции определяет верхнюю границу интервала, из которого необходимо выбирать случайные значения. Соответственно, три случайных значения выбираются из интервала [1, n], которые затем упаковывются конструктором данных Point в описание одной точки. Конечно, здесь надо было употребить аппликативный функтор, чтобы сделать определение намного проще и выразительнее: Point <$> randomRIO (1, n) <*> randomRIO (1, n) <*> randomRIO (1, n).

Ну и, наконец, последняя функция этого модуля:

save :: Show a => [Point a] -> Handle -> IO ()
save stars h = mapM_ (hPrint h) stars

Опять же, ничего необычного — берём список точек и выводим его в файл.

Таким образом, мы получили всё необходимое для генерации заданного количества точек и записи их координат в файл. Для этого необходимо выполнить выражение:

> main "stars.txt" 1000000 1000000

Это выражение запишет в файл stars.txt один миллион точек (второй параметр, первый миллион), каждая координата которых находится в интервале от 1 до 1000000 (третий параметр, второй миллион) включительно.

Подходы к решению

В целом есть два подхода к решению поставленной задачи, при этом в каждом подходе можно использовать различные алгоритмы и структуры данных для достижения всё более и более оптимального по времени решения. Ведь если попытаться сравнить все точки попарно, надо будет сделать 499 999 500 000 сравнений. Это печаль, поскольку если каждую миллисекунду делать одно сравнение, то на всё про всё уйдёт почти 16 лет. Так что надо думать.

Вот яйцеголовые учёные и придумали два метода:

1. Метод, основанный на динамическом программировании. Берутся две точки, между ними рассчитывается растояние. Оно становится текущим минимальным. Берётся следующая точка, из множества уже просмотренных точек выбираются такие точки, которые лежат от новой точки на расстоянии меньшем, чем текущее минимальное. Для всех них считается расстояние до новой точки, минимальное из которых становится новым текущим минимальным. Этот процесс повторяется для всех непросмотренных точек. Проблема заключается в нахождении эффективного алгоритма выборки уже просмотренных точек в окрестности новой точки.
2. Метод, основанный на рекурсии. Всё пространство делится плоскостью на два подпространства, в которых запускается тот же самый алгоритм. В подпространствах ищется пара точек с минимальным расстоянием между ними, после чего осуществляется ещё и поиск таких пар точек около границы разделения подпространств. Проблема заключается в нахождении эффективного алгоритма нахождения расстояний на границах.

Типовые алгоритмы имеют сложность O(n * log n), в то время как алгоритм попарного сравнения — O(n²).

В данной статье будет приведена реализация алгоритма по первому методу в двух вариантах.

Решение № 1

Первое решение будем делать «в лоб». Как вот написано в описании метода, так и запрограммируем. Ну и, как обычно, воспользуемся разработкой сверху вниз. Так что вот определение функции main (тут я уже нарочито пропущу декларации модуля и секцию импорта):

main :: FilePath -> IO ()
main fn = withFile fn ReadMode loadFromFile >>= minimalDistance >>= printSolution

Ну, опять двадцать пять! Но надо всегда помнить сигнатуру оператора (>>=), и тогда такие выражения не будут пугать. Ну а если кого уже не пугают — тому хорошо. В общем, всё, что делает эта функция, так это последовательно читает файл с перечнем координат точек, рассчитывает минимальное расстояние между всеми точками, выводит результат на экран. Собственно, теперь надо определить функции minimalDistance и printSolution (помните, функция loadFromFile определена в модуле Point).

minimalDistance :: [Point'] -> IO (Solution Double)
minimalDistance (p:q:ps) = return $ snd $ foldl countDistance ([q, p], (distance p q, (p, q))) ps

Как видно, эта функция принимает список, в котором, как минимум, два элемента (иначе как можно рассчитать минимальное расстояние?). Можно было бы определить ещё клоза для обработки ситуаций, когда в списке один или ноль элементов, но этого не сделано. Так что данное определение может быть источником потенциальных ошибок времени исполнения. Но не суть. Суть здесь в том, что она сворачивает список, который начинается с третьего элемента, причём делает это довольно забавным способом. Свёртка осуществляется в пару, первым элементом которой является список уже просмотренных точек (и тогда после полной свёртки здесь будет обращённый исходный список). А вторым элементом этой пары является текущее минимальное расстояние и пара точек, между которыми это минимальное расстояние и находится.

Здесь надо отметить, что в сигнатуре этой функции используются какие-то два типа, которых ещё не было. Не паниковать — это просто синонимы для упрощения сигнатур. Вот их определения:

type Point' = Point Integer

type Solution a = (a, (Point', Point'))

Итак, теперь определение функции для свёртки countDistance:

countDistance :: (Ord a, Floating a, RealFrac a) => ([Point'], Solution a) -> Point' -> ([Point'], Solution a)
countDistance (ss, v) s = (s:ss, minimumBy (comparing fst) (v : map (s' -> (distance s s', (s, s'))) [s' | s' <- ss, inArea s' s (fst v)]))

Тут, конечно, немного запутанно. Поскольку первым элементом пары, которую возвращает эта функция, является список уже просмотренных точек, то мы просто добавляем второй параметр (новую точку) в голову уже имеющегося там списка (потому-то список уже просмотренных точек в конце будет обращённым; и потому-то, если кто заметил, в определении функции minimalDistance в первоначальном формировании списка просмотренных точек две точки указаны в обратном порядке). Со вторым элементом пары сложнее. Здесь надо найти минимальный элемент среди всех пар, сравнивая их по первому элементу (minimumBy (comparing fst)), при этом пары формируются в виде расстояния между текущей точкой и всеми точками, находящимися в окрестности. За получение списка всех точек из окрестности отвечает генератор списка с предикатом охраны inArea. Этот предикат определяется следующим образом:

inArea :: (RealFrac a, Floating a) => Point' -> Point' -> a -> Bool
inArea (Point x' y' z') (Point x y z) md = abs (x - x') <= md' &&
                                           abs (y - y') <= md' &&
                                           abs (z - z') <= md'
  where md' = floor md

Тут нет ничего сложного — значение True возвращается тогда и только тогда, когда абсолютное значение разниц всех трёх координат не превышает текущее минимальное расстояние. И эта функция категорически неэффетивна, равно как и неэффективен генератор, в котором она используется. Из-за этого генератора программа проходит полностью по всему списку уже просмотренных точек. Сложность высокая, и программа работает более 15 суток на обычном ноутбуке (Windows 7, 4 Гб оперативной памяти) в режиме интерпретации. Что-то надо делать.

Но пока осталось реализовать функцию, которая выведет результат на экран. Она простая:

printSolution :: Show a => Solution a -> IO ()
printSolution (md, (p, q)) = putStrLn (prettyPrintPoint p ++ " - " ++
                                       prettyPrintPoint q ++ ": " ++ show md)

Вот так. Из-за неэффективности одного генератора программа не смогла бы участвовать в конкурсе, который длился всего трое суток. Думаем, как переделать…

Решение № 2

Начать рефакторинг надо с того, что желательно всё-таки определить нормальный АТД для представления решения. Пара, одним элементом которой является пара, одним элементом которой также является пара — это не язык Haskell, это чёрт знает что и с боку бантик. Так что начнём с типов:

type WPair a b = ([Point b], Solution a b)

data Solution a b = Empty
                  | Distance a (Point b) (Point b)

Синоним WPair — это опять же служебный тип для сокращения сигнатур. А вот АТД Solution представляет решение задачи. Если решение не пусто, то есть минимальное расстояние найдено, то первым элементом типа является это расстояние, а двумя другими — точки, между которыми это минимальное расстояние находится.

Теперь перейдём к главной функции модуля, которая получает на вход имя файла с описанием точек в трёхмерном пространстве и выводит на экран информацию о найденном решении. Она стандартно читает файл, потом сортирует список (это — важная часть решения), сворачивает уже отсортированный список при помощи функции addStar и выводит результат на экран. Вот тут-то собака и порылась. Мы отсортируем список по одной из координат, и это категорически сократит время поиска решения. Вот определение функции main:

main :: FilePath -> IO ()
main fn = do stars <- withFile fn ReadMode loadFromFile
             printSolution $ snd $ foldl addStar ([], Empty) $ sortBy (comparing getX) stars

Поскольку здесь выражения слишком длинны, резонно воспользоваться ключевым словом do, а не оператором связывания (>>=). Хотя и оператор смотрелся бы здесь вполне гармонично (как?).

Функцию printSolution придётся немного переделать по сравнению с предыдущим модулей, поскольку она должна принимать значение совершенно иного типа. Ну это просто:

printSolution :: (Show a, Show b) => Solution a b -> IO ()
printSolution  Empty           = putStrLn "No Solution."
printSolution (Distance d p q) = putStrLn (prettyPrintPoint p ++ " - " ++
                                           prettyPrintPoint q ++ ": " ++ show d)

Второй клоз функции полностью повторяет определение из ранее рассмотренного модуля.

Теперь функция addStar, функция, используемая для свёртки списка всех точек в решение. Для каждой ещё непросмотренной точки s она будет искать точки, лежащие в окрестности этой точки на расстоянии, не большем текущего минимального. Если список таких точек пуст, то она просто добавляет к списку уже просмотренных точек текущую, не меняя ранее найденного решения. Если же список не пуст, то она опять сворачивает его при помощи новой функции, которую мы напишем чуть позже, с частичным применением в виде текущей точки. Вот определение:

addStar :: (Ord a, Floating a, Integral b) => WPair a b -> Point b -> WPair a b
addStar (ss, slv) s = case findStars s slv ss of
                        []  -> wpair
                        ss' -> foldl (checkStar s) wpair ss'
  where
    wpair = (s:ss, slv)

Функция findStars здесь ключевая. Она возвращает список точек из окрестности заданной. Её определение категорически отличается от предиката inArea, который был определён в предыдущем решении задачи:

findStars :: (Ord a, Num a, Integral b) => Point b -> Solution a b -> [Point b] -> [Point b]
findStars  _             Empty            = id
findStars (Point x y z) (Distance md _ _) = filter    (p -> abs (fromIntegral $ z - getZ p) < md) .
                                            filter    (p -> abs (fromIntegral $ y - getY p) < md) .
                                            takeWhile (p -> fromIntegral (x - getX p) < md)

Этак функция самым категорическим образом использует тот факт, что список точек отсортирован по первой координате, из-за чего для ограничения окрестности по этой координате используется стандартная функция takeWhile, а не filter, как для других координат. В итоге по списку ещё нерассмотренных точек программа проходится только два раза (само собой разумеется, что функция takeWhile из-за своей лености останавливается сразу же, как только очередная точка выходит по расстоянию за окрестность; а из-за того, что текущее расстояние быстро уменьшается, функция takeWhile практически не проходится по списку). Ну и, к тому же, проходы по списку производятся последовательно, а не во вложенных циклах, как это было в предикате inArea и соответствующем генераторе.

В определении функции findStars используется бесточечная нотация, в связи с чем третий аргумент функции (список точек, среди которых ищутся точки в окрестности) не показан.

Наконец, функция checkStar, при помощи которой сворачивается список найденных в окрестности точек. Её определение таково:

checkStar :: (Ord a, Floating a, Integral b) => Point b -> WPair a b -> Point b -> WPair a b
checkStar s (ss, Empty) s' = (ss, Distance (distance s s') s s')
checkStar s (ss, slv@(Distance md _ _)) s' = if md' < md
                                               then (ss, Distance md' s s')
                                               else (ss, slv)
  where
    md' = distance s s'

Тут реализована та же самая логика — если находится точка, лежащая на меньшем расстоянии, чем текущее, то это расстояние заменяет собой текущее, равно как и пара точек становится новой. Если же для рассматриваемой точки расстояние больше текущего, то она просто игнорируется.

По сравнению с предыдущим решением, этот вариант работает на той же самой машине на порядки быстрее — минимальное расстояние среди миллиона точек было найдено чуть менее, чем за полтора часа в режиме интерпретации.

Приз зрительских симпатий

Великолепный dastapov написал краткую и очень выразительную программу на языке Haskell, которая и завоевала приз зрительских симпатий. Программа работает очень быстро, использует только стандартные средства работы с данными, а также модуль с разработанным алгоритмом k-d tree. Вот её текст полностью:

module Main where

import Control.Applicative ((<$>))
import Data.List (delete)
import Debug.Trace (trace)

import qualified Data.Trees.KdTree as KDT

main = do coords <- map (map read . words) . lines <$> getContents
          let points = map ([x, y, z] -> KDT.Point3d x y z) coords
              kdt = KDT.fromList points
              start_dist = 10^20
          print $ findBest kdt points (start_dist, start_dist, head points, head points)
  where
    findBest kdt [] (dist_sq, dist, p1, p2) = (dist, p1 ,p2)
    findBest kdt (p:rest) best@(dist_sq, dist, p1, p2) =
      case delete p (KDT.nearNeighbors kdt dist p) of
        []     -> findBest kdt rest best
        (p2:_) -> let d = KDT.dist2 p p2
                  in if d < dist_sq
                       then findBest kdt rest $
                              trace (show (length rest) ++ " to go, new best " ++ show (d, p, p2))
                                    (d, sqrt d, p, p2)
                       else findBest kdt rest best

Заключение

На этот раз я не стал рисовать диаграммы вызовов, чтобы не загромождать статью. Тем более, что пока ничего более конкретного, кроме самих диаграмм, я по нима сказать не могу.

В качестве незначительных замечаний, которые необходимо исправить для представления исходного кода модулей более «зрелым», необходимо отметить:

1. Исправление определения метода show для типа Point в одноимённом модуле.
2. Применение аппликативного функтора для определения функции generateStar в модуле StarsGenerator.
3. Сделать функцию minimalDistance безопасной.

Все желающие могут невозбранно скачать полные модули с приведёнными в данной статье программными сущностями по следующи адресам:

• Модуль Points
• Модуль StarsGenerator
• Модуль Distance (Решение № 1)
• Модуль Distance (Решение № 2)
• Модуль Distance (Решение на приз зрительских симпатий)

Как обычно, не могу не сказать о том, что у Фонда Поддержки Функционального Программирования ФП(ФП) имеется свой официальный форум, на котором всегда можно пообщаться с единомышленниками, узнать ответы на интересующие вас вопросы по ФП, предложить задачи для конкурсов, ну и во всяком случае себя показать и людей посмотреть. Всё это невозбранно можно сделать здесь.

Поскольку я занимаюсь организацией и проведением конкурсов, подведением итогов и описанием статей на абсолютно безвозмездной основе, идея «краудфандинга» мне особенно близка. Так что если вы считаете для себя возможным отблагодарить автора за хорошие материалы, но не знаете как это можно сделать, то прошу проследовать по данной ссылке.

Мои предыдущие статьи о конкурсах по ФП на Хаброхабре:

• Решение задачи о перегорающих лампочках на языке Haskell
• Фреймворк для решения задач о переправах на языке Haskell
• Конструирование и вычисление арифметических выражений на языке Haskell
• Расшифровка кода на языке Haskell (конкурс по ФП в январе 2012)
• Шахматные задачи на мат в один ход: решение на языке Haskell
• Измерение объёмов при помощи двух заданных сосудов: решение на языке Haskell
• Трансмутации слов друг в друга: решение на языке Haskell
• Решение арифметических задач — вероятностный подход против регулярных выражений

Автор: Darkus