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]
python类log()的实例源码
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)))
