boolean search()

in commons-statistics-inference/src/main/java/org/apache/commons/statistics/inference/BracketFinder.java [106:222]
91 lines of code
14 McCabe index (conditional complexity)

    boolean search(DoubleUnaryOperator func,
                   double a, double b,
                   double min, double max) {
        evaluations = 0;

        // Limit the range of x
        final DoubleUnaryOperator range;
        if (min < max) {
            // Limit: min <= x <= max
            range = x -> {
                if (x > min) {
                    return x < max ? x : max;
                }
                return min;
            };
        } else {
            range = DoubleUnaryOperator.identity();
        }

        double xA = range.applyAsDouble(a);
        double xB = range.applyAsDouble(b);
        double fA = value(func, xA);
        double fB = value(func, xB);
        // Ensure fB <= fA
        if (fA < fB) {
            double tmp = xA;
            xA = xB;
            xB = tmp;
            tmp = fA;
            fA = fB;
            fB = tmp;
        }

        double xC = range.applyAsDouble(xB + GOLD * (xB - xA));
        double fC = value(func, xC);

        // Note: When a [min, max] interval is provided and there is no minima then this
        // loop will terminate when B == C and both are at the min/max bound.
        while (fC < fB) {
            final double tmp1 = (xB - xA) * (fB - fC);
            final double tmp2 = (xB - xC) * (fB - fA);

            final double val = tmp2 - tmp1;
            // limit magnitude of val to a small value
            final double denom = 2 * Math.copySign(Math.max(Math.abs(val), EPS_MIN), val);

            double w = range.applyAsDouble(xB - ((xB - xC) * tmp2 - (xB - xA) * tmp1) / denom);
            final double wLim = range.applyAsDouble(xB + growLimit * (xC - xB));

            double fW;
            if ((w - xC) * (xB - w) > 0) {
                // xB < w < xC
                fW = value(func, w);
                if (fW < fC) {
                    // minimum in [xB, xC]
                    xA = xB;
                    xB = w;
                    fA = fB;
                    fB = fW;
                    break;
                } else if (fW > fB) {
                    // minimum in [xA, w]
                    xC = w;
                    fC = fW;
                    break;
                }
                // continue downhill
                w = range.applyAsDouble(xC + GOLD * (xC - xB));
                fW = value(func, w);
            } else if ((w - wLim) * (xC - w) > 0) {
                // xC < w < limit
                fW = value(func, w);
                if (fW < fC) {
                    // continue downhill
                    xB = xC;
                    xC = w;
                    w = range.applyAsDouble(xC + GOLD * (xC - xB));
                    fB = fC;
                    fC = fW;
                    fW = value(func, w);
                }
            } else if ((w - wLim) * (wLim - xC) >= 0) {
                // xC <= limit <= w
                w = wLim;
                fW = value(func, w);
            } else {
                // possibly w == xC; reject w and take a default step
                w = range.applyAsDouble(xC + GOLD * (xC - xB));
                fW = value(func, w);
            }

            xA = xB;
            fA = fB;
            xB = xC;
            fB = fC;
            xC = w;
            fC = fW;
        }

        mid = xB;
        fMid = fB;

        // Store the bracket: lo <= mid <= hi
        if (xC < xA) {
            lo = xC;
            fLo = fC;
            hi = xA;
            fHi = fA;
        } else {
            lo = xA;
            fLo = fA;
            hi = xC;
            fHi = fC;
        }

        return lo < mid && mid < hi;
    }