libitofin 0.5.0

A ground-up Rust port of QuantLib: quantitative-finance primitives for pricing, risk, and numerical methods.
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

//! Shared driver for line-search based optimization methods.
//!
//! Port of `ql/math/optimization/linesearchbasedmethod.{hpp,cpp}`. The C++
//! abstract base class' `minimize` becomes a free function taking the
//! `getUpdatedDirection` override as a closure; concrete methods
//! (conjugate gradient, steepest descent) wrap it. One deviation: on the
//! convergence exit the accepted point is stored in the problem, so the
//! reported minimum and its function value belong to the same point;
//! QuantLib returns there without `setCurrentValue`, leaving the problem's
//! current value one iterate behind.

use crate::errors::QlResult;
use crate::math::array::Array;
use crate::math::optimization::endcriteria::{EndCriteria, EndCriteriaType};
use crate::math::optimization::linesearch::LineSearch;
use crate::math::optimization::problem::Problem;
use crate::types::{Real, Size};

/// Minimizes `problem` by repeated line searches, with `updated_direction`
/// computing each new search direction from `(problem, gold2, line_search)`,
/// where `gold2` is the previous squared gradient norm.
pub(crate) fn line_search_based_minimize(
    problem: &mut Problem<'_>,
    end_criteria: &EndCriteria,
    line_search: &mut dyn LineSearch,
    mut updated_direction: impl FnMut(&Problem<'_>, Real, &dyn LineSearch) -> Array,
) -> QlResult<EndCriteriaType> {
    let ftol = end_criteria.function_epsilon();
    let mut max_stationary_state_iterations = end_criteria.max_stationary_state_iterations();
    let mut ec_type = EndCriteriaType::None;
    problem.reset();
    line_search.reset();
    let mut x = problem.current_value().clone();
    let mut iteration_number: Size = 0;
    // classical initial value for the line-search step
    let mut t = 1.0;

    let mut gradient = Array::with_size(x.size());
    let function_value = problem.value_and_gradient(&mut gradient, &x);
    problem.set_function_value(function_value);
    problem.set_gradient_norm_value(gradient.dot(&gradient));
    line_search.set_search_direction(-&gradient);

    loop {
        t = line_search.search(problem, &mut ec_type, end_criteria, t)?;
        // don't fail here: the search can stop just because maxIterations
        // was exceeded
        if !line_search.succeeded() {
            problem.set_current_value(x);
            return Ok(ec_type);
        }

        x = line_search.last_x().clone();
        let fold = problem.function_value();
        problem.set_function_value(line_search.last_function_value());
        // orthogonalization coefficient
        let gold2 = problem.gradient_norm_value();
        problem.set_gradient_norm_value(line_search.last_gradient_norm2());

        let direction = updated_direction(problem, gold2, line_search);
        line_search.set_search_direction(direction);

        // Numerical Recipes exit strategy on f(x) (NR in C++, p.423)
        let fnew = problem.function_value();
        let fdiff = 2.0 * (fnew - fold).abs() / (fnew.abs() + fold.abs() + Real::EPSILON);
        if fdiff < ftol || end_criteria.check_max_iterations(iteration_number, &mut ec_type) {
            end_criteria.check_stationary_function_value(
                0.0,
                0.0,
                &mut max_stationary_state_iterations,
                &mut ec_type,
            );
            end_criteria.check_max_iterations(iteration_number, &mut ec_type);
            problem.set_current_value(x);
            return Ok(ec_type);
        }
        problem.set_current_value(x.clone());
        iteration_number += 1;
    }
}