Ver 0.34.1

* Modify statistics of Weighted Uniform distribution
* Add more `Vec<f64>` algorithms

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "peroxide"
-version = "0.34.0"
+version = "0.34.1"
 authors = ["axect <[email protected]>"]
 edition = "2018"
 description = "Rust comprehensive scientific computation library contains linear algebra, numerical analysis, statistics and machine learning tools with farmiliar syntax"
@@ -29,7 +29,7 @@ matrixmultiply = { version = "0.3", features = ["threading"] }
 peroxide-ad = "0.3"
 #num-complex = "0.3"
 netcdf = { version = "0.7", optional = true, default_features = false }
-pyo3 = { version = "0.18", optional = true }
+pyo3 = { version = "0.19", optional = true }
 blas = { version = "0.22", optional = true }
 lapack = { version = "0.19", optional = true }
 serde = { version = "1.0", features = ["derive"], optional = true }

RELEASES.md

@@ -1,3 +1,14 @@
+# Release 0.34.1 (2023-08-03)
+
+## Modify Statistics of `WeightedUniform`
+
+* Modify `self.sum` to compatible with definition [Weighted Uniform Distribution](https://axect.github.io/posts/006_prs/#22-weighted-uniform-distribution)
+* Modify `mean` & `var` to compatible with definition
+
+## Add more `Vec<f64>` algorithms
+
+* Implement `arg_min`, `max`, `min`
+
 # Release 0.34.0 (2023-06-13)
 
 ## Change Gauss-Kronrod quadrature 

src/ml/reg.rs

@@ -4,7 +4,7 @@ use crate::structure::polynomial::Polynomial;
 ///
 /// # Type
 ///
-/// (Vec<f64>, Vec<f64>) -> Polynomial
+/// `(Vec<f64>, Vec<f64>) -> Polynomial`
 ///
 /// # Examples
 /// ```

src/numerical/ode.rs

@@ -49,7 +49,7 @@
 //!         * `ImMethod` : Implicit method
 //!             * `BDF` : Backward Euler 1st order (To be fixed)
 //!             * `GL4` : Gauss Legendre 4th order
-//!     * `Environment` : External environment (CubicSpline, Vec<f64>, Matrix or Another external table)
+//!     * `Environment` : External environment (`CubicSpline`, `Vec<f64>`, `Matrix` or Another external table)
 //!
 //!
 //! ### `State<T>` structure

src/numerical/utils.rs

@@ -9,7 +9,7 @@ use crate::util::non_macro::{cat, zeros};
 /// : Exact jacobian matrix using Automatic Differenitation
 ///
 /// # Type
-/// (Vector, F) -> Matrix where F: Fn(&Vec<AD>) -> Vec<AD>
+/// `(F, &Vec<f64>) -> Matrix where F: Fn(&Vec<AD>) -> Vec<AD>`
 ///
 /// # Examples
 /// ```

src/statistics/dist.rs

@@ -70,7 +70,7 @@
 //! ### Uniform Distribution
 //!
 //! * Definition
-//!     $$\text{Unif}(x | a, b) = \begin{cases} \frac{1}{b - a} & x \in [a,b]\\\ 0 & \text{otherwise} \end{cases}$$
+//!     $$\text{Unif}(x | a, b) = \begin{cases} \frac{1}{b - a} & x \in \[a,b\]\\\ 0 & \text{otherwise} \end{cases}$$
 //! * Representative value
 //!     * Mean: $\frac{a + b}{2}$
 //!     * Var : $\frac{1}{12}(b-a)^2$
@@ -126,6 +126,20 @@
 //!
 //! ### Binomial Distribution
 //!
+//! ### Student's t Distribution
+//!
+//! ### Weighted Uniform Distribution
+//!
+//! * Definition
+//!    $$\text{WUnif}(x | \mathbf{W}, \mathcal{I}) = \frac{1}{\sum_{j=1}^n w_j \mu(I_j)} \sum_{i=1}^n w_i
+//!    \mathbb{1}_{I_i}(x)$$
+//!    * $\mathbf{W} = (w_i)$: Weights
+//!    * $\mathcal{I} = \\{I_i\\}$: Intervals
+//!    * $\mu(I_i)$: Measure of $I_i$
+//!    * $\mathbb{1}_{I_i}(x)$: Indicator function
+//!
+//! * Reference
+//!     * [Piecewise Rejection Sampling](https://axect.github.io/posts/006_prs/#22-weighted-uniform-distribution)
 
 extern crate rand;
 extern crate rand_distr;
@@ -208,7 +222,9 @@ impl WeightedUniform<f64> {
             }
         }
 
