Julia Season of Contributions (JSoC '19 @ SpecialFunctions)

1. The Pledge

I am very happy to announce that this summer, I'll be working on developing the SpecialFunctions.jl package under The Julia Language organization.

My project is based on implementing new Special Functions which are lacking in SpecialFunctions.jl package which are of utmost importance to any researcher or student working in Statistical, physics, or function analysis domain. Currently, the package consists of primary functions like the error functions, bessel and airy functions, etc. I propose to add to this list the Incompete Gamma function, Incomplete Beta function, both of their inverse ones and also implement a Non-central Beta cumulative distribution function all in native Julia. All of these initially will also have a BigFloat implementation from the GNU MPFR library. But I'll be in constant search to find ways to implement them natively in Julia instead of calling them via wrappers.

2. Starting off with the work (The Turn)

The first day started with planning out the work for the week, which included wrapping up the work previously done by me during the proposal review period on Incomplete Gamma Function.

P(a,x) = 1/\Gamma (a) \int_{0}^{x} e^{-t}t^{a-1} dt

Mainly, the work that was left was to refactor the code to remove @goto's and convert those either into separate functions or control-statements. After some reviews made by my project mentor Simon Byrne (on my PR #146 ) I started working on those and successfully was able to remove a few big @goto blocks. Also I did some of the documentation for the newly created function definitions and provided DLMF links for the same as reference. The code consists of a fairly complex if-else structure due to which it will take a bit more time to finish this off, although I target to finish it within a day or two. Also I was reading through the paper to check for any flaws in implementation that I might have left. Testing was done to quite some extent before. Manually added one/two tests more to cover more cases. Now will work on evaluating accuracy plots in the domain the next day and continue with the refactoring.

After getting a taste of the immediate work to be completed, next day was very productive in the sense, that I completed refactoring the code. Then it was followed with addressing the reviews/comments made by my mentor on the PR to improve code structure further and clean up things a bit.

3. Picking up pace

Now with significant work being done on the Incomplete Gamma PR (as mentioned above) it stood at a good position to be merged. Finally, with a few more addition of extensive tests and cleaning up the code PR #146 finally got merged in SpecialFunctions.jl. With this the first objective of my proposal was complete well within the estimated time of one week.

Then immediately I moved on to implementing the Inverse Incomplete Gamma now.

Given a and :

finding the value of x.

Within two-three days I was done with writing a very crude version of the code which although wasn't tested but seemed to be working fine. It was based on this paper. But soon while doing extensive testing, I discovered a bug in the code since it wasn't converging for a < 0.1. Although I tried debugging the error but all in vain. To me it seemed a logical error somewhere and was difficult to find in the entire code. So without wasting much time I started implementing the function based on another famous paper. Since, I had gathered a bit more of Julia experience under my sleeves now so implementing it didn't take much time.

using BenchmarkTools
using SpecialFunctions
@btime gamma_inc(.5,.6,0)

5.3s

Julia

(0.726678, 0.273322)

#Abs. rel. error with Scipy.

13.5s

Julia

4. Well with the flow

Finally, this implementation of the Inverse Incomplete Gamma was working just fine. Finally I sent the PR, and waited for reviews by my mentor upon improving the code structure. This was already implemented using if-else structure instead of @goto and hence didn't take much time to clean up the code. Testing was also done quite fast simultaneously following up the reviews made by my mentor. Then documentation and docstrings alongwith equations were also added. Finally, PR #164 was ready to be merged too and all that in a week.

As it can be seen the maximum relative error for both the complementary functions are of order and which is good news.

gamma_inc_inv(0.5,.5,.5)

1.3s

Julia

0.227468

#Abs. rel. error plot with SciPy.

1.2s

Julia

Now, next comes up the Incomplete Beta Function the most challenging of all with lots of corner cases and evaluation domains to be dealt with which implies more bugs to fix too ;p

5. A bit relief

Finally I started off coding for the Incomplete Beta Function. I had to be careful the entire time while giving variable names, function names and had to care for every trivia because one mistake and I lose track of where I was ;p

The Incomplete Beta function is given by:

The implementation is based on this paper.

I spent the first few days understanding the paper, re-reading it 2-3 times and planning out how would I go about doing this. Then without wasting much time I started coding case-by-case. The code was full of @goto statements in the initial draft I wrote. Also there were a few bugs and few cases uncovered/missed. I carefully noted those down all at once and searched for the best way to integrate the missing cases into the already existing mess (code) ;p

I was lucky enough after this because all the tests passed in my 2nd attempt of correcting the code and making it clean. Then I submitted a PR for this and waited for reviews. Initially, most of the code reviews were on refactoring the code, i.e. removing the annoying @goto and writing them as if-else. Also we began modularizing our code by splitting every use-case into separate functions. After 2 weeks of this same process being repeated finally the PR was looking good enough to be merged. I added the documentation and LaTeX equations and the PR was good to go. The absolute relative error was bounded by as shown in the plot below. The function returns a tuple (x, 1.0-x) for the upper and lower incomplete beta functions.

using BenchmarkTools
using SpecialFunctions
@btime beta_inc(.5,.5,.01)

6.6s

Julia

(0.0637686, 0.936231)

#Relative error plot compared with SciPY

0.9s

Julia

6. Good to go

After completing the Incomplete Beta function, similar to it's Gamma counterpart, I started off with the Inverse Incomplete Beta. This felt relatively easier to deal with since I had by now gained more of Julia experience under my sleeves as evident from the very little code review required in this PR.

Given a, b and , find 'x'.

The implementation was based on this paper here. It involved choosing a proper initial estimate and then performing a modified version of the Newton-Raphson iteration to solve for 'x'. This took significantly less amount of time to finish. Although the initial paper I had implemented had some flaws so I had to rewrite the code at some places to remove all the bugs and make it look good enough to be merged. Also the abs. rel. error didn't exceed much over .

@btime beta_inc_inv(2.0,2.0,.216)

4.8s

Julia

(0.3, 0.7)

#Absolute relative error with SciPy

1.4s

Julia

7. Final one (The Prestige)

Now was the time for the last part of the project - CDF of Non-central Beta Function . (Enough with central ones ;p)

where q was the Poisson variable.

The first 2-3 days involved mainly reading papers and other relevant literature to find up some good implementation with better error bounds and overflow handling. Most of the papers stressed on one single thing - using backward recurrence instead of forward. So this was the paper that I found most suitable to implement since it started from and then performed both forward and backward recurrence to arrive at the result.

Bonus: Using this I could also implement the non-central F distribution , since it directly depends on the non-central beta given by:

Overall the entire work was finished within a week with some more code review and structure changes the PR was all set to be merged.

8. Conclusion

With this final PR being merged I justified my proposal well which I had submitted at the beginning of the program to serve as my rough timeline. All the works I did were completed well within the deadline without haste. This was a great experience for me. I would thank every member of the Julia community since they were always ready to help and were just one ping away on Slack. And not to forget my mentor Simon Byrne who helped me a lot giving advices and reviewing my not-so-readable code, pointing out silly mistakes and what not! It wouldn't have been possible without you.

I really look forward to what Julia achieves both as a computing langauge and as a community. Would surely keep in touch with everyone!

THANK YOU!