Python sympy 模块，log() 实例源码

def test_hyperexpand_bases():
    assert hyperexpand(hyper([2], [a], z)) == \
        a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)* \
        lowergamma(a - 1, z) - 1
    # TODO [a+1, a-S.Half], [2*a]
    assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
    assert hyperexpand(hyper([S.Half, 2], [S(3)/2], z)) == \
        -1/(2*z - 2) + atanh(sqrt(z))/sqrt(z)/2
    assert hyperexpand(hyper([S(1)/2, S(1)/2], [S(5)/2], z)) == \
        (-3*z + 3)/4/(z*sqrt(-z + 1)) \
        + (6*z - 3)*asin(sqrt(z))/(4*z**(S(3)/2))
    assert hyperexpand(hyper([1, 2], [S(3)/2], z)) == -1/(2*z - 2) \
        - asin(sqrt(z))/(sqrt(z)*(2*z - 2)*sqrt(-z + 1))
    assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == \
        sqrt(z)*(6*z/7 - S(6)/5)*atanh(sqrt(z)) \
        + (-30*z**2 + 32*z - 6)/35/z - 6*log(-z + 1)/(35*z**2)
    assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == \
        -4*log(sqrt(-z + 1)/2 + S(1)/2)/z
    # TODO hyperexpand(hyper([a], [2*a + 1], z))
    # TODO [S.Half, a], [S(3)/2, a+1]
    assert hyperexpand(hyper([2], [b, 1], z)) == \
        z**(-b/2 + S(1)/2)*besseli(b - 1, 2*sqrt(z))*gamma(b) \
        + z**(-b/2 + 1)*besseli(b, 2*sqrt(z))*gamma(b)
    # TODO [a], [a - S.Half, 2*a]

def _eval_power(self, other):
        from sympy.functions.elementary.exponential import log

        b, e = self.as_base_exp()
        b_nneg = b.is_nonnegative
        if b.is_real and not b_nneg and e.is_even:
            b = abs(b)
            b_nneg = True

        # Special case for when b is nan. See pull req 1714 for details
        if b is S.NaN:
            smallarg = abs(e).is_negative
        else:
            smallarg = (abs(e) - abs(S.Pi/log(b))).is_negative
        if (other.is_Rational and other.q == 2 and
                e.is_real is False and smallarg is False):
            return -self.func(b, e*other)
        if (other.is_integer or
            e.is_real and (b_nneg or (abs(e) < 1) == True) or
            e.is_real is False and smallarg is True or
                b.is_polar):
            return self.func(b, e*other)

def _eval_is_positive(self):
        if self.base.is_positive:
            if self.exp.is_real:
                return True
        elif self.base.is_negative:
            if self.exp.is_even:
                return True
            if self.exp.is_odd:
                return False
        elif self.base.is_nonpositive:
            if self.exp.is_odd:
                return False
        elif self.base.is_imaginary:
            if self.exp.is_integer:
                m = self.exp % 4
                if m.is_zero:
                    return True
                if m.is_integer and m.is_zero is False:
                    return False
            if self.exp.is_imaginary:
                return C.log(self.base).is_imaginary

def test__eis():
    assert _eis(z).diff(z) == -_eis(z) + 1/z

    assert _eis(1/z).series(z) == \
        z + z**2 + 2*z**3 + 6*z**4 + 24*z**5 + O(z**6)

    assert Ei(z).rewrite('tractable') == exp(z)*_eis(z)
    assert li(z).rewrite('tractable') == z*_eis(log(z))

    assert _eis(z).rewrite('intractable') == exp(-z)*Ei(z)

    assert expand(li(z).rewrite('tractable').diff(z).rewrite('intractable')) \
        == li(z).diff(z)

    assert expand(Ei(z).rewrite('tractable').diff(z).rewrite('intractable')) \
        == Ei(z).diff(z)

    assert _eis(z).series(z, n=3) == EulerGamma + log(z) + z*(-log(z) - \
        EulerGamma + 1) + z**2*(log(z)/2 - S(3)/4 + EulerGamma/2) + O(z**3*log(z))

def _eval_expand_func(self, **hints):
        from sympy import log, expand_mul, Dummy, exp_polar, I
        s, z = self.args
        if s == 1:
            return -log(1 + exp_polar(-I*pi)*z)
        if s.is_Integer and s <= 0:
            u = Dummy('u')
            start = u/(1 - u)
            for _ in range(-s):
                start = u*start.diff(u)
            return expand_mul(start).subs(u, z)
        return polylog(s, z)

