Merge 4ddacdfd21b6fe80d090491a10aed1a96e0b7383 into 11330e464ba3a7f08aaf73bc918281472c503b1d

Cargo.toml

@@ -14,8 +14,6 @@ crate-type = ["cdylib", "lib"]
 [dependencies]
 chrono = "^0.4.10"
 criterion = { version = "0.5", features = ["html_reports"], optional = true }
-
-[dev-dependencies]
 rand = "^0.9.1"
 
 [features]

src/compute/activations.rs

@@ -0,0 +1,60 @@
+use crate::matrix::{Matrix, SeriesOps};
+
+pub fn sigmoid(x: &Matrix<f64>) -> Matrix<f64> {
+    x.map(|v| 1.0 / (1.0 + (-v).exp()))
+}
+
+pub fn dsigmoid(y: &Matrix<f64>) -> Matrix<f64> {
+    // derivative w.r.t. pre-activation; takes y = sigmoid(x)
+    y.map(|v| v * (1.0 - v))
+}
+
+pub fn relu(x: &Matrix<f64>) -> Matrix<f64> {
+    x.map(|v| if v > 0.0 { v } else { 0.0 })
+}
+
+pub fn drelu(x: &Matrix<f64>) -> Matrix<f64> {
+    x.map(|v| if v > 0.0 { 1.0 } else { 0.0 })
+}
+
+
+mod tests {
+    use super::*;
+    
+    // Helper function to round all elements in a matrix to n decimal places
+    fn _round_matrix(mat: &Matrix<f64>, decimals: u32) -> Matrix<f64> {
+        let factor = 10f64.powi(decimals as i32);
+        let rounded: Vec<f64> = mat.to_vec().iter().map(|v| (v * factor).round() / factor).collect();
+        Matrix::from_vec(rounded, mat.rows(), mat.cols())
+    }
+
+    #[test]
+    fn test_sigmoid() {
+        let x = Matrix::from_vec(vec![-1.0, 0.0, 1.0], 3, 1);
+        let expected = Matrix::from_vec(vec![0.26894142, 0.5, 0.73105858], 3, 1);
+        let result = sigmoid(&x);
+        assert_eq!(_round_matrix(&result, 6), _round_matrix(&expected, 6));
+    }
+
+    #[test]
+    fn test_relu() {
+        let x = Matrix::from_vec(vec![-1.0, 0.0, 1.0], 3, 1);
+        let expected = Matrix::from_vec(vec![0.0, 0.0, 1.0], 3, 1);
+        assert_eq!(relu(&x), expected);
+    }
+
+    #[test]
+    fn test_dsigmoid() {
+        let y = Matrix::from_vec(vec![0.26894142, 0.5, 0.73105858], 3, 1);
+        let expected = Matrix::from_vec(vec![0.19661193, 0.25, 0.19661193], 3, 1);
+        let result = dsigmoid(&y);
+        assert_eq!(_round_matrix(&result, 6), _round_matrix(&expected, 6));
+    }
+    
+    #[test]
+    fn test_drelu() {
+        let x = Matrix::from_vec(vec![-1.0, 0.0, 1.0], 3, 1);
+        let expected = Matrix::from_vec(vec![0.0, 0.0, 1.0], 3, 1);
+        assert_eq!(drelu(&x), expected);    
+    }
+}

src/compute/mod.rs

@@ -0,0 +1,3 @@
+pub mod activations;
+
+pub mod models;

src/compute/models/dense_nn.rs

