diff --git a/Cargo.toml b/Cargo.toml
index 13d4d15..7addf66 100644
--- a/Cargo.toml
+++ b/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"] }
diff --git a/examples/benchmark_dmd.py b/examples/benchmark_dmd.py
index 6c40066..8c5f488 100644
--- a/examples/benchmark_dmd.py
+++ b/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()
diff --git a/readme.md b/readme.md
index 80c4ef1..946532c 100644
--- a/readme.md
+++ b/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
diff --git a/src/bin/hello_main.rs b/src/bin/hello_main.rs
deleted file mode 100644
index 2b58847..0000000
--- a/src/bin/hello_main.rs
+++ /dev/null
@@ -1,7 +0,0 @@
-use corrla_rs;
-
-fn main(){
-    let result_a = corrla_rs::my_add(4,8);
-    let result_b = corrla_rs::my_add(8,8);
-    println!("hello_world_from_bin: {}, {}", result_a, result_b);
-}
diff --git a/src/lib.rs b/src/lib.rs
index 06aa02d..0561581 100644
--- a/src/lib.rs
+++ b/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);
-    }
-}
diff --git a/src/lib_math_utils/active_subspaces.rs b/src/lib_math_utils/active_subspaces.rs
index e2511e4..3eb339d 100644
--- a/src/lib_math_utils/active_subspaces.rs
+++ b/src/lib_math_utils/active_subspaces.rs
@@ -6,14 +6,10 @@
 /// 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::*;
@@ -21,7 +17,6 @@ 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() {
diff --git a/src/lib_math_utils/dmd_rom.rs b/src/lib_math_utils/dmd_rom.rs
index 37a23e6..39d778a 100644
--- a/src/lib_math_utils/dmd_rom.rs
+++ b/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);
     }
 }
diff --git a/src/lib_math_utils/interp_utils.rs b/src/lib_math_utils/interp_utils.rs
index 8bf4fd2..b5dee7d 100644
--- a/src/lib_math_utils/interp_utils.rs
+++ b/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::*;
diff --git a/src/lib_math_utils/mat_utils.rs b/src/lib_math_utils/mat_utils.rs
index 416969c..a3553b1 100644
--- a/src/lib_math_utils/mat_utils.rs
+++ b/src/lib_math_utils/mat_utils.rs
@@ -4,11 +4,12 @@
 /// library, but are generally missing from rust libs
 ///
 use num_traits::Float;
-use faer::{prelude::*, IntoNdarray};
-use faer_core::{mat, Mat, MatRef, MatMut, Entity, AsMatRef, AsMatMut};
-use faer_core::ColRef;
-use faer_core::{c64, c32};
-use rand::{prelude::*};
+use faer::prelude::*;
+use faer::linalg::matmul;
+use faer::{mat, Mat, MatRef, ColRef, MatMut};
+use faer::diag::{DiagRef};
+use faer_ext::*;
+use rand::prelude::*;
 use rand_distr::{StandardNormal, Uniform};
 use rayon::prelude::*;
 use ndarray::prelude::*;
@@ -18,16 +19,16 @@ use ndarray::prelude::*;
 /// parallel control exposed
 pub fn par_matmul_helper<T>(res: MatMut<T>, lhs: MatRef<T>, rhs: MatRef<T>, beta: T, n_threads: usize)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
-    faer_core::mul::matmul(
+    matmul::matmul(
         res,
         lhs,
         rhs,
         None,
         beta,
-        // faer_core::Parallelism::Rayon(n_threads),
-        faer_core::get_global_parallelism()
+        // faer::Parallelism::Rayon(n_threads),
+        faer::get_global_parallelism()
         );
 }
 
