Решение задачи коммивояжера методом поиска ближайшего соседа на haskell

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

Итак, сформулируем варианты алгоритма ближайшего соседа, которые мы собираемся реализовать:

1. Берем первую точку из списка и находим ближайшую к ней из остатка списка. Получаем отрезок из двух точек. Для каждой из концевых точек отрезка ищем ближайшую точку из остатка списка, соединяем найденную точку с соответствующим концом отрезка - получаем кривую из трех точек. Проделываем то же самое с этой кривой - получаем кривую из четырех точек. И так далее, пока остаток списка точек не исчерпается и мы получим кривую, соединяющую все исходные точки - это и есть решение задачи. Заметим, что кривая, наращиваемая в процессе работы алгоритма имеет направление - это важно для механизма присоединения точек. Назовем первую концевую точку кривой головой, а последнюю - хвостом.
2. Этот вариант я назвал параллельным, хотя к параллельности вычисления это название не имеет отношения. Здесь примерно то же самое. Начинаем также с отрезка из двух точек, но затем не продолжаем его наращивать, а проделываем то же самое для остатка списка точек - получаем еще один отрезок и так далее, пока не превратим исходный список точек в список отрезков и, возможно (если полное число точек нечетное), оставшейся одиночной точки. Затем соединяем между собой отрезки, получая список кривых, состоящих по большей части из четырех точек (очевидно, из-за неравенства полного числа точек некоторой степени двойки, на этом этапе могут появиться 1-,2- и 3-точечные кривые). Проделываем то же самое, пока не соединим все кривые в одну. Это и есть решение задачи. Осталось показать, как искать и соединять ближайшие кривые. Итак, у нас есть кривая, для которой надо найти ближайшую. Берем две концевые точки заданной кривой, находим ближайшую концевую точку из остатка списка кривых (заметьте, в первом варианте алгоритма остатком был список точек, а не кривых) - если эта точка одиночная, то присоединяем ее так же, как в первом варианте алгоритма, если же это концевая точка кривой или отрезка, то возможны четыре варианта соединения этой кривой с заданной:
   1. Хвост-голова (TH). Ближайшими оказались хвост заданной кривой и голова ближайшей. В этом случае просто присоединяем к заданной кривой найденную.
   2. Голова-хвост (HT). Голова заданной кривой и хвост ближайшей. Присоединяем к найденной кривой заданную.
   3. Хвост-хвост (TT). Присоединяем к заданной кривой обращенную (reversed) найденную.
   4. Голова-голова (HH). Присоединяем к обращенной найденной кривой заданную.
   Возможны вырожденные случаи, когда одна из кривых (или обе) является одиночной точкой. Случай кривая-точка мы уже рассмотрели, формально тут можно выделить два варианта: TH и HT (так как вторая кривая вырождена, то ее единственную точку можно рассматривать и как хвост, и как голову, но нам будет удобнее использовать указанные два варианта, так как они проще реализуются). Если обе кривых - одиночные точки, то представим их соединение как TH. Если же первая кривая - невырожденная кривая, а вторая - одиночная точка, то мы опять можем соединить их двумя формальными способами TH и HT.

Давайте реализуем функцию findClosest, которая будет находить для заданной кривой ближайшую кривую из списка кривых, а также вычислять способ их соединения. Я приведу код, а затем его прокомментирую.

findClosest :: (Floating a, Ord a) =>
    [(a, a)] -> [[(a, a)]] -> ([(a, a)], Maybe Int)
