Matrices

use matrix.qweruse matrix_solve, matrix_inverse from matrix      // Quillo 0.6: only these

nalgebra = { version = "0.33", default-features = false, features = ["std"] }

// A bakery's numbers: what each loaf costs from three weeks of totals, the inverse of a transform,// and a trend line through monthly sales.////   quillo run matrix.qwe      then open http://localhost:8080/costsuse "matrix.qwer"// Three weeks of flour, butter and sugar used for bread, croissants and cakes (kg each),// and what they cost in total: the price of each ingredient per kg.route GET /costs {    let used = [[20, 2, 1], [10, 6, 3], [5, 3, 6]]    let paid = [41.5, 49.5, 36.0]    let p = matrix_solve(used, paid)?    { flour: p[0], butter: p[1], sugar: p[2] }}route GET /inverse = {    inverse: matrix_inverse([[4, 7], [2, 6]])?,    det: matrix_determinant([[4, 7], [2, 6]])?,    check: matrix_multiply([[4, 7], [2, 6]], matrix_inverse([[4, 7], [2, 6]])?)?,    singular: matrix_inverse([[1, 2], [2, 4]]).error,}route GET /trend {    let sales = [120, 135, 149, 160, 178, 190]    let fit = matrix_least_squares(range(len(sales)) | map((m) => [float(m), 1.0]), sales)?    { per_month: round(fit[0], 2), start: round(fit[1], 2), month_7: round(fit[0] * 6 + fit[1], 1) }}

//! Matrices as lists of rows: linear systems solved, inverses, determinants, products and transposes, and least-squares fits (a trend line through points). For pricing models, physics, statistics and graphics.//! version: 1.0.0//!//! [dependencies]//! nalgebra = { version = "0.33", default-features = false, features = ["std"] }use nalgebra::{DMatrix, DVector};const MAX: usize = 500;fn mat(rows: &[Vec<f64>], what: &str) -> Result<DMatrix<f64>, String> {    let n = rows.len();    let m = rows.first().map(Vec::len).unwrap_or(0);    if n == 0 || m == 0 {        return Err(format!("{} is empty", what));    }    if n > MAX || m > MAX {        return Err(format!("{} is at most {} by {}", what, MAX, MAX));    }    if let Some(i) = rows.iter().position(|r| r.len() != m) {        return Err(format!("row {} of {} has {} numbers, not {}", i, what, rows[i].len(), m));    }    if rows.iter().flatten().any(|x| !x.is_finite()) {        return Err(format!("{} holds a number that is not finite", what));    }    Ok(DMatrix::from_fn(n, m, |i, j| rows[i][j]))}/// Rounds away float dust: 1e-12 becomes 0, 2.0000000000004 becomes 2.fn clean(x: f64) -> f64 {    let y = (x * 1e10).round() / 1e10;    if y == 0.0 { 0.0 } else { y }}fn rows(m: &DMatrix<f64>) -> Vec<Vec<f64>> {    (0..m.nrows()).map(|i| (0..m.ncols()).map(|j| clean(m[(i, j)])).collect()).collect()}/// The `x` with `a · x = b`: `a` a square matrix (a list of rows), `b` a list as long as `a`. An error when `a` is singular (no single answer).pub fn matrix_solve(a: Vec<Vec<f64>>, b: Vec<f64>) -> Result<Vec<f64>, String> {    let a = mat(&a, "a")?;    if a.nrows() != a.ncols() {        return Err("a must be square: use matrix_least_squares for more equations than unknowns".into());    }    if b.len() != a.nrows() {        return Err(format!("b has {} numbers; a has {} rows", b.len(), a.nrows()));    }    let x = a.lu().solve(&DVector::from_vec(b)).ok_or("a is singular: the equations have no single solution")?;    Ok(x.iter().map(|v| clean(*v)).collect())}/// The inverse of a square matrix, or an error when it has none.pub fn matrix_inverse(a: Vec<Vec<f64>>) -> Result<Vec<Vec<f64>>, String> {    let a = mat(&a, "a")?;    if a.nrows() != a.ncols() {        return Err("only a square matrix has an inverse".into());    }    Ok(rows(&a.try_inverse().ok_or("the matrix is singular: it has no inverse")?))}/// The determinant of a square matrix (0 for a singular one).pub fn matrix_determinant(a: Vec<Vec<f64>>) -> Result<f64, String> {    let a = mat(&a, "a")?;    if a.nrows() != a.ncols() {        return Err("only a square matrix has a determinant".into());    }    Ok(clean(a.determinant()))}/// The product `a · b` (as many columns in `a` as rows in `b`).pub fn matrix_multiply(a: Vec<Vec<f64>>, b: Vec<Vec<f64>>) -> Result<Vec<Vec<f64>>, String> {    let (a, b) = (mat(&a, "a")?, mat(&b, "b")?);    if a.ncols() != b.nrows() {        return Err(format!("a has {} columns but b has {} rows", a.ncols(), b.nrows()));    }    Ok(rows(&(a * b)))}/// The matrix turned on its side: rows become columns.pub fn matrix_transpose(a: Vec<Vec<f64>>) -> Result<Vec<Vec<f64>>, String> {    Ok(rows(&mat(&a, "a")?.transpose()))}/// The `x` that makes `a · x` as close to `b` as can be (least squares), for more equations than unknowns: with rows `[x, 1]` it fits a line `y = slope·x + intercept`.pub fn matrix_least_squares(a: Vec<Vec<f64>>, b: Vec<f64>) -> Result<Vec<f64>, String> {    let a = mat(&a, "a")?;    if b.len() != a.nrows() {        return Err(format!("b has {} numbers; a has {} rows", b.len(), a.nrows()));    }    let x = a.svd(true, true).solve(&DVector::from_vec(b), 1e-12)?;    Ok(x.iter().map(|v| clean(*v)).collect())}