Repository: scipy/scipy
File: scipy/stats/_levy_stable/__init__.py (Python)

Code before:
    # following Nolan, we may now assume
    #   x0 > zeta when alpha != 1
    #   beta != 0 when alpha == 1

    # spare calculating integral on null set
    # use isclose as macos has fp differences
    if np.isclose(-xi, np.pi / 2, rtol=1e-014, atol=1e-014):
        return 0.0

    def integrand(theta):
        # limit any numerical issues leading to g_1 < 0 near theta limits
        g_1 = g(theta)
        if not np.isfinite(g_1) or g_1 < 0:
            g_1 = 0
        return g_1 * np.exp(-g_1)

    with np.errstate(all="ignore"):
        theta = np.linspace(-xi, np.pi/2, 1000, endpoint=True)
        # We should probably vectorize the function "g" in levyst.c 
        # so that this runs faster
        gvals = np.array([integrand(t) for t in theta])
        gmax = np.max(gvals)
        peak = theta[np.argmax(gvals)]
        try:
            tail_points = [optimize.bisect(lambda t: integrand(t) - gmax / exp_height, -xi, np.pi/2)
                            for exp_height in [1000, 500, 100, 10, 5]]
        except:
            tail_points = np.linspace(-xi, np.pi / 2, 50).tolist()
            
        points = tail_points + [0, peak, np.pi/2]
        
        intg, *ret = integrate.quad(
            integrand,
            -xi,
            np.pi / 2,
            points=points,
            limit=100,
            epsrel=quad_eps,
            epsabs=0,
            full_output=1,
        )

    return c2 * intg


def _cdf_single_value_piecewise_Z1(x, alpha, beta, **kwds):
    # convert from Nolan's S_1 (aka S) to S_0 (aka Zolaterev M)
    # parameterization

    zeta = -beta * np.tan(np.pi * alpha / 2.0)
    x0 = x + zeta if alpha != 1 else x

Review comment (bug):
Besides the addition of `1000` and `500` and the (bare?) `except` block, what has changed? It looks like maybe this is looking at *relative* heights compared to an estimate of the maximum of the integrand?

And when/why does bisection fail? Is the callable not monotonic? Is the root not guaranteed to be within the bounds? Why wasn't `except` needed before and is now? If it *is* needed now, what type of error is raised and can the `except` include that qualification?

Code after:
        # limit any numerical issues leading to g_1 < 0 near theta limits
        g_1 = g(theta)
        if not np.isfinite(g_1) or g_1 < 0:
            g_1 = 0
        return g_1 * np.exp(-g_1)

    with np.errstate(all="ignore"):
        peak = optimize.bisect(
            lambda t: g(t) - 1, -xi, np.pi / 2, xtol=quad_eps
        )

        # this integrand can be very peaked, so we need to force
        # QUADPACK to evaluate the function inside its support
        #

        # lastly, we add additional samples at
        #   ~exp(-100), ~exp(-10), ~exp(-5), ~exp(-1)
        # to improve QUADPACK's detection of rapidly descending tail behavior
        # (this choice is fairly ad hoc)
        tail_points = [
            optimize.bisect(lambda t: g(t) - exp_height, -xi, np.pi / 2)
            for exp_height in [100, 10, 5]
            # exp_height = 1 is handled by peak
        ]
        # Symmetry is used to ensure x > zeta as above. Integration 
        # becomes numerically unstable so let `quad` handle 
        # integration points
        if beta != -1:
            intg_points = [0, peak] + tail_points
        else:
            intg_points = None

        intg, *ret = integrate.quad(
            integrand,
            -xi,
            np.pi / 2,
            points=intg_points,
            limit=100,
            epsrel=quad_eps,
            epsabs=0,
            full_output=1,
        )

    return c2 * intg


def _cdf_single_value_piecewise_Z1(x, alpha, beta, **kwds):
    # convert from Nolan's S_1 (aka S) to S_0 (aka Zolaterev M)
    # parameterization

    zeta = -beta * np.tan(np.pi * alpha / 2.0)