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                        let ref_point = {
296                            //root + e_k
297                            let k = extension_elementary_basis_vectors[i];
298                            Vector::construct(ambient_space.clone(), |l| {
299                                ordered_field.add(
300                                    root.coordinate(l),
301                                    &match l == k {
302                                        false => ordered_field.zero(),
303                                        true => ordered_field.one(),
304                                    },
305                                )
306                            })
307                        };
308                        OrientedSimplex::new_with_positive_point(
309                            ambient_space.clone(),
310                            {
311                                let mut points = vec![root.clone()];
312                                for s in &span {
313                                    points.push(&root + s);
314                                }
315                                for j in 0..extension_elementary_basis_vectors.len() {
316                                    if i != j {
317                                        let k = extension_elementary_basis_vectors[j];
318                                        //push root + e_k
319                                        points.push(Vector::construct(ambient_space.clone(), |l| {
320                                            ordered_field.add(root.coordinate(l), &{
321                                                match l == k {
322                                                    false => ordered_field.zero(),
323                                                    true => ordered_field.one(),
324                                                }
325                                            })
326                                        }))
327                                    }
328                                }
329                                points
330                            },
331                            &ref_point,
332                        )
333                        .unwrap()
334                        .into_oriented_hyperplane()
335                    })
336                    .collect::<Vec<_>>();
337                debug_assert_eq!(hyperplanes.len(), dim_amb - dim_ss);
338                Some(hyperplanes)
339            }
340            None => None,
341        }
342    }
343}
344
345pub fn compose_affine_embeddings<
346    FS: OrderedRingStructure + FieldStructure,
347    SPA: Borrow<AffineSpace<FS>> + Clone,
348    SPB: Borrow<AffineSpace<FS>> + Clone,
349    SPC: Borrow<AffineSpace<FS>> + Clone,
350>(
351    _a_to_b: EmbeddedAffineSubspace<FS, SPB, SPA>,
352    _b_to_c: EmbeddedAffineSubspace<FS, SPC, SPB>,
353) -> EmbeddedAffineSubspace<FS, SPC, SPA> {
354    todo!() // call b_to_c.embed on the defining points of a_to_b
355}
356
357#[cfg(test)]
358mod tests {
359    use algebraeon_nzq::rational::*;
360    use algebraeon_sets::structure::*;
361
362    use super::*;
363
364    #[test]
365    fn make_affine_subspace() {
366        let space = AffineSpace::new_linear(Rational::structure(), 3);
367        let v1 = Vector::new(
368            &space,
369            vec![Rational::from(1), Rational::from(1), Rational::from(1)],
370        );
371        let v2 = Vector::new(
372            &space,
373            vec![Rational::from(1), Rational::from(0), Rational::from(0)],
374        );
375        let v3 = Vector::new(
376            &space,
377            vec![Rational::from(0), Rational::from(1), Rational::from(0)],
378        );
379        let s = EmbeddedAffineSubspace::new(&space, v1, vec![v2, v3]);
380        s.unwrap();
381
382        let space = AffineSpace::new_linear(Rational::structure(), 3);
383        let v1 = Vector::new(
384            &space,
385            vec![Rational::from(1), Rational::from(1), Rational::from(1)],
386        );
387        let v2 = Vector::new(
388            &space,
389            vec![Rational::from(1), Rational::from(2), Rational::from(0)],
390        );
391        let v3 = Vector::new(
392            &space,
393            vec![Rational::from(-2), Rational::from(-4), Rational::from(0)],
394        );
395        let s = EmbeddedAffineSubspace::new(&space, v1, vec![v2, v3]);
396        assert!(s.is_err());
397    }
398
399    #[test]
400    fn affine_subspace_embed_and_unembed() {
401        //1d embedded in 2d
402        {
403            let plane = AffineSpace::new_linear(Rational::structure(), 2);
404            //the line x + y = 2
405            let (line, _) = EmbeddedAffineSubspace::new(
406                &plane,
407                Vector::new(&plane, vec![Rational::from(1), Rational::from(1)]),
408                vec![Vector::new(
409                    &plane,
410                    vec![Rational::from(1), Rational::from(-1)],
411                )],
412            )
413            .unwrap();
414
415            assert_eq!(
416                line.embed_point(&Vector::new(
417                    line.embedded_space(),
418                    vec![Rational::from(-3)],
419                )),
420                Vector::new(&plane, vec![Rational::from(-2), Rational::from(4)])
421            );
422
423            assert_eq!(
424                line.unembed_point(&Vector::new(
425                    &plane,
426                    vec![Rational::from(-1), Rational::from(3)],
427                )),
428                Some(Vector::new(line.embedded_space(), vec![Rational::from(-2)],))
429            );
430
431            assert_eq!(
432                line.unembed_point(&Vector::new(
433                    &plane,
434                    vec![Rational::from(1), Rational::from(2)],
435                )),
436                None
437            );
438        }
439
440        //2d embedded in 3d
441        {
442            let space = AffineSpace::new_linear(Rational::structure(), 3);
443            let (plane, _) = EmbeddedAffineSubspace::new(
444                &space,
445                Vector::new(
446                    &space,
447                    vec![Rational::from(3), Rational::from(1), Rational::from(2)],
448                ),
449                vec![
450                    Vector::new(
451                        &space,
452                        vec![Rational::from(4), Rational::from(2), Rational::from(1)],
453                    ),
454                    Vector::new(
455                        &space,
456                        vec![Rational::from(1), Rational::from(-1), Rational::from(2)],
457                    ),
458                ],
459            )
460            .unwrap();
461
462            assert_eq!(
463                plane.embed_point(&Vector::new(
464                    plane.embedded_space(),
465                    vec![Rational::from(-3), Rational::from(2)],
466                )),
467                Vector::new(
468                    &space,
469                    vec![Rational::from(-7), Rational::from(-7), Rational::from(3)]
470                )
471            );
472
473            assert_eq!(
474                plane.unembed_point(&Vector::new(
475                    &space,
476                    vec![Rational::from(0), Rational::from(-2), Rational::from(3)],
477                )),
478                Some(Vector::new(
479                    plane.embedded_space(),
480                    vec![Rational::from(-1), Rational::from(1)],
481                ))
482            );
483
484            assert_eq!(
485                plane.unembed_point(&Vector::new(
486                    &space,
487                    vec![Rational::from(1), Rational::from(2), Rational::from(2)],
488                )),
489                None
490            );
491        }
492    }
493
494    #[test]
495    fn extend_by_point_embedding_composition() {
496        let space = AffineSpace::new_linear(Rational::structure(), 4);
497        let v1 = Vector::new(
498            &space,
499            vec![
500                Rational::from(1),
501                Rational::from(2),
502                Rational::from(1),
503                Rational::from(1),
504            ],
505        );
506        let v2 = Vector::new(
507            &space,
508            vec![
509                Rational::from(1),
510                Rational::from(-2),
511                Rational::from(2),
512                Rational::from(0),
513            ],
514        );
515        let v3 = Vector::new(
516            &space,
517            vec![
518                Rational::from(2),
519                Rational::from(1),
520                Rational::from(0),
521                Rational::from(2),
522            ],
523        );
524        let (h, _) = EmbeddedAffineSubspace::new(&space, v1, vec![v2, v3]).unwrap();
525        let v4 = Vector::new(
526            &space,
527            vec![
528                Rational::from(0),
529                Rational::from(3),
530                Rational::from(-2),
531                Rational::from(1),
532            ],
533        );
534        let (f, g, v4_inv) = h.extend_dimension_by_point_unsafe(v4.clone());
535        assert_eq!(g.embed_point(&v4_inv), v4);
536
537        let x = Vector::new(
538            h.embedded_space(),
539            vec![Rational::from(5), Rational::from(7)],
540        );
541        //check that g(f(x)) = h(x)
542        assert_eq!(g.embed_point(&f.embed_point(&x)), h.embed_point(&x));
543    }
544}
algebraeon_geometry/affine_subspace.rs

algebraeon_geometry/
affine_subspace.rs
