unthreaded/src/bisection_tree/aggregator.rs

/*!
Aggregation / summarisation of information in branches of bisection trees.
*/

use crate::types::*;
use crate::sets::Set;

/// Trait for aggregating information about a branch of a [bisection tree][super::BT].
///
/// Currently [`Bounds`] is the only provided aggregator.
/// It keeps track of upper and lower bounds of a function representeed by the `BT` by
/// summing [`Bounds`] produced by [`LocalAnalysis`][super::support::LocalAnalysis] of the
/// [`Support`][super::support::Support]s of the data stored in the tree.
/// For the `Bounds` aggregator:
///  * [`Self::aggregate`] sums input bounds to the current bound. This provides a conservative
///    estimate of the upper and lower bounds of a sum of functions.
///  * [`Self::summarise`] takes the maximum of the input bounds. This calculates the bounds
///    of a function on a greater domain from bounds on subdomains
///    (in practise [`Cube`][crate::sets::Cube]s).
///
pub trait Aggregator : Clone + std::fmt::Debug {
    /// Aggregate a new data to current state.
    fn aggregate<I>(&mut self, aggregates : I)
    where I : Iterator<Item=Self>;

    /// Summarise several other aggregators, resetting current state.
    fn summarise<'a, I>(&'a mut self, aggregates : I)
    where I : Iterator<Item=&'a Self>;

    /// Create a new “empty” aggregate data.
    fn new() -> Self;
}

/// An [`Aggregator`] that doesn't aggregate anything.
#[derive(Clone,Debug)]
pub struct NullAggregator;

impl Aggregator for NullAggregator {
    // TODO: these should probabably also take a Cube as argument to
    // allow integrating etc.
    fn aggregate<I>(&mut self, _aggregates : I)
    where I : Iterator<Item=Self> {}

    fn summarise<'a, I>(&'a mut self, _aggregates : I)
    where I : Iterator<Item=&'a Self> {}

    fn new() -> Self { NullAggregator }
}

/// Upper and lower bounds on an `F`-valued function.
#[derive(Copy,Clone,Debug)]
pub struct Bounds<F>(
    /// Lower bound
    pub F,
    /// Upper bound
    pub F
);

impl<F : Copy> Bounds<F> {
    /// Returns the lower bound
    #[inline]
    pub fn lower(&self) -> F { self.0 }

    /// Returns the upper bound
    #[inline]
    pub fn upper(&self) -> F { self.1 }
}

impl<F : Float> Bounds<F> {
    /// Returns a uniform bound.
    ///
    /// This is maximum over the absolute values of the upper and lower bound.
    #[inline]
    pub fn uniform(&self) -> F {
        let &Bounds(lower, upper) = self;
        lower.abs().max(upper.abs())
    }
}

impl<'a, F : Float> std::ops::Add<Self> for Bounds<F> {
    type Output = Self;
    #[inline]
    fn add(self, Bounds(l2, u2) : Self) -> Self::Output {
        let Bounds(l1, u1) = self;
        debug_assert!(l1 <= u1 && l2 <= u2);
        Bounds(l1 + l2, u1 + u2)
    }
}

impl<'a, F : Float> std::ops::Mul<Self> for Bounds<F> {
    type Output = Self;
    #[inline]
    fn mul(self, Bounds(l2, u2) : Self) -> Self::Output {
        let Bounds(l1, u1) = self;
        debug_assert!(l1 <= u1 && l2 <= u2);
        let a = l1 * l2;
        let b = u1 * u2;
        // The order may flip when negative numbers are involved, so need min/max
        Bounds(a.min(b), a.max(b))
    }
}

impl<F : Float> std::iter::Product for Bounds<F> {
    #[inline]
    fn product<I>(mut iter: I) -> Self
    where I: Iterator<Item = Self> {
        match iter.next() {
            None => Bounds(F::ZERO, F::ZERO),
            Some(init) => iter.fold(init, |a, b| a*b)
        }
    }
}

impl<F : Float> Set<F> for Bounds<F> {
    fn contains(&self, item : &F) -> bool {
        let &Bounds(l, u) = self;
        debug_assert!(l <= u);
        l <= *item && *item <= u
    }
}

impl<F : Float> Bounds<F> {
    /// Calculate a common bound (glb, lub) for two bounds.
    #[inline]
    pub fn common(&self, &Bounds(l2, u2) : &Self) -> Self {
        let &Bounds(l1, u1) = self;
        debug_assert!(l1 <= u1 && l2 <= u2);
        Bounds(l1.min(l2), u1.max(u2))
    }

    /// Indicates whether `Self` is a superset of the argument bound.
    #[inline]
    pub fn superset(&self, &Bounds(l2, u2) : &Self) -> bool {
        let &Bounds(l1, u1) = self;
        debug_assert!(l1 <= u1 && l2 <= u2);
        l1 <= l2 && u2 <= u1
    }

    /// Returns the greatest bound contained by both argument bounds, if one exists.
    #[inline]
    pub fn glb(&self, &Bounds(l2, u2) : &Self) -> Option<Self> {
        let &Bounds(l1, u1) = self;
        debug_assert!(l1 <= u1 && l2 <= u2);
        let l = l1.max(l2);
        let u = u1.min(u2);
        if l < u {
            Some(Bounds(l, u))
        } else {
            None
        }
    }
}

impl<F : Float> Aggregator for Bounds<F> {
    #[inline]
    fn aggregate<I>(&mut self, aggregates : I)
    where I : Iterator<Item=Self> {
        *self = aggregates.fold(*self, |a, b| a + b);
    }

    #[inline]
    fn summarise<'a, I>(&'a mut self, mut aggregates : I)
    where I : Iterator<Item=&'a Self> {
        *self = match aggregates.next() {
            None => Bounds(F::ZERO, F::ZERO), // No parts in this cube; the function is zero
            Some(&bounds) => {
                aggregates.fold(bounds, |a, b| a.common(b))
            }
        }
    }

    #[inline]
    fn new() -> Self {
        Bounds(F::ZERO, F::ZERO)
    }
}