Lecture 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A Short Tour of Julia\n",
    "\n",
    "- Calling Python and R code with ```PyCall.jl``` and ```RCall.jl```\n",
    "- Benchmarking with ```BenchmarkTools.jl```\n",
    "- Plotting: ```Plots.jl``` and others libraries\n",
    "- Managing data with ```DataFrames.jl```\n",
    "- Neural networks with ```Flux.jl```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Install all the packages used in this tutorial\n",
    "using Pkg\n",
    "Pkg.activate(\".\") #create a new virtual environment in the folder where the notebook is located\n",
    "Pkg.add(\"DataFrames\")\n",
    "Pkg.add(\"Plots\")\n",
    "Pkg.add(\"Statistics\")\n",
    "Pkg.add(\"PyCall\")\n",
    "Pkg.add(\"RCall\")\n",
    "Pkg.add(\"Flux\")\n",
    "Pkg.add(\"BenchmarkTools\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DataFrames\n",
    "using Plots\n",
    "using Statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calling Python code from Julia\n",
    "\n",
    "To call Python code (including ```numpy``` and any other library) fro Julia we can use the ```PyCall``` package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using PyCall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now import Python function (from libraries or builtin) and use them as standard Julia code.\n",
    "\n",
    "We will see:\n",
    "\n",
    "- how to import a builtin function\n",
    "- how to import a Python package\n",
    "- how to use inline Python code\n",
    "- type conversions and zero-copy arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us import the builtin Python functions ```sum``` and ```map``` using the ```pybuiltin``` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pysum = pybuiltin(\"sum\")\n",
    "pymap = pybuiltin(\"map\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use them as normal Julia functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pysum([1,2,3,4,5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pymap_obj = pymap(x -> x^2, 1:10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that python ```map``` returns a map object, we can iterate on it with the Julia function ```collect``` or use it in a ```for``` cycle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i ∈ pymap_obj\n",
    "    println(i)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Type conversion happens automatically and works most of the time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing a Python Module\n",
    "\n",
    "A python module (e.g., ```math```, ```matplotlib```, ```numpy```, etc.) can be imported using the ```pyimport``` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pymath = pyimport(\"math\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can access all the functions (and Objects) in the Python ```math``` module directly from the variable ```pymath```:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pymath.cos(2*pymath.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inline Python code\n",
    "\n",
    "We can also write directly Python code, enclosing it with ```py\" \"``` or ```py\"\"\" \"\"\"```:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = py\"[ i + 5 for i in range(0,10)]\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, Python code is executed and the result is automatically converted in the correct Julia type.\n",
    "\n",
    "We can avoid automatic type conversion by using ```py\" \"o``` or ```py\"\"\" \"\"\"o```, which return a ```PyObject```:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py\"3 + 5\"o"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also define entire Python functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py\"\n",
    "def a_python_function(func, x, key=None):\n",
    "    if key is not None:\n",
    "        return func(x, key)\n",
    "    else:\n",
    "        return func(x, x)\n",
    "\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function can then be called quite easily, including the use of keyword arguments"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py\"a_python_function\"((x,y) -> x^2 + y^2, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py\"a_python_function\"((x,y) -> x^2 + y^2, 4, key=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By assigning ```py\"a_python_function\"``` to a variable we can be more \"transparent\":"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_python_function = py\"a_python_function\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a_python_function(+, 4, key=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python/Numpy Arrays\n",
    "\n",
    "Pasing arrays to Python functions especting numpy arrays (```ndarray```) works as espected **but** a copy is performed:\n",
    "\n",
    "- Julia arrays are stored in _column-major order_\n",
    "- Numpy arrays are stored in _row-major order_\n",
    "\n",
    "The ```PyCall``` module has a ```PyArray``` type that can be used with Julia and, when passed to Python functions _no copy_ is made."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py\"import numpy as np\n",
    "\n",
    "xs = np.zeros([10,10])\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PyArray(py\"xs\"o) # No copy!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extras\n",
    "\n",
    "- It is also possible to call Julia code from Python (see ```PyCall.jl``` documentation)\n",
    "- If you want to use Matplotlib for plot the package ```PyPlots.jl``` provides a wrapper using ```PyCall.jl```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calling R from Julia\n",
    "\n",
    "R code can be called easily from Julia using the ```RCall.jl``` package. In case you want to use some R-only package or the very good ```ggplot``` library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using RCall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are multiple ways to move objects to and from R and to evaluate R code.\n",
    "\n",
    "- The ```@rput``` and ```@rget``` macros\n",
    "- Importing libraries with ```@rimport``` and ```@rlibrary```\n",
    "- Enclosing R code in ```R\" \"```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make a variable visible to R we use the ```@rput``` followed by the name of the variable, which will have the same name in the R environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1,2,3,4,5]\n",
    "\n",
    "@rput x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now use the variable ```x``` in the R code that we can write inside ```R\" \"```:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R\"y = max(x)\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can retrieve the variable ```y``` from the R environment using the macro ```@rget```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@rget y\n",
    "\n",
    "println(typeof(y))\n",
    "println(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we can use string interpolation inside ```R\" \"```:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R\"z <- $([1,2,3,4] .+ 9)\"\n",
    "\n",
    "@rget z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This also works for functions from Julia to R..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function function_for_R(x)\n",
    "    sum(x)\n",
    "end\n",
    "\n",
    "@rput function_for_R\n",
    "\n",
    "R\"w <- function_for_R(c(1,7,8,9))\"\n",
    "\n",
    "@rget w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and from R to Julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R\"r_function <- function(x, y) x + y\"\n",
    "\n",
    "@rget r_function\n",
    "\n",
    "r_function(3, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is also possible to convert R dataframes to Julia dataframes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R\"\"\"\n",
    "n <- c(2, 3, 4)\n",
    "s <- c(\"a\", \"b\", \"c\")\n",
    "df <- data.frame(s, n)\n",
    "\"\"\"\n",
    "\n",
    "@rget df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "typeof(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Managing Data with DataFrames\n",
    "\n",
    "DataFrames on Julia are similar to dataframes in R and Pandas dataframes in Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us start by creating an _empty_ dataframe. We will add/load new data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = DataFrame()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An empty dataframe is quite boring, let us generate some data:\n",
    "- $x$ positions from $1$ to $10$\n",
    "- a first random $y$ coordinate (uniform in $[0,1)$)\n",
    "- a second random $y$ coordinate ($N(0,1)$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = hcat(collect(1:10), rand(10, 1), randn(10, 1))\n",
    "df = DataFrame(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can rename the columns of the dataframe by passing a vector of Strings of Symbols to the ```rename!``` function (notice the ```!```)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "names = [:x, :y₁, :y₂]\n",
    "rename!(df, names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have also added the names during the creation of the dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = DataFrame(data, [:x, :y₁, :y₂])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the different columns of the dataframe by name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.y₁ # An array of Float64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[:, :x] # This looks more similar to array/dictionary access "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.\"x\" # We can even use strings..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add a new column by simply assigning a vector of suitable length ($10$ elements in this case)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.x₂ = 10*rand(10)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We might want to rename and reorder the columns.\n",
    "\n",
    "We can do this via the ```rename``` and ```select!``` functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rename!(df, :x => :x₁)\n",
    "select!(df, r\"x\", :) # group all columns matching the regexp \"x\" before all the rest (:)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us get some statistics on this data via the ```describe``` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "describe(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can select which statstics to get via additional arguments to the ```describe``` function.\n",
    "\n",
    "This has the form _column name_ ```=>``` _function to apply_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "describe(df, :mean => mean, :median => median, :sum => sum, :product => prod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we need a matrix, instead of a dataframe, simply using the ```Matrix()``` constructor works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Matrix(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before moving on, let us start with our first scatter plot.\n",
    "\n",
    "- ```scatter``` creates a new plot\n",
    "- ```scatter!``` adds to the existing plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter(df.x₁, df.y₁, label=\"data 1\")\n",
    "scatter!(df.x₂, df.y₂, label=\"data 2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A few notes on plotting in Julia\n",
    "\n",
    "There are multiple packages that can be used for plotting in Julia:\n",
    "\n",
    "- ```Plots.jl```: the \"main\" Julia plotting library with multiple backends (including in JavaScript)\n",
    "- ```PyPlots.jl```: wrapper for Python's matplotlib\n",
    "- ```Gadfly.jl```: promising package, inspired by ggplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Manipulation of Dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = DataFrame()\n",
    "\n",
    "for i ∈ 1:10^5\n",
    "    elem::Vector{Float64} = []\n",
    "    while sum(elem) ≤ 1\n",
    "        push!(elem, rand())\n",
    "    end\n",
    "    push!(df, (id = i, length = length(elem), elements = elem))\n",
    "end\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maybe we want to add the sum of all the elements in each list as an additional column in out dataframe.\n",
    "\n",
    "Notice that ```ByRow``` indicates that the function is applied to each row of the column, not to the entire column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform!(df, :elements => ByRow(sum))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```elements_sum``` is not a good name. Let us delete the column and create it again with a different name (without using ```rename!``` which would be better)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "select!(df, :id, :length, :elements)\n",
    "transform!(df, :elements => ByRow(sum) => :sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us find the average length of the list of elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean(df.length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to prove that the expecte value is the constant $e$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MathConstants.e"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have used ```combine``` to _combine_ all elements of a column in a single value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combine(df, :length => mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us explore how we can group the different rows of the dataframe using the ```groupby``` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "grouped_df = groupby(df, :length, sort=true) |> x -> combine(x, :length => length => :num_elems)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "histogram(df.length, yaxis = :log, bar_width = 0.75, title = \"number of sequences\", key=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building a Neural Network\n",
    "\n",
    "We are going to build a simple neural network from scratch, then we are going to use the facilities provided by ```Flux.jl``` to help us build and train neural networks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Flux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic gradient computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define a function of which we want to compute the derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f(x) = 3x^3 + 2x^2 + 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the derivative by using the ```gradient``` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "derivative_f(x) = gradient(f, x)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us plot both $f$ and its derivative:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_vals = -5:0.01:5\n",
    "\n",
    "plot(x_vals, f.(x_vals), label=\"f(x)\")\n",
    "plot!(x_vals, derivative_f.(x_vals), label=\"f'(x)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we expect $9x^2 + 4x$ as a derivative and a good automatic differentiation engine will actually write the code corresponding to it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@code_llvm derivative_f(3.0) # we expect 9x^2 + 4x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an introduction to automatic diffentiation the [wikipedia page](https://en.wikipedia.org/wiki/Automatic_differentiation) provides a good overview.\n",
    "\n",
    "For one of the automatic differentiation framework in Julia that is used in Flux see [Zigote.jl](https://github.com/FluxML/Zygote.jl) and the paper describing how automatic differentiation is performed on [arXiv](https://arxiv.org/abs/1810.07951)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Neural Networks from scratch\n",
    "\n",
    "For a general introduction to machine learning a quick read is [The hundred-page machine learning book](http://themlbook.com/wiki/doku.php) where all chapter are available online. For a more in-depth course on neural networks and deep learning, we refer to the [Deep Learning course](https://atcold.github.io/pytorch-Deep-Learning/) by Yann LeCun and Alfredo Canziani."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us build a simple fully connected layer (i.e., a simple linear function) with two inputs and one output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = rand(1, 2) .- 0.5;\n",
    "b = rand(1) .- 0.5;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of this linear function is $Wx + b$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_layer(x) = W*x .+ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The error between the expected outputs $y = (y_1, \\ldots, y_n)$ and the outputs $\\hat{y} = (\\hat{y}_1, \\ldots, \\hat{y}_n)$ given by the layer is $\\frac{1}{n} \\sum_{i=1}^n (y_i - \\hat{y}_i)^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function loss(x, y)\n",
    "    ŷ = simple_layer(x)\n",
    "    mean((y .- ŷ).^2)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the gradient by using the ```gradient``` function made available by Flux, and we can decide to derive whith respect to what parameters by using ```Flux.params```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_simple(x, y) = gradient(() -> loss(x,y), Flux.params(W, b))\n",
    "\n",
    "println(d_simple([2, 3], 4)[W])\n",
    "println(d_simple([2, 3], 4)[b])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How can we train this simple neural network/linear function? by using gradient descent. Notice that here we use ```global``` like in Python to modify a variable in the global scope"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function train!(X, Y; η=0.1)\n",
    "    grad = d_simple(X,Y)\n",
    "    W̃ = grad[W]\n",
    "    b̃ = grad[b]\n",
    "    global W = W - η*W̃\n",
    "    global b = b - η*b̃\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us generate a very simple training set as $100$ random point where the target value is actually a linear function of first component plus a gaussian noise:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simple_X = rand(2, 100)\n",
    "simple_y = [10*i + randn()*0.05 + 3 for i ∈ simple_X[1,:]];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us visualize the data in 3D:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function show_data()\n",
    "    scatter3d(simple_X[1,:], simple_X[2,:], simple_y, label=\"target\")\n",
    "    simple_ŷ = simple_layer(simple_X)\n",
    "    scatter3d!(simple_X[1,:], simple_X[2,:], reshape(simple_ŷ, (100,)), label=\"predicted\")\n",
    "end\n",
    "\n",
    "show_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now train for a few epochs the network, printing the loss before and after the training:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "println(\"Loss before training : $(loss(simple_X, simple_y))\")\n",
    "for _ in 1:100\n",
    "    train!(simple_X, simple_y)\n",
    "end\n",
    "println(\"Loss after training: $(loss(simple_X, simple_y))\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MNIST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us download the MNIST dataset, which contains $60,000$ images of handwritten digits as $28x28$ greyscale images.\n",
    "\n",
    "Flux downloads the MNIST dataset in a directory inside ```~/.julia/```. If you do not have access to it you can use ```@eval(Flux.Data.MNIST, dir=“.”)``` to download in the current directory or, if you download the dataset manually in the current directory, to make Flux find it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images = Flux.Data.MNIST.images()\n",
    "labels = Flux.Data.MNIST.labels();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us see one of the images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "images[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first image has label: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some standard preprocessing:\n",
    "\n",
    "- encoding the labels as one-hot vectors of $10$ elements\n",
    "- For this example we will use only $1,000$ images instead of $60,000$\n",
    "- change the type of the images as arrays of ```Float64``` and the shape of the input as $28 \\times 28 \\times \\textit{ num channels } \\times \\textit{ num samples}$\n",
    "- prepare the minibatches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_images = 1000\n",
    "\n",
    "Y = Flux.onehotbatch(labels[1:n_images], 0:9);\n",
    "\n",
    "X = Float64.(reshape(hcat(images[1:n_images]...), (28, 28, 1, n_images)))\n",
    "\n",
    "batches = Flux.Data.DataLoader((X, Y), batchsize=32);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now build our model as a convolutional neural network with:\n",
    "\n",
    "- convolutional layers, with a $3x3$ kernel, a padding of $1$ in all directions and with _input channels_ ```=>``` _output channels_\n",
    "- max pooling layers\n",
    "- a dense layer with $288$ inputs and $10$ outpus followed by a softmax layer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Chain(\n",
    "    Conv((3,3), 1 => 16, relu, pad=(1,1)),\n",
    "    MaxPool((2,2)),\n",
    "    Conv((3,3), 16 => 32, relu, pad=(1,1)),\n",
    "    MaxPool((2,2)),\n",
    "    Conv((3,3), 32 => 32, relu, pad=(1,1)),\n",
    "    MaxPool((2,2)),\n",
    "    Flux.flatten,\n",
    "    Dense(288, 10, identity),\n",
    "    softmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can define the loss as the _crossentropy_ loss. Notice that the model is now included in the definition of the loss function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function loss(x, y)\n",
    "    ŷ = model(x)\n",
    "    Flux.Losses.crossentropy(y, ŷ)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to the loss function, we are interested in the accuracy of the prediction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function accuracy(x, y)\n",
    "    mean(Flux.onecold(model(x)) .== Flux.onecold(y))\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We decide which optimizer to use (e.g., ADAM, ADAGrad, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "optim = Flux.ADAM()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also define a callback function to be called at most once every ```n_seconds``` during training to print the current value of the loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_seconds = 5\n",
    "cb = Flux.throttle(() -> println(\"Current loss: $(loss(X, Y))\"), n_seconds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How good is our untrained network?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now train the network. Notice that we also have a macro ```Flux.@epochs num_epochs code``` available"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i ∈ 1:5\n",
    "    println(\"Epoch $i\")\n",
    "    Flux.train!(loss, Flux.params(model), batches, optim, cb=cb)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that our accuracy has improved (but not by a lot, we had a very short learning phase):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracy(X, Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Benchmarking\n",
    "\n",
    "We have seen a few ways of exploring how much time a certain operation requires in Julia, using the ```@time``` or the ```@timed``` macros.\n",
    "\n",
    "Let us start by benchmarking this function ```my_sum``` with a $10^6$ vector of random elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function my_sum(v)\n",
    "    s = 0.0 # zero(eltype(v)) # would be better since it will use the \"correct\" zero\n",
    "    for x ∈ v\n",
    "        s += x\n",
    "    end\n",
    "    s\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rand_vec = rand(10^6);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us time the function using the ```@time``` macro:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@time my_sum(rand_vec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@time my_sum(rand_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since Julia is JIT-compiled the first execution includes the compilation and might not be representative of the successive execution. Furthermore, we need more than one execution to get some significant result!\n",
    "\n",
    "We can use the ```BenchmarkTools.jl``` package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using BenchmarkTools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have the ```@benchmark``` macro that executes the code multiple times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@benchmark my_sum(rand_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us compare this implementation with the Python implementation of ```sum```. First of all, we convert the array to a Python array before calling the function to avoid type conversion overhead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "py_rand_vec = PyVector(rand_vec);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@benchmark pysum(py_rand_vec) # recall that we defined pysum previously"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that there is still some overhead, since we are calling python code from Julia.\n",
    "\n",
    "But what about the Julia native sum implementation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@benchmark sum(rand_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The End"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter(randn(3,2000) .+ [-3,0,3], randn(3, 2000) .+ [-3,2,-3], \n",
    "        c=palette(:default)[2:4], key=:none, grid=false, showaxis=false,\n",
    "        ticks=false, size=(600,600), markerstrokewidth=0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.3",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}