Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Owen's T function #483

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open

Owen's T function #483

wants to merge 6 commits into from

Conversation

jbrea
Copy link

@jbrea jbrea commented Jan 24, 2025

Following up on #242 and JuliaStats/StatsFuns.jl#99.

This implementation should not have any license issues, because it is based on JuliaStats/StatsFuns.jl#99 (comment).

I compared the implementation to scipy with the following code. The Julia version is usually 1.5-3 times faster than scipy; outliers are due to shortcut implementations in the Julia version.

julia> using SpecialFunctions

julia> using PythonCall, DataFrames

julia> function timeit_py(a, h; number = 100_000, repeat = 5)
           pyexec("import timeit", Main)
           pyexec("import scipy", Main)
           times = pyeval(Vector{Float64}, "timeit.repeat(\"scipy.special.owens_t($a, $h)\", \"import scipy\", number=$number, repeat=$repeat)", Main)
           result = pyeval(Float64, "scipy.special.owens_t($a, $h)", Main)
           (; result, times)
       end
timeit_py (generic function with 1 method)

julia> function repeat_owent(a, h; number = 100_000)
           res = 0.
           for _ in 1:number
               res += owent(a, h)
           end
           res
       end
repeat_owent (generic function with 1 method)

julia> function timeit_jl(a, h; number = 100_000, repeat = 5)
           times = [@elapsed(repeat_owent(a, h; number)) for _ in 1:repeat]
           result = owent(a, h)
           (; result, times)
       end
timeit_jl (generic function with 1 method)

julia> result = DataFrame(a = [], h = [],
                          pyres = [], pytime = [],
                          jlres = [], jltime = []);

julia> for a in Any[0; 1; .9999999; randn(10); randn(Float32, 10)], h in Any[randn(10); randn(Float32, 10)]
           jlres = timeit_jl(a, h)
           pyres = timeit_py(a, h)
           push!(result, (a, h,
                          pyres.result, minimum(pyres.times),
                          jlres.result, minimum(jlres.times)))
       end

julia> result.differences = result.pyres - result.jlres;

julia> sort!(result, :differences)
460×7 DataFrame
 Row │ a           h          pyres        pytime     jlres        jltime       differences 
     │ Any         Any        Any          Any        Any          Any          Float64     
─────┼──────────────────────────────────────────────────────────────────────────────────────
   11           1.65593    0.0769836    0.0957145  0.0769836    0.0327059    -5.58371e-8
   21.8347      1.36702    0.016564     0.0848207  0.0165641    0.0415702    -5.40338e-8
   3-0.564564   2.21172    0.136236     0.07811    0.136236     0.0336247    -4.65142e-8
   42.29411     -1.34786   -0.00544201  0.0880784  -0.00544198  0.043255     -3.64403e-8
   50.280253    -0.597672  -0.0820831   0.0718078  -0.0820831   0.0329566    -3.28791e-8
   61.8347      2.15715    0.0166373    0.0827304  0.0166373    0.0421463    -3.28324e-8
   72.29411     -1.72526   -0.00544587  0.0840716  -0.00544584  0.0420629    -3.19978e-8
   82.29411     1.35615    0.00544229   0.0847247  0.00544232   0.0415061    -3.13535e-8
   90.135936    1.41673    0.150109     0.0735957  0.150109     0.0337166    -3.08821e-8
  101.8347      -1.319     -0.0165396   0.0839895  -0.0165396   0.0416502    -2.82707e-8
  110.280253    -0.631865  -0.0858245   0.0741887  -0.0858245   0.0316616    -1.8001e-8
  120.359144    1.1168     0.122817     0.074642   0.122817     0.0341947    -1.78076e-8
  130           -2.29284   -0.184545    0.0858694  -0.184545    0.000674779  -1.70548e-8
  140.135936    -0.840065  -0.109973    0.0744032  -0.109973    0.0320735    -1.52763e-8
  15-0.564564   -0.431946  -0.0548148   0.0714174  -0.0548148   0.0321709    -1.4652e-8
  160.00193627  1.25215    0.142745     0.073243   0.142745     0.0333754    -1.41235e-8
  17-0.564564   -0.285026  -0.0375227   0.0717169  -0.0375227   0.0319744    -1.36557e-8
  18-0.288648   -1.27571   -0.135996    0.0735466  -0.135996    0.0344008    -1.23922e-8
  190.359144    -0.698553  -0.0901374   0.0750477  -0.0901373   0.0318714    -1.22161e-8
  201.65177     -1.22102   -0.0241997   0.083956   -0.0241997   0.0390595    -1.12128e-8
  210.00193627  -0.610624  -0.0872478   0.0740678  -0.0872477   0.0321092    -1.0688e-8
  220.359144    -0.749227  -0.0949342   0.0734939  -0.0949342   0.0316966    -1.00754e-8
                                                                      
 4400.867207    0.402013   0.0409707    0.0751883  0.0409707    0.0315258     6.88819e-9
 4410.359144    0.174748   0.0257975    0.0694627  0.0257975    0.0318224     7.03036e-9
 4421           0.804224   0.0597782    0.0899204  0.0597782    0.0306946     7.04313e-9
 4430.359144    -0.959696  -0.11228     0.0717851  -0.11228     0.0320928     7.60505e-9
 4441           0.883784   0.0629212    0.0914889  0.0629212    0.0304683     7.90769e-9
 4451           0.290823   0.0269482    0.0882529  0.0269482    0.0333106     8.30915e-9
 4460.00193627  -0.678365  -0.0948653   0.0699623  -0.0948653   0.0321516     9.25489e-9
 4470.135936    0.594174   0.0844577    0.0732297  0.0844577    0.0316222     1.09388e-8
 448-0.288648   0.642862   0.0867688    0.0713109  0.0867688    0.0334977     1.17517e-8
 449-0.288648   1.06609    0.12322      0.0747424  0.12322      0.0341674     1.29416e-8
 4500.280253    1.18627    0.131323     0.0730108  0.131323     0.033275      1.6007e-8
 4511.65177     -1.17536   -0.0240919   0.0851801  -0.024092    0.0383333     1.78919e-8
 4521.65177     1.08506    0.0238082    0.0858568  0.0238082    0.0379493     1.90357e-8
 453-0.288648   0.514274   0.0722652    0.0720627  0.0722652    0.0318471     1.97314e-8
 4541           -1.38741   -0.0744888   0.0945712  -0.0744888   0.0332795     2.29497e-8
 4550.359144    -0.511267  -0.0701518   0.07315    -0.0701519   0.0334086     2.43759e-8
 4560.280253    -1.19764   -0.132002    0.0745334  -0.132003    0.0333158     2.58694e-8
 4571.8347      1.52049    0.0166094    0.0868804  0.0166094    0.0414594     3.30341e-8
 4580.00193627  1.21585    0.140455     0.0722593  0.140455     0.0333129     3.5323e-8
 4590.867207    -2.17019   -0.0952003   0.0843306  -0.0952004   0.0367136     5.17318e-8
 4600.280253    1.17843    0.130849     0.0745051  0.130849     0.0337025     5.89902e-8
                                                                            417 rows omitted

