From aa635248e13cb34cd62b906fc9b0c9ce765f9248 Mon Sep 17 00:00:00 2001 From: axect Date: Mon, 6 Jul 2026 21:04:57 +0800 Subject: [PATCH 01/11] FEAT: Make rand / rand_distr optional behind default-on rand feature (#88) Sandboxed targets (Typst plugins, wasm32) cannot use stateful RNG, so the random sampling stack is now optional: * rand / rand_distr become optional dependencies behind a new default-on `rand` feature; existing users see no change, and `default-features = false` drops the RNG stack entirely. * Gated behind the feature: statistics::dist, statistics::rand, util::wrapper, the rand() / rand_with_rng() / rand_with_dist() constructors, the runif! / rnorm! / dnorm! / pnorm! / rt! / dt! / pt! / rand! macros, the dist Printable impls, and the fuga / prelude re-exports of all of the above. * complex/matrix.rs imported num_traits through rand_distr's re-export; it now uses num-traits directly (new optional dep tied to the complex feature). * Examples that sample (dist, matmul, optim, clippy_verify) declare required-features = ["rand"]; rand-dependent tests and doctests are cfg-gated. * New CI job: cargo test --no-default-features plus a wasm32-unknown-unknown build, so the deterministic core (ODE, integration, splines, linear algebra) cannot silently regress. Verified locally: default 472 tests, complex 512 tests, no-default 444 tests, wasm32 build, clippy zero errors in both configurations. --- .github/workflows/Github.yml | 14 ++++++++++++++ Cargo.toml | 32 ++++++++++++++++++++++++++++---- README.md | 1 + src/complex/matrix.rs | 2 +- src/fuga/mod.rs | 12 ++++++++++-- src/lib.rs | 1 + src/macros/matlab_macro.rs | 1 + src/macros/r_macro.rs | 7 +++++++ src/numerical/optimize.rs | 2 ++ src/prelude/mod.rs | 13 +++++++++++-- src/statistics/mod.rs | 6 ++++-- src/statistics/stat.rs | 2 ++ src/structure/matrix.rs | 2 ++ src/util/mod.rs | 1 + src/util/non_macro.rs | 6 ++++++ src/util/print.rs | 8 +++++++- tests/dist.rs | 2 ++ tests/matrix.rs | 1 + tests/weighted_uniform.rs | 2 ++ 19 files changed, 103 insertions(+), 12 deletions(-) diff --git a/.github/workflows/Github.yml b/.github/workflows/Github.yml index 638636cb..dd231f1e 100644 --- a/.github/workflows/Github.yml +++ b/.github/workflows/Github.yml @@ -45,6 +45,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 diff --git a/Cargo.toml b/Cargo.toml index da0ad834..c96ff8c5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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"] } @@ -54,6 +54,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,8 +76,31 @@ 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"] @@ -88,5 +112,5 @@ 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 1bb249d9..84faceb6 100644 --- a/README.md +++ b/README.md @@ -363,6 +363,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/src/complex/matrix.rs b/src/complex/matrix.rs index 6162fcad..f126843a 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/fuga/mod.rs b/src/fuga/mod.rs index 892647f5..321e8c63 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 cab7066d..26985205 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -172,6 +172,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 c643cc88..c9c2c5e7 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/r_macro.rs b/src/macros/r_macro.rs index 6742e3e8..28bf8c1c 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/numerical/optimize.rs b/src/numerical/optimize.rs index b4eac36d..2e94cccb 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 6c46b349..823af049 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/mod.rs b/src/statistics/mod.rs index 3012a917..c128b652 100644 --- a/src/statistics/mod.rs +++ b/src/statistics/mod.rs @@ -1,11 +1,13 @@ //! Statistical Modules //! //! * 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 d2160b20..738b73a0 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 963d6c3d..1e45b6da 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/util/mod.rs b/src/util/mod.rs index a5bffe0b..b39f2963 100644 --- a/src/util/mod.rs +++ b/src/util/mod.rs @@ -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 edfd4a2f..9a49eda7 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 5ff93bf8..3ba6b9cf 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,12 +359,14 @@ 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); diff --git a/tests/dist.rs b/tests/dist.rs index 883a1d39..14dcad17 100644 --- a/tests/dist.rs +++ b/tests/dist.rs @@ -1,3 +1,5 @@ +#![cfg(feature = "rand")] + extern crate peroxide; use peroxide::fuga::*; diff --git a/tests/matrix.rs b/tests/matrix.rs index d36aed6d..6f0bb0c5 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 e7204995..e3c36c2d 100644 --- a/tests/weighted_uniform.rs +++ b/tests/weighted_uniform.rs @@ -1,3 +1,5 @@ +#![cfg(feature = "rand")] + extern crate peroxide; use peroxide::fuga::*; From ebf85c582c68d215d575f1b5936fe9677566ae06 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Sun, 19 Apr 2026 19:48:40 +0800 Subject: [PATCH 02/11] =?UTF-8?q?Implement=20Dirichlet(=CE=B1)=20distribut?= =?UTF-8?q?ion?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- TODO.md | 2 +- src/statistics/dist.rs | 206 ++++++++++++++++++++++++++++++++++++++++- tests/dist.rs | 18 ++++ 3 files changed, 223 insertions(+), 3 deletions(-) diff --git a/TODO.md b/TODO.md index 4842dcf4..b7289374 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/statistics/dist.rs b/src/statistics/dist.rs index 7e8e4d6b..2828c953 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,55 @@ //! * 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; +//! pub trait MVRNG { +//! /// Extract samples of multivariate distributions +//! fn sample(&self, n: usize) -> Vec>; +//! +//! /// Extract samples of distributions with specific rng +//! fn sample_with_rng(&self, rng: &mut R, n: usize) -> Vec>; +//! +//! /// 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; @@ -283,6 +334,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 +1060,145 @@ impl Statistics for WeightedUniform { vec![1f64] } } + +/// Multivariate Random Number Generator Trait +pub trait MVRNG { + /// Extract samples of multivariate distributions + fn sample(&self, n: usize) -> Vec> { + 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) -> Vec>; + + /// 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 = Vec>; + 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 { + // Element-wise standard deviation + 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_matrix = vec![vec![0f64; k]; k]; + + for i in 0..k { + for j in 0..k { + if i == j { + cov_matrix[i][j] = alpha[i] * (alpha0 - alpha[i]) / norm; + } else { + cov_matrix[i][j] = -alpha[i] * alpha[j] / norm; + } + } + } + cov_matrix + } + } + } + + fn cor(&self) -> Self::Array { + // Correlation matrix: Cor(X_i, X_j) = Cov(X_i, X_j) / (SD(X_i) * SD(X_j)) + let cov_matrix = self.cov(); + let sd_vec = self.sd(); + let k = sd_vec.len(); + let mut cor_matrix = vec![vec![0f64; k]; k]; + + for i in 0..k { + for j in 0..k { + cor_matrix[i][j] = cov_matrix[i][j] / (sd_vec[i] * sd_vec[j]); + } + } + cor_matrix + } +} + +impl> MVRNG for MVDist { + fn sample_with_rng(&self, rng: &mut R, n: usize) -> Vec> { + match self { + MVDist::Dirichlet(alpha_t) => { + let alpha: Vec = alpha_t.iter().map(|&a| a.into()).collect(); + let mut samples = Vec::with_capacity(n); + + for _ in 0..n { + let mut sum = 0f64; + let mut y = Vec::with_capacity(alpha.len()); + + for &a in &alpha { + let gamma_dist = rand_distr::Gamma::new(a, 1.0).unwrap(); + let val = gamma_dist.sample(rng); + sum += val; + y.push(val); + } + + samples.push(y.into_iter().map(|v| v / sum).collect()); + } + samples + } + } + } + + 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 += gamma(alpha_i).ln(); + sum_x += x_i; + alpha0 += alpha_i; + } + + assert!((sum_x - 1.0).abs() < 1e-4, "Arguments must sum up to 1"); + + term + gamma(alpha0).ln() - sum_alpha_ln_gamma + } + } + } +} diff --git a/tests/dist.rs b/tests/dist.rs index 14dcad17..b76567c4 100644 --- a/tests/dist.rs +++ b/tests/dist.rs @@ -10,3 +10,21 @@ 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]); + + 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)); +} From bd28135e067eb383de14b77931504fd7780f8ea4 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Sun, 19 Apr 2026 19:55:50 +0800 Subject: [PATCH 03/11] Implement `Print` --- src/util/print.rs | 265 ++++++++++++++++++++++++++++++++++++++++++++++ tests/dist.rs | 1 + 2 files changed, 266 insertions(+) diff --git a/src/util/print.rs b/src/util/print.rs index 3ba6b9cf..de6b461e 100644 --- a/src/util/print.rs +++ b/src/util/print.rs @@ -320,6 +320,265 @@ impl Printable for &Vec { } } +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for &Vec> { + fn print(&self) { + println!("{:?}", self); + } +} + +impl Printable for Vec> { + fn print(&self) { + if self.is_empty() { + println!("[]"); + return; + } + let mut result = String::new(); + result.push_str("[\n"); + for i in 0..self.len() { + result.push_str(" "); + result.push_str(&format_float_vec!(self[i])); + if i != self.len() - 1 { + result.push_str(",\n"); + } else { + result.push_str("\n"); + } + } + result.push_str("]"); + println!("{}", result); + } +} + +impl Printable for &Vec> { + fn print(&self) { + if self.is_empty() { + println!("[]"); + return; + } + let mut result = String::new(); + result.push_str("[\n"); + for i in 0..self.len() { + result.push_str(" "); + result.push_str(&format_float_vec!(self[i])); + if i != self.len() - 1 { + result.push_str(",\n"); + } else { + result.push_str("\n"); + } + } + result.push_str("]"); + println!("{}", result); + } +} + +impl Printable for Vec> { + fn print(&self) { + if self.is_empty() { + println!("[]"); + return; + } + let mut result = String::new(); + result.push_str("[\n"); + for i in 0..self.len() { + result.push_str(" "); + result.push_str(&format_float_vec!(self[i])); + if i != self.len() - 1 { + result.push_str(",\n"); + } else { + result.push_str("\n"); + } + } + result.push_str("]"); + println!("{}", result); + } +} + +impl Printable for &Vec> { + fn print(&self) { + if self.is_empty() { + println!("[]"); + return; + } + + let mut result = String::new(); + result.push_str("[\n"); + + for i in 0..self.len() { + result.push_str(" "); + result.push_str(&format_float_vec!(self[i])); + if i != self.len() - 1 { + result.push_str(",\n"); + } else { + result.push_str("\n"); + } + } + + result.push_str("]"); + println!("{}", result); + } +} + impl Printable for Matrix { fn print(&self) { println!("{}", self); @@ -373,6 +632,12 @@ impl> Printable for TPD } } +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 b76567c4..d153b29b 100644 --- a/tests/dist.rs +++ b/tests/dist.rs @@ -14,6 +14,7 @@ fn test_binomial() { #[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)); From f3269f839495451b6748ced3d29b3cd0dd23b81a Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Sun, 19 Apr 2026 21:22:13 +0800 Subject: [PATCH 04/11] use `Matrix` instead --- src/statistics/dist.rs | 56 +++++---- src/util/print.rs | 259 ----------------------------------------- 2 files changed, 31 insertions(+), 284 deletions(-) diff --git a/src/statistics/dist.rs b/src/statistics/dist.rs index 2828c953..de57923f 100644 --- a/src/statistics/dist.rs +++ b/src/statistics/dist.rs @@ -306,6 +306,7 @@ use self::WeightedUniformError::*; use crate::statistics::{ops::C, stat::Statistics}; use crate::util::non_macro::{linspace, seq}; use crate::util::useful::{auto_zip, find_interval}; +use crate::structure::matrix::{matrix, Matrix, Row}; use anyhow::{bail, Result}; use std::f64::consts::E; @@ -1063,14 +1064,14 @@ impl Statistics for WeightedUniform { /// Multivariate Random Number Generator Trait pub trait MVRNG { - /// Extract samples of multivariate distributions - fn sample(&self, n: usize) -> Vec> { + /// 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) -> Vec>; + fn sample_with_rng(&self, rng: &mut R, n: usize) -> Matrix; /// Probability Density Function fn pdf(&self, x: &[f64]) -> f64 { @@ -1082,7 +1083,7 @@ pub trait MVRNG { } impl> Statistics for MVDist { - type Array = Vec>; + type Array = Matrix; type Value = Vec; fn mean(&self) -> Self::Value { @@ -1107,7 +1108,6 @@ impl> Statistics for MVDist } fn sd(&self) -> Self::Value { - // Element-wise standard deviation self.var().into_iter().map(|v| v.sqrt()).collect() } @@ -1118,59 +1118,65 @@ impl> Statistics for MVDist let alpha0: f64 = alpha.iter().sum(); let k = alpha.len(); let norm = alpha0.powi(2) * (alpha0 + 1.0); - let mut cov_matrix = vec![vec![0f64; k]; k]; + 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_matrix[i][j] = alpha[i] * (alpha0 - alpha[i]) / norm; + cov_data[idx] = alpha[i] * (alpha0 - alpha[i]) / norm; } else { - cov_matrix[i][j] = -alpha[i] * alpha[j] / norm; + cov_data[idx] = -alpha[i] * alpha[j] / norm; } } } - cov_matrix + + matrix(cov_data, k, k, Row) } } } fn cor(&self) -> Self::Array { - // Correlation matrix: Cor(X_i, X_j) = Cov(X_i, X_j) / (SD(X_i) * SD(X_j)) let cov_matrix = self.cov(); let sd_vec = self.sd(); let k = sd_vec.len(); - let mut cor_matrix = vec![vec![0f64; k]; k]; + + let mut cor_data = vec![0f64; k * k]; for i in 0..k { for j in 0..k { - cor_matrix[i][j] = cov_matrix[i][j] / (sd_vec[i] * sd_vec[j]); + let idx = i * k + j; + cor_data[idx] = cov_matrix[(i, j)] / (sd_vec[i] * sd_vec[j]); } } - cor_matrix + matrix(cor_data, k, k, Row) } } impl> MVRNG for MVDist { - fn sample_with_rng(&self, rng: &mut R, n: usize) -> Vec> { + 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 mut samples = Vec::with_capacity(n); + let k = alpha.len(); + let mut sample_data = vec![0f64; n * k]; - for _ in 0..n { + for i in 0..n { let mut sum = 0f64; - let mut y = Vec::with_capacity(alpha.len()); + let mut y = vec![0f64; k]; - for &a in &alpha { - let gamma_dist = rand_distr::Gamma::new(a, 1.0).unwrap(); - let val = gamma_dist.sample(rng); - sum += val; - y.push(val); + 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; } - - samples.push(y.into_iter().map(|v| v / sum).collect()); } - samples + + matrix(sample_data, n, k, Row) } } } diff --git a/src/util/print.rs b/src/util/print.rs index de6b461e..ed64dc23 100644 --- a/src/util/print.rs +++ b/src/util/print.rs @@ -320,265 +320,6 @@ impl Printable for &Vec { } } -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for &Vec> { - fn print(&self) { - println!("{:?}", self); - } -} - -impl Printable for Vec> { - fn print(&self) { - if self.is_empty() { - println!("[]"); - return; - } - let mut result = String::new(); - result.push_str("[\n"); - for i in 0..self.len() { - result.push_str(" "); - result.push_str(&format_float_vec!(self[i])); - if i != self.len() - 1 { - result.push_str(",\n"); - } else { - result.push_str("\n"); - } - } - result.push_str("]"); - println!("{}", result); - } -} - -impl Printable for &Vec> { - fn print(&self) { - if self.is_empty() { - println!("[]"); - return; - } - let mut result = String::new(); - result.push_str("[\n"); - for i in 0..self.len() { - result.push_str(" "); - result.push_str(&format_float_vec!(self[i])); - if i != self.len() - 1 { - result.push_str(",\n"); - } else { - result.push_str("\n"); - } - } - result.push_str("]"); - println!("{}", result); - } -} - -impl Printable for Vec> { - fn print(&self) { - if self.is_empty() { - println!("[]"); - return; - } - let mut result = String::new(); - result.push_str("[\n"); - for i in 0..self.len() { - result.push_str(" "); - result.push_str(&format_float_vec!(self[i])); - if i != self.len() - 1 { - result.push_str(",\n"); - } else { - result.push_str("\n"); - } - } - result.push_str("]"); - println!("{}", result); - } -} - -impl Printable for &Vec> { - fn print(&self) { - if self.is_empty() { - println!("[]"); - return; - } - - let mut result = String::new(); - result.push_str("[\n"); - - for i in 0..self.len() { - result.push_str(" "); - result.push_str(&format_float_vec!(self[i])); - if i != self.len() - 1 { - result.push_str(",\n"); - } else { - result.push_str("\n"); - } - } - - result.push_str("]"); - println!("{}", result); - } -} - impl Printable for Matrix { fn print(&self) { println!("{}", self); From 63394424e7f1dd05a53b8031f7e697fabb09a3d3 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Sun, 19 Apr 2026 21:25:19 +0800 Subject: [PATCH 05/11] docs --- src/statistics/dist.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/statistics/dist.rs b/src/statistics/dist.rs index de57923f..96ea8113 100644 --- a/src/statistics/dist.rs +++ b/src/statistics/dist.rs @@ -252,10 +252,10 @@ //! use rand::Rng; //! pub trait MVRNG { //! /// Extract samples of multivariate distributions -//! fn sample(&self, n: usize) -> Vec>; +//! fn sample(&self, n: usize) -> Matrix; //! //! /// Extract samples of distributions with specific rng -//! fn sample_with_rng(&self, rng: &mut R, n: usize) -> Vec>; +//! fn sample_with_rng(&self, rng: &mut R, n: usize) -> Matrix; //! //! /// Probability Density Function //! fn pdf(&self, x: &[f64]) -> f64; From fbf5896254776509cbb2474b6b05c8636b4c73b9 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Tue, 7 Jul 2026 16:24:27 +0800 Subject: [PATCH 06/11] `cargo fmt`, unresolved imports, and `ln_gamma` --- src/statistics/dist.rs | 23 +++++++++++++++-------- tests/dist.rs | 2 +- 2 files changed, 16 insertions(+), 9 deletions(-) diff --git a/src/statistics/dist.rs b/src/statistics/dist.rs index 96ea8113..39afbfc9 100644 --- a/src/statistics/dist.rs +++ b/src/statistics/dist.rs @@ -250,6 +250,8 @@ //! * `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; @@ -297,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::*; @@ -304,9 +307,9 @@ 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 crate::structure::matrix::{matrix, Matrix, Row}; use anyhow::{bail, Result}; use std::f64::consts::E; @@ -1077,7 +1080,7 @@ pub trait MVRNG { fn pdf(&self, x: &[f64]) -> f64 { self.ln_pdf(x).exp() } - + /// Log Probability Density Function fn ln_pdf(&self, x: &[f64]) -> f64; } @@ -1140,7 +1143,7 @@ impl> Statistics for MVDist 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 { @@ -1164,7 +1167,7 @@ impl> MVRNG for MVDist { 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); @@ -1185,7 +1188,11 @@ impl> MVRNG for MVDist { 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."); + assert_eq!( + alpha.len(), + x.len(), + "Arguments must have correct dimensions." + ); let mut term = 0f64; let mut sum_x = 0f64; @@ -1194,16 +1201,16 @@ impl> MVRNG for MVDist { 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 += gamma(alpha_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 + gamma(alpha0).ln() - sum_alpha_ln_gamma + term + ln_gamma(alpha0) - sum_alpha_ln_gamma } } } diff --git a/tests/dist.rs b/tests/dist.rs index d153b29b..2549894d 100644 --- a/tests/dist.rs +++ b/tests/dist.rs @@ -25,7 +25,7 @@ fn test_dirichlet() { 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)); } From eda2336cb7b87299f2683e6800db0c5b3bb34981 Mon Sep 17 00:00:00 2001 From: jzeuzs <75403863+jzeuzs@users.noreply.github.com> Date: Tue, 7 Jul 2026 16:55:39 +0800 Subject: [PATCH 07/11] `#[cfg(feature = "rand")]` --- src/util/print.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/util/print.rs b/src/util/print.rs index ed64dc23..f9041e94 100644 --- a/src/util/print.rs +++ b/src/util/print.rs @@ -373,6 +373,7 @@ impl> Printable for TPD } } +#[cfg(feature = "rand")] impl> Printable for MVDist { fn print(&self) { println!("{:?}", self); From 6fdc95855f4280d273934da2538637add04212d0 Mon Sep 17 00:00:00 2001 From: axect Date: Fri, 10 Jul 2026 13:07:58 +0800 Subject: [PATCH 08/11] IMPL: Add O3-openblas-system feature and correct BLAS setup docs (JOSS jonaspleyer #98) - New O3-openblas-system feature links the system-installed OpenBLAS via pkg-config instead of compiling OpenBLAS from source; openblas-src is added as an optional dependency solely to forward its system feature. - README Pre-requisite section now states that O3-openblas builds OpenBLAS from source (C/Fortran toolchain, make, network access) and documents the rustls / native-tls requirement of openblas-src 0.10.16+ when used with default-features = false. - Drop the stale Peroxide_BLAS link from the README feature list; point to the Pre-requisite section instead. - CI: o3 job also builds and tests O3-openblas-system; cargo-hack excludes the bare openblas-src flag (no TLS backend on its own) and keeps the new composite flag out of the powerset. --- .github/workflows/Github.yml | 14 ++++++++++---- Cargo.toml | 6 ++++++ README.md | 36 ++++++++++++++++++++---------------- src/lib.rs | 6 +++++- 4 files changed, 41 insertions(+), 21 deletions(-) diff --git a/.github/workflows/Github.yml b/.github/workflows/Github.yml index dd231f1e..4719a0f0 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 @@ -75,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/Cargo.toml b/Cargo.toml index c96ff8c5..0839b5d2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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 } @@ -106,6 +109,9 @@ rand = ["dep:rand", "dep:rand_distr"] 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"] diff --git a/README.md b/README.md index 84faceb6..b3f7ad7e 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 | diff --git a/src/lib.rs b/src/lib.rs index 26985205..995dce94 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; From 737611ca576676743f03b702c628d5db2348306b Mon Sep 17 00:00:00 2001 From: axect Date: Fri, 10 Jul 2026 13:51:07 +0800 Subject: [PATCH 09/11] DOCS: Replace source layout table with pointer to docs.rs module index (JOSS jonaspleyer #99) The hand-maintained module table duplicated the module index that rustdoc generates on docs.rs and could only drift out of sync. What remains is the repo-only fact rustdoc cannot show: src/grave/ is retired code excluded from compilation. --- CONTRIBUTING.md | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index eee2e8fa..f62e2ea5 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 From ee6030d35b927757b53017ec019a055f2bae39f7 Mon Sep 17 00:00:00 2001 From: axect Date: Fri, 10 Jul 2026 13:57:37 +0800 Subject: [PATCH 10/11] DOCS: Move module descriptions from CONTRIBUTING table into module docs The removed source layout table carried descriptions that existed nowhere in rustdoc: complex and traits had no module docs at all, and several first lines (Useful macros, Statistical Modules, Main structures) were too vague to serve as an index. Each module now carries its one-line purpose in its own mod.rs, so the docs.rs module index is the single source for the source map and stays in sync with the code by construction. --- src/complex/mod.rs | 2 ++ src/macros/mod.rs | 2 +- src/ml/mod.rs | 2 +- src/numerical/mod.rs | 2 +- src/statistics/mod.rs | 2 +- src/structure/mod.rs | 6 +++--- src/traits/mod.rs | 2 ++ src/util/mod.rs | 2 +- 8 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/complex/mod.rs b/src/complex/mod.rs index b17fc18f..16a7e621 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/macros/mod.rs b/src/macros/mod.rs index 7b97f1f5..20ca7755 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/ml/mod.rs b/src/ml/mod.rs index e08e851d..d63c3a2d 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 0a21bfc3..e8bfe54d 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/statistics/mod.rs b/src/statistics/mod.rs index c128b652..8f2b93bd 100644 --- a/src/statistics/mod.rs +++ b/src/statistics/mod.rs @@ -1,4 +1,4 @@ -//! Statistical Modules +//! Probability distributions, random sampling, and statistical operations //! //! * Basic statistical tools - `stat.rs` //! * Popular distributions - `dist.rs` (`rand` feature) diff --git a/src/structure/mod.rs b/src/structure/mod.rs index 666bbeee..5951c364 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 67a0d15b..a3c5c159 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 b39f2963..97154d24 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; From c9b063eca3afc34278791272fb421014ba697470 Mon Sep 17 00:00:00 2001 From: axect Date: Sat, 11 Jul 2026 18:58:52 +0800 Subject: [PATCH 11/11] RELEASE: Version 0.43.0 --- Cargo.toml | 2 +- RELEASES.md | 25 +++++++++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 0839b5d2..c18f8d2a 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" diff --git a/RELEASES.md b/RELEASES.md index 2b3c166b..cdbd7f23 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