06_Composability_Code_Reuse.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Composable software from generic code\n",
    "\n",
    "In the previous notebook we produced the following sum function, which was basically en par in speed with C or numpy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function fastsum(v)\n",
    "    result = zero(eltype(v))\n",
    "    @simd for i in 1:length(v)    # @simd enforces vectorisation in the loop\n",
    "        @inbounds result += v[i]  # @inbounds suppresses bound checks\n",
    "    end\n",
    "    result\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But unlike a C implementation, we are not at all restricted to using a particular data type ... and this let's us do crazy things **even though the code is equally fast**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Composability\n",
    "\n",
    "... or the real power of Julia.\n",
    "\n",
    "\n",
    "(a) **Elevated precision** ... let's consider a nasty case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A vector that sums to one by construction\n",
    "function generate(N)\n",
    "    x = randn(N) .* exp.(10 .* randn(N))\n",
    "    x = [x; -x; 1.0]\n",
    "    x[sortperm(rand(length(x)))]\n",
    "end\n",
    "        \n",
    "v = generate(10000);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastsum(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... looks wrong, let's try elevated precision:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fastsum(big.(v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(b) **Tracking numerical error**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import IntervalArithmetic\n",
    "\n",
    "fastsum(IntervalArithmetic.Interval.(v))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... ups clearly something wrong here! But now we know."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(c) **Error propagation**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Measurements\n",
    "v = [1 ± 1e-10, 1.1 ± 1e-12, 1.2 ± 1e-10]\n",
    "fastsum(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(d) **Unlocking new features**\n",
    "\n",
    "Let us solve:\n",
    "\n",
    "$$\\frac{du(t)}{dt} = -cu(t)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using OrdinaryDiffEq\n",
    "\n",
    "# Half-life of Carbon-14 is 5730 years.\n",
    "c = 5.730\n",
    "\n",
    "# Setup\n",
    "u0 = 1.0\n",
    "tspan = (0.0, 1.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": [
    "using Plots\n",
    "\n",
    "plot(sol.t, sol.u, ylabel=\"u(t)\", xlabel=\"t\", lw=2, legend=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now what if we have measurement errors ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using OrdinaryDiffEq, Measurements, Plots\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, 1.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);\n",
    "\n",
    "plot(sol.t, sol.u, ylabel=\"u(t)\", xlabel=\"t\", lw=2, legend=false)"
   ]
  },
  {
   "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 collaboration on this.\n",
    "\n",
    "It **just works**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generic programming + multiple dispatch\n",
    "\n",
    "So how does this work?\n",
    "\n",
    "Julia's ecosystem is centred around convenient abstractions *by convention* (**duck typing**).\n",
    "\n",
    "- As a **user of a type** we can easily pick within the type hierarchy the most generic type, which satisfies all our needs. If we do not need the precise details of the `Array` implementation, we indicate our **algorithm** works with `AbstractArray`s. The resulting **generic implementation** immediately gives us access to a large range of datatypes.\n",
    "\n",
    "- As the **provider of a type** we only specify that our object behaves like an `AbstractArray`. Based on providing only a small set of functions in whatever optimised way we want our **new type works with a large range of algorithms**.\n",
    "\n",
    "In combination with **multiple dispatch** this leads to an ([unreasonable](https://www.youtube.com/watch?v=kc9HwsxE1OY)) amount of code reuse, for example:\n",
    "\n",
    "1. **Sharing types** (DataStructures.jl, OrderedCollections.jl, StaticArrays.jl, Colors.jl, Measurements.jl)\n",
    "2. **Sharing algorithms** (StatsBase.jl, SortingAlgorithms.jl, GenericLinearAlgebra.jl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code decoupling by multiple dispatch\n",
    "\n",
    "Now you might rightfully wonder why I emphasised *multiple dispatch* here. After all this sounds very much just like the generic programming available in other languages. Let me illustrate this with an example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We consider the case of a sparse compressed-stored column matrix type `SparseMatrixCSC`, which is available in the `SparseArrays` Julia standard library. Since we care about precision we want to use the `Double64` floating-point type (from `DoubleFloats`), which implements a fast datatype for higher-precision floating-point arithmetic.\n",
    "\n",
    "- To make these two datastructures interoperate, we need a set of methods. We consider here the scalar-matrix multiplication `*(factor::Double64, matrix::SparseMatrixCSC)`.\n",
    "\n",
    "- In many languages such as Python, implementing such a binary operator translates to a special operator function. In the Python case this is `__mul__`, which needs to be implemented for the LHS type, i.e. `Double64`. In other words which method to call for `*` is decided solely by considering the type of the *first* argument, hence the name *single dispatch*.\n",
    "\n",
    "- This is not always the most clever thing to do. In our case one could well argue that `SparseMatrixCSC` should be the place where the operator should be implemented, because the net effect will be to scale this matrix. Amongst other reasons this is why Python rolled out a second set of functions, called `__rmul__`, which dispatch on the *second* argument. Problem solved?\n",
    "\n",
    "- Well, not really. The issue is that the algorithm in either case ends up being closely coupled to one of the data structures. Moreover since these datastructures live in third-party packages (that know nothing from each other) one either needs to either\n",
    "   1. Implement `__rmul__` in `SparseMatrixCSC` and add a dependency `SparseArrays` -> `Double64`.\n",
    "   2. Implement `__mul__` in `Double64`: Same thing in reverse\n",
    "   3. Derive off e.g. `SparseMatrixCSC` in my code and implement there. Works, but now my code likely ends up in  close coupling to the internals of the `SparseArrays` package.\n",
    "\n",
    "The aforementioned problem is known as the *binary method problem* and is elegantly solved by multiple dispatch, because both arguments are treated on the same footing. So in Julia I just need to do\n",
    "```\n",
    "import Base: *\n",
    "*(n::Double64, str::SparseMatrixCSC) = ...\n",
    "```\n",
    "in my code, without changing any `Double64` or `SparseMatrixCSC` code. Notice how again algorithm and data types are well-separated and can be independently distributed as packages.\n",
    "\n",
    "- As an asided In C++ the very problem of implementing multiplication operators is not so prominent, because `operator*` can also be a free function. But: For other binary operations (like type conversion) the problem is exactly the same: As soon as two non-standard third-party types are involved and my code needs to connect them, I have to make one more special than the other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Workflow advantages by multiple dispatch\n",
    "\n",
    "Commonly in Julia for each algorithm there exists a generic fallback implementation, which is slow, but always works. An example is the `det`-function we saw earlier. For specific types, where more structure is known, appropriate specialisations are available (like for `StaticArray`s or `Tridiagonal` matrices),\n",
    "which more efficient.\n",
    "\n",
    "Let's assume now we program an algorithm, which involves a determinant, e.g. the simple normalisation function\n",
    "```\n",
    "normalise(A) = A / det(A)\n",
    "```\n",
    "(where of course we assume `A` to be nonsingular). We note:\n",
    "- As a user, if I know the structure of my matrix is special, I mark it e.g. with the `Tridiagonal` type. When running `normalise` on it, the `det`-method specialised for `Tridiagonal` will be used, so that the best possible performance results.\n",
    "- As a user, if I do not know the structure of the matrix, I leave it general. Then the slow `det` will still give me a sensible result.\n",
    "- As a programmer, writing the `normalise` algorithm I do not need to know what kind of matrix my users will pass me, since the dispatch of `det` will assure the best thing happens. Most notably even if I have standard matrices in mind when writing the code, my `normalise` will still work if a passed matrix `A` is non-standard, because the fallback `det` will do the job. \n",
    "- As a programmer, writing a custom matrix `A` which adheres to all conventions of the `AbstractMatrix` interface, I instantly have a working `normalise` function even if I could be more clever about it. This way I can rapidly prototype my fancy new matrix type and only later implement the `det` function once I see it's necessary to achieve better performance.\n",
    "- This allows for a gradual performance increase in codes by **making only local changes** where profiling flags up issues.\n",
    "\n",
    "- This leads to the general Julia mantra:\n",
    "\n",
    "<div style=\"text-align:center;\"><b>First get it to work, then get it to work <i>fast</i>.</b></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### More details\n",
    "- PhD Thesis from Jeff Bezanson\n",
    "- https://www.youtube.com/watch?v=kc9HwsxE1OY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Takeaways\n",
    "\n",
    "- Generic code in julia can be fast\n",
    "- Everyone can write generic code and everyone should do it.\n",
    "- Multiple dispatch: Code reuse and emergent features\n",
    "- Multiple dispatch allows gradual, local performance engineering."
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.7.2",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}