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

99.32

/shared_test_code/src/models.rs

use nalgebra::{DVector, Dyn, OMatrix, OVector, U1};
use varpro::model::SeparableModel;
use varpro::{model::errors::ModelError, prelude::*};

/// contains a double exponential model with constant offset that does not
/// perform caching, as preliminary results indicate that using caching can
/// lead to performance penalties for large models
pub mod uncached;

/// A separable model for double exponential decay
/// with a constant offset
/// f_j = c1*exp(-x_j/tau1) + c2*exp(-x_j/tau2) + c3
/// this is a handcrafted model which uses caching for
/// maximum performance.
///
/// This is an example of how to implement a separable model
/// using the trait directly without using the builder.
/// This allows us to use caching of the intermediate results
/// to calculate the derivatives more efficiently.
#[derive(Clone, Debug, PartialEq)]
pub struct DoubleExpModelWithConstantOffsetSepModel {
    /// the x vector associated with this model.
    /// We make this configurable to enable models to
    /// work with different x coordinates, but this is
    /// not an inherent requirement. We don't have to have a
    /// field like this. The only requirement given by the trait
    /// is that the model can be queried about the length
    /// of its output
    x_vector: DVector<f64>,
    /// current parameters [tau1,tau2] of the exponential
    /// functions
    params: OVector<f64, Dyn>,
    /// cached evaluation of the model
    /// the matrix columns contain the complete evaluation
    /// of the model. That is the first column contains the
    /// exponential exp(-x/tau), the second column contains
    /// exp(-x/tau) both evaluated on the x vector. The third
    /// column contains a column of straight ones for the constant
    /// offset.
    ///
    /// This value is calculated in the set_params method, which is
    /// the only method with mutable access to the model state.
    eval: OMatrix<f64, Dyn, Dyn>,
}

impl DoubleExpModelWithConstantOffsetSepModel {
    /// create a new model with the given x vector and initial guesses
    /// for the exponential decay constants tau_1 and tau_2
    pub fn new(x_vector: DVector<f64>, (tau1_guess, tau2_guess): (f64, f64)) -> Self {
        let x_len = x_vector.len();
        let mut ret = Self {
            x_vector,
            params: OVector::zeros_generic(Dyn(2), U1), //<-- will be overwritten by set_params
            eval: OMatrix::<f64, Dyn, Dyn>::zeros_generic(Dyn(x_len), Dyn(3)),
        };
        ret.set_params(OVector::from_column_slice_generic(
            Dyn(2),
            U1,
            &[tau1_guess, tau2_guess],
        ))
        .unwrap();
        ret
    }
}

impl SeparableNonlinearModel for DoubleExpModelWithConstantOffsetSepModel {
    /// we give our own mddel error here, but we
    /// could also have indicated that our calculations cannot
    /// fail by using [`std::convert::Infallible`].
    type Error = ModelError;
    /// the actual scalar type that our model uses for calculations
    type ScalarType = f64;

    #[inline]
    fn parameter_count(&self) -> usize {
        // regardless of the fact that we know at compile time
        // that the length is 2, we still have to return an instance
        // of that type
        2
    }

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

    // we use this method not only to set the parameters inside the
    // model but we also cache some calculations. The advantage is that
    // we don't have to recalculate the exponential terms for either
    // the evaluations or the derivatives for the same parameters.
    fn set_params(&mut self, parameters: OVector<f64, Dyn>) -> Result<(), Self::Error> {
        // even if it is not the only thing we do, we still
        // have to update the internal parameters of the model
        self.params = parameters;

        // parameters expected in this order
        // use unsafe to avoid bounds checks
        let tau1 = unsafe { self.params.get_unchecked(0) };
        let tau2 = unsafe { self.params.get_unchecked(1) };

        // the exponential exp(-x/tau1)
        let f1 = self.x_vector.map(|x| f64::exp(-x / tau1));
        // the exponential exp(-x/tau2)
        let f2 = self.x_vector.map(|x| f64::exp(-x / tau2));

        self.eval.set_column(0, &f1);
        self.eval.set_column(1, &f2);
        self.eval
            .set_column(2, &DVector::from_element(self.x_vector.len(), 1.));
        Ok(())
    }

