clifford_codegen/symbolic/

product.rs

//! Symbolic product computation.
//!
//! This module computes the symbolic output of product operations,
//! representing each output field as a symbolic expression of input fields.

use std::borrow::Cow;
use std::collections::HashMap;

use symbolica::atom::{Atom, DefaultNamespace};
use symbolica::parser::ParseSettings;

use crate::algebra::{Algebra, Blade, ProductTable};
use crate::spec::TypeSpec;

/// The kind of product to compute symbolically.
///
/// Product naming follows [Rigid Geometric Algebra](https://rigidgeometricalgebra.org/) conventions.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum ProductKind {
    /// Geometric product (full product).
    Geometric,
    /// Geometric antiproduct (complement(complement(a) × complement(b))).
    /// Used for versor composition with antisandwich-based transformations.
    Antigeometric,
    /// Wedge product (∧, exterior, grade-raising).
    Wedge,
    /// Left contraction (a ⌋ b, grade gb - ga when ga <= gb).
    LeftContraction,
    /// Right contraction (a ⌊ b, grade ga - gb when gb <= ga).
    RightContraction,
    /// Antiwedge product (∨, regressive/meet).
    Antiwedge,
    /// Bulk contraction (a ∨ b★).
    BulkContraction,
    /// Weight contraction (a ∨ b☆).
    WeightContraction,
    /// Bulk expansion (a ∧ b★).
    BulkExpansion,
    /// Weight expansion (a ∧ b☆).
    WeightExpansion,
    /// Dot product (• metric inner, same-grade only, returns scalar).
    Dot,
    /// Antidot product (⊚ metric antiproduct inner, same-antigrade only, returns scalar).
    Antidot,
    /// Scalar product (grade-0 projection of geometric product).
    Scalar,
    /// Projection: target ∨ (self ∧ target☆).
    Project,
    /// Antiprojection: target ∧ (self ∨ target☆).
    Antiproject,
}

/// A symbolic field expression.
#[derive(Debug, Clone)]
pub struct SymbolicField {
    /// Field name in the output type.
    pub name: String,
    /// Blade index this field represents.
    pub blade_index: usize,
    /// Symbolic expression for this field's value.
    pub expression: Atom,
}

/// Computes symbolic product outputs.
///
/// Given input types with symbolic field values, computes the symbolic
/// expressions for each output field.
pub struct SymbolicProduct {
    /// The product table.
    table: ProductTable,
}

impl SymbolicProduct {
    /// Creates a new symbolic product computer.
    pub fn new(algebra: &Algebra) -> Self {
        let table = ProductTable::new(algebra);
        Self { table }
    }

    /// Creates symbolic variables for a type's fields.
    ///
    /// Returns a map from field name to symbolic atom.
    pub fn create_field_symbols(&self, ty: &TypeSpec, prefix: &str) -> HashMap<String, Atom> {
        ty.fields
            .iter()
            .map(|f| {
                let symbol_name = format!("{}_{}", prefix, f.name);
                let input = DefaultNamespace {
                    namespace: Cow::Borrowed(env!("CARGO_CRATE_NAME")),
                    data: symbol_name.as_str(),
                    file: Cow::Borrowed(file!()),
                    line: line!() as usize,
                };
                let atom = Atom::parse(input, ParseSettings::symbolica()).unwrap();
                (f.name.clone(), atom)
            })
            .collect()
    }

    /// Computes the symbolic output of a product.
    ///
    /// # Arguments
    ///
    /// * `type_a` - Left operand type
    /// * `type_b` - Right operand type
    /// * `output_type` - Output type
    /// * `kind` - Product kind (geometric, outer, etc.)
    /// * `a_symbols` - Symbolic values for type_a fields
    /// * `b_symbols` - Symbolic values for type_b fields
    ///
    /// # Returns
    ///
    /// Symbolic expressions for each output field.
    pub fn compute(
        &self,
        type_a: &TypeSpec,
        type_b: &TypeSpec,
        output_type: &TypeSpec,
        kind: ProductKind,
        a_symbols: &HashMap<String, Atom>,
        b_symbols: &HashMap<String, Atom>,
    ) -> Vec<SymbolicField> {
        output_type
            .fields
            .iter()
            .map(|output_field| {
                let expr = self.compute_field(
                    type_a,
                    type_b,
                    output_field.blade_index,
                    kind,
                    a_symbols,
                    b_symbols,
                );
                // Apply output field sign for non-canonical blade ordering
                // (e.g., e31 = -e13, so sign = -1 and we negate the expression)
                let signed_expr = if output_field.sign < 0 { -expr } else { expr };
                SymbolicField {
                    name: output_field.name.clone(),
                    blade_index: output_field.blade_index,
                    expression: signed_expr,
                }
            })
            .collect()
    }

