algebraeon_geometry/
affine_subspace.rs

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

algebraeon_geometry/
affine_subspace.rs
