geosets_rs/sets/
vpolytope.rs

1#![allow(unused)]
2use crate::linalg_utils::{argmax, rank};
3use crate::qhull_wrapper::{convex_hull, convex_hull_vertices, qhull_volume};
4
5use super::*;
6use good_lp::{Expression, Solution, SolverModel, default_solver, variable, variables};
7use ndarray_rand::RandomExt;
8use ndarray_rand::rand_distr::{Exp1, Uniform};
9use plotly::common::Mode;
10use plotly::{Plot, Scatter};
11use thiserror::Error;
12
13#[derive(Clone, Debug)]
14pub struct VPolytope {
15    vertices: Array2<f64>,
16}
17
18#[derive(Error, Debug)]
19pub enum VPolytopeError {
20    #[error("Vertices must not be empty!")]
21    EmptyVertices,
22}
23
24impl VPolytope {
25    pub fn new(vertices: Array2<f64>) -> Result<VPolytope, VPolytopeError> {
26        if vertices.is_empty() {
27            return Err(VPolytopeError::EmptyVertices);
28        }
29        Ok(VPolytope { vertices })
30    }
31
32    pub fn from_random(dim: usize, n_vertices: usize) -> Result<VPolytope, VPolytopeError> {
33        let vertices = Array2::random((n_vertices, dim), Uniform::new(-1.0, 1.0));
34        VPolytope::new(vertices)
35    }
36
37    pub fn n_vertices(&self) -> usize {
38        self.vertices.nrows()
39    }
40
41    pub fn compact_(&mut self) -> Result<(), SetOperationError> {
42        self.vertices = convex_hull_vertices(&self.vertices)?;
43        Ok(())
44    }
45
46    pub fn compact(&self) -> Result<VPolytope, SetOperationError> {
47        let mut copy = self.clone();
48        copy.compact_()?;
49        Ok(copy)
50    }
51}
52
53impl GeoSet for VPolytope {
54    fn dim(&self) -> usize {
55        self.vertices.dim().1
56    }
57
58    fn empty(&self) -> Result<bool, SetOperationError> {
59        Ok(false)
60    }
61
62    fn from_unit_box(dim: usize) -> Self {
63        let vertices = Array2::from_shape_fn((1 << dim, dim), |(i, j)| {
64            if (i & (1 << j)) != 0 { 1.0 } else { -1.0 }
65        });
66        VPolytope::new(vertices).unwrap()
67    }
68
69    fn to_vertices(&self) -> Result<Array2<f64>, SetOperationError> {
70        Ok(self.compact()?.vertices)
71    }
72
73    fn center(&self) -> Result<Array1<f64>, SetOperationError> {
74        // Centroid. Chebyshev center requires halfspaces
75        let center = self.vertices.mean_axis(Axis(0)).unwrap();
76        Ok(center)
77    }
78
79    fn support_function(
80        &self,
81        direction: Array1<f64>,
82    ) -> Result<(Array1<f64>, f64), SetOperationError> {
83        self._check_operand_dim(direction.dim())?;
84
85        let dot_product = self.vertices.dot(&direction);
86        let max_index = argmax(&dot_product).unwrap();
87
88        let support_value = dot_product[max_index];
89        let support_vector = self.vertices.row(max_index).to_owned();
90
91        Ok((support_vector, support_value))
92    }
93
94    fn volume(&self) -> Result<f64, SetOperationError> {
95        if self.degenerate() {
96            return Ok(0.0);
97        }
98
99        let vertices = self.to_vertices()?;
100        let qh = convex_hull(&vertices, true)?;
101
102        Ok(qhull_volume(&qh, &vertices)?)
103    }
104
105    fn minkowski_sum_(&mut self, other: &Self) -> Result<(), SetOperationError> {
106        let mut vertices = Array2::zeros((self.n_vertices() * other.n_vertices(), self.dim()));
107
108        for (i, row_self) in self.vertices.outer_iter().enumerate() {
109            for (j, row_other) in other.vertices.outer_iter().enumerate() {
110                vertices
111                    .row_mut(i * other.n_vertices() + j)
112                    .assign(&(&row_self + &row_other));
113            }
114        }
115
116        self.vertices = convex_hull_vertices(&vertices)?;
117        Ok(())
118    }
119
120    fn matmul_(&mut self, mat: &Array2<f64>) -> Result<(), SetOperationError> {
121        self._check_operand_dim(mat.dim().0)?;
122        self.vertices = self.vertices.dot(&mat.t());
123        Ok(())
124    }
125
126    fn translate_(&mut self, vector: &Array1<f64>) -> Result<(), SetOperationError> {
127        self._check_operand_dim(vector.dim())?;
128        // Translate each vertex by the vector
129        self.vertices = &self.vertices + &vector.view().insert_axis(Axis(0));
130        Ok(())
131    }
132
133    fn degenerate(&self) -> bool {
134        if self.n_vertices() == 1 {
135            return true;
136        }
137        let mat = &self.vertices - self.vertices.mean_axis(Axis(0)).unwrap();
138        rank(&mat).unwrap() < self.dim()
139    }
140
141    /// Evaluates the feasibility of the optimization problem
142    /// $\min 0$ \
143    /// $\text{subject to } V \lambda = p; 1^\top \lambda = 1, \lambda \geq 0 b$ \
144    fn contains_point(&self, point: &Array1<f64>) -> Result<bool, SetOperationError> {
145        let mut vars = variables!();
146        let lambda: Vec<_> = (0..self.n_vertices())
147            .map(|_| vars.add(variable().min(0.0))) // \lambda \geq 0
148            .collect();
149
150        // Build the problem with dummy objective
151        let mut problem = vars.minimise(0.0).using(default_solver);
152
153        // Add constraint V \lambda = p
154        for i in 0..self.dim() {
155            let v = self.vertices.column(i);
156            let expr: Expression = v.iter().zip(&lambda).map(|(vi, li)| *vi * *li).sum();
157            problem = problem.with(expr.eq(point[i]));
158        }
159
160        // 1^\top \lambda = 1
161        let expr: Expression = lambda.iter().copied().sum();
162        problem = problem.with(expr.eq(1.0));
163
164        // Try solving
165        match problem.solve() {
166            Ok(_) => Ok(true),   // feasible → not empty
167            Err(_) => Ok(false), // infeasible → empty
168        }
169    }
170}
171
172#[cfg(test)]
173mod tests {
174    use super::*;
175
176    #[test]
177    fn test_polytope_new() {
178        let _ = VPolytope::new(Array::ones((2, 2))).unwrap();
179        let _ = VPolytope::new(Array::ones((2, 5))).unwrap();
180    }
181}
geosets_rs/sets/vpolytope.rs

geosets_rs/sets/
vpolytope.rs
