Не только лишь все

Мне всегда нравилось это высказывание Кличко, но не потому, что оно какое-то глупое или смишное, на что указывают многочисленные комментаторы, большинство из которых сами вполне достойны этих двух эпитетов. Просто есть в этой фразе ссылка на некоторую неинтуитивную теорию множеств. В самом деле, почему все должны принадлежать всем? Нет ли во множестве всех какого-то подмножества, всем не принадлежащего? Ничего не напоминает? Да это же шикарнейший парафраз парадокса Бертрана Рассела для наивной теории множеств Кантора, при преодолении которого были открыты целые новые направления в математике и формальной логике, в том числе лямбда-исчисление и, следовательно, основы современных функциональных языков программирования! К чему это я? Да к тому, что глупым людям что ни скажи, им всё или смешно или непонятно, а умные могут кучу полезных научных теорий наваять. P.S. Все животные все, но некоторые все́е.

Haskell, производительность, criterion

Стали с приятелем меряться скоростями на больших данных в задаче о монетах из предыдущей статьи. Он на C++, я на haskell. Само собой, даже самое быстрое решение из всех в ней представленных — с Map — очень, значительно уступает простому оптимизированному итеративному алгоритму на C. Из серьезных проблем решения с Map я вижу как минимум две. Во-первых, minNumCoinsMap генерирует огромный список от 0 до m, а затем берет его последний элемент. Список — классная структура данных в теории функционального программирования, вот только не очень быстрая, а для получения последнего элемента нужно пробежаться по всем элементам, начиная с первого. Все это очень нездо́рово, когда число m равно нескольким миллионам. Во-вторых, я использовал хотя и быструю, но неизменяемую (immutable) реализацию Map. Это значит, что на каждом шаге внутри функции step функция insert возвращает новый ассоциативный массив с уже накопленными данными и одним новым — очень расточительный подход! Чтобы хотя бы на порядок приблизиться к скорости алгоритма на C++, о списке и неизменяемом Map нужно забыть. Хотя на самом деле у меня было четыре итерации постепенного улучшения производительности за счет изменения деталей алгоритма и примитивов, его реализующих. На первой итерации я оставил immutable Map, но в функции step перестал копировать его полностью. В самом деле, алгоритм на шаге m проверяет уже полученные решения, записанные в массив, начиная со значения m - max c, где max c — максимальный номинал монеты в списке c. Соответственно, строки с определением функции step приняли следующий вид.

            where step a x = let cur = min (s ! (m - x) + 1) $ snd a
                             in (insert m cur (snd $ split (m - mc) s), cur)
                    where mc = maximum c

Функция split определена в модуле Data.IntMap.Strict и ее нужно импортировать. Это немного улучшило производительность. На второй итерации я заменил неизменяемый Map на изменяемый (mutable) Array из модуля Data.Array.IO. Новая функция minNumCoinsArray представлена ниже.

import           Data.Array.IO
import           Data.Array.Base             (unsafeRead, unsafeWrite)

-- ----------

minNumCoinsArray :: [Int] -> Int -> IO Int
minNumCoinsArray c m = do
    a <- newArray (0, m) 0 :: IO (IOUArray Int Int)
    mapM_ (mnc cs a) [1 .. m]
    unsafeRead a m
    where cs = sortBy (flip compare) c
          mnc :: [Int] -> IOUArray Int Int -> Int -> IO ()
          mnc c a m = foldM (step a) m (dropWhile (> m) c) >>= unsafeWrite a m
            where step a b x = liftM (min b . succ) $ unsafeRead a $ m - x

Этот массив работает с unboxed данными, расположенными непрерывно на участке памяти в куче от 0 до m — это очень эффективная структура данных (и уже не чистая). Кроме того, здесь оптимизирована свертка по списку монет c: сначала монеты сортируются в обратном порядке, а внутри свертки линейная функция filter заменена на dropWhile: в нормальной ситуации, когда самый старший номинал значительно меньше полной стоимости m, функция dropWhile будет почти всегда останавливаться на первом же элементе. Эта реализация при больших значениях m отставала от референсной программы на C++ на порядок. На третьей итерации я всего лишь заменил mutable unboxed Array на mutable unboxed Vector — чисто посмотреть, что из этого выйдет. Вот реализация функции minNumCoinsVector.

import qualified Data.Vector.Unboxed.Mutable as V

-- ----------

minNumCoinsVector :: [Int] -> Int -> IO Int
minNumCoinsVector c m = do
    v <- V.new $ m + 1
    V.unsafeWrite v 0 0
    mapM_ (mnc cs v) [1 .. m]
    V.unsafeRead v m
    where cs = sortBy (flip compare) c
          mnc c v m = foldM (step v) m (dropWhile (> m) c) >>= V.unsafeWrite v m
            where step v a x = liftM (min a . succ) $ V.unsafeRead v $ m - x