    fn params(&self) -> OVector<f64, Dyn> {
        self.params.clone()
    }

    // since we cached the model evaluation, we can just return
    // it here
    fn eval(&self) -> Result<OMatrix<f64, Dyn, Dyn>, Self::Error> {
        Ok(self.eval.clone())
    }

    // here we take advantage of the cached calculations
    // so that we do not have to recalculate the exponential
    // to calculate the derivative.
    fn eval_partial_deriv(
        &self,
        derivative_index: usize,
    ) -> Result<nalgebra::OMatrix<f64, Dyn, Dyn>, Self::Error> {
        let location = &self.x_vector;
        let parameters = &self.params;
        // 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];
        // the only nonzero derivative is the derivative of the exp(-x/tau) for
        // the corresponding tau at derivative_index
        // we can use the precalculated results so we don't have to use the
        // exponential function again
        let df = location
            .map(|x| x / (tau * tau))
            .component_mul(&self.eval.column(derivative_index));

        // two of the columns are always zero when we differentiate
        // with respect to tau_1 or tau_2. Remember the constant term
        // also occupies one column and will always be zero when differentiated
        // with respect to the nonlinear params of the model
        let mut derivatives = OMatrix::zeros_generic(Dyn(location.len()), Dyn(3));

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

    fn output_len(&self) -> usize {
        // this is how we give a length that is only known at runtime.
        // We wrap it in a `Dyn` instance.
        self.x_vector.len()
    }
}

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 nrows = self.x.len();
        let mut jacobian = Matrix::<f64, Dyn, U5, Self::JacobianStorage>::zeros(nrows);

        jacobian.set_column(
            0,
            &(c1 / (tau1 * tau1) * &(self.precalc_exp_tau1.component_mul(&self.x))),
        );
        jacobian.set_column(
            1,
            &(c2 / (tau2 * tau2) * &(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)
    }
}

/// example function from the matlab code that was published as part of the
/// O'Leary paper here: https://www.cs.umd.edu/~oleary/software/varpro/
/// The model function is
///
///   f(t) = c1 exp(-alpha2 t)*cos(alpha3 t)
///         + c2 exp(-alpha1 t)*cos(alpha2 t)
/// # performance
/// this model does reuse calculations but I think it could be optimized
/// further
#[derive(Clone, Debug, PartialEq)]
pub struct OLearyExampleModel {
    /// the t vector
    t: DVector<f64>,
    /// the current parameters (alpha1,alpha2,alpha3,alpha4)
    alpha: OVector<f64, Dyn>,
    /// the current evaluation of the model
    phi: OMatrix<f64, Dyn, Dyn>,
}

impl OLearyExampleModel {
    /// create a new model with the given t vector and initial guesses
    pub fn new(t: DVector<f64>, initial_guesses: OVector<f64, Dyn>) -> Self {
        let t_len = t.len();
        let mut ret = Self {
            t,
            alpha: initial_guesses.clone(),
            phi: OMatrix::<f64, Dyn, Dyn>::zeros_generic(Dyn(t_len), Dyn(2)),
        };
        ret.set_params(initial_guesses).unwrap();
        ret
    }
}

impl SeparableNonlinearModel for OLearyExampleModel {
    type ScalarType = f64;
    type Error = ModelError;

    fn parameter_count(&self) -> usize {
        3
    }

    fn base_function_count(&self) -> usize {
        2
    }

    fn output_len(&self) -> usize {
        self.t.len()
    }

    fn set_params(
        &mut self,
        parameters: OVector<Self::ScalarType, Dyn>,
    ) -> Result<(), Self::Error> {
        self.alpha = parameters;
        let alpha1 = self.alpha[0];
        let alpha2 = self.alpha[1];
        let alpha3 = self.alpha[2];

        let f1 = self.t.map(|t| f64::exp(-alpha2 * t) * f64::cos(alpha3 * t));
        let f2 = self.t.map(|t| f64::exp(-alpha1 * t) * f64::cos(alpha2 * t));

        self.phi = OMatrix::<f64, Dyn, Dyn>::from_columns(&[f1, f2]);
        Ok(())
    }

