amari_calculus/

lie.rs

//! Lie derivatives and Lie brackets
//!
//! This module provides Lie derivatives - derivatives along vector fields that measure
//! how fields change along the flow of a vector field.
//!
//! ## Mathematical Background
//!
//! The Lie derivative generalizes the directional derivative to arbitrary tensor fields:
//!
//! ### Lie Bracket of Vector Fields
//!
//! For vector fields X and Y, the Lie bracket is:
//! ```text
//! [X, Y] = ∇_X Y - ∇_Y X
//! ```
//!
//! In coordinates:
//! ```text
//! [X, Y]^i = X^j ∂_j Y^i - Y^j ∂_j X^i
//! ```
//!
//! ### Lie Derivative
//!
//! - **Scalar field**: L_X(f) = X(f) = X^i ∂_i f (directional derivative)
//! - **Vector field**: L_X(Y) = [X, Y]
//! - **Multivector field**: L_X(F) extends via Leibniz rule
//!
//! ## Examples
//!
//! ```
//! use amari_calculus::{VectorField, LieDerivative, vector_from_slice};
//!
//! // Define rotation field X = (-y, x, 0)
//! let x = VectorField::<3, 0, 0>::new(|coords| {
//!     vector_from_slice(&[-coords[1], coords[0], 0.0])
//! });
//!
//! // Define radial field Y = (x, y, 0)
//! let y = VectorField::<3, 0, 0>::new(|coords| {
//!     vector_from_slice(&[coords[0], coords[1], 0.0])
//! });
//!
//! let lie = LieDerivative::<3, 0, 0>::new();
//!
//! // Compute Lie bracket [X, Y] at origin
//! let bracket = lie.bracket(&x, &y, &[0.0, 0.0, 0.0]);
//! ```

use crate::fields::*;
use amari_core::Multivector;

/// Lie derivative operator
///
/// Computes Lie derivatives and Lie brackets of fields along vector fields.
pub struct LieDerivative<const P: usize, const Q: usize, const R: usize> {
    /// Step size for numerical differentiation
    h: f64,
}

impl<const P: usize, const Q: usize, const R: usize> LieDerivative<P, Q, R> {
    /// Create new Lie derivative operator
    ///
    /// # Examples
    ///
    /// ```
    /// use amari_calculus::LieDerivative;
    ///
    /// let lie = LieDerivative::<3, 0, 0>::new();
    /// ```
    pub fn new() -> Self {
        Self { h: 1e-5 }
    }

    /// Set step size for numerical differentiation
    pub fn with_step_size(mut self, h: f64) -> Self {
        self.h = h;
        self
    }

    /// Compute Lie derivative of scalar field along vector field: L_X(f)
    ///
    /// This is equivalent to the directional derivative X(f).
    ///
    /// # Arguments
    ///
    /// * `x` - Vector field along which to differentiate
    /// * `f` - Scalar field to differentiate
    /// * `coords` - Point at which to evaluate
    ///
    /// # Returns
    ///
    /// Scalar value L_X(f) = X^i ∂_i f
    ///
    /// # Examples
    ///
    /// ```
    /// use amari_calculus::{ScalarField, VectorField, LieDerivative, vector_from_slice};
    ///
    /// // f(x, y) = x² + y²
    /// let f = ScalarField::<3, 0, 0>::new(|coords| {
    ///     coords[0].powi(2) + coords[1].powi(2)
    /// });
    ///
    /// // X = (1, 0, 0) - unit vector in x direction
    /// let x = VectorField::<3, 0, 0>::new(|_coords| {
    ///     vector_from_slice(&[1.0, 0.0, 0.0])
    /// });
    ///
    /// let lie = LieDerivative::<3, 0, 0>::new();
    ///
    /// // L_X(f) at (3, 4) should be ∂f/∂x = 2x = 6
    /// let result = lie.of_scalar(&x, &f, &[3.0, 4.0, 0.0]);
    /// ```
    pub fn of_scalar(
        &self,
        x: &VectorField<P, Q, R>,
        f: &ScalarField<P, Q, R>,
        coords: &[f64],
    ) -> f64 {
        let dim = P + Q + R;
        let x_val = x.evaluate(coords);
        let mut result = 0.0;

        // L_X(f) = X^i ∂_i f
        for i in 0..dim {
            let x_i = x_val.vector_component(i);
            if let Ok(df_di) = f.partial_derivative(coords, i, self.h) {
                result += x_i * df_di;
            }
        }

        result
    }

