Add an implementation of Asymetric Numeral Systems (ANS).

Change-Id: Ie41bc72127e700887566dcc951da9d83a0b94891

configure

@@ -281,6 +281,7 @@ EXPERIMENT_LIST="
     ext_interp
     ext_refs
     supertx
+    ans
 "
 CONFIG_LIST="
     dependency_tracking

test/test.mk

@@ -165,6 +165,7 @@ ifeq ($(CONFIG_VP10),yes)
 
 LIBVPX_TEST_SRCS-yes                    += vp10_inv_txfm_test.cc
 LIBVPX_TEST_SRCS-$(CONFIG_VP10_ENCODER) += vp10_dct_test.cc
+LIBVPX_TEST_SRCS-$(CONFIG_ANS)          += vp10_ans_test.cc
 
 endif # VP10
 

test/vp10_ans_test.cc

+/*
+ *  Copyright (c) 2015 The WebM project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+#include <ctime>
+#include <utility>
+#include <vector>
+
+#include "third_party/googletest/src/include/gtest/gtest.h"
+
+#include "test/acm_random.h"
+#include "vp10/common/ans.h"
+#include "vp10/encoder/treewriter.h"
+#include "vpx_dsp/bitreader.h"
+#include "vpx_dsp/bitwriter.h"
+
+namespace {
+typedef std::vector<std::pair<uint8_t, bool> > PvVec;
+
+PvVec abs_encode_build_vals(int iters) {
+  PvVec ret;
+  libvpx_test::ACMRandom gen(0x30317076);
+  double entropy = 0;
+  for (int i = 0; i < iters; ++i) {
+    uint8_t p;
+    do {
+      p = gen.Rand8();
+    } while (p == 0);  // zero is not a valid coding probability
+    bool b = gen.Rand8() < p;
+    ret.push_back(std::make_pair(static_cast<uint8_t>(p), b));
+    double d = p / 256.;
+    entropy += -d * log2(d) - (1 - d) * log2(1 - d);
+  }
+  printf("entropy %f\n", entropy);
+  return ret;
+}
+
+bool check_rabs(const PvVec &pv_vec, uint8_t *buf) {
+  AnsCoder a;
+  ans_write_init(&a, buf);
+
+  std::clock_t start = std::clock();
+  for (PvVec::const_reverse_iterator it = pv_vec.rbegin(); it != pv_vec.rend();
+       ++it) {
+    rabs_write(&a, it->second, 256 - it->first);
+  }
+  std::clock_t enc_time = std::clock() - start;
+  int offset = ans_write_end(&a);
+  bool okay = true;
+  AnsDecoder d;
+  if (ans_read_init(&d, buf, offset)) return false;
+  start = std::clock();
+  for (PvVec::const_iterator it = pv_vec.begin(); it != pv_vec.end(); ++it) {
+    okay &= rabs_read(&d, 256 - it->first) == it->second;
+  }
+  std::clock_t dec_time = std::clock() - start;
+  if (!okay) return false;
+  printf("rABS size %d enc_time %f dec_time %f\n", offset,
+         static_cast<float>(enc_time) / CLOCKS_PER_SEC,
+         static_cast<float>(dec_time) / CLOCKS_PER_SEC);
+  return ans_read_end(&d);
+}
+
+bool check_rabs_asc(const PvVec &pv_vec, uint8_t *buf) {
+  AnsCoder a;
+  ans_write_init(&a, buf);
+
+  std::clock_t start = std::clock();
+  for (PvVec::const_reverse_iterator it = pv_vec.rbegin(); it != pv_vec.rend();
+       ++it) {
+    rabs_asc_write(&a, it->second, 256 - it->first);
+  }
+  std::clock_t enc_time = std::clock() - start;
+  int offset = ans_write_end(&a);
+  bool okay = true;
+  AnsDecoder d;
+  if (ans_read_init(&d, buf, offset)) return false;
+  start = std::clock();
+  for (PvVec::const_iterator it = pv_vec.begin(); it != pv_vec.end(); ++it) {
+    okay &= rabs_asc_read(&d, 256 - it->first) == it->second;
+  }
+  std::clock_t dec_time = std::clock() - start;
+  if (!okay) return false;
+  printf("rABS (asc) size %d enc_time %f dec_time %f\n", offset,
+         static_cast<float>(enc_time) / CLOCKS_PER_SEC,
+         static_cast<float>(dec_time) / CLOCKS_PER_SEC);
+  return ans_read_end(&d);
+}
+
+bool check_uabs(const PvVec &pv_vec, uint8_t *buf) {
+  AnsCoder a;
+  ans_write_init(&a, buf);
+
+  std::clock_t start = std::clock();
+  for (PvVec::const_reverse_iterator it = pv_vec.rbegin(); it != pv_vec.rend();
+       ++it) {
+    uabs_write(&a, it->second, 256 - it->first);
+  }
+  std::clock_t enc_time = std::clock() - start;
+  int offset = ans_write_end(&a);
+  bool okay = true;
+  AnsDecoder d;
+  if (ans_read_init(&d, buf, offset)) return false;
+  start = std::clock();
+  for (PvVec::const_iterator it = pv_vec.begin(); it != pv_vec.end(); ++it) {
+    okay &= uabs_read(&d, 256 - it->first) == it->second;
+  }
+  std::clock_t dec_time = std::clock() - start;
+  if (!okay) return false;
+  printf("uABS size %d enc_time %f dec_time %f\n", offset,
+         static_cast<float>(enc_time) / CLOCKS_PER_SEC,
+         static_cast<float>(dec_time) / CLOCKS_PER_SEC);
+  return ans_read_end(&d);
+}
+
+bool check_vpxbool(const PvVec &pv_vec, uint8_t *buf) {
+  vpx_writer w;
+  vpx_reader r;
+  vpx_start_encode(&w, buf);
+
+  std::clock_t start = std::clock();
+  for (PvVec::const_iterator it = pv_vec.begin(); it != pv_vec.end(); ++it) {
+    vpx_write(&w, it->second, 256 - it->first);
+  }
+  std::clock_t enc_time = std::clock() - start;
+  vpx_stop_encode(&w);
+  bool okay = true;
+  vpx_reader_init(&r, buf, w.pos, NULL, NULL);
+  start = std::clock();
+  for (PvVec::const_iterator it = pv_vec.begin(); it != pv_vec.end(); ++it) {
+    okay &= vpx_read(&r, 256 - it->first) == it->second;
+  }
+  std::clock_t dec_time = std::clock() - start;
+  printf("VPX size %d enc_time %f dec_time %f\n", w.pos,
+         static_cast<float>(enc_time) / CLOCKS_PER_SEC,
+         static_cast<float>(dec_time) / CLOCKS_PER_SEC);
+  return okay;
+}
+
+const rans_sym rans_sym_tab[] = {
+    {70, 186}, {70, 116}, {100, 16}, {16, 0},
+};
+const int kDistinctSyms = sizeof(rans_sym_tab) / sizeof(rans_sym_tab[0]);
+
+std::vector<int> ans_encode_build_vals(const rans_sym *tab, int iters) {
+  std::vector<int> p_to_sym;
+  int i = 0;
+  while (p_to_sym.size() < 256) {
+    p_to_sym.insert(p_to_sym.end(), tab[i].prob, i);
+    ++i;
+  }
+  assert(p_to_sym.size() == 256);
+  std::vector<int> ret;
+  libvpx_test::ACMRandom gen(18543637);
+  for (int i = 0; i < iters; ++i) {
+    int sym = p_to_sym[gen.Rand8()];
+    ret.push_back(sym);
+  }
+  return ret;
+}
+
+void rans_build_dec_tab(const struct rans_sym sym_tab[],
+                        rans_dec_lut dec_tab) {
+  int val = 0;
+  int i;
+  for (i = ans_p8_precision - 1; i >= 0; --i) {
+    dec_tab[i].val = val;
+    dec_tab[i].prob = sym_tab[val].prob;
+    dec_tab[i].cum_prob = sym_tab[val].cum_prob;
+    if (i == sym_tab[val].cum_prob) ++val;
+  }
+}
+
+bool check_rans(const std::vector<int> &sym_vec, const rans_sym *const tab,
+                uint8_t *buf) {
+  AnsCoder a;
+  ans_write_init(&a, buf);
+  rans_dec_lut dec_tab;
+  rans_build_dec_tab(tab, dec_tab);
+
+  std::clock_t start = std::clock();
+  for (std::vector<int>::const_reverse_iterator it = sym_vec.rbegin();
+       it != sym_vec.rend(); ++it) {
+    rans_write(&a, &tab[*it]);
+  }
+  std::clock_t enc_time = std::clock() - start;
+  int offset = ans_write_end(&a);
+  bool okay = true;
+  AnsDecoder d;
+  if (ans_read_init(&d, buf, offset)) return false;
+  start = std::clock();
+  for (std::vector<int>::const_iterator it = sym_vec.begin();
+       it != sym_vec.end(); ++it) {
+    okay &= rans_read(&d, dec_tab) == *it;
+  }
+  std::clock_t dec_time = std::clock() - start;
+  if (!okay) return false;
+  printf("rANS size %d enc_time %f dec_time %f\n", offset,
+         static_cast<float>(enc_time) / CLOCKS_PER_SEC,
+         static_cast<float>(dec_time) / CLOCKS_PER_SEC);
+  return ans_read_end(&d);
+}
+
+void build_tree(vpx_tree_index *tree, int num_syms) {
+  vpx_tree_index i;
+  int sym = 0;
+  for (i = 0; i < num_syms - 1; ++i) {
+    tree[2 * i] = sym--;
+    tree[2 * i + 1] = 2 * (i + 1);
+  }
+  tree[2 * i - 1] = sym;
+}
+
+// treep are the probabilites of tree nodes like:
+//          *
+//         / \
+//    -sym0  *
+//          / \
+//     -sym1  *
+//           / \
+//      -sym2  -sym3
+void tab2tree(const rans_sym *tab, int tab_size, vpx_prob *treep) {
+  const unsigned basep = 256;
+  unsigned pleft = basep;
+  for (int i = 0; i < tab_size - 1; ++i) {
+    unsigned prob = (tab[i].prob * basep + (basep / 2)) / pleft;
+    assert(prob > 0 && prob < 256);
+    treep[i] = prob;
+    pleft -= tab[i].prob;
+  }
+}
+
+struct sym_bools {
+  unsigned bits;
+  int len;
+};
+
+static void make_tree_bits_tab(sym_bools *tab, int num_syms) {
+  unsigned bits = 0;
+  int len = 0;
+  int i;
+  for (i = 0; i < num_syms - 1; ++i) {
+    bits *= 2;
+    ++len;
+    tab[i].bits = bits;
+    tab[i].len = len;
+    ++bits;
+  }
+  tab[i].bits = bits;
+  tab[i].len = len;
+}
+
+void build_tpb(vpx_prob probs[/*num_syms*/],
+               vpx_tree_index tree[/*2*num_syms*/],
+               sym_bools bit_len[/*num_syms*/],
+               const rans_sym sym_tab[/*num_syms*/], int num_syms) {
+  tab2tree(sym_tab, num_syms, probs);
+  build_tree(tree, num_syms);
+  make_tree_bits_tab(bit_len, num_syms);
+}
+
+bool check_vpxtree(const std::vector<int> &sym_vec, const rans_sym *sym_tab,
+                   uint8_t *buf) {
+  vpx_writer w;
+  vpx_reader r;
+  vpx_start_encode(&w, buf);
+
+  vpx_prob probs[kDistinctSyms];
+  vpx_tree_index tree[2 * kDistinctSyms];
+  sym_bools bit_len[kDistinctSyms];
+  build_tpb(probs, tree, bit_len, sym_tab, kDistinctSyms);
+
+  std::clock_t start = std::clock();
+  for (std::vector<int>::const_iterator it = sym_vec.begin();
+       it != sym_vec.end(); ++it) {
+    vp10_write_tree(&w, tree, probs, bit_len[*it].bits, bit_len[*it].len, 0);
+  }
+  std::clock_t enc_time = std::clock() - start;
+  vpx_stop_encode(&w);
+  vpx_reader_init(&r, buf, w.pos, NULL, NULL);
+  start = std::clock();
+  for (std::vector<int>::const_iterator it = sym_vec.begin();
+       it != sym_vec.end(); ++it) {
+    if (vpx_read_tree(&r, tree, probs) != *it) return false;
+  }
+  std::clock_t dec_time = std::clock() - start;
+  printf("VPXtree size %u enc_time %f dec_time %f\n", w.pos,
+         static_cast<float>(enc_time) / CLOCKS_PER_SEC,
+         static_cast<float>(dec_time) / CLOCKS_PER_SEC);
+  return true;
+}
+
+class Vp10AbsTest : public ::testing::Test {
+ protected:
+  static void SetUpTestCase() { pv_vec_ = abs_encode_build_vals(kNumBools); }
+  virtual void SetUp() { buf_ = new uint8_t[kNumBools / 8]; }
+  virtual void TearDown() { delete[] buf_; }
+  static const int kNumBools = 100000000;
+  static PvVec pv_vec_;
+  uint8_t *buf_;
+};
+PvVec Vp10AbsTest::pv_vec_;
+
+class Vp10AnsTest : public ::testing::Test {
+ protected:
+  static void SetUpTestCase() {
+    sym_vec_ = ans_encode_build_vals(rans_sym_tab, kNumSyms);
+  }
+  virtual void SetUp() { buf_ = new uint8_t[kNumSyms / 2]; }
+  virtual void TearDown() { delete[] buf_; }
+  static const int kNumSyms = 25000000;
+  static std::vector<int> sym_vec_;
+  uint8_t *buf_;
+};
+std::vector<int> Vp10AnsTest::sym_vec_;
+
+TEST_F(Vp10AbsTest, Vpxbool) { EXPECT_TRUE(check_vpxbool(pv_vec_, buf_)); }
+TEST_F(Vp10AbsTest, Rabs) { EXPECT_TRUE(check_rabs(pv_vec_, buf_)); }
+TEST_F(Vp10AbsTest, RabsAsc) { EXPECT_TRUE(check_rabs_asc(pv_vec_, buf_)); }
+TEST_F(Vp10AbsTest, Uabs) { EXPECT_TRUE(check_uabs(pv_vec_, buf_)); }
+
+TEST_F(Vp10AnsTest, Rans) {
+  EXPECT_TRUE(check_rans(sym_vec_, rans_sym_tab, buf_));
+}
+TEST_F(Vp10AnsTest, Vpxtree) {
+  EXPECT_TRUE(check_vpxtree(sym_vec_, rans_sym_tab, buf_));
+}
+}  // namespace

