bench_meijerint.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261
  1. # conceal the implicit import from the code quality tester
  2. from sympy.core.numbers import (oo, pi)
  3. from sympy.core.symbol import (Symbol, symbols)
  4. from sympy.functions.elementary.exponential import exp
  5. from sympy.functions.elementary.miscellaneous import sqrt
  6. from sympy.functions.special.bessel import besseli
  7. from sympy.functions.special.gamma_functions import gamma
  8. from sympy.integrals.integrals import integrate
  9. from sympy.integrals.transforms import (mellin_transform,
  10. inverse_fourier_transform, inverse_mellin_transform,
  11. laplace_transform, inverse_laplace_transform, fourier_transform)
  12. LT = laplace_transform
  13. FT = fourier_transform
  14. MT = mellin_transform
  15. IFT = inverse_fourier_transform
  16. ILT = inverse_laplace_transform
  17. IMT = inverse_mellin_transform
  18. from sympy.abc import x, y
  19. nu, beta, rho = symbols('nu beta rho')
  20. apos, bpos, cpos, dpos, posk, p = symbols('a b c d k p', positive=True)
  21. k = Symbol('k', real=True)
  22. negk = Symbol('k', negative=True)
  23. mu1, mu2 = symbols('mu1 mu2', real=True, nonzero=True, finite=True)
  24. sigma1, sigma2 = symbols('sigma1 sigma2', real=True, nonzero=True,
  25. finite=True, positive=True)
  26. rate = Symbol('lambda', positive=True)
  27. def normal(x, mu, sigma):
  28. return 1/sqrt(2*pi*sigma**2)*exp(-(x - mu)**2/2/sigma**2)
  29. def exponential(x, rate):
  30. return rate*exp(-rate*x)
  31. alpha, beta = symbols('alpha beta', positive=True)
  32. betadist = x**(alpha - 1)*(1 + x)**(-alpha - beta)*gamma(alpha + beta) \
  33. /gamma(alpha)/gamma(beta)
  34. kint = Symbol('k', integer=True, positive=True)
  35. chi = 2**(1 - kint/2)*x**(kint - 1)*exp(-x**2/2)/gamma(kint/2)
  36. chisquared = 2**(-k/2)/gamma(k/2)*x**(k/2 - 1)*exp(-x/2)
  37. dagum = apos*p/x*(x/bpos)**(apos*p)/(1 + x**apos/bpos**apos)**(p + 1)
  38. d1, d2 = symbols('d1 d2', positive=True)
  39. f = sqrt(((d1*x)**d1 * d2**d2)/(d1*x + d2)**(d1 + d2))/x \
  40. /gamma(d1/2)/gamma(d2/2)*gamma((d1 + d2)/2)
  41. nupos, sigmapos = symbols('nu sigma', positive=True)
  42. rice = x/sigmapos**2*exp(-(x**2 + nupos**2)/2/sigmapos**2)*besseli(0, x*
  43. nupos/sigmapos**2)
  44. mu = Symbol('mu', real=True)
  45. laplace = exp(-abs(x - mu)/bpos)/2/bpos
  46. u = Symbol('u', polar=True)
  47. tpos = Symbol('t', positive=True)
  48. def E(expr):
  49. integrate(expr*exponential(x, rate)*normal(y, mu1, sigma1),
  50. (x, 0, oo), (y, -oo, oo), meijerg=True)
  51. integrate(expr*exponential(x, rate)*normal(y, mu1, sigma1),
  52. (y, -oo, oo), (x, 0, oo), meijerg=True)
  53. bench = [
  54. 'MT(x**nu*Heaviside(x - 1), x, s)',
  55. 'MT(x**nu*Heaviside(1 - x), x, s)',
  56. 'MT((1-x)**(beta - 1)*Heaviside(1-x), x, s)',
  57. 'MT((x-1)**(beta - 1)*Heaviside(x-1), x, s)',
  58. 'MT((1+x)**(-rho), x, s)',
  59. 'MT(abs(1-x)**(-rho), x, s)',
  60. 'MT((1-x)**(beta-1)*Heaviside(1-x) + a*(x-1)**(beta-1)*Heaviside(x-1), x, s)',
  61. 'MT((x**a-b**a)/(x-b), x, s)',
  62. 'MT((x**a-bpos**a)/(x-bpos), x, s)',
  63. 'MT(exp(-x), x, s)',
  64. 'MT(exp(-1/x), x, s)',
  65. 'MT(log(x)**4*Heaviside(1-x), x, s)',
  66. 'MT(log(x)**3*Heaviside(x-1), x, s)',
  67. 'MT(log(x + 1), x, s)',
  68. 'MT(log(1/x + 1), x, s)',
  69. 'MT(log(abs(1 - x)), x, s)',
  70. 'MT(log(abs(1 - 1/x)), x, s)',
  71. 'MT(log(x)/(x+1), x, s)',
  72. 'MT(log(x)**2/(x+1), x, s)',
  73. 'MT(log(x)/(x+1)**2, x, s)',
  74. 'MT(erf(sqrt(x)), x, s)',
  75. 'MT(besselj(a, 2*sqrt(x)), x, s)',
  76. 'MT(sin(sqrt(x))*besselj(a, sqrt(x)), x, s)',
  77. 'MT(cos(sqrt(x))*besselj(a, sqrt(x)), x, s)',
  78. 'MT(besselj(a, sqrt(x))**2, x, s)',
  79. 'MT(besselj(a, sqrt(x))*besselj(-a, sqrt(x)), x, s)',
  80. 'MT(besselj(a - 1, sqrt(x))*besselj(a, sqrt(x)), x, s)',
  81. 'MT(besselj(a, sqrt(x))*besselj(b, sqrt(x)), x, s)',
  82. 'MT(besselj(a, sqrt(x))**2 + besselj(-a, sqrt(x))**2, x, s)',
  83. 'MT(bessely(a, 2*sqrt(x)), x, s)',
  84. 'MT(sin(sqrt(x))*bessely(a, sqrt(x)), x, s)',
  85. 'MT(cos(sqrt(x))*bessely(a, sqrt(x)), x, s)',
  86. 'MT(besselj(a, sqrt(x))*bessely(a, sqrt(x)), x, s)',
  87. 'MT(besselj(a, sqrt(x))*bessely(b, sqrt(x)), x, s)',
  88. 'MT(bessely(a, sqrt(x))**2, x, s)',
  89. 'MT(besselk(a, 2*sqrt(x)), x, s)',
  90. 'MT(besselj(a, 2*sqrt(2*sqrt(x)))*besselk(a, 2*sqrt(2*sqrt(x))), x, s)',
  91. 'MT(besseli(a, sqrt(x))*besselk(a, sqrt(x)), x, s)',
  92. 'MT(besseli(b, sqrt(x))*besselk(a, sqrt(x)), x, s)',
  93. 'MT(exp(-x/2)*besselk(a, x/2), x, s)',
  94. # later: ILT, IMT
  95. 'LT((t-apos)**bpos*exp(-cpos*(t-apos))*Heaviside(t-apos), t, s)',
  96. 'LT(t**apos, t, s)',
  97. 'LT(Heaviside(t), t, s)',
  98. 'LT(Heaviside(t - apos), t, s)',
  99. 'LT(1 - exp(-apos*t), t, s)',
  100. 'LT((exp(2*t)-1)*exp(-bpos - t)*Heaviside(t)/2, t, s, noconds=True)',
  101. 'LT(exp(t), t, s)',
  102. 'LT(exp(2*t), t, s)',
  103. 'LT(exp(apos*t), t, s)',
  104. 'LT(log(t/apos), t, s)',
  105. 'LT(erf(t), t, s)',
  106. 'LT(sin(apos*t), t, s)',
  107. 'LT(cos(apos*t), t, s)',
  108. 'LT(exp(-apos*t)*sin(bpos*t), t, s)',
  109. 'LT(exp(-apos*t)*cos(bpos*t), t, s)',
  110. 'LT(besselj(0, t), t, s, noconds=True)',
  111. 'LT(besselj(1, t), t, s, noconds=True)',
  112. 'FT(Heaviside(1 - abs(2*apos*x)), x, k)',
  113. 'FT(Heaviside(1-abs(apos*x))*(1-abs(apos*x)), x, k)',
  114. 'FT(exp(-apos*x)*Heaviside(x), x, k)',
  115. 'IFT(1/(apos + 2*pi*I*x), x, posk, noconds=False)',
  116. 'IFT(1/(apos + 2*pi*I*x), x, -posk, noconds=False)',
  117. 'IFT(1/(apos + 2*pi*I*x), x, negk)',
  118. 'FT(x*exp(-apos*x)*Heaviside(x), x, k)',
  119. 'FT(exp(-apos*x)*sin(bpos*x)*Heaviside(x), x, k)',
  120. 'FT(exp(-apos*x**2), x, k)',
  121. 'IFT(sqrt(pi/apos)*exp(-(pi*k)**2/apos), k, x)',
  122. 'FT(exp(-apos*abs(x)), x, k)',
  123. 'integrate(normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
  124. 'integrate(x*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
  125. 'integrate(x**2*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
  126. 'integrate(x**3*normal(x, mu1, sigma1), (x, -oo, oo), meijerg=True)',
  127. 'integrate(normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  128. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  129. 'integrate(x*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  130. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  131. 'integrate(y*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  132. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  133. 'integrate(x*y*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  134. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  135. 'integrate((x+y+1)*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  136. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  137. 'integrate((x+y-1)*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  138. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  139. 'integrate(x**2*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  140. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  141. 'integrate(y**2*normal(x, mu1, sigma1)*normal(y, mu2, sigma2),'
  142. ' (x, -oo, oo), (y, -oo, oo), meijerg=True)',
  143. 'integrate(exponential(x, rate), (x, 0, oo), meijerg=True)',
  144. 'integrate(x*exponential(x, rate), (x, 0, oo), meijerg=True)',
  145. 'integrate(x**2*exponential(x, rate), (x, 0, oo), meijerg=True)',
  146. 'E(1)',
  147. 'E(x*y)',
  148. 'E(x*y**2)',
  149. 'E((x+y+1)**2)',
  150. 'E(x+y+1)',
  151. 'E((x+y-1)**2)',
  152. 'integrate(betadist, (x, 0, oo), meijerg=True)',
  153. 'integrate(x*betadist, (x, 0, oo), meijerg=True)',
  154. 'integrate(x**2*betadist, (x, 0, oo), meijerg=True)',
  155. 'integrate(chi, (x, 0, oo), meijerg=True)',
  156. 'integrate(x*chi, (x, 0, oo), meijerg=True)',
  157. 'integrate(x**2*chi, (x, 0, oo), meijerg=True)',
  158. 'integrate(chisquared, (x, 0, oo), meijerg=True)',
  159. 'integrate(x*chisquared, (x, 0, oo), meijerg=True)',
  160. 'integrate(x**2*chisquared, (x, 0, oo), meijerg=True)',
  161. 'integrate(((x-k)/sqrt(2*k))**3*chisquared, (x, 0, oo), meijerg=True)',
  162. 'integrate(dagum, (x, 0, oo), meijerg=True)',
  163. 'integrate(x*dagum, (x, 0, oo), meijerg=True)',
  164. 'integrate(x**2*dagum, (x, 0, oo), meijerg=True)',
  165. 'integrate(f, (x, 0, oo), meijerg=True)',
  166. 'integrate(x*f, (x, 0, oo), meijerg=True)',
  167. 'integrate(x**2*f, (x, 0, oo), meijerg=True)',
  168. 'integrate(rice, (x, 0, oo), meijerg=True)',
  169. 'integrate(laplace, (x, -oo, oo), meijerg=True)',
  170. 'integrate(x*laplace, (x, -oo, oo), meijerg=True)',
  171. 'integrate(x**2*laplace, (x, -oo, oo), meijerg=True)',
  172. 'integrate(log(x) * x**(k-1) * exp(-x) / gamma(k), (x, 0, oo))',
  173. 'integrate(sin(z*x)*(x**2-1)**(-(y+S(1)/2)), (x, 1, oo), meijerg=True)',
  174. 'integrate(besselj(0,x)*besselj(1,x)*exp(-x**2), (x, 0, oo), meijerg=True)',
  175. 'integrate(besselj(0,x)*besselj(1,x)*besselk(0,x), (x, 0, oo), meijerg=True)',
  176. 'integrate(besselj(0,x)*besselj(1,x)*exp(-x**2), (x, 0, oo), meijerg=True)',
  177. 'integrate(besselj(a,x)*besselj(b,x)/x, (x,0,oo), meijerg=True)',
  178. 'hyperexpand(meijerg((-s - a/2 + 1, -s + a/2 + 1), (-a/2 - S(1)/2, -s + a/2 + S(3)/2), (a/2, -a/2), (-a/2 - S(1)/2, -s + a/2 + S(3)/2), 1))',
  179. "gammasimp(S('2**(2*s)*(-pi*gamma(-a + 1)*gamma(a + 1)*gamma(-a - s + 1)*gamma(-a + s - 1/2)*gamma(a - s + 3/2)*gamma(a + s + 1)/(a*(a + s)) - gamma(-a - 1/2)*gamma(-a + 1)*gamma(a + 1)*gamma(a + 3/2)*gamma(-s + 3/2)*gamma(s - 1/2)*gamma(-a + s + 1)*gamma(a - s + 1)/(a*(-a + s)))*gamma(-2*s + 1)*gamma(s + 1)/(pi*s*gamma(-a - 1/2)*gamma(a + 3/2)*gamma(-s + 1)*gamma(-s + 3/2)*gamma(s - 1/2)*gamma(-a - s + 1)*gamma(-a + s - 1/2)*gamma(a - s + 1)*gamma(a - s + 3/2))'))",
  180. 'mellin_transform(E1(x), x, s)',
  181. 'inverse_mellin_transform(gamma(s)/s, s, x, (0, oo))',
  182. 'mellin_transform(expint(a, x), x, s)',
  183. 'mellin_transform(Si(x), x, s)',
  184. 'inverse_mellin_transform(-2**s*sqrt(pi)*gamma((s + 1)/2)/(2*s*gamma(-s/2 + 1)), s, x, (-1, 0))',
  185. 'mellin_transform(Ci(sqrt(x)), x, s)',
  186. 'inverse_mellin_transform(-4**s*sqrt(pi)*gamma(s)/(2*s*gamma(-s + S(1)/2)),s, u, (0, 1))',
  187. 'laplace_transform(Ci(x), x, s)',
  188. 'laplace_transform(expint(a, x), x, s)',
  189. 'laplace_transform(expint(1, x), x, s)',
  190. 'laplace_transform(expint(2, x), x, s)',
  191. 'inverse_laplace_transform(-log(1 + s**2)/2/s, s, u)',
  192. 'inverse_laplace_transform(log(s + 1)/s, s, x)',
  193. 'inverse_laplace_transform((s - log(s + 1))/s**2, s, x)',
  194. 'laplace_transform(Chi(x), x, s)',
  195. 'laplace_transform(Shi(x), x, s)',
  196. 'integrate(exp(-z*x)/x, (x, 1, oo), meijerg=True, conds="none")',
  197. 'integrate(exp(-z*x)/x**2, (x, 1, oo), meijerg=True, conds="none")',
  198. 'integrate(exp(-z*x)/x**3, (x, 1, oo), meijerg=True,conds="none")',
  199. 'integrate(-cos(x)/x, (x, tpos, oo), meijerg=True)',
  200. 'integrate(-sin(x)/x, (x, tpos, oo), meijerg=True)',
  201. 'integrate(sin(x)/x, (x, 0, z), meijerg=True)',
  202. 'integrate(sinh(x)/x, (x, 0, z), meijerg=True)',
  203. 'integrate(exp(-x)/x, x, meijerg=True)',
  204. 'integrate(exp(-x)/x**2, x, meijerg=True)',
  205. 'integrate(cos(u)/u, u, meijerg=True)',
  206. 'integrate(cosh(u)/u, u, meijerg=True)',
  207. 'integrate(expint(1, x), x, meijerg=True)',
  208. 'integrate(expint(2, x), x, meijerg=True)',
  209. 'integrate(Si(x), x, meijerg=True)',
  210. 'integrate(Ci(u), u, meijerg=True)',
  211. 'integrate(Shi(x), x, meijerg=True)',
  212. 'integrate(Chi(u), u, meijerg=True)',
  213. 'integrate(Si(x)*exp(-x), (x, 0, oo), meijerg=True)',
  214. 'integrate(expint(1, x)*sin(x), (x, 0, oo), meijerg=True)'
  215. ]
  216. from time import time
  217. from sympy.core.cache import clear_cache
  218. import sys
  219. timings = []
  220. if __name__ == '__main__':
  221. for n, string in enumerate(bench):
  222. clear_cache()
  223. _t = time()
  224. exec(string)
  225. _t = time() - _t
  226. timings += [(_t, string)]
  227. sys.stdout.write('.')
  228. sys.stdout.flush()
  229. if n % (len(bench) // 10) == 0:
  230. sys.stdout.write('%s' % (10*n // len(bench)))
  231. print()
  232. timings.sort(key=lambda x: -x[0])
  233. for ti, string in timings:
  234. print('%.2fs %s' % (ti, string))