    /// Compute Lie bracket of two vector fields: [X, Y]
    ///
    /// The Lie bracket measures the non-commutativity of the flows generated by X and Y.
    ///
    /// # Arguments
    ///
    /// * `x` - First vector field
    /// * `y` - Second vector field
    /// * `coords` - Point at which to evaluate
    ///
    /// # Returns
    ///
    /// Vector field [X, Y]^i = X^j ∂_j Y^i - Y^j ∂_j X^i
    ///
    /// # Examples
    ///
    /// ```
    /// use amari_calculus::{VectorField, LieDerivative, vector_from_slice};
    ///
    /// // Coordinate vector fields
    /// let e1 = VectorField::<3, 0, 0>::new(|_coords| {
    ///     vector_from_slice(&[1.0, 0.0, 0.0])
    /// });
    /// let e2 = VectorField::<3, 0, 0>::new(|_coords| {
    ///     vector_from_slice(&[0.0, 1.0, 0.0])
    /// });
    ///
    /// let lie = LieDerivative::<3, 0, 0>::new();
    ///
    /// // [e1, e2] should be zero (coordinate fields commute)
    /// let bracket = lie.bracket(&e1, &e2, &[1.0, 1.0, 0.0]);
    /// ```
    pub fn bracket(
        &self,
        x: &VectorField<P, Q, R>,
        y: &VectorField<P, Q, R>,
        coords: &[f64],
    ) -> Multivector<P, Q, R> {
        let dim = P + Q + R;
        let x_val = x.evaluate(coords);
        let y_val = y.evaluate(coords);

        let mut bracket = Multivector::zero();

        // [X, Y]^i = X^j ∂_j Y^i - Y^j ∂_j X^i
        for i in 0..dim {
            let mut component = 0.0;

            // First term: X^j ∂_j Y^i
            for j in 0..dim {
                let x_j = x_val.vector_component(j);
                if let Ok(dy_dj) = y.partial_derivative(coords, j, self.h) {
                    component += x_j * dy_dj.vector_component(i);
                }
            }

            // Second term: -Y^j ∂_j X^i
            for j in 0..dim {
                let y_j = y_val.vector_component(j);
                if let Ok(dx_dj) = x.partial_derivative(coords, j, self.h) {
                    component -= y_j * dx_dj.vector_component(i);
                }
            }

            bracket.set_vector_component(i, component);
        }

        bracket
    }

    /// Compute Lie derivative of vector field along another vector field: L_X(Y) = [X, Y]
    ///
    /// This is an alias for the Lie bracket.
    ///
    /// # Arguments
    ///
    /// * `x` - Vector field along which to differentiate
    /// * `y` - Vector field to differentiate
    /// * `coords` - Point at which to evaluate
    ///
    /// # Returns
    ///
    /// Vector field L_X(Y) = [X, Y]
    pub fn of_vector(
        &self,
        x: &VectorField<P, Q, R>,
        y: &VectorField<P, Q, R>,
        coords: &[f64],
    ) -> Multivector<P, Q, R> {
        self.bracket(x, y, coords)
    }

    /// Compute Lie derivative of multivector field along vector field: L_X(F)
    ///
    /// For a general multivector field, the Lie derivative extends via the Leibniz rule.
    ///
    /// # Arguments
    ///
    /// * `x` - Vector field along which to differentiate
    /// * `f` - Multivector field to differentiate
    /// * `coords` - Point at which to evaluate
    ///
    /// # Returns
    ///
    /// Multivector field L_X(F)
    pub fn of_multivector(
        &self,
        x: &VectorField<P, Q, R>,
        f: &MultivectorField<P, Q, R>,
        coords: &[f64],
    ) -> Multivector<P, Q, R> {
        let dim = P + Q + R;
        let x_val = x.evaluate(coords);

        // L_X(F) = X^i ∂_i F (directional derivative for now)
        // Full implementation would account for geometric product structure
        let mut result = Multivector::zero();

        for i in 0..dim {
            let x_i = x_val.vector_component(i);
            if let Ok(df_di) = f.partial_derivative(coords, i, self.h) {
                result = result + df_di * x_i;
            }
        }

        result
    }
}

