feat: add syntactic sugar for signal operations

argus-core/Cargo.toml

@@ -5,6 +5,8 @@ edition = "2021"
 
 [dependencies]
 derive_more = "0.99.17"
+itertools = "0.10.5"
+paste = "1.0.12"
 num-traits = "0.2.15"
 thiserror = "1.0.39"
 

argus-core/src/signals.rs

@@ -8,8 +8,16 @@
 //! 2. [`ConstantSignal<T>`] is a signal that maintains a constant value throughtout
 //!    its domain, and thus, do not require interpolation and extrapolation. Moreover,
 //!    since they are defined over the entire time domain, they cannot be iterated over.
+pub mod bool_ops;
+pub mod cmp_ops;
 pub mod iter;
+pub mod num_ops;
 pub mod traits;
+mod utils;
+
+pub use bool_ops::*;
+pub use cmp_ops::*;
+pub use num_ops::*;
 
 use std::ops::{RangeFull, RangeInclusive};
 use std::time::Duration;

argus-core/src/signals/bool_ops.rs

@@ -0,0 +1,81 @@
+use crate::signals::utils::{apply1, apply2, apply2_const};
+use crate::signals::{ConstantSignal, Signal};
+
+impl core::ops::Not for &Signal<bool> {
+    type Output = Signal<bool>;
+
+    fn not(self) -> Self::Output {
+        apply1(self, |v| !v)
+    }
+}
+
+impl core::ops::BitAnd<Self> for &Signal<bool> {
+    type Output = Signal<bool>;
+
+    fn bitand(self, other: Self) -> Self::Output {
+        apply2(self, other, |lhs, rhs| lhs && rhs)
+    }
+}
+
+impl core::ops::BitAnd<&ConstantSignal<bool>> for &Signal<bool> {
+    type Output = Signal<bool>;
+
+    fn bitand(self, other: &ConstantSignal<bool>) -> Self::Output {
+        apply2_const(self, other, |lhs, rhs| lhs && rhs)
+    }
+}
+
+impl core::ops::BitOr<Self> for &Signal<bool> {
+    type Output = Signal<bool>;
+
+    fn bitor(self, other: Self) -> Self::Output {
+        apply2(self, other, |lhs, rhs| lhs && rhs)
+    }
+}
+
+impl core::ops::BitOr<&ConstantSignal<bool>> for &Signal<bool> {
+    type Output = Signal<bool>;
+
+    fn bitor(self, other: &ConstantSignal<bool>) -> Self::Output {
+        apply2_const(self, other, |lhs, rhs| lhs && rhs)
+    }
+}
+
+impl core::ops::Not for &ConstantSignal<bool> {
+    type Output = ConstantSignal<bool>;
+
+    fn not(self) -> Self::Output {
+        ConstantSignal::<bool>::new(!self.value)
+    }
+}
+
+impl core::ops::BitAnd<Self> for &ConstantSignal<bool> {
+    type Output = ConstantSignal<bool>;
+
+    fn bitand(self, rhs: Self) -> Self::Output {
+        ConstantSignal::<bool>::new(self.value && rhs.value)
+    }
+}
+
+impl core::ops::BitAnd<&Signal<bool>> for &ConstantSignal<bool> {
+    type Output = Signal<bool>;
+
+    fn bitand(self, rhs: &Signal<bool>) -> Self::Output {
+        rhs & self
+    }
+}
+impl core::ops::BitOr<Self> for &ConstantSignal<bool> {
+    type Output = ConstantSignal<bool>;
+
+    fn bitor(self, rhs: Self) -> Self::Output {
+        ConstantSignal::<bool>::new(self.value && rhs.value)
+    }
+}
+
+impl core::ops::BitOr<&Signal<bool>> for &ConstantSignal<bool> {
+    type Output = Signal<bool>;
+
+    fn bitor(self, rhs: &Signal<bool>) -> Self::Output {
+        rhs | self
+    }
+}

argus-core/src/signals/cmp_ops.rs

