Merge pull request #1 from wgurecky/faer_0190

Update to faer 0.19.0

Cargo.toml

@@ -16,20 +16,19 @@ assert_approx_eq = "*"
 itertools = "*"
 rand = {version=">=0.8.5", featurs=["alloc"]}
 statrs = "*"
-argmin = {version=">=0.9.0"}
+argmin = {version="0.10.0"}
 argmin-math = "*"
 finitediff = {version="0.1.4", features = ["ndarray"]}
 rand_distr = "*"
 num-traits = "*"
 num-complex = "*"
 ndarray = {version=">=0.15.6", features=["rayon"]}
-# nalgebra = "*"
-faer = {version = "0.16.0", features = ["ndarray"]}
-faer-core = "0.16.0"
-faer-svd = "0.16.0"
+faer = {version = "0.19.0"}
+faer-ext = {version="*", features=["ndarray"]}
 rayon = "*"
-pyo3 = {version = "0.20", features = ["extension-module"]}
-numpy = "0.20"
-# kiddo = "3.0.0"
+pyo3 = {version = "0.21.2", features = ["extension-module"]}
+numpy = "0.21"
 kdtree = "0.7.0"
+# nalgebra = "*"
+# kiddo = "3.0.0"
 # cudarc = { version = "*", default-features = false, optional = true, features = ["driver", "cublas", "nvrtc"] }

examples/benchmark_dmd.py

@@ -34,7 +34,7 @@ def u_fn(t):
     p_snapshots = np.asarray(p_snapshots).T
     plot_snapshots(x_points, t_points, p_snapshots)
 
-    return p_snapshots, u_seq
+    return x_points, p_snapshots, u_seq
 
 
 def plot_snapshots(x_points, t_points, p_snapshots):
@@ -77,7 +77,7 @@ def predict_dmdc(pydmdc_model: DMDc, x0, u_seq):
     return np.hstack(predicted_x)
 
 def fit_pydmd():
-    p_snapshots, u_seq = build_snapshots()
+    x_points, p_snapshots, u_seq = build_snapshots()
     n_modes = 12
     ti = time.time()
     pydmdc_model = DMDc(svd_rank=n_modes, svd_rank_omega=n_modes)
@@ -95,19 +95,26 @@ def fit_pydmd():
     #print("pydmdc modes: ", pydmdc_model.modes)
 
     # forecast
-    predicted = predict_dmdc(pydmdc_model,
-                             p_snapshots[:, 0:1],
-                             u_seq[:, 1:]
-                             )
+    x0 = p_snapshots[:, 0:1]
+    u_s = u_seq[:, 1:]
+    predicted = predict_dmdc(pydmdc_model, x0, u_s)
 
     print("pydmdc predicted: ", predicted[:, 20])
     print("pydmdc expected: ", p_snapshots[:, 20])
 
+    plt.figure()
+    for i in range(0, predicted.shape[1]-1):
+        plt.plot(x_points, predicted[:, i], label="step: %d, u: %0.2f" % (i, u_s[0, i]))
+    plt.legend()
+    plt.title("PyDMD DMDc prediction")
+    plt.savefig("pydmdc_pydmd_predict.png")
+    plt.close()
+
 
 def fit_corrla_dmd():
-    p_snapshots, u_seq = build_snapshots()
+    x_points, p_snapshots, u_seq = build_snapshots()
     n_modes = 12
-    n_iters = 10
+    n_iters = 20
     ti = time.time()
     rust_dmdc = PyDMDc(p_snapshots, u_seq, n_modes, n_iters)
     tf = time.time()
@@ -119,6 +126,14 @@ def fit_corrla_dmd():
     print("corrla dmdc predicted: ", predicted[:, 20])
     print("corrla dmdc expected: ", p_snapshots[:, 20])
 
+    plt.figure()
+    for i in range(0, predicted.shape[1]):
+        plt.plot(x_points, predicted[:, i], label="step: %d, u: %0.2f" % (i, u_s[0, i]))
+    plt.legend()
+    plt.title("CORRLA-RS DMDc prediction")
+    plt.savefig("pydmdc_corrla_predict.png")
+    plt.close()
+
 
 if __name__ == "__main__":
     fit_pydmd()

readme.md

@@ -1,6 +1,8 @@
 About
 ================
 