-        let sum = weights.iter().sum();
+        let sum = weights.iter()
+            .zip(intervals.iter())
+            .fold(0f64, |acc, (w, (a, b))| acc + w * (b - a));
 
         WeightedUniform {
             weights,
@@ -272,7 +288,9 @@ impl WeightedUniform<f64> {
             )
             .collect();
 
-        let sum = weights.iter().sum();
+        let sum = weights.iter()
+            .zip(intervals.iter())
+            .fold(0f64, |acc, (w, (x, y))| acc + w * (y - x));
 
         WeightedUniform {
             weights,
@@ -304,12 +322,17 @@ impl WeightedUniform<f64> {
     pub fn update_weights(&mut self, weights: Vec<f64>) {
         assert_eq!(self.intervals.len(), weights.len());
         self.weights = weights;
-        self.sum = self.weights.iter().sum();
+        self.sum = self.weights.iter()
+            .zip(self.intervals.iter())
+            .fold(0f64, |acc, (w, (a, b))| acc + w * (b - a));
     }
 
     pub fn update_intervals(&mut self, intervals: Vec<f64>) {
         assert_eq!(self.weights.len()+1, intervals.len());
         self.intervals = auto_zip(&intervals);
+        self.sum = self.weights.iter()
+            .zip(self.intervals.iter())
+            .fold(0f64, |acc, (w, (a, b))| acc + w * (b - a));
     }
 
     pub fn weight_at(&self, x: f64) -> f64 {
@@ -665,14 +688,26 @@ impl RNG for WeightedUniform<f64> {
 
     fn pdf<S: PartialOrd + SampleUniform + Copy + Into<f64>>(&self, x: S) -> f64 {
         let x: f64 = x.into();
+        if x < self.intervals[0].0 || x > self.intervals[self.intervals.len() - 1].1 {
+            return 0f64;
+        }
         let idx = find_interval(self.intervals(), x);
         self.weights[idx] / self.sum
     }
 
     fn cdf<S: PartialOrd + SampleUniform + Copy + Into<f64>>(&self, x: S) -> f64 {
         let x: f64 = x.into();
+        if x < self.intervals[0].0 {
+            return 0f64;
+        } else if x > self.intervals[self.intervals.len() - 1].1 {
+            return 1f64;
+        }
         let idx = find_interval(self.intervals(), x);
-        self.weights[0 ..=idx].iter().sum::<f64>() / self.sum
+        self.weights[0 ..=idx].iter()
+            .zip(self.intervals[0 ..=idx].iter())
+            .fold(0f64, |acc, (w, (a, b))| {
+                acc + w * (b - a)
+            }) / self.sum
     }
 }
 
@@ -769,14 +804,14 @@ impl Statistics for WeightedUniform<f64> {
 
     fn mean(&self) -> Self::Value {
         self.intervals().iter().zip(self.weights().iter())
-            .map(|((l, r), w)| (l+r)/2f64 * w)
+            .map(|((l, r), w)| (r.powi(2)-l.powi(2))/2f64 * w)
             .sum::<f64>() / self.sum
     }
 
     fn var(&self) -> Self::Value {
         let mean = self.mean();
         self.intervals().iter().zip(self.weights().iter())
-            .map(|((l, r), w)| w * (l*l + l*r + r*r) / 3f64)
+            .map(|((l, r), w)| w * (r.powi(3) - l.powi(3)) / 3f64)
             .sum::<f64>() / self.sum - mean * mean
     }
 
@@ -791,4 +826,4 @@ impl Statistics for WeightedUniform<f64> {
     fn cor(&self) -> Self::Array {
         vec![1f64]
     }
-}
+}

src/structure/matrix.rs

@@ -321,7 +321,7 @@
 //!     }
 //!     ```
 //!
-//! ## Conversion to Vec<f64>
+//! ## Conversion to `Vec<f64>`
 //!
 //! * Just use `row` or `col` method (I already showed at Basic method section).
 //!
@@ -1718,21 +1718,21 @@ impl MatrixProduct for Matrix {
 // =============================================================================
 // Common Properties of Matrix & Vec<f64>
 // =============================================================================
-/// Matrix to Vec<f64>
+/// `Matrix` to `Vec<f64>`
 impl Into<Vec<f64>> for Matrix {
     fn into(self) -> Vec<f64> {
         self.data
     }
 }
 
-/// &Matrix to &Vec<f64>
+/// `&Matrix` to `&Vec<f64>`
 impl<'a> Into<&'a Vec<f64>> for &'a Matrix {
     fn into(self) -> &'a Vec<f64> {
         &self.data
     }
 }
 
-/// Vec<f64> to Matrix
+/// `Vec<f64>` to `Matrix`
 impl Into<Matrix> for Vec<f64> {
     fn into(self) -> Matrix {
         let l = self.len();
@@ -2374,7 +2374,7 @@ impl<'a, 'b> Mul<&'b Vec<f64>> for &'a Matrix {
     }
 }
 
-/// Matrix multiplication for Vec<f64> vs Matrix
+/// Matrix multiplication for `Vec<f64>` vs `Matrix`
 ///
 /// # Examples
 /// ```

src/structure/vector.rs

@@ -39,7 +39,7 @@
 //!     }
 //!     ```
 //!
-//! ## From ranges to Vec<f64>
+//! ## From ranges to `Vec<f64>`
 //!
 //! * For `R`, there is `seq` to declare sequence.
 //!
@@ -63,7 +63,7 @@
 //!     }
 //!     ```
 //!
-//! ## Vec<f64> Operation
+//! ## `Vec<f64>` Operation
 //!
 //! There are two ways to do vector operations.
 //!
@@ -113,9 +113,9 @@
 //!
 //! ## Conversion to Matrix
 //!
-//! There are two ways to convert Vec<f64> to matrix.
+//! There are two ways to convert `Vec<f64>` to `Matrix`.
 //!
-//! * `into(self) -> Matrix` : Vec<f64> to column matrix
+//! * `into(self) -> Matrix` : `Vec<f64>` to column matrix
 //!
 //!     ```rust
 //!     #[macro_use]
@@ -135,7 +135,7 @@
 //!
 //! # Functional Programming {#functional}
 //!
-//! ## FP for Vec<f64>
+//! ## FP for `Vec<f64>`
 //!
 //! * There are some functional programming tools for `Vec<f64>`
 //!
@@ -282,7 +282,7 @@ use std::convert;
 impl FPVector for Vec<f64> {
     type Scalar = f64;
 
-    /// fmap for Vec<f64>
+    /// fmap for `Vec<f64>`
     ///
     /// # Examples
     /// ```
@@ -304,7 +304,7 @@ impl FPVector for Vec<f64> {
         v
     }
 
-    /// reduce for Vec<f64>
+    /// reduce for `Vec<f64>`
     ///
     /// # Examples
     /// ```
@@ -335,7 +335,7 @@ impl FPVector for Vec<f64> {
             .collect::<Vec<f64>>()
     }
 
-    /// Filter for Vec<f64>
+    /// Filter for `Vec<f64>`
     ///
     /// # Examples
     /// ```
@@ -359,7 +359,7 @@ impl FPVector for Vec<f64> {
             .collect::<Vec<f64>>()
     }
 
-    /// Take for Vec<f64>
+    /// Take for `Vec<f64>`
     ///
     /// # Examples
     /// ```
@@ -381,7 +381,7 @@ impl FPVector for Vec<f64> {
         return v;
     }
 
-    /// Skip for Vec<f64>
+    /// Skip for `Vec<f64>`
     ///
     /// # Examples
     /// ```
@@ -583,6 +583,75 @@ impl Algorithm for Vec<f64> {
         }
     }
 
+    /// arg min
+    ///
+    /// # Examples
+    /// ```
+    /// #[macro_use]
+    /// extern crate peroxide;
+    /// use peroxide::fuga::*;
+    ///
+    /// fn main() {
+    ///     let v = c!(1,3,2,4,3,7);
+    ///     assert_eq!(v.arg_min(),0);
+    ///
+    ///     let v2 = c!(4,3,2,5,1,6);
+    ///     assert_eq!(v2.arg_min(),4);
+    /// }
+    fn arg_min(&self) -> usize {
+        match () {
+            _ => {
+                self.into_iter()
+                    .enumerate()
+                    .fold((0usize, std::f64::MAX), |acc, (ics, &val)| {
+                        if acc.1 > val {
+                            (ics, val)
+                        } else {
+                            acc
+                        }
+                    })
+                    .0
+            }
+        }
+    }
+
+    fn max(&self) -> f64 {
+        match () {
+            #[cfg(feature = "O3")]
+            () => {
+                let n_i32 = self.len() as i32;
+                let idx: usize;
+                unsafe {
+                    idx = idamax(n_i32, self, 1) - 1;
+                }
+                self[idx]
+            }
+            _ => {
+                self.into_iter().fold(std::f64::MIN, |acc, &val| {
+                    if acc < val {
+                        val
+                    } else {
+                        acc
+                    }
+                })
+            }
+        }
+    }
+
+    fn min(&self) -> f64 {
+        match () {
+            _ => {
+                self.into_iter().fold(std::f64::MAX, |acc, &val| {
+                    if acc > val {
+                        val
+                    } else {
+                        acc
+                    }
+                })
+            }
+        }
+    }
+
     fn swap_with_perm(&mut self, p: &Vec<(usize, usize)>) {
         for (i, j) in p.iter() {
             self.swap(*i, *j);

src/traits/general.rs

@@ -3,5 +3,8 @@ pub trait Algorithm {
     fn rank(&self) -> Vec<usize>;
     fn sign(&self) -> f64;
     fn arg_max(&self) -> usize;
+    fn arg_min(&self) -> usize;
+    fn max(&self) -> f64;
+    fn min(&self) -> f64;
     fn swap_with_perm(&mut self, p: &Vec<(usize, usize)>);
 }