from __future__ import division, print_function, absolute_import

import sys
import math
import numpy as np
from numpy import sqrt, cos, sin, arctan, exp, log, pi, Inf
from numpy.testing import (assert_, TestCase, run_module_suite, dec,
        assert_allclose, assert_array_less, assert_almost_equal, assert_raises)
from scipy.integrate import quad, dblquad, tplquad, nquad
from scipy._lib.six import xrange
from scipy._lib._ccallback import LowLevelCallable

import ctypes
import ctypes.util
from scipy._lib._ccallback_c import sine_ctypes

import scipy.integrate._test_multivariate as clib_test


def assert_quad(value_and_err, tabled_value, errTol=1.5e-8):
    value, err = value_and_err
    assert_allclose(value, tabled_value, atol=err, rtol=0)
    if errTol is not None:
        assert_array_less(err, errTol)


class TestCtypesQuad(TestCase):
    def setUp(self):
        if sys.platform == 'win32':
            if sys.version_info < (3, 5):
                files = [ctypes.util.find_msvcrt()]
            else:
                files = ['api-ms-win-crt-math-l1-1-0.dll']
        elif sys.platform == 'darwin':
            files = ['libm.dylib']
        else:
            files = ['libm.so', 'libm.so.6']

        for file in files:
            try:
                self.lib = ctypes.CDLL(file)
                break
            except OSError:
                pass
        else:
            # This test doesn't work on some Linux platforms (Fedora for
            # example) that put an ld script in libm.so - see gh-5370
            self.skipTest("Ctypes can't import libm.so")

        restype = ctypes.c_double
        argtypes = (ctypes.c_double,)
        for name in ['sin', 'cos', 'tan']:
            func = getattr(self.lib, name)
            func.restype = restype
            func.argtypes = argtypes

    def test_typical(self):
        assert_quad(quad(self.lib.sin, 0, 5), quad(math.sin, 0, 5)[0])
        assert_quad(quad(self.lib.cos, 0, 5), quad(math.cos, 0, 5)[0])
        assert_quad(quad(self.lib.tan, 0, 1), quad(math.tan, 0, 1)[0])

    @dec.knownfailureif(True, msg="Unreliable test, see ticket 1684.")
    def test_improvement(self):
        import time
        start = time.time()
        for i in xrange(100):
            quad(self.lib.sin, 0, 100)
        fast = time.time() - start
        start = time.time()
        for i in xrange(100):
            quad(math.sin, 0, 100)
        slow = time.time() - start
        assert_(fast < 0.5*slow, (fast, slow))

    def test_ctypes_sine(self):
        quad(LowLevelCallable(sine_ctypes), 0, 1)

    def test_ctypes_variants(self):
        lib = ctypes.CDLL(clib_test.__file__)

        sin_0 = lib._sin_0
        sin_0.restype = ctypes.c_double
        sin_0.argtypes = [ctypes.c_double, ctypes.c_void_p]

        sin_1 = lib._sin_1
        sin_1.restype = ctypes.c_double
        sin_1.argtypes = [ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p]

        sin_2 = lib._sin_2
        sin_2.restype = ctypes.c_double
        sin_2.argtypes = [ctypes.c_double]

        sin_3 = lib._sin_3
        sin_3.restype = ctypes.c_double
        sin_3.argtypes = [ctypes.c_int, ctypes.POINTER(ctypes.c_double)]

        sin_4 = lib._sin_3
        sin_4.restype = ctypes.c_double
        sin_4.argtypes = [ctypes.c_int, ctypes.c_double]

        all_sigs = [sin_0, sin_1, sin_2, sin_3, sin_4]
        legacy_sigs = [sin_2, sin_4]
        legacy_only_sigs = [sin_4]

        # LowLevelCallables work for new signatures
        for j, func in enumerate(all_sigs):
            callback = LowLevelCallable(func)
            if func in legacy_only_sigs:
                assert_raises(ValueError, quad, callback, 0, pi)
            else:
                assert_allclose(quad(callback, 0, pi)[0], 2.0)

        # Plain ctypes items work only for legacy signatures
        for j, func in enumerate(legacy_sigs):
            if func in legacy_sigs:
                assert_allclose(quad(func, 0, pi)[0], 2.0)
            else:
                assert_raises(ValueError, quad, func, 0, pi)


