Skip to content

Commit f2adeec

Browse files
geo-antgeo-ant
authored
Feature/mrhs v2 (#26)
* get rid of OutputDim assoc type * start preparting removal of parameterdim * more preparations to get rid of parameter dim * make model dim dynamic * remove ModelDim and ParameterDim, doctest still borked * fix doctests * fix lints * simplify trait bounds * simplify more trait bounds * fix compile * add rayon, rewrite jacobian calc to functional style * remove rayon dependency for now but keep functional style jacobian * remove deprecated stuff * remove dead code * add matrix vectorization method * minor thing * multiple right hand sides, but not tested and not everything will work * get tests working * more tests * add test * start adding levmar crate mhrs test * finish mrhs levmar problem, untested * some cleanup * fix some bugs, but jacobian is not correct yet... * more fixes, still not correct * fix jacobian calculation * add test for conventional levmar mrhs * test for my mrhs model, but stuff still fails * fix mrhs calculation * start adding benches * benches, results are troubling: handrolled slower * some clippy fixes * benchmarking, handcrafted suspiciously slower * okay, it's a phantom * lints * more docs, better handling of observations * more consistent linear coefficient handling * update docs * remove unused imports * address lints * update changelog * add crates.io badge and link * address lints * more lints
1 parent 8d86def commit f2adeec

File tree

21 files changed

+1047
-795
lines changed

21 files changed

+1047
-795
lines changed

‎CHANGES.md

+11
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,17 @@
33
This is the changelog for the `varpro` library.
44
See also here for a [version history](https://crates.io/crates/varpro/versions).
55

6+
## 0.8.0 Multiple Right Hand Sides
7+
8+
- Observations can now be a matrix, in this case a global fit with multiple right
9+
hand sides is performed.
10+
- Linear coefficients are now in matrix form. For a single right hand side, this
11+
is still a matrix with one column, for multiple right hand sides these are now
12+
the coefficient vectors (columns) corresponding to the respective right hand
13+
side column of the observation matrix.
14+
- Removed lots of generics internally.
15+
- Deprecated the new function of the levmar solver.
16+
617
## 0.7.0 Fit Statistics
718

819
- Deprecate the `minimize` function of the LevMarSolver and

‎Cargo.toml

+8-2
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
[package]
22
name = "varpro"
3-
version = "0.7.1"
3+
version = "0.8.0"
44
authors = ["geo-ant"]
55
edition = "2021"
66
license = "MIT"
@@ -17,8 +17,9 @@ members = ["shared_test_code"]
1717
[dependencies]
1818
thiserror = "1"
1919
levenberg-marquardt = "0.13"
20-
nalgebra = "0.32"
20+
nalgebra = { version = "0.32" } #, features = ["rayon"]}
2121
num-traits = "0.2"
22+
# rayon = "1.6"
2223

2324
[dev-dependencies]
2425
approx = "0.5"
@@ -28,11 +29,16 @@ pprof = { version = "0.13", features = ["criterion", "flamegraph"] }
2829
shared_test_code = { path = "./shared_test_code" }
2930
assert_matches = "1.5"
3031
mockall = "0.11"
32+
rand = "0.8"
3133

3234
[[bench]]
3335
name = "double_exponential_without_noise"
3436
harness = false
3537

38+
[[bench]]
39+
name = "multiple_right_hand_sides"
40+
harness = false
41+
3642
[package.metadata.docs.rs]
3743
# To build locally use
3844
# RUSTDOCFLAGS="--html-in-header katex-header.html" cargo doc --no-deps --open

‎README.md

+14
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
![build](https://github.com/geo-ant/varpro/workflows/build/badge.svg?branch=main)
44
![tests](https://github.com/geo-ant/varpro/workflows/tests/badge.svg?branch=main)
55
![lints](https://github.com/geo-ant/varpro/workflows/lints/badge.svg?branch=main)
6+
[![crates](https://img.shields.io/crates/v/varpro)](https://crates.io/crates/varpro)
67
[![Coverage Status](https://coveralls.io/repos/github/geo-ant/varpro/badge.svg?branch=main)](https://coveralls.io/github/geo-ant/varpro?branch=main)
78
![maintenance-status](https://img.shields.io/badge/maintenance-actively--developed-brightgreen.svg)
89

@@ -106,6 +107,19 @@ Additionally to the `fit` member function, the `LevMarSolver` also provides a
106107
`fit_with_statistics` function that calculates some additional statistical
107108
information after the fit has finished.
108109

110+
### Global Fitting of Multiple Right Hand Sides
111+
112+
Before, we have passed a single column vector as the observations. It is also
113+
possible to pass a matrix, whose columns represent individual observations. We
114+
are now fitting a problem with multiple right hand sides. `vapro` will performa a _global fit_
115+
in which the nonlinear parameters are optimized across all right hand sides, but
116+
linear parameters of the fit are optimized for each right hand side individually.
117+
118+
This is an application where varpro really shines because it can take advantage
119+
of the separable nature of the problem. It allows us to perform a global fit over thousands
120+
or even tens of thousands of right hand sides in reasonable time (fractions of seconds to minutes),
121+
where conventional purely nonlinear solvers must perform much more work.
122+
109123
### Maximum Performance
110124

111125
The example code above will already run many times faster

‎benches/double_exponential_without_noise.rs

+16-76
Original file line numberDiff line numberDiff line change
@@ -3,12 +3,9 @@ use levenberg_marquardt::LeastSquaresProblem;
33
use levenberg_marquardt::LevenbergMarquardt;
44
use nalgebra::ComplexField;
55

6-
use nalgebra::Const;
76
use nalgebra::DefaultAllocator;
87

9-
use nalgebra::DimMin;
10-
use nalgebra::DimSub;
11-
8+
use nalgebra::Dyn;
129
use nalgebra::OVector;
1310
use nalgebra::RawStorageMut;
1411

@@ -19,7 +16,7 @@ use pprof::criterion::{Output, PProfProfiler};
1916
use shared_test_code::models::DoubleExpModelWithConstantOffsetSepModel;
2017
use shared_test_code::models::DoubleExponentialDecayFittingWithOffsetLevmar;
2118
use shared_test_code::*;
22-
use varpro::model::SeparableModel;
19+
2320
use varpro::prelude::SeparableNonlinearModel;
2421
use varpro::solvers::levmar::LevMarProblem;
2522
use varpro::solvers::levmar::LevMarProblemBuilder;
@@ -41,55 +38,12 @@ fn build_problem<Model>(
4138
) -> LevMarProblem<Model>
4239
where
4340
Model: SeparableNonlinearModel<ScalarType = f64>,
44-
DefaultAllocator: nalgebra::allocator::Allocator<f64, Model::ParameterDim>,
45-
DefaultAllocator: nalgebra::allocator::Allocator<f64, Model::ParameterDim, Model::OutputDim>,
46-
DefaultAllocator: nalgebra::allocator::Allocator<f64, Model::OutputDim, Model::ModelDim>,
47-
DefaultAllocator: nalgebra::allocator::Allocator<f64, Model::ModelDim>,
48-
DefaultAllocator: nalgebra::allocator::Allocator<f64, Model::OutputDim>,
49-
<DefaultAllocator as nalgebra::allocator::Allocator<f64, Model::OutputDim>>::Buffer:
50-
Storage<f64, Model::OutputDim>,
51-
<DefaultAllocator as nalgebra::allocator::Allocator<f64, Model::OutputDim>>::Buffer:
52-
RawStorageMut<f64, Model::OutputDim>,
53-
DefaultAllocator: nalgebra::allocator::Allocator<f64, Model::OutputDim, Model::ParameterDim>,
41+
DefaultAllocator: nalgebra::allocator::Allocator<f64, Dyn>,
42+
DefaultAllocator: nalgebra::allocator::Allocator<f64, Dyn, Dyn>,
43+
<DefaultAllocator as nalgebra::allocator::Allocator<f64, Dyn>>::Buffer: Storage<f64, Dyn>,
44+
<DefaultAllocator as nalgebra::allocator::Allocator<f64, Dyn>>::Buffer: RawStorageMut<f64, Dyn>,
5445
DefaultAllocator:
55-
nalgebra::allocator::Allocator<f64, <Model::OutputDim as DimMin<Model::ModelDim>>::Output>,
56-
DefaultAllocator: nalgebra::allocator::Allocator<
57-
(usize, usize),
58-
<Model::OutputDim as DimMin<Model::ModelDim>>::Output,
59-
>,
60-
DefaultAllocator: nalgebra::allocator::Allocator<
61-
f64,
62-
<Model::OutputDim as DimMin<Model::ModelDim>>::Output,
63-
Model::OutputDim,
64-
>,
65-
DefaultAllocator: nalgebra::allocator::Allocator<
66-
<f64 as ComplexField>::RealField,
67-
<<Model::OutputDim as DimMin<Model::ModelDim>>::Output as DimSub<Const<1>>>::Output,
68-
>,
69-
DefaultAllocator: nalgebra::allocator::Allocator<
70-
f64,
71-
<<Model::OutputDim as DimMin<Model::ModelDim>>::Output as DimSub<Const<1>>>::Output,
72-
>,
73-
DefaultAllocator: nalgebra::allocator::Allocator<
74-
(<f64 as ComplexField>::RealField, usize),
75-
<Model::OutputDim as DimMin<Model::ModelDim>>::Output,
76-
>,
77-
<Model::OutputDim as DimMin<Model::ModelDim>>::Output: DimSub<nalgebra::dimension::Const<1>>,
78-
Model::OutputDim: DimMin<Model::ModelDim>,
79-
DefaultAllocator: nalgebra::allocator::Allocator<
80-
f64,
81-
<Model::OutputDim as DimMin<Model::ModelDim>>::Output,
82-
Model::ModelDim,
83-
>,
84-
DefaultAllocator: nalgebra::allocator::Allocator<
85-
f64,
86-
Model::OutputDim,
87-
<Model::OutputDim as DimMin<Model::ModelDim>>::Output,
88-
>,
89-
DefaultAllocator: nalgebra::allocator::Allocator<
90-
<f64 as ComplexField>::RealField,
91-
<Model::OutputDim as DimMin<Model::ModelDim>>::Output,
92-
>,
46+
nalgebra::allocator::Allocator<(<f64 as ComplexField>::RealField, usize), Dyn>,
9347
{
9448
let DoubleExponentialParameters {
9549
tau1,
@@ -100,39 +54,25 @@ where
10054
} = true_parameters;
10155

10256
// save the initial guess so that we can reset the model to those
103-
let params = OVector::from_vec_generic(model.parameter_count(), U1, vec![tau1, tau2]);
57+
let params = OVector::from_vec_generic(Dyn(model.parameter_count()), U1, vec![tau1, tau2]);
10458

10559
let base_function_count = model.base_function_count();
10660
let y = evaluate_complete_model_at_params(
10761
&mut model,
10862
params,
109-
&OVector::from_vec_generic(base_function_count, U1, vec![c1, c2, c3]),
63+
&OVector::from_vec_generic(Dyn(base_function_count), U1, vec![c1, c2, c3]),
11064
);
11165
LevMarProblemBuilder::new(model)
11266
.observations(y)
11367
.build()
11468
.expect("Building valid problem should not panic")
11569
}
11670

117-
/// solve the double exponential fitting problem using a handrolled model
118-
fn run_minimization_for_handrolled_separable_model(
119-
problem: LevMarProblem<DoubleExpModelWithConstantOffsetSepModel>,
120-
) -> [f64; 5] {
121-
let result = LevMarSolver::new()
122-
.fit(problem)
123-
.expect("fitting must exit successfully");
124-
let params = result.nonlinear_parameters();
125-
let coeff = result.linear_coefficients().unwrap();
126-
[params[0], params[1], coeff[0], coeff[1], coeff[2]]
127-
}
128-
129-
/// solve the double exponential fitting problem using a separable model from the builder
130-
/// I should be able to unify this with the handrolled model, but I can't figure out how to do it
131-
/// because I cannot find the correct generic bounds to do it
132-
fn run_minimization_for_builder_separable_model(
133-
problem: LevMarProblem<SeparableModel<f64>>,
134-
) -> [f64; 5] {
135-
let result = LevMarSolver::new()
71+
fn run_minimization<Model>(problem: LevMarProblem<Model>) -> [f64; 5]
72+
where
73+
Model: SeparableNonlinearModel<ScalarType = f64> + std::fmt::Debug,
74+
{
75+
let result = LevMarSolver::default()
13676
.fit(problem)
13777
.expect("fitting must exit successfully");
13878
let params = result.nonlinear_parameters();
@@ -214,7 +154,7 @@ fn bench_double_exp_no_noise(c: &mut Criterion) {
214154
),
215155
)
216156
},
217-
run_minimization_for_builder_separable_model,
157+
run_minimization,
218158
criterion::BatchSize::SmallInput,
219159
)
220160
});
@@ -227,7 +167,7 @@ fn bench_double_exp_no_noise(c: &mut Criterion) {
227167
DoubleExpModelWithConstantOffsetSepModel::new(x.clone(), tau_guess),
228168
)
229169
},
230-
run_minimization_for_handrolled_separable_model,
170+
run_minimization,
231171
criterion::BatchSize::SmallInput,
232172
)
233173
});

‎benches/multiple_right_hand_sides.rs

+104
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
use criterion::{criterion_group, criterion_main, Criterion};
2+
use nalgebra::DMatrix;
3+
use nalgebra::DVector;
4+
use nalgebra::Dyn;
5+
use nalgebra::OVector;
6+
use nalgebra::U1;
7+
use pprof::criterion::{Output, PProfProfiler};
8+
use shared_test_code::models::DoubleExpModelWithConstantOffsetSepModel;
9+
use shared_test_code::*;
10+
use varpro::prelude::SeparableNonlinearModel;
11+
use varpro::solvers::levmar::LevMarProblem;
12+
use varpro::solvers::levmar::LevMarProblemBuilder;
13+
use varpro::solvers::levmar::LevMarSolver;
14+
15+
/// helper struct for the parameters of the double exponential
16+
#[derive(Clone, PartialEq, Debug)]
17+
struct DoubleExponentialParameters {
18+
tau1: f64,
19+
tau2: f64,
20+
coeffs: DMatrix<f64>,
21+
}
22+
23+
fn build_problem_mrhs<Model>(
24+
true_parameters: DoubleExponentialParameters,
25+
mut model: Model,
26+
) -> LevMarProblem<Model>
27+
where
28+
Model: SeparableNonlinearModel<ScalarType = f64>,
29+
{
30+
let DoubleExponentialParameters { tau1, tau2, coeffs } = true_parameters.clone();
31+
// save the initial guess so that we can reset the model to those
32+
let params = OVector::from_vec_generic(Dyn(model.parameter_count()), U1, vec![tau1, tau2]);
33+
let y = evaluate_complete_model_at_params_mrhs(&mut model, params, &coeffs);
34+
LevMarProblemBuilder::new(model)
35+
.observations(y)
36+
.build()
37+
.expect("Building valid problem should not panic")
38+
}
39+
40+
fn run_minimization<Model>(problem: LevMarProblem<Model>) -> (DVector<f64>, DMatrix<f64>)
41+
where
42+
Model: SeparableNonlinearModel<ScalarType = f64> + std::fmt::Debug,
43+
{
44+
let result = LevMarSolver::default()
45+
.fit(problem)
46+
.expect("fitting must exit successfully");
47+
let params = result.nonlinear_parameters();
48+
let coeff = result.linear_coefficients().unwrap();
49+
(params, coeff)
50+
}
51+
52+
fn bench_double_exp_no_noise_mrhs(c: &mut Criterion) {
53+
// see here on comparing functions
54+
// https://bheisler.github.io/criterion.rs/book/user_guide/comparing_functions.html
55+
let mut group = c.benchmark_group("Double Exponential Without Noise");
56+
// the support points for the model (moved into the closures separately)
57+
let x = linspace(0., 12.5, 1024);
58+
// initial guess for tau
59+
let tau_guess = (2., 6.5);
60+
61+
let dataset_size = 1000;
62+
let linear_coeffs = create_random_dmatrix((3, dataset_size), 2314093240213841123, 0.0..100.0);
63+
64+
let true_parameters = DoubleExponentialParameters {
65+
tau1: 1.,
66+
tau2: 3.,
67+
coeffs: linear_coeffs,
68+
};
69+
70+
group.bench_function("Handcrafted Model (MRHS)", |bencher| {
71+
bencher.iter_batched(
72+
|| {
73+
build_problem_mrhs(
74+
true_parameters.clone(),
75+
DoubleExpModelWithConstantOffsetSepModel::new(x.clone(), tau_guess),
76+
)
77+
},
78+
run_minimization,
79+
criterion::BatchSize::SmallInput,
80+
)
81+
});
82+
83+
group.bench_function("Using Model Builder (MRHS)", |bencher| {
84+
bencher.iter_batched(
85+
|| {
86+
build_problem_mrhs(
87+
true_parameters.clone(),
88+
get_double_exponential_model_with_constant_offset(
89+
x.clone(),
90+
vec![tau_guess.0, tau_guess.1],
91+
),
92+
)
93+
},
94+
run_minimization,
95+
criterion::BatchSize::SmallInput,
96+
)
97+
});
98+
}
99+
100+
criterion_group!(
101+
name = benches_mrhs;
102+
config = Criterion::default().with_profiler(PProfProfiler::new(100, Output::Flamegraph(None)));
103+
targets = bench_double_exp_no_noise_mrhs);
104+
criterion_main!(benches_mrhs);

‎shared_test_code/Cargo.toml

+2
Original file line numberDiff line numberDiff line change
@@ -10,3 +10,5 @@ varpro = {path = "../"}
1010
nalgebra = "0.32"
1111
num-traits = "0.2"
1212
levenberg-marquardt = "0.13"
13+
approx = "0.5"
14+
rand = "0.8"

0 commit comments

Comments
 (0)