19345778652

Committed 13 Nov 2025 09:04PM UTC coverage: 86.555% (+10.2%) from 76.356%

Build # 19345778652

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #53 from geo-ant/feature/nalgebra-lapack-qr-calculations-refactor

add QR and column-pivoted QR based logic varpro using nalgebra-lapack, resulting in much improved performance in the single RHS case.

Run Details

182 of 231 new or added lines in 10 files covered. (78.79%)

8 existing lines in 2 files now uncovered.

1030 of 1190 relevant lines covered (86.55%)

3692.73 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

63.46

/shared_test_code/src/lib.rs

#![warn(missing_docs)]
//! a helper crate which carries common code used by the benchtests and the
//! integration tests.
use std::ops::Range;

use levenberg_marquardt::LeastSquaresProblem;
use nalgebra::{
    ComplexField, Const, DMatrix, DVector, DefaultAllocator, Dyn, OMatrix, OVector, Owned, Scalar,
};
use num_traits::Float;
use rand::{Rng, SeedableRng};
use varpro::model::builder::SeparableModelBuilder;
use varpro::model::SeparableModel;
use varpro::prelude::SeparableNonlinearModel;
use varpro::problem::{RhsType, SeparableProblem};
use varpro::solvers::levmar::{LevMarProblem, LevMarSolver, LinearSolver};

/// multiple right hand sides for for levenberg marquardt
pub mod levmar_mrhs;
/// contains models both for the levmar crate as well as the
/// varpro crate
pub mod models;

/// create holding `count` the elements from range [first,last] with linear spacing. (equivalent to matlabs linspace)
pub fn linspace<ScalarType: Float + Scalar>(
    first: ScalarType,
    last: ScalarType,
    count: usize,
) -> DVector<ScalarType> {
    let n_minus_one = ScalarType::from(count - 1).expect("Could not convert usize to Float");
    let lin: Vec<ScalarType> = (0..count)
        .map(|n| {
            first
                + (first - last) / (n_minus_one)
                    * ScalarType::from(n).expect("Could not convert usize to Float")
        })
        .collect();
    DVector::from(lin)
}

/// create a random matrix with coefficients in the given range
pub fn create_random_dmatrix(
    (rows, cols): (usize, usize),
    seed: u64,
    range: Range<f64>,
) -> nalgebra::DMatrix<f64> {
    let mut rng = rand::rngs::StdRng::seed_from_u64(seed);
    nalgebra::DMatrix::from_fn(rows, cols, move |_, _| rng.gen_range(range.clone()))
}

/// evaluete the vector valued function of a model by evaluating the model at the given location
/// `x` with (nonlinear) parameters `params` and by calculating the linear superposition of the basisfunctions
/// with the given linear coefficients `linear_coeffs`.
pub fn evaluate_complete_model_at_params<Model>(
    model: &'_ mut Model,
    params: OVector<Model::ScalarType, Dyn>,
    linear_coeffs: &OVector<Model::ScalarType, Dyn>,
) -> OVector<Model::ScalarType, Dyn>
where
    Model::ScalarType: Scalar + ComplexField,
    Model: SeparableNonlinearModel,
    DefaultAllocator: nalgebra::allocator::Allocator<Dyn, Dyn>,
    DefaultAllocator: nalgebra::allocator::Allocator<Dyn>,
{
    let original_params = model.params();
    model
        .set_params(params)
        .expect("Setting params must not fail");
    let eval = (&model
        .eval()
        .expect("Evaluating model must not produce error"))
        * linear_coeffs;
    model
        .set_params(original_params)
        .expect("Setting params must not fail");
    eval
}

/// evaluate the model for multiple right hand sides
pub fn evaluate_complete_model_at_params_mrhs<Model>(
    model: &'_ mut Model,
    params: OVector<Model::ScalarType, Dyn>,
    linear_coeffs: &OMatrix<Model::ScalarType, Dyn, Dyn>,
) -> OMatrix<Model::ScalarType, Dyn, Dyn>
where
    Model::ScalarType: Scalar + ComplexField,
    Model: SeparableNonlinearModel,
    DefaultAllocator: nalgebra::allocator::Allocator<Dyn, Dyn>,
    DefaultAllocator: nalgebra::allocator::Allocator<Dyn>,
{
    let original_params = model.params();
    model
        .set_params(params)
        .expect("Setting params must not fail");
    let eval = (&model
        .eval()
        .expect("Evaluating model must not produce error"))
        * linear_coeffs;
    model
        .set_params(original_params)
        .expect("Setting params must not fail");
    eval
}

/// exponential decay f(t,tau) = exp(-t/tau)
pub fn exp_decay<ScalarType: Float + Scalar>(
    tvec: &DVector<ScalarType>,
    tau: ScalarType,
) -> DVector<ScalarType> {
    tvec.map(|t| (-t / tau).exp())
}

