symplex 0.42.0

//! Dimensional calculus: differentiation and integration that track
//! physical dimensions through the type system.
//!
//! Differentiation divides dimensions: `d(Qty<D1>)/d(Qty<D2>) → Qty<D1/D2>`
//! Integration multiplies dimensions: `∫ Qty<D1> d(Qty<D2>) → Qty<D1*D2>`
//!
//! Named newtypes use `.as_qty()` before calling these functions,
//! then convert back via `.into()`:
//!
//! ```ignore
//! let ctx = Context::new();
//! let x = Length::symbol(&ctx, "x");
//! let t = Time::symbol(&ctx, "t");
//! let v: Velocity = diff_qty(&x.as_qty(), &t.as_qty()).into();
//! ```

use std::ops;

use typenum::operator_aliases::{Diff, Sum};

use super::dim::Dim;
use super::qty::Qty;
use crate::prelude::Ex;

/// The symbol a quantity variable is measured in, with `dv/ds`, when the
/// variable is not itself a symbol: a unit conversion of one symbol
/// (`Time::minutes(τ)` is `60·τ`, `Temperature::from_celsius(ϑ)` is
/// `ϑ + 273.15`).  `None` for a bare symbol, a constant, or several
/// symbols.
fn measured_symbol(var: &Ex) -> Option<(Ex, Ex)> {
    if is_symbol(var) {
        return None;
    }
    let [s] = <[Ex; 1]>::try_from(var.free_symbols()).ok()?;
    let dv = var.diff(&s);
    (!dv.is_zero_structural()).then_some((s, dv))
}

fn is_symbol(e: &Ex) -> bool {
    matches!(
        e.inner.read().arena.node(e.raw_id()),
        crate::base::node::ExprNode::Symbol(_)
    )
}

/// `d expr / d var` for a quantity variable: the ordinary derivative for a
/// symbol, the chain rule `(d expr/ds) / (dv/ds)` for a variable that is a
/// function of one symbol `s`, else the unevaluated `Derivative(expr, var)`
/// (differentiating with respect to a non-symbol used to give `0`).
pub(crate) fn diff_wrt_quantity(expr: &Ex, var: &Ex) -> Ex {
    if is_symbol(var) {
        return expr.diff(var);
    }
    match measured_symbol(var) {
        Some((s, dv)) => expr.diff(&s) / &dv,
        None => expr.formal_diff(var),
    }
}