###############################################################################
###################### HURWITZ GENERALIZED ZETA FUNCTION ######################
###############################################################################

def _getunbranched(cls, ar):
        from sympy import exp_polar, log, polar_lift
        if ar.is_Mul:
            args = ar.args
        else:
            args = [ar]
        unbranched = 0
        for a in args:
            if not a.is_polar:
                unbranched += arg(a)
            elif a.func is exp_polar:
                unbranched += a.exp.as_real_imag()[1]
            elif a.is_Pow:
                re, im = a.exp.as_real_imag()
                unbranched += re*unbranched_argument(
                    a.base) + im*log(abs(a.base))
            elif a.func is polar_lift:
                unbranched += arg(a.args[0])
            else:
                return None
        return unbranched

def test_maxima_functions():
    assert parse_maxima('expand( (x+1)^2)') == x**2 + 2*x + 1
    assert parse_maxima('factor( x**2 + 2*x + 1)') == (x + 1)**2
    assert parse_maxima('2*cos(x)^2 + sin(x)^2') == 2*cos(x)**2 + sin(x)**2
    assert parse_maxima('trigexpand(sin(2*x)+cos(2*x))') == \
        -1 + 2*cos(x)**2 + 2*cos(x)*sin(x)
    assert parse_maxima('solve(x^2-4,x)') == [-2, 2]
    assert parse_maxima('limit((1+1/x)^x,x,inf)') == E
    assert parse_maxima('limit(sqrt(-x)/x,x,0,minus)') == -oo
    assert parse_maxima('diff(x^x, x)') == x**x*(1 + log(x))
    assert parse_maxima('sum(k, k, 1, n)', name_dict=dict(
        n=Symbol('n', integer=True),
        k=Symbol('k', integer=True)
    )) == (n**2 + n)/2
    assert parse_maxima('product(k, k, 1, n)', name_dict=dict(
        n=Symbol('n', integer=True),
        k=Symbol('k', integer=True)
    )) == factorial(n)
    assert parse_maxima('ratsimp((x^2-1)/(x+1))') == x - 1
    assert Abs( parse_maxima(
        'float(sec(%pi/3) + csc(%pi/3))') - 3.154700538379252) < 10**(-5)