@@ -0,0 +1,150 @@
+use std::{cmp::Ordering, time::Duration};
+
+use num_traits::NumCast;
+
+use crate::signals::{
+    utils::{find_intersection, Neighborhood},
+    InterpolationMethod, Sample,
+};
+
+use super::{
+    traits::{BaseSignal, LinearInterpolatable, SignalPartialOrd, SignalSyncPoints},
+    ConstantSignal, Signal,
+};
+
+fn sync_with_intersection<'a, T, Sig1, Sig2, F>(
+    sig1: &'a Sig1,
+    sig2: &'a Sig2,
+    sync_points: &[&'a Duration],
+    op: F,
+) -> Signal<bool>
+where
+    F: Fn(Ordering) -> bool,
+    T: PartialOrd + Copy + std::fmt::Debug + NumCast + LinearInterpolatable,
+    Sig1: BaseSignal<Value = T>,
+    Sig2: BaseSignal<Value = T>,
+{
+    // This has to be manually implemented and cannot use the apply2 functions.
+    // This is because if we have two signals that cross each other, then there is
+    // an intermediate point where the two signals are equal. This point must be
+    // added to the signal appropriately.
+    use Ordering::*;
+    // This will contain the new signal with an initial capacity of twice the input
+    // signals sample points (as that is the upper limit of the number of new points
+    // that will be added
+    let mut return_signal = Signal::<bool>::new_with_capacity(sync_points.len() * 2);
+    // this will contain the last sample point and ordering
+    let mut last_sample = None;
+    // We will now loop over the sync points, compare across signals and (if
+    // an intersection happens) we will have to compute the intersection point
+    for &t in sync_points {
+        let lhs = sig1.at(*t).expect("value must be present at given time");
+        let rhs = sig2.at(*t).expect("values must be present at given time");
+        let ord = lhs.partial_cmp(rhs).unwrap();
+
+        if let Some((tm1, last)) = last_sample {
+            // Check if the signals crossed, this will happen essentiall if the last
+            // and the current are opposites and were not Equal.
+            if let (Less, Greater) | (Greater, Less) = (last, ord) {
+                // Find the point of intersection between the points.
+                let a = Neighborhood {
+                    first: sig1.at(tm1).copied().map(|value| Sample { time: tm1, value }),
+                    second: sig1.at(*t).copied().map(|value| Sample { time: *t, value }),
+                };
+                let b = Neighborhood {
+                    first: sig2.at(tm1).copied().map(|value| Sample { time: tm1, value }),
+                    second: sig2.at(*t).copied().map(|value| Sample { time: *t, value }),
+                };
+                let intersect = find_intersection(&a, &b);
+                {
+                    let lhs = sig1
+                        .interpolate_at(intersect.time, InterpolationMethod::Linear)
+                        .unwrap();
+                    let rhs = sig2
+                        .interpolate_at(intersect.time, InterpolationMethod::Linear)
+                        .unwrap();
+                    assert_eq!(lhs, rhs);
+                }
+                return_signal
+                    .push(intersect.time, op(Equal))
+                    .expect("Signal should already be monotonic");
+            }
+        }
+        last_sample = Some((*t, ord));
+    }
+    return_signal.time_points.shrink_to_fit();
+    return_signal.values.shrink_to_fit();
+    return_signal
+}
+
+impl<T> SignalPartialOrd<Self> for Signal<T>
+where
+    T: PartialOrd + Copy + std::fmt::Debug + NumCast + LinearInterpolatable,
+{
+    type Output = Signal<bool>;
+
+    fn signal_cmp<F>(&self, other: &Self, op: F) -> Option<Self::Output>
+    where
+        F: Fn(Ordering) -> bool,
+    {
+        // This has to be manually implemented and cannot use the apply2 functions.
+        // This is because if we have two signals that cross each other, then there is
+        // an intermediate point where the two signals are equal. This point must be
+        // added to the signal appropriately.
+        // the union of the sample points in self and other
+        let sync_points = match self.synchronization_points(other) {
+            Some(points) => points,
+            None => return None,
+        };
+
+        Some(sync_with_intersection(self, other, &sync_points, op))
+    }
+}
+
+impl<T> SignalPartialOrd<ConstantSignal<T>> for Signal<T>
+where
+    T: PartialOrd + Copy + std::fmt::Debug + NumCast + LinearInterpolatable,
+{
+    type Output = Signal<bool>;
+
+    fn signal_cmp<F>(&self, other: &ConstantSignal<T>, op: F) -> Option<Self::Output>
+    where
+        F: Fn(Ordering) -> bool,
+    {
+        // the union of the sample points in self and other
+        let sync_points = match self.synchronization_points(other) {
+            Some(points) => points,
+            None => return None,
+        };
+
+        Some(sync_with_intersection(self, other, &sync_points, op))
+    }
+}
+
+impl<T> SignalPartialOrd<ConstantSignal<T>> for ConstantSignal<T>
+where
+    T: PartialOrd + Copy + std::fmt::Debug + NumCast,
+{
+    type Output = ConstantSignal<bool>;
+
+    fn signal_cmp<F>(&self, other: &ConstantSignal<T>, op: F) -> Option<Self::Output>
+    where
+        F: Fn(Ordering) -> bool,
+    {
+        self.value.partial_cmp(&other.value).map(op).map(ConstantSignal::new)
+    }
+}
+
+impl<T> SignalPartialOrd<Signal<T>> for ConstantSignal<T>
+where
+    T: PartialOrd + Copy + std::fmt::Debug + NumCast + LinearInterpolatable,
+{
+    type Output = Signal<bool>;
+
+    fn signal_cmp<F>(&self, other: &Signal<T>, op: F) -> Option<Self::Output>
+    where
+        F: Fn(Ordering) -> bool,
+    {
+        other.signal_cmp(self, op)
+    }
+}

