devela/num/geom/linear/

affine.rs

// devela::num::geom::linear::affine
//
//! Defines [`Affine`].
//
// > https://chatgpt.com/c/67bcf245-e818-8007-8453-eab0292bef86
//
// TOC
// -
// -
// - compute_minor_determinant
//
// TODO
// - constructor
// - validator of dimensions
// - validated constructors
//
// - macro to implement the D transform… for f32 or f64, …
//
// - up to 8x8 matrices (UNDERSTAND for affine) calculate a list of which ones…
//
// IDEA: could this file be also a script?

use crate::Matrix;

/// An affine transformation in `D`-dimensional space.
#[derive(Clone, Copy, Debug, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct Affine<T, const D: usize, const LEN: usize> {
    arr: [T; LEN],
}
impl<T, const D: usize, const LEN: usize> Affine<T, D, LEN> {
    ///
    pub fn new(arr: [T; LEN]) -> Self {
        Self::validate_d_len();
        Self { arr }
    }
    // MAYBE? PROBLEM:
    // pub fn new_2d(arr: [T; 2*3]) -> Affine<T, 2, {2*3}> { Self { arr } }

    const fn validate_d_len() {
        assert!(LEN == D * (D + 1), "Invalid array length for given dimension");
    }
}

impl<const D: usize, const LEN: usize> Affine<f64, D, LEN> {
    /// Computes the array index for a matrix element (row, col)
    #[inline(always)]
    const fn idx(row: usize, col: usize) -> usize {
        row * (D + 1) + col
    }
    /// Access a linear transformation element
    #[inline(always)]
    pub const fn at(&self, row: usize, col: usize) -> f64 {
        self.arr[Self::idx(row, col)]
    }
    /// Mutably access a linear transformation element
    #[inline(always)]
    pub fn at_mut(&mut self, row: usize, col: usize) -> &mut f64 {
        &mut self.arr[Self::idx(row, col)]
    }
    /// Access a translation component
    #[inline(always)]
    pub fn translation(&self, row: usize) -> f64 {
        self.arr[Self::idx(row, D)]
    }
    /// Mutably access a translation component
    #[inline(always)]
    pub fn translation_mut(&mut self, row: usize) -> &mut f64 {
        &mut self.arr[Self::idx(row, D)]
    }

    /// Computes the determinant of the Jacobian (D×D linear transformation).
    // RETHINK: duplicating Matrix::determinant_unchecked, but what about the data?
    #[rustfmt::skip]
    pub fn determinant(&self) -> f64 {
        let s = self;
        match D {
            1 => s.at(0, 0),
            2 => s.at(0, 0) * s.at(1, 1) - s.at(0, 1) * s.at(1, 0),
            3 => {
                  s.at(0, 0) * (s.at(1, 1) * s.at(2, 2) - s.at(1, 2) * s.at(2, 1))
                - s.at(0, 1) * (s.at(1, 0) * s.at(2, 2) - s.at(1, 2) * s.at(2, 0))
                + s.at(0, 2) * (s.at(1, 0) * s.at(2, 1) - s.at(1, 1) * s.at(2, 0))
            }
            _ => {
                let mut buf = [0.0; {8*8}];
                Matrix::<f64, D, D, LEN>::submatrix_determinant(D, &self.arr, &mut buf)
            }
        }
    }
}

// TEMP: STANDALONE
// ///
// const fn compute_minor_determinant<const LEN: usize>(
//     n: usize,
//     matrix: &[f64],
//     minor_buffer: &mut [f64],
// ) -> f64 {
//     iif![n == 1; return matrix[0]];
//     iif![n == 2; return matrix[0] * matrix[3] - matrix[1] * matrix[2]];
//     let (mut det, mut col) = (0.0, 0);
//     while col < n {
//         let minor_size = (n - 1) * (n - 1);
//         extract_minor(n, 0, col, matrix, Slice::range_to_mut(minor_buffer, minor_size));
//         let mut minor_copy = [0.0; LEN];
//         let mut i = 0;
//         while i < minor_size {
//             minor_copy[i] = minor_buffer[i];
//             i += 1;
//         }
//         let minor_det = compute_minor_determinant::<LEN>(
//             n - 1,
//             Slice::range_to(&minor_copy, minor_size),
//             minor_buffer,
//         );
//         det += sign(col) * matrix[col] * minor_det;
//         col += 1;
//     }
//     det
// }
//
// /// Extracts an (N-1)x(N-1) minor matrix into `buffer`, using a **flattened slice**.
// const fn extract_minor(n: usize, row: usize, col: usize, matrix: &[f64], buffer: &mut [f64]) {
//     let (mut idx, mut r) = (0, 0);
//     while r < n {
//         iif![r == row; { r += 1; continue }];
//         let mut c = 0;
//         while c < n {
//             iif![c == col; { c += 1; continue }];
//             buffer[idx] = matrix[r * n + c];
//             idx += 1;
//             c += 1;
//         }
//         r += 1;
//     }
// }
//
// /// Returns the sign (-1 or 1) for alternating determinant sums.
// const fn sign(i: usize) -> f64 {
//     iif![i % 2 == 0; 1.0; -1.0]
// }