    fn params(&self) -> OVector<Self::ScalarType, Dyn> {
        self.alpha.clone()
    }

    fn eval(&self) -> Result<OMatrix<Self::ScalarType, Dyn, Dyn>, Self::Error> {
        Ok(self.phi.clone())
    }

    fn eval_partial_deriv(
        &self,
        derivative_index: usize,
    ) -> Result<OMatrix<Self::ScalarType, Dyn, Dyn>, Self::Error> {
        let mut derivs = OMatrix::<f64, Dyn, Dyn>::zeros_generic(Dyn(self.t.len()), Dyn(2));
        let alpha1 = self.alpha[0];
        let alpha2 = self.alpha[1];
        let alpha3 = self.alpha[2];

        // from the original matlab impl:
        // % The nonzero partial derivatives of Phi with respect to alpha are
        // %              d Phi_1 / d alpha_2 ,
        // %              d Phi_1 / d alpha_3 ,
        // %              d Phi_2 / d alpha_1 ,
        // %              d Phi_2 / d alpha_2 ,
        // % and this determines Ind.
        // % The ordering of the columns of Ind is arbitrary but must match dPhi.

        // % Evaluate the four nonzero partial derivatives of Phi at each of
        // % the data points and store them in dPhi.

        // dPhi 1 / dalpha 2 = -t .* Phi(:,1);
        // dPhi 1 / dalpha 3 = -t .* exp(-alpha(2)*t).*sin(alpha(3)*t);
        // dPhi 2 / dalpha 1 = -t .* Phi(:,2);
        // dPhi 2 / dalpha 2 = -t .* exp(-alpha(1)*t).*sin(alpha(2)*t);

        match derivative_index {
            0 => {
                // d/d alpha1
                derivs.set_column(1, &self.phi.column(1).component_mul(&(-1. * &self.t)));
            }
            1 => {
                // d/d alpha2
                derivs.set_column(0, &self.phi.column(0).component_mul(&(-1. * &self.t)));
                derivs.set_column(
                    1,
                    &self
                        .t
                        .map(|t| -t * (-alpha1 * t).exp() * (alpha2 * t).sin()),
                );
            }
            2 => {
                // d/d alpha3
                derivs.set_column(
                    0,
                    &self
                        .t
                        .map(|t| -t * (-alpha2 * t).exp() * (alpha3 * t).sin()),
                );
            }
            _ => {
                unreachable!("derivative index must be 0,1,2");
            }
        }

        Ok(derivs)
    }
}

/// the oleary example model as above but this time it is generated using
/// the separable model builder
pub fn o_leary_example_model(t: DVector<f64>, initial_guesses: Vec<f64>) -> SeparableModel<f64> {
    fn phi1(t: &DVector<f64>, alpha2: f64, alpha3: f64) -> DVector<f64> {
        t.map(|t| f64::exp(-alpha2 * t) * f64::cos(alpha3 * t))
    }
    fn phi2(t: &DVector<f64>, alpha1: f64, alpha2: f64) -> DVector<f64> {
        t.map(|t| f64::exp(-alpha1 * t) * f64::cos(alpha2 * t))
    }

    SeparableModelBuilder::new(["alpha1", "alpha2", "alpha3"])
        .initial_parameters(initial_guesses)
        .independent_variable(t)
        // phi 1
        .function(["alpha2", "alpha3"], phi1)
        .partial_deriv("alpha2", |t: &DVector<f64>, alpha2: f64, alpha3: f64| {
            phi1(t, alpha2, alpha3).component_mul(&(-1. * t))
        })
        .partial_deriv("alpha3", |t: &DVector<f64>, alpha2: f64, alpha3: f64| {
            t.map(|t| -t * (-alpha2 * t).exp() * (alpha3 * t).sin())
        })
        .function(["alpha1", "alpha2"], phi2)
        .partial_deriv("alpha1", |t: &DVector<f64>, alpha1: f64, alpha2: f64| {
            phi2(t, alpha1, alpha2).component_mul(&(-1. * t))
        })
        .partial_deriv("alpha2", |t: &DVector<f64>, alpha1: f64, alpha2: f64| {
            t.map(|t| -t * (-alpha1 * t).exp() * (alpha2 * t).sin())
        })
        .build()
        .unwrap()
}

1	use nalgebra::{DVector, Dyn, OMatrix, OVector, U1};
2	use varpro::model::SeparableModel;
3	use varpro::{model::errors::ModelError, prelude::*};
4
5	/// contains a double exponential model with constant offset that does not
6	/// perform caching, as preliminary results indicate that using caching can
7	/// lead to performance penalties for large models
8	pub mod uncached;
9
10	/// A separable model for double exponential decay
11	/// with a constant offset
12	/// f_j = c1exp(-x_j/tau1) + c2exp(-x_j/tau2) + c3
13	/// this is a handcrafted model which uses caching for
14	/// maximum performance.
15	///
16	/// This is an example of how to implement a separable model
17	/// using the trait directly without using the builder.
18	/// This allows us to use caching of the intermediate results
19	/// to calculate the derivatives more efficiently.
20	#[derive(Clone, Debug, PartialEq)]
21	pub struct DoubleExpModelWithConstantOffsetSepModel {
22	/// the x vector associated with this model.
23	/// We make this configurable to enable models to
24	/// work with different x coordinates, but this is
25	/// not an inherent requirement. We don't have to have a
26	/// field like this. The only requirement given by the trait
27	/// is that the model can be queried about the length
28	/// of its output
29	x_vector: DVector<f64>,
30	/// current parameters [tau1,tau2] of the exponential
31	/// functions
32	params: OVector<f64, Dyn>,
33	/// cached evaluation of the model
34	/// the matrix columns contain the complete evaluation
35	/// of the model. That is the first column contains the
36	/// exponential exp(-x/tau), the second column contains
37	/// exp(-x/tau) both evaluated on the x vector. The third
38	/// column contains a column of straight ones for the constant
39	/// offset.
40	///
41	/// This value is calculated in the set_params method, which is
42	/// the only method with mutable access to the model state.
43	eval: OMatrix<f64, Dyn, Dyn>,
44	}
45
46	impl DoubleExpModelWithConstantOffsetSepModel {
47	/// create a new model with the given x vector and initial guesses
48	/// for the exponential decay constants tau_1 and tau_2
49	pub fn new(x_vector: DVector<f64>, (tau1_guess, tau2_guess): (f64, f64)) -> Self {	4✔
50	let x_len = x_vector.len();	4✔
51	let mut ret = Self {
52	x_vector,
53	params: OVector::zeros_generic(Dyn(2), U1), //<-- will be overwritten by set_params	4✔
54	eval: OMatrix::<f64, Dyn, Dyn>::zeros_generic(Dyn(x_len), Dyn(3)),	4✔
55	};
56	ret.set_params(OVector::from_column_slice_generic(	4✔
57	Dyn(2),	4✔
58	U1,	4✔
59	&[tau1_guess, tau2_guess],	4✔
60	))
61	.unwrap();
62	ret	4✔
63	}
64	}
65
66	impl SeparableNonlinearModel for DoubleExpModelWithConstantOffsetSepModel {
67	/// we give our own mddel error here, but we
68	/// could also have indicated that our calculations cannot
69	/// fail by using [`std::convert::Infallible`].
70	type Error = ModelError;
71	/// the actual scalar type that our model uses for calculations
72	type ScalarType = f64;
73
74	#[inline]
75	fn parameter_count(&self) -> usize {	59✔
76	// regardless of the fact that we know at compile time
77	// that the length is 2, we still have to return an instance
78	// of that type
79	2	59✔
80	}
81
82	#[inline]
83	fn base_function_count(&self) -> usize {	7✔
84	// same as above
85	3	7✔
86	}
87
88	// we use this method not only to set the parameters inside the
89	// model but we also cache some calculations. The advantage is that
90	// we don't have to recalculate the exponential terms for either
91	// the evaluations or the derivatives for the same parameters.
92	fn set_params(&mut self, parameters: OVector<f64, Dyn>) -> Result<(), Self::Error> {	80✔
93	// even if it is not the only thing we do, we still
94	// have to update the internal parameters of the model
95	self.params = parameters;	80✔
96
97	// parameters expected in this order
98	// use unsafe to avoid bounds checks
99	let tau1 = unsafe { self.params.get_unchecked(0) };	80✔
100	let tau2 = unsafe { self.params.get_unchecked(1) };	80✔
101
102	// the exponential exp(-x/tau1)
103	let f1 = self.x_vector.map(\|x\| f64::exp(-x / tau1));	54,348✔
104	// the exponential exp(-x/tau2)
105	let f2 = self.x_vector.map(\|x\| f64::exp(-x / tau2));	54,348✔
106
107	self.eval.set_column(0, &f1);	80✔
108	self.eval.set_column(1, &f2);	80✔
109	self.eval	80✔
110	.set_column(2, &DVector::from_element(self.x_vector.len(), 1.));	80✔
111	Ok(())	80✔
112	}
113
114	fn params(&self) -> OVector<f64, Dyn> {	21✔
115	self.params.clone()	21✔
116	}
117
118	// since we cached the model evaluation, we can just return
119	// it here
120	fn eval(&self) -> Result<OMatrix<f64, Dyn, Dyn>, Self::Error> {	83✔
121	Ok(self.eval.clone())	83✔
122	}
123
124	// here we take advantage of the cached calculations
125	// so that we do not have to recalculate the exponential
126	// to calculate the derivative.
127	fn eval_partial_deriv(	106✔
128	&self,
129	derivative_index: usize,
130	) -> Result<nalgebra::OMatrix<f64, Dyn, Dyn>, Self::Error> {
131	let location = &self.x_vector;	106✔
132	let parameters = &self.params;	106✔
133	// derivative index can be either 0,1 (corresponding to the linear parameters
134	// tau1, tau2). Since only one of the basis functions depends on
135	// tau_i, we can simplify calculations here
136
137	let tau = parameters[derivative_index];	106✔
138	// the only nonzero derivative is the derivative of the exp(-x/tau) for
139	// the corresponding tau at derivative_index
140	// we can use the precalculated results so we don't have to use the
141	// exponential function again
142	let df = location	106✔
143	.map(\|x\| x / (tau * tau))	73,660✔
144	.component_mul(&self.eval.column(derivative_index));	106✔
145
146	// two of the columns are always zero when we differentiate
147	// with respect to tau_1 or tau_2. Remember the constant term
148	// also occupies one column and will always be zero when differentiated
149	// with respect to the nonlinear params of the model
150	let mut derivatives = OMatrix::zeros_generic(Dyn(location.len()), Dyn(3));	106✔
151
152	derivatives.set_column(derivative_index, &df);	106✔
153	Ok(derivatives)	106✔
154	}
155
156	fn output_len(&self) -> usize {	63✔
157	// this is how we give a length that is only known at runtime.
158	// We wrap it in a `Dyn` instance.
159	self.x_vector.len()	63✔
160	}
161	}
162
163	use levenberg_marquardt::LeastSquaresProblem;
164	use nalgebra::storage::Owned;
165	use nalgebra::{Matrix, Vector, Vector5, U5};
166
167	/// Implementation of double exponential decay with constant offset
168	/// f(x) = c1exp(-x/tau1)+c2exp(-x/tau2)+c3
169	/// using the levenberg_marquardt crate
170	pub struct DoubleExponentialDecayFittingWithOffsetLevmar {
171	/// current problem paramters with layout (tau1,tau2,c1,c2,c3)
172	params: Vector5<f64>,
173	/// current independent variable
174	x: DVector<f64>,
175	/// the data (must have same length as x)
176	y: DVector<f64>,
177	/// precached calculations for the exponential terms so that the jacobian and residuals can
178	/// 1) exp(-x/tau1)
179	precalc_exp_tau1: DVector<f64>,
180	/// 2) exp(-x/tau2)
181	precalc_exp_tau2: DVector<f64>,
182	}
183
184	impl DoubleExponentialDecayFittingWithOffsetLevmar {
185	/// create new fitting problem with data and initial guesses
186	pub fn new(initial_guesses: &[f64], x: &DVector<f64>, y: &DVector<f64>) -> Self {	2✔
187	assert_eq!(	2✔
188	initial_guesses.len(),	2✔
189	5,
UNCOV 190	"Wrong parameter count. The 5 Parameters are: tau1, tau2, c1, c2, c3"	×
191	);
192	assert_eq!(y.len(), x.len(), "x and y must have same length");	2✔
193	let parameters = Vector5::from_iterator(initial_guesses.iter().copied());	2✔
194	let mut problem = Self {
195	params: parameters,
196	x: x.clone(),	2✔
197	y: y.clone(),	2✔
198	precalc_exp_tau1: DVector::from_element(x.len(), 0.), //will both be overwritten with correct vals	2✔
199	precalc_exp_tau2: DVector::from_element(x.len(), 0.), //by the following call to set params	2✔
200	};
201	problem.set_params(&parameters);	2✔
202	problem	2✔
203	}
204	}
205
206	impl LeastSquaresProblem<f64, Dyn, U5> for DoubleExponentialDecayFittingWithOffsetLevmar {
207	type ParameterStorage = Owned<f64, U5>;
208	type ResidualStorage = Owned<f64, Dyn>;
209	type JacobianStorage = Owned<f64, Dyn, U5>;
210
211	fn set_params(&mut self, params: &Vector<f64, U5, Self::ParameterStorage>) {	3,076✔
212	self.params = *params;	3,076✔
213	let tau1 = params[0];	3,076✔
214	let tau2 = params[1];	3,076✔
215
216	// and do precalculations
217	self.precalc_exp_tau1 = self.x.map(\|x: f64\| (-x / tau1).exp());	141,968✔
218	self.precalc_exp_tau2 = self.x.map(\|x: f64\| (-x / tau2).exp());	141,968✔
219	}
220
221	fn params(&self) -> Vector<f64, U5, Self::ParameterStorage> {	9✔
222	self.params	9✔
223	}
224
225	fn residuals(&self) -> Option<DVector<f64>> {	3,074✔
226	// get parameters from internal param storage
227	//let tau1 = self.params[0];
228	//let tau2 = self.params[1];
229	let c1 = self.params[2];	3,074✔
230	let c2 = self.params[3];	3,074✔
231	let c3 = self.params[4];	3,074✔
232
233	// function values at the parameters
234	let f = c1 * &self.precalc_exp_tau1	3,074✔
235	+ c2 * &self.precalc_exp_tau2	3,074✔
236	+ DVector::from_element(self.x.len(), c3);	3,074✔
237
238	// residuals: f-y
239	Some(f - &self.y)	3,074✔
240	}
241
242	fn jacobian(&self) -> Option<Matrix<f64, Dyn, U5, Self::JacobianStorage>> {	56✔
243	// get parameters from internal param storage
244	let tau1 = self.params[0];	56✔
245	let tau2 = self.params[1];	56✔
246	let c1 = self.params[2];	56✔
247	let c2 = self.params[3];	56✔
248	// populate jacobian
249	let nrows = self.x.len();	56✔
250	let mut jacobian = Matrix::<f64, Dyn, U5, Self::JacobianStorage>::zeros(nrows);	56✔
251
252	jacobian.set_column(	56✔
253	0,
254	&(c1 / (tau1 * tau1) * &(self.precalc_exp_tau1.component_mul(&self.x))),	56✔
255	);
256	jacobian.set_column(	56✔
257	1,
258	&(c2 / (tau2 * tau2) * &(self.precalc_exp_tau2.component_mul(&self.x))),	56✔
259	);
260	jacobian.set_column(2, &self.precalc_exp_tau1);	56✔
261	jacobian.set_column(3, &self.precalc_exp_tau2);	56✔
262	jacobian.set_column(4, &DVector::from_element(self.x.len(), 1.));	56✔
263
264	Some(jacobian)	56✔
265	}
266	}
267
268	/// example function from the matlab code that was published as part of the
269	/// O'Leary paper here: https://www.cs.umd.edu/~oleary/software/varpro/
270	/// The model function is
271	///
272	/// f(t) = c1 exp(-alpha2 t)*cos(alpha3 t)
273	/// + c2 exp(-alpha1 t)*cos(alpha2 t)
274	/// # performance
275	/// this model does reuse calculations but I think it could be optimized
276	/// further
277	#[derive(Clone, Debug, PartialEq)]
278	pub struct OLearyExampleModel {
279	/// the t vector
280	t: DVector<f64>,
281	/// the current parameters (alpha1,alpha2,alpha3,alpha4)
282	alpha: OVector<f64, Dyn>,
283	/// the current evaluation of the model
284	phi: OMatrix<f64, Dyn, Dyn>,
285	}
286
287	impl OLearyExampleModel {
288	/// create a new model with the given t vector and initial guesses
289	pub fn new(t: DVector<f64>, initial_guesses: OVector<f64, Dyn>) -> Self {	1✔
290	let t_len = t.len();	1✔
291	let mut ret = Self {
292	t,
293	alpha: initial_guesses.clone(),	1✔
294	phi: OMatrix::<f64, Dyn, Dyn>::zeros_generic(Dyn(t_len), Dyn(2)),	1✔
295	};
296	ret.set_params(initial_guesses).unwrap();	1✔
297	ret	1✔
298	}
299	}
300
301	impl SeparableNonlinearModel for OLearyExampleModel {
302	type ScalarType = f64;
303	type Error = ModelError;
304
305	fn parameter_count(&self) -> usize {	12✔
306	3	12✔
307	}
308
309	fn base_function_count(&self) -> usize {	2✔
310	2	2✔
311	}
312
313	fn output_len(&self) -> usize {	13✔
314	self.t.len()	13✔
315	}
316
317	fn set_params(	16✔
318	&mut self,
319	parameters: OVector<Self::ScalarType, Dyn>,
320	) -> Result<(), Self::Error> {
321	self.alpha = parameters;	16✔
322	let alpha1 = self.alpha[0];	16✔
323	let alpha2 = self.alpha[1];	16✔
324	let alpha3 = self.alpha[2];	16✔
325
326	let f1 = self.t.map(\|t\| f64::exp(-alpha2 * t) * f64::cos(alpha3 * t));	192✔
327	let f2 = self.t.map(\|t\| f64::exp(-alpha1 * t) * f64::cos(alpha2 * t));	192✔
328
329	self.phi = OMatrix::<f64, Dyn, Dyn>::from_columns(&[f1, f2]);	16✔
330	Ok(())	16✔
331	}
332
333	fn params(&self) -> OVector<Self::ScalarType, Dyn> {	3✔
334	self.alpha.clone()	3✔
335	}
336
337	fn eval(&self) -> Result<OMatrix<Self::ScalarType, Dyn, Dyn>, Self::Error> {	18✔
338	Ok(self.phi.clone())	18✔
339	}
340
341	fn eval_partial_deriv(	30✔
342	&self,
343	derivative_index: usize,
344	) -> Result<OMatrix<Self::ScalarType, Dyn, Dyn>, Self::Error> {
345	let mut derivs = OMatrix::<f64, Dyn, Dyn>::zeros_generic(Dyn(self.t.len()), Dyn(2));	30✔
346	let alpha1 = self.alpha[0];	30✔
347	let alpha2 = self.alpha[1];	30✔
348	let alpha3 = self.alpha[2];	30✔
349
350	// from the original matlab impl:
351	// % The nonzero partial derivatives of Phi with respect to alpha are
352	// % d Phi_1 / d alpha_2 ,
353	// % d Phi_1 / d alpha_3 ,
354	// % d Phi_2 / d alpha_1 ,
355	// % d Phi_2 / d alpha_2 ,
356	// % and this determines Ind.
357	// % The ordering of the columns of Ind is arbitrary but must match dPhi.
358
359	// % Evaluate the four nonzero partial derivatives of Phi at each of
360	// % the data points and store them in dPhi.
361
362	// dPhi 1 / dalpha 2 = -t .* Phi(:,1);
363	// dPhi 1 / dalpha 3 = -t .* exp(-alpha(2)t).sin(alpha(3)*t);
364	// dPhi 2 / dalpha 1 = -t .* Phi(:,2);
365	// dPhi 2 / dalpha 2 = -t .* exp(-alpha(1)t).sin(alpha(2)*t);
366
367	match derivative_index {	30✔
368	0 => {	10✔
369	// d/d alpha1
370	derivs.set_column(1, &self.phi.column(1).component_mul(&(-1. * &self.t)));	10✔
371	}
372	1 => {	10✔
373	// d/d alpha2
374	derivs.set_column(0, &self.phi.column(0).component_mul(&(-1. * &self.t)));	10✔
375	derivs.set_column(	10✔
376	1,
377	&self	10✔
378	.t	10✔
379	.map(\|t\| -t * (-alpha1 * t).exp() * (alpha2 * t).sin()),	120✔
380	);
381	}
382	2 => {	10✔
383	// d/d alpha3
384	derivs.set_column(	10✔
385	0,
386	&self	10✔
387	.t	10✔
388	.map(\|t\| -t * (-alpha2 * t).exp() * (alpha3 * t).sin()),	120✔
389	);
390	}
391	_ => {
392	unreachable!("derivative index must be 0,1,2");
393	}
394	}
395
396	Ok(derivs)	30✔
397	}
398	}
399
400	/// the oleary example model as above but this time it is generated using
401	/// the separable model builder
402	pub fn o_leary_example_model(t: DVector<f64>, initial_guesses: Vec<f64>) -> SeparableModel<f64> {	1✔
403	fn phi1(t: &DVector<f64>, alpha2: f64, alpha3: f64) -> DVector<f64> {	29✔
404	t.map(\|t\| f64::exp(-alpha2 * t) * f64::cos(alpha3 * t))	337✔
405	}
406	fn phi2(t: &DVector<f64>, alpha1: f64, alpha2: f64) -> DVector<f64> {	29✔
407	t.map(\|t\| f64::exp(-alpha1 * t) * f64::cos(alpha2 * t))	337✔
408	}
409
410	SeparableModelBuilder::new(["alpha1", "alpha2", "alpha3"])	1✔
411	.initial_parameters(initial_guesses)	1✔
412	.independent_variable(t)	1✔
413	// phi 1
414	.function(["alpha2", "alpha3"], phi1)	1✔
415	.partial_deriv("alpha2", \|t: &DVector<f64>, alpha2: f64, alpha3: f64\| {	11✔
416	phi1(t, alpha2, alpha3).component_mul(&(-1. * t))	10✔
417	})
418	.partial_deriv("alpha3", \|t: &DVector<f64>, alpha2: f64, alpha3: f64\| {	11✔
419	t.map(\|t\| -t * (-alpha2 * t).exp() * (alpha3 * t).sin())	120✔
420	})
421	.function(["alpha1", "alpha2"], phi2)	1✔
422	.partial_deriv("alpha1", \|t: &DVector<f64>, alpha1: f64, alpha2: f64\| {	11✔
423	phi2(t, alpha1, alpha2).component_mul(&(-1. * t))	10✔
424	})
425	.partial_deriv("alpha2", \|t: &DVector<f64>, alpha1: f64, alpha2: f64\| {	11✔
426	t.map(\|t\| -t * (-alpha1 * t).exp() * (alpha2 * t).sin())	120✔
427	})
428	.build()
429	.unwrap()
430	}

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