Feature/bump dependency versions 2 (#15)

* depend on newest levmar and nalgebra versions

* fix deprecation warnings

* add a test for diag matrix creation

* update changelog, bump version

Author: geo-ant

CHANGES.md

@@ -2,6 +2,10 @@
 
 This is the changelog for the `varpro` library. See also here for a [version history](https://crates.io/crates/varpro/versions).
 
+## 0.5.0
+- Upgrade nalgebra and levenberg_marquardt dependencies to current versions
+- Fix deprecation warnings and expand test coverage slightly
+
 ## 0.4.1
 - Remove snafu dependency and replace by `thiserror`
 - Use uninitialized matrices instead of zero initialisation for places where contents will be overwritten anyways

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "varpro"
-version = "0.4.1"
+version = "0.5.0"
 authors = ["geo-ant"]
 edition = "2021"
 license = "MIT"
@@ -15,12 +15,13 @@ keywords = ["nonlinear","least","squares","fitting","regression"]
 
 [dependencies]
 thiserror = "1"
-levenberg-marquardt = "0.12"
-nalgebra = "0.30"
+levenberg-marquardt = "0.13"
+nalgebra = "0.32"
 num-traits = "0.2"
 
 [dev-dependencies]
 approx = "0.5"
+num-complex = "0.4"
 
 [package.metadata.docs.rs]
 # To build locally use

src/lib.rs

@@ -1,5 +1,3 @@
-#![deny(missing_docs)]
-//! The varpro crate enables nonlinear least squares fitting for separable models using the Variable Projection (VarPro) algorithm.
 //!
 //! # Introduction
 //!

src/linalg_helpers/mod.rs

@@ -1,7 +1,7 @@
 #[cfg(test)]
 mod test;
 
-use nalgebra::{ClosedMul, ComplexField, DMatrix, DVector, Dynamic, Scalar};
+use nalgebra::{ClosedMul, ComplexField, DMatrix, DVector, Dyn, Scalar};
 use std::ops::Mul;
 
 /// A square diagonal matrix with dynamic dimension. Off-diagonal entries are assumed zero.
@@ -74,9 +74,8 @@ where
         // with valid values below. Ideally we would first call the
         // copy_from functionality below, but we can't because the Scalar
         // trait bounds will screw us. See https://github.com/dimforge/nalgebra/pull/949
-        let mut result_matrix = unsafe {
-            DMatrix::uninit(Dynamic::new(self.nrows()), Dynamic::new(rhs.ncols())).assume_init()
-        };
+        let mut result_matrix =
+            unsafe { DMatrix::uninit(Dyn(self.nrows()), Dyn(rhs.ncols())).assume_init() };
 
         for (col_idx, mut col) in result_matrix.column_iter_mut().enumerate() {
             col.copy_from(&self.diagonal.component_mul(&rhs.column(col_idx)));

src/linalg_helpers/test.rs

@@ -1,5 +1,6 @@
 use crate::linalg_helpers::DiagDMatrix;
-use nalgebra::{DMatrix, DVector};
+use approx::assert_relative_eq;
+use nalgebra::{ComplexField, DMatrix, DVector};
 
 #[test]
 #[allow(non_snake_case)]
@@ -74,3 +75,26 @@ fn diagonal_matrix_multiplication_must_panic_for_rhs_dimensions_too_large() {
     let A = DMatrix::from_element(ndiag - 1, 1, 1.);
     let _ = &D * &A;
 }
+
+#[test]
+#[allow(non_snake_case)]
+fn diagonal_matrix_test_from_real_field() {
+    use num_complex::Complex;
+    let diagonal = DVector::from(vec![1., 2., 3.]);
+    let D = DiagDMatrix::<Complex<f64>>::from_real_field(&diagonal);
+    let ones = DVector::from(vec![
+        Complex::new(1., 5.),
+        Complex::new(1., 5.),
+        Complex::new(1., 5.),
+    ]);
+    let result = &D * &ones;
+
+    assert_relative_eq!(
+        result.map(ComplexField::real),
+        DVector::from(vec![1., 2., 3.])
+    );
+    assert_relative_eq!(
+        result.map(ComplexField::imaginary),
+        DVector::from(vec![5., 10., 15.])
+    );
+}

src/model/mod.rs

@@ -1,7 +1,7 @@
 use crate::model::errors::ModelError;
 use crate::model::model_basis_function::ModelBasisFunction;
 use nalgebra::base::Scalar;
-use nalgebra::{DMatrix, DVector, Dynamic};
+use nalgebra::{DMatrix, DVector, Dyn};
 use num_traits::Zero;
 
 mod detail;
@@ -152,7 +152,7 @@ where
         // this pattern is not great, but the trait bounds in copy_from still
         // prevent us from doing something better
         let mut function_value_matrix =
-            unsafe { DMatrix::uninit(Dynamic::new(nrows), Dynamic::new(ncols)).assume_init() };
+            unsafe { DMatrix::uninit(Dyn(nrows), Dyn(ncols)).assume_init() };
 
         for (basefunc, mut column) in self
             .basefunctions

src/solvers/levmar/mod.rs

@@ -1,7 +1,7 @@
 use crate::model::SeparableModel;
 use levenberg_marquardt::LeastSquaresProblem;
 use nalgebra::storage::Owned;
-use nalgebra::{ComplexField, DMatrix, DVector, Dynamic, Matrix, Scalar, Vector, SVD};
+use nalgebra::{ComplexField, DMatrix, DVector, Dyn, Matrix, Scalar, Vector, SVD};
 
 mod builder;
 #[cfg(test)]
@@ -21,7 +21,7 @@ struct CachedCalculations<ScalarType: Scalar + ComplexField> {
     /// The current residual of model function values belonging to the current parameters
     current_residuals: DVector<ScalarType>,
     /// Singular value decomposition of the current function value matrix
-    current_svd: SVD<ScalarType, Dynamic, Dynamic>,
+    current_svd: SVD<ScalarType, Dyn, Dyn>,
     /// the linear coefficients `$\vec{c}$` providing the current best fit
     linear_coefficients: DVector<ScalarType>,
 }
@@ -80,22 +80,21 @@ impl<'a, ScalarType: Scalar + ComplexField + Copy> LevMarProblem<'a, ScalarType>
     }
 }
 
-impl<'a, ScalarType> LeastSquaresProblem<ScalarType, Dynamic, Dynamic>
-    for LevMarProblem<'a, ScalarType>
+impl<'a, ScalarType> LeastSquaresProblem<ScalarType, Dyn, Dyn> for LevMarProblem<'a, ScalarType>
 where
     ScalarType: Scalar + ComplexField + Copy,
     ScalarType::RealField: Mul<ScalarType, Output = ScalarType> + Float,
 {
-    type ResidualStorage = Owned<ScalarType, Dynamic>;
-    type JacobianStorage = Owned<ScalarType, Dynamic, Dynamic>;
-    type ParameterStorage = Owned<ScalarType, Dynamic>;
+    type ResidualStorage = Owned<ScalarType, Dyn>;
+    type JacobianStorage = Owned<ScalarType, Dyn, Dyn>;
+    type ParameterStorage = Owned<ScalarType, Dyn>;
 
     #[allow(non_snake_case)]
     /// Set the (nonlinear) model parameters `$\vec{\alpha}$` and update the internal state of the
     /// problem accordingly. The parameters are expected in the same order that the parameter
     /// names were provided in at model creation. So if we gave `&["tau","beta"]` as parameters at
     /// model creation, the function expects the layout of the parameter vector to be `$\vec{\alpha}=(\tau,\beta)^T$`.
-    fn set_params(&mut self, params: &Vector<ScalarType, Dynamic, Self::ParameterStorage>) {
+    fn set_params(&mut self, params: &Vector<ScalarType, Dyn, Self::ParameterStorage>) {
         self.model_parameters = params.iter().cloned().collect();
         // matrix of weighted model function values
         let Phi_w = self
@@ -135,7 +134,7 @@ where
     /// names given on model creation. E.g. if the parameters at model creation where given as
     /// `&["tau","beta"]`, then the returned vector is `$\vec{\alpha} = (\tau,\beta)^T$`, i.e.
     /// the value of parameter `$\tau$` is at index `0` and the value of `$\beta$` at index `1`.
-    fn params(&self) -> Vector<ScalarType, Dynamic, Self::ParameterStorage> {
+    fn params(&self) -> Vector<ScalarType, Dyn, Self::ParameterStorage> {
         DVector::from(self.model_parameters.clone())
     }
 
@@ -146,7 +145,7 @@ where
     /// algorithm calculates `$\vec{c}(\vec{\alpha})$` as the coefficients that provide the best linear least squares
     /// fit, given the current `$\vec{\alpha}$`. For more info on the math of VarPro, see
     /// e.g. [here](https://geo-ant.github.io/blog/2020/variable-projection-part-1-fundamentals/).
-    fn residuals(&self) -> Option<Vector<ScalarType, Dynamic, Self::ResidualStorage>> {
+    fn residuals(&self) -> Option<Vector<ScalarType, Dyn, Self::ResidualStorage>> {
         self.cached
             .as_ref()
             .map(|cached| cached.current_residuals.clone())
@@ -156,7 +155,7 @@ where
     /// Calculate the Jacobian matrix of the *weighted* residuals `$\vec{r}_w(\vec{\alpha})$`.
     /// For more info on how the Jacobian is calculated in the VarPro algorithm, see
     /// e.g. [here](https://geo-ant.github.io/blog/2020/variable-projection-part-1-fundamentals/).
-    fn jacobian(&self) -> Option<Matrix<ScalarType, Dynamic, Dynamic, Self::JacobianStorage>> {
+    fn jacobian(&self) -> Option<Matrix<ScalarType, Dyn, Dyn, Self::JacobianStorage>> {
         // TODO (Performance): make this more efficient by parallelizing
 
         if let Some(CachedCalculations {
@@ -168,11 +167,8 @@ where
             // this is not a great pattern, but the trait bounds on copy_from
             // as of now prevent us from doing something more idiomatic
             let mut jacobian_matrix = unsafe {
-                DMatrix::uninit(
-                    Dynamic::new(self.y_w.len()),
-                    Dynamic::new(self.model.parameter_count()),
-                )
-                .assume_init()
+                DMatrix::uninit(Dyn(self.y_w.len()), Dyn(self.model.parameter_count()))
+                    .assume_init()
             };
 
             let U = current_svd.u.as_ref()?; // will return None if this was not calculated

tests/integration_tests/levmar.rs

@@ -1,7 +1,7 @@
 use crate::common::linspace;
 use levenberg_marquardt::{differentiate_numerically, LeastSquaresProblem};
 use nalgebra::storage::Owned;
-use nalgebra::{DVector, Dynamic, Matrix, Vector, Vector5, U5};
+use nalgebra::{DVector, Dyn, Matrix, Vector, Vector5, U5};
 
 use approx::assert_relative_eq;
 
@@ -44,10 +44,10 @@ impl DoubleExponentialDecayFittingWithOffset {
     }
 }
 
-impl LeastSquaresProblem<f64, Dynamic, U5> for DoubleExponentialDecayFittingWithOffset {
+impl LeastSquaresProblem<f64, Dyn, U5> for DoubleExponentialDecayFittingWithOffset {
     type ParameterStorage = Owned<f64, U5>;
-    type ResidualStorage = Owned<f64, Dynamic>;
-    type JacobianStorage = Owned<f64, Dynamic, U5>;
+    type ResidualStorage = Owned<f64, Dyn>;
+    type JacobianStorage = Owned<f64, Dyn, U5>;
 
     fn set_params(&mut self, params: &Vector<f64, U5, Self::ParameterStorage>) {
         self.params = *params;
@@ -80,7 +80,7 @@ impl LeastSquaresProblem<f64, Dynamic, U5> for DoubleExponentialDecayFittingWith
         Some(&f - &self.y)
     }
 
-    fn jacobian(&self) -> Option<Matrix<f64, Dynamic, U5, Self::JacobianStorage>> {
+    fn jacobian(&self) -> Option<Matrix<f64, Dyn, U5, Self::JacobianStorage>> {
         // get parameters from internal param storage
         let tau1 = self.params[0];
         let tau2 = self.params[1];
@@ -89,7 +89,7 @@ impl LeastSquaresProblem<f64, Dynamic, U5> for DoubleExponentialDecayFittingWith
         // populate jacobian
         //let ncols = 5;
         let nrows = self.x.len();
-        let mut jacobian = Matrix::<f64, Dynamic, U5, Self::JacobianStorage>::zeros(nrows);
+        let mut jacobian = Matrix::<f64, Dyn, U5, Self::JacobianStorage>::zeros(nrows);
 
         jacobian.set_column(
             0,