def Simple_manifold_with_scalar_function_derivative():
    Print_Function()
    coords = (x, y, z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3', metric='[1,1,1]', coords=coords)
    # Define surface
    mfvar = (u, v) = symbols('u v')
    X = u*e1 + v*e2 + (u**2 + v**2)*e3
    print('\\f{X}{u,v} =', X)
    MF = Manifold(X, mfvar)
    (eu, ev) = MF.Basis()
    # Define field on the surface.
    g = (v + 1)*log(u)

    print('\\f{g}{u,v} =', g)

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g, u) + (MF.rbasis[1]/MF.E_sq)*diff(g, v)
    print('\\eval{\\nabla g}{u=1,v=0} =', VectorDerivative.subs({u: 1, v: 0}))

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print('\\eval{\\f{Grad}{g}}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
    dg = MF.grad*g
    print('\\eval{\\nabla g}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
    return

def main():
    x = Symbol("x")
    a = Symbol("a")
    h = Symbol("h")

    show( limit(sqrt(x**2 - 5*x + 6) - x, x, oo), -Rational(5)/2 )

    show( limit(x*(sqrt(x**2 + 1) - x), x, oo), Rational(1)/2 )

    show( limit(x - sqrt3(x**3 - 1), x, oo), Rational(0) )

    show( limit(log(1 + exp(x))/x, x, -oo), Rational(0) )

    show( limit(log(1 + exp(x))/x, x, oo), Rational(1) )

    show( limit(sin(3*x)/x, x, 0), Rational(3) )

    show( limit(sin(5*x)/sin(2*x), x, 0), Rational(5)/2 )

    show( limit(((x - 1)/(x + 1))**x, x, oo), exp(-2))

def main():
    x = sympy.Symbol('x')
    y = sympy.Symbol('y')

    e = 1/sympy.cos(x)
    print()
    pprint(e)
    print('\n')
    pprint(e.subs(sympy.cos(x), y))
    print('\n')
    pprint(e.subs(sympy.cos(x), y).subs(y, x**2))

    e = 1/sympy.log(x)
    e = e.subs(x, sympy.Float("2.71828"))
    print('\n')
    pprint(e)
    print('\n')
    pprint(e.evalf())
    print()

def real(self, nested_scope=None):
        """Return the correspond floating point number."""
        op = self.children[0].name
        expr = self.children[1]
        dispatch = {
            'sin': sympy.sin,
            'cos': sympy.cos,
            'tan': sympy.tan,
            'asin': sympy.asin,
            'acos': sympy.acos,
            'atan': sympy.atan,
            'exp': sympy.exp,
            'ln': sympy.log,
            'sqrt': sympy.sqrt
        }
        if op in dispatch:
            arg = expr.real(nested_scope)
            return dispatch[op](arg)
        else:
            raise NodeException("internal error: undefined external")

def sym(self, nested_scope=None):
        """Return the corresponding symbolic expression."""
        op = self.children[0].name
        expr = self.children[1]
        dispatch = {
            'sin': sympy.sin,
            'cos': sympy.cos,
            'tan': sympy.tan,
            'asin': sympy.asin,
            'acos': sympy.acos,
            'atan': sympy.atan,
            'exp': sympy.exp,
            'ln': sympy.log,
            'sqrt': sympy.sqrt
        }
        if op in dispatch:
            arg = expr.sym(nested_scope)
            return dispatch[op](arg)
        else:
            raise NodeException("internal error: undefined external")

def get_S95(b0, sigma):

    S95 = sy.Symbol('S95', positive = True, real = True)
    b = sy.Symbol('b', positive = True)
    chi21 = sy.Symbol('chi21')
    chi22 = sy.Symbol('chi22')

    chi2 = 3.84
    N = b0 

    replacements = [(b, (b0 - S95 - sigma**2)/2 + 1./2*((b0 - S95 - sigma**2)**2 + 4*(sigma**2*N - S95*sigma**2 + S95*b0))**0.5)]

    replacements2 = [(S95, 0.)]

    chi21 = -2*( N*sy.log(S95 + b) - (S95 + b) - ((b-b0)/sigma)**2)

    chi21 = chi21.subs(replacements)
    chi22 = chi21.subs(replacements2)

    eq = chi2 - chi21 + chi22

    S95_new = sy.nsolve(eq, S95, 1)

    return float(S95_new)

def test_hyperexpand_bases():
    assert hyperexpand(hyper([2], [a], z)) == \
        a + z**(-a + 1)*(-a**2 + 3*a + z*(a - 1) - 2)*exp(z)* \
        lowergamma(a - 1, z) - 1
    # TODO [a+1, a-S.Half], [2*a]
    assert hyperexpand(hyper([1, 2], [3], z)) == -2/z - 2*log(-z + 1)/z**2
    assert hyperexpand(hyper([S.Half, 2], [S(3)/2], z)) == \
        -1/(2*z - 2) + atanh(sqrt(z))/sqrt(z)/2
    assert hyperexpand(hyper([S(1)/2, S(1)/2], [S(5)/2], z)) == \
        (-3*z + 3)/4/(z*sqrt(-z + 1)) \
        + (6*z - 3)*asin(sqrt(z))/(4*z**(S(3)/2))
    assert hyperexpand(hyper([1, 2], [S(3)/2], z)) == -1/(2*z - 2) \
        - asin(sqrt(z))/(sqrt(z)*(2*z - 2)*sqrt(-z + 1))
    assert hyperexpand(hyper([-S.Half - 1, 1, 2], [S.Half, 3], z)) == \
        sqrt(z)*(6*z/7 - S(6)/5)*atanh(sqrt(z)) \
        + (-30*z**2 + 32*z - 6)/35/z - 6*log(-z + 1)/(35*z**2)
    assert hyperexpand(hyper([1 + S.Half, 1, 1], [2, 2], z)) == \
        -4*log(sqrt(-z + 1)/2 + S(1)/2)/z
    # TODO hyperexpand(hyper([a], [2*a + 1], z))
    # TODO [S.Half, a], [S(3)/2, a+1]
    assert hyperexpand(hyper([2], [b, 1], z)) == \
        z**(-b/2 + S(1)/2)*besseli(b - 1, 2*sqrt(z))*gamma(b) \
        + z**(-b/2 + 1)*besseli(b, 2*sqrt(z))*gamma(b)
    # TODO [a], [a - S.Half, 2*a]

def _eval_is_positive(self):
        from sympy import log
        if self.base == self.exp:
            if self.base.is_nonnegative:
                return True
        elif self.base.is_positive:
            if self.exp.is_real:
                return True
        elif self.base.is_negative:
            if self.exp.is_even:
                return True
            if self.exp.is_odd:
                return False
        elif self.base.is_nonpositive:
            if self.exp.is_odd:
                return False
        elif self.base.is_imaginary:
            if self.exp.is_integer:
                m = self.exp % 4
                if m.is_zero:
                    return True
                if m.is_integer and m.is_zero is False:
                    return False
            if self.exp.is_imaginary:
                return log(self.base).is_imaginary

def is_number(self):
        """Returns True if 'self' has no free symbols.
        It will be faster than `if not self.free_symbols`, however, since
        `is_number` will fail as soon as it hits a free symbol.

        Examples
        ========

        >>> from sympy import log, Integral
        >>> from sympy.abc import x

        >>> x.is_number
        False
        >>> (2*x).is_number
        False
        >>> (2 + log(2)).is_number
        True
        >>> (2 + Integral(2, x)).is_number
        False
        >>> (2 + Integral(2, (x, 1, 2))).is_number
        True

        """
        return all(obj.is_number for obj in self.args)

def compute_leading_term(self, x, logx=None):
        """
        as_leading_term is only allowed for results of .series()
        This is a wrapper to compute a series first.
        """
        from sympy import Dummy, log
        from sympy.series.gruntz import calculate_series

        if self.removeO() == 0:
            return self

        if logx is None:
            d = Dummy('logx')
            s = calculate_series(self, x, d).subs(d, log(x))
        else:
            s = calculate_series(self, x, logx)

        return s.as_leading_term(x)

def test_derivative_by_array():
    from sympy.abc import a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z

    bexpr = x*y**2*exp(z)*log(t)
    sexpr = sin(bexpr)
    cexpr = cos(bexpr)

    a = Array([sexpr])

    assert derive_by_array(sexpr, t) == x*y**2*exp(z)*cos(x*y**2*exp(z)*log(t))/t
    assert derive_by_array(sexpr, [x, y, z]) == Array([bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr, bexpr*cexpr])
    assert derive_by_array(a, [x, y, z]) == Array([[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr], [bexpr*cexpr]])

    assert derive_by_array(sexpr, [[x, y], [z, t]]) == Array([[bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr], [bexpr*cexpr, bexpr/log(t)/t*cexpr]])
    assert derive_by_array(a, [[x, y], [z, t]]) == Array([[[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr]], [[bexpr*cexpr], [bexpr/log(t)/t*cexpr]]])
    assert derive_by_array([[x, y], [z, t]], [x, y]) == Array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]])
    assert derive_by_array([[x, y], [z, t]], [[x, y], [z, t]]) == Array([[[[1, 0], [0, 0]], [[0, 1], [0, 0]]],
                                                                         [[[0, 0], [1, 0]], [[0, 0], [0, 1]]]])

