5. Custom types.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project: Defining new types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "log (generic function with 20 methods)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "include(\"dual.jl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's implement a very basic version of dual numbers (done much better [here](https://github.com/JuliaDiff/ForwardDiff.jl). A dual number is written as \n",
    "$$a + b\\epsilon$$\n",
    "where $a$ and $b$ are real numbers and \n",
    "$$\\epsilon^2 = 0$$\n",
    "Basically, a dual number encodes both a point in space and a local derivative.\n",
    "\n",
    "To start, we'll assume that dual numbers are only a generalization of `Float64`. \n",
    "\n",
    "We will also define (though it's not rigorously true) that these numbers are subtypes of `Real`. That way, any function that admits a `Real` will also admit a `Dual`."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "immutable Dual <: Real\n",
    "    val::Float64\n",
    "    adj::Float64\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes:\n",
    "- Again, we choose `Dual <: Real` for convenience, not math. We want `Dual`s to go where `Real`s are allowed.\n",
    "- Instead of `immutable` we could have used `type`. `immutable` means that we cannot change these numbers, but as a consequence, they can be stack-allocated and faster.\n",
    "- We need to give Julia the types of the *fields*. (If we don't, they're assumed to be `Any`.) In contrast to function arguments, we absolutely want these to be concrete, if possible, since abstract fields lead to a performance hit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0 + 2.0ϵ"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aa = Dual(1, 2)  # note implicit conversion to Float64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Symbol,1}:\n",
       " :val\n",
       " :adj"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fieldnames(aa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.0,2.0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aa.val, aa.adj"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Base.show(io::IO, x::Dual) = print(io, \"$(x.val) + $(x.adj)ϵ\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0 + 2.0ϵ"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "aa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making `Dual` work: conversion and promotion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make `Dual` \"just work\" we need to define how to upconvert other numbers to it. We can do this by extending `convert`."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Base.convert(::Type{Dual}, x::Real) = Dual(x, 0)  # first argument is not named, but is if type Type{Dual}\n",
    "Base.convert(::Type{Dual}, x::Dual) = x  # corner case, since Dual <: Real"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notes:\n",
    "- the first argument is not named (nothing before the `::`)\n",
    "- it is of `Type{Dual}`; `Type{T}` is a type of which the DataType `T` is the only member (that is, we're specifying tha that the call looked like `convert(Dual, 4.)` or similar\n",
    "- `x` can be any real number\n",
    "- we forward construction on to `Dual`; Real numbers just have `adj = 0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.0 + 0.0ϵ"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "convert(Dual, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to define a promotion rule. Promotion rules stipulate what happens when we need to find a type \"big enough\" to hold two things (for addition, comparison, etc.)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "Base.promote_rule(::Type{Float64}, ::Type{Dual}) = Dual\n",
    "Base.promote_rule(::Type{Int64}, ::Type{Dual}) = Dual\n",
    "Base.promote_rule(::Type{Float32}, ::Type{Dual}) = Dual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arithmetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With promotion and conversion defined, we *only* need to add arithmetic for `Dual`s. Other arithmetic operations automatically promote as needed."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "import Base: +, -, /, *, abs, <, ==, <="
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "+(x::Dual, y::Dual) = Dual(x.val + y.val, x.adj + y.adj)\n",
    "-(x::Dual) = Dual(-x.val, -x.adj)\n",
    "-(x::Dual, y::Dual) = x + -y\n",
    "*(x::Dual, y::Dual) = Dual(x.val * y.val, x.val * y.adj + x.adj + y.val)\n",
    "/(x::Dual, y::Dual) = Dual(x.val/y.val, x.adj/y.val - (x.val * y.adj)/y.val^2)\n",
    "abs(x::Dual) = x.val > 0 ? x : -x\n",
    "==(x::Dual, y::Dual) = x.val == y.val && x.adj == y.adj\n",
    "<(x::Dual, y::Dual) = x.val == y.val ? x.adj < y.adj : x.val < y.val\n",
    "<=(x::Dual, y::Dual) = x < y || x == y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.0 + 1.0ϵ"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dual(5, 2) - Dual(3, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.6666666666666667 + 0.11111111111111105ϵ"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dual(5, 2)/Dual(3, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now watch what we get for **free**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8.0 + 2.0ϵ,5.0 + 2.0ϵ)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "5 + Dual(3, 2), Dual(3, 2) + 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.5 + -4.5ϵ"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "4.5/Dual(1, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x5 Array{Dual,2}:\n",
       " 0.08117297108120636 + 0.0ϵ  …  0.7761537892483081 + 0.0ϵ \n",
       " 0.05062865704382702 + 0.0ϵ     0.3574831527572737 + 0.0ϵ \n",
       " 0.7042205664707128 + 0.0ϵ      0.7852508128888431 + 0.0ϵ \n",
       " 0.7750474721966927 + 0.0ϵ      0.45371424589626574 + 0.0ϵ\n",
       " 0.7305473566017255 + 0.0ϵ      0.7465167317384658 + 0.0ϵ "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Array{Dual}(rand(5, 5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5x5 Array{Dual,2}:\n",
       " 0.14552328892047406 + 1.5999828984728819ϵ  …  0.6566184836032978 + 1.28415435665137ϵ   \n",
       " 0.9110217627851693 + 1.1730441782246503ϵ      0.5056012271322503 + 1.1224056065104175ϵ \n",
       " 0.8297085162113442 + 1.5845091401000588ϵ      0.886713630546704 + 1.3861618580595891ϵ  \n",
       " 0.13998863884048807 + 1.210252817476958ϵ      0.21119856573946483 + 1.5117890146065958ϵ\n",
       " 0.8559721429778611 + 1.1646014444060433ϵ      0.6604833718120571 + 1.9889443641540232ϵ "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = rand(5, 5) + Dual(0, 1) * Array{Dual}(rand(5, 5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Dual,1}:\n",
       " 0.3467906904838485 + 8.156917884377215ϵ \n",
       " 0.5993701684355719 + 7.633324292347406ϵ \n",
       " 0.6716694785307615 + 8.483065746839014ϵ \n",
       " 0.24579398447062234 + 7.725375280944693ϵ\n",
       " 0.7679065688429951 + 8.526068040812186ϵ "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A * rand(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.03875850773520528 + 1.6630320021413ϵ,0.9633881400461819 + 1.1942231446682094ϵ)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minimum(A), maximum(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dual(1, 1) <= Dual(1, 0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Dual,1}:\n",
       "  0.11507610919563258 + -118.12697765176911ϵ\n",
       "  0.09480687542650598 + 956.749758164928ϵ   \n",
       "  0.8384680698058897 + -939.7213129761792ϵ  \n",
       "  0.45322362394265586 + -1038.554659367485ϵ \n",
       " -0.43926866649295715 + 216.8425258556168ϵ  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A \\ rand(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Automatic derivatives!"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "import Base.log\n",
    "log(x::Dual) = Dual(log(x.val), x.adj/x.val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "foo (generic function with 1 method)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function foo(x, y)\n",
    "    log(x + y) - y\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.9013877113318904 + 0.3333333333333333ϵ"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "foo(Dual(1, 1), 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.9013877113318904 + -0.6666666666666667ϵ"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "foo(1, Dual(2, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.6",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}