algebraeon_geometry/simplexes/
simplex.rs

1use itertools::Itertools;
2
3use super::*;
4
5#[derive(Clone)]
6pub struct Simplex<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> {
7    ambient_space: SP,
8    // points must be ordered w.r.t the ordering on vectors
9    // points must be non-degenerate
10    // points must belong to the ambient_space
11    points: Vec<Vector<FS, SP>>,
12}
13
14impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> std::fmt::Debug
15    for Simplex<FS, SP>
16{
17    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
18        f.debug_struct("Simplex")
19            .field("points", &self.points)
20            .finish()
21    }
22}
23
24impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> PartialEq
25    for Simplex<FS, SP>
26{
27    fn eq(&self, other: &Self) -> bool {
28        self.ambient_space.borrow() == other.ambient_space.borrow() && self.points == other.points
29    }
30}
31
32impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> Eq
33    for Simplex<FS, SP>
34{
35}
36
37impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> Hash
38    for Simplex<FS, SP>
39where
40    FS::Set: Hash,
41{
42    fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
43        self.points.hash(state);
44    }
45}
46
47impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> Simplex<FS, SP>
48where
49    SP: Clone,
50{
51    pub fn new(ambient_space: SP, mut points: Vec<Vector<FS, SP>>) -> Result<Self, &'static str> {
52        for point in &points {
53            assert_eq!(ambient_space.borrow(), point.ambient_space().borrow());
54        }
55        points.sort_unstable();
56        if !ambient_space
57            .borrow()
58            .are_points_affine_independent(points.iter().collect())
59        {
60            Err("Can't make a simplex using degenerate points")
61        } else {
62            Ok(Self {
63                ambient_space,
64                points,
65            })
66        }
67    }
68
69    pub fn ambient_space(&self) -> SP {
70        self.ambient_space.clone()
71    }
72
73    pub fn n(&self) -> usize {
74        self.points.len()
75    }
76
77    pub fn points(&self) -> &Vec<Vector<FS, SP>> {
78        &self.points
79    }
80
81    pub fn into_points(self) -> Vec<Vector<FS, SP>> {
82        self.points
83    }
84
85    pub fn point(&self, i: usize) -> &Vector<FS, SP> {
86        &self.points[i]
87    }
88
89    pub fn skeleton(&self, skel_n: isize) -> Vec<Self> {
90        if skel_n < 0 {
91            vec![]
92        } else {
93            let skel_n = skel_n as usize;
94            let mut parts = vec![];
95            for subset in (0..self.points.len()).combinations(skel_n) {
96                let part = Self::new(
97                    self.ambient_space.clone(),
98                    subset.into_iter().map(|i| self.points[i].clone()).collect(),
99                )
100                .unwrap();
101                parts.push(part);
102            }
103            parts
104        }
105    }
106
107    pub fn vertices(&self) -> Vec<Self> {
108        self.skeleton(1)
109    }
110
111    pub fn edges(&self) -> Vec<Self> {
112        self.skeleton(2)
113    }
114
115    pub fn faces(&self) -> Vec<Self> {
116        self.skeleton(3)
117    }
118
119    pub fn ridges(&self) -> Vec<Self> {
120        self.skeleton(self.points.len() as isize - 2)
121    }
122
123    pub fn facets(&self) -> Vec<Self> {
124        self.skeleton(self.points.len() as isize - 1)
125    }
126
127    pub fn facet(&self, k: usize) -> Self {
128        assert!(k <= self.points.len());
129        let facet = Self::new(self.ambient_space.clone(), {
130            let mut facet_points = self.points.clone();
131            facet_points.remove(k);
132            facet_points
133        })
134        .unwrap();
135        facet
136    }
137
138    pub fn sub_simplices(&self) -> Vec<Self> {
139        self.points()
140            .clone()
141            .into_iter()
142            .powerset()
143            .map(|sub_points| Self::new(self.ambient_space.clone(), sub_points).unwrap())
144            .collect()
145    }
146
147    pub fn sub_simplices_not_null(&self) -> Vec<Self> {
148        self.sub_simplices()
149            .into_iter()
150            .filter(|spx| spx.n() != 0)
151            .collect()
152    }
153
154    pub fn proper_sub_simplices_not_null(&self) -> Vec<Self> {
155        self.sub_simplices()
156            .into_iter()
157            .filter(|spx| spx.n() != 0 && spx.n() != self.n())
158            .collect()
159    }
160
161    pub fn sub_simplex(&self, pts: Vec<usize>) -> Self {
162        Self::new(
163            self.ambient_space(),
164            pts.iter().map(|idx| self.points[*idx].clone()).collect(),
165        )
166        .unwrap()
167    }
168
169    pub fn oriented_facet(&self, k: usize) -> OrientedSimplex<FS, SP> {
170        //return the oriented facet of self with negative side on the outside and positive side on the inside
171        assert!(k <= self.points.len());
172        let mut facet_points = self.points.clone();
173        let other_pt = facet_points.remove(k);
174        OrientedSimplex::new_with_positive_point(
175            self.ambient_space.clone(),
176            facet_points,
177            &other_pt,
178        )
179        .unwrap()
180    }
181
182    pub fn oriented_facets(&self) -> Vec<OrientedSimplex<FS, SP>> {
183        assert_eq!(self.ambient_space.borrow().affine_dimension(), self.n());
184        (0..self.n()).map(|k| self.oriented_facet(k)).collect()
185    }
186}
187
188#[derive(Debug, Clone, Copy, PartialEq, Eq)]
189pub enum OrientationSide {
190    Positive,
191    Neutral,
192    Negative,
193}
194
195/// a simplex spanning a hyperplane with a positive side and a negative side
196/// in 3d space it is a triangle
197/// in 2d space it is a line
198/// in 1d space it is a point
199/// in 0d space it is a null simplex. The orientation looses meaning but it is nice to still count this case.
200#[derive(Clone)]
201struct OrientedSimplexOrientation<
202    FS: OrderedRingStructure + FieldStructure,
203    SP: Borrow<AffineSpace<FS>> + Clone,
204> {
205    flip: bool,                     // flip the orientation if necessary
206    positive_point: Vector<FS, SP>, // a point on the positive side
207}
208
209impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> std::fmt::Debug
210    for OrientedSimplexOrientation<FS, SP>
211{
212    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
213        f.debug_struct("OrientedSimplexOrientation")
214            .field("flip", &self.flip)
215            .field("positive_point", &self.positive_point)
216            .finish()
217    }
218}
219
220#[derive(Clone)]
221pub struct OrientedSimplex<
222    FS: OrderedRingStructure + FieldStructure,
223    SP: Borrow<AffineSpace<FS>> + Clone,
224> {
225    simplex: Simplex<FS, SP>,
226    orientation: Option<OrientedSimplexOrientation<FS, SP>>, //None iff simplex is null
227}
228
229impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> std::fmt::Debug
230    for OrientedSimplex<FS, SP>
231{
232    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
233        f.debug_struct("OrientedSimplex")
234            .field("simplex", &self.simplex)
235            .field("orientation", &self.orientation)
236            .finish()
237    }
238}
239
240impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone>
241    OrientedSimplex<FS, SP>
242{
243    pub fn new_with_positive_point(
244        ambient_space: SP,
245        points: Vec<Vector<FS, SP>>,
246        ref_point: &Vector<FS, SP>,
247    ) -> Result<Self, &'static str> {
248        assert_eq!(ref_point.ambient_space().borrow(), ambient_space.borrow());
249        if points.len() != ambient_space.borrow().linear_dimension().unwrap() {
250            return Err("Oriented simplex must have dimension one less than the ambient space");
251        }
252        let n = points.len();
253        let simplex = Simplex::new(ambient_space, points)?;
254        if n == 0 {
255            Ok(Self {
256                simplex,
257                orientation: None,
258            })
259        } else {
260            let mut guess = Self {
261                simplex,
262                orientation: Some(OrientedSimplexOrientation {
263                    flip: false,
264                    positive_point: ref_point.clone(),
265                }),
266            };
267            match guess.classify_point(ref_point) {
268                OrientationSide::Positive => Ok(guess),
269                OrientationSide::Neutral => Err("The reference point lies on the hyperplane"),
270                OrientationSide::Negative => {
271                    let orientation = guess.orientation.as_mut().unwrap();
272                    orientation.flip = !orientation.flip;
273                    Ok(guess)
274                }
275            }
276        }
277    }
278
279    pub fn new_with_negative_point(
280        ambient_space: SP,
281        points: Vec<Vector<FS, SP>>,
282        ref_point: &Vector<FS, SP>,
283    ) -> Result<Self, &'static str> {
284        let mut ans = Self::new_with_positive_point(ambient_space, points, ref_point)?;
285        ans.flip();
286        Ok(ans)
287    }
288
289    pub fn positive_point(&self) -> Option<Vector<FS, SP>> {
290        Some(self.orientation.as_ref()?.positive_point.clone())
291    }
292
293    pub fn negative_point(&self) -> Option<Vector<FS, SP>> {
294        let positive_point = self.positive_point()?;
295        Some({
296            let pt: &Vector<FS, SP> = &self.simplex.points()[0];
297            &(pt + pt) - &positive_point
298        })
299    }
300
301    pub fn ambient_space(&self) -> SP {
302        self.simplex().ambient_space()
303    }
304
305    pub fn simplex(&self) -> &Simplex<FS, SP> {
306        &self.simplex
307    }
308
309    pub fn flip(&mut self) {
310        let negative_point = self.negative_point();
311        match &mut self.orientation {
312            Some(OrientedSimplexOrientation {
313                flip,
314                positive_point,
315            }) => {
316                (*flip, *positive_point) = (!*flip, negative_point.unwrap());
317            }
318            None => {}
319        }
320        match &self.orientation {
321            Some(OrientedSimplexOrientation {
322                flip: _flip,
323                positive_point,
324            }) => {
325                debug_assert_eq!(
326                    self.classify_point(positive_point),
327                    OrientationSide::Positive
328                );
329            }
330            None => {}
331        }
332    }
333
334    fn classify_point_quantitatively(&self, point: &Vector<FS, SP>) -> FS::Set {
335        let space = self.ambient_space();
336        let field = space.borrow().ordered_field();
337        match &self.orientation {
338            Some(OrientedSimplexOrientation {
339                flip,
340                positive_point: _,
341            }) => match space.borrow().linear_dimension().unwrap() {
342                0 => unreachable!(),
343                d => {
344                    let root = &self.simplex.points[0];
345                    let mut vecs = (1..d)
346                        .map(|i| &self.simplex.points[i] - root)
347                        .collect::<Vec<_>>();
348                    vecs.push(point - root);
349                    let det = space.borrow().determinant(vecs.iter().collect());
350                    match flip {
351                        true => field.neg(&det),
352                        false => det,
353                    }
354                }
355            },
356            None => field.zero(),
357        }
358    }
359
360    pub fn classify_point(&self, point: &Vector<FS, SP>) -> OrientationSide {
361        let space = self.ambient_space();
362        let field = space.borrow().ordered_field();
363        let value = self.classify_point_quantitatively(point);
364        match field.ring_cmp(&value, &field.zero()) {
365            std::cmp::Ordering::Less => OrientationSide::Negative,
366            std::cmp::Ordering::Equal => OrientationSide::Neutral,
367            std::cmp::Ordering::Greater => OrientationSide::Positive,
368        }
369    }
370
371    pub fn into_oriented_hyperplane(self) -> OrientedHyperplane<FS, SP> {
372        OrientedHyperplane {
373            oriented_simplex: self,
374        }
375    }
376}
377
378#[derive(Clone)]
379pub struct OrientedHyperplane<
380    FS: OrderedRingStructure + FieldStructure,
381    SP: Borrow<AffineSpace<FS>> + Clone,
382> {
383    oriented_simplex: OrientedSimplex<FS, SP>,
384}
385
386impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone> std::fmt::Debug
387    for OrientedHyperplane<FS, SP>
388{
389    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
390        f.debug_struct("OrientedHyperplane")
391            .field("oriented_simplex", &self.oriented_simplex)
392            .finish()
393    }
394}
395
396// impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone>
397//     From<OrientedSimplex<FS, SP>> for OrientedHyperplane<FS, SP>
398// {
399//     fn from(oriented_simplex: OrientedSimplex<FS, SP>) -> Self {
400//         Self { oriented_simplex }
401//     }
402// }
403
404pub enum OrientedHyperplaneIntersectLineSegmentResult<
405    FS: OrderedRingStructure + FieldStructure,
406    SP: Borrow<AffineSpace<FS>> + Clone,
407> {
408    PositivePositive,
409    NegativeNegative,
410    PositiveNeutral,
411    NegativeNeutral,
412    NeutralPositive,
413    NeutralNegative,
414    NeutralNeutral,
415    PositiveNegative { intersection_point: Vector<FS, SP> },
416    NegativePositive { intersection_point: Vector<FS, SP> },
417}
418
419impl<FS: OrderedRingStructure + FieldStructure, SP: Borrow<AffineSpace<FS>> + Clone>
420    OrientedHyperplane<FS, SP>
421{
422    pub fn ambient_space(&self) -> SP {
423        self.oriented_simplex.ambient_space()
424    }
425
426    pub fn classify_point(&self, point: &Vector<FS, SP>) -> OrientationSide {
427        self.oriented_simplex.classify_point(point)
428    }
429
430    pub fn intersect_line(
431        &self,
432        a: &Vector<FS, SP>,
433        b: &Vector<FS, SP>,
434    ) -> OrientedHyperplaneIntersectLineSegmentResult<FS, SP> {
435        let space = self.ambient_space();
436        let field = space.borrow().ordered_field();
437
438        let a_val = self.oriented_simplex.classify_point_quantitatively(a);
439        let b_val = self.oriented_simplex.classify_point_quantitatively(b);
440
441        match (
442            field.ring_cmp(&a_val, &field.zero()),
443            field.ring_cmp(&b_val, &field.zero()),
444        ) {
445            (std::cmp::Ordering::Less, std::cmp::Ordering::Greater) => {
446                let t = field
447                    .div(&a_val, &field.add(&a_val, &field.neg(&b_val)))
448                    .unwrap();
449                {
450                    debug_assert_eq!(
451                        self.oriented_simplex.classify_point(a),
452                        OrientationSide::Negative
453                    );
454                    debug_assert_eq!(
455                        self.oriented_simplex.classify_point(b),
456                        OrientationSide::Positive
457                    );
458                    OrientedHyperplaneIntersectLineSegmentResult::NegativePositive {
459                        intersection_point: a + &(b - a).scalar_mul(&t),
460                    }
461                }
462            }
463            (std::cmp::Ordering::Greater, std::cmp::Ordering::Less) => {
464                let t = field
465                    .div(&a_val, &field.add(&a_val, &field.neg(&b_val)))
466                    .unwrap();
467                {
468                    debug_assert_eq!(
469                        self.oriented_simplex.classify_point(a),
470                        OrientationSide::Positive
471                    );
472                    debug_assert_eq!(
473                        self.oriented_simplex.classify_point(b),
474                        OrientationSide::Negative
475                    );
476                    OrientedHyperplaneIntersectLineSegmentResult::PositiveNegative {
477                        intersection_point: a + &(b - a).scalar_mul(&t),
478                    }
479                }
480            }
481            (std::cmp::Ordering::Less, std::cmp::Ordering::Less) => {
482                OrientedHyperplaneIntersectLineSegmentResult::NegativeNegative
483            }
484            (std::cmp::Ordering::Less, std::cmp::Ordering::Equal) => {
485                OrientedHyperplaneIntersectLineSegmentResult::NegativeNeutral
486            }
487            (std::cmp::Ordering::Equal, std::cmp::Ordering::Less) => {
488                OrientedHyperplaneIntersectLineSegmentResult::NeutralNegative
489            }
490            (std::cmp::Ordering::Equal, std::cmp::Ordering::Equal) => {
491                OrientedHyperplaneIntersectLineSegmentResult::NeutralNeutral
492            }
493            (std::cmp::Ordering::Equal, std::cmp::Ordering::Greater) => {
494                OrientedHyperplaneIntersectLineSegmentResult::NeutralPositive
495            }
496            (std::cmp::Ordering::Greater, std::cmp::Ordering::Equal) => {
497                OrientedHyperplaneIntersectLineSegmentResult::PositiveNeutral
498            }
499            (std::cmp::Ordering::Greater, std::cmp::Ordering::Greater) => {
500                OrientedHyperplaneIntersectLineSegmentResult::PositivePositive
501            }
502        }
503    }
504
505    pub fn flip(&mut self) {
506        self.oriented_simplex.flip()
507    }
508}
509
510/*
511pub fn simplex_intersect_with_hyperplane<
512    FS: OrderedRingStructure + FieldStructure,
513    SP: Borrow<AffineSpace<FS>> + Clone,
514>(
515    simplex: &Simplex<FS, SP>,
516    hyperplane: &OrientedSimplex<FS, SP>,
517) -> Vec<Simplex<FS, SP>> {
518    todo!()
519}
520
521pub fn simplex_intersect_positive_side_hyperplane<
522    FS: OrderedRingStructure + FieldStructure,
523    SP: Borrow<AffineSpace<FS>> + Clone,
524>(
525    simplex: &Simplex<FS, SP>,
526    hyperplane: &OrientedSimplex<FS, SP>,
527) -> Vec<Simplex<FS, SP>> {
528    todo!()
529}
530
531pub fn simplex_intersect_negative_side_hyperplane<
532    FS: OrderedRingStructure + FieldStructure,
533    SP: Borrow<AffineSpace<FS>> + Clone,
534>(
535    simplex: &Simplex<FS, SP>,
536    hyperplane: &OrientedSimplex<FS, SP>,
537) -> Vec<Simplex<FS, SP>> {
538    simplex_intersect_positive_side_hyperplane(simplex, &hyperplane.clone().flip())
539}
540*/
541
542#[cfg(test)]
543mod tests {
544    use algebraeon_nzq::rational::*;
545    use algebraeon_sets::structure::*;
546
547    use super::*;
548
549    #[test]
550    fn make_simplex() {
551        let space = AffineSpace::new_linear(Rational::structure(), 2);
552        let v1 = Vector::new(&space, vec![Rational::from(1), Rational::from(1)]);
553        let v2 = Vector::new(&space, vec![Rational::from(1), Rational::from(0)]);
554        let v3 = Vector::new(&space, vec![Rational::from(0), Rational::from(1)]);
555        let s = Simplex::new(&space, vec![v1, v2, v3]);
556        assert!(s.is_ok());
557
558        let space = AffineSpace::new_linear(Rational::structure(), 2);
559        let v1 = Vector::new(&space, vec![Rational::from(0), Rational::from(0)]);
560        let v2 = Vector::new(&space, vec![Rational::from(1), Rational::from(0)]);
561        let v3 = Vector::new(&space, vec![Rational::from(2), Rational::from(0)]);
562        let s = Simplex::new(&space, vec![v1, v2, v3]);
563        assert!(s.is_err());
564    }
565
566    #[test]
567    fn simplex_skeleton() {
568        let space = AffineSpace::new_linear(Rational::structure(), 2);
569        let v1 = Vector::new(&space, vec![Rational::from(1), Rational::from(1)]);
570        let v2 = Vector::new(&space, vec![Rational::from(1), Rational::from(0)]);
571        let v3 = Vector::new(&space, vec![Rational::from(0), Rational::from(1)]);
572        let s = Simplex::new(&space, vec![v1, v2, v3]).unwrap();
573
574        assert_eq!(s.skeleton(-2).len(), 0);
575        assert_eq!(s.skeleton(-1).len(), 0);
576        assert_eq!(s.skeleton(0).len(), 1);
577        assert_eq!(s.vertices().len(), 3);
578        assert_eq!(s.edges().len(), 3);
579        assert_eq!(s.faces().len(), 1);
580        assert_eq!(s.skeleton(4).len(), 0);
581    }
582
583    #[test]
584    fn make_oriented_simplex() {
585        let space = AffineSpace::new_linear(Rational::structure(), 2);
586        let v1 = Vector::new(&space, vec![Rational::from(1), Rational::from(0)]);
587        let v2 = Vector::new(&space, vec![Rational::from(0), Rational::from(1)]);
588        let v3 = Vector::new(&space, vec![Rational::from(2), Rational::from(3)]);
589        let s_pos =
590            OrientedSimplex::new_with_positive_point(&space, vec![v1.clone(), v2.clone()], &v3)
591                .unwrap();
592        let s_neg = OrientedSimplex::new_with_negative_point(&space, vec![v1, v2], &v3).unwrap();
593
594        assert_ne!(
595            s_pos.orientation.unwrap().flip,
596            s_neg.orientation.unwrap().flip
597        );
598    }
599}
algebraeon_geometry/simplexes/simplex.rs

algebraeon_geometry/simplexes/
simplex.rs