argus-core/src/signals/num_ops.rs

@@ -0,0 +1,232 @@
+use num_traits::{Num, NumCast, Signed};
+
+use crate::signals::utils::{apply1, apply2, apply2_const};
+use crate::signals::{ConstantSignal, Signal};
+
+use super::traits::LinearInterpolatable;
+
+impl<T> core::ops::Neg for &Signal<T>
+where
+    T: Signed + Copy,
+{
+    type Output = Signal<T>;
+
+    /// Negate the signal at each time point
+    fn neg(self) -> Self::Output {
+        apply1(self, |v| -v)
+    }
+}
+
+impl<T> core::ops::Add for &Signal<T>
+where
+    T: core::ops::Add<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Add the given signal with another
+    fn add(self, rhs: Self) -> Self::Output {
+        apply2(self, rhs, |lhs, rhs| lhs + rhs)
+    }
+}
+
+impl<T> core::ops::Add<&ConstantSignal<T>> for &Signal<T>
+where
+    T: core::ops::Add<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Add the given signal with another
+    fn add(self, rhs: &ConstantSignal<T>) -> Self::Output {
+        apply2_const(self, rhs, |lhs, rhs| lhs + rhs)
+    }
+}
+
+impl<T> core::ops::Mul for &Signal<T>
+where
+    T: core::ops::Mul<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Multiply the given signal with another
+    fn mul(self, rhs: Self) -> Self::Output {
+        apply2(self, rhs, |lhs, rhs| lhs * rhs)
+    }
+}
+
+impl<T> core::ops::Mul<&ConstantSignal<T>> for &Signal<T>
+where
+    T: core::ops::Mul<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Multiply the given signal with another
+    fn mul(self, rhs: &ConstantSignal<T>) -> Self::Output {
+        apply2_const(self, rhs, |lhs, rhs| lhs * rhs)
+    }
+}
+
+impl<T> core::ops::Sub for &Signal<T>
+where
+    T: core::ops::Sub<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Subtract the given signal with another
+    fn sub(self, rhs: Self) -> Self::Output {
+        apply2(self, rhs, |lhs, rhs| lhs - rhs)
+    }
+}
+
+impl<T> core::ops::Sub<&ConstantSignal<T>> for &Signal<T>
+where
+    T: core::ops::Sub<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Subtiply the given signal with another
+    fn sub(self, rhs: &ConstantSignal<T>) -> Self::Output {
+        apply2_const(self, rhs, |lhs, rhs| lhs - rhs)
+    }
+}
+
+impl<T> core::ops::Div for &Signal<T>
+where
+    T: core::ops::Div<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Divide the given signal with another
+    fn div(self, rhs: Self) -> Self::Output {
+        apply2(self, rhs, |lhs, rhs| lhs / rhs)
+    }
+}
+impl<T> core::ops::Div<&ConstantSignal<T>> for &Signal<T>
+where
+    T: core::ops::Div<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Divide the given signal with another
+    fn div(self, rhs: &ConstantSignal<T>) -> Self::Output {
+        apply2_const(self, rhs, |lhs, rhs| lhs / rhs)
+    }
+}
+
+impl<T> num_traits::Pow<Self> for &Signal<T>
+where
+    T: num_traits::Pow<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Returns the values in `self` to the power of the values in `other`
+    fn pow(self, other: Self) -> Self::Output {
+        apply2(self, other, |lhs, rhs| lhs.pow(rhs))
+    }
+}
+
+impl<T> core::ops::Neg for &ConstantSignal<T>
+where
+    T: Signed + Copy,
+{
+    type Output = ConstantSignal<T>;
+
+    /// Negate the signal at each time point
+    fn neg(self) -> Self::Output {
+        ConstantSignal::new(self.value.neg())
+    }
+}
+
+impl<T> core::ops::Add for &ConstantSignal<T>
+where
+    T: core::ops::Add<T, Output = T> + Num + Copy,
+{
+    type Output = ConstantSignal<T>;
+
+    /// Add the given signal with another
+    fn add(self, rhs: Self) -> Self::Output {
+        ConstantSignal::<T>::new(self.value + rhs.value)
+    }
+}
+
+impl<T> core::ops::Add<&Signal<T>> for &ConstantSignal<T>
+where
+    T: core::ops::Add<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Add the given signal with another
+    fn add(self, rhs: &Signal<T>) -> Self::Output {
+        rhs + self
+    }
+}
+
+impl<T> core::ops::Mul for &ConstantSignal<T>
+where
+    T: core::ops::Mul<T, Output = T> + Num + Copy,
+{
+    type Output = ConstantSignal<T>;
+
+    /// Multiply the given signal with another
+    fn mul(self, rhs: Self) -> Self::Output {
+        ConstantSignal::<T>::new(self.value * rhs.value)
+    }
+}
+
+impl<T> core::ops::Mul<&Signal<T>> for &ConstantSignal<T>
+where
+    T: core::ops::Mul<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Multiply the given signal with another
+    fn mul(self, rhs: &Signal<T>) -> Self::Output {
+        rhs * self
+    }
+}
+
+impl<T> core::ops::Sub for &ConstantSignal<T>
+where
+    T: core::ops::Sub<T, Output = T> + Num + Copy,
+{
+    type Output = ConstantSignal<T>;
+
+    /// Subtract the given signal with another
+    fn sub(self, rhs: Self) -> Self::Output {
+        ConstantSignal::<T>::new(self.value - rhs.value)
+    }
+}
+impl<T> core::ops::Sub<&Signal<T>> for &ConstantSignal<T>
+where
+    T: core::ops::Sub<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Subtract the given signal with another
+    fn sub(self, rhs: &Signal<T>) -> Self::Output {
+        apply2_const(rhs, self, |rhs, lhs| lhs - rhs)
+    }
+}
+
+impl<T> core::ops::Div for &ConstantSignal<T>
+where
+    T: core::ops::Div<T, Output = T> + Num + Copy,
+{
+    type Output = ConstantSignal<T>;
+
+    /// Divide the given signal with another
+    fn div(self, rhs: Self) -> Self::Output {
+        ConstantSignal::<T>::new(self.value / rhs.value)
+    }
+}
+
+impl<T> core::ops::Div<&Signal<T>> for &ConstantSignal<T>
+where
+    T: core::ops::Div<T, Output = T> + Num + NumCast + Copy + LinearInterpolatable,
+{
+    type Output = Signal<T>;
+
+    /// Divide the given signal with another
+    fn div(self, rhs: &Signal<T>) -> Self::Output {
+        apply2_const(rhs, self, |rhs, lhs| lhs / rhs)
+    }
+}

