peroxide/numerical/
utils.rs

1use crate::structure::ad::AD::*;
2use crate::structure::ad::*;
3use crate::structure::matrix::*;
4use crate::traits::matrix::MatrixTrait;
5use crate::util::non_macro::{cat, zeros};
6
7/// Jacobian Matrix
8///
9/// # Description
10/// : Exact jacobian matrix using Automatic Differenitation
11///
12/// # Type
13/// `(F, &Vec<f64>) -> Matrix where F: Fn(&Vec<AD>) -> Vec<AD>`
14///
15/// # Examples
16/// ```
17/// #[macro_use]
18/// extern crate peroxide;
19/// use peroxide::fuga::*;
20///
21/// fn main() {
22///     let x = c!(1, 1);
23///     let j = jacobian(f, &x);
24///     j.print();
25///
26///     //      c[0] c[1]
27///     // r[0]    1   -1
28///     // r[1]    1    2
29/// }
30///
31/// fn f(xs: &Vec<AD>) -> Vec<AD> {
32///     let x = xs[0];
33///     let y = xs[1];
34///
35///     vec![
36///        x - y,
37///        x + 2.*y,
38///    ]
39/// }
40/// ```
41#[allow(non_snake_case)]
42pub fn jacobian<F: Fn(&Vec<AD>) -> Vec<AD>>(f: F, x: &Vec<f64>) -> Matrix {
43    let l = x.len();
44    let mut x_ad: Vec<AD> = x.iter().map(|&x| AD1(x, 0f64)).collect();
45    let l2 = f(&x_ad).len();
46
47    let mut J = zeros(l2, l);
48
49    for i in 0..l {
50        x_ad[i][1] = 1f64;
51        let slopes: Vec<f64> = f(&x_ad).iter().map(|ad| ad.dx()).collect();
52        J.subs_col(i, &slopes);
53        x_ad[i][1] = 0f64;
54    }
55    J
56}
57
58///// Hessian Matrix
59//#[allow(non_snake_case)]
60//pub fn hessian<F: Fn(&Vec<AD>) -> AD>(f: F, x: &Vec<f64>) -> Matrix {
61//    let l = x.len();
62//    let mut x_ad: Vec<AD> = x.iter().map(|&x| AD2(x, 0f64, 0f64)).collect();
63//
64//    let mut H = zeros(l, l);
65//
66//    for i in 0 .. l {
67//        for j in 0 .. l {
68//        }
69//    }
70//
71//    unimplemented!()
72//}
73
74//#[allow(non_snake_case)]
75//pub fn jacobian_ad<F: Fn(&Vec<AD>) -> Vec<AD>>(f: F, x: &Vec<AD>) -> Vec<Vec<AD>> {
76//    let l = x.len();
77//    let mut x_ad: Vec<AD> = x.clone().into_iter().map(|mut t| {
78//        t.iter_mut().skip(1).for_each(|k| *k = 0f64);
79//        t
80//    }).collect();
81//    let l2 = f(&x_ad).len();
82//
83//    let mut JT: Vec<Vec<AD>> = vec![vec![AD0(0f64); l2]; l];
84//
85//    for i in 0 .. l {
86//        x_ad[i][1] = 1f64;
87//        let ads = f(&x_ad);
88//        JT[i] = ads;
89//        x_ad[i][1] = 0f64;
90//    }
91//    JT
92//}
93
94/// TriDiagonal Matrix Algorithm (TDMA)
95///
96/// # Description
97///
98/// Solve tri-diagonal matrix system efficiently (O(n))
99/// ```bash
100/// |b0 c0         | |x0|   |y0|
101/// |a1 b1 c1      | |x1|   |y1|
102/// |   a2 b2 c2   | |x2| = |y2|
103/// |      ...     | |..|   |..|
104/// |         am bm| |xm|   |ym|
105/// ```
106///
107/// # Caution
108///
109/// You should apply boundary condition yourself
110pub fn tdma(a_input: Vec<f64>, b_input: Vec<f64>, c_input: Vec<f64>, y_input: Vec<f64>) -> Matrix {
111    let n = b_input.len();
112    assert_eq!(a_input.len(), n - 1);
113    assert_eq!(c_input.len(), n - 1);
114    assert_eq!(y_input.len(), n);
115
116    let a = cat(0f64, &a_input);
117    let mut b = b_input.clone();
118    let c = {
119        let mut c_temp = c_input.clone();
120        c_temp.push(0f64);
121        c_temp.clone()
122    };
123    let mut y = y_input.clone();
124
125    // Forward substitution
126    let mut w = vec![0f64; n];
127    for i in 1..n {
128        w[i] = a[i] / b[i - 1];
129        b[i] = b[i] - w[i] * c[i - 1];
130        y[i] = y[i] - w[i] * y[i - 1];
131    }
132
133    // Backward substitution
134    let mut x = vec![0f64; n];
135    x[n - 1] = y[n - 1] / b[n - 1];
136    for i in (0..n - 1).rev() {
137        x[i] = (y[i] - c[i] * x[i + 1]) / b[i];
138    }
139    x.into()
140}
peroxide/numerical/utils.rs

peroxide/numerical/
utils.rs