@@ -0,0 +1,65 @@
+use crate::matrix::{Matrix, SeriesOps};
+use crate::compute::activations::{relu, sigmoid, drelu};  
+use rand::Rng;
+
+pub struct DenseNN {
+    w1: Matrix<f64>, // (n_in, n_hidden)
+    b1: Matrix<f64>, // (1,    n_hidden)
+    w2: Matrix<f64>, // (n_hidden, n_out)
+    b2: Matrix<f64>, // (1,    n_out)
+}
+
+impl DenseNN {
+    pub fn new(n_in: usize, n_hidden: usize, n_out: usize) -> Self {
+        let mut rng = rand::rng();
+        let mut init = |rows, cols| {
+            let data = (0..rows * cols)
+                .map(|_| rng.random_range(-1.0..1.0))
+                .collect::<Vec<_>>();
+            Matrix::from_vec(data, rows, cols)
+        };
+        Self {
+            w1: init(n_in, n_hidden),
+            b1: Matrix::zeros(1, n_hidden),
+            w2: init(n_hidden, n_out),
+            b2: Matrix::zeros(1, n_out),
+        }
+    }
+
+    pub fn forward(&self, x: &Matrix<f64>) -> (Matrix<f64>, Matrix<f64>, Matrix<f64>) {
+        // z1 = X·W1 + b1 ; a1 = ReLU(z1)
+        let z1 = x.dot(&self.w1) + &self.b1;
+        let a1 = relu(&z1);
+        // z2 = a1·W2 + b2 ; a2 = softmax(z2) (here binary => sigmoid)
+        let z2 = a1.dot(&self.w2) + &self.b2;
+        let a2 = sigmoid(&z2); // binary output
+        (a1, z2, a2)           // keep intermediates for back-prop
+    }
+
+    pub fn train(&mut self, x: &Matrix<f64>, y: &Matrix<f64>, lr: f64, epochs: usize) {
+        let m = x.rows() as f64;
+        for _ in 0..epochs {
+            let (a1, _z2, y_hat) = self.forward(x);
+
+            // -------- backwards ----------
+            // dL/da2 = y_hat - y  (BCE derivative)
+            let dz2 = &y_hat - y;                           // (m, n_out)
+            let dw2 = a1.transpose().dot(&dz2) / m;         // (n_h, n_out)
+            // let db2 = dz2.sum_vertical() * (1.0 / m);       // broadcast ok
+            let db2 = Matrix::from_vec(dz2.sum_vertical(), 1, dz2.cols()) * (1.0 / m); // (1, n_out)
+            let da1 = dz2.dot(&self.w2.transpose());        // (m,n_h)
+            let dz1 = da1.zip(&a1, |g, act| g * drelu(&Matrix::from_cols(vec![vec![act]])).data()[0]); // (m,n_h)
+            
+            // real code: drelu returns Matrix, broadcasting needed; you can optimise.
+
+            let dw1 = x.transpose().dot(&dz1) / m;          // (n_in,n_h)
+            let db1 = Matrix::from_vec(dz1.sum_vertical(), 1, dz1.cols()) * (1.0 / m); // (1, n_h)
+
+            // -------- update ----------
+            self.w2 = &self.w2 - &(dw2 * lr);
+            self.b2 = &self.b2 - &(db2 * lr);
+            self.w1 = &self.w1 - &(dw1 * lr);
+            self.b1 = &self.b1 - &(db1 * lr);
+        }
+    }
+}

src/compute/models/gaussian_nb.rs

