From ea107b87b25e2d6e55c5027c0fdf99dc118c8550 Mon Sep 17 00:00:00 2001 From: Geo <54497890+geo-ant@users.noreply.github.com> Date: Tue, 21 May 2024 23:18:11 +0200 Subject: [PATCH] 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 | 2 +- .github/workflows/coverage.yml | 2 +- .github/workflows/lints.yml | 2 +- .github/workflows/tests.yml | 2 +- CHANGES.md | 12 + Cargo.toml | 4 +- README.md | 130 +++--- Todo.md | 10 + benches/double_exponential_without_noise.rs | 7 +- benches/multiple_right_hand_sides.rs | 12 +- python/.gitignore | 1 + python/multiexp_decay.py | 43 ++ python/weighted_multiexp_decay.py | 43 ++ src/lib.rs | 81 +++- src/model/builder/error.rs | 4 + src/model/builder/mod.rs | 384 +++++++----------- .../builder/modelfunction_builder/mod.rs | 2 +- src/model/builder/test.rs | 26 ++ src/model/model_basis_function.rs | 3 +- src/readme.rs | 5 + src/solvers/levmar/builder.rs | 103 +++-- src/solvers/levmar/mod.rs | 221 ++++++++-- src/solvers/levmar/test.rs | 9 - src/statistics/mod.rs | 294 ++++++++++++-- src/statistics/numeric_traits/mod.rs | 45 ++ src/statistics/test.rs | 4 +- src/util/mod.rs | 9 - .../multiexp_decay/conf_1000_64bit.raw | Bin 0 -> 8000 bytes .../multiexp_decay/covmat_5x5_64bit.raw | 1 + .../multiexp_decay/xdata_1000_64bit.raw | Bin 0 -> 8000 bytes .../multiexp_decay/ydata_1000_64bit.raw | Bin 0 -> 8000 bytes .../conf_1000_64bit.raw | Bin 0 -> 8000 bytes .../covmat_5x5_64bit.raw | Bin 0 -> 200 bytes .../xdata_1000_64bit.raw | Bin 0 -> 8000 bytes .../ydata_1000_64bit.raw | Bin 0 -> 8000 bytes tests/integration_tests/main.rs | 212 ++++++++-- 36 files changed, 1187 insertions(+), 486 deletions(-) create mode 100644 Todo.md create mode 100644 python/.gitignore create mode 100644 python/multiexp_decay.py create mode 100644 python/weighted_multiexp_decay.py create mode 100644 src/readme.rs create mode 100644 src/statistics/numeric_traits/mod.rs create mode 100644 test_assets/multiexp_decay/conf_1000_64bit.raw create mode 100644 test_assets/multiexp_decay/covmat_5x5_64bit.raw create mode 100644 test_assets/multiexp_decay/xdata_1000_64bit.raw create mode 100644 test_assets/multiexp_decay/ydata_1000_64bit.raw create mode 100644 test_assets/weighted_multiexp_decay/conf_1000_64bit.raw create mode 100644 test_assets/weighted_multiexp_decay/covmat_5x5_64bit.raw create mode 100644 test_assets/weighted_multiexp_decay/xdata_1000_64bit.raw create mode 100644 test_assets/weighted_multiexp_decay/ydata_1000_64bit.raw diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index d6c2d74..0533c22 100644 --- a/.github/workflows/build.yml +++ b/.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 diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index d11cdaa..2e3c179 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -4,7 +4,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main, dmz ] + branches: [ main, dev ] env: RUST_BACKTRACE: 1 diff --git a/.github/workflows/lints.yml b/.github/workflows/lints.yml index f7268cd..e7ba568 100644 --- a/.github/workflows/lints.yml +++ b/.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 diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 5ade0dd..a630c78 100644 --- a/.github/workflows/tests.yml +++ b/.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 diff --git a/CHANGES.md b/CHANGES.md index 3b54682..c62a73d 100644 --- a/CHANGES.md +++ b/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 diff --git a/Cargo.toml b/Cargo.toml index 5dc406c..ffe480a 100644 --- a/Cargo.toml +++ b/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" diff --git a/README.md b/README.md index ecaf681..83fcb4d 100644 --- a/README.md +++ b/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, tau : f64) -> DVector { - 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, tau : f64) + -> DVector { + t.map(|t|(-t/tau).exp()) } -fn exp_decay_dtau(tvec: &DVector,tau: f64) -> DVector { - 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,tau: f64) + -> DVector { + 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::::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 diff --git a/Todo.md b/Todo.md new file mode 100644 index 0000000..c7eb0cf --- /dev/null +++ b/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. diff --git a/benches/double_exponential_without_noise.rs b/benches/double_exponential_without_noise.rs index 4ff874e..71c3a40 100644 --- a/benches/double_exponential_without_noise.rs +++ b/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( true_parameters: DoubleExponentialParameters, mut model: Model, -) -> LevMarProblem +) -> LevMarProblem where Model: SeparableNonlinearModel, DefaultAllocator: nalgebra::allocator::Allocator, @@ -68,7 +69,7 @@ where .expect("Building valid problem should not panic") } -fn run_minimization(problem: LevMarProblem) -> [f64; 5] +fn run_minimization(problem: LevMarProblem) -> (DVector, DVector) where Model: SeparableNonlinearModel + 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 diff --git a/benches/multiple_right_hand_sides.rs b/benches/multiple_right_hand_sides.rs index cbfd43f..af06bd9 100644 --- a/benches/multiple_right_hand_sides.rs +++ b/benches/multiple_right_hand_sides.rs @@ -23,7 +23,7 @@ struct DoubleExponentialParameters { fn build_problem_mrhs( true_parameters: DoubleExponentialParameters, mut model: Model, -) -> LevMarProblem +) -> LevMarProblem where Model: SeparableNonlinearModel, { @@ -31,13 +31,13 @@ where // 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(problem: LevMarProblem) -> (DVector, DMatrix) +fn run_minimization_mrhs(problem: LevMarProblem) -> (DVector, DMatrix) where Model: SeparableNonlinearModel + 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, ) }); diff --git a/python/.gitignore b/python/.gitignore new file mode 100644 index 0000000..e51c42e --- /dev/null +++ b/python/.gitignore @@ -0,0 +1 @@ +*.raw diff --git a/python/multiexp_decay.py b/python/multiexp_decay.py new file mode 100644 index 0000000..3b87551 --- /dev/null +++ b/python/multiexp_decay.py @@ -0,0 +1,43 @@ +import lmfit as lm +import numpy as np + +# Defines the double exponential model with constant offset +# y = c1 exp(-t/tau1) + c2 exp(-t/tau2) + c3 +def multiexp_decay(t,c1,c2,c3,tau1,tau2): + return c1*np.exp(-t/tau1)+c2*np.exp(-t/tau2)+c3 + +# here are the true parameters of the model for +# generating synthetic data +c1 = 2.2 +c2 = 6.8 +c3 = 1.6 +tau1 = 2.4 +tau2 = 6.0 +ndata = 1000 + +# to have reproducible results +np.random.seed(0xdeadbeef) + +tdata = np.linspace(0,20,ndata) +ydata = multiexp_decay(tdata,c1,c2,c3,tau1,tau2)+np.random.normal(size=ndata,scale=0.01); + +model = lm.Model(multiexp_decay) + +# guess initial parameters and fit model +params = model.make_params(c1 = 1.0,c2 = 5, c3 = 0.3,tau1 = 1.,tau2 = 7.) +result = model.fit(ydata,params,t=tdata) + +print(result.fit_report()) +confidence_radius = result.eval_uncertainty(sigma=0.88,dscale = 0.0000001); + +# print(cb) +# print(result.ci_report()) + +# now write the data and the results to disk +# as raw float 64 values +tdata.astype(np.float64).tofile(f"xdata_{ndata}_64bit.raw"); +ydata.astype(np.float64).tofile(f"ydata_{ndata}_64bit.raw"); +confidence_radius.astype(np.float64).tofile(f"conf_{ndata}_64bit.raw"); + +cov = result.covar.astype(np.float64); +cov.tofile(f"covmat_{cov.shape[0]}x{cov.shape[1]}_64bit.raw"); diff --git a/python/weighted_multiexp_decay.py b/python/weighted_multiexp_decay.py new file mode 100644 index 0000000..d8f2579 --- /dev/null +++ b/python/weighted_multiexp_decay.py @@ -0,0 +1,43 @@ +import lmfit as lm +import numpy as np + +# Defines the double exponential model with constant offset +# y = c1 exp(-t/tau1) + c2 exp(-t/tau2) + c3 +def multiexp_decay(t,c1,c2,c3,tau1,tau2): + return c1*np.exp(-t/tau1)+c2*np.exp(-t/tau2)+c3 + +# here are the true parameters of the model for +# generating synthetic data +c1 = 2.2 +c2 = 6.8 +c3 = 1.6 +tau1 = 2.4 +tau2 = 6.0 +ndata = 1000 + +# to have reproducible results +np.random.seed(0xdeadbeef) + +tdata = np.linspace(0,20,ndata) +ydata = multiexp_decay(tdata,c1,c2,c3,tau1,tau2)+np.random.normal(size=ndata,scale=0.01); + +model = lm.Model(multiexp_decay) + +# guess initial parameters and fit model +params = model.make_params(c1 = 1.0,c2 = 5, c3 = 0.3,tau1 = 1.,tau2 = 7.) +result = model.fit(ydata,params,t=tdata, weights = 1/np.sqrt(ydata)) + +print(result.fit_report()) +confidence_radius = result.eval_uncertainty(sigma=0.88,dscale = 0.0000001); + +# print(cb) +# print(result.ci_report()) + +# now write the data and the results to disk +# as raw float 64 values +tdata.astype(np.float64).tofile(f"xdata_{ndata}_64bit.raw"); +ydata.astype(np.float64).tofile(f"ydata_{ndata}_64bit.raw"); +confidence_radius.astype(np.float64).tofile(f"conf_{ndata}_64bit.raw"); + +cov = result.covar.astype(np.float64); +cov.tofile(f"covmat_{cov.shape[0]}x{cov.shape[1]}_64bit.raw"); diff --git a/src/lib.rs b/src/lib.rs index be1c59b..7fc84b6 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -138,10 +138,11 @@ //! //! ```no_run //! # use varpro::model::*; -//! # let problem : varpro::solvers::levmar::LevMarProblem> = unimplemented!(); +//! # let problem : varpro::solvers::levmar::LevMarProblem,false> = unimplemented!(); //! use varpro::solvers::levmar::LevMarSolver; -//! let fit_result = LevMarSolver::new().fit(problem).unwrap(); +//! let fit_result = LevMarSolver::default().fit(problem).unwrap(); //! ``` +//! //! Finally, check the minimization report and, if successful, retrieve the nonlinear parameters `$\alpha$` //! using the [LevMarProblem::params](levenberg_marquardt::LeastSquaresProblem::params) and the linear //! coefficients `$\vec{c}$` using [LevMarProblem::linear_coefficients](crate::solvers::levmar::LevMarProblem::linear_coefficients) @@ -153,9 +154,9 @@ //! ```no_run //! # use varpro::model::SeparableModel; //! # use varpro::prelude::*; -//! # let problem : varpro::solvers::levmar::LevMarProblem> = unimplemented!(); +//! # let problem : varpro::solvers::levmar::LevMarProblem,false> = unimplemented!(); //! # use varpro::solvers::levmar::LevMarSolver; -//! # let fit_result = LevMarSolver::new().fit(problem).unwrap(); +//! # let fit_result = LevMarSolver::default().fit(problem).unwrap(); //! let alpha = fit_result.nonlinear_parameters(); //! let coeff = fit_result.linear_coefficients().unwrap(); //! ``` @@ -227,7 +228,7 @@ //! # let x = DVector::from_vec(vec![0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.]); //! # let y = DVector::from_vec(vec![1.0,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2,0.1,0.01]); //! -//! //1. create the model by giving only the nonlinear parameter names it depends on +//! // 1. create the model by giving only the nonlinear parameter names it depends on //! let model = SeparableModelBuilder::::new(&["tau1", "tau2"]) //! // add the first exponential decay and its partial derivative to the model //! // give all parameter names that the function depends on @@ -252,7 +253,7 @@ //! .build() //! .unwrap(); //! // 4. Solve using the fitting problem -//! let fit_result = LevMarSolver::new() +//! let fit_result = LevMarSolver::default() //! .fit(problem) //! .expect("fit must succeed"); //! // the nonlinear parameters after fitting @@ -332,7 +333,7 @@ //! .unwrap(); //! //! // fit the data -//! let fit_result = LevMarSolver::new() +//! let fit_result = LevMarSolver::default() //! .fit(problem) //! .expect("fit must succeed"); //! // the nonlinear parameters @@ -342,25 +343,71 @@ //! ``` //! # Global Fitting with Multiple Right Hand Sides //! -//! This is an application that lets varpro really shine. It solves the following -//! problem: +//! Instead of fitting a single data vector (i.e. a single _right hand side_), +//! this library can also solve a related, but slightly different problem. This +//! is the problem of global fitting for _multiple right hand sides_. The problem +//! statement is the following: +//! //! * We have not only one observation but a set `$\{\vec{y}_s\}$`, `$s=1,...,S$` of //! observations. //! * We want to fit the separable nonlinear function `$\vec{f}(\vec{\alpha},\vec{c})$` //! to all vectors of observations, but in such a way that the linear parameters -//! are allowed to be different for each vector, but the nonlinear parameters +//! are allowed to vary with `$s$`, but the nonlinear parameters //! are the same for the whole dataset. //! //! This is called global fitting with multiple right hand sides, //! because the nonlinear parameters are not //! allowed to change and are optimized for the complete dataset, whereas the linear -//! parameters are allowed to vary with each vector of observations. For a more -//! in depth explanation [see my article here](https://geo-ant.github.io/blog/2024/variable-projection-part-2-multiple-right-hand-sides/). +//! parameters are allowed to vary with each vector of observations. This is an application +//! where varpro really shines. Note that it is not the same as fitting the data +//! vectors independently. For a more in depth explanation +//! [see my article here](https://geo-ant.github.io/blog/2024/variable-projection-part-2-multiple-right-hand-sides/). //! //! To take advantage of global fitting we don't need to change anything about the -//! model or the building process. The only difference is that we now pass a matrix -//! as the observations rather than a single vector. The columns in the matrix -//! correspond to the observations. +//! model, we just have to make a slight modification to the way we build a problem. +//! +//! ## Example +//! ```no_run +//! # use varpro::prelude::*; +//! # use varpro::solvers::levmar::{LevMarProblemBuilder, LevMarSolver}; +//! # let model : varpro::model::SeparableModel:: = unimplemented!(); +//! # let Y : nalgebra::DMatrix:: = unimplemented!(); +//! # let w : nalgebra::DVector:: = unimplemented!(); +//! +//! // use the model as before but now invoke the `mrhs` +//! // constructor for the fitting problem +//! let problem = LevMarProblemBuilder::mrhs(model) +//! // the observations is a matrix where each column vector represents +//! // a single observation +//! .observations(Y) +//! .weights(w) +//! .build() +//! .unwrap(); +//! +//! // fit the data +//! let fit_result = LevMarSolver::default() +//! .fit(problem) +//! .expect("fit must succeed"); +//! +//! // the nonlinear parameters +//! // these parameters are fitted globally, meaning they +//! // are the same for each observation. Hence alpha is a single +//! // vector of parameters. +//! let alpha = fit_result.nonlinear_parameters(); +//! +//! // the linear coefficients +//! // those are a matrix for global fitting with multiple right hand sides. +//! // Each colum corresponds to the linear coefficients for the same column +//! // in the matrix of observations Y. +//! # #[allow(non_snake_case)] +//! let C = fit_result.linear_coefficients().unwrap(); +//! ``` +//! +//! The main difference to fitting problems with a single right hand side is that +//! the observations now become a matrix. Each column of this matrix is an +//! observation. Since the linear coefficients are allowed to vary, they now +//! also become a matrix instead of a single vector. Each column corresponds to +//! the best fit linear coefficients of the observations in the same matrix column. //! //! # 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) @@ -384,5 +431,9 @@ pub mod statistics; /// implemented in the nalgebra crate pub mod util; +/// a helper module for doctesting code in my readme +#[cfg(doctest)] +mod readme; + #[cfg(any(test, doctest))] pub mod test_helpers; diff --git a/src/model/builder/error.rs b/src/model/builder/error.rs index 3f706d8..a983d14 100644 --- a/src/model/builder/error.rs +++ b/src/model/builder/error.rs @@ -108,4 +108,8 @@ pub enum ModelBuildError { #[error("Missing initial guesses for model parameters")] /// missing initial guesses for parameters MissingInitialParameters, + + #[error("Illegal call to 'partial_deriv': a call to this function can only follow a call to 'function' or another call to 'partial_deriv'")] + /// an illegal call to the partial_deriv function + IllegalCallToPartialDeriv, } diff --git a/src/model/builder/mod.rs b/src/model/builder/mod.rs index 55afabe..f032cad 100644 --- a/src/model/builder/mod.rs +++ b/src/model/builder/mod.rs @@ -249,16 +249,32 @@ pub mod error; /// method. This method returns a valid model or an error variant doing a pretty good job of /// explaning why the model is invalid. #[must_use = "The builder should be transformed into a model using the build() method"] -pub struct SeparableModelBuilder +pub enum SeparableModelBuilder where ScalarType: Scalar, { - model_result: Result, ModelBuildError>, + #[doc(hidden)] + /// builder is in an error state + Error(ModelBuildError), + /// builder is in normal state (i.e. NOT in the process of building a function and not in an error state) + #[doc(hidden)] + Normal(UnfinishedModel), + /// builder is in the state of building a function + #[doc(hidden)] + FunctionBuilding { + /// the currently held model + model: UnfinishedModel, + /// the function we are currently building. It is added + /// to the model once a builder function is called that + /// indicates that building this function is over + function_builder: ModelBasisFunctionBuilder, + }, } /// a helper structure that represents an unfinished separable model #[derive(Default)] -struct UnfinishedModel { +#[doc(hidden)] +pub struct UnfinishedModel { /// the parameter names parameter_names: Vec, /// the base functions @@ -275,9 +291,17 @@ where ScalarType: Scalar, { fn from(err: ModelBuildError) -> Self { - Self { - model_result: Err(err), - } + Self::Error(err) + } +} + +impl From> for SeparableModelBuilder +where + ScalarType: Scalar, +{ + #[inline] + fn from(model: UnfinishedModel) -> Self { + Self::Normal(model) } } @@ -287,8 +311,12 @@ impl From, ModelBuildError>> where ScalarType: Scalar, { + #[inline] fn from(model_result: Result, ModelBuildError>) -> Self { - Self { model_result } + match model_result { + Ok(model) => Self::Normal(model), + Err(err) => Self::Error(err), + } } } @@ -318,17 +346,15 @@ where .collect(); if let Err(parameter_error) = check_parameter_names(¶meter_names) { - Self { - model_result: Err(parameter_error), - } + Self::Error(parameter_error) } else { - let model_result = Ok(UnfinishedModel { + let model_result = UnfinishedModel { parameter_names, basefunctions: Vec::new(), x_vector: None, initial_parameters: None, - }); - Self { model_result } + }; + Self::from(model_result) } } @@ -337,16 +363,30 @@ where /// # Usage /// For usage see the documentation of the [SeparableModelBuilder](crate::model::builder::SeparableModelBuilder) /// struct documentation. - pub fn invariant_function(mut self, function: F) -> Self + pub fn invariant_function(self, function: F) -> Self where F: Fn(&DVector) -> DVector + 'static, { - if let Ok(model) = self.model_result.as_mut() { - model - .basefunctions - .push(ModelBasisFunction::parameter_independent(function)); + match self { + SeparableModelBuilder::Error(err) => Self::from(err), + SeparableModelBuilder::Normal(mut model) => { + model + .basefunctions + .push(ModelBasisFunction::parameter_independent(function)); + Self::from(model) + } + SeparableModelBuilder::FunctionBuilding { + model, + function_builder, + } => + // this will not be an infinite recursion because it goes to either the normal + // or error state. This template can be used pretty much everywhere, because + // we just extend the model (finalize the function builder) and then invoke the + // called function again + { + Self::from(extend_model(model, function_builder)).invariant_function(function) + } } - self } /// Add a function `$\vec{f}(\vec{x},\alpha_1,...,\alpha_n)$` to the model that depends on the @@ -358,41 +398,104 @@ where self, function_params: StrCollection, function: F, - ) -> SeparableModelBuilderProxyWithDerivatives + ) -> Self where F: BasisFunction + 'static, StrCollection: IntoIterator, StrCollection::Item: AsRef, { - SeparableModelBuilderProxyWithDerivatives::new(self.model_result, function_params, function) + match self { + SeparableModelBuilder::Error(err) => Self::from(err), + SeparableModelBuilder::Normal(model) => { + let function_builder = ModelBasisFunctionBuilder::new( + model.parameter_names.clone(), + function_params, + function, + ); + Self::FunctionBuilding { + model, + function_builder, + } + } + SeparableModelBuilder::FunctionBuilding { + model, + function_builder, + } => Self::from(extend_model(model, function_builder)) + .function(function_params, function), + } + } + + /// Add a partial derivative to a function, see the example in the documentation + /// to this structure. + /// A call to this function must only occur when it follows a call to + /// `function` or another call to `partial_derivative`. Other cases will + /// lead to errors when building the model. + pub fn partial_deriv, F, ArgList>( + self, + parameter: StrType, + derivative: F, + ) -> Self + where + F: BasisFunction + 'static, + { + match self { + SeparableModelBuilder::Error(err) => Self::from(err), + SeparableModelBuilder::Normal(_model) => { + // only when we are in the process of building a function, may + // we call this function + Self::from(Err(ModelBuildError::IllegalCallToPartialDeriv)) + } + SeparableModelBuilder::FunctionBuilding { + model, + function_builder, + } => Self::FunctionBuilding { + model, + function_builder: function_builder.partial_deriv(parameter.as_ref(), derivative), + }, + } } /// Set the independent variable `$x$` which will be used when evaluating the model. /// Also see the struct documentation of [SeparableModelBuilder](crate::model::builder::SeparableModelBuilder) /// for information on how to use this method. - pub fn independent_variable(mut self, x: DVector) -> Self { - if let Ok(model) = self.model_result.as_mut() { - model.x_vector = Some(x); + pub fn independent_variable(self, x: DVector) -> Self { + match self { + SeparableModelBuilder::Error(err) => Self::from(err), + SeparableModelBuilder::Normal(mut model) => { + model.x_vector = Some(x); + Self::from(model) + } + SeparableModelBuilder::FunctionBuilding { + model, + function_builder, + } => Self::from(extend_model(model, function_builder)).independent_variable(x), } - self } /// Set the initial values for the model parameters `$\vec{\alpha}$`. /// Also see the struct documentation of [SeparableModelBuilder](crate::model::builder::SeparableModelBuilder) /// for information on how to use this method. - pub fn initial_parameters(mut self, initial_parameters: Vec) -> Self { - if let Ok(model) = self.model_result.as_mut() { - let expected = model.parameter_names.len(); - if expected != initial_parameters.len() { - self.model_result = Err(ModelBuildError::IncorrectParameterCount { - expected, - actual: initial_parameters.len(), - }); - } else { - model.initial_parameters = Some(initial_parameters); + pub fn initial_parameters(self, initial_parameters: Vec) -> Self { + match self { + SeparableModelBuilder::Error(err) => Self::from(err), + SeparableModelBuilder::Normal(mut model) => { + let expected = model.parameter_names.len(); + if expected != initial_parameters.len() { + Self::from(Err(ModelBuildError::IncorrectParameterCount { + expected, + actual: initial_parameters.len(), + })) + } else { + model.initial_parameters = Some(initial_parameters); + Self::from(model) + } } + SeparableModelBuilder::FunctionBuilding { + model, + function_builder, + } => Self::from(extend_model(model, function_builder)) + .initial_parameters(initial_parameters), } - self } /// Build a separable model from the contents of this builder. @@ -411,7 +514,14 @@ where /// where provided during the builder stage. That means the first basis functions gets index `0` in /// the model, the second gets index `1` and so on. pub fn build(self) -> Result, ModelBuildError> { - self.model_result.and_then(TryInto::try_into) + match self { + SeparableModelBuilder::Error(err) => Err(err), + SeparableModelBuilder::Normal(model) => model.try_into(), + SeparableModelBuilder::FunctionBuilding { + model, + function_builder, + } => extend_model(model, function_builder).and_then(TryInto::try_into), + } } } @@ -459,211 +569,15 @@ impl TryInto> for UnfinishedModel } } } -/// helper struct that contains a seperable model as well as a model function builder -/// used inside the SeparableModelBuilderProxyWithDerivatives. -struct ModelAndModelBasisFunctionBuilderPair -where - ScalarType: Scalar, -{ - model: UnfinishedModel, - builder: ModelBasisFunctionBuilder, -} - -impl ModelAndModelBasisFunctionBuilderPair -where - ScalarType: Scalar, -{ - fn new( - model: UnfinishedModel, - builder: ModelBasisFunctionBuilder, - ) -> Self { - Self { model, builder } - } -} - -/// This is just a proxy that does need to be used directly. For constructing a model from -/// a builder see the documentation for [SeparableModelBuilder](self::SeparableModelBuilder). -/// **Sidenote** This structure will hopefully be made more elegant using some metaprogramming techniques -/// in the future. Right now this exists to make sure that partial derivatives cannot accidentally -/// be added to invariant functions. The compiler simply will forbid it. In future the library aims -/// to make more such checks at compile time and reduce the need for runtime errors caused by invalid model -/// construction. -#[must_use = "This is meant as a transient expression proxy. Use build() to build a model."] -pub struct SeparableModelBuilderProxyWithDerivatives -where - ScalarType: Scalar, -{ - current_result: Result, ModelBuildError>, -} - -impl From for SeparableModelBuilderProxyWithDerivatives -where - ScalarType: Scalar, -{ - fn from(err: ModelBuildError) -> Self { - Self { - current_result: Err(err), - } - } -} - -impl SeparableModelBuilderProxyWithDerivatives -where - ScalarType: Scalar, -{ - /// Construct an instance. This is invoked when adding a function that depends on model - /// parameters to the [SeparableModelBuilder](self::SeparableModelBuilder) - /// # Arguments - /// * `model_result`: the current model or an error - /// * `function_parameters`: the given list of nonlinear function parameters - /// * `function`: the function - fn new( - model_result: Result, ModelBuildError>, - function_parameters: StrCollection, - function: F, - ) -> Self - where - F: BasisFunction + 'static, - StrCollection: IntoIterator, - StrCollection::Item: AsRef, - { - match model_result { - Ok(model) => { - let model_parameters = model.parameter_names.clone(); - Self { - current_result: Ok(ModelAndModelBasisFunctionBuilderPair::new( - model, - ModelBasisFunctionBuilder::new( - model_parameters, - function_parameters, - function, - ), - )), - } - } - Err(err) => Self { - current_result: Err(err), - }, - } - } - - /// Add a partial derivative to a function. This function is documented as part of the - /// [SeparableModelBuilder](self::SeparableModelBuilder) documentation. - /// *Info*: the reason it appears in this proxy only is to make sure that partial derivatives - /// can only be added to functions that depend on model parameters. - pub fn partial_deriv, F, ArgList>( - self, - parameter: StrType, - derivative: F, - ) -> Self - where - F: BasisFunction + 'static, - { - match self.current_result { - Ok(result) => Self { - current_result: Ok(ModelAndModelBasisFunctionBuilderPair { - model: result.model, - builder: result.builder.partial_deriv(parameter.as_ref(), derivative), - }), - }, - Err(err) => Self::from(err), - } - } - - /// Add a function `\vec{f}(\vec{x})` that does not depend on the model parameters. This is documented - /// as part of the [SeparableModelBuilder](self::SeparableModelBuilder) documentation. - pub fn invariant_function(self, function: F) -> SeparableModelBuilder - where - F: Fn(&DVector) -> DVector + 'static, - { - match self.current_result { - Ok(pair) => { - let model_result = extend_model(pair.model, pair.builder); - SeparableModelBuilder::from(model_result).invariant_function(function) - } - Err(err) => SeparableModelBuilder::from(err), - } - } - - /// Add a function `\vec{f}(\vec{x},\alpha_j,...,\alpha_k)` that depends on a subset of the - /// model parameters. This functionality is documented as part of the [SeparableModelBuilder](self::SeparableModelBuilder) - /// documentation. - pub fn function( - self, - function_params: StrCollection, - function: F, - ) -> SeparableModelBuilderProxyWithDerivatives - where - F: BasisFunction + 'static, - StrCollection: IntoIterator, - StrCollection::Item: AsRef, - { - match self.current_result { - Ok(pair) => { - let model_result = extend_model(pair.model, pair.builder); - Self::new(model_result, function_params, function) - } - Err(err) => SeparableModelBuilderProxyWithDerivatives::from(err), - } - } - - /// Set the independent variable `\vec{x}` - /// - /// # Usage - /// For usage see the documentation of the [SeparableModelBuilder](crate::model::builder::SeparableModelBuilder) - /// struct documentation. - pub fn independent_variable(self, x: DVector) -> SeparableModelBuilder { - match self.current_result { - Ok(pair) => { - let model_result = extend_model(pair.model, pair.builder); - SeparableModelBuilder::from(model_result).independent_variable(x) - } - Err(err) => SeparableModelBuilder::from(err), - } - } - - /// Set the initial value for the nonlinear parameters `\vec{\alpha}` - /// - /// # Usage - /// For usage see the documentation of the [SeparableModelBuilder](crate::model::builder::SeparableModelBuilder) - /// struct documentation. - pub fn initial_parameters( - self, - initial_parameters: Vec, - ) -> SeparableModelBuilder { - match self.current_result { - Ok(pair) => { - let model_result = extend_model(pair.model, pair.builder); - SeparableModelBuilder::from(model_result).initial_parameters(initial_parameters) - } - Err(err) => SeparableModelBuilder::from(err), - } - } - - /// Finalized the building process and build a separable model. - /// This functionality is documented as part of the [SeparableModelBuilder](self::SeparableModelBuilder) - /// documentation. - pub fn build(self) -> Result, ModelBuildError> { - // this method converts the internal results into a separable model and uses its - // facilities to check for completion and the like - match self.current_result { - Ok(pair) => { - let model_result = extend_model(pair.model, pair.builder); - SeparableModelBuilder::from(model_result).build() - } - Err(err) => SeparableModelBuilder::from(err).build(), - } - } -} /// try and extend a model with the given function in the builder /// if building the function in the builder fails, an error is returned, /// otherwise the extended model is returned fn extend_model( mut model: UnfinishedModel, - builder: ModelBasisFunctionBuilder, + function_builder: ModelBasisFunctionBuilder, ) -> Result, ModelBuildError> { - let function = builder.build()?; + let function = function_builder.build()?; model.basefunctions.push(function); Ok(model) } diff --git a/src/model/builder/modelfunction_builder/mod.rs b/src/model/builder/modelfunction_builder/mod.rs index e0f65c6..0c92c74 100644 --- a/src/model/builder/modelfunction_builder/mod.rs +++ b/src/model/builder/modelfunction_builder/mod.rs @@ -19,7 +19,7 @@ use crate::model::model_basis_function::ModelBasisFunction; /// (all other derivatives are zero, because the function does not depends on other params) /// #[doc(hidden)] -pub(crate) struct ModelBasisFunctionBuilder +pub struct ModelBasisFunctionBuilder where ScalarType: Scalar, { diff --git a/src/model/builder/test.rs b/src/model/builder/test.rs index 9f98c15..a492efb 100644 --- a/src/model/builder/test.rs +++ b/src/model/builder/test.rs @@ -258,6 +258,32 @@ fn test_model_builder_fails_when_initial_parameters_have_incorrect_parameter_cou assert_matches!(result, Err(ModelBuildError::IncorrectParameterCount { .. })); } +#[test] +// this is just to make sure that this code compiles when used like that. +fn model_builder_can_be_used_programmatically_in_a_loop() { + let param_names = &["param1", "param2", "param3"]; + let mut builder = SeparableModelBuilder::::new(param_names) + .independent_variable(nalgebra::dvector![1., 2., 3., 4.]) + .initial_parameters(vec![1., 2., 3.]); + + fn dummy_function(_x: &DVector, _param: f64) -> DVector { + todo!() + } + fn dummy_deriv(_x: &DVector, _param: f64) -> DVector { + todo!() + } + + for param in param_names { + builder = builder + .function(&[*param], dummy_function) + .partial_deriv(*param, dummy_deriv); + } + + // we just need this to build the model successfully + // all evaluations would panic + let _ = builder.build().unwrap(); +} + /// a function that calculates exp( -(t-t0)/tau)) for every location t /// **ATTENTION** for this kind of exponential function the shift will /// just be a linear multiplier exp(t0/tau), so it might not a good idea to include it in the fitting diff --git a/src/model/model_basis_function.rs b/src/model/model_basis_function.rs index 61d89e6..483f839 100644 --- a/src/model/model_basis_function.rs +++ b/src/model/model_basis_function.rs @@ -13,7 +13,8 @@ type BaseFuncType = /// An internal type that is used to store basefunctions whose interface has been wrapped in /// such a way that they can accept the location and the *complete model parameters as arguments*. -pub(crate) struct ModelBasisFunction +#[doc(hidden)] +pub struct ModelBasisFunction where ScalarType: Scalar, { diff --git a/src/readme.rs b/src/readme.rs new file mode 100644 index 0000000..7c7f3fb --- /dev/null +++ b/src/readme.rs @@ -0,0 +1,5 @@ +// this is a way of running doctests on my readme so I know the code is still valid +// https://github.com/rust-lang/cargo/issues/383. +// That means errors that occur here need to be fixed in the readme. +#[doc = include_str!("../README.md")] +struct _Readme {} diff --git a/src/solvers/levmar/builder.rs b/src/solvers/levmar/builder.rs index c24ad9e..1758bab 100644 --- a/src/solvers/levmar/builder.rs +++ b/src/solvers/levmar/builder.rs @@ -2,7 +2,7 @@ use crate::prelude::*; use crate::solvers::levmar::LevMarProblem; use crate::util::Weights; use levenberg_marquardt::LeastSquaresProblem; -use nalgebra::{ComplexField, DMatrix, DVector, Dyn, OVector, Scalar}; +use nalgebra::{ComplexField, DMatrix, Dyn, OMatrix, OVector, Scalar}; use num_traits::{Float, Zero}; use std::ops::Mul; use thiserror::Error as ThisError; @@ -11,7 +11,7 @@ use thiserror::Error as ThisError; #[derive(Debug, Clone, ThisError, PartialEq, Eq)] pub enum LevMarBuilderError { /// the data for y variable was not given to the builder - #[error("Data for vector y not provided.")] + #[error("Right hand side(s) not provided")] YDataMissing, /// x and y vector have different lengths @@ -38,31 +38,6 @@ pub enum LevMarBuilderError { InvalidLengthOfWeights, } -/// a type describing the observation `$\vec{y}$` or `$\boldsymbol{Y}$` that we want to fit -/// with a model. The observation can either be create from a single column vector, corresponding -/// to the common case of fitting a single right hand side with a model. It can -/// also be create from a matrix, which corresponds to a global fitting problem -/// with multiple right hand sides. The observations in this case are a matrix -/// where each _column_ represents an observation. -pub struct Observation { - data: DMatrix, -} - -// construct the obseration from a single vector or matrix -impl From> for Observation { - fn from(data: DMatrix) -> Self { - Self { data } - } -} - -impl From> for Observation { - fn from(data: DVector) -> Self { - Self { - data: DMatrix::from_column_slice(data.nrows(), 1, data.as_slice()), - } - } -} - /// A builder structure to create a [LevMarProblem](super::LevMarProblem), which can be used for /// fitting a separable model to data. /// # Example @@ -81,15 +56,23 @@ impl From> for Observation { /// .unwrap(); /// # } /// ``` +/// /// # Building a Model -/// A new builder is constructed with the [new](LevMarProblemBuilder::new) method. It must be filled with +/// +/// A new builder is constructed with the [new](LevMarProblemBuilder::new) constructor. It must be filled with /// content using the methods described in the following. After all mandatory fields have been filled, /// the [build](LevMarProblemBuilder::build) method can be called. This returns a [Result](std::result::Result) /// type that contains the finished model iff all mandatory fields have been set with valid values. Otherwise /// it contains an error variant. +/// +/// ## Multiple Right Hand Sides +/// +/// We can also construct a problem with multiple right hand sides, using the +/// [mrhs](LevMarProblemBuilder::mrhs) constructor, see [LevMarProblem] for +/// additional details. #[derive(Clone)] #[allow(non_snake_case)] -pub struct LevMarProblemBuilder +pub struct LevMarProblemBuilder where Model::ScalarType: Scalar + ComplexField + Copy, ::RealField: Float, @@ -111,7 +94,7 @@ where weights: Weights, } -impl LevMarProblemBuilder +impl LevMarProblemBuilder where Model::ScalarType: Scalar + ComplexField + Zero + Copy, ::RealField: Float, @@ -119,7 +102,9 @@ where Mul + Float, Model: SeparableNonlinearModel, { - /// Create a new builder based on the given model + /// Create a new builder based on the given model for a problem + /// with a **single right hand side**. This is the standard use case, + /// where the data is a vector that is fitted by the model. pub fn new(model: Model) -> Self { Self { Y: None, @@ -128,7 +113,48 @@ where weights: Weights::default(), } } + /// **Mandatory**: Set the data which we want to fit: since this + /// is called on a model builder for problems with **single right hand sides**, + /// this is a column vector `$\vec{y}=\vec{y}(\vec{x})$` containing the + /// values we want to fit with the model. + /// + /// The length of `$\vec{y}` and the output dimension of the model must be + /// the same. + pub fn observations(self, observed: OVector) -> Self { + let nrows = observed.nrows(); + Self { + Y: Some(observed.reshape_generic(Dyn(nrows), Dyn(1))), + ..self + } + } +} +impl LevMarProblemBuilder +where + Model::ScalarType: Scalar + ComplexField + Zero + Copy, + ::RealField: Float, + <::ScalarType as ComplexField>::RealField: + Mul + Float, + Model: SeparableNonlinearModel, +{ + /// Create a new builder based on the given model + /// for a problem with **multiple right hand sides** and perform a global + /// fit. + /// + /// That means the observations are expected to be a matrix, where the + /// columns correspond to the individual observations. The nonlinear + /// parameters will be optimized across all the observations, but the + /// linear coefficients are calculated for each observation individually. + /// Hence, they also become a matrix where the columns correspond to the + /// linear coefficients of the observation in the same column. + pub fn mrhs(model: Model) -> Self { + Self { + Y: None, + separable_model: model, + epsilon: None, + weights: Weights::default(), + } + } /// **Mandatory**: Set the data which we want to fit: This is either a single vector /// `$\vec{y}=\vec{y}(\vec{x})$` or a matrix `$\boldsymbol{Y}$` of multiple /// vectors. In the former case this corresponds to fitting a single right hand side, @@ -136,13 +162,22 @@ where /// multiple right hand sides. /// The length of `$\vec{x}$` and the number of _rows_ in the data must /// be the same. - pub fn observations(self, observed: impl Into>) -> Self { + pub fn observations(self, observed: OMatrix) -> Self { Self { - Y: Some(observed.into().data), + Y: Some(observed), ..self } } +} +impl LevMarProblemBuilder +where + Model::ScalarType: Scalar + ComplexField + Zero + Copy, + ::RealField: Float, + <::ScalarType as ComplexField>::RealField: + Mul + Float, + Model: SeparableNonlinearModel, +{ /// **Optional** This value is relevant for the solver, because it uses singular value decomposition /// internally. This method sets a value `\epsilon` for which smaller (i.e. absolute - wise) singular /// values are considered zero. In essence this gives a truncation of the SVD. This might be @@ -185,7 +220,7 @@ where /// If all prerequisites are fulfilled, returns a [LevMarProblem](super::LevMarProblem) with the given /// content and the parameters set to the initial guess. Otherwise returns an error variant. #[allow(non_snake_case)] - pub fn build(self) -> Result, LevMarBuilderError> { + pub fn build(self) -> Result, LevMarBuilderError> { // and assign the defaults to the values we don't have let Y = self.Y.ok_or(LevMarBuilderError::YDataMissing)?; let model = self.separable_model; diff --git a/src/solvers/levmar/mod.rs b/src/solvers/levmar/mod.rs index 49db639..21dff4d 100644 --- a/src/solvers/levmar/mod.rs +++ b/src/solvers/levmar/mod.rs @@ -3,8 +3,8 @@ use crate::statistics::FitStatistics; use levenberg_marquardt::{LeastSquaresProblem, MinimizationReport}; use nalgebra::storage::Owned; use nalgebra::{ - ComplexField, DMatrix, DefaultAllocator, Dim, DimMin, Dyn, Matrix, OVector, RawStorageMut, - RealField, Scalar, UninitMatrix, Vector, SVD, U1, + ComplexField, DMatrix, DefaultAllocator, Dim, DimMin, Dyn, Matrix, MatrixView, OMatrix, + OVector, RawStorageMut, RealField, Scalar, UninitMatrix, Vector, VectorView, SVD, U1, }; mod builder; @@ -23,7 +23,7 @@ use std::ops::Mul; /// [`LevenbergMarquardt`](https://docs.rs/levenberg-marquardt/latest/levenberg_marquardt/struct.LevenbergMarquardt.html) /// solver from the `levenberg_marquardt` crate. The core benefit of this /// wrapper is that we can also use it to calculate statistics. -pub struct LevMarSolver +pub struct LevMarSolver where Model: SeparableNonlinearModel, { @@ -33,14 +33,14 @@ where /// A helper type that contains the fitting problem after the /// minimization, as well as a report and some convenience functions #[derive(Debug)] -pub struct FitResult +pub struct FitResult where Model: SeparableNonlinearModel, Model::ScalarType: RealField + Scalar + Float, { /// the final state of the fitting problem after the /// minimization finished (regardless of whether fitting was successful or not) - pub problem: LevMarProblem, + pub problem: LevMarProblem, /// the minimization report of the underlying solver. /// It contains information about the minimization process @@ -49,7 +49,68 @@ where pub minimization_report: MinimizationReport, } -impl FitResult +impl FitResult +// take trait bounds from above: +where + Model: SeparableNonlinearModel, + Model::ScalarType: RealField + Scalar + Float, +{ + /// **Note** This implementation is for fitting problems with multiple right hand sides. + /// + /// Convenience function to get the linear coefficients after the fit has + /// finished. Will return None if there was an error during fitting. + /// + /// The coefficients vectors for the individual + /// members of the datasets are the colums of the returned matrix. That means + /// one coefficient vector for each right hand side. + pub fn linear_coefficients(&self) -> Option> { + Some(self.problem.cached.as_ref()?.linear_coefficients.as_view()) + } + + /// **Note** This implementation is for fitting problems with a single right hand side. + /// + /// Calculate the values of the model at the best fit parameters. + /// Returns None if there was an error during fitting. + /// Since this is for single right hand sides, the output is + /// a column vector. + pub fn best_fit(&self) -> Option> { + let coeff = self.linear_coefficients()?; + let eval = self.problem.model().eval().ok()?; + Some(eval * coeff) + } +} + +impl FitResult +// take trait bounds from above: +where + Model: SeparableNonlinearModel, + Model::ScalarType: RealField + Scalar + Float, +{ + /// **Note** This implementation is for fitting problems with a single right hand side. + /// + /// Convenience function to get the linear coefficients after the fit has + /// finished. Will return None if there was an error during fitting. + /// The coefficients are given as a single vector. + pub fn linear_coefficients(&self) -> Option> { + let coeff = &self.problem.cached.as_ref()?.linear_coefficients; + debug_assert_eq!(coeff.ncols(),1, + "Coefficient matrix must have exactly one colum for problem with single right hand side. This indicates a programming error inside this library!"); + Some(self.problem.cached.as_ref()?.linear_coefficients.column(0)) + } + /// **Note** This implementation is for fitting problems with multiple right hand sides + /// + /// Calculate the values of the model at the best fit parameters. + /// Returns None if there was an error during fitting. + /// Since this is for single right hand sides, the output is + /// a matrix containing composed of column vectors for the right hand sides. + pub fn best_fit(&self) -> Option> { + let coeff = self.linear_coefficients()?; + let eval = self.problem.model().eval().ok()?; + Some(eval * coeff) + } +} + +impl FitResult // take trait bounds from above: where Model: SeparableNonlinearModel, @@ -57,7 +118,7 @@ where { /// internal helper for constructing an instance fn new( - problem: LevMarProblem, + problem: LevMarProblem, minimization_report: MinimizationReport, ) -> Self { Self { @@ -72,17 +133,6 @@ where self.problem.model().params() } - /// convenience function to get the linear coefficients after the fit has - /// finished - /// - /// in case of multiple datasets, the coefficients vectors for the individual - /// members of the datasets are the colums of the matrix. If the problem - /// only has a single right hand side, then the matrix will have one column - /// only. - pub fn linear_coefficients(&self) -> Option> { - self.problem.linear_coefficients().cloned() - } - /// whether the fit was deemeed successful. The fit might still be not /// be optimal for numerical reasons, but the minimization process /// terminated successfully. @@ -91,7 +141,7 @@ where } } -impl LevMarSolver +impl LevMarSolver where Model: SeparableNonlinearModel, { @@ -118,11 +168,14 @@ where #[deprecated(since = "0.7.0", note = "use the fit(...) method instead")] pub fn minimize( &self, - problem: LevMarProblem, - ) -> (LevMarProblem, MinimizationReport) + problem: LevMarProblem, + ) -> ( + LevMarProblem, + MinimizationReport, + ) where Model: SeparableNonlinearModel, - LevMarProblem: LeastSquaresProblem, + LevMarProblem: LeastSquaresProblem, Model::ScalarType: Scalar + ComplexField + Copy + RealField + Float + FromPrimitive, { self.solver.minimize(problem) @@ -141,10 +194,13 @@ where /// On failure (when the minimization was not deemeed successful), returns /// an error with the same information as in the success case. #[allow(clippy::result_large_err)] - pub fn fit(&self, problem: LevMarProblem) -> Result, FitResult> + pub fn fit( + &self, + problem: LevMarProblem, + ) -> Result, FitResult> where Model: SeparableNonlinearModel, - LevMarProblem: LeastSquaresProblem, + LevMarProblem: LeastSquaresProblem, Model::ScalarType: Scalar + ComplexField + RealField + Float + FromPrimitive, { #[allow(deprecated)] @@ -156,7 +212,12 @@ where Err(result) } } +} +impl LevMarSolver +where + Model: SeparableNonlinearModel, +{ /// Same as the [LevMarProblem::fit] function but also calculates some additional /// statistical information about the fit, if the fit was successful. /// @@ -164,14 +225,20 @@ where /// /// See also the [LevMarProblem::fit] function, but on success also returns statistical /// information about the fit. + /// + /// ## Problems with Multiple Right Hand Sides + /// + /// **Note** For now, fitting with statistics is only supported for problems + /// with a single right hand side. If this function is invoked on a problem + /// with multiple right hand sides, an error is returned. #[allow(clippy::result_large_err)] pub fn fit_with_statistics( &self, - problem: LevMarProblem, - ) -> Result<(FitResult, FitStatistics), FitResult> + problem: LevMarProblem, + ) -> Result<(FitResult, FitStatistics), FitResult> where Model: SeparableNonlinearModel, - LevMarProblem: LeastSquaresProblem, + LevMarProblem: LeastSquaresProblem, Model::ScalarType: Scalar + ComplexField + RealField + Float, { let FitResult { @@ -198,8 +265,7 @@ where } } } - -impl Default for LevMarSolver +impl Default for LevMarSolver where Model: SeparableNonlinearModel, Model::ScalarType: RealField + Float, @@ -247,24 +313,41 @@ where /// This is a the problem of fitting the separable model to data in a form that the /// [levenberg_marquardt](https://crates.io/crates/levenberg-marquardt) crate can use it to /// perform the least squares fit. +/// /// # Construction +/// /// Use the [LevMarProblemBuilder](self::builder::LevMarProblemBuilder) to create an instance of a /// levmar problem. +/// /// # Usage +/// /// After obtaining an instance of `LevMarProblem` we can pass it to the [LevenbergMarquardt](levenberg_marquardt::LevenbergMarquardt) /// structure of the levenberg_marquardt crate for minimization. Refer to the documentation of the /// [levenberg_marquardt](https://crates.io/crates/levenberg-marquardt) for an overview. A usage example /// is provided in this crate documentation as well. The [LevenbergMarquardt](levenberg_marquardt::LevenbergMarquardt) /// solver is reexported by this module as [LevMarSolver](self::LevMarSolver) for naming consistency. +/// +/// # `MRHS`: Multiple Right Hand Sides +/// +/// The problem generic on the boolean `MRHS` which indicates whether the +/// problem fits a single (`MRHS == false`) or multiple (`MRHS == true`) right +/// hand sides. This is decided during the building process. The underlying +/// math does not change, but the interface changes to use vectors for coefficients +/// and data in case of a single right hand side. For multiple right hand sides, +/// the coefficients and the data are matrices corresponding to columns of +/// coefficient vectors and data vectors respectively. #[derive(Clone)] #[allow(non_snake_case)] -pub struct LevMarProblem +pub struct LevMarProblem where Model: SeparableNonlinearModel, Model::ScalarType: Scalar + ComplexField + Copy, { - /// the *weighted* data vector to which to fit the model `$\vec{y}_w$` - /// **Attention** the data vector is weighted with the weights if some weights + /// the *weighted* data matrix to which to fit the model `$\boldsymbol{Y}_w$`. + /// It is a matrix so it can accomodate multiple right hand sides. If + /// the problem has only a single right hand side (MRHS = false), this is just + /// a matrix with one column. The underlying math does not change in either case. + /// **Attention** the data matrix is weighted with the weights if some weights /// where provided (otherwise it is unweighted) Y_w: DMatrix, /// a reference to the separable model we are trying to fit to the data @@ -282,7 +365,7 @@ where cached: Option>, } -impl std::fmt::Debug for LevMarProblem +impl std::fmt::Debug for LevMarProblem where Model: SeparableNonlinearModel, Model::ScalarType: Scalar + ComplexField + Copy, @@ -298,7 +381,40 @@ where } } -impl LevMarProblem +impl LevMarProblem +where + Model: SeparableNonlinearModel, + Model::ScalarType: Scalar + ComplexField + Copy, +{ + /// Get the linear coefficients for the current problem. After a successful pass of the solver, + /// this contains a value with the best fitting linear coefficients + /// + /// # Returns + /// + /// Either the current best estimate coefficients or None, if none were calculated or the solver + /// encountered an error. After the solver finished, this is the least squares best estimate + /// for the linear coefficients of the base functions. + /// + /// Since this method is for fitting a single right hand side, the coefficients + /// are a single column vector. + pub fn linear_coefficients(&self) -> Option> { + self.cached + .as_ref() + .map(|cache| cache.linear_coefficients.as_view()) + } + + /// the weighted data matrix`$\boldsymbol{Y}_w$` to which to fit the model. Note + /// that the weights are already applied to the data matrix and this + /// is not the original data vector. + /// + /// This method is for fitting a multiple right hand sides, hence the data + /// matrix is a matrix that contains the right hand sides as columns. + pub fn weighted_data(&self) -> MatrixView { + self.Y_w.as_view() + } +} + +impl LevMarProblem where Model: SeparableNonlinearModel, Model::ScalarType: Scalar + ComplexField + Copy, @@ -313,10 +429,33 @@ where /// The linear coefficients are column vectors that are ordered /// into a matrix, where the column at index $$s$$ are the best linear /// coefficients for the member at index $$s$$ of the dataset. - pub fn linear_coefficients(&self) -> Option<&DMatrix> { - self.cached.as_ref().map(|cache| &cache.linear_coefficients) + pub fn linear_coefficients(&self) -> Option> { + self.cached + .as_ref() + .map(|cache|{ + debug_assert_eq!(cache.linear_coefficients.ncols(),1, + "coefficient matrix must have exactly one column for single right hand side. This indicates a programming error in the library."); + cache.linear_coefficients.as_view() + }) } + /// the weighted data vector `$\vec{y}_w$` to which to fit the model. Note + /// that the weights are already applied to the data vector and this + /// is not the original data vector. + /// + /// This method is for fitting a single right hand side, hence the data + /// is a single column vector. + pub fn weighted_data(&self) -> VectorView { + debug_assert_eq!(self.Y_w.ncols(),1, + "data matrix must have exactly one column for single right hand side. This indicates a programming error in the library."); + self.Y_w.as_view() + } +} +impl LevMarProblem +where + Model: SeparableNonlinearModel, + Model::ScalarType: Scalar + ComplexField + Copy, +{ /// access the contained model immutably pub fn model(&self) -> &Model { &self.model @@ -326,16 +465,10 @@ where pub fn weights(&self) -> &Weights { &self.weights } - - /// the weighted data vector `$\vec{y}_w$` to which to fit the model. Note - /// that the weights are already applied to the data vector and this - /// is not the original data vector. - pub fn weighted_data(&self) -> &DMatrix { - &self.Y_w - } } -impl LeastSquaresProblem for LevMarProblem +impl LeastSquaresProblem + for LevMarProblem where Model::ScalarType: Scalar + ComplexField + Copy, <::ScalarType as ComplexField>::RealField: diff --git a/src/solvers/levmar/test.rs b/src/solvers/levmar/test.rs index 09a2679..fb38cb6 100644 --- a/src/solvers/levmar/test.rs +++ b/src/solvers/levmar/test.rs @@ -3,7 +3,6 @@ use crate::model::test::MockSeparableNonlinearModel; use crate::test_helpers::differentiation::numerical_derivative; use crate::test_helpers::get_double_exponential_model_with_constant_offset; #[cfg(test)] -use crate::util::to_matrix; use approx::assert_relative_eq; use levenberg_marquardt::differentiate_numerically; use nalgebra::DVector; @@ -241,14 +240,6 @@ fn matrix_to_vector_works() { assert_eq!(vec, nalgebra::dvector![1., 2., 3., 4., 5., 6.]); } -#[test] -fn vector_to_matrix_works() { - let vec = nalgebra::dvector![1., 2., 3., 4., 5., 6.]; - let expected = nalgebra::dmatrix![1.; 2.; 3.; 4.; 5.; 6.]; - assert_eq!(expected.ncols(), 1); - assert_eq!(to_matrix(vec), expected); -} - #[test] fn vector_matrix_copy_works() { let mut mat = DMatrix::::zeros(4, 2); diff --git a/src/statistics/mod.rs b/src/statistics/mod.rs index 6f0256e..962ada1 100644 --- a/src/statistics/mod.rs +++ b/src/statistics/mod.rs @@ -1,17 +1,19 @@ -use crate::{ - prelude::SeparableNonlinearModel, - util::{to_vector, Weights}, -}; +use self::numeric_traits::CastF64; +use crate::{prelude::SeparableNonlinearModel, util::Weights}; use nalgebra::{ - allocator::Allocator, ComplexField, DVector, DefaultAllocator, Dim, DimAdd, DimMin, DimSub, - Dyn, Matrix, OMatrix, OVector, RealField, Scalar, U0, U1, + allocator::Allocator, ComplexField, DefaultAllocator, Dim, DimAdd, DimMin, DimSub, Dyn, Matrix, + OMatrix, OVector, RealField, Scalar, VectorView, U0, U1, }; -use num_traits::{Float, FromPrimitive, Zero}; +use num_traits::{Float, FromPrimitive, One, Zero}; use thiserror::Error as ThisError; #[cfg(any(test, doctest))] mod test; +/// contains some helper traits for numeric conversions that +/// are too specialized for the num_traits crate +pub mod numeric_traits; + /// Information about an error that can occur during calculation of the /// of the fit statistics. #[derive(Debug, Clone, ThisError)] @@ -24,7 +26,7 @@ pub(crate) enum Error { Underdetermined, #[error("Floating point unable to capture integral value {}", .0)] /// the floating point type was unable to capture an integral value - FloatToIntConversion(usize), + IntegerToFloatConversion(usize), /// Failed to calculate the inverse of a matrix #[error("Matrix inversion error")] MatrixInversion, @@ -70,26 +72,32 @@ where /// # References /// See [O'Leary and Rust 2012](https://www.nist.gov/publications/variable-projection-nonlinear-least-squares-problems) /// for reference. - #[allow(clippy::type_complexity)] covariance_matrix: OMatrix, - /// the correlation matrix, ordered the same way as the covariance matrix. - #[allow(clippy::type_complexity)] - correlation_matrix: OMatrix, - /// the weighted residuals `$\vec{r_w}$ = W * (\vec{y} - \vec{f}(vec{\alpha},\vec{c}))$`, /// where `$\vec{y}$` is the data, `$\vec{f}$` is the model function and `$W$` is the /// weights - weighted_residuals: OMatrix, + weighted_residuals: OVector, - /// the _weighted residual mean square_ or _regression standard error_. - sigma: Model::ScalarType, + // reduced chi squared + reduced_chi2: Model::ScalarType, /// the number of linear coefficients linear_coefficient_count: usize, + /// the number of degrees of freedom + degrees_of_freedom: usize, + /// the number of nonlinear parameters nonlinear_parameter_count: usize, + + /// a helper variable that stores the calculated unscaled + /// standard deviation that we can use in calculation for the + /// confidence band. + /// calculated according to eq (97) in + /// [this example](https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#confidence-and-prediction-intervals) + /// with the difference that the square root has already been applied + unscaled_confidence_sigma: OVector, } impl FitStatistics @@ -117,15 +125,30 @@ where /// * `$C_{12}$` is the covariance between `$c_1$`and `$c_2$`, /// * `$C_{13}$` is the covariance between `$c_1$` and `$\alpha_1$`, /// * and so on. - #[allow(clippy::type_complexity)] pub fn covariance_matrix(&self) -> &OMatrix { &self.covariance_matrix } - /// the correlation matrix, ordered the same way as the covariance matrix. - #[allow(clippy::type_complexity)] - pub fn correlation_matrix(&self) -> &OMatrix { - &self.correlation_matrix + #[deprecated(note = "Use the method calc_correlation_matrix.", since = "0.9.0")] + /// calculate the correlation matrix. **Deprecated**, use the `calc_correlation_matrix`` + /// function instead. + pub fn correlation_matrix(&self) -> OMatrix + where + Model::ScalarType: Float, + { + self.calculate_correlation_matrix().clone() + } + + /// The correlation matrix, ordered the same way as the covariance matrix. + /// + /// **Note** The correlation matrix is calculated on the fly from the + /// covariance matrix when this function is called. It is suggested to + /// store this matrix somewhere to avoid having to recalculate it. + pub fn calculate_correlation_matrix(&self) -> OMatrix + where + Model::ScalarType: Float, + { + calc_correlation_matrix(&self.covariance_matrix) } /// the weighted residuals @@ -136,8 +159,8 @@ where /// /// In case of a dataset with multiple members, the residuals are /// written one after the other into the vector - pub fn weighted_residuals(&self) -> DVector { - to_vector(self.weighted_residuals.clone()) + pub fn weighted_residuals(&self) -> OVector { + self.weighted_residuals.clone() } /// the _regression standard error_ (also called _weighted residual mean square_, or _sigma_). @@ -148,8 +171,17 @@ where /// where `$N_{data}$` is the number of data points (observations), `$N_{params}$` is the number of nonlinear /// parameters, and `$N_{basis}$` is the number of linear parameters (i.e. the number /// of basis functions). - pub fn regression_standard_error(&self) -> Model::ScalarType { - self.sigma.clone() + pub fn regression_standard_error(&self) -> Model::ScalarType + where + Model::ScalarType: Float, + { + Float::sqrt(self.reduced_chi2) + } + + /// the reduced chi-squared after the fit, i.e. `$\chi^2/\nu$`, where + /// `$\nu$` is the number of degrees of freedom. + pub fn reduced_chi2(&self) -> Model::ScalarType { + self.reduced_chi2.clone() } /// helper function to extract the estimated _variance_ @@ -179,6 +211,143 @@ where let diagonal = self.covariance_matrix.diagonal(); extract_range(&diagonal, U0, Dyn(self.linear_coefficient_count)) } + + /// Calculate the radius (i.e. half of the width) of the confidence band + /// of the fitted model function at the best fit parameters evaluated at + /// the support points. + /// + /// This function was inspired by the function + /// [eval_uncertainty](https://lmfit.github.io/lmfit-py/model.html#lmfit.model.ModelResult.eval_uncertainty) + /// of the python package [lmfit](https://lmfit.github.io/lmfit-py/intro.html). + /// It will produce an identical result (within numerical accuracy) for the + /// same fit, given the same probability level. + /// + /// # Arguments + /// + /// * `probability` the confidence level of the confidence interval, expressed + /// as a probability value between 0 and 1. Note that it can't be 0 or 1, but + /// anything in between. Since least squares fitting assumes Gaussian error + /// distributions, we can use the quantiles of the normal distribution to + /// relate this to often used "number of sigmas". For example the probability + /// `$68.3 \% = 0.683$` corresponds to (approximately) `$1\sigma$`. + /// + /// # Returns + /// + /// The half-width of the confidence band with the given probability. To + /// calculate the confidence band, we first need to calculate the fit result at + /// the best fit parameters. Then the upper bound of the confidence interval + /// is given by adding this radius and the lower bound is given by subtracting + /// this radius. + /// + /// # Example + /// Fit model `model` to observations `y` and calculate the best fit and the + /// 1-sigma confidence band upper and lower bounds. + /// + /// ```no_run + /// # let model : varpro::model::SeparableModel = unimplemented!(); + /// # let y = nalgebra::DVector::from_vec(vec![0.0, 10.0]); + /// use varpro::solvers::levmar::LevMarProblemBuilder; + /// use varpro::solvers::levmar::LevMarSolver; + /// let problem = LevMarProblemBuilder::new(model) + /// .observations(y) + /// .build() + /// .unwrap(); + /// + /// let (fit_result,fit_statistics) = LevMarSolver::default() + /// .fit_with_statistics(problem) + /// .unwrap(); + /// + /// let best_fit = fit_result.best_fit().unwrap(); + /// let cb_radius = fit_statistics.confidence_band_radius(0.683); + /// // upper and lower bounds of the confidence band + /// let cb_upper = best_fit + cb_radius; + /// let cb_lower = best_fit - cb_radius; + /// ``` + /// # An Important Note on the Confidence Band Values + /// + /// This library chooses to implement the confidence interval such that it + /// gives the same results as the python library [`lmfit`](https://lmfit.github.io/lmfit-py/). + /// That means that the confidence interval is given _as if_ during the + /// fitting process the weights had been uniformly scaled such that the + /// reduced `$\chi^2$` after fitting is equal to unity: `$\chi_/nu^2 = \chi^2/\nu = 1$`, + /// where `$\nu$` is the number of degrees of freedom (i.e. the number of data + /// points minus the number of total fit parameters). + /// + /// Scaling the weights does not influence the best fit parameters themselves, + /// but it does influence the resulting uncertainties. Read [here](https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#reduced-chi-squared) + /// for an in-depth explanation. Briefly, the expected value for + /// `$\chi^2_\nu$` for a fit with `$\nu$` degrees of freedom is one. + /// Therefore it can be reasonable to apply the scaling such that we force + /// `$chi^2_\nu$` to take its expected value. This will correct an overestimation + /// or underestimation of the errors and is often a reasonable choice to make. + /// + /// However, this will make this confidence band inconsistent with the other + /// information from the fit, such as the standard errors from the covariance matrix + /// or the reduced `$\chi^2$`. Luckily, it's easy to obtain the confidence + /// bands with the errors exactly as given. Just multiply the result of this + /// function with the _reduced_ `$\chi^2$` of this fit. + /// + /// ```no_run + /// # let model : varpro::model::SeparableModel = unimplemented!(); + /// # let y = nalgebra::DVector::from_vec(vec![0.0, 10.0]); + /// # use varpro::solvers::levmar::LevMarProblemBuilder; + /// # use varpro::solvers::levmar::LevMarSolver; + /// # let problem = LevMarProblemBuilder::new(model) + /// # .observations(y) + /// # .build() + /// # .unwrap(); + /// + /// let (fit_result,fit_statistics) = LevMarSolver::default() + /// .fit_with_statistics(problem) + /// .unwrap(); + /// + /// let best_fit = fit_result.best_fit().unwrap(); + /// // get the unscaled confidence band by multiplying with the + /// // reduced chi2 + /// let cb_radius_unscaled = fit_statistics.reduced_chi2() * fit_statistics.confidence_band_radius(0.683); + /// // upper and lower bounds of the confidence band + /// let cb_upper = best_fit + cb_radius_unscaled; + /// let cb_lower = best_fit - cb_radius_unscaled; + /// ``` + /// + /// # Panics + /// + /// Panics if `probability` is not within the half-open interval + /// `$(0,1)$`. + pub fn confidence_band_radius( + &self, + probability: Model::ScalarType, + ) -> OVector + where + Model::ScalarType: num_traits::Float + One + Zero + CastF64, + { + assert!( + probability.is_finite() + && probability > Model::ScalarType::ZERO + && probability < Model::ScalarType::ONE, + "probability must be in open interval (0.,1.)" + ); + + let t_scale = distrs::StudentsT::ppf( + (probability.into_f64() + 1.) / 2., + f64::from_usize(self.degrees_of_freedom).expect("failed int to float conversion"), + ); + + // this is a bit cumbersome due to the method being generic + // and t_scale only being available as f64. We are just + // multiplying t_scale to the unscaled confidence sigma element wise + let mut confidence_radius = + OVector::::zeros(self.unscaled_confidence_sigma.nrows()); + + confidence_radius + .iter_mut() + .zip(self.unscaled_confidence_sigma.iter()) + .for_each(|(cb, sigma)| { + *cb = CastF64::from_f64(t_scale * sigma.into_f64()); + }); + + confidence_radius + } } /// extrant the half open range `[Start,End)` from the given vector. @@ -228,9 +397,9 @@ where #[allow(non_snake_case)] pub(crate) fn try_calculate( model: &Model, - weighted_data: &OMatrix, + weighted_data: VectorView, weights: &Weights, - linear_coefficients: &OMatrix, + linear_coefficients: VectorView, ) -> Result> where Model::ScalarType: Scalar + ComplexField + Float + Zero + FromPrimitive, @@ -239,39 +408,80 @@ where DefaultAllocator: Allocator, Model::ScalarType: Scalar + ComplexField + Copy + RealField + Float, { + // this is a sanity check and should never actually fail + // it just makes sure I have been using this correctly inside of this crate + // if this fails, this indicates a programming error + debug_assert_eq!( + weighted_data.ncols(), + linear_coefficients.ncols(), + "Data dims and linear coefficient dims don't match. Indicates logic error in library!" + ); + debug_assert_eq!(weighted_data.nrows(), + model.output_len(), + "model output dimensions and data dimensions do not match. Indicates a programming error in this library!" + ); + let output_len = model.output_len(); + let J = model_function_jacobian(model, linear_coefficients)?; + // see the OLeary and Rust Paper for reference // the names are taken from the paper - - let H = weights * model_function_jacobian(model, linear_coefficients)?; - let output_len = model.output_len(); + let H = weights * J.clone(); let weighted_residuals = weighted_data - weights * model.eval()? * linear_coefficients; - let degrees_of_freedom = model.parameter_count() + model.base_function_count(); - if output_len <= degrees_of_freedom { + let total_parameter_count = model.parameter_count() + model.base_function_count(); + let degrees_of_freedom = output_len - total_parameter_count; + if output_len <= total_parameter_count { return Err(Error::Underdetermined); } - let sigma: Model::ScalarType = weighted_residuals.norm() - / Float::sqrt( - Model::ScalarType::from_usize(output_len - degrees_of_freedom) - .ok_or(Error::FloatToIntConversion(output_len - degrees_of_freedom))?, - ); + let reduced_chi2 = weighted_residuals.norm_squared() + / Model::ScalarType::from_usize(degrees_of_freedom) + .ok_or(Error::IntegerToFloatConversion(degrees_of_freedom))?; + + let sigma: Model::ScalarType = Float::sqrt(reduced_chi2); let HTH_inv = (H.transpose() * H) .try_inverse() .ok_or(Error::MatrixInversion)?; let covariance_matrix = HTH_inv * sigma * sigma; - let correlation_matrix = calc_correlation_matrix(&covariance_matrix); // we don't calculate R^2, see the notes on the documentation // of this struct + // what follows is the calculation for the confidence bands. + // this is the correspoding github issue #29: https://github.com/geo-ant/varpro/issues/29 + // also see documentation for the python lmfit library here + // https://lmfit.github.io/lmfit-py/model.html#lmfit.model.ModelResult.eval_uncertainty + // which references the formula here: + // https://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#confidence-and-prediction-intervals + + //@todo(georgios) this logic assumes that unscaled_sigma + // is a column vector. That means iter_mut will just iterate + // over the elements. + let mut unscaled_confidence_sigma = OVector::::zeros(output_len); + + unscaled_confidence_sigma + .iter_mut() + // we have to iterate over the columns of J^T (transpose) + // so we iterate over the rows of J. + // @todo(georgios) check how to make this more efficient (performance) + .zip(J.row_iter()) + .for_each(|(sig, j)| { + let j = j.transpose(); + // in the linked docs, the code is given as a double sum. However, + // this can be simplified into a quadratic form b^T A b, where b are the + // rows of the Jacobian and A is the covariance matrix. + *sig = j.dot(&(&covariance_matrix * &j)); + *sig = Float::sqrt(*sig); + }); + Ok(Self { covariance_matrix, - correlation_matrix, + reduced_chi2, weighted_residuals, - sigma, linear_coefficient_count: model.base_function_count(), + degrees_of_freedom, nonlinear_parameter_count: model.parameter_count(), + unscaled_confidence_sigma, }) } } @@ -315,7 +525,7 @@ where #[allow(non_snake_case)] fn model_function_jacobian( model: &Model, - C: &OMatrix, + c: VectorView, ) -> Result, Model::Error> where Model: SeparableNonlinearModel, @@ -336,7 +546,7 @@ where { //@todo this is not very efficient, make this better //but this does not happen repeatedly, so it might not be as bad - col.copy_from(&(to_vector(model.eval_partial_deriv(idx)? * C))); + col.copy_from(&(model.eval_partial_deriv(idx)? * c)); } Ok(concat_colwise( diff --git a/src/statistics/numeric_traits/mod.rs b/src/statistics/numeric_traits/mod.rs new file mode 100644 index 0000000..faf8c43 --- /dev/null +++ b/src/statistics/numeric_traits/mod.rs @@ -0,0 +1,45 @@ +/// a helper trait for floating point numbers that can be cast from +/// f64. This is only implemented for f32 and f64. Casting f64 into +/// f32 is typically associated with a loss of precision. +pub trait CastF64 { + /// helper for the constant 0 (zero) + const ZERO: Self; + /// helper for the constant 1 (one) + const ONE: Self; + + /// make an f64 into a value of this type + fn from_f64(value: f64) -> Self; + + /// make a value of this type into an f64 + fn into_f64(self) -> f64; +} + +impl CastF64 for f64 { + const ZERO: Self = 0.; + const ONE: Self = 1.; + + #[inline] + fn from_f64(value: f64) -> Self { + value + } + + #[inline] + fn into_f64(self) -> Self { + self + } +} + +impl CastF64 for f32 { + const ZERO: Self = 0.; + const ONE: Self = 1.; + + #[inline] + fn from_f64(value: f64) -> Self { + value as _ + } + + #[inline] + fn into_f64(self) -> f64 { + self as _ + } +} diff --git a/src/statistics/test.rs b/src/statistics/test.rs index deb6d40..8a44298 100644 --- a/src/statistics/test.rs +++ b/src/statistics/test.rs @@ -1,7 +1,5 @@ use std::convert::Infallible; -#[cfg(test)] -use crate::util::to_matrix; use crate::{ prelude::SeparableNonlinearModel, statistics::{extract_range, model_function_jacobian}, @@ -114,7 +112,7 @@ fn model_function_jacobian_is_calculated_correctly() { model.eval().unwrap(), DMatrix::from_column_slice(3, 2, &[12., 12., 12., 13., 13., 13.]), ); - let calculated_jac = model_function_jacobian(&model, &to_matrix(c)); + let calculated_jac = model_function_jacobian(&model, c.as_view()); assert_relative_eq!(expected_jac, calculated_jac.unwrap()); } diff --git a/src/util/mod.rs b/src/util/mod.rs index d52673b..9ca7196 100644 --- a/src/util/mod.rs +++ b/src/util/mod.rs @@ -105,12 +105,3 @@ pub(crate) fn to_vector( let new_rows = Dyn(mat.nrows() * mat.ncols()); mat.reshape_generic(new_rows, U1) } - -#[inline] -#[cfg(test)] -pub(crate) fn to_matrix( - vec: DVector, -) -> nalgebra::DMatrix { - let nrows = Dyn(vec.nrows()); - vec.reshape_generic(nrows, Dyn(1)) -} diff --git a/test_assets/multiexp_decay/conf_1000_64bit.raw b/test_assets/multiexp_decay/conf_1000_64bit.raw new file mode 100644 index 0000000000000000000000000000000000000000..d63444e274401b6d00cbc8def01fd45f7b69eebf GIT binary patch literal 8000 zcmWkzcRZEfA5ThVNQe+oic%C2dY6QIpXa)R>)M$~gOW-_%4|s~%BrlAq?GhcXhidjD z_vu|GKj&;h!9cQ5d9*J=r`2XxjrgGRv-B_dZ9b6ca^8Qc%NuJNms~U6?2W*drKP{W zd*Mt=S;(UhF9>#X#tw^m!L%(Ua7~6M^f8&8Qzx_nmAs(0;=3_rg z-2;c6^MAZMyAfem;`HPGx}(4POPY|cJESEqURHVMhEdg`4ViXsxba=X%FpTAUFaF79l6%o)vh-NQXBobf2V zKs#^J30eOs`&Qg{f=cC~?T+zID1M@OV+G#{Cud8I{95P)(SaTr-w8)VsO@=b^v)4w zp99QgN*$4x7&%#R&JiX?6&8&L9AP~xXg&NL;pyCAmcn-g>$I2D4NFJtADp|GW8?_Q zAU{tX9Y?762eqlGJE9?O;M8YDM<}YuORtwDbY~=Alb3YFAEkr8mWh&QuOd00>5lLn zP^*a(c7)u+(Y;Yq9dX9|L7XDNb${CqtRd*Vxp!GD!5Jm~|4PVzI<2R4%cnUa)S*b~ z2U(vrnr=3Kh9mT;KcfbOZa|cEo|%{YC0VRhaHmg?-4xH zWN~0M!Ch_&Q&{BxSt@!T`UoA-S_AJfLf5p(%Q~9ylg?@8#Y`b|tUqZFjSF!9z4(>q zzXhmy>Y^<%EI?O-RbxTF0BN$p_5XDXaB>G*NxegW1krP^yPE~*_RBaV+bF<^$2~47 zH3IB$n(7w&Kmf6~bu(U;3h*bp%qOx?fOS1rlTxz;NY>TfAeSOQz}DOk<3|M8*fF8T zj}Sm(iEA^}d;LbEoerqJap>?~pgtY}&AL_fj zSU~`hcM|V@P8VSIyLOfQk>&yiFvfEl zZ=ZEQ&WFn{Ok*5yA>`li_MHy6w(Ux+tD6HXPrJ!a*f_xAb!#qvr2`J@zH;W>p$eeM@;08#RVJqsLG2p_Fl$E$|cn-1tmUY_iYh>Je7}o*}o$8e6z<= z%bD#fn(VQorD5Y(r9I}Re3EUE%^1t0wE@i@qCwIa_R zCr;no9BO8dFppJJCKlQw=lsO?YqIuO9cJNqeB2IYVS*>AU3S>fIv&&f$__O?9U}Y+ zJ2?19o)Wua2WlwLXXK0>=4SaV(mrGd>lcqUo(-^r?9^8rmWv%2j^oufHg-_Ung24; z&<;(z=3CF#w1XfaFmkG_9W0&Re8_gy6E$VHC#8V#mB7eifI0j(S^Ok7)(3s_uKb{7R#FuCwPKjtoBL*z`hHOSn+ zMdpPUjoa6Akrk?V<^)Yq+Q=t1qiE<08n_`rtWJVvn?cn=)2& zp%S>OrOAW~*FP%_Ob9->99C3i%tgU9uj(LzW!p6jXo8~b7GD-Y9hF|~tpsn@o0ME8 z_>jffF_o-)W%rkLfMAv6X^Yu}uHF~dIgbcF8k4VzBK*t>kEpv5{>o=t_cG17@EYGS z*M`W|Qa`!bfynu`FP6DuEf+qM2P@Nv{QDK`m)2Skz4!NT6t?7|Hb%XQOY{uc?JauF ziVF*k>SsOHTqwrJiJLOG*xmR-?*y@n(zld`c8ZH(w<{h7#9q5o?kL6DauJ~-)zHM| zqH(fi!(1*GD3eU%^SC(Es&9G5j*BUIm5NV^y^l}(Qu&?FMKUi}RaC&mt&DAv%8p!g zgn#>`=ETK~tPxpxXJU`riE(2Dza9-QswIB==^Gh;)RiDNm~G?6g}azcNKox}lir*yU?!T!B*Eo0<;_!bW5 z13_)ccUGqf&d%i}ttD9di!YY`Z;AXfHF5ns(Z{N9rwL8;o6t#`{6h2%JTf;(K=ePaf4w`C*hNAt(Y}Mk zU2 z^4U2~0c*Hu{k}f?1o8j!GkqH_ka&@JulOt8kmS>4O}H9x5xTAIuE;VjO4s(}r0NoV zM=yBGEh6!_F?PD278iGgS2D8IxtPIJz5hss`04cS{ZHm_Vg5CYacwriFGEX1WQafL z2gVwbjNQ^zBCXFbC$Myrgdf96Z13 z;oaWL!5!PW*8e&=7%XpkvAdmvbJeaaRQg4yLApe|I=YtgPEIUd(|^i)noKRSx{5Jl!AVanRYA=MZ^;gDm#+{$?YFV~TE3_)_n zY7U+awhp?i;K1tU=zFmx94Hpoe!Z#1!Ou-OckEO-c(W+_O_w|e8}w=x2F&CjPfu>m zpa=)lyMEnjhDT)zl77fWZBo;thsA8j6jguO zn$1RkPE)k-X*QG|OM9G*W+U)y+(YBtY=p^PiYW4BW0v-%^HvUQTpOt~e`>|X(2m)s z9afR`>RHp?X|u84&a##_hmDGKn%yVLMtzFGg6KcCP>#4RG}&Vdw_eS%v}Rj;4UW9I z>bWg=w_Ft7mf2#v-PN)aS8Wl`ms-6(!xjgFjb2Y2w}r9x!0@XGTSzvyvTg<1BG2>r zkrWSGY(CyOoW!+-=(bfY=PhheG}P?*-wIo}THkN((z1m@MQ!ME1zQ|HFxMe=hAqZb zQuhx3W?|-~@TmR0EUaEDYB;-v1-Y%q*w5=&`0W+%pIXU6(Yg`Os2ePtlqj)0l+D8E z&$8B{Gb}8f5tP{z$3pdUZHf+MK}9dT{K+;J=(OR&EgM;oGW?~a$z$PQR`kcI8_55A zbKXxkW+B(7w_~X;3jwp1_ynu5uz9WJie6b3UIYja@1M?s$WQJbqdzn%4a+r@zSF2N zwJ_BDL}QK6#XT->XgE%M{BrXdjlQMNFRZSh;a)w_)LTqr{XQG_=a*r)bFcgSonQ#f}_`eAi6g}0igV`C2z{$KZIhD1=f^z7^XvtblYymwXX z387FrZTFn5dni;~$}Cq7rVuTaA2Gg@!l^Y2+{JfLu=c+y!U&+yvwo%G&28jaW4;_~ zD}|N0dC78HDCkPmtBLqiV4Zb5thJef{yy&doqiP7tib-RO+>zf7HKgA*K|Dp$|0Cu zSfpu9aDLI)Ztwqp7t(uji6G_czEY8_yV#r2e1%~9m$_9wgx=4Z{_RXcf6s9D7B1l@ zd|A;xobWF@7nfN>rog@zlP=(2qj#M>>V%I>EirO17Gl-NbNe98u8Vy7qSHJZ<(2tAP{bI-+6I2_kt z5pkHporF8hhmKJ=nR0YVCb9QCiRF^7Pf-Yer_iUCOyQw}H4dI5c1tumK1uwP?_ute zMEvF$&zoG8Lm@U|-_9SGC?uC0$SEtJaAQ`CdwLO}r_(5T;wFXqr+O_(r4+^k=e;c; zaWQ|)dBclJ3Xet4TobOM5PYL`$GWE^zCI0pI$290#K%0VSws%oDpE%uY`Zo$Cb$5T~4N~BCg>*WMka(z)kQV<- zVN1lOsTCx?CC;C!2^FT1w)rXk$9)bfZP4MXV7-rAfra;bczDOfibeqYch?+}r|v&bOWkMF5Lpu{x0B>=>e0kD1t%JN7av)7%Z-N8 zs>m#RFOr{0LJwYTqT%+aX49svG|bz-Z*2*rG4Sni7%P~@GIrCQn<1pGC~5^rhttqf z>s;t}kVep>3bX4mG@2${Y`+|)LEkUV6-%J;w(jJ(E~zhSZHugol4;ny?T%QTO2e~g zTxLNgjcpg66ia2%c$FQJ*PKg3JWb6$H=oqC{4dh%LK<(H1XKQAr|~H-#wh$2jZJk+ z{~Ic$v3s+nkZTzYebd~cJEXo9EN%HZ^nk|vuELkvk7-;spLo9EDXDAB1PYVfn`$qg6|JHI* zKaCylKVP{%Kx2=({>4QKt#xcdwrz+=ZJU$VWaVC{U%i7snf1II_ z6*Q8ebCQOfR#+|j7>)1A1uI-g-@G60e`v!Y;+J%@Ky{KA%Y_bTeBMj^@!qoU0O_N8 z3w>+#p2jOr(yv#KtBAYOD2vikHguqI(a)#0dq>R)J!kvc+*KqF&RBO^>ytdKUU2uB4vi5X^$d_Y@iC=uVX`W*V^n)c zrvi<(WxDcnWoSrcy-_s~qp|kg_iBr&G)mXaO<6uh>Q+@ny4)~@MOTX}yGR`itaO@I z*h!v^g#G_t?}9k)Jv~jNuDdboch{3XvD?tA_X+7cYSFG-(r?SZ>4@jup)liDm|{Qa zvj?hfjLpmce;w7ivf&~H!9`U#llr|;cJZRkCn>lawg!8Wy1&0EXw07U1;fdGTUYF( zP*e2kt|aL%_uDHUydwQcGxpDx6CM;s7dE6(juf_*wm1G@lRlrKC3eD^!kocWiN&BW zyWpr!p%I0=18#g{J<^Xf3Xhg9ApJ0K;e+|AWIu3bbH$n26n0F%Z8s@FK{`aG-&2^v zO~#?W_y02Sqf*gi?hq3th z2$~i9`^uU4I#jl9NeL6U)x9j@Dich>C#8qkWW8?IGOVH zw`pYBK_;dh4Kd}1Fp)W9cE*iBCZy`sBo+Lauq*t&DAb*a^L-_LeS9X?{b`W$VliQ* z;N#J3$;2n4aHf+P6Q3XMavNO1gbn-D(*#{6N=kEnF41J7{$c64c4a1Fce)y$oW(@5 z*Vcy~;!McA3Mn?2$^_-&*QfZG0TG?bK=~mC^ullRq z&UZ(=VjxAwDD+tk1FMJI<+UC$5JR)8Q_C2ba&OCHwOb4%Sjlg^ca?$sf#R{fxeVk4 zZd;5@2AW$U?G@CG>YO3m1ypjR+Rl;uB%NQ8*;?1nqVW4A5VZGEsBLC6Q z_+SkNW+cWRlTc;g_LA&-wR0HQ`Xu;i$!sE56tB5XhJoeN7b;6jF%WQM*4$7r296X~ zPgfFQfMtFs_Ulvz%Dz4`ZvAJ2=CIOS;RzdL7TsIs_SXhmx((QUqc+e`?)!OZ#0L7+ z_2K7#+Mqu8#o>EDY%scfN3{H)4fbT7`kvKqgYZMR;rrbNRbwLh-rsDn6zMLf`)u%M zfw}(QUh@C@f446ZlrcP}9ZJwsd#-UJ!O+Bjk~V^@f#?O>$p38t55^71dcW4>Q-*|I zvJFRO523HyoL@RRU;_i^FMn^ffXCnYHy1&aH{5q^CL+Hjvp&`l$B-RdCc^4Ir0o- z&M&KpP-b9ddv;;!JO-#Qr#b($7>K@WWs*F+|Mf@aPBo)OYd>ZeJYB@7daKYxYtQP}kk7(W5u!RAYdtEh;gBaK) zA(3*N#L<21YC0dyz*N88C;uI1z`1VW$;%}ET+h2^Wn3h2-LfLRrGTIy!OHJ8@!#hi zZfhPeV3))H=~>Ibe3`;kFW(Zsmj~ZG`k8?(X@L<*0}NQ3Cy8~9G4OZE38C@FItBBY7`dRaUVk|gw?BQ<6JEo_om{OE4VsD96pv1DX5zPZ*_RFeB>&3y z9$=DuH}zm$NQ-0Q;9iBtCMis?JLC3fUSi_i*KhjHrA+kfl~!wc!bH=(h22+3-I=v4 zHSjU1M^{ey{?i>}!sYtK?=6xPOf?%Xe^8^KUY=sFsZYV;$@oa=dQv~z?r5brko~8U zSW_k0*XIp%9x)6jb@af`k{^j=pZfgZr!lf$K07pLko6zgS0B;u>n7)?v9h}pb!0zV z(!?~H7^RT0eTH5z*_XEO?wHiqq~VogP-JXEo}Wh6#gqM`{9%n;b+SK{Uh=ePf}9hl ztL^sFI8DQ)FHiJ%0gbH_-t){!{~dZUHQ#{rah=ri-axWnc&Kq$YEXiOHT_}9x3yT9 zA{T2_Zq9;u{{^OmfQ6EG*YeVLuyBk{H`sZEh1uVw{Nr+1IDYA^_s0qr_I>6HhrVaQ zuu%W6{V)sLJ;x<>O4~xLHnqD>#}?&0m1=)0TkKk2v&+iM7O&Vv3@>v2mus8aR+egu zl)v@YxTUrbcdy_tdToo5(BwI5e%K->H~#IhnQZj^8thY8!iG*nTW~vr4UYN71z-Kx z*p(w3uNljRfyxxG^c*%OVxMecKO*yys>n3+4mPZHT@^y;c} z@av(|pIOVv`vv|^a%(ww=r_-11)0}~dCLuM_26LU;X^+wc5txzKieX`gB+YSTiks6 zB$?NYSwuQzk-5^R%07jg95`K-QUCskg9G8`F1Ebp;GbaApH4Ej+PX7uljtZ1ZKInW ztrz9u-X6}3WF;<4^@Qh5C37|H7{!m#Wd3^KgxI}#WS-Y3ojTA+=7`y?y>#YIG7o(r zF%}WUg$Uc&dCwU#CvG`AJLVGkzOAXO;6E-_M?YWoP28OGul{_-HI zy-)#-p^ZQyV%VRtIS^8%kH!z=WZnpsTe!_ z*{q#E=bRnZcW#wgRcMF$3jxO75A2}qQQ>;=jU9R&D~6?9l#+j(xNmBY51}_DTsZcS80GM$`PxJ4l*hRT`|QcNVQqeFf<4%~q2iKF z?q^;qT3Xz($35S>nFe+CNU1+QwxWa3JDRy=?WjGXrK|g$CCNPyN?*rl@UhlbkAK&Q zkDcoSqQ)t54^(@-n)2kMt(LEUc@H1oU#=0CJ;ukc3+;Z<7x^d;k!lnv<)d0_Uf8Kx za^Dnodg+o5K3+Y3rr0pb2eT(d`GAxIif4%0nQA&fqA@O0$ixAD(-xg*uysJ^pJ!XI z`Z(aC?reqAp$-WDlQt(N(E;boWYZ7kIpEg&!$#5f9KgPpp_TN?0WXwA+6uZI5LKD; zzJAgHUC|}O6SD=-pWh{Bpew+lrLVPu){*Bs?Urx^6,m"Ҳ$c?d?UrC)r?+?d%c/L>[cx^6,+?4}5>80ϾzHm"Ҳ$c?d%c80 lU?W0XiT?d?/L>[cϾzHW0XiT?t1GU? \ No newline at end of file diff --git a/test_assets/multiexp_decay/xdata_1000_64bit.raw b/test_assets/multiexp_decay/xdata_1000_64bit.raw new file mode 100644 index 0000000000000000000000000000000000000000..9b66f05e4a0786c582b430ad1efc7c1441e48dbb GIT binary patch literal 8000 zcmXY$1(a6R7sZKTf|&07y5ofkQf~_gDs9k8TC_@ofXH7Z#2}?Xq)R|hKuV+yLPT1c z0qGv%fA4#H7HhEvX1?#fbI;lPx6fCql=W5Sh0?#-hK@6=zvmcbFMp(~o0jtXRYsY~ z&z17e{E7d*!6;MIer5i(&HVqZMw$0!opU8`=g;jl%Kq|iPFUH)pa0z$R{pEGx%2~+mfAMwxHirH0y_=pEhxmH`7#}Z9sT90$n6G=(`1tSd zJ<)#0`1;3&xCyg!X_ZR0i_?PeZuTjgt^QVGwr})07jV(34*;aJ)G~fS> zQS0S*_uE}(d5*Kjmbacg>Fj!r=Q(fG8kq2(d)0ZK>w>Z6i+}EWA70@3E*iDxEQk%% zDdstgjbDE3pB5T@iRZm+)c);7TjaoHp8G$;EKmLwp8twb`_9E@wN6+09M_Bvxn|5x`qiKJ;2NLnx^YTWH-GwdKHqi2;Tr!;L9H8n&KpLD$fN3o=Wg(MZyK|! zyjiDJM;;STF_*KqXt>s@z;yR6e)qr=b|Jw37aSg(7=?2+H~@h-f_y4^EQjjQybuhD(h z@4n%j(fecnjr*+Q1Eb@DV|+>tRY zJnSLs{LtvQ@565+o=2?rBjeknXBR|wJ!0J-8>i1TT&{F~%=$kzoY#L?t)@QVJv=cw zmK6M~t$xD$D207xrzvitfzaTME`am4cSqwo9I~G^{=> zjhC~=UdiiS8rHp)#{2o(-^<@#8rB?@#@F_eTlqO>)u-B~*p)*Qvc9Ix@Qu=*wz?{|;;&UHQ(mTtu1>kq2=+|S0r z>Z3TU=+|SBdsG~(xr>8&oyX!}^;JA>43h#74@)QFk@{JMr@bG=!|JnmXrCKDdiTe} zn$vi^Y`#GPtiDUY`xEDF@byiArE3WoJN3wO{+$W1`Y-`2XUTyp4{NT=WB=URF9f=k zht-$madT1M-GPneVd+wNqJ7B&@vI_?}uV36?G;VgJ342Wp=u!J3OC+hfTmLtCY%=F38M`XHv3>ITWX@4CS|#3#pK>;t z^OTI@WZ9u9oU0TJ5b8Rmath}w1?4lh{Whgt3g;{Zle60;O&yWKc}qbc??IEPD^fUj zDOg)D?(3;XQaFDpcuA=9)R4N2vEreb7Z z!r^K2Q#q%p$d&TEHUej#uwA z&T|@yg*(icp2oRO!+`30zMHWzs7XdpczGXVyyR+^1u0 z-KQ4M?2^v;Psd9y{{G+0ap}~7beQe*b2{}P9i1D=Nlm9Nq(fG7R#FD_Ap;{{dFIxv zDjC#?4CFTZ^Xb{`GpH9CSlFU--`OKGs2drm+B)W^*()=s9~qFgnZgL3o@xUnFw^NI(*K) zOzKW1*7i8?>ztBI>Q5G4G6gA%I+TUOz2aZ|u3;ASC=2rR-}T9&E@k0Ce|Z5})TbDhSh^(^XDHWq#)0Vtcgl?{2CxwW&YU)k98NzL?mU9+iU*=RN5 z&};L?XH(CzQ9Md2RyK7l8w18v{Ab?zZ0cJ!OuS0Yq0Z%C@>sN)|7;HRE(d0NcF3Xb z{%Ax+{puUv`6$^v=*g0e`PAopjQTzE zxuqBKsnhw$+q-<&(v$+~bpaObe@VPU0d>0o)efF1x2$6U^}7JO4-KxjY;*y2ya26_ zWRF<(Ljm=?0GE!f5ROnlT`$1ElMN~^k1M3U7ovhG?1j|%LQFY3^o!*!3#s>okaV(q zSRr-45bKI%K?|w>g{c2uqgTH_P)Hvr#F48P#()2)kbY2vE~eoY(HDv!X=O#@BKkuS zM%`WAYDK>y`a}^-xt&!+zbL|@hsELvis&0fF!AJO5&gr4bMm-}Hu{JS5fdzI^b;HU z2s&(8*+yTnVZ-n{sfX*^=r1Fh6$g&-n-QZ8-2%y z^PqWc~Z19*5V51M&(0%mQzYlM-(T{9c{rOk-S{<>`muz@E=8c>qS8ens8|n%= zZ5?B$PuVd>{&=>Xer3nOucg7;>05RbO8a@vPXDswt?_MtZr#jIAG2erpxYyz?esG{ z?o7P@xb+Y_ea()VGQNzp)8Fj)#2k1#ea?>Gr?hK-biJK^XGf-><2HNk^gTP?kaBz0 zPXDuG?u_lfw0U5s4?1vV=D0IQ6CLzJ2c9v(&Ou*vV6dR;W0f8BM+dgeDXIH<0|$N5 zfyB8x-agjOLBDjM*}U=M1|0NF2c}DT8R4LRI&erw;U4yi*W@ zgTCs(28sVi9Q0QQVwZJn_{KE{eb#|Sg6@yUIO(@eO!)qhxG^Vv*NNkDoWf4}uM?hC z6aRdpmXkj0ME4&$T|eH;Nk4XCwMYkTyEy5~PCSHT`=*ey^Fr@f~2jJyf$C}xDk7(?AslNy2%-C zH2UYyPfv|=lQ-Oubk<>ko7~}s7%rjcZt{m4qL~F^xXB@Ih~^QJ>L!o4Vd8v=n_S|7 zq_xvY9`cC?;xHl*d&nssj5%FJ{PB=iJTT?Hp@-b!L81A;9`cI^lHNM@@sMLYSbCvr z%9)WK@{9-OdSHf!T;oB_OR)m!J>(k?KDoUAgHAg<wA@6vQd9|D9Xb-u^ zgEy{K6-V-tf4rD`T?Se&ImnAEH)O!@l83yIboXX$FS*Ez!MD4&Jlot$KJsGQooWKf zyyPS=One{eB`x+F}@{|{S9w!LG z^^&W+Fvq3DOTO|!(%)OjK5~{1Ca$`C@9`03lLeB?AAOg!7+Bd_@&>G17iK60B6g-J4j z`N(fRyp=4k-cOG6!(4ae`pI*CNVj_@>L=IvA!+eKZ9n?yj%%vesZ86S7jXB;3p6IQMo{F0{rAcKZckO;3psY zv0b>)yAS>3!~l{6O^>FGwy8;)ziIOVJ?tG>CdKPSpvLQ-k>8>HeE9 zwFr_|gZNv@Rkt9yHHe}rX(DL974u4UvyS*e-E-bcmcBLXsIML*(TUn%9>K6e2f=Fhi*4 zm7O8-a|q|4g_Xvb}iG^u| z93O#+GgTwx`3U5;v1j85xjq6j?strk?;{Z6v#Dpl2su9j(|$)s$omn9;SvEpLhg@1 zZYQt*5F!6ZAl2l(T@mJh2u!;;9$_AcKom!8 z6y_HNtK@pSwZa^upoM9u3iFJDhk~xx_E4B>6wH*k-d|z9QBX&Y?{I}VN5O?1Uw$}X zjKaL5VDyW3KVLgZVeV0&mxvcTCIwC8dij>Z{G{NPTsJ;Zn4?ro z5;Q(AL1ms&QN8J}FRx2enX6Qs6x-F30+snn#c*@Mt1@S)@RKLt+ij1^oTuWQW9`IBN@dK)Nn>@r>hrg%!wL4GZ*a|^P+~3xj5FC8#VkX z*V)xI=0^>^I_;=gqn^ebsUcfzw?AyEF;8k(FXK~djk!`oI};c+=1UDRrhsb9nHuJw zmjj?NZ)#{Dw&On!*O)ssT$AgKF&guyhH;{rHNMf9Lp4+}7mgb9sD`65E-cZQOEnA? z+x16lH0DzcuC9_WHRe^Nmc~4*;i1IiCmM6Dj+r7YJWkM=Z*|mpPi_!(=3E^Y#CCshq0YRkW3>5!htAxq zLopZlI`glN10p@ttf4ap>*#lDm)J<@%)>hJB@Q;#nTvI7x-AzHI`grPP9~u1%*i_9 zO~cZemvtfsGXMWaktKX#h#uqwsw2nz~-SLgiJguX;x$xDQt96`| z?YTr}zSc2Zd^%zLs558l@X7W37M*!p#~!(^+oLmg>*y)Vc~EEm){$-oTAevu$B#1J XU(%V!b+nOoe_LlR*I|8X{uKTX*48QV literal 0 HcmV?d00001 diff --git a/test_assets/multiexp_decay/ydata_1000_64bit.raw b/test_assets/multiexp_decay/ydata_1000_64bit.raw new file mode 100644 index 0000000000000000000000000000000000000000..27c23d15a611b55e36104e5eac8e337056d99fe3 GIT binary patch literal 8000 zcmWMqcQ}=O7(R=GV;}R_t3e4VMR{)-5tWscq>^MdlmmAHexiViQ+PC_-Gnm^3Br5qsl|6>r=Jy-S@OK_k=d=o@i5> zJ+u+~l5#=aQ5)ySS6q8zt&RVRtL80RsEvo$ONN}JwGr((X0T3B8~a|xw!a_OLX(CA z4EnWT7N#_GvR(^EFGv;k-PeM&<6TQyz7}k5(78QlwJ;jX{*@Q1g`%ZKVTX5WA!2d7 zv6!nC9Inc=E0LJIu+lVlr52bgTyN*8YC(QRt35zg3)Qc)c)O`uFfDzYb!k`=3Q~(b zhz?Espx!YbeXEI8MXgsROEjU}upBnmG_h{@*v^?WP0Ujbxh@f<32s*b$K6*GNA51# zY~ZJf_AZ~~7+Xy|J6U>&L!#Zl=xEJCO}y{%Uq7U+3F87@&J0r%pVy3L)RHLKTkyAK zPy>^G{arTyX+UZ7ukNY08t@5rkMgh9fd6iV87DdMrx=`(o+p9cML=QD+gc zIzj`=`*b9F{53EyURh1TLj&2=#GOyqY2eyOhs^vn8W=zLG)Hf)2E5*cXnE*pU_-K{ z-a%PX&!yhwOsONIbaJ%|iK)#dvQIy&W9nOZdB8h$44iKJ6Z%{o{BOZ3?&a$E+`q}N zGgBQ0F2rv7n5d4Pq)(&UW7Xj%Q4-u7qKd6@nv}? zSKUk<6H%Fa(hSvM_40)eOG_Prn}!qjvej{9>f}*1L9)L7pO*&T)!;I{An{PM8n&!= ztrsp-!|IQDoa*aph>p6d+>xM$)1Q=;KJ8RP{8A$`wzV1pIYBoU7^R9v>Kia z_{cr@tBRK2>!m&VRMB^_dEuissyN(#weZPJRScB0#MWi1;y{sHI4@8Y-&(X?2P0MC zdMI#isD~=73JNn%8LGm3N3MK@vMK}<5AQoUrUI!lH{EliD%jWSfBSEb3VO=N7i3kb zK+N7NZT(k;PxbC_}=do*bu1gnld<%A|Cn5%J^lM_Dk@W666zfWqezeu-+-S zM(nu~GD6Ijbv;tT_E*2M{=2M%pv#eAXU`}>X??9|rLPhKjnk7KZ&ZT)@tNnAP(tH? zUuwCo5=zXCHb9r({GVN&s!96f3;A2?^8uA8It>f3Pt39II-(vl_EY%%jP;NgA}l0)v4`UJQXk+^*x)vQ~}4=tty~vDxl_U3{6jh#1}d2XIy!N zNYwdS|B***wsoom(@GxM$35n`tt8_MgKtFW$zxC3d9Id{ zJcM>E9+}|D!Sc<{Eq6uavF50=UUiQgibC{vz3-HRc+9uj-ZD9CF7O+(O_M|K2mdRR zSLJa1i$YCIgdAu_SshyEdxWl^|wgWumsS-hOBYgz3qi~Qr9 zvII|AuuQ&4%X`SeUFo!`ptUU4d}u58R+q(H{WJaoMzZj)+qt`&C5zZpix?laEH1LI ztQ)1s!ejnhABkZZJYjcD3$@E2b-ZKGe@|u5-*&^gxLF2HvrA1vZ_9wcV|SVJEg5`! zQ`Iw*Ap<#A`L^6cG6;RoE}}=uU}6PV{e74WvNG>0+w79Tzmyt<6Psjk!|{%&yMYWe z6}sYAI>(8%AnZw3=?RVWqPE^%W9*Z&~xgVo0C& zj_jO7HtbjnjvVzP^_x{rt8HZ?*f7X#ojn_iE~rEfJFxM`cyPAFnvMFw6ZUG$*x)^% ztZma{Bkt{}u>T@9((|3tw~DZ_N3y?Xg2u*{M6+k#howPvFt7LcDvkI#M!GlOOXJ0} zMx#y5(xCbJ4GX>^`B;8YVYM_ux`)K=A4p^Qw?&!!yV8)|7Zz=OP8yA2sYlrH(%^`^ zx}I}b8n;eY|Jiy}8d8H<4WPs_M#IeQ&Je| zSvY=vKnfwBu4>=wm4b7{O3fFQQg|IbEYbX23c2qBwjR4Hg;t-K`wy>6LB^~*zUhh- zrn>5?N;0L8a`gCk{0S+nyL8^U^t=>w*n-LrBc-r2I`QR$P$|4%EV#4ILkdTee*Wv+ zA%*i#bVagvNny#7qwDF*q)=I~W7yS13dWOk`Jj1H_%Y#p+g5{|=ho)eG!-el_h0Kc zE+++!(fY0*G%3`dP1_;Bmco;H&QhKfDI`Z$9C!aBiI6S1`OH2^cqX+=C-q8VUz|*W zS%)N`KJev8nIz`e74;>)l|)kW=(i;$lBlH^n(ru(L_clhdTXvE-b+G$MS>*sZ?Dnc ze^U}>{{+@GpOM7L3XKNggOZR{HO?>fm&CV}?2(CJN$|oN`fqKNgwu{ES5uuOp=m$d z9!a8pv2$3hlO#rNNA+B@mBe)$^>6b9-Gqj6C&qS-F;&% zMG}$^j}DvslYrx8M%uCo3F!aVv|Xk}0$MQ|LDo|eNOWKSRlZ&V1D$g#RBuZl#pjgr zz*z~bI`+59BuoN2zx*Y}10`UuY<7l4NI+%g+xg2L66j}~yYynC1lk*$N^Q*~5WLE$ zI{ddd9!gwV*!4pkf0|d(bwu{#P7zqT%MjZ^eJMAF`>sb%5_ z7T*18$V41Yr+F*Z>4@XleWQ7*d@(Gt**!IPQVd7-95w0~6GM3Thd2JcV)(7kV=n$I zhBNQWH`?73!v;^=aqCnue66cq{wH4y>>mZ(mK-sxl@)TlA1Vgds_YG}C&Un6^dMOG zG>ID%=IeNg;f?)msU2=&_^2owrRE`q80UC@2U9WZyBfL5Vuct+D#w&7&BZXB6Sy{< zBZg)Zb(!B|qL4fgvgYScQP@>2`y<{fij=UWPeeLJ5w(IId8|YfZ>tQ}ez_tFHd>{(xw)1?aJ)&rO$M`(6Qxv9ZN`c=jMd3Siz)<6$D89UXRW!6t6g>Wm zTV6{6W2f#)+QajRyv~a?vV{q`pkmHZBg;~X%@yaX12~SM6vPbCRX@33kw61 z8&)*1P;Z-eDj|=B)={cWXfq3WbB_m!b+Pc_YI}!50SoTAQwht-{5}1-r?>gAaQyW3 zAC*U0=z4E?c3~t7*|YhHIgTtO2kbH5XvG43Xm4z$DGSC6GMyxiS%?e?tZiS!fLgT!(&3 zwoL@zW7u9fg(A53tbc*uYY|jh*NTS5ilE$caT1F}gTJ}XE)gQIYo{oA9v6YtpV=V` z4-r_!- z!YE?fiPt<5#^pN3(7Zds5TD$d9B^A0&MA)WpHqafyX4Bt_*7x&J-ybuB3>BEL&43j z4hy4rW1P?C)57@kHUE#HuQ1*>t^JkaDvVJjS9-jzFlIhoi16_cM$K!BKhn#E5jU{+ zildb zgaWTg-+3|HUF+aZLm($w4QC4`_ec~B}%L zVio702yo8+HOa)mExw0MpEKdnP#9YOjEVFb@2r## zCU_%nw(YyY#MGKC<_2UuMkeNNe(LMHl?j0&SLf^2Ow78Om>O?qqV{vjV?QTS z->>t-Vht05T}tuyRGHY+WN)Q7hlwQV%O6@yN#9+@hFUfg^Itr5UoFOjlj#1Etv>~^ zAkD;M>sLXjmwV5B)F}v_RXQYUlSYfSyMAgY*CrptN-@g6$*^erSszWLk;a9J#ffcvMvj*;<)w0siq!b1Q>iAH7j zNPaEizvxN_0oXDecidksfF;x!!~0}hIcfOohehNb8;mD~778F;aL#c0I5YKz*4DA_nXLkfwSvi%mM*?7Z(kLga8s^-1;6a69DDoy{s2%t)WgqQD<|l_&B_?|_aLn^( z+xC47Nauvl>Gflvw$D^e(S`x)?f#eDmJIwInVxgmhXKm1cTb}|8IasyH@Mk=^vPU) znxVo#=$g(gm2wQѪO)yYh*8FwzA_jB}ySr_y8SwtK-prRp?x(zO%a3U~E?hkG zJMk|a|MBldR){cAPg&Q#LV$sxU(S09hUgf*mO8y}f{u)N>xUZq>Db5AkXS)uV9SP) zYx9{5MEuV4TJevLmAv`kUO(yZR(rJ0vW1Q_`afTlHPV5>qm6MjbeO&9S%?%mo;)cl z9Z#Uc>h725?z?noOu7!PPp9LeMCPIQ6?AaME*XuUrQ?a&Ea_+H^v{mdphINm9}uKt($cG_V@kvtb7{w zuzHowH_$L)Ks)BzLWBRg;ZK&=A{M zd9&~^SvO^`x|9tK_m{O=RBR-l$B*NW_M6jSC*o)0VMT-7lf#U8_B5QV3)Njh=0!XQ zcCA}NL)q8)+uE1Wu*oV_trsSQupD4-rr0^vhTBk zB_$jh%AZ|sW=YdvQqOUYP^00(yRs^`IV67Beq=&{hW*|5ZR}>KNTeyY?q$+2ru&G$ zLzae@@d@>#(^M#KSbivXkcwCliNT0}RP?`enPVeB!%F+=w`+!}nAbEN*;Y?Qn*O~8 z3R&mN-D59x`l!(SPuG3TJ1U}n{wk~;p~9f4^59Sv71z=yhn*i&v9@tp_sU8ttadae z9wBv&p(Q6*l~K|AGjQc_Ar+DSF4In}RO~D)C;CdMV7>kt^z0!O36a%qKDVjpNM(Gt zdP+rHreh&9hYDBKu*G7e|7dR0wxBaq$oafm;E+Ux%dXpF%Vgi`5Em zGO4(^Ml1K-9xA%O&zQZ9p+cchWkqBZ74tVYb7q36C`dHWHa$xtddpN*sKEU7R&I;i%-hRj=bWz3K$ z>Kd;+IG2ha2DX1X!l)>m%;jsjQqeZ2QS#P$Dtt>j?mU;J;(B~y^D%uY;-eoZo0(G4 z_xbUgbQ3DvE26}9FCgbB%_gs)R0st;KVZ3miod}F#x-;*Qog*;;VF^7hu3<057EfI z18Nf*gs2#|X`5iEQ4w4xcdlQW)D_EVET>Y*d(nzpJt9;r7}~^Pv#IcMJ6nH*)R}ur z2F5cyA zO%yC?O?eU9NkK$2Tl4%d1#zqdMfFz{NOUgtNqbGfIj4~KmLz}4XmM=(PJxTuZoB)n z6iBDrNrv}RV0!wCgL@4Hi#K_%n?Fdw{NrIBQ;#XITm5!$nS`wLb-;#eO#0B==Z)ICff$JdeVM^P0>k*xt2aZXfB}YMZU} z*_nb%H;;}m^eA|!Jhh}yo&ve?hfCC#Q?RFS@3WLe!CuSiiab57@;*dU(>JaBEHm>>l!HP70plX(q3jpkp%C$`$?*gg;@ z&+YQ~^><7PqO+N&QkwWgrOdfQ^|O3pB&6Yb0fjtY?*3OJKlsF**P#P9a{0u8^hG*l z6?~%Sn3=zN7oUh)A+Rl)PQfJS+m$De_(cBJxzDyv@`(WlKRdNnK9N+|ckyHzpE&Pu zM`)&)PrMSIL;N7~I~tn8gCFw=IbA~&hXy`T)A!6MKa)@7vRZs%68J>U?BJw}H=j@$ z+Wo=zDxWy0F&mXa`r~_l6r2p@6QSk(Ym)c!iH+V7%SS@^#F4w-?g+;5i9HHmj0b)B z#Qsw!JFH9jgr3vy{S9V(;)&q6e^U^jxOgeLrO%Jtb2j(+_ho$I(Yu7X!47<4b-bVX z+M9f0*Ke1rJ%{*&XO*9a@d7?^@5P=E2eSA?(cAgIt9xPKf7;gj7Kc6%FaA&&L?!I%N;D{^ND0b@h2j3d}2B>)U1)iCpI_!b?PDO4fN*a za6a*f&$a)>Ew14a9uG`9ezEw3W~EYs`+Xj<^7QXg;m158#^}_$C@zmU`Kh$l>^qMr z4$965|HvaW9H^Ymt2~0)bNu6t*E}MjbIG<1l{_MiWwgqyfk)_-{_qwd=P@@ceoFmE z>I{?2UWV}qxr^-Y2XlGEcxg`Oedmiy@X0+naejc&?AIHvckVn*dwvJ_Q=Mj6Af@`!hc|>-AY(+82HUECi z9(3UmZRHxjC%t%tyZNkzKgpBVl$^|{<`Ey34h=sKyz=Q%rggx#&0;kGS2;uO<9&SE8xa6YX2;D7^w1z$Jc)dOVy_asu3AjCaEYY{ z7c{gNafzP8+E)%%afwcqT7@+axdgj=;af32msn84^1t(#OIRv3Y91u}i>L|doX_AA z;oa;^qjoNlOw^lpPjHFRiVHSf#av>u>E+>yTrTnDZ+hqcYg}S#?|hYpTU=syNlWM@ zk~`{krT4^hiTAse(@JBxgh;u5st@VQO)23;-r*7rGO72r{l_Jmma?|QG?7?vUuwq+ zE-^RVb@PfKE}?rZFTc8+OB{+h?^APw+=p*8dL)WVw1v}mhaKS(9hP73iRE#L^~~0U Qv))`n-eaJ`>=BpvAC47ng#Z8m literal 0 HcmV?d00001 diff --git a/test_assets/weighted_multiexp_decay/conf_1000_64bit.raw b/test_assets/weighted_multiexp_decay/conf_1000_64bit.raw new file mode 100644 index 0000000000000000000000000000000000000000..44f82964dc5e55bc8dfa6e9b4b169f7b9428ae22 GIT binary patch literal 8000 zcmWldcRbbKAIBppqmok5S4JhJw2ZP|k-I*)?y}dt$Ot8?B1vgV3L&AQ(n4hv+9HY; z$=(eSMT#_j=l9q1an3#C{eHdA=iKvt_g%T-8a`F>ZX$JTH0m%_mPJ$0c_F$_zWKcCqcABw#)ZmT{>g<>SA_w#{6A;8X+%dSrgL7qkL zMVq)_{M7o|6sQmk>8@KEBj($DO ze8H4Bx&cXWQcyBUnN&lSRgH)FTf>woIno6-2t>T3OfCl>Uqoa2}83Dcl^ z7A2vcP*S~FlLb$h<*L{Ik@dt*3xk^TZ5{}4nc?Vv(*ud$J=0eo_5kbig4?e}9#Gah zQKW9-fwmGiCuwC5Xq|6ru>I=}{!5PKr!IGB#oUxU|I{59b|^dX^W9;&?RVAtRCfp& z8x8eh2KOKFqclanH}I2SAkh!x@Y zUAk9aj0gcIqjU83h`=bg{fDzdgbgE^H*^9;aCP8mHTsB<_NK|lR3yS~F{fT12N8z+ zp4d@b5u!jFf6o+Q>%Yy<3^$5!EpuZJQ(puLcKDB}+9F(5cKEPku?UNv?{>^m5<$Vv zzHFDA2uIXBnSDas|0=|FpWMjyULhu36nn9uNr=R+Wrnm`AylGY zJd%1OgrrT0$KFCAjBjk8b01qAMAg9;PS-Sr$l4^f%TZAX$G3}4HcJVy z%d<|A{@o3sV!3`H^=@cM%++HSx}oyk-h#43H*DIGJf$_r4RYhVV+y!#cyAK;N_Vvz z%+Ec~51HYH8CQ>%dJMQC^Qxzz@-tU>pU}}s%W(y-HgBjS#ua);r`CNJx#C|=&g^xsX+_`~+_lq+_ILV@fW?F463*qhYZY<0%+=K4h^ zDQEazx}Vml;fy5<$G2{n?2Im#1l^JjC+MjDZMMGage~srCNobvA&D6q`aRSMhS|>U z|9DR5j#ZXfx7rD`?r!6=)1A<=*Kjd`45o*G` zSDhY?IFVjD_nMg_Obc98@2WY%rYmc=oR}kGKMWfcb~-@*#H8*E4;^5ztgfW@oCDUp zcHbEm?EpiQsQLEp4(K1dE1$B-0rtD}ZfGtc*s9^2i;_jAJRp>r-KZNH;E=xbjrTL639;w{x!u-G1gC-z77N!g<^ zF;Sj3WQTj-70<=JwZnwutj%lg+Trz~5z0K<4khM;W)F_q;qHpaqnCs1pe>*;v_pHU)Nl$+HVV$VOh3RgDs{q_cRJ0 z*`hdDq_`y47Sld`l5#k03wiqSLq8ALVy|2GRPkV2_E+yVq@_Ws3_Ayh))-wkUZj-!&>_3p?qEj`ndIT+cuAY~^PgG-f3`bT!*RHFV|m zpU-XZZr0kAZTD;tUZkBu<=bG|vTw8YoVS5(bJzH+Bpc+8)@>0xWP{d|w{Ba7+u+*D znx?@38ywxrDCrZ~;O?PYmke!f;2IXWxt${W+-0tR3mfp{BZ?GtZ6LaEtvYYH4a$q` z!V*<&Kx2n$KbdO-k4%fly3=i7@8TLbI@ty)-Fr@t{t}?tenezEEP&}k&FJTS0({V1 zHyYn1K&yVXd1{jYD``66jjshz`h2h0_L%^K`@0ln$^|(1u66ilnE<;pS3Q|iB!Fwx z#m%010z~w^9O=yw;CqHmC$EB)Pk*5S$m4C_mR=faS$@z6tjtY>N zb6G4pMgShoT2nKMjPE|;{eG7Ks&etGDnkTl+thKZdYb^u^j#;qw+KMxm2Fe<5?~v@ zFC$1K0AqNl>7BCxb-G61Z0rSy-L;0-#}^>@Wu*FHYXP>K&(P4f5Mc5h7dJ_U0NTGV z1->*B;I8`CIYlM{{C7MdGs8fDAlI++59^WkkZkwX@?oNVNw54DA1;Z_hUb6q;W(^zIAWL&@et)x zR)c&fNu0SV{h5z%KdT1rzUO1CN4j0q#mA3;IAuubdaz51#RHD5+ROqLL5a{c+BQkN6m>Xt#E~%g5&1=~w)3^HH8xkSDsuhxF}( zb%I#z+499hdJL^07+6Q@cio56drE z;#Dj8kl3n!^tJ{c^D572hpF&!Vd!G#qWOHJ+HiN2$nnv>`~3nFDL&pED%*QqjE^9h z)Q{uiJXA4L*RLApLCdf2rP*g5F3wu7A?W5o>Qw0q(`FuuoTU!Uf5U^$>Kn<8)jUYY znh!)i;^DDSfi|mz2Tsb>XGd@FU{W*t(m)OmK|ISDdS`hszo~R=TQUzNR%I!PaXggY z>-5Wv;o({2^7lt0cyMcAwrmUH;mD<~IR?HwP(7ZCV{SamyoaJ>TOKkQmzS@!;$f9f z3+Kuv9uC+lC(SV6!8M0^>$RGPv5gM)*EM)xE4H=t{>MYBN6i`W**w^4tUEYeiU*qu z--qV?re`f<17p3mJSFhn>aJ2Bl(??v~ z;~T2@-R44k|AvD@*SR?3^+|eb78gAm59{=%aA6Al2kba5s%_K9a`tiIeI0)IM|Y=%K)Tb_WOH zk;x^K>p3uzvyuE<#ewAB4&j@798j;99lCRigBIn)eP?qxs8D(zwkMr~R1K2@UMId^TBH5f8<+0K)-ay?tJOBw8Xgf(N1mlvL#$v@tJ6Mf$f)@@N(NZNEPLPBLpy7X z%oSggY-Wwt2j%hmS6aiqtoQW3`PSHExX?6vvNeX}#|+;Mvau^ZtYAem8+tS|x6BGQ zg4_x#RBp2I($MzDi*z>XIoiGnhuIjmGFTrS!bamrRs2CA8<*$Y5Z$$8qdZh6e_}lw zl`AE3wk%=eyIJPMWO+6od$dHA{E+)F znqLVs*?J$kEucfo`qKE;0&$NQ$xg1eKuU2}qsLtf?9kF*Gn#LKyz66^A7xv>;Izu) zhba~a_|aAKGu8snYQA)N?zMna)61LEp%$QIr^xmAS|DNat#@BtNqd)HfPugQ_nO~? z7qKiLKlkG$KUkniBYZMPm+ZU9r z&lC%UKZ;$Y`G>-io4aiff2U9(5gVvDKtWYJiQnBrVSm7AYf}pa!TZ5J>3RzD`KH4G zH5B$P*|1%_g2I8L;&IRKQD8jx?yo4O;3nlAGLlardwJqU!4(Qey~o*ISro4DzR6!q zr=U$$?xpZtsqXKL2nzB6xliYW zQ3$Q0x3ITUNWXqC?Up~;fBrhIsSgE{w_ncx_Mo7dwTa*9M!`ZbIb_g@0<~zin5I32 zUmK5w9~4lix10ZZ7MH??BVC~tR>a<<&(3KSg@0S zkTk(AP`DpF*JGM71%v6c8Xp@_um9st$!G@*NQ$*HZYa z6;g{er2oI_>$g`?nAW*7$U>XKi_M{flU7n#D&4xcUW-C3d$+;;09deETwQT`KQ~E1_c-YS1O@|PuMjzGYFO5?mP8p359RHHwsa^IO2R#wenLB9PseAAvp)$Zn1Z|vV(+VM-YY0#|LVK*0wst) za#NRGt0c_H?UhJgLBV^lX`wsuSMYU7>&E~4FW%drNc1@9t(7Len!;I*@y32aoyV1O z?L@EIeS1smiGEHar=CApN0EE#Q%D)nb?J)QU&Xp44!o0@xkT^YsUh__`V>N1d2h2w z92zq|=wxi5@TlmXQ7VZ~tMlmML__j><-o)qBMJ`Afv5aP{Pr$>5$4@R_~80G z=R-6Kep*tcCUgqIkSi0dggP>I+Bix^2 zmBk|L?K-)wFU>OyD<^STf7D{LIEmBvoT5+bN&Mb~ zRA00ram}3~s-=^7XXVbix1OAbNS%1jUvgeHzdbtn^>PXerI#>VG>J~TTBm;^=d?1< zm?d09&Xaj%T#+gT<%?a(6Ur3o=(`8i7f?7d`t{J-`Q#kGoBv_C02H6pBBT z7BtQv_u1OTkq2Z*oR60+pF54*51oZO5+umI@ja@mZxXo|?(L-;{b3L3eV*>@9`>}8=PAhq`6E@J1%UB~^wECeh#Y$O-J!VA%8YP1&%aSiVe zx4W^h>TJ9MNYO|U5|x+vmG8o z+AJ(un^SsdDGMeiSM66*CHrnpopxqE3;W|U2mZ*hux_F)@MBba!w^C8bGh>2NtPmekHGBNc`Wdc*kM1OyBqlO(5Mup?Ll58eyXUs;q zITImgry55ZGSU3)RMyh9Ok~fo+VNbIi6)^H!%2mSVtI{-E(IoJWD6>{$TIO)rQN(= zoC(FbnUuvi1HRK^lT7ir+BMts5>jR>eTb`-`vi z9x$+ML9R(~5ot%)`Da~Y!2I>d^7``(d~8q8_?^PQyB_!8CC3@yN3RUGiecdL3C8X% z;SAU<`(qay$bjQw-F^9942b=i@w48Ufp3%kwSVL>@L%dhO&JygwhQg%ECB;Q&P7~g z=rZ6CZYcvF);G-KgVbD7_hJPKKE-n1AEUEF_%j)kaO_KcAsB#2%7_@ z<$R?hO{25=#|JtZdR&7UZFKw#(u~QhBV2Jd=5I9}&gZ5_Y<)*Xb zbLkkf4!!gvn~q^GZP`<4bS#qRIc_*f$DV4Lfw4n$^ekR=zHAR2fh`}VrG(P4oUg4F z8Z*;G#pmjy72ob z4LNDY-fkGAL2%CV$E6Q6$Xu~7)99dqeebNxi+42K-zCG3e?^0FS({6AH4T^gC2t;i zOoJM=_|DBT8VWy6{yn9LhQudEWjpd{I8uAud(LGVRvQGmHlL?K^*c-I@fjMz=Jr`O zoT7nx_*-wm2^#i=bia!@N<(M0>Wn!DXc#(S`|L#&4M#g)ZM_pt!zt6XJ56Dh&gG(26mVEqG28u&7jdl(cA5||ex&Ls0DwbYH#Xt-XyFi+c@ zM(Q%XlTyrRa60j7#Ke?_r6)6!7H_2C-;6OebI|au(8}VX360cEY2qxxiK>l~T7>g( zB$FbXU3__7F5%CIG8t?#e)r$_nQO`XdV{rjLb6_#`scQ4vcCU#vGg9YUsGG>v~W7n zH#*X_kl6FA_+_Sn1<}3FY5r6;=`W2stHdSiCQM_63TRlQ7NV#`{Ipr|lNB8}e?Z(Qsa2^2@bZ#9qM@bIGe@-K_~p>bGdvYxinl zLMaVf>^+US*Ci^YBQ(q!e7e@* zFAXI&4}>f!I&$*1&xxHw$BEkZB4t%N)?HPO5?evX6V7=xz4de?Dklw>nbR>jXuKea zdr*0+p%AM*np3w~adzx<;nZ!Wu_Eup} z4g&?w`@SZZF-!!CjCVFC?j2*`MwR(yI4h%yf=I~GKF$CGzXQXwUEoxImxnVKh7k-RCVve@B0c?Z&0 z#5jpt;NWHJFK1LN@bidXi?Xo=O4i>ht9K-M;PD;{%D~#=$^Wc7q739op9*_JXdHY!xNqbc`3iB@9|6<05TN-m! z!DcpmZ)*(akh;Oc3+yRX*Vst&Nm5XH$%e}7cZ}p=HfqEyXPeKp#*fW`>lYeWgZ*~) zX?0g?lt}$$T9bN4XSibF^-I=x7Ier==eafF<~v;L8nTA_qqiG!R*QfeQg2`t_pfz-tgvNj)T;!Q{kj;lU{0kJdT6&;)ZT%=Qxlu z4jj2%z`;?8Ffe>QKf#~dyijLH}EFXckE zA^XZF$}{s^g~ zp4!>)p?D`BBRWcFzr~UMebZF*viVTowB=yT9X_6EQjwjtd~9M`dog?XFzs2ot8|R? ij}J^)PwLYHD|cDmTO2Uop`5oJh^&7W7v(qa2b!no8!hY%d$)9)Y sD8S8!>4(X`h+6buNkyFfv~+mfAMwxHirH0y_=pEhxmH`7#}Z9sT90$n6G=(`1tSd zJ<)#0`1;3&xCyg!X_ZR0i_?PeZuTjgt^QVGwr})07jV(34*;aJ)G~fS> zQS0S*_uE}(d5*Kjmbacg>Fj!r=Q(fG8kq2(d)0ZK>w>Z6i+}EWA70@3E*iDxEQk%% zDdstgjbDE3pB5T@iRZm+)c);7TjaoHp8G$;EKmLwp8twb`_9E@wN6+09M_Bvxn|5x`qiKJ;2NLnx^YTWH-GwdKHqi2;Tr!;L9H8n&KpLD$fN3o=Wg(MZyK|! zyjiDJM;;STF_*KqXt>s@z;yR6e)qr=b|Jw37aSg(7=?2+H~@h-f_y4^EQjjQybuhD(h z@4n%j(fecnjr*+Q1Eb@DV|+>tRY zJnSLs{LtvQ@565+o=2?rBjeknXBR|wJ!0J-8>i1TT&{F~%=$kzoY#L?t)@QVJv=cw zmK6M~t$xD$D207xrzvitfzaTME`am4cSqwo9I~G^{=> zjhC~=UdiiS8rHp)#{2o(-^<@#8rB?@#@F_eTlqO>)u-B~*p)*Qvc9Ix@Qu=*wz?{|;;&UHQ(mTtu1>kq2=+|S0r z>Z3TU=+|SBdsG~(xr>8&oyX!}^;JA>43h#74@)QFk@{JMr@bG=!|JnmXrCKDdiTe} zn$vi^Y`#GPtiDUY`xEDF@byiArE3WoJN3wO{+$W1`Y-`2XUTyp4{NT=WB=URF9f=k zht-$madT1M-GPneVd+wNqJ7B&@vI_?}uV36?G;VgJ342Wp=u!J3OC+hfTmLtCY%=F38M`XHv3>ITWX@4CS|#3#pK>;t z^OTI@WZ9u9oU0TJ5b8Rmath}w1?4lh{Whgt3g;{Zle60;O&yWKc}qbc??IEPD^fUj zDOg)D?(3;XQaFDpcuA=9)R4N2vEreb7Z z!r^K2Q#q%p$d&TEHUej#uwA z&T|@yg*(icp2oRO!+`30zMHWzs7XdpczGXVyyR+^1u0 z-KQ4M?2^v;Psd9y{{G+0ap}~7beQe*b2{}P9i1D=Nlm9Nq(fG7R#FD_Ap;{{dFIxv zDjC#?4CFTZ^Xb{`GpH9CSlFU--`OKGs2drm+B)W^*()=s9~qFgnZgL3o@xUnFw^NI(*K) zOzKW1*7i8?>ztBI>Q5G4G6gA%I+TUOz2aZ|u3;ASC=2rR-}T9&E@k0Ce|Z5})TbDhSh^(^XDHWq#)0Vtcgl?{2CxwW&YU)k98NzL?mU9+iU*=RN5 z&};L?XH(CzQ9Md2RyK7l8w18v{Ab?zZ0cJ!OuS0Yq0Z%C@>sN)|7;HRE(d0NcF3Xb z{%Ax+{puUv`6$^v=*g0e`PAopjQTzE zxuqBKsnhw$+q-<&(v$+~bpaObe@VPU0d>0o)efF1x2$6U^}7JO4-KxjY;*y2ya26_ zWRF<(Ljm=?0GE!f5ROnlT`$1ElMN~^k1M3U7ovhG?1j|%LQFY3^o!*!3#s>okaV(q zSRr-45bKI%K?|w>g{c2uqgTH_P)Hvr#F48P#()2)kbY2vE~eoY(HDv!X=O#@BKkuS zM%`WAYDK>y`a}^-xt&!+zbL|@hsELvis&0fF!AJO5&gr4bMm-}Hu{JS5fdzI^b;HU z2s&(8*+yTnVZ-n{sfX*^=r1Fh6$g&-n-QZ8-2%y z^PqWc~Z19*5V51M&(0%mQzYlM-(T{9c{rOk-S{<>`muz@E=8c>qS8ens8|n%= zZ5?B$PuVd>{&=>Xer3nOucg7;>05RbO8a@vPXDswt?_MtZr#jIAG2erpxYyz?esG{ z?o7P@xb+Y_ea()VGQNzp)8Fj)#2k1#ea?>Gr?hK-biJK^XGf-><2HNk^gTP?kaBz0 zPXDuG?u_lfw0U5s4?1vV=D0IQ6CLzJ2c9v(&Ou*vV6dR;W0f8BM+dgeDXIH<0|$N5 zfyB8x-agjOLBDjM*}U=M1|0NF2c}DT8R4LRI&erw;U4yi*W@ zgTCs(28sVi9Q0QQVwZJn_{KE{eb#|Sg6@yUIO(@eO!)qhxG^Vv*NNkDoWf4}uM?hC z6aRdpmXkj0ME4&$T|eH;Nk4XCwMYkTyEy5~PCSHT`=*ey^Fr@f~2jJyf$C}xDk7(?AslNy2%-C zH2UYyPfv|=lQ-Oubk<>ko7~}s7%rjcZt{m4qL~F^xXB@Ih~^QJ>L!o4Vd8v=n_S|7 zq_xvY9`cC?;xHl*d&nssj5%FJ{PB=iJTT?Hp@-b!L81A;9`cI^lHNM@@sMLYSbCvr z%9)WK@{9-OdSHf!T;oB_OR)m!J>(k?KDoUAgHAg<wA@6vQd9|D9Xb-u^ zgEy{K6-V-tf4rD`T?Se&ImnAEH)O!@l83yIboXX$FS*Ez!MD4&Jlot$KJsGQooWKf zyyPS=One{eB`x+F}@{|{S9w!LG z^^&W+Fvq3DOTO|!(%)OjK5~{1Ca$`C@9`03lLeB?AAOg!7+Bd_@&>G17iK60B6g-J4j z`N(fRyp=4k-cOG6!(4ae`pI*CNVj_@>L=IvA!+eKZ9n?yj%%vesZ86S7jXB;3p6IQMo{F0{rAcKZckO;3psY zv0b>)yAS>3!~l{6O^>FGwy8;)ziIOVJ?tG>CdKPSpvLQ-k>8>HeE9 zwFr_|gZNv@Rkt9yHHe}rX(DL974u4UvyS*e-E-bcmcBLXsIML*(TUn%9>K6e2f=Fhi*4 zm7O8-a|q|4g_Xvb}iG^u| z93O#+GgTwx`3U5;v1j85xjq6j?strk?;{Z6v#Dpl2su9j(|$)s$omn9;SvEpLhg@1 zZYQt*5F!6ZAl2l(T@mJh2u!;;9$_AcKom!8 z6y_HNtK@pSwZa^upoM9u3iFJDhk~xx_E4B>6wH*k-d|z9QBX&Y?{I}VN5O?1Uw$}X zjKaL5VDyW3KVLgZVeV0&mxvcTCIwC8dij>Z{G{NPTsJ;Zn4?ro z5;Q(AL1ms&QN8J}FRx2enX6Qs6x-F30+snn#c*@Mt1@S)@RKLt+ij1^oTuWQW9`IBN@dK)Nn>@r>hrg%!wL4GZ*a|^P+~3xj5FC8#VkX z*V)xI=0^>^I_;=gqn^ebsUcfzw?AyEF;8k(FXK~djk!`oI};c+=1UDRrhsb9nHuJw zmjj?NZ)#{Dw&On!*O)ssT$AgKF&guyhH;{rHNMf9Lp4+}7mgb9sD`65E-cZQOEnA? z+x16lH0DzcuC9_WHRe^Nmc~4*;i1IiCmM6Dj+r7YJWkM=Z*|mpPi_!(=3E^Y#CCshq0YRkW3>5!htAxq zLopZlI`glN10p@ttf4ap>*#lDm)J<@%)>hJB@Q;#nTvI7x-AzHI`grPP9~u1%*i_9 zO~cZemvtfsGXMWaktKX#h#uqwsw2nz~-SLgiJguX;x$xDQt96`| z?YTr}zSc2Zd^%zLs558l@X7W37M*!p#~!(^+oLmg>*y)Vc~EEm){$-oTAevu$B#1J XU(%V!b+nOoe_LlR*I|8X{uKTX*48QV literal 0 HcmV?d00001 diff --git a/test_assets/weighted_multiexp_decay/ydata_1000_64bit.raw b/test_assets/weighted_multiexp_decay/ydata_1000_64bit.raw new file mode 100644 index 0000000000000000000000000000000000000000..27c23d15a611b55e36104e5eac8e337056d99fe3 GIT binary patch literal 8000 zcmWMqcQ}=O7(R=GV;}R_t3e4VMR{)-5tWscq>^MdlmmAHexiViQ+PC_-Gnm^3Br5qsl|6>r=Jy-S@OK_k=d=o@i5> zJ+u+~l5#=aQ5)ySS6q8zt&RVRtL80RsEvo$ONN}JwGr((X0T3B8~a|xw!a_OLX(CA z4EnWT7N#_GvR(^EFGv;k-PeM&<6TQyz7}k5(78QlwJ;jX{*@Q1g`%ZKVTX5WA!2d7 zv6!nC9Inc=E0LJIu+lVlr52bgTyN*8YC(QRt35zg3)Qc)c)O`uFfDzYb!k`=3Q~(b zhz?Espx!YbeXEI8MXgsROEjU}upBnmG_h{@*v^?WP0Ujbxh@f<32s*b$K6*GNA51# zY~ZJf_AZ~~7+Xy|J6U>&L!#Zl=xEJCO}y{%Uq7U+3F87@&J0r%pVy3L)RHLKTkyAK zPy>^G{arTyX+UZ7ukNY08t@5rkMgh9fd6iV87DdMrx=`(o+p9cML=QD+gc zIzj`=`*b9F{53EyURh1TLj&2=#GOyqY2eyOhs^vn8W=zLG)Hf)2E5*cXnE*pU_-K{ z-a%PX&!yhwOsONIbaJ%|iK)#dvQIy&W9nOZdB8h$44iKJ6Z%{o{BOZ3?&a$E+`q}N zGgBQ0F2rv7n5d4Pq)(&UW7Xj%Q4-u7qKd6@nv}? zSKUk<6H%Fa(hSvM_40)eOG_Prn}!qjvej{9>f}*1L9)L7pO*&T)!;I{An{PM8n&!= ztrsp-!|IQDoa*aph>p6d+>xM$)1Q=;KJ8RP{8A$`wzV1pIYBoU7^R9v>Kia z_{cr@tBRK2>!m&VRMB^_dEuissyN(#weZPJRScB0#MWi1;y{sHI4@8Y-&(X?2P0MC zdMI#isD~=73JNn%8LGm3N3MK@vMK}<5AQoUrUI!lH{EliD%jWSfBSEb3VO=N7i3kb zK+N7NZT(k;PxbC_}=do*bu1gnld<%A|Cn5%J^lM_Dk@W666zfWqezeu-+-S zM(nu~GD6Ijbv;tT_E*2M{=2M%pv#eAXU`}>X??9|rLPhKjnk7KZ&ZT)@tNnAP(tH? zUuwCo5=zXCHb9r({GVN&s!96f3;A2?^8uA8It>f3Pt39II-(vl_EY%%jP;NgA}l0)v4`UJQXk+^*x)vQ~}4=tty~vDxl_U3{6jh#1}d2XIy!N zNYwdS|B***wsoom(@GxM$35n`tt8_MgKtFW$zxC3d9Id{ zJcM>E9+}|D!Sc<{Eq6uavF50=UUiQgibC{vz3-HRc+9uj-ZD9CF7O+(O_M|K2mdRR zSLJa1i$YCIgdAu_SshyEdxWl^|wgWumsS-hOBYgz3qi~Qr9 zvII|AuuQ&4%X`SeUFo!`ptUU4d}u58R+q(H{WJaoMzZj)+qt`&C5zZpix?laEH1LI ztQ)1s!ejnhABkZZJYjcD3$@E2b-ZKGe@|u5-*&^gxLF2HvrA1vZ_9wcV|SVJEg5`! zQ`Iw*Ap<#A`L^6cG6;RoE}}=uU}6PV{e74WvNG>0+w79Tzmyt<6Psjk!|{%&yMYWe z6}sYAI>(8%AnZw3=?RVWqPE^%W9*Z&~xgVo0C& zj_jO7HtbjnjvVzP^_x{rt8HZ?*f7X#ojn_iE~rEfJFxM`cyPAFnvMFw6ZUG$*x)^% ztZma{Bkt{}u>T@9((|3tw~DZ_N3y?Xg2u*{M6+k#howPvFt7LcDvkI#M!GlOOXJ0} zMx#y5(xCbJ4GX>^`B;8YVYM_ux`)K=A4p^Qw?&!!yV8)|7Zz=OP8yA2sYlrH(%^`^ zx}I}b8n;eY|Jiy}8d8H<4WPs_M#IeQ&Je| zSvY=vKnfwBu4>=wm4b7{O3fFQQg|IbEYbX23c2qBwjR4Hg;t-K`wy>6LB^~*zUhh- zrn>5?N;0L8a`gCk{0S+nyL8^U^t=>w*n-LrBc-r2I`QR$P$|4%EV#4ILkdTee*Wv+ zA%*i#bVagvNny#7qwDF*q)=I~W7yS13dWOk`Jj1H_%Y#p+g5{|=ho)eG!-el_h0Kc zE+++!(fY0*G%3`dP1_;Bmco;H&QhKfDI`Z$9C!aBiI6S1`OH2^cqX+=C-q8VUz|*W zS%)N`KJev8nIz`e74;>)l|)kW=(i;$lBlH^n(ru(L_clhdTXvE-b+G$MS>*sZ?Dnc ze^U}>{{+@GpOM7L3XKNggOZR{HO?>fm&CV}?2(CJN$|oN`fqKNgwu{ES5uuOp=m$d z9!a8pv2$3hlO#rNNA+B@mBe)$^>6b9-Gqj6C&qS-F;&% zMG}$^j}DvslYrx8M%uCo3F!aVv|Xk}0$MQ|LDo|eNOWKSRlZ&V1D$g#RBuZl#pjgr zz*z~bI`+59BuoN2zx*Y}10`UuY<7l4NI+%g+xg2L66j}~yYynC1lk*$N^Q*~5WLE$ zI{ddd9!gwV*!4pkf0|d(bwu{#P7zqT%MjZ^eJMAF`>sb%5_ z7T*18$V41Yr+F*Z>4@XleWQ7*d@(Gt**!IPQVd7-95w0~6GM3Thd2JcV)(7kV=n$I zhBNQWH`?73!v;^=aqCnue66cq{wH4y>>mZ(mK-sxl@)TlA1Vgds_YG}C&Un6^dMOG zG>ID%=IeNg;f?)msU2=&_^2owrRE`q80UC@2U9WZyBfL5Vuct+D#w&7&BZXB6Sy{< zBZg)Zb(!B|qL4fgvgYScQP@>2`y<{fij=UWPeeLJ5w(IId8|YfZ>tQ}ez_tFHd>{(xw)1?aJ)&rO$M`(6Qxv9ZN`c=jMd3Siz)<6$D89UXRW!6t6g>Wm zTV6{6W2f#)+QajRyv~a?vV{q`pkmHZBg;~X%@yaX12~SM6vPbCRX@33kw61 z8&)*1P;Z-eDj|=B)={cWXfq3WbB_m!b+Pc_YI}!50SoTAQwht-{5}1-r?>gAaQyW3 zAC*U0=z4E?c3~t7*|YhHIgTtO2kbH5XvG43Xm4z$DGSC6GMyxiS%?e?tZiS!fLgT!(&3 zwoL@zW7u9fg(A53tbc*uYY|jh*NTS5ilE$caT1F}gTJ}XE)gQIYo{oA9v6YtpV=V` z4-r_!- z!YE?fiPt<5#^pN3(7Zds5TD$d9B^A0&MA)WpHqafyX4Bt_*7x&J-ybuB3>BEL&43j z4hy4rW1P?C)57@kHUE#HuQ1*>t^JkaDvVJjS9-jzFlIhoi16_cM$K!BKhn#E5jU{+ zildb zgaWTg-+3|HUF+aZLm($w4QC4`_ec~B}%L zVio702yo8+HOa)mExw0MpEKdnP#9YOjEVFb@2r## zCU_%nw(YyY#MGKC<_2UuMkeNNe(LMHl?j0&SLf^2Ow78Om>O?qqV{vjV?QTS z->>t-Vht05T}tuyRGHY+WN)Q7hlwQV%O6@yN#9+@hFUfg^Itr5UoFOjlj#1Etv>~^ zAkD;M>sLXjmwV5B)F}v_RXQYUlSYfSyMAgY*CrptN-@g6$*^erSszWLk;a9J#ffcvMvj*;<)w0siq!b1Q>iAH7j zNPaEizvxN_0oXDecidksfF;x!!~0}hIcfOohehNb8;mD~778F;aL#c0I5YKz*4DA_nXLkfwSvi%mM*?7Z(kLga8s^-1;6a69DDoy{s2%t)WgqQD<|l_&B_?|_aLn^( z+xC47Nauvl>Gflvw$D^e(S`x)?f#eDmJIwInVxgmhXKm1cTb}|8IasyH@Mk=^vPU) znxVo#=$g(gm2wQѪO)yYh*8FwzA_jB}ySr_y8SwtK-prRp?x(zO%a3U~E?hkG zJMk|a|MBldR){cAPg&Q#LV$sxU(S09hUgf*mO8y}f{u)N>xUZq>Db5AkXS)uV9SP) zYx9{5MEuV4TJevLmAv`kUO(yZR(rJ0vW1Q_`afTlHPV5>qm6MjbeO&9S%?%mo;)cl z9Z#Uc>h725?z?noOu7!PPp9LeMCPIQ6?AaME*XuUrQ?a&Ea_+H^v{mdphINm9}uKt($cG_V@kvtb7{w zuzHowH_$L)Ks)BzLWBRg;ZK&=A{M zd9&~^SvO^`x|9tK_m{O=RBR-l$B*NW_M6jSC*o)0VMT-7lf#U8_B5QV3)Njh=0!XQ zcCA}NL)q8)+uE1Wu*oV_trsSQupD4-rr0^vhTBk zB_$jh%AZ|sW=YdvQqOUYP^00(yRs^`IV67Beq=&{hW*|5ZR}>KNTeyY?q$+2ru&G$ zLzae@@d@>#(^M#KSbivXkcwCliNT0}RP?`enPVeB!%F+=w`+!}nAbEN*;Y?Qn*O~8 z3R&mN-D59x`l!(SPuG3TJ1U}n{wk~;p~9f4^59Sv71z=yhn*i&v9@tp_sU8ttadae z9wBv&p(Q6*l~K|AGjQc_Ar+DSF4In}RO~D)C;CdMV7>kt^z0!O36a%qKDVjpNM(Gt zdP+rHreh&9hYDBKu*G7e|7dR0wxBaq$oafm;E+Ux%dXpF%Vgi`5Em zGO4(^Ml1K-9xA%O&zQZ9p+cchWkqBZ74tVYb7q36C`dHWHa$xtddpN*sKEU7R&I;i%-hRj=bWz3K$ z>Kd;+IG2ha2DX1X!l)>m%;jsjQqeZ2QS#P$Dtt>j?mU;J;(B~y^D%uY;-eoZo0(G4 z_xbUgbQ3DvE26}9FCgbB%_gs)R0st;KVZ3miod}F#x-;*Qog*;;VF^7hu3<057EfI z18Nf*gs2#|X`5iEQ4w4xcdlQW)D_EVET>Y*d(nzpJt9;r7}~^Pv#IcMJ6nH*)R}ur z2F5cyA zO%yC?O?eU9NkK$2Tl4%d1#zqdMfFz{NOUgtNqbGfIj4~KmLz}4XmM=(PJxTuZoB)n z6iBDrNrv}RV0!wCgL@4Hi#K_%n?Fdw{NrIBQ;#XITm5!$nS`wLb-;#eO#0B==Z)ICff$JdeVM^P0>k*xt2aZXfB}YMZU} z*_nb%H;;}m^eA|!Jhh}yo&ve?hfCC#Q?RFS@3WLe!CuSiiab57@;*dU(>JaBEHm>>l!HP70plX(q3jpkp%C$`$?*gg;@ z&+YQ~^><7PqO+N&QkwWgrOdfQ^|O3pB&6Yb0fjtY?*3OJKlsF**P#P9a{0u8^hG*l z6?~%Sn3=zN7oUh)A+Rl)PQfJS+m$De_(cBJxzDyv@`(WlKRdNnK9N+|ckyHzpE&Pu zM`)&)PrMSIL;N7~I~tn8gCFw=IbA~&hXy`T)A!6MKa)@7vRZs%68J>U?BJw}H=j@$ z+Wo=zDxWy0F&mXa`r~_l6r2p@6QSk(Ym)c!iH+V7%SS@^#F4w-?g+;5i9HHmj0b)B z#Qsw!JFH9jgr3vy{S9V(;)&q6e^U^jxOgeLrO%Jtb2j(+_ho$I(Yu7X!47<4b-bVX z+M9f0*Ke1rJ%{*&XO*9a@d7?^@5P=E2eSA?(cAgIt9xPKf7;gj7Kc6%FaA&&L?!I%N;D{^ND0b@h2j3d}2B>)U1)iCpI_!b?PDO4fN*a za6a*f&$a)>Ew14a9uG`9ezEw3W~EYs`+Xj<^7QXg;m158#^}_$C@zmU`Kh$l>^qMr z4$965|HvaW9H^Ymt2~0)bNu6t*E}MjbIG<1l{_MiWwgqyfk)_-{_qwd=P@@ceoFmE z>I{?2UWV}qxr^-Y2XlGEcxg`Oedmiy@X0+naejc&?AIHvckVn*dwvJ_Q=Mj6Af@`!hc|>-AY(+82HUECi z9(3UmZRHxjC%t%tyZNkzKgpBVl$^|{<`Ey34h=sKyz=Q%rggx#&0;kGS2;uO<9&SE8xa6YX2;D7^w1z$Jc)dOVy_asu3AjCaEYY{ z7c{gNafzP8+E)%%afwcqT7@+axdgj=;af32msn84^1t(#OIRv3Y91u}i>L|doX_AA z;oa;^qjoNlOw^lpPjHFRiVHSf#av>u>E+>yTrTnDZ+hqcYg}S#?|hYpTU=syNlWM@ zk~`{krT4^hiTAse(@JBxgh;u5st@VQO)23;-r*7rGO72r{l_Jmma?|QG?7?vUuwq+ zE-^RVb@PfKE}?rZFTc8+OB{+h?^APw+=p*8dL)WVw1v}mhaKS(9hP73iRE#L^~~0U Qv))`n-eaJ`>=BpvAC47ng#Z8m literal 0 HcmV?d00001 diff --git a/tests/integration_tests/main.rs b/tests/integration_tests/main.rs index 927510b..5df6094 100644 --- a/tests/integration_tests/main.rs +++ b/tests/integration_tests/main.rs @@ -4,7 +4,6 @@ use nalgebra::DMatrix; use nalgebra::DVector; use nalgebra::Dyn; use nalgebra::OVector; -use nalgebra::Scalar; use nalgebra::U1; use shared_test_code::evaluate_complete_model_at_params; use shared_test_code::get_double_exponential_model_with_constant_offset; @@ -17,11 +16,6 @@ use shared_test_code::models::OLearyExampleModel; use varpro::prelude::*; use varpro::solvers::levmar::*; -fn to_vector(mat: DMatrix) -> DVector { - let new_rows = Dyn(mat.nrows() * mat.ncols()); - mat.reshape_generic(new_rows, U1) -} - #[test] // sanity check my calculations above fn sanity_check_jacobian_of_levenberg_marquardt_problem_is_correct() { @@ -117,7 +111,7 @@ fn double_exponential_fitting_without_noise_produces_accurate_results() { ); let problem = LevMarProblemBuilder::new(model) - .observations(y) + .observations(y.clone()) .build() .expect("Building valid problem should not panic"); _ = format!("{:?}", problem); @@ -129,6 +123,7 @@ fn double_exponential_fitting_without_noise_produces_accurate_results() { fit_result.minimization_report.termination.was_successful(), "Levenberg Marquardt did not converge" ); + assert_relative_eq!(fit_result.best_fit().unwrap(), y, epsilon = 1e-5); assert_relative_eq!( fit_result.problem.residuals().unwrap(), statistics.weighted_residuals(), @@ -184,7 +179,7 @@ fn double_exponential_fitting_without_noise_produces_accurate_results_with_handr ); let problem = LevMarProblemBuilder::new(model) - .observations(y) + .observations(y.clone()) .build() .expect("Building valid problem should not panic"); @@ -198,6 +193,8 @@ fn double_exponential_fitting_without_noise_produces_accurate_results_with_handr epsilon = 1e-5 ); + assert_relative_eq!(fit_result.best_fit().unwrap(), y, epsilon = 1e-5); + // extract the calculated paramters, because tau1 and tau2 might switch places here let (tau1_index, tau2_index) = if fit_result.nonlinear_parameters()[0] < fit_result.nonlinear_parameters()[1] { @@ -423,26 +420,27 @@ fn double_exponential_model_with_handrolled_model_mrhs_produces_accurate_results ); let model = DoubleExpModelWithConstantOffsetSepModel::new(x, (tau1_guess, tau2_guess)); - let problem = LevMarProblemBuilder::new(model) - .observations(Y) + let problem = LevMarProblemBuilder::mrhs(model) + .observations(Y.clone()) .build() .expect("building the lev mar problem must not fail"); - let (levenberg_marquardt_solution, report) = LevenbergMarquardt::new().minimize(problem); - assert!( - report.termination.was_successful(), - "Levenberg Marquardt did not converge" - ); + let fit_result = LevMarSolver::default() + .fit(problem) + .expect("fitting must not fail"); + + assert_relative_eq!(fit_result.best_fit().unwrap(), Y, epsilon = 1e-5); + // extract the calculated paramters, because tau1 and tau2 might switch places here let (tau1_index, tau2_index) = - if levenberg_marquardt_solution.params()[0] < levenberg_marquardt_solution.params()[1] { + if fit_result.nonlinear_parameters()[0] < fit_result.nonlinear_parameters()[1] { (0, 1) } else { (1, 0) }; - let tau1_calc = levenberg_marquardt_solution.params()[tau1_index]; - let tau2_calc = levenberg_marquardt_solution.params()[tau2_index]; - let coeff = levenberg_marquardt_solution + let tau1_calc = fit_result.nonlinear_parameters()[tau1_index]; + let tau2_calc = fit_result.nonlinear_parameters()[tau2_index]; + let coeff = fit_result .linear_coefficients() .expect("linear coefficients must exist"); let a1_calc = coeff[tau1_index]; @@ -462,6 +460,164 @@ fn double_exponential_model_with_handrolled_model_mrhs_produces_accurate_results assert_relative_eq!(tau2, tau2_calc, epsilon = 1e-8); } +#[test] +fn double_exponential_model_with_noise_gives_same_confidence_interval_as_lmfit() { + // I have python scripts using the lmfit package that allow me to test + // my results. + // this tests against the file python/multiexp_decay.py + // see there for more details. The parameters are taken from there. + + let x = read_vec_f64( + "test_assets/multiexp_decay/xdata_1000_64bit.raw", + Some(1000), + ); + let y = read_vec_f64( + "test_assets/multiexp_decay/ydata_1000_64bit.raw", + Some(1000), + ); + let conf_radius = read_vec_f64("test_assets/multiexp_decay/conf_1000_64bit.raw", Some(1000)); + let covmat = read_vec_f64("test_assets/multiexp_decay/covmat_5x5_64bit.raw", Some(25)); + let model = DoubleExpModelWithConstantOffsetSepModel::new(DVector::from_vec(x), (1., 7.)); + let problem = LevMarProblemBuilder::new(model) + .observations(DVector::from_vec(y)) + .build() + .expect("building the lev mar problem must not fail"); + + let (fit_result, fit_stat) = LevMarSolver::default() + .fit_with_statistics(problem) + .expect("fitting must not fail"); + + // extract the calculated paramters, because tau1 and tau2 might switch places here + let tau1_calc = fit_result.nonlinear_parameters()[0]; + let tau2_calc = fit_result.nonlinear_parameters()[1]; + let coeff = fit_result + .linear_coefficients() + .expect("linear coefficients must exist"); + let a1_calc = coeff[0]; + let a2_calc = coeff[1]; + let a3_calc = coeff[2]; + + // parameters are taken from the python/multiexp_decay + // script in the root dir of this library. We compare against the + // fit results of the python lmfit library + // run the script to see the output + assert_relative_eq!(2.19344628, a1_calc, epsilon = 1e-5); + assert_relative_eq!(6.80462652, a2_calc, epsilon = 1e-5); + assert_relative_eq!(1.59995673, a3_calc, epsilon = 1e-5); + assert_relative_eq!(2.40392137, tau1_calc, epsilon = 1e-5); + assert_relative_eq!(5.99571068, tau2_calc, epsilon = 1e-5); + + assert_relative_eq!(fit_stat.reduced_chi2(), 1.0109e-4, epsilon = 1e-8); + + // covariance matrix is correct + let expected_covmat = DMatrix::from_row_slice(5, 5, &covmat); + let calculated_covmat = fit_stat.covariance_matrix(); + assert_relative_eq!(expected_covmat, calculated_covmat, epsilon = 1e-6); + + // now for the confidence intervals + assert_relative_eq!( + DVector::from_vec(conf_radius), + fit_stat.confidence_band_radius(0.88), + epsilon = 1e-6, + ); +} + +#[test] +fn weighted_double_exponential_model_with_noise_gives_same_confidence_interval_as_lmfit() { + // I have python scripts using the lmfit package that allow me to test + // my results. + // this tests against the file python/weighted_multiexp_decay.py + // see there for more details. The parameters are taken from there. + + let x = read_vec_f64( + "test_assets/weighted_multiexp_decay/xdata_1000_64bit.raw", + Some(1000), + ); + let y = DVector::from_vec(read_vec_f64( + "test_assets/weighted_multiexp_decay/ydata_1000_64bit.raw", + Some(1000), + )); + let conf_radius = read_vec_f64( + "test_assets/weighted_multiexp_decay/conf_1000_64bit.raw", + Some(1000), + ); + let covmat = read_vec_f64( + "test_assets/weighted_multiexp_decay/covmat_5x5_64bit.raw", + Some(25), + ); + let model = DoubleExpModelWithConstantOffsetSepModel::new(DVector::from_vec(x), (1., 7.)); + let problem = LevMarProblemBuilder::new(model) + .observations(y.clone()) + // in the python script we also apply these weights + .weights(y.map(|v| 1. / v.sqrt())) + .build() + .expect("building the lev mar problem must not fail"); + + let (fit_result, fit_stat) = LevMarSolver::default() + .fit_with_statistics(problem) + .expect("fitting must not fail"); + + // extract the calculated paramters, because tau1 and tau2 might switch places here + let tau1_calc = fit_result.nonlinear_parameters()[0]; + let tau2_calc = fit_result.nonlinear_parameters()[1]; + let coeff = fit_result + .linear_coefficients() + .expect("linear coefficients must exist"); + let a1_calc = coeff[0]; + let a2_calc = coeff[1]; + let a3_calc = coeff[2]; + + // parameters are taken from the python/weighted_multiexp_decay + // script in the root dir of this library. We compare against the + // fit results of the python lmfit library + // run the script to see the output + assert_relative_eq!(2.24275841, a1_calc, epsilon = 1e-5); + assert_relative_eq!(6.75609070, a2_calc, epsilon = 1e-5); + assert_relative_eq!(1.59790007, a3_calc, epsilon = 1e-5); + assert_relative_eq!(2.43119160, tau1_calc, epsilon = 1e-5); + assert_relative_eq!(6.02052311, tau2_calc, epsilon = 1e-5); + + assert_relative_eq!(fit_stat.reduced_chi2(), 3.2117e-5, epsilon = 1e-8); + assert_relative_eq!( + fit_stat.reduced_chi2().sqrt(), + fit_stat.regression_standard_error(), + epsilon = 1e-8 + ); + + // covariance matrix is correct + let expected_covmat = DMatrix::from_row_slice(5, 5, &covmat); + let calculated_covmat = fit_stat.covariance_matrix(); + assert_relative_eq!(expected_covmat, calculated_covmat, epsilon = 1e-6); + + // now for the confidence intervals + assert_relative_eq!( + DVector::from_vec(conf_radius), + fit_stat.confidence_band_radius(0.88), + epsilon = 1e-6, + ); +} + +// helper function to read a vector of f64 from a file +fn read_vec_f64(path: impl AsRef, size_hint: Option) -> Vec { + use byteorder::{LittleEndian, ReadBytesExt}; + let mut vect = Vec::with_capacity(size_hint.unwrap_or(1024)); + + let f = std::fs::File::open(path).expect("error opening file"); + let mut r = std::io::BufReader::new(f); + + loop { + match r.read_f64::() { + Ok(val) => { + vect.push(val); + } + Err(e) if e.kind() == std::io::ErrorKind::UnexpectedEof => break, + Err(_e) => panic!("error parsing file"), + } + } + + vect +} + #[test] // this also tests the correct application of weights fn oleary_example_with_handrolled_model_produces_correct_results() { @@ -481,7 +637,7 @@ fn oleary_example_with_handrolled_model_produces_correct_results() { let model = OLearyExampleModel::new(t, initial_guess); let problem = LevMarProblemBuilder::new(model) - .observations(y) + .observations(y.clone()) .weights(w) .build() .unwrap(); @@ -493,6 +649,9 @@ fn oleary_example_with_handrolled_model_produces_correct_results() { fit_result.minimization_report.termination.was_successful(), "fitting did not terminate successfully" ); + + assert_relative_eq!(fit_result.best_fit().unwrap(), y, epsilon = 1e-2); + let alpha_fit = fit_result.nonlinear_parameters(); let c_fit = fit_result .linear_coefficients() @@ -503,7 +662,7 @@ fn oleary_example_with_handrolled_model_produces_correct_results() { OVector::::from_vec(vec![1.0132255e+00, 2.4968675e+00, 4.0625148e+00]); let c_true = OVector::::from_vec(vec![5.8416357e+00, 1.1436854e+00]); assert_relative_eq!(alpha_fit, alpha_true, epsilon = 1e-5); - assert_relative_eq!(to_vector(c_fit), c_true, epsilon = 1e-5); + assert_relative_eq!(c_fit.into_owned(), c_true, epsilon = 1e-5); let expected_weighted_residuals = DVector::from_column_slice(&[ -1.1211e-03, @@ -568,7 +727,7 @@ fn oleary_example_with_handrolled_model_produces_correct_results() { -0.9914, 0.9918, 0.1103, 0.7729, 1.0000; ]; assert_relative_eq!( - statistics.correlation_matrix(), + statistics.calculate_correlation_matrix(), &expected_correlation_matrix, epsilon = 1e-4 ); @@ -593,7 +752,7 @@ fn test_oleary_example_with_separable_model() { let model = o_leary_example_model(t, initial_guess); let problem = LevMarProblemBuilder::new(model) - .observations(y) + .observations(y.clone()) .weights(w) .build() .unwrap(); @@ -605,6 +764,9 @@ fn test_oleary_example_with_separable_model() { fit_result.minimization_report.termination.was_successful(), "fitting did not terminate successfully" ); + + assert_relative_eq!(fit_result.best_fit().unwrap(), y, epsilon = 1e-2); + let alpha_fit = fit_result.nonlinear_parameters(); let c_fit = fit_result .linear_coefficients() @@ -614,7 +776,7 @@ fn test_oleary_example_with_separable_model() { let alpha_true = DVector::from_vec(vec![1.0132255e+00, 2.4968675e+00, 4.0625148e+00]); let c_true = DVector::from_vec(vec![5.8416357e+00, 1.1436854e+00]); assert_relative_eq!(alpha_fit, alpha_true, epsilon = 1e-5); - assert_relative_eq!(to_vector(c_fit), &c_true, epsilon = 1e-5); + assert_relative_eq!(c_fit.into_owned(), &c_true, epsilon = 1e-5); let expected_weighted_residuals = DVector::from_column_slice(&[ -1.1211e-03, @@ -704,7 +866,7 @@ fn test_oleary_example_with_separable_model() { ], ); assert_relative_eq!( - statistics.correlation_matrix(), + statistics.calculate_correlation_matrix(), &expected_correlation_matrix, epsilon = 1e-4 );