argus-core/src/signals/traits.rs

@@ -1,9 +1,13 @@
-use std::ops::RangeBounds;
+use paste::paste;
+use std::cmp::Ordering;
 use std::time::Duration;
+use std::{iter::Empty, ops::RangeBounds};
 
+use itertools::Itertools;
 use num_traits::Num;
 
-use super::{InterpolationMethod, Sample};
+use super::{ConstantSignal, InterpolationMethod, Sample, Signal};
+use crate::signals::utils::intersect_bounds;
 use crate::ArgusResult;
 
 /// A general Signal trait
@@ -157,3 +161,168 @@ interpolate_for_num!(u32);
 interpolate_for_num!(u64);
 interpolate_for_num!(f32);
 interpolate_for_num!(f64);
+
+pub trait SignalSamplePoints {
+    type Output<'a>: IntoIterator<Item = &'a Duration>
+    where
+        Self: 'a;
+
+    /// Get the time points where the signal is sampled.
+    fn time_points(&'_ self) -> Option<Self::Output<'_>>;
+}
+
+pub trait SignalSyncPoints<Rhs = Self> {
+    type Output<'a>: IntoIterator<Item = &'a Duration>
+    where
+        Self: 'a,
+        Rhs: 'a;
+
+    /// Return the union list of time points where each of the given signals is sampled.
+    fn synchronization_points<'a>(&'a self, other: &'a Rhs) -> Option<Self::Output<'a>>;
+}
+
+impl<T> SignalSamplePoints for Signal<T>
+where
+    Signal<T>: BaseSignal,
+    T: Copy,
+{
+    type Output<'a> = Vec<&'a Duration>
+    where
+        Self: 'a;
+
+    fn time_points(&'_ self) -> Option<Self::Output<'_>> {
+        if self.is_empty() {
+            None
+        } else {
+            self.time_points.iter().collect_vec().into()
+        }
+    }
+}
+
+impl<T> SignalSamplePoints for ConstantSignal<T>
+where
+    T: Copy,
+{
+    type Output<'a> = Empty<&'a Duration>
+    where
+        Self: 'a;
+
+    fn time_points(&'_ self) -> Option<Self::Output<'_>> {
+        if self.is_empty() {
+            None
+        } else {
+            core::iter::empty().into()
+        }
+    }
+}
+
+impl<T> SignalSyncPoints<Self> for Signal<T>
+where
+    T: Copy,
+    Self: BaseSignal<Value = T>,
+{
+    type Output<'a> = Vec<&'a Duration>
+    where
+        Self: 'a,
+        Self: 'a;
+
+    fn synchronization_points<'a>(&'a self, other: &'a Self) -> Option<Self::Output<'a>> {
+        use core::ops::Bound::*;
+        if self.is_empty() || other.is_empty() {
+            return None;
+        }
+
+        let bounds = match intersect_bounds(&self.bounds(), &other.bounds()) {
+            (Included(start), Included(end)) => start..=end,
+            (..) => unreachable!(),
+        };
+
+        self.time_points
+            .iter()
+            .merge(other.time_points.iter())
+            .filter(|time| bounds.contains(time))
+            .dedup()
+            .collect_vec()
+            .into()
+    }
+}
+
+impl<T> SignalSyncPoints<ConstantSignal<T>> for Signal<T>
+where
+    T: Copy,
+    Self: BaseSignal<Value = T>,
+{
+    type Output<'a> = Vec<&'a Duration>
+    where
+        Self: 'a,
+        Self: 'a;
+
+    fn synchronization_points<'a>(&'a self, other: &'a ConstantSignal<T>) -> Option<Self::Output<'a>> {
+        if self.is_empty() || other.is_empty() {
+            return None;
+        }
+
+        self.time_points.iter().collect_vec().into()
+    }
+}
+
+impl<T> SignalSyncPoints<ConstantSignal<T>> for ConstantSignal<T>
+where
+    T: Copy,
+    Self: BaseSignal<Value = T>,
+{
+    type Output<'a> = Empty<&'a Duration>
+    where
+        Self: 'a,
+        Self: 'a;
+
+    fn synchronization_points<'a>(&'a self, _other: &'a ConstantSignal<T>) -> Option<Self::Output<'a>> {
+        Some(core::iter::empty())
+    }
+}
+
+impl<T> SignalSyncPoints<Signal<T>> for ConstantSignal<T>
+where
+    T: Copy,
+    Self: BaseSignal<Value = T>,
+{
+    type Output<'a> = Vec<&'a Duration>
+    where
+        Self: 'a,
+        Self: 'a;
+
+    fn synchronization_points<'a>(&'a self, other: &'a Signal<T>) -> Option<Self::Output<'a>> {
+        other.synchronization_points(self)
+    }
+}
+
+macro_rules! impl_signal_cmp {
+    ($cmp:ident) => {
+        paste! {
+            fn [<signal_ $cmp>](&self, other: &Rhs) -> Option<Self::Output> {
+                self.signal_cmp(other, |ord| ord.[<is_ $cmp>]())
+            }
+        }
+    };
+}
+
+/// A time-wise partial ordering defined for signals
+pub trait SignalPartialOrd<Rhs = Self>: BaseSignal {
+    type Output: BaseSignal<Value = bool>;
+
+    /// Compare two signals within each of their domains (using [`PartialOrd`]) and
+    /// apply the given function `op` to the ordering to create a signal.
+    ///
+    /// This function returns `None` if the comparison isn't possible, namely, when
+    /// either of the signals are empty.
+    fn signal_cmp<F>(&self, other: &Rhs, op: F) -> Option<Self::Output>
+    where
+        F: Fn(Ordering) -> bool;
+
+    impl_signal_cmp!(lt);
+    impl_signal_cmp!(le);
+    impl_signal_cmp!(gt);
+    impl_signal_cmp!(ge);
+    impl_signal_cmp!(eq);
+    impl_signal_cmp!(ne);
+}

argus-core/src/signals/utils.rs

@@ -0,0 +1,345 @@
+//! A bunch of utility code for argus
+//!
+//! - The implementation for Range intersection is based on the library
+//! [`range_ext`](https://github.com/AnickaBurova/range-ext), but adapted for my use a
+//! bit.
+
+use core::ops::{Bound, RangeBounds};
+use core::time::Duration;
+
+use num_traits::NumCast;
+
+use super::traits::{LinearInterpolatable, SignalSyncPoints};
+use super::{BaseSignal, ConstantSignal, InterpolationMethod, Sample, Signal};
+
+/// The neighborhood around a signal such that the time `at` is between the `first` and
+/// `second` samples.
+///
+/// The values of `first` and `second` are `None` if and only if `at` lies outside the
+/// domain over which the signal is defined.
+///
+/// This can be used to interpolate the value at the given `at` time using strategies
+/// like constant previous, constant following, and linear interpolation.
+#[derive(Copy, Clone, Debug)]
+pub struct Neighborhood<T: ?Sized + Copy> {
+    pub first: Option<Sample<T>>,
+    pub second: Option<Sample<T>>,
+}
+
+/// Given two signals with two sample points each, find the intersection of the two
+/// lines.
+pub fn find_intersection<T>(a: &Neighborhood<T>, b: &Neighborhood<T>) -> Sample<T>
+where
+    T: Copy + std::fmt::Debug + NumCast,
+{
+    // https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line
+    use num_traits::cast;
+
+    let Sample { time: t1, value: y1 } = a.first.unwrap();
+    let Sample { time: t2, value: y2 } = a.second.unwrap();
+    let Sample { time: t3, value: y3 } = b.first.unwrap();
+    let Sample { time: t4, value: y4 } = b.second.unwrap();
+
+    let t1 = t1.as_secs_f64();
+    let t2 = t2.as_secs_f64();
+    let t3 = t3.as_secs_f64();
+    let t4 = t4.as_secs_f64();
+
+    let y1: f64 = cast(y1).unwrap();
+    let y2: f64 = cast(y2).unwrap();
+    let y3: f64 = cast(y3).unwrap();
+    let y4: f64 = cast(y4).unwrap();
+
+    let denom = ((t1 - t2) * (y3 - y4)) - ((y1 - y2) * (t3 - t4));
+
+    let t_top = (((t1 * y2) - (y1 * t2)) * (t3 - t4)) - ((t1 - t2) * (t3 * y4 - y3 * t4));
+    let y_top = (((t1 * y2) - (y1 * t2)) * (y3 - y4)) - ((y1 - y2) * (t3 * y4 - y3 * t4));
+
+    let t = Duration::from_secs_f64(t_top / denom);
+    let y: T = cast(y_top / denom).unwrap();
+    Sample { time: t, value: y }
+}
+
+pub fn apply1<T, F>(signal: &Signal<T>, op: F) -> Signal<T>
+where
+    T: Copy,
+    F: Fn(T) -> T,
+    Signal<T>: std::iter::FromIterator<(Duration, T)>,
+{
+    signal.iter().map(|(t, v)| (*t, op(*v))).collect()
+}
+
+pub fn apply2<'a, T, U, F>(lhs: &'a Signal<T>, rhs: &'a Signal<T>, op: F) -> Signal<U>
+where
+    T: Copy + LinearInterpolatable,
+    U: Copy,
+    F: Fn(T, T) -> U,
+{
+    // If either of the signals are empty, we return an empty signal.
+    if lhs.is_empty() || rhs.is_empty() {
+        // Intersection with empty signal should yield an empty signal
+        return Signal::<U>::new();
+    }
+    // We determine the range of the signal (as the output signal can only be
+    // defined in the domain where both signals are defined).
+    let time_points = lhs.synchronization_points(rhs).unwrap();
+    // Now, at each of the merged time points, we sample each signal and operate on
+    // them
+    time_points
+        .into_iter()
+        .map(|t| {
+            let v1 = lhs.interpolate_at(*t, InterpolationMethod::Linear).unwrap();
+            let v2 = rhs.interpolate_at(*t, InterpolationMethod::Linear).unwrap();
+            (*t, op(v1, v2))
+        })
+        .collect()
+}
+
+pub fn apply2_const<'a, T, U, F>(lhs: &'a Signal<T>, rhs: &'a ConstantSignal<T>, op: F) -> Signal<U>
+where
+    T: Copy + LinearInterpolatable,
+    U: Copy,
+    F: Fn(T, T) -> U,
+{
+    // If either of the signals are empty, we return an empty signal.
+    if lhs.is_empty() {
+        // Intersection with empty signal should yield an empty signal
+        return Signal::<U>::new();
+    }
+    lhs.time_points
+        .iter()
+        .map(|&t| {
+            let v1 = lhs.interpolate_at(t, InterpolationMethod::Linear).unwrap();
+            let v2 = rhs.interpolate_at(t, InterpolationMethod::Linear).unwrap();
+            (t, op(v1, v2))
+        })
+        .collect()
+}
+
+fn partial_min<T>(a: T, b: T) -> Option<T>
+where
+    T: PartialOrd,
+{
+    a.partial_cmp(&b).map(|ord| if ord.is_lt() { a } else { b })
+}
+
+fn partial_max<T>(a: T, b: T) -> Option<T>
+where
+    T: PartialOrd,
+{
+    a.partial_cmp(&b).map(|ord| if ord.is_gt() { a } else { b })
+}
+
+/// Compute the intersection of two ranges
+pub fn intersect_bounds<T>(lhs: &impl RangeBounds<T>, rhs: &impl RangeBounds<T>) -> (Bound<T>, Bound<T>)
+where
+    T: PartialOrd + Copy,
+{
+    use core::ops::Bound::*;
+
+    let start = match (lhs.start_bound(), rhs.start_bound()) {
+        (Included(&l), Included(&r)) => Included(partial_max(l, r).unwrap()),
+        (Excluded(&l), Excluded(&r)) => Excluded(partial_max(l, r).unwrap()),
+
+        (Included(l), Excluded(r)) | (Excluded(r), Included(l)) => {
+            if l > r {
+                Included(*l)
+            } else {
+                Excluded(*r)
+            }
+        }
+
+        (Unbounded, Included(&l)) | (Included(&l), Unbounded) => Included(l),
+        (Unbounded, Excluded(&l)) | (Excluded(&l), Unbounded) => Excluded(l),
+        (Unbounded, Unbounded) => Unbounded,
+    };
+
+    let end = match (lhs.end_bound(), rhs.end_bound()) {
+        (Included(&l), Included(&r)) => Included(partial_min(l, r).unwrap()),
+        (Excluded(&l), Excluded(&r)) => Excluded(partial_min(l, r).unwrap()),
+
+        (Included(l), Excluded(r)) | (Excluded(r), Included(l)) => {
+            if l < r {
+                Included(*l)
+            } else {
+                Excluded(*r)
+            }
+        }
+
+        (Unbounded, Included(&l)) | (Included(&l), Unbounded) => Included(l),
+        (Unbounded, Excluded(&l)) | (Excluded(&l), Unbounded) => Excluded(l),
+        (Unbounded, Unbounded) => Unbounded,
+    };
+
+    (start, end)
+}
+
+/// More precise intersection of two ranges from point of the first range
+#[derive(Debug, PartialEq)]
+pub enum Intersection {
+    /// The self is below the other
+    Below,
+    /// The self is below but overlaping
+    BelowOverlap,
+    /// The self is within the other
+    Within,
+    /// The self is same as the other
+    Same,
+    /// The self is over the other, the other is within the self
+    Over,
+    /// The self is above but overlaping
+    AboveOverlap,
+    /// The self is above the other
+    Above,
+}
+
+impl Intersection {
+    /// Test if there is any intersection
+    pub fn is_any(&self) -> bool {
+        !matches!(self, Intersection::Below | Intersection::Above)
+    }
+    /// Test if the range is fully within the other
+    pub fn is_within(&self) -> bool {
+        matches!(self, Intersection::Within | Intersection::Same)
+    }
+
+    /// Test if the range is fully over the other
+    pub fn is_over(&self) -> bool {
+        matches!(self, Intersection::Over | Intersection::Same)
+    }
+}
+
+pub trait Intersect<T, U>
+where
+    T: PartialOrd,
+    U: RangeBounds<T>,
+{
+    /// Test two ranges for an intersection
+    fn check_intersect(&self, other: &U) -> Intersection;
+}
+
+impl<T, U, R> Intersect<T, U> for R
+where
+    T: PartialOrd + PartialEq,
+    U: RangeBounds<T>,
+    R: RangeBounds<T>,
+{
+    fn check_intersect(&self, other: &U) -> Intersection {
+        use core::cmp::Ordering::*;
+        use core::ops::Bound::*;
+
+        // We find where the start of self is with respect to that of other
+        let (left_rel_pos, me_start) = match (self.start_bound(), other.start_bound()) {
+            (Included(me), Excluded(them)) if me == them => (Less, Some(me)), // [a, _} left of (a, }
+            //
+            (Excluded(me), Included(them)) if me == them => (Greater, Some(me)), // (a, _} right of [a, }
+
+            // If both are consistently open or close, or they are not equal then we
+            // just compare them
+            (Included(me), Excluded(them))
+            | (Excluded(me), Included(them))
+            | (Included(me), Included(them))
+            | (Excluded(me), Excluded(them)) => (me.partial_cmp(them).unwrap(), Some(me)),
+
+            // start of self > start of other
+            (Included(me), Unbounded) | (Excluded(me), Unbounded) => (Greater, Some(me)),
+
+            (Unbounded, Unbounded) => (Equal, None), // unbounded start
+
+            (Unbounded, _) => (Less, None), // start of self < start of other
+        };
+
+        // We find where the end of self is with respect to that of other
+        let (right_rel_pos, me_end) = match (self.end_bound(), other.end_bound()) {
+            (Included(me), Excluded(them)) if me == them => (Greater, Some(me)), // {_, a] right of {_, a)
+
+            (Excluded(me), Included(them)) if me == them => (Less, Some(me)), // {_, a) right of {_, a]
+
+            // If both are consistently open or close, or they are not equal then we just compare them
+            (Included(me), Excluded(them))
+            | (Excluded(me), Included(them))
+            | (Included(me), Included(them))
+            | (Excluded(me), Excluded(them)) => (me.partial_cmp(them).unwrap(), Some(me)),
+
+            (Included(me), Unbounded) | (Excluded(me), Unbounded) => (Less, Some(me)), // end of self < end of other
+
+            (Unbounded, Unbounded) => (Equal, None), // unbounded end
+
+            (Unbounded, _) => (Greater, None), // end of self > end of other
+        };
+
+        // We have gotten the relative position of the ends. But we need to check if one
+        // of the ends are contained within the bounds of the other.
+
+        match (left_rel_pos, right_rel_pos) {
+            (Less, Less) => {
+                // Check if the end of self is contained within other's domain
+                // NOTE: Since right is less than, me_end must not be None
+                assert!(me_end.is_some());
+                if other.contains(me_end.unwrap()) {
+                    // self is below but overlaps
+                    Intersection::BelowOverlap
+                } else {
+                    // self is strictly below
+                    Intersection::Below
+                }
+            }
+            (Greater, Greater) => {
+                // Check if the start of self is contained within other's domain
+                // NOTE: Since left is greater than, me_start must not be None
+                assert!(me_start.is_some());
+                if other.contains(me_start.unwrap()) {
+                    // self is to the right of but overlaps other
+                    Intersection::AboveOverlap
+                } else {
+                    // self is strictly above
+                    Intersection::Above
+                }
+            }
+            (Less, Greater) | (Equal, Greater) | (Less, Equal) => Intersection::Over, // self contains other
+            (Equal, Less) | (Greater, Equal) | (Greater, Less) => Intersection::Within, // self within other
+            (Equal, Equal) => Intersection::Same,                                     // The ranges are equal
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use std::ops::Range;
+
+    use proptest::prelude::*;
+
+    use super::*;
+
+    proptest! {
+        #[test]
+        fn range_range_intersection(r1 in any::<Range<i32>>(), r2 in any::<Range<i32>>()) {
+            use Intersection::*;
+            let intersect_p = r1.check_intersect(&r2);
+            match intersect_p {
+                Below => assert!(r1.end < r2.start, "Expected strict below"),
+                BelowOverlap => {
+                    assert!(r1.start < r2.start, "Expected below with overlap");
+                    assert!(r2.contains(&r1.end), "Expected below with overlap");
+                },
+                Within => {
+                    assert!(r2.contains(&r1.end), "Expected to be contained");
+                    assert!(r2.contains(&r1.start), "Expected to be contained");
+                }
+                Same => {
+                    assert!(r1.start == r2.start, "Expected to be same");
+                    assert!(r1.end == r2.end, "Expected to be same");
+                }
+                Over => {
+                    assert!(r1.contains(&r2.start), "Expected to cover");
+                    assert!(r1.contains(&r2.end), "Expected to cover");
+                }
+                AboveOverlap => {
+                    assert!(r2.contains(&r1.start), "Expected above with overlap");
+                    assert!(r1.end > r2.end, "Expected above with overlap");
+                }
+                Above => assert!(r1.start > r2.end, "Expected strict above"),
+            }
+        }
+    }
+}