+![status](https://github.com/wgurecky/CORRLA_RS/actions/workflows/rust.yml/badge.svg)
+
 CORRLA-RS is a package for correlation, regression and random linear algebra methods in Rust.
 
 Methods

src/lib.rs

@@ -1,45 +1,2 @@
-// root level lib
-use faer::{prelude::*};
-use faer_core::{mat, Mat, MatRef};
-
-// other modules
 pub mod lib_math_utils;
 pub mod lib_math_utils_py;
-
-// example root level routines in the lib
-pub fn my_add(left: usize, right: usize) -> usize {
-    left + right
-}
-
-
-pub fn svd_random_faer<'a>(my_mat: Mat<f64>) -> (Mat<f64>, Mat<f64>, Mat<f64>) {
-    print!("test");
-    print!("test2");
-    let my_svd = my_mat.svd();
-    (my_svd.u().to_owned(), my_svd.s_diagonal().to_owned(), my_svd.v().to_owned())
-}
-
-
-// example simple within-lib test
-#[cfg(test)]
-mod lib_unit_tests {
-    use super::*;
-
-    #[test]
-    fn it_works() {
-        let result = my_add(2, 2);
-        assert_eq!(result, 4);
-    }
-
-    #[test]
-    fn tst_svd(){
-        let a = mat![
-            [1., 2.],
-            [3., 4.],
-            ];
-        // multiple assignment with types
-        let (b, c, d): (Mat<f64>,Mat<f64>,Mat<f64>) = svd_random_faer(a);
-        print!("svd test \n");
-        print!("{:?} \n", b);
-    }
-}

src/lib_math_utils/active_subspaces.rs

@@ -6,22 +6,17 @@
 /// practice. https://arxiv.org/pdf/1304.2070.pdf
 ///
 use std::cmp;
-use assert_approx_eq::assert_approx_eq;
-use faer::{prelude::*, IntoFaer};
-use faer_core::{mat, Mat, MatRef, MatMut, Entity, AsMatRef, AsMatMut};
-use faer_core::{ComplexField, RealField, c32, c64};
-use faer::solvers::{Eigendecomposition};
-// use kiddo::{KdTree, SquaredEuclidean};
+use faer::prelude::*;
+use faer_ext::*;
+use faer::solvers::Eigendecomposition;
 use kdtree::KdTree;
-use kdtree::ErrorKind;
 use kdtree::distance::squared_euclidean;
 use ndarray::prelude::*;
 use rayon::prelude::*;
 
 use crate::lib_math_utils::mat_utils::*;
 use crate::lib_math_utils::random_svd::*;
 use crate::lib_math_utils::stats_corr::*;
-use crate::lib_math_utils::pca_rsvd::{ApplyTransform};
 
 /// Stores owned data for gradient estimation
 /// over point cloud datasets
