def test_zero_der_nz_dp():
    """Test secant method with a non-zero dp, but an infinite newton step"""
    # pick a symmetrical functions and choose a point on the side that with dx
    # makes a secant that is a flat line with zero slope, EG: f = (x - 100)**2,
    # which has a root at x = 100 and is symmetrical around the line x = 100
    # we have to pick a really big number so that it is consistently true
    # now find a point on each side so that the secant has a zero slope
    dx = np.finfo(float).eps ** 0.33
    # 100 - p0 = p1 - 100 = p0 * (1 + dx) + dx - 100
    # -> 200 = p0 * (2 + dx) + dx
    p0 = (200.0 - dx) / (2.0 + dx)
    with suppress_warnings() as sup:
        sup.filter(RuntimeWarning, "RMS of")
        x = zeros.newton(lambda y: (y - 100.0)**2, x0=[p0] * 10)
    assert_allclose(x, [100] * 10)
    # test scalar cases too
    p0 = (2.0 - 1e-4) / (2.0 + 1e-4)
    with suppress_warnings() as sup:
        sup.filter(RuntimeWarning, "Tolerance of")
        x = zeros.newton(lambda y: (y - 1.0) ** 2, x0=p0)
    assert_allclose(x, 1)
    p0 = (-2.0 + 1e-4) / (2.0 + 1e-4)
    with suppress_warnings() as sup:
        sup.filter(RuntimeWarning, "Tolerance of")
        x = zeros.newton(lambda y: (y + 1.0) ** 2, x0=p0)
    assert_allclose(x, -1)

def test_moments(distname, arg, normalization_ok, higher_ok):
    try:
        distfn = getattr(stats, distname)
    except TypeError:
        distfn = distname
        distname = 'rv_histogram_instance'

    with suppress_warnings() as sup:
        sup.filter(IntegrationWarning, "The integral is probably divergent, or slowly convergent.")
        m, v, s, k = distfn.stats(*arg, moments='mvsk')

    if normalization_ok:
        check_normalization(distfn, arg, distname)

    if higher_ok:
        check_mean_expect(distfn, arg, m, distname)
        with suppress_warnings() as sup:
            sup.filter(IntegrationWarning,
                       "The integral is probably divergent, or slowly convergent.")
            check_skew_expect(distfn, arg, m, v, s, distname)
            check_var_expect(distfn, arg, m, v, distname)
            check_kurt_expect(distfn, arg, m, v, k, distname)

    check_loc_scale(distfn, arg, m, v, distname)
    check_moment(distfn, arg, m, v, distname)

    def test_imsave(self):
        picdir = os.path.join(datapath, "data")
        for png in glob.iglob(picdir + "/*.png"):
            with suppress_warnings() as sup:
                # PIL causes a Py3k ResourceWarning
                sup.filter(message="unclosed file")
                img = misc.imread(png)
            tmpdir = tempfile.mkdtemp()
            try:
                fn1 = os.path.join(tmpdir, 'test.png')
                fn2 = os.path.join(tmpdir, 'testimg')
                with suppress_warnings() as sup:
                    # PIL causes a Py3k ResourceWarning
                    sup.filter(message="unclosed file")
                    misc.imsave(fn1, img)
                    misc.imsave(fn2, img, 'PNG')

                with suppress_warnings() as sup:
                    # PIL causes a Py3k ResourceWarning
                    sup.filter(message="unclosed file")
                    data1 = misc.imread(fn1)
                    data2 = misc.imread(fn2)
                assert_allclose(data1, img)
                assert_allclose(data2, img)
                assert_equal(data1.shape, img.shape)
                assert_equal(data2.shape, img.shape)
            finally:
                shutil.rmtree(tmpdir)

    def test_callback(self):

        def store_residual(r, rvec):
            rvec[rvec.nonzero()[0].max()+1] = r

        # Define, A,b
        A = csr_matrix(array([[-2,1,0,0,0,0],[1,-2,1,0,0,0],[0,1,-2,1,0,0],[0,0,1,-2,1,0],[0,0,0,1,-2,1],[0,0,0,0,1,-2]]))
        b = ones((A.shape[0],))
        maxiter = 1
        rvec = zeros(maxiter+1)
        rvec[0] = 1.0
        callback = lambda r:store_residual(r, rvec)
        with suppress_warnings() as sup:
            sup.filter(DeprecationWarning, ".*called without specifying.*")
            x,flag = gmres(A, b, x0=zeros(A.shape[0]), tol=1e-16, maxiter=maxiter, callback=callback)

        # Expected output from Scipy 1.0.0
        assert_allclose(rvec, array([1.0, 0.81649658092772603]), rtol=1e-10)

        # Test preconditioned callback
        M = 1e-3 * np.eye(A.shape[0])
        rvec = zeros(maxiter+1)
        rvec[0] = 1.0
        with suppress_warnings() as sup:
            sup.filter(DeprecationWarning, ".*called without specifying.*")
            x, flag = gmres(A, b, M=M, tol=1e-16, maxiter=maxiter, callback=callback)

        # Expected output from Scipy 1.0.0 (callback has preconditioned residual!)
        assert_allclose(rvec, array([1.0, 1e-3 * 0.81649658092772603]), rtol=1e-10)

