1 // Copyright 2015 Emilie Gillet.
2 //
3 // Author: Emilie Gillet (emilie.o.gillet@gmail.com)
4 //
5 // Permission is hereby granted, free of charge, to any person obtaining a copy
6 // of this software and associated documentation files (the "Software"), to deal
7 // in the Software without restriction, including without limitation the rights
8 // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9 // copies of the Software, and to permit persons to whom the Software is
10 // furnished to do so, subject to the following conditions:
11 //
12 // The above copyright notice and this permission notice shall be included in
13 // all copies or substantial portions of the Software.
14 //
15 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18 // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20 // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
21 // THE SOFTWARE.
22 //
23 // See http://creativecommons.org/licenses/MIT/ for more information.
24 //
25 // -----------------------------------------------------------------------------
26 //
27 // Generates samples from various kinds of random distributions.
28
29 #ifndef MARBLES_RANDOM_DISTRIBUTIONS_H_
30 #define MARBLES_RANDOM_DISTRIBUTIONS_H_
31
32 #include "stmlib/stmlib.h"
33
34 #include <algorithm>
35
36 #include "stmlib/dsp/dsp.h"
37
38 #include "marbles/resources.h"
39
40 namespace marbles {
41
42 const size_t kNumBiasValues = 5;
43 const size_t kNumRangeValues = 9;
44 const float kIcdfTableSize = 128.0f;
45
46 // Generates samples from beta distribution, from uniformly distributed samples.
47 // For higher throughput, uses pre-computed tables of inverse cdfs.
BetaDistributionSample(float uniform,float spread,float bias)48 inline float BetaDistributionSample(float uniform, float spread, float bias) {
49 // Tables are pre-computed only for bias <= 0.5. For values above 0.5,
50 // symmetry is used.
51 bool flip_result = bias > 0.5f;
52 if (flip_result) {
53 uniform = 1.0f - uniform;
54 bias = 1.0f - bias;
55 }
56
57 bias *= (static_cast<float>(kNumBiasValues) - 1.0f) * 2.0f;
58 spread *= (static_cast<float>(kNumRangeValues) - 1.0f);
59
60 MAKE_INTEGRAL_FRACTIONAL(bias);
61 MAKE_INTEGRAL_FRACTIONAL(spread);
62
63 size_t cell = bias_integral * (kNumRangeValues + 1) + spread_integral;
64
65 // Lower 5% and 95% percentiles use a different table with higher resolution.
66 size_t offset = 0;
67 if (uniform <= 0.05f) {
68 offset = kIcdfTableSize + 1;
69 uniform *= 20.0f;
70 } else if (uniform >= 0.95f) {
71 offset = 2 * (kIcdfTableSize + 1);
72 uniform = (uniform - 0.95f) * 20.0f;
73 }
74
75 float x1y1 = stmlib::Interpolate(
76 distributions_table[cell] + offset,
77 uniform,
78 kIcdfTableSize);
79 float x2y1 = stmlib::Interpolate(
80 distributions_table[cell + 1] + offset,
81 uniform,
82 kIcdfTableSize);
83 float x1y2 = stmlib::Interpolate(
84 distributions_table[cell + kNumRangeValues + 1] + offset,
85 uniform,
86 kIcdfTableSize);
87 float x2y2 = stmlib::Interpolate(
88 distributions_table[cell + kNumRangeValues + 2] + offset,
89 uniform,
90 kIcdfTableSize);
91
92 float y1 = x1y1 + (x2y1 - x1y1) * spread_fractional;
93 float y2 = x1y2 + (x2y2 - x1y2) * spread_fractional;
94 float y = y1 + (y2 - y1) * bias_fractional;
95
96 if (flip_result) {
97 y = 1.0f - y;
98 }
99 return y;
100 }
101
102 // Pre-computed beta(3, 3) with a fatter tail.
FastBetaDistributionSample(float uniform)103 inline float FastBetaDistributionSample(float uniform) {
104 return stmlib::Interpolate(dist_icdf_4_3, uniform, kIcdfTableSize);
105 }
106
107 // Draws samples from a discrete distribution. Used for the quantizer.
108 // Example:
109 // * 1 with probability 0.2
110 // * 20 with probability 0.7
111 // * 666 with probability 0.1
112 //
113 // DiscreteDistribution d;
114 // d.Init();
115 // d.AddToken(1, 0.2);
116 // d.AddToken(20, 0.7);
117 // d.AddToken(666, 0.1);
118 // d.NoMoreTokens();
119 // Result r = d.Sample(u);
120 // cout << r.token_id;
121 //
122 // Weights do not have to add to 1.0f - the class handles normalization.
123 //
124 template<size_t size>
125 class DiscreteDistribution {
126 public:
DiscreteDistribution()127 DiscreteDistribution() { }
~DiscreteDistribution()128 ~DiscreteDistribution() { }
129
Init()130 void Init() {
131 sum_ = 0.0f;
132 num_tokens_ = 1;
133
134 cdf_[0] = 0.0f;
135 token_ids_[0] = 0;
136 }
137
AddToken(int token_id,float weight)138 void AddToken(int token_id, float weight) {
139 if (weight <= 0.0f) {
140 return;
141 }
142 sum_ += weight;
143 token_ids_[num_tokens_] = token_id;
144 cdf_[num_tokens_] = sum_;
145 ++num_tokens_;
146 }
147
NoMoreTokens()148 void NoMoreTokens() {
149 token_ids_[num_tokens_] = token_ids_[num_tokens_ - 1];
150 cdf_[num_tokens_] = sum_ + 1.0f;
151 }
152
153 struct Result {
154 int token_id;
155 float fraction;
156 float start;
157 float width;
158 };
159
Sample(float u)160 inline Result Sample(float u) const {
161 Result r;
162 u *= sum_;
163 int n = std::upper_bound(&cdf_[1], &cdf_[num_tokens_ + 1], u) - &cdf_[0];
164 float norm = 1.0f / sum_;
165 r.token_id = token_ids_[n];
166 r.width = (cdf_[n] - cdf_[n - 1]) * norm;
167 r.start = (cdf_[n - 1]) * norm;
168 r.fraction = (u - cdf_[n - 1]) / (cdf_[n] - cdf_[n - 1]);
169 return r;
170 }
171
172 float sum_;
173 float cdf_[size + 2];
174 int token_ids_[size + 2];
175 int num_tokens_;
176
177 DISALLOW_COPY_AND_ASSIGN(DiscreteDistribution);
178 };
179
180 } // namespace marbles
181
182 #endif // MARBLES_RANDOM_DISTRIBUTIONS_H_
183