def test__eis():
    assert _eis(z).diff(z) == -_eis(z) + 1/z

    assert _eis(1/z).series(z) == \
        z + z**2 + 2*z**3 + 6*z**4 + 24*z**5 + O(z**6)

    assert Ei(z).rewrite('tractable') == exp(z)*_eis(z)
    assert li(z).rewrite('tractable') == z*_eis(log(z))

    assert _eis(z).rewrite('intractable') == exp(-z)*Ei(z)

    assert expand(li(z).rewrite('tractable').diff(z).rewrite('intractable')) \
        == li(z).diff(z)

    assert expand(Ei(z).rewrite('tractable').diff(z).rewrite('intractable')) \
        == Ei(z).diff(z)

    assert _eis(z).series(z, n=3) == EulerGamma + log(z) + z*(-log(z) - \
        EulerGamma + 1) + z**2*(log(z)/2 - S(3)/4 + EulerGamma/2) + O(z**3*log(z))

def _getunbranched(cls, ar):
        from sympy import exp_polar, log, polar_lift
        if ar.is_Mul:
            args = ar.args
        else:
            args = [ar]
        unbranched = 0
        for a in args:
            if not a.is_polar:
                unbranched += arg(a)
            elif a.func is exp_polar:
                unbranched += a.exp.as_real_imag()[1]
            elif a.is_Pow:
                re, im = a.exp.as_real_imag()
                unbranched += re*unbranched_argument(
                    a.base) + im*log(abs(a.base))
            elif a.func is polar_lift:
                unbranched += arg(a.args[0])
            else:
                return None
        return unbranched

