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_derivative::derivator::DerivatorMultiVariable;
use crate::numerical_derivative::hessian::Hessian;
use num_complex::ComplexFloat;
#[derive(Debug)]
pub struct QuadraticApproximationResult<T: ComplexFloat, const NUM_VARS: usize>
{
pub intercept: T,
pub linear_coefficients: [T; NUM_VARS],
pub quadratic_coefficients: [[T; NUM_VARS]; NUM_VARS]
}
#[derive(Debug)]
pub struct QuadraticApproximationPredictionMetrics<T: ComplexFloat>
{
pub mean_absolute_error: T::Real,
pub mean_squared_error: T::Real,
pub root_mean_squared_error: T::Real,
pub r_squared: T::Real,
pub adjusted_r_squared: T::Real
}
///Helper functions if you don't care about the details and just want the predictor directly
impl<T: ComplexFloat, const NUM_VARS: usize> QuadraticApproximationResult<T, NUM_VARS>
{
pub fn get_prediction_value(&self, args: &[T; NUM_VARS]) -> T
{
let mut result = self.intercept;
for (i, arg) in args.iter().enumerate().take(NUM_VARS)
{
result = result + self.linear_coefficients[i]**arg;
}
for i in 0..NUM_VARS
{
for j in 1..NUM_VARS
{
result = result + self.quadratic_coefficients[i][j]*args[i]*args[j];
}
}
return result;
}
//get prediction metrics by feeding a list of points and the original function
pub fn get_prediction_metrics<const NUM_POINTS: usize>(&self, points: &[[T; NUM_VARS]; NUM_POINTS], original_function: &dyn Fn(&[T; NUM_VARS]) -> T) -> QuadraticApproximationPredictionMetrics<T>
{
//let num_points = points.len() as f64;
let mut mae = T::zero();
let mut mse = T::zero();
for point in points.iter().take(NUM_POINTS)
{
let predicted_y = self.get_prediction_value(point);
mae = mae + (predicted_y - original_function(point));
mse = mse + num_complex::ComplexFloat::powi(predicted_y - original_function(point), 2);
}
mae = mae/T::from(NUM_POINTS).unwrap();
mse = mse/T::from(NUM_POINTS).unwrap();
let rmse = mse.sqrt().abs();
let mut r2_numerator = T::zero();
let mut r2_denominator = T::zero();
for point in points.iter().take(NUM_POINTS)
{
let predicted_y = self.get_prediction_value(point);
r2_numerator = r2_numerator + num_complex::ComplexFloat::powi(predicted_y - original_function(point), 2);
r2_denominator = r2_numerator + num_complex::ComplexFloat::powi(mae - original_function(point), 2);
}
let r2 = T::one() - (r2_numerator/r2_denominator);
let r2_adj = T::one() - (T::one() - r2)*(T::from(NUM_POINTS).unwrap())/(T::from(NUM_POINTS).unwrap() - T::from(2.0).unwrap());
return QuadraticApproximationPredictionMetrics
{
mean_absolute_error: mae.abs(),
mean_squared_error: mse.abs(),
root_mean_squared_error: rmse,
r_squared: r2.abs(),
adjusted_r_squared: r2_adj.abs()
};
}
}
pub struct QuadraticApproximator<D: DerivatorMultiVariable>
{
derivator: D,
}
impl<D: DerivatorMultiVariable> Default for QuadraticApproximator<D>
{
fn default() -> Self
{
return QuadraticApproximator { derivator: D::default()};
}
}
impl<D: DerivatorMultiVariable> QuadraticApproximator<D>
{
pub fn from_derivator(derivator: D) -> Self
{
return QuadraticApproximator {derivator};
}
/// For an n-dimensional approximation, the equation is approximated as I + L + Q, where:
/// I = intercept
/// L = linear_coefficients[0]*var_1 + linear_coefficients[1]*var_2 + ... + linear_coefficients[n-1]*var_n
/// Q = quadratic_coefficients[0][0]*var_1*var_1 + quadratic_coefficients[0][1]*var_1*var_2 + ... + quadratic_coefficients[n-1][n-1]*var_n*var_n
///
/// NOTE: Returns a Result<T, &'static str>
/// Possible &'static str are:
/// NumberOfStepsCannotBeZero -> if the derivative step size is zero
///
///example function is e^(x/2) + sin(y) + 2.0*z, which we want to approximate. First define the function:
///```
///use multicalc::approximation::quadratic_approximation::*;
///use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;
///
///let function_to_approximate = | args: &[f64; 3] | -> f64
///{
/// return f64::exp(args[0]/2.0) + f64::sin(args[1]) + 2.0*args[2];
///};
///
///let point = [0.0, 1.57, 10.0]; //the point we want to approximate around
///
///let approximator = QuadraticApproximator::<MultiVariableSolver>::default();
///let result = approximator.get(&function_to_approximate, &point).unwrap();
///
///assert!(f64::abs(function_to_approximate(&point) - result.get_prediction_value(&point)) < 1e-9);
/// ```
/// you can also inspect the results of the approximation. For an n-dimensional approximation, the equation is linearized as
///
/// [`QuadraticApproximationResult::intercept`] gives you the required intercept
/// [`QuadraticApproximationResult::linear_coefficients`] gives you the required linear coefficients in order
/// [`QuadraticApproximationResult::quadratic_coefficients`] gives you the required quadratic coefficients as a matrix
///
/// if you don't care about the results and want the predictor directly, use [`QuadraticApproximationResult::get_prediction_value()`]
/// you can also inspect the prediction metrics by providing list of points, use [`QuadraticApproximationResult::get_prediction_metrics()`]
///
/// to see how the [QuadraticApproximationResult::quadratic_coefficients] matrix should be used, refer to [`QuadraticApproximationResult::get_prediction_metrics()`]
/// or refer to its tests.
///
pub fn get<T: ComplexFloat, const NUM_VARS: usize>(&self, function: &dyn Fn(&[T; NUM_VARS]) -> T, point: &[T; NUM_VARS]) -> Result<QuadraticApproximationResult<T, NUM_VARS>, &'static str>
{
let mut intercept_ = function(point);
let mut linear_coeffs_ = [T::zero(); NUM_VARS];
let hessian_matrix = Hessian::from_derivator(self.derivator).get(function, point)?;
for iter in 0..NUM_VARS
{
linear_coeffs_[iter] = self.derivator.get(1, function, &[iter], point)?;
intercept_ = intercept_ - self.derivator.get(1, function, &[iter], point)?*point[iter];
}
let mut quad_coeff = [[T::zero(); NUM_VARS]; NUM_VARS];
for row in 0..NUM_VARS
{
for col in row..NUM_VARS
{
quad_coeff[row][col] = hessian_matrix[row][col];
intercept_ = intercept_ + hessian_matrix[row][col]*point[row]*point[row];
if row == col
{
linear_coeffs_[row] = linear_coeffs_[row] - T::from(2.0).unwrap()*hessian_matrix[row][col]*point[row]
}
else
{
linear_coeffs_[row] = linear_coeffs_[row] - hessian_matrix[row][col]*point[col];
linear_coeffs_[col] = linear_coeffs_[col] - hessian_matrix[row][col]*point[row];
}
}
}
return Ok(QuadraticApproximationResult
{
intercept: intercept_,
linear_coefficients: linear_coeffs_,
quadratic_coefficients: quad_coeff
});
}
}