Динамическое сжатие Маркова.

6.11.97, 30.03.1993

Данное документ составлен на основе первоначального описания алгоритма, опубликованного в [1]. Кратко описывается алгоритм и приводятся результаты сравнения этого алгоритма с другими методами сжатия. Более подробно о работе алгоритма, в частности, о тонких местах в выборе параметров и совместном использовании GUAZZO CODING ALGORITHM, можно прочитать в [1]. В работе [2] показывается, что модель Динамического сжатия Маркова эквивалентна автомату с конечным контекстом (Finite Context Automata, FCA). Т.к. класс FCA является подмножеством класса автоматов с конечным числом состояний (Finite State Automata, FSA), то делается предположение, что т.к. фрагменты текстов на английском языке могут быть распознаны FSA, но не могут распознаваться Марковской моделью переменного порядка (variable-order Markov model), то использование модели FSA потенциально должно давать лучшие результаты сжатия.

1. G.V. Cormack, R.N. Horspool.
   Data Compression Using Dynamic Markov Modelling.
   // The Computer Journal, vol.30,No.6,1987, pp.541-550
2. T. Bell, A. Moffat.
   A note on the DMC data compression scheme.
   // The Computer Journal, vol.32,No.1,1989, pp.16-20
------------------------

На практике, сообщения почти всегда имеют некоторую степень зависимости между соседними символами, что соответствует источнику сообщений с памятью. В частности, таким источником может быть Марковская модель (в каждом состоянии она помнит только об одном предыдущем своем состоянии).

Для примера используем простую Марковскую модель. Будем считать, что источник порождает строки, состоящие из двоичных цифр такие, что:

     за 0 следует 0 с вероятностью 2/3                                        
     за 0 следует 1 с вероятностью 1/3                                        
     за 1 следует 0 с вероятностью 2/5                                        
     за 1 следует 1 с вероятностью 3/5                                        
и первая цифра может быть 0 или 1 равновероятно (1/2).

Сообщения, порождаемые этой моделью имеют общее свойство: за 0 следует чаще другой 0, чем 1. Аналогично - за 1 чаще следует 1, чем 0. Марковская модель может быть описана с помощью ориентированного графа, каждой дуге которого приписываются соответствующие вероятности перехода. Диаграмма состояний для этой модели выглядит следующим образом:

                       ---                                                    
                      | А |                                                   
                     / ---  \
                    /        \                  
            0(50%) /          \ 1(50%)                                                     
                  /            \
                |/   / 0 (40%)  \|            
            \  ---  ---------  ---  /                                         
          --- | B |           | C | ----                                      
         |     ---  ---------  ---      |                                     
  0 (67%)|      |   1 (33%) /   |       |1 (60%)                              
          ------                 --------                                     

При автоматическом создании такой модели возникает две проблемы: проблема определения вероятностей перехода, приписываемых дугам графа, и проблема в определении самой структуры графа. Эти проблемы рассматриваются по отдельности.

Определение вероятностей перехода

Пусть уже известна структура графа. По мере прочитывания двоичных цифр сообщения Марковская модель будет переходить из одного своего состояния в другое в соответствии с прочитанными символами. Необходимо подсчитать сколько раз модель будет переходить по каждой дуге (для каждой дуги заводятся счетчики переходов по этой дуге). С помощью этих счетчиков можно получить оценки вероятностей переходов. Если переход из состояния А по дуге для символа 0 проходил N0 раз, а по дуге для 1 - N1, то оценки вероятностей будут выглядеть следующим образом:

     Prob{digit=0 | current state=A} = N0/(N0 + N1)
     Prob{digit=1 | current state=A} = N1/(N0 + N1)
Чем чаще мы будем попадать в состояние А, тем точнее будут эти оценки.

При адаптивном кодировании заранее определить значения N0 и N1 невозможно. В этом случае эти оценки могут выглядеть следующим образом:

     Prob{digit=0 | current state=A} = (N0 + c)/(N0 + N1 + 2c)
     Prob{digit=1 | current state=A} = (N1 + c)/(N0 + N1 + 2c)
