9181626839

Committed 21 May 2024 09:18PM UTC coverage: 81.737% (-1.2%) from 82.953%

Build # 9181626839

Build Type

push

github

Committed by

web-flow

Commit Message

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... (continued)

Run Details

135 of 178 new or added lines in 7 files covered. (75.84%)

9 existing lines in 3 files now uncovered.

819 of 1002 relevant lines covered (81.74%)

1259.22 hits per line

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

89.55

/shared_test_code/src/levmar_mrhs.rs

#![allow(non_snake_case)]
use levenberg_marquardt::LeastSquaresProblem;
use nalgebra::DMatrix;
use nalgebra::DVector;
use nalgebra::Dyn;
use nalgebra::Matrix;
use nalgebra::OMatrix;
use nalgebra::Owned;
use nalgebra::Vector;
use nalgebra::U1;
use nalgebra::U3;

/// double exponential model with constant offset where the nonlinear
/// parameters are shared across the dataset and the linear parameters are not
/// f_s(x) = c_{s,1} * exp(-x/alpha_1) + c_{s,2} * exp(-x/alpha_2) + c_{s,3},
/// where s is the index of the element in the dataset
pub struct DoubleExponentialModelWithConstantOffsetLevmarMrhs {
    /// the independent variable. Also defines the output length of the model.
    x: DVector<f64>,
    /// the nonlinear parameters [alpha1,alpha2]
    alpha: [f64; 2],
    /// the matrix of datasets organized as columns
    ///       |       |
    /// Y = (y_1,...,y_s,...)
    ///       |       |
    Y: DMatrix<f64>,
    /// the matrix of linear coefficients organized as columns
    ///       |       |
    /// C = (c_1,...,c_s,...)
    ///       |       |
    C: OMatrix<f64, U3, Dyn>,
    /// the precalculated function matrix
    /// where the first exponential comes first (alpha1),
    /// then the second exponential and finally
    /// a column of all ones for the constant offset
    Phi: OMatrix<f64, Dyn, U3>,
}

impl DoubleExponentialModelWithConstantOffsetLevmarMrhs {
    /// create a new fitting problem with the given data and initial parameters
    /// See the internal struct documentation for the parameter layout
    pub fn new(initial_params: impl AsRef<[f64]>, x: DVector<f64>, data: DMatrix<f64>) -> Self {
        assert_eq!(
            x.len(),
            data.nrows(),
            "x vector must have same number of rows as data"
        );
        assert!(
            data.ncols() > 0,
            "data matrix must have at least one column"
        );
        assert!(!x.is_empty(), "x vector must have at least one element");
        assert_eq!(
            data.ncols() * 3 + 2,
            initial_params.as_ref().len(),
            "we need  exactly 3 linear parameters per dataset"
        );
        let mut Phi = OMatrix::zeros_generic(Dyn(x.len()), U3);
        // the third column of only ones. We never touch this one again
        Phi.column_mut(2).copy_from_slice(&vec![1f64; x.len()]);
        let mut me = Self {
            C: OMatrix::zeros_generic(U3, Dyn(data.ncols())),
            alpha: [0., 0.], //<- this will be overwritten by set params
            x,
            Y: data,
            Phi,
        };
        me.set_params(&DVector::from_column_slice(initial_params.as_ref()));
        me
    }
}

impl DoubleExponentialModelWithConstantOffsetLevmarMrhs {
    /// get the linear parameter matrix
    pub fn lin_param_matrix(&self) -> &OMatrix<f64, U3, Dyn> {
        &self.C
    }
}

impl LeastSquaresProblem<f64, Dyn, Dyn> for DoubleExponentialModelWithConstantOffsetLevmarMrhs {
    type ParameterStorage = Owned<f64, Dyn>;
    type ResidualStorage = Owned<f64, Dyn>;
    type JacobianStorage = Owned<f64, Dyn, Dyn>;