@@ -35,9 +36,10 @@ pub fn par_matmul_helper<T>(res: MatMut<T>, lhs: MatRef<T>, rhs: MatRef<T>, beta
 /// Compute the Moore-Penrose inverse
 pub fn mat_pinv<T>(x: MatRef<T>) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let x_svd = x.svd();
+    // x_svd.pseudoinverse()
     let u = x_svd.u();
     let eps = T::from(1.0e-14).unwrap();
     let s_vec = x_svd.s_diagonal();
@@ -45,7 +47,7 @@ pub fn mat_pinv<T>(x: MatRef<T>) -> Mat<T>
     let mut s_inv_mat = faer::Mat::zeros(x.ncols(), x.nrows());
     for i in 0..s_vec.nrows(){
         s_inv_mat.write(i, i, T::from(1.).unwrap() /
-                        (s_vec.read(i, 0) + eps));
+                        (s_vec.read(i) + eps));
     }
     v * s_inv_mat * u.transpose()
 }
@@ -54,6 +56,7 @@ pub fn mat_pinv<T>(x: MatRef<T>) -> Mat<T>
 pub fn mat_pinv_comp(x: MatRef<c64>) -> Mat<c64>
 {
     let x_svd = x.svd();
+    // x_svd.pseudoinverse()
     let u = x_svd.u();
     let eps = c64::new(1.0e-16, 1.0e-16);
     let comp_one = c64::new(1.0, 0.0);
@@ -62,7 +65,7 @@ pub fn mat_pinv_comp(x: MatRef<c64>) -> Mat<c64>
     let mut s_inv_mat: Mat<c64> = faer::Mat::zeros(x.ncols(), x.nrows());
     for i in 0..s_vec.nrows(){
         s_inv_mat.write(i, i, comp_one /
-                        (s_vec.read(i, 0) + eps));
+                        (s_vec.read(i) + eps));
     }
     v * s_inv_mat * u.adjoint()
 }
@@ -70,11 +73,11 @@ pub fn mat_pinv_comp(x: MatRef<c64>) -> Mat<c64>
 /// Truncated SVD
 pub fn mat_truncated_svd<T>(my_mat: MatRef<T>, rank: usize) -> (Mat<T>, Mat<T>, Mat<T>)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let my_svd = my_mat.svd();
     let u_r = my_svd.u().get(.., 0..rank).to_owned();
-    let s_r = my_svd.s_diagonal().get(0..rank, ..).to_owned();
+    let s_r = my_svd.s_diagonal().get(0..rank).as_2d().to_owned();
     let v_r = my_svd.v().get(.., 0..rank).to_owned();
     (u_r, s_r, v_r)
 }
