diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 7282ed4..ecf3cbb 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -1,22 +1,29 @@ -name: tests +name: ci on: push: - branches: [ master ] + branches: [master] pull_request: - branches: [ master ] + branches: [master] env: CARGO_TERM_COLOR: always jobs: - build: - + test: runs-on: ubuntu-latest - + strategy: + matrix: + features: ["", "rand,nalgebra", "io,rand,nalgebra,metadata"] steps: - - uses: actions/checkout@v2 - - name: Build - run: cargo build --verbose --all-features - - name: Run tests - run: cargo test --verbose --all-features \ No newline at end of file + - uses: actions/checkout@v4 + - uses: actions-rs/toolchain@v1 + with: + toolchain: stable + override: true + components: rustfmt + - uses: Swatinem/rust-cache@v2 + - name: Format + run: cargo fmt --all -- --check + - name: Test + run: cargo test --no-default-features --features "${{ matrix.features }}" diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..77c9227 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,26 @@ +# Changelog + +## 0.4.0 — Unreleased + +### Numerical workflows + +- Add `row`, `col`, and `mat` constructors over ordinary Rhai arrays with explicit orientation and numeric shape validation. +- Add `dot` as a real scalar inner product of equal-length lists, rows, or columns. Use `mtimes` for matrix multiplication. +- Accept numeric row and column vectors consistently in statistics, moving and cumulative operations, differences, interpolation, and trapezoidal integration. Mixed INT/FLOAT samples are supported; sequence results remain flat lists. +- Validate complete matrix shapes and preserve row/column orientation through transpose, concatenation, and related operations. +- Preserve integer precision in integer statistics, compare extrema numerically, and report malformed numerical input as script errors in the revised paths. +- Return the fitted `intercept` from `regress`, alongside the existing predictor coefficients, p-values, and standard errors. Validate response shape and length before fitting. +- Add examples for local CSV calibration and residual diagnostics, matrix inversion, projectile motion, and explicit XOR backpropagation. + +### Compatibility and migration + +- `regress` fits an intercept automatically. Do not add a column of ones; calculate predictions as `intercept + X * parameters`. Existing result fields retain their meanings, with `intercept` added as a separate field. +- The Rhai names `transpose`, `horzcat`, `vertcat`, and `mtimes` remain available. Rust callers of matrix helpers should use the new `RhaiMatrix` wrappers or the corresponding `*_from_array` functions for array inputs. +- New constructors reject empty or ragged shapes. The flat-list identities for `sum`, `prod`, `diff`, and `unique` remain available on empty arrays. +- `dot` supports real vectors only; it does not implement matrix-axis overloads or complex conjugation. No implicit broadcasting is introduced. + +### Build and validation + +- Restore minimal-feature builds without the optional linear-algebra backend. +- Limit Polars to the supported 0.45 API and enable only its CSV feature. +- Test minimal, numerical, and full/documentation feature configurations in CI. diff --git a/Cargo.toml b/Cargo.toml index 70d6dd9..d2dee8d 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "rhai-sci" -version = "0.3.0" +version = "0.4.0" edition = "2021" authors = ["Chris McComb "] description = "Scientific computing in the Rhai scripting language" @@ -23,7 +23,7 @@ rand = ["randlib"] [dependencies] rhai = "1.8.0" nalgebralib = { version = ">=0.33.2,<1", optional = true, package = "nalgebra" } -polars = { version = ">=0.45.1,<1", optional = true } +polars = { version = "0.45.1", default-features = false, features = ["csv"], optional = true } url = { version = "2.0.0", optional = true } temp-file = { version = "0.2.0", optional = true } csv-sniffer = { version = "0.3.1", optional = true } @@ -35,7 +35,7 @@ linregress = { version = "0.5.0", optional = true } [build-dependencies] rhai = "1.8.0" nalgebralib = { version = ">=0.33.2,<1", optional = true, package = "nalgebra" } -polars = { version = ">=0.45.1,<1", optional = true } +polars = { version = "0.45.1", default-features = false, features = ["csv"], optional = true } url = { version = "2.0.0", optional = true } temp-file = { version = "0.2.0", optional = true } csv-sniffer = { version = "0.3.1", optional = true } diff --git a/README.md b/README.md index 5fc4d97..06f653c 100644 --- a/README.md +++ b/README.md @@ -2,52 +2,71 @@ [![Crates.io](https://img.shields.io/crates/v/rhai-sci.svg)](https://crates.io/crates/rhai-sci) [![docs.rs](https://img.shields.io/docsrs/rhai-sci/latest?logo=rust)](https://docs.rs/rhai-sci) -# About `rhai-sci` +# rhai-sci -This crate provides some basic scientific computing utilities for the [`Rhai`](https://rhai.rs/) scripting language, -inspired by languages like MATLAB, Octave, and R. For a complete API reference, -check [the docs](https://docs.rs/rhai-sci). +Scientific computing for the [Rhai](https://rhai.rs/) scripting language, inspired +by MATLAB, Octave, and R. Includes statistics, linear algebra, interpolation, +integration, and regression. -# Install +## Quickstart -To use the latest released version of `rhai-sci`, add this to your `Cargo.toml`: +Add the crate to your `Cargo.toml`: ```toml -rhai-sci = "0.2.3" +rhai-sci = "0.4.0" ``` -# Usage - -Using this crate is pretty simple! If you just want to evaluate a single line of [`Rhai`](https://rhai.rs/), then you -only need: +Evaluate a Rhai expression: ```rust use rhai::INT; use rhai_sci::eval; + let result = eval::("argmin([43, 42, -500])").unwrap(); +assert_eq!(result, 2); ``` -If you need to use `rhai-sci` as part of a persistent [`Rhai`](https://rhai.rs/) scripting engine, then do this instead: +For a persistent engine, register `SciPackage` as shown in the +[Rust host example](examples/regression_workflow.rs). -```rust -use rhai::{Engine, packages::Package, INT}; -use rhai_sci::SciPackage; +## Numerical workflows -// Create a new Rhai engine -let mut engine = Engine::new(); +Use `row`, `col`, and `mat` to construct vectors and matrices from ordinary Rhai +arrays. Statistics accept lists, rows, or columns; sequence results are flat lists. +Use `mtimes` for matrix multiplication and `dot` for a scalar vector inner product. -// Add the rhai-sci package to the new engine -engine.register_global_module(SciPackage::new().as_shared_module()); +`regress(X, y)` fits an intercept automatically and returns it separately from the +predictor coefficients. See the [workflow guide](docs/numerical-workflows.md) for +shape conventions, predictions, and input validation. -// Now run your code -let value = engine.eval::("argmin([43, 42, -500])").unwrap(); +Run the bundled CSV example to fit a model and summarize its residuals: + +```bash +cargo run --example regression_workflow ``` -# Features +More examples: [matrix inversion](examples/matrix_inversion.rhai), +[projectile motion](examples/projectile_motion.rhai), and +[XOR backpropagation](examples/neural_network_backprop.rhai). + +## Features + +| Feature | Default | Enables | +| --- | --- | --- | +| `io` | On | CSV loading with `read_matrix` | +| `nalgebra` | On | Matrix operations and regression | +| `rand` | On | Random values and matrices | +| `metadata` | Off | Function metadata and Rhai documentation tests | + +Disable default features and select only what you need for smaller builds; +CSV support brings in Polars. + +## Reference + +[API documentation](https://docs.rs/rhai-sci) · +[Changelog](CHANGELOG.md) · +[Development checks](docs/numerical-workflows.md#development) + +## License -| Feature | Default | Description | -|------------|----------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------| -| `metadata` | Disabled | Enables exporting function metadata and is ___necessary for running doc-tests on Rhai examples___. | -| `io` | Enabled | Enables the [`read_matrix`](#read_matrixfile_path-string---array) function but pulls in several additional dependencies (`polars`, `url`, `temp-file`, `csv-sniffer`, `minreq`). | -| `nalgebra` | Enabled | Enables several functions ([`regress`](#regressx-array-y-array---map), [`inv`](#invmatrix-array---array), [`mtimes`](#mtimesmatrix1-array-matrix2-array---array), [`horzcat`](#horzcatmatrix1-array-matrix2-array---array), [`vertcat`](#vertcatmatrix1-array-matrix2-array---array), [`repmat`](#repmatmatrix-array-nx-i64-ny-i64---array), [`svd`](#svdmatrix-array---map), [`hessenberg`](#hessenbergmatrix-array---map), and [`qr`](#qrmatrix-array---map)) but brings in the `nalgebra` and `linregress` crates. | -| `rand` | Enabled | Enables the [`rand`](#rand) function for generating random FLOAT values and random matrices, but brings in the `rand` crate. | +Licensed under [MIT](LICENSE-MIT.txt) or [Apache-2.0](LICENSE-APACHE.txt), at your option. diff --git a/build.rs b/build.rs index dde45cc..a44f009 100644 --- a/build.rs +++ b/build.rs @@ -225,3 +225,8 @@ mod functions { #[cfg(feature = "metadata")] pub use functions::*; + +#[cfg(feature = "metadata")] +pub mod matrix { + include!("src/matrix/mod.rs"); +} diff --git a/docs/numerical-workflows.md b/docs/numerical-workflows.md new file mode 100644 index 0000000..0c3a629 --- /dev/null +++ b/docs/numerical-workflows.md @@ -0,0 +1,95 @@ +# Numerical workflows + +[Back to the README](../README.md) + +## Matrix and vector conventions + +Matrices use ordinary Rhai arrays of rows. Constructors make orientation explicit +and validate numeric values while preserving their INT/FLOAT types: + +```typescript +let values = [1, 2, 3]; // plain Rhai list +let c = col(values); // N by 1: [[1], [2], [3]] +let r = row(values); // 1 by N: [[1, 2, 3]] +let A = mat([[1, 2], [3, 4]]); // rectangular numeric matrix +``` + +`row` and `col` also convert between vector orientations. Constructors reject empty, +ragged, or nonnumeric inputs. The returned arrays remain editable; matrix operations +check their inputs again. General matrix arithmetic does not implicitly broadcast +or turn a flat list into a row or column. + +Use `mtimes` for matrix multiplication and `dot` for a real scalar inner product: + +```typescript +let A = mat([[1, 2], [3, 4]]); +let x = col([5, 6]); +let prediction = mtimes(A, x); // [[17.0], [39.0]] +let energy = dot(x, x); // 61.0; also accepts lists or row vectors +let At = transpose(A); +let augmented = horzcat(A, x); +let extended = vertcat(A, row([7, 8])); +``` + +`dot` accepts equal-length, nonempty vectors in any combination of orientations and +returns FLOAT. It does not implement MATLAB's matrix/axis overloads or complex +conjugation. `mtimes` requires matching inner dimensions and returns a matrix, +including a 1 by 1 matrix for a row times a column. + +Statistics, moving/cumulative operations, differences, interpolation, and trapezoidal +integration accept numeric lists, rows, and columns. Mixed INT/FLOAT values are +supported. Scalar statistics return scalars; sequence operations return flat lists, +so use `col` or `row` when feeding those results back into matrix operations. +Inputs retain their original shape: + +```typescript +let samples = col([1, 2.0, 3]); +let average = mean(samples); // 2.0 +let smoothed = movmean(samples, 3); // [1.5, 2.0, 2.5] +let area = trapz(row([0, 1, 2]), samples); // 4.0 +``` + +Empty numeric samples produce a script error where a value is required. The +flat-list identities `sum([]) == 0`, `prod([]) == 1`, `diff([]) == []`, and +`unique([]) == []` remain available. + +## Regression and predictions + +`regress(X, y)` treats rows as observations and fits an intercept automatically. +Do not add a column of ones. The returned `parameters`, `pvalues`, and +`standard_errors` correspond to predictor columns in order. The fitted `intercept` +is returned separately: + +```typescript +let X = col([0, 1, 2]); +let fit = regress(X, col([1.1, 2.8, 5.1])); +let linear_part = mtimes(X, col(fit.parameters)); +let first_prediction = fit.intercept + linear_part[0][0]; +``` + +Earlier releases omitted the fitted intercept from the result. Existing result +fields are retained; use the new `intercept` field when calculating predictions. + +## CSV calibration example + +Run from the repository root: + +```bash +cargo run --example regression_workflow +``` + +This loads the bundled calibration CSV, constructs two predictor columns, fits a +linear model, predicts responses, and summarizes residuals. The sample has an +intercept of 1, slopes of 2 and 0.5, and an RMSE of approximately 0.1414. +It uses local data and requires the `io` and `nalgebra` features (both on by default). +Rust hosts that already have data can pass an `observations` array to the same +[script](../examples/regression_workflow.rhai) without enabling `io`. + +## Development + +```bash +cargo fmt --all -- --check +cargo test --no-default-features +cargo test --no-default-features --features rand,nalgebra +cargo test --all-features +``` diff --git a/examples/data/calibration.csv b/examples/data/calibration.csv new file mode 100644 index 0000000..0583805 --- /dev/null +++ b/examples/data/calibration.csv @@ -0,0 +1,7 @@ +temperature,load,response +0,0,1.1 +1,1,3.3 +2,0,5.1 +3,1,7.6 +4,0,8.8 +5,1,11.6 diff --git a/examples/download_and_regress.rhai b/examples/download_and_regress.rhai index fc6301f..97e552f 100644 --- a/examples/download_and_regress.rhai +++ b/examples/download_and_regress.rhai @@ -5,8 +5,8 @@ let x = read_matrix(url).transpose(); // Massage data let L = x.len; let y = x.drain(|v, i| i == (L-1)); -let x = ones(1, size(x)[1]) + x; +// regress fits the intercept automatically; no column of ones is needed. // Do regression and report let b = regress(x.transpose(), y.transpose()); -print(b); \ No newline at end of file +print(b); diff --git a/examples/matrix_inversion.rhai b/examples/matrix_inversion.rhai new file mode 100644 index 0000000..d70740e --- /dev/null +++ b/examples/matrix_inversion.rhai @@ -0,0 +1,4 @@ +let m = [[1, 2], [3, 4]]; +let inv_m = inv(m); +print(inv_m); +inv_m diff --git a/examples/matrix_inversion.rs b/examples/matrix_inversion.rs new file mode 100644 index 0000000..422d275 --- /dev/null +++ b/examples/matrix_inversion.rs @@ -0,0 +1,21 @@ +//! Demonstrates computing the inverse of a matrix using rhai-sci. + +fn main() { + #[cfg(feature = "nalgebra")] + { + use rhai::{packages::Package, Engine}; + use rhai_sci::SciPackage; + + // Create a new Rhai engine + let mut engine = Engine::new(); + + // Add the rhai-sci package to the engine + engine.register_global_module(SciPackage::new().as_shared_module()); + + // Run the script that inverts a matrix + let result = engine + .run_file("examples/matrix_inversion.rhai".into()) + .expect("script should run"); + println!("{:?}", result); + } +} diff --git a/examples/neural_network_backprop.rhai b/examples/neural_network_backprop.rhai new file mode 100644 index 0000000..161066b --- /dev/null +++ b/examples/neural_network_backprop.rhai @@ -0,0 +1,187 @@ +// Train a tiny 2-2-1 neural network on XOR with explicit backpropagation. + +fn sigmoid(x) { + 1.0 / (1.0 + exp(0.0 - x)) +} + +fn sigmoid_matrix(A) { + let out = []; + + for row in A { + let out_row = []; + for value in row { + out_row.push(sigmoid(value)); + } + out.push(out_row); + } + + out +} + +fn sigmoid_prime_from_activation(A) { + let out = []; + + for row in A { + let out_row = []; + for value in row { + out_row.push(value * (1.0 - value)); + } + out.push(out_row); + } + + out +} + +fn sub_matrix(A, B) { + let out = []; + + for i in 0..A.len() { + let a_row = A[i]; + let b_row = B[i]; + let row = []; + for j in 0..a_row.len() { + row.push(a_row[j] - b_row[j]); + } + out.push(row); + } + + out +} + +fn scale_matrix(A, scale) { + let out = []; + + for row in A { + let out_row = []; + for value in row { + out_row.push(value * scale); + } + out.push(out_row); + } + + out +} + +fn hadamard(A, B) { + let out = []; + + for i in 0..A.len() { + let a_row = A[i]; + let b_row = B[i]; + let row = []; + for j in 0..a_row.len() { + row.push(a_row[j] * b_row[j]); + } + out.push(row); + } + + out +} + +fn sum_squares(A) { + let total = 0.0; + + for row in A { + for value in row { + total += value * value; + } + } + + total +} + +fn forward(W1, W2, x) { + let x_with_bias = vertcat(x, col([1.0])); + let hidden = sigmoid_matrix(mtimes(W1, x_with_bias)); + let hidden_with_bias = vertcat(hidden, col([1.0])); + let output = sigmoid_matrix(mtimes(W2, hidden_with_bias)); + + #{ + "x_with_bias": x_with_bias, + "hidden": hidden, + "hidden_with_bias": hidden_with_bias, + "output": output + } +} + +fn total_loss(W1, W2, inputs, targets) { + let loss = 0.0; + + for i in 0..inputs.len() { + let x = inputs[i]; + let target_value = targets[i]; + let step = forward(W1, W2, x); + let prediction = step.output; + let target = col([target_value]); + let residual = sub_matrix(prediction, target); + loss += 0.5 * sum_squares(residual); + } + + loss +} + +fn predictions(W1, W2, inputs) { + let out = []; + + for x in inputs { + let prediction = forward(W1, W2, x).output; + let row = prediction[0]; + out.push(row[0]); + } + + out +} + +let inputs = [ + col([0.0, 0.0]), + col([0.0, 1.0]), + col([1.0, 0.0]), + col([1.0, 1.0]) +]; +let targets = [0.0, 1.0, 1.0, 0.0]; + +// Fixed starting weights keep the example deterministic while still learning. +let W1 = mat([[0.6888, 0.5159, -0.1589], [-0.4822, 0.0225, -0.1901]]); +let W2 = row([0.5676, -0.3934, -0.0468]); +let learning_rate = 1.5; +let epochs = 800; +let initial_loss = total_loss(W1, W2, inputs, targets); + +for epoch in 0..epochs { + for i in 0..inputs.len() { + let x = inputs[i]; + let target_value = targets[i]; + let step = forward(W1, W2, x); + let target = col([target_value]); + + let output_error = sub_matrix(step.output, target); + let output_slope = sigmoid_prime_from_activation(step.output); + let output_delta = hadamard(output_error, output_slope); + let output_weight_row = W2[0]; + let output_weights_without_bias = transpose(row([output_weight_row[0], output_weight_row[1]])); + + let hidden_error = mtimes(output_weights_without_bias, output_delta); + let hidden_slope = sigmoid_prime_from_activation(step.hidden); + let hidden_delta = hadamard(hidden_error, hidden_slope); + + let grad_W2 = mtimes(output_delta, transpose(step.hidden_with_bias)); + let grad_W1 = mtimes(hidden_delta, transpose(step.x_with_bias)); + + W2 = sub_matrix(W2, scale_matrix(grad_W2, learning_rate)); + W1 = sub_matrix(W1, scale_matrix(grad_W1, learning_rate)); + } +} + +let final_loss = total_loss(W1, W2, inputs, targets); +let final_predictions = predictions(W1, W2, inputs); + +let result = #{ + "initial_loss": initial_loss, + "final_loss": final_loss, + "predictions": final_predictions, + "W1": W1, + "W2": W2 +}; + +print(result); +result diff --git a/examples/neural_network_backprop.rs b/examples/neural_network_backprop.rs new file mode 100644 index 0000000..6045d45 --- /dev/null +++ b/examples/neural_network_backprop.rs @@ -0,0 +1,17 @@ +//! Trains a tiny neural network with backpropagation using rhai-sci matrices. + +fn main() { + #[cfg(feature = "nalgebra")] + { + use rhai::{packages::Package, Engine, Map}; + use rhai_sci::SciPackage; + + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + + let result: Map = engine + .eval_file("examples/neural_network_backprop.rhai".into()) + .expect("script should run"); + println!("{result:?}"); + } +} diff --git a/examples/projectile_motion.rhai b/examples/projectile_motion.rhai new file mode 100644 index 0000000..3b49c64 --- /dev/null +++ b/examples/projectile_motion.rhai @@ -0,0 +1,33 @@ +// Simulate simple projectile motion using rhai-sci functions +let g = 9.81; // gravitational acceleration (m/s^2) +let v0 = 25.0; // launch speed (m/s) +let angle = 45.0; // launch angle (degrees) + +// Generate a time vector from 0 to total flight time +let t_flight = 2.0 * v0 * sind(angle) / g; +let times = linspace(0.0, t_flight, 50); + +let vx = v0 * cosd(angle); +let vy0 = v0 * sind(angle); + +let x = []; +let y = []; + +for t in times { + x.push(vx * t); + y.push(vy0 * t - 0.5 * g * t * t); +} + +let max_y = max(y); +let idx = argmax(y); +let peak_time = times[idx]; +let range = max(x); + +let result = #{ + "max_height": max_y, + "time_of_flight": t_flight, + "peak_time": peak_time, + "range": range +}; +print(result); +result diff --git a/examples/projectile_motion.rs b/examples/projectile_motion.rs new file mode 100644 index 0000000..9e1f8bd --- /dev/null +++ b/examples/projectile_motion.rs @@ -0,0 +1,14 @@ +//! Simulates projectile motion using rhai-sci. + +fn main() { + use rhai::{packages::Package, Engine}; + use rhai_sci::SciPackage; + + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + + let result: rhai::Map = engine + .eval_file("examples/projectile_motion.rhai".into()) + .expect("script should run"); + println!("{result:?}"); +} diff --git a/examples/regression_workflow.rhai b/examples/regression_workflow.rhai new file mode 100644 index 0000000..46fc053 --- /dev/null +++ b/examples/regression_workflow.rhai @@ -0,0 +1,29 @@ +// The host loads a numeric CSV into `observations` (one observation per row). +// Columns are temperature, load, and measured response; see regression_workflow.rs. +let columns = transpose(mat(observations)); +let predictors = horzcat(col(columns[0]), col(columns[1])); +let response = col(columns[2]); + +// regress fits an intercept automatically; do not add a column of ones. +let fit = regress(predictors, response); +let linear_part = mtimes(predictors, col(fit.parameters)); +let predictions = []; +let residual_values = []; +for i in 0..observations.len() { + let prediction = fit.intercept + linear_part[i][0]; + predictions.push(prediction); + residual_values.push(response[i][0] - prediction); +} + +// Statistics consume the same column-vector representation as matrix operations. +let residuals = col(residual_values); +let result = #{ + "intercept": fit.intercept, + "parameters": fit.parameters, + "predictions": predictions, + "residual_mean": mean(residuals), + "residual_std": std(residuals), + "rmse": sqrt(dot(residuals, residuals) / residual_values.len().to_float()) +}; +print(result); +result diff --git a/examples/regression_workflow.rs b/examples/regression_workflow.rs new file mode 100644 index 0000000..0b5b22b --- /dev/null +++ b/examples/regression_workflow.rs @@ -0,0 +1,17 @@ +#[cfg(all(feature = "io", feature = "nalgebra"))] +fn main() -> Result<(), Box> { + use rhai::{packages::Package, Array, Engine, Scope}; + use rhai_sci::SciPackage; + + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + let observations = engine.eval::(r#"read_matrix("examples/data/calibration.csv")"#)?; + let mut scope = Scope::new(); + scope.push("observations", observations); + engine.run_with_scope(&mut scope, include_str!("regression_workflow.rhai")) +} + +#[cfg(not(all(feature = "io", feature = "nalgebra")))] +fn main() { + eprintln!("This CSV example requires the io and nalgebra features."); +} diff --git a/src/constants.rs b/src/constants.rs index de4c276..8641571 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -7,15 +7,15 @@ pub mod constant_definitions { // The ratio of a circle's circumference to its diameter. #[allow(non_upper_case_globals)] - pub const pi: FLOAT = 3.14159265358979323846264338327950288; + pub const pi: FLOAT = 3.141_592_653_589_793_238_462_643_383_279_502_88; //Speed of light in meters per second (m/s). #[allow(non_upper_case_globals)] - pub const c: FLOAT = 299792458.0; + pub const c: FLOAT = 299_792_458.0; // Euler's number. #[allow(non_upper_case_globals)] - pub const e: FLOAT = 2.71828182845904523536028747135266250; + pub const e: FLOAT = 2.718_281_828_459_045_235_360_287_471_352_662_50; // Acceleration due to gravity on Earth in meters per second per second (m/s^2). #[allow(non_upper_case_globals)] @@ -23,14 +23,14 @@ pub mod constant_definitions { // The Planck constant in Joules per Hertz (J/Hz) #[allow(non_upper_case_globals)] - pub const h: FLOAT = 6.62607015e-34; + pub const h: FLOAT = 6.626_070_15e-34; // The golden ratio #[allow(non_upper_case_globals)] - pub const phi: FLOAT = 1.61803398874989484820; + pub const phi: FLOAT = 1.618_033_988_749_894_848_20; // Newtonian gravitational constant - pub const G: FLOAT = 6.6743015e-11; + pub const G: FLOAT = 6.674_301_5e-11; /// Physical constants useful for science. /// ### `pi: FLOAT` diff --git a/src/cumulative.rs b/src/cumulative.rs index 602e792..5134011 100644 --- a/src/cumulative.rs +++ b/src/cumulative.rs @@ -72,25 +72,24 @@ pub mod cum_functions { /// ``` #[rhai_fn(name = "cumtrapz", return_raw)] pub fn cumtrapz(x: Array, y: Array) -> Result> { - if x.len() != y.len() { - Err(EvalAltResult::ErrorArithmetic( - "The arrays must have the same length".to_string(), - Position::NONE, - ) - .into()) - } else { - if_list_convert_to_vec_float_and_do(&mut y.clone(), |yf| { - if_list_convert_to_vec_float_and_do(&mut x.clone(), |xf| { - let mut trapsum = 0.0; - let mut cumtrapsum = vec![Dynamic::FLOAT_ZERO]; - for i in 1..x.len() { - trapsum += (yf[i] + yf[i - 1]) * (xf[i] - xf[i - 1]) / 2.0; - cumtrapsum.push(Dynamic::from_float(trapsum)); - } - Ok(cumtrapsum) - }) + if_list_convert_to_vec_float_and_do(&mut y.clone(), |yf| { + if_list_convert_to_vec_float_and_do(&mut x.clone(), |xf| { + if xf.len() != yf.len() { + return Err(EvalAltResult::ErrorArithmetic( + "The arrays must have the same length".into(), + Position::NONE, + ) + .into()); + } + let mut trapsum = 0.0; + let mut cumulative = vec![Dynamic::FLOAT_ZERO]; + for i in 1..xf.len() { + trapsum += (yf[i] + yf[i - 1]) * (xf[i] - xf[i - 1]) / 2.0; + cumulative.push(Dynamic::from_float(trapsum)); + } + Ok(cumulative) }) - } + }) } /// Returns the cumulative approximate integral of the curve defined by Y and x using the diff --git a/src/integration_and_differentiation.rs b/src/integration_and_differentiation.rs index ac48ac6..9823e3b 100644 --- a/src/integration_and_differentiation.rs +++ b/src/integration_and_differentiation.rs @@ -20,23 +20,23 @@ pub mod int_and_diff { /// ``` #[rhai_fn(name = "trapz", return_raw)] pub fn trapz(x: Array, y: Array) -> Result> { - if x.len() != y.len() { - Err(EvalAltResult::ErrorArithmetic( - "The arrays must have the same length".to_string(), - Position::NONE, - ) - .into()) - } else { - if_list_convert_to_vec_float_and_do(&mut y.clone(), |yf| { - if_list_convert_to_vec_float_and_do(&mut x.clone(), |xf| { - let mut trapsum = 0.0; - for i in 1..x.len() { - trapsum += (yf[i] + yf[i - 1]) * (xf[i] - xf[i - 1]) / 2.0; - } - Ok(Dynamic::from_float(trapsum)) - }) + if_list_convert_to_vec_float_and_do(&mut y.clone(), |yf| { + if_list_convert_to_vec_float_and_do(&mut x.clone(), |xf| { + if xf.len() != yf.len() { + return Err(EvalAltResult::ErrorArithmetic( + "The arrays must have the same length".into(), + Position::NONE, + ) + .into()); + } + let mut trapsum = 0.0; + + for i in 1..xf.len() { + trapsum += (yf[i] + yf[i - 1]) * (xf[i] - xf[i - 1]) / 2.0; + } + Ok(Dynamic::from_float(trapsum)) }) - } + }) } /// Returns the approximate integral of the curve defined by `y` using the trapezoidal method. @@ -70,6 +70,9 @@ pub mod int_and_diff { /// ``` #[rhai_fn(name = "diff", return_raw, pure)] pub fn diff(arr: &mut Array) -> Result> { + if arr.is_empty() { + return Ok(Array::new()); + } crate::if_list_do_int_or_do_float( arr, |arr| { diff --git a/src/lib.rs b/src/lib.rs index 48228b3..63098e6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -7,6 +7,8 @@ mod patterns; pub use patterns::*; +pub mod matrix; +pub use matrix::{RhaiMatrix, RhaiVector}; use rhai::{def_package, packages::Package, plugin::*, Engine, EvalAltResult}; mod matrices_and_arrays; pub use matrices_and_arrays::matrix_functions; diff --git a/src/matrices_and_arrays.rs b/src/matrices_and_arrays.rs index a03e61f..69968bc 100644 --- a/src/matrices_and_arrays.rs +++ b/src/matrices_and_arrays.rs @@ -2,20 +2,87 @@ use rhai::plugin::*; #[export_module] pub mod matrix_functions { + #[cfg(feature = "nalgebra")] + use crate::matrix::RhaiVector; + use crate::matrix::{numeric_vector_data, RhaiMatrix}; + use crate::validation_functions::{is_column_vector, is_row_vector}; use crate::{ array_to_vec_float, if_int_convert_to_float_and_do, if_int_do_else_if_array_do, if_list_do, if_matrix_convert_to_vec_array_and_do, }; #[cfg(feature = "nalgebra")] - use crate::{ - if_matrices_and_compatible_convert_to_vec_array_and_do, if_matrix_do, - omatrix_to_vec_dynamic, ovector_to_vec_dynamic, FOIL, - }; + use crate::{if_matrices_and_compatible_convert_to_vec_array_and_do, FOIL}; #[cfg(feature = "nalgebra")] use nalgebralib::DMatrix; use rhai::{Array, Dynamic, EvalAltResult, Map, Position, FLOAT, INT}; use std::collections::BTreeMap; + /// Construct a numeric row vector (1 by N) from a list or vector. + /// ```typescript + /// assert_eq(row([1, 2, 3]), [[1, 2, 3]]); + /// ``` + #[rhai_fn(name = "row", return_raw)] + pub fn row_from_array(values: Array) -> Result> { + Ok(RhaiMatrix::row_vector(numeric_vector_data(&values)?).to_array()) + } + + /// Construct a numeric column vector (N by 1) from a list or vector. + /// ```typescript + /// assert_eq(col([1, 2, 3]), [[1], [2], [3]]); + /// ``` + #[rhai_fn(name = "col", return_raw)] + pub fn col_from_array(values: Array) -> Result> { + Ok(RhaiMatrix::column_vector(numeric_vector_data(&values)?).to_array()) + } + + /// Validate a nonempty rectangular numeric matrix, preserving INT and FLOAT values. + /// The result is an ordinary Rhai array; later edits are validated by each operation. + /// ```typescript + /// assert_eq(mat([[1, 2.5], [3, 4]]), [[1, 2.5], [3, 4]]); + /// ``` + #[rhai_fn(name = "mat", return_raw)] + pub fn mat_from_array(values: Array) -> Result> { + if crate::matrix::matrix_dimensions(&values).is_none() { + return Err(EvalAltResult::ErrorArithmetic( + "mat expects nonempty row arrays of equal length".into(), + Position::NONE, + ) + .into()); + } + for row in &values { + numeric_vector_data(&row.clone().into_array().unwrap())?; + } + Ok(values) + } + + /// Compute a real scalar inner product of equal-length numeric vectors. + /// Accepts lists, rows, and columns in any combination and returns FLOAT. + /// Use `mtimes` for matrix multiplication; matrices with multiple rows and columns + /// are not accepted by this vector-only function. + /// ```typescript + /// assert_eq(dot(row([1, 2]), col([3, 4])), 11.0); + /// ``` + /// ```typescript + /// assert_eq(dot(col([1, 2]), col([3, 4])), 11.0); + /// ``` + #[rhai_fn(name = "dot", return_raw)] + pub fn dot(left: Array, right: Array) -> Result> { + let mut left = numeric_vector_data(&left)?; + let mut right = numeric_vector_data(&right)?; + if left.len() != right.len() { + return Err(EvalAltResult::ErrorArithmetic( + "dot expects vectors of the same length".into(), + Position::NONE, + ) + .into()); + } + Ok(array_to_vec_float(&mut left) + .iter() + .zip(array_to_vec_float(&mut right)) + .map(|(a, b)| a * b) + .sum()) + } + /// Calculates the inverse of a matrix. Fails if the matrix if not invertible, or if the /// elements of the matrix aren't FLOAT or INT. /// ```typescript @@ -39,26 +106,16 @@ pub mod matrix_functions { #[cfg(feature = "nalgebra")] #[rhai_fn(name = "inv", return_raw, pure)] pub fn invert_matrix(matrix: &mut Array) -> Result> { - if_matrix_convert_to_vec_array_and_do(matrix, |matrix_as_vec| { - let dm = DMatrix::from_fn(matrix_as_vec.len(), matrix_as_vec[0].len(), |i, j| { - if matrix_as_vec[0][0].is_float() { - matrix_as_vec[i][j].as_float().unwrap() - } else { - matrix_as_vec[i][j].as_int().unwrap() as FLOAT - } - }); - - // Try to invert - let dm = dm.try_inverse(); - - dm.map(omatrix_to_vec_dynamic).ok_or_else(|| { + let dm = RhaiMatrix::from_array(matrix.clone()).to_dmatrix()?; + dm.try_inverse() + .map(|m| RhaiMatrix::from_dmatrix(&m).to_array()) + .ok_or_else(|| { EvalAltResult::ErrorArithmetic( "Matrix cannot be inverted".to_string(), Position::NONE, ) .into() }) - }) } /// Calculate the eigenvalues and eigenvectors for a matrix. Specifically, the output is an @@ -78,73 +135,58 @@ pub mod matrix_functions { #[cfg(feature = "nalgebra")] #[rhai_fn(name = "eigs", return_raw, pure)] pub fn matrix_eigs_alt(matrix: &mut Array) -> Result> { - if_matrix_convert_to_vec_array_and_do(matrix, |matrix_as_vec| { - // Convert vec_array to omatrix - let dm = DMatrix::from_fn(matrix_as_vec.len(), matrix_as_vec[0].len(), |i, j| { - if matrix_as_vec[0][0].is_float() { - matrix_as_vec[i][j].as_float().unwrap() - } else { - matrix_as_vec[i][j].as_int().unwrap() as FLOAT - } - }); - - // Grab shape for later - let dms = dm.shape().1; - - // Get teh eigenvalues - let eigenvalues = dm.complex_eigenvalues(); - - // Iterate through eigenvalues to get eigenvectors - let mut imaginary_values = vec![Dynamic::from_float(1.0); 0]; - let mut real_values = vec![Dynamic::from_float(1.0); 0]; - let mut residuals = vec![Dynamic::from_float(1.0); 0]; - let mut eigenvectors = DMatrix::from_element(dms, 0, 0.0); - for (idx, ev) in eigenvalues.iter().enumerate() { - // Eigenvalue components - imaginary_values.push(Dynamic::from_float(ev.im)); - real_values.push(Dynamic::from_float(ev.re)); - - // Get eigenvector - let mut A = dm.clone() - DMatrix::from_diagonal_element(dms, dms, ev.re); - A = A.insert_column(0, 0.0); - A = A.insert_row(0, 0.0); - A[(0, idx + 1)] = 1.0; - let mut b = DMatrix::from_element(dms + 1, 1, 0.0); - b[(0, 0)] = 1.0; - let eigenvector = A - .svd(true, true) - .solve(&b, 1e-10) - .unwrap() - .remove_rows(0, 1) - .normalize(); - - // Verify solution - residuals.push(Dynamic::from_float( - (dm.clone() * eigenvector.clone() - ev.re * eigenvector.clone()).amax(), - )); - - eigenvectors.extend(eigenvector.column_iter()); - } + let dm = RhaiMatrix::from_array(matrix.clone()).to_dmatrix()?; + + let dms = dm.shape().1; + + let eigenvalues = dm.complex_eigenvalues(); + + let mut imaginary_values = vec![Dynamic::from_float(1.0); 0]; + let mut real_values = vec![Dynamic::from_float(1.0); 0]; + let mut residuals = vec![Dynamic::from_float(1.0); 0]; + let mut eigenvectors = DMatrix::from_element(dms, 0, 0.0); + for (idx, ev) in eigenvalues.iter().enumerate() { + imaginary_values.push(Dynamic::from_float(ev.im)); + real_values.push(Dynamic::from_float(ev.re)); + + let mut a = dm.clone() - DMatrix::from_diagonal_element(dms, dms, ev.re); + a = a.insert_column(0, 0.0); + a = a.insert_row(0, 0.0); + a[(0, idx + 1)] = 1.0; + let mut b = DMatrix::from_element(dms + 1, 1, 0.0); + b[(0, 0)] = 1.0; + let eigenvector = a + .svd(true, true) + .solve(&b, 1e-10) + .unwrap() + .remove_rows(0, 1) + .normalize(); + + residuals.push(Dynamic::from_float( + (dm.clone() * eigenvector.clone() - ev.re * eigenvector.clone()).amax(), + )); + + eigenvectors.extend(eigenvector.column_iter()); + } - let mut result = BTreeMap::new(); - let mut vid = smartstring::SmartString::new(); - vid.push_str("eigenvectors"); - result.insert( - vid, - Dynamic::from_array(omatrix_to_vec_dynamic(eigenvectors)), - ); - let mut did = smartstring::SmartString::new(); - did.push_str("real_eigenvalues"); - result.insert(did, Dynamic::from_array(real_values)); - let mut eid = smartstring::SmartString::new(); - eid.push_str("imaginary_eigenvalues"); - result.insert(eid, Dynamic::from_array(imaginary_values)); - let mut rid = smartstring::SmartString::new(); - rid.push_str("residuals"); - result.insert(rid, Dynamic::from_array(residuals)); - - Ok(result) - }) + let mut result = BTreeMap::new(); + let mut vid = smartstring::SmartString::new(); + vid.push_str("eigenvectors"); + result.insert( + vid, + Dynamic::from_array(RhaiMatrix::from_dmatrix(&eigenvectors).to_array()), + ); + let mut did = smartstring::SmartString::new(); + did.push_str("real_eigenvalues"); + result.insert(did, Dynamic::from_array(real_values)); + let mut eid = smartstring::SmartString::new(); + eid.push_str("imaginary_eigenvalues"); + result.insert(eid, Dynamic::from_array(imaginary_values)); + let mut rid = smartstring::SmartString::new(); + rid.push_str("residuals"); + result.insert(rid, Dynamic::from_array(residuals)); + + Ok(result) } /// Calculates the singular value decomposition of a matrix @@ -156,54 +198,50 @@ pub mod matrix_functions { #[cfg(feature = "nalgebra")] #[rhai_fn(name = "svd", return_raw, pure)] pub fn svd_decomp(matrix: &mut Array) -> Result> { - if_matrix_convert_to_vec_array_and_do(matrix, |matrix_as_vec| { - let dm = DMatrix::from_fn(matrix_as_vec.len(), matrix_as_vec[0].len(), |i, j| { - if matrix_as_vec[0][0].is::() { - matrix_as_vec[i][j].as_float().unwrap() - } else { - matrix_as_vec[i][j].as_int().unwrap() as FLOAT - } - }); - - // Try ot invert - let svd = nalgebralib::linalg::SVD::new(dm, true, true); - - let mut result = BTreeMap::new(); - let mut uid = smartstring::SmartString::new(); - uid.push_str("u"); - match svd.u { - Some(u) => result.insert(uid, Dynamic::from_array(omatrix_to_vec_dynamic(u))), - None => { - return Err(EvalAltResult::ErrorArithmetic( - format!("SVD decomposition cannot be computed for this matrix."), - Position::NONE, - ) - .into()); - } - }; - - let mut vid = smartstring::SmartString::new(); - vid.push_str("v"); - match svd.v_t { - Some(v) => result.insert(vid, Dynamic::from_array(omatrix_to_vec_dynamic(v))), - None => { - return Err(EvalAltResult::ErrorArithmetic( - format!("SVD decomposition cannot be computed for this matrix."), - Position::NONE, - ) - .into()); - } - }; + let dm = RhaiMatrix::from_array(matrix.clone()).to_dmatrix()?; + let svd = nalgebralib::linalg::SVD::new(dm, true, true); + + let mut result = BTreeMap::new(); + let mut u_key = smartstring::SmartString::new(); + u_key.push_str("u"); + match svd.u { + Some(u) => result.insert( + u_key, + Dynamic::from_array(RhaiMatrix::from_dmatrix(&u).to_array()), + ), + None => { + return Err(EvalAltResult::ErrorArithmetic( + "SVD decomposition cannot be computed for this matrix.".to_string(), + Position::NONE, + ) + .into()); + } + }; + + let mut v_key = smartstring::SmartString::new(); + v_key.push_str("v"); + match svd.v_t { + Some(v) => result.insert( + v_key, + Dynamic::from_array(RhaiMatrix::from_dmatrix(&v).to_array()), + ), + None => { + return Err(EvalAltResult::ErrorArithmetic( + "SVD decomposition cannot be computed for this matrix.".to_string(), + Position::NONE, + ) + .into()); + } + }; - let mut sid = smartstring::SmartString::new(); - sid.push_str("s"); - result.insert( - sid, - Dynamic::from_array(ovector_to_vec_dynamic(svd.singular_values)), - ); + let mut s_key = smartstring::SmartString::new(); + s_key.push_str("s"); + result.insert( + s_key, + Dynamic::from_array(RhaiVector::from_dvector(&svd.singular_values).to_array()), + ); - Ok(result) - }) + Ok(result) } /// Calculates the QR decomposition of a matrix @@ -215,29 +253,25 @@ pub mod matrix_functions { #[cfg(feature = "nalgebra")] #[rhai_fn(name = "qr", return_raw, pure)] pub fn qr_decomp(matrix: &mut Array) -> Result> { - if_matrix_convert_to_vec_array_and_do(matrix, |matrix_as_vec| { - let dm = DMatrix::from_fn(matrix_as_vec.len(), matrix_as_vec[0].len(), |i, j| { - if matrix_as_vec[0][0].is::() { - matrix_as_vec[i][j].as_float().unwrap() - } else { - matrix_as_vec[i][j].as_int().unwrap() as FLOAT - } - }); - - // Try ot invert - let qr = nalgebralib::linalg::QR::new(dm); - - let mut result = BTreeMap::new(); - let mut qid = smartstring::SmartString::new(); - qid.push_str("q"); - result.insert(qid, Dynamic::from_array(omatrix_to_vec_dynamic(qr.q()))); - - let mut rid = smartstring::SmartString::new(); - rid.push_str("r"); - result.insert(rid, Dynamic::from_array(omatrix_to_vec_dynamic(qr.r()))); - - Ok(result) - }) + let dm = RhaiMatrix::from_array(matrix.clone()).to_dmatrix()?; + let qr = nalgebralib::linalg::QR::new(dm); + + let mut result = BTreeMap::new(); + let mut qid = smartstring::SmartString::new(); + qid.push_str("q"); + result.insert( + qid, + Dynamic::from_array(RhaiMatrix::from_dmatrix(&qr.q()).to_array()), + ); + + let mut rid = smartstring::SmartString::new(); + rid.push_str("r"); + result.insert( + rid, + Dynamic::from_array(RhaiMatrix::from_dmatrix(&qr.r()).to_array()), + ); + + Ok(result) } /// Calculates the QR decomposition of a matrix @@ -249,58 +283,56 @@ pub mod matrix_functions { #[cfg(feature = "nalgebra")] #[rhai_fn(name = "hessenberg", return_raw, pure)] pub fn hessenberg(matrix: &mut Array) -> Result> { - if_matrix_convert_to_vec_array_and_do(matrix, |matrix_as_vec| { - let dm = DMatrix::from_fn(matrix_as_vec.len(), matrix_as_vec[0].len(), |i, j| { - if matrix_as_vec[0][0].is::() { - matrix_as_vec[i][j].as_float().unwrap() - } else { - matrix_as_vec[i][j].as_int().unwrap() as FLOAT - } - }); - - // Try ot invert - let h = nalgebralib::linalg::Hessenberg::new(dm); - - let mut result = BTreeMap::new(); - let mut hid = smartstring::SmartString::new(); - hid.push_str("h"); - result.insert(hid, Dynamic::from_array(omatrix_to_vec_dynamic(h.h()))); - - let mut qid = smartstring::SmartString::new(); - qid.push_str("q"); - result.insert(qid, Dynamic::from_array(omatrix_to_vec_dynamic(h.q()))); - - Ok(result) - }) + let dm = RhaiMatrix::from_array(matrix.clone()).to_dmatrix()?; + let h = nalgebralib::linalg::Hessenberg::new(dm); + + let mut result = BTreeMap::new(); + let mut hid = smartstring::SmartString::new(); + hid.push_str("h"); + result.insert( + hid, + Dynamic::from_array(RhaiMatrix::from_dmatrix(&h.h()).to_array()), + ); + + let mut qid = smartstring::SmartString::new(); + qid.push_str("q"); + result.insert( + qid, + Dynamic::from_array(RhaiMatrix::from_dmatrix(&h.q()).to_array()), + ); + + Ok(result) } /// Transposes a matrix. /// ```typescript /// let row = [[1, 2, 3, 4]]; /// let column = transpose(row); - /// assert_eq(column, [[1], - /// [2], - /// [3], - /// [4]]); + /// assert_eq(column, [[1.0], + /// [2.0], + /// [3.0], + /// [4.0]]); /// ``` /// ```typescript /// let matrix = transpose(eye(3)); /// assert_eq(matrix, eye(3)); /// ``` - #[rhai_fn(name = "transpose", pure, return_raw)] - pub fn transpose(matrix: &mut Array) -> Result> { - if_matrix_convert_to_vec_array_and_do(matrix, |matrix_as_vec| { - // Turn into Array - let mut out = vec![]; - for idx in 0..matrix_as_vec[0].len() { - let mut new_row = vec![]; - for jdx in 0..matrix_as_vec.len() { - new_row.push(matrix_as_vec[jdx][idx].clone()); - } - out.push(Dynamic::from_array(new_row)); - } - Ok(out) - }) + #[cfg(feature = "nalgebra")] + #[rhai_fn(name = "transpose", return_raw)] + pub fn transpose(matrix: RhaiMatrix) -> Result> { + let mut raw = matrix.clone().to_array(); + if !raw.is_empty() && matrix_size_by_reference(&mut raw).len() == 1 { + return RhaiMatrix::row_vector(raw).transpose(); + } + + matrix.transpose() + } + + /// Transpose an array by first converting it to a [`RhaiMatrix`]. + #[cfg(feature = "nalgebra")] + #[rhai_fn(name = "transpose", return_raw)] + pub fn transpose_from_array(matrix: Array) -> Result> { + transpose(RhaiMatrix::from_array(matrix)).map(RhaiMatrix::to_array) } /// Returns an array indicating the size of the matrix along each dimension, passed by reference. @@ -317,9 +349,9 @@ pub mod matrix_functions { let mut new_matrix = matrix.clone(); let mut shape = vec![Dynamic::from_int(new_matrix.len() as INT)]; - loop { - if new_matrix[0].is_array() { - new_matrix = new_matrix[0].clone().into_array().unwrap(); + while let Some(first) = new_matrix.first() { + if first.is_array() { + new_matrix = first.clone().into_array().unwrap(); shape.push(Dynamic::from_int(new_matrix.len() as INT)); } else { break; @@ -375,16 +407,15 @@ pub mod matrix_functions { .count() as INT } - #[cfg(all(feature = "io"))] + #[cfg(feature = "io")] pub mod read_write { use polars::prelude::{CsvReadOptions, DataType, SerReader}; use rhai::{Array, Dynamic, EvalAltResult, ImmutableString, FLOAT}; - /// Reads a numeric csv file from a url + /// Reads a numeric CSV file from the filesystem /// ```typescript - /// let url = "https://raw.githubusercontent.com/plotly/datasets/master/diabetes.csv"; - /// let x = read_matrix(url); - /// assert_eq(size(x), [768, 9]); + /// let x = read_matrix("tests/fixtures/sample_matrix.csv"); + /// assert_eq(x, [[1.0, 2.0], [3.0, 4.0]]); /// ``` #[rhai_fn(name = "read_matrix", return_raw)] pub fn read_matrix(file_path: ImmutableString) -> Result> { @@ -448,7 +479,7 @@ pub mod matrix_functions { // Convert into vec of vec let mut final_output = vec![]; - for series in x.columns() { + for series in x.get_columns() { let col: Vec = series .cast(&DataType::Float64) .map_err(|err| { @@ -715,6 +746,9 @@ pub mod matrix_functions { /// Returns an identity matrix. If argument is a single number, then the output is /// a square matrix. The argument can also be an array specifying the dimensions separately. + /// Passing `[n]` is equivalent to `eye(n)` (a square `n x n` matrix) while `[rows, cols]` + /// creates a rectangular matrix with the provided row and column counts. Any other shape is + /// rejected. /// ```typescript /// let matrix = eye(3); /// assert_eq(matrix, [[1.0, 0.0, 0.0], @@ -729,28 +763,36 @@ pub mod matrix_functions { /// ``` #[rhai_fn(name = "eye", return_raw)] pub fn eye_single_input(n: Dynamic) -> Result> { + fn parse_eye_dimension(value: &Dynamic) -> Result> { + value.as_int().map_err(|_| { + EvalAltResult::ErrorMismatchDataType( + "Size vector for eye must contain integers".to_string(), + String::new(), + Position::NONE, + ) + .into() + }) + } + if_int_do_else_if_array_do( n, |n| Ok(eye_double_input(n, n)), - |m| { - if m.len() == 1 { - Ok(eye_double_input(1, m[0].as_int().unwrap())[0] - .clone() - .into_array() - .unwrap()) - } else if m.len() == 2 { - Ok(eye_double_input( - m[0].as_int().unwrap(), - m[1].as_int().unwrap(), - )) - } else { - Err(EvalAltResult::ErrorMismatchDataType( - format!("Cannot create an identity matrix with more than 2 dimensions"), - format!(""), - Position::NONE, - ) - .into()) + |m| match m.len() { + 1 => { + let size = parse_eye_dimension(&m[0])?; + Ok(eye_double_input(size, size)) + } + 2 => { + let rows = parse_eye_dimension(&m[0])?; + let cols = parse_eye_dimension(&m[1])?; + Ok(eye_double_input(rows, cols)) } + _ => Err(EvalAltResult::ErrorMismatchDataType( + "Cannot create an identity matrix with more than 2 dimensions".to_string(), + String::new(), + Position::NONE, + ) + .into()), }, ) } @@ -921,157 +963,72 @@ pub mod matrix_functions { &mut matrix1.clone(), &mut matrix2.clone(), |matrix_as_vec1, matrix_as_vec2| { - let dm1 = - DMatrix::from_fn(matrix_as_vec1.len(), matrix_as_vec1[0].len(), |i, j| { - if matrix_as_vec1[0][0].is_float() { - matrix_as_vec1[i][j].as_float().unwrap() - } else { - matrix_as_vec1[i][j].as_int().unwrap() as FLOAT - } - }); - - let dm2 = - DMatrix::from_fn(matrix_as_vec2.len(), matrix_as_vec2[0].len(), |i, j| { - if matrix_as_vec2[0][0].is_float() { - matrix_as_vec2[i][j].as_float().unwrap() - } else { - matrix_as_vec2[i][j].as_int().unwrap() as FLOAT - } - }); - - // Try to multiply + let arr1: Array = matrix_as_vec1 + .into_iter() + .map(Dynamic::from_array) + .collect(); + let arr2: Array = matrix_as_vec2 + .into_iter() + .map(Dynamic::from_array) + .collect(); + let dm1 = RhaiMatrix::from_array(arr1).to_dmatrix()?; + let dm2 = RhaiMatrix::from_array(arr2).to_dmatrix()?; let mat = dm1 * dm2; - - // Turn into Array - let mut out = vec![]; - for idx in 0..mat.shape().0 { - let mut new_row = vec![]; - for jdx in 0..mat.shape().1 { - new_row.push(Dynamic::from_float(mat[(idx, jdx)])); - } - out.push(Dynamic::from_array(new_row)); - } - Ok(out) + Ok(RhaiMatrix::from_dmatrix(&mat).to_array()) }, ) } /// Concatenate two arrays horizontally. /// ```typescript - /// let arr1 = eye(3); - /// let arr2 = eye(3); - /// let combined = horzcat(arr1, arr2); - /// assert_eq(size(combined), [3, 6]); + /// let left = [[1, 2]]; + /// let right = [[3, 4]]; + /// let row = horzcat(left, right); + /// assert_eq(row, [[1.0, 2.0, 3.0, 4.0]]); /// ``` #[cfg(feature = "nalgebra")] #[rhai_fn(name = "horzcat", return_raw)] - pub fn horzcat(matrix1: Array, matrix2: Array) -> Result> { - if_matrices_and_compatible_convert_to_vec_array_and_do( - FOIL::First, - &mut matrix1.clone(), - &mut matrix2.clone(), - |matrix_as_vec1, matrix_as_vec2| { - let dm1 = - DMatrix::from_fn(matrix_as_vec1.len(), matrix_as_vec1[0].len(), |i, j| { - if matrix_as_vec1[0][0].is_float() { - matrix_as_vec1[i][j].as_float().unwrap() - } else { - matrix_as_vec1[i][j].as_int().unwrap() as FLOAT - } - }); - - let dm2 = - DMatrix::from_fn(matrix_as_vec2.len(), matrix_as_vec2[0].len(), |i, j| { - if matrix_as_vec2[0][0].is_float() { - matrix_as_vec2[i][j].as_float().unwrap() - } else { - matrix_as_vec2[i][j].as_int().unwrap() as FLOAT - } - }); - - // Try to multiple - let w0 = dm1.shape().1; - let w = dm1.shape().1 + dm2.shape().1; - let h = dm1.shape().0; - let mat = DMatrix::from_fn(h, w, |i, j| { - if j >= w0 { - dm2[(i, j - w0)] - } else { - dm1[(i, j)] - } - }); - - // Turn into Array - let mut out = vec![]; - for idx in 0..h { - let mut new_row = vec![]; - for jdx in 0..w { - new_row.push(Dynamic::from_float(mat[(idx, jdx)])); - } - out.push(Dynamic::from_array(new_row)); - } - Ok(out) - }, + pub fn horzcat( + matrix1: RhaiMatrix, + matrix2: RhaiMatrix, + ) -> Result> { + matrix1.concat_h(&matrix2) + } + + #[cfg(feature = "nalgebra")] + #[rhai_fn(name = "horzcat", return_raw)] + pub fn horzcat_from_array(matrix1: Array, matrix2: Array) -> Result> { + horzcat( + RhaiMatrix::from_array(matrix1), + RhaiMatrix::from_array(matrix2), ) + .map(RhaiMatrix::to_array) } /// Concatenates two array vertically. /// ```typescript - /// let arr1 = eye(3); - /// let arr2 = eye(3); - /// let combined = vertcat(arr1, arr2); - /// assert_eq(size(combined), [6, 3]); + /// let top = [[1], [2]]; + /// let bottom = [[3], [4]]; + /// let column = vertcat(top, bottom); + /// assert_eq(column, [[1.0], [2.0], [3.0], [4.0]]); /// ``` #[cfg(feature = "nalgebra")] #[rhai_fn(name = "vertcat", return_raw)] - pub fn vertcat(matrix1: Array, matrix2: Array) -> Result> { - if_matrices_and_compatible_convert_to_vec_array_and_do( - FOIL::Last, - &mut matrix1.clone(), - &mut matrix2.clone(), - |matrix_as_vec1, matrix_as_vec2| { - let dm1 = - DMatrix::from_fn(matrix_as_vec1.len(), matrix_as_vec1[0].len(), |i, j| { - if matrix_as_vec1[0][0].is_float() { - matrix_as_vec1[i][j].as_float().unwrap() - } else { - matrix_as_vec1[i][j].as_int().unwrap() as FLOAT - } - }); - - let dm2 = - DMatrix::from_fn(matrix_as_vec2.len(), matrix_as_vec2[0].len(), |i, j| { - if matrix_as_vec2[0][0].is_float() { - matrix_as_vec2[i][j].as_float().unwrap() - } else { - matrix_as_vec2[i][j].as_int().unwrap() as FLOAT - } - }); - - // Try to multiple - let h0 = dm1.shape().0; - let w = dm1.shape().1; - let h = dm1.shape().0 + dm2.shape().0; - let mat = DMatrix::from_fn(h, w, |i, j| { - if i >= h0 { - dm2[(i - h0, j)] - } else { - dm1[(i, j)] - } - }); - - // Turn into Array - let mut out = vec![]; - for idx in 0..h { - let mut new_row = vec![]; - for jdx in 0..w { - new_row.push(Dynamic::from_float(mat[(idx, jdx)])); - } - out.push(Dynamic::from_array(new_row)); - } - Ok(out) - }, + pub fn vertcat( + matrix1: RhaiMatrix, + matrix2: RhaiMatrix, + ) -> Result> { + matrix1.concat_v(&matrix2) + } + + #[cfg(feature = "nalgebra")] + #[rhai_fn(name = "vertcat", return_raw)] + pub fn vertcat_from_array(matrix1: Array, matrix2: Array) -> Result> { + vertcat( + RhaiMatrix::from_array(matrix1), + RhaiMatrix::from_array(matrix2), ) + .map(RhaiMatrix::to_array) } /// This function can be used in two distinct ways. @@ -1094,44 +1051,60 @@ pub mod matrix_functions { /// ``` #[rhai_fn(name = "diag", return_raw)] pub fn diag(matrix: Array) -> Result> { - if ndims_by_reference(&mut matrix.clone()) == 2 { - // Turn into Vec + let dims = ndims_by_reference(&mut matrix.clone()); + if dims == 2 { + let mut candidate_for_row = matrix.clone(); + let mut candidate_for_col = matrix.clone(); + if is_row_vector(&mut candidate_for_row) || is_column_vector(&mut candidate_for_col) { + let mut flattened_vector = matrix.clone(); + let vector = flatten(&mut flattened_vector); + return Ok(diagonal_matrix_from_vector(vector)); + } + let matrix_as_vec = matrix .into_iter() .map(|x| x.into_array().unwrap()) .collect::>(); - let mut out = vec![]; - for i in 0..matrix_as_vec.len() { - out.push(matrix_as_vec[i][i].clone()); + if matrix_as_vec.is_empty() { + return Ok(vec![]); } - Ok(out) - } else if ndims_by_reference(&mut matrix.clone()) == 1 { + let cols = matrix_as_vec[0].len(); + let diag_len = matrix_as_vec.len().min(cols); let mut out = vec![]; - for idx in 0..matrix.len() { - let mut new_row = vec![]; - for jdx in 0..matrix.len() { - if idx == jdx { - new_row.push(matrix[idx].clone()); - } else { - if matrix[idx].is_int() { - new_row.push(Dynamic::ZERO); - } else { - new_row.push(Dynamic::FLOAT_ZERO); - } - } - } - out.push(Dynamic::from_array(new_row)); + for i in 0..diag_len { + out.push(matrix_as_vec[i][i].clone()); } + Ok(out) + } else if dims == 1 { + Ok(diagonal_matrix_from_vector(matrix)) } else { - return Err(EvalAltResult::ErrorArithmetic( + Err(EvalAltResult::ErrorArithmetic( "Argument must be a 2-D matrix (to extract the diagonal) or a 1-D array (to create a matrix with that diagonal".to_string(), Position::NONE, ) - .into()); + .into()) + } + } + + fn diagonal_matrix_from_vector(vector: Array) -> Array { + let mut out = vec![]; + for idx in 0..vector.len() { + let mut new_row = vec![]; + for jdx in 0..vector.len() { + if idx == jdx { + new_row.push(vector[idx].clone()); + } else if vector[idx].is_int() { + new_row.push(Dynamic::ZERO); + } else { + new_row.push(Dynamic::FLOAT_ZERO); + } + } + out.push(Dynamic::from_array(new_row)); } + out } /// Repeats copies of a matrix @@ -1142,18 +1115,24 @@ pub mod matrix_functions { /// ``` #[cfg(feature = "nalgebra")] #[rhai_fn(name = "repmat", return_raw)] - pub fn repmat(matrix: &mut Array, nx: INT, ny: INT) -> Result> { - if_matrix_do(matrix, |matrix| { - let mut row_matrix = matrix.clone(); - for _ in 1..ny { - row_matrix = horzcat(row_matrix, matrix.clone())?; - } - let mut new_matrix = row_matrix.clone(); - for _ in 1..nx { - new_matrix = vertcat(new_matrix, row_matrix.clone())?; - } - Ok(new_matrix) - }) + pub fn repmat(matrix: RhaiMatrix, nx: INT, ny: INT) -> Result> { + let oriented = matrix + .as_column() + .or_else(|| matrix.as_row()) + .unwrap_or(matrix); + let dm = oriented.to_dmatrix()?; + let nx = if nx < 1 { 1 } else { nx as usize }; + let ny = if ny < 1 { 1 } else { ny as usize }; + let mat = DMatrix::from_fn(dm.nrows() * nx, dm.ncols() * ny, |i, j| { + dm[(i % dm.nrows(), j % dm.ncols())] + }); + Ok(RhaiMatrix::from_dmatrix(&mat)) + } + + #[cfg(feature = "nalgebra")] + #[rhai_fn(name = "repmat", return_raw)] + pub fn repmat_from_array(matrix: Array, nx: INT, ny: INT) -> Result> { + repmat(RhaiMatrix::from_array(matrix), nx, ny).map(RhaiMatrix::to_array) } /// Returns an object map containing 2-D grid coordinates based on the uni-axial coordinates @@ -1166,6 +1145,12 @@ pub mod matrix_functions { /// [1, 2]], /// "y": [[3, 3], /// [4, 4]]}); + /// + /// let x = [0, 1, 2]; + /// let y = [10]; + /// let g = meshgrid(x, y); + /// assert_eq(g, #{"x": [[0, 1, 2]], + /// "y": [[10, 10, 10]]}); /// ``` #[rhai_fn(name = "meshgrid", return_raw)] pub fn meshgrid(x: Array, y: Array) -> Result> { @@ -1173,8 +1158,14 @@ pub mod matrix_functions { if_list_do(&mut y.clone(), |y| { let nx = x.len(); let ny = y.len(); - let x_dyn: Array = vec![Dynamic::from_array(x.to_vec()); nx]; - let mut y_dyn: Array = vec![Dynamic::from_array(y.to_vec()); ny]; + let x_dyn: Array = (0..ny).map(|_| Dynamic::from_array(x.to_vec())).collect(); + let y_dyn: Array = y + .iter() + .map(|value| { + let y_row: Array = (0..nx).map(|_| value.clone()).collect(); + Dynamic::from_array(y_row) + }) + .collect(); let mut result = BTreeMap::new(); let mut xid = smartstring::SmartString::new(); @@ -1182,7 +1173,7 @@ pub mod matrix_functions { let mut yid = smartstring::SmartString::new(); yid.push_str("y"); result.insert(xid, Dynamic::from_array(x_dyn)); - result.insert(yid, Dynamic::from_array(transpose(&mut y_dyn).unwrap())); + result.insert(yid, Dynamic::from_array(y_dyn)); Ok(result) }) }) diff --git a/src/matrix/mod.rs b/src/matrix/mod.rs new file mode 100644 index 0000000..871680d --- /dev/null +++ b/src/matrix/mod.rs @@ -0,0 +1,392 @@ +#[cfg(feature = "nalgebra")] +use nalgebralib::{DMatrix, DVector}; +#[cfg(feature = "nalgebra")] +use rhai::FLOAT; +use rhai::{Array, Dynamic, EvalAltResult, Position}; + +/// Inspect every row rather than inferring a rectangular shape from the first row. +pub(crate) fn matrix_dimensions(values: &Array) -> Option<(usize, usize)> { + let first = values.first()?.clone().into_array().ok()?; + let cols = first.len(); + if cols == 0 { + return None; + } + for value in values { + let row = value.clone().into_array().ok()?; + if row.len() != cols || row.iter().any(Dynamic::is_array) { + return None; + } + } + Some((values.len(), cols)) +} + +/// Normalize a nonempty numeric list, row, or column without changing its values. +pub(crate) fn numeric_vector_data(values: &Array) -> Result> { + let numeric = |value: &Dynamic| value.is_int() || value.is_float(); + let data = if values.iter().all(|value| !value.is_array()) { + values.clone() + } else { + match matrix_dimensions(values) { + Some((1, _)) => values[0].clone().into_array().unwrap(), + Some((_, 1)) => values + .iter() + .map(|value| value.clone().into_array().unwrap().remove(0)) + .collect(), + _ => { + return Err(EvalAltResult::ErrorArithmetic( + "Expected a numeric list, row vector, or column vector".into(), + Position::NONE, + ) + .into()) + } + } + }; + if data.is_empty() { + return Err(EvalAltResult::ErrorArithmetic( + "Numeric vectors must contain at least one value".into(), + Position::NONE, + ) + .into()); + } + if !data.iter().all(numeric) { + return Err(EvalAltResult::ErrorArithmetic( + "Vector elements must be INT or FLOAT".into(), + Position::NONE, + ) + .into()); + } + Ok(data) +} + +/// Wrapper around [`rhai::Array`] representing a matrix. +/// +/// This type provides conversions between Rhai arrays and +/// `nalgebra::DMatrix`. +/// +/// # Examples +/// ``` +/// use rhai::{Array, Dynamic}; +/// use rhai_sci::matrix::RhaiMatrix; +/// let raw: Array = vec![ +/// Dynamic::from_array(vec![Dynamic::from_float(1.0), Dynamic::from_float(2.0)]), +/// Dynamic::from_array(vec![Dynamic::from_float(3.0), Dynamic::from_float(4.0)]), +/// ]; +/// let matrix = RhaiMatrix::from_array(raw.clone()); +/// assert_eq!(matrix.to_array().len(), raw.len()); +/// ``` +#[derive(Clone, Debug)] +pub struct RhaiMatrix(Array); + +impl RhaiMatrix { + /// Construct a [`RhaiMatrix`] from a [`rhai::Array`]. + #[must_use] + pub fn from_array(arr: Array) -> Self { + Self(arr) + } + + /// Construct a [`RhaiMatrix`] representing a row vector (`1×N`). + /// + /// # Examples + /// ``` + /// use rhai::{Array, Dynamic}; + /// use rhai_sci::matrix::RhaiMatrix; + /// let data: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + /// let row = RhaiMatrix::row_vector(data.clone()); + /// assert!(row.as_row().is_some()); + /// ``` + #[must_use] + pub fn row_vector(data: Array) -> Self { + Self(vec![Dynamic::from_array(data)]) + } + + /// Construct a [`RhaiMatrix`] representing a column vector (`N×1`). + /// + /// # Examples + /// ``` + /// use rhai::{Array, Dynamic}; + /// use rhai_sci::matrix::RhaiMatrix; + /// let data: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + /// let column = RhaiMatrix::column_vector(data.clone()); + /// assert!(column.as_column().is_some()); + /// ``` + #[must_use] + pub fn column_vector(data: Array) -> Self { + let rows = data + .into_iter() + .map(|v| Dynamic::from_array(vec![v])) + .collect(); + Self(rows) + } + + /// Return the matrix as a row vector (`1×N`), reshaping a column vector if necessary. + /// + /// Returns `None` if the matrix is not `1×N` or `N×1`. + /// + /// # Examples + /// ``` + /// use rhai::{Array, Dynamic}; + /// use rhai_sci::matrix::RhaiMatrix; + /// let data: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + /// let column = RhaiMatrix::column_vector(data.clone()); + /// let row = column.as_row().unwrap(); + /// assert!(row.as_row().is_some()); + /// ``` + #[must_use] + pub fn as_row(&self) -> Option { + match matrix_dimensions(&self.0) { + Some((1, _)) => Some(self.clone()), + Some((_, 1)) => Some(Self::row_vector(crate::matrix_functions::flatten( + &mut self.0.clone(), + ))), + _ => None, + } + } + + /// Return the matrix as a column vector (`N×1`), reshaping a row vector if necessary. + /// + /// Returns `None` if the matrix is not `1×N` or `N×1`. + /// + /// # Examples + /// ``` + /// use rhai::{Array, Dynamic}; + /// use rhai_sci::matrix::RhaiMatrix; + /// let data: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + /// let row = RhaiMatrix::row_vector(data.clone()); + /// let column = row.as_column().unwrap(); + /// assert!(column.as_column().is_some()); + /// ``` + #[must_use] + pub fn as_column(&self) -> Option { + match matrix_dimensions(&self.0) { + Some((_, 1)) => Some(self.clone()), + Some((1, _)) => Some(Self::column_vector(self.0[0].clone().into_array().unwrap())), + _ => None, + } + } + + /// Convert the matrix back into a [`rhai::Array`]. + #[must_use] + pub fn to_array(self) -> Array { + self.0 + } + + /// Convert the matrix into a `nalgebra::DMatrix`. + /// + /// # Errors + /// Returns an error if any element is non-numeric or rows have differing lengths. + /// + /// # Panics + /// Panics if an integer value cannot be represented as `FLOAT`. + #[cfg(feature = "nalgebra")] + #[allow(clippy::cast_precision_loss)] + pub fn to_dmatrix(&self) -> Result, Box> { + if self.0.is_empty() { + return Ok(DMatrix::from_element(0, 0, 0.0)); + } + let rows = self.0.len(); + let first_row = self.0[0].clone().into_array().map_err(|_| { + EvalAltResult::ErrorArithmetic( + "Matrix must contain row arrays".to_string(), + Position::NONE, + ) + })?; + let cols = first_row.len(); + let mut dm = DMatrix::zeros(rows, cols); + for (i, row_dyn) in self.0.iter().enumerate() { + let row = row_dyn.clone().into_array().map_err(|_| { + EvalAltResult::ErrorArithmetic( + "Matrix must contain row arrays".to_string(), + Position::NONE, + ) + })?; + if row.len() != cols { + return Err(EvalAltResult::ErrorArithmetic( + "Matrix rows must have equal length".to_string(), + Position::NONE, + ) + .into()); + } + for (j, val) in row.iter().enumerate() { + dm[(i, j)] = if val.is_float() { + val.as_float().unwrap() + } else if val.is_int() { + val.as_int().unwrap() as FLOAT + } else { + return Err(EvalAltResult::ErrorArithmetic( + "Matrix elements must be INT or FLOAT".to_string(), + Position::NONE, + ) + .into()); + }; + } + } + Ok(dm) + } + + /// Create a [`RhaiMatrix`] from a `nalgebra::DMatrix`. + #[cfg(feature = "nalgebra")] + #[must_use] + pub fn from_dmatrix(mat: &DMatrix) -> Self { + let mut rows = Vec::with_capacity(mat.nrows()); + for i in 0..mat.nrows() { + let mut row = Vec::with_capacity(mat.ncols()); + for j in 0..mat.ncols() { + row.push(Dynamic::from_float(mat[(i, j)])); + } + rows.push(Dynamic::from_array(row)); + } + Self(rows) + } + + /// Transpose the matrix. + /// + /// # Errors + /// Returns an error if the matrix contains non-numeric values or rows of + /// unequal length. + #[cfg(feature = "nalgebra")] + pub fn transpose(&self) -> Result> { + let dm = self.to_dmatrix()?; + Ok(Self::from_dmatrix(&dm.transpose())) + } + + /// Horizontally concatenate two matrices. + /// + /// # Examples + /// ``` + /// use rhai::{Array, Dynamic}; + /// use rhai_sci::{matrix::RhaiMatrix, validation_functions::is_row_vector}; + /// let left = RhaiMatrix::row_vector(vec![Dynamic::from_int(1), Dynamic::from_int(2)]); + /// let right = RhaiMatrix::row_vector(vec![Dynamic::from_int(3), Dynamic::from_int(4)]); + /// let combined = left.concat_h(&right).unwrap(); + /// let mut arr = combined.to_array(); + /// assert!(is_row_vector(&mut arr)); + /// ``` + /// # Errors + /// Returns an error if the matrices have differing row counts or contain + /// non-numeric values. + #[cfg(feature = "nalgebra")] + pub fn concat_h(&self, other: &Self) -> Result> { + let left = self.to_dmatrix()?; + let right = other.to_dmatrix()?; + if left.nrows() != right.nrows() { + return Err(EvalAltResult::ErrorArithmetic( + "Matrices must have the same number of rows".to_string(), + Position::NONE, + ) + .into()); + } + let cols = left.ncols() + right.ncols(); + let rows = left.nrows(); + let mat = DMatrix::from_fn(rows, cols, |i, j| { + if j < left.ncols() { + left[(i, j)] + } else { + right[(i, j - left.ncols())] + } + }); + Ok(Self::from_dmatrix(&mat)) + } + + /// Vertically concatenate two matrices. + /// + /// # Examples + /// ``` + /// use rhai::{Array, Dynamic}; + /// use rhai_sci::{matrix::RhaiMatrix, validation_functions::is_column_vector}; + /// let top = RhaiMatrix::column_vector(vec![Dynamic::from_int(1), Dynamic::from_int(2)]); + /// let bottom = RhaiMatrix::column_vector(vec![Dynamic::from_int(3), Dynamic::from_int(4)]); + /// let combined = top.concat_v(&bottom).unwrap(); + /// let mut arr = combined.to_array(); + /// assert!(is_column_vector(&mut arr)); + /// ``` + /// # Errors + /// Returns an error if the matrices have differing column counts or contain + /// non-numeric values. + #[cfg(feature = "nalgebra")] + pub fn concat_v(&self, other: &Self) -> Result> { + let top = self.to_dmatrix()?; + let bottom = other.to_dmatrix()?; + if top.ncols() != bottom.ncols() { + return Err(EvalAltResult::ErrorArithmetic( + "Matrices must have the same number of columns".to_string(), + Position::NONE, + ) + .into()); + } + let rows = top.nrows() + bottom.nrows(); + let cols = top.ncols(); + let mat = DMatrix::from_fn(rows, cols, |i, j| { + if i < top.nrows() { + top[(i, j)] + } else { + bottom[(i - top.nrows(), j)] + } + }); + Ok(Self::from_dmatrix(&mat)) + } +} + +/// Wrapper around [`rhai::Array`] representing a vector. +/// +/// # Examples +/// ``` +/// use rhai::{Array, Dynamic}; +/// use rhai_sci::matrix::RhaiVector; +/// let raw: Array = vec![Dynamic::from_float(1.0), Dynamic::from_float(2.0)]; +/// let vector = RhaiVector::from_array(raw.clone()); +/// assert_eq!(vector.to_array().len(), raw.len()); +/// ``` +#[derive(Clone, Debug)] +pub struct RhaiVector(Array); + +impl RhaiVector { + /// Construct a [`RhaiVector`] from a [`rhai::Array`]. + #[must_use] + pub fn from_array(arr: Array) -> Self { + Self(arr) + } + + /// Convert the vector back into a [`rhai::Array`]. + #[must_use] + pub fn to_array(self) -> Array { + self.0 + } + + /// Convert the vector into a `nalgebra::DVector`. + /// + /// # Errors + /// Returns an error if any element is non-numeric. + /// + /// # Panics + /// Panics if an integer value cannot be represented as `FLOAT`. + #[cfg(feature = "nalgebra")] + #[allow(clippy::cast_precision_loss)] + pub fn to_dvector(&self) -> Result, Box> { + let mut dv = DVector::zeros(self.0.len()); + for (i, val) in self.0.iter().enumerate() { + dv[i] = if val.is_float() { + val.as_float().unwrap() + } else if val.is_int() { + val.as_int().unwrap() as FLOAT + } else { + return Err(EvalAltResult::ErrorArithmetic( + "Vector elements must be INT or FLOAT".to_string(), + Position::NONE, + ) + .into()); + }; + } + Ok(dv) + } + + /// Create a [`RhaiVector`] from a `nalgebra::DVector`. + #[cfg(feature = "nalgebra")] + #[must_use] + pub fn from_dvector(vec: &DVector) -> Self { + let mut data = Vec::with_capacity(vec.len()); + for i in 0..vec.len() { + data.push(Dynamic::from_float(vec[i])); + } + Self(data) + } +} diff --git a/src/misc.rs b/src/misc.rs index 37c2b6a..2e0a8e8 100644 --- a/src/misc.rs +++ b/src/misc.rs @@ -28,6 +28,9 @@ pub mod misc_functions { /// ``` #[rhai_fn(name = "unique", return_raw, pure)] pub fn unique(arr: &mut Array) -> Result> { + if arr.is_empty() { + return Ok(Array::new()); + } if_list_do_int_or_do_float( arr, |arr| { @@ -71,25 +74,25 @@ pub mod misc_functions { .into()); }; - if x.len() < 2 { - return Err(EvalAltResult::ErrorArithmetic( - "The arrays must have at least 2 elements".to_string(), - Position::NONE, - ) - .into()); - } - if x.len() != y.len() { - return Err(EvalAltResult::ErrorArithmetic( - "The arrays must have the same length".to_string(), - Position::NONE, - ) - .into()); - } - let mut y = y; if_list_convert_to_vec_float_and_do(&mut y, |new_y| { if_list_convert_to_vec_float_and_do(x, |new_x| { + if new_x.len() < 2 { + return Err(EvalAltResult::ErrorArithmetic( + "The arrays must have at least 2 elements".to_string(), + Position::NONE, + ) + .into()); + } + if new_x.len() != new_y.len() { + return Err(EvalAltResult::ErrorArithmetic( + "The arrays must have the same length".to_string(), + Position::NONE, + ) + .into()); + } + if new_xq >= *new_x.last().unwrap() { return Ok(*new_y.last().unwrap()); } else if new_xq <= *new_x.first().unwrap() { diff --git a/src/patterns.rs b/src/patterns.rs index baef622..2d4344c 100644 --- a/src/patterns.rs +++ b/src/patterns.rs @@ -1,3 +1,5 @@ +#[cfg(feature = "nalgebra")] +use crate::matrix::{RhaiMatrix, RhaiVector}; use rhai::{Array, Dynamic, EvalAltResult, Position, FLOAT, INT}; /// Matrix compatibility conditions @@ -36,6 +38,8 @@ where FA: FnMut(&mut Array) -> Result>, FB: FnMut(&mut Array) -> Result>, { + let mut normalized = crate::matrix::numeric_vector_data(arr)?; + let arr = &mut normalized; let (int, float, total) = int_and_float_totals(arr); if int == total { f_int(arr) @@ -69,13 +73,8 @@ pub fn if_list_do(arr: &mut Array, mut f: F) -> Result Result>, { - crate::validation_functions::is_numeric_list(arr) - .then(|| f(arr)) - .unwrap_or(Err(EvalAltResult::ErrorArithmetic( - format!("The elements of the input array must either be INT or FLOAT."), - Position::NONE, - ) - .into())) + let mut normalized = crate::matrix::numeric_vector_data(arr)?; + f(&mut normalized) } pub fn if_list_convert_to_vec_float_and_do( @@ -189,11 +188,11 @@ pub fn if_matrix_convert_to_vec_array_and_do( where F: FnMut(Vec) -> Result>, { - let matrix_as_vec = matrix - .into_iter() - .map(|x| x.clone().into_array().unwrap()) - .collect::>(); if crate::validation_functions::is_matrix(matrix) { + let matrix_as_vec = matrix + .iter() + .map(|x| x.clone().into_array().unwrap()) + .collect::>(); f(matrix_as_vec) } else { Err(EvalAltResult::ErrorArithmetic( @@ -228,36 +227,32 @@ where pub fn array_to_vec_int(arr: &mut Array) -> Vec { arr.iter() - .map(|el| el.as_int().unwrap()) - .collect::>() + .map(|value| { + value.as_int().unwrap_or_else(|_| { + value.as_float().expect("Array elements must be numeric") as INT + }) + }) + .collect() } pub fn array_to_vec_float(arr: &mut Array) -> Vec { - arr.into_iter() - .map(|el| el.as_float().unwrap()) - .collect::>() + arr.iter() + .map(|value| { + value.as_float().unwrap_or_else(|_| { + value.as_int().expect("Array elements must be numeric") as FLOAT + }) + }) + .collect() } #[cfg(feature = "nalgebra")] pub fn omatrix_to_vec_dynamic( mat: nalgebralib::OMatrix, ) -> Vec { - let mut out = vec![]; - for idx in 0..mat.shape().0 { - let mut new_row = vec![]; - for jdx in 0..mat.shape().1 { - new_row.push(Dynamic::from_float(mat[(idx, jdx)])); - } - out.push(Dynamic::from_array(new_row)); - } - out + RhaiMatrix::from_dmatrix(&mat).to_array() } #[cfg(feature = "nalgebra")] pub fn ovector_to_vec_dynamic(mat: nalgebralib::OVector) -> Vec { - let mut out = vec![]; - for idx in 0..mat.shape().0 { - out.push(Dynamic::from_float(mat[idx])); - } - out + RhaiVector::from_dvector(&mat).to_array() } diff --git a/src/statistics.rs b/src/statistics.rs index 4228785..992c6d9 100644 --- a/src/statistics.rs +++ b/src/statistics.rs @@ -1,9 +1,31 @@ use rhai::plugin::*; +fn extremum_index( + values: &[T], + direction: std::cmp::Ordering, +) -> Result> { + let mut best = 0; + for (index, value) in values.iter().enumerate() { + let ordering = value.partial_cmp(&values[best]).ok_or_else(|| { + rhai::EvalAltResult::ErrorArithmetic( + "Cannot compare NaN values".into(), + rhai::Position::NONE, + ) + })?; + if ordering == direction { + best = index; + } + } + Ok(rhai::Dynamic::from_int(best as rhai::INT)) +} + #[export_module] pub mod stats { + use super::extremum_index; + #[cfg(feature = "nalgebra")] + use crate::matrix::RhaiMatrix; use crate::{ - array_to_vec_float, array_to_vec_int, if_list_convert_to_vec_float_and_do, if_list_do, + array_to_vec_float, array_to_vec_int, if_list_convert_to_vec_float_and_do, if_list_do_int_or_do_float, }; #[cfg(feature = "nalgebra")] @@ -208,6 +230,9 @@ pub mod stats { /// ``` #[rhai_fn(name = "sum", return_raw, pure)] pub fn sum(arr: &mut Array) -> Result> { + if arr.is_empty() { + return Ok(Dynamic::from_int(0)); + } if_list_do_int_or_do_float( arr, |arr| { @@ -230,13 +255,15 @@ pub mod stats { /// ``` #[rhai_fn(name = "mean", return_raw, pure)] pub fn mean(arr: &mut Array) -> Result> { - let l = arr.len() as FLOAT; if_list_do_int_or_do_float( arr, |arr: &mut Array| { - sum(arr).map(|s| Dynamic::from_float(s.as_int().unwrap() as FLOAT / l)) + sum(arr) + .map(|s| Dynamic::from_float(s.as_int().unwrap() as FLOAT / arr.len() as FLOAT)) + }, + |arr: &mut Array| { + sum(arr).map(|s| Dynamic::from_float(s.as_float().unwrap() / arr.len() as FLOAT)) }, - |arr: &mut Array| sum(arr).map(|s| Dynamic::from_float(s.as_float().unwrap() / l)), ) } @@ -249,15 +276,11 @@ pub mod stats { /// ``` #[rhai_fn(name = "argmax", return_raw, pure)] pub fn argmax(arr: &mut Array) -> Result> { - if_list_do(arr, |arr| { - array_max(arr).map(|m| { - Dynamic::from_int( - arr.iter() - .position(|r| format!("{r}") == format!("{m}")) - .unwrap() as INT, - ) - }) - }) + if_list_do_int_or_do_float( + arr, + |values| extremum_index(&array_to_vec_int(values), std::cmp::Ordering::Greater), + |values| extremum_index(&array_to_vec_float(values), std::cmp::Ordering::Greater), + ) } /// Return the index of the smallest array element. Fails if the input is not an array, or if @@ -269,15 +292,11 @@ pub mod stats { /// ``` #[rhai_fn(name = "argmin", return_raw, pure)] pub fn argmin(arr: &mut Array) -> Result> { - if_list_do(arr, |arr| { - array_min(arr).map(|m| { - Dynamic::from_int( - arr.iter() - .position(|r| format!("{r}") == format!("{m}")) - .unwrap() as INT, - ) - }) - }) + if_list_do_int_or_do_float( + arr, + |values| extremum_index(&array_to_vec_int(values), std::cmp::Ordering::Less), + |values| extremum_index(&array_to_vec_float(values), std::cmp::Ordering::Less), + ) } /// Compute the product of an array. Fails if the input is not an array, or if @@ -294,6 +313,9 @@ pub mod stats { /// ``` #[rhai_fn(name = "prod", return_raw, pure)] pub fn prod(arr: &mut Array) -> Result> { + if arr.is_empty() { + return Ok(Dynamic::from_int(1)); + } if_list_do_int_or_do_float( arr, |arr| { @@ -533,27 +555,38 @@ pub mod stats { ) } - /// Performs ordinary least squares regression and provides a statistical assessment. + /// Performs ordinary least squares regression with an automatically fitted intercept. + /// Rows of `x` are observations; columns are predictors. Do not add a column of ones. + /// `y` may be a list, row vector, or column vector with one value per observation. + /// `parameters`, `pvalues`, and `standard_errors` describe the predictor columns in order; + /// `intercept` is returned separately. Predict with intercept + X * parameters. /// ```typescript - /// let x = [[1.0, 0.0], - /// [1.0, 1.0], - /// [1.0, 2.0]]; + /// let x = col([0.0, 1.0, 2.0]); /// let y = [[0.1], /// [0.8], /// [2.1]]; /// let b = regress(x, y); - /// assert_eq(b, #{"parameters": [-2.220446049250313e-16, 1.0000000000000002], - /// "pvalues": [1.0, 0.10918255350924745], - /// "standard_errors": [0.11180339887498947, 0.17320508075688767]}); + /// assert_approx_eq(b.intercept, 0.0); + /// assert_approx_eq(b.parameters, [1.0]); /// ``` #[cfg(feature = "nalgebra")] #[rhai_fn(name = "regress", return_raw, pure)] pub fn regress(x: &mut Array, y: Array) -> Result> { use linregress::{FormulaRegressionBuilder, RegressionDataBuilder}; - let x_transposed = crate::matrix_functions::transpose(x)?; + crate::matrix_functions::mat_from_array(x.clone())?; + let mut response = crate::matrix::numeric_vector_data(&y)?; + if response.len() != x.len() { + return Err(EvalAltResult::ErrorArithmetic( + "regress expects one response per predictor row".into(), + Position::NONE, + ) + .into()); + } + let x_transposed = crate::matrix_functions::transpose(RhaiMatrix::from_array(x.clone()))?; + let x_arr = x_transposed.to_array(); let mut data: Vec<(String, Vec)> = vec![]; let mut vars = vec![]; - for (iter, column) in x_transposed.iter().enumerate() { + for (iter, column) in x_arr.iter().enumerate() { let var_name = format!("x_{iter}"); vars.push(var_name.clone()); data.push(( @@ -561,12 +594,11 @@ pub mod stats { array_to_vec_float(&mut column.clone().into_array().unwrap()), )); } - data.push(( - "y".to_string(), - array_to_vec_float(&mut crate::matrix_functions::flatten(&mut y.clone())), - )); + data.push(("y".to_string(), array_to_vec_float(&mut response))); - let regress_data = RegressionDataBuilder::new().build_from(data).unwrap(); + let regress_data = RegressionDataBuilder::new() + .build_from(data) + .map_err(|e| EvalAltResult::ErrorArithmetic(e.to_string(), Position::NONE))?; let model = FormulaRegressionBuilder::new() .data(®ress_data) @@ -594,6 +626,10 @@ pub mod stats { ); let mut result = BTreeMap::new(); + result.insert( + "intercept".into(), + Dynamic::from_float(model.parameters()[0]), + ); let mut params = smartstring::SmartString::new(); params.push_str("parameters"); result.insert(params, parameters); diff --git a/src/validate.rs b/src/validate.rs index 5c74c96..c5df933 100644 --- a/src/validate.rs +++ b/src/validate.rs @@ -2,9 +2,10 @@ use rhai::plugin::*; #[export_module] pub mod validation_functions { + use crate::matrix::{matrix_dimensions, numeric_vector_data}; use rhai::Array; - /// Tests whether the input in a simple list array + /// Tests whether the input is a simple list or a numeric vector. /// ```typescript /// let x = [1, 2, 3, 4]; /// assert_eq(is_list(x), true); @@ -15,11 +16,7 @@ pub mod validation_functions { /// ``` #[rhai_fn(name = "is_list", pure)] pub fn is_list(arr: &mut Array) -> bool { - if crate::matrix_functions::matrix_size_by_reference(arr).len() == 1 { - true - } else { - false - } + arr.iter().all(|value| !value.is_array()) || numeric_vector_data(arr).is_ok() } /// Determines if the entire array is numeric (ints or floats). @@ -79,7 +76,8 @@ pub mod validation_functions { }; } - /// Tests whether the input in a simple list array composed of either floating point or integer values. + /// Tests whether the input in a simple list array composed of either floating point or integer values + /// or a numeric row/column vector. /// ```typescript /// let x = [1.0, 2.0, 3.0, 4.0]; /// assert_eq(is_numeric_list(x), true) @@ -94,50 +92,35 @@ pub mod validation_functions { /// ``` #[rhai_fn(name = "is_numeric_list", pure)] pub fn is_numeric_list(arr: &mut Array) -> bool { - let (int, float, total) = crate::int_and_float_totals(arr); - if (int == total || float == total) && is_list(arr) { - true - } else { - false - } + numeric_vector_data(arr).is_ok() } - /// Tests whether the input is a row vector + /// Tests whether the input is a row vector. /// ```typescript - /// let x = ones([1, 5]); - /// assert_eq(is_row_vector(x), true) + /// let row = [[1, 2, 3]]; + /// assert_eq(is_row_vector(row), true); /// ``` /// ```typescript - /// let x = ones([5, 5]); - /// assert_eq(is_row_vector(x), false) + /// let column = [[1], [2], [3]]; + /// assert_eq(is_row_vector(column), false); /// ``` #[rhai_fn(name = "is_row_vector", pure)] pub fn is_row_vector(arr: &mut Array) -> bool { - let s = crate::matrix_functions::matrix_size_by_reference(arr); - if s.len() == 2 && s[0].as_int().unwrap() == 1 { - true - } else { - false - } + matches!(matrix_dimensions(arr), Some((1, _))) } - /// Tests whether the input is a column vector + /// Tests whether the input is a column vector. /// ```typescript - /// let x = ones([5, 1]); - /// assert_eq(is_column_vector(x), true) + /// let column = [[1], [2], [3]]; + /// assert_eq(is_column_vector(column), true); /// ``` /// ```typescript - /// let x = ones([5, 5]); - /// assert_eq(is_column_vector(x), false) + /// let row = [[1, 2, 3]]; + /// assert_eq(is_column_vector(row), false); /// ``` #[rhai_fn(name = "is_column_vector", pure)] pub fn is_column_vector(arr: &mut Array) -> bool { - let s = crate::matrix_functions::matrix_size_by_reference(arr); - if s.len() == 2 && s[1].as_int().unwrap() == 1 { - true - } else { - false - } + matches!(matrix_dimensions(arr), Some((_, 1))) } /// Tests whether the input is a matrix @@ -151,19 +134,6 @@ pub mod validation_functions { /// ``` #[rhai_fn(name = "is_matrix", pure)] pub fn is_matrix(arr: &mut Array) -> bool { - if crate::matrix_functions::matrix_size_by_reference(arr).len() != 2 { - false - } else { - if crate::stats::prod(&mut crate::matrix_functions::matrix_size_by_reference(arr)) - .unwrap() - .as_int() - .unwrap() - == crate::matrix_functions::numel_by_reference(arr) - { - true - } else { - false - } - } + matrix_dimensions(arr).is_some() } } diff --git a/tests/fixtures/sample_matrix.csv b/tests/fixtures/sample_matrix.csv new file mode 100644 index 0000000..ebb6763 --- /dev/null +++ b/tests/fixtures/sample_matrix.csv @@ -0,0 +1,2 @@ +1,2 +3,4 diff --git a/tests/list_like_inputs.rs b/tests/list_like_inputs.rs new file mode 100644 index 0000000..0f0c68e --- /dev/null +++ b/tests/list_like_inputs.rs @@ -0,0 +1,51 @@ +use rhai::{Array, Dynamic, FLOAT}; +use rhai_sci::moving_functions::movmean; +use rhai_sci::stats::argmax; + +fn row_vector(values: &[i64]) -> Array { + let row: Array = values + .iter() + .map(|value| Dynamic::from_int(*value)) + .collect(); + vec![Dynamic::from_array(row)] +} + +fn column_vector(values: &[i64]) -> Array { + values + .iter() + .map(|value| Dynamic::from_array(vec![Dynamic::from_int(*value)])) + .collect() +} + +fn array_to_floats(values: Array) -> Vec { + values + .into_iter() + .map(|value| value.as_float().unwrap()) + .collect() +} + +#[test] +fn movmean_accepts_row_and_column_vectors() { + let mut row = row_vector(&[1, 2, 3, 4]); + let mut column = column_vector(&[1, 2, 3, 4]); + + let expected = vec![1.5, 2.0, 3.0, 3.5]; + + let row_result = movmean(&mut row, 3).unwrap(); + assert_eq!(array_to_floats(row_result), expected); + + let column_result = movmean(&mut column, 3).unwrap(); + assert_eq!(array_to_floats(column_result), expected); +} + +#[test] +fn argmax_accepts_row_and_column_vectors() { + let mut row = row_vector(&[1, 9, 3]); + let mut column = column_vector(&[1, 9, 3]); + + let row_index = argmax(&mut row).unwrap(); + let column_index = argmax(&mut column).unwrap(); + + assert_eq!(row_index.as_int().unwrap(), 1); + assert_eq!(column_index.as_int().unwrap(), 1); +} diff --git a/tests/matrix_conventions.rs b/tests/matrix_conventions.rs new file mode 100644 index 0000000..a352383 --- /dev/null +++ b/tests/matrix_conventions.rs @@ -0,0 +1,166 @@ +use rhai::{packages::Package, Array, Dynamic, Engine, EvalAltResult}; +use rhai_sci::SciPackage; + +#[test] +fn constructors_preserve_values_and_convert_vector_orientation() { + for script in ["row([1, 2.5, 3])", "row(col([1, 2.5, 3]))"] { + assert_matrix_eq(eval_array(script).unwrap(), &[&[1.0, 2.5, 3.0]]); + } + for script in ["col([1, 2.5, 3])", "col(row([1, 2.5, 3]))"] { + assert_matrix_eq(eval_array(script).unwrap(), &[&[1.0], &[2.5], &[3.0]]); + } + let values = eval_array("mat([[1, 2.5]])").unwrap()[0] + .clone() + .into_array() + .unwrap(); + assert!(values[0].is_int()); + assert!(values[1].is_float()); +} + +#[test] +fn dot_is_a_scalar_inner_product_independent_of_orientation() { + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + for left in ["[1, 2]", "row([1, 2])", "col([1, 2])"] { + for right in ["[3, 4.0]", "row([3, 4.0])", "col([3, 4.0])"] { + let value = engine + .eval::(&format!("dot({left}, {right})")) + .unwrap(); + assert_eq!(value, 11.0); + } + } + assert_eq!( + engine.eval::("col([1, 2]).dot(row([3, 4]))").unwrap(), + 11.0 + ); + for script in [ + "dot([1], [1, 2])", + "dot([], [])", + "dot([[1, 2], [3, 4]], [1, 2])", + "dot([1, \"x\"], [1, 2])", + ] { + assert!(engine.eval::(script).is_err(), "{script}"); + } +} + +#[cfg(feature = "nalgebra")] +#[test] +fn matrix_products_and_concatenation_compose_with_constructors() { + assert_matrix_eq( + eval_array("mtimes(mat([[1, 2], [3, 4]]), col([5, 6]))").unwrap(), + &[&[17.0], &[39.0]], + ); + assert_matrix_eq( + eval_array("transpose(row([1, 2]))").unwrap(), + &[&[1.0], &[2.0]], + ); + assert_matrix_eq( + eval_array("transpose(col([1, 2]))").unwrap(), + &[&[1.0, 2.0]], + ); + assert_matrix_eq( + eval_array("transpose(transpose(col([1, 2])))").unwrap(), + &[&[1.0], &[2.0]], + ); + assert_matrix_eq( + eval_array("mtimes(row([1, 2]), col([3, 4]))").unwrap(), + &[&[11.0]], + ); + assert_matrix_eq( + eval_array("horzcat(mat([[1, 2], [3, 4]]), col([5, 6]))").unwrap(), + &[&[1.0, 2.0, 5.0], &[3.0, 4.0, 6.0]], + ); + assert_matrix_eq( + eval_array("vertcat(mat([[1, 2], [3, 4]]), row([5, 6]))").unwrap(), + &[&[1.0, 2.0], &[3.0, 4.0], &[5.0, 6.0]], + ); + assert_error_contains("horzcat(row([1, 2]), col([3, 4]))", "same number of rows"); + assert_error_contains( + "vertcat(row([1, 2]), col([3, 4]))", + "same number of columns", + ); + assert_error_contains("mtimes(col([1, 2]), col([3, 4]))", "not compatible"); + assert_error_contains( + "let A = mat([[1, 2], [3, 4]]); A[1] = [3]; mtimes(A, col([1, 2]))", + "matrix", + ); +} + +#[test] +fn constructors_reject_invalid_shapes_and_values() { + assert_error_contains("mat([[1, 2], [3]])", "equal length"); + assert_error_contains("mat([[1, \"x\"]])", "INT or FLOAT"); + assert_error_contains("mat([[]])", "nonempty"); + assert_error_contains("mat([])", "nonempty"); + assert_error_contains("row([])", "at least one value"); + assert_error_contains("col([])", "at least one value"); + for constructor in ["row", "col"] { + assert_error_contains(&format!("{constructor}([[1], [2, 3]])"), "vector"); + assert_error_contains(&format!("{constructor}([[1, 2], [3, 4]])"), "vector"); + } +} + +#[test] +fn shape_predicates_check_every_row() { + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + engine + .run( + r#" + assert_eq(is_matrix([[1], [2, 3]]), false); + assert_eq(is_column_vector([[1], [2, 3]]), false); + assert_eq(is_matrix([[1, 2], [3], [4, 5, 6]]), false); + assert_eq(is_row_vector([[[1, 2]]]), false); + assert_eq(is_numeric_list(row([1, 2.5])), true); + assert_eq(is_numeric_list(col([1, 2.5])), true); + "#, + ) + .unwrap(); +} + +fn eval_array(script: &str) -> Result> { + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + engine.eval::(script) +} + +fn assert_error_contains(script: &str, expected: &str) { + let err = eval_array(script).unwrap_err(); + match err.as_ref() { + EvalAltResult::ErrorArithmetic(message, _) => { + assert!( + message.contains(expected), + "expected error message `{message}` to contain `{expected}`" + ); + } + other => panic!("unexpected error: {other:?}"), + } +} + +fn assert_matrix_eq(actual: Array, expected: &[&[f64]]) { + let actual = numeric_matrix(actual); + let expected = expected + .iter() + .map(|row| row.to_vec()) + .collect::>>(); + assert_eq!(actual, expected); +} + +fn numeric_matrix(matrix: Array) -> Vec> { + matrix + .into_iter() + .map(|row| { + row.into_array() + .expect("matrix rows should be arrays") + .into_iter() + .map(|value| { + if value.is_float() { + value.as_float().expect("value should be FLOAT") + } else { + value.as_int().expect("value should be INT") as f64 + } + }) + .collect() + }) + .collect() +} diff --git a/tests/matrix_inverse_example.rs b/tests/matrix_inverse_example.rs new file mode 100644 index 0000000..593c8a6 --- /dev/null +++ b/tests/matrix_inverse_example.rs @@ -0,0 +1,24 @@ +#![cfg(feature = "nalgebra")] +use rhai::{packages::Package, Array, Engine}; +use rhai_sci::SciPackage; + +#[test] +fn matrix_inverse_example_produces_expected_result() { + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + + let result: Array = engine + .eval("inv([[1, 2], [3, 4]])") + .expect("script evaluation should succeed"); + + let first_row: Array = result[0].clone().cast::(); + let second_row: Array = result[1].clone().cast::(); + + let r0: Vec = first_row.into_iter().map(|v| v.cast::()).collect(); + let r1: Vec = second_row.into_iter().map(|v| v.cast::()).collect(); + + assert!((r0[0] + 2.0).abs() < f64::EPSILON); + assert!((r0[1] - 1.0).abs() < f64::EPSILON); + assert!((r1[0] - 1.5).abs() < f64::EPSILON); + assert!((r1[1] + 0.5).abs() < f64::EPSILON); +} diff --git a/tests/matrix_ops.rs b/tests/matrix_ops.rs new file mode 100644 index 0000000..f349a78 --- /dev/null +++ b/tests/matrix_ops.rs @@ -0,0 +1,247 @@ +#![cfg(feature = "nalgebra")] + +use rhai::{Array, Dynamic, EvalAltResult, FLOAT, INT}; +use rhai_sci::matrix::RhaiMatrix; +use rhai_sci::matrix_functions::{ + horzcat, matrix_size_by_reference, meshgrid, repmat, transpose, vertcat, +}; +use rhai_sci::validation_functions::{is_column_vector, is_row_vector}; + +#[test] +fn transpose_orients_row_vector() { + let data: Array = vec![ + Dynamic::from_int(1), + Dynamic::from_int(2), + Dynamic::from_int(3), + ]; + let row = RhaiMatrix::row_vector(data); + let mut result = transpose(row).unwrap().to_array(); + assert!(is_column_vector(&mut result)); +} + +#[test] +fn transpose_orients_column_vector() { + let data: Array = vec![ + Dynamic::from_int(1), + Dynamic::from_int(2), + Dynamic::from_int(3), + ]; + let column = RhaiMatrix::column_vector(data); + let mut result = transpose(column).unwrap().to_array(); + assert!(is_row_vector(&mut result)); + + let row = result[0].clone().into_array().unwrap(); + let values: Vec = row.into_iter().map(|d| d.as_float().unwrap()).collect(); + assert_eq!(values, vec![1.0, 2.0, 3.0]); +} + +#[test] +fn horzcat_concatenates_rows() { + let a: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + let b: Array = vec![Dynamic::from_int(3), Dynamic::from_int(4)]; + let m1 = RhaiMatrix::row_vector(a); + let m2 = RhaiMatrix::row_vector(b); + let mut result = horzcat(m1, m2).unwrap().to_array(); + assert!(is_row_vector(&mut result)); + let row = result[0].clone().into_array().unwrap(); + let values: Vec = row.into_iter().map(|d| d.as_float().unwrap()).collect(); + assert_eq!(values, vec![1.0, 2.0, 3.0, 4.0]); +} + +#[test] +fn horzcat_column_vectors_result_in_matrix_with_two_columns() { + let a = RhaiMatrix::column_vector(vec![Dynamic::from_int(1), Dynamic::from_int(2)]); + let b = RhaiMatrix::column_vector(vec![Dynamic::from_int(3), Dynamic::from_int(4)]); + let mut result = horzcat(a, b).unwrap().to_array(); + let shape = matrix_size_by_reference(&mut result); + let dims: Vec = shape.into_iter().map(|d| d.as_int().unwrap()).collect(); + assert_eq!(dims, vec![2, 2]); +} + +#[test] +fn horzcat_mixed_shapes_error_out() { + let row = RhaiMatrix::row_vector(vec![Dynamic::from_int(1), Dynamic::from_int(2)]); + let column = RhaiMatrix::column_vector(vec![Dynamic::from_int(3), Dynamic::from_int(4)]); + let err = horzcat(row, column).unwrap_err(); + match err.as_ref() { + EvalAltResult::ErrorArithmetic(message, _) => { + assert!(message.contains("same number of rows")); + } + other => panic!("unexpected error: {:?}", other), + } +} + +#[test] +fn vertcat_concatenates_columns() { + let a: Array = vec![ + Dynamic::from_array(vec![Dynamic::from_int(1)]), + Dynamic::from_array(vec![Dynamic::from_int(2)]), + ]; + let b: Array = vec![ + Dynamic::from_array(vec![Dynamic::from_int(3)]), + Dynamic::from_array(vec![Dynamic::from_int(4)]), + ]; + let m1 = RhaiMatrix::from_array(a); + let m2 = RhaiMatrix::from_array(b); + let mut result = vertcat(m1, m2).unwrap().to_array(); + assert!(is_column_vector(&mut result)); +} + +#[test] +fn vertcat_row_vectors_result_in_matrix_with_two_rows() { + let m1 = RhaiMatrix::row_vector(vec![Dynamic::from_int(1), Dynamic::from_int(2)]); + let m2 = RhaiMatrix::row_vector(vec![Dynamic::from_int(3), Dynamic::from_int(4)]); + let mut result = vertcat(m1, m2).unwrap().to_array(); + let shape = matrix_size_by_reference(&mut result); + let dims: Vec = shape.into_iter().map(|d| d.as_int().unwrap()).collect(); + assert_eq!(dims, vec![2, 2]); +} + +#[test] +fn vertcat_mixed_shapes_error_out() { + let column = RhaiMatrix::column_vector(vec![Dynamic::from_int(1), Dynamic::from_int(2)]); + let row = RhaiMatrix::row_vector(vec![Dynamic::from_int(3), Dynamic::from_int(4)]); + let err = vertcat(column, row).unwrap_err(); + match err.as_ref() { + EvalAltResult::ErrorArithmetic(message, _) => { + assert!(message.contains("same number of columns")); + } + other => panic!("unexpected error: {:?}", other), + } +} + +#[test] +fn repmat_replicates_matrix() { + let data: Array = vec![ + Dynamic::from_array(vec![Dynamic::from_int(1), Dynamic::from_int(2)]), + Dynamic::from_array(vec![Dynamic::from_int(3), Dynamic::from_int(4)]), + ]; + let m = RhaiMatrix::from_array(data); + let mut result = repmat(m, 2, 2).unwrap().to_array(); + let shape = matrix_size_by_reference(&mut result); + let dims: Vec = shape.into_iter().map(|d| d.as_int().unwrap()).collect(); + assert_eq!(dims, vec![4, 4]); +} + +#[test] +fn meshgrid_matches_matlab_shape_for_mismatched_lengths() { + let x: Array = vec![ + Dynamic::from_int(1), + Dynamic::from_int(2), + Dynamic::from_int(3), + ]; + let y: Array = vec![Dynamic::from_int(4), Dynamic::from_int(5)]; + let grid = meshgrid(x.clone(), y.clone()).unwrap(); + + let x_grid = grid.get("x").unwrap().clone().into_array().unwrap(); + let y_grid = grid.get("y").unwrap().clone().into_array().unwrap(); + + let mut x_grid_for_size = x_grid.clone(); + let x_shape = matrix_size_by_reference(&mut x_grid_for_size); + let x_dims: Vec = x_shape.into_iter().map(|d| d.as_int().unwrap()).collect(); + assert_eq!(x_dims, vec![2, 3]); + + let mut y_grid_for_size = y_grid.clone(); + let y_shape = matrix_size_by_reference(&mut y_grid_for_size); + let y_dims: Vec = y_shape.into_iter().map(|d| d.as_int().unwrap()).collect(); + assert_eq!(y_dims, vec![2, 3]); + + for row in x_grid.into_iter() { + let row_values: Vec = row + .into_array() + .unwrap() + .into_iter() + .map(|d| d.as_int().unwrap()) + .collect(); + assert_eq!(row_values, vec![1, 2, 3]); + } + + let y_rows: Vec> = y_grid + .into_iter() + .map(|row| { + row.into_array() + .unwrap() + .into_iter() + .map(|d| d.as_int().unwrap()) + .collect() + }) + .collect(); + assert_eq!(y_rows, vec![vec![4, 4, 4], vec![5, 5, 5]]); +} + +#[test] +fn meshgrid_accepts_row_vector_input() { + let row: Array = vec![Dynamic::from_array(vec![ + Dynamic::from_int(0), + Dynamic::from_int(1), + Dynamic::from_int(2), + ])]; + let y: Array = vec![Dynamic::from_int(3), Dynamic::from_int(4)]; + + let grid = meshgrid(row, y).unwrap(); + + let x_grid = grid.get("x").unwrap().clone().into_array().unwrap(); + let y_grid = grid.get("y").unwrap().clone().into_array().unwrap(); + + for row in x_grid.into_iter() { + let values: Vec = row + .into_array() + .unwrap() + .into_iter() + .map(|d| d.as_int().unwrap()) + .collect(); + assert_eq!(values, vec![0, 1, 2]); + } + + let y_rows: Vec> = y_grid + .into_iter() + .map(|row| { + row.into_array() + .unwrap() + .into_iter() + .map(|d| d.as_int().unwrap()) + .collect() + }) + .collect(); + assert_eq!(y_rows, vec![vec![3, 3, 3], vec![4, 4, 4]]); +} + +#[test] +fn meshgrid_accepts_column_vector_inputs() { + let column_x: Array = vec![ + Dynamic::from_array(vec![Dynamic::from_int(0)]), + Dynamic::from_array(vec![Dynamic::from_int(1)]), + Dynamic::from_array(vec![Dynamic::from_int(2)]), + ]; + let column_y: Array = vec![ + Dynamic::from_array(vec![Dynamic::from_int(3)]), + Dynamic::from_array(vec![Dynamic::from_int(4)]), + ]; + + let grid = meshgrid(column_x, column_y).unwrap(); + + let x_grid = grid.get("x").unwrap().clone().into_array().unwrap(); + let y_grid = grid.get("y").unwrap().clone().into_array().unwrap(); + + for row in x_grid.into_iter() { + let values: Vec = row + .into_array() + .unwrap() + .into_iter() + .map(|d| d.as_int().unwrap()) + .collect(); + assert_eq!(values, vec![0, 1, 2]); + } + + let y_rows: Vec> = y_grid + .into_iter() + .map(|row| { + row.into_array() + .unwrap() + .into_iter() + .map(|d| d.as_int().unwrap()) + .collect() + }) + .collect(); + assert_eq!(y_rows, vec![vec![3, 3, 3], vec![4, 4, 4]]); +} diff --git a/tests/matrix_vectors.rs b/tests/matrix_vectors.rs new file mode 100644 index 0000000..cac0f33 --- /dev/null +++ b/tests/matrix_vectors.rs @@ -0,0 +1,81 @@ +#![cfg(feature = "nalgebra")] + +use rhai::{Array, Dynamic}; +use rhai_sci::matrix::RhaiMatrix; +use rhai_sci::matrix_functions::{horzcat, matrix_size_by_reference, transpose, vertcat}; +use rhai_sci::validation_functions::{is_column_vector, is_row_vector}; + +#[test] +fn constructors_create_properly_oriented_vectors() { + // Row vector constructor produces 1xN matrix + let row_data: Array = vec![ + Dynamic::from_int(1), + Dynamic::from_int(2), + Dynamic::from_int(3), + ]; + let mut row = RhaiMatrix::row_vector(row_data).to_array(); + assert!(is_row_vector(&mut row)); + + // Column vector constructor produces Nx1 matrix + let column_data: Array = vec![ + Dynamic::from_int(4), + Dynamic::from_int(5), + Dynamic::from_int(6), + ]; + let mut column = RhaiMatrix::column_vector(column_data).to_array(); + assert!(is_column_vector(&mut column)); +} + +#[test] +fn as_column_converts_row_to_column() { + let data: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + let row = RhaiMatrix::row_vector(data); + let mut column = row.as_column().unwrap().to_array(); + assert!(is_column_vector(&mut column)); +} + +#[test] +fn as_row_converts_column_to_row() { + let data: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + let column = RhaiMatrix::column_vector(data); + let mut row = column.as_row().unwrap().to_array(); + assert!(is_row_vector(&mut row)); +} + +#[test] +fn transpose_flips_vector_orientation() { + let data: Array = vec![ + Dynamic::from_int(1), + Dynamic::from_int(2), + Dynamic::from_int(3), + ]; + let row = RhaiMatrix::row_vector(data); + let mut transposed = transpose(row).unwrap().to_array(); + assert!(is_column_vector(&mut transposed)); +} + +#[test] +fn horzcat_produces_row_vector() { + let left: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + let right: Array = vec![Dynamic::from_int(3), Dynamic::from_int(4)]; + let m1 = RhaiMatrix::row_vector(left); + let m2 = RhaiMatrix::row_vector(right); + let mut result = horzcat(m1, m2).unwrap().to_array(); + assert!(is_row_vector(&mut result)); + let dims = matrix_size_by_reference(&mut result); + assert_eq!(dims[0].as_int().unwrap(), 1); + assert_eq!(dims[1].as_int().unwrap(), 4); +} + +#[test] +fn vertcat_produces_column_vector() { + let top: Array = vec![Dynamic::from_int(1), Dynamic::from_int(2)]; + let bottom: Array = vec![Dynamic::from_int(3), Dynamic::from_int(4)]; + let m1 = RhaiMatrix::column_vector(top); + let m2 = RhaiMatrix::column_vector(bottom); + let mut result = vertcat(m1, m2).unwrap().to_array(); + assert!(is_column_vector(&mut result)); + let dims = matrix_size_by_reference(&mut result); + assert_eq!(dims[0].as_int().unwrap(), 4); + assert_eq!(dims[1].as_int().unwrap(), 1); +} diff --git a/tests/neural_network_backprop_example.rs b/tests/neural_network_backprop_example.rs new file mode 100644 index 0000000..82e2dfa --- /dev/null +++ b/tests/neural_network_backprop_example.rs @@ -0,0 +1,30 @@ +#![cfg(feature = "nalgebra")] + +use rhai::{packages::Package, Array, Engine, Map}; +use rhai_sci::SciPackage; + +#[test] +fn neural_network_backprop_example_learns_xor() { + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + + let result: Map = engine + .eval_file("examples/neural_network_backprop.rhai".into()) + .expect("script evaluation should succeed"); + + let initial_loss = result["initial_loss"].clone().cast::(); + let final_loss = result["final_loss"].clone().cast::(); + let predictions = result["predictions"].clone().cast::(); + let predictions = predictions + .into_iter() + .map(|value| value.cast::()) + .collect::>(); + + assert!(initial_loss > 0.45); + assert!(final_loss < 0.02); + assert!(final_loss < initial_loss); + assert!(predictions[0] < 0.15); + assert!(predictions[1] > 0.85); + assert!(predictions[2] > 0.85); + assert!(predictions[3] < 0.15); +} diff --git a/tests/orientation.rs b/tests/orientation.rs new file mode 100644 index 0000000..eeb3bd4 --- /dev/null +++ b/tests/orientation.rs @@ -0,0 +1,84 @@ +use rhai::{Array, Dynamic, FLOAT}; +use rhai_sci::matrix::RhaiMatrix; +use rhai_sci::matrix_functions; +use rhai_sci::validation_functions::{is_column_vector, is_row_vector}; + +#[test] +fn row_column_constructors_and_orientation() { + let data: Array = vec![ + Dynamic::from_int(1), + Dynamic::from_int(2), + Dynamic::from_int(3), + ]; + let row = RhaiMatrix::row_vector(data.clone()); + let column = RhaiMatrix::column_vector(data.clone()); + let mut row_to_col = row.as_column().unwrap().to_array(); + assert!(is_column_vector(&mut row_to_col)); + let mut col_to_row = column.as_row().unwrap().to_array(); + assert!(is_row_vector(&mut col_to_row)); +} + +#[test] +fn validate_orientation_helpers() { + let mut row: Array = vec![Dynamic::from_array(vec![ + Dynamic::from_int(1), + Dynamic::from_int(2), + ])]; + let column: Array = vec![ + Dynamic::from_array(vec![Dynamic::from_int(1)]), + Dynamic::from_array(vec![Dynamic::from_int(2)]), + ]; + assert!(is_row_vector(&mut row.clone())); + assert!(!is_row_vector(&mut column.clone())); + assert!(is_column_vector(&mut column.clone())); + assert!(!is_column_vector(&mut row)); +} + +#[test] +fn eye_vector_size_matches_scalar_size() { + for size in [1, 2, 5] { + let scalar = matrix_functions::eye_single_input(Dynamic::from_int(size)) + .expect("scalar eye should succeed"); + let vector = + matrix_functions::eye_single_input(Dynamic::from_array(vec![Dynamic::from_int(size)])) + .expect("vector eye should succeed"); + assert_eq!( + normalize_matrix(&scalar), + normalize_matrix(&vector), + "size {size} should match", + ); + } +} + +#[test] +fn eye_vector_two_dimensions_remains_rectangular() { + let rectangular = matrix_functions::eye_single_input(Dynamic::from_array(vec![ + Dynamic::from_int(2), + Dynamic::from_int(3), + ])) + .expect("rectangular eye should succeed"); + assert_eq!( + normalize_matrix(&rectangular), + normalize_matrix(&matrix_functions::eye_double_input(2, 3)), + ); +} + +fn normalize_matrix(matrix: &Array) -> Vec> { + matrix + .iter() + .map(|row| { + row.clone() + .into_array() + .expect("matrix rows should be arrays") + .into_iter() + .map(|value| { + if value.is_float() { + value.as_float().expect("value is float") + } else { + value.as_int().expect("value is int") as FLOAT + } + }) + .collect() + }) + .collect() +} diff --git a/tests/projectile_motion_example.rs b/tests/projectile_motion_example.rs new file mode 100644 index 0000000..64f49b9 --- /dev/null +++ b/tests/projectile_motion_example.rs @@ -0,0 +1,28 @@ +use rhai::{packages::Package, Engine, Map}; +use rhai_sci::SciPackage; + +#[test] +fn projectile_motion_example_produces_expected_result() { + // Arrange: set up engine with rhai-sci package + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + + // Act: evaluate projectile motion script + let result: Map = engine + .eval_file("examples/projectile_motion.rhai".into()) + .expect("script evaluation should succeed"); + + // Assert: compare against analytical solution + let max_height = result["max_height"].clone().cast::(); + let time_of_flight = result["time_of_flight"].clone().cast::(); + let range = result["range"].clone().cast::(); + + let expected_max_height = + (25.0_f64.powi(2) * (45_f64.to_radians().sin().powi(2))) / (2.0 * 9.81); + let expected_time_of_flight = 2.0 * 25.0 * 45_f64.to_radians().sin() / 9.81; + let expected_range = (25.0_f64.powi(2) * (2.0 * 45_f64.to_radians()).sin()) / 9.81; + + assert!((max_height - expected_max_height).abs() < 1e-2); + assert!((time_of_flight - expected_time_of_flight).abs() < 1e-6); + assert!((range - expected_range).abs() < 1e-6); +} diff --git a/tests/regression_workflow.rs b/tests/regression_workflow.rs new file mode 100644 index 0000000..7de0d11 --- /dev/null +++ b/tests/regression_workflow.rs @@ -0,0 +1,79 @@ +#![cfg(feature = "nalgebra")] + +use rhai::{packages::Package, Array, Dynamic, Engine, Map, Scope}; +use rhai_sci::SciPackage; + +fn engine() -> Engine { + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + engine +} + +fn check_calibration(engine: &Engine, observations: Array) { + let mut scope = Scope::new(); + scope.push("observations", observations); + let result: Map = engine + .eval_with_scope( + &mut scope, + include_str!("../examples/regression_workflow.rhai"), + ) + .unwrap(); + let parameters = result["parameters"].clone().cast::(); + // The chosen residuals are orthogonal to the intercept and both predictors. + assert!((result["intercept"].as_float().unwrap() - 1.0).abs() < 1e-10); + assert!((parameters[0].as_float().unwrap() - 2.0).abs() < 1e-10); + assert!((parameters[1].as_float().unwrap() - 0.5).abs() < 1e-10); + assert!(result["residual_mean"].as_float().unwrap().abs() < 1e-10); + assert!((result["rmse"].as_float().unwrap() - 0.02_f64.sqrt()).abs() < 1e-10); + assert!((result["residual_std"].as_float().unwrap() - 0.024_f64.sqrt()).abs() < 1e-10); + let predictions = result["predictions"].clone().cast::(); + for (prediction, expected) in predictions.iter().zip([1.0, 3.5, 5.0, 7.5, 9.0, 11.5]) { + assert!((prediction.as_float().unwrap() - expected).abs() < 1e-10); + } +} + +#[test] +fn embedded_regression_predicts_and_summarizes_residuals() { + let engine = engine(); + let observations: Array = engine + .eval("[[0, 0, 1.1], [1, 1, 3.3], [2, 0, 5.1], [3, 1, 7.6], [4, 0, 8.8], [5, 1, 11.6]]") + .unwrap(); + check_calibration(&engine, observations); +} + +#[cfg(feature = "io")] +#[test] +fn csv_regression_workflow_loads_the_bundled_data() { + let engine = engine(); + let observations = engine + .eval::(r#"read_matrix("examples/data/calibration.csv")"#) + .unwrap(); + check_calibration(&engine, observations); +} + +#[test] +fn regression_accepts_vector_responses_and_reports_invalid_data() { + let engine = engine(); + for response in [ + "[1.1, 2.8, 5.1]", + "row([1.1, 2.8, 5.1])", + "col([1.1, 2.8, 5.1])", + ] { + let fit: Map = engine + .eval(&format!("regress(col([0, 1, 2]), {response})")) + .unwrap(); + assert!((fit["intercept"].as_float().unwrap() - 1.0).abs() < 1e-10); + let parameters = fit["parameters"].clone().cast::(); + assert!((parameters[0].as_float().unwrap() - 2.0).abs() < 1e-10); + } + for script in [ + "regress([], [])", + "regress([[1], [2, 3]], [1, 2])", + "regress(col([1, 2, 3]), [1, 2])", + "regress(col([1, 2, 3]), [[1, 2], [3, 4]])", + "regress(col([1, 2, 3]), [1, \"x\", 3])", + "regress(col([1, 2, 3]), [inf, inf, inf])", + ] { + assert!(engine.eval::(script).is_err(), "{script}"); + } +} diff --git a/tests/rhai-sci-tests.rs b/tests/rhai-sci-tests.rs index f833da3..f5f4a0c 100644 --- a/tests/rhai-sci-tests.rs +++ b/tests/rhai-sci-tests.rs @@ -1,2 +1,2 @@ #[cfg(feature = "metadata")] -include!(concat!(env!("OUT_DIR"), "/rhai-sci-tests.rs")); \ No newline at end of file +include!(concat!(env!("OUT_DIR"), "/rhai-sci-tests.rs")); diff --git a/tests/vector_workflows.rs b/tests/vector_workflows.rs new file mode 100644 index 0000000..4bf38e6 --- /dev/null +++ b/tests/vector_workflows.rs @@ -0,0 +1,111 @@ +use rhai::{packages::Package, Engine}; +use rhai_sci::SciPackage; + +fn engine() -> Engine { + let mut engine = Engine::new(); + engine.register_global_module(SciPackage::new().as_shared_module()); + engine +} + +#[test] +fn statistics_and_signal_operations_accept_every_vector_orientation() { + let engine = engine(); + for values in ["[1, 2, 3]", "[1, 2.0, 3]"] { + for vector in [ + values.to_string(), + format!("row({values})"), + format!("col({values})"), + ] { + let script = format!( + r#" + let samples = {vector}; + let before = samples; + assert_eq(mean(samples), 2.0); + assert_eq(std(samples), 1.0); + assert_eq(median(samples), 2.0); + assert_approx_eq(min(samples).to_float(), 1.0); + assert_approx_eq(max(samples).to_float(), 3.0); + assert_eq(argmax(samples), 2); + assert_eq(argmin(samples), 0); + assert_approx_eq(movmean(samples, 3), [1.5, 2.0, 2.5]); + assert_approx_eq(cumsum(samples), [1.0, 3.0, 6.0]); + assert_approx_eq(diff(samples), [1.0, 1.0]); + assert_eq(samples, before); + "# + ); + engine + .run(&script) + .unwrap_or_else(|error| panic!("{vector}: {error}")); + } + } +} + +#[test] +fn paired_sample_operations_compare_vector_lengths_after_normalizing() { + engine() + .run( + r#" + let x = row([0, 1, 2]); + let y = col([1, 3.0, 5]); + assert_eq(trapz(x, y), 6.0); + assert_approx_eq(cumtrapz(x, y), [0.0, 2.0, 6.0]); + assert_eq(interp1(x, y, 0.5), 2.0); + "#, + ) + .unwrap(); +} + +#[test] +fn malformed_vectors_return_script_errors() { + let engine = engine(); + for data in [ + "[[1], [2, 3]]", + "[[1, 2], [3, 4]]", + "[1, [2]]", + "[[1, \"x\"]]", + "[]", + "[[]]", + ] { + for function in ["mean", "median", "std", "max", "argmax"] { + let script = format!("{function}({data})"); + assert!(engine.eval::(&script).is_err(), "{script}"); + } + } + for script in [ + "trapz(row([1, 2]), col([1, 2, 3]))", + "cumtrapz(row([1, 2]), row([1, 2, 3]))", + "interp1(row([1, 2]), row([1, 2, 3]), 1.5)", + ] { + assert!(engine.eval::(script).is_err(), "{script}"); + } +} + +#[test] +fn integer_statistics_preserve_values_above_float_precision() { + engine() + .run( + r#" + let n = 9007199254740993; + assert_eq(max(col([n, n - 1])), n); + assert_eq(argmax(col([n - 1, n])), 1); + assert_eq(sum(row([n, 0])), n); + "#, + ) + .unwrap(); +} + +#[test] +fn empty_list_identities_remain_available() { + engine() + .run( + r#" + assert_eq(sum([]), 0); + assert_eq(prod([]), 1); + assert_eq(diff([]), []); + assert_eq(unique([]), []); + assert_eq(size([]), [0]); + assert_eq(size([[]]), [1, 0]); + "#, + ) + .unwrap(); +}