    /// parameters are organized as
    /// alpha_1,alpha_2, c_{1,1},c_{1,2},c_{1,3},...,c_{s,1},c_{s,2},c_{s,3},...
    fn set_params(&mut self, params: &Vector<f64, Dyn, Self::ParameterStorage>) {
        debug_assert_eq!(
            params.len(),
            self.Y.ncols() * 3 + 2,
            "number of parameters does not match the dataset and the decay model"
        );
        //tau1
        self.alpha[0] = params[0];
        //tau2
        self.alpha[1] = params[1];

        // and do precalculations
        self.Phi
            .column_mut(0)
            .copy_from(&self.x.map(|x: f64| (-x / self.alpha[0]).exp()));
        self.Phi
            .column_mut(1)
            .copy_from(&self.x.map(|x: f64| (-x / self.alpha[1]).exp()));
        // and load the matrix of linear coefficients
        // I have a test down below that shows this should actually work
        self.C.copy_from_slice(&params.as_slice()[2..]);
    }

    fn params(&self) -> Vector<f64, Dyn, Self::ParameterStorage> {
        let param_count = self.C.len() + 2;
        let params = DVector::from_iterator(
            param_count,
            self.alpha
                .iter()
                .cloned()
                .chain(self.C.as_slice().iter().cloned()),
        );
        params
    }

    fn residuals(&self) -> Option<DVector<f64>> {
        let R = &self.Y - &self.Phi * &self.C;
        let new_nrows = Dyn(R.nrows() * R.ncols());
        Some(R.reshape_generic(new_nrows, U1))
    }

    fn jacobian(&self) -> Option<Matrix<f64, Dyn, Dyn, Self::JacobianStorage>> {
        // let mut jacobian = DMatrix::<f64>::zeros(self.
        let mut dPhi_dalpha1 = OMatrix::<_, Dyn, U3>::zeros_generic(Dyn(self.Phi.nrows()), U3);
        dPhi_dalpha1.set_column(
            0,
            &(-1. / (self.alpha[0] * self.alpha[0]) * self.Phi.column(0).component_mul(&self.x)),
        );
        let mut dPhi_dalpha2 = OMatrix::<_, Dyn, U3>::zeros_generic(Dyn(self.Phi.nrows()), U3);
        dPhi_dalpha2.set_column(
            1,
            &(-1. / (self.alpha[1] * self.alpha[1]) * self.Phi.column(1).component_mul(&self.x)),
        );
        let mut jac_block = DMatrix::from_element(self.Phi.nrows(), 3, -1.);
        jac_block.set_column(0, &-self.Phi.column(0));
        jac_block.set_column(1, &-self.Phi.column(1));

        let x_len = self.Y.nrows();
        let total_data_len = x_len * self.Y.ncols();
        let mut jacobian = DMatrix::zeros(self.Y.nrows() * self.Y.ncols(), self.C.len() + 2);
        debug_assert_eq!(
            self.Y.ncols() * 3,
            self.C.len(),
            "we need exactly 3 linear parameters per dataset"
        );

        let to_vector = |mat: OMatrix<f64, Dyn, Dyn>| {
            let new_nrows = Dyn(mat.nrows() * mat.ncols());
            mat.reshape_generic(new_nrows, U1)
        };
        jacobian.set_column(0, &to_vector(dPhi_dalpha1 * &self.C));
        jacobian.set_column(1, &to_vector(dPhi_dalpha2 * &self.C));

        for (col, row) in (2..jacobian.nrows())
            .step_by(3)
            .zip((0..total_data_len).step_by(x_len))
        {
            let mut view = jacobian.view_mut((row, col), (x_len, 3));
            view.copy_from(&jac_block);
        }
        Some(jacobian)
    }
}

#[test]
// make sure that the way I copy the parameters into a matrix works as
// I think it does
fn matrix_is_colum_major_layout() {
    let data = &[1., 2., 3., 4., 5., 6.];
    let mut mat = OMatrix::zeros_generic(Dyn(2), U3);
    mat.copy_from_slice(data);
    assert_eq!(mat.column(0).as_slice(), &[1., 2.]);
    assert_eq!(mat.column(1).as_slice(), &[3., 4.]);
    assert_eq!(mat.column(2).as_slice(), &[5., 6.]);
    // also assert that the reverse is true
    assert_eq!(mat.as_slice(), data);
}

