4141916717

Build # 4141916717

Build Type

github

Committed by Unknown Committer

Commit Message

Unknown Commit Message

Pull Request Pull Request #19: Feature/revamp model function experiment

Run Details

98 of 98 new or added lines in 5 files covered. (100.0%)

330 of 488 relevant lines covered (67.62%)

3.13 hits per line

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

0.0

/shared_test_code/src/models.rs

use nalgebra::{DMatrix, DVector, Dyn};
use varpro::{model::errors::ModelError, prelude::*};

#[derive(Default, Clone)]
/// A separable model for double exponential decay
/// with a constant offset
/// f_j = c1*exp(-x_j/tau1) + c2*exp(-x_j/tau2) + c3
pub struct DoubleExpModelWithConstantOffsetSepModel {}

impl SeparableNonlinearModel<f64> for DoubleExpModelWithConstantOffsetSepModel {
    type Error = ModelError;

    #[inline]
    fn parameter_count(&self) -> usize {
        2
    }

    #[inline]
    fn base_function_count(&self) -> usize {
        3
    }

    fn eval(
        &self,
        location: &nalgebra::DVector<f64>,
        parameters: &[f64],
    ) -> Result<nalgebra::DMatrix<f64>, Self::Error> {
        // parameters expected in this order
        // use unsafe to avoid bounds checks
        let tau1 = unsafe { parameters.get_unchecked(0) };
        let tau2 = unsafe { parameters.get_unchecked(1) };

        let f1 = location.map(|x| f64::exp(-x / tau1));
        let f2 = location.map(|x| f64::exp(-x / tau2));

        let mut basefuncs = unsafe {
            nalgebra::DMatrix::uninit(Dyn(location.len()), nalgebra::Dyn(3)).assume_init()
        };

        basefuncs.set_column(0, &f1);
        basefuncs.set_column(1, &f2);
        basefuncs.set_column(2, &DVector::from_element(location.len(), 1.));
        Ok(basefuncs)
    }

    fn eval_partial_deriv(
        &self,
        location: &nalgebra::DVector<f64>,
        parameters: &[f64],
        derivative_index: usize,
    ) -> Result<nalgebra::DMatrix<f64>, Self::Error> {
        // derivative index can be either 0,1 (corresponding to the linear parameters
        // tau1, tau2). Since only one of the basis functions depends on
        // tau_i, we can simplify calculations here

        let tau = parameters[derivative_index];
        let df = location.map(|x| x / (tau * tau) * f64::exp(-x / tau));

        let mut basefuncs = DMatrix::zeros(location.len(), 3);

        basefuncs.set_column(derivative_index, &df);
        Ok(basefuncs)
    }
}

use levenberg_marquardt::LeastSquaresProblem;
use nalgebra::storage::Owned;
use nalgebra::{Matrix, Vector, Vector5, U5};

/// Implementation of double exponential decay with constant offset
/// f(x) = c1*exp(-x/tau1)+c2*exp(-x/tau2)+c3
/// using the levenberg_marquardt crate
pub struct DoubleExponentialDecayFittingWithOffsetLevmar {
    /// current problem paramters with layout (tau1,tau2,c1,c2,c3)
    params: Vector5<f64>,
    /// current independent variable
    x: DVector<f64>,
    /// the data (must have same length as x)
    y: DVector<f64>,
    /// precached calculations for the exponential terms so that the jacobian and residuals can
    /// 1) exp(-x/tau1)
    precalc_exp_tau1: DVector<f64>,
    /// 2) exp(-x/tau2)
    precalc_exp_tau2: DVector<f64>,
}

