peroxide/numerical/
utils.rs

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

peroxide/numerical/
utils.rs