// TEMP moved to matrix
// // -----
//
// /// Row-major storage matrix
// pub struct Matrix<T, const ROWS: usize, const COLS: usize, const LEN: usize> {
//     data: [T; LEN],
// }
//
// // Generic struct only requires `T: Copy`
// impl<T: Copy, const ROWS: usize, const COLS: usize, const LEN: usize> Matrix<T, ROWS, COLS, LEN> {
//     /// Returns an immutable reference to the element at (row, col)
//     #[inline(always)]
//     pub const fn at(&self, row: usize, col: usize) -> T {
//         self.data[row * COLS + col]
//     }
//
//     /// Returns a mutable reference to the element at (row, col)
//     #[inline(always)]
//     pub fn at_mut(&mut self, row: usize, col: usize) -> &mut T {
//         &mut self.data[row * COLS + col]
//     }
// }
//
// // Specialized implementation for `f64`
// impl<const D: usize, const LEN: usize> Matrix<f64, D, D, LEN> {
//     /// Computes the determinant of the square matrix.
//     pub const fn determinant(&self) -> f64 {
//         match D {
//             1 => self.at(0, 0),
//             2 => self.at(0, 0) * self.at(1, 1) - self.at(0, 1) * self.at(1, 0),
//             3 => {
//                 self.at(0, 0) * (self.at(1, 1) * self.at(2, 2) - self.at(1, 2) * self.at(2, 1))
//                     - self.at(0, 1)
//                         * (self.at(1, 0) * self.at(2, 2) - self.at(1, 2) * self.at(2, 0))
//                     + self.at(0, 2)
//                         * (self.at(1, 0) * self.at(2, 1) - self.at(1, 1) * self.at(2, 0))
//             }
//             _ => {
//                 let mut minor_matrix = [0.0; 64]; // Supports up to 8×8 matrices
//                 Self::minor_determinant(D, &self.data, &mut minor_matrix)
//             }
//         }
//     }
//
//     /// Computes the determinant of a square matrix using Laplace expansion.
//     ///
//     /// This method is **only valid for square matrices** of size `dim × dim`
//     /// and is intended for cases where `dim > 3`.
//     ///
//     /// ### How It Works:
//     /// - Expands along the first row using minors (submatrices).
//     /// - Recursively computes determinants of `(dim-1)×(dim-1)` matrices.
//     /// - Uses `minor_buffer` to avoid extra allocations.
//     ///
//     /// ### Buffer Requirements:
//     /// - `minor_buffer` must have **at least** `(dim-1)² + (dim-2)²` elements.
//     /// - This ensures enough space for the current minor and recursive calls.
//     /// - A **debug assertion** checks this before execution.
//     ///
//     /// ### Complexity:
//     /// - **O(dim!)** (factorial time complexity) due to recursive expansion.
//     /// - Should only be used for **small matrices (`dim ≤ 8`)**.
//     ///
//     /// # Panics
//     /// Panics in debug mode if `minor_buffer` is too small.
//     pub const fn minor_determinant(dim: usize, matrix: &[f64], minor_buffer: &mut [f64]) -> f64 {
//         iif![dim == 1; return matrix[0]];
//         iif![dim == 2; return matrix[0] * matrix[3] - matrix[1] * matrix[2]];
//
//         let required_size = (dim - 1) * (dim - 1) + (dim - 2) * (dim - 2);
//         debug_assert!(minor_buffer.len() >= required_size, "minor_buffer is too small");
//
//         let (mut determinant, mut col_idx) = (0.0, 0);
//         while col_idx < dim {
//             let minor_len = (dim - 1) * (dim - 1);
//             // 1. Split buffer into current level & next (for recursion)
//             let (minor_matrix, next_minor_buffer) = minor_buffer.split_at_mut(minor_len);
//             // 2. Extract the minor submatrix
//             Self::minor_matrix(dim, 0, col_idx, matrix, minor_matrix);
//             // 3. Compute determinant recursively
//             let sub_det = Self::minor_determinant(dim - 1, minor_matrix, next_minor_buffer);
//             determinant += Self::parity_sign(col_idx) * matrix[col_idx] * sub_det;
//             col_idx += 1;
//         }
//         determinant
//     }
//
//     /// Extracts a `(D-1)x(D-1)` minor matrix by removing the given row and column.
//     ///
//     /// - `n`: Matrix dimension.
//     /// - `row`, `col`: The row and column to exclude.
//     /// - `matrix`: Source matrix in row-major order.
//     /// - `buffer`: Target buffer to store the minor matrix.
//     ///
//     /// Used internally by determinant computation.
//     pub const fn minor_matrix(
//         n: usize,
//         row: usize,
//         col: usize,
//         matrix: &[f64],
//         buffer: &mut [f64],
//     ) {
//         let (mut idx, mut r) = (0, 0);
//         while r < n {
//             iif![r == row; { r += 1; continue; }];
//             let mut c = 0;
//             while c < n {
//                 iif![c == col; { c += 1; continue; }];
//                 buffer[idx] = matrix[r * n + c];
//                 idx += 1;
//                 c += 1;
//             }
//             r += 1;
//         }
//     }
//
//     /// Returns alternating ±1 based on the column index for determinant expansion.
//     ///
//     /// - Returns `1.0` for even indices.
//     /// - Returns `-1.0` for odd indices.
//     ///
//     /// Used for Laplace expansion in determinant calculation.
//     pub const fn parity_sign(i: usize) -> f64 {
//         if i % 2 == 0 {
//             1.0
//         } else {
//             -1.0
//         }
//     }
// }