Rewrite `gamma_inc_taylor` and `gamma_inc_asym` more concisely #443

BioTurboNick · 2023-05-04T05:09:40Z

Many functions were translated directly from Fortran. This PR rewrites two of them using more concise Julia forms.

~~Ideally would like to eliminate the array allocation as well.~~ (See #442) And now does!

codecov · 2023-05-04T05:17:11Z

Codecov Report

Patch coverage: 100.00% and project coverage change: +0.01 🎉

Comparison is base (ae35d10) 94.32% compared to head (8d37a34) 94.34%.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #443      +/-   ##
==========================================
+ Coverage   94.32%   94.34%   +0.01%     
==========================================
  Files          14       14              
  Lines        2909     2881      -28     
==========================================
- Hits         2744     2718      -26     
+ Misses        165      163       -2

Flag	Coverage Δ
unittests	`94.34% <100.00%> (+0.01%)`	⬆️

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

Impacted Files	Coverage Δ
src/gamma_inc.jl	`93.49% <100.00%> (+0.03%)`	⬆️

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

src/gamma_inc.jl

BioTurboNick · 2023-05-04T12:09:17Z

Oops, left an nmc off that commit. Apologies to that person.

heltonmc · 2023-05-04T12:21:38Z

I know I originally suggested it but it would still be nice to know why it was implemented in this way to begin with as it is quite unusual and we are really just relying on the strength of the original test set here.

My last comment is just about the convergence check. My original suggestion (#442 (comment)) checked for relative tolerance where this just checks for absolute tolerance. It probably doesn't matter but if the final result is much smaller than 1 they could be different. (edit: relative check will be slower so not saying just to do it instead)

BioTurboNick · 2023-05-05T12:30:27Z

Alright, thanks! :-) I reverted to the original version, with your smaller suggestions. I kept the other one on another branch, gamma-fix2, just in case.

heltonmc · 2023-05-05T12:43:33Z

Sounds good 😊 I’ll let other people comment on what they think of the two versions.

This version will allocate though which might not solve your initial problem. If you want to sum in reverse order your probably best bet is to to use @nexprs to remove the array then can sum in any order you want. I’m not for sure it’s really needed here though but it would be nice since we are putting in the effort to make this allocation free !

Co-authored-by: David Widmann <[email protected]>

BioTurboNick · 2023-05-05T13:36:34Z

Ah, found an allocation-free way! Before I commit it:

# compute first 21 terms
ts = cumprod(ntuple(i -> x / (a + i), Val(21)))
last_t = findfirst(<(1.0e-3), ts)
last_t = last_t === nothing ? 20 : last_t - 1
next_t = last_t + 1

BioTurboNick · 2023-05-09T23:19:25Z

Ah, cumprod for tuples doesn't exist in Julia 1.3... any chance of dropping support for Julia before 1.6? (It was introduced in 1.5)

~~EDIT: Or can just define it here.~~ Nevermind, too much to bring in.

BioTurboNick · 2023-06-07T14:20:47Z

Bumping - any more work needed on this?

BioTurboNick · 2023-07-31T14:10:27Z

What do I need to do to get this in? I have a couple bits of code that would really benefit from this change, and I'd like not to have to disable precompilation.

ViralBShah · 2024-01-15T08:03:27Z

@devmotion Is this good to merge?

devmotion · 2024-01-15T13:50:18Z

Project.toml

@@ -22,7 +22,7 @@ IrrationalConstants = "0.1, 0.2"
 LogExpFunctions = "0.3.2"
 OpenLibm_jll = "0.7, 0.8"
 OpenSpecFun_jll = "0.5"
-julia = "1.3"
+julia = "1.5"


Due to the use of cumprod(::Tuple) IIRC?

devmotion · 2024-01-15T13:53:02Z

src/gamma_inc.jl

-           break
-       end
+
+    # compute first 21 terms


I wonder if summation in reverse order is necessary (something like #442 (comment) would be much simpler implementation-wise)? I wonder if there is any argument/comment regarding this choice in the original Fortran code - does anyone know where the current implementation was copied from?

In theory better for numerical stability and precision to do it this way, because otherwise you can lose the accumulated effect of the small components when they're very much smaller than the large components, but I don't know in practice if/when it affects results of this function.

I've found a couple possible sources, but I think I'm gathering that the impetus was to support arbitrarily low precision. Maybe with Float64 it works pretty much fine, which this package forces. But with lower precision, maybe the errors would add up too much if floating point math wasn't carefully managed.

Looks like @Sumegh-git and @simonbyrne were writing/reviewing the PR that introduced them, if they can comment :-)

does anyone know where the current implementation was copied from?

ACM TOMS (FREE ACCESS): https://dl.acm.org/doi/abs/10.1145/22721.23109
Fig. 2 in the paper shows the flowchart of the algorithm as a whole.

gamma_inc_taylor: Equ. 15

gamma_inc_asym: Equ. 16

The original paper is mentioned in the code comments.
Perhaps those comments should be moved to the function documentation.

SpecialFunctions.jl/src/gamma_inc.jl

Lines 831 to 834 in fdf95ab

# Reference : 'Computation of the incomplete gamma function ratios and their inverse' by Armido R DiDonato, Alfred H Morris.

# Published in Journal: ACM Transactions on Mathematical Software (TOMS)

# Volume 12 Issue 4, Dec. 1986 Pages 377-393

# doi>10.1145/22721.23109

Thanks! In that case maybe this is the FORTRAN code itself? https://jblevins.org/mirror/amiller/inc_gam.f90

BioTurboNick added 2 commits May 4, 2023 01:07

Rewrite more concisely

eb71c65

Remove old code

2363f8f

devmotion reviewed May 4, 2023

View reviewed changes

src/gamma_inc.jl Outdated Show resolved Hide resolved

src/gamma_inc.jl Outdated Show resolved Hide resolved

Remove sum (does pairwise addition?)

0b0a61a

devmotion reviewed May 4, 2023

View reviewed changes

src/gamma_inc.jl Outdated Show resolved Hide resolved

src/gamma_inc.jl Outdated Show resolved Hide resolved

heltonmc reviewed May 4, 2023

View reviewed changes

src/gamma_inc.jl Outdated Show resolved Hide resolved

devmotion reviewed May 4, 2023

View reviewed changes

src/gamma_inc.jl Outdated Show resolved Hide resolved

Remove view

46a7b39

BioTurboNick force-pushed the gamma-fix branch from ade38cb to 46a7b39 Compare May 5, 2023 12:29

Missed a couple

4e2edc6

Co-authored-by: David Widmann <[email protected]>

Move to the stack!

7ac105e

BioTurboNick added 7 commits May 10, 2023 10:02

Try to add cumprod for compat with 1.3

c095329

Add Base ref

c058a2a

Add Base ref

ee41ab7

Revert compat attempts

83afebb

Bump CI

34229ef

Change minimum supported Julia from 1.3 to 1.5

482516c

Change minimum supported Julia to 1.5

8d37a34

devmotion reviewed Jan 15, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Rewrite `gamma_inc_taylor` and `gamma_inc_asym` more concisely #443

Rewrite `gamma_inc_taylor` and `gamma_inc_asym` more concisely #443

BioTurboNick commented May 4, 2023 •

edited

Loading

codecov bot commented May 4, 2023 •

edited

Loading

BioTurboNick commented May 4, 2023 •

edited

Loading

heltonmc commented May 4, 2023 •

edited

Loading

BioTurboNick commented May 5, 2023 •

edited

Loading

heltonmc commented May 5, 2023

BioTurboNick commented May 5, 2023

BioTurboNick commented May 9, 2023 •

edited

Loading

BioTurboNick commented Jun 7, 2023

BioTurboNick commented Jul 31, 2023

ViralBShah commented Jan 15, 2024

devmotion Jan 15, 2024

BioTurboNick Jan 15, 2024

devmotion Jan 15, 2024

BioTurboNick Jan 15, 2024

BioTurboNick Feb 14, 2024 •

edited

Loading

BioTurboNick Feb 14, 2024

inkydragon Feb 19, 2024

BioTurboNick Feb 19, 2024 •

edited

Loading

	# Reference : 'Computation of the incomplete gamma function ratios and their inverse' by Armido R DiDonato, Alfred H Morris.
	# Published in Journal: ACM Transactions on Mathematical Software (TOMS)
	# Volume 12 Issue 4, Dec. 1986 Pages 377-393
	# doi>10.1145/22721.23109

Rewrite gamma_inc_taylor and gamma_inc_asym more concisely #443

Are you sure you want to change the base?

Rewrite gamma_inc_taylor and gamma_inc_asym more concisely #443

Conversation

BioTurboNick commented May 4, 2023 • edited Loading

codecov bot commented May 4, 2023 • edited Loading

Codecov Report

BioTurboNick commented May 4, 2023 • edited Loading

heltonmc commented May 4, 2023 • edited Loading

BioTurboNick commented May 5, 2023 • edited Loading

heltonmc commented May 5, 2023

BioTurboNick commented May 5, 2023

BioTurboNick commented May 9, 2023 • edited Loading

BioTurboNick commented Jun 7, 2023

BioTurboNick commented Jul 31, 2023

ViralBShah commented Jan 15, 2024

devmotion Jan 15, 2024

Choose a reason for hiding this comment

BioTurboNick Jan 15, 2024

Choose a reason for hiding this comment

devmotion Jan 15, 2024

Choose a reason for hiding this comment

BioTurboNick Jan 15, 2024

Choose a reason for hiding this comment

BioTurboNick Feb 14, 2024 • edited Loading

Choose a reason for hiding this comment

BioTurboNick Feb 14, 2024

Choose a reason for hiding this comment

inkydragon Feb 19, 2024

Choose a reason for hiding this comment

BioTurboNick Feb 19, 2024 • edited Loading

Choose a reason for hiding this comment

Rewrite `gamma_inc_taylor` and `gamma_inc_asym` more concisely #443

Rewrite `gamma_inc_taylor` and `gamma_inc_asym` more concisely #443

BioTurboNick commented May 4, 2023 •

edited

Loading

codecov bot commented May 4, 2023 •

edited

Loading

BioTurboNick commented May 4, 2023 •

edited

Loading

heltonmc commented May 4, 2023 •

edited

Loading

BioTurboNick commented May 5, 2023 •

edited

Loading

BioTurboNick commented May 9, 2023 •

edited

Loading

BioTurboNick Feb 14, 2024 •

edited

Loading

BioTurboNick Feb 19, 2024 •

edited

Loading