@@ -0,0 +1,124 @@
+use crate::matrix::{Matrix};
+use std::collections::HashMap;
+
+pub struct GaussianNB {
+    classes: Vec<f64>, // distinct labels
+    priors: Vec<f64>,  // P(class)
+    means: Vec<Matrix<f64>>,
+    variances: Vec<Matrix<f64>>,
+    eps: f64, // var-smoothing
+}
+
+impl GaussianNB {
+    pub fn new(var_smoothing: f64) -> Self {
+        Self {
+            classes: vec![],
+            priors: vec![],
+            means: vec![],
+            variances: vec![],
+            eps: var_smoothing,
+        }
+    }
+
+    pub fn fit(&mut self, x: &Matrix<f64>, y: &Matrix<f64>) {
+        let m = x.rows();
+        let n = x.cols();
+        assert_eq!(y.rows(), m);
+        assert_eq!(y.cols(), 1);
+        if m == 0 || n == 0 {
+            panic!("Input matrix x or y is empty");
+        }
+
+        // ----- group samples by label -----
+        let mut groups: HashMap<u64, Vec<usize>> = HashMap::new();
+        for i in 0..m {
+            let label_bits = y[(i, 0)].to_bits();
+            groups.entry(label_bits).or_default().push(i);
+        }
+        if groups.is_empty() {
+            panic!("No class labels found in y");
+        }
+
+        self.classes = groups
+            .keys()
+            .cloned()
+            .map(f64::from_bits)
+            .collect::<Vec<_>>();
+        // Note: If NaN is present in class labels, this may panic. Ensure labels are valid floats.
+        self.classes.sort_by(|a, b| a.partial_cmp(b).unwrap());
+
+        self.priors.clear();
+        self.means.clear();
+        self.variances.clear();
+
+        for &c in &self.classes {
+            let label_bits = c.to_bits();
+            let idx = &groups[&label_bits];
+            let count = idx.len();
+            if count == 0 {
+                panic!("Class group for label {c} is empty");
+            }
+            self.priors.push(count as f64 / m as f64);
+
+            let mut mean = Matrix::zeros(1, n);
+            let mut var = Matrix::zeros(1, n);
+
+            // mean
+            for &i in idx {
+                for j in 0..n {
+                    mean[(0, j)] += x[(i, j)];
+                }
+            }
+            for j in 0..n {
+                mean[(0, j)] /= count as f64;
+            }
+
+            // variance
+            for &i in idx {
+                for j in 0..n {
+                    let d = x[(i, j)] - mean[(0, j)];
+                    var[(0, j)] += d * d;
+                }
+            }
+            for j in 0..n {
+                var[(0, j)] = var[(0, j)] / count as f64 + self.eps; // always add eps after division
+                if var[(0, j)] <= 0.0 {
+                    var[(0, j)] = self.eps; // ensure strictly positive variance
+                }
+            }
+
+            self.means.push(mean);
+            self.variances.push(var);
+        }
+    }
+
+    /// Return class labels (shape m×1) for samples in X.
+    pub fn predict(&self, x: &Matrix<f64>) -> Matrix<f64> {
+        let m = x.rows();
+        let k = self.classes.len();
+        let n = x.cols();
+        let mut preds = Matrix::zeros(m, 1);
+        let ln_2pi = (2.0 * std::f64::consts::PI).ln();
+
+        for i in 0..m {
+            let mut best_class = 0usize;
+            let mut best_log_prob = f64::NEG_INFINITY;
+            for c in 0..k {
+                // log P(y=c) + Σ log N(x_j | μ, σ²)
+                let mut log_prob = self.priors[c].ln();
+                for j in 0..n {
+                    let mean = self.means[c][(0, j)];
+                    let var = self.variances[c][(0, j)];
+                    let diff = x[(i, j)] - mean;
+                    log_prob += -0.5 * (diff * diff / var + var.ln() + ln_2pi);
+                }
+                if log_prob > best_log_prob {
+                    best_log_prob = log_prob;
+                    best_class = c;
+                }
+            }
+            preds[(i, 0)] = self.classes[best_class];
+        }
+        preds
+    }
+}

src/compute/models/k_means.rs