def power_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if symbol not in exp.free_symbols and isinstance(base, sympy.Symbol):
        if sympy.simplify(exp + 1) == 0:
            return ReciprocalRule(base, integrand, symbol)
        return PowerRule(base, exp, integrand, symbol)
    elif symbol not in base.free_symbols and isinstance(exp, sympy.Symbol):
        rule = ExpRule(base, exp, integrand, symbol)

        if fuzzy_not(sympy.log(base).is_zero):
            return rule
        elif sympy.log(base).is_zero:
            return ConstantRule(1, 1, symbol)

        return PiecewiseRule([
            (ConstantRule(1, 1, symbol), sympy.Eq(sympy.log(base), 0)),
            (rule, True)
        ], integrand, symbol)

def test_maxima_functions():
    assert parse_maxima('expand( (x+1)^2)') == x**2 + 2*x + 1
    assert parse_maxima('factor( x**2 + 2*x + 1)') == (x + 1)**2
    assert parse_maxima('2*cos(x)^2 + sin(x)^2') == 2*cos(x)**2 + sin(x)**2
    assert parse_maxima('trigexpand(sin(2*x)+cos(2*x))') == \
        -1 + 2*cos(x)**2 + 2*cos(x)*sin(x)
    assert parse_maxima('solve(x^2-4,x)') == [-2, 2]
    assert parse_maxima('limit((1+1/x)^x,x,inf)') == E
    assert parse_maxima('limit(sqrt(-x)/x,x,0,minus)') == -oo
    assert parse_maxima('diff(x^x, x)') == x**x*(1 + log(x))
    assert parse_maxima('sum(k, k, 1, n)', name_dict=dict(
        n=Symbol('n', integer=True),
        k=Symbol('k', integer=True)
    )) == (n**2 + n)/2
    assert parse_maxima('product(k, k, 1, n)', name_dict=dict(
        n=Symbol('n', integer=True),
        k=Symbol('k', integer=True)
    )) == factorial(n)
    assert parse_maxima('ratsimp((x^2-1)/(x+1))') == x - 1
    assert Abs( parse_maxima(
        'float(sec(%pi/3) + csc(%pi/3))') - 3.154700538379252) < 10**(-5)

def fdebug(self):

        theta = self.N_theta_right
        phi = self.N_phi
        dma = self.frontend.data_ma

        _theta, _phi = self._reduce_latent()
        print(_theta.sum())
        print(_phi.sum())

        p_ij = _theta.dot(_phi).dot(_theta.T)
        pij = self.data_A * p_ij + self.data_B

        ll = - np.log(pij).sum() / self._len['nnz']

        print(dma.sum())
        print(p_ij.sum())
        print(pij.sum())
        print(ll)

def recursive_line(self, new_line=5246):
        stir = self.load()
        J = stir.shape[0]
        K = stir.shape[1]
        for x in range(new_line):
            n = J + x
            new_l =  np.ones((1, K)) * np.inf
            print(n)
            for m in range(1,K):
                if m > n:
                    continue
                elif m == n:
                    new_l[0, m] = 0
                elif m == 1:
                    new_l[0, 1] = logsumexp( [  np.log(n-1) + stir[n-1, m] ] )
                else:
                    new_l[0, m] = logsumexp( [ stir[n-1, m-1] , np.log(n-1) + stir[n-1, m] ] )
            stir = np.vstack((stir, new_l))

        #np.save('stirling.npy', stir)
        #np.load('stirling.npy')
        return stir

def recursive_row(self, new_row=''):
        stir = np.load('stirling.npy')
        J = stir.shape[0]
        K = stir.shape[1]
        x = 0
        while stir.shape[0] != stir.shape[1]:
            m = K + x
            new_c =  np.ones((J, 1)) * np.inf
            stir = np.hstack((stir, new_c))
            print(m)
            for n in range(K,J):
                if m > n:
                    continue
                elif m == n:
                    stir[n, m] = 0

                else:
                    stir[n,m] = logsumexp( [ stir[n-1, m-1] , np.log(n-1) + stir[n-1, m] ] )
            x += 1

        #np.save('stirling.npy', stir)
        #np.load('stirling.npy',)
        return stir

def NFW_mass_profile(r="r", rho_s="rho_s", r_s="r_s"):
    """
    An NFW mass profile (Navarro, J.F., Frenk, C.S.,
    & White, S.D.M. 1996, ApJ, 462, 563).

    Parameters
    ----------
    rho_s : string
        The symbol for the scale density of the profile.
    r_s : string
        The symbol for the scale radius.
    """
    r, rho_s, r_s = symbols((r, rho_s, r_s))
    x = r/r_s
    profile = 4*pi*rho_s*r_s**3*(log(1+x)-x/(1+x))
    return Formula1D(profile, r, [rho_s, r_s])

