Merge pull request #427 from Zip753/master

Add backward reference research tools to brotli repository.

.gitmodules

@@ -1,3 +1,6 @@
 [submodule "terryfy"]
 	path = terryfy
 	url = https://github.com/MacPython/terryfy.git
+[submodule "research/esaxx"]
+	path = research/esaxx
+	url = https://github.com/hillbig/esaxx

research/Makefile

@@ -0,0 +1,17 @@
+CC = g++
+CFLAGS += -O2
+CPPFLAGS += -std=c++11
+SOURCES = $(wildcard *.cc)
+EXECUTABLES = $(SOURCES:.cc=)
+BINDIR = bin
+
+all: $(EXECUTABLES)
+
+$(BINDIR):
+	mkdir -p $@
+
+$(EXECUTABLES): $(BINDIR)
+	$(CC) $(CFLAGS) $(CPPFLAGS) $(addsuffix .cc, $@) -o $(BINDIR)/$@
+
+clean:
+	rm -rf $(BINDIR)

research/README.md

@@ -0,0 +1,67 @@
+## Introduction
+
+In this directory we publish simple tools to analyze backward reference distance distributions in LZ77 compression. We developed these tools to be able to make more efficient encoding of distances in large-window brotli. In large-window compression the average cost of a backward reference distance is higher, and this may allow for more advanced encoding strategies, such as delta coding or an increase in context size, to bring significant compression density improvements. Our tools visualize the backward references as histogram images, i.e., one pixel in the image shows how many distances of a certain range exist at a certain locality in the data. The human visual system is excellent at pattern detection, so we tried to roughly identify patterns visually before going into more quantitative analysis. These tools can turn out to be useful in development of  other LZ77-based compressors and we hope you try them out.
+
+
+## Tools
+### find\_opt\_references
+
+This tool generates optimal (match-length-wise) backward references for every position in the input files and stores them in `*.dist` file described below.
+
+Example usage:
+
+    find_opt_references input.txt output.dist
+
+### draw\_histogram
+
+This tool generates a visualization of the distribution of backward references stored in `*.dist` file. The original file size has to be specified as a second parameter. The output is a grayscale PGM (binary) image.
+
+Example usage:
+
+    draw_histogram input.dist 65536 output.pgm
+
+Here's an example of resulting image:
+
+![](img/enwik9_brotli.png)
+
+### draw\_diff
+
+This tool generates a diff PPM (binary) image between two input 8-bit PGM (binary) images. Input images must be of same size. Useful for comparing different backward references distributions for same input file. Normally used for comparison of output images from `draw_histogram` tool.
+
+Example usage:
+
+    draw_diff image1.pgm image2.pgm diff.ppm
+
+For example the diff of this image
+
+![](img/enwik9_brotli.png)
+
+and this image
+
+![](img/enwik9_opt.png)
+
+looks like this:
+
+![](img/enwik9_diff.png)
+
+
+## Backward distance file format
+
+The format of `*.dist` files is as follows:
+
+    [[     0| match length][     1|position|distance]...]
+     [1 byte|      4 bytes][1 byte| 4 bytes| 4 bytes]
+
+More verbose explanation: for each backward reference there is a position-distance pair, also a copy length may be specified. Copy length is prefixed with flag byte 0, position-distance pair is prefixed with flag byte 1. Each number is a 32-bit integer. Copy length always comes before position-distance pair. Standalone copy length is allowed, in this case it is ignored.
+
+Here's an example of how to read from `*.dist` file:
+
+```c++
+#include "read_dist.h"
+
+FILE* f;
+int copy, pos, dist;
+while (ReadBackwardReference(fin, &copy, &pos, &dist)) {
+   ...
+}
+```

research/draw_diff.cc

