Привет! Я Рома, бэкендер-питонист в KTS.
Это вторая статья в моем цикле об алгоритме архивации bzip2. Первую можно прочитать здесь, но для понимания сегодняшней темы она необязательна. Ниже я разберу преобразование Барроуза-Уилера — ключевой этап сжатия bzip2.
Немного контекста: в своем цикле статей я разбираю механизм работы алгоритма архивации bzip2 с примерами на Python. Всю предысторию такого сомнительного решения я объясняю в своей первой публикации. Вкратце — цикл преследует образовательную цель. Он вырос из небольшого пет-проекта, в котором я писал собственный bzip2-архиватор на Python (медленный и бинарно несовместимый со стандартным bzip2, но работающий).
Должен предупредить, что эта статья — весьма увесистый кусок «научного гранита». Вполне вероятно, что на практике эти знания ни мне, ни большинству из читателей в жизни не пригодятся. И все же я приложил некоторые усилия, чтобы разбить этот гранит на фракции удобоваримого размера.
Моя идея — сделать так, чтобы читатель, осиливший весь цикл, мог написать bzip2-архиватор на любом языке программирования (которым он владеет), пользуясь только информацией из моих статей.
Еще одно важное замечание: в имплементации намеренно не используется ничего кроме стандартной библиотеки Python3. Считаю, что такой подход увеличивает образовательную ценность статьи.
Оглавление
Обзор алгоритма bzip2
Словом «bzip2» в разных контекстах может обозначаться как алгоритм, так и реализующая его одноименная программа. Здесь и далее (если не указано иное) под словом «bzip2» мы будем подразумевать именно алгоритм.
Для простоты введем упрощенное определение процесса архивации. Архивация — это преобразование данных, которое характеризуется двумя свойствами:
сжатие: размер данных на выходе как правило меньше, чем размер данных на входе;
обратимость: можно «разжать» то, что было сжато, и получить в точности исходные данные.
Алгоритм bzip2 состоит из цепочки преобразований. Каждое из них обратимо, и некоторые из них сжимают:
RLE -> BWT -> MTF -> RLE -> HFC
При архивации большой файл разбивается на относительно небольшие блоки фиксированной длины (условно, 100 KiB), и каждый такой блок проходит через вышеописанный конвейер трансформаций слева направо. При распаковке (деархивации) данные будут идти по узлам конвейера в обратном порядке. Это возможно благодаря обратимости всех преобразований в конвейере.
Преобразования RLE (Run-Length Encoding) и HFC (Huffman coding) отвечают за сжатие информации. Преобразования BWT (Burrows-Wheeler transform) и MTF (Move-to-front transform) выполняют вспомогательную функцию: не сжимают данные, но реструктурируют их таким образом, что другим алгоритмам становится удобнее с ними работать.
Предыдущая статья была посвящена алгоритму RLE, который встречается в цепочке целых два раза. RLE (aka Run-Length Encoding aka кодирование длин серий) — это алгоритм, который строку "goooogle" превращает в "g4ogle" (на самом деле, всё гораздо сложнее и интереснее, но в первом приближении он работает именно так). RLE в bzip2 берет на себя львиную долю сжатия, но он хорошо сжимает только специальном образом подготовленные данные.
BWT расшифровывается как Burrows-Wheeler transform (Преобразование Барроуза-Уилера). Он является ключевым звеном цепи, и, пожалуй, самым сложным для понимания компонентом алгоритма bzip2.
На мой взгляд, в русскоязычном интернете сейчас мало информации о BWT. По крайней мере, такое впечатление у меня сложилось два года назад, когда я работал над своим пет-проектом. Особенно сложно было разобраться с оптимизациями, а без них алгоритм вообще непригоден для практических задач. Собственно, поэтому я и решил написать эту статью. По моей задумке, она очень пригодится тем полутора землекопам, которые захотят повторить мой путь.
Сначала мы разберемся с механизмом работы преобразования Барроуза-Уилера и напишем наивную реализацию. Затем мы постараемся понять, какой смысл стоит за действиями, которые алгоритм выполняет, и почему это реально работает. И наконец, в финальной части мы добавим ряд оптимизаций, которые превратят исходный наивный алгоритм во что-то более жизнеспособное и менее прожорливое.
Что такое BWT и как он работает?
Сейчас для упрощения мы будем рассматривать только вариант алгоритма на строках. В общем случае BWT работает на байтах, но в эту версию алгоритма мы погрузимся в последнем разделе, который будет посвящен оптимизации.
Как я уже отметил выше, BWT является алгоритмом обратимой трансформации. Это значит, что помимо метода «encode» (прямая трансформация) будет еще метод «decode» (обратная трансформация), который возвращает данные в исходное состояние.
Кроме того, BWT является перестановкой. То есть это такое преобразование, которое меняет позиции символов в строке, но не меняет сами символы.
Вот небольшая демонстрация работы BWT в формате чёрного ящика:
bwt.encode("habrahabr") -> "rhhraaa$bb"
bwt.decode("rhhraaa$bb") -> "habrahabr"
bwt.encode("abracadabra") -> "ard$rcaaaabb"
bwt.decode("ard$rcaaaabb") -> "abracadabra"Обратите внимание, как в кодированной строке образовались группы из повторяющихся символов (кроме того, появился дополнительный символ "$", но о нем чуть позже). Примечательное свойство BWT заключается в том, что на не случайных данных со скрытыми закономерностями он часто группирует символы в длинные повторяющиеся серии (а такие серии как раз хорошо сжимаются через RLE).
Вот еще для примера кусок трансформированных данных из англоязычной художественной книги:
eeeeeaoeeaioeaeeeaeeV7ZgfKHOOLXRArEs/BHQ+rujieeeEnyeReeeeeoooeoeooeeaaao
oooooooooaaiiiiieeinaaaaaaaaaaeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
eeeeeeeejeeeeeeeeeeeeeepfw\nAf874n24BZKGCufni/u/bbbI?/x5t97v9IIOoNB0pr/f
0RfQQFaeedef4Gke/CnWNeeeeeeeeaeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
eeeeeeeeeeeeeeeeeeeeeeEeeeEeeeeeeEEeeeeeeeeeeeeeeeeeeNseeeeeeeeeeeeeeeee
eeeeeEeeeeeeeeeeeeeeeeeeeeEeeeeeeeeeeeeeEeeeeeeeeeeWpeeeeeeeeeeeFPf+yMPd
6P2rteOptbY9vYVfPtqQy4Yeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
eeeeeeeeeeeLE2eeiieEeeeeeeeeEeeeeeeeeeeeeeeeeeeeeeeeeiiiieeeeeeeeeEe8bee
eeeEeeeeeeeeEEeeeeeeeeeeeeeeeeeeeiiiiiiiiiee9FFtSp26Lpaeeuees/VKDB1Taaaa
eQTDIQdbPooooooOoo9KoowazlAPB+Kea3rFMCaalrblzrbllltabtadmmrrtrltcttMlaar
lhollralbeatheeedlrrlhlelleeleerllmdllelnnldmlemalBeeelhnamrleeelelleael
dbllallellhhllrlklllblbmblhlltttrlnddlbnMMnmMMmMmmmMrlhlllhlhhhahhlrlhll
ldplhugghlmlbbmbrbbllbtttbbahhhlllrrhdhaeerktkkhblrllllbblbaabblblllllll
lbablllblllblbllbbarbrbblblrallaarzraalprrrrrprrlmnrgellllcalllrendeelll
nlmllrllllrllrlltalrdlblalallMktlttlbblaeBlebnbaaatlaaeelleshlllbbllbllm
rtdlmglrradrrxtrtlranrdaltrekltrrtrzrrarttrcrgrtkgaaghblrrrtrtbttgrrstll
rllltbrarbmlrhlbflbkhbaolaatmmlkllndmteeeeeeeleleenehheeeeeeehebtmmlmmmr
kmnnallllllalalMmelrrgvlrtdlllellllslltlbarkvlglalzrlalnmmmllllllllabggt
allcrllllrllllngenaallllMalladlaallrmrllaadllaallaaammmamdllllndnnnnraee
laaalaaalllllllllddllrlctakrddllolblaelaldrmmmmrrmmmelnalaalamdlmrrrnnbn
lrmlbllltmmMnmmMmnrMtehoemmmmmmmMmtmrrMammverllllelmplheehclreeldttlrlrl
larreeureeeedleeeeeedeteeeeehhleeeeeehhhllrlrmlrntmthnnBllrknmlmrnlllnba
errellneettklneadeermmmdbelmmbrrrllllrrntdlnlglllnleemlnrlllbllelglrrelk
eeeeeedeeeeeemeeeehleeeeeeeellrMnelrrdrrbllrmtkdnaeeeeterrrlnnssrttallmrОбратите внимание на частоту и длину повторов в получившемся тексте. Сжатие через RLE покажет высокую производительность на таких входных данных.
Пора заглянуть внутрь черного ящика. Давайте посмотрим на не оптимизированную (зато наглядную) реализацию.
Encoding
На вход BWT принимает строку символов. Пусть это будет строка "habrahabr".
Далее BWT будет работать с ней как с кольцом (не путать с одноименным термином из общей алгебры): последний символ строки сшивается с первым, и алгоритму становится неважно, где исходная строка начиналась, а где заканчивалась. BWT не видит различий между "habrahabr" и "abrahabrh".
Но мы-то видим эти различия. Поэтому, чтобы не потерять информацию о месте конца строки, мы помечаем его специальным служебным символом "$". Назовём его «терминатор».
"habrahabr" => "habrahabr$"Необходимо, чтобы терминатор был уникален в пределах кодируемой строки. Эта необходимость сразу же поднимает ряд дополнительных вопросов, но мы пока оставим эти вопросы за скобками.
Теперь нам нужно построить таблицу всех циклических сдвигов. Она (в полном соответствии со своим названием) включает в себя все циклические смещения исходной строки:
"habrahabr$"
"abrahabr$h"
"brahabr$ha"
"rahabr$hab"
"ahabr$habr"
"habr$habra"
"abr$habrah"
"br$habraha"
"r$habrahab"
"$habrahabr"На следующем шаге мы сортируем строки таблицы в алфавитном порядке:
"$habrahabr"
"abr$habrah"
"abrahabr$h"
"ahabr$habr"
"br$habraha"
"brahabr$ha"
"habr$habra"
"habrahabr$"
"r$habrahab"
"rahabr$hab"Далее мы берем последний столбец таблицы (т.е. последнюю букву из каждой строки). Это и есть результат прямой трансформации:
"$habrahab" + "r"
"abr$habra" + "h"
"abrahabr$" + "h"
"ahabr$hab" + "r"
"br$habrah" + "a"
"brahabr$h" + "a"
"habr$habr" + "a"
"habrahabr" + "$"
"r$habraha" + "b"
"rahabr$ha" + "b"
╰→ "rhhraaa$bb"На что стоит обратить внимание
Символ конца строки будет ставить нам палки в колеса. Мы выбираем его из какого-то алфавита, а нам может попасться текст, содержащий все символы этого алфавита, и выбрать уникальный терминатор будет невозможно.
Очевидно, но все же: каждая строка таблицы перестановок содержит все символы исходной строки.
Несколько менее очевидно: каждый столбец таблицы перестановок содержит все символы исходной строки.
Более того, если бы мы как-то могли различать одинаковые символы (например, пометили бы их числами:
h₁ a₁ b₁ r₁ a₂ h₂ a₃ b₂ r₂), мы бы обнаружили, что это правило все еще выполняется. Это происходит по двум причинам:пока мы не отсортировали строки таблицы, она представляет собой симметричную матрицу, т.е. утверждения об уникальности символов в строках автоматически применимы и к столбцам (ведь из симметричности следует, что для каждой строки есть посимвольно совпадающий столбец);
при сортировке строк символы в столбцах меняются местами, но состав самих столбцов не меняется.
Использование таблицы перестановок даст нам сложность по памяти
, и это не очень хорошо.
Decoding
Декодирование устроено несколько сложнее. Основная идея следующая:
мы шаг за шагом восстанавливаем таблицу перестановок по ее последнему столбцу;
получив восстановленную таблицу, мы берем из нее строку, оканчивающуюся на терминатор. Эта строка (за вычетом терминатора) и будет исходной строкой.
Вот так выглядит начальное состояние таблицы перестановок:
".........r"
".........h"
".........h"
".........r"
".........a"
".........a"
".........a"
".........$"
".........b"
".........b"Здесь точками обозначены неизвестные нам символы. Мы знаем только последний столбец.
Кроме того, у нас есть такие предпосылки:
строки таблицы отсортированы в алфавитном порядке;
все столбцы таблицы составлены из одинакового набора символов (см. Encoding. На что стоит обратить внимание).
Из этого сразу следует, что мы можем восстановить первый столбец таблицы. Первый столбец – это отсортированный последний столбец:
"$........r"
"a........h"
"a........h"
"a........r"
"b........a"
"b........a"
"h........a"
"h........$"
"r........b"
"r........b"В наших дальнейших рассуждениях весомую роль сыграет тот факт, что строки таблицы — это все циклические перестановки одной и той же строки.
К примеру, мы могли бы провернуть такую цепочку рассуждений:
В таблице есть строка
"h........$".Значит, в зацикленной исходной строке после
"$"будет идти"h":"....$h...."Но тогда в строке
"$........r"мы можем восстановить ещё один символ:"$h.......r"
Но мы можем это сделать только потому, что терминатор присутствует в строке в единственном экземпляре. К сожалению, такие рассуждения дают сбой на неуникальных символах. После "a" может идти "b" (строка 6), а может "h" (строка 7).
Первый и последний столбцы таблицы дают нам информацию обо всех парах идущих подряд символов в строке (потому что строки закольцованы):
"$........r" // ....r$....
"a........h" // ....ha....
"a........h" // ....ha....
"a........r" // ....ra....
"b........a" // ....ab....
"b........a" // ....ab....
"h........a" // ....ah....
"h........$" // ....$h....
"r........b" // ....br....
"r........b" // ....br....И теперь, используя эту информацию, мы можем восстановить второй столбец таблицы:
// 1. non-sorted 2. sorted 3. with last column
"r$........" "$h........" "$h.......r"
"ha........" "ab........" "ab.......h"
"ha........" "ab........" "ab.......h"
"ra........" "ah........" "ah.......r"
"ab........" ==> "br........" ==> "br.......a"
"ab........" "br........" "br.......a"
"ah........" "ha........" "ha.......a"
"$h........" "ha........" "ha.......$"
"br........" "r$........" "r$.......b"
"br........" "ra........" "ra.......b"Пояснение:
(non-sorted) Для всех строк в таблице мы делаем сдвиг влево на один символ. К примеру, строка
"a........h"превратится в"ha........". Важно:так как изначально строки были разные, и мы одновременно применили к ним это преобразование, то мы гарантированно получим весь набор строк таблицы без повторов;
сортировка при этой операции теряется.
(sorted) Таблица перестановок должна быть отсортирована (по определению алгоритма) — сортируем строки в алфавитном порядке.
(with last column) В отсортированной таблице перестановок мы знаем, каким будет последний столбец, поэтому смело добавляем его.
Но теперь у нас на руках уже все возможные подстроки длиной 3. И мы можем повторить предыдущие шаги, пока не достроим всю таблицу:
"$habrahabr"
"abr$habrah"
"abrahabr$h"
"ahabr$habr"
"br$habraha"
"brahabr$ha"
"habr$habra"
"habrahabr$" // <- "habrahabr"
"r$habrahab"
"rahabr$hab"Теперь мы берем строку, оканчивающуюся на "$", отбрасываем последний символ и получаем исходный текст.
На что стоит обратить внимание
Теперь становится ясно, зачем нам нужен символ конца строки: без него невозможно понять, который из циклических сдвигов является «нулевым». Как следствие, трансформированный текст будет длиннее на один символ.
Для декодинга мы используем таблицу перестановок (
памяти — нехорошо).Кроме того, мы
раз сортируем строки длиной в 
; в общем случае это займет у нас 
времени. Довольно печальный результат, который очень хочется оптимизировать (спойлер: у нас получится).Столбцы таблицы сдвигов несут тем больше информации об исходном тексте, чем правее они расположены:
самый левый столбец содержит только информацию о составе и количестве символов, но никакой информации о порядке (такая информация уничтожается процессом сортировки);
столбец из середины (пусть он расположен на позиции
) уже позволяет восстановить все куски исходных данных длиной 
. Для большинства текстов мы даже могли бы восстановить исходный текст, не проходя алгоритм до конца: если мы имеем "abc","bcd","cda"и"dab", то единственный способ склеить из них слово из четырех букв будет"abcd"(точнее, все циклические сдвиги такой строки);и, наконец, последний столбец содержит информацию обо всех других столбцах таблицы. Это наводит на мысль, что таблица избыточна в плане использования памяти.
Имплементация наивного алгоритма
"""
Алгоритм без претензий на оптимальность:
- O(n^3 log n) по времени
- O(n^2) по памяти
"""
END = "$"
DBG = True
def encode(text: str) -> str:
assert END not in text, "входная строка не должна содержать символ конца строки"
# Формируем таблицу циклических перестановок:
characters = list(text)
characters.append(END)
table = [characters.copy()]
for _ in range(len(characters) - 1):
first_char = characters.pop(0)
characters.append(first_char)
table.append(characters.copy())
table.sort()
if DBG:
print()
print(*["".join(r) for r in table], sep="\n")
print()
# Извлекаем последний столбец:
encoded = [row[-1] for row in table]
return "".join(encoded)
def decode(encoded: str) -> str:
n = len(encoded)
last_column = list(encoded)
table = ["" for _ in range(n)]
# Восстанавливаем таблицу:
for _ in range(n):
# Переносим символ из последнего столбца в начало строки:
for i in range(n):
table[i] = last_column[i] + table[i]
table.sort()
# Находим строку с нулевым циклическим сдвигом:
for row in table:
if row.endswith(END):
return row[:-1]
raise ValueError
if __name__ == "__main__":
input_str = input() # "habrahabr"
encoded = encode(input_str)
if DBG:
print(encoded)
decoded = decode(encoded)
if DBG:
print(decoded)
assert input_str == decoded
Почему BWT работает
Только что мы в лучших традициях современной системы образования вбивали себе в голову многоступенчатый алгоритм, действуя по принципу «заткнись-и-считай». А сейчас мы предпримем амбициозную попытку осознать глубинный смысл сделанных нами манипуляций.
Почему в трансформированном тексте появляются повторы? В общем случае речь идет не о тексте, а скорее о потоке байт. Но для наглядности пока разберемся на примере текста.
Допустим, мы захотели применить преобразование Барроуза-Уилера на тексте настоящей статьи. Она написана на русском с латинскими вставками. Одной из таких характерных вставок будет слово "BWT" — оно встречается в тексте много раз. Когда мы будем формировать таблицу смещений, то получим нечто такое:
... // какие-то другие строки
"WT.......очень_много_текста.......B"
"WT.......очень_много_текста.......B"
"WT.......очень_много_текста.......B"
... // какие-то другие строкиСреди всех возможных циклических перестановок будут и такие, которые разрезают текст по промежутку между "B" и "WT" в слове "BWT". Ввиду того, что строки таблицы отсортированы, эти перестановки будут сгруппированы. А значит, в последнем столбце мы получим длинную серию из латинских "...BBBB...".
Если какому-то неизвестному символу в тексте этой статьи предшествует строка WT, то с большой долей достоверности мы можем утверждать, что ...?WT... - это ...BWT... (в другой конфигурации эти латинские буквы в тексте статьи почти не встречаются).
Кстати, из-за того, что абзацем выше я поставил знак вопроса вместо первой буквы BWT, в таблице сдвигов появится вот такой фрагмент:
... // какие-то другие строки
"WT.......очень_много_текста.......B"
"WT.......очень_много_текста.......?"
"WT.......очень_много_текста.......B"
... // какие-то другие строкиПозиция знака ? в серии ...BBBBB... определяется тем, как отсортируется текст после него. Этот знак вопроса может появиться где-то посередине последовательности из «бэшек», в будто бы случайном (а на самом деле, строго детерминированном) месте. В одной из версий текста статьи кусок последовательности с вопросительным знаком выглядел вот так: BBBBBBBB"BBB BBBBBBB?""""""/--eeLAAAAAA.
Другими словами, одни и те же символы (буквы или байты) в BWT-кодированных данных собираются в группы потому, что лексикографическая сортировка группирует вместе их суффиксы.
Альтернативное объяснение через цепи Маркова
Еще один способ интерпретации BWT — цепь Маркова для генерации случайного текста. Цепи Маркова можно рассматривать как расширение конечных автоматов. В них переход в новое состояние происходит с определенной вероятностью, в отличие от конечных автоматов, где переходы детерминированы.
На основе цепей Маркова можно делать языковые модели для генерации текста (правда, по всем параметрам они будут уступать языковым моделям на нейронных сетях). Состояние модели определяется префиксом в тексте. Например, таким состоянием может быть четырехсимвольный префикс "и.т.". Переходами из этого состояния могут быть, например, переход в ".т.д" с вероятностью 83 % и переход в ".т.п" с вероятностью 17 %. Вероятности переходов определяются на основе анализа большого количества текстов. Перемещаясь из состояния в состояние, можно генерировать текст произвольной длины, в котором каждая новая буква выбирается случайно на основе 
предыдущих букв.
Количество состояний модели растет экспоненциально от длины суффикса, и если у нас есть русский алфавит без заглавных букв, но с пробелами и минимальными знаками препинания, то суффикс длиной в 4 символа может дать нам 
— полтора миллиона состояний (это оценка сверху, т.к. не все возможные сочетания символов будут встречаться в обучающих текстах).
Архивация через BWT похожа на очень странную попытку построить такую Марковскую цепь.
Каждая строка таблицы сдвигов — это как бы предсказание символа по его суффиксу. Если бы мы убрали в таблице все столбцы, оставив только 
первых и один последний, и сгруппировали бы эти данные, то получили бы цепь Маркова, которая угадывает следующую букву текста по 
предыдущим.
"$h.......r" "$h": [("r", 1.0)] # перед "$h" идёт "r" с вероятностью 1.0
"ab.......h" "ab": [("h", 1.0)]
"ab.......h"
"ah.......r" "ah": [("r", 1.0)]
"br.......a" ==> "br": [("a", 1.0)]
"br.......a"
"ha.......a" "ha": [("a", 0.5), ("$", 0.5)]
"ha.......$"
"r$.......b" "r$": [("b", 1.0)]
"ra.......b" "ra": [("b", 1.0)]Если же мы не удаляем ни одного столбца (т.е. берем суффикс размером с весь текст), то цепь Маркова вырождается в конечный автомат, который детерминированно выдает следующий символ.
Примечательно, что процесс восстановления текста происходит задом наперед: мы вычисляем букву по суффиксу, который за ней следует. В естественных языках мы формируем фразы, слова и буквы в порядке, противоположном работе BWT.
Можно предположить, что человеческий язык эволюционировал таким образом, чтобы слушатели и читатели могли предугадывать слова, когда читают/слышат их. И тогда резонно было бы ожидать, что если бы BWT строил таблицу по тексту в обратном порядке, то это дало бы нам прирост в эффективности (мы получили бы больше длинных серий из повторяющихся символов).
Возможные контраргументы к этому доводу:
Утверждение опирается на неверную предпосылку, что процесс обработки текста алгоритмом аналогичен человеческому чтению. Но процесс предугадывания человеком слов при чтении устроен не так прямолинейно: мы охватываем все слово одним взглядом, и наш мозг угадывает его по кортежу
(первая_буква, последняя_буква, примерная_длина, состав_букв_посередине).Достаточно долго существующая письменность постепенно перегружается атавизмами, и это привносит недетерминированность. Написание слов становится малосопоставимо с их произношением (посмотрите, к примеру, на английский язык с его ough-словами).
Я решил проверить это предположение на практике и прогонял BWT в прямом и обратном порядке на разных русских и английских текстах. В качестве метрики эффективности использовал суммарную длину всех серий из повторов длиннее двух. Но мои опыты не выявили какого-то явного преимущества в порядке обработки.
Оптимизация
Как мы выяснили выше, ситуация с потреблением ресурсов у текущей версии алгоритма очень печальная. Если попытаться применить код наивной имплементации на текст умеренного размера (например, текст данной статьи без картинок), то в ожидании результата можно не только успеть попить чай, но и даже поразмышлять над некоторыми фундаментальными вопросами теории вычислимости. А владельцы менее оснащенных машин даже рискуют словить OutOfMemoryError.
Злопыхатели скажут, что всему виной выбор языка. Дескать, чего еще ожидать от питона. Этот язык непригоден для написания кода архиваторов.
Но на самом деле гораздо большую роль играет сама реализация. Ее достоинство — наглядность. Ее недостаток — она, мягко скажем, не вполне оптимизирована.
Во-первых, таблица циклических сдвигов занимает 
памяти, и это очень плохо. Квадратичная память гораздо хуже квадратичного времени. Сейчас у нас вошло в привычку пренебрегать затратами RAM в угоду скорости разработки, но конкретно в этом случае такой подход нам выйдет боком. Квадрат по памяти — это угроза любому размеру RAM. Если исходный блок данных имеет размер 128 KiB (
), то его таблица перестановок будет уже иметь размер как минимум 16 GiB (
) (в данном случае предполагается, что 1 символ = 1 байт).
Во-вторых, сортировка 
строк длиной в 
— это 
. А при декодировании эту сортировку нужно делать 
раз, и тогда уже совсем тушите свет: 
. То есть, сравнивая строки друг с другом, в некоторых случаях нам будет достаточно проверить два первых символа, а в некоторых — зайти гораздо глубже. И это дает нам сложность на одну операцию сравнения 
, где параметр 
— это коэффициент глубины захода. Но O-нотации нет дела до константных множителей.
Сейчас мы попробуем оптимизировать наш алгоритм, и попутно перепишем его со строк на байты.
Оптимизация потребления памяти при кодировании
Как при кодировании, так и при декодировании мы работаем с таблицей сдвигов, и памяти она занимает изрядно. При этом все строки в ней — это де-факто одна и та же строка, но с разными значениями размера сдвига:
"habrahabr$" 0
"abrahabr$h" 1
"brahabr$ha" 2
"rahabr$hab" 3
"ahabr$habr" = "habrahabr$" + 4
"habr$habra" 5
"abr$habrah" 6
"br$habraha" 7
"r$habrahab" 8
"$habrahabr" 9Нам нужно создать свою структуру данных для таких строк, которая поддерживала бы операции сравнения (ради сортировки) и обращения по индексу (чтобы достать последнюю колонку). Вот пример реализации такой структуры данных:
from typing import Self
def encode(text: str) -> str:
"""Версия BWT-энкодера, которая потребляет O(N) памяти"""
array = list(text)
array.append("$")
array_size = len(array)
class ComparableRotation: # динамически объявляемый класс
def __init__(self, shift: int) -> None:
self.shift = shift
# перегрузка оператора `[]`
def __getitem__(self, index: int) -> str:
return array[(index + self.shift) % array_size]
# перегрузка оператора `<`
def __lt__(self, other: Self) -> bool:
for i in range(array_size):
if self[i] < other[i]:
return True
elif self[i] > other[i]:
return False
elif self[i] == other[i]:
continue
return False
table = [ComparableRotation(shift) for shift in range(array_size)]
table.sort()
# Извлекаем последний столбец:
encoded = [row[-1] for row in table]
return "".join(encoded)
Избавляемся от терминатора
Задача символа-терминатора — быть индикатором конца строки. Ведь декодирование BWT дает нам не исходный текст, а его закольцованную версию. И мы понятия не имеем, в каком месте это кольцо нужно разорвать. Точнее, в восстановленной таблице мы не сможем определить позицию строки с нулевым смещением (которая и будет являться исходным текстом). С символом-терминатором мы просто брали строку, которая на него заканчивается.
Но используя нашу структуру данных ComparableRotation мы и так можем определить позицию символа, который был последним в тексте: это будет позиция объекта ComparableRotation со смещением 0.
Таким образом, в качестве кодированных данных мы будем возвращать два значения: кодированную строку и индекс последнего элемента в ней.
Финальный код будет выглядеть вот так:
from typing import Self
def encode(text: str) -> tuple[str, int]:
"""Версия BWT-энкодера, которая потребляет O(N) памяти"""
array = list(text)
array.append("$")
array_size = len(array)
class ComparableRotation: # динамически объявляемый класс
def __init__(self, shift: int) -> None:
self.shift = shift
# перегрузка оператора `[]`
def __getitem__(self, index: int) -> str:
return array[(index + self.shift) % array_size]
# перегрузка оператора `<`
def __lt__(self, other: Self) -> bool:
for i in range(array_size):
if self[i] < other[i]:
return True
elif self[i] > other[i]:
return False
elif self[i] == other[i]:
continue
return False
table = [ComparableRotation(shift) for shift in range(array_size)]
table.sort()
# Извлекаем последний столбец:
encoded = []
origin = -1
for i, row in enumerate(table):
if row.shift == 0:
origin = i # индекс строки с нулевым сдвигом
encoded.append(row[-1])
return ("".join(encoded), origin)
Этот код уже с минимальными изменениями можно использовать для работы с байтами (полная имплементация будет в конце статьи).
Оптимизация декодирования
Вернемся к процессу декодирования таблицы перестановок из неоптимизированного алгоритма и посмотрим на процесс повнимательнее. Будем использовать версию с терминатором, потому что с ним будет удобнее ориентироваться.
Здесь стрелками отмечено, в какое место попал символ из последнего столбца после сортировки на первом шаге. В слове "habrahabr" есть несколько повторяющихся символов, поэтому для детерминированности будем использовать стабильный алгоритм сортировки (т.е. такой, который сохраняет взаимный порядок одинаковых символов).
Получился такой своеобразный двудольный граф, который описывает перестановку символов. А теперь самое интересное: перестановка остается неизменной вне зависимости от шага трансформации. Вот, к примеру, двудольный граф для второго шага:
Вообще здесь совсем не очевидно, почему граф не меняется…
…ведь если в данных встречаются повторяющиеся символы, то откуда мы знаем, в каком порядке их расположить друг относительно друга? На первом шаге порядок был аналогичен стабильной сортировке, и на втором тоже. Но где гарантия, что на шаге 
порядок одинаковых букв не поменяется?
Оказывается, эта гарантия следует из того, что до добавления последнего столбца строки были отсортированы. Вот небольшой демонстрационный пример:
1. 2. 3.
... ... ...
"aaa...d" "daaa..." "daaa..."
"aab...e" => "eaab..." => "daac..."
"aac...d" "daac..." "eaab..."
... ... ... Пояснение:
На очередном шаге декодирования (восстановления таблицы сдвигов) мы имеем две буквы
dв последнем столбце (таблица в данный момент отсортирована).Переносим последний столбец влево (теряем сортировку).
Сортируем строки лексикографически. Буквы
dсохранили свой порядок друг относительно друга, т.к."aaa"лексикографически меньше, чем"aac".
Используя выше описанный двудольный граф, мы можем очень легко восстановить исходную строку. Вот как это делается:
Пояснение:
Находим терминатор в последнем столбце (в варианте без терминатора находим символ с нулевым смещением по предоставленному индексу). Соответствующая ему буква в первом столбце будет идти сразу за ним (строка
"h........$"третья снизу).Согласно графу переходов, буква
hв начале третей снизу строки — это букваhв конце третьей сверху строки (которая"a........h").За ней будет следовать буква
a. Получаем"ha...".Согласно графу переходов,
aиз начала третьей строки — этоaиз конца шестой строки ("b........a"). За ней будет следовать букваb("hab...").Повторяем предыдущий шаг еще 6 раз.
Для построения двудольного графа переходов нам нужны только первый и последний столбцы таблицы. Последний столбец — это наши входные данные. Первый столбец мы получаем сортировкой из последнего (
). Также для хранения вспомогательных структур данных используется 
дополнительной памяти (что не критично, т.к. входные данные уже занимают 
).
Имплементация:
def decode(encoded: str, origin: int) -> str:
n = len(encoded)
permutation = list(range(n))
permutation.sort(key=lambda i: encoded[i]) # получаем перестановку
current = origin
decoded = []
for _ in range(n):
current = permutation[current] # получаем индекс следующего символа
decoded.append(encoded[current])
return "".join(decoded)
Оптимизация времени при кодировании
Следующее «бутылочное горлышко», в которое мы упремся — это время, занимаемое операцией кодирования. Если на вход подается блок данных объемом 100 килобайт (вполне адекватный размер блока данных для архивации), то его кодирование может измеряться в секундах (на моём скромном AMD Ryzen 5 PRO 3500U — около полутора секунд на одном ядре). Тому есть несколько причин.
Во-первых, при кодировании мы сортируем действительно большой объем данных. В текущей реализации мы работаем над байтовым алфавитом, а это значит, что данные в 100 килобайт будут порождать 100 000 циклических перестановок. И список из них нужно отсортировать. Максимальная эффективность алгоритмов сортировки: 
.
Во-вторых, мы сортируем массивы, а не числа. Операция сравнения двух массивов (или строк) занимает 
времени. Работа с массивами также не позволяет нам эффективно применить алгоритмы сортировки за линейное время типа Radix Sort. Совокупная сложность операции сортировки таблицы циклических сдвигов таким образом становится 
. Тут, правда, стоит заметить, что практически никогда для сравнения двух строк нам не придется проходить по всем их символам до конца (если только мы не подсунем какие-то специально подготовленные данные). В моих экспериментах на блоках из 100 килобайт английского текста нужно было сделать максимум 23 сравнения (23 первых символа обеих сравниваемых строк совпадали), и это были единичные случаи. То есть, это будет 
. Но в терминологии О-нотации это все равно считается за 
.
В-третьих, мы используем Python — интерпретируемый язык программирования. И по скорости он сильно уступает практически любому компилируемому языку. Различные источники из интернета оценивают Python как язык в 10-100 раз более медленный, чем C.
Что же делать? Тратить десять секунд на архивацию одного мегабайта данных кажется не очень практичным, как, впрочем, и сама идея писать архиватор на Python.
Можно поменять интерпретатор (использование pypy позволяет улучшить производительность ценой бóльших расходов памяти; в моих тестах pypy работает в 2.5 раз быстрее). Вычисление отдельных блоков можно легко распараллелить на несколько ядер и получить кратное ускорение. Но это уже будет история не совсем про алгоритмы.
Существуют хитрые алгоритмы сортировки циклических сдвигов с использованием суффиксных деревьев, которые выполняются со сложностью 
по времени и 
по памяти. Такое возможно потому, что мы сортируем не случайные строки, а сдвиги одной и той же строки. Но мне кажется, тут это будет оверкилл. На небольших блоках данных время работы растет примерно как 
(да, я только что отмахнулся от своего же утверждения несколькими параграфами выше, что мы не можем игнорировать константу 
в 
).
Если проанализировать выполнение алгоритма кодирования при помощи входящего в стандартную библиотеку инструмента cProfile, то видно, что большая часть времени приходится на операции сравнения (в этом примере: 1.471 из 2.316 секунд):
Все сравнения выполняются через метод __lt__ класса ComparableRotation, который мы придумали, чтобы оптимизировать хранение памяти. Исходно этот метод выглядит вот так:
def __lt__(self, other: Self) -> bool:
for i in range(array_size):
if self[i] < other[i]:
return True
elif self[i] > other[i]:
return False
elif self[i] == other[i]:
continue
return FalseНо если переписать вот так:
def __lt__(self, other: Self) -> bool:
si = self.shift
oi = other.shift
for i in range(array_size):
s = array[si]
o = array[oi]
if s < o:
return True
elif s > o:
return False
si += 1
oi += 1
if si == array_size:
si = 0
if oi == array_size:
oi = 0
return False
то сравнение будет работать примерно в 3 раза быстрее. Полагаю, это происходит потому, что мы используем обращение по индексу объекта array (который является экземпляром стандартного класса list), а не объектов self и other, чьи методы обращения к индексу мы заменили на самописные.
Еще одна оптимизация: Bucket sort (блочная сортировка / корзинная сортировка) — один из видов сортировки за линейное время, о которых я писал несколькими параграфами выше, что они не будут работать.
Чистая блочная сортировка действительно показывает результаты худшие, чем стандартный питоновский sort. Но если их скомбинировать, то можно получить двукратное ускорение. Идея в том, чтобы сделать несколько первых итераций блочной сортировки, а затем отсортировать содержимое получившихся «корзин» уже стандартным sort. Такой подход позволяет добиться двукратного (x2) увеличения скорости сортировки.
MAX_DEPTH = 6 # MAX_DEPTH подбирается экспериментально
def _sort( rows: list[Any], depth: int = 0) -> list[Any]:
if depth == MAX_DEPTH: # условие выхода из рекурсии
rows.sort() # стандартный алгоритм сортировки
return rows
buckets: dict[Any, list[Any]] = {}
for row in rows:
key = row[depth]
if key not in buckets:
buckets[key] = []
buckets[key].append(row)
ans: list[Any] = []
for key in sorted(buckets.keys()):
# рекурсивно сортируем каждую из корзин
buckets[key] = _sort(buckets[key], depth + 1)
ans.extend(buckets[key])
return ansТаким образом, получаем в совокупности ускорение в 3*2=6 раз. Т.е. теперь на кодирование блока в 100 килобайт будем тратить не 1,5 секунды, а 0,25 секунды.
Плюс, небольшое ускорение примерно в 5–15 % даст переход со строк на байты.
Конечный исходный код со всеми оптимизациями
import time
from typing import Any, Self
MAX_DEPTH = 6 # MAX_DEPTH подбирается экспериментально
class Bwt:
def _sort(self, rows: list[Any], depth: int = 0) -> list[Any]:
# Оптимизация быстродействия: частичная блочная сортировка
if depth == MAX_DEPTH: # условие выхода из рекурсии
rows.sort() # стандартный алгоритм сортировки
return rows
buckets: dict[Any, Any] = {}
for row in rows:
key = row[depth]
if key not in buckets:
buckets[key] = []
buckets[key].append(row)
ans = []
for key in sorted(buckets.keys()):
# рекурсивно сортируем каждую из корзин
buckets[key] = self._sort(buckets[key], depth + 1)
ans.extend(buckets[key])
return ans
def encode(self, data: bytes) -> tuple[bytes, int]:
array = list(data)
array_size = len(array)
class ComparableRotation:
__slots__ = ("shift",)
# Оптимизация памяти: структура данных для сравнения циклических перестановок
def __init__(self, shift: int) -> None:
self.shift = shift
def __getitem__(self, index: int) -> int:
return array[(index + self.shift) % array_size]
def __lt__(self, other: Self) -> bool:
# Оптимизация быстродействия: более эффективное сравнение
si = self.shift
oi = other.shift
for _ in range(array_size):
s = array[si]
o = array[oi]
if s < o:
return True
elif s > o:
return False
si += 1
oi += 1
if si == array_size:
si = 0
if oi == array_size:
oi = 0
return False
table = [ComparableRotation(shift) for shift in range(array_size)]
table = self._sort(table)
encoded = []
origin = -1
for i, rot in enumerate(table):
if rot.shift == 0:
origin = i
encoded.append(rot[-1])
return (bytes(encoded), origin)
def decode(self, data: tuple[bytes, int]) -> bytes:
encoded, origin = data
n = len(encoded)
# формируем двудольный граф переходов
transmissions_graph = list(range(n))
transmissions_graph.sort(key=lambda i: encoded[i])
current = origin
decoded = bytearray(n)
for i in range(n):
current = transmissions_graph[current]
decoded[i] = encoded[current]
return bytes(decoded)
if __name__ == "__main__":
bwt = Bwt()
input_bytes = None
block_size = -1 # 100 * 2**10 # 100 KiB
with open("./Watts Peter. Blindsight.fb2", "rb") as f:
input_bytes = f.read(block_size)
start = time.process_time()
encoded = bwt.encode(input_bytes)
decoded = bwt.decode(encoded)
print(f"The input size is {len(input_bytes)} bytes.")
print(f"Encoding and decoding took {time.process_time() - start:.4f} seconds.")
assert input_bytes == decoded
В сухом остатке
В отличие от других компонентов алгоритма bzip2 (RLE -> BWT -> MTF -> RLE -> HFC), преобразование Барроуза-Уилера — вещь очень нишевая и малоизвестная. На мой взгляд, это еще и самый сложный компонент всего алгоритма. К тому же, в моей реализации оно является самым узким местом в плане потребляемых ресурсов, несмотря на все примененные оптимизации. Что же касается других звеньев цепи…
RLE — обманчиво прост, но в целом не представляет серьезного челленджа (и у меня имеется статья по нему).
MTF — мне кажется, его эффективную имплементацию можно написать с нуля минут за 40 (и в это время я закладываю чтение Вики-статьи о нем).
HFC — я бы по сложности расположил где-то между RLE и BWT. Кодирование Хаффмана имеет множество практических приложений, и по нему есть огромное количество информации в интернете с наглядными визуализациями.
Итого: чтобы полностью закрыть тему bzip2, нам осталось разобрать Move-To-Front, кодирование Хаффмана, и, возможно, некоторые общие особенности компоновки данных. Я планирую написать статью и об этом, но практика показывает, что такие планы очень медленно воплощаются в жизнь. Если вы смогли осилить этот материал и хотя бы примерно понять, что тут происходит, то:
во-первых, спасибо вам, что разделили со мной этот опыт;
во-вторых, вы, вероятно, справитесь со всей оставшейся цепочкой и без моей помощи.
А тем читателям, которым не хватило технической духоты, я рекомендую другие статьи из нашего блога. Там мои коллеги разбирают чуть менее наукоемкие, но чуть более применимые на практике нюансы работы с Python:
