1d newton

[ryanelandt]

Closes #808

This PR adds Newton-aided bisection (tries both directions) as a fallback if local minimization fails. This provides reasonable protection against bad initial guesses. If a root is bracketed, it provides the guaranteed robustness of rtsafe. It also adds the ability to find roots that are difficult to or cannot be bracketed (e.g., cusps). Solver specialization is done without branching and without the need for additional function calls.

The solver makes three attempts to solve find a root. As summarized in the code comment below:

   // Makes three attempts to find a root:
   //   1. Local minimization starting with the initial x
   //   2. Bracketing the root in the direction of the initial dx if there is a sign flip
   //      between f(x) and f(x_bound in direction of dx)
   //   3. Bracketing the root in the opposite direction of the initial dx if there is a
   //      sign flip between f(x) and f(x_bound in opposite direction of dx).
   //
   // If all three attempts fail, an error is thrown.

The initial local minimization attempt transitions between up to three solvers depending on particular sets of circumstances. The final two attempts use the Case 2B solver (assuming a root is bracketed).

   ////// 1D Newton root finder explanation. //////
   //
   // Consider the following three root finding problems:
   //
   //   Case 1:  Find a root of the function f(x) given bounds [x_l_orig, x_h_orig] and
   //            an initial evaluation (x_0, f(x_0), f'(x_0)).
   //   Case 2U: Find a root of the function f(x) given evaluations at a low point
   //            (x_l, f(x_l), f'(x_l)) and a high point (x_h, f(x_h), f'(x_h)), where
   //            and f(x_l) and f(x_h) have the same sign, but have slopes that indicate
   //            that the function is sloping towards zero between x_l and x_h, forming
   //            a U shape (hence the name). The U may or may not cross zero.
   //   Case 2B: Find a root of the function f(x) given evaluations at a low point
   //            (x_l, f(x_l), f'(x_l)) and a high point (x_h, f(x_h), f'(x_h)), where
   //            and f(x_l) and f(x_h) have opposite signs that BRACKET a root.
   //
   // All three cases operate on similar data, but require different algorithms to
   // solve. These cases can also transitions between each other. Specifically...
   //
   //   Case 1:  The problem will remain in Case 1 as long as dx does not change sign.
   //            If dx changes sign, the problem will transition to Case 2U or Case 2B.
   //   Case 2U: The problem will transition to Case 2B if the criteria for Case 2B is
   //            met, i.e., f(x_l) and f(x_h) have opposite signs.
   //   Case 2B: Does not transition to other cases.
   //
   //
   ////// Code organization. //////
   //
   // The code for the root finding problem is divided into two parts: a state that
   // contains data and solvers that are dataless. State is stored in the Root1D_State
   // class. The following solver classes operate on the state:
   //      - SolverCase1
   //      - SolverCase2U
   //      - SolverCase2B
   // The process of finding a root begins by iterating with the solver for Case 1.
   // Iterations continue until one of four following things happen:
   //   - Returns success (f(x) is zero, ... etc)
   //   - Transition to Case 2U
   //   - Transition to Case 2B
   //   - Evaluated limit without being able to transition to Case 2U or Case 2B
   // Similarly, when iterating in Case 2U, iterations will continue until one of the
   // following happens:
   //   - Transition to Case 2B
   //   - Returns success (f(x) is zero, ... etc)
   //   - Returns failure (converging to non-zero extrema)
   // When iterating in Case 2B, iterations will continue until:
   //    - Returns success (f(x) is zero, ... etc)

[ryanelandt]

The remaining failures seems unrelated to this PR. Let me know if there are any changes I should make.

[NAThompson]

@ryanelandt : Will try to get to this this weekend! First impressions are good though.

[ryanelandt]

@NAThompson Thanks! I look forward to your feedback.

[ryanelandt]

@NAThompson Any updates on possible changes?

[NAThompson]

@ryanelandt : I'm so sorry! I'm really trying to get to this!

(FYI: I hate it when my PRs don't get reviewed-I feel your pain.)

[NAThompson]

Approved the CI run, but this is quite a complicated PR. I think @jzmaddock (who wrote the rootfinder) may need to take a look.

[ryanelandt]

Thanks for looking at this @NAThompson. I look forward to what @jzmaddock has to say.

[mborland]

@ryanelandt can you please rebase/pull develop into this now that the other chunks have been merged?

[ryanelandt]

@mborland after doing this I'm getting test failures. These failures are the result of the 1D Newton solver running out of iterations here:

math/test/test_roots.cpp

Lines 657 to 659 in 50ef83a

	// bug reports:
	boost::math::skew_normal_distribution<> dist(2.0, 1.0, -2.5);
	BOOST_CHECK(boost::math::isfinite(quantile(dist, 0.075)));

Running out of iterations in this test is a bug introduced by the removal of large step protection in #1002. This bug isn't obvious on the develop branch because running out of iterations is treated as a success criteria. #1008 describes this issue in greater detail.

I'm trying to figure out the best way to address this for this PR and will update when I figure something out.

[ryanelandt]

@mborland Thanks for your patience. I found a solution to #1008. I made this another chunk for this PR. Please see #1012.