bio/data_structures/
interpolation_table.rs

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

bio/data_structures/
interpolation_table.rs
