Adds CMA-ES algorithm

.gitignore

@@ -4,3 +4,6 @@ target/
 target/*
 *.log
 justfile
+.vscode
+.DS_Store
+checkpoints/

argmin-math/Cargo.toml

@@ -37,6 +37,7 @@ num-complex_0_2 = { package = "num-complex", version = "0.2", optional = true, d
 num-traits = { version = "0.2" }
 num-integer = { version = "0.1" }
 rand = { version = "0.8.3" }
+rand_distr = { version = "0.4.3" }
 anyhow = { version = "1.0" }
 thiserror = { version = "1.0" }
 

argmin-math/src/lib.rs

@@ -161,8 +161,6 @@ pub use crate::nalgebra_m::*;
 
 #[cfg(feature = "vec")]
 mod vec;
-#[cfg(feature = "vec")]
-pub use crate::vec::*;
 
 use anyhow::Error;
 
@@ -172,6 +170,12 @@ pub trait ArgminDot<T, U> {
     fn dot(&self, other: &T) -> U;
 }
 
+/// Dot/scalar product of `T`.transpose() and `self`
+pub trait ArgminTDot<T, U> {
+    /// Dot/scalar product of `T`.transpose() and `self`
+    fn tdot(&self, other: &T) -> U;
+}
+
 /// Dot/scalar product of `T` and `self` weighted by W (p^TWv)
 pub trait ArgminWeightedDot<T, U, V> {
     /// Dot/scalar product of `T` and `self`
@@ -275,3 +279,108 @@ pub trait ArgminMinMax {
     /// Select piecewise maximum
     fn max(x: &Self, y: &Self) -> Self;
 }
+
+/// Defines associated types
+pub trait ArgminTransition {
+    /// Associated type for 1-D array
+    type Array1D;
+    /// Associated type for 2-D array
+    type Array2D;
+}
+
+/// Create a random array
+pub trait ArgminRandomMatrix {
+    /// Draw samples from a standard Normal distribution (mean=0, stdev=1)
+    fn standard_normal(nrows: usize, ncols: usize) -> Self;
+}
+
+/// Get dimensions of an array
+pub trait ArgminSize<U> {
+    /// Get array shape
+    fn shape(&self) -> U;
+}
+
+/// Eigendecomposition for symmetric matrices
+pub trait ArgminEigenSym<U> {
+    /// Get the eigenvalues and eigenvectors of a complex Hermitian (conjugate symmetric) or a real symmetric matrix.
+    /// Returns a 1-D array containing the eigenvalues and a 2-D square array with eigenvectors (in columns).
+    fn eig_sym(&self) -> (U, Self);
+}
+
+/// Eigendecomposition
+pub trait ArgminEigen<U> {
+    /// Get the eigenvalues and eigenvectors of any matrix.
+    fn eig(&self) -> (U, U, Self);
+}
+
+/// Returns the indices that would sort an array.
+pub trait ArgminArgsort {
+    /// Perform sort and return an array of indices of the same shape.
+    fn argsort(&self) -> Vec<usize>;
+}
+
+/// Take elements from an array along an axis.
+pub trait ArgminTake<I> {
+    /// Take elements from indices along a given axis.
+    /// Returns as new array of the same size as indices
+    fn take(&self, indices: &[I], axis: u8) -> Self;
+}
+
+/// Defines iterator
+pub trait ArgminIter<T> {
+    /// Get immutable, read-only iterator
+    fn iterator<'a>(&'a self) -> Box<dyn Iterator<Item = &'a T> + 'a>;
+}
+
+/// Defines mutable iterator
+pub trait ArgminMutIter<T> {
+    /// Get mutable iterator
+    fn iterator_mut<'a>(&'a mut self) -> Box<dyn Iterator<Item = &'a mut T> + 'a>;
+}
+
+/// Defines axis iterator
+pub trait ArgminAxisIter<T> {
+    /// Get immutable, read-only row iterator
+    fn row_iterator<'a>(&'a self) -> Box<dyn Iterator<Item = T> + 'a>;
+    /// Get immutable, read-only column iterator
+    fn column_iterator<'a>(&'a self) -> Box<dyn Iterator<Item = T> + 'a>;
+}
+
+/// Defines methods that can be used to access an element at a given location
+pub trait ArgminGet<S, T> {
+    /// Get element at position
+    fn get(&self, pos: S) -> &T;
+}
+
+/// Defines methods that can be used to set an element at a given location
+pub trait ArgminSet<S, T> {
+    /// Set element at position to a new value
+    fn set(&mut self, pos: S, x: T);
+}
+
+/// Create new instance of an array
+pub trait ArgminFrom<T, S> {
+    /// Create new array from an iterator
+    fn from_iterator<I: Iterator<Item = T>>(size: S, iter: I) -> Self;
+}
+
+/// Compute the outer product of two vectors.
+pub trait ArgminOuter<T, U> {
+    /// Get an outer product of self with array other
+    fn outer(&self, other: &T) -> U;
+}
+
+/// Broadcast operation.
+pub trait ArgminBroadcast<T, U> {
+    /// Broadcast add `T` to `self`
+    fn broadcast_add(&self, other: &T) -> U;
+
+    /// Broadcast subtract `T` from `self`
+    fn broadcast_sub(&self, other: &T) -> U;
+
+    /// Broadcast multiply a `T` by `self`
+    fn broadcast_mul(&self, other: &T) -> U;
+
+    /// Broadcast divide a `T` by `self`
+    fn broadcast_div(&self, other: &T) -> U;
+}

argmin-math/src/nalgebra_m/argsort.rs

@@ -0,0 +1,47 @@
+// Copyright 2018-2022 argmin developers
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// http://apache.org/licenses/LICENSE-2.0> or the MIT license <LICENSE-MIT or
+// http://opensource.org/licenses/MIT>, at your option. This file may not be
+// copied, modified, or distributed except according to those terms.
+
+use crate::ArgminArgsort;
+use nalgebra::{
+    base::{dimension::Dim, storage::Storage},
+    Matrix,
+};
+
+impl<N, R, C, S> ArgminArgsort for Matrix<N, R, C, S>
+where
+    N: PartialOrd,
+    R: Dim,
+    C: Dim,
+    S: Storage<N, R, C>,
+{
+    #[inline]
+    fn argsort(&self) -> Vec<usize> {
+        let (n, m) = self.shape();
+        let l = if n > m { n } else { m };
+        let mut indices = (0..l).collect::<Vec<_>>();
+        indices.sort_by(|&i, &j| self[i].partial_cmp(&self[j]).unwrap());
+        indices
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use nalgebra::matrix;
+
+    #[test]
+    fn test_argsort() {
+        assert_eq!(
+            ArgminArgsort::argsort(&matrix![2, 4, 5, 2, 5, 0, 1]),
+            vec![5, 6, 0, 3, 1, 2, 4]
+        );
+        assert_eq!(
+            ArgminArgsort::argsort(&matrix![2., 4., 5., 2., 5., 0., 1.]),
+            vec![5, 6, 0, 3, 1, 2, 4]
+        );
+    }
+}

argmin-math/src/nalgebra_m/broadcast.rs

@@ -0,0 +1,202 @@
+// Copyright 2018-2022 argmin developers
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// http://apache.org/licenses/LICENSE-2.0> or the MIT license <LICENSE-MIT or
+// http://opensource.org/licenses/MIT>, at your option. This file may not be
+// copied, modified, or distributed except according to those terms.
+
+use crate::ArgminBroadcast;
+
+use nalgebra::{
+    base::{allocator::Allocator, dimension::Dim, storage::Storage, Scalar},
+    ClosedAdd, ClosedDiv, ClosedMul, ClosedSub, DefaultAllocator, Matrix, OMatrix,
+};
+
+impl<N, R1, C1, R2, C2, SV, SM> ArgminBroadcast<Matrix<N, R2, C2, SV>, OMatrix<N, R1, C1>>
+    for Matrix<N, R1, C1, SM>
+where
+    N: Copy + Scalar + ClosedAdd + ClosedMul + ClosedSub + ClosedDiv,
+    R1: Dim,
+    C1: Dim,
+    R2: Dim,
+    C2: Dim,
+    SM: Storage<N, R1, C1>,
+    SV: Storage<N, R2, C2>,
+    DefaultAllocator: Allocator<N, R1, C1>,
+{
+    #[inline]
+    fn broadcast_add(&self, other: &Matrix<N, R2, C2, SV>) -> OMatrix<N, R1, C1> {
+        let (n, m) = self.shape();
+        match other.shape() {
+            (1, _) | (_, 1) => OMatrix::<N, R1, C1>::from_iterator_generic(
+                R1::from_usize(n),
+                C1::from_usize(m),
+                (0..m).flat_map(move |i| (0..n).map(move |j| self[(j, i)] + other[i])),
+            ),
+            _ => panic!(
+                "Can't broadcast matrix {}x{} to matrix {}x{}, yet",
+                other.nrows(),
+                other.ncols(),
+                self.nrows(),
+                self.ncols()
+            ),
+        }
+    }
+
+    #[inline]
+    fn broadcast_sub(&self, other: &Matrix<N, R2, C2, SV>) -> OMatrix<N, R1, C1> {
+        let (n, m) = self.shape();
+        match other.shape() {
+            (1, _) | (_, 1) => OMatrix::<N, R1, C1>::from_iterator_generic(
+                R1::from_usize(n),
+                C1::from_usize(m),
+                (0..m).flat_map(move |i| (0..n).map(move |j| self[(j, i)] - other[i])),
+            ),
+            _ => panic!(
+                "Can't broadcast matrix {}x{} to matrix {}x{}, yet",
+                other.nrows(),
+                other.ncols(),
+                self.nrows(),
+                self.ncols()
+            ),
+        }
+    }
+
+    #[inline]
+    fn broadcast_mul(&self, other: &Matrix<N, R2, C2, SV>) -> OMatrix<N, R1, C1> {
+        let (n, m) = self.shape();
+        match other.shape() {
+            (1, _) | (_, 1) => OMatrix::<N, R1, C1>::from_iterator_generic(
+                R1::from_usize(n),
+                C1::from_usize(m),
+                (0..m).flat_map(move |i| (0..n).map(move |j| self[(j, i)] * other[i])),
+            ),
+            _ => panic!(
+                "Can't broadcast matrix {}x{} to matrix {}x{}, yet",
+                other.nrows(),
+                other.ncols(),
+                self.nrows(),
+                self.ncols()
+            ),
+        }
+    }
+
+    #[inline]
+    fn broadcast_div(&self, other: &Matrix<N, R2, C2, SV>) -> OMatrix<N, R1, C1> {
+        let (n, m) = self.shape();
+        match other.shape() {
+            (1, _) | (_, 1) => OMatrix::<N, R1, C1>::from_iterator_generic(
+                R1::from_usize(n),
+                C1::from_usize(m),
+                (0..m).flat_map(move |i| (0..n).map(move |j| self[(j, i)] / other[i])),
+            ),
+            _ => panic!(
+                "Can't broadcast matrix {}x{} to matrix {}x{}, yet",
+                other.nrows(),
+                other.ncols(),
+                self.nrows(),
+                self.ncols()
+            ),
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use nalgebra::{Matrix2x3, RowVector3, Vector3};
+    use paste::item;
+
+    macro_rules! make_test {
+        ($t:ty) => {
+            item! {
+                #[test]
+                fn [<test_broadcast_add_sub_mul_div_row_vec_ $t>]() {
+                    let a = Matrix2x3::new(
+                        2 as $t, 4 as $t, 6 as $t,
+                        1 as $t, 2 as $t, 3 as $t
+                    );
+                    let b = RowVector3::new(1 as $t, 2 as $t, 3 as $t);
+                    let expected_add = Matrix2x3::new(
+                        3 as $t, 6 as $t, 9 as $t,
+                        2 as $t, 4 as $t, 6 as $t
+                    );
+                    let expected_sub = Matrix2x3::new(
+                        1 as $t, 2 as $t, 3 as $t,
+                        0 as $t, 0 as $t, 0 as $t
+                    );
+                    let expected_mul = Matrix2x3::new(
+                        2 as $t, 8 as $t, 18 as $t,
+                        1 as $t, 4 as $t, 9 as $t
+                    );
+                    let expected_div = Matrix2x3::new(
+                        2 as $t, 2 as $t, 2 as $t,
+                        1 as $t, 1 as $t, 1 as $t
+                    );
+                    let res_add = a.broadcast_add(&b);
+                    let res_sub = a.broadcast_sub(&b);
+                    let res_mul = a.broadcast_mul(&b);
+                    let res_div = a.broadcast_div(&b);
+                    for i in 0..2 {
+                        for j in 0..3 {
+                            assert!((((res_add[(i, j)] - expected_add[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                            assert!((((res_sub[(i, j)] - expected_sub[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                            assert!((((res_mul[(i, j)] - expected_mul[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                            assert!((((res_div[(i, j)] - expected_div[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                        }
+                    }
+                }
+            }
+
+            item! {
+                #[test]
+                fn [<test_broadcast_add_sub_mul_div_vec_ $t>]() {
+                    let a = Matrix2x3::new(
+                        2 as $t, 4 as $t, 6 as $t,
+                        1 as $t, 2 as $t, 3 as $t
+                    );
+                    let b = Vector3::new(1 as $t, 2 as $t, 3 as $t);
+                    let expected_add = Matrix2x3::new(
+                        3 as $t, 6 as $t, 9 as $t,
+                        2 as $t, 4 as $t, 6 as $t
+                    );
+                    let expected_sub = Matrix2x3::new(
+                        1 as $t, 2 as $t, 3 as $t,
+                        0 as $t, 0 as $t, 0 as $t
+                    );
+                    let expected_mul = Matrix2x3::new(
+                        2 as $t, 8 as $t, 18 as $t,
+                        1 as $t, 4 as $t, 9 as $t
+                    );
+                    let expected_div = Matrix2x3::new(
+                        2 as $t, 2 as $t, 2 as $t,
+                        1 as $t, 1 as $t, 1 as $t
+                    );
+                    let res_add = a.broadcast_add(&b);
+                    let res_sub = a.broadcast_sub(&b);
+                    let res_mul = a.broadcast_mul(&b);
+                    let res_div = a.broadcast_div(&b);
+                    for i in 0..2 {
+                        for j in 0..3 {
+                            assert!((((res_add[(i, j)] - expected_add[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                            assert!((((res_sub[(i, j)] - expected_sub[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                            assert!((((res_mul[(i, j)] - expected_mul[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                            assert!((((res_div[(i, j)] - expected_div[(i, j)]) as f64).abs()) < std::f64::EPSILON);
+                        }
+                    }
+                }
+            }
+        };
+    }
+
+    make_test!(i8);
+    make_test!(u8);
+    make_test!(i16);
+    make_test!(u16);
+    make_test!(i32);
+    make_test!(u32);
+    make_test!(i64);
+    make_test!(u64);
+    make_test!(f32);
+    make_test!(f64);
+}