Smoothing Functions For Actors

Cargo.toml

@@ -33,6 +33,8 @@ members = [
     "utils/test_utils",
     "utils/commcid",
     "utils/json_utils",
+    "utils/smoothing",
+    "utils/math",
     "types",
     "key_management",
 ]

utils/math/Cargo.toml

@@ -0,0 +1,8 @@
+[package]
+name = "forest_math"
+version = "0.1.0"
+authors = ["ChainSafe Systems <info@chainsafe.io>"]
+edition = "2018"
+
+[dependencies]
+num-bigint = { path = "../../utils/bigint", package = "forest_bigint" }

utils/math/src/lib.rs

@@ -0,0 +1,22 @@
+// Copyright 2020 ChainSafe Systems
+// SPDX-License-Identifier: Apache-2.0, MIT
+
+use num_bigint::{BigInt, ParseBigIntError};
+
+pub const PRECISION: u64 = 128;
+
+/// polyval evaluates a polynomial given by coefficients `p` in Q.128 format
+/// at point `x` in Q.128 format. Output is in Q.128.
+/// Coefficients should be ordered from the highest order coefficient to the lowest.
+pub fn poly_val(poly: &[BigInt], x: &BigInt) -> BigInt {
+    let mut res = BigInt::default();
+
+    for coeff in poly {
+        res = ((res * x) >> PRECISION) + coeff;
+    }
+    res
+}
+
+pub fn parse(coefs: &[&str]) -> Result<Vec<BigInt>, ParseBigIntError> {
+    coefs.iter().map(|c| c.parse()).collect()
+}

utils/smoothing/Cargo.toml

@@ -0,0 +1,16 @@
+[package]
+name = "forest_smoothing"
+version = "0.1.0"
+authors = ["ChainSafe Systems <info@chainsafe.io>"]
+edition = "2018"
+
+[dependencies]
+num-bigint = { path = "../../utils/bigint", package = "forest_bigint" }
+clock = { path = "../../node/clock" }
+forest_math = {path = "../math"}
+lazy_static = "1.4"
+num-traits = "0.2"
+actor = {path = "../../vm/actor", package = "actor"}
+fil_types = { path = "../../types" }
+encoding = { package = "forest_encoding", path = "../../encoding" }
+serde = { version = "1.0", features = ["derive"] }

utils/smoothing/src/alpha_beta_filter.rs

@@ -0,0 +1,69 @@
+// Copyright 2020 ChainSafe Systems
+// SPDX-License-Identifier: Apache-2.0, MIT
+
+use clock::ChainEpoch;
+use encoding::Cbor;
+use forest_math::PRECISION;
+use num_bigint::{bigint_ser, BigInt};
+use serde::{Deserialize, Serialize};
+
+#[derive(Default, Serialize, Deserialize)]
+pub struct FilterEstimate {
+    #[serde(with = "bigint_ser")]
+    pub pos: BigInt,
+    #[serde(with = "bigint_ser")]
+    pub velo: BigInt,
+}
+
+impl FilterEstimate {
+    pub fn new(pos: BigInt, velo: BigInt) -> Self {
+        FilterEstimate {
+            pos: pos << PRECISION,
+            velo: velo << PRECISION,
+        }
+    }
+
+    pub fn estimate(&self) -> BigInt {
+        &self.pos >> PRECISION
+    }
+
+    pub fn extrapolate(&self, delta: ChainEpoch) -> BigInt {
+        let delta_t = BigInt::from(delta) << PRECISION;
+        let pos = &self.pos << PRECISION;
+        (&self.velo * delta_t) + pos
+    }
+}
+
+impl Cbor for FilterEstimate {}
+
+#[derive(Default)]
+pub struct AlphaBetaFilter {
+    alpha: BigInt,
+    beta: BigInt,
+    prev_est: FilterEstimate,
+}
+
+impl AlphaBetaFilter {
+    pub fn load_filter(prev_est: FilterEstimate, alpha: BigInt, beta: BigInt) -> Self {
+        AlphaBetaFilter {
+            alpha,
+            beta,
+            prev_est,
+        }
+    }
+
+    pub fn next_estimate(&self, obs: BigInt, epoch_delta: ChainEpoch) -> FilterEstimate {
+        let delta_t = BigInt::from(epoch_delta) << PRECISION;
+        let delta_x = (&delta_t * &self.prev_est.velo) >> PRECISION;
+        let mut pos = delta_x + &self.prev_est.pos;
+
+        let obs = obs << PRECISION;
+        let residual = obs - &pos;
+        let revision_x = (&self.alpha * &residual) >> PRECISION;
+        pos += &revision_x;
+
+        let revision_v = (residual * &self.beta) / delta_t;
+        let velo = revision_v + &self.prev_est.velo;
+        FilterEstimate { pos, velo }
+    }
+}

