1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
// Copyright 2018 Johannes Köster.
// Licensed under the MIT license (http://opensource.org/licenses/MIT)
// This file may not be copied, modified, or distributed
// except according to those terms.
#[inline]
pub fn interpolate(a: f64, b: f64, fraction: f64) -> f64 {
a * (1.0 - fraction) + b * fraction
}
/// Fast lookup table for arbitrary floating point functions.
/// This can be used to e.g., provide fast lookups of distribution values.
/// Input values are sampled with a given precision and results are stored in a vector.
/// During lookup, infimum and supremum of a given value are calculated and the result is
/// interpolated.
pub struct InterpolationTable<F: Fn(f64) -> f64> {
inner: Vec<f64>,
func: F,
offset: usize,
min_x: f64,
max_x: f64,
shift: f64,
}
impl<F: Fn(f64) -> f64> InterpolationTable<F> {
/// Create a new `InterpolationTable`.
///
/// # Arguments
///
/// * `min_x` - minimum sample value
/// * `max_x` - maximum sample value
/// * `frac_digits` - number of fraction digits to store in sample
/// * `func` - Function to emulate.
///
/// If given value is outside of min_x and max_x, the lookup falls back to applying the
/// function itself.
/// The table size grows with the number of fraction digits.
pub fn new(min_x: f64, max_x: f64, frac_digits: i32, func: F) -> Self {
let shift = 10.0_f64.powi(frac_digits);
let offset = (min_x * shift) as usize;
let mut table = InterpolationTable {
inner: Vec::new(),
func,
min_x,
max_x,
shift,
offset,
};
let mut i = table.index(min_x);
while i < table.index(max_x) {
let x = i as f64 / shift;
table.inner.push((table.func)(x));
i += 1;
}
table
}
/// Return lower bound index for given f64.
#[inline]
fn index(&self, x: f64) -> usize {
(x * self.shift) as usize - self.offset
}
/// Lookup given value in table, and interpolate the result between the sampled values if
/// necessary. This provides an approximation that is better the more fraction digits are
/// used to generate this table.
pub fn get(&self, x: f64) -> f64 {
if x < self.min_x || x >= self.max_x {
(self.func)(x)
} else {
let i = self.index(x);
// interpolate
let fraction = (x * self.shift - i as f64) / self.shift;
interpolate(self.inner[i], self.inner[i + 1], fraction)
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_interpolation_table() {
let table = InterpolationTable::new(0.0, 10.0, 5, |x| x.ln_1p());
for &x in &[0.02, 0.04, 0.45678686, 0.23875, 1.45345e-6] {
assert_relative_eq!(table.get(x), x.ln_1p(), epsilon = 0.00001);
}
}
}