Как я упоминал в посте пару дней назад, я писал алгоритм на основе ближайших пар в Haskell, и хотя версия с перебором работает для небольшого числа пар, она начинает распадаться по мере увеличения числа пар:
time ./closest_pairs 100 bf ./closest_pairs 100 bf 0.01s user 0.00s system 87% cpu 0.016 total time ./closest_pairs 1000 bf ./closest_pairs 1000 bf 3.59s user 0.01s system 99% cpu 3.597 total time ./closest_pairs 5000 bf ./closest_pairs 5000 554.09s user 0.36s system 99% cpu 9:14.46 total
К счастью, мы можем использовать алгоритм «разделяй и властвуй», который сокращает время работы от O (n 2
Этот алгоритм определен так :
- Сортировка точек по оси X
- Разбейте множество точек на два равных по размеру подмножества вертикальной линией x = x mid
- Решите проблему рекурсивно в левом и правом подмножествах. Это даст минимальные расстояния слева и справа d Lmin и d Rmin соответственно.
- Найдите минимальное расстояние d LRmin между парой точек, в которой одна точка лежит слева от делительной вертикали, а вторая — справа (расщепленная пара).
- Окончательный ответ является минимальным среди d Lmin , d Rmin и d LRmin .
На шаге 4, в котором мы будем искать ближайшую пару, где значения x и y находятся в противоположных подмножествах, нам не нужно снова рассматривать весь список точек, поскольку мы уже знаем, что ближайшая пара точек больше не отделена чем dist = min (d Lmin , d Rmin ).
Поэтому мы можем отфильтровать много значений , лишь сохраняя ценности , которые находятся в пределах диста среднего значения х.
Если точка находится дальше , чем мы уже знаем , что расстояние между ним и любой точкой на другой стороне будет больше диста поэтому нет смысла рассматривать его.
Я использовал версию C # из Rosetta Code в качестве шаблона, и вот чем я закончил:
import Data.Maybe
import Data.List
import Data.List.Split
import Data.Function
dcClosest :: (Ord a, Floating a) => [Point a] -> (Point a, Point a)
dcClosest pairs
if length pairs <= 3 then = fromJust $ bfClosest pairs
else
foldl (\closest (p1:p2:_) -> if distance (p1, p2) < distance closest then (p1, p2) else closest)
closestPair
(windowed 2 pairsWithinMinimumDelta)
where sortedByX = sortBy compare pairs
(leftByX:rightByX:_) = chunk (length sortedByX `div` 2) sortedByX
closestPair = if distance closestLeftPair < distance closestRightPair then closestLeftPair else closestRightPair
where closestLeftPair = dcClosest leftByX
closestRightPair = dcClosest rightByX
pairsWithinMinimumDelta = sortBy (compare `on` snd) $ filter withinMinimumDelta sortedByX
where withinMinimumDelta (x, _) = abs (xMidPoint - x) <= distance closestPair
where (xMidPoint, _) = last leftByX
Если количество пар равно 3 или меньше, мы можем просто использовать алгоритм грубой силы, поскольку разделение списка пополам не так хорошо работает со списком из менее чем 4 элементов.
Основная функция — это сгиб, который проходит по всем парам точек, где мы определили, что может быть ближайшая пара с одной точкой с левой стороны от нашего «отсортировано по x списку», а другой с правой стороны. ,
Мы используем « оконную » функцию для объединения точек, которая в этом случае делает что-то вроде этого:
> windowed 2 [(0,0), (1,1), (2,2), (1.1, 1.2)] [[(0.0,0.0),(1.0,1.0)],[(1.0,1.0),(2.0,2.0)],[(2.0,2.0),(1.1,1.2)]]
Мы передаем ближайшую пару, которую мы нашли из пар точек слева и справа от списка, в качестве начального значения и проверяем, ближе ли какая-либо из расщепленных пар, чем эта. К тому времени, когда фолд закончится, у нас будет самая близкая пара в списке.
Код в разделе where используется для разбиения списка на две половины, сортировки по x и определения, какая сторона списка имеет ближайшую пару. Все это делается рекурсивно, и затем последний бит кода определяет, какие значения нам нужно учитывать для части алгоритма с расщепленными парами.
Я узнал о функции « on » при написании этого алгоритма, который позволяет очень легко выбрать функцию для сортировки коллекции, например, «сравнить` on` snd »позволяет нам сортировать список в порядке возрастания на основе второго значения в кортеже ,
Одна из проблем, связанных с тем, как я написал этот алгоритм, заключается в том, что он уделяет большое внимание биту кода с разделенной парой, хотя на самом деле большая часть того, почему он работает, состоит в том, что мы сортируем его в порядке, который облегчает его выполнение. работать с, а затем делить задачу на 2 каждый раз.
В настоящее время я думаю, что, возможно, мне следовало бы иметь этот код в основной части функции, а не прятать его в секции where.
На самом деле нет необходимости выполнять сгибание на всем протяжении, так как мы часто будем знать, что не будет разделенной пары, прежде чем мы пройдемся по всем парам. Я думал, что это читается довольно хорошо, поскольку это думается, и это бежит довольно быстро, поскольку это так или иначе!
Используя тот же набор данных, что и раньше (и получая те же ответы!):
> time ./closest_pairs 1000 dc ./closest_pairs 1000 dc 0.02s user 0.00s system 91% cpu 0.024 total > time ./closest_pairs 5000 dc ./closest_pairs 5000 dc 0.11s user 0.01s system 97% cpu 0.118 total
Я описал код для генерации случайных чисел в предыдущем посте в блоге, но я включу его снова, чтобы показать, как работает алгоритм:
import Control.Monad.State (State, evalState, get, put) import System.Random (StdGen, mkStdGen, random) import System type R a = State StdGen a rand :: R Double rand = do gen <- get let (r, gen') = random gen put gen' return r randPair :: R (Double, Double) randPair = do x <- rand y <- rand return (x,y) runRandom :: R a -> Int -> a runRandom action seed = evalState action $ mkStdGen seed normals :: R [(Double, Double)] normals = mapM (\_ -> randPair) $ repeat () main = do args <- getArgs let numberOfPairs = read (head args) :: Int if length args > 1 && args !! 1 == "bf" then putStrLn $ show ( (bfClosest $ take numberOfPairs $ runRandom normals 42)) else putStrLn $ show ( (dcClosest $ take numberOfPairs $ runRandom normals 42))