rpm -Uvh --nodeps.
воскресенье, 21 декабря 2014 г.
Падение xneur в Fedora 21
Все, кого волнует данная тема, могут скачать xneur src.rpm отсюда.
Кстати, лично для меня установка Fedora 21 оказалась наиболее косячной в сравнении с версиями 19 и 20. Хотя косяки эти в общем-то мелкие, но неприятные: кроме падучести xneur, это отсутствие шрифтов infinality, черный экран MATE после пробуждения из длительного сна, произвольное вычищение истории команд оболочки, внезапное подключение репозитория rawhide после обновления, из-за чего я чуть было не обновил все пакеты до следующей Fedora 22 (возможно, это стало следствием не до конца отработавшего скрипта rfremix-upgrade от Russian Fedora), еще что-то там — уже не помню. Надеюсь, все это будет со временем исправлено.
Наиболее жесткий косяк был с установкой Fedora 21 на новый ноутбук MSI GE60 2PL Apache. В нем установлены две видеокарты — встроенная Intel и более мощная NVIDIA GeForce GTX 850M. Переключение карт в BIOS не предусмотрено. При первой же загрузке подключение драйвера Nouveau сломало напрочь весь процесс с самого начала. Дальнейшие перезагрузки тоже приводили к разнообразным ошибкам памяти в самом начале загрузки. Пару раз мне все-таки удалось загрузиться (сказалась неопределенность поведения), и в один из этих светлых моментов я установил драйвер bumblebee: это исправило проблему, и даже optirun теперь переключает видеокарты как следует.
Update. Собрал пропатченный билд xneur на федориной копре. Репозиторий здесь. Установить пакеты с подключенным репозиторием через yum не удастся, если в системе установлены графические оболочки типа gxneur, которые зависят от оригинального пакета xneur. Поэтому можно загрузить подходящие rpm пакеты непосредственно отсюда и установить их командой
пятница, 28 ноября 2014 г.
Перезапуск падучего приложения из самого приложения (спасение утопающих ...)
Представьте себе ситуацию, когда у вас есть некоторое серверное приложение, которое потенциально может падать (я не буду рассуждать о причинах: допустим, это новое, написанное вами приложение, которое вы хотите считать абсолютно надежным, но это не так). Предположим, что вам ни в коем случае нельзя допускать простоя службы, реализуемой этим приложением. Есть разные способы следить за работой программ и перезапускать их в случае надобности извне (например, cron). Но …
Спасение утопающих — дело рук самих утопающих. Именно этот подход я и хочу продемонстрировать в самых общих чертах. Разумеется, вы должны обладать исходным кодом программы. В этом случае достаточно выделить исходный полезный код в отдельный дочерний процесс (worker), а исходный процесс превратить в надсмотрщика (master или supervisor), который будет перезапускать полезный дочерний процесс в случае его падения.
Вот код, а пояснения ниже.
#include <unistd.h> #include <sys/wait.h> #include <string.h> #include <iostream> #ifdef MAXCYCLES #define LOOPSTOPCOND i < MAXCYCLES #else #define LOOPSTOPCOND #endif int main( int argc, char ** argv ) { int pid( 0 ); std::cout << "Master: " << getpid() << std::endl; for ( int i( 0 ); LOOPSTOPCOND; ++i ) { if ( ( pid = fork() ) == 0 ) /* Worker process */ { std::cout << "(cycle " << i << ") Worker: " << getpid() << std::endl; /* do segv if option -s was specified in command line */ if ( argc > 1 && ! strcmp( argv[ 1 ], "-s" ) ) int a( *( ( int * )0 ) ); break; } else /* Master process */ { if ( pid < 0 ) { std::cout << "Failed to fork a worker process, exiting" << std::endl; break; } int status( 0 ); bool respawn( false ); waitpid( pid, &status, 0 ); if ( ! WIFSIGNALED( status ) ) break; int sig( WTERMSIG( status ) ); std::cout << "Worker process was signaled " << sig << std::endl; switch ( sig ) { case SIGSEGV: respawn = true; default: break; } if ( ! respawn ) break; } } return 0; }Прежде всего нужно отметить, что код написан на C++, хотя ничто не мешает переписать его на C, если понадобится. Макросы MAXCYCLES и LOOPSTOPCOND нужны здесь только в демонстрационных целях: дабы ограничить число перезапусков воркеров величиной, заданной во время компиляции (MAXCYCLES). В релизной версии приложения цикл for скорее всего примет вид
for ( ; ; ), так что эти макросы окажутся не нужны.
Итак, цикл for нужен для перезапуска воркер-процессов в случае их падения. Внутри цикла for с помощью вызова fork() производится расщепление основного процесса, связанного с приложением, на две части — новый процесс (воркер) и старый (мастер). Как известно, fork() возвращает 0 в новом процессе, и какое-либо значение большее нуля — в старом. В новом процессе (блок после if ( ( pid = fork() ) == 0 )) мы выводим сообщение о старте воркера, вызываем его полезный код (в этом примере полезный код присутствует только в случае, когда в программу была передана опция командной строки -s: он приводит к посылке ядром сигнала SIGSEGV из-за разыменования нулевого указателя), и в конце обязательно ставим break;, поскольку воркер по-прежнему находится в цикле, и без его прерывания продолжит запускать новые воркеры, притворившись мастером!
Внутри кода, относящегося к мастер-процессу (блок после else), прежде всего проверяется, что вызов функции fork() завершился успешно. Затем мастер ожидает завершения воркера с помощью вызова waitpid() и проверяет, каким образом тот завершился. Макрос WIFSIGNALED() позволяет установить, был ли воркер остановлен сигналом ядра или вышел самостоятельно. Во втором случае цикл for прерывается и мастер завершает свою работу. В случае, если воркер был остановлен сигналом, мы хотим узнать каким именно. Для этого предназначен макрос WTERMSIG(). Если это был сигнал SIGSEGV, то локальной переменной respawn присваивается истинное значение. Если ее значение остается ложным (когда воркер был прерван любым другим сигналом), то цикл завершается.
Давайте посмотрим, как это работает (исходный файл я назвал main.cpp). Компилируем.
g++ -g -DMAXCYCLES=4 -o test main.cpp
Я ограничил число перезапусков четырьмя. Запускаем программу с нормальным завершением воркера.
./test Master: 6061 (cycle 0) Worker: 6062Завершился воркер — завершился мастер. А теперь запустим программу с падучим воркером.
./test -s Master: 6085 (cycle 0) Worker: 6086 Worker process was signaled 11 (cycle 1) Worker: 6152 Worker process was signaled 11 (cycle 2) Worker: 6154 Worker process was signaled 11 (cycle 3) Worker: 6156 Worker process was signaled 11Мастер-процесс четыре раза перезапустил упавшие воркер-процессы, как и ожидалось. Я не стал касаться вопросов наследования ресурсов мастера при инициализации воркера и корректного завершения воркера при получении мастером сигнала прерывания — это отдельные интересные темы.
понедельник, 10 ноября 2014 г.
Ручное распараллеливание вычислений в haskell
Возьмем самую быструю реализацию задачи из предыдущей статьи.
Вот такая картинка появится на экране. Что ж, здесь достаточно много информации как на графике в центральном окне, так и во вкладках, расположенных в нижнем окне. Верхняя часть графика — Activity — отражает степень загрузки потоков выполнения (у нас их четыре) в течении всего времени жизни программы. Четыре плотных зеленых слоя, начинающиеся на 25 мс означают, что все четыре потока работают на полную катушку: это очень хорошо. Начиная с 0.13 с программа стала работать в одном потоке: это финальная сортировка по длинам в функции grpSeqsPar (вы можете убедиться в этом, удалив
Ага, переполнение кучи. Опция
На этот раз продуктивность программы составила 99% против 93% в предыдущем запуске, сборщик мусора практически не задействован. Однако время выполнения почти не изменилось. Помните, я говорил, что это странно. Так вот, никакой странности здесь оказывается нет: просто время сборки мусора эффективно заменилось временем, затраченным на аллокацию и управление гораздо бо́льшими объемами памяти. В данном случае значительное время было потрачено в фазе инициализации (INIT): это хорошо видно на графике.
А что будет, если использовать параметры аллокации кучи для сборщика мусора по умолчанию?
Интересный результат. Время выполнения немного выросло. Время сборщика мусора (GC time) теперь равно времени полезной нагрузки (Mutator time), и это видно по соотношению зеленого и оранжевого цветов в потоках. Загрузка просела до двух (чуть позже трех) потоков одновременно. Несмотря на это, время выполнения, как я уже сказал, выросло ненамного: видимо это связано с более эффективной стратегией аллокации и меньшим потреблением памяти.
Update 2. Функцию concat’ из grpSeqsPar можно переписать с помощью стандартной функции ceiling.
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% от всего времени выполнения программы, за ней идет функция nearSeqs — 20.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.
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 покажет такую картинку.
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
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 не поддерживаются.
среда, 5 ноября 2014 г.
Решение задачи о поиске наиболее часто встречающейся подпоследовательности с заданным максимальным числом мутаций на haskell
Эта задача имеет большое значение в биологии и часто формулируется как поиск наиболее часто встречающегося k-мера с максимально допустимым числом мутаций n внутри заданной генетической последовательности. Генетическая последовательность обычно моделируется строкой произвольной длины, составленной из символов некоторого алфавита. Соответственно, k-мер — это подстрока длиной k внутри генетической последовательности. В постановке задачи искомые k-меры могут накладываться друг на друга. Так, если общая длина исходной последовательности равна 10, то в ней определены один 10-мер, два 9-мера, три 8-мера и т.д. Благодаря условию о допустимых мутациях искомая наиболее часто повторяющаяся подпоследовательность может вообще ни разу не встретиться в исходной последовательности!
Для решения задачи предложим такой алгоритм: для каждого k-мера в заданной последовательности находим все его мутации, включая исходное значение, складываем их вместе и находим наиболее часто встречающийся элемент: он и будет ответом на поставленную задачу.
Я приведу решение, а ниже его прокомментирую.
import qualified Data.List as L 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 = map head . L.group . L.sort . 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 $ nearSeqs n main = do let seqs = ["ACGTTGCATGTCGCATGATGCATGAGAGCT", "CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGC\ \GGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCC\ \GGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCG\ \CCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTA\ \CCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACAC\ \ACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGG\ \GCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC"] {-print $ nearSeqs 2 "TTTT"-} {-mapM_ print' $ take 10 $ grpSeqs 1 $ subSeqs 4 $ head seqs-} mapM_ print' $ take 10 $ grpSeqs 2 $ subSeqs 10 $ seqs !! 1 where print' = print . (head &&& length)В функции main происходит тестирование алгоритма. Функция alphabet задает алфавит последовательности. Как видим, он состоит всего из четырех символов. Я намеренно задал тип alphabet как
[Char], а не String, дабы подчеркнуть, что это алфавит, а не простая строка. Для простоты я не стал каким-либо образом верифицировать исходные последовательности в seqs, объявленной в main. Функция subSeqs возвращает все k-меры из заданной последовательности (ее первый формальный параметр n соответствует числу k).
Функция nearSeqs возвращает список уникальных мутаций (включая исходное значение) некоторой последовательности, ее формальный параметр n задает максимально допустимое количество мутаций. Реализация этой функции основана на рекурсивном применении генератора списков по всевозможным перестановкам общего числа n позиций внутри заданной последовательности. Для каждой перестановки генерируются всевозможные мутации на основе алфавита alphabet: при таком подходе результирующий список может содержать повторяющиеся элементы, поэтому он фильтруется функцией map head . L.group . L.sort. Как выяснилось впоследствии, правильный выбор этой функции очень сильно влияет на скорость программы. Так, в моем оригинальном варианте L.nub . L.sort программа работала в два раза медленнее! Просто L.nub работает еще медленнее (все-таки это функция класса O(n^2)). Текущий вариант я подсмотрел здесь. Позже я нашел самый быстрый вариант: Data.List.Ordered.nubSort, но все же решил на нем не останавливаться, поскольку он требует импорта еще одного модуля.
К слову, функцию nearSeqs можно определить без необходимости финальной сортировки списка, заранее вынимая те элементы из alpahabet, которые приведут к дубликатам.
nearSeqs'' :: Int -> String -> [String] nearSeqs'' n s = s : concatMap (`nearSeqs'` s) [1 .. n] where nearSeqs' 0 s = [s] nearSeqs' n s = concatMap (\(a, b) -> map (a++) b) $ concatMap (\m -> [(z ++ [x], nearSeqs' (n - 1) $ tail z') | let (z, z') = splitAt m s, x <- L.delete (head z') alphabet]) [0 .. length s - 1]Однако, эта функция работает медленнее, поэтому я ее забраковал. Чтобы посмотреть на функцию nearSeqs в деле, раскомментируйте первую закомментированную строку в main, скомпилируйте программу и запустите ее. Функция grpSeqs вызывает nearSeqs, сортирует и группирует возвращенный ею список, и далее сортирует сгруппированный список по длине его подсписков в порядке убывания. В функции main из этого списка берется 10 первых подсписков и выводятся их первый элемент (все элементы подсписков должны совпадать) и длина. Для компиляции всех примеров в этой статье я использовал флаг -O2, имя исходника — gen.hs. Вывод программы с замером времени выполнения:
time ./gen ("GCACACAGAC",20) ("GCGCACACAC",20) ("ACACACACAC",19) ("CCCGCACACA",19) ("CACACACACG",18) ("CGCGCACACA",18) ("CGGCACACAC",18) ("GCACACACCC",18) ("GCACACACGC",18) ("GCCCGCACAC",18) real 0m0.639s user 0m0.612s sys 0m0.022sЧуть больше, чем полсекунды. Неплохо. Однако, этого я достиг не сразу и поэтому перепробовал множество подходов. В частности, я пробовал использовать для сборки результатов nearSeqs хэш-таблицу Data.HashTable.IO. В список импортируемых модулей добавляем
import qualified Data.HashTable.IO as H import Data.Maybe (fromMaybe)В определения записываем
type HashTable k v = H.CuckooHashTable k v grpSeqsHt :: Int -> Int -> [String] -> IO (HashTable String Int) grpSeqsHt m n s = do h <- H.newSized m let insert' h k = H.lookup h k >>= H.insert h k . (+1) . fromMaybe 0 mapM_ (mapM_ (insert' h) . nearSeqs n) s return h mostOftenSeqsHt :: HashTable String Int -> IO [(String, Int)] mostOftenSeqsHt = H.foldM findMax [] where findMax [] a = return [a] findMax c@((_, v') : _) a@(_, v) = return $ case v `compare` v' of GT -> [a]; EQ -> a : c; LT -> cВ функцию main записываем строку
grpSeqsHt 70000 2 (subSeqs 10 $ seqs !! 1) >>= mostOftenSeqsHt >>= printФункция grpSeqsHt создает хэш-таблицу исходного размера m (величина 70000 подобрана эмпирическим путем: примерно столько уникальных k-меров генерируется из
seqs !! 1). Ключами таблицы являются значения k-меров, значениями — их общее количество в исходных данных. Если при добавлении новой ячейки ее еще не существует в таблице (это проверяется функцией H.lookup), то она создается со значением 1, иначе соответствующее ей значение увеличивается на единицу. Функция mostOftenSeqsHt проходится по уже заполненной таблице и возвращает список наиболее часто встречающихся k-меров (заметьте отличие: первый вариант решения задачи формировал список всех k-меров в порядке уменьшения их частоты в исходной последовательности).
Второй вариант с использованием хэш-таблицы учитывает тот факт, что вычисление списка наиболее часто встречающихся k-меров можно произвести прямо во время заполнения таблицы внутри функции grpSeqsHt, и функция mostOftenSeqsHt становится ненужной.
Итак, добавляем в список импортируемых модулей
import Control.Monad (foldM), в определения
grpSeqsHt' :: Int -> Int -> [String] -> IO (HashTable String Int, [(String, Int)]) grpSeqsHt' m n s = do h <- H.newSized m let insert' h k = H.lookup h k >>= (\v -> H.insert h k v >> return (k, v)) . (+1) . fromMaybe 0 findMax [] a = return [a] findMax c@((_, v') : _) a@(_, v) = return $ case v `compare` v' of GT -> [a]; EQ -> a : c; LT -> c c <- foldM (\c s -> foldM (\c k -> insert' h k >>= findMax c) c $ nearSeqs n s) [] s return (h, c), и в функцию main
grpSeqsHt' 70000 2 (subSeqs 10 $ seqs !! 1) >>= print . sndВ этом случае использование хэш-таблицы становится побочным артефактом реализации, поскольку за пределами функции grpSeqsHt’ она не используется, несмотря на то, что grpSeqsHt’ возвращает ее как первый элемент кортежа. Как ни странно, варианты с хэш-таблицей оказались медленнее. Так, первый вариант выдал
time ./gen [("GCACACAGAC",20),("GCGCACACAC",20)] real 0m0.656s user 0m0.643s sys 0m0.011s, а второй
time ./gen [("GCGCACACAC",20),("GCACACAGAC",20)] real 0m0.676s user 0m0.664s sys 0m0.009sНу и, собственно, чемпион! Возьмем исходное решение и чуточку изменим алгоритм. Предположим, что список, возвращенный из функции subSeqs, может содержать дубликаты (как оказалось,
seqs !! 1 действительно содержит 98 дубликатов k-меров при их общем числе 352). Для дубликатов не нужно получать список мутаций каждый раз заново: достаточно получить список единожды и реплицировать его по числу дубликатов. Давайте просто перепишем функцию grpSeqs в предположении, что она по-прежнему будет получать данные из subSeqs.
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Вывод программы:
time ./gen ("GCACACAGAC",20) ("GCGCACACAC",20) ("ACACACACAC",19) ("CCCGCACACA",19) ("CACACACACG",18) ("CGCGCACACA",18) ("CGGCACACAC",18) ("GCACACACCC",18) ("GCACACACGC",18) ("GCCCGCACAC",18) real 0m0.546s user 0m0.525s sys 0m0.018s
Подписаться на:
Сообщения (Atom)