где с - положительная константа. Для файлов большого размера (на много больше значения c) выбор конкретного значения c существенной роли не играет. Для файлов малого размера при малом значении c модель быстрее 'выучит' характеристики источника, однако, при этом увеличивается вероятность неправильного прогноза. (Т.е. типичная проблема адаптивных алгоритмов).

Построение диаграммы состояний

Алгоритм динамического построения диаграммы состояний поясняется на простом примере: пусть уже построена модель, содержащая состояния A,B,...,E:

    ---    0   \  ---    0   \  ---                                           
   | A |-------- | C |-------- | D |                                          
    ---         / --- \         ---                                           
               /|      \
              /         \
           1 /           \ 1                                                                           
            /             \ 
           /               \ 
    ---   /                 \|  ---                                                               
   | B | /------------------ \ | E |                                          
    ---                         ---                                           

Из диаграммы видно, что в состояние C ведут две дуги из A и B и выходит две дуги: в D и E. Всякий раз, когда модель переходит в состояние C, теряется некоторая информация о контексте, т.е. забывается откуда мы пришли - из A или B. Однако вполне может быть, что выбор следующего состояния (D или E) зависит от предыдущего состояния (A или B). Чтобы учесть эту зависимость прибегнем к следующему приему: разделим состояние C на два состояния. Измененная модель будет иметь вид:

    ---    0   \  ---   0    \  ---                                           
   | A |-------- | C |-------- | D |                                          
    ---           --- \ 1   /|   ---                                          
                       \   /                                                  
                        \ /                                                   
                         /                                                    
                      0 / \
                       /   \
    ---    1   \  --- /     \|  ---                                                             
   | B |-------- | С'|-------  | E |                                          
    ---           ---   1   /   ---                                           

Значения счетчиков переходов из C в D или E делится между счетчиками C и C' в новой модели пропорционально числу переходов из A в C и из B в C для старой модели. Теперь модель отражает степень зависимости между состояниями A,B и D,E. Для этой модели возможно выбор следующего состояния (D или E) не зависит от предыдущего состояния (A или B), но зависит от непосредственного предшественника A или B. В этом случае делим состояния A,B,C' и еще раз делим состояние C. Т.е. модель теперь будет учитывать эти зависимости.

В основном, чем больше мы делим состояний, тем больше увеличивает размах зависимостей, используемых моделью для прогноза. Само собой процесс деления необходимо производить только в том случае, когда мы уверены в существовании соответствующих зависимостей, иначе мы потеряем информацию о контексте. Отсюда желательно, чтобы деление состояния C проходило тогда, когда оба счетчика переходов AC и BC достаточно велики. На основе этого можно вывести критерий деления состояния: Выбранное состояние делится в том и только в том случае, если число переходов из текущего состояния в выбранное состояние больше, чем MIN_CNT1, и число переходов из всех других состояний, кроме текущего, в выбранное состояние больше, чем MIN_CNT2.

Начальная модель

Для полного определения метода динамического сжатия Маркова осталось определить стартовую модель. Самая простая модель выглядит следующим образом:

                       -----------------
                      |/                |
                     ---                |
          ----------|   |---------------
         |           ---                                                      
         |            | \
          -------------                              
Самая простая модель для байт-ориентированного варианта имеет 255 состояний, представляется в виде двоичного дерева у которого переход из каждого листа ведет в корень дерева.

В целом, древовидная начальная модель работает хорошо, но можно построить модель, которая будет работать немного лучше. В общем случае, состояния на k-м уровне дерева учитывает зависимость от предыдущих k-1 бит. Размер контекста, от которого зависит текущее состояние увеличивается от 0 до 7 бит, затем, при переходе к корню дерева, эта зависимость теряется. Лучшей по сравнению с этой будет модель, сохраняющая постоянным размер контекста. Т.е. такая модель должна учитывать зависимость между последними битами предыдущего байта и первыми битами текущего бита. К сожалению, такую модель трудно представить в виде диаграммы, но алгоритм построения этой модели приводится в конце данной статьи.

