asd
This commit is contained in:
@ -0,0 +1,17 @@
|
||||
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)
|
||||
@ -0,0 +1,54 @@
|
||||
"""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(f"#define nA {len(A)}\n")
|
||||
for k, Ak in enumerate(A):
|
||||
', '.join([str(x.evalf(18)) for x in Ak.coeffs()])
|
||||
f.write(f"static const double A{k}[] = {{tmp}};\n")
|
||||
", ".join([f"A{k}" for k in range(K + 1)])
|
||||
f.write("static const double *A[] = {{tmp}};\n")
|
||||
", ".join([str(Ak.degree()) for Ak in A])
|
||||
f.write("static const int Adegs[] = {{tmp}};\n")
|
||||
os.rename(fn + '.new', fn)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@ -0,0 +1,116 @@
|
||||
"""
|
||||
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()
|
||||
@ -0,0 +1,124 @@
|
||||
"""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',
|
||||
f'{func.__name__}.txt')
|
||||
np.savetxt(filename, dataset)
|
||||
|
||||
print(f"{(time() - t0)/60} minutes elapsed")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@ -0,0 +1,484 @@
|
||||
"""This script evaluates scipy's implementation of hyp2f1 against mpmath's.
|
||||
|
||||
Author: Albert Steppi
|
||||
|
||||
This script is long running and generates a large output file. With default
|
||||
arguments, the generated file is roughly 700MB in size and it takes around
|
||||
40 minutes using an Intel(R) Core(TM) i5-8250U CPU with n_jobs set to 8
|
||||
(full utilization). There are optional arguments which can be used to restrict
|
||||
(or enlarge) the computations performed. These are described below.
|
||||
The output of this script can be analyzed to identify suitable test cases and
|
||||
to find parameter and argument regions where hyp2f1 needs to be improved.
|
||||
|
||||
The script has one mandatory positional argument for specifying the path to
|
||||
the location where the output file is to be placed, and 4 optional arguments
|
||||
--n_jobs, --grid_size, --regions, and --parameter_groups. --n_jobs specifies
|
||||
the number of processes to use if running in parallel. The default value is 1.
|
||||
The other optional arguments are explained below.
|
||||
|
||||
Produces a tab separated values file with 11 columns. The first four columns
|
||||
contain the parameters a, b, c and the argument z. The next two contain |z| and
|
||||
a region code for which region of the complex plane belongs to. The regions are
|
||||
|
||||
0) z == 1
|
||||
1) |z| < 0.9 and real(z) >= 0
|
||||
2) |z| <= 1 and real(z) < 0
|
||||
3) 0.9 <= |z| <= 1 and |1 - z| < 0.9:
|
||||
4) 0.9 <= |z| <= 1 and |1 - z| >= 0.9 and real(z) >= 0:
|
||||
5) 1 < |z| < 1.1 and |1 - z| >= 0.9 and real(z) >= 0
|
||||
6) |z| > 1 and not in 5)
|
||||
|
||||
The --regions optional argument allows the user to specify a list of regions
|
||||
to which computation will be restricted.
|
||||
|
||||
Parameters a, b, c are taken from a 10 * 10 * 10 grid with values at
|
||||
|
||||
-16, -8, -4, -2, -1, 1, 2, 4, 8, 16
|
||||
|
||||
with random perturbations applied.
|
||||
|
||||
There are 9 parameter groups handling the following cases.
|
||||
|
||||
1) A, B, C, B - A, C - A, C - B, C - A - B all non-integral.
|
||||
2) B - A integral
|
||||
3) C - A integral
|
||||
4) C - B integral
|
||||
5) C - A - B integral
|
||||
6) A integral
|
||||
7) B integral
|
||||
8) C integral
|
||||
9) Wider range with c - a - b > 0.
|
||||
|
||||
The seventh column of the output file is an integer between 1 and 8 specifying
|
||||
the parameter group as above.
|
||||
|
||||
The --parameter_groups optional argument allows the user to specify a list of
|
||||
parameter groups to which computation will be restricted.
|
||||
|
||||
The argument z is taken from a grid in the box
|
||||
-box_size <= real(z) <= box_size, -box_size <= imag(z) <= box_size.
|
||||
with grid size specified using the optional command line argument --grid_size,
|
||||
and box_size specified with the command line argument --box_size.
|
||||
The default value of grid_size is 20 and the default value of box_size is 2.0,
|
||||
yielding a 20 * 20 grid in the box with corners -2-2j, -2+2j, 2-2j, 2+2j.
|
||||
|
||||
The final four columns have the expected value of hyp2f1 for the given
|
||||
parameters and argument as calculated with mpmath, the observed value
|
||||
calculated with scipy's hyp2f1, the relative error, and the absolute error.
|
||||
|
||||
As special cases of hyp2f1 are moved from the original Fortran implementation
|
||||
into Cython, this script can be used to ensure that no regressions occur and
|
||||
to point out where improvements are needed.
|
||||
"""
|
||||
|
||||
|
||||
import os
|
||||
import csv
|
||||
import argparse
|
||||
import numpy as np
|
||||
from itertools import product
|
||||
from multiprocessing import Pool
|
||||
|
||||
|
||||
from scipy.special import hyp2f1
|
||||
from scipy.special.tests.test_hyp2f1 import mp_hyp2f1
|
||||
|
||||
|
||||
def get_region(z):
|
||||
"""Assign numbers for regions where hyp2f1 must be handled differently."""
|
||||
if z == 1 + 0j:
|
||||
return 0
|
||||
elif abs(z) < 0.9 and z.real >= 0:
|
||||
return 1
|
||||
elif abs(z) <= 1 and z.real < 0:
|
||||
return 2
|
||||
elif 0.9 <= abs(z) <= 1 and abs(1 - z) < 0.9:
|
||||
return 3
|
||||
elif 0.9 <= abs(z) <= 1 and abs(1 - z) >= 0.9:
|
||||
return 4
|
||||
elif 1 < abs(z) < 1.1 and abs(1 - z) >= 0.9 and z.real >= 0:
|
||||
return 5
|
||||
else:
|
||||
return 6
|
||||
|
||||
|
||||
def get_result(a, b, c, z, group):
|
||||
"""Get results for given parameter and value combination."""
|
||||
expected, observed = mp_hyp2f1(a, b, c, z), hyp2f1(a, b, c, z)
|
||||
if (
|
||||
np.isnan(observed) and np.isnan(expected) or
|
||||
expected == observed
|
||||
):
|
||||
relative_error = 0.0
|
||||
absolute_error = 0.0
|
||||
elif np.isnan(observed):
|
||||
# Set error to infinity if result is nan when not expected to be.
|
||||
# Makes results easier to interpret.
|
||||
relative_error = float("inf")
|
||||
absolute_error = float("inf")
|
||||
else:
|
||||
absolute_error = abs(expected - observed)
|
||||
relative_error = absolute_error / abs(expected)
|
||||
|
||||
return (
|
||||
a,
|
||||
b,
|
||||
c,
|
||||
z,
|
||||
abs(z),
|
||||
get_region(z),
|
||||
group,
|
||||
expected,
|
||||
observed,
|
||||
relative_error,
|
||||
absolute_error,
|
||||
)
|
||||
|
||||
|
||||
def get_result_no_mp(a, b, c, z, group):
|
||||
"""Get results for given parameter and value combination."""
|
||||
expected, observed = complex('nan'), hyp2f1(a, b, c, z)
|
||||
relative_error, absolute_error = float('nan'), float('nan')
|
||||
return (
|
||||
a,
|
||||
b,
|
||||
c,
|
||||
z,
|
||||
abs(z),
|
||||
get_region(z),
|
||||
group,
|
||||
expected,
|
||||
observed,
|
||||
relative_error,
|
||||
absolute_error,
|
||||
)
|
||||
|
||||
|
||||
def get_results(params, Z, n_jobs=1, compute_mp=True):
|
||||
"""Batch compute results for multiple parameter and argument values.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
params : iterable
|
||||
iterable of tuples of floats (a, b, c) specifying parameter values
|
||||
a, b, c for hyp2f1
|
||||
Z : iterable of complex
|
||||
Arguments at which to evaluate hyp2f1
|
||||
n_jobs : Optional[int]
|
||||
Number of jobs for parallel execution.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list
|
||||
List of tuples of results values. See return value in source code
|
||||
of `get_result`.
|
||||
"""
|
||||
input_ = (
|
||||
(a, b, c, z, group) for (a, b, c, group), z in product(params, Z)
|
||||
)
|
||||
|
||||
with Pool(n_jobs) as pool:
|
||||
rows = pool.starmap(
|
||||
get_result if compute_mp else get_result_no_mp,
|
||||
input_
|
||||
)
|
||||
return rows
|
||||
|
||||
|
||||
def _make_hyp2f1_test_case(a, b, c, z, rtol):
|
||||
"""Generate string for single test case as used in test_hyp2f1.py."""
|
||||
expected = mp_hyp2f1(a, b, c, z)
|
||||
return (
|
||||
" pytest.param(\n"
|
||||
" Hyp2f1TestCase(\n"
|
||||
f" a={a},\n"
|
||||
f" b={b},\n"
|
||||
f" c={c},\n"
|
||||
f" z={z},\n"
|
||||
f" expected={expected},\n"
|
||||
f" rtol={rtol},\n"
|
||||
" ),\n"
|
||||
" ),"
|
||||
)
|
||||
|
||||
|
||||
def make_hyp2f1_test_cases(rows):
|
||||
"""Generate string for a list of test cases for test_hyp2f1.py.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rows : list
|
||||
List of lists of the form [a, b, c, z, rtol] where a, b, c, z are
|
||||
parameters and the argument for hyp2f1 and rtol is an expected
|
||||
relative error for the associated test case.
|
||||
|
||||
Returns
|
||||
-------
|
||||
str
|
||||
String for a list of test cases. The output string can be printed
|
||||
or saved to a file and then copied into an argument for
|
||||
`pytest.mark.parameterize` within `scipy.special.tests.test_hyp2f1.py`.
|
||||
"""
|
||||
result = "[\n"
|
||||
result += '\n'.join(
|
||||
_make_hyp2f1_test_case(a, b, c, z, rtol)
|
||||
for a, b, c, z, rtol in rows
|
||||
)
|
||||
result += "\n]"
|
||||
return result
|
||||
|
||||
|
||||
def main(
|
||||
outpath,
|
||||
n_jobs=1,
|
||||
box_size=2.0,
|
||||
grid_size=20,
|
||||
regions=None,
|
||||
parameter_groups=None,
|
||||
compute_mp=True,
|
||||
):
|
||||
outpath = os.path.realpath(os.path.expanduser(outpath))
|
||||
|
||||
random_state = np.random.RandomState(1234)
|
||||
# Parameters a, b, c selected near these values.
|
||||
root_params = np.array(
|
||||
[-16, -8, -4, -2, -1, 1, 2, 4, 8, 16]
|
||||
)
|
||||
# Perturbations to apply to root values.
|
||||
perturbations = 0.1 * random_state.random_sample(
|
||||
size=(3, len(root_params))
|
||||
)
|
||||
|
||||
params = []
|
||||
# Parameter group 1
|
||||
# -----------------
|
||||
# No integer differences. This has been confirmed for the above seed.
|
||||
A = root_params + perturbations[0, :]
|
||||
B = root_params + perturbations[1, :]
|
||||
C = root_params + perturbations[2, :]
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 1) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 2
|
||||
# -----------------
|
||||
# B - A an integer
|
||||
A = root_params + 0.5
|
||||
B = root_params + 0.5
|
||||
C = root_params + perturbations[1, :]
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 2) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 3
|
||||
# -----------------
|
||||
# C - A an integer
|
||||
A = root_params + 0.5
|
||||
B = root_params + perturbations[1, :]
|
||||
C = root_params + 0.5
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 3) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 4
|
||||
# -----------------
|
||||
# C - B an integer
|
||||
A = root_params + perturbations[0, :]
|
||||
B = root_params + 0.5
|
||||
C = root_params + 0.5
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 4) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 5
|
||||
# -----------------
|
||||
# C - A - B an integer
|
||||
A = root_params + 0.25
|
||||
B = root_params + 0.25
|
||||
C = root_params + 0.5
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 5) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 6
|
||||
# -----------------
|
||||
# A an integer
|
||||
A = root_params
|
||||
B = root_params + perturbations[0, :]
|
||||
C = root_params + perturbations[1, :]
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 6) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 7
|
||||
# -----------------
|
||||
# B an integer
|
||||
A = root_params + perturbations[0, :]
|
||||
B = root_params
|
||||
C = root_params + perturbations[1, :]
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 7) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 8
|
||||
# -----------------
|
||||
# C an integer
|
||||
A = root_params + perturbations[0, :]
|
||||
B = root_params + perturbations[1, :]
|
||||
C = root_params
|
||||
params.extend(
|
||||
sorted(
|
||||
((a, b, c, 8) for a, b, c in product(A, B, C)),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
)
|
||||
|
||||
# Parameter group 9
|
||||
# -----------------
|
||||
# Wide range of magnitudes, c - a - b > 0.
|
||||
phi = (1 + np.sqrt(5))/2
|
||||
P = phi**np.arange(16)
|
||||
P = np.hstack([-P, P])
|
||||
group_9_params = sorted(
|
||||
(
|
||||
(a, b, c, 9) for a, b, c in product(P, P, P) if c - a - b > 0
|
||||
),
|
||||
key=lambda x: max(abs(x[0]), abs(x[1])),
|
||||
)
|
||||
|
||||
if parameter_groups is not None:
|
||||
# Group 9 params only used if specified in arguments.
|
||||
params.extend(group_9_params)
|
||||
params = [
|
||||
(a, b, c, group) for a, b, c, group in params
|
||||
if group in parameter_groups
|
||||
]
|
||||
|
||||
# grid_size * grid_size grid in box with corners
|
||||
# -2 - 2j, -2 + 2j, 2 - 2j, 2 + 2j
|
||||
X, Y = np.meshgrid(
|
||||
np.linspace(-box_size, box_size, grid_size),
|
||||
np.linspace(-box_size, box_size, grid_size)
|
||||
)
|
||||
Z = X + Y * 1j
|
||||
Z = Z.flatten().tolist()
|
||||
# Add z = 1 + 0j (region 0).
|
||||
Z.append(1 + 0j)
|
||||
if regions is not None:
|
||||
Z = [z for z in Z if get_region(z) in regions]
|
||||
|
||||
# Evaluate scipy and mpmath's hyp2f1 for all parameter combinations
|
||||
# above against all arguments in the grid Z
|
||||
rows = get_results(params, Z, n_jobs=n_jobs, compute_mp=compute_mp)
|
||||
|
||||
with open(outpath, "w", newline="") as f:
|
||||
writer = csv.writer(f, delimiter="\t")
|
||||
writer.writerow(
|
||||
[
|
||||
"a",
|
||||
"b",
|
||||
"c",
|
||||
"z",
|
||||
"|z|",
|
||||
"region",
|
||||
"parameter_group",
|
||||
"expected", # mpmath's hyp2f1
|
||||
"observed", # scipy's hyp2f1
|
||||
"relative_error",
|
||||
"absolute_error",
|
||||
]
|
||||
)
|
||||
for row in rows:
|
||||
writer.writerow(row)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
parser = argparse.ArgumentParser(
|
||||
description="Test scipy's hyp2f1 against mpmath's on a grid in the"
|
||||
" complex plane over a grid of parameter values. Saves output to file"
|
||||
" specified in positional argument \"outpath\"."
|
||||
" Caution: With default arguments, the generated output file is"
|
||||
" roughly 700MB in size. Script may take several hours to finish if"
|
||||
" \"--n_jobs\" is set to 1."
|
||||
)
|
||||
parser.add_argument(
|
||||
"outpath", type=str, help="Path to output tsv file."
|
||||
)
|
||||
parser.add_argument(
|
||||
"--n_jobs",
|
||||
type=int,
|
||||
default=1,
|
||||
help="Number of jobs for multiprocessing.",
|
||||
)
|
||||
parser.add_argument(
|
||||
"--box_size",
|
||||
type=float,
|
||||
default=2.0,
|
||||
help="hyp2f1 is evaluated in box of side_length 2*box_size centered"
|
||||
" at the origin."
|
||||
)
|
||||
parser.add_argument(
|
||||
"--grid_size",
|
||||
type=int,
|
||||
default=20,
|
||||
help="hyp2f1 is evaluated on grid_size * grid_size grid in box of side"
|
||||
" length 2*box_size centered at the origin."
|
||||
)
|
||||
parser.add_argument(
|
||||
"--parameter_groups",
|
||||
type=int,
|
||||
nargs='+',
|
||||
default=None,
|
||||
help="Restrict to supplied parameter groups. See the Docstring for"
|
||||
" this module for more info on parameter groups. Calculate for all"
|
||||
" parameter groups by default."
|
||||
)
|
||||
parser.add_argument(
|
||||
"--regions",
|
||||
type=int,
|
||||
nargs='+',
|
||||
default=None,
|
||||
help="Restrict to argument z only within the supplied regions. See"
|
||||
" the Docstring for this module for more info on regions. Calculate"
|
||||
" for all regions by default."
|
||||
)
|
||||
parser.add_argument(
|
||||
"--no_mp",
|
||||
action='store_true',
|
||||
help="If this flag is set, do not compute results with mpmath. Saves"
|
||||
" time if results have already been computed elsewhere. Fills in"
|
||||
" \"expected\" column with None values."
|
||||
)
|
||||
args = parser.parse_args()
|
||||
compute_mp = not args.no_mp
|
||||
print(args.parameter_groups)
|
||||
main(
|
||||
args.outpath,
|
||||
n_jobs=args.n_jobs,
|
||||
box_size=args.box_size,
|
||||
grid_size=args.grid_size,
|
||||
parameter_groups=args.parameter_groups,
|
||||
regions=args.regions,
|
||||
compute_mp=compute_mp,
|
||||
)
|
||||
@ -0,0 +1,68 @@
|
||||
"""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(f"p = {p}")
|
||||
print(f"q = {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()
|
||||
@ -0,0 +1,43 @@
|
||||
"""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()
|
||||
@ -0,0 +1,131 @@
|
||||
"""
|
||||
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:
|
||||
def sh(v, z):
|
||||
return float(mpmath.struveh(mpmath.mpf(v), mpmath.mpf(z)))
|
||||
else:
|
||||
def sh(v, z):
|
||||
return 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()
|
||||
@ -0,0 +1,38 @@
|
||||
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
|
||||
@ -0,0 +1,342 @@
|
||||
"""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 convention: 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 convention: 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] alone."
|
||||
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 arc length of a function f from 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([f'{x:.5g}' 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(f"\n{(time() - t0)/60:.1f} minutes elapsed.\n")
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
@ -0,0 +1,152 @@
|
||||
"""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(f"{(time() - t0)/60:.1f} minutes elapsed")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@ -0,0 +1,41 @@
|
||||
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()
|
||||
@ -0,0 +1,27 @@
|
||||
"""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