def solved_vc_inequality(probability, error, approximate_datapoints_needed):
    n = Symbol('n')
    return nsolve(log(4) + 10 * log(2*n) - 1/8 * error ** 2 * n - log(probability), n, approximate_datapoints_needed)

def original_vc_bound(n, delta, growth_function):
    return sqrt(8/n * log(4*growth_function(2*n)/delta))

def rademacher_penalty_bound(n, delta, growth_function):
    return sqrt(2 * log(2 * n * growth_function(n)) / n) + sqrt(2/n * log(1/delta)) + 1/n

def parrondo_van_den_broek_right(error, n, delta, growth_function):
    return sqrt(1/n * (2 * error + log(6 * growth_function(2*n)/delta)))

def devroye(error, n, delta, growth_function):
    return sqrt(1/(2*Decimal(n)) * (4 * error * (1 + error) + log(4 * growth_function(Decimal(n**2))/Decimal(delta))))

def _fit_parameters(feature_matrix, data_quantities, feature_tuple):
    """Solve Ax = b, where 'feature_matrix' is A and 'data_quantities' is b.

    Parameters
    ----------
    feature_matrix : ndarray
        (M*N) regressor matrix.
    data_quantities : ndarray
        (M,) response vector
    feature_tuple : (float
        Polynomial coefficient corresponding to each column of 'feature_matrix'

    Returns
    -------
    OrderedDict
        {featured_tuple: fitted_parameter}. Maps 'feature_tuple'
        to fitted parameter value. If a coefficient is not used, it maps to zero.

    """
    # Now generate candidate models; add parameters one at a time
    model_scores = []
    results = np.zeros((len(feature_tuple), len(feature_tuple)))
    clf = LinearRegression(fit_intercept=False, normalize=True)
    for num_params in range(1, feature_matrix.shape[-1] + 1):
        current_matrix = feature_matrix[:, :num_params]
        clf.fit(current_matrix, data_quantities)
        # This may not exactly be the correct form for the likelihood
        # We're missing the "ridge" contribution here which could become relevant for sparse data
        rss = np.square(np.dot(current_matrix, clf.coef_) - data_quantities.astype(np.float)).sum()
        # Compute Aikaike Information Criterion
        # Form valid under assumption all sample variances are equal and unknown
        score = 2*num_params + current_matrix.shape[-2] * np.log(rss)
        model_scores.append(score)
        results[num_params - 1, :num_params] = clf.coef_
        logging.debug('{} rss: {}, AIC: {}'.format(feature_tuple[:num_params], rss, score))
    return OrderedDict(zip(feature_tuple, results[np.argmin(model_scores), :]))

def test_bug():
    assert residue(2**(z)*(s + z)*(1 - s - z)/z**2, z, 0) == \
        1 + s*log(2) - s**2*log(2) - 2*s

def test_hyperexpand():
    # Luke, Y. L. (1969), The Special Functions and Their Approximations,
    # Volume 1, section 6.2

    assert hyperexpand(hyper([], [], z)) == exp(z)
    assert hyperexpand(hyper([1, 1], [2], -z)*z) == log(1 + z)
    assert hyperexpand(hyper([], [S.Half], -z**2/4)) == cos(z)
    assert hyperexpand(z*hyper([], [S('3/2')], -z**2/4)) == sin(z)
    assert hyperexpand(hyper([S('1/2'), S('1/2')], [S('3/2')], z**2)*z) \
        == asin(z)

def can_do_meijer(a1, a2, b1, b2, numeric=True):
    """
    This helper function tries to hyperexpand() the meijer g-function
    corresponding to the parameters a1, a2, b1, b2.
    It returns False if this expansion still contains g-functions.
    If numeric is True, it also tests the so-obtained formula numerically
    (at random values) and returns False if the test fails.
    Else it returns True.
    """
    from sympy import unpolarify, expand
    r = hyperexpand(meijerg(a1, a2, b1, b2, z))
    if r.has(meijerg):
        return False
    # NOTE hyperexpand() returns a truly branched function, whereas numerical
    #      evaluation only works on the main branch. Since we are evaluating on
    #      the main branch, this should not be a problem, but expressions like
    #      exp_polar(I*pi/2*x)**a are evaluated incorrectly. We thus have to get
    #      rid of them. The expand heuristically does this...
    r = unpolarify(expand(r, force=True, power_base=True, power_exp=False,
                          mul=False, log=False, multinomial=False, basic=False))

    if not numeric:
        return True

    repl = {}
    for n, a in enumerate(meijerg(a1, a2, b1, b2, z).free_symbols - set([z])):
        repl[a] = randcplx(n)
    return tn(meijerg(a1, a2, b1, b2, z).subs(repl), r.subs(repl), z)