@@ -83,7 +86,7 @@ pub fn mat_truncated_svd<T>(my_mat: MatRef<T>, rank: usize) -> (Mat<T>, Mat<T>,
 /// Computes mean along axis
 pub fn mat_mean<T>(a_mat: MatRef<T>, axis: usize) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut acc: T = T::from(0.0).unwrap();
     let a_ncols: usize = a_mat.ncols();
@@ -118,7 +121,7 @@ pub fn mat_mean<T>(a_mat: MatRef<T>, axis: usize) -> Mat<T>
 /// Computes std dev along axis
 pub fn mat_std<T>(a_mat: MatRef<T>, axis: usize) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut acc: T = T::from(0.0).unwrap();
     let a_ncols: usize = a_mat.ncols();
@@ -158,7 +161,7 @@ pub fn mat_std<T>(a_mat: MatRef<T>, axis: usize) -> Mat<T>
 pub fn random_mat_normal<T>(n_rows: usize, n_cols: usize)
     -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let omega: Mat<T> = Mat::from_fn(
         n_rows,
@@ -175,7 +178,7 @@ pub fn random_mat_normal<T>(n_rows: usize, n_cols: usize)
 pub fn random_mat_uniform<T>(n_rows: usize, n_cols: usize, lb: f64, ub: f64)
     -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let uni_dist = Uniform::new(lb, ub);
     let omega: Mat<T> = Mat::from_fn(
@@ -256,7 +259,7 @@ pub fn mat_vec_add(x: MatRef<f64>, pv: MatRef<f64>, axis: usize)
 /// Raise all elements in a matrix to a power
 pub fn mat_ele_pow<T>(a_mat: MatRef<T>, pwr: T) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut out_mat = a_mat.to_owned();
     for i in 0..out_mat.nrows() {
@@ -270,7 +273,7 @@ pub fn mat_ele_pow<T>(a_mat: MatRef<T>, pwr: T) -> Mat<T>
 /// Element by element matrix multiplication
 pub fn mat_mat_ele_mul<T>(a: MatRef<T>, b: MatRef<T>) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut out_mat = a.to_owned();
     for i in 0..out_mat.nrows() {
@@ -284,7 +287,7 @@ pub fn mat_mat_ele_mul<T>(a: MatRef<T>, b: MatRef<T>) -> Mat<T>
 /// Peforms matrix scalar addition in-place
 pub fn mat_scalar_add<T>(mut out_mat: MatMut<T>, b: T)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     // consume a, modify in-place
     for i in 0..out_mat.nrows() {
@@ -297,7 +300,7 @@ pub fn mat_scalar_add<T>(mut out_mat: MatMut<T>, b: T)
 /// Modifies matrix row in-place
 pub fn mat_row_mod<T>(mut out_mat: MatMut<T>, row: usize, vec: MatRef<T>)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     for j in 0..vec.ncols() {
         out_mat.write(row, j, vec.read(0, j));
@@ -307,7 +310,7 @@ pub fn mat_row_mod<T>(mut out_mat: MatMut<T>, row: usize, vec: MatRef<T>)
 /// Modifies matrix col in-place
 pub fn mat_col_mod<T>(mut out_mat: MatMut<T>, col: usize, vec: MatRef<T>)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     for i in 0..vec.nrows() {
         out_mat.write(i, col, vec.read(i, 0));
@@ -342,7 +345,7 @@ pub fn mat_parts_from_complex(a_comp: MatRef<c64>) -> (Mat<f64>, Mat<f64>) {
 /// converts a col vector to a diagnoal matrix
 pub fn mat_colvec_to_diag<T>(vec: MatRef<T>) -> Mat<T>
     where
-    T: faer_core::ComplexField
+    T: faer::ComplexField
 {
     let mut out_mat = faer::Mat::zeros(vec.nrows(), vec.nrows());
     for i in 0..vec.nrows() {
@@ -354,7 +357,7 @@ pub fn mat_colvec_to_diag<T>(vec: MatRef<T>) -> Mat<T>
 /// converts a col mat to a diagnoal matrix
 pub fn mat_colmat_to_diag<T>(vec: ColRef<T>) -> Mat<T>
     where
-    T: faer_core::ComplexField
+    T: faer::ComplexField
 {
     let mut out_mat = faer::Mat::zeros(vec.nrows(), vec.nrows());
     for i in 0..vec.nrows() {
@@ -367,7 +370,7 @@ pub fn mat_colmat_to_diag<T>(vec: ColRef<T>) -> Mat<T>
 /// converts a row vector to a diagnoal matrix
 pub fn mat_rowvec_to_diag<T>(vec: MatRef<T>) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut out_mat = faer::Mat::zeros(vec.ncols(), vec.ncols());
     for i in 0..vec.ncols() {
@@ -376,10 +379,18 @@ pub fn mat_rowvec_to_diag<T>(vec: MatRef<T>) -> Mat<T>
     out_mat
 }
 
+/// converts a DiagRef to a diagnoal in 2d matrix
+pub fn mat_diagref_to_2d<T>(diag: DiagRef<T>) -> Mat<T>
+    where
+    T: faer::ComplexField
+{
+    mat_colmat_to_diag(diag.column_vector()).to_owned()
+}
+
 /// pseudo inv of diag matrix
 pub fn mat_pinv_diag<T>(diag_mat: MatRef<T>) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let eps = T::from(1.0e-20).unwrap();
     let mut out_mat = faer::Mat::zeros(diag_mat.nrows(), diag_mat.ncols());
@@ -398,7 +409,7 @@ pub fn mat_pinv_diag<T>(diag_mat: MatRef<T>) -> Mat<T>
 /// create a owned Vec from row, results in data copy
 pub fn mat_row_to_vec<T>(a_mat: MatRef<T>, row: usize) -> Vec<T>
     where
-    T: faer_core::SimpleEntity
+    T: faer::SimpleEntity
 {
     let tmp_ndarray = a_mat.get(row..row+1, ..).into_ndarray();
     let tmp_1darray: Array1<T> = tmp_ndarray.remove_axis(Axis(0)).to_owned();
@@ -408,7 +419,7 @@ pub fn mat_row_to_vec<T>(a_mat: MatRef<T>, row: usize) -> Vec<T>
 /// create a owned Vec from col, results in data copy
 pub fn mat_col_to_vec<T>(a_mat: MatRef<T>, col: usize) -> Vec<T>
     where
-    T: faer_core::SimpleEntity
+    T: faer::SimpleEntity
 {
     let tmp_ndarray = a_mat.get(.., col..col+1).into_ndarray();
     let tmp_1darray: Array1<T> = tmp_ndarray.remove_axis(Axis(1)).to_owned();
@@ -418,7 +429,7 @@ pub fn mat_col_to_vec<T>(a_mat: MatRef<T>, col: usize) -> Vec<T>
 /// create a owned Vec from mat diag, results in data copy
 pub fn mat_diag_to_vec<T>(a_mat: MatRef<T>) -> Vec<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut diag_vec: Vec<T> = Vec::new();
     for i in 0..a_mat.ncols() {
@@ -431,7 +442,7 @@ pub fn mat_diag_to_vec<T>(a_mat: MatRef<T>) -> Vec<T>
 /// TODO: make this work for real and imag vecs
 pub fn argsort_float<T>(a: &[T]) -> Vec<usize>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut idx = (0..a.len()).collect::<Vec<_>>();
     idx.sort_by(|&i, &j| a[i].partial_cmp(&a[j]).unwrap());
@@ -440,7 +451,7 @@ pub fn argsort_float<T>(a: &[T]) -> Vec<usize>
 /// Returns indicies that would sort a vec of floats, high to lowest
 pub fn argsort_float_rev<T>(a: &[T]) -> Vec<usize>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut idx = (0..a.len()).collect::<Vec<_>>();
     idx.sort_by(|&i, &j| a[j].partial_cmp(&a[i]).unwrap());
@@ -452,7 +463,7 @@ pub fn argsort_float_rev<T>(a: &[T]) -> Vec<usize>
 /// Only works for real matricies. TODO: fix for both real and complex mats
 pub fn sort_evd<T>(eigs: MatRef<T>, eigvs: MatRef<T>) -> (Mat<T>, Mat<T>)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let dim = eigs.ncols();
     // compute ordering
@@ -529,7 +540,7 @@ pub fn zcenter_mat_col(a_mat: MatRef<f64>) -> Mat<f64>
 // Helper function to ensure two matrix are almost equal
 pub fn mat_mat_approx_eq<T>(a: MatRef<T>, b: MatRef<T>, tol: T)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     use assert_approx_eq::assert_approx_eq;
     assert_eq!(a.ncols(), b.ncols());
@@ -545,7 +556,7 @@ pub fn mat_mat_approx_eq<T>(a: MatRef<T>, b: MatRef<T>, tol: T)
 // Helper function to ensure all matrix ele are close to const
 pub fn mat_scale_approx_eq<T>(a: MatRef<T>, b: T, tol: T)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     use assert_approx_eq::assert_approx_eq;
     for j in 0..a.ncols() {
@@ -559,7 +570,7 @@ pub fn mat_scale_approx_eq<T>(a: MatRef<T>, b: T, tol: T)
 // Stack two matrix together, horizontally
 pub fn mat_hstack<T>(a: MatRef<T>, b: MatRef<T>) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     assert!(a.nrows() == b.nrows());
     let o_nrows = a.nrows();
@@ -583,7 +594,7 @@ pub fn mat_hstack<T>(a: MatRef<T>, b: MatRef<T>) -> Mat<T>
 // Stack two matrix together, vertically
 pub fn mat_vstack<T>(a: MatRef<T>, b: MatRef<T>) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     assert!(a.ncols() == b.ncols());
     let o_nrows = a.nrows() + b.nrows();
@@ -606,7 +617,7 @@ pub fn mat_vstack<T>(a: MatRef<T>, b: MatRef<T>) -> Mat<T>
 /// Produces a single column vector of elements evenly spaced
 pub fn mat_linspace<T>(start: T, end: T, n_steps: usize) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut o_mat = faer::Mat::zeros(n_steps, 1);
     let delta = (end - start) / (T::from(n_steps).unwrap());
@@ -620,7 +631,7 @@ pub fn mat_linspace<T>(start: T, end: T, n_steps: usize) -> Mat<T>
 /// Modifies the input matrix in-place
 pub fn mat_set_col<T>(mut a_mat: MatMut<T>, col: usize, col_mat: MatRef<T>)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     for i in 0..col_mat.nrows() {
         a_mat.write(i, col, col_mat.read(i, 0) );
@@ -630,7 +641,7 @@ pub fn mat_set_col<T>(mut a_mat: MatMut<T>, col: usize, col_mat: MatRef<T>)
 /// Converts a rust Vec into a faer Mat.  Done by a data copy (can be expensive)
 pub fn mat_from_vec<T>(in_vec: &Vec<T>) -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     let mut o_mat: Mat<T> = faer::Mat::zeros(in_vec.len(), 1);
     for (i, ele) in in_vec.into_iter().enumerate() {
diff --git a/src/lib_math_utils/pca_rsvd.rs b/src/lib_math_utils/pca_rsvd.rs
index 620c2bf..961a9ac 100644
--- a/src/lib_math_utils/pca_rsvd.rs
+++ b/src/lib_math_utils/pca_rsvd.rs
@@ -1,8 +1,8 @@
 // Principle component analysis (PCA) using
 // randomized svd methods.
 use std::cmp;
-use faer::{prelude::*};
-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::*;
diff --git a/src/lib_math_utils/pod_rom.rs b/src/lib_math_utils/pod_rom.rs
index 2fcc396..bc84608 100644
--- a/src/lib_math_utils/pod_rom.rs
+++ b/src/lib_math_utils/pod_rom.rs
@@ -7,8 +7,8 @@
 /// Ex: estimate the pressure field as a function of angle of attack
 /// using a few pre-computed CFD snapshots.  The POD model interpolates
 /// the CFD data to obtain an estimated pressure field at any attack angle.
-use faer::{prelude::*};
-use faer_core::{mat, Mat, MatRef, MatMut, Entity, AsMatRef, AsMatMut};
+use faer::prelude::*;
+use faer::{mat, Mat, MatRef, MatMut};
 use num_traits::Float;
 
 // internal imports
diff --git a/src/lib_math_utils/random_svd.rs b/src/lib_math_utils/random_svd.rs
index 672ebb0..7646cc3 100644
--- a/src/lib_math_utils/random_svd.rs
+++ b/src/lib_math_utils/random_svd.rs
@@ -1,11 +1,10 @@
 // Impl random SVD from:
 use std::cmp;
 use faer::{prelude::*};
-use faer_core::{mat, Mat, MatRef, MatMut, Entity, AsMatRef, AsMatMut};
+use faer::{mat, Mat, MatRef, MatMut};
 use rand::{prelude::*};
 use rand_distr::{StandardNormal, Uniform};
 use num_traits::Float;
-use std::time::SystemTime;
 
 // Internal imports
 use crate::lib_math_utils::mat_utils::*;
@@ -16,7 +15,7 @@ use crate::lib_math_utils::mat_utils::*;
 pub fn power_iter<T>(a_mat: MatRef<T>, omega_rank: usize, n_iter: usize)
     -> Mat<T>
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     // Generate random matrix, omega
     let a_ncols = a_mat.ncols();
@@ -64,7 +63,7 @@ pub fn power_iter<T>(a_mat: MatRef<T>, omega_rank: usize, n_iter: usize)
 pub fn random_svd<T>(a_mat: MatRef<T>, omega_rank: usize, n_iter: usize, n_oversamples: usize)
     -> (Mat<T>, Mat<T>, Mat<T>)
     where
-    T: faer_core::RealField + Float
+    T: faer::RealField + Float
 {
     // if matrix is fat, make thin
     let mut aa_mat = a_mat;
@@ -97,14 +96,14 @@ pub fn random_svd<T>(a_mat: MatRef<T>, omega_rank: usize, n_iter: usize, n_overs
     if fat == true {
         (
          my_rsvd.v().get(.., 0..omega_rank).to_owned(),
-         my_rsvd.s_diagonal().get(0..omega_rank, ..).to_owned(),
+         my_rsvd.s_diagonal().get(0..omega_rank).as_2d().to_owned(),
          u.transpose().get(0..omega_rank, ..).to_owned(),
         )
     }
     else {
         (
          u.get(.., 0..omega_rank).to_owned(),
-         my_rsvd.s_diagonal().get(0..omega_rank, ..).to_owned(),
+         my_rsvd.s_diagonal().get(0..omega_rank).as_2d().to_owned(),
          my_rsvd.v().get(.., 0..omega_rank).transpose().to_owned()
         )
     }
@@ -115,11 +114,12 @@ pub fn random_svd<T>(a_mat: MatRef<T>, omega_rank: usize, n_iter: usize, n_overs
 mod rsvd_unit_tests {
     // bring everything from above (parent) module into scope
     use super::*;
+    use std::time::SystemTime;
 
     #[test]
     fn test_rsvd_f64() {
     // set global parallelism
-    faer_core::set_global_parallelism(faer_core::Parallelism::Rayon(8));
+    faer::set_global_parallelism(faer::Parallelism::Rayon(8));
 
     // create random matrix
     let sys_timer = SystemTime::now();
@@ -132,7 +132,7 @@ mod rsvd_unit_tests {
     print!("done rng...{:?} s \n", (tf - ti).as_secs_f64());
 
     // print!("sigular values: {:?}", test_a);
-    let (ur, sr, vr) = random_svd(test_a.as_ref(), 4, 8, 10);
+    let (_ur, sr, _vr) = random_svd(test_a.as_ref(), 4, 8, 10);
     let tf2 = sys_timer.elapsed().unwrap();
     print!("done rsvd...{:?} s \n", (tf2 - tf).as_secs_f64());
     print!("sigular values: {:?}", sr);
diff --git a/src/lib_math_utils/space_samplers.rs b/src/lib_math_utils/space_samplers.rs
index 78337bc..29b78b2 100644
--- a/src/lib_math_utils/space_samplers.rs
+++ b/src/lib_math_utils/space_samplers.rs
@@ -487,11 +487,11 @@ mod space_samplers_unit_tests {
         let mut tst_mcmc_sampler = DeMcSampler::new(tst_ln_like_prior, tst_chains, 1, 0.8, 1.0e-10);
 
         // draw samples
-        let n_samples: usize = 1000;
+        let n_samples: usize = 5000;
         tst_mcmc_sampler.sample_mcmc(n_samples);
 
         // Check the estimated mean and var against known]
-        let tst_samples = tst_mcmc_sampler.get_samples(500);
+        let tst_samples = tst_mcmc_sampler.get_samples(2000);
         let ar = tst_mcmc_sampler.accept_ratio();
         println!("MCMC Samples: {:?}", tst_samples);
         println!("Accept ratio: {:?}", ar);
@@ -554,13 +554,13 @@ mod space_samplers_unit_tests {
         tst_mcmc_sampler.set_prop_fixup(proposal_fix_fn);
 
         // draw samples
-        let n_samples: usize = 1000;
+        let n_samples: usize = 5000;
         tst_mcmc_sampler.sample_mcmc_par(n_samples);
         let ar = tst_mcmc_sampler.accept_ratio();
         println!("Accept ratio: {:?}", ar);
 
         // TODO: check samples
-        let tst_samples = tst_mcmc_sampler.get_samples(500/8);
+        let tst_samples = tst_mcmc_sampler.get_samples(2000/8);
         println!("MCMC Diriclet Samples: {:?}", tst_samples);
         // assert_eq!(tst_samples.nrows(), n_samples);
         for (_i, sample) in tst_samples.axis_iter(Axis(0)).enumerate() {
diff --git a/src/lib_math_utils/stats_corr.rs b/src/lib_math_utils/stats_corr.rs
index a5085df..d6e18c6 100644
--- a/src/lib_math_utils/stats_corr.rs
+++ b/src/lib_math_utils/stats_corr.rs
@@ -1,8 +1,8 @@
 /// Calculate correlation coefficients and
 /// compute and manipulate covariance matricies
 ///
-use faer::{prelude::*};
-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::*;
diff --git a/src/lib_math_utils/univariate_rv.rs b/src/lib_math_utils/univariate_rv.rs
index 2b10819..5f40123 100644
--- a/src/lib_math_utils/univariate_rv.rs
+++ b/src/lib_math_utils/univariate_rv.rs
@@ -5,11 +5,10 @@
 use assert_approx_eq::assert_approx_eq;
 use rand_distr::{Normal, Distribution, Beta, Exp};
 use rand::Rng;
-use rand::seq::index;
 use rayon::prelude::*;
 use itertools::Itertools;
 use ndarray::prelude::*;
-use statrs::function::{erf, gamma::{gamma, self}, beta::beta_reg};
+use statrs::function::{erf, gamma::gamma, beta::beta_reg};
 use std::f64::consts::PI;
 use finitediff::FiniteDiff;
 use argmin::{
diff --git a/src/lib_math_utils_py.rs b/src/lib_math_utils_py.rs
index 68adc20..b03b92d 100644
--- a/src/lib_math_utils_py.rs
+++ b/src/lib_math_utils_py.rs
@@ -3,8 +3,9 @@ use numpy::{IntoPyArray, PyArray1, PyArray2, PyReadonlyArray2, PyReadonlyArray1,
 use pyo3::{exceptions::PyRuntimeError, pyclass, pymodule, types::PyModule, PyResult, Python};
 use pyo3::prelude::*;
 
-use faer::{mat, Mat, MatRef, IntoFaer, IntoNdarray};
-use faer::{prelude::*};
+use faer::prelude::*;
+use faer::Mat;
+use faer_ext::*;
 use crate::lib_math_utils::dmd_rom::*;
 use crate::lib_math_utils::space_samplers::*;
 use crate::lib_math_utils::pod_rom::*;
@@ -14,12 +15,12 @@ use crate::lib_math_utils::pca_rsvd::{PcaRsvd};
 use crate::lib_math_utils::active_subspaces::{ActiveSsRsvd, PolyGradientEstimator};
 
 #[pymodule]
-fn corrla_rs<'py>(_py: Python<'py>, m: &'py PyModule)
+fn corrla_rs<'py>(_py: Python<'py>, m: &Bound<'py, PyModule>)
     -> PyResult<()>
 {
     #[pyfn(m)]
     fn rsvd<'py>(py: Python<'py>, a_mat: PyReadonlyArray2<'py, f64>, n_rank: usize, n_iters: usize, n_oversamples: usize)
-        -> (&'py PyArray2<f64>, &'py PyArray2<f64>, &'py PyArray2<f64>)
+        -> (Bound<'py, PyArray2<f64>>, Bound<'py, PyArray2<f64>>, Bound<'py, PyArray2<f64>>)
     {
         // convert numpy array into rust ndarray and
         // convert to faer-rs matrix
@@ -31,7 +32,7 @@ fn corrla_rs<'py>(_py: Python<'py>, m: &'py PyModule)
         let ndarray_sr: Array2<f64> = sr.as_ref().into_ndarray().to_owned();
         let ndarray_ur: Array2<f64> = ur.as_ref().into_ndarray().to_owned();
         let ndarray_vr: Array2<f64> = vr.as_ref().into_ndarray().to_owned();
-        (ndarray_ur.into_pyarray(py), ndarray_sr.into_pyarray(py), ndarray_vr.into_pyarray(py))
+        (ndarray_ur.into_pyarray_bound(py), ndarray_sr.into_pyarray_bound(py), ndarray_vr.into_pyarray_bound(py))
     }
 
     #[pyfn(m)]
diff --git a/tests/test_a.rs b/tests/test_a.rs
index d7e707e..7d21115 100644
--- a/tests/test_a.rs
+++ b/tests/test_a.rs
@@ -6,7 +6,5 @@ mod integration_tests {
 
     #[test]
     fn it_works() {
-        let result = my_add(2, 2);
-        assert_eq!(result, 4);
     }
 }