class TestMultivariateCtypesQuad(TestCase):
    def setUp(self):
        self.lib = ctypes.CDLL(clib_test.__file__)
        restype = ctypes.c_double
        argtypes = (ctypes.c_int, ctypes.c_double)
        for name in ['_multivariate_typical', '_multivariate_indefinite',
                     '_multivariate_sin']:
            func = getattr(self.lib, name)
            func.restype = restype
            func.argtypes = argtypes

    def test_typical(self):
        # 1) Typical function with two extra arguments:
        assert_quad(quad(self.lib._multivariate_typical, 0, pi, (2, 1.8)),
                    0.30614353532540296487)

    def test_indefinite(self):
        # 2) Infinite integration limits --- Euler's constant
        assert_quad(quad(self.lib._multivariate_indefinite, 0, Inf),
                    0.577215664901532860606512)

    def test_threadsafety(self):
        # Ensure multivariate ctypes are threadsafe
        def threadsafety(y):
            return y + quad(self.lib._multivariate_sin, 0, 1)[0]
        assert_quad(quad(threadsafety, 0, 1), 0.9596976941318602)

    def test_improvement(self):
        def myfunc(x):           # Euler's constant integrand
            return -exp(-x)*log(x)
        import time
        start = time.time()
        for i in xrange(20):
            quad(self.lib._multivariate_indefinite, 0, 100)
        fast = time.time() - start
        start = time.time()
        for i in xrange(20):
            quad(myfunc, 0, 100)
        slow = time.time() - start
        # 2+ times faster speeds generated by nontrivial ctypes
        # function (single variable)
        assert_(fast < 0.5*slow, (fast, slow))


class TestQuad(TestCase):
    def test_typical(self):
        # 1) Typical function with two extra arguments:
        def myfunc(x, n, z):       # Bessel function integrand
            return cos(n*x-z*sin(x))/pi
        assert_quad(quad(myfunc, 0, pi, (2, 1.8)), 0.30614353532540296487)

    def test_indefinite(self):
        # 2) Infinite integration limits --- Euler's constant
        def myfunc(x):           # Euler's constant integrand
            return -exp(-x)*log(x)
        assert_quad(quad(myfunc, 0, Inf), 0.577215664901532860606512)

    def test_singular(self):
        # 3) Singular points in region of integration.
        def myfunc(x):
            if 0 < x < 2.5:
                return sin(x)
            elif 2.5 <= x <= 5.0:
                return exp(-x)
            else:
                return 0.0

        assert_quad(quad(myfunc, 0, 10, points=[2.5, 5.0]),
                    1 - cos(2.5) + exp(-2.5) - exp(-5.0))

    def test_sine_weighted_finite(self):
        # 4) Sine weighted integral (finite limits)
        def myfunc(x, a):
            return exp(a*(x-1))

        ome = 2.0**3.4
        assert_quad(quad(myfunc, 0, 1, args=20, weight='sin', wvar=ome),
                    (20*sin(ome)-ome*cos(ome)+ome*exp(-20))/(20**2 + ome**2))

    def test_sine_weighted_infinite(self):
        # 5) Sine weighted integral (infinite limits)
        def myfunc(x, a):
            return exp(-x*a)

        a = 4.0
        ome = 3.0
        assert_quad(quad(myfunc, 0, Inf, args=a, weight='sin', wvar=ome),
                    ome/(a**2 + ome**2))

    def test_cosine_weighted_infinite(self):
        # 6) Cosine weighted integral (negative infinite limits)
        def myfunc(x, a):
            return exp(x*a)

        a = 2.5
        ome = 2.3
        assert_quad(quad(myfunc, -Inf, 0, args=a, weight='cos', wvar=ome),
                    a/(a**2 + ome**2))

    def test_algebraic_log_weight(self):
        # 6) Algebraic-logarithmic weight.
        def myfunc(x, a):
            return 1/(1+x+2**(-a))

        a = 1.5
        assert_quad(quad(myfunc, -1, 1, args=a, weight='alg',
                         wvar=(-0.5, -0.5)),
                    pi/sqrt((1+2**(-a))**2 - 1))

    def test_cauchypv_weight(self):
        # 7) Cauchy prinicpal value weighting w(x) = 1/(x-c)
        def myfunc(x, a):
            return 2.0**(-a)/((x-1)**2+4.0**(-a))

        a = 0.4
        tabledValue = ((2.0**(-0.4)*log(1.5) -
                        2.0**(-1.4)*log((4.0**(-a)+16) / (4.0**(-a)+1)) -
                        arctan(2.0**(a+2)) -
                        arctan(2.0**a)) /
                       (4.0**(-a) + 1))
        assert_quad(quad(myfunc, 0, 5, args=0.4, weight='cauchy', wvar=2.0),
                    tabledValue, errTol=1.9e-8)

    def test_double_integral(self):
        # 8) Double Integral test
        def simpfunc(y, x):       # Note order of arguments.
            return x+y

        a, b = 1.0, 2.0
        assert_quad(dblquad(simpfunc, a, b, lambda x: x, lambda x: 2*x),
                    5/6.0 * (b**3.0-a**3.0))

    def test_double_integral2(self):
        def func(x0, x1, t0, t1):
            return x0 + x1 + t0 + t1
        g = lambda x: x
        h = lambda x: 2 * x
        args = 1, 2
        assert_quad(dblquad(func, 1, 2, g, h, args=args),35./6 + 9*.5)

    def test_triple_integral(self):
        # 9) Triple Integral test
        def simpfunc(z, y, x, t):      # Note order of arguments.