findClosest x [] = ([], Just 0)
findClosest x xs = foldr (\x' y' -> findCloser' x' y' x) ([], Just 0) xs
    where findCloser' x y z =
              let zx = z `distance'` x
                  zy = z `distance'` (fst y)
              in if fst zx < fst zy
                  then (x, snd zx)
                  else (fst y, snd zy)
                  where distance' []     _      = (infiniteDistance, Just 0)
                        distance' _      []     = (infiniteDistance, Just 0)
                        distance' [x]    [y]    = (distance (x, y),  Just 0)
                        distance' (x:xs) [y]    = (d, i)
                            where x' = last xs
                                  zs = [(x', y), (x, y)]
                                  z  = minimumBy findCloser zs
                                  d  = distance z
                                  i  = elemIndex z zs
                        distance' [x]    (y:ys) = (d, i)
                            where y' = last ys
                                  zs = [(x, y), (x, y')]
                                  z  = minimumBy findCloser zs
                                  d  = distance z
                                  i  = elemIndex z zs
                        distance' (x:xs) (y:ys) = (d, i)
                            where x' = last xs
                                  y' = last ys
                                  zs = [(x', y), (x, y'), (x', y'), (x, y)]
                                  z  = minimumBy findCloser zs
                                  d  = distance z
                                  i  = elemIndex z zs

Вспомогательные функции:

distance :: Floating a =>
    ((a, a), (a, a)) -> a
distance ((x0, y0), (x1, y1)) = sqrt $ (x0 - x1) ^ 2 + (y0 - y1) ^ 2

findCloser :: (Num a, Ord a) =>
    ((a, a), (a, a)) -> ((a, a), (a, a)) -> Ordering
findCloser ((x0, y0), (x1, y1)) ((x2, y2), (x3, y3)) = compare d1 d2
    where d1 = (x0 - x1) ^ 2 + (y0 - y1) ^ 2
          d2 = (x2 - x3) ^ 2 + (y2 - y3) ^ 2

infiniteDistance :: Floating a => a
infiniteDistance = 1e15

Для использования функций minimumBy и elemIndex в теле findClosest необходимо включить модуль Data.List.

Начнем со вспомогательных функций. Функция distance рассчитывает расстояние между двумя точками на плоскости, ее реализация элементарна. В функцию findCloser передаются два отрезка (расстояния между различными концевыми точками двух кривых), она используется в качестве предиката сравнения для minimumBy в теле findClosest, ее реализация тоже элементарна. Функция infiniteDistance просто объявляет какую-то большую константу (1e15) - она будет использоваться в качестве значения расстояния для вырожденных случаев, соответствующих отсутствию одной из двух кривых, в функции distance', объявленной внутри тела findClosest.

Что же делает функция findClosest? Как видим, она принимает в качестве первого параметра список точек: он будет соответствовать заданной кривой из описания алгоритмов ближайшего соседа. Второй параметр findClosest - список списков точек - это остаток списка одиночных точек из первого варианта алгоритма, или остаток списка кривых из второго варианта. Возвращаемым значением является пара, первый элемент которой - найденная ближайшая кривая (или одиночная точка) из остатка списка кривых (точек). Второй элемент возвращаемой пары - один из четырех способов соединения - TH, HT, TT или HH, обернутый в Maybe Int. Значения Just 0 .. Just 3 соответствуют указанным вариантам соединения в порядке их перечисления. Значение Nothing не может возникнуть и (как впрочем и само появление обертки Maybe в возвращаемом значении findClosest) является артефактом применения функции elemIndex в distance'.

Итак, функция findClosest пробегает по списку кривых (второй аргумент функции) с целью поиска ближайшей к заданной (первый аргумент) с помощью правой свертки foldr. Исходным значением аккумулятора свертки является ([], Just 0), а комбинирующей функцией - локально объявленная функция findCloser', которая проверяет, находится ли текущая кривая из списка кривых ближе к заданной, чем та, которая на данном этапе записана в первый элемент аккумулятора. Если это так, то findCloser' возвращает новую ближайшую кривую и способ соединения ее с заданной (и, конечно, эти значения запоминаются в аккумуляторе). Расстояние между кривыми и способ их соединения рассчитываются внутри функции distance', для которой определены аж 6 шаблонов. Первые два соответствуют вырожденным случаям, когда одна из кривых не задана - в этом случае расстояние должно быть равно бесконечности (в нашем случае бесконечности соответствует infiniteDistance), а способ соединения - неважно какой (в нашем случае это Just 0, то есть TH). Третий шаблон соответствует расстоянию между двумя одиночными точками, четвертый и пятый - когда одна из кривых вырождена в одиночную точку, шестой - когда обе кривые состоят из более чем одной точки. В трех последних случаях составляются комбинации zs из концевых точек кривых, находится элемент z, соответствующий минимальному расстоянию между парами точек, и, с помощью функции elemIndex, определяется его индекс внутри списка zs. Поскольку список zs составлен хитрым способом, то по значению этого индекса можно определить, каким способом следует соединить найденные кривые.
 
Проверим, как работает findClosest (при условии, что вы сохранили эту функцию в файле и загрузили его в ghci, либо загрузили готовый модуль, доступный по ссылке в конце этой статьи).

*ShortestPath> findClosest [(0,1),(0,6),(10,23)] [[(0,2),(0,5)],[(5,6)]]
([(0.0,2.0),(0.0,5.0)],Just 3)
*ShortestPath>

Здесь задана кривая, ограниченная точками (0,1) и (10,23), нужно найти ближайшую к ней кривую из списка кривых [[(0,2),(0,5)],[(5,6)]]. Очевидно, это дожна быть кривая [(0,2),(0,5)], поскольку точка (0,2) ближе к точке (0,1) заданной кривой по сравнению с другими возможными вариантами. Эта кривая выведена в первом элементе ответа. Как должны быть соединены заданная кривая [(0,1),(0,6),(10,23)] с ближайшей к ней кривой [(0,2),(0,5)]? Точка (0,1) из заданной кривой - это ее голова, точка (0,2) из найденной ближайшей кривой - тоже голова, соответственно реализуется соединение НН, то есть Just 3, что мы и видим в ответе. Имеет смысл поиграться еще, но я больше не буду этого делать, функция findClosest работает правильно!

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

moveClosest :: (Floating a, Ord a) =>
    [(a, a)] -> [[(a, a)]] -> ([(a, a)], [[(a, a)]])
moveClosest x [] = (x, [])
moveClosest x xs
    | i == Just 0 = (x ++ y, delete y xs)
    | i == Just 1 = (y ++ x, delete y xs)
    | i == Just 2 = (x ++ reverse y, delete y xs)
    | otherwise   = (reverse y ++ x, delete y xs)
    where (y, i) = findClosest x xs

Эта функция реализована с помощью функции findClosest и использует информацию, возвращаемую последней, для определения способа соединения заданной кривой с найденной ближайшей. Реализация способа соединения соответствует методу из описания второго варианта алгоритма ближайшего соседа. Проверяем.

*ShortestPath> moveClosest [(0,1),(0,6),(10,23)] [[(0,2),(0,5)],[(5,6)]]
([(0.0,5.0),(0.0,2.0),(0.0,1.0),(0.0,6.0),(10.0,23.0)],[[(5.0,6.0)]])
*ShortestPath>

Всё верно. Следующий шаг - реализация обоих вариантов алгоритма ближайшего соседа.

Первый вариант.

shortestPath :: (Floating a, Ord a) =>
     [[(a, a)]] -> [(a, a)]
shortestPath []     = []
shortestPath [x]    = x
shortestPath (x:xs) = fst $ shortestPath' $ moveClosest x xs
    where shortestPath' (x, []) = (x, [])
          shortestPath' (x, y)  = shortestPath' $ moveClosest x y

Очень просто! Функция shortestPath (хотя правильнее было бы назвать ее shortPath) принимает список кривых (для задачи коммивояжера список кривых вырождается в список одиночных точек) и возвращает короткий (не кратчайший!) незамкнутый путь, соединяющий их все. Внутри тела shortestPath определена вспомогательная функция shortestPath', которая рекурсивно вызывает moveClosest до тех пор, пока все кривые из исходного списка не будут соединены.

Второй вариант алгоритма чуть-чуть посложнее.

shortestPathP :: (Floating a, Ord a) =>
     [[(a, a)]] -> [(a, a)]
shortestPathP []  = []
shortestPathP [x] = x
shortestPathP xs  = shortestPath' $ combineClosest xs
    where shortestPath' []  = []
          shortestPath' [x] = x
          shortestPath' xs  = shortestPath' $ combineClosest xs'
              where xs' = combineClosest xs

Последняя буква P в названии shortestPathP обозначает Parallel. Вспомогательная функция combineClosest рекурсивно соединяет каждую первую кривую (это будет заданная кривая) из остатка кривых с найденной ближайшей к ней кривой до полного исчерпания остатка кривых. Затем список уже соединенных кривых подвергается такой же обработке до тех пор, пока не будет получена единственная кривая, являющаяся решением задачи. Вот реализация combineClosest:

combineClosest :: (Floating a, Ord a) =>
    [[(a, a)]] -> [[(a, a)]]
combineClosest (x:xs) = moveClosest' x xs []
    where moveClosest' x [] t       = x : t
          moveClosest' x xs@[_] t   = (fst $ moveClosest x xs) : t
          moveClosest' x xs@(_:_) t = moveClosest' x' y' t'
              where z  = moveClosest x xs
                    x' = head $ snd z
                    y' = tail $ snd z
                    t' = fst z : t

Список t в реализации вспомогательной функции moveClosest' нужен для хранения уже соединенных кривых. Исключений, связанных с вызовом функций head и tail для пустых списков в определении moveClosest' x xs@(_:_), возникнуть не может, поскольку для данного шаблона функции будет вызвана функция moveClosest, у которой второй аргумент (список кривых) состоит как минимум из двух элементов. В этом случае второй элемент возвращенного moveClosest значения (то есть список кривых после соединения ближайшей найденной кривой с заданной) будет состоять как минимум из одного элемента, поэтому на вход функций head и tail в moveClosest' никогда не попадет пустой список. Проверим, как работает функция combineClosest.

*ShortestPath> combineClosest [[(0,1),(0,6),(10,23)],[(0,2),(0,5)],[(5,6)]]
[[(5.0,6.0)],[(0.0,5.0),(0.0,2.0),(0.0,1.0),(0.0,6.0),(10.0,23.0)]]
*ShortestPath>

Замечательно! Исходный список состоит всего из трех кривых, две из них были соединены правильно, а третья осталась в одиночестве. Проверьте combineClosest на большем количестве кривых, результаты должны быть верными.

Ну и, собственно, функция, рассчитывающая длину кривой.

pathDistance :: (Floating a, Ord a) =>
    [(a, a)] -> a
pathDistance []        = 0
pathDistance [_]       = 0
pathDistance zs@(_:xs) = foldr (\x' y' -> distance x' + y') 0 (zip zs $ xs)

Здесь все просто. С помощью zip zs $ xs составляем список отрезков из соседних точек в заданной кривой и затем сворачиваем его, используя функцию, находящую расстояние между двумя точками элемента из списка отрезков. Проверяем.

*ShortestPath> pathDistance [(0,0),(1,1),(2,2),(3,3),(4,4)]
5.656854249492381
*ShortestPath>

Ответ верный, это число равно квадратному корню из 32.

А теперь давайте поиграем с большим количеством точек. У меня есть файл, который я назвал tsp.dat,  с такими данными:

20833.3333 17100.0000
20900.0000 17066.6667
21300.0000 13016.6667
21600.0000 14150.0000
21600.0000 14966.6667
21600.0000 16500.0000
22183.3333 13133.3333
22583.3333 14300.0000
22683.3333 12716.6667
23616.6667 15866.6667
23700.0000 15933.3333
23883.3333 14533.3333
24166.6667 13250.0000
25149.1667 12365.8333
26133.3333 14500.0000
26150.0000 10550.0000
26283.3333 12766.6667
26433.3333 13433.3333
26550.0000 13850.0000
26733.3333 11683.3333
27026.1111 13051.9444
27096.1111 13415.8333
27153.6111 13203.3333
27166.6667 9833.3333
27233.3333 10450.0000

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

calcFromFile' :: String -> Bool -> Bool -> IO ()
calcFromFile' file p r = do
    content <- readFile file
    let fileLines        = if r
                               then reverse $ lines content
                               else lines content
        fileData         = map readFloat fileLines
        path             = if p
                               then shortestPathP fileData
                               else shortestPath fileData
        distance         = pathDistance path
        printPair (x, y) = printf "(%10.4f, %10.4f)\n" x y
    putStrLn "*** Original set ***"
    mapM_ (mapM_ printPair) fileData
    putStrLn "\n*** Short path ***"
    mapM_ printPair path
    putStrLn "\n*** Distance ***"
    print distance
    putStrLn ""

readFloat :: String -> [(Float, Float)]
readFloat x = [(read y :: Float, read z :: Float)]
    where y:z:_ = words x

calcFromFile :: String -> IO ()
calcFromFile file = calcFromFile' file False False

calcFromFileR :: String -> IO ()
calcFromFileR file = calcFromFile' file False True

calcFromFileP :: String -> IO ()
calcFromFileP file = calcFromFile' file True False

calcFromFilePR :: String -> IO ()
calcFromFilePR file = calcFromFile' file True True

Функция calcFromFile' принимает три аргумента: название файла, из которого будут читаться данные, логическое значение p, которое определяет, какой из двух вариантов алгоритма ближайшего соседа будет использован (если p равно True, то будет использован параллельный алгоритм), и еще одно логическое значение r, в случае равенства которого True, исходный список будет предварительно обращен (reversed). Соответственно, определены четыре более высокоуровневые функции calcFromFile, calcFromFileR, calcFromFileP и calcFromFilePR, принимающие только один аргумент - имя файла - и комбинирующие аргументы p и r функции calcFromFile' четырьмя возможными способами. Преобразование строки, состоящей из двух чисел с плавающей точкой в список, состоящий из единственной точки происходит во вспомогательной функции readFloat. Надо отметить, что функция calcFromFile' в своей реализации использует функцию printf, а это значит, что нужно подключить еще один модуль Text.Printf.

Рассчитаем короткий путь между точками из файла tsp.dat.

*ShortestPath> calcFromFile "tsp.dat"
*** Original set ***
(20833.3340, 17100.0000)
(20900.0000, 17066.6660)
(21300.0000, 13016.6670)
(21600.0000, 14150.0000)
(21600.0000, 14966.6670)
(21600.0000, 16500.0000)
(22183.3340, 13133.3330)
(22583.3340, 14300.0000)
(22683.3340, 12716.6670)
(23616.6660, 15866.6670)
(23700.0000, 15933.3330)
(23883.3340, 14533.3330)
(24166.6660, 13250.0000)
(25149.1660, 12365.8330)
(26133.3340, 14500.0000)
(26150.0000, 10550.0000)
(26283.3340, 12766.6670)
(26433.3340, 13433.3330)
(26550.0000, 13850.0000)
(26733.3340, 11683.3330)
(27026.1110, 13051.9440)
(27096.1110, 13415.8330)
(27153.6110, 13203.3330)
(27166.6660,  9833.3330)
(27233.3340, 10450.0000)

*** Short path ***
(20833.3340, 17100.0000)
(20900.0000, 17066.6660)
(21600.0000, 16500.0000)
(21600.0000, 14966.6670)
(21600.0000, 14150.0000)
(22583.3340, 14300.0000)
(22183.3340, 13133.3330)
(22683.3340, 12716.6670)
(21300.0000, 13016.6670)
(24166.6660, 13250.0000)
(23883.3340, 14533.3330)
(23616.6660, 15866.6670)
(23700.0000, 15933.3330)
(26133.3340, 14500.0000)
(26550.0000, 13850.0000)
(26433.3340, 13433.3330)
(27096.1110, 13415.8330)
(27153.6110, 13203.3330)
(27026.1110, 13051.9440)
(26283.3340, 12766.6670)
(26733.3340, 11683.3330)
(26150.0000, 10550.0000)
(27233.3340, 10450.0000)
(27166.6660,  9833.3330)
(25149.1660, 12365.8330)

*** Distance ***
26575.813

Давайте посмотрим, что мы посчитали. Для этого сначала изобразим на графике исходные точки.

А теперь построим полученный короткий путь.

Проделаем то же самое с тремя оставшимися комбинациями p и r. Я не стану показывать вывод программы дабы не загромождать статью. Приведу лишь результаты рассчитанных длин: (p = False, r = True): 25318.4;  (p = True, r = False): 26018.463; (p = True, r = True): 29142.959, и картинки.

Самым коротким оказался путь, рассчитанный с помощью первого варианта алгоритма с предварительно обращенным списком точек. Только это ничего не значит! Во-первых, мы, скорее всего, вообще не получили кратчайшего пути, а во-вторых, те короткие пути, которые мы получили, зависят от порядка перечисления исходных точек в оригинальном списке! Единственное, что здесь можно гарантировать - это то, что полученные пути являются короткими (насколько бы неформально это ни звучало), и это видно на картинках.

Исходный код примера можно взять здесь. В тарболе, кроме него, находятся файл tsp.dat, результаты расчетов коротких путей (другие файлы с расширением .dat) и скрипт tsp.m  для octave, который можно использовать для построения приведенных изображений.

Update. Функции ввода-вывода можно переписать следующим образом:

type Transform = ([[(Float, Float)]] -> [[(Float, Float)]])
type Algorithm = ([[(Float, Float)]] -> [(Float, Float)])

calcFromFile' :: String -> Transform -> Algorithm -> IO ()
calcFromFile' file tr alg = do
    content <- readFile file
    let fileData         = map readFloat (lines content)
        path             = alg . tr $ fileData
        distance         = pathDistance path
        printPair (x, y) = printf "(%10.4f, %10.4f)\n" x y
    putStrLn "*** Original set ***"
    mapM_ (mapM_ printPair) fileData
    putStrLn "\n*** Short path ***"
    mapM_ printPair path
    putStrLn "\n*** Distance ***"
    print distance
    putStrLn ""

readFloat :: String -> [(Float, Float)]
readFloat x = [(read y :: Float, read z :: Float)]
    where y:z:_ = words x

calcFromFile :: String -> IO ()
calcFromFile file = calcFromFile' file id shortestPath

calcFromFileR :: String -> IO ()
calcFromFileR file = calcFromFile' file reverse shortestPath

calcFromFileP :: String -> IO ()
calcFromFileP file = calcFromFile' file id shortestPathP

calcFromFilePR :: String -> IO ()
calcFromFilePR file = calcFromFile' file reverse shortestPathP

Маловыразительные флаги-объекты p и r, которые передавались в функцию calcFromFile' в предыдущей версии, теперь заменены на аргументы, являющиеся функциями типа Algorithm и Transform, предоставляющими конкретные реализации алгоритма ближайшего соседа и трансформации исходного списка кривых соответственно. Это решение намного более гибкое, чем предыдущее, в котором конкретные алгоритмы были зашиты в тело функции calcFromFile'.

Я намеренно оставил предыдущую реализацию функций ввода-вывода, так как различие между ней и новой реализацией, на мой взгляд, прекрасно иллюстрирует эволюцию мышления, зашоренного императивными и ООП шаблонами, при изучении функционального подхода в программировании. В самом деле, объекты типа p и r из предыдущей реализации calcFromFile', являющиеся неотъемлемой частью модели императивного мышления, в функциональном подходе оказываются объектами-паразитами. Они не приносят ничего нового, но при этом усложняют понимание кода и лишают его гибкости без всякой на то необходимости.

Хочу отметить еще один интересный момент. Мы могли бы объявить тип Transform как ([String] -> [String]), а локальные определения внутри let в функции calcFromFile' записать так:

    let fileData         = map readFloat (tr . lines $ content)
        path             = alg fileData
        distance         = pathDistance path
        printPair (x, y) = printf "(%10.4f, %10.4f)\n" x y

В этом случае трансформации подвергался бы не список кривых типа [(Float, Float)], а строки в исходном файле. А в качестве конкретных функций-трансформеров можно было бы по-прежнему использовать id и reverse. Это возможно, несмотря на строгую статическую типизацию, благодаря мощной реализации параметрического полиморфизма в haskell. В самом деле, типы id и reverse, a -> a и [a] -> [a] соответственно, полиморфны, и оба варианта определения типа Transform прекрасно согласуются с  этими шаблонами.