nexsys_math/
mvcalc.rs

1use super::{NxN, Variable};
2use meval::{Context, eval_str_with_context};
3use std::{collections::HashMap, ops::{Add, Sub, Div}, fmt::Display, hash::Hash};
4
5/// Takes a mathematical expression given as a string and returns a function.
6pub fn functionify<S>(text: S) -> impl Fn(&HashMap<S, Variable>) -> f64
7where
8    S: Copy + AsRef<str> + Display + Into<String>
9{
10    let func = move |v:&HashMap<S, Variable>| -> f64 {
11        let mut ctx = Context::new();
12        
13        for k in v {
14            ctx.var(*k.0, k.1.as_f64());
15        }
16
17        eval_str_with_context(text, ctx).expect(&format!("failed to evaluate expression: {}", text))
18    };
19    func
20}
21
22/// Returns the derivative of a function at a point.
23pub fn d_dx<T>(mut func: impl FnMut(T) -> T, x: T) -> T 
24where
25    T: Copy + Add<T, Output = T> + Add<f64, Output = T> + Sub<T, Output = T> + Div<f64, Output = T>
26{
27    let dx = 1e-7;
28    ( func(x + dx) - func(x) ) / dx
29}
30
31/// Returns the partial derivative of a function w.r.t. the `target` variable.
32/// # Example
33/// ```
34/// use nexsys_math::partial_d_dx;
35/// use nexsys_math::Variable;
36/// use std::collections::HashMap;
37/// let expr = "x^2 + y - z";
38/// 
39/// let X = HashMap::from([
40///     ("x", Variable::new(1_f64, None)),
41///     ("y", Variable::new(1_f64, None)),
42///     ("z", Variable::new(1_f64, None))
43/// ]);
44/// 
45/// let dFdx = partial_d_dx(expr, &X, "x");
46/// assert_eq!(dFdx.round(), 2_f64);
47/// ```
48pub fn partial_d_dx<S>(expr: S, guess: &HashMap<S, Variable>, target: S) -> f64 
49where 
50    S: Copy + AsRef<str> + Display + Into<String> + Eq + Hash
51{
52    // copy the guess vector
53    let mut temp = guess.clone();
54
55    // create an actual function from the given expression
56    let func = functionify(expr);
57
58    // create a partial function of the target variable
59    let partial = move |x:f64| -> f64 {
60        if let Some(v) = temp.get_mut(&target) {
61            v.change(x);
62        }
63        func(&temp)
64    };
65
66    // take the derivative of the partial function
67    d_dx(partial, guess[&target].as_f64())
68}
69
70/// Returns a tuple of `Vec`s that contain the keys and values of the original HashMap. 
71/// The index of the key will be the same as its corresponding value's index.
72/// 
73/// This function only exists for use in `pub fn jacobian()`.
74fn split_hm<K, V>(hm: HashMap<K, V>) -> (Vec<K>, Vec<V>) {
75    let mut keys = Vec::new();
76    let mut vals = Vec::new();
77
78    for i in hm {
79        keys.push(i.0);
80        vals.push(i.1);
81    }
82
83    (keys, vals)
84}
85
86/// Returns the (numerical) `NxN` Jacobian matrix of a given system of equations at the vector given by `guess`.
87/// 
88/// Note that the resulting matrix's columns will be in a random order, so extra care is needed to identify which
89/// variable occupies which column by checking the ordering of `self.vars`.
90/// # Example
91/// ```
92/// use nexsys_math::jacobian;
93/// use nexsys_math::Variable;
94/// use std::collections::HashMap;
95/// 
96/// let my_sys = vec![
97///     "x^2 + y",
98///     "y   - x"
99/// ];
100/// let guess = HashMap::from([
101///     ("x", Variable::new(1.0, None)),
102///     ("y", Variable::new(1.0, None))
103/// ]);
104/// 
105/// let j = jacobian(&my_sys, &guess);
106/// 
107/// // j.to_vec() will return roughly:
108/// // vec![
109/// //      vec![2.0, -1.0],
110/// //      vec![1.0, 1.0]
111/// // ];
112/// ```
113pub fn jacobian(system: &Vec<&str>, guess: &HashMap<&str, Variable>) -> Result<NxN, &'static str> {
114    if system.len() != guess.keys().len() { 
115        panic!("ERR: System is not properly constrained!") // guard clause against invalid problems
116    } 
117
118    let size = system.len();
119    let mut mat = Vec::new();
120    let vec = split_hm(guess.clone());
121
122    for c in 0..size {
123        let col = Vec::from_iter(
124            system.iter().map(
125                |&i| partial_d_dx(i, guess, vec.0[c])
126            )
127        );
128        mat.push(col);
129    };
130
131    NxN::from_cols( mat, Some(vec.0) )
132}
nexsys_math/mvcalc.rs

nexsys_math/
mvcalc.rs