@@ -0,0 +1,105 @@
+use crate::matrix::Matrix;
+use std::collections::HashMap;
+
+pub struct GaussianNB {
+    classes:      Vec<f64>,        // distinct labels
+    priors:       Vec<f64>,        // P(class)
+    means:        Vec<Matrix<f64>>,
+    variances:    Vec<Matrix<f64>>,
+    eps:          f64,             // var-smoothing
+}
+
+impl GaussianNB {
+    pub fn new(var_smoothing: f64) -> Self {
+        Self {
+            classes: vec![],
+            priors: vec![],
+            means: vec![],
+            variances: vec![],
+            eps: var_smoothing,
+        }
+    }
+
+    pub fn fit(&mut self, x: &Matrix<f64>, y: &Matrix<f64>) {
+        let m = x.rows();
+        let n = x.cols();
+        assert_eq!(y.rows(), m);
+        assert_eq!(y.cols(), 1);
+
+        // ----- group samples by label -----
+        let mut groups: HashMap<i64, Vec<usize>> = HashMap::new();
+        for i in 0..m {
+            groups.entry(y[(i, 0)] as i64).or_default().push(i);
+        }
+
+        self.classes = groups.keys().cloned().map(|v| v as f64).collect::<Vec<_>>();
+        self.classes.sort_by(|a, b| a.partial_cmp(b).unwrap());
+
+        self.priors.clear();
+        self.means.clear();
+        self.variances.clear();
+
+        for &c in &self.classes {
+            let idx = &groups[&(c as i64)];
+            let count = idx.len();
+            self.priors.push(count as f64 / m as f64);
+
+            let mut mean = Matrix::zeros(1, n);
+            let mut var  = Matrix::zeros(1, n);
+
+            // mean
+            for &i in idx {
+                for j in 0..n {
+                    mean[(0, j)] += x[(i, j)];
+                }
+            }
+            for j in 0..n {
+                mean[(0, j)] /= count as f64;
+            }
+
+            // variance
+            for &i in idx {
+                for j in 0..n {
+                    let d = x[(i, j)] - mean[(0, j)];
+                    var[(0, j)] += d * d;
+                }
+            }
+            for j in 0..n {
+                var[(0, j)] = var[(0, j)] / count as f64 + self.eps;
+            }
+
+            self.means.push(mean);
+            self.variances.push(var);
+        }
+    }
+
+    /// Return class labels (shape m×1) for samples in X.
+    pub fn predict(&self, x: &Matrix<f64>) -> Matrix<f64> {
+        let m = x.rows();
+        let k = self.classes.len();
+        let n = x.cols();
+        let mut preds = Matrix::zeros(m, 1);
+        let ln_2pi = (2.0 * std::f64::consts::PI).ln();
+
+        for i in 0..m {
+            let mut best_class = 0usize;
+            let mut best_log_prob = f64::NEG_INFINITY;
+            for c in 0..k {
+                // log P(y=c) + Σ log N(x_j | μ, σ²)
+                let mut log_prob = self.priors[c].ln();
+                for j in 0..n {
+                    let mean = self.means[c][(0, j)];
+                    let var  = self.variances[c][(0, j)];
+                    let diff = x[(i, j)] - mean;
+                    log_prob += -0.5 * (diff * diff / var + var.ln() + ln_2pi);
+                }
+                if log_prob > best_log_prob {
+                    best_log_prob = log_prob;
+                    best_class = c;
+                }
+            }
+            preds[(i, 0)] = self.classes[best_class];
+        }
+        preds
+    }
+}

src/compute/models/linreg.rs

@@ -0,0 +1,54 @@
+use crate::matrix::{Matrix, SeriesOps};
+
+pub struct LinReg {
+    w: Matrix<f64>, // shape (n_features, 1)
+    b: f64,
+}
+
+impl LinReg {
+    pub fn new(n_features: usize) -> Self {
+        Self {
+            w: Matrix::from_vec(vec![0.0; n_features], n_features, 1),
+            b: 0.0,
+        }
+    }
+
+    pub fn predict(&self, x: &Matrix<f64>) -> Matrix<f64> {
+        // X.dot(w) + b
+        x.dot(&self.w) + self.b
+    }
+
+    pub fn fit(&mut self, x: &Matrix<f64>, y: &Matrix<f64>, lr: f64, epochs: usize) {
+        let m = x.rows() as f64;
+        for _ in 0..epochs {
+            let y_hat = self.predict(x);
+            let err = &y_hat - y; // shape (m,1)
+
+            // grads
+            let grad_w = x.transpose().dot(&err) * (2.0 / m); // (n,1)
+            let grad_b = (2.0 / m) * err.sum_vertical().iter().sum::<f64>();
+            // update
+            self.w = &self.w - &(grad_w * lr);
+            self.b -= lr * grad_b;
+        }
+    }
+}
+
+mod tests {
+
+    use super::LinReg;
+    use crate::matrix::{Matrix};
+
+    #[test]
+    fn test_linreg_fit_predict() {
+        let x = Matrix::from_vec(vec![1.0, 2.0, 3.0, 4.0], 4, 1);
+        let y = Matrix::from_vec(vec![2.0, 3.0, 4.0, 5.0], 4, 1);
+        let mut model = LinReg::new(1);
+        model.fit(&x, &y, 0.01, 10000);
+        let preds = model.predict(&x);
+        assert!((preds[(0, 0)] - 2.0).abs() < 1e-2);
+        assert!((preds[(1, 0)] - 3.0).abs() < 1e-2);
+        assert!((preds[(2, 0)] - 4.0).abs() < 1e-2);
+        assert!((preds[(3, 0)] - 5.0).abs() < 1e-2);
+    }
+}

src/compute/models/logreg.rs