impl DoubleExponentialDecayFittingWithOffsetLevmar {
    /// create new fitting problem with data and initial guesses
    pub fn new(initial_guesses: &[f64], x: &DVector<f64>, y: &DVector<f64>) -> Self {
        assert_eq!(
            initial_guesses.len(),
            5,
            "Wrong parameter count. The 5 Parameters are: tau1, tau2, c1, c2, c3"
        );
        assert_eq!(y.len(), x.len(), "x and y must have same length");
        let parameters = Vector5::from_iterator(initial_guesses.iter().copied());
        let mut problem = Self {
            params: parameters,
            x: x.clone(),
            y: y.clone(),
            precalc_exp_tau1: DVector::from_element(x.len(), 0.), //will both be overwritten with correct vals
            precalc_exp_tau2: DVector::from_element(x.len(), 0.), //by the following call to set params
        };
        problem.set_params(&parameters);
        problem
    }
}

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

    fn set_params(&mut self, params: &Vector<f64, U5, Self::ParameterStorage>) {
        self.params = *params;
        let tau1 = params[0];
        let tau2 = params[1];

        // and do precalculations
        self.precalc_exp_tau1 = self.x.map(|x: f64| (-x / tau1).exp());
        self.precalc_exp_tau2 = self.x.map(|x: f64| (-x / tau2).exp());
    }

    fn params(&self) -> Vector<f64, U5, Self::ParameterStorage> {
        self.params
    }

    fn residuals(&self) -> Option<DVector<f64>> {
        // get parameters from internal param storage
        //let tau1 = self.params[0];
        //let tau2 = self.params[1];
        let c1 = self.params[2];
        let c2 = self.params[3];
        let c3 = self.params[4];

        // function values at the parameters
        let f = c1 * &self.precalc_exp_tau1
            + c2 * &self.precalc_exp_tau2
            + &DVector::from_element(self.x.len(), c3);

        // residuals: f-y
        Some(&f - &self.y)
    }

    fn jacobian(&self) -> Option<Matrix<f64, Dyn, U5, Self::JacobianStorage>> {
        // get parameters from internal param storage
        let tau1 = self.params[0];
        let tau2 = self.params[1];
        let c1 = self.params[2];
        let c2 = self.params[3];
        // populate jacobian
        //let ncols = 5;
        let nrows = self.x.len();
        let mut jacobian = Matrix::<f64, Dyn, U5, Self::JacobianStorage>::zeros(nrows);

        jacobian.set_column(
            0,
            &(c1 / (tau1.powi(2)) * &(self.precalc_exp_tau1.component_mul(&self.x))),
        );
        jacobian.set_column(
            1,
            &(c2 / (tau2.powi(2)) * &(self.precalc_exp_tau2.component_mul(&self.x))),
        );
        jacobian.set_column(2, &self.precalc_exp_tau1);
        jacobian.set_column(3, &self.precalc_exp_tau2);
        jacobian.set_column(4, &DVector::from_element(self.x.len(), 1.));

        Some(jacobian)
    }
}

