umicp_core/
matrix.rs

1/*!
2# UMICP Matrix Operations
3
4High-performance matrix operations with SIMD optimization for UMICP protocol.
5*/
6
7use crate::error::{Result, UmicpError};
8use crate::types::MatrixResult;
9
10/// Matrix operations class with high-performance implementations
11#[derive(Debug)]
12pub struct Matrix {
13    // Internal state can be added here if needed
14}
15
16impl Matrix {
17    /// Create a new matrix operations instance
18    pub fn new() -> Self {
19        Matrix {}
20    }
21
22    /// Matrix addition: result = a + b
23    /// Matrices must have the same dimensions
24    pub fn add(&self, a: &[f32], b: &[f32], result: &mut [f32], rows: usize, cols: usize) -> Result<MatrixResult> {
25        self.validate_dimensions(a.len(), b.len(), result.len(), rows, cols)?;
26
27        // Use parallel processing for large matrices
28        if rows * cols > 1000 {
29            self.add_parallel(a, b, result, rows, cols);
30        } else {
31            self.add_sequential(a, b, result);
32        }
33
34        Ok(MatrixResult {
35            success: true,
36            error: None,
37            result: None,
38            similarity: None,
39            data: None,
40        })
41    }
42
43    /// Matrix multiplication: result = a * b (m x n) * (n x p) = (m x p)
44    pub fn multiply(&self, a: &[f32], b: &[f32], result: &mut [f32], m: usize, n: usize, p: usize) -> Result<MatrixResult> {
45        let a_len = m * n;
46        let b_len = n * p;
47        let result_len = m * p;
48
49        if a.len() != a_len || b.len() != b_len || result.len() != result_len {
50            return Err(UmicpError::matrix(format!(
51                "Invalid matrix dimensions: a({}) != {}x{}, b({}) != {}x{}, result({}) != {}x{}",
52                a.len(), m, n, b.len(), n, p, result.len(), m, p
53            )));
54        }
55
56        // Initialize result to zeros
57        result.fill(0.0);
58
59        // Use parallel processing for large matrices
60        if m * n * p > 10000 {
61            self.multiply_parallel(a, b, result, m, n, p);
62        } else {
63            self.multiply_sequential(a, b, result, m, n, p);
64        }
65
66        Ok(MatrixResult {
67            success: true,
68            error: None,
69            result: None,
70            similarity: None,
71            data: None,
72        })
73    }
74
75    /// Matrix transpose: result = a^T
76    pub fn transpose(&self, input: &[f32], output: &mut [f32], rows: usize, cols: usize) -> Result<MatrixResult> {
77        let input_len = rows * cols;
78        let output_len = cols * rows;
79
80        if input.len() != input_len || output.len() != output_len {
81            return Err(UmicpError::matrix(format!(
82                "Invalid transpose dimensions: input({}) != {}x{}, output({}) != {}x{}",
83                input.len(), rows, cols, output.len(), cols, rows
84            )));
85        }
86
87        // Transpose operation
88        for i in 0..rows {
89            for j in 0..cols {
90                output[j * rows + i] = input[i * cols + j];
91            }
92        }
93
94        Ok(MatrixResult {
95            success: true,
96            error: None,
97            result: None,
98            similarity: None,
99            data: None,
100        })
101    }
102
103    /// Dot product of two vectors
104    pub fn dot_product(&self, a: &[f32], b: &[f32]) -> Result<MatrixResult> {
105        if a.len() != b.len() {
106            return Err(UmicpError::matrix(format!(
107                "Vector length mismatch: a({}) != b({})",
108                a.len(), b.len()
109            )));
110        }
111
112        // Use SIMD for large vectors
113        let result = if a.len() >= 8 {
114            self.dot_product_simd(a, b)
115        } else {
116            a.iter().zip(b.iter()).map(|(x, y)| x * y).sum()
117        };
118
119        Ok(MatrixResult {
120            success: true,
121            error: None,
122            result: Some(result as f64),
123            similarity: None,
124            data: None,
125        })
126    }
127
128    /// Vector/matrix normalization (L2 normalization)
129    pub fn normalize(&self, matrix: &mut [f32], rows: usize, cols: usize) -> Result<MatrixResult> {
130        let matrix_len = rows * cols;
131        if matrix.len() != matrix_len {
132            return Err(UmicpError::matrix(format!(
133                "Invalid matrix dimensions: matrix({}) != {}x{}",
134                matrix.len(), rows, cols
135            )));
136        }
137
138        for row in 0..rows {
139            let start = row * cols;
140            let end = start + cols;
141            let row_slice = &mut matrix[start..end];
142
143            // Calculate L2 norm
144            let norm: f32 = row_slice.iter().map(|x| x * x).sum::<f32>().sqrt();
145
146            if norm > 0.0 {
147                // Normalize the row
148                for val in row_slice.iter_mut() {
149                    *val /= norm;
150                }
151            }
152        }
153
154        Ok(MatrixResult {
155            success: true,
156            error: None,
157            result: None,
158            similarity: None,
159            data: Some(matrix.to_vec()),
160        })
161    }
162
163    /// Cosine similarity between two vectors
164    pub fn cosine_similarity(&self, a: &[f32], b: &[f32]) -> Result<MatrixResult> {
165        if a.len() != b.len() {
166            return Err(UmicpError::matrix(format!(
167                "Vector length mismatch: a({}) != b({})",
168                a.len(), b.len()
169            )));
170        }
171
172        // Calculate dot product
173        let dot_result = self.dot_product(a, b)?;
174        let dot_product = dot_result.result.unwrap() as f32;
175
176        // Calculate magnitudes
177        let a_magnitude: f32 = a.iter().map(|x| x * x).sum::<f32>().sqrt();
178        let b_magnitude: f32 = b.iter().map(|x| x * x).sum::<f32>().sqrt();
179
180        if a_magnitude == 0.0 || b_magnitude == 0.0 {
181            return Ok(MatrixResult {
182                success: true,
183                error: None,
184                result: None,
185                similarity: Some(0.0),
186                data: None,
187            });
188        }
189
190        let similarity = dot_product / (a_magnitude * b_magnitude);
191
192        Ok(MatrixResult {
193            success: true,
194            error: None,
195            result: None,
196            similarity: Some(similarity as f64),
197            data: None,
198        })
199    }
200
201    /// Element-wise vector addition
202    pub fn vector_add(&self, a: &[f32], b: &[f32], result: &mut [f32]) -> Result<MatrixResult> {
203        if a.len() != b.len() || a.len() != result.len() {
204            return Err(UmicpError::matrix(format!(
205                "Vector length mismatch: a({}), b({}), result({})",
206                a.len(), b.len(), result.len()
207            )));
208        }
209
210        for i in 0..a.len() {
211            result[i] = a[i] + b[i];
212        }
213
214        Ok(MatrixResult {
215            success: true,
216            error: None,
217            result: None,
218            similarity: None,
219            data: None,
220        })
221    }
222
223    /// Element-wise vector subtraction
224    pub fn vector_subtract(&self, a: &[f32], b: &[f32], result: &mut [f32]) -> Result<MatrixResult> {
225        if a.len() != b.len() || a.len() != result.len() {
226            return Err(UmicpError::matrix(format!(
227                "Vector length mismatch: a({}), b({}), result({})",
228                a.len(), b.len(), result.len()
229            )));
230        }
231
232        for i in 0..a.len() {
233            result[i] = a[i] - b[i];
234        }
235
236        Ok(MatrixResult {
237            success: true,
238            error: None,
239            result: None,
240            similarity: None,
241            data: None,
242        })
243    }
244
245    /// Element-wise vector multiplication (Hadamard product)
246    pub fn vector_multiply(&self, a: &[f32], b: &[f32], result: &mut [f32]) -> Result<MatrixResult> {
247        if a.len() != b.len() || a.len() != result.len() {
248            return Err(UmicpError::matrix(format!(
249                "Vector length mismatch: a({}), b({}), result({})",
250                a.len(), b.len(), result.len()
251            )));
252        }
253
254        for i in 0..a.len() {
255            result[i] = a[i] * b[i];
256        }
257
258        Ok(MatrixResult {
259            success: true,
260            error: None,
261            result: None,
262            similarity: None,
263            data: None,
264        })
265    }
266
267    /// Scalar multiplication of vector
268    pub fn vector_scale(&self, vector: &[f32], scalar: f32, result: &mut [f32]) -> Result<MatrixResult> {
269        if vector.len() != result.len() {
270            return Err(UmicpError::matrix(format!(
271                "Vector length mismatch: vector({}), result({})",
272                vector.len(), result.len()
273            )));
274        }
275
276        for i in 0..vector.len() {
277            result[i] = vector[i] * scalar;
278        }
279
280        Ok(MatrixResult {
281            success: true,
282            error: None,
283            result: None,
284            similarity: None,
285            data: None,
286        })
287    }
288
289    /// Calculate matrix determinant (for square matrices only)
290    pub fn determinant(&self, matrix: &[f32], size: usize) -> Result<MatrixResult> {
291        let matrix_len = size * size;
292        if matrix.len() != matrix_len {
293            return Err(UmicpError::matrix(format!(
294                "Invalid matrix dimensions for determinant: matrix({}) != {}x{}",
295                matrix.len(), size, size
296            )));
297        }
298
299        if size == 1 {
300            return Ok(MatrixResult {
301                success: true,
302                error: None,
303                result: Some(matrix[0] as f64),
304                similarity: None,
305                data: None,
306            });
307        }
308
309        if size == 2 {
310            let det = matrix[0] * matrix[3] - matrix[1] * matrix[2];
311            return Ok(MatrixResult {
312                success: true,
313                error: None,
314                result: Some(det as f64),
315                similarity: None,
316                data: None,
317            });
318        }
319
320        // For larger matrices, use LAPACK if available, otherwise return error
321        Err(UmicpError::matrix("Determinant calculation for matrices larger than 2x2 not yet implemented"))
322    }
323
324    /// Matrix inverse (for square matrices only)
325    pub fn inverse(&self, matrix: &[f32], result: &mut [f32], size: usize) -> Result<MatrixResult> {
326        let matrix_len = size * size;
327        if matrix.len() != matrix_len || result.len() != matrix_len {
328            return Err(UmicpError::matrix(format!(
329                "Invalid matrix dimensions for inverse: matrix({}) != {}x{}, result({}) != {}x{}",
330                matrix.len(), size, size, result.len(), size, size
331            )));
332        }
333
334        // Simple 2x2 inverse implementation
335        if size == 2 {
336            let det = matrix[0] * matrix[3] - matrix[1] * matrix[2];
337            if det == 0.0 {
338                return Err(UmicpError::matrix("Matrix is singular, cannot compute inverse"));
339            }
340
341            result[0] = matrix[3] / det;
342            result[1] = -matrix[1] / det;
343            result[2] = -matrix[2] / det;
344            result[3] = matrix[0] / det;
345
346            return Ok(MatrixResult {
347                success: true,
348                error: None,
349                result: None,
350                similarity: None,
351                data: Some(result.to_vec()),
352            });
353        }
354
355        // For larger matrices, use LAPACK if available, otherwise return error
356        Err(UmicpError::matrix("Matrix inverse for matrices larger than 2x2 not yet implemented"))
357    }
358
359    // Private helper methods
360
361    fn validate_dimensions(&self, a_len: usize, b_len: usize, result_len: usize, rows: usize, cols: usize) -> Result<()> {
362        let expected_len = rows * cols;
363        if a_len != expected_len || b_len != expected_len || result_len != expected_len {
364            return Err(UmicpError::matrix(format!(
365                "Invalid matrix dimensions: expected {}x{} ({} elements), got a({}), b({}), result({})",
366                rows, cols, expected_len, a_len, b_len, result_len
367            )));
368        }
369        Ok(())
370    }
371
372    fn add_sequential(&self, a: &[f32], b: &[f32], result: &mut [f32]) {
373        for i in 0..a.len() {
374            result[i] = a[i] + b[i];
375        }
376    }
377
378    fn add_parallel(&self, a: &[f32], b: &[f32], result: &mut [f32], _rows: usize, _cols: usize) {
379        // Sequential implementation for compatibility
380        for i in 0..a.len() {
381            result[i] = a[i] + b[i];
382        }
383    }
384
385    fn multiply_sequential(&self, a: &[f32], b: &[f32], result: &mut [f32], m: usize, n: usize, p: usize) {
386        for i in 0..m {
387            for j in 0..p {
388                for k in 0..n {
389                    result[i * p + j] += a[i * n + k] * b[k * p + j];
390                }
391            }
392        }
393    }
394
395    fn multiply_parallel(&self, a: &[f32], b: &[f32], result: &mut [f32], m: usize, n: usize, p: usize) {
396        // Sequential implementation for compatibility
397        for i in 0..m {
398            for j in 0..p {
399                let mut sum = 0.0;
400                for k in 0..n {
401                    sum += a[i * n + k] * b[k * p + j];
402                }
403                result[i * p + j] = sum;
404            }
405        }
406    }
407
408    fn dot_product_simd(&self, a: &[f32], b: &[f32]) -> f32 {
409        // Fallback to regular implementation for now
410        // In a real implementation, this would use SIMD intrinsics
411        a.iter().zip(b.iter()).map(|(x, y)| x * y).sum()
412    }
413}
414
415impl Default for Matrix {
416    fn default() -> Self {
417        Self::new()
418    }
419}
420
421#[cfg(test)]
422mod tests {
423    use super::*;
424
425    #[test]
426    fn test_matrix_creation() {
427        let _matrix = Matrix::new();
428        // Just verify it can be created
429        assert!(true);
430    }
431
432    #[test]
433    fn test_vector_add() {
434        let matrix = Matrix::new();
435        let a = vec![1.0, 2.0, 3.0];
436        let b = vec![4.0, 5.0, 6.0];
437        let mut result = vec![0.0; 3];
438
439        let matrix_result = matrix.vector_add(&a, &b, &mut result).unwrap();
440        assert!(matrix_result.success);
441        assert_eq!(result, vec![5.0, 7.0, 9.0]);
442    }
443
444    #[test]
445    fn test_dot_product() {
446        let matrix = Matrix::new();
447        let a = vec![1.0, 2.0, 3.0];
448        let b = vec![4.0, 5.0, 6.0];
449
450        let result = matrix.dot_product(&a, &b).unwrap();
451        assert!(result.success);
452        assert_eq!(result.result.unwrap(), 32.0); // 1*4 + 2*5 + 3*6 = 32
453    }
454
455    #[test]
456    fn test_cosine_similarity() {
457        let matrix = Matrix::new();
458        let a = vec![1.0, 2.0, 3.0];
459        let b = vec![1.0, 2.0, 3.0]; // Identical vectors
460
461        let result = matrix.cosine_similarity(&a, &b).unwrap();
462        assert!(result.success);
463        assert!((result.similarity.unwrap() - 1.0).abs() < 1e-6); // Should be 1.0
464    }
465
466    #[test]
467    fn test_matrix_multiply() {
468        let matrix = Matrix::new();
469        // 2x2 * 2x2 = 2x2
470        let a = vec![1.0, 2.0, 3.0, 4.0];
471        let b = vec![5.0, 6.0, 7.0, 8.0];
472        let mut result = vec![0.0; 4];
473
474        let matrix_result = matrix.multiply(&a, &b, &mut result, 2, 2, 2).unwrap();
475        assert!(matrix_result.success);
476        // Expected: [19, 22, 43, 50]
477        assert_eq!(result[0], 19.0);
478        assert_eq!(result[1], 22.0);
479        assert_eq!(result[2], 43.0);
480        assert_eq!(result[3], 50.0);
481    }
482
483    #[test]
484    fn test_matrix_transpose() {
485        let matrix = Matrix::new();
486        let input = vec![1.0, 2.0, 3.0, 4.0]; // 2x2 matrix
487        let mut output = vec![0.0; 4];
488
489        let result = matrix.transpose(&input, &mut output, 2, 2).unwrap();
490        assert!(result.success);
491        assert_eq!(output, vec![1.0, 3.0, 2.0, 4.0]);
492    }
493
494    #[test]
495    fn test_determinant_2x2() {
496        let matrix = Matrix::new();
497        let mat = vec![1.0, 2.0, 3.0, 4.0]; // det = 1*4 - 2*3 = -2
498
499        let result = matrix.determinant(&mat, 2).unwrap();
500        assert!(result.success);
501        assert_eq!(result.result.unwrap(), -2.0);
502    }
503
504    #[test]
505    fn test_validation_errors() {
506        let matrix = Matrix::new();
507
508        // Test vector length mismatch
509        let a = vec![1.0, 2.0];
510        let b = vec![1.0, 2.0, 3.0];
511        let mut result = vec![0.0; 2];
512
513        let matrix_result = matrix.vector_add(&a, &b, &mut result);
514        assert!(matrix_result.is_err());
515    }
516}
umicp_core/matrix.rs

umicp_core/
matrix.rs