@@ -0,0 +1,36 @@
+use crate::matrix::{Matrix, SeriesOps};
+use crate::compute::activations::sigmoid;
+
+pub struct LogReg {
+    w: Matrix<f64>,
+    b: f64,
+}
+
+impl LogReg {
+    pub fn new(n_features: usize) -> Self {
+        Self {
+            w: Matrix::zeros(n_features, 1),
+            b: 0.0,
+        }
+    }
+
+    pub fn predict_proba(&self, x: &Matrix<f64>) -> Matrix<f64> {
+        sigmoid(&(x.dot(&self.w) + self.b))          // σ(Xw + b)
+    }
+
+    pub fn fit(&mut self, x: &Matrix<f64>, y: &Matrix<f64>, lr: f64, epochs: usize) {
+        let m = x.rows() as f64;
+        for _ in 0..epochs {
+            let p  = self.predict_proba(x);          // shape (m,1)
+            let err = &p - y;                        // derivative of BCE wrt pre-sigmoid
+            let grad_w = x.transpose().dot(&err) / m;
+            let grad_b = err.sum_vertical().iter().sum::<f64>() / m;
+            self.w = &self.w - &(grad_w * lr);
+            self.b -= lr * grad_b;
+        }
+    }
+
+    pub fn predict(&self, x: &Matrix<f64>) -> Matrix<f64> {
+        self.predict_proba(x).map(|p| if p >= 0.5 { 1.0 } else { 0.0 })
+    }
+}

src/compute/models/mod.rs

@@ -0,0 +1,6 @@
+pub mod linreg;
+pub mod logreg;
+pub mod dense_nn;
+pub mod k_means;    
+pub mod pca;
+pub mod gaussian_nb;

src/compute/models/pca.rs

@@ -0,0 +1,85 @@
+use crate::matrix::{Matrix, SeriesOps};
+use rand;
+
+/// Returns the `n_components` principal axes (rows) and the centred data’s mean.
+pub struct PCA {
+    pub components: Matrix<f64>,   // (n_components, n_features)
+    pub mean:       Matrix<f64>,   // (1, n_features)
+}
+
+impl PCA {
+    pub fn fit(x: &Matrix<f64>, n_components: usize, iters: usize) -> Self {
+        let m = x.rows();
+        let n = x.cols();
+        assert!(n_components <= n);
+
+        // ----- centre data -----
+        let mean_vec = {
+            let mut v = Matrix::zeros(1, n);
+            for j in 0..n {
+                let mut s = 0.0;
+                for i in 0..m {
+                    s += x[(i, j)];
+                }
+                v[(0, j)] = s / m as f64;
+            }
+            v
+        };
+        let x_centered = x - &mean_vec;
+
+        // ----- covariance matrix C = Xᵀ·X / (m-1) -----
+        let cov = x_centered.transpose().dot(&x_centered) * (1.0 / (m as f64 - 1.0));
+
+        // ----- power iteration to find top eigenvectors -----
+        let mut comp = Matrix::zeros(n_components, n);
+        let mut b = Matrix::zeros(1, n);            // current vector
+        for c in 0..n_components {
+            // random initial vector
+            for j in 0..n {
+                b[(0, j)] = rand::random::<f64>() - 0.5;
+            }
+            // subtract projections on previously found components
+            for prev in 0..c {
+                // let proj = b.dot(Matrix::from_vec(data, rows, cols).transpose())[(0, 0)];
+                // let proj = b.dot(&comp.row(prev).transpose())[(0, 0)];
+                let proj = b.dot(&Matrix::from_vec(comp.row(prev).to_vec(), 1, n).transpose())[(0, 0)];
+                // subtract projection to maintain orthogonality
+                for j in 0..n {
+                    b[(0, j)] -= proj * comp[(prev, j)];
+                }
+            }
+            // iterate
+            for _ in 0..iters {
+                // b = C·bᵀ
+                let mut nb = cov.dot(&b.transpose()).transpose();
+                // subtract projections again to maintain orthogonality
+                for prev in 0..c {
+                    let proj = nb.dot(&Matrix::from_vec(comp.row(prev).to_vec(), 1, n).transpose())[(0, 0)];
+                    for j in 0..n {
+                        nb[(0, j)] -= proj * comp[(prev, j)];
+                    }
+                }
+                // normalise
+                let norm = nb.data().iter().map(|v| v * v).sum::<f64>().sqrt();
+                for j in 0..n {
+                    nb[(0, j)] /= norm;
+                }
+                b = nb;
+            }
+            // store component
+            for j in 0..n {
+                comp[(c, j)] = b[(0, j)];
+            }
+        }
+        Self {
+            components: comp,
+            mean: mean_vec,
+        }
+    }
+
+    /// Project new data on the learned axes.
+    pub fn transform(&self, x: &Matrix<f64>) -> Matrix<f64> {
+        let x_centered = x - &self.mean;
+        x_centered.dot(&self.components.transpose())
+    }
+}

