Confidence bands and Usability Improvements (#33)

* add dev branch checks

* Feature/remove typestate pattern from model builder (#30)

* break everything, but now we have a state machine

* more work on removing typestate

* refactor should be finished

* fix clippy lints

* Feature/confidence bands (#32)

* start with confidence bands

* update comment

* don't store correlation matrix anymore, instead calculate it on the fly

* disable fit statistics for mrhs because that was not working correctly

* start adding calculations

* minor changes to calculations

* finish draft for confidence bands

* add generics for mrhs vs single rhs

* compiling, but doctests are failing

* offer different APIs for single and multiple rhs

* single vector api in fit statistics

* compile and tests working, doctests still fail

* remove obsolete code

* add best fit method and start testing it

* add more tests for best fit

* more tests for best fit

* add docs for confidence bands

* fix doctests

* start changing docs to better reflect mrhs

* start with python script using lmfit for comparison

* fiddle with parameters until fit works, add random noise for ci calculation

* minor cleanups in script

* write results

* start with tests for confidence band

* add x data to output

* more test assets

* test and fix bugs in confidence band

* move some test assets around

* add weighted decay

* test fitting with weights, found problem with covmat

* smaller refactor

* use correct cov matrix for weighted problem

* shorten todo comment

* use correct conf interval, fix test

* doctest readme, fix problems

* increment version

* fmt

* doc fixes

* add reduced chi2 and add comment about scaling

* test reduced chi2

* update readme

* add todo list

* update changelog, add todo list

* more documentation

* add test for the remove typestate feature

* more doc

* overhaul readme again and add todos

* more corrections in readme, append todo

.github/workflows/build.yml

@@ -4,7 +4,7 @@ on:
   push:
     branches: [ main ]
   pull_request:
-    branches: [ main, dmz ]
+    branches: [ main, dev ]
 
 env:
   CARGO_TERM_COLOR: always

.github/workflows/coverage.yml

@@ -4,7 +4,7 @@ on:
   push:
     branches: [ main ]
   pull_request:
-    branches: [ main, dmz ]
+    branches: [ main, dev ]
 
 env:
   RUST_BACKTRACE: 1

.github/workflows/lints.yml

@@ -4,7 +4,7 @@ on:
   push:
     branches: [ main ]
   pull_request:
-    branches: [ main, dmz ]
+    branches: [ main, dev]
 
 env:
   CARGO_TERM_COLOR: always

.github/workflows/tests.yml

@@ -4,7 +4,7 @@ on:
   push:
     branches: [ main ]
   pull_request:
-    branches: [ main, dmz ]
+    branches: [ main, dev ]
 
 env:
   CARGO_TERM_COLOR: always

CHANGES.md

@@ -3,6 +3,18 @@
 This is the changelog for the `varpro` library.
 See also here for a [version history](https://crates.io/crates/varpro/versions).
 
+## 0.9.0
+
+- We can now calculate confidence bands using via `FitStatistics`. Only available
+for problems with a single RHS, though. 
+- The API for calculating the correlation matrix changed to on-the-fly
+calculations, deprecated (a drop-in replacement) for the old API.
+- Cleaner separation (API changes) for problems with single and multiple right
+hand sides.
+- Removed some overly clever (but ultimately bad) generics from the `SeparableModelBuilder`
+that were making it hard to use to add functions in a loop. State machine is now
+a runtime thing.
+
 ## 0.8.0 Multiple Right Hand Sides
 
 - Observations can now be a matrix, in this case a global fit with multiple right

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "varpro"
-version = "0.8.0"
+version = "0.9.0"
 authors = ["geo-ant"]
 edition = "2021"
 license = "MIT"
@@ -19,6 +19,7 @@ thiserror = "1"
 levenberg-marquardt = "0.13"
 nalgebra = { version = "0.32" } #, features = ["rayon"]}
 num-traits = "0.2"
+distrs = "0.2"
 # rayon = "1.6"
 
 [dev-dependencies]
@@ -30,6 +31,7 @@ shared_test_code = { path = "./shared_test_code" }
 assert_matches = "1.5"
 mockall = "0.11"
 rand = "0.8"
+byteorder = "1.5"
 
 [[bench]]
 name = "double_exponential_without_noise"

README.md

@@ -11,77 +11,94 @@ Nonlinear function fitting made simple. This library provides robust and fast
 least-squares fitting of a wide class of model functions to data.
 It uses the VarPro algorithm to achieve this, hence the name.
 
-## Brief Introduction
+## Introduction
 
-This crate implements a powerful algorithm
-to fit model functions to data, but it is restricted to so called _separable_
-models. See the next section for an explanation. 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.
+This crate implements a powerful algorithm to fit model functions to data, 
+but it is restricted to so called _separable_ models. The lack of formulas on 
+this  site makes it hard to go into detail, but a brief overview is provided in
+the next sections. [Refer to the documentation](https://docs.rs/varpro/) for all
+the meaty details including the math.
 
 ### What are Separable Models?
 
-Put simply, separable models are nonlinear functions that can be 
+Put simply, separable models are nonlinear functions which 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 given fitting problem can be separated into linear and truly nonlinear parameters.
-The linear parameters are eliminated using some clever linear algebra
-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,
+Variable Projection (VarPro) is an algorithm that exploits that its fitting
+problem can be separated into linear and nonlinear parameters.
+First, the linear parameters are eliminated using some clever linear algebra. Then,
+the fitting problem is rewritten so that it depends only on the nonlinear parameters.
+Finally, this reduced problem is solved by using a general purpose nonlinear minimization algorithm,
 such as Levenberg-Marquardt (LM).
 
 ### When Should You Give it a Try?
 
 VarPro can dramatically increase the robustness and speed of the fitting process
-compared to using a "normal" nonlinear least squares fitting algorithm. When
+compared to using a general purpose nonlinear least squares fitting algorithm. When
 
-* the model function you want to fit is a linear combination of nonlinear functions
+* 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.
+_then_ you should give it a whirl. Also consider the section on global fitting below,
+which provides another great use case for this crate.
 
 ## Example Usage
 
-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`. 
+The following example shows, how to use this crate to fit a double exponential decay
+with constant offset to a data vector `y` obtained at time points `t`. 
 [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())
+use varpro::prelude::*;
+use varpro::solvers::levmar::{LevMarProblemBuilder, LevMarSolver};
+use nalgebra::{dvector,DVector};
+
+// Define the exponential decay e^(-t/tau).
+// Both of the nonlinear basis functions in this example
+// are exponential decays.
+fn exp_decay(t :&DVector<f64>, tau : f64) 
+  -> DVector<f64> {
+  t.map(|t|(-t/tau).exp())
 }
 
-fn exp_decay_dtau(tvec: &DVector<f64>,tau: f64) -> DVector<f64> {
-  tvec.map(|t| (-t / tau).exp() * t / tau.powi(2))
+// the partial derivative of the exponential
+// decay with respect to the nonlinear parameter tau.
+// d/dtau e^(-t/tau) = e^(-t/tau)*t/tau^2
+fn exp_decay_dtau(t: &DVector<f64>,tau: f64) 
+  -> DVector<f64> {
+  t.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};
-
-let x = /*time or spatial coordinates of the observations*/;
-let y = /*the observed data we want to fit*/;
+// temporal (or spatial) coordintates of the observations
+let t = dvector![0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.];
+// the observations we want to fit
+let y = dvector![6.0,4.8,4.0,3.3,2.8,2.5,2.2,1.9,1.7,1.6,1.5];
 
 // 1. create the model by giving only the nonlinear parameter names it depends on
 let model = SeparableModelBuilder::<f64>::new(&["tau1", "tau2"])
+  // provide the nonlinear basis functions and their derivatives.
+  // In general, base functions can depend on more than just one parameter.
+  // first function:
   .function(&["tau1"], exp_decay)
   .partial_deriv("tau1", exp_decay_dtau)
+  // second function and derivatives with respect to all parameters
+  // that it depends on (just one in this case)
   .function(&["tau2"], exp_decay)
   .partial_deriv("tau2", exp_decay_dtau)
-  .invariant_function(|x|DVector::from_element(x.len(),1.))
-  .independent_variable(x)
-  .initial_parameters(initial_params)
+  // a constant offset is added as an invariant basefunction
+  // as a vector of ones. It is multiplied with its own linear coefficient,
+  // creating a fittable offset
+  .invariant_function(|v|DVector::from_element(v.len(),1.))
+  // give the coordinates of the problem
+  .independent_variable(t)
+  // provide guesses only for the nonlinear parameters in the
+  // order that they were given on construction.
+  .initial_parameters(vec![2.5,5.5])
   .build()
   .unwrap();
 // 2. Cast the fitting problem as a nonlinear least squares minimization problem
@@ -90,7 +107,7 @@ let problem = LevMarProblemBuilder::new(model)
   .build()
   .unwrap();
 // 3. Solve the fitting problem
-let fit_result = LevMarSolver::new()
+let fit_result = LevMarSolver::default()
     .fit(problem)
     .expect("fit must exit successfully");
 // 4. obtain the nonlinear parameters after fitting
@@ -99,37 +116,42 @@ let alpha = fit_result.nonlinear_parameters();
 let c = fit_result.linear_coefficients().unwrap();
 ```
 
-For more examples please refer to the crate documentation.
+For more in depth examples, please refer to the crate documentation.
 
 ### Fit Statistics
 
-Additionally to the `fit` member function, the `LevMarSolver` also provides a 
-`fit_with_statistics` function that calculates some additional statistical
-information after the fit has finished.
+Additionally to the `fit` member function, the `LevMarSolver` provides a 
+`fit_with_statistics` function that calculates quite a bit of useful additional statistical
+information.
 
 ### Global Fitting of Multiple Right Hand Sides
 
-Before, we have passed a single column vector as the observations. It is also
-possible to pass a matrix, whose columns represent individual observations. We
-are now fitting a problem with multiple right hand sides. `vapro` will performa a _global fit_
-in which the nonlinear parameters are optimized across all right hand sides, but
-linear parameters of the fit are optimized for each right hand side individually.
+In the example above, we have passed a single column vector as the observations.
+The library also allows fitting multiple right hand sides, by constructing a
+problem via `LevMarProblemBuilder::mrhs`. When fitting multiple right hand sides,
+`vapro` will performa a _global fit_, in which the nonlinear parameters are optimized
+across all right hand sides, but linear coefficients of the fit are optimized for
+each right hand side individually.
 
-This is an application where varpro really shines because it can take advantage
-of the separable nature of the problem. It allows us to perform a global fit over thousands
+This is another application where varpro really shines, since it can take advantage
+of the separable nature of the problem. It allows us to perform a global fit over thousands,
 or even tens of thousands of right hand sides in reasonable time (fractions of seconds to minutes),
 where conventional purely nonlinear solvers must perform much more work.
 
-### Maximum Performance
+### Maximum Performance and Advanced Use Cases
 
 The example code above will already run many times faster
 than just using a nonlinear solver without the magic of varpro.
-But this crate offers an additional way to eek out the last bits of  performance.
+But this crate offers an additional way to eek out the last bits of performance.
+
+The `SeparableNonlinearModel` trait can be manually implemented to describe a
+model function. This often allows us to shave of the last hundreds of microseconds
+from the computation, e.g. by caching intermediate calculations. The crate documentation
+contains detailed examples.
 
-The `SeparableNonlinearModel` trait can be manually
-implemented to describe a model function. This often allows us to shave of the 
-last hundreds of microseconds from the computation e.g. by caching intermediate
-calculations. The crate documentation contains detailed examples.
+This is not only useful for performance, but also for use cases that are difficult
+or impossible to accomodate using only the `SeparableModelBuilder`. The builder
+was created for ease of use _and_ performance, but it has some limitations by design.
 
 ## Acknowledgements
 

Todo.md

@@ -0,0 +1,10 @@
+# ToDo List
+
+[ ] Better interface for fit result: successful fits should not have optional
+linear coefficients. Use const generic bool like I did for the MRHS case?
+[ ] Fit statistics (and confidence bands, but also correlation matrix etc) for
+problems with multiple RHS
+[ ] Provide a more convenient way to add fittable offsets (plus some safeguards, such
+that offsets cannot be added twice). Also think of a better name than fittable offset,
+but make it clear that it is not just a constant offset, but
+one that comes with an extra parameter.

benches/double_exponential_without_noise.rs

@@ -3,6 +3,7 @@ use levenberg_marquardt::LeastSquaresProblem;
 use levenberg_marquardt::LevenbergMarquardt;
 use nalgebra::ComplexField;
 
+use nalgebra::DVector;
 use nalgebra::DefaultAllocator;
 
 use nalgebra::Dyn;
@@ -35,7 +36,7 @@ struct DoubleExponentialParameters {
 fn build_problem<Model>(
     true_parameters: DoubleExponentialParameters,
     mut model: Model,
-) -> LevMarProblem<Model>
+) -> LevMarProblem<Model, false>
 where
     Model: SeparableNonlinearModel<ScalarType = f64>,
     DefaultAllocator: nalgebra::allocator::Allocator<f64, Dyn>,
@@ -68,7 +69,7 @@ where
         .expect("Building valid problem should not panic")
 }
 
-fn run_minimization<Model>(problem: LevMarProblem<Model>) -> [f64; 5]
+fn run_minimization<Model>(problem: LevMarProblem<Model, false>) -> (DVector<f64>, DVector<f64>)
 where
     Model: SeparableNonlinearModel<ScalarType = f64> + std::fmt::Debug,
 {
@@ -77,7 +78,7 @@ where
         .expect("fitting must exit successfully");
     let params = result.nonlinear_parameters();
     let coeff = result.linear_coefficients().unwrap();
-    [params[0], params[1], coeff[0], coeff[1], coeff[2]]
+    (params, coeff.into_owned())
 }
 
 /// solve the problem by using nonlinear least squares with levenberg marquardt

benches/multiple_right_hand_sides.rs

@@ -23,21 +23,21 @@ struct DoubleExponentialParameters {
 fn build_problem_mrhs<Model>(
     true_parameters: DoubleExponentialParameters,
     mut model: Model,
-) -> LevMarProblem<Model>
+) -> LevMarProblem<Model, true>
 where
     Model: SeparableNonlinearModel<ScalarType = f64>,
 {
     let DoubleExponentialParameters { tau1, tau2, coeffs } = true_parameters.clone();
     // save the initial guess so that we can reset the model to those
     let params = OVector::from_vec_generic(Dyn(model.parameter_count()), U1, vec![tau1, tau2]);
     let y = evaluate_complete_model_at_params_mrhs(&mut model, params, &coeffs);
-    LevMarProblemBuilder::new(model)
+    LevMarProblemBuilder::mrhs(model)
         .observations(y)
         .build()
         .expect("Building valid problem should not panic")
 }
 
-fn run_minimization<Model>(problem: LevMarProblem<Model>) -> (DVector<f64>, DMatrix<f64>)
+fn run_minimization_mrhs<Model>(problem: LevMarProblem<Model, true>) -> (DVector<f64>, DMatrix<f64>)
 where
     Model: SeparableNonlinearModel<ScalarType = f64> + std::fmt::Debug,
 {
@@ -46,7 +46,7 @@ where
         .expect("fitting must exit successfully");
     let params = result.nonlinear_parameters();
     let coeff = result.linear_coefficients().unwrap();
-    (params, coeff)
+    (params, coeff.into_owned())
 }
 
 fn bench_double_exp_no_noise_mrhs(c: &mut Criterion) {
@@ -75,7 +75,7 @@ fn bench_double_exp_no_noise_mrhs(c: &mut Criterion) {
                     DoubleExpModelWithConstantOffsetSepModel::new(x.clone(), tau_guess),
                 )
             },
-            run_minimization,
+            run_minimization_mrhs,
             criterion::BatchSize::SmallInput,
         )
     });
@@ -91,7 +91,7 @@ fn bench_double_exp_no_noise_mrhs(c: &mut Criterion) {
                     ),
                 )
             },
-            run_minimization,
+            run_minimization_mrhs,
             criterion::BatchSize::SmallInput,
         )
     });

python/.gitignore

@@ -0,0 +1 @@
+*.raw