def test_errprint():
    with suppress_warnings() as sup:
        sup.filter(DeprecationWarning, "`errprint` is deprecated!")
        flag = sc.errprint(True)

    try:
        assert_(isinstance(flag, bool))
        with pytest.warns(sc.SpecialFunctionWarning):
            sc.loggamma(0)
    finally:
        with suppress_warnings() as sup:
            sup.filter(DeprecationWarning, "`errprint` is deprecated!")
            sc.errprint(flag)

 def test_bytescale_cscale_lowhigh(self):
     a = np.arange(10)
     with suppress_warnings() as sup:
         sup.filter(DeprecationWarning)
         actual = misc.bytescale(a, cmin=3, cmax=6, low=100, high=200)
     expected = [100, 100, 100, 100, 133, 167, 200, 200, 200, 200]
     assert_equal(actual, expected)

    def test_triangularity_perturbation(self):
        # Experiment (1) of
        # Awad H. Al-Mohy and Nicholas J. Higham (2012)
        # Improved Inverse Scaling and Squaring Algorithms
        # for the Matrix Logarithm.
        A = np.array([
            [3.2346e-1, 3e4, 3e4, 3e4],
            [0, 3.0089e-1, 3e4, 3e4],
            [0, 0, 3.221e-1, 3e4],
            [0, 0, 0, 3.0744e-1]],
            dtype=float)
        A_logm = np.array([
            [-1.12867982029050462e+00, 9.61418377142025565e+04,
             -4.52485573953179264e+09, 2.92496941103871812e+14],
            [0.00000000000000000e+00, -1.20101052953082288e+00,
             9.63469687211303099e+04, -4.68104828911105442e+09],
            [0.00000000000000000e+00, 0.00000000000000000e+00,
             -1.13289322264498393e+00, 9.53249183094775653e+04],
            [0.00000000000000000e+00, 0.00000000000000000e+00,
             0.00000000000000000e+00, -1.17947533272554850e+00]],
            dtype=float)
        assert_allclose(expm(A_logm), A, rtol=1e-4)

        # Perturb the upper triangular matrix by tiny amounts,
        # so that it becomes technically not upper triangular.
        random.seed(1234)
        tiny = 1e-17
        A_logm_perturbed = A_logm.copy()
        A_logm_perturbed[1, 0] = tiny
        with suppress_warnings() as sup:
            sup.filter(RuntimeWarning, "Ill-conditioned.*")
            A_expm_logm_perturbed = expm(A_logm_perturbed)
        rtol = 1e-4
        atol = 100 * tiny
        assert_(not np.allclose(A_expm_logm_perturbed, A, rtol=rtol, atol=atol))

    def test_multiple_constraint_objects(self):
        fun = lambda x: (x[0] - 1)**2 + (x[1] - 2.5)**2 + (x[2] - 0.75)**2
        x0 = [2, 0, 1]
        coni = []  # only inequality constraints (can use cobyla)
        methods = ["slsqp", "cobyla", "trust-constr"]

        # mixed old and new
        coni.append([{'type': 'ineq', 'fun': lambda x: x[0] - 2 * x[1] + 2},
                     NonlinearConstraint(lambda x: x[0] - x[1], -1, 1)])

        coni.append([LinearConstraint([1, -2, 0], -2, np.inf),
                     NonlinearConstraint(lambda x: x[0] - x[1], -1, 1)])

        coni.append([NonlinearConstraint(lambda x: x[0] - 2 * x[1] + 2, 0, np.inf),
                     NonlinearConstraint(lambda x: x[0] - x[1], -1, 1)])

        for con in coni:
            funs = {}
            for method in methods:
                with suppress_warnings() as sup:
                    sup.filter(UserWarning)
                    result = minimize(fun, x0, method=method, constraints=con)
                    funs[method] = result.fun
            assert_allclose(funs['slsqp'], funs['trust-constr'], rtol=1e-4)
            assert_allclose(funs['cobyla'], funs['trust-constr'], rtol=1e-4)

    def test_bug_6690(self):
        # https://github.com/scipy/scipy/issues/6690
        A_eq = np.array([[0., 0., 0., 0.93, 0., 0.65, 0., 0., 0.83, 0.]])
        b_eq = np.array([0.9626])
        A_ub = np.array([[0., 0., 0., 1.18, 0., 0., 0., -0.2, 0.,
                          -0.22],
                         [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
                         [0., 0., 0., 0.43, 0., 0., 0., 0., 0., 0.],
                         [0., -1.22, -0.25, 0., 0., 0., -2.06, 0., 0.,
                          1.37],
                         [0., 0., 0., 0., 0., 0., 0., -0.25, 0., 0.]])
        b_ub = np.array([0.615, 0., 0.172, -0.869, -0.022])
        bounds = np.array(
            [[-0.84, -0.97, 0.34, 0.4, -0.33, -0.74, 0.47, 0.09, -1.45, -0.73],
             [0.37, 0.02, 2.86, 0.86, 1.18, 0.5, 1.76, 0.17, 0.32, -0.15]]).T
        c = np.array([-1.64, 0.7, 1.8, -1.06, -1.16,
                      0.26, 2.13, 1.53, 0.66, 0.28])

        with suppress_warnings() as sup:
            sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll...")
            sup.filter(OptimizeWarning, "Solving system with option...")
            sol = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq,
                          bounds=bounds, method=self.method,
                          options=self.options)
        _assert_success(sol, desired_fun=-1.191)

    def test_network_flow_limited_capacity(self):
        # A network flow problem with supply and demand at nodes
        # and with costs and capacities along directed edges.
        # http://blog.sommer-forst.de/2013/04/10/
        cost = [2, 2, 1, 3, 1]
        bounds = [
            [0, 4],
            [0, 2],
            [0, 2],
            [0, 3],
            [0, 5]]
        n, p = -1, 1
        A_eq = [
            [n, n, 0, 0, 0],
            [p, 0, n, n, 0],
            [0, p, p, 0, n],
            [0, 0, 0, p, p]]
        b_eq = [-4, 0, 0, 4]

        if self.method == "simplex":
            # Including the callback here ensures the solution can be
            # calculated correctly, even when phase 1 terminated
            # with some of the artificial variables as pivots
            # (i.e. basis[:m] contains elements corresponding to
            # the artificial variables)
            res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
                          method=self.method, options=self.options,
                          callback=lambda x, **kwargs: None)
        else:
            with suppress_warnings() as sup:
                sup.filter(RuntimeWarning, "scipy.linalg.solve\nIll...")
                sup.filter(OptimizeWarning, "A_eq does not appear...")
                res = linprog(c=cost, A_eq=A_eq, b_eq=b_eq, bounds=bounds,
                              method=self.method, options=self.options)
        _assert_success(res, desired_fun=14)

