1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
use crate::numerical::optimization::LM_optimization::MinimizationReport;
use crate::numerical::optimization::varpro::model::SeparableNonlinearModel;
use crate::numerical::optimization::varpro::problem::{
MultiRhs, RhsType, SeparableProblem, SingleRhs,
};
use nalgebra::{DMatrix, DMatrixView, DVectorView, Dyn, OMatrix, OVector};
use std::collections::HashMap;
/// A helper type that contains the fitting problem after the
/// minimization, as well as a report and some convenience functions.
///
/// This structure is returned by the [`LevMarSolver::solve`](crate::solvers::levmar::LevMarSolver::solve)
/// and other solve methods of [`LevMarSolver`](crate::solvers::levmar::LevMarSolver).
#[derive(Debug)]
pub struct FitResult<Model, Rhs>
where
Model: SeparableNonlinearModel<ScalarType = f64>,
Rhs: RhsType,
{
/// The final state of the fitting problem after the
/// minimization finished (regardless of whether fitting was successful or not).
pub problem: SeparableProblem<Model, Rhs>,
/// the linear coefficient matrix (in case of single rhs: matrix has
/// only one column) at the solution.
pub(crate) linear_coefficients: Option<DMatrix<Model::ScalarType>>,
/// The minimization report of the underlying solver.
/// It contains information about the minimization process
/// and should be queried to see whether the minimization
/// was considered successful.
pub minimization_report: MinimizationReport,
}
impl<Model, Rhs> FitResult<Model, Rhs>
where
Model: SeparableNonlinearModel<ScalarType = f64>,
Rhs: RhsType,
{
/// this returns the linear coefficients in matrix form. For a single
/// right hand side, this is a matrix with one column. Typically,
/// the [`linear_coefficients`](Self::linear_coefficients) function should be called to get the
/// appropriate view for single or multiple right hand sides out of
/// the box. This function exists to overcome some limitations on how
/// generics work.
pub fn linear_coefficients_generic(&self) -> Option<DMatrixView<'_, Model::ScalarType>> {
self.linear_coefficients.as_ref().map(|c| c.as_view())
}
}
impl<Model> FitResult<Model, MultiRhs>
// take trait bounds from above:
where
Model: SeparableNonlinearModel<ScalarType = f64>,
{
/// **Note** This implementation is for fitting problems with multiple right hand sides.
///
/// Convenience function to get the linear coefficients after the fit has
/// finished. Will return None if there was an error during fitting.
///
/// The coefficients vectors for the individual
/// members of the datasets are the colums of the returned matrix. That means
/// one coefficient vector for each right hand side.
pub fn linear_coefficients(&self) -> Option<DMatrixView<'_, Model::ScalarType>> {
self.linear_coefficients.as_ref().map(|c| c.as_view())
}
/// **Note** This implementation is for fitting problems with multiple right hand sides.
///
/// Calculate the values of the model at the best fit parameters.
/// Returns None if there was an error during fitting.
/// Since this is for multiple right hand sides, the output is a matrix
/// where each column corresponds to one right hand side.
pub fn best_fit(&self) -> Option<OMatrix<Model::ScalarType, Dyn, Dyn>> {
let coeff = self.linear_coefficients()?;
let eval = self.problem.model().eval().ok()?;
Some(eval * coeff)
}
}
impl<Model> FitResult<Model, SingleRhs>
// take trait bounds from above:
where
Model: SeparableNonlinearModel<ScalarType = f64>,
{
/// **Note** This implementation is for fitting problems with a single right hand side.
///
/// Convenience function to get the linear coefficients after the fit has
/// finished. Will return None if there was an error during fitting.
/// The coefficients are given as a single vector.
pub fn linear_coefficients(&self) -> Option<DVectorView<'_, Model::ScalarType>> {
let coeff = self.linear_coefficients.as_ref()?;
debug_assert_eq!(
coeff.ncols(),
1,
"Coefficient matrix must have exactly one colum for problem with single right hand side. This indicates a programming error inside this library!"
);
Some(coeff.column(0))
}
/// **Note** This implementation is for fitting problems with multiple right hand sides
///
/// Calculate the values of the model at the best fit parameters.
/// Returns None if there was an error during fitting.
/// Since this is for a single right hand side, the output is a vector.
pub fn best_fit(&self) -> Option<OVector<Model::ScalarType, Dyn>> {
let coeff = self.linear_coefficients()?;
let eval = self.problem.model().eval().ok()?;
Some(eval * coeff)
}
/// Return all fitted parameters in a single map.
///
/// Nonlinear parameters keep their original symbolic names. Linear
/// coefficients are exposed as `c0`, `c1`, ... in basis-function order.
pub fn solution_map(&self) -> Option<HashMap<String, f64>> {
let nonlinear_names = self.problem.model().parameter_names();
let nonlinear_values = self.nonlinear_parameters();
let linear_coefficients = self.linear_coefficients()?;
let mut map = HashMap::with_capacity(nonlinear_names.len() + linear_coefficients.len());
for (name, value) in nonlinear_names.iter().zip(nonlinear_values.iter()) {
map.insert(name.clone(), *value);
}
for (index, value) in linear_coefficients.iter().enumerate() {
map.insert(format!("c{index}"), *value);
}
Some(map)
}
/// Return the coefficient of determination for the fitted data.
///
/// This uses the problem data stored by the solver. For unweighted
/// problems this matches the ordinary `r^2` value.
pub fn r_squared(&self) -> Option<f64> {
let observed: Vec<f64> = self.problem.weighted_data().iter().copied().collect();
let predicted: Vec<f64> = self.best_fit()?.iter().copied().collect();
if observed.len() != predicted.len() || observed.is_empty() {
return None;
}
let mean = observed.iter().sum::<f64>() / observed.len() as f64;
let ss_tot = observed.iter().map(|y| (y - mean).powi(2)).sum::<f64>();
if ss_tot == 0.0 {
return None;
}
let ss_res = observed
.iter()
.zip(predicted.iter())
.map(|(y, y_pred)| (y - y_pred).powi(2))
.sum::<f64>();
Some(1.0 - ss_res / ss_tot)
}
/// Short alias for [`r_squared`](Self::r_squared).
pub fn r2(&self) -> Option<f64> {
self.r_squared()
}
}
impl<Model, Rhs: RhsType> FitResult<Model, Rhs>
// take trait bounds from above:
where
Model: SeparableNonlinearModel<ScalarType = f64>,
{
/// internal helper for constructing an instance
pub(crate) fn new(
problem: SeparableProblem<Model, Rhs>,
linear_coefficients: Option<DMatrix<Model::ScalarType>>,
minimization_report: MinimizationReport,
) -> Self {
Self {
problem,
minimization_report,
linear_coefficients,
}
}
/// convenience function to get the nonlinear parameters of the model after
/// the fitting process has finished.
pub fn nonlinear_parameters(&self) -> OVector<Model::ScalarType, Dyn> {
self.problem.model().params()
}
/// whether the fit was deemeed successful. The fit might still be not
/// be optimal for numerical reasons, but the minimization process
/// terminated successfully.
pub fn was_successful(&self) -> bool {
self.minimization_report.termination.was_successful()
}
}