    /// Computes the symbolic output of a sandwich product: v × x × rev(v).
    ///
    /// # Arguments
    ///
    /// * `versor_type` - Versor type (Motor, Rotor, Flector)
    /// * `operand_type` - Operand type being transformed
    /// * `versor_symbols` - Symbolic values for versor fields
    /// * `operand_symbols` - Symbolic values for operand fields
    /// * `use_antiproduct` - If true, uses antiproduct and antireverse (for antisandwich)
    ///
    /// # Returns
    ///
    /// Symbolic expressions for each output field (output has same type as operand).
    pub fn compute_sandwich(
        &self,
        versor_type: &TypeSpec,
        operand_type: &TypeSpec,
        versor_symbols: &HashMap<String, Atom>,
        operand_symbols: &HashMap<String, Atom>,
        use_antiproduct: bool,
    ) -> Vec<SymbolicField> {
        operand_type
            .fields
            .iter()
            .map(|output_field| {
                let expr = self.compute_sandwich_field(
                    versor_type,
                    operand_type,
                    output_field,
                    versor_symbols,
                    operand_symbols,
                    use_antiproduct,
                );
                // Apply output field sign for non-canonical blade ordering
                let signed_expr = if output_field.sign < 0 { -expr } else { expr };
                SymbolicField {
                    name: output_field.name.clone(),
                    blade_index: output_field.blade_index,
                    expression: signed_expr,
                }
            })
            .collect()
    }

    /// Computes the symbolic expression for a single sandwich output field.
    ///
    /// Computes: Σ_{i,j,k} v_i × x_j × rev(v_k) for the target blade.
    fn compute_sandwich_field(
        &self,
        versor_type: &TypeSpec,
        operand_type: &TypeSpec,
        output_field: &crate::spec::FieldSpec,
        versor_symbols: &HashMap<String, Atom>,
        operand_symbols: &HashMap<String, Atom>,
        use_antiproduct: bool,
    ) -> Atom {
        let dim = self.table.dim();
        let result_blade = output_field.blade_index;
        let mut terms: Vec<Atom> = Vec::new();

        for field_v1 in &versor_type.fields {
            for field_x in &operand_type.fields {
                for field_v2 in &versor_type.fields {
                    let v1_blade = field_v1.blade_index;
                    let x_blade = field_x.blade_index;
                    let v2_blade = field_v2.blade_index;

                    // Compute v_i × x_j (or v_i ⊛ x_j)
                    let (sign_vx, vx) = if use_antiproduct {
                        self.table.antiproduct(v1_blade, x_blade)
                    } else {
                        self.table.geometric(v1_blade, x_blade)
                    };
                    if sign_vx == 0 {
                        continue;
                    }

                    // Compute the reverse/antireverse sign for v2
                    let v2_grade = Blade::from_index(v2_blade).grade();
                    let rev_sign: i8 = if use_antiproduct {
                        // Antireverse sign: (-1)^((n-k)(n-k-1)/2)
                        let antigrade = dim - v2_grade;
                        if (antigrade * antigrade.saturating_sub(1) / 2).is_multiple_of(2) {
                            1
                        } else {
                            -1
                        }
                    } else {
                        // Reverse sign: (-1)^(k(k-1)/2)
                        if (v2_grade * v2_grade.saturating_sub(1) / 2).is_multiple_of(2) {
                            1
                        } else {
                            -1
                        }
                    };

                    // Compute (v_i × x_j) × rev(v_k) (or (v_i ⊛ x_j) ⊛ antirev(v_k))
                    let (sign_vxr, result) = if use_antiproduct {
                        self.table.antiproduct(vx, v2_blade)
                    } else {
                        self.table.geometric(vx, v2_blade)
                    };
                    if sign_vxr == 0 {
                        continue;
                    }

                    if result != result_blade {
                        continue;
                    }

                    // Apply input field signs for non-canonical blade orderings.
                    // v1 and v2 are both from versor_type, x is from operand_type.
                    let input_sign = field_v1.sign * field_x.sign * field_v2.sign;
                    let final_sign = sign_vx * sign_vxr * rev_sign * input_sign;
                    if final_sign == 0 {
                        continue;
                    }

                    let v1_sym = versor_symbols.get(&field_v1.name).unwrap();
                    let x_sym = operand_symbols.get(&field_x.name).unwrap();
                    let v2_sym = versor_symbols.get(&field_v2.name).unwrap();

                    // Create term: sign * v1 * x * v2
                    let product = v1_sym * x_sym * v2_sym;
                    let term = if final_sign > 0 { product } else { -product };
                    terms.push(term);
                }
            }
        }

        if terms.is_empty() {
            Atom::num(0)
        } else {
            // Sort terms by string representation for deterministic output
            terms.sort_by_cached_key(|t| t.to_string());
            terms.into_iter().reduce(|acc, t| acc + t).unwrap()
        }
    }

