add padding function to imgScalePadding()
This commit is contained in:
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -1,18 +0,0 @@
|
||||
|
||||
import mpmath
|
||||
|
||||
|
||||
def f(x):
|
||||
return (mpmath.pi + x + mpmath.sin(x)) / (2*mpmath.pi)
|
||||
|
||||
|
||||
# Note: 40 digits might be overkill; a few more digits than the default
|
||||
# might be sufficient.
|
||||
mpmath.mp.dps = 40
|
||||
ts = mpmath.taylor(f, -mpmath.pi, 20)
|
||||
p, q = mpmath.pade(ts, 9, 10)
|
||||
|
||||
p = [float(c) for c in p]
|
||||
q = [float(c) for c in q]
|
||||
print('p =', p)
|
||||
print('q =', q)
|
||||
@@ -1,54 +0,0 @@
|
||||
"""Precompute the polynomials for the asymptotic expansion of the
|
||||
generalized exponential integral.
|
||||
|
||||
Sources
|
||||
-------
|
||||
[1] NIST, Digital Library of Mathematical Functions,
|
||||
https://dlmf.nist.gov/8.20#ii
|
||||
|
||||
"""
|
||||
import os
|
||||
|
||||
try:
|
||||
import sympy
|
||||
from sympy import Poly
|
||||
x = sympy.symbols('x')
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def generate_A(K):
|
||||
A = [Poly(1, x)]
|
||||
for k in range(K):
|
||||
A.append(Poly(1 - 2*k*x, x)*A[k] + Poly(x*(x + 1))*A[k].diff())
|
||||
return A
|
||||
|
||||
|
||||
WARNING = """\
|
||||
/* This file was automatically generated by _precompute/expn_asy.py.
|
||||
* Do not edit it manually!
|
||||
*/
|
||||
"""
|
||||
|
||||
|
||||
def main():
|
||||
print(__doc__)
|
||||
fn = os.path.join('..', 'cephes', 'expn.h')
|
||||
|
||||
K = 12
|
||||
A = generate_A(K)
|
||||
with open(fn + '.new', 'w') as f:
|
||||
f.write(WARNING)
|
||||
f.write("#define nA {}\n".format(len(A)))
|
||||
for k, Ak in enumerate(A):
|
||||
tmp = ', '.join([str(x.evalf(18)) for x in Ak.coeffs()])
|
||||
f.write("static const double A{}[] = {{{}}};\n".format(k, tmp))
|
||||
tmp = ", ".join(["A{}".format(k) for k in range(K + 1)])
|
||||
f.write("static const double *A[] = {{{}}};\n".format(tmp))
|
||||
tmp = ", ".join([str(Ak.degree()) for Ak in A])
|
||||
f.write("static const int Adegs[] = {{{}}};\n".format(tmp))
|
||||
os.rename(fn + '.new', fn)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -1,116 +0,0 @@
|
||||
"""
|
||||
Precompute coefficients of Temme's asymptotic expansion for gammainc.
|
||||
|
||||
This takes about 8 hours to run on a 2.3 GHz Macbook Pro with 4GB ram.
|
||||
|
||||
Sources:
|
||||
[1] NIST, "Digital Library of Mathematical Functions",
|
||||
https://dlmf.nist.gov/
|
||||
|
||||
"""
|
||||
import os
|
||||
from scipy.special._precompute.utils import lagrange_inversion
|
||||
|
||||
try:
|
||||
import mpmath as mp
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def compute_a(n):
|
||||
"""a_k from DLMF 5.11.6"""
|
||||
a = [mp.sqrt(2)/2]
|
||||
for k in range(1, n):
|
||||
ak = a[-1]/k
|
||||
for j in range(1, len(a)):
|
||||
ak -= a[j]*a[-j]/(j + 1)
|
||||
ak /= a[0]*(1 + mp.mpf(1)/(k + 1))
|
||||
a.append(ak)
|
||||
return a
|
||||
|
||||
|
||||
def compute_g(n):
|
||||
"""g_k from DLMF 5.11.3/5.11.5"""
|
||||
a = compute_a(2*n)
|
||||
g = [mp.sqrt(2)*mp.rf(0.5, k)*a[2*k] for k in range(n)]
|
||||
return g
|
||||
|
||||
|
||||
def eta(lam):
|
||||
"""Function from DLMF 8.12.1 shifted to be centered at 0."""
|
||||
if lam > 0:
|
||||
return mp.sqrt(2*(lam - mp.log(lam + 1)))
|
||||
elif lam < 0:
|
||||
return -mp.sqrt(2*(lam - mp.log(lam + 1)))
|
||||
else:
|
||||
return 0
|
||||
|
||||
|
||||
def compute_alpha(n):
|
||||
"""alpha_n from DLMF 8.12.13"""
|
||||
coeffs = mp.taylor(eta, 0, n - 1)
|
||||
return lagrange_inversion(coeffs)
|
||||
|
||||
|
||||
def compute_d(K, N):
|
||||
"""d_{k, n} from DLMF 8.12.12"""
|
||||
M = N + 2*K
|
||||
d0 = [-mp.mpf(1)/3]
|
||||
alpha = compute_alpha(M + 2)
|
||||
for n in range(1, M):
|
||||
d0.append((n + 2)*alpha[n+2])
|
||||
d = [d0]
|
||||
g = compute_g(K)
|
||||
for k in range(1, K):
|
||||
dk = []
|
||||
for n in range(M - 2*k):
|
||||
dk.append((-1)**k*g[k]*d[0][n] + (n + 2)*d[k-1][n+2])
|
||||
d.append(dk)
|
||||
for k in range(K):
|
||||
d[k] = d[k][:N]
|
||||
return d
|
||||
|
||||
|
||||
header = \
|
||||
r"""/* This file was automatically generated by _precomp/gammainc.py.
|
||||
* Do not edit it manually!
|
||||
*/
|
||||
|
||||
#ifndef IGAM_H
|
||||
#define IGAM_H
|
||||
|
||||
#define K {}
|
||||
#define N {}
|
||||
|
||||
static const double d[K][N] =
|
||||
{{"""
|
||||
|
||||
footer = \
|
||||
r"""
|
||||
#endif
|
||||
"""
|
||||
|
||||
|
||||
def main():
|
||||
print(__doc__)
|
||||
K = 25
|
||||
N = 25
|
||||
with mp.workdps(50):
|
||||
d = compute_d(K, N)
|
||||
fn = os.path.join(os.path.dirname(__file__), '..', 'cephes', 'igam.h')
|
||||
with open(fn + '.new', 'w') as f:
|
||||
f.write(header.format(K, N))
|
||||
for k, row in enumerate(d):
|
||||
row = [mp.nstr(x, 17, min_fixed=0, max_fixed=0) for x in row]
|
||||
f.write('{')
|
||||
f.write(", ".join(row))
|
||||
if k < K - 1:
|
||||
f.write('},\n')
|
||||
else:
|
||||
f.write('}};\n')
|
||||
f.write(footer)
|
||||
os.rename(fn + '.new', fn)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -1,124 +0,0 @@
|
||||
"""Compute gammainc and gammaincc for large arguments and parameters
|
||||
and save the values to data files for use in tests. We can't just
|
||||
compare to mpmath's gammainc in test_mpmath.TestSystematic because it
|
||||
would take too long.
|
||||
|
||||
Note that mpmath's gammainc is computed using hypercomb, but since it
|
||||
doesn't allow the user to increase the maximum number of terms used in
|
||||
the series it doesn't converge for many arguments. To get around this
|
||||
we copy the mpmath implementation but use more terms.
|
||||
|
||||
This takes about 17 minutes to run on a 2.3 GHz Macbook Pro with 4GB
|
||||
ram.
|
||||
|
||||
Sources:
|
||||
[1] Fredrik Johansson and others. mpmath: a Python library for
|
||||
arbitrary-precision floating-point arithmetic (version 0.19),
|
||||
December 2013. http://mpmath.org/.
|
||||
|
||||
"""
|
||||
import os
|
||||
from time import time
|
||||
import numpy as np
|
||||
from numpy import pi
|
||||
|
||||
from scipy.special._mptestutils import mpf2float
|
||||
|
||||
try:
|
||||
import mpmath as mp
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def gammainc(a, x, dps=50, maxterms=10**8):
|
||||
"""Compute gammainc exactly like mpmath does but allow for more
|
||||
summands in hypercomb. See
|
||||
|
||||
mpmath/functions/expintegrals.py#L134
|
||||
|
||||
in the mpmath github repository.
|
||||
|
||||
"""
|
||||
with mp.workdps(dps):
|
||||
z, a, b = mp.mpf(a), mp.mpf(x), mp.mpf(x)
|
||||
G = [z]
|
||||
negb = mp.fneg(b, exact=True)
|
||||
|
||||
def h(z):
|
||||
T1 = [mp.exp(negb), b, z], [1, z, -1], [], G, [1], [1+z], b
|
||||
return (T1,)
|
||||
|
||||
res = mp.hypercomb(h, [z], maxterms=maxterms)
|
||||
return mpf2float(res)
|
||||
|
||||
|
||||
def gammaincc(a, x, dps=50, maxterms=10**8):
|
||||
"""Compute gammaincc exactly like mpmath does but allow for more
|
||||
terms in hypercomb. See
|
||||
|
||||
mpmath/functions/expintegrals.py#L187
|
||||
|
||||
in the mpmath github repository.
|
||||
|
||||
"""
|
||||
with mp.workdps(dps):
|
||||
z, a = a, x
|
||||
|
||||
if mp.isint(z):
|
||||
try:
|
||||
# mpmath has a fast integer path
|
||||
return mpf2float(mp.gammainc(z, a=a, regularized=True))
|
||||
except mp.libmp.NoConvergence:
|
||||
pass
|
||||
nega = mp.fneg(a, exact=True)
|
||||
G = [z]
|
||||
# Use 2F0 series when possible; fall back to lower gamma representation
|
||||
try:
|
||||
def h(z):
|
||||
r = z-1
|
||||
return [([mp.exp(nega), a], [1, r], [], G, [1, -r], [], 1/nega)]
|
||||
return mpf2float(mp.hypercomb(h, [z], force_series=True))
|
||||
except mp.libmp.NoConvergence:
|
||||
def h(z):
|
||||
T1 = [], [1, z-1], [z], G, [], [], 0
|
||||
T2 = [-mp.exp(nega), a, z], [1, z, -1], [], G, [1], [1+z], a
|
||||
return T1, T2
|
||||
return mpf2float(mp.hypercomb(h, [z], maxterms=maxterms))
|
||||
|
||||
|
||||
def main():
|
||||
t0 = time()
|
||||
# It would be nice to have data for larger values, but either this
|
||||
# requires prohibitively large precision (dps > 800) or mpmath has
|
||||
# a bug. For example, gammainc(1e20, 1e20, dps=800) returns a
|
||||
# value around 0.03, while the true value should be close to 0.5
|
||||
# (DLMF 8.12.15).
|
||||
print(__doc__)
|
||||
pwd = os.path.dirname(__file__)
|
||||
r = np.logspace(4, 14, 30)
|
||||
ltheta = np.logspace(np.log10(pi/4), np.log10(np.arctan(0.6)), 30)
|
||||
utheta = np.logspace(np.log10(pi/4), np.log10(np.arctan(1.4)), 30)
|
||||
|
||||
regimes = [(gammainc, ltheta), (gammaincc, utheta)]
|
||||
for func, theta in regimes:
|
||||
rg, thetag = np.meshgrid(r, theta)
|
||||
a, x = rg*np.cos(thetag), rg*np.sin(thetag)
|
||||
a, x = a.flatten(), x.flatten()
|
||||
dataset = []
|
||||
for i, (a0, x0) in enumerate(zip(a, x)):
|
||||
if func == gammaincc:
|
||||
# Exploit the fast integer path in gammaincc whenever
|
||||
# possible so that the computation doesn't take too
|
||||
# long
|
||||
a0, x0 = np.floor(a0), np.floor(x0)
|
||||
dataset.append((a0, x0, func(a0, x0)))
|
||||
dataset = np.array(dataset)
|
||||
filename = os.path.join(pwd, '..', 'tests', 'data', 'local',
|
||||
'{}.txt'.format(func.__name__))
|
||||
np.savetxt(filename, dataset)
|
||||
|
||||
print("{} minutes elapsed".format((time() - t0)/60))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -1,68 +0,0 @@
|
||||
"""Compute a Pade approximation for the principal branch of the
|
||||
Lambert W function around 0 and compare it to various other
|
||||
approximations.
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
|
||||
try:
|
||||
import mpmath
|
||||
import matplotlib.pyplot as plt
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def lambertw_pade():
|
||||
derivs = [mpmath.diff(mpmath.lambertw, 0, n=n) for n in range(6)]
|
||||
p, q = mpmath.pade(derivs, 3, 2)
|
||||
return p, q
|
||||
|
||||
|
||||
def main():
|
||||
print(__doc__)
|
||||
with mpmath.workdps(50):
|
||||
p, q = lambertw_pade()
|
||||
p, q = p[::-1], q[::-1]
|
||||
print("p = {}".format(p))
|
||||
print("q = {}".format(q))
|
||||
|
||||
x, y = np.linspace(-1.5, 1.5, 75), np.linspace(-1.5, 1.5, 75)
|
||||
x, y = np.meshgrid(x, y)
|
||||
z = x + 1j*y
|
||||
lambertw_std = []
|
||||
for z0 in z.flatten():
|
||||
lambertw_std.append(complex(mpmath.lambertw(z0)))
|
||||
lambertw_std = np.array(lambertw_std).reshape(x.shape)
|
||||
|
||||
fig, axes = plt.subplots(nrows=3, ncols=1)
|
||||
# Compare Pade approximation to true result
|
||||
p = np.array([float(p0) for p0 in p])
|
||||
q = np.array([float(q0) for q0 in q])
|
||||
pade_approx = np.polyval(p, z)/np.polyval(q, z)
|
||||
pade_err = abs(pade_approx - lambertw_std)
|
||||
axes[0].pcolormesh(x, y, pade_err)
|
||||
# Compare two terms of asymptotic series to true result
|
||||
asy_approx = np.log(z) - np.log(np.log(z))
|
||||
asy_err = abs(asy_approx - lambertw_std)
|
||||
axes[1].pcolormesh(x, y, asy_err)
|
||||
# Compare two terms of the series around the branch point to the
|
||||
# true result
|
||||
p = np.sqrt(2*(np.exp(1)*z + 1))
|
||||
series_approx = -1 + p - p**2/3
|
||||
series_err = abs(series_approx - lambertw_std)
|
||||
im = axes[2].pcolormesh(x, y, series_err)
|
||||
|
||||
fig.colorbar(im, ax=axes.ravel().tolist())
|
||||
plt.show()
|
||||
|
||||
fig, ax = plt.subplots(nrows=1, ncols=1)
|
||||
pade_better = pade_err < asy_err
|
||||
im = ax.pcolormesh(x, y, pade_better)
|
||||
t = np.linspace(-0.3, 0.3)
|
||||
ax.plot(-2.5*abs(t) - 0.2, t, 'r')
|
||||
fig.colorbar(im, ax=ax)
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@@ -1,43 +0,0 @@
|
||||
"""Precompute series coefficients for log-Gamma."""
|
||||
|
||||
try:
|
||||
import mpmath
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def stirling_series(N):
|
||||
with mpmath.workdps(100):
|
||||
coeffs = [mpmath.bernoulli(2*n)/(2*n*(2*n - 1))
|
||||
for n in range(1, N + 1)]
|
||||
return coeffs
|
||||
|
||||
|
||||
def taylor_series_at_1(N):
|
||||
coeffs = []
|
||||
with mpmath.workdps(100):
|
||||
coeffs.append(-mpmath.euler)
|
||||
for n in range(2, N + 1):
|
||||
coeffs.append((-1)**n*mpmath.zeta(n)/n)
|
||||
return coeffs
|
||||
|
||||
|
||||
def main():
|
||||
print(__doc__)
|
||||
print()
|
||||
stirling_coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0)
|
||||
for x in stirling_series(8)[::-1]]
|
||||
taylor_coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0)
|
||||
for x in taylor_series_at_1(23)[::-1]]
|
||||
print("Stirling series coefficients")
|
||||
print("----------------------------")
|
||||
print("\n".join(stirling_coeffs))
|
||||
print()
|
||||
print("Taylor series coefficients")
|
||||
print("--------------------------")
|
||||
print("\n".join(taylor_coeffs))
|
||||
print()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@@ -1,120 +0,0 @@
|
||||
"""
|
||||
Convergence regions of the expansions used in ``struve.c``
|
||||
|
||||
Note that for v >> z both functions tend rapidly to 0,
|
||||
and for v << -z, they tend to infinity.
|
||||
|
||||
The floating-point functions over/underflow in the lower left and right
|
||||
corners of the figure.
|
||||
|
||||
|
||||
Figure legend
|
||||
=============
|
||||
|
||||
Red region
|
||||
Power series is close (1e-12) to the mpmath result
|
||||
|
||||
Blue region
|
||||
Asymptotic series is close to the mpmath result
|
||||
|
||||
Green region
|
||||
Bessel series is close to the mpmath result
|
||||
|
||||
Dotted colored lines
|
||||
Boundaries of the regions
|
||||
|
||||
Solid colored lines
|
||||
Boundaries estimated by the routine itself. These will be used
|
||||
for determining which of the results to use.
|
||||
|
||||
Black dashed line
|
||||
The line z = 0.7*|v| + 12
|
||||
|
||||
"""
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
import mpmath
|
||||
|
||||
|
||||
def err_metric(a, b, atol=1e-290):
|
||||
m = abs(a - b) / (atol + abs(b))
|
||||
m[np.isinf(b) & (a == b)] = 0
|
||||
return m
|
||||
|
||||
|
||||
def do_plot(is_h=True):
|
||||
from scipy.special._ufuncs import (_struve_power_series,
|
||||
_struve_asymp_large_z,
|
||||
_struve_bessel_series)
|
||||
|
||||
vs = np.linspace(-1000, 1000, 91)
|
||||
zs = np.sort(np.r_[1e-5, 1.0, np.linspace(0, 700, 91)[1:]])
|
||||
|
||||
rp = _struve_power_series(vs[:,None], zs[None,:], is_h)
|
||||
ra = _struve_asymp_large_z(vs[:,None], zs[None,:], is_h)
|
||||
rb = _struve_bessel_series(vs[:,None], zs[None,:], is_h)
|
||||
|
||||
mpmath.mp.dps = 50
|
||||
if is_h:
|
||||
sh = lambda v, z: float(mpmath.struveh(mpmath.mpf(v), mpmath.mpf(z)))
|
||||
else:
|
||||
sh = lambda v, z: float(mpmath.struvel(mpmath.mpf(v), mpmath.mpf(z)))
|
||||
ex = np.vectorize(sh, otypes='d')(vs[:,None], zs[None,:])
|
||||
|
||||
err_a = err_metric(ra[0], ex) + 1e-300
|
||||
err_p = err_metric(rp[0], ex) + 1e-300
|
||||
err_b = err_metric(rb[0], ex) + 1e-300
|
||||
|
||||
err_est_a = abs(ra[1]/ra[0])
|
||||
err_est_p = abs(rp[1]/rp[0])
|
||||
err_est_b = abs(rb[1]/rb[0])
|
||||
|
||||
z_cutoff = 0.7*abs(vs) + 12
|
||||
|
||||
levels = [-1000, -12]
|
||||
|
||||
plt.cla()
|
||||
|
||||
plt.hold(1)
|
||||
plt.contourf(vs, zs, np.log10(err_p).T, levels=levels, colors=['r', 'r'], alpha=0.1)
|
||||
plt.contourf(vs, zs, np.log10(err_a).T, levels=levels, colors=['b', 'b'], alpha=0.1)
|
||||
plt.contourf(vs, zs, np.log10(err_b).T, levels=levels, colors=['g', 'g'], alpha=0.1)
|
||||
|
||||
plt.contour(vs, zs, np.log10(err_p).T, levels=levels, colors=['r', 'r'], linestyles=[':', ':'])
|
||||
plt.contour(vs, zs, np.log10(err_a).T, levels=levels, colors=['b', 'b'], linestyles=[':', ':'])
|
||||
plt.contour(vs, zs, np.log10(err_b).T, levels=levels, colors=['g', 'g'], linestyles=[':', ':'])
|
||||
|
||||
lp = plt.contour(vs, zs, np.log10(err_est_p).T, levels=levels, colors=['r', 'r'], linestyles=['-', '-'])
|
||||
la = plt.contour(vs, zs, np.log10(err_est_a).T, levels=levels, colors=['b', 'b'], linestyles=['-', '-'])
|
||||
lb = plt.contour(vs, zs, np.log10(err_est_b).T, levels=levels, colors=['g', 'g'], linestyles=['-', '-'])
|
||||
|
||||
plt.clabel(lp, fmt={-1000: 'P', -12: 'P'})
|
||||
plt.clabel(la, fmt={-1000: 'A', -12: 'A'})
|
||||
plt.clabel(lb, fmt={-1000: 'B', -12: 'B'})
|
||||
|
||||
plt.plot(vs, z_cutoff, 'k--')
|
||||
|
||||
plt.xlim(vs.min(), vs.max())
|
||||
plt.ylim(zs.min(), zs.max())
|
||||
|
||||
plt.xlabel('v')
|
||||
plt.ylabel('z')
|
||||
|
||||
|
||||
def main():
|
||||
plt.clf()
|
||||
plt.subplot(121)
|
||||
do_plot(True)
|
||||
plt.title('Struve H')
|
||||
|
||||
plt.subplot(122)
|
||||
do_plot(False)
|
||||
plt.title('Struve L')
|
||||
|
||||
plt.savefig('struve_convergence.png')
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -1,38 +0,0 @@
|
||||
try:
|
||||
import mpmath as mp
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
try:
|
||||
from sympy.abc import x
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def lagrange_inversion(a):
|
||||
"""Given a series
|
||||
|
||||
f(x) = a[1]*x + a[2]*x**2 + ... + a[n-1]*x**(n - 1),
|
||||
|
||||
use the Lagrange inversion formula to compute a series
|
||||
|
||||
g(x) = b[1]*x + b[2]*x**2 + ... + b[n-1]*x**(n - 1)
|
||||
|
||||
so that f(g(x)) = g(f(x)) = x mod x**n. We must have a[0] = 0, so
|
||||
necessarily b[0] = 0 too.
|
||||
|
||||
The algorithm is naive and could be improved, but speed isn't an
|
||||
issue here and it's easy to read.
|
||||
|
||||
"""
|
||||
n = len(a)
|
||||
f = sum(a[i]*x**i for i in range(n))
|
||||
h = (x/f).series(x, 0, n).removeO()
|
||||
hpower = [h**0]
|
||||
for k in range(n):
|
||||
hpower.append((hpower[-1]*h).expand())
|
||||
b = [mp.mpf(0)]
|
||||
for k in range(1, n):
|
||||
b.append(hpower[k].coeff(x, k - 1)/k)
|
||||
b = [mp.mpf(x) for x in b]
|
||||
return b
|
||||
@@ -1,342 +0,0 @@
|
||||
"""Precompute coefficients of several series expansions
|
||||
of Wright's generalized Bessel function Phi(a, b, x).
|
||||
|
||||
See https://dlmf.nist.gov/10.46.E1 with rho=a, beta=b, z=x.
|
||||
"""
|
||||
from argparse import ArgumentParser, RawTextHelpFormatter
|
||||
import numpy as np
|
||||
from scipy.integrate import quad
|
||||
from scipy.optimize import minimize_scalar, curve_fit
|
||||
from time import time
|
||||
|
||||
try:
|
||||
import sympy
|
||||
from sympy import EulerGamma, Rational, S, Sum, \
|
||||
factorial, gamma, gammasimp, pi, polygamma, symbols, zeta
|
||||
from sympy.polys.polyfuncs import horner
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def series_small_a():
|
||||
"""Tylor series expansion of Phi(a, b, x) in a=0 up to order 5.
|
||||
"""
|
||||
order = 5
|
||||
a, b, x, k = symbols("a b x k")
|
||||
A = [] # terms with a
|
||||
X = [] # terms with x
|
||||
B = [] # terms with b (polygammas)
|
||||
# Phi(a, b, x) = exp(x)/gamma(b) * sum(A[i] * X[i] * B[i])
|
||||
expression = Sum(x**k/factorial(k)/gamma(a*k+b), (k, 0, S.Infinity))
|
||||
expression = gamma(b)/sympy.exp(x) * expression
|
||||
|
||||
# nth term of taylor series in a=0: a^n/n! * (d^n Phi(a, b, x)/da^n at a=0)
|
||||
for n in range(0, order+1):
|
||||
term = expression.diff(a, n).subs(a, 0).simplify().doit()
|
||||
# set the whole bracket involving polygammas to 1
|
||||
x_part = (term.subs(polygamma(0, b), 1)
|
||||
.replace(polygamma, lambda *args: 0))
|
||||
# sign convetion: x part always positive
|
||||
x_part *= (-1)**n
|
||||
|
||||
A.append(a**n/factorial(n))
|
||||
X.append(horner(x_part))
|
||||
B.append(horner((term/x_part).simplify()))
|
||||
|
||||
s = "Tylor series expansion of Phi(a, b, x) in a=0 up to order 5.\n"
|
||||
s += "Phi(a, b, x) = exp(x)/gamma(b) * sum(A[i] * X[i] * B[i], i=0..5)\n"
|
||||
for name, c in zip(['A', 'X', 'B'], [A, X, B]):
|
||||
for i in range(len(c)):
|
||||
s += f"\n{name}[{i}] = " + str(c[i])
|
||||
return s
|
||||
|
||||
|
||||
# expansion of digamma
|
||||
def dg_series(z, n):
|
||||
"""Symbolic expansion of digamma(z) in z=0 to order n.
|
||||
|
||||
See https://dlmf.nist.gov/5.7.E4 and with https://dlmf.nist.gov/5.5.E2
|
||||
"""
|
||||
k = symbols("k")
|
||||
return -1/z - EulerGamma + \
|
||||
sympy.summation((-1)**k * zeta(k) * z**(k-1), (k, 2, n+1))
|
||||
|
||||
|
||||
def pg_series(k, z, n):
|
||||
"""Symbolic expansion of polygamma(k, z) in z=0 to order n."""
|
||||
return sympy.diff(dg_series(z, n+k), z, k)
|
||||
|
||||
|
||||
def series_small_a_small_b():
|
||||
"""Tylor series expansion of Phi(a, b, x) in a=0 and b=0 up to order 5.
|
||||
|
||||
Be aware of cancellation of poles in b=0 of digamma(b)/Gamma(b) and
|
||||
polygamma functions.
|
||||
|
||||
digamma(b)/Gamma(b) = -1 - 2*M_EG*b + O(b^2)
|
||||
digamma(b)^2/Gamma(b) = 1/b + 3*M_EG + b*(-5/12*PI^2+7/2*M_EG^2) + O(b^2)
|
||||
polygamma(1, b)/Gamma(b) = 1/b + M_EG + b*(1/12*PI^2 + 1/2*M_EG^2) + O(b^2)
|
||||
and so on.
|
||||
"""
|
||||
order = 5
|
||||
a, b, x, k = symbols("a b x k")
|
||||
M_PI, M_EG, M_Z3 = symbols("M_PI M_EG M_Z3")
|
||||
c_subs = {pi: M_PI, EulerGamma: M_EG, zeta(3): M_Z3}
|
||||
A = [] # terms with a
|
||||
X = [] # terms with x
|
||||
B = [] # terms with b (polygammas expanded)
|
||||
C = [] # terms that generate B
|
||||
# Phi(a, b, x) = exp(x) * sum(A[i] * X[i] * B[i])
|
||||
# B[0] = 1
|
||||
# B[k] = sum(C[k] * b**k/k!, k=0..)
|
||||
# Note: C[k] can be obtained from a series expansion of 1/gamma(b).
|
||||
expression = gamma(b)/sympy.exp(x) * \
|
||||
Sum(x**k/factorial(k)/gamma(a*k+b), (k, 0, S.Infinity))
|
||||
|
||||
# nth term of taylor series in a=0: a^n/n! * (d^n Phi(a, b, x)/da^n at a=0)
|
||||
for n in range(0, order+1):
|
||||
term = expression.diff(a, n).subs(a, 0).simplify().doit()
|
||||
# set the whole bracket involving polygammas to 1
|
||||
x_part = (term.subs(polygamma(0, b), 1)
|
||||
.replace(polygamma, lambda *args: 0))
|
||||
# sign convetion: x part always positive
|
||||
x_part *= (-1)**n
|
||||
# expansion of polygamma part with 1/gamma(b)
|
||||
pg_part = term/x_part/gamma(b)
|
||||
if n >= 1:
|
||||
# Note: highest term is digamma^n
|
||||
pg_part = pg_part.replace(polygamma,
|
||||
lambda k, x: pg_series(k, x, order+1+n))
|
||||
pg_part = (pg_part.series(b, 0, n=order+1-n)
|
||||
.removeO()
|
||||
.subs(polygamma(2, 1), -2*zeta(3))
|
||||
.simplify()
|
||||
)
|
||||
|
||||
A.append(a**n/factorial(n))
|
||||
X.append(horner(x_part))
|
||||
B.append(pg_part)
|
||||
|
||||
# Calculate C and put in the k!
|
||||
C = sympy.Poly(B[1].subs(c_subs), b).coeffs()
|
||||
C.reverse()
|
||||
for i in range(len(C)):
|
||||
C[i] = (C[i] * factorial(i)).simplify()
|
||||
|
||||
s = "Tylor series expansion of Phi(a, b, x) in a=0 and b=0 up to order 5."
|
||||
s += "\nPhi(a, b, x) = exp(x) * sum(A[i] * X[i] * B[i], i=0..5)\n"
|
||||
s += "B[0] = 1\n"
|
||||
s += "B[i] = sum(C[k+i-1] * b**k/k!, k=0..)\n"
|
||||
s += "\nM_PI = pi"
|
||||
s += "\nM_EG = EulerGamma"
|
||||
s += "\nM_Z3 = zeta(3)"
|
||||
for name, c in zip(['A', 'X'], [A, X]):
|
||||
for i in range(len(c)):
|
||||
s += f"\n{name}[{i}] = "
|
||||
s += str(c[i])
|
||||
# For C, do also compute the values numerically
|
||||
for i in range(len(C)):
|
||||
s += f"\n# C[{i}] = "
|
||||
s += str(C[i])
|
||||
s += f"\nC[{i}] = "
|
||||
s += str(C[i].subs({M_EG: EulerGamma, M_PI: pi, M_Z3: zeta(3)})
|
||||
.evalf(17))
|
||||
|
||||
# Does B have the assumed structure?
|
||||
s += "\n\nTest if B[i] does have the assumed structure."
|
||||
s += "\nC[i] are derived from B[1] allone."
|
||||
s += "\nTest B[2] == C[1] + b*C[2] + b^2/2*C[3] + b^3/6*C[4] + .."
|
||||
test = sum([b**k/factorial(k) * C[k+1] for k in range(order-1)])
|
||||
test = (test - B[2].subs(c_subs)).simplify()
|
||||
s += f"\ntest successful = {test==S(0)}"
|
||||
s += "\nTest B[3] == C[2] + b*C[3] + b^2/2*C[4] + .."
|
||||
test = sum([b**k/factorial(k) * C[k+2] for k in range(order-2)])
|
||||
test = (test - B[3].subs(c_subs)).simplify()
|
||||
s += f"\ntest successful = {test==S(0)}"
|
||||
return s
|
||||
|
||||
|
||||
def asymptotic_series():
|
||||
"""Asymptotic expansion for large x.
|
||||
|
||||
Phi(a, b, x) ~ Z^(1/2-b) * exp((1+a)/a * Z) * sum_k (-1)^k * C_k / Z^k
|
||||
Z = (a*x)^(1/(1+a))
|
||||
|
||||
Wright (1935) lists the coefficients C_0 and C_1 (he calls them a_0 and
|
||||
a_1). With slightly different notation, Paris (2017) lists coefficients
|
||||
c_k up to order k=3.
|
||||
Paris (2017) uses ZP = (1+a)/a * Z (ZP = Z of Paris) and
|
||||
C_k = C_0 * (-a/(1+a))^k * c_k
|
||||
"""
|
||||
order = 8
|
||||
|
||||
class g(sympy.Function):
|
||||
"""Helper function g according to Wright (1935)
|
||||
|
||||
g(n, rho, v) = (1 + (rho+2)/3 * v + (rho+2)*(rho+3)/(2*3) * v^2 + ...)
|
||||
|
||||
Note: Wright (1935) uses square root of above definition.
|
||||
"""
|
||||
nargs = 3
|
||||
|
||||
@classmethod
|
||||
def eval(cls, n, rho, v):
|
||||
if not n >= 0:
|
||||
raise ValueError("must have n >= 0")
|
||||
elif n == 0:
|
||||
return 1
|
||||
else:
|
||||
return g(n-1, rho, v) \
|
||||
+ gammasimp(gamma(rho+2+n)/gamma(rho+2)) \
|
||||
/ gammasimp(gamma(3+n)/gamma(3))*v**n
|
||||
|
||||
class coef_C(sympy.Function):
|
||||
"""Calculate coefficients C_m for integer m.
|
||||
|
||||
C_m is the coefficient of v^(2*m) in the Taylor expansion in v=0 of
|
||||
Gamma(m+1/2)/(2*pi) * (2/(rho+1))^(m+1/2) * (1-v)^(-b)
|
||||
* g(rho, v)^(-m-1/2)
|
||||
"""
|
||||
nargs = 3
|
||||
|
||||
@classmethod
|
||||
def eval(cls, m, rho, beta):
|
||||
if not m >= 0:
|
||||
raise ValueError("must have m >= 0")
|
||||
|
||||
v = symbols("v")
|
||||
expression = (1-v)**(-beta) * g(2*m, rho, v)**(-m-Rational(1, 2))
|
||||
res = expression.diff(v, 2*m).subs(v, 0) / factorial(2*m)
|
||||
res = res * (gamma(m + Rational(1, 2)) / (2*pi)
|
||||
* (2/(rho+1))**(m + Rational(1, 2)))
|
||||
return res
|
||||
|
||||
# in order to have nice ordering/sorting of expressions, we set a = xa.
|
||||
xa, b, xap1 = symbols("xa b xap1")
|
||||
C0 = coef_C(0, xa, b)
|
||||
# a1 = a(1, rho, beta)
|
||||
s = "Asymptotic expansion for large x\n"
|
||||
s += "Phi(a, b, x) = Z**(1/2-b) * exp((1+a)/a * Z) \n"
|
||||
s += " * sum((-1)**k * C[k]/Z**k, k=0..6)\n\n"
|
||||
s += "Z = pow(a * x, 1/(1+a))\n"
|
||||
s += "A[k] = pow(a, k)\n"
|
||||
s += "B[k] = pow(b, k)\n"
|
||||
s += "Ap1[k] = pow(1+a, k)\n\n"
|
||||
s += "C[0] = 1./sqrt(2. * M_PI * Ap1[1])\n"
|
||||
for i in range(1, order+1):
|
||||
expr = (coef_C(i, xa, b) / (C0/(1+xa)**i)).simplify()
|
||||
factor = [x.denominator() for x in sympy.Poly(expr).coeffs()]
|
||||
factor = sympy.lcm(factor)
|
||||
expr = (expr * factor).simplify().collect(b, sympy.factor)
|
||||
expr = expr.xreplace({xa+1: xap1})
|
||||
s += f"C[{i}] = C[0] / ({factor} * Ap1[{i}])\n"
|
||||
s += f"C[{i}] *= {str(expr)}\n\n"
|
||||
import re
|
||||
re_a = re.compile(r'xa\*\*(\d+)')
|
||||
s = re_a.sub(r'A[\1]', s)
|
||||
re_b = re.compile(r'b\*\*(\d+)')
|
||||
s = re_b.sub(r'B[\1]', s)
|
||||
s = s.replace('xap1', 'Ap1[1]')
|
||||
s = s.replace('xa', 'a')
|
||||
# max integer = 2^31-1 = 2,147,483,647. Solution: Put a point after 10
|
||||
# or more digits.
|
||||
re_digits = re.compile(r'(\d{10,})')
|
||||
s = re_digits.sub(r'\1.', s)
|
||||
return s
|
||||
|
||||
|
||||
def optimal_epsilon_integral():
|
||||
"""Fit optimal choice of epsilon for integral representation.
|
||||
|
||||
The integrand of
|
||||
int_0^pi P(eps, a, b, x, phi) * dphi
|
||||
can exhibit oscillatory behaviour. It stems from the cosine of P and can be
|
||||
minimized by minimizing the arc length of the argument
|
||||
f(phi) = eps * sin(phi) - x * eps^(-a) * sin(a * phi) + (1 - b) * phi
|
||||
of cos(f(phi)).
|
||||
We minimize the arc length in eps for a grid of values (a, b, x) and fit a
|
||||
parametric function to it.
|
||||
"""
|
||||
def fp(eps, a, b, x, phi):
|
||||
"""Derivative of f w.r.t. phi."""
|
||||
eps_a = np.power(1. * eps, -a)
|
||||
return eps * np.cos(phi) - a * x * eps_a * np.cos(a * phi) + 1 - b
|
||||
|
||||
def arclength(eps, a, b, x, epsrel=1e-2, limit=100):
|
||||
"""Compute Arc length of f.
|
||||
|
||||
Note that the arg length of a function f fro t0 to t1 is given by
|
||||
int_t0^t1 sqrt(1 + f'(t)^2) dt
|
||||
"""
|
||||
return quad(lambda phi: np.sqrt(1 + fp(eps, a, b, x, phi)**2),
|
||||
0, np.pi,
|
||||
epsrel=epsrel, limit=100)[0]
|
||||
|
||||
# grid of minimal arc length values
|
||||
data_a = [1e-3, 0.1, 0.5, 0.9, 1, 2, 4, 5, 6, 8]
|
||||
data_b = [0, 1, 4, 7, 10]
|
||||
data_x = [1, 1.5, 2, 4, 10, 20, 50, 100, 200, 500, 1e3, 5e3, 1e4]
|
||||
data_a, data_b, data_x = np.meshgrid(data_a, data_b, data_x)
|
||||
data_a, data_b, data_x = (data_a.flatten(), data_b.flatten(),
|
||||
data_x.flatten())
|
||||
best_eps = []
|
||||
for i in range(data_x.size):
|
||||
best_eps.append(
|
||||
minimize_scalar(lambda eps: arclength(eps, data_a[i], data_b[i],
|
||||
data_x[i]),
|
||||
bounds=(1e-3, 1000),
|
||||
method='Bounded', options={'xatol': 1e-3}).x
|
||||
)
|
||||
best_eps = np.array(best_eps)
|
||||
# pandas would be nice, but here a dictionary is enough
|
||||
df = {'a': data_a,
|
||||
'b': data_b,
|
||||
'x': data_x,
|
||||
'eps': best_eps,
|
||||
}
|
||||
|
||||
def func(data, A0, A1, A2, A3, A4, A5):
|
||||
"""Compute parametric function to fit."""
|
||||
a = data['a']
|
||||
b = data['b']
|
||||
x = data['x']
|
||||
return (A0 * b * np.exp(-0.5 * a)
|
||||
+ np.exp(A1 + 1 / (1 + a) * np.log(x) - A2 * np.exp(-A3 * a)
|
||||
+ A4 / (1 + np.exp(A5 * a))))
|
||||
|
||||
func_params = list(curve_fit(func, df, df['eps'], method='trf')[0])
|
||||
|
||||
s = "Fit optimal eps for integrand P via minimal arc length\n"
|
||||
s += "with parametric function:\n"
|
||||
s += "optimal_eps = (A0 * b * exp(-a/2) + exp(A1 + 1 / (1 + a) * log(x)\n"
|
||||
s += " - A2 * exp(-A3 * a) + A4 / (1 + exp(A5 * a)))\n\n"
|
||||
s += "Fitted parameters A0 to A5 are:\n"
|
||||
s += ', '.join(['{:.5g}'.format(x) for x in func_params])
|
||||
return s
|
||||
|
||||
|
||||
def main():
|
||||
t0 = time()
|
||||
parser = ArgumentParser(description=__doc__,
|
||||
formatter_class=RawTextHelpFormatter)
|
||||
parser.add_argument('action', type=int, choices=[1, 2, 3, 4],
|
||||
help='chose what expansion to precompute\n'
|
||||
'1 : Series for small a\n'
|
||||
'2 : Series for small a and small b\n'
|
||||
'3 : Asymptotic series for large x\n'
|
||||
' This may take some time (>4h).\n'
|
||||
'4 : Fit optimal eps for integral representation.'
|
||||
)
|
||||
args = parser.parse_args()
|
||||
|
||||
switch = {1: lambda: print(series_small_a()),
|
||||
2: lambda: print(series_small_a_small_b()),
|
||||
3: lambda: print(asymptotic_series()),
|
||||
4: lambda: print(optimal_epsilon_integral())
|
||||
}
|
||||
switch.get(args.action, lambda: print("Invalid input."))()
|
||||
print("\n{:.1f} minutes elapsed.\n".format((time() - t0)/60))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@@ -1,152 +0,0 @@
|
||||
"""Compute a grid of values for Wright's generalized Bessel function
|
||||
and save the values to data files for use in tests. Using mpmath directly in
|
||||
tests would take too long.
|
||||
|
||||
This takes about 10 minutes to run on a 2.7 GHz i7 Macbook Pro.
|
||||
"""
|
||||
from functools import lru_cache
|
||||
import os
|
||||
from time import time
|
||||
|
||||
import numpy as np
|
||||
from scipy.special._mptestutils import mpf2float
|
||||
|
||||
try:
|
||||
import mpmath as mp
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
# exp_inf: smallest value x for which exp(x) == inf
|
||||
exp_inf = 709.78271289338403
|
||||
|
||||
|
||||
# 64 Byte per value
|
||||
@lru_cache(maxsize=100_000)
|
||||
def rgamma_cached(x, dps):
|
||||
with mp.workdps(dps):
|
||||
return mp.rgamma(x)
|
||||
|
||||
|
||||
def mp_wright_bessel(a, b, x, dps=50, maxterms=2000):
|
||||
"""Compute Wright's generalized Bessel function as Series with mpmath.
|
||||
"""
|
||||
with mp.workdps(dps):
|
||||
a, b, x = mp.mpf(a), mp.mpf(b), mp.mpf(x)
|
||||
res = mp.nsum(lambda k: x**k / mp.fac(k)
|
||||
* rgamma_cached(a * k + b, dps=dps),
|
||||
[0, mp.inf],
|
||||
tol=dps, method='s', steps=[maxterms]
|
||||
)
|
||||
return mpf2float(res)
|
||||
|
||||
|
||||
def main():
|
||||
t0 = time()
|
||||
print(__doc__)
|
||||
pwd = os.path.dirname(__file__)
|
||||
eps = np.finfo(float).eps * 100
|
||||
|
||||
a_range = np.array([eps,
|
||||
1e-4 * (1 - eps), 1e-4, 1e-4 * (1 + eps),
|
||||
1e-3 * (1 - eps), 1e-3, 1e-3 * (1 + eps),
|
||||
0.1, 0.5,
|
||||
1 * (1 - eps), 1, 1 * (1 + eps),
|
||||
1.5, 2, 4.999, 5, 10])
|
||||
b_range = np.array([0, eps, 1e-10, 1e-5, 0.1, 1, 2, 10, 20, 100])
|
||||
x_range = np.array([0, eps, 1 - eps, 1, 1 + eps,
|
||||
1.5,
|
||||
2 - eps, 2, 2 + eps,
|
||||
9 - eps, 9, 9 + eps,
|
||||
10 * (1 - eps), 10, 10 * (1 + eps),
|
||||
100 * (1 - eps), 100, 100 * (1 + eps),
|
||||
500, exp_inf, 1e3, 1e5, 1e10, 1e20])
|
||||
|
||||
a_range, b_range, x_range = np.meshgrid(a_range, b_range, x_range,
|
||||
indexing='ij')
|
||||
a_range = a_range.flatten()
|
||||
b_range = b_range.flatten()
|
||||
x_range = x_range.flatten()
|
||||
|
||||
# filter out some values, especially too large x
|
||||
bool_filter = ~((a_range < 5e-3) & (x_range >= exp_inf))
|
||||
bool_filter = bool_filter & ~((a_range < 0.2) & (x_range > exp_inf))
|
||||
bool_filter = bool_filter & ~((a_range < 0.5) & (x_range > 1e3))
|
||||
bool_filter = bool_filter & ~((a_range < 0.56) & (x_range > 5e3))
|
||||
bool_filter = bool_filter & ~((a_range < 1) & (x_range > 1e4))
|
||||
bool_filter = bool_filter & ~((a_range < 1.4) & (x_range > 1e5))
|
||||
bool_filter = bool_filter & ~((a_range < 1.8) & (x_range > 1e6))
|
||||
bool_filter = bool_filter & ~((a_range < 2.2) & (x_range > 1e7))
|
||||
bool_filter = bool_filter & ~((a_range < 2.5) & (x_range > 1e8))
|
||||
bool_filter = bool_filter & ~((a_range < 2.9) & (x_range > 1e9))
|
||||
bool_filter = bool_filter & ~((a_range < 3.3) & (x_range > 1e10))
|
||||
bool_filter = bool_filter & ~((a_range < 3.7) & (x_range > 1e11))
|
||||
bool_filter = bool_filter & ~((a_range < 4) & (x_range > 1e12))
|
||||
bool_filter = bool_filter & ~((a_range < 4.4) & (x_range > 1e13))
|
||||
bool_filter = bool_filter & ~((a_range < 4.7) & (x_range > 1e14))
|
||||
bool_filter = bool_filter & ~((a_range < 5.1) & (x_range > 1e15))
|
||||
bool_filter = bool_filter & ~((a_range < 5.4) & (x_range > 1e16))
|
||||
bool_filter = bool_filter & ~((a_range < 5.8) & (x_range > 1e17))
|
||||
bool_filter = bool_filter & ~((a_range < 6.2) & (x_range > 1e18))
|
||||
bool_filter = bool_filter & ~((a_range < 6.2) & (x_range > 1e18))
|
||||
bool_filter = bool_filter & ~((a_range < 6.5) & (x_range > 1e19))
|
||||
bool_filter = bool_filter & ~((a_range < 6.9) & (x_range > 1e20))
|
||||
|
||||
# filter out known values that do not meet the required numerical accuracy
|
||||
# see test test_wright_data_grid_failures
|
||||
failing = np.array([
|
||||
[0.1, 100, 709.7827128933841],
|
||||
[0.5, 10, 709.7827128933841],
|
||||
[0.5, 10, 1000],
|
||||
[0.5, 100, 1000],
|
||||
[1, 20, 100000],
|
||||
[1, 100, 100000],
|
||||
[1.0000000000000222, 20, 100000],
|
||||
[1.0000000000000222, 100, 100000],
|
||||
[1.5, 0, 500],
|
||||
[1.5, 2.220446049250313e-14, 500],
|
||||
[1.5, 1.e-10, 500],
|
||||
[1.5, 1.e-05, 500],
|
||||
[1.5, 0.1, 500],
|
||||
[1.5, 20, 100000],
|
||||
[1.5, 100, 100000],
|
||||
]).tolist()
|
||||
|
||||
does_fail = np.full_like(a_range, False, dtype=bool)
|
||||
for i in range(x_range.size):
|
||||
if [a_range[i], b_range[i], x_range[i]] in failing:
|
||||
does_fail[i] = True
|
||||
|
||||
# filter and flatten
|
||||
a_range = a_range[bool_filter]
|
||||
b_range = b_range[bool_filter]
|
||||
x_range = x_range[bool_filter]
|
||||
does_fail = does_fail[bool_filter]
|
||||
|
||||
dataset = []
|
||||
print(f"Computing {x_range.size} single points.")
|
||||
print("Tests will fail for the following data points:")
|
||||
for i in range(x_range.size):
|
||||
a = a_range[i]
|
||||
b = b_range[i]
|
||||
x = x_range[i]
|
||||
# take care of difficult corner cases
|
||||
maxterms = 1000
|
||||
if a < 1e-6 and x >= exp_inf/10:
|
||||
maxterms = 2000
|
||||
f = mp_wright_bessel(a, b, x, maxterms=maxterms)
|
||||
if does_fail[i]:
|
||||
print("failing data point a, b, x, value = "
|
||||
f"[{a}, {b}, {x}, {f}]")
|
||||
else:
|
||||
dataset.append((a, b, x, f))
|
||||
dataset = np.array(dataset)
|
||||
|
||||
filename = os.path.join(pwd, '..', 'tests', 'data', 'local',
|
||||
'wright_bessel.txt')
|
||||
np.savetxt(filename, dataset)
|
||||
|
||||
print("{:.1f} minutes elapsed".format((time() - t0)/60))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -1,41 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
try:
|
||||
import mpmath
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def mpmath_wrightomega(x):
|
||||
return mpmath.lambertw(mpmath.exp(x), mpmath.mpf('-0.5'))
|
||||
|
||||
|
||||
def wrightomega_series_error(x):
|
||||
series = x
|
||||
desired = mpmath_wrightomega(x)
|
||||
return abs(series - desired) / desired
|
||||
|
||||
|
||||
def wrightomega_exp_error(x):
|
||||
exponential_approx = mpmath.exp(x)
|
||||
desired = mpmath_wrightomega(x)
|
||||
return abs(exponential_approx - desired) / desired
|
||||
|
||||
|
||||
def main():
|
||||
desired_error = 2 * np.finfo(float).eps
|
||||
print('Series Error')
|
||||
for x in [1e5, 1e10, 1e15, 1e20]:
|
||||
with mpmath.workdps(100):
|
||||
error = wrightomega_series_error(x)
|
||||
print(x, error, error < desired_error)
|
||||
|
||||
print('Exp error')
|
||||
for x in [-10, -25, -50, -100, -200, -400, -700, -740]:
|
||||
with mpmath.workdps(100):
|
||||
error = wrightomega_exp_error(x)
|
||||
print(x, error, error < desired_error)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@@ -1,27 +0,0 @@
|
||||
"""Compute the Taylor series for zeta(x) - 1 around x = 0."""
|
||||
try:
|
||||
import mpmath
|
||||
except ImportError:
|
||||
pass
|
||||
|
||||
|
||||
def zetac_series(N):
|
||||
coeffs = []
|
||||
with mpmath.workdps(100):
|
||||
coeffs.append(-1.5)
|
||||
for n in range(1, N):
|
||||
coeff = mpmath.diff(mpmath.zeta, 0, n)/mpmath.factorial(n)
|
||||
coeffs.append(coeff)
|
||||
return coeffs
|
||||
|
||||
|
||||
def main():
|
||||
print(__doc__)
|
||||
coeffs = zetac_series(10)
|
||||
coeffs = [mpmath.nstr(x, 20, min_fixed=0, max_fixed=0)
|
||||
for x in coeffs]
|
||||
print("\n".join(coeffs[::-1]))
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
Reference in New Issue
Block a user