<jupyter_start><jupyter_text>**Authors:** <NAME>, <NAME>   *[Faculty of Science](https://www.upjs.sk/en/faculty-of-science/?prefferedLang=EN), P. J. Šafárik University in Košice, Slovakia*  emails: [<EMAIL>](mailto:<EMAIL>)***  Numerical inversion for $\mathcal{GDD}$ pdf calculations Python-based computational tools:   **Python, NumPy, Cython, Numba**    Trapezoidal rule vs. DE quadratureTrapezoidal (complex functions - Python), DE (complex, real functions - Python, Cython, Numba)<jupyter_code>N = 4 #test up to 10^N values of pdf

import numpy as np
import pandas as pd
from numpy import array as v
import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')

from time import time
import math, cmath
import scipy

# accuracy in Python
def abs_errs(df1, df2):
    N = len(df1)
    errors = [abs(df1[i]-df2[i]) for i in range(N)]
    return errors<jupyter_output><empty_output><jupyter_text>$\mathcal{GDD}$ pdf  $-$ analytic expression$f(z)=\dfrac{\beta_{1}^{\alpha_{1}} \beta_{2}^{\alpha_{2}}}{\beta^{\alpha-1}}\begin{cases}{\dfrac{e^{z \beta_{2}}}{\Gamma\left(\alpha_{2}\right)} U\left(1-\alpha_{2}, 2-\alpha,-z \beta\right),} & {z<0} \\[12pt]\begin{array}{cc} \frac{\Gamma(\alpha-1)}{\beta^{a / 2-1} \Gamma\left(\alpha_{1}\right) \Gamma\left(\alpha_{2}\right)}, & \scriptstyle 1<\alpha \\ \infty,  & \scriptstyle 0<\alpha \leq 1 \end{array}, & z=0 \\[12pt] {\dfrac{e^{-z \beta_{1}}}{\Gamma\left(\alpha_{1}\right)} U\left(1-\alpha_{1}, 2-\alpha, z \beta\right),} & {z>0}\end{cases}$$\alpha=\alpha_{1}+\alpha_{2}, \quad \beta=\beta_{1}+\beta_{2}$<jupyter_code># pre-calulated analytic results from mpmath, Pari GP and Arb
dpari = {str(10**(n+1)):np.loadtxt('data/Pari_Sage_pdf'+str(10**(n+1))+'.txt', delimiter=',',
                                   dtype=np.longdouble) for n in range(N)}
dx = {str(10**(n+1)):np.loadtxt('data/xvalues'+str(10**(n+1))+'.txt', dtype=np.longdouble) for n in range(N)}<jupyter_output><empty_output><jupyter_text>Pdf as a numerical inversion integral$ \newcommand{\a}{\alpha} \newcommand{\b}{\beta} \newcommand{\cft}{\varphi(t)}  \newcommand{\Rv}[1]{\mathscr{R}^{1}}$$$\large f(x)=  \frac{1}{\pi} \int_0^\infty \Re \left[e^{-itx}\cft\right]dt $$$$ \normalsize\cft=  \left(1-\dfrac{it}{\b_1}\right)^{-\a_1}\left(1+\dfrac{it}{\b_2}\right)^{-\a_2}, t\in \Rv{}$$ ** alternative form**$$\large f(x)= \int_{0}^{\infty} g(x,t) dt $$$$ \normalsize g(x,t) = r(t)\cos(xt-\phi(t)), $$$$\normalsize r(t)= \frac{\b_1^{\a_1} \b_2^{\a_2}}{\pi\left(\b_1^2+t^2\right)^{\a_1/2}\left(\b_2^2+t^2\right)^{\a_2 / 2}}, \quad \phi(t) = \a_1 \arctan \frac{t}{\b_1}-\a_2 \arctan \frac{t}{\b_2}$$  DE, TR<jupyter_code># Ooura DE quadrature implementation in Python, Numba; folder TRandDE
from TRandDE.DEQuadratureNumba import *
from TRandDE.DEQuadrature import *
scipy.__version__, np.__version__<jupyter_output>C:\Users\jozef\.conda\envs\py38\lib\site-packages\numba\core\cpu.py:78: UserWarning: Numba extension module 'numba_scipy' failed to load due to 'ValueError(No function '__pyx_fuse_0pdtr' found in __pyx_capi__ of 'scipy.special.cython_special')'.
  numba.core.entrypoints.init_all()<jupyter_text>Integrands and parameters<jupyter_code># setting parameters and integrands
eps, epss, epsss = 1e-3, 1e-10, 1e-15
omg = 1

a1, b1, a2, b2 = 0.5, 1.0, 8.5, 93.0
d = b1**a1*b2**a2/math.pi
# integrand Python, complex char function
def gcpy(t, x):
    phit = (1-1j*t/b1)**(-a1)*(1+1j*t/b2)**(-a2)
    return (1/math.pi*cmath.exp(-1j*t*x)*phit).real
# integrand Numba, complex and real char function
@njit(fastmath=True)
def gcn(t, x):
    phit = (1-1j*t/b1)**(-a1)*(1+1j*t/b2)**(-a2)
    return (1/math.pi*cmath.exp(-1j*t*x)*phit).real

@njit(fastmath=True)
def gn(t,x):
    phi = a1*math.atan(t/b1) - a2*math.atan(t/b2)
    r = 1/((b1*b1 + t*t)**(a1/2)*(b2*b2 + t*t)**(a2/2))
    return d*r*math.cos(x*t-phi)
%load_ext Cython
%%cython
cimport scipy.special.cython_special as cs
cimport libc.math as lc

def Gcy(double t, double x, double a1, double b1, double a2, double b2):
    cdef:
        double phi, r, d
    
    phi = a1*lc.atan(t/b1)-a2*lc.atan(t/b2)
    r = 1/(lc.pow(b1*b1 + t*t,a1/2)*lc.pow(b2*b2 + t*t,a2/2))
    d = lc.pow(b1,a1)*lc.pow(b2,a2)/lc.pi
    return d*r*lc.cos(x*t-phi)
# integrad Cython
def gcy(t,x):
    return Gcy(t,x, a1,b1,a2,b2)<jupyter_output><empty_output><jupyter_text>Comparing integrands<jupyter_code>t,x = 0.32, 0.8
# cross-check values
gint = [gcpy, gcn, gn, gcy]
[gcpy(t,x), gcn(t,x), gn(t,x), gcy(t,x)]
# speed of integrands
runtimes = []
for fun in gint:
    print(fun.__name__, '('+str(t)+','+str(x)+') =', fun(t, x))
    %timeit -r 3 -n 10 fun(t, x)
    rt = %timeit -r 1 -oq fun(t, x)
    runtimes += [rt.average] 
    print('')
(runtimes[0]/v(runtimes)).round(2)<jupyter_output><empty_output><jupyter_text>DE implementation<jupyter_code>#python, complex
fde3cpy = lambda x: intdeo(lambda t: gcpy(t,x), 0, x, eps)[0]
fde10cpy = lambda x: intdeo(lambda t: gcpy(t,x), 0, x, epss)[0]
fde15cpy = lambda x: intdeo(lambda t: gcpy(t,x), 0, x, epsss)[0]

# numpy vectorization
fde3cpy_np = np.vectorize(fde3cpy)
fde10cpy_np = np.vectorize(fde10cpy)
fde15cpy_np = np.vectorize(fde15cpy)
intdeogcn = make_intdeo(gcn)
intdeogn = make_intdeo(gn)

# numba, complex
fde3cn = lambda x: intdeogcn(0, x, eps)[0]
fde10cn = lambda x: intdeogcn(0, x, epss)[0]
fde15cn = lambda x: intdeogcn(0, x, epsss)[0]

fde3cn_np = np.vectorize(fde3cn)
fde10cn_np = np.vectorize(fde10cn)
fde15cn_np = np.vectorize(fde15cn)

# numba, float
fde3n = lambda x: intdeogn(0, x, eps)[0]
fde10n = lambda x: intdeogn(0, x, epss)[0]
fde15n = lambda x: intdeogn(0, x, epsss)[0]

fde3n_np = np.vectorize(fde3n)
fde10n_np = np.vectorize(fde10n)
fde15n_np = np.vectorize(fde15n)
# cross-check
u = 0.45
[fde3cpy(u), fde10cpy(u), fde15cpy(u), #python complex
 fde3cn(u), fde10cn(u), fde15cn(u),    #numba, complex
 fde3n(u), fde10n(u), fde15n(u)        #numba, float 
]
# cross-check
[fde3cpy_np(dx['10'][3:7]), fde10cpy_np(dx['10'][3:7]), fde15cpy_np(dx['10'][3:7]), #python complex
 fde3cn_np(dx['10'][3:7]), fde10cn_np(dx['10'][3:7]), fde15cn_np(dx['10'][3:7]),    #numba, complex
 fde3n_np(dx['10'][3:7]), fde10n_np(dx['10'][3:7]), fde15n_np(dx['10'][3:7])        #numba, float 
]<jupyter_output><empty_output><jupyter_text>TR implementation<jupyter_code># Witkovsky TR implementation in Python; folder TRandDE
import TRandDE.cf2DistGP as ninv
# python
def cf1(t, a1, b1, a2, b2):
    result = ((1-1j*t/b1)**(-a1))*((1+1j*t/b2)**(-a2))
    return(result)

def cf(t): 
    result = cf1(t=t, a1=a1, b1=b1, a2=a2, b2=b2)
    return(result)

fwpy_np = lambda x: ninv.cf2DistGP(cf, x)['pdf']
fwpy = np.vectorize(lambda t: ninv.cf2DistGP(cf, x=t)['pdf'][0])
fwpy_np(dx['10'])
fwpy(dx['10'])
%timeit fwpy_np(dx['10000'])
%timeit fwpy(dx['10000'])<jupyter_output>548 ms ± 10.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
3.8 s ± 57.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>Benchmark<jupyter_code># values x over interval <-3,4>
x = dx
#x = lambda n: np.linspace(-3, 4, n)
methods = ['W py', 'W pynp', 'DEcpy 1e-15', 'DEcpy 1e-10', 'DEcpy 1e-3', 
           'DEcn 1e-15', 'DEcn 1e-10', 'DEcn 1e-3', 'DEn 1e-15', 'DEn 1e-10', 'DEn 1e-3']
benchmark = pd.DataFrame(index =['10', '100', '1000', '10000'], 
            columns = methods)
functions = [fwpy, fwpy_np, fde15cpy_np, fde10cpy_np, fde3cpy_np, 
             fde15cn_np, fde10cn_np, fde3cn_np, 
             fde15n_np, fde10n_np, fde3n_np]
pdf = {m:f for m,f in zip(methods, functions)}
# start time
tic = time()
for n in range(N):
    points = x[str(10**(n+1))]
    for method, fun in pdf.items():
        rt = %timeit -r 3 -n 10 -oq fun(points)
        print(10**(n+1), method, ', runs = ', rt.repeat, ', loops =', rt.loops)
        dm = fun(points)
        dp = dpari[str(10**(n+1))]
        ch_abs = abs_errs(dm,dp)
        #ch_rel = rel_errs(dm,dp)
        benchmark.loc[str(10**(n+1)),method] = [rt.average, rt.stdev, min(ch_abs), max(ch_abs)]
    #end time
    toc = time()-tic; 
    print(10**(n+1),' runtime =',toc,'s')
    print(25*'*')<jupyter_output>10 W py , runs =  3 , loops = 10
10 W pynp , runs =  3 , loops = 10
10 DEcpy 1e-15 , runs =  3 , loops = 10
10 DEcpy 1e-10 , runs =  3 , loops = 10
10 DEcpy 1e-3 , runs =  3 , loops = 10
10 DEcn 1e-15 , runs =  3 , loops = 10
10 DEcn 1e-10 , runs =  3 , loops = 10
10 DEcn 1e-3 , runs =  3 , loops = 10
10 DEn 1e-15 , runs =  3 , loops = 10
10 DEn 1e-10 , runs =  3 , loops = 10
10 DEn 1e-3 , runs =  3 , loops = 10
10  runtime = 1.0413157939910889 s
*************************
100 W py , runs =  3 , loops = 10
100 W pynp , runs =  3 , loops = 10
100 DEcpy 1e-15 , runs =  3 , loops = 10
100 DEcpy 1e-10 , runs =  3 , loops = 10
100 DEcpy 1e-3 , runs =  3 , loops = 10
100 DEcn 1e-15 , runs =  3 , loops = 10
100 DEcn 1e-10 , runs =  3 , loops = 10
100 DEcn 1e-3 , runs =  3 , loops = 10
100 DEn 1e-15 , runs =  3 , loops = 10
100 DEn 1e-10 , runs =  3 , loops = 10
100 DEn 1e-3 , runs =  3 , loops = 10
100  runtime = 10.280387163162231 s
*************************
1000 W py , runs =  3 , loops = 10[...]<jupyter_text>Results<jupyter_code>scn = lambda series, prec: [np.format_float_scientific(item, precision=prec) for item in v(series)]
benchres = benchmark
show_allrowscols(benchres, fullcolwidth=True)
benchres.to_excel('PdfNumericTRDEcomparison.xlsx')<jupyter_output><empty_output><jupyter_text>Run Times<jupyter_code># runtimes for methods and tools
times = benchres.applymap(lambda x: x[0])
times.loc[:,times.columns.str.contains('W')] = times.loc[:,times.columns.str.contains('W')]/2
# average runtimes
times
# relative stdev = variation coefficients in %
delta = benchres.applymap(lambda x:x[1]/x[0])
delta.multiply(100)
tmax = times.max()
tmax
times.loc[:,times.columns.str.contains('15') | times.columns.str.contains('W')].apply(lambda x: scn(x,2))<jupyter_output><empty_output><jupyter_text>Accelerations<jupyter_code>Accel = times.applymap(lambda x: times.loc['10000','W py']/x).round(2)
Accel
Accel.loc[:,Accel.columns.str.contains('15') | Accel.columns.str.contains('W') | Accel.columns.str.contains('DEn')].loc[['10000']]
# reference runtime for pdf 
tWpy = 3.7066538833333325/2 # times.loc['10000','W py']
tWpy
Accel = times.applymap(lambda x: tWpy/x).round(2)
Accel<jupyter_output><empty_output><jupyter_text>Errors<jupyter_code>errs = benchres.applymap(lambda x: scn(x[2:4], 0))
errs.loc[:,errs.columns.str.contains('15') | errs.columns.str.contains('W') | errs.columns.str.contains('DEn')]<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>xh = dx['1000']
for method in pdf.keys():
    yh = [pdf[method](x) for x in xh]
    plt.plot(xh,yh)
    plt.title(method)
    plt.show();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19964499252615844, 'fraction_numerical': 0.06707772795216742, 'mean_word_length': 4.425747592498733, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 56, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '84461', 'license': 'MIT', 'n_tokens_mistral': 5026, 'n_tokens_neox': 4655, 'n_words': 959, 'path': 'Jupyter/3.3 fx pdf NumInv - Trapezoidal rule, DE quadrature - Python tools.ipynb', 'repo_name': 'fdslrm/GDD', 'size': '10749'}