git.sesse.net Git - narabu/blob - ryg_rans/renormalize.cpp

   1 // Copyright (c) 2017, Steinar H. Gunderson
   2 // All rights reserved.
   3 //
   4 // Redistribution and use in source and binary forms, with or without
   5 // modification, are permitted.
   6 //
   7 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   8 // “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   9 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
  10 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
  11 // HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
  12 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
  13 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
  14 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
  15 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
  16 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  17 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  18
  19 #include "renormalize.h"
  20
  21 #include <assert.h>
  22 #include <math.h>
  23
  24 #include <unordered_map>
  25 #include <map>
  26 #include <memory>
  27 #include <utility>
  28
  29 using std::equal_to;
  30 using std::hash;
  31 using std::max;
  32 using std::min;
  33 using std::make_pair;
  34 using std::pair;
  35 using std::unique_ptr;
  36 using std::unordered_map;
  37
  38 namespace {
  39
  40 struct OptimalChoice {
  41         double cost;  // In bits.
  42         uint32_t chosen_freq;
  43 };
  44 struct CacheKey {
  45         int num_syms;
  46         int available_slots;
  47
  48         bool operator== (const CacheKey &other) const
  49         {
  50                 return num_syms == other.num_syms && available_slots == other.available_slots;
  51         }
  52 };
  53 struct HashCacheKey {
  54         size_t operator() (const CacheKey &key) const
  55         {
  56                 return hash<int64_t>()((uint64_t(key.available_slots) << 32) | key.num_syms);
  57         }
  58 };
  59 using CacheMap = unordered_map<CacheKey, OptimalChoice, HashCacheKey>;
  60
  61 // Find, recursively, the optimal cost of encoding the symbols [0, num_syms),
  62 // assuming an optimal distribution of those symbols to "available_slots".
  63 // The cache is used for memoization, and also to remember the best choice.
  64 // No frequency can be zero.
  65 //
  66 // Returns HUGE_VAL if there's no legal mapping.
  67 double FindOptimalCost(uint32_t *cum_freqs, int num_syms, int available_slots, const double *log2cache, CacheMap *cache)
  68 {
  69         if (num_syms == 0) {
  70                 // Encoding zero symbols needs zero bits.
  71                 return 0.0;
  72         }
  73         if (num_syms > available_slots) {
  74                 // Every (non-zero-frequency) symbol needs at least one slot.
  75                 return HUGE_VAL;
  76         }
  77         if (num_syms == 1) {
  78                 return cum_freqs[1] * log2cache[available_slots];
  79         }
  80
  81         CacheKey cache_key{num_syms, available_slots};
  82         auto insert_result = cache->insert(make_pair(cache_key, OptimalChoice()));
  83         if (!insert_result.second) {
  84                 // There was already an item in the cache, so return it.
  85                 return insert_result.first->second.cost;
  86         }
  87
  88         // Minimize the number of total bits spent as a function of how many slots
  89         // we assign to this symbol.
  90         //
  91         // The cost function is convex (at least in practice; I suppose also in
  92         // theory because it's the sum of an increasing and a decreasing function?).
  93         // Find a reasonable guess and see in what direction the function is decreasing,
  94         // then iterate until we either hit the end or we start increasing again.
  95         //
  96         // Since the function is a sum of log() terms, it is differentiable, and we
  97         // could in theory use this; however, it doesn't seem to be worth the complexity.
  98         uint32_t freq = cum_freqs[num_syms] - cum_freqs[num_syms - 1];
  99         assert(freq > 0);
 100         double guess = lrint(available_slots * double(freq) / cum_freqs[num_syms]);
 101
 102         int x1 = max<int>(floor(guess), 1);
 103         int x2 = x1 + 1;
 104
 105         double cost1 = freq * log2cache[x1] + FindOptimalCost(cum_freqs, num_syms - 1, available_slots - x1, log2cache, cache);
 106         double cost2 = freq * log2cache[x2] + FindOptimalCost(cum_freqs, num_syms - 1, available_slots - x2, log2cache, cache);
 107
 108         int x;
 109         int direction;  // -1 or +1.
 110         double best_cost;
 111         if (isinf(cost1) && isinf(cost2)) {
 112                 // The cost isn't infinite due to the first term, so we need to go downwards
 113                 // to give the second term more room to breathe.
 114                 x = x1;
 115                 best_cost = cost1;
 116                 direction = -1;
 117         } else if (cost1 < cost2) {
 118                 x = x1;
 119                 best_cost = cost1;
 120                 direction = -1;
 121         } else {
 122                 x = x2;
 123                 best_cost = cost2;
 124                 direction = 1;
 125         }
 126         int best_choice = x;
 127
 128         for ( ;; ) {
 129                 x += direction;
 130                 if (x == 0 || x > available_slots) {
 131                         // We hit the end; we can't assign zero slots to this symbol,
 132                         // and we can't assign more slots than we have. This extreme
 133                         // is the best choice.
 134                         break;
 135                 }
 136                 double cost = freq * log2cache[x] + FindOptimalCost(cum_freqs, num_syms - 1, available_slots - x, log2cache, cache);
 137                 if (cost > best_cost) {
 138                         // The cost started increasing again, so we've found the optimal choice.
 139                         break;
 140                 }
 141                 best_choice = x;
 142                 best_cost = cost;
 143         }
 144         insert_result.first->second.cost = best_cost;
 145         insert_result.first->second.chosen_freq = best_choice;
 146         return best_cost;
 147 }
 148
 149 }  // namespace
 150
 151 void OptimalRenormalize(uint32_t *cum_freqs, uint32_t num_syms, uint32_t target_total)
 152 {
 153         // First remove all symbols that have a zero frequency; they tend to
 154         // complicate the analysis. We'll put them back afterwards.
 155         unique_ptr<uint32_t[]> remapped_cum_freqs(new uint32_t[num_syms + 1]);
 156         unique_ptr<uint32_t[]> mapping(new uint32_t[num_syms + 1]);
 157
 158         uint32_t new_num_syms = 0;
 159         remapped_cum_freqs[0] = 0;
 160         for (uint32_t i = 0; i < num_syms; ++i) {
 161                 if (cum_freqs[i + 1] == cum_freqs[i]) {
 162                         continue;
 163                 }
 164                 mapping[new_num_syms] = i;
 165                 remapped_cum_freqs[new_num_syms + 1] = cum_freqs[i + 1];
 166                 new_num_syms++;
 167         }
 168
 169         // Calculate the cost of encoding a symbol with frequency f/target_total.
 170         // We call log2() quite a lot, so it's best to cache it once at the start.
 171         unique_ptr<double[]> log2cache(new double[target_total + 1]);
 172         for (uint32_t i = 0; i <= target_total; ++i) {
 173                 log2cache[i] = -log2(i * (1.0 / target_total));
 174         }
 175
 176         CacheMap cache;
 177         FindOptimalCost(remapped_cum_freqs.get(), new_num_syms, target_total, log2cache.get(), &cache);
 178
 179         for (uint32_t i = 0; i <= num_syms; ++i) {
 180                 cum_freqs[i] = 0;
 181         }
 182
 183         // Reconstruct the optimal choices from the cache. Note that during this,
 184         // cum_freq contains frequencies, _not_ cumulative frequencies.
 185         int available_slots = target_total;
 186         for (int symbol_idx = new_num_syms; symbol_idx --> 0; ) {  // :-)
 187                 uint32_t freq;
 188                 if (symbol_idx == 0) {
 189                         // Last symbol isn't in the cache, but it's obvious what the answer is.
 190                         freq = available_slots;
 191                 } else {
 192                         CacheKey cache_key{symbol_idx + 1, available_slots};
 193                         assert(cache.count(cache_key));
 194                         freq = cache[cache_key].chosen_freq;
 195                 }
 196                 cum_freqs[mapping[symbol_idx]] = freq;
 197                 assert(available_slots >= 0 && unsigned(available_slots) >= freq);
 198                 available_slots -= freq;
 199         }
 200
 201         // Convert the frequencies back to cumulative frequencies.
 202         uint32_t total = 0;
 203         for (uint32_t i = 0; i <= num_syms; ++i) {
 204                 uint32_t freq = cum_freqs[i];
 205                 cum_freqs[i] = total;
 206                 total += freq;
 207         }
 208 }