1	use nalgebra::{DMatrix, DVector, Dyn};
2	use varpro::{model::errors::ModelError, prelude::*};
3
4	#[derive(Default, Clone)]
5	/// A separable model for double exponential decay
6	/// with a constant offset
7	/// f_j = c1exp(-x_j/tau1) + c2exp(-x_j/tau2) + c3
8	pub struct DoubleExpModelWithConstantOffsetSepModel {}
9
10	impl SeparableNonlinearModel<f64> for DoubleExpModelWithConstantOffsetSepModel {
11	type Error = ModelError;
12
13	#[inline]
14	fn parameter_count(&self) -> usize {	×
15	2	×
16	}
17
18	#[inline]
19	fn base_function_count(&self) -> usize {	×
20	3	×
21	}
22
23	fn eval(	×
24	&self,
25	location: &nalgebra::DVector<f64>,
26	parameters: &[f64],
27	) -> Result<nalgebra::DMatrix<f64>, Self::Error> {
28	// parameters expected in this order
29	// use unsafe to avoid bounds checks
30	let tau1 = unsafe { parameters.get_unchecked(0) };	×
31	let tau2 = unsafe { parameters.get_unchecked(1) };	×
32
33	let f1 = location.map(\|x\| f64::exp(-x / tau1));	×
34	let f2 = location.map(\|x\| f64::exp(-x / tau2));	×
35
36	let mut basefuncs = unsafe {
37	nalgebra::DMatrix::uninit(Dyn(location.len()), nalgebra::Dyn(3)).assume_init()	×
38	};
39
40	basefuncs.set_column(0, &f1);	×
41	basefuncs.set_column(1, &f2);	×
42	basefuncs.set_column(2, &DVector::from_element(location.len(), 1.));	×
43	Ok(basefuncs)	×
44	}
45
46	fn eval_partial_deriv(	×
47	&self,
48	location: &nalgebra::DVector<f64>,
49	parameters: &[f64],
50	derivative_index: usize,
51	) -> Result<nalgebra::DMatrix<f64>, Self::Error> {
52	// derivative index can be either 0,1 (corresponding to the linear parameters
53	// tau1, tau2). Since only one of the basis functions depends on
54	// tau_i, we can simplify calculations here
55
56	let tau = parameters[derivative_index];	×
57	let df = location.map(\|x\| x / (tau * tau) * f64::exp(-x / tau));	×
58
59	let mut basefuncs = DMatrix::zeros(location.len(), 3);	×
60
61	basefuncs.set_column(derivative_index, &df);	×
62	Ok(basefuncs)	×
63	}
64	}
65
66	use levenberg_marquardt::LeastSquaresProblem;
67	use nalgebra::storage::Owned;
68	use nalgebra::{Matrix, Vector, Vector5, U5};
69
70	/// Implementation of double exponential decay with constant offset
71	/// f(x) = c1exp(-x/tau1)+c2exp(-x/tau2)+c3
72	/// using the levenberg_marquardt crate
73	pub struct DoubleExponentialDecayFittingWithOffsetLevmar {
74	/// current problem paramters with layout (tau1,tau2,c1,c2,c3)
75	params: Vector5<f64>,
76	/// current independent variable
77	x: DVector<f64>,
78	/// the data (must have same length as x)
79	y: DVector<f64>,
80	/// precached calculations for the exponential terms so that the jacobian and residuals can
81	/// 1) exp(-x/tau1)
82	precalc_exp_tau1: DVector<f64>,
83	/// 2) exp(-x/tau2)
84	precalc_exp_tau2: DVector<f64>,
85	}
86
87	impl DoubleExponentialDecayFittingWithOffsetLevmar {
88	/// create new fitting problem with data and initial guesses
89	pub fn new(initial_guesses: &[f64], x: &DVector<f64>, y: &DVector<f64>) -> Self {	×
90	assert_eq!(	×
91	initial_guesses.len(),	×
92	5,
93	"Wrong parameter count. The 5 Parameters are: tau1, tau2, c1, c2, c3"
94	);
95	assert_eq!(y.len(), x.len(), "x and y must have same length");	×
96	let parameters = Vector5::from_iterator(initial_guesses.iter().copied());	×
97	let mut problem = Self {
98	params: parameters,
99	x: x.clone(),	×
100	y: y.clone(),	×
101	precalc_exp_tau1: DVector::from_element(x.len(), 0.), //will both be overwritten with correct vals	×
102	precalc_exp_tau2: DVector::from_element(x.len(), 0.), //by the following call to set params	×
103	};
104	problem.set_params(&parameters);	×
105	problem
106	}
107	}
108
109	impl LeastSquaresProblem<f64, Dyn, U5> for DoubleExponentialDecayFittingWithOffsetLevmar {
110	type ParameterStorage = Owned<f64, U5>;
111	type ResidualStorage = Owned<f64, Dyn>;
112	type JacobianStorage = Owned<f64, Dyn, U5>;
113
114	fn set_params(&mut self, params: &Vector<f64, U5, Self::ParameterStorage>) {	×
115	self.params = *params;	×
116	let tau1 = params[0];	×
117	let tau2 = params[1];	×
118
119	// and do precalculations
120	self.precalc_exp_tau1 = self.x.map(\|x: f64\| (-x / tau1).exp());	×
121	self.precalc_exp_tau2 = self.x.map(\|x: f64\| (-x / tau2).exp());	×
122	}
123
124	fn params(&self) -> Vector<f64, U5, Self::ParameterStorage> {	×
125	self.params	×
126	}
127
128	fn residuals(&self) -> Option<DVector<f64>> {	×
129	// get parameters from internal param storage
130	//let tau1 = self.params[0];
131	//let tau2 = self.params[1];
132	let c1 = self.params[2];	×
133	let c2 = self.params[3];	×
134	let c3 = self.params[4];	×
135
136	// function values at the parameters
137	let f = c1 * &self.precalc_exp_tau1	×
138	+ c2 * &self.precalc_exp_tau2	×
139	+ &DVector::from_element(self.x.len(), c3);	×
140
141	// residuals: f-y
142	Some(&f - &self.y)	×
143	}
144
145	fn jacobian(&self) -> Option<Matrix<f64, Dyn, U5, Self::JacobianStorage>> {	×
146	// get parameters from internal param storage
147	let tau1 = self.params[0];	×
148	let tau2 = self.params[1];	×
149	let c1 = self.params[2];	×
150	let c2 = self.params[3];	×
151	// populate jacobian
152	//let ncols = 5;
153	let nrows = self.x.len();	×
154	let mut jacobian = Matrix::<f64, Dyn, U5, Self::JacobianStorage>::zeros(nrows);	×
155
156	jacobian.set_column(	×
157	0,
158	&(c1 / (tau1.powi(2)) * &(self.precalc_exp_tau1.component_mul(&self.x))),	×
159	);
160	jacobian.set_column(	×
161	1,
162	&(c2 / (tau2.powi(2)) * &(self.precalc_exp_tau2.component_mul(&self.x))),	×
163	);
164	jacobian.set_column(2, &self.precalc_exp_tau1);	×
165	jacobian.set_column(3, &self.precalc_exp_tau2);	×
166	jacobian.set_column(4, &DVector::from_element(self.x.len(), 1.));	×
167
168	Some(jacobian)	×
169	}
170	}

geo-ant / varpro / 4141916717

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