/// derivative of exp decay with respect to tau
pub fn exp_decay_dtau<ScalarType: Scalar + Float>(
    tvec: &DVector<ScalarType>,
    tau: ScalarType,
) -> DVector<ScalarType> {
    tvec.map(|t| (-t / tau).exp() * t / (tau * tau))
}

/// A helper function that returns a double exponential decay model
/// f(x,tau1,tau2) = c1*exp(-x/tau1)+c2*exp(-x/tau2)+c3
/// Model parameters are: tau1, tau2
pub fn get_double_exponential_model_with_constant_offset(
    x: DVector<f64>,
    initial_params: Vec<f64>,
) -> SeparableModel<f64> {
    let ones = |t: &DVector<_>| DVector::from_element(t.len(), 1.);

    SeparableModelBuilder::new(&["tau1", "tau2"])
        .initial_parameters(initial_params)
        .function(&["tau1"], exp_decay)
        .partial_deriv("tau1", exp_decay_dtau)
        .function(&["tau2"], exp_decay)
        .partial_deriv("tau2", exp_decay_dtau)
        .invariant_function(ones)
        .independent_variable(x)
        .build()
        .expect("double exponential model builder should produce a valid model")
}

/// run the mimimization problem in a way that's pretty generic and can be
/// used in the benchmark code.
pub fn run_minimization_generic<Model, Rhs, Solver>(
    problem: SeparableProblem<Model, Rhs>,
) -> (DVector<f64>, DMatrix<f64>)
where
    Model: SeparableNonlinearModel<ScalarType = f64> + std::fmt::Debug,
    Solver: LinearSolver<ScalarType = f64>,
    LevMarProblem<Model, Rhs, Solver>: LeastSquaresProblem<
        Model::ScalarType,
        Dyn,
        Dyn,
        ParameterStorage = Owned<Model::ScalarType, Dyn, Const<1>>,
    >,
    Rhs: RhsType,
{
    let problem = LevMarProblem::<_, _, Solver>::from(problem);
    let result = LevMarSolver::default()
        .solve_generic::<_, Solver>(problem)
        .expect("fitting must exit successfully");
    let params = result.nonlinear_parameters();
    let coeff = result.linear_coefficients_generic().unwrap();
    (params, coeff.into_owned())
}