    /// Computes the symbolic expression for a single output field.
    fn compute_field(
        &self,
        type_a: &TypeSpec,
        type_b: &TypeSpec,
        result_blade: usize,
        kind: ProductKind,
        a_symbols: &HashMap<String, Atom>,
        b_symbols: &HashMap<String, Atom>,
    ) -> Atom {
        let mut terms: Vec<Atom> = Vec::new();

        for field_a in &type_a.fields {
            for field_b in &type_b.fields {
                let a_blade = field_a.blade_index;
                let b_blade = field_b.blade_index;

                // Compute the product based on kind using ProductTable methods
                let (sign, result) = match kind {
                    ProductKind::Geometric => self.table.geometric(a_blade, b_blade),
                    ProductKind::Antigeometric => self.table.antiproduct(a_blade, b_blade),
                    ProductKind::Wedge => self.table.exterior(a_blade, b_blade),
                    ProductKind::LeftContraction => self.table.left_contraction(a_blade, b_blade),
                    ProductKind::RightContraction => self.table.right_contraction(a_blade, b_blade),
                    ProductKind::Antiwedge => self.table.regressive(a_blade, b_blade),
                    ProductKind::BulkContraction => self.table.bulk_contraction(a_blade, b_blade),
                    ProductKind::WeightContraction => {
                        self.table.weight_contraction(a_blade, b_blade)
                    }
                    ProductKind::BulkExpansion => self.table.bulk_expansion(a_blade, b_blade),
                    ProductKind::WeightExpansion => self.table.weight_expansion(a_blade, b_blade),
                    ProductKind::Dot => self.table.dot(a_blade, b_blade),
                    ProductKind::Antidot => self.table.antidot(a_blade, b_blade),
                    ProductKind::Scalar => {
                        // Scalar product: grade-0 projection of geometric product
                        let (s, r) = self.table.geometric(a_blade, b_blade);
                        // Only include if result is grade 0 (scalar blade index = 0)
                        let result_grade = Blade::from_index(r).grade();
                        if result_grade == 0 { (s, r) } else { (0, 0) }
                    }
                    ProductKind::Project => self.table.project(a_blade, b_blade),
                    ProductKind::Antiproject => self.table.antiproject(a_blade, b_blade),
                };

                if result != result_blade || sign == 0 {
                    continue;
                }

                let a_sym = a_symbols.get(&field_a.name).unwrap();
                let b_sym = b_symbols.get(&field_b.name).unwrap();

                // Apply input field signs for non-canonical blade orderings.
                // If a field is specified as e.g., e31 (sign = -1), the stored value
                // represents the negated canonical blade value.
                let input_sign = sign * field_a.sign * field_b.sign;

                // Create term: total_sign * a_field * b_field
                let product = a_sym * b_sym;
                let term = if input_sign > 0 { product } else { -product };
                terms.push(term);
            }
        }

        if terms.is_empty() {
            // Return symbolic zero
            Atom::num(0)
        } else {
            // Sort terms by string representation for deterministic output
            terms.sort_by_cached_key(|t| t.to_string());
            // Sum all terms
            terms.into_iter().reduce(|acc, t| acc + t).unwrap()
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::spec::parse_spec;
    use std::sync::Mutex;

    // Symbolica uses global state that conflicts when tests run in parallel.
    // Tests prefixed with `symbolica_` are configured to run serially via nextest.
    // The mutex provides a fallback for `cargo test` users.
    static SYMBOLICA_LOCK: Mutex<()> = Mutex::new(());

    #[test]
    fn symbolica_create_symbols_for_vector() {
        let _guard = SYMBOLICA_LOCK.lock().unwrap();
        let spec = parse_spec(
            r#"
            [algebra]
            name = "test"
            complete = false

            [signature]
            positive = ["e1", "e2", "e3"]

            [types.Vector]
            grades = [1]
            field_map = [
                { name = "x", blade = "e1" },
                { name = "y", blade = "e2" },
                { name = "z", blade = "e3" }
            ]
            "#,
        )
        .unwrap();

        let algebra = Algebra::euclidean(3);
        let product = SymbolicProduct::new(&algebra);
        let vector = spec.types.iter().find(|t| t.name == "Vector").unwrap();

        let symbols = product.create_field_symbols(vector, "a");

        assert!(symbols.contains_key("x"));
        assert!(symbols.contains_key("y"));
        assert!(symbols.contains_key("z"));
    }
}