Результат не сильно отличался от предыдущего. И тогда я решил выжать из этого решения все что можно до последней капли (ну, наверное, все-таки не до последней). Никаких функциональных штучек вроде списков и сверток, никаких min, succ, dropWhile и тому подобного. Берем тупой итеративный алгоритм как в C и применяем его к тупой C-подобной структуре данных mutable unboxed Vector. Список монет реализуем как immutable unboxed Vector — ведь он не изменяется.

import qualified Data.Vector.Unboxed         as VI

-- ----------

minNumCoinsVectorBasic :: VI.Vector Int -> Int -> IO Int
minNumCoinsVectorBasic c m = do
    v <- V.new $ m + 1
    V.unsafeWrite v 0 0
    mapM_ (mnc c v) [1 .. m]
    V.unsafeRead v m
    where mnc c v m = VI.foldM_ (step v) m c
            where step v a x =
                    if x > m
                        then return a
                        else do y <- V.unsafeRead v $ m - x
                                let cur = y + 1
                                if cur < a
                                    then do V.unsafeWrite v m cur
                                            return cur
                                    else return a

Список монет теперь имеет тип VI.Vector Int и должен формироваться перед вызовом minNumCoinsVectorBasic, например так.

          coins    = [24, 18, 14, 11, 10, 8, 5, 3, 1]
          vcoins   = VI.fromList $ sortBy (flip compare) coins

Предельно просто и уныло, бейсик какой-то, но … Вы будете смеяться, но это заработало на порядок быстрее и догнало программу на C++. Репутация haskell была восстановлена! Для измерения производительности я использовал библиотеку criterion. Функция main приняла вид

import           Criterion.Main
import           System.Environment

-- ----------

main = do
    args <- getArgs
    if benchArg `elem` args
        then withArgs (filter (/= benchArg) args) $ defaultMain
             [
                 bench ("map " ++ show m2)          $ whnf
                     (minNumCoinsMap        coins) m2,
                 bench ("array " ++ show m2)        $ whnfIO $
                     minNumCoinsArray       coins  m2,
                 bench ("vector " ++ show m2)       $ whnfIO $
                     minNumCoinsVector      coins  m2,
                 bench ("vector basic " ++ show m2) $ whnfIO $
                     minNumCoinsVectorBasic vcoins m2
             ]
        else {-print $ minNumCoinsPlain [1, 2, 5] 20   -- m = 50 will hang it!-}
             {-print $ minNumCoins_     [1, 2, 5] 20   -- m = 50 will hang it!-}
             {-print $ minNumCoins      [1, 2, 5] 10002-}
             {-print $ minNumCoinsMemo  [1, 2, 5] 10001-}
             {-print $ minNumCoinsMap   [1, 2, 5] 10001-}
             {-print $ minNumCoins      coins m1-}
             {-print $ minNumCoinsMap   coins m1-}
             {-print $ minNumCoinsMap   coins m2-}
             {-minNumCoinsArray         coins m2 >>= print-}
             minNumCoinsVectorBasic   vcoins m2 >>= print
    where benchArg = "-b"
          coins    = [24, 18, 14, 11, 10, 8, 5, 3, 1]
          vcoins   = VI.fromList $ sortBy (flip compare) coins
          m1       = 16659
          m2       = 1665900

То есть, для запуска criterion в эту программу нужно предать опцию -b. Скомпилировав программу с оптимизацией под LLVM,

ghc --make -O2 -fllvm -optlo-O3 minNumCoins.hs

и запустив ее с опцией -b,

./minNumCoins -b -o minNumCoins.html
benchmarking map 1665900
time                 3.304 s    (3.158 s .. 3.478 s)
                     1.000 R²   (0.999 R² .. 1.000 R²)
mean                 3.394 s    (3.370 s .. 3.411 s)
std dev              25.30 ms   (0.0 s .. 29.13 ms)
variance introduced by outliers: 19% (moderately inflated)

benchmarking array 1665900
time                 145.4 ms   (143.3 ms .. 147.7 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 144.1 ms   (143.6 ms .. 145.0 ms)
std dev              890.1 μs   (399.6 μs .. 1.244 ms)
variance introduced by outliers: 12% (moderately inflated)

benchmarking vector 1665900
time                 136.2 ms   (135.4 ms .. 137.5 ms)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 137.7 ms   (137.0 ms .. 139.8 ms)
std dev              1.518 ms   (389.3 μs .. 2.314 ms)
variance introduced by outliers: 11% (moderately inflated)

benchmarking vector basic 1665900
time                 20.87 ms   (20.48 ms .. 21.20 ms)
                     0.996 R²   (0.988 R² .. 0.999 R²)
mean                 21.22 ms   (20.88 ms .. 21.96 ms)
std dev              1.090 ms   (452.6 μs .. 1.646 ms)
variance introduced by outliers: 18% (moderately inflated)

я получил вот такой замечательный графический отчет. Все картинки внутри отчета интерактивные — для их просмотра нужно разрешить javascript. Я не хочу делать никаких выводов, хотя они напрашиваются сами собой. Исходный код программы по-прежнему доступен отсюда.