Feature/better docu and usability (#9)

* minor cleanups in code and major additions to docs

* more docs

* document many things

* new readme

* amend readme

* start documenting separable model builder

* document separable model builder

* correct documentation

* bump version to 0.1.0

* mention nalgebra in readme

* apply fmt

Cargo.toml

@@ -1,12 +1,12 @@
 [package]
 name = "varpro"
-version = "0.0.0"
+version = "0.1.0"
 authors = ["geo-ant"]
 edition = "2018"
 license = "MIT"
 homepage = "https://github.com/geo-ant/varpro"
 repository = "https://github.com/geo-ant/varpro"
-description = "An easy to use nonlinear least-squares fitting library which makes use of the Variable Projection algorithm."
+description = "A straightforward nonlinear least-squares fitting library which uses the Variable Projection algorithm."
 readme = "README.md"
 categories = ["mathematics","science","algorithms"]
 keywords = ["nonlinear","least","squares","fitting","regression"]
@@ -24,5 +24,5 @@ approx = "0.4"
 
 [package.metadata.docs.rs]
 # To build locally use
-#   RUSTDOCFLAGS="--html-in-header katex-header.html" cargo doc --no-deps --document-private-items --open
+#   RUSTDOCFLAGS="--html-in-header katex-header.html" cargo doc --no-deps --open
 rustdoc-args = [ "--html-in-header", "katex-header.html" ]

README.md

@@ -5,6 +5,82 @@
 ![lints](https://github.com/geo-ant/varpro/workflows/lints/badge.svg?branch=main)
 [![coverage](https://codecov.io/gh/geo-ant/varpro/branch/main/graph/badge.svg?token=1L2PEJFMXP)](https://codecov.io/gh/geo-ant/varpro)
 
-_Work in Progress_
+This library enables robust and fast least-squares fitting of nonlinear, separable model functions to observations. It uses the VarPro algorithm to achieve this.
 
-A work in progress striving to bring separable nonlinear model fitting using the Variable Projection algorithm to Rust. Aims to combine an efficient implementation with a clean and usable interface.
+## Brief Introduction
+The lack of formulas on this site makes it hard to get into the depth of the what and how of this crate at this point. [Refer to the documentation](https://docs.rs/varpro/) for all the meaty details including the math.
+
+### What are Separable Model Functions?
+Put simply, separable model functions are nonlinear functions that can be written as a *linear combination* of some *nonlinear* basis functions. A common use case for VarPro is e.g. fitting sums of exponentials, which is a notoriously ill-conditioned problem.
+
+### What is VarPro?
+Variable Projection (VarPro) is an algorithm that takes advantage of the fact that the fitting problem can be separated into linear and truly nonlinear parameters. The linear parameters are eliminated and the fitting problem is cast into a problem that only depends on the nonlinear parameters. This reduced problem is then solved by using a common nonlinear fitting algorithm, such as Levenberg-Marquardt (LM). The varpro library uses the [levenberg-marquardt](https://crates.io/crates/levenberg-marquardt) crate as a backend for the LM minimizer. This crate further uses [nalgebra](https://crates.io/crates/nalgebra) for vector and matrix calculations.
+
+### When Should You Give it a Try?
+VarPro can dramatically increase the robustness and speed of the fitting process compared to using the nonlinear approach directly. When
+* the model function you want to fit is a linear combination of nonlinear functions
+* *and* you know the analytical derivatives of all those functions
+
+*then* you should give it a whirl.
+
+## Example
+The following example shows how to use varpro to fit a double exponential decay with constant offset to a data vector `y` obtained at grid points `x`. [Refer to the documentation](https://docs.rs/varpro/) for a more in-depth guide.
+
+The exponential decay and it's derivative are given as:
+
+```rust
+use nalgebra::DVector;
+fn exp_decay(x :&DVector<f64>, tau : f64) -> DVector<f64> {
+  x.map(|x|(-x/tau).exp())
+}
+
+fn exp_decay_dtau(tvec: &DVector<f64>,tau: f64) -> DVector<f64> {
+  tvec.map(|t| (-t / tau).exp() * t / tau.powi(2))
+}
+```
+
+The steps to perform the fitting are:
+
+```rust
+use varpro::prelude::*;
+use varpro::solvers::levmar::{LevMarProblemBuilder, LevMarSolver};
+// 1. create the model by giving only the nonlinear parameter names it depends on
+let model = SeparableModelBuilder::<f64>::new(&["tau1", "tau2"])
+  // add the first exponential decay and its partial derivative to the model
+// give all parameter names that the function depends on
+// and subsequently provide the partial derivative for each parameter
+  .function(&["tau1"], exp_decay)
+  .partial_deriv("tau1", exp_decay_dtau)
+  // add the second exponential decay and its partial derivative to the model
+  .function(&["tau2"], exp_decay)
+  .partial_deriv("tau2", exp_decay_dtau)
+  // add the constant as a vector of ones as an invariant function
+  .invariant_function(|x|DVector::from_element(x.len(),1.))
+  // build the model
+  .build()
+  .unwrap();
+// 2. Cast the fitting problem as a nonlinear least squares minimization problem
+let problem = LevMarProblemBuilder::new()
+  .model(&model)
+  .x(x)
+  .y(y)
+  .initial_guess(&[1., 2.])
+  .build()
+  .expect("Building valid problem should not panic");
+// 3. Solve the fitting problem
+let (solved_problem, report) = LevMarSolver::new().minimize(problem);
+assert!(report.termination.was_successful());
+// 4. obtain  the nonlinear parameters after fitting
+// they are in the same order as the parameter names given to the model
+let alpha = solved_problem.params();
+// the linear coefficients after fitting
+// the are in the same order as the basis functions that were added to the model
+let c = solved_problem.linear_coefficients().unwrap();
+```
+
+## References and Further Reading
+(O'Leary2013) O’Leary, D.P., Rust, B.W. Variable projection for nonlinear least squares problems. *Comput Optim Appl* **54**, 579–593 (2013). DOI: [10.1007/s10589-012-9492-9](https://doi.org/10.1007/s10589-012-9492-9)
+
+**attention**: the O'Leary paper contains errors that are fixed (so I hope) in [this blog article](https://geo-ant.github.io/blog/2020/variable-projection-part-1-fundamentals/) of mine.
+
+(Golub2003) Golub, G. , Pereyra, V Separable nonlinear least squares: the variable projection method and its applications. Inverse Problems **19** R1 (2003) [https://iopscience.iop.org/article/10.1088/0266-5611/19/2/201](https://iopscience.iop.org/article/10.1088/0266-5611/19/2/201)

src/basis_function/detail.rs

@@ -16,11 +16,11 @@ where
     ScalarType: Scalar,
     Func: Fn(&DVector<ScalarType>, ScalarType) -> DVector<ScalarType>,
 {
-    fn eval(&self, location: &DVector<ScalarType>, params: &[ScalarType]) -> DVector<ScalarType> {
+    fn eval(&self, x: &DVector<ScalarType>, params: &[ScalarType]) -> DVector<ScalarType> {
         if params.len() != Self::ARGUMENT_COUNT {
             panic!("Basisfunction expected {} arguments but the provided parameter slice has {} elements.",Self::ARGUMENT_COUNT,params.len());
         }
-        (&self)(location, params[0].clone())
+        (&self)(x, params[0].clone())
     }
 
     const ARGUMENT_COUNT: usize = 1;
@@ -35,7 +35,8 @@ where
 // i.e. these example expressions evaluate to true:
 // count_args(1,2,3)==count_args(1,2,3,) == 3usize, as well as
 // count_args(a,b)==count_args(a,b,)==2usize
-//TODO: see also daniel keep post (see above): recursion is not the most efficient way to implement this!
+//TODO (Performance): maybe see also daniel keep post (see above): recursion is not the most efficient way to implement this!
+// however, the alternative solutions can't be used in constant expressions, iirc.
 macro_rules! count_args {
     () => {0usize};
     ($_head:tt, $($tail:tt),*) => {1usize + count_args!($($tail,)*)}; //these overloads are there to allow with and w/o trailing comma syntax
@@ -111,4 +112,4 @@ basefunction_impl_helper!(
     (8, T),
     (9, T)
 ); //10 parameter arguments
-   //todo: if more are implemented, add tests as well
+   //if more are implemented, add tests as well

src/basis_function/mod.rs

@@ -5,8 +5,8 @@ mod detail;
 
 use nalgebra::{DVector, Scalar};
 
-/// This trait allows us to pass vector valued functions `$f(\vec{x},\alpha_1,...\alpha_N$)` in a generic fashion,
-/// where
+/// This trait allows us to pass vector valued functions `$\vec{f}(\vec{x},\alpha_1,...\alpha_N$)` in a generic fashion, where
+///
 /// * `$\vec{x}$` is an independent variable, like time, spatial coordinates and so on,
 /// * `$\alpha_j$` are scalar parameter of this function
 /// The functions must have at least one parameter argument additional ot the location parameter, i.e. a function
@@ -15,12 +15,11 @@ use nalgebra::{DVector, Scalar};
 /// add this function to a model this is actually a requirement that will lead to errors, even panics
 /// if violated).
 /// Unfortunately, requirement on the length of the output vector cannot be enforced by the type system.
-/// If it is violated, then calculations using the basis function will fail in the [SeparableModel].
+/// If it is violated, then calculations using the basis function will return errors in the [SeparableModel](crate::model::SeparableModel).
 ///
 /// # Variadic Functions
 /// Since Rust does not have variadic functions or generics, this trait is implemented for all functions up to
-/// a maximum number of arguments. This maximum number of arguments can be found out by checking
-/// the blanket implementations.
+/// a maximum number of arguments. Right now this is implemented for up to 10 arguments.
 ///
 /// # Generic Parameters
 /// ## Scalar Type
@@ -37,21 +36,26 @@ pub trait BasisFunction<ScalarType, ArgList>
 where
     ScalarType: Scalar + Clone,
 {
-    /// A common calling interface to evaluate this function by passing a slice of scalar types
-    /// that is dispatched to the arguments in order. The slice must have the same
-    /// length as the parameter argument list.
+    /// A common calling interface to evaluate this function `$\vec{f}(\vec{x},\alpha_1,...\alpha_N)$` by passing a slice of scalar types
+    /// that is dispatched to the arguments in order.
+    ///
+    /// # Arguments
+    /// * `x`: the vector `$\vec{x}$`
+    /// * `params`: The parameters `$(\alpha_1,...,\alpha_N)$` as a slice. The slice must have
+    /// the correct number of arguments for calling the function (no more, no less).
+    ///
     /// # Effect and Panics
     /// If the slice has fewer elements than the parameter argument list. If the slice has more element,
     /// only the first `$N$` elements are used and dispatched to the parameters in order, i.e.
     /// `$\alpha_1$`=`param[0]`, ..., `$\alpha_N$`=`param[N-1]`. Calling eval will result in evaluating
     /// the implementing callable at the location and arguments given.
     ///
-    /// **Attention** The library takes care that no panics can be caused by calls to `eval` (from
+    /// **Attention** The `varpro` library takes care that no panics can be caused by calls to `eval` (from
     /// within this library) by making sure the parameter slice has the correct length. Therefore
     /// it is not recommended (and not necessary) to use this function in other code than inside this library.
-    fn eval(&self, location: &DVector<ScalarType>, params: &[ScalarType]) -> DVector<ScalarType>;
+    fn eval(&self, x: &DVector<ScalarType>, params: &[ScalarType]) -> DVector<ScalarType>;
 
     /// This gives the number of parameter arguments to the callable. So for a function `$f(\vec{x},\alpha_1,...\alpha_N)$`
-    /// this is give `N`.
+    /// this is equal to `$N$`.
     const ARGUMENT_COUNT: usize;
 }

src/basis_function/test.rs

@@ -1,11 +1,11 @@
 use super::*;
 
-/// a function taking a location vector and 1 parameter argument
+/// a function taking a location vector x and 1 parameter argument
 fn callable1(x: &DVector<f64>, alpha: f64) -> DVector<f64> {
     alpha * x.clone()
 }
 
-/// a function taking a location vector and 2 parameter arguments
+/// a function taking a location vector y  and 2 parameter arguments
 fn callable2(x: &DVector<f64>, alpha1: f64, alpha2: f64) -> DVector<f64> {
     alpha1 * alpha2 * x.clone()
 }