понедельник, 10 ноября 2014 г.

Ручное распараллеливание вычислений в haskell

Возьмем самую быструю реализацию задачи из предыдущей статьи.
import qualified Data.List as L
import qualified Data.List.Ordered as LO
import Data.Ord (comparing)
import Control.Arrow

alphabet :: [Char]
alphabet = "ACGT"

subSeqs :: Int -> String -> [String]
subSeqs n = takeWhile ((==n) . length) . map (take n) . L.tails

nearSeqs :: Int -> String -> [String]
nearSeqs n = LO.nubSort . nearSeqs' n
    where nearSeqs' 0 s = [s]
          nearSeqs' n s =
            concatMap (\(a, b) -> map (a++) b) $
                concatMap (\m -> [(z ++ [x], nearSeqs' (n - 1) $ tail z') |
                                  x <- alphabet, let (z, z') = splitAt m s])
                [0 .. length s - 1]

grpSeqs :: Int -> [String] -> [[String]]
grpSeqs n = L.sortBy (flip $ comparing length) . L.group . L.sort . seqs
    where seqs = concatMap (\(s, m) -> concat . replicate m . nearSeqs n $ s) .
                    map (head &&& length) . L.group . L.sort

main = do
    let seqs = ["ACGTTGCATGTCGCATGATGCATGAGAGCT",
                "CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGC\
                \GGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCC\
                \GGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCG\
                \CCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTA\
                \CCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACAC\
                \ACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGG\
                \GCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC"]
    mapM_ print' $ take 10 $ grpSeqs 2 $ subSeqs 10 $ seqs !! 1
    where print' = print . (head &&& length)
Я все же решил использовать Data.List.Ordered.nubSort для удаления дубликатов в nearSeqs, поэтому программа должна работать еще чуточку быстрее. В среднем новая версия работает порядка 0.51 секунды, то есть действительно быстрее на ∼0.03 секунды. Отлично. Давайте отпрофилируем нашу программу и посмотрим, что можно сделать, чтобы она работала еще быстрее. Компилируем с поддержкой профилирования.
ghc --make -rtsopts -prof -auto-all -caf-all -fforce-recomp gen
Запускаем программу с записью данных профилирования в файл.
./gen +RTS -p
Смотрим результаты профилирования в файле gen.prof.
    Mon Nov 10 16:38 2014 Time and Allocation Profiling Report  (Final)

       gen +RTS -p -RTS

    total time  =        0.77 secs   (766 ticks @ 1000 us, 1 processor)
    total alloc = 593,276,296 bytes  (excludes profiling overheads)

COST CENTRE                MODULE  %time %alloc

grpSeqs                    Main     52.0   24.6
nearSeqs                   Main     20.6   10.7
nearSeqs.nearSeqs'.\       Main      8.4   19.9
nearSeqs.nearSeqs'.\       Main      7.8   19.5
nearSeqs.nearSeqs'         Main      4.6    9.9
nearSeqs.nearSeqs'.\.(...) Main      4.3   12.3
grpSeqs.seqs               Main      0.9    1.5
grpSeqs.seqs.\             Main      0.4    1.5


                                                                              individual     inherited
COST CENTRE                         MODULE                  no.     entries  %time %alloc   %time %alloc

MAIN                                MAIN                     51           0    0.0    0.0   100.0  100.0
 CAF:main                           :Main                   101           0    0.0    0.0     0.0    0.0
  main                              Main                    116           0    0.0    0.0     0.0    0.0
   main.print'                      Main                    118           0    0.0    0.0     0.0    0.0
 CAF:main                           Main                    100           0    0.0    0.0   100.0  100.0
  main                              Main                    102           1    0.0    0.0   100.0  100.0
   main.print'                      Main                    117           1    0.0    0.0     0.0    0.0
   main.seqs                        Main                    106           1    0.0    0.0     0.0    0.0
   subSeqs                          Main                    105           1    0.0    0.0     0.0    0.0
   grpSeqs                          Main                    103           1   52.0   24.6   100.0   99.9
    grpSeqs.seqs                    Main                    104           1    0.9    1.5    48.0   75.3
     grpSeqs.seqs.\                 Main                    107         254    0.4    1.5    47.1   73.8
      nearSeqs                      Main                    108         254   20.6   10.7    46.7   72.3
       nearSeqs.nearSeqs'           Main                    109      193294    4.6    9.9    26.1   61.7
        nearSeqs.nearSeqs'.\        Main                    112      193040    8.4   19.9     8.4   19.9
        nearSeqs.nearSeqs'.\        Main                    110       48260    7.8   19.5    13.2   31.8
         nearSeqs.nearSeqs'.\.z     Main                    115      192024    0.3    0.0     0.3    0.0
         nearSeqs.nearSeqs'.\.(...) Main                    114      193040    4.3   12.3     4.3   12.3
         nearSeqs.nearSeqs'.\.z'    Main                    113      183951    0.8    0.0     0.8    0.0
 CAF:$dShow_r1F9                    Main                     99           0    0.0    0.0     0.0    0.0
 CAF:$dOrd_r1F8                     Main                     98           0    0.0    0.0     0.0    0.0
 CAF:$dEq_r1D9                      Main                     97           0    0.0    0.0     0.0    0.0
 CAF:alphabet_ryE                   Main                     96           0    0.0    0.0     0.0    0.0
  alphabet                          Main                    111           1    0.0    0.0     0.0    0.0
 CAF                                GHC.Conc.Signal          92           0    0.0    0.0     0.0    0.0
 CAF                                GHC.IO.Handle.FD         84           0    0.0    0.0     0.0    0.0
 CAF                                GHC.IO.Encoding          79           0    0.0    0.0     0.0    0.0
 CAF                                GHC.IO.Encoding.Iconv    68           0    0.0    0.0     0.0    0.0
Видим, что абсолютным чемпионом по времязатратам является функция grpSeqs с результатом 52% от всего времени выполнения программы, за ней идет функция nearSeqs20.6%. Функция grpSeqs состоит сплошь из сортировок и группировок списков — эти типы вычислений не имеют значительных перспектив для распараллеливания. С другой стороны, списки, возвращаемые функцией nearSeqs для 254 уникальных подпоследовательностей исходной последовательности seqs !! 1, могут быть вычислены независимо, а это может означать потенциально ожидаемый 15% выигрыш по времени для 4-ядерной машины (20.6% - (20.6% / 4) ≈ 15%). Перепишем функцию grpSeqs с параллельным вычислением nearSeqs. В список импортируемых модулей добавляем
import Control.Parallel.Stategies
Определяем новую функцию grpSeqsPar.
grpSeqsPar :: Int -> [String] -> [[String]]
grpSeqsPar n = L.sortBy (flip $ comparing length) . L.group . L.sort . seqs
    where seqs = concat . parMap rdeepseq
                    (\(s, m) -> concat . replicate m . nearSeqs n $ s) .
                        map (head &&& length) . L.group . L.sort
Не забываем заменить вызов grpSeqs на grpSeqsPar в функции main.
    mapM_ print' $ take 10 $ grpSeqsPar 2 $ subSeqs 10 $ seqs !! 1
Компилируем с поддержкой распараллеливания.
ghc --make -O2 -rtsopts -threaded -feager-blackholing -fforce-recomp gen
На моей машине установлен 4-ядерный процессор Intel Core i5. Запускаем с распараллеливанием на 4 ядра.
time ./gen +RTS -N4
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)

real    0m0.521s
user    0m1.049s
sys     0m0.339s
Что-то никакого толку. Давайте запустим программу с опцией -s, которая выводит разную полезную информацию на экран.
time ./gen +RTS -N4 -s
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)
     316,334,968 bytes allocated in the heap
     161,177,512 bytes copied during GC
      18,120,968 bytes maximum residency (10 sample(s))
         310,528 bytes maximum slop
              53 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0       365 colls,   365 par    0.54s    0.14s     0.0004s    0.0147s
  Gen  1        10 colls,     9 par    0.30s    0.08s     0.0084s    0.0159s

  Parallel GC work balance: 38.30% (serial 0%, perfect 100%)

  TASKS: 10 (1 bound, 9 peak workers (9 total), using -N4)

  SPARKS: 254 (254 converted, 0 overflowed, 0 dud, 0 GC'd, 0 fizzled)

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.60s  (  0.30s elapsed)
  GC      time    0.84s  (  0.23s elapsed)
  EXIT    time    0.00s  (  0.00s elapsed)
  Total   time    1.44s  (  0.53s elapsed)

  Alloc rate    529,187,979 bytes per MUT second

  Productivity  41.9% of total user, 113.3% of total elapsed

gc_alloc_block_sync: 67687
whitehole_spin: 0
gen[0].sync: 22
gen[1].sync: 10747

real    0m0.536s
user    0m1.025s
sys     0m0.417s
Продуктивность программы составила 41.9%, остальное время ушло на сборку мусора (GC time 0.84s)! Обычно с этим борются увеличением предполагаемого размера кучи (suggested heap size) для сборщика мусора (см. документацию здесь). Давайте установим эту величину в 1 гигабайт.
time ./gen +RTS -N4 -H1G -s
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)
     315,385,264 bytes allocated in the heap
      17,826,760 bytes copied during GC
      15,081,784 bytes maximum residency (2 sample(s))
         118,472 bytes maximum slop
             970 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0         1 colls,     1 par    0.03s    0.01s     0.0110s    0.0110s
  Gen  1         2 colls,     1 par    0.12s    0.04s     0.0178s    0.0321s

  Parallel GC work balance: 83.58% (serial 0%, perfect 100%)

  TASKS: 10 (1 bound, 9 peak workers (9 total), using -N4)

  SPARKS: 254 (254 converted, 0 overflowed, 0 dud, 0 GC'd, 0 fizzled)

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.56s  (  0.31s elapsed)
  GC      time    0.15s  (  0.05s elapsed)
  EXIT    time    0.00s  (  0.00s elapsed)
  Total   time    0.72s  (  0.37s elapsed)

  Alloc rate    567,298,189 bytes per MUT second

  Productivity  78.5% of total user, 153.2% of total elapsed

gc_alloc_block_sync: 15012
whitehole_spin: 0
gen[0].sync: 2831
gen[1].sync: 0

real    0m0.391s
user    0m0.600s
sys     0m0.137s
А вот это уже здорово! Время выполнения программы сократилось до 0.39 секунды, продуктивность увеличилась до 78.5%. Кстати, в выводе программы есть строка, озаглавленная SPARKS: в ней находится статистика запусков параллельных вычислений, которые называются искрами (sparks). Всего было запущено 254 искры, что соответствует числу уникальных 10-меров в исходной последовательности. А что будет, если в оригинальной, нераспараллеленной версии выделить 1 гигабайт для кучи? Может быть она тоже ускорится? Как показала практика — нет. Продуктивность программы действительно возрастает, но время выполнения практически не изменяется, как это ни странно. И, наконец, следующая версия функции grpSeqsPar работает еще чуточку быстрее (но незначительно).
grpSeqsPar :: Int -> Int -> [String] -> [[String]]
grpSeqsPar m n = L.sortBy (flip $ comparing length) . L.group . sortPar . seqs
    where sortPar = foldr LO.merge [] . parMap rdeepseq L.sort
          seqs = concat' m . parMap rdeepseq
                    (\(s, m) -> concat . replicate m . nearSeqs n $ s) .
                        map (head &&& length) . L.group . L.sort
          concat' m s = let l = length s
                        in map concat $ chunksOf
                            (l `div` m + if l `mod` m == 0 then 0 else 1) s
Вот только я совсем не уверен в ее оптимальности (именно в оптимальности, а не в корректности): в самом деле, вложенные стратегии rdeepseq скорее всего не самый лучший вариант. В этой версии затратная сортировка списка, возвращаемого функцией seqs, заменена на серию последовательных слияний LO.merge его по отдельности отсортированных частей. Количество частей, на которые разбивается исходный список, задается новым параметром m функции grpSeqsPar. На мой взгляд, он должен быть равен количеству ядер процессора. Соответственно, вызов в функции main принимает вид
    mapM_ print' $ take 10 $ grpSeqsPar 4 2 $ subSeqs 10 $ seqs !! 1
Для использования функции chunksOf в список импортируемых модулей следует добавить
import Data.List.Split (chunksOf)
Если скомпилировать новую программу и запустить ее на выполнение, то в строке SPARKS появятся несколько десятков погасших искр (fizzled sparks). Искры гаснут в том случае, когда танки (thunks), то есть элементы для вычисления, на которые они ссылаются, уже были вычислены в некотором параллельном потоке программы. Такая ситуация может свидетельствовать о неоптимальном выборе стратегии параллелизации. Опция компиляции -feager-blackholing позволяет предотвратить многократные вычисления одних и тех же танков, поэтому само по себе появление погасших искр не должно влиять на производительность программы. Что еще можно распараллелить в этой программе? Ну, например, функцию subSeqs. Хотя в данном случае это не даст никаких преимуществ, поскольку время выполнения этой функции по результатам профилирования составляет 0.0%. Тем не менее, для полноты картины я приведу параллельную версию subSeqs.
subSeqsPar :: Int -> Int -> String -> [String]
subSeqsPar m n s = subSeqs n s `using` parListChunk (length s `div` m) rdeepseq
Здесь я решил применить parListChunk со стратегией rdeepseq. Размер чанка вычисляется на основании значения первого параметра m функции subSeqsPar путем разделения исходной последовательности s на m равных частей: это упрощенный аналог алгоритма из функции concat’. Update. Вложенная стратегия parMap rdeepseq для вычисления всех мутаций с помощью nearSeqs в функции grpSeqsPar действительно не нужна: внешняя стратегия из функции sortPar должна прекрасно справиться с распараллеливанием сама. По этой же самой причине subSeqsPar тоже не нужна в принципе. Таким образом, функция grpSeqsPar переписывается так:
grpSeqsPar :: Int -> Int -> [String] -> [[String]]
grpSeqsPar m n = L.sortBy (flip $ comparing length) . L.group . sortPar . seqs
    where sortPar = foldr LO.merge [] . parMap rdeepseq L.sort
          seqs = concat' m .
                    map (\(s, m) -> concat . replicate m . nearSeqs n $ s) .
                        map (head &&& length) . L.group . L.sort
          concat' m s = let l = length s
                        in map concat $ chunksOf
                            (l `div` m + if l `mod` m == 0 then 0 else 1) s
Новая стратегия должна сократить число искр до четырех и потенциально ускорить выполнение программы. Кстати, для визуализации выполнения параллельной программы на haskell существует прекраснейший инструмент — программа threadscope. Для того, чтобы воспользоваться threadscope, нужно скомпилировать нашу программу с поддержкой журнала событий.
ghc --make -O2 -rtsopts -threaded -eventlog -feager-blackholing -fforce-recomp gen
Давайте запустим программу с опцией -H1G, как мы это делали раньше.
time ./gen +RTS -N4 -H1G -qa -lf
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)

real    0m0.282s
user    0m0.525s
sys     0m0.118s
Не обращайте внимание на новый рекорд скорости (0.282 сек.): просто я запустил программу на более быстрой машине. Зато обратите внимание на новые опции: -qa и -lf. Первая опция экспериментальная и в документации haskell говорится, что она может ускорить выполнение параллельной программы, а может и не ускорить. Для меня она действительно работает. Вторая опция приводит к записи событий в файл gen.eventlog. Давайте откроем этот файл в программе threadscope. Вот такая картинка появится на экране. Что ж, здесь достаточно много информации как на графике в центральном окне, так и во вкладках, расположенных в нижнем окне. Верхняя часть графика — Activity — отражает степень загрузки потоков выполнения (у нас их четыре) в течении всего времени жизни программы. Четыре плотных зеленых слоя, начинающиеся на 25 мс означают, что все четыре потока работают на полную катушку: это очень хорошо. Начиная с 0.13 с программа стала работать в одном потоке: это финальная сортировка по длинам в функции grpSeqsPar (вы можете убедиться в этом, удалив L.sortBy (flip $ comparing length) . из определения grpSeqsPar и запустив программу заново). А что же это за провал, длившийся 20 мс в самом начале программы, когда активность была нулевой и все потоки занимались сборкой мусора (оранжевые полосы)? К счастью, график является интерактивным, и мы можем узнать, что произошло в определенной точке выполнения программы, кликнув мышкой по графику в этой точке и просмотрев события во вкладке Raw events. Ага, переполнение кучи. Опция -H1G не влияет на первый запуск сборщика мусора: для этого нужно настроить опцию -A. Давайте запустим программу с этой опцией.
time ./gen +RTS -N4 -A128M -qa -lf
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)