Последний важный момент

Этим последним моментом является вопрос когда прекращать производить деление состояния модели. Если этого не производить, то через некоторое время модель может занять всю имеющуюся память и потребовать еще. Одно из возможных решений - определить предел числа состояний модели по достижении которого модель разрушается и начинает создаваться заново из начальной модели. Последствия такой операции можно облегчить, если запомнить перед разрушением модели предыдущие n байт текста и затем по этому фрагменту улучшить начальную модель перед продолжением обработки текста.

   Зависимость коэффициента сжатия от параметров деления                      
                   (MIN_CNT1, MIN_CNT2)                                       
-----------+---------------------+-----------------------                     
Значение   |  Число вершин графа |  Коэффициент сжатия                        
параметров |                     |          (%)                               
-----------+---------------------+-----------------------                     
 ( 1,1 )*  |      > 194000       |          34.7                              
 ( 2,2 )   |        150901       |          33.8                              
 ( 4,4 )   |         84090       |          35.8                              
 ( 8,8 )   |         44296       |          38.9                              
 (16,16)   |         23889       |          42.7                              
 (32,32)   |         12089       |          46.5                              
 (64,64)   |          6347       |          50.6                              
(128,128)  |          3211       |          54.6                              
(256,256)  |          1711       |          58.6                              
-----------+---------------------+-----------------------                     
Файл, на котором проводились исследования, содержал 97393 символов ASCII текста.
* Сжав 90% исходного файла, программа исчерпала ресурс памяти для узлов графа и была прервана.
              Сравнительные результаты работы                                 
------------+-------------------------------------------------+               
            |       Исследуемые файлы                         |               
            +--------------+----------------+---------+-------+               
  Метод     |Форматированый|Неформатированый|Объектный|С-текст|               
 сжатия     |    текст     |     текст      |   код   |       |               
------------+--------------+----------------+---------+-------+               
  Adaptive  |     59.7     |      61.6      |   79.6  |  62.9 |               
  Huffman   |              |                |         |       |               
------------+--------------+----------------+---------+-------+               
   Normal   |              |                |         |       |               
 Ziv-Lempel |     38.2     |      42.9      |   98.4  |  40.8 |               
   (LZW)    |              |                |         |       |               
------------+--------------+----------------+---------+-------+               
Bit-oriented|     74.2     |      83.6      |   91.3  |  86.7 |               
   (LZ-2)   |              |                |         |       |               
------------+--------------+----------------+---------+-------+               
 Cleary and |              |                |         |       |               
   Witten   |     26.5     |      30.2      |   69.4  |  26.3 |               
    (CW)    |              |                |         |       |               
------------+--------------+----------------+---------+-------+               
  Dynamic   |              |                |         |       |               
  Markov    |     27.2     |      31.8      |   54.8  |  27.5 |               
   (DMC)    |              |                |         |       |               
------------+--------------+----------------+---------+-------+               
  File size |     74510    |      61536     |   68608 | 31479 |               
------------+--------------+----------------+---------+-------+               

Алгоритм деления состояния модели

{ NEXT_STATE[S,D] = state reached from S after trsnsition on digit D          
  TRANS_CNT[S,D]  = number of observations of input D when in state S         
  STATE           = number of current state                                   
  LAST_STATE      = largest state number used so far                          
  MIN_CNT1        = minimum number of transitions from the current            
                    state to state S before S is eligible for cloning         
  MIN_CNT2        = minimum number of visits to a state S from all            
                    predecessors of S other than the current state            
                    before S is eligible for cloning                }         