impl<const P: usize, const Q: usize, const R: usize> Default for LieDerivative<P, Q, R> {
    fn default() -> Self {
        Self::new()
    }
}

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

    #[test]
    fn test_lie_derivative_of_scalar() {
        // f(x, y) = x² + y²
        let f = ScalarField::<3, 0, 0>::new(|coords| coords[0].powi(2) + coords[1].powi(2));

        // X = (1, 0, 0)
        let x = VectorField::<3, 0, 0>::new(|_coords| crate::vector_from_slice(&[1.0, 0.0, 0.0]));

        let lie = LieDerivative::<3, 0, 0>::new();
        let result = lie.of_scalar(&x, &f, &[3.0, 4.0, 0.0]);

        // L_X(f) = X · ∇f = (1,0,0) · (2x, 2y, 0) = 2x = 6
        assert!(
            (result - 6.0).abs() < 1e-4,
            "Lie derivative should be 6.0, got {}",
            result
        );
    }

    #[test]
    fn test_lie_bracket_coordinate_fields() {
        // Coordinate vector fields commute, so bracket should be zero
        let e1 = VectorField::<3, 0, 0>::new(|_coords| crate::vector_from_slice(&[1.0, 0.0, 0.0]));
        let e2 = VectorField::<3, 0, 0>::new(|_coords| crate::vector_from_slice(&[0.0, 1.0, 0.0]));

        let lie = LieDerivative::<3, 0, 0>::new();
        let bracket = lie.bracket(&e1, &e2, &[1.0, 1.0, 0.0]);

        // All components should be approximately zero
        for i in 0..3 {
            assert!(
                bracket.vector_component(i).abs() < 1e-4,
                "Coordinate fields should commute, got component {} = {}",
                i,
                bracket.vector_component(i)
            );
        }
    }

    #[test]
    fn test_lie_bracket_noncommuting_fields() {
        // X = (y, 0, 0)
        let x =
            VectorField::<3, 0, 0>::new(|coords| crate::vector_from_slice(&[coords[1], 0.0, 0.0]));

        // Y = (0, x, 0)
        let y =
            VectorField::<3, 0, 0>::new(|coords| crate::vector_from_slice(&[0.0, coords[0], 0.0]));

        let lie = LieDerivative::<3, 0, 0>::new();

        // ∂_x Y = (0, 1, 0), ∂_y Y = (0, 0, 0)
        // ∂_x X = (0, 0, 0), ∂_y X = (1, 0, 0)
        // At point (2, 3, 0): X = (3, 0, 0), Y = (0, 2, 0)
        // [X, Y]^x = X^x ∂_x Y^x + X^y ∂_y Y^x - (Y^x ∂_x X^x + Y^y ∂_y X^x)
        //          = 3*0 + 0*0 - (0*0 + 2*1)
        //          = 0 - 2 = -2
        // [X, Y]^y = X^x ∂_x Y^y + X^y ∂_y Y^y - (Y^x ∂_x X^y + Y^y ∂_y X^y)
        //          = 3*1 + 0*0 - (0*0 + 2*0)
        //          = 3 - 0 = 3
        let bracket = lie.bracket(&x, &y, &[2.0, 3.0, 0.0]);

        assert!(
            (bracket.vector_component(0) + 2.0).abs() < 1e-4,
            "x-component should be -2, got {}",
            bracket.vector_component(0)
        );
        assert!(
            (bracket.vector_component(1) - 3.0).abs() < 1e-4,
            "y-component should be 3, got {}",
            bracket.vector_component(1)
        );
    }

    #[test]
    fn test_lie_bracket_antisymmetry() {
        // [X, Y] = -[Y, X]
        let x = VectorField::<3, 0, 0>::new(|coords| {
            crate::vector_from_slice(&[-coords[1], coords[0], 0.0])
        });
        let y = VectorField::<3, 0, 0>::new(|coords| {
            crate::vector_from_slice(&[coords[0], coords[1], 0.0])
        });

        let lie = LieDerivative::<3, 0, 0>::new();
        let coords = [1.0, 2.0, 0.0];

        let bracket_xy = lie.bracket(&x, &y, &coords);
        let bracket_yx = lie.bracket(&y, &x, &coords);

        for i in 0..3 {
            let sum = bracket_xy.vector_component(i) + bracket_yx.vector_component(i);
            assert!(
                sum.abs() < 1e-4,
                "Antisymmetry violated: [X,Y]_{} + [Y,X]_{} = {}",
                i,
                i,
                sum
            );
        }
    }
}