utils/smoothing/src/constants.rs

@@ -0,0 +1,30 @@
+// Copyright 2020 ChainSafe Systems
+// SPDX-License-Identifier: Apache-2.0, MIT
+
+use forest_math::parse;
+use num_bigint::BigInt;
+
+lazy_static! {
+    pub static ref NUM: Vec<BigInt> = parse(&[
+        "261417938209272870992496419296200268025",
+        "7266615505142943436908456158054846846897",
+        "32458783941900493142649393804518050491988",
+        "17078670566130897220338060387082146864806",
+        "-35150353308172866634071793531642638290419",
+        "-20351202052858059355702509232125230498980",
+        "-1563932590352680681114104005183375350999",
+    ])
+    .unwrap();
+    pub static ref DENOM: Vec<BigInt> = parse(&[
+        "49928077726659937662124949977867279384",
+        "2508163877009111928787629628566491583994",
+        "21757751789594546643737445330202599887121",
+        "53400635271583923415775576342898617051826",
+        "41248834748603606604000911015235164348839",
+        "9015227820322455780436733526367238305537",
+        "340282366920938463463374607431768211456",
+    ])
+    .unwrap();
+    pub static ref LN_2: BigInt = "235865763225513294137944142764154484399".parse().unwrap();
+    pub static ref EPSILON: BigInt = "302231454903657293676544".parse().unwrap();
+}

utils/smoothing/src/lib.rs

@@ -0,0 +1,81 @@
+// Copyright 2020 ChainSafe Systems
+// SPDX-License-Identifier: Apache-2.0, MIT
+
+#[macro_use]
+extern crate lazy_static;
+
+use num_bigint::BigInt;
+
+mod alpha_beta_filter;
+mod constants;
+pub use alpha_beta_filter::*;
+use clock::ChainEpoch;
+use constants::*;
+use forest_math::{poly_val, PRECISION};
+use num_traits::sign::Signed;
+
+fn get_bit_len(z: &BigInt) -> usize {
+    z.abs().to_radix_le(2).1.len()
+}
+
+pub fn extrapolated_cum_sum_of_ratio(
+    delta: ChainEpoch,
+    relative_start: ChainEpoch,
+    est_num: &FilterEstimate,
+    est_denom: &FilterEstimate,
+) -> BigInt {
+    let delta_t = BigInt::from(delta) << PRECISION;
+    let t0 = BigInt::from(relative_start) << PRECISION;
+
+    let pos_1 = &est_num.pos;
+    let pos_2 = &est_denom.pos;
+    let velo_1 = &est_num.velo;
+    let velo_2 = &est_denom.velo;
+
+    let squared_velo_2 = (velo_2 * velo_2) >> PRECISION;
+
+    if squared_velo_2 >= EPSILON.clone() {
+        let mut x2a = ((velo_2 * t0) >> PRECISION) + pos_2;
+        let mut x2b = ((velo_2 * &delta_t) >> PRECISION) + &x2a;
+        x2a = ln(&x2a);
+        x2b = ln(&x2b);
+
+        let m1 = ((&x2b - &x2a) * pos_1 * velo_2) >> PRECISION;
+
+        let m2_l = (&x2a - &x2b) * pos_2;
+        let m2_r = velo_2 * &delta_t;
+        let m2 = ((m2_l + m2_r) * velo_1) >> PRECISION;
+
+        return (m2 + m1) / squared_velo_2;
+    }
+
+    let half_delta = &delta_t >> 1;
+    let mut x1m = velo_1 * (t0 + half_delta);
+    x1m = (x1m >> PRECISION) + pos_1;
+
+    (x1m * delta_t) / pos_2
+}
+
+pub fn ln(z: &BigInt) -> BigInt {
+    let k: i64 = get_bit_len(z) as i64 - 1 - PRECISION as i64;
+
+    let x: BigInt = if k > 0 { z >> k } else { z << k.abs() };
+
+    BigInt::from(k) * LN_2.clone() + ln_between_one_and_two(x)
+}
+
+fn ln_between_one_and_two(x: BigInt) -> BigInt {
+    let num = poly_val(&NUM, &x) << PRECISION;
+    let denom = poly_val(&DENOM, &x);
+    num / denom
+}
+
+// Returns an estimate with position val and velocity 0
+pub fn testing_constant_estimate(val: BigInt) -> FilterEstimate {
+    FilterEstimate::new(val, BigInt::from(0u8))
+}
+
+// Returns and estimate with postion x and velocity v
+pub fn testing_estimate(x: BigInt, v: BigInt) -> FilterEstimate {
+    FilterEstimate::new(x, v)
+}