static void minimize()

in include/maths/common/CSolvers.h [109:207]
83 lines of code
22 McCabe index (conditional complexity)

    static void minimize(double a,
                         double b,
                         double fa,
                         double fb,
                         const F& f,
                         double tolerance,
                         std::size_t& maxIterations,
                         double lb,
                         double& x,
                         double& fx) {
        tolerance = std::max(tolerance, std::sqrt(std::numeric_limits<double>::epsilon()));
        const double golden = 0.3819660;

        if (fa < fb) {
            x = a;
            fx = fa;
        } else {
            x = b;
            fx = fb;
        }

        double w = x;
        double v = x;
        double fw = fx;
        double fv = fx;
        double s = 0.0;
        double sLast = 0.0;

        std::size_t n = std::max(maxIterations, std::size_t{1});

        do {
            double xm = bisect(a, b);

            double t1 = tolerance * (std::fabs(x) + 0.25);
            double t2 = 2.0 * t1;
            if (fx <= lb || std::fabs(x - xm) <= (t2 - (b - a) / 2.0)) {
                break;
            }

            double sign = (x >= xm) ? -1.0 : 1.0;

            if (sLast > t1) {
                // Fit parabola to abscissa.
                double r = (x - w) * (fx - fv);
                double q = (x - v) * (fx - fw);
                double p = (x - v) * q - (x - w) * r;
                q = 2 * (q - r);
                if (q > 0) {
                    p = -p;
                }
                q = std::fabs(q);

                double td = sLast;
                sLast = std::fabs(s);

                if (std::fabs(p) >= q * td / 2.0 || p <= q * (a - x) || p >= q * (b - x)) {
                    // Minimum not in range or converging too slowly.
                    sLast = (x >= xm) ? x - a : b - x;
                    s = sign * std::max(golden * sLast, t1);
                } else {
                    s = p / q;
                    double u = x + s;
                    if ((u - a) < t2 || (b - u) < t2) {
                        s = sign * t1;
                    }
                }
            } else {
                // Don't have a suitable abscissa so just use golden
                // section in to the larger of [a, x] and [x, b].
                sLast = (x >= xm) ? x - a : b - x;
                s = sign * std::max(golden * sLast, t1);
            }

            double u = x + s;
            double fu = f(u);
            LOG_TRACE(<< "s = " << s << ", u = " << u);
            LOG_TRACE(<< "f(u) = " << fu << ", f(x) = " << fx);

            if (fu <= fx) {
                u >= x ? a = x : b = x;
                shift(v, w, x, u);
                shift(fv, fw, fx, fu);
            } else {
                u < x ? a = u : b = u;
                if (fu <= fw || w == x) {
                    shift(v, w, u);
                    shift(fv, fw, fu);
                } else if (fu <= fv || v == x || v == w) {
                    v = u;
                    fv = fu;
                }
            }
            LOG_TRACE(<< "a = " << a << ", b = " << b);
            LOG_TRACE(<< "x = " << x << ", v = " << v << ", w = " << w);
            LOG_TRACE(<< "f(x) = " << fx << ", f(v) = " << fv << ", f(w) = " << fw);
        } while (--n > 0);

        maxIterations -= n;
    }