while not eof do                                                              
   begin                                                                      
      resd( B );    { read one input bit }                                    
      TRANS_CNT[STATE,B] := TRANS_CNT[STATE,B] + 1;                           
      NXT                := NEXT_STATE[STATE,B];                              
      NXT_CNT            := TRANS_CNT[NXT,0] + TRANS_CNT[NXT,1];              
      if (TRANS_CNT[STATE,B] > MIN_CNT1) and                                  
         ((TRANS_CNT - TRANS_CNT[STATE,B]) > MIN_CNT2) then                   
          begin                                                               
             LAST_STATE := LAST_STATE + 1;                                    
             NEW        := LAST_STATE;    { Obtain new state number }         
             NEXT_STATE[STATE,B] := NEW;                                      
             RATIO := TRANS_CNT[STATE,B] / NXT_CNT;                           
             for B:=0 to 1 do                                                 
              begin                                                           
               NEXT_STATE[NEW,B] := NEXT_STATE[NXT,B];                        
               TRANS_CNT[NEW,B]  := RATIO * TRANS_CNT[NXT,B];                 
               TRABS_CNT[NXT,B]  := TRANS_CNT[NXT,B] - TRANS_CNT[NEW,B]       
              end;                                                            
             NXT := NEW                                                       
          end;                                                                
       STATE := NXT                                                           
   end;                                                                       

Алгоритм генерации начального состояния

const                                                                         
   NBITS = 8;        { Number of bits per byte }                              
   STRANDS = 256;    { 2 ** NBITS }                                           
                                                                              
for I := 0 to NBITS-1 do                                                      
  for J :=0 to STRANDS-1 do                                                   
      begin                                                                   
         STATE := I + NBITS * J;                                              
         K := (I + 1) mod NBITS;                                              
         NEXT_STATE[STATE,0] := K + ((2 * j) mod STRANDS) * NBITS;            
         NEXT_STATE[STATE,1] := K + ((2 * j + 1) mod STRANDS) * NBITS;        
         TRANS_CNT[STATE,0]  := 1;                                            
         TRANS_CNT[STATE,0]  := 1                                             
      end;                                                                    
LAST_STATE := NBITS * STRANDS -1;                                             

Учебный вариант алгоритма Динамического сжатия Маркова

Ниже приводится реализация алгоритма Динамического сжатия Маркова, написанный автором алгоритма для учебных целей. Этот алгоритм уже публиковался в этой конференции и приводится здесь для полноты.


/*   Dynamic Markov Compression (DMC)    Version 0.0.0                        
                                                                              
     Copyright 1993, 1987                                                     
    Gordon V. Cormack                                                         
    University of Waterloo                                                    
    cormack@uwaterloo.ca                                                      
                                                                              
                                                                              
    All rights reserved.                                                      
                                                                              
    This code and the algorithms herein are the property of Gordon V. Cormack.
                                                                              
    Neither the code nor any algorithm herein may be included in any software,
    device, or process which is sold, exchanged for profit, or for which a    
    licence or royalty fee is charged.                                        
                                                                              
    Permission is granted to use this code for educational, research, or      
    commercial purposes, provided this notice is included, and provided this  
    code is not used as described in the above paragraph.                     
                                                                              
*/                                                                            
/*    This program implements DMC as described in                             
      "Data Compression using Dynamic Markov Modelling",                      
      by Gordon Cormack and Nigel Horspool                                    
      in Computer Journal 30:6 (December 1987)                                
                                                                              
      It uses floating point so it isn't fast.  Converting to fixed point     
      isn't too difficult.                                                    
                                                                              
      comp() and exp() implement Guazzo's version of arithmetic coding.       
                                                                              
      pinit(), predict(), and pupdate() are the DMC predictor.                
                                                                              
      pflush() reinitializes the DMC table and reclaims space                 
                                                                              
      preset() starts the DMC predictor at its start state, but doesn't       
               reinitialize the tables.  This is used for packetized          
               communications, but not here.                                  
*/                                                                            
#include <stdio.h>
                                                                              
float predict();                                                              
int pinit();                                                                  
int pupdate();                                                                
                                                                              