@@ -263,11 +258,11 @@ impl ActiveSsRsvd {
         assert!(grad_mat_sc.ncols() == x_mat.ncols());
         // all eigenvalues of a real sym mat are real
         let evd: Eigendecomposition<c64> = grad_mat_sc.eigendecomposition();
-        let eigs = evd.s_diagonal();
+        let eigs = evd.s();
         let eig_vs = evd.u();
         // split real from complex part, discard imag parts (zero)
         let (real_eigs, _imag_eigs) =
-            mat_parts_from_complex(mat_colmat_to_diag(eigs).as_ref());
+            mat_parts_from_complex(mat_diagref_to_2d(eigs).as_ref());
         let (real_eig_vs, _imag_eig_vs) =
             mat_parts_from_complex(eig_vs);
         // sort eigenvalues and eigenvectors from largest to smallest
@@ -287,6 +282,7 @@ impl ActiveSsRsvd {
 mod active_subspace_unit_tests {
     // bring everything from above (parent) module into scope
     use super::*;
+    use assert_approx_eq::assert_approx_eq;
 
     #[test]
     fn test_grad_est() {

src/lib_math_utils/dmd_rom.rs

@@ -9,16 +9,13 @@
 /// Where u_t is the control or forcing term
 /// and x_t is the state vector.
 use faer::{prelude::*};
-use faer_core::{mat, Mat, MatRef, MatMut, Entity, AsMatRef, AsMatMut, ColMut};
-use faer_core::{ComplexField, RealField, c32, c64};
-use faer::solvers::{Eigendecomposition};
-use num_traits::Float;
+use faer::{mat, Mat, MatRef, MatMut};
+use faer::solvers::Eigendecomposition;
 use std::marker;
 
 // internal imports
 use crate::lib_math_utils::mat_utils::*;
 use crate::lib_math_utils::random_svd::*;
-use crate::lib_math_utils::interp_utils::*;
 
 pub struct DMDc <'a> {
     // number of supplied data snapshots
@@ -72,7 +69,7 @@ impl <'a> DMDc <'a> {
         // let (u_til, s_til, v_til_) = mat_truncated_svd(self._X(omega.as_ref()), self.n_modes);
         // let v_til = v_til_.to_owned();
         // println!("u_til_exact: {:?}, v_til_exact: {:?}", u_til, v_til_);
-        let (u_til, s_til, v_til_) = random_svd(self._X(omega.as_ref()), self.n_modes, n_iters, 8);
+        let (u_til, s_til, v_til_) = random_svd(self._X(omega.as_ref()), self.n_modes, n_iters, 12);
         let v_til = v_til_.transpose().to_owned();
 
         let u_til_1 = u_til.as_ref().submatrix(
@@ -117,7 +114,7 @@ impl <'a> DMDc <'a> {
         //    self._A.as_ref().unwrap().as_ref());
         let ev: Eigendecomposition<c64> = (self._A.as_ref()).unwrap().eigendecomposition();
         let a_til_eigenvectors = ev.u();
-        let a_til_eigenvalues = ev.s_diagonal().as_2d();
+        let a_til_eigenvalues = mat_diagref_to_2d(ev.s());
         // convert to real and imag components
         let (a_til_eigenvectors_re, a_til_eigenvectors_im) =
             mat_parts_from_complex(a_til_eigenvectors);
@@ -130,8 +127,7 @@ impl <'a> DMDc <'a> {
     /// Computes DMD modes
     fn _calc_modes(&mut self, omega: MatRef<f64>, v_til: MatRef<f64>, s_til: MatRef<f64>, u_til_1: MatRef<f64>, u_hat: MatRef<f64>) {
         let (lambdas, w_re, w_im) = self._calc_eigs();
-        let lambdas_diag: Mat<c64> = mat_colvec_to_diag(lambdas.as_ref());
-        self.lambdas = Some(lambdas_diag);
+        self.lambdas = Some(lambdas);
         // from eq 36 in Proctor. et. al DMDc
         // BUT we only need the real part of the modes, since
         // when we recombine with
@@ -275,8 +271,8 @@ mod dmd_unit_tests {
         println!("u_data shape: {:?}, {:?}", u_mat.nrows(), u_mat.ncols());
 
         // build DMDc model
-        let n_modes = 8;
-        let dmdc_model = DMDc::new(p_snapshots.as_ref(), u_mat.as_ref(), 1.0, n_modes, 10);
+        let n_modes = 12;
+        let dmdc_model = DMDc::new(p_snapshots.as_ref(), u_mat.as_ref(), 1.0, n_modes, 30);
 
         // test the DMDc model
         let estimated_a_op = dmdc_model.est_a_til();
@@ -309,6 +305,6 @@ mod dmd_unit_tests {
         println!("Predicted: {:?}", p20_predicted);
         println!("DMDc Eigs: {:?}", dmdc_model.lambdas.as_ref());
         assert_eq!(dmdc_model.lambdas.unwrap().nrows(), n_modes);
-        mat_mat_approx_eq(p20_predicted.as_ref(), p20_expected.as_ref(), 1e-3);
+        mat_mat_approx_eq(p20_predicted.as_ref(), p20_expected.as_ref(), 5e-3);
     }
 }

src/lib_math_utils/interp_utils.rs

@@ -1,6 +1,6 @@
 /// Interpolation utility methods
-use faer::{prelude::*, IntoNdarray};
-use faer_core::{mat, Mat, MatRef, MatMut, Entity, AsMatRef, AsMatMut};
+use faer::prelude::*;
+use faer::{mat, Mat, MatRef, MatMut};
 
 // internal imports
 use crate::lib_math_utils::mat_utils::*;