solvr 0.2.0

Advanced computing library for real-world problem solving - optimization, differential equations, interpolation, statistics, and more
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

//! Delaunay triangulation trait.
//!
//! Computes the Delaunay triangulation of a point set. The Delaunay triangulation
//! maximizes the minimum angle of all triangles and has important properties
//! for interpolation and mesh generation.
use crate::DType;

use numr::error::Result;
use numr::runtime::Runtime;
use numr::tensor::Tensor;

/// Delaunay triangulation result.
#[derive(Debug, Clone)]
pub struct Delaunay<R: Runtime<DType = DType>> {
    /// Original points `[n, d]`.
    pub points: Tensor<R>,

    /// Simplices of the triangulation `[n_simplices, d+1]` (I64 dtype).
    /// For 2D: triangles as triples of vertex indices.
    /// For 3D: tetrahedra as 4-tuples of vertex indices.
    pub simplices: Tensor<R>,

    /// Neighboring simplex indices `[n_simplices, d+1]` (I64 dtype).
    /// `neighbors[i, j]` is the simplex sharing the facet opposite vertex j.
    /// -1 indicates no neighbor (boundary).
    pub neighbors: Tensor<R>,

    /// For each point, index of one simplex containing it `[n]` (I64 dtype).
    /// Used for efficient point location queries.
    pub vertex_to_simplex: Tensor<R>,

    /// Vertex indices of the convex hull `[n_hull]` (I64 dtype).
    pub convex_hull: Tensor<R>,
}

/// Algorithmic contract for Delaunay triangulation operations.
///
/// All backends implementing Delaunay algorithms MUST implement this trait.
pub trait DelaunayAlgorithms<R: Runtime<DType = DType>> {
    /// Compute the Delaunay triangulation of a point set.
    ///
    /// # Arguments
    ///
    /// * `points` - Point set with shape (n, d) where n >= d+1
    ///
    /// # Returns
    ///
    /// Delaunay structure containing simplices and neighbor information.
    ///
    /// # Algorithm
    ///
    /// Uses the Bowyer-Watson incremental algorithm:
    /// 1. Create a super-triangle/tetrahedron containing all points
    /// 2. Insert points one by one
    /// 3. For each point, find all simplices whose circumsphere contains it
    /// 4. Remove those simplices and re-triangulate the cavity
    /// 5. Remove simplices connected to the super-simplex
    ///
    /// Complexity: O(n^(d/2+1)) worst case, O(n log n) for well-distributed points.
    ///
    /// # Examples
    ///
    /// ```
    /// # use numr::runtime::cpu::{CpuClient, CpuDevice};
    /// # use numr::tensor::Tensor;
    /// use solvr::spatial::DelaunayAlgorithms;
    /// # let device = CpuDevice::new();
    /// # let client = CpuClient::new(device.clone());
    /// // 2D points
    /// # let points = Tensor::from_slice(&[
    /// #     0.0, 0.0,
    /// #     1.0, 0.0,
    /// #     0.5, 1.0,
    /// #     0.5, 0.3,
    /// # ], &[4, 2], &device);
    /// let tri = client.delaunay(&points)?;
    /// # Ok::<(), Box<dyn std::error::Error>>(())
    /// ```
    fn delaunay(&self, points: &Tensor<R>) -> Result<Delaunay<R>>;

    /// Find the simplex containing each query point.
    ///
    /// # Arguments
    ///
    /// * `tri` - The Delaunay triangulation
    /// * `query` - Query points with shape (m, d)
    ///
    /// # Returns
    ///
    /// Tensor `[m]` (I64 dtype) with simplex indices. -1 for points outside the hull.
    fn delaunay_find_simplex(&self, tri: &Delaunay<R>, query: &Tensor<R>) -> Result<Tensor<R>>;

    /// Get the vertex neighbors of each point.
    ///
    /// # Arguments
    ///
    /// * `tri` - The Delaunay triangulation
    ///
    /// # Returns
    ///
    /// For each vertex, the indices of neighboring vertices.
    /// Returns (indices, indptr) in CSR-like format:
    /// - indices: concatenated neighbor indices
    /// - indptr: start/end positions for each vertex [n+1]
    fn delaunay_vertex_neighbors(&self, tri: &Delaunay<R>) -> Result<(Tensor<R>, Tensor<R>)>;
}

#[cfg(test)]
mod tests {
    // Tests will be in the implementation files
}