@@ -0,0 +1,104 @@
+/* Copyright 2016 Google Inc. All Rights Reserved.
+   Author: vanickulin@google.com (Ivan Nikulin)
+
+   Distributed under MIT license.
+   See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
+*/
+
+/* Tool for drawing diff PPM images between two input PGM images. Normally used
+   with backward reference histogram drawing tool. */
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <cstdint>
+#include <cstdio>
+#include <vector>
+
+void ReadPGM(FILE* f, uint8_t*** image, size_t* height, size_t* width) {
+  int colors;
+  assert(fscanf(f, "P5\n%lu %lu\n%d\n", width, height, &colors) == 3);
+  assert(colors == 255);
+  *image = new uint8_t*[*height];
+  for (int i = *height - 1; i >= 0; --i) {
+    (*image)[i] = new uint8_t[*width];
+    assert(fread((*image)[i], 1, *width, f) == *width);
+  }
+}
+
+void CalculateDiff(int** diff, uint8_t** image1, uint8_t** image2,
+                   size_t height, size_t width) {
+  for (size_t i = 0; i < height; ++i) {
+    for (size_t j = 0; j < width; ++j) {
+      diff[i][j] = static_cast<int>(image1[i][j]) - image2[i][j];
+    }
+  }
+}
+
+void DrawDiff(int** diff, uint8_t** image1, uint8_t** image2,
+              size_t height, size_t width, FILE* f) {
+  int max = -1234;
+  int min = +1234;
+  for (size_t i = 0; i < height; ++i) {
+    for (size_t j = 0; j < width; ++j) {
+      if (max < diff[i][j]) max = diff[i][j];
+      if (min > diff[i][j]) min = diff[i][j];
+      int img_min = std::min(255 - image1[i][j], 255 - image2[i][j]);
+      if (max < img_min) max = img_min;
+    }
+  }
+
+  int abs_max = -min;
+  if (abs_max < max) abs_max = max;
+
+  fprintf(f, "P6\n%lu %lu\n%d\n", width, height, abs_max);
+
+  uint8_t* row = new uint8_t[3 * width];
+  for (int i = height - 1; i >= 0; --i) {
+    for (int j = 0; j < width; ++j) {
+      int min_val = std::min(255 - image1[i][j], 255 - image2[i][j]);
+      int max_val = std::max(min_val, abs(diff[i][j]));
+      if (diff[i][j] > 0) { /* red */
+        row[3 * j + 0] = abs_max - max_val + diff[i][j];
+        row[3 * j + 1] = abs_max - max_val;
+        row[3 * j + 2] = abs_max - max_val + min_val;
+      } else { /* green */
+        row[3 * j + 0] = abs_max - max_val;
+        row[3 * j + 1] = abs_max - max_val - diff[i][j];
+        row[3 * j + 2] = abs_max - max_val + min_val;
+      }
+    }
+    fwrite(row, 1, 3 * width, f);
+  }
+  delete[] row;
+}
+
+int main(int argc, char* argv[]) {
+  if (argc != 4) {
+    printf("usage: %s pgm1 pgm2 diff_ppm_path\n", argv[0]);
+    return 1;
+  }
+
+  FILE* fimage1 = fopen(argv[1], "rb");
+  FILE* fimage2 = fopen(argv[2], "rb");
+  FILE* fdiff = fopen(argv[3], "wb");
+
+  uint8_t **image1, **image2;
+  size_t h1, w1, h2, w2;
+  ReadPGM(fimage1, &image1, &h1, &w1);
+  ReadPGM(fimage2, &image2, &h2, &w2);
+  fclose(fimage1);
+  fclose(fimage2);
+  if (!(h1 == h2 && w1 == w2)) {
+    printf("Images must have the same size.\n");
+    return 1;
+  }
+
+  int** diff = new int*[h1];
+  for (size_t i = 0; i < h1; ++i) diff[i] = new int[w1];
+  CalculateDiff(diff, image1, image2, h1, w1);
+
+  DrawDiff(diff, image1, image2, h1, w1, fdiff);
+  fclose(fdiff);
+  return 0;
+}

research/draw_histogram.cc

