algebraeon_geometry/
affine_subspace.rs

1use super::*;
2use crate::{
3    ambient_space::AffineSpace,
4    oriented_simplex::{OrientedHyperplane, OrientedSimplex},
5    simplex::Simplex,
6    vector::Vector,
7};
8use algebraeon_rings::matrix::{Matrix, MatrixStructure};
9
10#[derive(Debug, Clone)]
11pub struct EmbeddedAffineSubspace<'f, FS: OrderedRingSignature + FieldSignature> {
12    // The ordered_field of ambient_space and subspace must match
13    ambient_space: AffineSpace<'f, FS>,
14    embedded_space: AffineSpace<'f, FS>,
15    /*
16    these vectors must be equal in length to the affine dimension of subspace
17    they define the embedding of subspace into ambient space
18        if they are empty then they define the empty embedding
19        if there is one vector then it defines the location of the embedded point
20        if there are vectors [v0, v1, v2, ..., vn] then the embedding sends (a1, a2, ..., an) in subspace to v0 + a1*v1, v0 + a2*v2, ..., v0 + an*vn in ambient space
21    */
22    embedding_points: Vec<Vector<'f, FS>>,
23}
24
25impl<'f, FS: OrderedRingSignature + FieldSignature> EmbeddedAffineSubspace<'f, FS> {
26    pub(crate) fn new_affine_independent_span(
27        ambient_space: AffineSpace<'f, FS>,
28        points: Vec<Vector<'f, FS>>,
29    ) -> Result<(Self, Vec<Vector<'f, FS>>), &'static str> {
30        for point in &points {
31            debug_assert_eq!(point.ambient_space(), ambient_space);
32        }
33        if !ambient_space
34            .borrow()
35            .are_points_affine_independent(points.iter().collect())
36        {
37            return Err("Affine embedding points must be affine independent");
38        }
39        let field = ambient_space.borrow().field();
40        let embedded_space: AffineSpace<'f, FS> = AffineSpace::new_affine(field, points.len());
41        let n = points.len();
42        let embedded_pts = (0..n)
43            .map(|i| {
44                Vector::construct(embedded_space, |j| {
45                    if i == 0 {
46                        field.borrow().zero()
47                    } else if i == j + 1 {
48                        field.borrow().one()
49                    } else {
50                        field.borrow().zero()
51                    }
52                })
53            })
54            .collect();
55        Ok((
56            Self {
57                ambient_space,
58                embedded_space,
59                embedding_points: points,
60            },
61            embedded_pts,
62        ))
63    }
64
65    pub(crate) fn new_empty(ambient_space: AffineSpace<'f, FS>) -> Self {
66        Self::new_affine_independent_span(ambient_space, vec![])
67            .unwrap()
68            .0
69    }
70}
71
72impl<'f, FS: OrderedRingSignature + FieldSignature> EmbeddedAffineSubspace<'f, FS> {
73    pub(crate) fn new_root_and_linear_independent_span(
74        ambient_space: AffineSpace<'f, FS>,
75        root: &Vector<'f, FS>,
76        span: Vec<&Vector<'f, FS>>,
77    ) -> Result<(Self, Vec<Vector<'f, FS>>), &'static str> {
78        let mut points = vec![root.clone()];
79        points.extend(span.into_iter().map(|vec| root + vec));
80        Self::new_affine_independent_span(ambient_space, points)
81    }
82
83    pub(crate) fn new_root_and_linear_span(
84        ambient_space: AffineSpace<'f, FS>,
85        root: &Vector<'f, FS>,
86        span: Vec<&Vector<'f, FS>>,
87    ) -> Self {
88        let mut points = vec![root.clone()];
89        points.extend(span.into_iter().map(|vec| root + vec));
90        Self::new_affine_span(ambient_space, points.iter().collect())
91    }
92
93    pub(crate) fn new_affine_span(
94        ambient_space: AffineSpace<'f, FS>,
95        points: Vec<&Vector<'f, FS>>,
96    ) -> Self {
97        if points.is_empty() {
98            Self::new_empty(ambient_space)
99        } else {
100            let dim = ambient_space.borrow().linear_dimension().unwrap();
101            let field = ambient_space.borrow().field();
102            let mut points = points.into_iter();
103            let root = points.next().unwrap();
104            let span = points.map(|pt| pt - root).collect::<Vec<_>>();
105            //matrix whose columns are pt - root for every other pt in points
106            let mat = Matrix::construct(dim, span.len(), |r, c| span[c].coordinate(r).clone());
107            let (_, _, _, pivs) = MatrixStructure::new(field.clone()).row_hermite_algorithm(mat);
108            Self::new_root_and_linear_independent_span(
109                ambient_space,
110                root,
111                pivs.into_iter().map(|i| &span[i]).collect(),
112            )
113            .unwrap()
114            .0
115        }
116    }
117}
118
119impl<'f, FS: OrderedRingSignature + FieldSignature> AffineSpace<'f, FS> {
120    pub fn affine_subspace_from_affine_span(
121        &self,
122        points: Vec<&Vector<'f, FS>>,
123    ) -> EmbeddedAffineSubspace<'f, FS> {
124        EmbeddedAffineSubspace::new_affine_span(*self, points)
125    }
126
127    pub fn affine_subspace_from_root_and_linear_span(
128        &self,
129        root: &Vector<'f, FS>,
130        span: Vec<&Vector<'f, FS>>,
131    ) -> EmbeddedAffineSubspace<'f, FS> {
132        EmbeddedAffineSubspace::new_root_and_linear_span(*self, root, span)
133    }
134
135    pub fn affine_subspace_from_affine_independent_span(
136        &self,
137        points: Vec<Vector<'f, FS>>,
138    ) -> Result<(EmbeddedAffineSubspace<'f, FS>, Vec<Vector<'f, FS>>), &'static str> {
139        EmbeddedAffineSubspace::new_affine_independent_span(*self, points)
140    }
141
142    pub fn affine_subspace_from_root_and_linear_independent_span(
143        &self,
144        root: &Vector<'f, FS>,
145        span: Vec<&Vector<'f, FS>>,
146    ) -> Result<(EmbeddedAffineSubspace<'f, FS>, Vec<Vector<'f, FS>>), &'static str> {
147        EmbeddedAffineSubspace::new_root_and_linear_independent_span(*self, root, span)
148    }
149
150    pub fn empty_affine_subspace(&self) -> EmbeddedAffineSubspace<'f, FS> {
151        EmbeddedAffineSubspace::new_empty(*self)
152    }
153}
154
155impl<'f, FS: OrderedRingSignature + FieldSignature> EmbeddedAffineSubspace<'f, FS> {
156    pub fn field(&self) -> &'f FS {
157        self.ambient_space.borrow().field()
158    }
159
160    pub fn ambient_space(&self) -> AffineSpace<'f, FS> {
161        self.ambient_space
162    }
163
164    pub fn embedded_space(&self) -> AffineSpace<'f, FS> {
165        self.embedded_space
166    }
167
168    //Let A be the affine subspace and let S be its ambient space
169    //Find an affine subspace B of S obtained by linearly extending by pt
170    //Return the embeddings (f, g, pt) where
171    //  f : A -> B
172    //  g : B -> S
173    //  pt is pt in B
174    pub(crate) fn extend_dimension_by_point_unsafe(
175        &self,
176        pt: Vector<'f, FS>,
177    ) -> (
178        EmbeddedAffineSubspace<'f, FS>,
179        EmbeddedAffineSubspace<'f, FS>,
180        Vector<'f, FS>,
181    ) {
182        debug_assert_eq!(self.ambient_space, pt.ambient_space());
183        debug_assert!(self.unembed_point(&pt).is_none());
184        let field = self.field();
185
186        let n = self.embedded_space.borrow().affine_dimension();
187        let extended_embedded_space = AffineSpace::new_affine(field, n + 1);
188
189        (
190            EmbeddedAffineSubspace {
191                ambient_space: extended_embedded_space,
192                embedded_space: self.embedded_space,
193                // 0, e_1, e_2, ..., e_(n-1)
194                embedding_points: {
195                    (0..n)
196                        .map(|k| {
197                            Vector::construct(extended_embedded_space, |i| {
198                                if k == 0 {
199                                    field.zero()
200                                } else {
201                                    let j = k - 1;
202                                    if i == j { field.one() } else { field.zero() }
203                                }
204                            })
205                        })
206                        .collect()
207                },
208            },
209            EmbeddedAffineSubspace {
210                ambient_space: self.ambient_space,
211                embedded_space: extended_embedded_space,
212                embedding_points: {
213                    let mut pts = self.embedding_points.clone();
214                    pts.push(pt);
215                    pts
216                },
217            },
218            Vector::construct(extended_embedded_space, |i| {
219                if i + 1 == n {
220                    self.field().one()
221                } else {
222                    self.field().zero()
223                }
224            }),
225        )
226    }
227
228    pub fn get_root_and_span(&self) -> Option<(Vector<'f, FS>, Vec<Vector<'f, FS>>)> {
229        let mut points = self.embedding_points.iter();
230        let root = points.next()?;
231        let span = points.map(|pt| pt - root).collect::<Vec<_>>();
232        Some((root.clone(), span))
233    }
234
235    pub fn get_embedding_points(&self) -> &Vec<Vector<'f, FS>> {
236        &self.embedding_points
237    }
238
239    pub fn embed_point(&self, pt: &Vector<'f, FS>) -> Vector<'f, FS> {
240        assert_eq!(pt.ambient_space(), self.embedded_space);
241        let (root, span) = self.get_root_and_span().unwrap(); //pt exists in the embedded space, so the embedded space is non-empty, so has a root and span
242        let mut total = root.clone();
243        for (i, vec) in span.iter().enumerate() {
244            total += &vec.scalar_mul(pt.coordinate(i));
245        }
246        total
247    }
248
249    pub fn embed_simplex(&self, spx: &Simplex<'f, FS>) -> Simplex<'f, FS> {
250        self.ambient_space()
251            .simplex(spx.points().iter().map(|p| self.embed_point(p)).collect())
252            .unwrap()
253    }
254
255    pub fn unembed_point(&self, pt: &Vector<'f, FS>) -> Option<Vector<'f, FS>> {
256        assert_eq!(pt.ambient_space(), self.ambient_space);
257        match self.get_root_and_span() {
258            Some((root, span)) => {
259                //solve root + x * basis = v for x
260                let y = (pt - &root).into_coordinates();
261                let basis_matrix = self
262                    .ambient_space
263                    .borrow()
264                    .cols_from_vectors(span.iter().collect());
265                let x = MatrixStructure::new(self.ambient_space.borrow().field().clone())
266                    .col_solve(basis_matrix, &y);
267                Some(self.embedded_space().vector(x?))
268            }
269            None => None,
270        }
271    }
272
273    pub fn unembed_simplex(&self, spx: &Simplex<'f, FS>) -> Option<Simplex<'f, FS>> {
274        let mut pts = vec![];
275        for embedded_pt in spx.points() {
276            match self.unembed_point(embedded_pt) {
277                Some(pt) => {
278                    pts.push(pt);
279                }
280                None => {
281                    return None;
282                }
283            }
284        }
285        Some(self.embedded_space().simplex(pts).unwrap())
286    }
287
288    pub fn to_oriented_hyperplane_intersection(&self) -> Option<Vec<OrientedHyperplane<'f, FS>>> {
289        let ambient_space = self.ambient_space();
290        let field = ambient_space.field();
291        match self.get_root_and_span() {
292            Some((root, span)) => {
293                let dim_amb = ambient_space.linear_dimension().unwrap();
294                let dim_ss = span.len();
295                // step 1: extend span to a basis by appending a subset of the elementary basis vectors
296                //first columns = span, remaining columns = identity matrix
297                let mat = Matrix::construct(dim_amb, dim_ss + dim_amb, |r, c| {
298                    if c < dim_ss {
299                        span[c].coordinate(r).clone()
300                    } else {
301                        let c = c - dim_ss;
302                        if r == c { field.one() } else { field.zero() }
303                    }
304                });
305                let (_, _, _, pivs) =
306                    MatrixStructure::new(field.clone()).row_hermite_algorithm(mat);
307                debug_assert_eq!(pivs.len(), dim_amb);
308                #[allow(clippy::needless_range_loop)]
309                for i in 0..dim_ss {
310                    debug_assert_eq!(pivs[i], i); //span is linearly independent so we expect this
311                }
312                let extension_elementary_basis_vectors = (dim_ss..dim_amb)
313                    .map(|i| pivs[i] - dim_ss)
314                    .collect::<Vec<_>>();
315
316                // step 2: take the hyperplanes formed by removing exactly one of each of the added elementary basis vectors at a time
317                let hyperplanes = (0..extension_elementary_basis_vectors.len())
318                    .map(|i| {
319                        let ref_point = {
320                            //root + e_k
321                            let k = extension_elementary_basis_vectors[i];
322                            Vector::construct(ambient_space, |l| {
323                                field.add(
324                                    root.coordinate(l),
325                                    &if l == k { field.one() } else { field.zero() },
326                                )
327                            })
328                        };
329                        OrientedSimplex::new_with_positive_point(
330                            ambient_space,
331                            {
332                                let mut points = vec![root.clone()];
333                                for s in &span {
334                                    points.push(&root + s);
335                                }
336                                for (j, &k) in extension_elementary_basis_vectors.iter().enumerate()
337                                {
338                                    if i != j {
339                                        //push root + e_k
340                                        points.push(Vector::construct(ambient_space, |l| {
341                                            field.add(root.coordinate(l), &{
342                                                if l == k { field.one() } else { field.zero() }
343                                            })
344                                        }));
345                                    }
346                                }
347                                points
348                            },
349                            &ref_point,
350                        )
351                        .unwrap()
352                        .into_oriented_hyperplane()
353                    })
354                    .collect::<Vec<_>>();
355                debug_assert_eq!(hyperplanes.len(), dim_amb - dim_ss);
356                Some(hyperplanes)
357            }
358            None => None,
359        }
360    }
361
362    pub fn to_oriented_hyperplane(&self) -> Option<OrientedHyperplane<'f, FS>> {
363        let mut oriented_hyperplanes = self.to_oriented_hyperplane_intersection()?;
364        if oriented_hyperplanes.len() == 1 {
365            Some(oriented_hyperplanes.pop().unwrap())
366        } else {
367            None
368        }
369    }
370}
371
372pub fn compose_affine_embeddings<'f, FS: OrderedRingSignature + FieldSignature>(
373    _a_to_b: EmbeddedAffineSubspace<'f, FS>,
374    _b_to_c: EmbeddedAffineSubspace<'f, FS>,
375) -> EmbeddedAffineSubspace<'f, FS> {
376    todo!() // call b_to_c.embed on the defining points of a_to_b
377}
378
379#[cfg(test)]
380mod tests {
381    use algebraeon_nzq::Rational;
382
383    use super::*;
384
385    #[test]
386    fn make_affine_subspace() {
387        let space = AffineSpace::new_linear(Rational::structure_ref(), 3);
388        let v1 = space.vector([1, 1, 1]);
389        let v2 = space.vector([1, 0, 0]);
390        let v3 = space.vector([0, 1, 0]);
391        let s = EmbeddedAffineSubspace::new_root_and_linear_independent_span(
392            space,
393            &v1,
394            vec![&v2, &v3],
395        );
396        s.unwrap();
397
398        let space = AffineSpace::new_linear(Rational::structure_ref(), 3);
399        let v1 = space.vector([1, 1, 1]);
400        let v2 = space.vector([1, 2, 0]);
401        let v3 = space.vector([-2, -4, 0]);
402        let s = EmbeddedAffineSubspace::new_root_and_linear_independent_span(
403            space,
404            &v1,
405            vec![&v2, &v3],
406        );
407        assert!(s.is_err());
408    }
409
410    #[test]
411    fn affine_subspace_embed_and_unembed() {
412        //1d embedded in 2d
413        {
414            let plane = AffineSpace::new_linear(Rational::structure_ref(), 2);
415            //the line x + y = 2
416            let (line, _) = EmbeddedAffineSubspace::new_root_and_linear_independent_span(
417                plane,
418                &plane.vector([1, 1]),
419                vec![&plane.vector([1, -1])],
420            )
421            .unwrap();
422
423            assert_eq!(
424                line.embed_point(&line.embedded_space().vector([-3])),
425                plane.vector([-2, 4])
426            );
427
428            assert_eq!(
429                line.unembed_point(&plane.vector([-1, 3])),
430                Some(line.embedded_space().vector([-2]))
431            );
432
433            assert_eq!(line.unembed_point(&plane.vector([1, 2])), None);
434        }
435
436        //2d embedded in 3d
437        {
438            let space = AffineSpace::new_linear(Rational::structure_ref(), 3);
439            let (plane, _) = EmbeddedAffineSubspace::new_root_and_linear_independent_span(
440                space,
441                &space.vector([3, 1, 2]),
442                vec![&space.vector([4, 2, 1]), &space.vector([1, -1, 2])],
443            )
444            .unwrap();
445
446            assert_eq!(
447                plane.embed_point(&plane.embedded_space().vector([-3, 2])),
448                space.vector([-7, -7, 3])
449            );
450
451            assert_eq!(
452                plane.unembed_point(&space.vector([0, -2, 3],)),
453                Some(plane.embedded_space().vector([-1, 1]))
454            );
455
456            assert_eq!(plane.unembed_point(&space.vector([1, 2, 1],)), None);
457        }
458    }
459
460    #[test]
461    fn extend_by_point_embedding_composition() {
462        let space = AffineSpace::new_linear(Rational::structure_ref(), 4);
463        let v1 = space.vector([1, 2, 1, 1]);
464        let v2 = space.vector([1, -2, 2, 0]);
465        let v3 = space.vector([2, 1, 0, 2]);
466        let (h, _) = EmbeddedAffineSubspace::new_root_and_linear_independent_span(
467            space,
468            &v1,
469            vec![&v2, &v3],
470        )
471        .unwrap();
472        let v4 = space.vector([0, 3, -2, 1]);
473        let (f, g, v4_inv) = h.extend_dimension_by_point_unsafe(v4.clone());
474        assert_eq!(g.embed_point(&v4_inv), v4);
475
476        let x = h
477            .embedded_space()
478            .vector([Rational::from(5), Rational::from(7)]);
479        //check that g(f(x)) = h(x)
480        assert_eq!(g.embed_point(&f.embed_point(&x)), h.embed_point(&x));
481    }
482}
algebraeon_geometry/affine_subspace.rs

algebraeon_geometry/
affine_subspace.rs
