qualia_core_db/solvers/calculus/
manifold.rs

1//! Finite-dimensional chart and Riemannian metric primitives.
2
3use super::analysis::{AnalysisError, Interval, LinearMap, Vector};
4
5#[repr(C)]
6#[derive(Debug, Clone, Copy, PartialEq)]
7pub struct Chart<const N: usize> {
8    pub chart_id: u64,
9    pub coordinate_domain: [Interval; N],
10}
11
12impl<const N: usize> Chart<N> {
13    pub fn contains(&self, point: Vector<N>) -> bool {
14        self.coordinate_domain
15            .iter()
16            .zip(point.data)
17            .all(|(interval, value)| interval.contains(value))
18    }
19}
20
21#[repr(C)]
22#[derive(Debug, Clone, Copy, PartialEq)]
23pub struct TransitionEvidence {
24    pub samples_checked: usize,
25    pub maximum_round_trip_error: f64,
26    pub tolerance: f64,
27}
28
29pub fn verify_transition<const N: usize, F, G>(
30    source: &Chart<N>,
31    target: &Chart<N>,
32    samples: &[Vector<N>],
33    forward: F,
34    inverse: G,
35    tolerance: f64,
36) -> Result<TransitionEvidence, AnalysisError>
37where
38    F: Fn(Vector<N>) -> Vector<N>,
39    G: Fn(Vector<N>) -> Vector<N>,
40{
41    if samples.is_empty() || !tolerance.is_finite() || tolerance <= 0.0 {
42        return Err(AnalysisError::InvalidDomain);
43    }
44    let mut maximum_error = 0.0_f64;
45    for sample in samples {
46        if !source.contains(*sample) {
47            return Err(AnalysisError::InvalidDomain);
48        }
49        let mapped = forward(*sample);
50        mapped.validate()?;
51        if !target.contains(mapped) {
52            return Err(AnalysisError::InvalidDomain);
53        }
54        let recovered = inverse(mapped);
55        let error = recovered.distance(*sample)?;
56        maximum_error = maximum_error.max(error);
57    }
58    if maximum_error > tolerance {
59        return Err(AnalysisError::NotCertified);
60    }
61    Ok(TransitionEvidence {
62        samples_checked: samples.len(),
63        maximum_round_trip_error: maximum_error,
64        tolerance,
65    })
66}
67
68#[repr(C)]
69#[derive(Debug, Clone, Copy, PartialEq)]
70pub struct RiemannMetric<const N: usize> {
71    pub coefficients: LinearMap<N, N>,
72}
73
74impl<const N: usize> RiemannMetric<N> {
75    pub fn new(coefficients: LinearMap<N, N>) -> Result<Self, AnalysisError> {
76        coefficients.validate()?;
77        for row in 0..N {
78            for column in 0..N {
79                let scale = coefficients.coefficients[row][column]
80                    .abs()
81                    .max(coefficients.coefficients[column][row].abs())
82                    .max(1.0);
83                if (coefficients.coefficients[row][column] - coefficients.coefficients[column][row])
84                    .abs()
85                    > 32.0 * f64::EPSILON * scale
86                {
87                    return Err(AnalysisError::InvalidDomain);
88                }
89            }
90        }
91        // Sylvester's criterion via leading principal determinants.
92        for order in 1..=N {
93            let mut block = [[0.0; N]; N];
94            for row in 0..order {
95                for column in 0..order {
96                    block[row][column] = coefficients.coefficients[row][column];
97                }
98            }
99            // Fill unused diagonal entries so the full determinant equals the
100            // leading-principal determinant.
101            for index in order..N {
102                block[index][index] = 1.0;
103            }
104            if LinearMap::new(block).determinant()? <= 0.0 {
105                return Err(AnalysisError::InvalidDomain);
106            }
107        }
108        Ok(Self { coefficients })
109    }
110
111    pub fn lower(&self, tangent: Vector<N>) -> Result<Vector<N>, AnalysisError> {
112        self.coefficients.apply(tangent)
113    }
114
115    pub fn raise(&self, covector: Vector<N>) -> Result<Vector<N>, AnalysisError> {
116        self.coefficients.solve(covector)
117    }
118
119    pub fn inner(&self, left: Vector<N>, right: Vector<N>) -> Result<f64, AnalysisError> {
120        left.dot(self.lower(right)?)
121    }
122
123    pub fn volume_density(&self) -> Result<f64, AnalysisError> {
124        Ok(self.coefficients.determinant()?.sqrt())
125    }
126}
127
128pub fn orientation_sign<const N: usize>(
129    transition_derivative: &LinearMap<N, N>,
130) -> Result<i8, AnalysisError> {
131    let determinant = transition_derivative.determinant()?;
132    if determinant > 0.0 {
133        Ok(1)
134    } else if determinant < 0.0 {
135        Ok(-1)
136    } else {
137        Err(AnalysisError::Singular)
138    }
139}
140
141#[cfg(test)]
142mod tests {
143    use super::*;
144
145    #[test]
146    fn chart_transition_round_trip_is_certified() {
147        let domain = [
148            Interval::new(-10.0, 10.0).unwrap(),
149            Interval::new(-10.0, 10.0).unwrap(),
150        ];
151        let source = Chart {
152            chart_id: 1,
153            coordinate_domain: domain,
154        };
155        let target = Chart {
156            chart_id: 2,
157            coordinate_domain: domain,
158        };
159        let samples = [
160            Vector::new([0.0, 0.0]),
161            Vector::new([1.0, -2.0]),
162            Vector::new([-3.0, 4.0]),
163        ];
164        let evidence = verify_transition(
165            &source,
166            &target,
167            &samples,
168            |x| Vector::new([x.data[0] + x.data[1], x.data[0] - x.data[1]]),
169            |y| Vector::new([0.5 * (y.data[0] + y.data[1]), 0.5 * (y.data[0] - y.data[1])]),
170            1e-14,
171        )
172        .unwrap();
173        assert_eq!(evidence.samples_checked, 3);
174    }
175
176    #[test]
177    fn musical_maps_and_volume_density_are_consistent() {
178        let metric = RiemannMetric::new(LinearMap::new([[4.0, 1.0], [1.0, 3.0]])).unwrap();
179        let vector = Vector::new([2.0, -1.0]);
180        let covector = metric.lower(vector).unwrap();
181        let recovered = metric.raise(covector).unwrap();
182        assert!(recovered.distance(vector).unwrap() < 1e-14);
183        assert!((metric.volume_density().unwrap() - 11.0_f64.sqrt()).abs() < 1e-14);
184        assert_eq!(
185            orientation_sign(&LinearMap::new([[0.0, 1.0], [1.0, 0.0]])),
186            Ok(-1)
187        );
188    }
189}
qualia_core_db/solvers/calculus/manifold.rs

qualia_core_db/solvers/calculus/
manifold.rs