@@ -0,0 +1,197 @@
+/* Copyright 2016 Google Inc. All Rights Reserved.
+   Author: vanickulin@google.com (Ivan Nikulin)
+
+   Distributed under MIT license.
+   See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
+*/
+
+/* Backward reference visualization tool. Accepts file with backward references
+   as an input and produces PGM image with histogram of those references. */
+
+#include <algorithm> /* min */
+#include <cassert>
+#include <cstring> /* memset */
+#include <cmath> /* log, round */
+#include <cstdio> /* fscanf, fprintf */
+#include <cstdint>
+
+#include "./read_dist.h"
+
+const int FLAGS_height = 1000;  // Height of the resulting histogam.
+const int FLAGS_width = 1000;  // Width of the resulting histogam.
+int FLAGS_size;  // Size of the compressed file.
+const int FLAGS_brotli_window = 0;  // Size of brotli window in bits.
+const uint64_t FLAGS_min_distance = 0;  // Minimum distance.
+uint64_t FLAGS_max_distance = 0;  // Maximum distance.
+const bool FLAGS_with_copies = false;  // True if input contains copy length.
+const bool FLAGS_simple = false;  // True if using only black and white pixels.
+const bool FLAGS_linear = true;  // True if using linear distance mapping.
+const uint64_t FLAGS_skip = 0;  // Number of bytes to skip.
+
+inline double DistanceTransform(double x) {
+  static bool linear = FLAGS_linear;
+  if (linear) {
+    return x;
+  } else {
+    /* Using log^2 scale because log scale produces big white gap at the bottom
+       of image. */
+    return log(x) * log(x);
+  }
+}
+
+/* Mapping pixel density on arc function to increase contrast. */
+inline double DensityTransform(double x) {
+  double z = 255 - x;
+  return sqrt(255 * 255 - z * z);
+}
+
+inline int GetMaxDistance() {
+  return FLAGS_max_distance;
+}
+
+void AdjustPosition(int* pos) {
+  static uint32_t offset = 0;
+  static int last = 0;
+  static uint32_t window_size = (1 << FLAGS_brotli_window);
+  assert(*pos >= 0 && *pos < window_size);
+  if (*pos < last) {
+    offset += window_size;
+  }
+  last = *pos;
+  *pos += offset;
+}
+
+void BuildHistogram(FILE* fin, int** histo) {
+  int height = FLAGS_height;
+  int width = FLAGS_width;
+  int skip = FLAGS_skip;
+  size_t min_distance = FLAGS_min_distance;
+
+  printf("height = %d, width = %d\n", height, width);
+
+  for (int i = 0; i < height; i++) {
+    for (int j = 0; j < width; j++) {
+      histo[i][j] = 0;
+    }
+  }
+
+  int max_pos = FLAGS_size - skip;
+  double min_dist = min_distance > 0 ? DistanceTransform(min_distance) : 0;
+  double max_dist = DistanceTransform(GetMaxDistance()) - min_dist;
+  int copy, pos, distance, x, y;
+  double dist;
+  while (ReadBackwardReference(fin, &copy, &pos, &distance)) {
+    if (pos == -1) continue;  // In case when only insert is present.
+    if (distance < min_distance || distance >= GetMaxDistance()) continue;
+    if (FLAGS_brotli_window != 0) {
+      AdjustPosition(&pos);
+    }
+    if (pos >= skip && distance <= pos) {
+      pos -= skip;
+      if (pos >= max_pos) break;
+      dist = DistanceTransform(static_cast<double>(distance)) - min_dist;
+
+      x = std::min(static_cast<int>(round(dist / max_dist * height)),
+                   height - 1);
+      y = 1ul * pos * width / max_pos;
+      if (!(y >= 0 && y < width)) {
+        printf("pos = %d, max_pos = %d, y = %d\n", pos, max_pos, y);
+        assert(y >= 0 && y < width);
+      }
+
+      if (FLAGS_with_copies) {
+        int right = 1ul * (pos + copy - 1) * width / max_pos;
+        if (right < 0) {
+          printf("pos = %d, distance = %d, copy = %d, y = %d, right = %d\n",
+                  pos, distance, copy, y, right);
+          assert(right >= 0);
+        }
+        if (y == right) {
+          histo[x][y] += copy;
+        } else {
+          int pos2 = static_cast<int>(ceil(1.0 * (y + 1) * max_pos / width));
+          histo[x][y] += pos2 - pos;
+          for (int i = y + 1; i < right && i < width; ++i) {
+            histo[x][i] += max_pos / width;  // Sometimes 1 more, but who cares.
+          }
+          // Make sure the match doesn't go beyond the image.
+          if (right < width) {
+            pos2 = static_cast<int>(ceil(1.0 * right * max_pos / width));
+            histo[x][right] += pos + copy - 1 - pos2 + 1;
+          }
+        }
+      } else {
+        histo[x][y]++;
+      }
+    }
+  }
+}
+
+void ConvertToPixels(int** histo, uint8_t** pixel) {
+  int height = FLAGS_height;
+  int width = FLAGS_width;
+
+  int maxs = 0;
+  for (int i = 0; i < height; i++) {
+    for (int j = 0; j < width; j++) {
+      if (maxs < histo[i][j]) maxs = histo[i][j];
+    }
+  }
+
+  bool simple = FLAGS_simple;
+  double max_histo = static_cast<double>(maxs);
+  for (int i = 0; i < height; i++) {
+    for (int j = 0; j < width; j++) {
+      if (simple) {
+        pixel[i][j] = histo[i][j] > 0 ? 0 : 255;
+      } else {
+        pixel[i][j] = static_cast<uint8_t>(
+            255 - DensityTransform(histo[i][j] / max_histo * 255));
+      }
+    }
+  }
+}
+
+void DrawPixels(uint8_t** pixel, FILE* fout) {
+  int height = FLAGS_height;
+  int width = FLAGS_width;
+
+  fprintf(fout, "P5\n%d %d\n255\n", width, height);
+  for (int i = height - 1; i >= 0; i--) {
+    fwrite(pixel[i], 1, width, fout);
+  }
+}
+
+int main(int argc, char* argv[]) {
+  if (argc != 4) {
+    printf("usage: draw_histogram.cc dist_file input_size output_file\n");
+    return 1;
+  }
+
+  int height = FLAGS_height;
+  int width = FLAGS_width;
+
+  FILE* fin = fopen(argv[1], "r");
+  FILE* fout = fopen(argv[3], "wb");
+
+  FLAGS_size = atoi(argv[2]);
+
+  if (FLAGS_max_distance == 0) FLAGS_max_distance = FLAGS_size;
+
+  uint8_t** pixel = new uint8_t*[height];
+  int** histo = new int*[height];
+  for (int i = 0; i < height; i++) {
+    pixel[i] = new uint8_t[width];
+    histo[i] = new int[width];
+  }
+
+  BuildHistogram(fin, histo);
+  fclose(fin);
+
+  ConvertToPixels(histo, pixel);
+
+  DrawPixels(pixel, fout);
+  fclose(fout);
+
+  return 0;
+}

