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
use crate::numerical_derivative::derivator::DerivatorMultiVariable;
use crate::utils::error_codes::*;
use num_complex::ComplexFloat;
#[cfg(feature = "heap")]
use std::{boxed::Box, vec::Vec};
pub struct Jacobian<D: DerivatorMultiVariable>
{
derivator: D
}
impl<D: DerivatorMultiVariable> Default for Jacobian<D>
{
fn default() -> Self
{
return Jacobian { derivator: D::default() };
}
}
impl<D: DerivatorMultiVariable> Jacobian<D>
{
///custom constructor, optimal for fine tuning
/// You can create a custom multivariable derivator from this crate
/// or supply your own by implementing the base traits yourself
pub fn from_derivator(derivator: D) -> Self
{
return Jacobian {derivator}
}
/// Returns the jacobian matrix for a given vector of functions
/// Can handle multivariable functions of any order or complexity
///
/// The 2-D matrix returned has the structure [[df1/dvar1, df1/dvar2, ... , df1/dvarN],
/// [ ... ],
/// [dfM/dvar1, dfM/dvar2, ... , dfM/dvarN]]
///
/// where 'N' is the total number of variables, and 'M' is the total number of functions
///
/// NOTE: Returns a Result<T, &'static str>
/// Possible &'static str are:
/// VECTOR_OF_FUNCTIONS_CANNOT_BE_EMPTY -> if function_matrix argument is an empty array
/// NUMBER_OF_DERIVATIVE_STEPS_CANNOT_BE_ZERO -> if the derivative step size is zero
///
/// assume our function vector is (x*y*z , x^2 + y^2). First define both the functions
/// ```
/// use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;
/// use multicalc::numerical_derivative::jacobian::Jacobian;
/// let my_func1 = | args: &[f64; 3] | -> f64
/// {
/// return args[0]*args[1]*args[2]; //x*y*z
/// };
///
/// let my_func2 = | args: &[f64; 3] | -> f64
/// {
/// return args[0].powf(2.0) + args[1].powf(2.0); //x^2 + y^2
/// };
///
/// //define the function vector
/// let function_matrix: [&dyn Fn(&[f64; 3]) -> f64; 2] = [&my_func1, &my_func2];
/// let points = [1.0, 2.0, 3.0]; //the point around which we want the jacobian matrix
///
/// let jacobian = Jacobian::<MultiVariableSolver>::default();
/// let result = jacobian.get(&function_matrix, &points).unwrap();
/// ```
///
pub fn get<T: ComplexFloat, const NUM_FUNCS: usize, const NUM_VARS: usize>(&self, function_matrix: &[&dyn Fn(&[T; NUM_VARS]) -> T; NUM_FUNCS], vector_of_points: &[T; NUM_VARS]) -> Result<[[T; NUM_VARS]; NUM_FUNCS], &'static str>
{
if function_matrix.is_empty()
{
return Err(VECTOR_OF_FUNCTIONS_CANNOT_BE_EMPTY);
}
let mut result = [[T::zero(); NUM_VARS]; NUM_FUNCS];
for row_index in 0..NUM_FUNCS
{
for col_index in 0..NUM_VARS
{
result[row_index][col_index] = self.derivator.get_single_partial(&function_matrix[row_index], col_index, vector_of_points)?;
}
}
return Ok(result);
}
/// same as [Jacobian::get] but uses heap-allocated vectors to generate the jacobian matrix.
/// Useful when handling large datasets to avoid a stack overflow. To use, turn on the feature "heap" (false by default)
/// Can handle multivariable functions of any order or complexity
///
/// The 2-D matrix returned has the structure [[df1/dvar1, df1/dvar2, ... , df1/dvarN],
/// [ ... ],
/// [dfM/dvar1, dfM/dvar2, ... , dfM/dvarN]]
///
/// where 'N' is the total number of variables, and 'M' is the total number of functions
///
/// NOTE: Returns a Result<T, &'static str>
/// Possible &'static str are:
/// VectorOfFunctionsCannotBeEmpty -> if function_matrix argument is an empty array
/// NumberOfStepsCannotBeZero -> if the derivative step size is zero
#[cfg(feature = "heap")]
pub fn get_on_heap<T: ComplexFloat, const NUM_VARS: usize>(&self, function_matrix: &Vec<Box<dyn Fn(&[T; NUM_VARS]) -> T>>, vector_of_points: &[T; NUM_VARS]) -> Result<Vec<Vec<T>>, &'static str>
{
if function_matrix.is_empty()
{
return Err(VECTOR_OF_FUNCTIONS_CANNOT_BE_EMPTY);
}
let num_funcs = function_matrix.len();
let mut result: Vec<Vec<T>> = Vec::new();
for row_index in 0..num_funcs
{
let mut cur_row: Vec<T> = Vec::new();
for col_index in 0..NUM_VARS
{
cur_row.push(self.derivator.get_single_partial(&function_matrix[row_index], col_index, vector_of_points)?);
}
result.push(cur_row);
}
return Ok(result);
}
}