vp10/common/ans.h

+/*
+ *  Copyright (c) 2015 The WebM project authors. All Rights Reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license
+ *  that can be found in the LICENSE file in the root of the source
+ *  tree. An additional intellectual property rights grant can be found
+ *  in the file PATENTS.  All contributing project authors may
+ *  be found in the AUTHORS file in the root of the source tree.
+ */
+
+#ifndef VP10_COMMON_ANS_H_
+#define VP10_COMMON_ANS_H_
+// An implementation of Asymmetric Numeral Systems
+// http://arxiv.org/abs/1311.2540v2
+
+#include "./vpx_config.h"
+#include "vpx/vpx_integer.h"
+#include "vpx_ports/mem_ops.h"
+
+#define ANS_DIVIDE_BY_MULTIPLY 1
+#if ANS_DIVIDE_BY_MULTIPLY
+#include "vp10/common/divide.h"
+#define ANS_DIVREM(quotient, remainder, dividend, divisor) \
+  do { \
+    quotient = fastdiv(dividend, divisor); \
+    remainder = dividend - quotient * divisor; \
+  } while (0)
+#define ANS_DIV(dividend, divisor) \
+  fastdiv(dividend, divisor)
+#else
+#define ANS_DIVREM(quotient, remainder, dividend, divisor) \
+  do { \
+    quotient = dividend / divisor; \
+    remainder = dividend % divisor; \
+  } while (0)
+#define ANS_DIV(dividend, divisor) \
+    ((dividend) / (divisor))
+#endif
+
+#ifdef __cplusplus
+extern "C" {
+#endif  // __cplusplus
+
+struct AnsCoder {
+  uint8_t *buf;
+  int buf_offset;
+  uint32_t state;
+};
+
+struct AnsDecoder {
+  const uint8_t *buf;
+  int buf_offset;
+  uint32_t state;
+};
+
+typedef uint8_t AnsP8;
+#define ans_p8_precision 256u
+#define ans_p8_shift 8
+#define l_base (ans_p8_precision * 4)  // l_base % precision must be 0
+#define io_base 256
+// Range I = { l_base, l_base + 1, ..., l_base * io_base - 1 }
+
+static INLINE void ans_write_init(struct AnsCoder *const ans,
+                                  uint8_t *const buf) {
+  ans->buf = buf;
+  ans->buf_offset = 0;
+  ans->state = l_base;
+}
+
+static INLINE int ans_write_end(struct AnsCoder *const ans) {
+  mem_put_le24(ans->buf + ans->buf_offset, ans->state);
+  return ans->buf_offset + 3;
+}
+
+// rABS with descending spread
+// p or p0 takes the place of l_s from the paper
+// ans_p8_precision is m
+static INLINE void rabs_desc_write(struct AnsCoder *ans, int val, AnsP8 p0) {
+  const AnsP8 p = ans_p8_precision - p0;
+  const unsigned l_s = val ? p : p0;
+  unsigned quot, rem;
+  if (ans->state >= l_base / ans_p8_precision * io_base * l_s) {
+    ans->buf[ans->buf_offset++] = ans->state % io_base;
+    ans->state /= io_base;
+  }
+  ANS_DIVREM(quot, rem, ans->state, l_s);
+  ans->state = quot * ans_p8_precision + rem + (val ? 0 : p);
+}
+
+#define ANS_IMPL1 0
+#define UNPREDICTABLE(x) x
+static INLINE int rabs_desc_read(struct AnsDecoder *ans, AnsP8 p0) {
+  int val;
+#if ANS_IMPL1
+  unsigned l_s;
+#else
+  unsigned quot, rem, x, xn;
+#endif
+  const AnsP8 p = ans_p8_precision - p0;
+  if (ans->state < l_base) {
+    ans->state = ans->state * io_base + ans->buf[--ans->buf_offset];
+  }
+#if ANS_IMPL1
+  val = ans->state % ans_p8_precision < p;
+  l_s = val ? p : p0;
+  ans->state = (ans->state / ans_p8_precision) * l_s +
+               ans->state % ans_p8_precision - (!val * p);
+#else
+  x = ans->state;
+  quot = x / ans_p8_precision;
+  rem = x % ans_p8_precision;
+  xn = quot * p;
+  val = rem < p;
+  if (UNPREDICTABLE(val)) {
+    ans->state = xn + rem;
+  } else {
+    // ans->state = quot * p0 + rem - p;
+    ans->state = x - xn - p;
+  }
+#endif
+  return val;
+}
+
+// rABS with ascending spread
+// p or p0 takes the place of l_s from the paper
+// ans_p8_precision is m
+static INLINE void rabs_asc_write(struct AnsCoder *ans, int val, AnsP8 p0) {
+  const AnsP8 p = ans_p8_precision - p0;
+  const unsigned l_s = val ? p : p0;
+  unsigned quot, rem;
+  if (ans->state >= l_base / ans_p8_precision * io_base * l_s) {
+    ans->buf[ans->buf_offset++] = ans->state % io_base;
+    ans->state /= io_base;
+  }
+  ANS_DIVREM(quot, rem, ans->state, l_s);
+  ans->state = quot * ans_p8_precision + rem + (val ? p0 : 0);
+}
+
+static INLINE int rabs_asc_read(struct AnsDecoder *ans, AnsP8 p0) {
+  int val;
+#if ANS_IMPL1
+  unsigned l_s;
+#else
+  unsigned quot, rem, x, xn;
+#endif
+  const AnsP8 p = ans_p8_precision - p0;
+  if (ans->state < l_base) {
+    ans->state = ans->state * io_base + ans->buf[--ans->buf_offset];
+  }
+#if ANS_IMPL1
+  val = ans->state % ans_p8_precision < p;
+  l_s = val ? p : p0;
+  ans->state = (ans->state / ans_p8_precision) * l_s +
+               ans->state % ans_p8_precision - (!val * p);
+#else
+  x = ans->state;
+  quot = x / ans_p8_precision;
+  rem = x % ans_p8_precision;
+  xn = quot * p;
+  val = rem >= p0;
+  if (UNPREDICTABLE(val)) {
+    ans->state = xn + rem - p0;
+  } else {
+    // ans->state = quot * p0 + rem - p0;
+    ans->state = x - xn;
+  }
+#endif
+  return val;
+}
+
+#define rabs_read rabs_desc_read
+#define rabs_write rabs_desc_write
+
+// uABS with normalization
+static INLINE void uabs_write(struct AnsCoder *ans, int val, AnsP8 p0) {
+  AnsP8 p = ans_p8_precision - p0;
+  const unsigned l_s = val ? p : p0;
+  if (ans->state >= l_base / ans_p8_precision * io_base * l_s) {
+    ans->buf[ans->buf_offset++] = ans->state % io_base;
+    ans->state /= io_base;
+  }
+  if (!val)
+    ans->state = ANS_DIV(ans->state * ans_p8_precision, p0);
+  else
+    ans->state = ANS_DIV((ans->state + 1) * ans_p8_precision + p - 1, p) - 1;
+}
+
+static INLINE int uabs_read(struct AnsDecoder *ans, AnsP8 p0) {
+  AnsP8 p = ans_p8_precision - p0;
+  int s;
+  // unsigned int xp1;
+  unsigned xp, sp;
+  unsigned state = ans->state;
+  if (state < l_base && ans->buf_offset > 0) {
+    state = state * io_base + ans->buf[--ans->buf_offset];
+  }
+  sp = state * p;
+  // xp1 = (sp + p) / ans_p8_precision;
+  xp = sp / ans_p8_precision;
+  // s = xp1 - xp;
+  s = (sp & 0xFF) >= p0;
+  if (UNPREDICTABLE(s))
+    ans->state = xp;
+  else
+    ans->state = state - xp;
+  return s;
+}
+
+static INLINE int uabs_read_bit(struct AnsDecoder *ans) {
+  int s;
+  unsigned state = ans->state;
+  if (state < l_base && ans->buf_offset > 0) {
+    state = state * io_base + ans->buf[--ans->buf_offset];
+  }
+  s = (int)(state & 1);
+  ans->state = state >> 1;
+  return s;
+}
+
+struct rans_sym {
+  AnsP8 prob;
+  AnsP8 cum_prob;  // not-inclusive
+};
+
+struct rans_dec_sym {
+  uint8_t val;
+  AnsP8 prob;
+  AnsP8 cum_prob;  // not-inclusive
+};
+
+typedef struct rans_dec_sym rans_dec_lut[ans_p8_precision];
+
+static INLINE void rans_build_dec_tab(const AnsP8 token_probs[],
+                                      rans_dec_lut dec_tab) {
+  int val = 0;
+  int cum_prob = 0;
+  int sym_end = token_probs[0];
+  int i;
+  for (i = 0; i < 256; ++i) {
+    if (i == sym_end) {
+      ++val;
+      cum_prob = sym_end;
+      sym_end += token_probs[val];
+    }
+    dec_tab[i].val = val;
+    dec_tab[i].prob = token_probs[val];
+    dec_tab[i].cum_prob = cum_prob;
+  }
+}
+
+// rANS with normalization