julia> result.speedup = result.pytime ./ result.jltime;

julia> sort!(result, :speedup)
460×8 DataFrame
 Row │ a           h           pyres        pytime     jlres        jltime       differences   speedup   
     │ Any         Any         Any          Any        Any          Any          Float64       Float64   
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────
   10.91965     -1.08832    -0.076221    0.0834341  -0.076221    0.0552344     6.25891e-11    1.51055
   20.91965     0.326923    0.0324702    0.0692716  0.0324702    0.0448327    -1.92725e-10    1.54511
   30.153105    -0.855227   -0.111022    0.0712328  -0.111022    0.043876     -1.38778e-17    1.6235
   40.0395757   -0.251856   -0.0392361   0.0733503  -0.0392361   0.0444967    -4.10545e-10    1.64844
   5-1.87946    0.953713    0.014473     0.0805918  0.014473     0.0488578     3.46945e-18    1.64952
   60.91965     -0.620778   -0.0551835   0.0739141  -0.0551835   0.0442225    -1.38778e-17    1.67142
   70.177284    -0.272255   -0.0416299   0.0711275  -0.0416299   0.0422422     1.86172e-10    1.6838
   80.177284    0.351235    0.0528872    0.0722759  0.0528872    0.0416287     1.38778e-17    1.7362
   9-1.87946    -1.7683     -0.0150419   0.0831276  -0.0150419   0.0465378    -2.60209e-17    1.78624
  100.177284    -2.0328     -0.172248    0.074197   -0.172248    0.040407      5.55112e-17    1.83624
  11-0.700541   0.0122647   0.00152716   0.0704533  0.00152716   0.0383499     0.0            1.83712
  12-0.288648   0.0577138   0.00880047   0.0706493  0.00880047   0.0382655     1.8624e-11     1.84629
  13-1.87946    1.16593     0.0148605    0.0894758  0.0148605    0.0483338    -5.86829e-12    1.85121
  140.00193627  0.426092    0.0641067    0.0701449  0.0641067    0.0376274    -2.90046e-15    1.8642
  151.8347      2.9939      0.0166375    0.0823891  0.0166375    0.0436107    -5.77745e-10    1.88919
  160.564926    -0.259144   -0.0342837   0.072101   -0.0342837   0.0380073    -2.06629e-10    1.89703
  170.135936    0.238109    0.0368549    0.0710897  0.0368549    0.0374547    -4.7046e-12     1.89802
  180.517077    0.46296     0.0598289    0.0707332  0.0598289    0.0370159     0.0            1.91089
  19-1.87946    1.08073     0.0147494    0.0877235  0.0147494    0.0454492     9.31177e-11    1.93015
  200.91965     -0.142303   -0.0146974   0.0689571  -0.0146974   0.0356151    -3.46945e-18    1.93618
  21-0.700541   -0.0815047  -0.0101214   0.0811347  -0.0101214   0.041866     -2.2112e-11     1.93796
  22-1.87946    -0.105722   -0.00284766  0.0704058  -0.00284766  0.0362593    -5.3258e-11     1.94173
                                                                                   
 4401           -0.664301   -0.0530681   0.0937582  -0.0530681   0.0308511     2.07798e-9     3.03906
 4410           -2.09427    -0.179099    0.0851074  -0.179099    0.000908579   2.77556e-17   93.6709
 4420           1.09876     0.132483     0.0852197  0.132483     0.000901666  -2.77556e-17   94.5137
 4430           1.37593     0.149975     0.085271   0.149975     0.000892317  -2.77556e-17   95.5614
 4440           -2.2709     -0.183982    0.0862938  -0.183982    0.00089851    0.0           96.041
 4450           1.87557     0.172041     0.0853501  0.172041     0.00088792   -2.77556e-17   96.1237
 4460           -1.1233     -0.134232    0.0856927  -0.134232    0.000875295   0.0           97.9015
 4470           0.869305    0.113891     0.085465   0.113891     0.000870727   0.0           98.1536
 4480           -0.58978    -0.084809    0.0846976  -0.084809    0.000727338   9.15938e-10  116.449
 4490           -0.441735   -0.0662021   0.0852595  -0.0662021   0.000713922  -3.92105e-9   119.424
 4500           0.589119    0.0847309    0.0838691  0.0847309    0.000676782   6.61299e-10  123.923
 4510           -2.29284    -0.184545    0.0858694  -0.184545    0.000674779  -1.70548e-8   127.256
 4520           -0.546127   -0.079556    0.0840459  -0.079556    0.000632089  -5.87656e-9   132.965
 4530           -0.723372   -0.0996693   0.0847082  -0.0996693   0.000628521  -9.42135e-9   134.774
 4540           -0.727786   -0.100129    0.0847101  -0.100129    0.000624714  -5.86455e-9   135.598
 4550           -0.988578   -0.124086    0.0838366  -0.124086    0.000616509   5.9324e-10   135.986
 4560           -0.980688   -0.123448    0.0860391  -0.123448    0.000624714  -7.53333e-9   137.726
 4570           0.414568    0.0625481    0.0845706  0.0625481    0.000577215   0.0          146.515
 4580           -0.317979   -0.0489988   0.0847555  -0.0489988   0.00053174    0.0          159.393
 4590           -0.406442   -0.0614413   0.0848829  -0.0614413   0.000511402   0.0          165.981
 4600           -0.0449544  -0.0071499   0.0841268  -0.0071499   0.000413357  -2.75596e-10  203.521
                                                                                         417 rows omitted