int memsize = 0x1000000;                                                      
                                                                              
main(argc,argv)                                                               
   int argc;                                                                  
   char *argv[];                                                              
{                                                                             
   if (argc == 3 && isdigit(*argv[2])) sscanf(argv[2],"%d",&memsize);         
   if (argc >= 2 && *argv[1] == 'c') comp();                                  
   else if (argc >= 2 && *argv[1] == 'e') exp();                              
   else {                                                                     
      fprintf(stderr,"usage:  dmc [ce] memsize outfile\n");          
      exit(1);                                                                
   }                                                                          
   return 0;                                                                  
}                                                                             
                                                                              
comp(){                                                                       
   int max = 0x1000000,                                                       
       min = 0,                                                               
       mid,                                                                   
       c,i,                                                                   
       inbytes = 0,                                                           
       outbytes =3,                                                           
       pout = 3,                                                              
       bit;                                                                   
                                                                              
   pinit(memsize);                                                            
                                                                              
   for(;;){                                                                   
      c = getchar();                                                          
      if (c == EOF) {                                                         
         min = max-1;                                                         
        fprintf(stderr,"compress done: bytes in %d, bytes out %d, ratio %f\n",
                 inbytes,outbytes,(float)outbytes/inbytes);                   
         break;                                                               
      }                                                                       
      for (i=0;i<8;i++){                                                      
         bit = (c << i) & 0x80;                                               
         mid = min + (max-min-1) * predict();                                 
         pupdate(bit != 0);                                                   
         if (mid == min) mid++;                                               
         if (mid == (max-1)) mid--;                                           
                                                                              
         if (bit) {                                                           
            min = mid;                                                        
         } else {                                                             
            max = mid;                                                        
         }                                                                    
         while ((max-min) < 256) {                                            
            if(bit)max--;                                                     
            putchar(min >> 16);                                               
            outbytes++;                                                       
            min = (min << 8) & 0xffff00;                                      
            max = ((max << 8) & 0xffff00 ) ;                                  
            if (min >= max) max = 0x1000000;                                  
         }                                                                    
      }                                                                       
      if(!(++inbytes & 0xff)){                                                
         if(!(inbytes & 0xffff)){                                             
               fprintf(stderr,                                                
                       "compressing... bytes in %d, bytes out %d, ratio %f\r",
                       inbytes,outbytes,(float)outbytes/inbytes);             
         }                                                                    
         if (outbytes - pout > 256) { /* compression failing */               
            pflush();                                                         
         }                                                                    
         pout = outbytes;                                                     
      }                                                                       
   }                                                                          
   putchar(min>>16);                                                          
   putchar((min>>8) & 0xff);                                                  
   putchar(min & 0x00ff);                                                     
}                                                                             
                                                                              
                                                                              
exp(){                                                                        
   int max = 0x1000000,                                                       
       min = 0,                                                               
       mid,                                                                   
       val,                                                                   
       i,                                                                     
       inbytes=3,                                                             
       pin=3,                                                                 
       outbytes=0,                                                            
       bit,                                                                   
       c;                                                                     
                                                                              
   pinit(memsize);                                                            
                                                                              
   val = getchar()<<16;                                                       
   val += getchar()<<8;                                                       
   val += getchar();                                                          
   while(1) {                                                                 
      c = 0;                                                                  
      if (val == (max-1)) {                                                   
         fprintf(stderr,"expand: input %d output %d\n",inbytes,outbytes);     
         break;                                                               
      }                                                                       
      for (i=0;i<8;i++){                                                      
         mid = min + (max-min-1)*predict();                                   
         if (mid == min) mid++;                                               
         if (mid == (max-1)) mid--;                                           
         if (val >= mid) {                                                    
            bit = 1;                                                          
            min = mid;                                                        
         } else {                                                             
            bit = 0;                                                          
            max = mid;                                                        
         }                                                                    
         pupdate(bit != 0);                                                   
         c = c + c + bit;                                                     
         while ((max-min) < 256) {                                            
            if(bit)max--;                                                     
            inbytes++;                                                        
            val = (val << 8) & 0xffff00 | (getchar()& 0xff);                  
            min = (min << 8) & 0xffff00;                                      
            max = ((max << 8) & 0xffff00 ) ;                                  
            if (min >= max) max = 0x1000000;                                  
         }                                                                    
      }                                                                       
      putchar(c);                                                             
      if(!(++outbytes & 0xff)){                                               
         if (inbytes - pin > 256) { /* compression was failing */             
            pflush();                                                         
         }                                                                    
         pin = inbytes;                                                       
      }                                                                       
   }                                                                          
}                                                                             
                                                                              