def __new__(cls, b, e, evaluate=None):
        if evaluate is None:
            evaluate = global_evaluate[0]
        from sympy.functions.elementary.exponential import exp_polar

        # don't optimize "if e==0; return 1" here; it's better to handle that
        # in the calling routine so this doesn't get called
        b = _sympify(b)
        e = _sympify(e)
        if evaluate:
            if e is S.Zero:
                return S.One
            elif e is S.One:
                return b
            elif b is S.One:
                if e in (S.NaN, S.Infinity, -S.Infinity):
                    return S.NaN
                return S.One
            elif S.NaN in (b, e):
                return S.NaN
            else:
                # recognize base as E
                if not e.is_Atom and b is not S.Exp1 and b.func is not exp_polar:
                    from sympy import numer, denom, log, sign, im, factor_terms
                    c, ex = factor_terms(e, sign=False).as_coeff_Mul()
                    den = denom(ex)
                    if den.func is log and den.args[0] == b:
                        return S.Exp1**(c*numer(ex))
                    elif den.is_Add:
                        s = sign(im(b))
                        if s.is_Number and s and den == \
                                log(-factor_terms(b, sign=False)) + s*S.ImaginaryUnit*S.Pi:
                            return S.Exp1**(c*numer(ex))

                obj = b._eval_power(e)
                if obj is not None:
                    return obj
        obj = Expr.__new__(cls, b, e)
        obj.is_commutative = (b.is_commutative and e.is_commutative)
        return obj

def _eval_derivative(self, s):
        dbase = self.base.diff(s)
        dexp = self.exp.diff(s)
        return self * (dexp * C.log(self.base) + dbase * self.exp/self.base)

def _eval_as_leading_term(self, x):
        if not self.exp.has(x):
            return self.func(self.base.as_leading_term(x), self.exp)
        return C.exp(self.exp * C.log(self.base)).as_leading_term(x)

def test_special_is_rational():
    i = Symbol('i', integer=True)
    r = Symbol('r', rational=True)
    x = Symbol('x')
    assert sqrt(3).is_rational is False
    assert (3 + sqrt(3)).is_rational is False
    assert (3*sqrt(3)).is_rational is False
    assert exp(3).is_rational is False
    assert exp(i).is_rational is False
    assert exp(r).is_rational is False
    assert exp(x).is_rational is None
    assert exp(log(3), evaluate=False).is_rational is True
    assert log(exp(3), evaluate=False).is_rational is True
    assert log(3).is_rational is False
    assert log(i).is_rational is False
    assert log(r).is_rational is False
    assert log(x).is_rational is None
    assert (sqrt(3) + sqrt(5)).is_rational is None
    assert (sqrt(3) + S.Pi).is_rational is None
    assert (x**i).is_rational is None
    assert (i**i).is_rational is True
    assert (r**i).is_rational is True
    assert (r**r).is_rational is None
    assert (r**x).is_rational is None
    assert sin(1).is_rational is False
    assert sin(i).is_rational is False
    assert sin(r).is_rational is False
    assert sin(x).is_rational is None
    assert asin(r).is_rational is False
    assert sin(asin(3), evaluate=False).is_rational is True

def test_singularities_non_rational():
    x = Symbol('x', real=True)

    assert singularities(exp(1/x), x) == (0)
    assert singularities(log((x - 2)**2), x) == (2)

