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

use structure::matrix::*;
use structure::dual::*;
use util::non_macro::zeros;

/// Jacobian Matrix
/// 
/// # Description
/// : Exact jacobian matrix using Automatic Differenitation
/// 
/// # Type
/// (Vector, F) -> Matrix where F: Fn(Vec<Dual>) -> Vec<Dual>
/// 
/// # Examples
/// ```
/// extern crate peroxide;
/// use peroxide::*;
/// 
/// let x = c!(1, 1);
/// let j = jacobian(x, f);
/// j.print();
/// 
/// //      c[0] c[1]
/// // r[0]    1    1
/// // r[1]   -1    2
/// 
/// fn f(xs: Vec<Dual>) -> Vec<Dual> {
///     let x = xs[0];
///     let y = xs[1];
/// 
///     vec![
///        x - y,
///        x + 2.*y,
///    ]
/// }
/// ```
#[allow(non_snake_case)]
pub fn jacobian<F>(x: Vec<f64>, f: F) -> Matrix
    where F: Fn(Vec<Dual>) -> Vec<Dual>
{
    let l = x.len();

    let x_var: Vec<Dual> = merge_dual(x.clone(), vec![1f64; l]);
    let x_const = x.clone().conv_dual();

    let l2 = f(x_const.clone()).len();

    let mut J = zeros(l2, l);

    let mut x_temp = x_const.clone();

    for i in 0 .. l {
        x_temp[i] = x_var[i];
        let dual_temp = f(x_temp.clone());
        let slope_temp = dual_temp.slopes();
        for j in 0 .. l2 {
            J[(j, i)] = slope_temp[j];
        }
        x_temp = x_const.clone();
    }
    J
}