v0.0.1 - MWE with Laguerre and blog writeup

First commit. Be sure to read the thingsilearned.md, tests, and the README.md on how this was made!

README.md

@@ -5,3 +5,118 @@
 [![Coverage Status](https://coveralls.io/repos/miguelraz/OrthogonalPolynomials.jl/badge.svg?branch=master&service=github)](https://coveralls.io/github/miguelraz/OrthogonalPolynomials.jl?branch=master)
 
 [![codecov.io](http://codecov.io/github/miguelraz/OrthogonalPolynomials.jl/coverage.svg?branch=master)](http://codecov.io/github/miguelraz/OrthogonalPolynomials.jl?branch=master)
+
+## Introduction
+
+Hello and welcome to OrthogonalPolynomials v0.0.1!
+
+This is an open source project with a few simple goals:
+
+1. Provide an easy interface for generating arbitrary order generalized Orthogonal Polynomials and efficient evaluation and using them in Julia while being robustly tested.
+2. Show the process of producing a working software that is shared for others for academic and applied purposes, whilst building a compendium of handy software practices and principles.
+
+In essence, you rarely get to see anything but the highlights of other's work. This project is an ode to all the pulled hairs and silly mistakes and drowning frustration to getting. this. darn. code. to run properly.
+
+This goal will be met by [video tutorials](https://www.youtube.com/channel/UC840v4b_71e78fmPHiCPQVg) and [livestreaming](https://www.twitch.tv/brainrpg) where you can join in on the action if you so desire, and observe the entirety of the software design process.
+
+## Code and Roadmap
+
+This is how we well try to meet goal 1 above, in about ~100 lines of Julia code.
+We will use some:
+
+1. Metaproramming in Julia á la [@horner's macro](https://youtu.be/dK3zRXhrFZY) and some generated functions
+
+2. A very convenient formula that we found in [Abramowitz and Stegun](http://people.math.sfu.ca/~cbm/aands/page_789.htm) that lends itself to multiple dispatch and generating the coefficients we want for our `@horner` formula at compile time.
+
+3. Multiple dispatch design.
+
+In short, the user dream is this:
+
+```
+julia> Pkg.add("OrthogonalPolynomials");
+
+julia> using OrthogonalPolynomials, BenchmarkTools;
+
+julia> L₁₂(x) = a(Laguerre(12));
+
+julia> @btime L₁₂(.5)
+  1.449 ns (0 allocations: 0 bytes)
+-0.23164963886602852
+```
+
+## Motivation and related work
+
+Orthogonal Polynomials are polynomials that appear almost everywhere in applied mathematics and sciences due to their centuries-old connections established to differential equations, group theory, numerical analysis, etc.
+
+People find them very interesting for many different reasons, and you can find them in some well established software like `MATLAB` or `Mathematica` ( if you want to pay for it) or in Julialand through packages like [ApproxFun.jl](https://github.com/JuliaApproximation/ApproxFun.jl/tree/master/src/Spaces). If you're interested in more thorough manipulation of Orthogonal Polynomial spaces, they are probably your best bet right now.
+
+So why another package?
+
+First, I wanted to build my own package, and do not mind upstreaming it later on into SpecialFunctions.jl or ApproxFun.jl, if it so happens to be found a happy home.
+
+Second of all, I wanted to learn by myself and have a lightweight standalone no-dependency repository where I could simply point to others if they wanted a package that does one thing, and does it very well. Most importantly, no other packages, as known to the author at time of writing, offer the full optimizations possible for user chosen parameters, e.g., offering Legendre 17th order with α = 3, β = 9.
+
+Third, I wanted to compare my approach to established literature and other working code and see how I fared, all of this in parallel to producing a living document of videos and notes so that others may learn from my mistakes. As a comparison, here are some benchmarks for the 12th order Laguerre polynomial with different implementations:
+
+```
+# Hand-Crafted implementation of L₁₂ - writing out the entire polynomial
+julia> @btime L12(.5)
+  478.579 ns (0 allocations: 0 bytes)
+-0.23164963886551918
+
+# Approx Fun implementation
+julia> @btime h12(.5)
+  416.910 ns (4 allocations: 64 bytes)
+-0.2316496388655192
+
+# Naive recurrence formula
+julia> @btime laguerre(12,.5)
+  430.729 ns (52 allocations: 832 bytes)
+-0.23164963886551906
+
+# Recurrence formula guarding the 0th case
+julia> @btime laguerre_corrected(12,.5)
+  442.377 ns (52 allocations: 832 bytes)
+-0.23164963886551906
+
+# Written Horner method case for degree 12
+julia> @btime goal12(.5)
+  1.449 ns (0 allocations: 0 bytes)
+-0.23164963886602852
+```
+
+The key takeaway is that we can generate efficient code (since all our polynomials' coefficients can be computed at compile time) and thus cut our computation time by almost 2 orders of magnitude. Note as well that because we are reducing dramatically the number of floating point operations, we are losing less precision than other methods.
+It is also worth noting that the Julia JIT compiler performs quite well for simply writing the most simplistic code possible, as in the `laguerre` code. (Credits to `@yingboma`, who contributed that function in slack.) It doesn't scale into the `@horner` land of ridiculous speed due to not being spoonfed the appropriate structure to SIMD the generated code.
+
+## Approach
+
+This is how we will solve the problem as was hinted above in the roadmap.
+There is a couple of things you need to learn about first.
+
+1. How the `@horner` macro works, explained in the video link above.
+
+2. How the formula found matches this approach, and how we our code will build up incrementally towards it.
+
+3. Why using `@generated` functions is necessary at all.
+
+The `@horner` macro works by taking a given list of coefficients to compute a polynomials and rearranging it algebraically so that it a) faster b) more precise, a classical result found in Knuth, TAOCP II. We will have to modify the macro just a tinybit so that it fits with our master formula.
+
+The master formula was found by ~procrastination~ leafing through the book. After reading a discourse post by [Steven G Johnson](https://discourse.julialang.org/t/simple-metaprogramming-exercises-challenges/731/6) and looking at the formula, it is a clear match if you just rearrange it a bit. To wit:
+
+$$ a\_{m-1}(x) = 1 - \frac{b\_m}{c\_m}f(x) a\_m(x)$$
+
+Fits just right in the `muladd` patter if we rearrange the `1` to be the addition. This is why I found the idea of this package exciting: we can take classical orthogonal polynomials and give them a reach and speed not yet taken care of in Julia.
+
+From here, we want to be able to produce functions for the coefficients that we want for a given polynomial. This means that, based on the polynomial we want, we will have to have different behavior for each coefficient, depending on the type of polynomial (this screams multiple dispatch! Do you see how?)
+
+Eventually, we will want to produce these coefficients at compile time so that a new function can optimized by the JIT compiler when we actually call it. This will require the use of `@generated` functions, which can only see types as arguments. This means that any and all computation we want to do in our generated functions like cooking up the desired polynomial coefficients will have to be crammed into the types and the functions have to work on the types alone.
+
+Lastly, The codebase in this package is very much alpha (development stage). It is meant to showcase the basic building blocks and should not be used in production.
+
+Our starting point is to have a Minimal Working Example (MWE) of the vanilla Laguerre polynomial - even if it is a bit wonky.
+
+First Lesson: build a basic example where your code runs, then build on top of it.
+
+The code is all stashed in `src/OrthogonalPolynomials.jl` so go and have a look at the basic machinery before moving on to `v0.0.2`.
+
+Be sure to look at the tests in `tests/runtests.jl` where we verified our workflow along the way.

docs/.gitignore

@@ -0,0 +1,2 @@
+build/
+site/

docs/make.jl

@@ -0,0 +1,13 @@
+using Documenter
+using OrthogonalPolynomials
+
+makedocs(
+    modules = [OrthogonalPolynomials]
+)
+
+# Documenter can also automatically deploy documentation to gh-pages.
+# See "Hosting Documentation" and deploydocs() in the Documenter manual
+# for more information.
+#=deploydocs(
+    repo = "<repository url>"
+)=#

docs/mkdocs.yml

@@ -0,0 +1,27 @@
+# See the mkdocs user guide for more information on these settings.
+#   http://www.mkdocs.org/user-guide/configuration/
+
+site_name:        OrthogonalPolynomials.jl
+#repo_url:         https://github.com/USER_NAME/PACKAGE_NAME.jl
+#site_description: Description...
+#site_author:      USER_NAME
+
+theme: readthedocs
+
+extra_css:
+  - assets/Documenter.css
+
+extra_javascript:
+  - https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML
+  - assets/mathjaxhelper.js
+
+markdown_extensions:
+  - extra
+  - tables
+  - fenced_code
+  - mdx_math
+
+docs_dir: 'build'
+
+pages:
+  - Home: index.md

docs/src/index.md

@@ -0,0 +1,3 @@
+# OrthogonalPolynomials.jl
+
+Documentation for OrthogonalPolynomials.jl

src/OrthogonalPolynomials.jl

@@ -1,5 +1,25 @@
+__precompile__()
+
 module OrthogonalPolynomials
 
-# package code goes here
+# 1. First design - Laguerre only
+d(n, α=0) = binomial(n+α, n)
+b(n, m) = n-m+1
+c(m, α) = m*(α+m)
+k(n, α=0) = [-b(n,i)*inv(c(i, α)) * d(n, α) for i in 1:n]
+f(x) = x
+
+a(x, n, ks = k(n,α), i=0) = i == n ? :(1) : return :(muladd( $(ks[i+1]*f(x)), $((a)(x, n, ks, i+1)) , 1))
+
+# 2. Macro function as a for loop
+macro a(x,n,ks=k(n,α))
+    ex = 1
+    for i in n-1 : -1 : 0
+        ex = :(muladd( $(ks[i+1]*f(x)), $ex, 1))
+    end
+    return :($ex)
+end
+
+export d,b,c,k,f,a
 
 end # module

test/runtests.jl

@@ -1,5 +1,103 @@
 using OrthogonalPolynomials
-using Base.Test
+using Base.Test, Base.Math
 
 # write your own tests here
-@test 1 == 2
+
+# Hand Implementation
+L0(x) = 1.0
+L1(x) = -x + 1.0
+L2(x) = .5 * (x^2 - 4x + 2)
+L3(x) = (1/6)*(-x^3 +9x^2 -18x + 6)
+L4(x) = (1/24)*(x^4 - 16x^3 + 72x^2 - 96x + 24)
+L5(x) = (1/120)*(-x^5 + 25x^4 -200x^3 + 600x^2 - 600x + 120)
+L6(x) = (1/720)*(x^6 - 36x^5 + 450x^4 - 2400x^3 + 5400x^2 - 4320x + 720)
+L7(x) = (1/5040)*(-x^7 + 49x^6 - 882x^5 + 7350x^4 - 29_400x^3 + 52_920x^2 - 35_280x + 5040)
+L8(x) = (1/40320)*(x^8 - 64x^7 + 1568x^6 - 18_816x^5 + 117_600x^4 - 376_320x^3 + 564_480x^2 - 322_560x + 40320)
+L9(x) = (1/362880)*(-x^9 + 81x^8 - 2592x^7 + 42_336x^6 - 381_024x^5 + 1_905_120x^4 - 5_080_320x^3 + 6_531_840x^2 - 3_265_920x + 362_880)
+L10(x) = (1/3628800)*(x^10 - 100x^9 + 4_050x^8 - 86_400x^7 + 1_058_400x^6 - 7_620_480x^5 + 31_752_000x^4 - 72_576_000x^3  + 81_648_000x^2 - 36_288_000x + 3_628_800)
+L11(x) = (1/39916800)*(-x^11 + 121x^10 - 6_050x^9 + 163_350x^8 - 2_613_600x^7 + 25_613_280x^6 - 153_679_680x^5 + 548_856_000x^4 - 1_097_712_000x^3 + 1_097_712_000x^2 - 439_084_800x + 39_916_800)
+L12(x) = (1/479001600)*(x^12 - 144x^11 + 8_712x^10 - 290_400x^9 + 5_880_600x^8 - 75_271_680x^7 + 614_718_720x^6 - 3_161_410_560x^5 + 9_879_408_000x^4 - 17_563_392_000x^3 + 15_807_052_800x^2 - 5_748_019_200x + 479_001_600)
+
+Ls = [L0, L1, L2, L3, L4, L5, L6, L7, L8, L9, L10, L11, L12]
+
+# Naive implemenatation
+"Credit to @yingboma"
+function laguerre(n, x)
+
+    l0, l1 = 1, 1-x
+    for k in 1:n-1
+        l0, l1 = l1, ( (2k+1-x)*l1 - k*l0 )/(k+1)
+    end
+    l1
+end
+
+# Correct 0 implemenatation
+"Credit to @yingboma"
+function laguerre_corrected(n, x)
+
+    n == 0 ? (return 1) : nothing
+
+    l0, l1 = 1, 1-x
+    for k in 1:n-1
+        l0, l1 = l1, ( (2k+1-x)*l1 - k*l0 )/(k+1)
+    end
+    l1
+end
+
+goal12(x) = @evalpoly x 1.0 -12.0 33.0 -36.666666666666664 20.625 -6.6 1.2833333333333334 -0.15714285714285714 0.012276785714285714 -0.0006062610229276896 1.8187830687830687e-5 -3.0062530062530064e-7
+
+
+# Abramowitz and Stegun pg. 799
+LaguerreTable = [
+    1 1 1 2 6 24 120 720 5040 40320 362880 3628800 39916800 479001600;
+    1 1 1 -4 -18 -96 -600 -4320 -35280 -322560 -3265920 -36288000 -439084800 -6748019200;
+    2 2 -4 2 18 144 1200 10800 105840 1128960 13063680 163296000 2195424000 31614106600;
+    6 6 -18 9 6 -96 -1200 -14400 -176400 -2257920 -30481920 -435456000 -6586272000 -105380352000;
+    24 24 -96 72 -16 24 600 10800 176400 2822400 45722880 762048000 13172544000 237105792000;
+    120 120 -600 600 -200 25 120 -4320 -106840 -2257920 -45722880 -914457600 -18441561600 -379369267200;
+    720 720 -4320 5400 -2400 450 -36 720 35280 1128960 30481920 762048000 18441561600 442597478400;
+    5040 5040 -35280 52920 -29400 7350 -882 49 5040 -322560 -13063680 -435456000 -13172544000 -379369267200;
+    40320 40320 -322560 564480 -376320 117600 -18816 1568 -64 40320 3269520 163296000 6586272000 237105792000;
+    362880 362880 -3265920 6531840 -5080320 1905120 -381024 42336 -2592 81 362880 -36288000 -2195424000 -105380352000;
+    3623800 3628800 -36288000 81648000 -72576000 31752000 -7620480 1058400 -86400 4050 -100 3628800 439084800 31614106600;
+    39916800 39916800 -439084800 1097712000 -1097712000 648806000 -103679680 25613280 -2613600 163350 -6050 121 39916800 -5748019200;
+    479001600 479001600 -5748019200 15807052800 -17563392000 9879408000 -3161410560 614718720 -75271680 5880600 -290400 8712 -144 479001600]
+
+#n\z is the first column
+# Abramowitz and Stegun p. 800
+testvals = [0. .5 1 3 5 10];
+TestLaguerreTable = [
+    0 1.0000000000 1.0000000000 1.0000000000 1.0000000000 1.0000000000;
+    1 0.5000000000 0.0000000000 -2.0000000000 -4.0000000000 -9.0000000000;
+    2 0.1250000000 -0.5000000000 -0.5000000000 3.5000000000 31.0000000000;
+    3 -0.1458333333 -0.6666666667 1.0000000000 2.6666666667 -45.6666666667;
+    4 -0.3307291667 -0.6250000000 1.3750000000 -1.2916666667 11.0000000000;
+    5 -0.4455729167 -0.4666666667 0.8500000000 -3.1666666667 34.3333333333;
+    6 -0.5041449653 -0.2569444444 -0.0125000000 -2.0902777778 -3.4444444444;
+    7 -0.5183392237 -0.0404761905 -0.7464285714 0.3253968254 -30.9047619048;
+    8 -0.4983629984 0.1539930556 -1.1087053571 2.2357390873 -16.3015873016;
+    9 -0.4529195204 0.3097442681 -1.0611607143 2.6917438272 14.7918871252;
+    10 -0.3893744141 0.4189459325 -0.7000223214 1.7562761795 27.9841269841;
+    11 -0.3139072988 0.4801341791 -0.1807995130 0.1075436909 14.5369568703;
+    12 -0.2316496389 0.4962122235 0.3403546063 -1.4486042948 -9.9037464593];
+
+@testset "Base Table" begin
+    @testset "Hand Implemenation" begin
+        @test all(f(0) ≈ 1.0 for f in Ls)
+        @test all(Ls[i](testvals[j]) ≈ TestLaguerreTable[i,j] for i in 1:13, j in 2:6)
+    end
+    @testset "Naive Implementation" begin
+        @test all(laguerre(i,0.0) ≈ 1.0 for i in 0:12)
+        @test_broken all([laguerre(i-1,testvals[j]) ≈ TestLaguerreTable[i,j] for i in 1:13, j in 2:6])
+    end
+    @testset "Naive Corrected" begin
+        @test all(laguerre_corrected(i,0.0) ≈ 1.0 for i in 0:12)
+        @test all([laguerre_corrected(i-1, testvals[j]) ≈ TestLaguerreTable[i,j] for i in 1:13, j in 2:6])
+    end
+    @testset "Stubborn Way" begin
+        @test all(eval.(a(0,i)) .≈ 1 for i in 1:12)
+        @test all(eval.([a(testvals[j],i) for i in 0:12, j in 2:6]) .≈ TestLaguerreTable[:,2:6])
+    end
+end
+
+goal12(x) = @horner x 1.0 -12.0 33.0 -36.666666666666664 20.625 -6.6 1.2833333333333334 -0.15714285714285714 0.012276785714285714 -0.0006062610229276896 1.8187830687830687e-5 -3.0062530062530064e-7

thingsilearned.md

@@ -0,0 +1,17 @@
+1. Testing for all functions in an array:
+@test all(f(theta)==1 for f in fs)
+Credit to @chrisrackauckas
+
+2. You an also `0 .|> fs` for "reverse broadcasting" a value to an array of functions.
+Credit to @MichaelBorregaard
+
+3. Generated functions ONLY SEE TYPE INFO. You gotta use dispatch for all the goodies.
+
+4. SIMD needs a for loop - that will make it preferable to use the for loop approach over the recursive one!
+
+5. Ternary operator requires a `()` around a return statemnt for it to work.
+`x == zero(x) ? (return 1) : nothing
+
+6. Use `@test_broken to know that your tests are broken, and you won't forget later!`
+
+7. 

utils.jl

@@ -0,0 +1,35 @@
+
+
+#Laguerre breaks at 16
+for i in 12:20
+    println("n = $i ", all(L_rec(i) .< 1/eps()))
+end
+
+
+
+# Laguere Coefficient Calculator
+# Credit to [Peter Luschny Apr. 11, 2015](https://oeis.org/search?q=laguerre+coefficient+triangle&sort=&language=&go=Search)
+function r(n,m)
+    if 0 == n == m
+        1
+    elseif m == -1
+        0
+    elseif n < m
+        0
+    elseif n >= 1
+        (n+m+1) * r(n-1, m) - r(n-1, m-1)
+    end
+end
+
+L_rec(n) = [r(j,i) for i in 0:n,j in 0:n]
+
+# Credit to [Paul Barry, Aug 28, 2005.](https://oeis.org/search?q=coefficient+triangle+hermite&sort=&language=&go=Search)
+# Hermite Coefficient Calculator
+function T(n,k)
+
+    if iseven(n-k) && (n-k >= 0)
+        ( (-1)^((n-k)/2) * (2^k) * factorial(n)) / (factorial(k) * factorial((n-k)/2))
+    else
+        0
+    end
+end