geosets_rs/cddlib_rs/
mod.rs

1mod cdd_bindings;
2
3use crate::sets::errors::SetOperationError;
4use cdd_bindings::*;
5use ndarray::{Array1, Array2};
6use once_cell::sync::Lazy;
7use std::os::raw::c_long;
8use std::sync::Mutex;
9
10static CDD_MUTEX: Lazy<Mutex<()>> = Lazy::new(|| Mutex::new(()));
11
12/// Computes the vertices of a polytope defined by the inequality Ax ≤ b.
13///
14/// This function uses rust bindings of the cddlib library to convert the H-representation
15/// (halfspace representation) to V-representation (vertex representation).
16pub fn compute_polytope_vertices(
17    a: &Array2<f64>,
18    b: &Array1<f64>,
19) -> Result<Array2<f64>, SetOperationError> {
20    // Check dimension compatibility
21    if a.nrows() != b.len() {
22        return Err(SetOperationError::DimensionMismatch {
23            expected: a.nrows(),
24            got: b.len(),
25        });
26    }
27
28    // Due to the global_constants calls, we need to ensure that only one thread
29    // is executing this at a time.
30    let _guard = CDD_MUTEX.lock().unwrap();
31
32    unsafe {
33        dd_set_global_constants();
34        let result = compute_polytope_vertices_inner(a, b);
35        dd_free_global_constants();
36        result
37    }
38}
39
40unsafe fn compute_polytope_vertices_inner(
41    a: &Array2<f64>,
42    b: &Array1<f64>,
43) -> Result<Array2<f64>, SetOperationError> {
44    let m = a.nrows() as c_long;
45    let n = a.ncols() as c_long;
46    // Initialize cddlib
47    // Create matrix [b | -A]
48    unsafe {
49        let mat = dd_CreateMatrix(m, n + 1);
50        if mat.is_null() {
51            dd_free_global_constants();
52            return Err(SetOperationError::DataConversionError {
53                source: "Failed to create cddlib matrix".into(),
54            });
55        }
56
57        (*mat).representation = dd_RepresentationType::dd_Inequality;
58
59        for i in 0..m {
60            let i_usize = i as usize;
61            let i_isize = i as isize;
62
63            // RHS - get row pointer, then column pointer, then access the [f64; 1] array
64            let row_ptr = (*mat).matrix.offset(i_isize);
65            let rhs_ptr = (*row_ptr).offset(0);
66            (*rhs_ptr)[0] = b[i_usize];
67
68            // Coefficients (-A)
69            for j in 0..n {
70                let j_usize = j as usize;
71                let j_isize = (j + 1) as isize;
72                let coeff_ptr = (*row_ptr).offset(j_isize);
73                (*coeff_ptr)[0] = -a[[i_usize, j_usize]];
74            }
75        }
76
77        // Build polyhedron
78        let mut err: dd_ErrorType = dd_ErrorType::dd_NoError;
79        let poly = dd_DDMatrix2Poly(mat, &mut err);
80        if err != dd_ErrorType::dd_NoError || poly.is_null() {
81            dd_FreeMatrix(mat);
82            dd_free_global_constants();
83            return Err(SetOperationError::InfeasibleOptimization {
84                source: format!("cddlib error: {:?}", err).into(),
85            });
86        }
87
88        // Extract generators (vertices + rays)
89        let gens = dd_CopyGenerators(poly);
90        if gens.is_null() {
91            dd_FreePolyhedra(poly);
92            dd_FreeMatrix(mat);
93            dd_free_global_constants();
94            return Err(SetOperationError::DataConversionError {
95                source: "Failed to extract generators from polyhedron".into(),
96            });
97        }
98
99        let mut vertices_data = Vec::new();
100        let mut vertex_count = 0;
101
102        for i in 0..(*gens).rowsize {
103            let i_isize = i as isize;
104            let gen_row_ptr = (*gens).matrix.offset(i_isize);
105            let kind_ptr = (*gen_row_ptr).offset(0);
106            let kind = (*kind_ptr)[0];
107
108            if (kind - 1.0).abs() < 1e-9 {
109                // It's a vertex
110                for j in 1..=n {
111                    let j_isize = j as isize;
112                    let vertex_ptr = (*gen_row_ptr).offset(j_isize);
113                    vertices_data.push((*vertex_ptr)[0]);
114                }
115                vertex_count += 1;
116            }
117        }
118
119        // Cleanup
120        dd_FreePolyhedra(poly);
121        dd_FreeMatrix(mat);
122        dd_FreeMatrix(gens);
123        dd_free_global_constants();
124
125        // Convert to Array2
126        if vertex_count == 0 {
127            Ok(Array2::zeros((0, n as usize)))
128        } else {
129            Array2::from_shape_vec((vertex_count, n as usize), vertices_data).map_err(|e| {
130                SetOperationError::DataConversionError {
131                    source: format!("Failed to create Array2 from vertices data: {}", e).into(),
132                }
133            })
134        }
135    }
136}
geosets_rs/cddlib_rs/mod.rs

geosets_rs/cddlib_rs/
mod.rs