1	#![warn(missing_docs)]
2	//! a helper crate which carries common code used by the benchtests and the
3	//! integration tests.
4	use std::ops::Range;
5
6	use levenberg_marquardt::LeastSquaresProblem;
7	use nalgebra::{
8	ComplexField, Const, DMatrix, DVector, DefaultAllocator, Dyn, OMatrix, OVector, Owned, Scalar,
9	};
10	use num_traits::Float;
11	use rand::{Rng, SeedableRng};
12	use varpro::model::builder::SeparableModelBuilder;
13	use varpro::model::SeparableModel;
14	use varpro::prelude::SeparableNonlinearModel;
15	use varpro::problem::{RhsType, SeparableProblem};
16	use varpro::solvers::levmar::{LevMarProblem, LevMarSolver, LinearSolver};
17
18	/// multiple right hand sides for for levenberg marquardt
19	pub mod levmar_mrhs;
20	/// contains models both for the levmar crate as well as the
21	/// varpro crate
22	pub mod models;
23
24	/// create holding `count` the elements from range [first,last] with linear spacing. (equivalent to matlabs linspace)
25	pub fn linspace<ScalarType: Float + Scalar>(	17✔
26	first: ScalarType,
27	last: ScalarType,
28	count: usize,
29	) -> DVector<ScalarType> {
30	let n_minus_one = ScalarType::from(count - 1).expect("Could not convert usize to Float");	85✔
31	let lin: Vec<ScalarType> = (0..count)	51✔
32	.map(\|n\| {	8,389✔
33	first	16,744✔
34	+ (first - last) / (n_minus_one)	16,744✔
35	* ScalarType::from(n).expect("Could not convert usize to Float")	25,116✔
36	})
37	.collect();
38	DVector::from(lin)	34✔
39	}
40
41	/// create a random matrix with coefficients in the given range
42	pub fn create_random_dmatrix(	×
43	(rows, cols): (usize, usize),
44	seed: u64,
45	range: Range<f64>,
46	) -> nalgebra::DMatrix<f64> {
47	let mut rng = rand::rngs::StdRng::seed_from_u64(seed);	×
48	nalgebra::DMatrix::from_fn(rows, cols, move \|_, _\| rng.gen_range(range.clone()))	×
49	}
50
51	/// evaluete the vector valued function of a model by evaluating the model at the given location
52	/// `x` with (nonlinear) parameters `params` and by calculating the linear superposition of the basisfunctions
53	/// with the given linear coefficients `linear_coeffs`.
54	pub fn evaluate_complete_model_at_params<Model>(	8✔
55	model: &'_ mut Model,
56	params: OVector<Model::ScalarType, Dyn>,
57	linear_coeffs: &OVector<Model::ScalarType, Dyn>,
58	) -> OVector<Model::ScalarType, Dyn>
59	where
60	Model::ScalarType: Scalar + ComplexField,
61	Model: SeparableNonlinearModel,
62	DefaultAllocator: nalgebra::allocator::Allocator<Dyn, Dyn>,
63	DefaultAllocator: nalgebra::allocator::Allocator<Dyn>,
64	{
65	let original_params = model.params();	24✔
66	model	8✔
67	.set_params(params)	16✔
68	.expect("Setting params must not fail");
69	let eval = (&model	24✔
70	.eval()	16✔
71	.expect("Evaluating model must not produce error"))	8✔
72	* linear_coeffs;	8✔
73	model	8✔
74	.set_params(original_params)	16✔
75	.expect("Setting params must not fail");
76	eval	8✔
77	}
78
79	/// evaluate the model for multiple right hand sides
80	pub fn evaluate_complete_model_at_params_mrhs<Model>(
81	model: &'_ mut Model,
82	params: OVector<Model::ScalarType, Dyn>,
83	linear_coeffs: &OMatrix<Model::ScalarType, Dyn, Dyn>,
84	) -> OMatrix<Model::ScalarType, Dyn, Dyn>
85	where
86	Model::ScalarType: Scalar + ComplexField,
87	Model: SeparableNonlinearModel,
88	DefaultAllocator: nalgebra::allocator::Allocator<Dyn, Dyn>,
89	DefaultAllocator: nalgebra::allocator::Allocator<Dyn>,
90	{
91	let original_params = model.params();	×
92	model	×
93	.set_params(params)	×
94	.expect("Setting params must not fail");
95	let eval = (&model	×
96	.eval()	×
97	.expect("Evaluating model must not produce error"))	×
98	* linear_coeffs;	×
99	model	×
100	.set_params(original_params)	×
101	.expect("Setting params must not fail");
102	eval	×
103	}
104
105	/// exponential decay f(t,tau) = exp(-t/tau)
106	pub fn exp_decay<ScalarType: Float + Scalar>(	176✔
107	tvec: &DVector<ScalarType>,
108	tau: ScalarType,
109	) -> DVector<ScalarType> {
110	tvec.map(\|t\| (-t / tau).exp())	360,800✔
111	}
112
113	/// derivative of exp decay with respect to tau
114	pub fn exp_decay_dtau<ScalarType: Scalar + Float>(	98✔
115	tvec: &DVector<ScalarType>,
116	tau: ScalarType,
117	) -> DVector<ScalarType> {
118	tvec.map(\|t\| (-t / tau).exp() * t / (tau * tau))	301,252✔
119	}
120
121	/// A helper function that returns a double exponential decay model
122	/// f(x,tau1,tau2) = c1exp(-x/tau1)+c2exp(-x/tau2)+c3
123	/// Model parameters are: tau1, tau2
124	pub fn get_double_exponential_model_with_constant_offset(	5✔
125	x: DVector<f64>,
126	initial_params: Vec<f64>,
127	) -> SeparableModel<f64> {
128	let ones = \|t: &DVector<_>\| DVector::from_element(t.len(), 1.);	269✔
129
130	SeparableModelBuilder::new(&["tau1", "tau2"])	10✔
131	.initial_parameters(initial_params)	10✔
132	.function(&["tau1"], exp_decay)	10✔
133	.partial_deriv("tau1", exp_decay_dtau)	5✔
134	.function(&["tau2"], exp_decay)	10✔
135	.partial_deriv("tau2", exp_decay_dtau)	5✔
136	.invariant_function(ones)	10✔
137	.independent_variable(x)	10✔
138	.build()
139	.expect("double exponential model builder should produce a valid model")
140	}
141
142	/// run the mimimization problem in a way that's pretty generic and can be
143	/// used in the benchmark code.
144	pub fn run_minimization_generic<Model, Rhs, Solver>(
145	problem: SeparableProblem<Model, Rhs>,
146	) -> (DVector<f64>, DMatrix<f64>)
147	where
148	Model: SeparableNonlinearModel<ScalarType = f64> + std::fmt::Debug,
149	Solver: LinearSolver<ScalarType = f64>,
150	LevMarProblem<Model, Rhs, Solver>: LeastSquaresProblem<
151	Model::ScalarType,
152	Dyn,
153	Dyn,
154	ParameterStorage = Owned<Model::ScalarType, Dyn, Const<1>>,
155	>,
156	Rhs: RhsType,
157	{
NEW 158	let problem = LevMarProblem::<_, _, Solver>::from(problem);	×
NEW 159	let result = LevMarSolver::default()	×
NEW 160	.solve_generic::<_, Solver>(problem)	×
161	.expect("fitting must exit successfully");
NEW 162	let params = result.nonlinear_parameters();	×
NEW 163	let coeff = result.linear_coefficients_generic().unwrap();	×
NEW 164	(params, coeff.into_owned())	×
165	}

geo-ant / varpro / 19345778652

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous