algebraeon_geometry/
affine_subspace.rs

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

algebraeon_geometry/
affine_subspace.rs