src/lib.rs

@@ -8,3 +8,6 @@ pub mod frame;
 
 /// Documentation for the [`crate::utils`] module.
 pub mod utils;
+
+/// Documentation for the [`crate::compute`] module.
+pub mod compute;

src/matrix/mat.rs

@@ -179,6 +179,21 @@ impl<T: Clone> Matrix<T> {
         self.cols -= 1;
     }
 
+    #[inline]
+    pub fn row(&self, r: usize) -> Vec<T> {
+        assert!(
+            r < self.rows,
+            "row index {} out of bounds for {} rows",
+            r,
+            self.rows
+        );
+        let mut row_data = Vec::with_capacity(self.cols);
+        for c in 0..self.cols {
+            row_data.push(self[(r, c)].clone()); // Clone each element
+        }
+        row_data
+    }
+
     /// Deletes a row from the matrix. Panics on out-of-bounds.
     /// This is O(N) where N is the number of elements, as it rebuilds the data vec.
     pub fn delete_row(&mut self, row: usize) {
@@ -310,6 +325,26 @@ impl<T: Clone> Matrix<T> {
     }
 }
 
+impl Matrix<f64> {
+    /// Creates a new matrix filled with a specific value of the specified size.
+    pub fn filled(rows: usize, cols: usize, value: f64) -> Self {
+        Matrix {
+            rows,
+            cols,
+            data: vec![value; rows * cols], // Fill with the specified value
+        }
+    }
+    /// Creates a new matrix filled with zeros of the specified size.
+    pub fn zeros(rows: usize, cols: usize) -> Self {
+        Matrix::filled(rows, cols, 0.0)
+    }
+
+    /// Creates a new matrix filled with ones of the specified size.
+    pub fn ones(rows: usize, cols: usize) -> Self {
+        Matrix::filled(rows, cols, 1.0)
+    }
+}
+
 impl<T> Index<(usize, usize)> for Matrix<T> {
     type Output = T;
 
@@ -1110,6 +1145,21 @@ mod tests {
         matrix[(0, 3)] = 99;
     }
 
+    #[test]
+    fn test_row() {
+        let ma = static_test_matrix();
+        assert_eq!(ma.row(0), &[1, 4, 7]);
+        assert_eq!(ma.row(1), &[2, 5, 8]);
+        assert_eq!(ma.row(2), &[3, 6, 9]);
+    }
+
+    #[test]
+    #[should_panic(expected = "row index 3 out of bounds for 3 rows")]
+    fn test_row_out_of_bounds() {
+        let ma = static_test_matrix();
+        ma.row(3);
+    }
+
     #[test]
     fn test_column() {
         let matrix = static_test_matrix_2x4();
@@ -1794,4 +1844,25 @@ mod tests {
             }
         }
     }
+
+    #[test]
+    fn test_matrix_zeros_ones_filled() {
+        // Test zeros
+        let m = Matrix::<f64>::zeros(2, 3);
+        assert_eq!(m.rows(), 2);
+        assert_eq!(m.cols(), 3);
+        assert_eq!(m.data(), &[0.0, 0.0, 0.0, 0.0, 0.0, 0.0]);
+
+        // Test ones
+        let m = Matrix::<f64>::ones(3, 2);
+        assert_eq!(m.rows(), 3);
+        assert_eq!(m.cols(), 2);
+        assert_eq!(m.data(), &[1.0, 1.0, 1.0, 1.0, 1.0, 1.0]);
+
+        // Test filled
+        let m = Matrix::<f64>::filled(2, 2, 42.5);
+        assert_eq!(m.rows(), 2);
+        assert_eq!(m.cols(), 2);
+        assert_eq!(m.data(), &[42.5, 42.5, 42.5, 42.5]);
+    }
 }