1	#![allow(non_snake_case)]
2	use levenberg_marquardt::LeastSquaresProblem;
3	use nalgebra::DMatrix;
4	use nalgebra::DVector;
5	use nalgebra::Dyn;
6	use nalgebra::Matrix;
7	use nalgebra::OMatrix;
8	use nalgebra::Owned;
9	use nalgebra::Vector;
10	use nalgebra::U1;
11	use nalgebra::U3;
12
13	/// double exponential model with constant offset where the nonlinear
14	/// parameters are shared across the dataset and the linear parameters are not
15	/// f_s(x) = c_{s,1} * exp(-x/alpha_1) + c_{s,2} * exp(-x/alpha_2) + c_{s,3},
16	/// where s is the index of the element in the dataset
17	pub struct DoubleExponentialModelWithConstantOffsetLevmarMrhs {
18	/// the independent variable. Also defines the output length of the model.
19	x: DVector<f64>,
20	/// the nonlinear parameters [alpha1,alpha2]
21	alpha: [f64; 2],
22	/// the matrix of datasets organized as columns
23	/// \| \|
24	/// Y = (y_1,...,y_s,...)
25	/// \| \|
26	Y: DMatrix<f64>,
27	/// the matrix of linear coefficients organized as columns
28	/// \| \|
29	/// C = (c_1,...,c_s,...)
30	/// \| \|
31	C: OMatrix<f64, U3, Dyn>,
32	/// the precalculated function matrix
33	/// where the first exponential comes first (alpha1),
34	/// then the second exponential and finally
35	/// a column of all ones for the constant offset
36	Phi: OMatrix<f64, Dyn, U3>,
37	}
38
39	impl DoubleExponentialModelWithConstantOffsetLevmarMrhs {
40	/// create a new fitting problem with the given data and initial parameters
41	/// See the internal struct documentation for the parameter layout
42	pub fn new(initial_params: impl AsRef<[f64]>, x: DVector<f64>, data: DMatrix<f64>) -> Self {	2✔
43	assert_eq!(	2✔
44	x.len(),	2✔
45	data.nrows(),	2✔
UNCOV 46	"x vector must have same number of rows as data"	×
47	);
48	assert!(	2✔
49	data.ncols() > 0,	2✔
UNCOV 50	"data matrix must have at least one column"	×
51	);
52	assert!(!x.is_empty(), "x vector must have at least one element");	2✔
53	assert_eq!(	2✔
54	data.ncols() * 3 + 2,	2✔
55	initial_params.as_ref().len(),	2✔
UNCOV 56	"we need exactly 3 linear parameters per dataset"	×
57	);
58	let mut Phi = OMatrix::zeros_generic(Dyn(x.len()), U3);	2✔
59	// the third column of only ones. We never touch this one again
60	Phi.column_mut(2).copy_from_slice(&vec![1f64; x.len()]);	2✔
61	let mut me = Self {
62	C: OMatrix::zeros_generic(U3, Dyn(data.ncols())),	2✔
63	alpha: [0., 0.], //<- this will be overwritten by set params	2✔
64	x,
65	Y: data,
66	Phi,
67	};
68	me.set_params(&DVector::from_column_slice(initial_params.as_ref()));	2✔
69	me	2✔
70	}
71	}
72
73	impl DoubleExponentialModelWithConstantOffsetLevmarMrhs {
74	/// get the linear parameter matrix
75	pub fn lin_param_matrix(&self) -> &OMatrix<f64, U3, Dyn> {	×
76	&self.C	×
77	}
78	}
79
80	impl LeastSquaresProblem<f64, Dyn, Dyn> for DoubleExponentialModelWithConstantOffsetLevmarMrhs {
81	type ParameterStorage = Owned<f64, Dyn>;
82	type ResidualStorage = Owned<f64, Dyn>;
83	type JacobianStorage = Owned<f64, Dyn, Dyn>;
84
85	/// parameters are organized as
86	/// alpha_1,alpha_2, c_{1,1},c_{1,2},c_{1,3},...,c_{s,1},c_{s,2},c_{s,3},...
87	fn set_params(&mut self, params: &Vector<f64, Dyn, Self::ParameterStorage>) {	9,667✔
88	debug_assert_eq!(	9,667✔
89	params.len(),	9,667✔
90	self.Y.ncols() * 3 + 2,	9,667✔
UNCOV 91	"number of parameters does not match the dataset and the decay model"	×
92	);
93	//tau1
94	self.alpha[0] = params[0];	9,667✔
95	//tau2
96	self.alpha[1] = params[1];	9,667✔
97
98	// and do precalculations
99	self.Phi	9,667✔
100	.column_mut(0)
101	.copy_from(&self.x.map(\|x: f64\| (-x / self.alpha[0]).exp()));	203,007✔
102	self.Phi
103	.column_mut(1)
104	.copy_from(&self.x.map(\|x: f64\| (-x / self.alpha[1]).exp()));	193,340✔
105	// and load the matrix of linear coefficients
106	// I have a test down below that shows this should actually work
107	self.C.copy_from_slice(&params.as_slice()[2..]);
108	}
109
110	fn params(&self) -> Vector<f64, Dyn, Self::ParameterStorage> {	13✔
111	let param_count = self.C.len() + 2;	13✔
112	let params = DVector::from_iterator(
113	param_count,	13✔
114	self.alpha	13✔
115	.iter()	13✔
116	.cloned()	13✔
117	.chain(self.C.as_slice().iter().cloned()),	13✔
118	);
119	params	13✔
120	}
121
122	fn residuals(&self) -> Option<DVector<f64>> {	9,666✔
123	let R = &self.Y - &self.Phi * &self.C;	9,666✔
124	let new_nrows = Dyn(R.nrows() * R.ncols());	9,666✔
125	Some(R.reshape_generic(new_nrows, U1))	9,666✔
126	}
127
128	fn jacobian(&self) -> Option<Matrix<f64, Dyn, Dyn, Self::JacobianStorage>> {	50✔
129	// let mut jacobian = DMatrix::<f64>::zeros(self.
130	let mut dPhi_dalpha1 = OMatrix::<_, Dyn, U3>::zeros_generic(Dyn(self.Phi.nrows()), U3);	50✔
131	dPhi_dalpha1.set_column(	50✔
132	0,
133	&(-1. / (self.alpha[0] * self.alpha[0]) * self.Phi.column(0).component_mul(&self.x)),	50✔
134	);
135	let mut dPhi_dalpha2 = OMatrix::<_, Dyn, U3>::zeros_generic(Dyn(self.Phi.nrows()), U3);	50✔
136	dPhi_dalpha2.set_column(	50✔
137	1,
138	&(-1. / (self.alpha[1] * self.alpha[1]) * self.Phi.column(1).component_mul(&self.x)),	50✔
139	);
140	let mut jac_block = DMatrix::from_element(self.Phi.nrows(), 3, -1.);	50✔
141	jac_block.set_column(0, &-self.Phi.column(0));	50✔
142	jac_block.set_column(1, &-self.Phi.column(1));	50✔
143
144	let x_len = self.Y.nrows();	50✔
145	let total_data_len = x_len * self.Y.ncols();	50✔
146	let mut jacobian = DMatrix::zeros(self.Y.nrows() * self.Y.ncols(), self.C.len() + 2);	50✔
147	debug_assert_eq!(	50✔
148	self.Y.ncols() * 3,	50✔
149	self.C.len(),	50✔
UNCOV 150	"we need exactly 3 linear parameters per dataset"	×
151	);
152
153	let to_vector = \|mat: OMatrix<f64, Dyn, Dyn>\| {	150✔
154	let new_nrows = Dyn(mat.nrows() * mat.ncols());	100✔
155	mat.reshape_generic(new_nrows, U1)	100✔
156	};
157	jacobian.set_column(0, &to_vector(dPhi_dalpha1 * &self.C));
158	jacobian.set_column(1, &to_vector(dPhi_dalpha2 * &self.C));
159
160	for (col, row) in (2..jacobian.nrows())	100✔
161	.step_by(3)
162	.zip((0..total_data_len).step_by(x_len))
163	{
164	let mut view = jacobian.view_mut((row, col), (x_len, 3));	100✔
165	view.copy_from(&jac_block);	100✔
166	}
167	Some(jacobian)	50✔
168	}
169	}
170
171	#[test]
172	// make sure that the way I copy the parameters into a matrix works as
173	// I think it does
174	fn matrix_is_colum_major_layout() {
175	let data = &[1., 2., 3., 4., 5., 6.];
176	let mut mat = OMatrix::zeros_generic(Dyn(2), U3);
177	mat.copy_from_slice(data);
178	assert_eq!(mat.column(0).as_slice(), &[1., 2.]);
179	assert_eq!(mat.column(1).as_slice(), &[3., 4.]);
180	assert_eq!(mat.column(2).as_slice(), &[5., 6.]);
181	// also assert that the reverse is true
182	assert_eq!(mat.as_slice(), data);
183	}

geo-ant / varpro / 9181626839

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