rosa 0.1.0

use crate::Matrix;

impl Matrix {
    /// Transpose: (m×n) → (n×m). Returns a new Matrix with copied data.
    pub fn transpose(&self) -> Matrix {
        let mut data = vec![0.0; self.rows() * self.cols()];
        for r in 0..self.rows() {
            for c in 0..self.cols() {
                data[c * self.rows() + r] = self.get(r, c);
            }
        }
        Matrix::from_vec(data, self.cols(), self.rows())
    }

    /// Matrix multiply: `self (m×k) @ other (k×n) → result (m×n)`.
    /// Panics if `self.cols() != other.rows()`.
    pub fn matmul(&self, other: &Matrix) -> Matrix {
        let m = self.rows();
        let k = self.cols();
        let n = other.cols();
        assert_eq!(
            k,
            other.rows(),
            "matmul dimension mismatch: {}x{} @ {}x{}",
            m,
            k,
            other.rows(),
            n
        );

        let mut result = vec![0.0; m * n];
        unsafe {
            matrixmultiply::dgemm(
                m,
                k,
                n,
                1.0,
                self.as_slice().as_ptr(),
                k as isize,
                1,
                other.as_slice().as_ptr(),
                n as isize,
                1,
                0.0,
                result.as_mut_ptr(),
                n as isize,
                1,
            );
        }
        Matrix::from_vec(result, m, n)
    }

    /// Matrix-vector multiply: `self (m×k) @ vec (k,) → result (m,)`.
    /// Panics if `self.cols() != vec.len()`.
    pub fn vecmul(&self, vec: &[f64]) -> Vec<f64> {
        let m = self.rows();
        let k = self.cols();
        assert_eq!(
            k,
            vec.len(),
            "vecmul dimension mismatch: {}x{} @ {}",
            m,
            k,
            vec.len()
        );

        let mut result = vec![0.0; m];
        // Use dgemm with n=1 for vector multiply
        unsafe {
            matrixmultiply::dgemm(
                m,
                k,
                1,
                1.0,
                self.as_slice().as_ptr(),
                k as isize,
                1,
                vec.as_ptr(),
                1,
                1,
                0.0,
                result.as_mut_ptr(),
                1,
                1,
            );
        }
        result
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use approx::assert_relative_eq;

    #[test]
    fn test_transpose() {
        let m = Matrix::from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]);
        let t = m.transpose();
        assert_eq!(t.shape(), (3, 2));
        assert_eq!(t.row(0), &[1.0, 4.0]);
        assert_eq!(t.row(1), &[2.0, 5.0]);
        assert_eq!(t.row(2), &[3.0, 6.0]);
    }

    #[test]
    fn test_transpose_roundtrip() {
        let m = Matrix::from_fn(3, 4, |r, c| (r * 4 + c) as f64);
        let tt = m.transpose().transpose();
        assert_eq!(m, tt);
    }

    #[test]
    fn test_transpose_shape_swap() {
        let m = Matrix::zeros(3, 7);
        let t = m.transpose();
        assert_eq!(t.shape(), (7, 3));
    }

    #[test]
    fn test_matmul_identity() {
        let a = Matrix::from_rows(&[&[1.0, 2.0], &[3.0, 4.0]]);
        let eye = Matrix::from_fn(2, 2, |r, c| if r == c { 1.0 } else { 0.0 });
        let result = a.matmul(&eye);
        assert_eq!(result, a);
    }

    #[test]
    fn test_matmul_2x3_3x2() {
        let a = Matrix::from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]);
        let b = Matrix::from_rows(&[&[7.0, 8.0], &[9.0, 10.0], &[11.0, 12.0]]);
        let c = a.matmul(&b);
        assert_eq!(c.shape(), (2, 2));
        // [1*7+2*9+3*11, 1*8+2*10+3*12] = [58, 64]
        // [4*7+5*9+6*11, 4*8+5*10+6*12] = [139, 154]
        assert_relative_eq!(c.get(0, 0), 58.0);
        assert_relative_eq!(c.get(0, 1), 64.0);
        assert_relative_eq!(c.get(1, 0), 139.0);
        assert_relative_eq!(c.get(1, 1), 154.0);
    }

    #[test]
    #[should_panic(expected = "matmul dimension mismatch")]
    fn test_matmul_dimension_mismatch() {
        let a = Matrix::zeros(2, 3);
        let b = Matrix::zeros(4, 2);
        a.matmul(&b);
    }

    #[test]
    fn test_matmul_transpose_property() {
        // (A @ B)^T == B^T @ A^T
        let a = Matrix::from_fn(2, 3, |r, c| (r * 3 + c + 1) as f64);
        let b = Matrix::from_fn(3, 4, |r, c| (r * 4 + c + 1) as f64);
        let ab_t = a.matmul(&b).transpose();
        let bt_at = b.transpose().matmul(&a.transpose());
        for r in 0..ab_t.rows() {
            for c in 0..ab_t.cols() {
                assert_relative_eq!(ab_t.get(r, c), bt_at.get(r, c), epsilon = 1e-10);
            }
        }
    }

    #[test]
    fn test_vecmul() {
        let m = Matrix::from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]);
        let v = vec![1.0, 0.0, 1.0];
        let result = m.vecmul(&v);
        assert_eq!(result.len(), 2);
        assert_relative_eq!(result[0], 4.0); // 1*1 + 2*0 + 3*1
        assert_relative_eq!(result[1], 10.0); // 4*1 + 5*0 + 6*1
    }

    #[test]
    fn test_vecmul_filterbank() {
        // Simulate filterbank-like operation: (128×1025) @ (1025,) → (128,)
        let filterbank = Matrix::from_fn(128, 1025, |r, c| {
            if c >= r * 8 && c < (r + 1) * 8 {
                1.0
            } else {
                0.0
            }
        });
        let spectrum: Vec<f64> = (0..1025).map(|i| i as f64).collect();
        let result = filterbank.vecmul(&spectrum);
        assert_eq!(result.len(), 128);
        // First mel bin sums spectrum[0..8] = 0+1+2+3+4+5+6+7 = 28
        assert_relative_eq!(result[0], 28.0);
    }

    #[test]
    #[should_panic(expected = "vecmul dimension mismatch")]
    fn test_vecmul_dimension_mismatch() {
        let m = Matrix::zeros(2, 3);
        let v = vec![1.0, 2.0];
        m.vecmul(&v);
    }
}