flavio 0.5.0

flavio welcomes you
Documentation 
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

#[cfg(test)]
mod test;

use super::{
    super::{Hessian, Tensor, TensorRank0},
    Dirichlet, Neumann, OptimizeError, SecondOrder,
};
use crate::ABS_TOL;
use std::ops::Div;

/// The Newton-Raphson method.
#[derive(Debug)]
pub struct NewtonRaphson {
    /// Absolute error tolerance.
    pub abs_tol: TensorRank0,
    /// Whether to check if solution is minimum.
    pub check_minimum: bool,
    /// Maximum number of steps.
    pub max_steps: usize,
}

impl Default for NewtonRaphson {
    fn default() -> Self {
        Self {
            abs_tol: ABS_TOL,
            check_minimum: true,
            max_steps: 250,
        }
    }
}

impl<H: Hessian, J: Tensor, X: Tensor> SecondOrder<H, J, X> for NewtonRaphson
where
    J: Div<H, Output = X>,
{
    fn minimize(
        &self,
        jacobian: impl Fn(&X) -> Result<J, OptimizeError>,
        hessian: impl Fn(&X) -> Result<H, OptimizeError>,
        initial_guess: X,
        _dirichlet: Option<Dirichlet>,
        _neumann: Option<Neumann>,
    ) -> Result<X, OptimizeError> {
        //
        // consider having separate implementations for constrained and unconstrained optimization
        // mostly because what you are planning for constrained optimization isn't really NewtonRaphson anymore
        // so would have unconstrained GD, unconstrained NR, and the constrained solver for fem
        // maybe you can pass Option<SecondOrder> into that solver?
        //
        let mut residual;
        let mut solution = initial_guess;
        // if let Some(ref bc) = dirichlet {
        //     bc.places
        //         .iter()
        //         .zip(bc.values.iter())
        //         .for_each(|(place, value)| *solution.get_at_mut(place) = *value)
        // }
        let mut tangent;
        for _ in 0..self.max_steps {
            residual = jacobian(&solution)?;
            // if let Some(ref bc) = neumann {
            //     bc.places
            //         .iter()
            //         .zip(bc.values.iter())
            //         .for_each(|(place, value)| *residual.get_at_mut(place) -= value)
            // }
            // if let Some(ref bc) = dirichlet {
            //     bc.places
            //         .iter()
            //         .for_each(|place| *residual.get_at_mut(place) = 0.0)
            // }
            tangent = hessian(&solution)?;
            if residual.norm() < self.abs_tol {
                if self.check_minimum && !tangent.is_positive_definite() {
                    return Err(OptimizeError::NotMinimum(
                        format!("{}", solution),
                        format!("{:?}", &self),
                    ));
                } else {
                    return Ok(solution);
                }
            } else {
                solution -= residual / tangent;
            }
        }
        Err(OptimizeError::MaximumStepsReached(
            self.max_steps,
            format!("{:?}", &self),
        ))
    }
}