alice/groups/
matrix_group.rs

1use alloc::vec::Vec;
2use faer::Mat;
3use faer::prelude::Solve;
4use crate::maps::faer_bridge;
5
6#[derive(Debug, Clone, PartialEq, Eq)]
7pub struct MatrixGroup {
8    data: Vec<OrderedF64>,
9    n: usize,
10}
11
12#[derive(Debug, Clone, Copy, PartialEq)]
13pub struct OrderedF64(pub f64);
14
15impl Eq for OrderedF64 {}
16
17#[derive(Debug, Clone, PartialEq, Eq)]
18pub enum GroupError {
19    DimensionMismatch,
20    NotInvertible,
21    ConstraintViolated(&'static str),
22}
23
24impl core::fmt::Display for GroupError {
25    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
26        match self {
27            GroupError::DimensionMismatch => write!(f, "matrix dimension mismatch"),
28            GroupError::NotInvertible => write!(f, "matrix is not invertible"),
29            GroupError::ConstraintViolated(msg) => write!(f, "group constraint violated: {}", msg),
30        }
31    }
32}
33
34impl MatrixGroup {
35    pub fn new(data: Vec<f64>, n: usize) -> Result<Self, GroupError> {
36        if data.len() != n * n {
37            return Err(GroupError::DimensionMismatch);
38        }
39        Ok(Self {
40            data: data.into_iter().map(OrderedF64).collect(),
41            n,
42        })
43    }
44
45    pub fn identity(n: usize) -> Self {
46        let data = (0..n * n)
47            .map(|k| OrderedF64(if k / n == k % n { 1.0 } else { 0.0 }))
48            .collect();
49        Self { data, n }
50    }
51
52    pub fn n(&self) -> usize {
53        self.n
54    }
55
56    pub fn data(&self) -> Vec<f64> {
57        self.data.iter().map(|x| x.0).collect()
58    }
59
60    pub fn data_ref(&self) -> &[OrderedF64] {
61        &self.data
62    }
63
64    pub fn get(&self, i: usize, j: usize) -> f64 {
65        self.data[i * self.n + j].0
66    }
67
68    pub fn mul(&self, other: &Self) -> Result<Self, GroupError> {
69        if self.n != other.n {
70            return Err(GroupError::DimensionMismatch);
71        }
72        let a = to_faer_inner(&self.data(), self.n);
73        let b = to_faer_inner(&other.data(), self.n);
74        let c = a * b;
75        Ok(Self {
76            data: from_faer_inner(&c).into_iter().map(OrderedF64).collect(),
77            n: self.n,
78        })
79    }
80
81    pub fn inverse(&self) -> Result<Self, GroupError> {
82        let a = to_faer_inner(&self.data(), self.n);
83        let det = a.determinant();
84        if libm::fabs(det) < 1e-14 {
85            return Err(GroupError::NotInvertible);
86        }
87        let lu = a.partial_piv_lu();
88        let identity = Mat::<f64>::identity(self.n, self.n);
89        let inv = lu.solve(identity);
90        Ok(Self {
91            data: from_faer_inner(&inv).into_iter().map(OrderedF64).collect(),
92            n: self.n,
93        })
94    }
95
96    pub fn det(&self) -> f64 {
97        let a = to_faer_inner(&self.data(), self.n);
98        a.determinant()
99    }
100
101    pub fn transpose(&self) -> Self {
102        let mut data = alloc::vec![OrderedF64(0.0); self.n * self.n];
103        for i in 0..self.n {
104            for j in 0..self.n {
105                data[j * self.n + i] = OrderedF64(self.data[i * self.n + j].0);
106            }
107        }
108        Self { data, n: self.n }
109    }
110
111    pub fn exp(&self) -> Self {
112        let result = faer_bridge::matrix_exp(&self.data(), self.n);
113        Self {
114            data: result.into_iter().map(OrderedF64).collect(),
115            n: self.n,
116        }
117    }
118
119    pub fn log(&self) -> Option<Self> {
120        let result = faer_bridge::matrix_log(&self.data(), self.n)?;
121        Some(Self {
122            data: result.into_iter().map(OrderedF64).collect(),
123            n: self.n,
124        })
125    }
126}
127
128pub(crate) fn to_faer_inner(data: &[f64], n: usize) -> Mat<f64> {
129    Mat::from_fn(n, n, |i, j| data[i * n + j])
130}
131
132pub(crate) fn from_faer_inner(m: &Mat<f64>) -> Vec<f64> {
133    let n = m.nrows();
134    let mut out = alloc::vec![0.0f64; n * n];
135    for i in 0..n {
136        for j in 0..n {
137            out[i * n + j] = m[(i, j)];
138        }
139    }
140    out
141}
alice/groups/matrix_group.rs

alice/groups/
matrix_group.rs
