Repository: scipy/scipy
File: scipy/optimize/tests/test_bracket.py (Python)

Code before:
        f.count = 0
        res = _bracket_root(f, xp.asarray(5.), xp.asarray(10.), 
                            factor=2)

        assert res.nfev == 4
        xp_assert_close(res.xl, xp.asarray(0.), atol=1e-15)
        xp_assert_close(res.xr, xp.asarray(5.), atol=1e-15)

        # 3. bracket limit hits root exactly
        with np.errstate(over='ignore'):
            res = _bracket_root(f, xp.asarray(5.), xp.asarray(10.), 
                                xmin=0)
        xp_assert_close(res.xl, xp.asarray(0.), atol=1e-15)

        with np.errstate(over='ignore'):
            res = _bracket_root(f, xp.asarray(-10.), xp.asarray(-5.), 
                                xmax=0)
        xp_assert_close(res.xr, xp.asarray(0.), atol=1e-15)

        # 4. bracket not within min, max
        with np.errstate(over='ignore'):
            res = _bracket_root(f, xp.asarray(5.), xp.asarray(10.),
                                xmin=1)
        assert not res.success

        # 5. Edge case. Multiple searches in same direction terminate
        # simultaneously in some iteration after the corresponding searches
        # in the other direction terminated without finding a root.

        # Example based on that in
        # https://github.com/scipy/scipy/pull/22560#discussion_r1962853839
        def f(x, p):
            return np.exp(x) - p

        p = np.asarray([0.29, 0.35])
        res = _bracket_root(f, xl0=-1, xmin=-np.inf, xmax=0, args=(p, ))

        # 6. Edge case: https://github.com/scipy/scipy/pull/22560/files#r1962952517
        def f(x, p, c):
            return np.exp(x*c) - p

        p = [0.32061201, 0.39175242, 0.40047535, 0.50527218, 0.55654373,
             0.11911647, 0.37507896, 0.66554191]
        c = [1., -1., 1., 1., -1., 1., 1., 1.]
        xl0 = [-7.63108551,  3.27840947, -8.36968526, -1.78124372,
               0.92201295, -2.48930123, -0.66733533, -0.44606749]
        xr0 = [-6.63108551,  4.27840947, -7.36968526, -0.78124372,
               1.92201295, -1.48930123, 0., 0.]
        xmin = [-np.inf, 0., -np.inf, -np.inf, 0., -np.inf, -np.inf,
                -np.inf]
        xmax = [0., np.inf, 0., 0., np.inf, 0., 0., 0.]

Review comment (bug):
I wrote these tests before I understood what was going on. I guess these aren't even edge cases, just a straight up bug we missed. I'll separate these out into a different test since they're not special cases.

Code after:
        f.count = 0
        res = _bracket_root(f, xp.asarray(5.), xp.asarray(10.), 
                            factor=2)

        assert res.nfev == 4
        xp_assert_close(res.xl, xp.asarray(0.), atol=1e-15)
        xp_assert_close(res.xr, xp.asarray(5.), atol=1e-15)

        # 3. bracket limit hits root exactly
        with np.errstate(over='ignore'):
            res = _bracket_root(f, xp.asarray(5.), xp.asarray(10.), 
                                xmin=0)
        xp_assert_close(res.xl, xp.asarray(0.), atol=1e-15)

        with np.errstate(over='ignore'):
            res = _bracket_root(f, xp.asarray(-10.), xp.asarray(-5.), 
                                xmax=0)
        xp_assert_close(res.xr, xp.asarray(0.), atol=1e-15)

        # 4. bracket not within min, max
        with np.errstate(over='ignore'):
            res = _bracket_root(f, xp.asarray(5.), xp.asarray(10.),
                                xmin=1)
        assert not res.success

    def test_bug_fixes(self):
        # 1. Bug in double sided bracket search.
        # Happened in some cases where there are terminations on one side
        # after corresponding searches on other side failed due to reaching the
        # boundary.

        # https://github.com/scipy/scipy/pull/22560#discussion_r1962853839
        def f(x, p):
            return np.exp(x) - p

        p = np.asarray([0.29, 0.35])
        res = _bracket_root(f, xl0=-1, xmin=-np.inf, xmax=0, args=(p, ))

        # https://github.com/scipy/scipy/pull/22560/files#r1962952517
        def f(x, p, c):
            return np.exp(x*c) - p

        p = [0.32061201, 0.39175242, 0.40047535, 0.50527218, 0.55654373,
             0.11911647, 0.37507896, 0.66554191]
        c = [1., -1., 1., 1., -1., 1., 1., 1.]
        xl0 = [-7.63108551,  3.27840947, -8.36968526, -1.78124372,
               0.92201295, -2.48930123, -0.66733533, -0.44606749]
        xr0 = [-6.63108551,  4.27840947, -7.36968526, -0.78124372,
               1.92201295, -1.48930123, 0., 0.]
        xmin = [-np.inf, 0., -np.inf, -np.inf, 0., -np.inf, -np.inf,
                -np.inf]