surge_sparse/
csc.rs

1// SPDX-License-Identifier: LicenseRef-PolyForm-Noncommercial-1.0.0
2//! Canonical validated CSC sparse matrix representation.
3
4use std::ops::AddAssign;
5
6use crate::error::{SparseError, SparseResult};
7
8/// Sparse triplet entry (row, col, value).
9#[derive(Debug, Clone, Copy, PartialEq)]
10pub struct Triplet<T> {
11    pub row: usize,
12    pub col: usize,
13    pub val: T,
14}
15
16/// Sparse matrix in validated compressed sparse column (CSC) form.
17#[derive(Debug, Clone, PartialEq)]
18pub struct CscMatrix<T> {
19    nrows: usize,
20    ncols: usize,
21    col_ptrs: Vec<usize>,
22    row_indices: Vec<usize>,
23    values: Vec<T>,
24}
25
26impl<T> CscMatrix<T> {
27    /// Build a validated CSC matrix from owned parts.
28    pub fn try_new(
29        nrows: usize,
30        ncols: usize,
31        col_ptrs: Vec<usize>,
32        row_indices: Vec<usize>,
33        values: Vec<T>,
34    ) -> SparseResult<Self> {
35        validate_csc_pattern(nrows, ncols, &col_ptrs, &row_indices)?;
36        if row_indices.len() != values.len() {
37            return Err(SparseError::ValueCountMismatch {
38                expected: row_indices.len(),
39                found: values.len(),
40            });
41        }
42        Ok(Self {
43            nrows,
44            ncols,
45            col_ptrs,
46            row_indices,
47            values,
48        })
49    }
50
51    pub fn nrows(&self) -> usize {
52        self.nrows
53    }
54
55    pub fn ncols(&self) -> usize {
56        self.ncols
57    }
58
59    pub fn nnz(&self) -> usize {
60        self.values.len()
61    }
62
63    pub fn is_square(&self) -> bool {
64        self.nrows == self.ncols
65    }
66
67    pub fn col_ptrs(&self) -> &[usize] {
68        &self.col_ptrs
69    }
70
71    pub fn row_indices(&self) -> &[usize] {
72        &self.row_indices
73    }
74
75    pub fn values(&self) -> &[T] {
76        &self.values
77    }
78
79    pub fn values_mut(&mut self) -> &mut [T] {
80        &mut self.values
81    }
82
83    pub fn into_parts(self) -> (usize, usize, Vec<usize>, Vec<usize>, Vec<T>) {
84        (
85            self.nrows,
86            self.ncols,
87            self.col_ptrs,
88            self.row_indices,
89            self.values,
90        )
91    }
92
93    pub fn has_same_pattern<U>(&self, other: &CscMatrix<U>) -> bool {
94        self.nrows == other.nrows
95            && self.ncols == other.ncols
96            && self.col_ptrs == other.col_ptrs
97            && self.row_indices == other.row_indices
98    }
99
100    pub fn has_same_pattern_slices(&self, col_ptrs: &[usize], row_indices: &[usize]) -> bool {
101        self.col_ptrs == col_ptrs && self.row_indices == row_indices
102    }
103}
104
105impl<T> CscMatrix<T>
106where
107    T: Clone,
108{
109    /// Convert CSC indices to `i32` with checked bounds.
110    pub fn try_to_i32(&self) -> SparseResult<(Vec<i32>, Vec<i32>, Vec<T>)> {
111        let col_ptrs = self
112            .col_ptrs
113            .iter()
114            .copied()
115            .map(|value| try_usize_to_i32("column pointer", value))
116            .collect::<SparseResult<Vec<_>>>()?;
117        let row_indices = self
118            .row_indices
119            .iter()
120            .copied()
121            .map(|value| try_usize_to_i32("row index", value))
122            .collect::<SparseResult<Vec<_>>>()?;
123        Ok((col_ptrs, row_indices, self.values.clone()))
124    }
125}
126
127impl<T> CscMatrix<T>
128where
129    T: Copy + AddAssign,
130{
131    /// Convert triplets to canonical CSC, summing duplicate entries.
132    pub fn try_from_triplets(
133        nrows: usize,
134        ncols: usize,
135        triplets: &[Triplet<T>],
136    ) -> SparseResult<Self> {
137        let tuples = triplets
138            .iter()
139            .map(|triplet| (triplet.row, triplet.col, triplet.val));
140        build_from_triplets(nrows, ncols, tuples)
141    }
142}
143
144fn build_from_triplets<T, I>(nrows: usize, ncols: usize, triplets: I) -> SparseResult<CscMatrix<T>>
145where
146    T: Copy + AddAssign,
147    I: IntoIterator<Item = (usize, usize, T)>,
148{
149    let mut sorted = Vec::new();
150    for (row, col, value) in triplets {
151        if row >= nrows || col >= ncols {
152            return Err(SparseError::TripletOutOfBounds {
153                row,
154                col,
155                nrows,
156                ncols,
157            });
158        }
159        sorted.push((row, col, value));
160    }
161
162    sorted.sort_by(|lhs, rhs| lhs.1.cmp(&rhs.1).then(lhs.0.cmp(&rhs.0)));
163
164    let mut row_indices = Vec::with_capacity(sorted.len());
165    let mut values = Vec::with_capacity(sorted.len());
166    let mut entry_cols = Vec::with_capacity(sorted.len());
167
168    let mut previous = None;
169    for (row, col, value) in sorted {
170        if previous == Some((row, col)) {
171            *values
172                .last_mut()
173                .expect("duplicate entry requires an existing merged value") += value;
174        } else {
175            row_indices.push(row);
176            values.push(value);
177            entry_cols.push(col);
178            previous = Some((row, col));
179        }
180    }
181
182    let mut col_ptrs = vec![0usize; ncols + 1];
183    for col in entry_cols {
184        col_ptrs[col + 1] += 1;
185    }
186    for col in 0..ncols {
187        col_ptrs[col + 1] += col_ptrs[col];
188    }
189
190    CscMatrix::try_new(nrows, ncols, col_ptrs, row_indices, values)
191}
192
193pub(crate) fn validate_csc_pattern(
194    nrows: usize,
195    ncols: usize,
196    col_ptrs: &[usize],
197    row_indices: &[usize],
198) -> SparseResult<()> {
199    let expected = ncols + 1;
200    if col_ptrs.len() != expected {
201        return Err(SparseError::InvalidColPtrLen {
202            expected,
203            found: col_ptrs.len(),
204        });
205    }
206
207    let Some(&first) = col_ptrs.first() else {
208        return Err(SparseError::InvalidColPtrLen { expected, found: 0 });
209    };
210    if first != 0 {
211        return Err(SparseError::ColPtrMustStartAtZero { found: first });
212    }
213
214    for index in 1..col_ptrs.len() {
215        if col_ptrs[index] < col_ptrs[index - 1] {
216            return Err(SparseError::ColPtrNotMonotonic {
217                index,
218                current: col_ptrs[index],
219                prev_index: index - 1,
220                previous: col_ptrs[index - 1],
221            });
222        }
223    }
224
225    let last = *col_ptrs.last().expect("validated non-empty col_ptrs");
226    if last != row_indices.len() {
227        return Err(SparseError::InvalidColPtrEnd {
228            expected: row_indices.len(),
229            found: last,
230        });
231    }
232
233    for (index, &row) in row_indices.iter().enumerate() {
234        if row >= nrows {
235            return Err(SparseError::RowIndexOutOfBounds { index, row, nrows });
236        }
237    }
238
239    for col in 0..ncols {
240        let start = col_ptrs[col];
241        let end = col_ptrs[col + 1];
242        for index in start + 1..end {
243            if row_indices[index] <= row_indices[index - 1] {
244                return Err(SparseError::RowIndicesNotStrictlyIncreasing {
245                    col,
246                    previous: row_indices[index - 1],
247                    current: row_indices[index],
248                });
249            }
250        }
251    }
252
253    Ok(())
254}
255
256pub(crate) fn try_usize_to_i32(what: &'static str, value: usize) -> SparseResult<i32> {
257    i32::try_from(value).map_err(|_| SparseError::IndexTooLarge { what, value })
258}
259
260#[cfg(test)]
261mod tests {
262    use super::*;
263
264    #[test]
265    fn triplets_merge_duplicates_into_sorted_csc() {
266        let matrix = CscMatrix::try_from_triplets(
267            3,
268            3,
269            &[
270                Triplet {
271                    row: 2,
272                    col: 2,
273                    val: 1.0,
274                },
275                Triplet {
276                    row: 0,
277                    col: 0,
278                    val: 1.0,
279                },
280                Triplet {
281                    row: 1,
282                    col: 0,
283                    val: 4.0,
284                },
285                Triplet {
286                    row: 0,
287                    col: 0,
288                    val: 2.0,
289                },
290            ],
291        )
292        .expect("triplets should produce a valid CSC matrix");
293
294        assert_eq!(matrix.col_ptrs(), &[0, 2, 2, 3]);
295        assert_eq!(matrix.row_indices(), &[0, 1, 2]);
296        assert_eq!(matrix.values(), &[3.0, 4.0, 1.0]);
297    }
298
299    #[test]
300    fn rejects_out_of_bounds_triplet() {
301        let error = CscMatrix::try_from_triplets(
302            2,
303            2,
304            &[
305                Triplet {
306                    row: 0,
307                    col: 0,
308                    val: 1.0,
309                },
310                Triplet {
311                    row: 2,
312                    col: 1,
313                    val: 5.0,
314                },
315            ],
316        )
317        .expect_err("triplets with invalid coordinates must be rejected");
318        assert!(matches!(
319            error,
320            SparseError::TripletOutOfBounds {
321                row: 2,
322                col: 1,
323                nrows: 2,
324                ncols: 2
325            }
326        ));
327    }
328
329    #[test]
330    fn rejects_duplicate_rows_in_existing_csc() {
331        let error = CscMatrix::try_new(2, 1, vec![0, 2], vec![0, 0], vec![1.0, 2.0])
332            .expect_err("existing CSC with duplicate row indices must be rejected");
333        assert!(matches!(
334            error,
335            SparseError::RowIndicesNotStrictlyIncreasing { .. }
336        ));
337    }
338
339    #[test]
340    fn checked_i32_conversion_rejects_large_indices() {
341        let matrix =
342            CscMatrix::try_new(1, 1, vec![0, 1], vec![0], vec![2.0]).expect("valid matrix");
343        let overflow = CscMatrix::try_new(1, 1, vec![0, i32::MAX as usize + 1], vec![0], vec![2.0])
344            .expect_err("oversized CSC structure must be rejected by validation");
345        assert!(matches!(overflow, SparseError::InvalidColPtrEnd { .. }));
346
347        let (col_ptrs, row_indices, values) =
348            matrix.try_to_i32().expect("small indices should convert");
349        assert_eq!(col_ptrs, vec![0, 1]);
350        assert_eq!(row_indices, vec![0]);
351        assert_eq!(values, vec![2.0]);
352    }
353
354    /// Empty matrix: all entries are zero (no triplets) for a 3x3 matrix.
355    #[test]
356    fn empty_matrix_all_zeros() {
357        let mat = CscMatrix::<f64>::try_from_triplets(3, 3, &[]).unwrap();
358        assert_eq!(mat.nrows(), 3);
359        assert_eq!(mat.ncols(), 3);
360        assert_eq!(mat.nnz(), 0);
361        assert_eq!(mat.col_ptrs(), &[0, 0, 0, 0]);
362        assert!(mat.row_indices().is_empty());
363        assert!(mat.values().is_empty());
364    }
365
366    /// Out-of-order triplets that still merge correctly.
367    #[test]
368    fn out_of_order_triplets_merge_correctly() {
369        // Triplets submitted in reverse order: (1,1), (0,1), (1,0), (0,0)
370        let triplets = vec![
371            Triplet {
372                row: 1,
373                col: 1,
374                val: 4.0,
375            },
376            Triplet {
377                row: 0,
378                col: 1,
379                val: 2.0,
380            },
381            Triplet {
382                row: 1,
383                col: 0,
384                val: 3.0,
385            },
386            Triplet {
387                row: 0,
388                col: 0,
389                val: 1.0,
390            },
391        ];
392        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
393        // CSC should be column-sorted with rows in order.
394        assert_eq!(mat.col_ptrs(), &[0, 2, 4]);
395        assert_eq!(mat.row_indices(), &[0, 1, 0, 1]);
396        assert_eq!(mat.values(), &[1.0, 3.0, 2.0, 4.0]);
397    }
398
399    /// Dense matrix: every entry present in a 3x3 matrix.
400    #[test]
401    fn dense_matrix_all_entries() {
402        let mut triplets = Vec::new();
403        for row in 0..3 {
404            for col in 0..3 {
405                triplets.push(Triplet {
406                    row,
407                    col,
408                    val: (row * 3 + col + 1) as f64,
409                });
410            }
411        }
412        let mat = CscMatrix::try_from_triplets(3, 3, &triplets).unwrap();
413        assert_eq!(mat.nnz(), 9);
414        assert_eq!(mat.col_ptrs(), &[0, 3, 6, 9]);
415        // Column 0: rows 0,1,2 with values 1,4,7
416        assert_eq!(mat.row_indices(), &[0, 1, 2, 0, 1, 2, 0, 1, 2]);
417        assert_eq!(mat.values(), &[1.0, 4.0, 7.0, 2.0, 5.0, 8.0, 3.0, 6.0, 9.0]);
418    }
419
420    /// Duplicates across multiple triplets at the same position are summed.
421    #[test]
422    fn multiple_duplicates_summed() {
423        let triplets = vec![
424            Triplet {
425                row: 0,
426                col: 0,
427                val: 1.0,
428            },
429            Triplet {
430                row: 0,
431                col: 0,
432                val: 2.0,
433            },
434            Triplet {
435                row: 0,
436                col: 0,
437                val: 3.0,
438            },
439            Triplet {
440                row: 1,
441                col: 1,
442                val: 10.0,
443            },
444        ];
445        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
446        assert_eq!(mat.nnz(), 2);
447        assert!((mat.values()[0] - 6.0_f64).abs() < 1e-14); // 1 + 2 + 3
448        assert!((mat.values()[1] - 10.0_f64).abs() < 1e-14);
449    }
450
451    /// 0x0 matrix can be constructed.
452    #[test]
453    fn zero_dimension_matrix() {
454        let mat = CscMatrix::<f64>::try_new(0, 0, vec![0], vec![], vec![]).unwrap();
455        assert_eq!(mat.nrows(), 0);
456        assert_eq!(mat.ncols(), 0);
457        assert_eq!(mat.nnz(), 0);
458        assert!(mat.is_square());
459    }
460
461    /// has_same_pattern correctly detects matching and non-matching patterns.
462    #[test]
463    fn has_same_pattern_check() {
464        let a = CscMatrix::try_new(2, 2, vec![0, 1, 2], vec![0, 1], vec![1.0, 2.0]).unwrap();
465        let b = CscMatrix::try_new(2, 2, vec![0, 1, 2], vec![0, 1], vec![9.0, 8.0]).unwrap();
466        let c = CscMatrix::try_new(
467            2,
468            2,
469            vec![0, 2, 4],
470            vec![0, 1, 0, 1],
471            vec![1.0, 2.0, 3.0, 4.0],
472        )
473        .unwrap();
474
475        assert!(a.has_same_pattern(&b), "same pattern, different values");
476        assert!(!a.has_same_pattern(&c), "different patterns");
477    }
478
479    /// into_parts decomposes the matrix correctly.
480    #[test]
481    fn into_parts_roundtrip() {
482        let mat = CscMatrix::try_new(2, 2, vec![0, 1, 2], vec![0, 1], vec![5.0, 6.0]).unwrap();
483        let (nrows, ncols, cp, ri, vals) = mat.into_parts();
484        assert_eq!(nrows, 2);
485        assert_eq!(ncols, 2);
486        assert_eq!(cp, vec![0, 1, 2]);
487        assert_eq!(ri, vec![0, 1]);
488        assert_eq!(vals, vec![5.0, 6.0]);
489    }
490
491    /// Single-column matrix (tall).
492    #[test]
493    fn single_column_matrix() {
494        let triplets = vec![
495            Triplet {
496                row: 0,
497                col: 0,
498                val: 1.0,
499            },
500            Triplet {
501                row: 2,
502                col: 0,
503                val: 3.0,
504            },
505        ];
506        let mat = CscMatrix::try_from_triplets(4, 1, &triplets).unwrap();
507        assert_eq!(mat.nrows(), 4);
508        assert_eq!(mat.ncols(), 1);
509        assert!(!mat.is_square());
510        assert_eq!(mat.nnz(), 2);
511        assert_eq!(mat.col_ptrs(), &[0, 2]);
512        assert_eq!(mat.row_indices(), &[0, 2]);
513    }
514
515    /// Rejects non-zero starting column pointer.
516    #[test]
517    fn rejects_nonzero_start_col_ptr() {
518        let err = CscMatrix::try_new(2, 2, vec![1, 1, 2], vec![1], vec![1.0]).unwrap_err();
519        assert!(matches!(err, SparseError::ColPtrMustStartAtZero { .. }));
520    }
521
522    /// Rejects non-monotonic column pointers.
523    #[test]
524    fn rejects_non_monotonic_col_ptrs() {
525        let err = CscMatrix::try_new(2, 2, vec![0, 2, 1], vec![0, 1], vec![1.0, 2.0]).unwrap_err();
526        assert!(matches!(err, SparseError::ColPtrNotMonotonic { .. }));
527    }
528
529    /// values_mut allows in-place modification.
530    #[test]
531    fn values_mut_modification() {
532        let mut mat = CscMatrix::try_new(2, 2, vec![0, 1, 2], vec![0, 1], vec![1.0, 2.0]).unwrap();
533        mat.values_mut()[0] = 99.0;
534        assert_eq!(mat.values(), &[99.0, 2.0]);
535    }
536}
surge_sparse/csc.rs

surge_sparse/
csc.rs