research/esaxx

@@ -0,0 +1 @@
+Subproject commit ca7cb332011ec37a8436487f210f396b84bd8273

research/find_opt_references.cc

@@ -0,0 +1,147 @@
+/* Copyright 2016 Google Inc. All Rights Reserved.
+   Author: vanickulin@google.com (Ivan Nikulin)
+
+   Distributed under MIT license.
+   See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
+*/
+
+/* Tool for generating optimal backward references for the input file. Uses
+   sais-lite library for building suffix array. */
+
+#include <algorithm>
+#include <cassert>
+#include <cstdio>
+#include <cstring>
+
+#include "./esaxx/sais.hxx"
+
+const int FLAGS_min_length = 1;  // Minimal length of found backward references.
+const int FLAGS_skip = 1;  // Number of bytes to skip.
+
+const size_t kFileBufferSize = (1 << 16);  // 64KB
+
+typedef int sarray_type;  // Can't make it unsigned because of templates :(
+typedef uint8_t input_type;
+typedef uint32_t lcp_type;
+
+void ReadInput(FILE* fin, input_type* storage, size_t input_size) {
+  size_t last_pos = 0;
+  size_t available_in;
+  fseek(fin, 0, SEEK_SET);
+  do {
+    available_in = fread(storage + last_pos, 1, kFileBufferSize, fin);
+    last_pos += available_in;
+  } while (available_in != 0);
+  assert(last_pos == input_size);
+}
+
+void BuildLCP(input_type* storage, sarray_type* sarray, lcp_type* lcp,
+              size_t size, uint32_t* pos) {
+  for (int i = 0; i < size; ++i) {
+    pos[sarray[i]] = i;
+  }
+  uint32_t k = 0;
+  lcp[size - 1] = 0;
+  for (int i = 0; i < size; ++i) {
+    if (pos[i] == size - 1) {
+      k = 0;
+      continue;
+    }
+    uint32_t j = sarray[pos[i] + 1];  // Suffix which follow i-th suffix in SA.
+    while (i + k < size && j + k < size && storage[i + k] == storage[j + k]) {
+      ++k;
+    }
+    lcp[pos[i]] = k;
+    if (k > 0) --k;
+  }
+}
+
+void ProcessReferences(input_type* storage, sarray_type* sarray, lcp_type* lcp,
+                       size_t size, uint32_t* pos, FILE* fout) {
+  int min_length = FLAGS_min_length;
+  for (int idx = FLAGS_skip; idx < size; ++idx) {
+    int max_lcp = -1;
+    int max_lcp_ix;
+    int left_lcp = -1;
+    int left_ix;
+    for (left_ix = pos[idx] - 1; left_ix >= 0; --left_ix) {
+      if (left_lcp == -1 || left_lcp > lcp[left_ix]) {
+        left_lcp = lcp[left_ix];
+      }
+      if (left_lcp == 0) break;
+      if (sarray[left_ix] < idx) break;
+    }
+    if (left_ix >= 0) {
+      max_lcp = left_lcp;
+      max_lcp_ix = left_ix;
+    }
+
+    int right_lcp = -1;
+    int right_ix;
+    for (right_ix = pos[idx]; right_ix < size - 1; ++right_ix) {
+      if (right_lcp == -1 || right_lcp > lcp[right_ix]) {
+        right_lcp = lcp[right_ix];
+      }
+      // Stop if we have better result from the left side already.
+      if (right_lcp < max_lcp) break;
+      if (right_lcp == 0) break;
+      if (sarray[right_ix] < idx) break;
+    }
+    if (right_lcp > max_lcp && right_ix < size - 1) {
+      max_lcp = right_lcp;
+      max_lcp_ix = right_ix;
+    }
+
+    if (max_lcp >= min_length) {
+      int dist = idx - sarray[max_lcp_ix];
+      if (dist <= 0) {
+        printf("idx = %d, pos[idx] = %u\n", idx, pos[idx]);
+        printf("left_ix = %d, right_ix = %d\n",
+                left_ix, right_ix);
+        printf("left_lcp = %d, right_lcp = %d\n",
+                left_lcp, right_lcp);
+        printf("sarray[left_ix] = %d, sarray[right_ix] = %d\n",
+                sarray[left_ix], sarray[right_ix]);
+        assert(dist > 0);
+      }
+      fputc(1, fout);
+      fwrite(&idx, sizeof(int), 1, fout);  // Position in input.
+      fwrite(&dist, sizeof(int), 1, fout);  // Backward distance.
+    }
+  }
+}
+
+int main(int argc, char* argv[]) {
+  if (argc != 3) {
+    printf("usage: %s input_file output_file\n", argv[0]);
+    return 1;
+  }
+
+  FILE* fin = fopen(argv[1], "rb");
+  FILE* fout = fopen(argv[2], "w");
+
+  fseek(fin, 0, SEEK_END);
+  int input_size = ftell(fin);
+  fseek(fin, 0, SEEK_SET);
+  printf("The file size is %u bytes\n", input_size);
+
+  input_type* storage = new input_type[input_size];
+
+  ReadInput(fin, storage, input_size);
+  fclose(fin);
+
+  sarray_type* sarray = new sarray_type[input_size];
+  saisxx(storage, sarray, input_size);
+  printf("Suffix array calculated.\n");
+
+  // Inverse suffix array.
+  uint32_t* pos = new uint32_t[input_size];
+
+  lcp_type* lcp = new lcp_type[input_size];
+  BuildLCP(storage, sarray, lcp, input_size, pos);
+  printf("LCP array constructed.\n");
+
+  ProcessReferences(storage, sarray, lcp, input_size, pos, fout);
+  fclose(fout);
+  return 0;
+}

