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

use super::{
    functionify, 
    split_hm, 
    linalg::NxN
};
use std::collections::HashMap;

#[derive(Clone)]
#[derive(Debug)]
pub struct Variable {
    value: f64,
    domain: Option<[f64; 2]>
}
impl Variable {
    /// Instantiates a new `Variable` struct with a specified value and domain
    pub fn new(value: f64, domain:Option<[f64; 2]>) -> Variable {
        Variable {
            value, 
            domain,
        }
    }

    /// Allows the ability to mutate `self.value` if the new value is on `self.domain`
    pub fn change(&mut self, qty: f64) {
        match self.domain {
            Some(bounds) => {
                if bounds[0] < qty && qty < bounds[1] {
                    self.value = qty;
                } else if bounds[0] > qty { // if qty is o.o.b., then move self.value to the bound 
                    self.value = bounds[0];
                } else {
                    self.value = bounds[1];
                }
            }
            None => {
                // This comment is here exclusively to commemorate the STUPIDEST bug I have ever written:
                // self.value += qty; <- note how the variable's value is increased instead of changed
                //            ~~         if no domain is specified. 
                self.value = qty;
            }
        }
    }

    /// returns `self.value` as `f64`
    pub fn as_f64(&self) -> f64 {
        self.value
    }
}

/// Returns the derivative of a function at a point
fn d_dx(mut func: impl FnMut(f64) -> f64, x: f64) -> f64 {
    let dx = 1e-7;
    ( func(x + dx) - func(x) ) / dx
}

/// Returns the partial derivative of a function w.r.t. the `target` variable
/// # Example
/// ```
/// use nexsys_core::mvcalc::partial_d_dx;
/// use nexsys_core::mvcalc::Variable;
/// use std::collections::HashMap;
/// let expr = "x^2 + y - z";
/// 
/// let X = HashMap::from([
///     ("x", Variable::new(1_f64, None)),
///     ("y", Variable::new(1_f64, None)),
///     ("z", Variable::new(1_f64, None))
/// ]);
/// 
/// let dFdx = partial_d_dx(expr, &X, "x");
/// assert_eq!(dFdx.round(), 2_f64);
/// ```
pub fn partial_d_dx<'a>(
    expr: &str, 
    guess: &HashMap<&'a str, Variable>, 
    target: &'a str 
) -> f64 {

    // copy the guess vector
    let mut temp = guess.clone();

    // create an actual function from the given expression
    let func = functionify(expr);

    // create a partial function of the target variable
    let partial = move |x:f64| -> f64 {
        if let Some(v) = temp.get_mut(target) {
            v.change(x);
        }
        func(&temp)
    };

    // take the derivative of the partial function
    d_dx(partial, guess[target].as_f64())
}

/// Returns the NxN Jacobian matrix of  
pub fn jacobian<'a>(system: &Vec<&str>, guess: &HashMap<&str, Variable>) -> NxN {

    if system.len() != guess.keys().len() {
        // guard clause against invalid problems
        panic!("ERR: System is not properly constrained!")
    } 
    let size = system.len();
    
    let mut mat = Vec::new();
    let vec = split_hm(guess.clone());

    for c in 0..size {
        let col = Vec::from_iter(
            system.iter().map(
                |i| partial_d_dx(i, guess, vec.0[c])
            )
        );
        mat.push(col);
    };

    NxN::from_cols( size, mat, Some(vec.0) )
}