real    0m0.270s
user    0m0.471s
sys     0m0.115s
Я подобрал размер аллокации (128M) таким образом, чтобы единожды выделенной кучи хватило на все время выполнения программы. Теперь threadscope покажет такую картинку. На этот раз продуктивность программы составила 99% против 93% в предыдущем запуске, сборщик мусора практически не задействован. Однако время выполнения почти не изменилось. Помните, я говорил, что это странно. Так вот, никакой странности здесь оказывается нет: просто время сборки мусора эффективно заменилось временем, затраченным на аллокацию и управление гораздо бо́льшими объемами памяти. В данном случае значительное время было потрачено в фазе инициализации (INIT): это хорошо видно на графике. А что будет, если использовать параметры аллокации кучи для сборщика мусора по умолчанию?
time ./gen +RTS -N4 -qa -lf
("GCACACAGAC",20)
("GCGCACACAC",20)
("ACACACACAC",19)
("CCCGCACACA",19)
("CACACACACG",18)
("CGCGCACACA",18)
("CGGCACACAC",18)
("GCACACACCC",18)
("GCACACACGC",18)
("GCCCGCACAC",18)

real    0m0.312s
user    0m0.877s
sys     0m0.169s
Интересный результат. Время выполнения немного выросло. Время сборщика мусора (GC time) теперь равно времени полезной нагрузки (Mutator time), и это видно по соотношению зеленого и оранжевого цветов в потоках. Загрузка просела до двух (чуть позже трех) потоков одновременно. Несмотря на это, время выполнения, как я уже сказал, выросло ненамного: видимо это связано с более эффективной стратегией аллокации и меньшим потреблением памяти. Update 2. Функцию concat’ из grpSeqsPar можно переписать с помощью стандартной функции ceiling.
          concat' m s = let l = let (/.) = (/) `on` fromIntegral
                                in ceiling $ length s /. m
                        in map concat $ chunksOf l s
Для функции on нужен
import Data.Function (on)
Определение функции (/.) я взял отсюда. Она нужна для преобразования типов (или коэрции) аргументов операции вещественного деления (/) из Int в абстрактный тип Num a => a. Без этого преобразования типов программа не скомпилируется: неявные преобразования типов в стиле C и C++ в haskell не поддерживаются.

Комментариев нет:

Отправить комментарий