research/read_dist.h

@@ -0,0 +1,40 @@
+/* Copyright 2016 Google Inc. All Rights Reserved.
+   Author: vanickulin@google.com (Ivan Nikulin)
+
+   Distributed under MIT license.
+   See file LICENSE for detail or copy at https://opensource.org/licenses/MIT
+*/
+
+/* API for reading distances from *.dist file. */
+
+#include <cassert>
+#include <cstdio>
+
+#ifndef BROTLI_RESEARCH_READ_DIST_H_
+#define BROTLI_RESEARCH_READ_DIST_H_
+
+/* Reads backwards reference from .dist file. Sets all missing fields to -1.
+   Returns false when EOF is met or input is corrupt. */
+bool ReadBackwardReference(FILE* fin, int* copy, int* pos, int* dist) {
+  int c = getc(fin);
+  if (c == EOF) return false;
+  if (c == 0) {
+    assert(fread(copy, sizeof(int), 1, fin) == 1);
+    if ((c = getc(fin)) != 1) {
+      ungetc(c, fin);
+      *pos = *dist = -1;
+    } else {
+      assert(fread(pos, sizeof(int), 1, fin) == 1);
+      assert(fread(dist, sizeof(int), 1, fin) == 1);
+    }
+  } else if (c != 1) {
+    return false;
+  } else {
+    assert(fread(pos, sizeof(int), 1, fin) == 1);
+    assert(fread(dist, sizeof(int), 1, fin) == 1);
+    *copy = -1;
+  }
+  return true;
+}
+
+#endif  /* BROTLI_RESEARCH_READ_DIST_H_ */