def test_zero_rhs(solver):
    np.random.seed(1234)
    A = np.random.rand(10, 10)
    A = A.dot(A.T) + 10 * np.eye(10)

    b = np.zeros(10)
    tols = np.r_[np.logspace(np.log10(1e-10), np.log10(1e2), 7)]

    for tol in tols:
        with suppress_warnings() as sup:
            sup.filter(DeprecationWarning, ".*called without specifying.*")

            x, info = solver(A, b, tol=tol)
            assert_equal(info, 0)
            assert_allclose(x, 0, atol=1e-15)

            x, info = solver(A, b, tol=tol, x0=ones(10))
            assert_equal(info, 0)
            assert_allclose(x, 0, atol=tol)

            if solver is not minres:
                x, info = solver(A, b, tol=tol, atol=0, x0=ones(10))
                if info == 0:
                    assert_allclose(x, 0)

                x, info = solver(A, b, tol=tol, atol=tol)
                assert_equal(info, 0)
                assert_allclose(x, 0, atol=1e-300)

                x, info = solver(A, b, tol=tol, atol=0)
                assert_equal(info, 0)
                assert_allclose(x, 0, atol=1e-300)

def test_imread_indexed_png():
    # The file `foo3x5x4indexed.png` was created with this array
    # (3x5 is (height)x(width)):
    data = np.array([[[127, 0, 255, 255],
                      [127, 0, 255, 255],
                      [127, 0, 255, 255],
                      [127, 0, 255, 255],
                      [127, 0, 255, 255]],
                     [[192, 192, 255, 0],
                      [192, 192, 255, 0],
                      [0, 0, 255, 0],
                      [0, 0, 255, 0],
                      [0, 0, 255, 0]],
                     [[0, 31, 255, 255],
                      [0, 31, 255, 255],
                      [0, 31, 255, 255],
                      [0, 31, 255, 255],
                      [0, 31, 255, 255]]], dtype=np.uint8)

    filename = os.path.join(datapath, 'data', 'foo3x5x4indexed.png')
    with open(filename, 'rb') as f:
        with suppress_warnings() as sup:
            sup.filter(DeprecationWarning)
            im = misc.imread(f)
    assert_array_equal(im, data)

 def test_bytescale(self):
     x = np.array([0, 1, 2], np.uint8)
     y = np.array([0, 1, 2])
     with suppress_warnings() as sup:
         sup.filter(DeprecationWarning)
         assert_equal(misc.bytescale(x), x)
         assert_equal(misc.bytescale(y), [0, 128, 255])

 def test_bytescale_low_equals_high(self):
     a = np.arange(3)
     with suppress_warnings() as sup:
         sup.filter(DeprecationWarning)
         actual = misc.bytescale(a, low=10, high=10)
     expected = [10, 10, 10]
     assert_equal(actual, expected)

 def test_bytescale_rounding(self):
     a = np.array([-0.5, 0.5, 1.5, 2.5, 3.5])
     with suppress_warnings() as sup:
         sup.filter(DeprecationWarning)
         actual = misc.bytescale(a, cmin=0, cmax=10, low=0, high=10)
     expected = [0, 1, 2, 3, 4]
     assert_equal(actual, expected)

 def test_spherical_jn_inf_complex(self):
     # https://dlmf.nist.gov/10.52.E3
     n = 7
     x = np.array([-inf + 0j, inf + 0j, inf*(1+1j)])
     with suppress_warnings() as sup:
         sup.filter(RuntimeWarning, "invalid value encountered in multiply")
         assert_allclose(spherical_jn(n, x), np.array([0, 0, inf*(1+1j)]))

    def test_integral(self):
        x = [1,1,1,2,2,2,4,4,4]
        y = [1,2,3,1,2,3,1,2,3]
        z = array([0,7,8,3,4,7,1,3,4])

        with suppress_warnings() as sup:
            # This seems to fail (ier=1, see ticket 1642).
            sup.filter(UserWarning, "\nThe required storage space")
            lut = SmoothBivariateSpline(x, y, z, kx=1, ky=1, s=0)

        tx = [1,2,4]
        ty = [1,2,3]

        tz = lut(tx, ty)
        trpz = .25*(diff(tx)[:,None]*diff(ty)[None,:]
                    * (tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum()
        assert_almost_equal(lut.integral(tx[0], tx[-1], ty[0], ty[-1]), trpz)

        lut2 = SmoothBivariateSpline(x, y, z, kx=2, ky=2, s=0)
        assert_almost_equal(lut2.integral(tx[0], tx[-1], ty[0], ty[-1]), trpz,
                            decimal=0)  # the quadratures give 23.75 and 23.85

        tz = lut(tx[:-1], ty[:-1])
        trpz = .25*(diff(tx[:-1])[:,None]*diff(ty[:-1])[None,:]
                    * (tz[:-1,:-1]+tz[1:,:-1]+tz[:-1,1:]+tz[1:,1:])).sum()
        assert_almost_equal(lut.integral(tx[0], tx[-2], ty[0], ty[-2]), trpz)

    def test_breakdown_underdetermined(self):
        # Should find LSQ solution in the Krylov span in one inner
        # iteration, despite solver breakdown from nilpotent A.
        A = np.array([[0, 1, 1, 1],
                      [0, 0, 1, 1],
                      [0, 0, 0, 1],
                      [0, 0, 0, 0]], dtype=float)

        bs = [
            np.array([1, 1, 1, 1]),
            np.array([1, 1, 1, 0]),
            np.array([1, 1, 0, 0]),
            np.array([1, 0, 0, 0]),
        ]

        for b in bs:
            with suppress_warnings() as sup:
                sup.filter(DeprecationWarning, ".*called without specifying.*")
                xp, info = lgmres(A, b, maxiter=1)
            resp = np.linalg.norm(A.dot(xp) - b)

            K = np.c_[b, A.dot(b), A.dot(A.dot(b)), A.dot(A.dot(A.dot(b)))]
            y, _, _, _ = np.linalg.lstsq(A.dot(K), b, rcond=-1)
            x = K.dot(y)
            res = np.linalg.norm(A.dot(x) - b)

            assert_allclose(resp, res, err_msg=repr(b))

def test_reentrancy():
    non_reentrant = [cg, cgs, bicg, bicgstab, gmres, qmr]
    reentrant = [lgmres, minres, gcrotmk]
    for solver in reentrant + non_reentrant:
        with suppress_warnings() as sup:
            sup.filter(DeprecationWarning, ".*called without specifying.*")
            _check_reentrancy(solver, solver in reentrant)

    def test_cornercase(self):
        np.random.seed(1234)

        # Rounding error may prevent convergence with tol=0 --- ensure
        # that the return values in this case are correct, and no
        # exceptions are raised

        for n in [3, 5, 10, 100]:
            A = 2*eye(n)

            with suppress_warnings() as sup:
                sup.filter(DeprecationWarning, ".*called without specifying.*")
                b = np.ones(n)
                x, info = gcrotmk(A, b, maxiter=10)
                assert_equal(info, 0)
                assert_allclose(A.dot(x) - b, 0, atol=1e-14)

                x, info = gcrotmk(A, b, tol=0, maxiter=10)
                if info == 0:
                    assert_allclose(A.dot(x) - b, 0, atol=1e-14)

                b = np.random.rand(n)
                x, info = gcrotmk(A, b, maxiter=10)
                assert_equal(info, 0)
                assert_allclose(A.dot(x) - b, 0, atol=1e-14)

                x, info = gcrotmk(A, b, tol=0, maxiter=10)
                if info == 0:
                    assert_allclose(A.dot(x) - b, 0, atol=1e-14)