def test_gosper_sum():
    assert gosper_sum(1, (k, 0, n)) == 1 + n
    assert gosper_sum(k, (k, 0, n)) == n*(1 + n)/2
    assert gosper_sum(k**2, (k, 0, n)) == n*(1 + n)*(1 + 2*n)/6
    assert gosper_sum(k**3, (k, 0, n)) == n**2*(1 + n)**2/4

    assert gosper_sum(2**k, (k, 0, n)) == 2*2**n - 1

    assert gosper_sum(factorial(k), (k, 0, n)) is None
    assert gosper_sum(binomial(n, k), (k, 0, n)) is None

    assert gosper_sum(factorial(k)/k**2, (k, 0, n)) is None
    assert gosper_sum((k - 3)*factorial(k), (k, 0, n)) is None

    assert gosper_sum(k*factorial(k), k) == factorial(k)
    assert gosper_sum(
        k*factorial(k), (k, 0, n)) == n*factorial(n) + factorial(n) - 1

    assert gosper_sum((-1)**k*binomial(n, k), (k, 0, n)) == 0
    assert gosper_sum((
        -1)**k*binomial(n, k), (k, 0, m)) == -(-1)**m*(m - n)*binomial(n, m)/n

    assert gosper_sum((4*k + 1)*factorial(k)/factorial(2*k + 1), (k, 0, n)) == \
        (2*factorial(2*n + 1) - factorial(n))/factorial(2*n + 1)

    # issue 6033:
    assert gosper_sum(
        n*(n + a + b)*a**n*b**n/(factorial(n + a)*factorial(n + b)), \
        (n, 0, m)) == -a*b*(exp(m*log(a))*exp(m*log(b))*factorial(a)* \
        factorial(b) - factorial(a + m)*factorial(b + m))/(factorial(a)* \
        factorial(b)*factorial(a + m)*factorial(b + m))

def test_entropy():
    up = JzKet(S(1)/2, S(1)/2)
    down = JzKet(S(1)/2, -S(1)/2)
    d = Density((up, 0.5), (down, 0.5))

    # test for density object
    ent = entropy(d)
    assert entropy(d) == 0.5*log(2)
    assert d.entropy() == 0.5*log(2)

    np = import_module('numpy', min_module_version='1.4.0')
    if np:
        #do this test only if 'numpy' is available on test machine
        np_mat = represent(d, format='numpy')
        ent = entropy(np_mat)
        assert isinstance(np_mat, np.matrixlib.defmatrix.matrix)
        assert ent.real == 0.69314718055994529
        assert ent.imag == 0

    scipy = import_module('scipy', __import__kwargs={'fromlist': ['sparse']})
    if scipy and np:
        #do this test only if numpy and scipy are available
        mat = represent(d, format="scipy.sparse")
        assert isinstance(mat, scipy_sparse_matrix)
        assert ent.real == 0.69314718055994529
        assert ent.imag == 0

def test_polylog_expansion():
    from sympy import factor, log
    assert polylog(s, 0) == 0
    assert polylog(s, 1) == zeta(s)
    assert polylog(s, -1) == dirichlet_eta(s)

    assert myexpand(polylog(1, z), -log(1 + exp_polar(-I*pi)*z))
    assert myexpand(polylog(0, z), z/(1 - z))
    assert myexpand(polylog(-1, z), z**2/(1 - z)**2 + z/(1 - z))
    assert myexpand(polylog(-5, z), None)

def test_ei():
    pos = Symbol('p', positive=True)
    neg = Symbol('n', negative=True)
    assert Ei(-pos) == Ei(polar_lift(-1)*pos) - I*pi
    assert Ei(neg) == Ei(polar_lift(neg)) - I*pi
    assert tn_branch(Ei)
    assert mytd(Ei(x), exp(x)/x, x)
    assert mytn(Ei(x), Ei(x).rewrite(uppergamma),
                -uppergamma(0, x*polar_lift(-1)) - I*pi, x)
    assert mytn(Ei(x), Ei(x).rewrite(expint),
                -expint(1, x*polar_lift(-1)) - I*pi, x)
    assert Ei(x).rewrite(expint).rewrite(Ei) == Ei(x)
    assert Ei(x*exp_polar(2*I*pi)) == Ei(x) + 2*I*pi
    assert Ei(x*exp_polar(-2*I*pi)) == Ei(x) - 2*I*pi

    assert mytn(Ei(x), Ei(x).rewrite(Shi), Chi(x) + Shi(x), x)
    assert mytn(Ei(x*polar_lift(I)), Ei(x*polar_lift(I)).rewrite(Si),
                Ci(x) + I*Si(x) + I*pi/2, x)

    assert Ei(log(x)).rewrite(li) == li(x)
    assert Ei(2*log(x)).rewrite(li) == li(x**2)

    assert gruntz(Ei(x+exp(-x))*exp(-x)*x, x, oo) == 1

    assert Ei(x).series(x) == EulerGamma + log(x) + x + x**2/4 + \
        x**3/18 + x**4/96 + x**5/600 + O(x**6)