/// `∫ expr d var` for a quantity variable: substitution `∫ expr·(dv/ds) ds`
/// for a variable that is a function of one symbol `s` (integrating with
/// respect to a non-symbol used to stay unevaluated).
pub(crate) fn integrate_wrt_quantity(expr: &Ex, var: &Ex) -> Ex {
    match measured_symbol(var) {
        Some((s, dv)) => (expr * &dv).integrate(&s),
        None => expr.integrate(var),
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// diff_qty — symbolic differentiation with dimension tracking
// ═══════════════════════════════════════════════════════════════════════════

/// Differentiate a dimensioned expression with respect to a dimensioned variable.
///
/// The output dimension is the quotient of the two input dimensions:
/// `D_expr / D_var` (exponents are subtracted).
///
/// The variable may be a unit conversion of one symbol — `Time::minutes(τ)`
/// is `60·τ` seconds — and the derivative is then taken by the chain rule,
/// `(d expr/dτ) / 60`; a variable with several symbols gives the
/// unevaluated `Derivative(expr, var)`.
///
/// # Type-level mechanics
///
/// Given `expr: Qty<Dim<L1, M1, T1, …>>` and `var: Qty<Dim<L2, M2, T2, …>>`,
/// the result has dimension `Dim<L1−L2, M1−M2, T1−T2, …>`.
///
/// # Examples
///
/// ```ignore
/// use symplex::units::*;
///
/// let ctx = Context::new();
/// let x: Qty<LengthDim> = Qty::from_ex(ctx.symbol("x"));
/// let t: Qty<TimeDim>   = Qty::from_ex(ctx.symbol("t"));
/// let v: Qty<VelocityDim> = diff_qty(&x, &t);   // Length / Time = Velocity
/// ```
#[allow(clippy::type_complexity)]
pub fn diff_qty<L1, M1, T1, I1, Th1, N1x, J1, L2, M2, T2, I2, Th2, N2x, J2>(
    expr: &Qty<Dim<L1, M1, T1, I1, Th1, N1x, J1>>,
    var: &Qty<Dim<L2, M2, T2, I2, Th2, N2x, J2>>,
) -> Qty<
    Dim<
        Diff<L1, L2>,
        Diff<M1, M2>,
        Diff<T1, T2>,
        Diff<I1, I2>,
        Diff<Th1, Th2>,
        Diff<N1x, N2x>,
        Diff<J1, J2>,
    >,
>
where
    L1: ops::Sub<L2>,
    M1: ops::Sub<M2>,
    T1: ops::Sub<T2>,
    I1: ops::Sub<I2>,
    Th1: ops::Sub<Th2>,
    N1x: ops::Sub<N2x>,
    J1: ops::Sub<J2>,
{
    Qty::from_ex(diff_wrt_quantity(expr.inner(), var.inner()))
}

// ═══════════════════════════════════════════════════════════════════════════
// integrate_qty — symbolic integration with dimension tracking
// ═══════════════════════════════════════════════════════════════════════════

/// Integrate a dimensioned expression with respect to a dimensioned variable.
///
/// The output dimension is the product of the two input dimensions:
/// `D_expr × D_var` (exponents are added).
///
/// A variable that is a unit conversion of one symbol (`Length::kilometers(ξ)`
/// is `1000·ξ` metres) is integrated by substitution: `∫ expr·1000 dξ`.
///
/// # Type-level mechanics
///
/// Given `expr: Qty<Dim<L1, M1, T1, …>>` and `var: Qty<Dim<L2, M2, T2, …>>`,
/// the result has dimension `Dim<L1+L2, M1+M2, T1+T2, …>`.
///
/// # Examples
///
/// ```ignore
/// use symplex::units::*;
///
/// let ctx = Context::new();
/// let f: Qty<ForceDim>  = Qty::from_ex(ctx.symbol("F"));
/// let x: Qty<LengthDim> = Qty::from_ex(ctx.symbol("x"));
/// let w: Qty<EnergyDim> = integrate_qty(&f, &x);  // Force × Length = Energy
/// ```
#[allow(clippy::type_complexity)]
pub fn integrate_qty<L1, M1, T1, I1, Th1, N1x, J1, L2, M2, T2, I2, Th2, N2x, J2>(
    expr: &Qty<Dim<L1, M1, T1, I1, Th1, N1x, J1>>,
    var: &Qty<Dim<L2, M2, T2, I2, Th2, N2x, J2>>,
) -> Qty<
    Dim<
        Sum<L1, L2>,
        Sum<M1, M2>,
        Sum<T1, T2>,
        Sum<I1, I2>,
        Sum<Th1, Th2>,
        Sum<N1x, N2x>,
        Sum<J1, J2>,
    >,
>
where
    L1: ops::Add<L2>,
    M1: ops::Add<M2>,
    T1: ops::Add<T2>,
    I1: ops::Add<I2>,
    Th1: ops::Add<Th2>,
    N1x: ops::Add<N2x>,
    J1: ops::Add<J2>,
{
    Qty::from_ex(integrate_wrt_quantity(expr.inner(), var.inner()))
}

// ═══════════════════════════════════════════════════════════════════════════
// Tests
// ═══════════════════════════════════════════════════════════════════════════

#[cfg(test)]
mod tests {
    use super::super::dim::*;
    use super::super::si::*;
    use super::*;

    // --- Differentiation tests ---

    /// d(Length)/d(Time) = Velocity
    #[test]
    fn diff_length_by_time_is_velocity() {
        let ctx = crate::api::context::Context::new();
        let x: Qty<LengthDim> = Qty::from_ex(ctx.symbol("x"));
        let t: Qty<TimeDim> = Qty::from_ex(ctx.symbol("t"));

        let result = diff_qty(&x, &t);

        // If dimensions are wrong, this `.into()` won't compile.
        let _v: Velocity = result.into();
    }

    /// d(Velocity)/d(Time) = Acceleration
    #[test]
    fn diff_velocity_by_time_is_acceleration() {
        let ctx = crate::api::context::Context::new();
        let v: Qty<VelocityDim> = Qty::from_ex(ctx.symbol("v"));
        let t: Qty<TimeDim> = Qty::from_ex(ctx.symbol("t"));

        let result = diff_qty(&v, &t);

        let _a: Acceleration = result.into();
    }

    /// d(Energy)/d(Length) = Force
    #[test]
    fn diff_energy_by_length_is_force() {
        let ctx = crate::api::context::Context::new();
        let e: Qty<EnergyDim> = Qty::from_ex(ctx.symbol("E"));
        let x: Qty<LengthDim> = Qty::from_ex(ctx.symbol("x"));

        let result = diff_qty(&e, &x);

        let _f: Force = result.into();
    }

    // --- Integration tests ---

    /// ∫ Force d(Length) = Energy
    #[test]
    fn integrate_force_over_length_is_energy() {
        let ctx = crate::api::context::Context::new();
        let f: Qty<ForceDim> = Qty::from_ex(ctx.symbol("F"));
        let x: Qty<LengthDim> = Qty::from_ex(ctx.symbol("x"));

        let result = integrate_qty(&f, &x);

        let _w: Energy = result.into();
    }

    /// ∫ Velocity d(Time) = Length
    #[test]
    fn integrate_velocity_over_time_is_length() {
        let ctx = crate::api::context::Context::new();
        let v: Qty<VelocityDim> = Qty::from_ex(ctx.symbol("v"));
        let t: Qty<TimeDim> = Qty::from_ex(ctx.symbol("t"));

        let result = integrate_qty(&v, &t);

        let _x: Length = result.into();
    }

    // --- Fundamental Theorem of Calculus roundtrip ---

    /// diff(integrate(v, t), t) must have Velocity dimension.
    ///
    /// This verifies that the type-level arithmetic is consistent:
    /// Velocity×Time = Length, then Length/Time = Velocity.
    #[test]
    fn ftc_roundtrip_velocity() {
        let ctx = crate::api::context::Context::new();
        let v: Qty<VelocityDim> = Qty::from_ex(ctx.symbol("v"));
        let t: Qty<TimeDim> = Qty::from_ex(ctx.symbol("t"));

        // integrate: Velocity × Time = Length
        let integrated = integrate_qty(&v, &t);

        // differentiate: Length / Time = Velocity
        let t2: Qty<TimeDim> = Qty::from_ex(ctx.symbol("t"));
        let roundtrip = diff_qty(&integrated, &t2);

        // Must compile as Velocity — proves FTC dimensional consistency.
        let _v2: Velocity = roundtrip.into();
    }

    // --- Named-type workflow via as_qty ---

    /// Demonstrate the full named → Qty → calculus → named workflow.
    #[test]
    fn named_type_workflow() {
        let ctx = crate::api::context::Context::new();
        let x = Length::symbol(&ctx, "x");
        let t = Time::symbol(&ctx, "t");

        let v: Velocity = diff_qty(&x.as_qty(), &t.as_qty()).into();

        // The result is a proper Velocity; we can use named-type methods.
        assert_eq!(Velocity::dim_name_str(), "Velocity");
        assert_eq!(Velocity::dim_symbol_str(), "m/s");
        let _ = v;
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Typed calculus: DiffWrt and IntWrt traits
// ═══════════════════════════════════════════════════════════════════════════

/// Typed differentiation. Differentiating `Self` with respect to `Var`
/// produces `Self::Output` with the correct physical dimension.
pub trait DiffWrt<Var> {
    /// The output quantity type after differentiation.
    type Output;
    /// Differentiate this quantity with respect to the given variable.
    fn diff_wrt(&self, var: &Var) -> Self::Output;
}

/// Typed integration. Integrating `Self` with respect to `Var`
/// produces `Self::Output` with the correct physical dimension.
pub trait IntWrt<Var> {
    /// The output quantity type after integration.
    type Output;
    /// Integrate this quantity with respect to the given variable.
    fn integrate_wrt(&self, var: &Var) -> Self::Output;
}

macro_rules! impl_diff_wrt {
    ($Expr:ty, $Var:ty => $Out:ty) => {
        impl DiffWrt<$Var> for $Expr {
            type Output = $Out;
            fn diff_wrt(&self, var: &$Var) -> $Out {
                <$Out>::from_ex(self.inner().diff(var.inner()))
            }
        }
    };
}

macro_rules! impl_int_wrt {
    ($Expr:ty, $Var:ty => $Out:ty) => {
        impl IntWrt<$Var> for $Expr {
            type Output = $Out;
            fn integrate_wrt(&self, var: &$Var) -> $Out {
                <$Out>::from_ex(self.inner().integrate(var.inner()))
            }
        }
    };
}

use super::si::*;

// ── Kinematics: d/dt ──
impl_diff_wrt!(Length, Time => Velocity);
impl_diff_wrt!(Velocity, Time => Acceleration);
impl_diff_wrt!(Angle, Time => AngularVelocity);
impl_diff_wrt!(AngularVelocity, Time => AngularAcceleration);

// ── Energy/Power: d/dt ──
impl_diff_wrt!(Energy, Time => Power);
impl_diff_wrt!(Momentum, Time => Force);
impl_diff_wrt!(AngularMomentum, Time => Torque);
impl_diff_wrt!(Charge, Time => Current);
impl_diff_wrt!(MagneticFlux, Time => Voltage);

// ── Spatial derivatives ──
impl_diff_wrt!(Energy, Length => Force);
impl_diff_wrt!(Energy, Angle => Torque);
impl_diff_wrt!(Momentum, Length => Stiffness); // dp/dx in wave context
impl_diff_wrt!(Force, Length => Stiffness);

// ── Energy w.r.t. generalized velocities ──
impl_diff_wrt!(Energy, Velocity => Momentum);
impl_diff_wrt!(Energy, AngularVelocity => AngularMomentum);

// ── Electrical ──
impl_diff_wrt!(Power, Current => Voltage);
impl_diff_wrt!(Power, Voltage => Current);

// ── Integration (reverse of differentiation) ──
impl_int_wrt!(Velocity, Time => Length);
impl_int_wrt!(Acceleration, Time => Velocity);
impl_int_wrt!(AngularVelocity, Time => Angle);
impl_int_wrt!(AngularAcceleration, Time => AngularVelocity);
impl_int_wrt!(Power, Time => Energy);
impl_int_wrt!(Force, Time => Momentum);
impl_int_wrt!(Torque, Time => AngularMomentum);
impl_int_wrt!(Current, Time => Charge);
impl_int_wrt!(Voltage, Time => MagneticFlux);
impl_int_wrt!(Force, Length => Energy);
impl_int_wrt!(Stiffness, Length => Force);
impl_int_wrt!(Momentum, Velocity => Energy); // ∫p dv = KE
impl_int_wrt!(AngularMomentum, AngularVelocity => Energy);

// ── Angle-based derivatives (for robotics Jacobians) ──
// Differentiating w.r.t. Angle preserves dimension (Angle is dimensionless in SI)
impl_diff_wrt!(Length, Angle => Length);
impl_diff_wrt!(Velocity, Angle => Velocity);
impl_diff_wrt!(Acceleration, Angle => Acceleration);
impl_diff_wrt!(Force, Angle => Force);
impl_diff_wrt!(Torque, Angle => Torque);
impl_diff_wrt!(AngularMomentum, Angle => AngularMomentum);
impl_diff_wrt!(MomentOfInertia, Angle => MomentOfInertia);

// Angle-based integration
impl_int_wrt!(Torque, Angle => Energy); // ∫τ dθ = work done
impl_int_wrt!(Length, Angle => Length); // Arc length integration
impl_int_wrt!(Force, Angle => Force);
impl_int_wrt!(AngularMomentum, Angle => AngularMomentum);