1_One/4_generic_programming.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generic programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From [Wikipedia](https://en.wikipedia.org/wiki/Generic_programming):\n",
    "> **Generic programming** is a style of computer programming in which algorithms are written in terms of types *to-be-specified-later* that are then *instantiated* when needed for specific types provided as parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Vandermonde matrix (revisited)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Vandermonde matrix:](https://en.wikipedia.org/wiki/Vandermonde_matrix)\n",
    "\\begin{align}V=\\begin{bmatrix}1&\\alpha _{1}&\\alpha _{1}^{2}&\\dots &\\alpha _{1}^{n-1}\\\\1&\\alpha _{2}&\\alpha _{2}^{2}&\\dots &\\alpha _{2}^{n-1}\\\\1&\\alpha _{3}&\\alpha _{3}^{2}&\\dots &\\alpha _{3}^{n-1}\\\\\\vdots &\\vdots &\\vdots &\\ddots &\\vdots \\\\1&\\alpha _{m}&\\alpha _{m}^{2}&\\dots &\\alpha _{m}^{n-1}\\end{bmatrix}\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Keep function argument types generic if possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function vander_naive(x::Vector)\n",
    "    m = length(x)\n",
    "    V = Matrix{Float64}(undef, m, m)\n",
    "    for j = 1:m\n",
    "        V[j,1] = 1.0\n",
    "    end\n",
    "    for i= 2:m\n",
    "        for j = 1:m\n",
    "            V[j,i] = x[j] * V[j,i-1]\n",
    "            end\n",
    "        end\n",
    "    return V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rand(3)\n",
    "vander_naive(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_naive(1:3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The annotation `::Vector` in the the function signature is unnecessarily specific."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function vander_naive(x::AbstractVector)\n",
    "    m = length(x)\n",
    "    V = Matrix{Float64}(undef, m, m)\n",
    "    for j = 1:m\n",
    "        V[j,1] = 1.0\n",
    "    end\n",
    "    for i= 2:m\n",
    "        for j = 1:m\n",
    "            V[j,i] = x[j] * V[j,i-1]\n",
    "            end\n",
    "        end\n",
    "    return V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_naive(1:3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Avoid explicit typing if possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_naive([1,2,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is the result a matrix of floating point numbers....?\n",
    "\n",
    "Even worse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_naive(rand(ComplexF64, 3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily cover those cases as well by only slightly modifying our code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function vander_almost_generic(x::AbstractVector)\n",
    "    T = eltype(x)\n",
    "    m = length(x)\n",
    "    V = Matrix{T}(undef, m, m)\n",
    "    for j = 1:m\n",
    "        V[j,1] = 1.0\n",
    "    end\n",
    "    for i= 2:m\n",
    "        for j = 1:m\n",
    "            V[j,i] = x[j] * V[j,i-1]\n",
    "            end\n",
    "        end\n",
    "    return V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_almost_generic([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_almost_generic(rand(ComplexF64, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_almost_generic([\"Stadt\", \"Land\", \"Fluss\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function vander_generic(x::AbstractVector{T}) where T # this is the same as just x::AbstractVector\n",
    "    m = length(x)\n",
    "    V = Matrix{T}(undef, m, m)\n",
    "    for j = 1:m\n",
    "        V[j,1] = one(T)\n",
    "    end\n",
    "    for i= 2:m\n",
    "        for j = 1:m\n",
    "            V[j,i] = x[j] * V[j,i-1]\n",
    "            end\n",
    "        end\n",
    "    return V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_generic([\"Stadt\", \"Land\", \"Fluss\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One more level of generality, just because we can. :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_generic([3, \"Stadt\", 4 + 5im])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function vander_supergeneric(x::AbstractVector{T}) where T\n",
    "    m = length(x)\n",
    "    V = Matrix{T}(undef, m, m)\n",
    "    for j = 1:m\n",
    "        V[j,1] = one(x[j])\n",
    "    end\n",
    "    for i= 2:m\n",
    "        for j = 1:m\n",
    "            V[j,i] = x[j] * V[j,i-1]\n",
    "            end\n",
    "        end\n",
    "    return V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vander_supergeneric([3, \"Stadt\", 4 + 5im])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### And all of this comes at no performance penality"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using BenchmarkTools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rand(Float64, 100);\n",
    "@btime vander_naive($x);\n",
    "@btime vander_supergeneric($x);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, for this specific example **our generic code is faster** in a few cases inasmuch as type conversions are unnecessary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rand(Int, 100);\n",
    "@btime vander_naive($x);\n",
    "@btime vander_supergeneric($x);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rand(Bool, 100);\n",
    "@btime vander_naive($x);\n",
    "@btime vander_supergeneric($x);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, sometimes it is worth converting to a different type to dispatch to a faster method or to utilize magic like compiler optimizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = rand(Float32, 100);\n",
    "@btime vander_naive($x);\n",
    "@btime vander_supergeneric($x);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arbitrary precision computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say you have implemented the following crazily complex physics code as part of your thesis project which takes in a number and spits out the answer to life, the universe, and everything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function answer_to_life_universe_and_everything(x::Integer)\n",
    "    m = sin((2*x)^100)\n",
    "    c = 42\n",
    "    E = m*c^2\n",
    "    a = sqrt(abs(E))\n",
    "    b = atan(m)\n",
    "    c = a^2 + b^2\n",
    "    answer = sqrt(1764)/(1+exp(-c))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "answer_to_life_universe_and_everything(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alright, apparently the answer is 21!\n",
    "\n",
    "The author: \"I checked the code multiple times, it is correct. So let's publish.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Without changing a line of code we can check the *correctness* (in the numerical sense) of our result using [arbitrary precision arithmetics](https://docs.julialang.org/en/latest/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic-1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(big(2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "answer_to_life_universe_and_everything(big(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's get a bit more fancy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Interact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@manipulate for n in 1:20\n",
    "    [i*j for i in 1:n, j in 1:n]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function insert_block(A::AbstractMatrix, i, j, what=7)\n",
    "    B = copy(A)\n",
    "    B[i:i+2, j:j+2] .= what\n",
    "    B\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = fill(0, 9, 9)\n",
    "insert_block(A, 3, 5)  # this returns the new matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = fill(0, 10, 10)\n",
    "n = size(A, 1)\n",
    "\n",
    "@manipulate for i in 1:n-2, j in 1:n-2\n",
    "    insert_block(A, i, j)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's add some color!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our function `insert_block` is generic. Since the first argument `A isa AbstractArray`, we can index into it and set new values. Pretty much every value type is fine!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@manipulate for n in 1:80\n",
    "    distinguishable_colors(n)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = distinguishable_colors(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = fill(colors[1], 10, 10)\n",
    "n = size(A, 1)\n",
    "\n",
    "@manipulate for i in 1:n-2, j in 1:n-2\n",
    "    insert_block(A, i, j, colors[4])\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generic Programming + Multiple Dispatch + JIT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The possibility to write generic algorithms that compile to fast machine code in combination with multiple dispatch leads to an ([unreasonable](https://pretalx.com/juliacon2019/talk/BCYWZJ/)) amount of code reuse. This sharing of code comes in two forms:\n",
    "1. **Sharing types** among a wide variety of packages implementing different algorithms;\n",
    "2. **Sharing generic algorithms** that work for different package-defined types implementing common abstractions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <img src=\"revdeps.png\" alt=\"drawing\" width=\"800\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. **Sharing types:** DataStructures.jl, OrderedCollections.jl, StaticArrays.jl, Colors.jl, Measurements.jl ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As of the time of this writing, **804 packages** depend on the data types provided in [DataStructures.jl](https://juliacollections.github.io/DataStructures.jl/latest/).\n",
    "\n",
    "**851 packages** reuse type implementations in [OrderedCollections.jl](https://github.com/JuliaCollections/OrderedCollections.jl).\n",
    "\n",
    "**That's about every third package!**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. **Sharing generic algorithms:** StatsBase.jl, SortingAlgorithms.jl, GenericLinearAlgebra.jl, ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Emergent features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Measurement` type from [Measurements.jl]() and differential equation solver from [OrdinaryDiffEq.jl](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) (i.e. [DifferentialEquations.jl](https://github.com/JuliaDiffEq/DifferentialEquations.jl))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using OrdinaryDiffEq, Measurements, PyPlot\n",
    "\n",
    "#Half-life of Carbon-14 is 5730 years.\n",
    "c = 5.730 ± 2\n",
    "\n",
    "#Setup\n",
    "u0 = 1.0 ± 0.1\n",
    "tspan = (0.0 ± 0.0, 1.0 ± 0.0)\n",
    "\n",
    "#Define the problem\n",
    "radioactivedecay(u,p,t) = -c*u\n",
    "\n",
    "#Pass to solver\n",
    "prob = ODEProblem(radioactivedecay,u0,tspan)\n",
    "sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# analytic solution\n",
    "u = u0 .* exp.(-c .* sol.t);\n",
    "\n",
    "# plot solution\n",
    "ts = getfield.(sol.t, :val)\n",
    "solvals = getfield.(sol, :val)\n",
    "solerrs = getfield.(sol, :err);\n",
    "\n",
    "errorbar(ts, solvals, yerr=solerrs)\n",
    "plot(ts, getfield.(u, :val), color=\"red\", lw=2)\n",
    "ylabel(\"u(t)\")\n",
    "xlabel(\"t\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that, in some sense, **Julia implemented that feature by itself**.\n",
    "\n",
    "The authors of Measurements.jl and DifferentialEquations.jl never had any collabration on this.\n",
    "\n",
    "It **just works**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Core messages of this Notebook\n",
    "\n",
    "* It is simple to write **type-generic code** in Julia and you should do it.\n",
    "* Generally, **generic code is just as fast as specific code**.\n",
    "* Generic Programming + Multiple Dispatch + JIT = **lots of code sharing and emergent features**"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}