julia> versioninfo()
Julia Version 1.11.3
Commit d63adeda50d (2025-01-21 19:42 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 16 × AMD Ryzen 7 PRO 5850U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 16 virtual cores)

Copy link

codecov bot commented Jan 24, 2025

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 94.17%. Comparing base (f8fd782) to head (e4b203b).

Additional details and impacted files
@@            Coverage Diff             @@
##           master     #483      +/-   ##
==========================================
+ Coverage   94.11%   94.17%   +0.06%     
==========================================
  Files          14       15       +1     
  Lines        2905     2935      +30     
==========================================
+ Hits         2734     2764      +30     
  Misses        171      171              
Flag Coverage Δ
unittests 94.17% <100.00%> (+0.06%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.


Worst case accuracy is about 2e-16.
"""
function owent(h::T, a::T) where {T <: Real}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems like this implementation is specific to Float64.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed. I adapted the integration.

src/owent.jl Outdated
Comment on lines 123 to 126
towen = zero(a)
@inbounds for i in eachindex(w)
towen += w[i] * t2(h,a,x[i])
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since these are tuples, you should be able to just do:

Suggested change
towen = zero(a)
@inbounds for i in eachindex(w)
towen += w[i] * t2(h,a,x[i])
end
towen = sum(w .* t2.(h, a, x))

and it will be allocation-free and unrolled, no?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, thanks.

@stevengj
Copy link
Member

In general, I think we're hoping to move in the direction of breaking different special functions into their own packages, rather than increasing the size of this package. Maybe this should just be its own package?

@jbrea
Copy link
Author

jbrea commented Jan 27, 2025

In general, I think we're hoping to move in the direction of breaking different special functions into their own packages, rather than increasing the size of this package. Maybe this should just be its own package?

A single package just for this function? Or should it be a package devoted to Owen's table?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants