diff --git a/.github/workflows/Github.yml b/.github/workflows/Github.yml index 638636c..4719a0f 100644 --- a/.github/workflows/Github.yml +++ b/.github/workflows/Github.yml @@ -30,10 +30,14 @@ jobs: - uses: actions/checkout@v4 - name: Install OpenBLAS / LAPACK run: sudo apt-get update && sudo apt-get install -y libopenblas-dev liblapack-dev gfortran - - name: Build + - name: Build (OpenBLAS from source) run: cargo build --verbose --features O3-openblas - - name: Run tests + - name: Run tests (OpenBLAS from source) run: cargo test --verbose --features O3-openblas + - name: Build (system OpenBLAS via pkg-config) + run: cargo build --verbose --features O3-openblas-system + - name: Run tests (system OpenBLAS via pkg-config) + run: cargo test --verbose --features O3-openblas-system pure-rust: name: Pure-Rust feature set @@ -45,6 +49,20 @@ jobs: - name: Run tests run: cargo test --verbose --features "arrow complex csv indexmap json num-complex parallel parquet rayon rkyv serde" + # The deterministic core (no `rand` stack) must keep building for sandboxed + # targets such as Typst / wasm plugins (#88). + no-rand-wasm: + name: No-rand core (wasm32) + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Add wasm32 target + run: rustup target add wasm32-unknown-unknown + - name: Test without the rand stack + run: cargo test --verbose --no-default-features + - name: Build for wasm32 + run: cargo build --verbose --no-default-features --target wasm32-unknown-unknown + feature-combinations: name: Feature combinations (cargo-hack) runs-on: ubuntu-latest @@ -61,17 +79,19 @@ jobs: # HDF5 features (nc / netcdf) are excluded because they need an external # toolchain (vendor BLAS, HDF5 headers) or a non-Linux platform rather than # being a code problem; plot / pyo3 are exercised by the dedicated plot job. + # openblas-src alone (default-features = false) has no TLS backend for its + # source download, so it is only meaningful through O3-openblas-system. - name: Each feature builds on its own run: > cargo hack build --each-feature --optional-deps --keep-going - --exclude-features O3-accelerate,O3-mkl,O3-netlib,nc,netcdf,plot,pyo3,blas-src,lapack-src + --exclude-features O3-accelerate,O3-mkl,O3-netlib,nc,netcdf,plot,pyo3,blas-src,lapack-src,openblas-src # Pairwise (depth 2) coverage over the pure-Rust / numerical features. The # heavy compile units (arrow / parquet) and the backend / HDF5 / Python # features are excluded so the target dir does not blow up on a CI runner. - name: Pairwise feature powerset (pure-Rust / numerical) run: > cargo hack build --feature-powerset --depth 2 --optional-deps --keep-going - --exclude-features O3,O3-openblas,O3-accelerate,O3-mkl,O3-netlib,nc,netcdf,plot,pyo3,blas-src,lapack-src,arrow,parquet + --exclude-features O3,O3-openblas,O3-openblas-system,O3-accelerate,O3-mkl,O3-netlib,nc,netcdf,plot,pyo3,blas-src,lapack-src,openblas-src,arrow,parquet plot: name: plot (pyo3 + matplotlib) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index eee2e8f..f62e2ea 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -41,21 +41,8 @@ Peroxide follows the [Gitflow workflow]. A few practical rules: ## Source layout -A high-level map of `src/`; see each module's `mod.rs` and the [API docs](https://docs.rs/peroxide) for details. - -| Module | Purpose | -| ------------------ | ------------------------------------------------------------------ | -| [`structure`](src/structure) | Core data structures: `Matrix`, `Vec` extensions, `DataFrame`, `Polynomial`, `Jet` forward AD | -| [`numerical`](src/numerical) | Numerical algorithms: ODE solvers, integration, interpolation, splines, root finding, optimization, eigenvalues | -| [`statistics`](src/statistics) | Probability distributions, RNG wrappers, ordered statistics | -| [`complex`](src/complex) | Complex vectors, matrices, and integrals (`complex` feature) | -| [`special`](src/special) | Special functions (wrapper of `puruspe`) | -| [`traits`](src/traits) | Shared trait definitions (math, functional programming, pointers) | -| [`macros`](src/macros) | R / MATLAB / Julia style macros | -| [`fuga`](src/fuga), [`prelude`](src/prelude) | The two user-facing import styles (explicit vs simple) | -| [`util`](src/util) | Constructors, printing, plotting, low-level helpers | -| [`ml`](src/ml) | Basic machine learning tools (beta) | -| [`grave`](src/grave) | Retired implementations kept for reference; not compiled | +The directories under `src/` map one-to-one to the public modules, so the module index in the [API docs](https://docs.rs/peroxide) doubles as the source map; start from a module's documentation and its `mod.rs`. +The only exception is `src/grave/`, which holds retired implementations kept for reference and is excluded from compilation. ## Code of conduct diff --git a/Cargo.toml b/Cargo.toml index da0ad83..c18f8d2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "peroxide" -version = "0.42.0" +version = "0.43.0" authors = ["axect "] edition = "2018" description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with familiar syntax" @@ -30,8 +30,8 @@ criterion = { version = "0.5.1", features = ["html_reports"] } [dependencies] csv = { version = "1.3", optional = true, default-features = false } -rand = { version = "0.9", features = ["small_rng"] } -rand_distr = "0.5" +rand = { version = "0.9", features = ["small_rng"], optional = true } +rand_distr = { version = "0.5", optional = true } order-stat = "0.1" puruspe = "0.4" matrixmultiply = { version = "0.3", features = ["threading"] } @@ -47,6 +47,9 @@ blas = { version = "0.22", optional = true } lapack = { version = "0.19", optional = true } blas-src = { version = "0.14", optional = true, default-features = false } lapack-src = { version = "0.13", optional = true, default-features = false } +# Only used to forward the `system` feature to the `openblas-src` instance +# that `blas-src` / `lapack-src` already pull in (see `O3-openblas-system`). +openblas-src = { version = "0.10", optional = true, default-features = false } serde = { version = "1.0", features = ["derive"], optional = true } rkyv = { version = "0.8", optional = true } json = { version = "0.12", optional = true } @@ -54,6 +57,7 @@ parquet = { version = "55", features = ["arrow", "snap"], optional = true } arrow = { version = "55", optional = true } indexmap = { version = "1", optional = true } num-complex = { version = "0.4", optional = true } +num-traits = { version = "0.2", optional = true } rayon = { version = "1.10", optional = true } [package.metadata.docs.rs] @@ -75,18 +79,44 @@ features = [ "serde", ] +# Examples that exercise the random sampling stack are skipped when the +# `rand` feature is disabled. +[[example]] +name = "dist" +required-features = ["rand"] + +[[example]] +name = "matmul" +required-features = ["rand"] + +[[example]] +name = "optim" +required-features = ["rand"] + +[[example]] +name = "clippy_verify" +required-features = ["rand"] + [features] -default = [] +default = ["rand"] +# Random sampling stack (`rand` / `rand_distr`): probability distributions, +# RNG wrappers, and random constructors. Enabled by default; disable with +# `default-features = false` for the deterministic core (ODE, integration, +# splines, linear algebra), e.g. for wasm32 sandboxes without IO. +rand = ["dep:rand", "dep:rand_distr"] # Bare BLAS / LAPACK FFI; downstream binary must also pull in a # `blas-src` / `lapack-src` backend to resolve the link symbols. O3 = ["blas", "lapack"] # Convenience flags that bundle a backend choice for `O3`. O3-openblas = ["O3", "blas-src/openblas", "lapack-src/openblas"] +# Same as `O3-openblas`, but links the OpenBLAS already installed on the +# system (found via pkg-config) instead of compiling OpenBLAS from source. +O3-openblas-system = ["O3-openblas", "openblas-src/system"] O3-accelerate = ["O3", "blas-src/accelerate", "lapack-src/accelerate"] O3-mkl = ["O3", "blas-src/intel-mkl-dynamic-parallel", "lapack-src/intel-mkl-dynamic-parallel"] O3-netlib = ["O3", "blas-src/netlib", "lapack-src/netlib"] plot = ["pyo3"] nc = ["netcdf"] parquet = ["dep:parquet", "arrow", "indexmap"] -complex = ["num-complex", "matrixmultiply/cgemm"] +complex = ["num-complex", "num-traits", "matrixmultiply/cgemm"] parallel = ["rayon"] diff --git a/README.md b/README.md index 1bb249d..b3f7ad7 100644 --- a/README.md +++ b/README.md @@ -66,7 +66,7 @@ For accelerated linear algebra, plotting, or DataFrame I/O, enable the matching Peroxide provides various features. - `default` - Pure Rust (No dependencies of architecture - Perfect cross compilation) -- `O3` - BLAS & LAPACK (Perfect performance but little bit hard to set-up - Strongly recommend to look [Peroxide with BLAS](https://github.com/Axect/Peroxide_BLAS)) +- `O3-openblas` / `O3-openblas-system` / `O3-accelerate` / `O3-mkl` / `O3-netlib` - BLAS & LAPACK accelerated linear algebra; pick one backend flag (see [Pre-requisite](#pre-requisite)) - `plot` - With matplotlib of python, we can draw any plots. - `complex` - With complex numbers (vector, matrix and integral) - `parallel` - With some parallel functions @@ -76,7 +76,7 @@ Peroxide provides various features. - `serde` - serialization with [Serde](https://serde.rs/). - `rkyv` - serialization with [rkyv](https://rkyv.org). -If you want to do high performance computation and more linear algebra, then choose `O3` feature. +If you want to do high performance computation and more linear algebra, then choose one of the `O3-*` backend features. If you don't want to depend C/C++ or Fortran libraries, then choose `default` feature. If you want to draw plot with some great templates, then choose `plot` feature. @@ -267,25 +267,27 @@ The three groups below depend on external libraries or runtimes; install the rel Those crates only provide function signatures, so the link backend that supplies the actual `dgemv_` / `dpotrf_` / ... symbols must be selected separately. The simplest path is to enable one of the convenience flags below; each pulls in [`blas-src`](https://crates.io/crates/blas-src) and [`lapack-src`](https://crates.io/crates/lapack-src) with the matching backend. -| Convenience flag | Backend | Typical platform / use case | -| ----------------- | ------------------ | --------------------------------------- | -| `O3-openblas` | OpenBLAS | Linux, Windows, macOS via Homebrew | -| `O3-accelerate` | Apple Accelerate | macOS (no extra system install) | -| `O3-mkl` | Intel MKL | Intel CPUs, vendor-tuned performance | -| `O3-netlib` | Netlib reference | Portability, lowest performance | +| Convenience flag | Backend | Build-time requirements | +| -------------------- | ----------------------------- | ------------------------------------------------ | +| `O3-openblas` | OpenBLAS, compiled from source | C + Fortran toolchain, `make`, network access | +| `O3-openblas-system` | System-installed OpenBLAS | `pkg-config` + the OpenBLAS system package | +| `O3-accelerate` | Apple Accelerate | macOS only (no extra system install) | +| `O3-mkl` | Intel MKL | Intel's redistributable (fetched automatically) | +| `O3-netlib` | Netlib reference, compiled from source | `cmake` + Fortran toolchain (lowest performance) | -If you need a backend not in the list above (for example BLIS or R's BLAS), enable the bare `O3` flag and add `blas-src` / `lapack-src` to your downstream binary's `Cargo.toml` with the appropriate features yourself. +`O3-openblas` does **not** use a system-installed OpenBLAS: the [`openblas-src`](https://crates.io/crates/openblas-src) crate downloads the OpenBLAS source tarball during the cargo build and compiles it, so it needs `gcc`, `gfortran`, `make`, and network access, but no BLAS system package. +The download happens over HTTPS through `openblas-src`'s default `rustls` TLS backend; if you depend on `openblas-src` directly with `default-features = false` (as some older guides suggest), you must re-enable one of its `rustls` / `native-tls` features yourself or the build will fail. -System libraries still need to be present on the host for `O3-openblas` and `O3-netlib`; install them with: +`O3-openblas-system` skips the source build and links the OpenBLAS already installed on the host, discovered via `pkg-config`: | Platform | Install | | --------------------- | ---------------------------------------------------- | -| Debian / Ubuntu | `sudo apt install libopenblas-dev liblapack-dev` | -| Fedora / RHEL | `sudo dnf install openblas-devel lapack-devel` | -| Arch Linux | `sudo pacman -S openblas lapack` | -| macOS (Homebrew) | `brew install openblas lapack` | +| Debian / Ubuntu | `sudo apt install libopenblas-dev` | +| Fedora / RHEL | `sudo dnf install openblas-devel` | +| Arch Linux | `sudo pacman -S openblas` | +| macOS (Homebrew) | `brew install openblas` | -`O3-accelerate` and `O3-mkl` ship their own backend (Apple's framework and Intel's redistributable, respectively), so they need no further system packages. +If you need a backend not in the list above (for example BLIS or R's BLAS), enable the bare `O3` flag and add `blas-src` / `lapack-src` to your downstream binary's `Cargo.toml` with the appropriate features yourself. > **Note:** `O3-accelerate` only builds on Apple targets. Enabling it on Linux or Windows fails while compiling `accelerate-src` with ``error: library kind `framework` is only supported on Apple targets``; pick `O3-openblas`, `O3-mkl`, or `O3-netlib` instead. For the same reason, exclude `O3-accelerate` (and `O3-mkl` / `O3-netlib` unless their toolchains are installed) when running tools like `cargo hack --each-feature` on Linux. @@ -335,6 +337,7 @@ cargo add peroxide --features "" # opt-in features | Goal | Command | | ------------------------------------------------- | ------------------------------------------------------------------------- | | Linear algebra on Linux / Windows | `cargo add peroxide --features O3-openblas` | +| Linear algebra with system OpenBLAS | `cargo add peroxide --features O3-openblas-system` | | Linear algebra on macOS | `cargo add peroxide --features O3-accelerate` | | Plotting via Python / matplotlib | `cargo add peroxide --features plot` | | DataFrame + Parquet I/O | `cargo add peroxide --features parquet` | @@ -350,7 +353,8 @@ The remaining single-crate flags exist so advanced users can pull in just one op | Flag | Requires | Purpose | | ---------------- | ----------------------- | ------------------------------------------------------------- | -| `O3-openblas` | OpenBLAS | BLAS / LAPACK accelerated linear algebra (Linux / Windows) | +| `O3-openblas` | OpenBLAS (from source) | BLAS / LAPACK accelerated linear algebra (Linux / Windows) | +| `O3-openblas-system` | OpenBLAS (system) | Same, linking the system-installed OpenBLAS via pkg-config | | `O3-accelerate` | Apple Accelerate | Same, using the Accelerate framework on macOS | | `O3-mkl` | Intel MKL | Same, using Intel MKL | | `O3-netlib` | Netlib | Same, using the reference Netlib BLAS | @@ -363,6 +367,7 @@ The remaining single-crate flags exist so advanced users can pull in just one op | `json` | (pure Rust) | JSON I/O for `DataFrame` | | `serde` | (pure Rust) | `serde` (de)serialization | | `rkyv` | (pure Rust) | `rkyv` zero-copy (de)serialization | +| `rand` | (pure Rust) | Random sampling stack: distributions, RNG wrappers, `rand()` constructors. **On by default**; disable with `default-features = false` to get the deterministic core (ODE, integration, splines, linear algebra) for sandboxed targets like wasm32 |
Advanced: single-crate flags diff --git a/RELEASES.md b/RELEASES.md index 2b3c166..cdbd7f2 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -1,3 +1,28 @@ +# Release 0.43.0 (2026-07-11) + +## Breaking changes + +- Make the `rand` / `rand_distr` sampling stack optional behind the default-on `rand` feature ([#88](https://github.com/Axect/Peroxide/issues/88), [#104](https://github.com/Axect/Peroxide/pull/104)) + - Existing default-feature users are unchanged. + - `default-features = false` now provides a deterministic core without RNG dependencies; sampling APIs require the `rand` feature. + - This is a breaking change for users who already disabled default features and relied on sampling APIs. + +## New features + +- Add the `Dirichlet(α)` probability distribution ([#95](https://github.com/Axect/Peroxide/pull/95)) +- Add `O3-openblas-system` for linking a system-installed OpenBLAS through `pkg-config` ([#98](https://github.com/Axect/Peroxide/issues/98), [#107](https://github.com/Axect/Peroxide/pull/107)) + +## Documentation + +- Clarify BLAS/LAPACK backend selection, OpenBLAS source versus system builds, TLS prerequisites, and HDF5 constraints ([#98](https://github.com/Axect/Peroxide/issues/98)) +- Remove the stale `Peroxide_BLAS` setup link from the main README; the archived repository is retained for historical reference +- Replace the hand-maintained source-layout table with module-level docs.rs pointers and improve module descriptions ([#99](https://github.com/Axect/Peroxide/issues/99), [#108](https://github.com/Axect/Peroxide/pull/108)) + +## CI / Lint + +- Add cargo-hack coverage for individual features and pairwise pure-Rust feature combinations ([#98](https://github.com/Axect/Peroxide/issues/98)) +- Add dedicated CI coverage for system/source OpenBLAS, no-rand `wasm32`, plotting, formatting, and clippy + # Release 0.42.0 (2026-07-06) ## Breaking changes diff --git a/TODO.md b/TODO.md index 4842dcf..b728937 100644 --- a/TODO.md +++ b/TODO.md @@ -21,7 +21,7 @@ - [x] Bernoulli - [x] Beta - [x] Binomial - - [ ] Dirichlet + - [x] Dirichlet - [x] Gamma - [x] Student's t - [x] Uniform diff --git a/src/complex/matrix.rs b/src/complex/matrix.rs index 6162fca..f126843 100644 --- a/src/complex/matrix.rs +++ b/src/complex/matrix.rs @@ -7,8 +7,8 @@ use std::{ use anyhow::{bail, Result}; use matrixmultiply::CGemmOption; use num_complex::Complex; +use num_traits::{One, Zero}; use peroxide_num::{ExpLogOps, PowOps, TrigOps}; -use rand_distr::num_traits::{One, Zero}; use crate::{ complex::C64, diff --git a/src/complex/mod.rs b/src/complex/mod.rs index b17fc18..16a7e62 100644 --- a/src/complex/mod.rs +++ b/src/complex/mod.rs @@ -1,3 +1,5 @@ +//! Complex vectors, matrices, and integrals, enabled by the `complex` feature + use num_complex::Complex; pub type C64 = Complex; diff --git a/src/fuga/mod.rs b/src/fuga/mod.rs index 892647f..321e8c6 100644 --- a/src/fuga/mod.rs +++ b/src/fuga/mod.rs @@ -183,10 +183,17 @@ pub use crate::complex::{integral::*, matrix::*, vector::*, C64}; #[allow(ambiguous_glob_reexports)] pub use crate::structure::{ad::*, dataframe::*, matrix::*, polynomial::*, vector::*}; -pub use crate::util::{api::*, low_level::*, non_macro::*, print::*, useful::*, wrapper::*}; +pub use crate::util::{api::*, low_level::*, non_macro::*, print::*, useful::*}; +#[cfg(feature = "rand")] +pub use crate::util::wrapper::*; + +#[allow(unused_imports)] +pub use crate::statistics::{ops::*, stat::*}; + +#[cfg(feature = "rand")] #[allow(unused_imports)] -pub use crate::statistics::{dist::*, ops::*, rand::*, stat::*}; +pub use crate::statistics::{dist::*, rand::*}; #[allow(unused_imports)] pub use crate::special::function::*; @@ -205,6 +212,7 @@ pub use crate::util::plot::*; pub use anyhow; pub use paste; +#[cfg(feature = "rand")] pub use rand::prelude::*; // ============================================================================= diff --git a/src/lib.rs b/src/lib.rs index cab7066..995dce9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -140,7 +140,7 @@ //! //! ## Useful tips for features //! -//! * If you want to use _QR_, _SVD_, or _Cholesky Decomposition_, you should use the `O3` feature. These decompositions are not implemented in the `default` feature. +//! * If you want to use _QR_, _SVD_, or _Cholesky Decomposition_, enable one of the BLAS backend features: `O3-openblas` (compiles OpenBLAS from source), `O3-openblas-system` (links the system-installed OpenBLAS), `O3-accelerate`, `O3-mkl`, or `O3-netlib`. These decompositions are not implemented in the `default` feature. The bare `O3` flag is an advanced escape hatch that expects you to supply a `blas-src` / `lapack-src` backend yourself. //! //! * If you want to save your numerical results, consider using the `parquet` or `nc` features, which correspond to the `parquet` and `netcdf` file formats, respectively. These formats are much more efficient than `csv` and `json`. //! @@ -165,6 +165,10 @@ extern crate lapack; extern crate blas_src as _; #[cfg(feature = "lapack-src")] extern crate lapack_src as _; +// `O3-openblas-system` re-exposes `openblas-src` directly so its `system` +// feature can be forwarded; link the crate here for the same reason. +#[cfg(feature = "openblas-src")] +extern crate openblas_src as _; #[cfg(feature = "plot")] extern crate pyo3; @@ -172,6 +176,7 @@ extern crate pyo3; #[cfg(feature = "serde")] extern crate serde; +#[cfg(feature = "rand")] extern crate rand; // extern crate json; diff --git a/src/macros/matlab_macro.rs b/src/macros/matlab_macro.rs index c643cc8..c9c2c5e 100644 --- a/src/macros/matlab_macro.rs +++ b/src/macros/matlab_macro.rs @@ -39,6 +39,7 @@ macro_rules! zeros { /// println!("{}", a); // 2 x 2 random matrix (0 ~ 1) /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! rand { () => {{ diff --git a/src/macros/mod.rs b/src/macros/mod.rs index 7b97f1f..20ca775 100644 --- a/src/macros/mod.rs +++ b/src/macros/mod.rs @@ -1,4 +1,4 @@ -//! Useful macros +//! R / MATLAB / Julia style macros pub mod julia_macro; pub mod matlab_macro; diff --git a/src/macros/r_macro.rs b/src/macros/r_macro.rs index 6742e3e..28bf8c1 100644 --- a/src/macros/r_macro.rs +++ b/src/macros/r_macro.rs @@ -269,6 +269,7 @@ macro_rules! rbind { /// println!("{:?}", b); /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! runif { ( $x0:expr, $start:expr, $end:expr ) => {{ @@ -304,6 +305,7 @@ macro_rules! runif { /// println!("{:?}", b); /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! rnorm { ( $n:expr, $mean:expr, $sd:expr ) => {{ @@ -333,6 +335,7 @@ macro_rules! rnorm { /// println!("{:?}", b); /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! dnorm { ( $x:expr, $mean: expr, $sd:expr ) => {{ @@ -364,6 +367,7 @@ macro_rules! dnorm { /// println!("{:?}", b); /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! pnorm { ( $x:expr, $mean:expr, $sd:expr ) => {{ @@ -392,6 +396,7 @@ macro_rules! pnorm { /// println!("{:?}", a); /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! rt { ( $n:expr, $df:expr ) => {{ @@ -414,6 +419,7 @@ macro_rules! rt { /// println!("{:?}", a); /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! dt { ( $x:expr, $df:expr ) => {{ @@ -437,6 +443,7 @@ macro_rules! dt { /// println!("{:?}", a); // 0.5 /// } /// ``` +#[cfg(feature = "rand")] #[macro_export] macro_rules! pt { ( $x:expr, $df:expr ) => {{ diff --git a/src/ml/mod.rs b/src/ml/mod.rs index e08e851..d63c3a2 100644 --- a/src/ml/mod.rs +++ b/src/ml/mod.rs @@ -1,3 +1,3 @@ -//! Machine learning tools +//! Basic machine learning tools (beta) pub mod reg; diff --git a/src/numerical/mod.rs b/src/numerical/mod.rs index 0a21bfc..e8bfe54 100644 --- a/src/numerical/mod.rs +++ b/src/numerical/mod.rs @@ -1,4 +1,4 @@ -//! Differential equations & Numerical Analysis tools +//! Numerical algorithms: ODE solvers, quadrature, interpolation, splines, root finding, optimization, and eigenvalue computation pub mod eigen; pub mod integral; diff --git a/src/numerical/optimize.rs b/src/numerical/optimize.rs index b4eac36..2e94ccc 100644 --- a/src/numerical/optimize.rs +++ b/src/numerical/optimize.rs @@ -53,6 +53,7 @@ //! use peroxide::fuga::*; //! //! fn main() { +//! # #[cfg(feature = "rand")] { //! // To prepare noise //! let normal = Normal(0f64, 0.1f64); //! let normal2 = Normal(0f64, 100f64); @@ -94,6 +95,7 @@ //! // .set_marker(vec![Point, Line]) //! // .savefig().expect("Can't draw a plot"); //! //} +//! # } //! } //! //! fn quad(x: &Vec, n: Vec) -> Option> { diff --git a/src/prelude/mod.rs b/src/prelude/mod.rs index 6c46b34..823af04 100644 --- a/src/prelude/mod.rs +++ b/src/prelude/mod.rs @@ -191,10 +191,18 @@ pub use crate::complex::{integral::*, matrix::*, vector::*, C64}; pub use simpler::{solve, SimplerLinearAlgebra}; #[allow(unused_imports)] -pub use crate::util::{api::*, low_level::*, non_macro::*, print::*, useful::*, wrapper::*}; +pub use crate::util::{api::*, low_level::*, non_macro::*, print::*, useful::*}; +#[cfg(feature = "rand")] #[allow(unused_imports)] -pub use crate::statistics::{dist::*, ops::*, rand::*, stat::*}; +pub use crate::util::wrapper::*; + +#[allow(unused_imports)] +pub use crate::statistics::{ops::*, stat::*}; + +#[cfg(feature = "rand")] +#[allow(unused_imports)] +pub use crate::statistics::{dist::*, rand::*}; #[allow(unused_imports)] pub use crate::special::function::{ @@ -228,4 +236,5 @@ pub use crate::util::plot::*; pub use anyhow; pub use paste; +#[cfg(feature = "rand")] pub use rand::prelude::*; diff --git a/src/statistics/dist.rs b/src/statistics/dist.rs index 7e8e4d6..39afbfc 100644 --- a/src/statistics/dist.rs +++ b/src/statistics/dist.rs @@ -13,12 +13,14 @@ //! * Uniform //! * Weighted Uniform //! * Log Normal -//! * There are two enums to represent probability distribution +//! * There are three enums to represent probability distribution //! * `OPDist` : One parameter distribution (Bernoulli) //! * `TPDist` : Two parameter distribution (Uniform, Normal, Beta, Gamma) +//! * `MVDist` : Multivariate distribution (Dirichlet) //! * `T: PartialOrd + SampleUniform + Copy + Into` //! * There are some traits for pdf -//! * `RNG` trait - extract sample & calculate pdf +//! * `RNG` trait - extract sample & calculate pdf for 1D distributions +//! * `MVRNG` trait - extract sample & calculate pdf for multivariate distributions //! * `Statistics` trait - already shown above //! //! ### `RNG` trait @@ -239,6 +241,57 @@ //! * Mean: $e^{\mu + \frac{\sigma^2}{2}}$ //! * Var: $(e^{\sigma^2} - 1)e^{2\mu + \sigma^2}$ //! * To generate log-normal random samples, Peroxide uses the `rand_distr::LogNormal` distribution from the `rand_distr` crate. +//! ### `MVRNG` trait +//! +//! * `MVRNG` trait is composed of four fields +//! * `sample`: Extract samples +//! * `sample_with_rng`: Extract samples with specific rng +//! * `pdf` : Calculate pdf value at specific point +//! * `ln_pdf` : Calculate log-pdf value at specific point +//! ```no_run +//! use rand::Rng; +//! use peroxide::fuga::*; +//! +//! pub trait MVRNG { +//! /// Extract samples of multivariate distributions +//! fn sample(&self, n: usize) -> Matrix; +//! +//! /// Extract samples of distributions with specific rng +//! fn sample_with_rng(&self, rng: &mut R, n: usize) -> Matrix; +//! +//! /// Probability Density Function +//! fn pdf(&self, x: &[f64]) -> f64; +//! +//! /// Log Probability Density Function +//! fn ln_pdf(&self, x: &[f64]) -> f64; +//! } +//! ``` +//! +//! ### Dirichlet Distribution +//! +//! * Definition +//! $$ \text{Dir}(\mathbf{x} | \boldsymbol{\alpha}) = \frac{1}{\text{B}(\boldsymbol{\alpha})} \prod_{i=1}^K x_i^{\alpha_i - 1} $$ +//! where $\text{B}(\boldsymbol{\alpha}) = \frac{\prod_{i=1}^K \Gamma(\alpha_i)}{\Gamma(\sum_{i=1}^K \alpha_i)}$ +//! * Representative value +//! * Mean: $\frac{\alpha_i}{\alpha_0}$ +//! * Var : $\frac{\alpha_i(\alpha_0 - \alpha_i)}{\alpha_0^2(\alpha_0 + 1)}$ +//! * To generate Dirichlet random samples, Peroxide generates $K$ independent Gamma samples and normalizes them. +//! * **Caution**: `MVDist` utilizes the existing `Statistics` trait but outputs vectors and matrices. +//! +//! ```rust +//! use peroxide::fuga::*; +//! +//! fn main() { +//! let mut rng = smallrng_from_seed(42); +//! let a = Dirichlet(vec![1.0, 2.0, 3.0]); // Dir(x | 1.0, 2.0, 3.0) +//! a.sample(100).print(); // Generate 100 samples +//! a.sample_with_rng(&mut rng, 100).print(); // Generate 100 samples with specific rng +//! a.pdf(&[0.16, 0.33, 0.51]).print(); // Probability density +//! a.mean().print(); // Mean vector +//! a.var().print(); // Variance vector +//! a.cov().print(); // Covariance matrix +//! } +//! ``` extern crate rand; extern crate rand_distr; @@ -246,6 +299,7 @@ use rand_distr::weighted::WeightedAliasIndex; use self::rand::prelude::*; use self::rand_distr::uniform::SampleUniform; +pub use self::MVDist::*; pub use self::OPDist::*; pub use self::TPDist::*; use crate::special::function::*; @@ -253,6 +307,7 @@ use crate::traits::fp::FPVector; //use statistics::rand::ziggurat; use self::WeightedUniformError::*; use crate::statistics::{ops::C, stat::Statistics}; +use crate::structure::matrix::{matrix, Matrix, Row}; use crate::util::non_macro::{linspace, seq}; use crate::util::useful::{auto_zip, find_interval}; use anyhow::{bail, Result}; @@ -283,6 +338,15 @@ pub enum TPDist> { LogNormal(T, T), } +/// Multivariate Distribution +/// +/// # Distributions +/// * `Dirichlet(alpha)`: Dirichlet distribution +#[derive(Debug, Clone)] +pub enum MVDist> { + Dirichlet(Vec), +} + pub struct WeightedUniform> { weights: Vec, sum: T, @@ -1000,3 +1064,154 @@ impl Statistics for WeightedUniform { vec![1f64] } } + +/// Multivariate Random Number Generator Trait +pub trait MVRNG { + /// Extract samples of multivariate distributions (Returns an n x k Matrix) + fn sample(&self, n: usize) -> Matrix { + let mut rng = rand::rng(); + self.sample_with_rng(&mut rng, n) + } + + /// Extract samples of distributions with specific rng + fn sample_with_rng(&self, rng: &mut R, n: usize) -> Matrix; + + /// Probability Density Function + fn pdf(&self, x: &[f64]) -> f64 { + self.ln_pdf(x).exp() + } + + /// Log Probability Density Function + fn ln_pdf(&self, x: &[f64]) -> f64; +} + +impl> Statistics for MVDist { + type Array = Matrix; + type Value = Vec; + + fn mean(&self) -> Self::Value { + match self { + MVDist::Dirichlet(alpha_t) => { + let alpha: Vec = alpha_t.iter().map(|&a| a.into()).collect(); + let alpha0: f64 = alpha.iter().sum(); + alpha.iter().map(|&a| a / alpha0).collect() + } + } + } + + fn var(&self) -> Self::Value { + match self { + MVDist::Dirichlet(alpha_t) => { + let alpha: Vec = alpha_t.iter().map(|&a| a.into()).collect(); + let alpha0: f64 = alpha.iter().sum(); + let norm = alpha0.powi(2) * (alpha0 + 1.0); + alpha.iter().map(|&a| a * (alpha0 - a) / norm).collect() + } + } + } + + fn sd(&self) -> Self::Value { + self.var().into_iter().map(|v| v.sqrt()).collect() + } + + fn cov(&self) -> Self::Array { + match self { + MVDist::Dirichlet(alpha_t) => { + let alpha: Vec = alpha_t.iter().map(|&a| a.into()).collect(); + let alpha0: f64 = alpha.iter().sum(); + let k = alpha.len(); + let norm = alpha0.powi(2) * (alpha0 + 1.0); + let mut cov_data = vec![0f64; k * k]; + + for i in 0..k { + for j in 0..k { + let idx = i * k + j; + if i == j { + cov_data[idx] = alpha[i] * (alpha0 - alpha[i]) / norm; + } else { + cov_data[idx] = -alpha[i] * alpha[j] / norm; + } + } + } + + matrix(cov_data, k, k, Row) + } + } + } + + fn cor(&self) -> Self::Array { + let cov_matrix = self.cov(); + let sd_vec = self.sd(); + let k = sd_vec.len(); + + let mut cor_data = vec![0f64; k * k]; + + for i in 0..k { + for j in 0..k { + let idx = i * k + j; + cor_data[idx] = cov_matrix[(i, j)] / (sd_vec[i] * sd_vec[j]); + } + } + matrix(cor_data, k, k, Row) + } +} + +impl> MVRNG for MVDist { + fn sample_with_rng(&self, rng: &mut R, n: usize) -> Matrix { + match self { + MVDist::Dirichlet(alpha_t) => { + let alpha: Vec = alpha_t.iter().map(|&a| a.into()).collect(); + let k = alpha.len(); + let mut sample_data = vec![0f64; n * k]; + + for i in 0..n { + let mut sum = 0f64; + let mut y = vec![0f64; k]; + + for j in 0..k { + let gamma_dist = rand_distr::Gamma::new(alpha[j], 1.0).unwrap(); + y[j] = gamma_dist.sample(rng); + sum += y[j]; + } + + for j in 0..k { + sample_data[i * k + j] = y[j] / sum; + } + } + + matrix(sample_data, n, k, Row) + } + } + } + + fn ln_pdf(&self, x: &[f64]) -> f64 { + match self { + MVDist::Dirichlet(alpha_t) => { + let alpha: Vec = alpha_t.iter().map(|&a| a.into()).collect(); + assert_eq!( + alpha.len(), + x.len(), + "Arguments must have correct dimensions." + ); + + let mut term = 0f64; + let mut sum_x = 0f64; + let mut sum_alpha_ln_gamma = 0f64; + let mut alpha0 = 0f64; + + for (&x_i, &alpha_i) in x.iter().zip(alpha.iter()) { + assert!(x_i > 0f64 && x_i < 1f64, "Arguments must be in (0, 1)"); + + term += (alpha_i - 1.0) * x_i.ln(); + sum_alpha_ln_gamma += ln_gamma(alpha_i); + sum_x += x_i; + alpha0 += alpha_i; + } + + assert!((sum_x - 1.0).abs() < 1e-4, "Arguments must sum up to 1"); + + term + ln_gamma(alpha0) - sum_alpha_ln_gamma + } + } + } +} diff --git a/src/statistics/mod.rs b/src/statistics/mod.rs index 3012a91..8f2b93b 100644 --- a/src/statistics/mod.rs +++ b/src/statistics/mod.rs @@ -1,11 +1,13 @@ -//! Statistical Modules +//! Probability distributions, random sampling, and statistical operations //! //! * Basic statistical tools - `stat.rs` -//! * Popular distributions - `dist.rs` -//! * Simple Random Number Generator - `rand.rs` +//! * Popular distributions - `dist.rs` (`rand` feature) +//! * Simple Random Number Generator - `rand.rs` (`rand` feature) //! * Basic probabilistic operations - `ops.rs` +#[cfg(feature = "rand")] pub mod dist; pub mod ops; +#[cfg(feature = "rand")] pub mod rand; pub mod stat; diff --git a/src/statistics/stat.rs b/src/statistics/stat.rs index d2160b2..738b73a 100644 --- a/src/statistics/stat.rs +++ b/src/statistics/stat.rs @@ -552,6 +552,7 @@ pub fn cor(v1: &Vec, v2: &Vec) -> f64 { /// use peroxide::fuga::*; /// /// fn main() { +/// # #[cfg(feature = "rand")] { /// let a: Matrix = c!(1,2,3,4,5).into(); /// let noise: Matrix = Normal(0,1).sample(5).into(); /// let b: Matrix = &a + &noise; @@ -560,6 +561,7 @@ pub fn cor(v1: &Vec, v2: &Vec) -> f64 { /// // c[0] /// // r[0] 0.7219 /// // r[1] 0.8058 +/// # } /// } /// ``` pub fn lm(input: &Matrix, target: &Matrix) -> Matrix { diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs index 963d6c3..1e45b6d 100644 --- a/src/structure/matrix.rs +++ b/src/structure/matrix.rs @@ -351,8 +351,10 @@ //! let b = eye(2); //! assert_eq!(b, ml_matrix("1 0;0 1")); //! +//! # #[cfg(feature = "rand")] { //! let c = rand(2, 2); //! c.print(); // Random 2x2 matrix +//! # } //! } //! ``` //! # Linear Algebra diff --git a/src/structure/mod.rs b/src/structure/mod.rs index 666bbee..5951c36 100644 --- a/src/structure/mod.rs +++ b/src/structure/mod.rs @@ -1,8 +1,8 @@ -//! Main structures for peroxide +//! Core data structures: `Matrix`, `Vec` extensions, `DataFrame`, `Polynomial`, and const-generic forward-mode AD (`Jet`) //! -//! * Matrix +//! * Matrix (dense & sparse) //! * Vector -//! * Automatic derivatives +//! * Automatic differentiation (`Jet`) //! * Polynomial //! * DataFrame //! * Multinomial (not yet implemented) diff --git a/src/traits/mod.rs b/src/traits/mod.rs index 67a0d15..a3c5c15 100644 --- a/src/traits/mod.rs +++ b/src/traits/mod.rs @@ -1,3 +1,5 @@ +//! Shared trait definitions: mathematics, functional programming, mutability, pointers, and numeric abstractions + pub mod float; pub mod fp; pub mod general; diff --git a/src/util/mod.rs b/src/util/mod.rs index a5bffe0..97154d2 100644 --- a/src/util/mod.rs +++ b/src/util/mod.rs @@ -1,4 +1,4 @@ -//! Utility - plot, print, pickle and etc. +//! Utilities: constructors, printing, plotting, and low-level helpers pub mod api; pub mod non_macro; @@ -9,5 +9,6 @@ pub mod plot; pub mod low_level; pub mod print; pub mod useful; +#[cfg(feature = "rand")] pub mod wrapper; pub mod writer; diff --git a/src/util/non_macro.rs b/src/util/non_macro.rs index edfd4a2..9a49eda 100644 --- a/src/util/non_macro.rs +++ b/src/util/non_macro.rs @@ -30,7 +30,9 @@ //! - concat //! - cat +#[cfg(feature = "rand")] extern crate rand; +#[cfg(feature = "rand")] use self::rand::prelude::*; use crate::structure::{ matrix::Shape::{Col, Row}, @@ -39,6 +41,7 @@ use crate::structure::{ use crate::traits::float::FloatWithPrecision; use crate::traits::matrix::MatrixTrait; use anyhow::{bail, Result}; +#[cfg(feature = "rand")] use rand_distr::{Distribution, Uniform}; #[derive(Debug, Copy, Clone)] @@ -322,6 +325,7 @@ where /// # Description /// /// Range = from 0 to 1 +#[cfg(feature = "rand")] pub fn rand(r: usize, c: usize) -> Matrix { let mut rng = rand::rng(); rand_with_rng(r, c, &mut rng) @@ -332,6 +336,7 @@ pub fn rand(r: usize, c: usize) -> Matrix { /// # Description /// /// Range = from 0 to 1 +#[cfg(feature = "rand")] pub fn rand_with_rng(r: usize, c: usize, rng: &mut R) -> Matrix { let uniform = Uniform::new_inclusive(0f64, 1f64).unwrap(); rand_with_dist(r, c, rng, uniform) @@ -342,6 +347,7 @@ pub fn rand_with_rng(r: usize, c: usize, rng: &mut R) -> Matrix { /// # Description /// /// Any range +#[cfg(feature = "rand")] pub fn rand_with_dist, R: Rng, D: Distribution>( r: usize, c: usize, diff --git a/src/util/print.rs b/src/util/print.rs index 5ff93bf..f9041e9 100644 --- a/src/util/print.rs +++ b/src/util/print.rs @@ -1,5 +1,6 @@ //! Easy to print any structures +#[cfg(feature = "rand")] use crate::statistics::dist::*; use crate::statistics::stat::ConfusionMatrix; #[allow(unused_imports)] @@ -10,8 +11,11 @@ use crate::structure::{ multinomial::Multinomial, polynomial::Polynomial, }; +#[cfg(feature = "rand")] use rand_distr::uniform::SampleUniform; -use std::fmt::{Debug, LowerExp, UpperExp}; +#[cfg(feature = "rand")] +use std::fmt::Debug; +use std::fmt::{LowerExp, UpperExp}; pub trait Printable { fn print(&self); @@ -355,18 +359,27 @@ impl Printable for Multinomial { // } //} +#[cfg(feature = "rand")] impl> Printable for OPDist { fn print(&self) { println!("{:?}", self); } } +#[cfg(feature = "rand")] impl> Printable for TPDist { fn print(&self) { println!("{:?}", self); } } +#[cfg(feature = "rand")] +impl> Printable for MVDist { + fn print(&self) { + println!("{:?}", self); + } +} + //impl Printable for Number { // fn print(&self) { // println!("{:?}", self) diff --git a/tests/dist.rs b/tests/dist.rs index 883a1d3..2549894 100644 --- a/tests/dist.rs +++ b/tests/dist.rs @@ -1,3 +1,5 @@ +#![cfg(feature = "rand")] + extern crate peroxide; use peroxide::fuga::*; @@ -8,3 +10,22 @@ fn test_binomial() { assert!(nearly_eq(b.mean(), 80f64)); assert!(nearly_eq(b.var(), 16f64)); } + +#[test] +fn test_dirichlet() { + let dir = MVDist::Dirichlet(vec![1.0, 2.0, 3.0]); + dir.sample(10).print(); + + let m = dir.mean(); + assert!(nearly_eq(m[0], 1.0 / 6.0)); + assert!(nearly_eq(m[1], 1.0 / 3.0)); + assert!(nearly_eq(m[2], 0.5)); + + let v = dir.var(); + assert!(nearly_eq(v[0], 5.0 / 252.0)); // 1 * 5 / (36 * 7) + assert!(nearly_eq(v[1], 8.0 / 252.0)); // 2 * 4 / (36 * 7) + assert!(nearly_eq(v[2], 9.0 / 252.0)); // 3 * 3 / (36 * 7) + + let pdf_val = dir.pdf(&[0.33333, 0.33333, 0.33333]); + assert!(nearly_eq(pdf_val, 2.222155556222205)); +} diff --git a/tests/matrix.rs b/tests/matrix.rs index d36aed6..6f0bb0c 100644 --- a/tests/matrix.rs +++ b/tests/matrix.rs @@ -54,6 +54,7 @@ fn test_row() { assert_eq!(a.row(0), c!(1, 2)); } +#[cfg(feature = "rand")] #[test] fn test_print() { let op = Bernoulli(0); diff --git a/tests/weighted_uniform.rs b/tests/weighted_uniform.rs index e720499..e3c36c2 100644 --- a/tests/weighted_uniform.rs +++ b/tests/weighted_uniform.rs @@ -1,3 +1,5 @@ +#![cfg(feature = "rand")] + extern crate peroxide; use peroxide::fuga::*;