## Coverage Report

#### Created: 2026-01-25 15:05

next uncovered line (L), next uncovered region (R), next uncovered branch (B)

    /home/noah/src/trueno/src/vector/ops/reductions.rs

    Line
    Count
    Source
    1
    //! Reduction operations for Vector<f32>
    2
    //!
    3
    //! This module provides reduction operations that aggregate vector elements:
    4
    //! - Basic: `sum`, `dot`, `max`, `min`
    5
    //! - Index-finding: `argmax`, `argmin`
    6
    //! - Statistical: `mean`, `variance`, `stddev`, `covariance`, `correlation`
    7
    //! - Numerically stable: `sum_kahan`, `sum_of_squares`
    8
    9
    #[cfg(target_arch = "x86_64")]
    10
    use crate::backends::avx2::Avx2Backend;
    11
    #[cfg(target_arch = "x86_64")]
    12
    use crate::backends::avx512::Avx512Backend;
    13
    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    14
    use crate::backends::neon::NeonBackend;
    15
    use crate::backends::scalar::ScalarBackend;
    16
    #[cfg(target_arch = "x86_64")]
    17
    use crate::backends::sse2::Sse2Backend;
    18
    #[cfg(target_arch = "wasm32")]
    19
    use crate::backends::wasm::WasmBackend;
    20
    use crate::backends::VectorBackend;
    21
    use crate::vector::Vector;
    22
    use crate::{dispatch_reduction, Backend, Result, TruenoError};
    23
    24
    impl Vector<f32> {
    25
        /// Dot product
    26
        ///
    27
        /// # Examples
    28
        ///
    29
        /// ```
    30
        /// use trueno::Vector;
    31
        ///
    32
        /// let a = Vector::from_slice(&[1.0, 2.0, 3.0]);
    33
        /// let b = Vector::from_slice(&[4.0, 5.0, 6.0]);
    34
        /// let result = a.dot(&b)?;
    35
        ///
    36
        /// assert_eq!(result, 32.0); // 1*4 + 2*5 + 3*6 = 4 + 10 + 18 = 32
    37
        /// # Ok::<(), trueno::TruenoError>(())
    38
        /// ```
    39
    226k
        pub fn dot(&self, other: &Self) -> Result<f32> {
    40
    226k
            if self.len() != other.len() {
    41
    0
                return Err(TruenoError::SizeMismatch {
    42
    0
                    expected: self.len(),
    43
    0
                    actual: other.len(),
    44
    0
                });
    45
    226k
            }
    46
    47
            // SAFETY: Unsafe block delegates to backend implementation which maintains safety invariants
    48
    226k
            let result = unsafe {
    49
    226k
                match self.backend {
    50
    0
                    Backend::Scalar => ScalarBackend::dot(&self.data, &other.data),
    51
                    #[cfg(target_arch = "x86_64")]
    52
    0
                    Backend::SSE2 | Backend::AVX => Sse2Backend::dot(&self.data, &other.data),
    53
                    #[cfg(target_arch = "x86_64")]
    54
    226k
                    Backend::AVX2 => Avx2Backend::dot(&self.data, &other.data),
    55
                    #[cfg(target_arch = "x86_64")]
    56
    0
                    Backend::AVX512 => Avx512Backend::dot(&self.data, &other.data),
    57
                    #[cfg(not(target_arch = "x86_64"))]
    58
                    Backend::SSE2 | Backend::AVX | Backend::AVX2 | Backend::AVX512 => {
    59
                        ScalarBackend::dot(&self.data, &other.data)
    60
                    }
    61
                    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    62
                    Backend::NEON => NeonBackend::dot(&self.data, &other.data),
    63
                    #[cfg(not(any(target_arch = "aarch64", target_arch = "arm")))]
    64
    0
                    Backend::NEON => ScalarBackend::dot(&self.data, &other.data),
    65
                    #[cfg(target_arch = "wasm32")]
    66
                    Backend::WasmSIMD => WasmBackend::dot(&self.data, &other.data),
    67
                    #[cfg(not(target_arch = "wasm32"))]
    68
    0
                    Backend::WasmSIMD => ScalarBackend::dot(&self.data, &other.data),
    69
    0
                    Backend::GPU | Backend::Auto => ScalarBackend::dot(&self.data, &other.data),
    70
                }
    71
            };
    72
    73
    226k
            Ok(result)
    74
    226k
        }
    75
    76
        /// Sum all elements
    77
        ///
    78
        /// # Examples
    79
        ///
    80
        /// ```
    81
        /// use trueno::Vector;
    82
        ///
    83
        /// let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0]);
    84
        /// assert_eq!(v.sum()?, 10.0);
    85
        /// # Ok::<(), trueno::TruenoError>(())
    86
        /// ```
    87
    118
        pub fn sum(&self) -> Result<f32> {
    88
    118
            Ok(dispatch_reduction!0(self.backend, sum, &self.data0))
    89
    118
        }
    90
    91
        /// Find maximum element
    92
        ///
    93
        /// # Examples
    94
        ///
    95
        /// ```
    96
        /// use trueno::Vector;
    97
        ///
    98
        /// let v = Vector::from_slice(&[1.0, 5.0, 3.0, 2.0]);
    99
        /// assert_eq!(v.max()?, 5.0);
    100
        /// # Ok::<(), trueno::TruenoError>(())
    101
        /// ```
    102
        ///
    103
        /// # Errors
    104
        ///
    105
        /// Returns [`TruenoError::InvalidInput`] if vector is empty.
    106
    10.4k
        pub fn max(&self) -> Result<f32> {
    107
    10.4k
            if self.data.is_empty() {
    108
    0
                return Err(TruenoError::InvalidInput("Empty vector".to_string()));
    109
    10.4k
            }
    110
    111
            // SAFETY: Unsafe block delegates to backend implementation which maintains safety invariants
    112
    10.4k
            let result = unsafe {
    113
    10.4k
                match self.backend {
    114
    0
                    Backend::Scalar => ScalarBackend::max(&self.data),
    115
                    #[cfg(target_arch = "x86_64")]
    116
    0
                    Backend::SSE2 | Backend::AVX => Sse2Backend::max(&self.data),
    117
                    #[cfg(target_arch = "x86_64")]
    118
    10.4k
                    Backend::AVX2 | Backend::AVX512 => Avx2Backend::max(&self.data),
    119
                    #[cfg(not(target_arch = "x86_64"))]
    120
                    Backend::SSE2 | Backend::AVX | Backend::AVX2 | Backend::AVX512 => {
    121
                        ScalarBackend::max(&self.data)
    122
                    }
    123
                    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    124
                    Backend::NEON => NeonBackend::max(&self.data),
    125
                    #[cfg(not(any(target_arch = "aarch64", target_arch = "arm")))]
    126
    0
                    Backend::NEON => ScalarBackend::max(&self.data),
    127
                    #[cfg(target_arch = "wasm32")]
    128
                    Backend::WasmSIMD => WasmBackend::max(&self.data),
    129
                    #[cfg(not(target_arch = "wasm32"))]
    130
    0
                    Backend::WasmSIMD => ScalarBackend::max(&self.data),
    131
    0
                    Backend::GPU | Backend::Auto => ScalarBackend::max(&self.data),
    132
                }
    133
            };
    134
    135
    10.4k
            Ok(result)
    136
    10.4k
        }
    137
    138
        /// Find minimum value in the vector
    139
        ///
    140
        /// Returns the smallest element in the vector using SIMD optimization.
    141
        ///
    142
        /// # Examples
    143
        ///
    144
        /// ```
    145
        /// use trueno::Vector;
    146
        ///
    147
        /// let v = Vector::from_slice(&[1.0, 5.0, 3.0, 2.0]);
    148
        /// assert_eq!(v.min()?, 1.0);
    149
        /// # Ok::<(), trueno::TruenoError>(())
    150
        /// ```
    151
        ///
    152
        /// # Errors
    153
        ///
    154
        /// Returns [`TruenoError::InvalidInput`] if vector is empty.
    155
    0
        pub fn min(&self) -> Result<f32> {
    156
    0
            if self.data.is_empty() {
    157
    0
                return Err(TruenoError::InvalidInput("Empty vector".to_string()));
    158
    0
            }
    159
    160
            // SAFETY: Unsafe block delegates to backend implementation which maintains safety invariants
    161
    0
            let result = unsafe {
    162
    0
                match self.backend {
    163
    0
                    Backend::Scalar => ScalarBackend::min(&self.data),
    164
                    #[cfg(target_arch = "x86_64")]
    165
    0
                    Backend::SSE2 | Backend::AVX => Sse2Backend::min(&self.data),
    166
                    #[cfg(target_arch = "x86_64")]
    167
    0
                    Backend::AVX2 | Backend::AVX512 => Avx2Backend::min(&self.data),
    168
                    #[cfg(not(target_arch = "x86_64"))]
    169
                    Backend::SSE2 | Backend::AVX | Backend::AVX2 | Backend::AVX512 => {
    170
                        ScalarBackend::min(&self.data)
    171
                    }
    172
                    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    173
                    Backend::NEON => NeonBackend::min(&self.data),
    174
                    #[cfg(not(any(target_arch = "aarch64", target_arch = "arm")))]
    175
    0
                    Backend::NEON => ScalarBackend::min(&self.data),
    176
                    #[cfg(target_arch = "wasm32")]
    177
                    Backend::WasmSIMD => WasmBackend::min(&self.data),
    178
                    #[cfg(not(target_arch = "wasm32"))]
    179
    0
                    Backend::WasmSIMD => ScalarBackend::min(&self.data),
    180
    0
                    Backend::GPU | Backend::Auto => ScalarBackend::min(&self.data),
    181
                }
    182
            };
    183
    184
    0
            Ok(result)
    185
    0
        }
    186
    187
        /// Find index of maximum value in the vector
    188
        ///
    189
        /// Returns the index of the first occurrence of the maximum value using SIMD optimization.
    190
        ///
    191
        /// # Examples
    192
        ///
    193
        /// ```
    194
        /// use trueno::Vector;
    195
        ///
    196
        /// let v = Vector::from_slice(&[1.0, 5.0, 3.0, 2.0]);
    197
        /// assert_eq!(v.argmax()?, 1); // max value 5.0 is at index 1
    198
        /// # Ok::<(), trueno::TruenoError>(())
    199
        /// ```
    200
        ///
    201
        /// # Errors
    202
        ///
    203
        /// Returns [`TruenoError::InvalidInput`] if vector is empty.
    204
    0
        pub fn argmax(&self) -> Result<usize> {
    205
    0
            if self.data.is_empty() {
    206
    0
                return Err(TruenoError::InvalidInput("Empty vector".to_string()));
    207
    0
            }
    208
    209
            // SAFETY: Unsafe block delegates to backend implementation which maintains safety invariants
    210
    0
            let result = unsafe {
    211
    0
                match self.backend {
    212
    0
                    Backend::Scalar => ScalarBackend::argmax(&self.data),
    213
                    #[cfg(target_arch = "x86_64")]
    214
    0
                    Backend::SSE2 | Backend::AVX => Sse2Backend::argmax(&self.data),
    215
                    #[cfg(target_arch = "x86_64")]
    216
    0
                    Backend::AVX2 | Backend::AVX512 => Avx2Backend::argmax(&self.data),
    217
                    #[cfg(not(target_arch = "x86_64"))]
    218
                    Backend::SSE2 | Backend::AVX | Backend::AVX2 | Backend::AVX512 => {
    219
                        ScalarBackend::argmax(&self.data)
    220
                    }
    221
                    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    222
                    Backend::NEON => NeonBackend::argmax(&self.data),
    223
                    #[cfg(not(any(target_arch = "aarch64", target_arch = "arm")))]
    224
    0
                    Backend::NEON => ScalarBackend::argmax(&self.data),
    225
                    #[cfg(target_arch = "wasm32")]
    226
                    Backend::WasmSIMD => WasmBackend::argmax(&self.data),
    227
                    #[cfg(not(target_arch = "wasm32"))]
    228
    0
                    Backend::WasmSIMD => ScalarBackend::argmax(&self.data),
    229
    0
                    Backend::GPU | Backend::Auto => ScalarBackend::argmax(&self.data),
    230
                }
    231
            };
    232
    233
    0
            Ok(result)
    234
    0
        }
    235
    236
        /// Find index of minimum value in the vector
    237
        ///
    238
        /// Returns the index of the first occurrence of the minimum value using SIMD optimization.
    239
        ///
    240
        /// # Examples
    241
        ///
    242
        /// ```
    243
        /// use trueno::Vector;
    244
        ///
    245
        /// let v = Vector::from_slice(&[1.0, 5.0, 3.0, 2.0]);
    246
        /// assert_eq!(v.argmin()?, 0); // min value 1.0 is at index 0
    247
        /// # Ok::<(), trueno::TruenoError>(())
    248
        /// ```
    249
        ///
    250
        /// # Errors
    251
        ///
    252
        /// Returns [`TruenoError::InvalidInput`] if vector is empty.
    253
    0
        pub fn argmin(&self) -> Result<usize> {
    254
    0
            if self.data.is_empty() {
    255
    0
                return Err(TruenoError::InvalidInput("Empty vector".to_string()));
    256
    0
            }
    257
    258
            // SAFETY: Unsafe block delegates to backend implementation which maintains safety invariants
    259
    0
            let result = unsafe {
    260
    0
                match self.backend {
    261
    0
                    Backend::Scalar => ScalarBackend::argmin(&self.data),
    262
                    #[cfg(target_arch = "x86_64")]
    263
    0
                    Backend::SSE2 | Backend::AVX => Sse2Backend::argmin(&self.data),
    264
                    #[cfg(target_arch = "x86_64")]
    265
    0
                    Backend::AVX2 | Backend::AVX512 => Avx2Backend::argmin(&self.data),
    266
                    #[cfg(not(target_arch = "x86_64"))]
    267
                    Backend::SSE2 | Backend::AVX | Backend::AVX2 | Backend::AVX512 => {
    268
                        ScalarBackend::argmin(&self.data)
    269
                    }
    270
                    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    271
                    Backend::NEON => NeonBackend::argmin(&self.data),
    272
                    #[cfg(not(any(target_arch = "aarch64", target_arch = "arm")))]
    273
    0
                    Backend::NEON => ScalarBackend::argmin(&self.data),
    274
                    #[cfg(target_arch = "wasm32")]
    275
                    Backend::WasmSIMD => WasmBackend::argmin(&self.data),
    276
                    #[cfg(not(target_arch = "wasm32"))]
    277
    0
                    Backend::WasmSIMD => ScalarBackend::argmin(&self.data),
    278
    0
                    Backend::GPU | Backend::Auto => ScalarBackend::argmin(&self.data),
    279
                }
    280
            };
    281
    282
    0
            Ok(result)
    283
    0
        }
    284
    285
        /// Kahan summation (numerically stable sum)
    286
        ///
    287
        /// Uses the Kahan summation algorithm to reduce floating-point rounding errors
    288
        /// when summing many numbers. This is more accurate than the standard sum() method
    289
        /// for vectors with many elements or elements of vastly different magnitudes.
    290
        ///
    291
        /// # Performance
    292
        ///
    293
        /// Note: Kahan summation is inherently sequential and cannot be effectively
    294
        /// parallelized with SIMD. All backends use the scalar implementation.
    295
        ///
    296
        /// # Examples
    297
        ///
    298
        /// ```
    299
        /// use trueno::Vector;
    300
        ///
    301
        /// let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0]);
    302
        /// assert_eq!(v.sum_kahan()?, 10.0);
    303
        /// # Ok::<(), trueno::TruenoError>(())
    304
        /// ```
    305
    0
        pub fn sum_kahan(&self) -> Result<f32> {
    306
    0
            if self.data.is_empty() {
    307
    0
                return Ok(0.0);
    308
    0
            }
    309
    310
            // SAFETY: Unsafe block delegates to backend implementation which maintains safety invariants
    311
    0
            let result = unsafe {
    312
    0
                match self.backend {
    313
    0
                    Backend::Scalar => ScalarBackend::sum_kahan(&self.data),
    314
                    #[cfg(target_arch = "x86_64")]
    315
    0
                    Backend::SSE2 | Backend::AVX => Sse2Backend::sum_kahan(&self.data),
    316
                    #[cfg(target_arch = "x86_64")]
    317
    0
                    Backend::AVX2 | Backend::AVX512 => Avx2Backend::sum_kahan(&self.data),
    318
                    #[cfg(not(target_arch = "x86_64"))]
    319
                    Backend::SSE2 | Backend::AVX | Backend::AVX2 | Backend::AVX512 => {
    320
                        ScalarBackend::sum_kahan(&self.data)
    321
                    }
    322
                    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    323
                    Backend::NEON => NeonBackend::sum_kahan(&self.data),
    324
                    #[cfg(not(any(target_arch = "aarch64", target_arch = "arm")))]
    325
    0
                    Backend::NEON => ScalarBackend::sum_kahan(&self.data),
    326
                    #[cfg(target_arch = "wasm32")]
    327
                    Backend::WasmSIMD => WasmBackend::sum_kahan(&self.data),
    328
                    #[cfg(not(target_arch = "wasm32"))]
    329
    0
                    Backend::WasmSIMD => ScalarBackend::sum_kahan(&self.data),
    330
    0
                    Backend::GPU | Backend::Auto => ScalarBackend::sum_kahan(&self.data),
    331
                }
    332
            };
    333
    334
    0
            Ok(result)
    335
    0
        }
    336
    337
        /// Sum of squared elements
    338
        ///
    339
        /// Computes the sum of squares: sum(a\[i\]^2).
    340
        /// This is the building block for computing L2 norm and variance.
    341
        ///
    342
        /// # Examples
    343
        ///
    344
        /// ```
    345
        /// use trueno::Vector;
    346
        ///
    347
        /// let v = Vector::from_slice(&[1.0, 2.0, 3.0]);
    348
        /// let sum_sq = v.sum_of_squares()?;
    349
        /// assert_eq!(sum_sq, 14.0); // 1^2 + 2^2 + 3^2 = 1 + 4 + 9 = 14
    350
        /// # Ok::<(), trueno::TruenoError>(())
    351
        /// ```
    352
        ///
    353
        /// # Empty vectors
    354
        ///
    355
        /// Returns 0.0 for empty vectors.
    356
    1
        pub fn sum_of_squares(&self) -> Result<f32> {
    357
    1
            if self.data.is_empty() {
    358
    0
                return Ok(0.0);
    359
    1
            }
    360
    361
            // Use dot product with self: dot(self, self) = sum(a[i]^2)
    362
    1
            self.dot(self)
    363
    1
        }
    364
    365
        /// Arithmetic mean (average)
    366
        ///
    367
        /// Computes the arithmetic mean of all elements: sum(a\[i\]) / n.
    368
        ///
    369
        /// # Performance
    370
        ///
    371
        /// Uses optimized SIMD sum() implementation, then divides by length.
    372
        ///
    373
        /// # Examples
    374
        ///
    375
        /// ```
    376
        /// use trueno::Vector;
    377
        ///
    378
        /// let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0]);
    379
        /// let avg = v.mean()?;
    380
        /// assert!((avg - 2.5).abs() < 1e-5); // (1+2+3+4)/4 = 2.5
    381
        /// # Ok::<(), trueno::TruenoError>(())
    382
        /// ```
    383
        ///
    384
        /// # Empty vectors
    385
        ///
    386
        /// Returns an error for empty vectors (division by zero).
    387
        ///
    388
        /// ```
    389
        /// use trueno::{Vector, TruenoError};
    390
        ///
    391
        /// let v: Vector<f32> = Vector::from_slice(&[]);
    392
        /// assert!(matches!(v.mean(), Err(TruenoError::EmptyVector)));
    393
        /// ```
    394
    0
        pub fn mean(&self) -> Result<f32> {
    395
    0
            if self.data.is_empty() {
    396
    0
                return Err(TruenoError::EmptyVector);
    397
    0
            }
    398
    399
    0
            let total = self.sum()?;
    400
    0
            Ok(total / self.len() as f32)
    401
    0
        }
    402
    403
        /// Population variance
    404
        ///
    405
        /// Computes the population variance: Var(X) = E\[(X - μ)²\] = E\[X²\] - μ²
    406
        /// Uses the computational formula to avoid two passes over the data.
    407
        ///
    408
        /// # Performance
    409
        ///
    410
        /// Uses optimized SIMD implementations via sum_of_squares() and mean().
    411
        ///
    412
        /// # Examples
    413
        ///
    414
        /// ```
    415
        /// use trueno::Vector;
    416
        ///
    417
        /// let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0, 5.0]);
    418
        /// let var = v.variance()?;
    419
        /// assert!((var - 2.0).abs() < 1e-5); // Population variance
    420
        /// # Ok::<(), trueno::TruenoError>(())
    421
        /// ```
    422
        ///
    423
        /// # Empty vectors
    424
        ///
    425
        /// Returns an error for empty vectors.
    426
        ///
    427
        /// ```
    428
        /// use trueno::{Vector, TruenoError};
    429
        ///
    430
        /// let v: Vector<f32> = Vector::from_slice(&[]);
    431
        /// assert!(matches!(v.variance(), Err(TruenoError::EmptyVector)));
    432
        /// ```
    433
    0
        pub fn variance(&self) -> Result<f32> {
    434
    0
            if self.data.is_empty() {
    435
    0
                return Err(TruenoError::EmptyVector);
    436
    0
            }
    437
    438
    0
            let mean_val = self.mean()?;
    439
    0
            let sum_sq = self.sum_of_squares()?;
    440
    0
            let mean_sq = sum_sq / self.len() as f32;
    441
    442
            // Var(X) = E[X²] - μ²
    443
    0
            Ok(mean_sq - mean_val * mean_val)
    444
    0
        }
    445
    446
        /// Population standard deviation
    447
        ///
    448
        /// Computes the population standard deviation: σ = sqrt(Var(X)).
    449
        /// This is the square root of the variance.
    450
        ///
    451
        /// # Performance
    452
        ///
    453
        /// Uses optimized SIMD implementations via variance().
    454
        ///
    455
        /// # Examples
    456
        ///
    457
        /// ```
    458
        /// use trueno::Vector;
    459
        ///
    460
        /// let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0, 5.0]);
    461
        /// let sd = v.stddev()?;
    462
        /// assert!((sd - 1.4142135).abs() < 1e-5); // sqrt(2) ≈ 1.414
    463
        /// # Ok::<(), trueno::TruenoError>(())
    464
        /// ```
    465
        ///
    466
        /// # Empty vectors
    467
        ///
    468
        /// Returns an error for empty vectors.
    469
        ///
    470
        /// ```
    471
        /// use trueno::{Vector, TruenoError};
    472
        ///
    473
        /// let v: Vector<f32> = Vector::from_slice(&[]);
    474
        /// assert!(matches!(v.stddev(), Err(TruenoError::EmptyVector)));
    475
        /// ```
    476
    0
        pub fn stddev(&self) -> Result<f32> {
    477
    0
            let var = self.variance()?;
    478
    0
            Ok(var.sqrt())
    479
    0
        }
    480
    481
        /// Population covariance between two vectors
    482
        ///
    483
        /// Computes the population covariance: Cov(X,Y) = E[(X - μx)(Y - μy)]
    484
        /// Uses the computational formula: Cov(X,Y) = E\[XY\] - μx·μy
    485
        ///
    486
        /// # Performance
    487
        ///
    488
        /// Uses optimized SIMD implementations via dot() and mean().
    489
        ///
    490
        /// # Examples
    491
        ///
    492
        /// ```
    493
        /// use trueno::Vector;
    494
        ///
    495
        /// let x = Vector::from_slice(&[1.0, 2.0, 3.0]);
    496
        /// let y = Vector::from_slice(&[2.0, 4.0, 6.0]);
    497
        /// let cov = x.covariance(&y)?;
    498
        /// assert!((cov - 1.333).abs() < 0.01); // Perfect positive covariance
    499
        /// # Ok::<(), trueno::TruenoError>(())
    500
        /// ```
    501
        ///
    502
        /// # Size mismatch
    503
        ///
    504
        /// Returns an error if vectors have different lengths.
    505
        ///
    506
        /// ```
    507
        /// use trueno::{Vector, TruenoError};
    508
        ///
    509
        /// let x = Vector::from_slice(&[1.0, 2.0]);
    510
        /// let y = Vector::from_slice(&[1.0, 2.0, 3.0]);
    511
        /// assert!(matches!(x.covariance(&y), Err(TruenoError::SizeMismatch { .. })));
    512
        /// ```
    513
        ///
    514
        /// # Empty vectors
    515
        ///
    516
        /// Returns an error for empty vectors.
    517
        ///
    518
        /// ```
    519
        /// use trueno::{Vector, TruenoError};
    520
        ///
    521
        /// let x: Vector<f32> = Vector::from_slice(&[]);
    522
        /// let y: Vector<f32> = Vector::from_slice(&[]);
    523
        /// assert!(matches!(x.covariance(&y), Err(TruenoError::EmptyVector)));
    524
        /// ```
    525
    0
        pub fn covariance(&self, other: &Self) -> Result<f32> {
    526
    0
            if self.data.is_empty() {
    527
    0
                return Err(TruenoError::EmptyVector);
    528
    0
            }
    529
    0
            if self.len() != other.len() {
    530
    0
                return Err(TruenoError::SizeMismatch {
    531
    0
                    expected: self.len(),
    532
    0
                    actual: other.len(),
    533
    0
                });
    534
    0
            }
    535
    536
    0
            let mean_x = self.mean()?;
    537
    0
            let mean_y = other.mean()?;
    538
    0
            let dot_xy = self.dot(other)?;
    539
    0
            let mean_xy = dot_xy / self.len() as f32;
    540
    541
            // Cov(X,Y) = E[XY] - μx·μy
    542
    0
            Ok(mean_xy - mean_x * mean_y)
    543
    0
        }
    544
    545
        /// Pearson correlation coefficient
    546
        ///
    547
        /// Computes the Pearson correlation coefficient: ρ(X,Y) = Cov(X,Y) / (σx·σy)
    548
        /// Normalized covariance in range [-1, 1].
    549
        ///
    550
        /// # Performance
    551
        ///
    552
        /// Uses optimized SIMD implementations via covariance() and stddev().
    553
        ///
    554
        /// # Examples
    555
        ///
    556
        /// ```
    557
        /// use trueno::Vector;
    558
        ///
    559
        /// let x = Vector::from_slice(&[1.0, 2.0, 3.0]);
    560
        /// let y = Vector::from_slice(&[2.0, 4.0, 6.0]);
    561
        /// let corr = x.correlation(&y)?;
    562
        /// assert!((corr - 1.0).abs() < 1e-5); // Perfect positive correlation
    563
        /// # Ok::<(), trueno::TruenoError>(())
    564
        /// ```
    565
        ///
    566
        /// # Size mismatch
    567
        ///
    568
        /// Returns an error if vectors have different lengths.
    569
        ///
    570
        /// # Division by zero
    571
        ///
    572
        /// Returns DivisionByZero error if either vector has zero standard deviation
    573
        /// (i.e., is constant).
    574
        ///
    575
        /// ```
    576
        /// use trueno::{Vector, TruenoError};
    577
        ///
    578
        /// let x = Vector::from_slice(&[5.0, 5.0, 5.0]); // Constant
    579
        /// let y = Vector::from_slice(&[1.0, 2.0, 3.0]);
    580
        /// assert!(matches!(x.correlation(&y), Err(TruenoError::DivisionByZero)));
    581
        /// ```
    582
    0
        pub fn correlation(&self, other: &Self) -> Result<f32> {
    583
    0
            let cov = self.covariance(other)?;
    584
    0
            let std_x = self.stddev()?;
    585
    0
            let std_y = other.stddev()?;
    586
    587
            // Check for zero standard deviation (constant vectors)
    588
    0
            if std_x.abs() < 1e-10 || std_y.abs() < 1e-10 {
    589
    0
                return Err(TruenoError::DivisionByZero);
    590
    0
            }
    591
    592
            // ρ(X,Y) = Cov(X,Y) / (σx·σy)
    593
            // Clamp to [-1, 1] to handle floating-point precision errors
    594
    0
            let corr = cov / (std_x * std_y);
    595
    0
            Ok(corr.clamp(-1.0, 1.0))
    596
    0
        }
    597
    }
    598
    599
    #[cfg(test)]
    600
    mod tests {
    601
        use super::*;
    602
        use crate::TruenoError;
    603
    604
        // ========== Basic Reductions ==========
    605
    606
        #[test]
    607
        fn test_dot_basic() {
    608
            let a = Vector::from_slice(&[1.0, 2.0, 3.0]);
    609
            let b = Vector::from_slice(&[4.0, 5.0, 6.0]);
    610
            let result = a.dot(&b).unwrap();
    611
            assert!((result - 32.0).abs() < 1e-6); // 1*4 + 2*5 + 3*6 = 32
    612
        }
    613
    614
        #[test]
    615
        fn test_dot_size_mismatch() {
    616
            let a = Vector::from_slice(&[1.0, 2.0, 3.0]);
    617
            let b = Vector::from_slice(&[4.0, 5.0]);
    618
            assert!(matches!(a.dot(&b), Err(TruenoError::SizeMismatch { .. })));
    619
        }
    620
    621
        #[test]
    622
        fn test_dot_empty() {
    623
            let a = Vector::<f32>::from_slice(&[]);
    624
            let b = Vector::<f32>::from_slice(&[]);
    625
            let result = a.dot(&b).unwrap();
    626
            assert!((result - 0.0).abs() < 1e-6);
    627
        }
    628
    629
        #[test]
    630
        fn test_dot_single() {
    631
            let a = Vector::from_slice(&[3.0]);
    632
            let b = Vector::from_slice(&[4.0]);
    633
            let result = a.dot(&b).unwrap();
    634
            assert!((result - 12.0).abs() < 1e-6);
    635
        }
    636
    637
        #[test]
    638
        fn test_dot_large_aligned() {
    639
            // Test SIMD path with aligned size
    640
            let a = Vector::from_slice(&[1.0; 256]);
    641
            let b = Vector::from_slice(&[2.0; 256]);
    642
            let result = a.dot(&b).unwrap();
    643
            assert!((result - 512.0).abs() < 1e-3); // 256 * 1 * 2 = 512
    644
        }
    645
    646
        #[test]
    647
        fn test_dot_large_unaligned() {
    648
            // Test SIMD path with unaligned size
    649
            let a = Vector::from_slice(&[1.0; 259]);
    650
            let b = Vector::from_slice(&[2.0; 259]);
    651
            let result = a.dot(&b).unwrap();
    652
            assert!((result - 518.0).abs() < 1e-3); // 259 * 1 * 2 = 518
    653
        }
    654
    655
        #[test]
    656
        fn test_sum_basic() {
    657
            let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0]);
    658
            assert!((v.sum().unwrap() - 10.0).abs() < 1e-6);
    659
        }
    660
    661
        #[test]
    662
        fn test_sum_empty() {
    663
            let v = Vector::<f32>::from_slice(&[]);
    664
            assert!((v.sum().unwrap() - 0.0).abs() < 1e-6);
    665
        }
    666
    667
        #[test]
    668
        fn test_sum_single() {
    669
            let v = Vector::from_slice(&[42.0]);
    670
            assert!((v.sum().unwrap() - 42.0).abs() < 1e-6);
    671
        }
    672
    673
        #[test]
    674
        fn test_sum_negatives() {
    675
            let v = Vector::from_slice(&[-1.0, -2.0, 3.0, 4.0]);
    676
            assert!((v.sum().unwrap() - 4.0).abs() < 1e-6);
    677
        }
    678
    679
        #[test]
    680
        fn test_max_basic() {
    681
            let v = Vector::from_slice(&[1.0, 5.0, 3.0, 2.0]);
    682
            assert!((v.max().unwrap() - 5.0).abs() < 1e-6);
    683
        }
    684
    685
        #[test]
    686
        fn test_max_empty() {
    687
            let v = Vector::<f32>::from_slice(&[]);
    688
            assert!(matches!(v.max(), Err(TruenoError::InvalidInput(_))));
    689
        }
    690
    691
        #[test]
    692
        fn test_max_single() {
    693
            let v = Vector::from_slice(&[42.0]);
    694
            assert!((v.max().unwrap() - 42.0).abs() < 1e-6);
    695
        }
    696
    697
        #[test]
    698
        fn test_max_all_negative() {
    699
            let v = Vector::from_slice(&[-5.0, -1.0, -3.0, -2.0]);
    700
            assert!((v.max().unwrap() - (-1.0)).abs() < 1e-6);
    701
        }
    702
    703
        #[test]
    704
        fn test_min_basic() {
    705
            let v = Vector::from_slice(&[1.0, 5.0, 3.0, 2.0]);
    706
            assert!((v.min().unwrap() - 1.0).abs() < 1e-6);
    707
        }
    708
    709
        #[test]
    710
        fn test_min_empty() {
    711
            let v = Vector::<f32>::from_slice(&[]);
    712
            assert!(matches!(v.min(), Err(TruenoError::InvalidInput(_))));
    713
        }
    714
    715
        #[test]
    716
        fn test_min_single() {
    717
            let v = Vector::from_slice(&[42.0]);
    718
            assert!((v.min().unwrap() - 42.0).abs() < 1e-6);
    719
        }
    720
    721
        #[test]
    722
        fn test_min_all_negative() {
    723
            let v = Vector::from_slice(&[-5.0, -1.0, -3.0, -2.0]);
    724
            assert!((v.min().unwrap() - (-5.0)).abs() < 1e-6);
    725
        }
    726
    727
        // ========== Index-finding ==========
    728
    729
        #[test]
    730
        fn test_argmax_basic() {
    731
            let v = Vector::from_slice(&[1.0, 5.0, 3.0, 2.0]);
    732
            assert_eq!(v.argmax().unwrap(), 1);
    733
        }
    734
    735
        #[test]
    736
        fn test_argmax_empty() {
    737
            let v = Vector::<f32>::from_slice(&[]);
    738
            assert!(matches!(v.argmax(), Err(TruenoError::InvalidInput(_))));
    739
        }
    740
    741
        #[test]
    742
        fn test_argmax_single() {
    743
            let v = Vector::from_slice(&[42.0]);
    744
            assert_eq!(v.argmax().unwrap(), 0);
    745
        }
    746
    747
        #[test]
    748
        fn test_argmax_duplicate_max() {
    749
            let v = Vector::from_slice(&[1.0, 5.0, 5.0, 2.0]);
    750
            // Should return first occurrence
    751
            assert_eq!(v.argmax().unwrap(), 1);
    752
        }
    753
    754
        #[test]
    755
        fn test_argmin_basic() {
    756
            let v = Vector::from_slice(&[3.0, 1.0, 5.0, 2.0]);
    757
            assert_eq!(v.argmin().unwrap(), 1);
    758
        }
    759
    760
        #[test]
    761
        fn test_argmin_empty() {
    762
            let v = Vector::<f32>::from_slice(&[]);
    763
            assert!(matches!(v.argmin(), Err(TruenoError::InvalidInput(_))));
    764
        }
    765
    766
        #[test]
    767
        fn test_argmin_single() {
    768
            let v = Vector::from_slice(&[42.0]);
    769
            assert_eq!(v.argmin().unwrap(), 0);
    770
        }
    771
    772
        // ========== Numerically Stable ==========
    773
    774
        #[test]
    775
        fn test_sum_kahan_basic() {
    776
            let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0]);
    777
            assert!((v.sum_kahan().unwrap() - 10.0).abs() < 1e-6);
    778
        }
    779
    780
        #[test]
    781
        fn test_sum_kahan_empty() {
    782
            let v = Vector::<f32>::from_slice(&[]);
    783
            assert!((v.sum_kahan().unwrap() - 0.0).abs() < 1e-6);
    784
        }
    785
    786
        #[test]
    787
        fn test_sum_kahan_precision() {
    788
            // Kahan summation provides better precision for certain scenarios
    789
            // but f32 limits mean 1e10 + 1 = 1e10 in float representation
    790
            // Test with values that demonstrate the benefit of Kahan
    791
            let v = Vector::from_slice(&[1.0, 1e-8, 1e-8, 1e-8, 1e-8]);
    792
            let result = v.sum_kahan().unwrap();
    793
            // Should be close to 1.0 + 4e-8
    794
            assert!((result - 1.00000004).abs() < 1e-6);
    795
        }
    796
    797
        #[test]
    798
        fn test_sum_of_squares_basic() {
    799
            let v = Vector::from_slice(&[1.0, 2.0, 3.0]);
    800
            // 1 + 4 + 9 = 14
    801
            assert!((v.sum_of_squares().unwrap() - 14.0).abs() < 1e-6);
    802
        }
    803
    804
        #[test]
    805
        fn test_sum_of_squares_empty() {
    806
            let v = Vector::<f32>::from_slice(&[]);
    807
            assert!((v.sum_of_squares().unwrap() - 0.0).abs() < 1e-6);
    808
        }
    809
    810
        // ========== Statistical ==========
    811
    812
        #[test]
    813
        fn test_mean_basic() {
    814
            let v = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0, 5.0]);
    815
            assert!((v.mean().unwrap() - 3.0).abs() < 1e-6);
    816
        }
    817
    818
        #[test]
    819
        fn test_mean_empty() {
    820
            let v = Vector::<f32>::from_slice(&[]);
    821
            assert!(matches!(v.mean(), Err(TruenoError::EmptyVector)));
    822
        }
    823
    824
        #[test]
    825
        fn test_mean_single() {
    826
            let v = Vector::from_slice(&[42.0]);
    827
            assert!((v.mean().unwrap() - 42.0).abs() < 1e-6);
    828
        }
    829
    830
        #[test]
    831
        fn test_variance_basic() {
    832
            let v = Vector::from_slice(&[2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0]);
    833
            // Mean = 5, Variance = 4
    834
            let var = v.variance().unwrap();
    835
            assert!((var - 4.0).abs() < 1e-3);
    836
        }
    837
    838
        #[test]
    839
        fn test_variance_empty() {
    840
            let v = Vector::<f32>::from_slice(&[]);
    841
            assert!(matches!(v.variance(), Err(TruenoError::EmptyVector)));
    842
        }
    843
    844
        #[test]
    845
        fn test_variance_constant() {
    846
            let v = Vector::from_slice(&[5.0, 5.0, 5.0, 5.0]);
    847
            assert!((v.variance().unwrap() - 0.0).abs() < 1e-6);
    848
        }
    849
    850
        #[test]
    851
        fn test_stddev_basic() {
    852
            let v = Vector::from_slice(&[2.0, 4.0, 4.0, 4.0, 5.0, 5.0, 7.0, 9.0]);
    853
            // Stddev = sqrt(4) = 2
    854
            let std = v.stddev().unwrap();
    855
            assert!((std - 2.0).abs() < 1e-3);
    856
        }
    857
    858
        #[test]
    859
        fn test_stddev_empty() {
    860
            let v = Vector::<f32>::from_slice(&[]);
    861
            assert!(matches!(v.stddev(), Err(TruenoError::EmptyVector)));
    862
        }
    863
    864
        #[test]
    865
        fn test_covariance_basic() {
    866
            let x = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0, 5.0]);
    867
            let y = Vector::from_slice(&[2.0, 4.0, 6.0, 8.0, 10.0]); // y = 2x
    868
            let cov = x.covariance(&y).unwrap();
    869
            // Cov(X, 2X) = 2 * Var(X) = 2 * 2 = 4
    870
            assert!((cov - 4.0).abs() < 1e-3);
    871
        }
    872
    873
        #[test]
    874
        fn test_covariance_size_mismatch() {
    875
            let x = Vector::from_slice(&[1.0, 2.0, 3.0]);
    876
            let y = Vector::from_slice(&[1.0, 2.0]);
    877
            assert!(matches!(x.covariance(&y), Err(TruenoError::SizeMismatch { .. })));
    878
        }
    879
    880
        #[test]
    881
        fn test_covariance_empty() {
    882
            let x = Vector::<f32>::from_slice(&[]);
    883
            let y = Vector::<f32>::from_slice(&[]);
    884
            assert!(matches!(x.covariance(&y), Err(TruenoError::EmptyVector)));
    885
        }
    886
    887
        #[test]
    888
        fn test_correlation_positive() {
    889
            let x = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0, 5.0]);
    890
            let y = Vector::from_slice(&[2.0, 4.0, 6.0, 8.0, 10.0]); // y = 2x
    891
            let corr = x.correlation(&y).unwrap();
    892
            // Perfect positive correlation
    893
            assert!((corr - 1.0).abs() < 1e-3);
    894
        }
    895
    896
        #[test]
    897
        fn test_correlation_negative() {
    898
            let x = Vector::from_slice(&[1.0, 2.0, 3.0, 4.0, 5.0]);
    899
            let y = Vector::from_slice(&[10.0, 8.0, 6.0, 4.0, 2.0]); // y = -2x + 12
    900
            let corr = x.correlation(&y).unwrap();
    901
            // Perfect negative correlation
    902
            assert!((corr - (-1.0)).abs() < 1e-3);
    903
        }
    904
    905
        #[test]
    906
        fn test_correlation_constant_x() {
    907
            let x = Vector::from_slice(&[5.0, 5.0, 5.0]);
    908
            let y = Vector::from_slice(&[1.0, 2.0, 3.0]);
    909
            assert!(matches!(x.correlation(&y), Err(TruenoError::DivisionByZero)));
    910
        }
    911
    912
        #[test]
    913
        fn test_correlation_constant_y() {
    914
            let x = Vector::from_slice(&[1.0, 2.0, 3.0]);
    915
            let y = Vector::from_slice(&[5.0, 5.0, 5.0]);
    916
            assert!(matches!(x.correlation(&y), Err(TruenoError::DivisionByZero)));
    917
        }
    918
    919
        // ========== Backend Tests ==========
    920
    921
        #[test]
    922
        fn test_dot_scalar_backend() {
    923
            let a = Vector::from_slice_with_backend(&[1.0, 2.0, 3.0], Backend::Scalar);
    924
            let b = Vector::from_slice_with_backend(&[4.0, 5.0, 6.0], Backend::Scalar);
    925
            let result = a.dot(&b).unwrap();
    926
            assert!((result - 32.0).abs() < 1e-6);
    927
        }
    928
    929
        #[test]
    930
        #[cfg(target_arch = "x86_64")]
    931
        fn test_dot_sse2_backend() {
    932
            let a = Vector::from_slice_with_backend(
    933
                &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0],
    934
                Backend::SSE2,
    935
            );
    936
            let b = Vector::from_slice_with_backend(
    937
                &[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
    938
                Backend::SSE2,
    939
            );
    940
            let result = a.dot(&b).unwrap();
    941
            assert!((result - 36.0).abs() < 1e-6); // sum 1..8 = 36
    942
        }
    943
    944
        #[test]
    945
        #[cfg(target_arch = "x86_64")]
    946
        fn test_dot_avx2_backend() {
    947
            if !is_x86_feature_detected!("avx2") {
    948
                return;
    949
            }
    950
            let a = Vector::from_slice_with_backend(&[1.0; 32], Backend::AVX2);
    951
            let b = Vector::from_slice_with_backend(&[2.0; 32], Backend::AVX2);
    952
            let result = a.dot(&b).unwrap();
    953
            assert!((result - 64.0).abs() < 1e-4);
    954
        }
    955
    956
        #[test]
    957
        fn test_sum_scalar_backend() {
    958
            let v = Vector::from_slice_with_backend(&[1.0, 2.0, 3.0, 4.0], Backend::Scalar);
    959
            assert!((v.sum().unwrap() - 10.0).abs() < 1e-6);
    960
        }
    961
    962
        #[test]
    963
        fn test_max_scalar_backend() {
    964
            let v = Vector::from_slice_with_backend(&[1.0, 5.0, 3.0, 2.0], Backend::Scalar);
    965
            assert!((v.max().unwrap() - 5.0).abs() < 1e-6);
    966
        }
    967
    968
        #[test]
    969
        fn test_min_scalar_backend() {
    970
            let v = Vector::from_slice_with_backend(&[1.0, 5.0, 3.0, 2.0], Backend::Scalar);
    971
            assert!((v.min().unwrap() - 1.0).abs() < 1e-6);
    972
        }
    973
    }