1 // Copyright (c) 2017, Steinar H. Gunderson
2 // All rights reserved.
4 // Redistribution and use in source and binary forms, with or without
5 // modification, are permitted.
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.
19 #include "renormalize.h"
24 #include <unordered_map>
35 using std::unique_ptr;
36 using std::unordered_map;
40 struct OptimalChoice {
41 double cost; // In bits.
48 bool operator== (const CacheKey &other) const
50 return num_syms == other.num_syms && available_slots == other.available_slots;
54 size_t operator() (const CacheKey &key) const
56 return hash<int64_t>()((uint64_t(key.available_slots) << 32) | key.num_syms);
59 using CacheMap = unordered_map<CacheKey, OptimalChoice, HashCacheKey>;
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.
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)
70 // Encoding zero symbols needs zero bits.
73 if (num_syms > available_slots) {
74 // Every (non-zero-frequency) symbol needs at least one slot.
78 return cum_freqs[1] * log2cache[available_slots];
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;
88 // Minimize the number of total bits spent as a function of how many slots
89 // we assign to this symbol.
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.
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];
100 double guess = lrint(available_slots * double(freq) / cum_freqs[num_syms]);
102 int x1 = max<int>(floor(guess), 1);
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);
109 int direction; // -1 or +1.
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.
117 } else if (cost1 < cost2) {
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.
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.
144 insert_result.first->second.cost = best_cost;
145 insert_result.first->second.chosen_freq = best_choice;
151 void OptimalRenormalize(uint32_t *cum_freqs, uint32_t num_syms, uint32_t target_total)
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]);
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]) {
164 mapping[new_num_syms] = i;
165 remapped_cum_freqs[new_num_syms + 1] = cum_freqs[i + 1];
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));
177 FindOptimalCost(remapped_cum_freqs.get(), new_num_syms, target_total, log2cache.get(), &cache);
179 for (uint32_t i = 0; i <= num_syms; ++i) {
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; ) { // :-)
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;
192 CacheKey cache_key{symbol_idx + 1, available_slots};
193 assert(cache.count(cache_key));
194 freq = cache[cache_key].chosen_freq;
196 cum_freqs[mapping[symbol_idx]] = freq;
197 assert(available_slots >= 0 && unsigned(available_slots) >= freq);
198 available_slots -= freq;
201 // Convert the frequencies back to cumulative frequencies.
203 for (uint32_t i = 0; i <= num_syms; ++i) {
204 uint32_t freq = cum_freqs[i];
205 cum_freqs[i] = total;