typedef struct nnn {                                                          
           float count[2];                                                    
           struct nnn    *next[2];                                            
} node;                                                                       
                                                                              
static int threshold = 2, bigthresh = 2;                                      
                                                                              
static node *p, *new, nodes[256][256];                                        
                                                                              
static node *nodebuf;                                                         
static node *nodemaxp;                                                        
static node *nodesptr;                                                        
                                                                              
#include <malloc.h>
                                                                              
pinit(memsize)                                                                
   int memsize;                                                               
{                                                                             
   fprintf(stderr,"using %d bytes of predictor memory\n",memsize);            
   nodebuf = (node *) malloc (memsize);                                       
   if (nodebuf == (node *) NULL) {                                            
      fprintf(stderr,"memory alloc failed; try smaller predictor memory\n");  
      exit(1);                                                                
   }                                                                          
   nodemaxp = nodebuf + (memsize/sizeof(node)) - 20;                          
   pflush();                                                                  
}                                                                             
                                                                              
pflush(){                                                                     
   int i,j;                                                                   
   for (j=0;j<256;j++){                                                       
      for (i=0;i<127;i++) {                                                   
         nodes[j][i].count[0] = 0.2;                                          
         nodes[j][i].count[1] = 0.2;                                          
         nodes[j][i].next[0] = &nodes[j][2*i + 1];                            
         nodes[j][i].next[1] = &nodes[j][2*i + 2];                            
      }                                                                       
      for (i=127;i<255;i++) {                                                 
         nodes[j][i].count[0] = 0.2;                                          
         nodes[j][i].count[1] = 0.2;                                          
         nodes[j][i].next[0] = &nodes[i+1][0];                                
         nodes[j][i].next[1] = &nodes[i-127][0];                              
      }                                                                       
   }                                                                          
   nodesptr = nodebuf;                                                        
   preset();                                                                  
}                                                                             
                                                                              
preset(){                                                                     
   p = &nodes[0][0];                                                          
}                                                                             
                                                                              
float predict(){                                                              
   return   p->count[0] / (p->count[0] + p->count[1]);                        
}                                                                             
                                                                              
pupdate(b)                                                                    
   int b;                                                                     
{                                                                             
   float r;                                                                   
   if (p->count[b] >= threshold &&                                            
      p->next[b]->count[0]+p->next[b]->count[1]                               
       >= bigthresh + p->count[b]){                                           
      new = nodesptr++;                                                       
      p->next[b]->count[0] -= new->count[0] =                                 
         p->next[b]->count[0] *                                               
         (r = p->count[b]/(p->next[b]->count[1]+p->next[b]->count[0]));       
      p->next[b]->count[1] -= new->count[1] =                                 
         p->next[b]->count[1] * r;                                            
      new->next[0] = p->next[b]->next[0];                                     
      new->next[1] = p->next[b]->next[1];                                     
      p->next[b] = new;                                                       
   }                                                                          
   p->count[b]++;                                                             
   p = p->next[b];                                                            
   if (nodesptr > nodemaxp){                                                  
      fprintf(stderr,"flushing ...\n");                                       
      pflush();                                                               
   }                                                                          
}                                                                             

Изменена 19.03.2011 06:45 MSK Яндекс цитирования Рейтинг@Mail.ru   quill