algebraeon_geometry/simplexes/
simplex.rs

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

algebraeon_geometry/simplexes/
simplex.rs
