{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Introduction to Types and Generic Programming"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contents\n",
"\n",
"- [Introduction to Types and Generic Programming](#Introduction-to-Types-and-Generic-Programming) \n",
" - [Overview](#Overview) \n",
" - [Finding and Interpreting Types](#Finding-and-Interpreting-Types) \n",
" - [The Type Hierarchy](#The-Type-Hierarchy) \n",
" - [Deducing and Declaring Types](#Deducing-and-Declaring-Types) \n",
" - [Creating New Types](#Creating-New-Types) \n",
" - [Introduction to Multiple Dispatch](#Introduction-to-Multiple-Dispatch) \n",
" - [Exercises](#Exercises) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Overview\n",
"\n",
"In Julia, arrays and tuples are the most important data type for working with numerical data.\n",
"\n",
"In this lecture we give more details on\n",
"\n",
"- declaring types \n",
"- abstract types \n",
"- motivation for generic programming \n",
"- multiple dispatch \n",
"- building user-defined types "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"hide-output": true
},
"outputs": [],
"source": [
"using InstantiateFromURL\n",
"# optionally add arguments to force installation: instantiate = true, precompile = true\n",
"github_project(\"QuantEcon/quantecon-notebooks-julia\", version = \"0.8.0\")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"using LinearAlgebra, Statistics"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finding and Interpreting Types"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Finding The Type\n",
"\n",
"As we have seen in the previous lectures, in Julia all values have a type, which can be queried using the `typeof` function"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(1) = Int64\n",
"typeof(1.0) = Float64\n"
]
}
],
"source": [
"@show typeof(1)\n",
"@show typeof(1.0);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The hard-coded values `1` and `1.0` are called literals in a programming\n",
"language, and the compiler deduces their types (`Int64` and `Float64` respectively in the example above).\n",
"\n",
"You can also query the type of a value"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"Int64"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 1\n",
"typeof(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The name `x` binds to the value `1`, created as a literal."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Parametric Types\n",
"\n",
"(See [parametric types documentation](https://docs.julialang.org/en/v1/manual/types/#Parametric-Types-1)).\n",
"\n",
"The next two types use curly bracket notation to express the fact that they are *parametric*"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(1.0 + 1im) = Complex{Float64}\n",
"typeof(ones(2, 2)) = Array{Float64,2}\n"
]
}
],
"source": [
"@show typeof(1.0 + 1im)\n",
"@show typeof(ones(2, 2));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will learn more details about [generic programming](../more_julia/generic_programming.html) later, but the key is to interpret the curly brackets as swappable parameters for a given type.\n",
"\n",
"For example, `Array{Float64, 2}` can be read as\n",
"\n",
"1. `Array` is a parametric type representing a dense array, where the first parameter is the type stored, and the second is the number of dimensions. \n",
"1. `Float64` is a concrete type declaring that the data stored will be a particular size of floating point. \n",
"1. `2` is the number of dimensions of that array. \n",
"\n",
"\n",
"A concrete type is one where values can be created by the compiler (equivalently, one which can be the result of `typeof(x)` for some object `x`).\n",
"\n",
"Values of a **parametric type** cannot be concretely constructed unless all of the parameters are given (themselves with concrete types).\n",
"\n",
"In the case of `Complex{Float64}`\n",
"\n",
"1. `Complex` is an abstract complex number type. \n",
"1. `Float64` is a concrete type declaring what the type of the real and imaginary parts of the value should store. \n",
"\n",
"\n",
"Another type to consider is the `Tuple` and `NamedTuple`"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x) = Tuple{Int64,Float64,String}\n"
]
},
{
"data": {
"text/plain": [
"Tuple{Int64,Float64,String}"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = (1, 2.0, \"test\")\n",
"@show typeof(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, `Tuple` is the parametric type, and the three parameters are a list of the types of each value.\n",
"\n",
"For a named tuple"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x) = NamedTuple{(:a, :b, :c),Tuple{Int64,Float64,String}}\n"
]
},
{
"data": {
"text/plain": [
"NamedTuple{(:a, :b, :c),Tuple{Int64,Float64,String}}"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = (a = 1, b = 2.0, c = \"test\")\n",
"@show typeof(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The parametric `NamedTuple` type contains two parameters: first a list of names for each field of the tuple, and second the underlying `Tuple` type to store the values.\n",
"\n",
"Anytime a value is prefixed by a colon, as in the `:a` above, the type is `Symbol` – a special kind of string used by the compiler."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"Symbol"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"typeof(:a)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Remark:** Note that, by convention, type names use CamelCase – `Array`, `AbstractArray`, etc."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Variables, Types, and Values\n",
"\n",
"Since variables and functions are lower case by convention, this can be used to easily identify types when reading code and output.\n",
"\n",
"After assigning a variable name to a value, we can query the type of the\n",
"value via the name."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x) = Int64\n"
]
}
],
"source": [
"x = 42\n",
"@show typeof(x);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Thus, `x` is just a symbol bound to a value of type `Int64`.\n",
"\n",
"We can *rebind* the symbol `x` to any other value, of the same type or otherwise."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"42.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = 42.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now `x` “points to” another value, of type `Float64`"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"Float64"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"typeof(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, beyond a few notable exceptions (e.g. `nothing` used for [error handling](fundamental_types.html#error-handling)),\n",
"changing types is usually a symptom of poorly organized code, and makes\n",
"[type inference](#type-inference) more difficult for the compiler."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Type Hierarchy\n",
"\n",
"Let’s discuss how types are organized."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Abstract vs Concrete Types\n",
"\n",
"(See [abstract types documentation](https://docs.julialang.org/en/v1/manual/types/#man-abstract-types-1))\n",
"\n",
"Up to this point, most of the types we have worked with (e.g., `Float64, Int64`) are examples of **concrete types**.\n",
"\n",
"Concrete types are types that we can *instantiate* – i.e., pair with data in memory.\n",
"\n",
"We will now examine **abstract types** that cannot be instantiated (e.g., `Real`, `AbstractFloat`).\n",
"\n",
"For example, while you will never have a `Real` number directly in memory, the abstract types\n",
"help us organize and work with related concrete types."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Subtypes and Supertypes\n",
"\n",
"How exactly do abstract types organize or relate different concrete types?\n",
"\n",
"In the Julia language specification, the types form a hierarchy.\n",
"\n",
"You can check if a type is a subtype of another with the `<:` operator."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Float64 <: Real = true\n",
"Int64 <: Real = true\n",
"Complex{Float64} <: Real = false\n",
"Array <: Real = false\n"
]
}
],
"source": [
"@show Float64 <: Real\n",
"@show Int64 <: Real\n",
"@show Complex{Float64} <: Real\n",
"@show Array <: Real;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the above, both `Float64` and `Int64` are **subtypes** of `Real`, whereas the `Complex` numbers are not.\n",
"\n",
"They are, however, all subtypes of `Number`"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Real <: Number = true\n",
"Float64 <: Number = true\n",
"Int64 <: Number = true\n",
"Complex{Float64} <: Number = true\n"
]
}
],
"source": [
"@show Real <: Number\n",
"@show Float64 <: Number\n",
"@show Int64 <: Number\n",
"@show Complex{Float64} <: Number;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`Number` in turn is a subtype of `Any`, which is a parent of all types."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Number <: Any"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In particular, the type tree is organized with `Any` at the top and the concrete types at the bottom.\n",
"\n",
"We never actually see *instances* of abstract types (i.e., `typeof(x)` never returns an abstract type).\n",
"\n",
"The point of abstract types is to categorize the concrete types, as well as other abstract types that sit below them in the hierarchy.\n",
"\n",
"There are some further functions to help you explore the type hierarchy, such as `show_supertypes` which walks up the tree of types to `Any` for a given type."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Int64 <: Signed <: Integer <: Real <: Number <: Any"
]
}
],
"source": [
"using Base: show_supertypes # import the function from the `Base` package\n",
"\n",
"show_supertypes(Int64)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the `subtypes` which gives a list of the available subtypes for any packages or code currently loaded"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"subtypes(Real) = Any[AbstractFloat, AbstractIrrational, Integer, Rational]\n",
"subtypes(AbstractFloat) = Any[BigFloat, Float16, Float32, Float64]\n"
]
}
],
"source": [
"@show subtypes(Real)\n",
"@show subtypes(AbstractFloat);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Deducing and Declaring Types\n",
"\n",
"We will discuss this in detail in [generic programming](../more_julia/generic_programming.html),\n",
"but much of Julia’s performance gains and generality of notation comes from its type system.\n",
"\n",
"For example"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x1) = Array{Int64,1}\n",
"typeof(x2) = Array{Float64,1}\n"
]
},
{
"data": {
"text/plain": [
"Array{Float64,1}"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x1 = [1, 2, 3]\n",
"x2 = [1.0, 2.0, 3.0]\n",
"\n",
"@show typeof(x1)\n",
"@show typeof(x2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These return `Array{Int64,1}` and `Array{Float64,1}` respectively, which the compiler is able to infer from the right hand side of the expressions.\n",
"\n",
"Given the information on the type, the compiler can work through the sequence of expressions to infer other types."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Int64,1}:\n",
" 2\n",
" 4\n",
" 6"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"f(y) = 2y # define some function\n",
"\n",
"\n",
"x = [1, 2, 3]\n",
"z = f(x) # call with an integer array - compiler deduces type"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Good Practices for Functions and Variable Types\n",
"\n",
"In order to keep many of the benefits of Julia, you will sometimes want to ensure\n",
"the compiler can always deduce a single type from any function or expression.\n",
"\n",
"An example of bad practice is to use an array to hold unrelated types"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Any,1}:\n",
" 1.0\n",
" \"test\"\n",
" 1"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = [1.0, \"test\", 1] # typically poor style"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The type of this array is `Array{Any,1}`, where `Any` means the compiler has determined that any valid Julia type can be added to the array.\n",
"\n",
"While occasionally useful, this is to be avoided whenever possible in performance sensitive code.\n",
"\n",
"The other place this can come up is in the declaration of functions.\n",
"\n",
"As an example, consider a function which returns different types depending on the arguments."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f(1) = 1.0\n",
"f(-1) = 0\n"
]
}
],
"source": [
"function f(x)\n",
" if x > 0\n",
" return 1.0\n",
" else\n",
" return 0 # probably meant `0.0`\n",
" end\n",
"end\n",
"\n",
"@show f(1)\n",
"@show f(-1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The issue here is relatively subtle: `1.0` is a floating point, while `0` is an integer.\n",
"\n",
"Consequently, given the type of `x`, the compiler cannot in general determine what type the function will return.\n",
"\n",
"This issue, called **type stability**, is at the heart of most Julia performance considerations.\n",
"\n",
"Luckily, trying to ensure that functions return the same types is also generally consistent with simple, clear code."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Manually Declaring Function and Variable Types\n",
"\n",
"(See [type declarations documentation](https://docs.julialang.org/en/v1/manual/types/#Type-Declarations-1))\n",
"\n",
"You will notice that in the lecture notes we have never directly declared any types.\n",
"\n",
"This is intentional both for exposition and as a best practice for using packages (as opposed to writing new packages, where declaring these types is very important).\n",
"\n",
"It is also in contrast to some of the sample code you will see in other Julia sources, which you will need to be able to read.\n",
"\n",
"To give an example of the declaration of types, the following are equivalent"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 9.1\n",
" 14.3"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f(x, A)\n",
" b = [5.0, 6.0]\n",
" return A * x .+ b\n",
"end\n",
"\n",
"val = f([0.1, 2.0], [1.0 2.0; 3.0 4.0])"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 9.1\n",
" 14.3"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f2(x::Vector{Float64}, A::Matrix{Float64})::Vector{Float64}\n",
" # argument and return types\n",
" b::Vector{Float64} = [5.0, 6.0]\n",
" return A * x .+ b\n",
"end\n",
"\n",
"val = f2([0.1; 2.0], [1.0 2.0; 3.0 4.0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While declaring the types may be verbose, would it ever generate faster code?\n",
"\n",
"The answer is almost never.\n",
"\n",
"Furthermore, it can lead to confusion and inefficiencies since many things that behave like vectors and matrices are not `Matrix{Float64}` and `Vector{Float64}`.\n",
"\n",
"Here, the first line works and the second line fails"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f([0.1; 2.0], [1 2; 3 4]) = [9.1, 14.3]\n",
"f([0.1; 2.0], Diagonal([1.0, 2.0])) = "
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[5.1, 10.0]\n"
]
},
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 5.1\n",
" 10.0"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@show f([0.1; 2.0], [1 2; 3 4])\n",
"@show f([0.1; 2.0], Diagonal([1.0, 2.0]))\n",
"\n",
"# f2([0.1; 2.0], [1 2; 3 4]) # not a `Float64`\n",
"# f2([0.1; 2.0], Diagonal([1.0, 2.0])) # not a `Matrix{Float64}`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating New Types\n",
"\n",
"(See [type declarations documentation](https://docs.julialang.org/en/v1/manual/types/#Type-Declarations-1))\n",
"\n",
"Up until now, we have used `NamedTuple` to collect sets of parameters for our models and examples.\n",
"\n",
"These are useful for maintaining values for model parameters,\n",
"but you will eventually need to be able to use code that creates its own types."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Syntax for Creating Concrete Types\n",
"\n",
"(See [composite types documentation](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1))\n",
"\n",
"While other sorts of types exist, we almost always use the `struct` keyword, which is for creation of composite data types\n",
"\n",
"- “Composite” refers to the fact that the data types in question can be used as collection of named fields. \n",
"- The `struct` terminology is used in a number of programming languages to refer to composite data types. \n",
"\n",
"\n",
"Let’s start with a trivial example where the `struct` we build has fields named `a, b, c`, are not typed"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"struct FooNotTyped # immutable by default, use `mutable struct` otherwise\n",
" a # BAD! not typed\n",
" b\n",
" c\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And another where the types of the fields are chosen"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"struct Foo\n",
" a::Float64\n",
" b::Int64\n",
" c::Vector{Float64}\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In either case, the compiler generates a function to create new values of the data type, called a “constructor”.\n",
"\n",
"It has the same name as the data type but uses function call notion"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(foo) = Foo\n",
"foo.a = 2.0\n",
"foo.b = 3\n",
"foo.c = [1.0, 2.0, 3.0]\n"
]
}
],
"source": [
"foo_nt = FooNotTyped(2.0, 3, [1.0, 2.0, 3.0]) # new `FooNotTyped`\n",
"foo = Foo(2.0, 3, [1.0, 2.0, 3.0]) # creates a new `Foo`\n",
"\n",
"@show typeof(foo)\n",
"@show foo.a # get the value for a field\n",
"@show foo.b\n",
"@show foo.c;\n",
"\n",
"# foo.a = 2.0 # fails since it is immutable"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You will notice two differences above for the creation of a `struct` compared to our use of `NamedTuple`.\n",
"\n",
"- Types are declared for the fields, rather than inferred by the compiler. \n",
"- The construction of a new instance has no named parameters to prevent accidental misuse if the wrong order is chosen. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Issues with Type Declarations\n",
"\n",
"Was it necessary to manually declare the types `a::Float64` in the above struct?\n",
"\n",
"The answer, in practice, is usually yes.\n",
"\n",
"Without a declaration of the type, the compiler is unable to generate efficient code, and the use of a `struct` declared without types could drop performance by orders of magnitude.\n",
"\n",
"Moreover, it is very easy to use the wrong type, or unnecessarily constrain the types.\n",
"\n",
"The first example, which is usually just as low-performance as no declaration of types at all, is to accidentally declare it with an abstract type"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"hide-output": false
},
"outputs": [],
"source": [
"struct Foo2\n",
" a::Float64\n",
" b::Integer # BAD! Not a concrete type\n",
" c::Vector{Real} # BAD! Not a concrete type\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second issue is that by choosing a type (as in the `Foo` above), you may\n",
"be unnecessarily constraining what is allowed"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f(foo) = 11.0\n",
"f(foo_nt) = 11.0\n"
]
}
],
"source": [
"f(x) = x.a + x.b + sum(x.c) # use the type\n",
"a = 2.0\n",
"b = 3\n",
"c = [1.0, 2.0, 3.0]\n",
"foo = Foo(a, b, c)\n",
"@show f(foo) # call with the foo, no problem\n",
"\n",
"# some other typed for the values\n",
"a = 2 # not a floating point but `f()` would work\n",
"b = 3\n",
"c = [1.0, 2.0, 3.0]' # transpose is not a `Vector` but `f()` would work\n",
"# foo = Foo(a, b, c) # fails to compile\n",
"\n",
"# works with `NotTyped` version, but low performance\n",
"foo_nt = FooNotTyped(a, b, c)\n",
"@show f(foo_nt);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Declaring Parametric Types (Advanced)\n",
"\n",
"(See [type parametric types documentation](https://docs.julialang.org/en/v1/manual/types/#Parametric-Types-1))\n",
"\n",
"Motivated by the above, we can create a type which can adapt to holding fields of different types."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(foo) = Foo3{Int64,Int64,Adjoint{Float64,Array{Float64,1}}}\n"
]
},
{
"data": {
"text/plain": [
"11.0"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"struct Foo3{T1, T2, T3}\n",
" a::T1 # could be any type\n",
" b::T2\n",
" c::T3\n",
"end\n",
"\n",
"# works fine\n",
"a = 2\n",
"b = 3\n",
"c = [1.0, 2.0, 3.0]' # transpose is not a `Vector` but `f()` would work\n",
"foo = Foo3(a, b, c)\n",
"@show typeof(foo)\n",
"f(foo)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, this is probably too flexible, and the `f` function might not work on an arbitrary set of `a, b, c`.\n",
"\n",
"You could constrain the types based on the abstract parent type using the `<:` operator"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(foo) = Foo4{Int64,Int64,Adjoint{Float64,Array{Float64,1}}}\n"
]
},
{
"data": {
"text/plain": [
"11.0"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"struct Foo4{T1 <: Real, T2 <: Real, T3 <: AbstractVecOrMat{<:Real}}\n",
" a::T1\n",
" b::T2\n",
" c::T3 # should check dimensions as well\n",
"end\n",
"foo = Foo4(a, b, c) # no problem, and high performance\n",
"@show typeof(foo)\n",
"f(foo)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This ensures that\n",
"\n",
"- `a` and `b` are a subtype of `Real`, and `+` in the definition of `f` works \n",
"- `c` is a one dimensional abstract array of `Real` values \n",
"\n",
"\n",
"The code works, and is equivalent in performance to a `NamedTuple`, but is more verbose and error prone."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Keyword Argument Constructors (Advanced)\n",
"\n",
"There is no way to avoid learning parametric types to achieve high performance code.\n",
"\n",
"However, the other issue where constructor arguments are error-prone, can be remedied with the `Parameters.jl` library."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"foo = Foo5\n",
" a: Float64 0.1\n",
" b: Int64 2\n",
" c: Array{Float64}((3,)) [1.0, 2.0, 3.0]\n",
"\n",
"foo2 = Foo5\n",
" a: Float64 2.0\n",
" b: Int64 2\n",
" c: Array{Float64}((3,)) [1.0, 2.0, 3.0]\n",
"\n"
]
},
{
"data": {
"text/plain": [
"8.1"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Parameters\n",
"\n",
"@with_kw struct Foo5\n",
" a::Float64 = 2.0 # adds default value\n",
" b::Int64\n",
" c::Vector{Float64}\n",
"end\n",
"\n",
"foo = Foo5(a = 0.1, b = 2, c = [1.0, 2.0, 3.0])\n",
"foo2 = Foo5(c = [1.0, 2.0, 3.0], b = 2) # rearrange order, uses default values\n",
"\n",
"@show foo\n",
"@show foo2\n",
"\n",
"function f(x)\n",
" @unpack a, b, c = x # can use `@unpack` on any struct\n",
" return a + b + sum(c)\n",
"end\n",
"\n",
"f(foo)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Tips and Tricks for Writing Generic Functions\n",
"\n",
"As discussed in the previous sections, there is major advantage to never declaring a type unless it is absolutely necessary.\n",
"\n",
"The main place where it is necessary is designing code around [multiple dispatch](#intro-multiple-dispatch).\n",
"\n",
"If you are careful to write code that doesn’t unnecessarily assume types,\n",
"you will both achieve higher performance and allow seamless use of a\n",
"number of powerful libraries such as\n",
"[auto-differentiation](https://github.com/JuliaDiff/ForwardDiff.jl),\n",
"[static arrays](https://github.com/JuliaArrays/StaticArrays.jl),\n",
"[GPUs](https://github.com/JuliaGPU/CuArrays.jl),\n",
"[interval arithmetic and root finding](https://github.com/JuliaIntervals/IntervalRootFinding.jl),\n",
"[arbitrary precision numbers](https://docs.julialang.org/en/v1/manual/integers-and-floating-point-numbers/index.html#Arbitrary-Precision-Arithmetic-1),\n",
"and many more packages – including ones that have not even been written yet.\n",
"\n",
"A few simple programming patterns ensure that this is possible\n",
"\n",
"- Do not declare types when declaring variables or functions unless necessary. "
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 6.0\n",
" 7.0\n",
" 3.1"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# BAD\n",
"x = [5.0, 6.0, 2.1]\n",
"\n",
"function g(x::Array{Float64, 1}) # not generic!\n",
" y = zeros(length(x)) # not generic, hidden float!\n",
" z = Diagonal(ones(length(x))) # not generic, hidden float!\n",
" q = ones(length(x))\n",
" y .= z * x + q\n",
" return y\n",
"end\n",
"\n",
"g(x)\n",
"\n",
"# GOOD\n",
"function g2(x) # or `x::AbstractVector`\n",
" y = similar(x)\n",
" z = I\n",
" q = ones(eltype(x), length(x)) # or `fill(one(x), length(x))`\n",
" y .= z * x + q\n",
" return y\n",
"end\n",
"\n",
"g2(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"- Preallocate related vectors with `similar` where possible, and use `eltype` or `typeof`. This is important when using Multiple Dispatch given the different input types the function can call "
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{BigInt,1}:\n",
" 1\n",
" 4"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function g(x)\n",
" y = similar(x)\n",
" for i in eachindex(x)\n",
" y[i] = x[i]^2 # could broadcast\n",
" end\n",
" return y\n",
"end\n",
"\n",
"g([BigInt(1), BigInt(2)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"- Use `typeof` or `eltype` to declare a type "
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof([1.0, 2.0, 3.0]) = Array{Float64,1}\n",
"eltype([1.0, 2.0, 3.0]) = Float64\n"
]
}
],
"source": [
"@show typeof([1.0, 2.0, 3.0])\n",
"@show eltype([1.0, 2.0, 3.0]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"- Beware of hidden floating points "
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(ones(3)) = Array{Float64,1}\n",
"typeof(ones(Int64, 3)) = Array{Int64,1}\n",
"typeof(zeros(3)) = Array{Float64,1}\n",
"typeof(zeros(Int64, 3)) = Array{Int64,1}\n"
]
}
],
"source": [
"@show typeof(ones(3))\n",
"@show typeof(ones(Int64, 3))\n",
"@show typeof(zeros(3))\n",
"@show typeof(zeros(Int64, 3));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"- Use `one` and `zero` to write generic code "
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(1) = Int64\n",
"typeof(1.0) = Float64\n",
"typeof(BigFloat(1.0)) = BigFloat\n",
"typeof(one(BigFloat)) = BigFloat\n",
"typeof(zero(BigFloat)) = BigFloat\n",
"typeof(one(x)) = BigFloat\n",
"typeof(zero(x)) = BigFloat\n"
]
}
],
"source": [
"@show typeof(1)\n",
"@show typeof(1.0)\n",
"@show typeof(BigFloat(1.0))\n",
"@show typeof(one(BigFloat)) # gets multiplicative identity, passing in type\n",
"@show typeof(zero(BigFloat))\n",
"\n",
"x = BigFloat(2)\n",
"\n",
"@show typeof(one(x)) # can call with a variable for convenience\n",
"@show typeof(zero(x));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"This last example is a subtle, because of something called [type promotion](https://docs.julialang.org/en/v1/manual/conversion-and-promotion/#Promotion-1)\n",
"\n",
"- Assume reasonable type promotion exists for numeric types "
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(g(x)) = BigFloat\n"
]
}
],
"source": [
"# ACCEPTABLE\n",
"function g(x::AbstractFloat)\n",
" return x + 1.0 # assumes `1.0` can be converted to something compatible with `typeof(x)`\n",
"end\n",
"\n",
"x = BigFloat(1.0)\n",
"\n",
"@show typeof(g(x)); # this has \"promoted\" the `1.0` to a `BigFloat`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But sometimes assuming promotion is not enough "
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(g2(x)) = BigFloat\n",
"typeof(g2(x2)) = Nothing\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(g3(x)) = BigFloat\n",
"typeof(g3(x2)) = Nothing\n"
]
}
],
"source": [
"# BAD\n",
"function g2(x::AbstractFloat)\n",
" if x > 0.0 # can't efficiently call with `x::Integer`\n",
" return x + 1.0 # OK - assumes you can promote `Float64` to `AbstractFloat`\n",
" otherwise\n",
" return 0 # BAD! Returns a `Int64`\n",
" end\n",
"end\n",
"\n",
"x = BigFloat(1.0)\n",
"x2 = BigFloat(-1.0)\n",
"\n",
"@show typeof(g2(x))\n",
"@show typeof(g2(x2)) # type unstable\n",
"\n",
"# GOOD\n",
"function g3(x) #\n",
" if x > zero(x) # any type with an additive identity\n",
" return x + one(x) # more general but less important of a change\n",
" otherwise\n",
" return zero(x)\n",
" end\n",
"end\n",
"\n",
"@show typeof(g3(x))\n",
"@show typeof(g3(x2)); # type stable"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"These patterns are relatively straightforward, but generic programming can be thought of\n",
"as a Leontief production function: if *any* of the functions you write or call are not\n",
"precise enough, then it may break the chain.\n",
"\n",
"This is all the more reason to exploit carefully designed packages rather than “do-it-yourself”."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A Digression on Style and Naming\n",
"\n",
"The previous section helps to establish some of the reasoning behind the style\n",
"choices in these lectures: “be aware of types, but avoid declaring them”.\n",
"\n",
"The purpose of this is threefold:\n",
"\n",
"- Provide easy to read code with minimal “syntactic noise” and a clear correspondence to the math. \n",
"- Ensure that code is sufficiently generic to exploit other packages and types. \n",
"- Avoid common mistakes and unnecessary performance degradations. \n",
"\n",
"\n",
"This is just one of many decisions and patterns to ensure that your code is consistent and clear.\n",
"\n",
"The best resource is to carefully read other peoples code, but a few sources to review are\n",
"\n",
"- [Julia Style Guide](https://docs.julialang.org/en/v1/manual/style-guide/). \n",
"- [Invenia Blue Style Guide](https://github.com/invenia/BlueStyle). \n",
"- [Julia Praxis Naming Guides](https://github.com/JuliaPraxis/Naming/tree/master/guides). \n",
"- [QuantEcon Style Guide](https://github.com/QuantEcon/lecture-source-jl/blob/master/style.md) used in these lectures. \n",
"\n",
"\n",
"Now why would we emphasize naming and style as a crucial part of the lectures?\n",
"\n",
"Because it is an essential tool for creating research that is\n",
"**reproducible** and [**correct**](https://en.wikipedia.org/wiki/Correctness_%28computer_science%29).\n",
"\n",
"Some helpful ways to think about this are\n",
"\n",
"- **Clearly written code is easier to review for errors**: The first-order\n",
" concern of any code is that it correctly implements the whiteboard math. \n",
"- **Code is read many more times than it is written**: Saving a few keystrokes\n",
" in typing a variable name is never worth it, nor is a divergence from the\n",
" mathematical notation where a single symbol for a variable name would map better to the model. \n",
"- **Write code to be read in the future, not today**: If you are not sure\n",
" anyone else will read the code, then write it for an ignorant future version\n",
" of yourself who may have forgotten everything, and is likely to misuse the code. \n",
"- **Maintain the correspondence between the whiteboard math and the code**:\n",
" For example, if you change notation in your model, then immediately update\n",
" all variables in the code to reflect it. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Commenting Code\n",
"\n",
"One common mistake people make when trying to apply these goals is to add in a large number of comments.\n",
"\n",
"Over the years, developers have found that excess comments in code (and *especially* big comment headers used before every function declaration) can make code *harder* to read.\n",
"\n",
"The issue is one of syntactic noise: if most of the comments are redundant given clear variable and function names, then the comments make it more difficult to mentally parse and read the code.\n",
"\n",
"If you examine Julia code in packages and the core language, you will see a great amount of care taken in function and variable names, and comments are only added where helpful.\n",
"\n",
"For creating packages that you intend others to use, instead of a comment header, you should use [docstrings](https://docs.julialang.org/en/v1/manual/documentation/index.html#Syntax-Guide-1).\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Introduction to Multiple Dispatch\n",
"\n",
"One of the defining features of Julia is **multiple dispatch**, whereby the same function name can do different things depending on the underlying types.\n",
"\n",
"Without realizing it, in nearly every function call within packages or the standard library you have used this feature.\n",
"\n",
"To see this in action, consider the absolute value function `abs`"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"abs(-1) = 1\n",
"abs(-1.0) = 1.0\n",
"abs(0.0 - 1.0im) = 1.0\n"
]
}
],
"source": [
"@show abs(-1) # `Int64`\n",
"@show abs(-1.0) # `Float64`\n",
"@show abs(0.0 - 1.0im); # `Complex{Float64}`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In all of these cases, the `abs` function has specialized code depending on the type passed in.\n",
"\n",
"To do this, a function specifies different **methods** which operate on a particular set of types.\n",
"\n",
"Unlike most cases we have seen before, this requires a type annotation.\n",
"\n",
"To rewrite the `abs` function"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ourabs(-1) = 1\n",
"ourabs(-1.0) = 1.0\n",
"ourabs(1.0 - 2.0im) = 2.23606797749979\n"
]
}
],
"source": [
"function ourabs(x::Real)\n",
" if x > zero(x) # note, not 0!\n",
" return x\n",
" else\n",
" return -x\n",
" end\n",
"end\n",
"\n",
"function ourabs(x::Complex)\n",
" sqrt(real(x)^2 + imag(x)^2)\n",
"end\n",
"\n",
"@show ourabs(-1) # `Int64`\n",
"@show ourabs(-1.0) # `Float64`\n",
"@show ourabs(1.0 - 2.0im); # `Complex{Float64}`"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that in the above, `x` works for any type of `Real`, including `Int64`, `Float64`, and ones you may not have realized exist"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x) = Rational{Int64}\n",
"ourabs(x) = "
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"2//3\n"
]
}
],
"source": [
"x = -2//3 # `Rational` number, -2/3\n",
"@show typeof(x)\n",
"@show ourabs(x);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You will also note that we used an abstract type, `Real`, and an incomplete\n",
"parametric type, `Complex`, when defining the above functions.\n",
"\n",
"Unlike the creation of `struct` fields, there is no penalty in using abstract\n",
"types when you define function parameters, as they are used purely to determine which version of a function to use."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Multiple Dispatch in Algorithms (Advanced)\n",
"\n",
"If you want an algorithm to have specialized versions when given different input types, you need to declare the types for the function inputs.\n",
"\n",
"As an example where this could come up, assume that we have some grid `x` of values, the results of a function `f` applied at those values, and want to calculate an approximate derivative using forward differences.\n",
"\n",
"In that case, given $ x_n, x_{n+1}, f(x_n) $ and $ f(x_{n+1}) $, the forward-difference approximation of the derivative is\n",
"\n",
"$$\n",
"f'(x_n) \\approx \\frac{f(x_{n+1}) - f(x_n)}{x_{n+1} - x_n}\n",
"$$\n",
"\n",
"To implement this calculation for a vector of inputs, we notice that there is a specialized implementation if the grid is uniform.\n",
"\n",
"The uniform grid can be implemented using an `AbstractRange`, which we can analyze with\n",
"`typeof`, `supertype` and `show_supertypes`."
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x) = StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}\n",
"typeof(x_2) = StepRange{Int64,Int64}\n",
"supertype(typeof(x)) = AbstractRange{Float64}\n"
]
},
{
"data": {
"text/plain": [
"AbstractRange{Float64}"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = range(0.0, 1.0, length = 20)\n",
"x_2 = 1:1:20 # if integers\n",
"\n",
"@show typeof(x)\n",
"@show typeof(x_2)\n",
"@show supertype(typeof(x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To see the entire tree about a particular type, use `show_supertypes`."
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}} <: AbstractRange{Float64} <: AbstractArray{Float64,1} <: Any"
]
}
],
"source": [
"show_supertypes(typeof(x)) # or typeof(x) |> show_supertypes"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"StepRange{Int64,Int64} <: OrdinalRange{Int64,Int64} <: AbstractRange{Int64} <: AbstractArray{Int64,1} <: Any"
]
}
],
"source": [
"show_supertypes(typeof(x_2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The types of the range objects can be very complicated, but are both subtypes of `AbstractRange`."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x) <: AbstractRange = true\n",
"typeof(x_2) <: AbstractRange = true\n"
]
}
],
"source": [
"@show typeof(x) <: AbstractRange\n",
"@show typeof(x_2) <: AbstractRange;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"While you may not know the exact concrete type, any `AbstractRange` has an informal set of operations that are available."
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"minimum(x) = 0.0\n",
"maximum(x) = 1.0\n",
"length(x) = 20\n",
"step(x) = 0.05263157894736842\n"
]
}
],
"source": [
"@show minimum(x)\n",
"@show maximum(x)\n",
"@show length(x)\n",
"@show step(x);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly, there are a number of operations available for any `AbstractVector`, such as `length`."
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(f_x) = Array{Float64,1}\n",
"supertype(typeof(f_x)) = DenseArray{Float64,1}\n",
"supertype(supertype(typeof(f_x))) = AbstractArray{Float64,1}\n",
"length(f_x) = 20\n"
]
}
],
"source": [
"f(x) = x^2\n",
"f_x = f.(x) # calculating at the range values\n",
"\n",
"@show typeof(f_x)\n",
"@show supertype(typeof(f_x))\n",
"@show supertype(supertype(typeof(f_x))) # walk up tree again!\n",
"@show length(f_x); # and many more"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Array{Float64,1} <: DenseArray{Float64,1} <: AbstractArray{Float64,1} <: Any"
]
}
],
"source": [
"show_supertypes(typeof(f_x))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are also many functions that can use any `AbstractArray`, such as `diff`."
]
},
{
"cell_type": "markdown",
"metadata": {
"hide-output": false
},
"source": [
"```julia\n",
"?diff\n",
"\n",
"search: diff symdiff setdiff symdiff! setdiff! Cptrdiff_t\n",
"\n",
"diff(A::AbstractVector) # finite difference operator of matrix or vector A\n",
"\n",
"# if A is a matrix, specify the dimension over which to operate with the dims keyword argument\n",
"diff(A::AbstractMatrix; dims::Integer)\n",
"```\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hence, we can call this function for anything of type `AbstractVector`.\n",
"\n",
"Finally, we can make a high performance specialization for any `AbstractVector` and `AbstractRange`."
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"slopes (generic function with 1 method)"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"slopes(f_x::AbstractVector, x::AbstractRange) = diff(f_x) / step(x)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can use auto-differentiation to compare the results."
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdd2ATZR8H8Oe5uyRNm66kLRQotGXvvSmjLbSssoeCLEFQEGSDCGULvCiCiIJMRZC9N2VTNkhB9p7deyW58f4RPCqglq6Mfj9/tb9ckqfi5Zvn7hlUkiQCAABQWDHmbgAAAIA5IQgBAKBQQxACAEChhiAEAIBCDUEIAACFGoIQAAAKNQQhAAAUaghCAAAo1BCEAABQqCEIAQCgUDNbEN64ccNoNGbzYEmSRFHM1/ZA9gmCYO4mwCuSJGGVRMuBU8NyiKKY/VPDbEEYFBQUGRmZzYMFQdDr9fnaHsi+9PR0czcBXjEajQaDwdytgFdwalgOg8HA83w2D8alUQAAKNQQhAAAUKghCAEAoFBDEAIAQKHG5ebJgiDcvn371q1blStXLl++/BuPnj59+syZMz4+Ph07dmRZNjdvBAAAkE9y1SMMCAgIDAwcNGjQjh073nho8eLFPXr0SEhImDdvXs+ePXPzLgAAAPknV0G4Y8eOFy9eNGnS5I26wWCYMWPG2rVrZ82adfjw4UOHDl29ejU3bwQAAJBPchWEzs7O76z/8ccfer2+adOmhBBHR8fmzZsfPHgwN28EAACQT3J1j/CfvHjxomjRogzzKmU9PT1fvHjxxjFGo3HevHlOTk5Zi0OGDPHw8Hj7BXmeX3shknd0c7C3c1VRFyVxUUouSuqskNT58hfAv9Hr9Uql0tytAEIIMRgMWFnGcuDUsBx6vZ5lWVEUlUolpfTfD86XGKGUZj05JUl6ZztYls3+IJra4Ut08fcMjCpGXeSF0uOSqugD1v0uU/S5skiyg7uDWuWskFxUxF1FqrqSGjpS3VVywf+QAADwX/IlCD09PaOjo0VRNHUKIyMjGzdu/MYxCoVi9OjRXl5e2XlBlmUrj5ypVqvF1ESfuEg+PkqIi+TjnwjxF4yxkcL9aMnO0eBSNMOxaKyD58WUKnOelL+YqHCzozV1tIaO1tTRmjpS3OE/vhRANhkMBpVKZe5WACF/fenEP4eFwKlhOSRJYllWoVBk5+A8DsLk5OTY2NgaNWoolcqTJ082a9YsJSXl+PHj06ZNy5PXZzQuSo2LslSFv1UlSUiO5+MihbhIr6gnle798sHLRwqvsmku1W7ZVT9lKLf4BnslTpIkUlNHa7rRmjpaS0fLOiMXAQAgd0G4dOnSsLCwixcvPnr06OLFi0OHDn369OnUqVPv3bs3efLkXr169e3b9/Dhw4GBgdWrV8+rFr8DpayzjnXWEd/KpoJkyDQ8umn34M86EeuqP7nDeXqrfCunlqpx1bnq5Xhm4wNp7DlRyZJOpWgnb6aBB2WQiQAAhRXNzZ32S5cuPXjwQP61bt26HMfdunUrMDCQEBIeHn7mzBlfX9+QkJC37wV6eXmFh4dn89Ioz/NGo1GtVuegkWJmmuH+tcw7V/X3IoS4SKVvJVWZ6qqy1e9qyux+Iu16It5OklqXYLr50qASjBIr7WRDSkqKo6OjuVsBhPw1WAaX4ywETg3LkZmZmf1Lo7kKwtwosCDMSkxL1t+/pr8XkXnzIiHEoV5L+7qBT1ndjsfS7qfi+WipqSft5sN09GacsvVfr5DC2W45EIQWBaeG5UAQZovh6d30C4fTLx9XFPNxqBugru4XJ6r2PhU3PRRPvJTqedB2XkzP0kyRvHxPG4Gz3XIgCC0KTg3LgSB8DxJvzLx+Nu3CYcOjm+rqTezrBKh8KycayJ4n4vbH0qHnon8x5osqTNOiuIv4Gs52y4EgtCg4NSzHewVhYZ+OTjmFuoafuoafkBSXfjEsYf23hGEd6gX2rNuyVxnXdJ799Z445JSg5sgXVZgevriJCABgawp7j/BNkqR/+Gf6uYMZ18JVvlXs67eyq1SPstypSGluhHAxRupbjhlemSlmX6g7iPjaaznQI7QoODUsBy6N5gFJn5F+9VT62f1CUrxjq572dQIoy91NkhbfEH+9J7bzYsZWY6pqC2kc4my3HAhCi2I6NX777bf58+ebuy02qEWLFt9++202D8al0TxAVWqHei0d6rU0PLyRfHhD8oHfHJt1LtO47cKGium12dV3xLYHhJIaMqIK09mbYQtpIALAO0RHR1euXHn06NHmbohNOXbs2PHjx/PpxRGE/0HpU8lt0DRTHKYc3+rYrLNT47Yjqig+rcRsfCDOvSpOviiOqML0Kcs44L8lABBCCHF3d69Zs6a5W2FTHj9+nH9BiLEf2WKKQ91HEzLvXImc/XHq8e0K0di7DHOxI/ezH3v4uVRuI//zLVHANgAAANYGQfgeXsVhn4lyHEq80a8o3RLI7mrFrrsvVtnCb3oomruZAADwHhCE703pXfHtOKzlRo+25b5vyE6/LLbcx1+NR98QAMA6IAhz6HUc3rwQNXtg+vlDRJICi9MrnbhuPkzwPr7PMeFlurlbCQAA/wVBmCtK74puQ2Zp+0xIPb0n5vuxxsjHHEM+qcDc7a7wdSJVtxgnXBBSjOZuJQAA/DMEYR5Qelf0+GKBQ8Pg2CUTErf+KOkzNAoytRb7R2cuQU8qbeaXYRwNAIClQhDmEUrt6wYWmbBM4o2RswemXzhMCCnhQJc2Ybe1ZNfdF6tt4fc8RRgCAFgcBGFeYuwdXbsP1/aZkBK2KfbnUD4+ihBSx40ea8vNqMOMPCu0P8g/S0McAgBYEARh3lOVruox9ge7cjWjvx2evH+txBsJIZ29metduAYeTO3t/Ko7mGIBAOYkSdL06dOzf3zVqlVfvHhx69atpKQkU2XlypVjx479z3eZOnXq7t275Up6evqECRMmTJgwY8aMX3755dy5c+Za5jMrBGG+oCynadaxyJgfjC8fRc0dor9zhRCiZMikGsyRNtySG2LwfnQNAcBsRFGcNWtW9o/v16+fRqPp37//6dOnTZWkpKSYmJh/f9apU6e++eab0aNHy2mXkZExd+5cjuM4jjt//nyPHj3q16//8uXLnP0VeQVBmI9YFzdd/69cOg6K/31B3OpZYmoiIaSyKz0TwrXwZGpv55fdQtcQAArCrVu3FixYsGbNmqioqP379799wMmTJ588eUII0ev1mzZtio6OJoQkJSWZ+nMlS5a8ceNGXFzcqVOnNm3adPPmTdOzIiIi5syZ8/vvv4viOz7NVq5cOW7cOJ7nw8PDs9aHDRs2ceLExYsX37p1y9HR8ZNPPsnzv/e9IAjznV3lBkXHL+W0RSLnDEk7u58QwjFkfHUmrA237JbY5gC6hgCQv86ePevn55eUlHT79u1WrVqNGDGCUtqmTZusx+zevfvHH38khJw6dapnz57r1q0jhOzZs2fJkiWEkI8//vjhw4d6vT4yMvLBgweJiYmEkHPnzk2aNIll2Xnz5o0ZM+aNN01NTd26dWvv3r179eq1atWqdzbMzs4uNDR03759KSkp+fGHZxMWii4IVKV2DhloXzcwYd03mX+ec+05knFwquJKz4Zw31wTa2/nZ9RmP6mALyUANqvXUSGNL6D3WteCtf/7R/uMGTO++uqrESNGEEKUSuX69esZhtm2bVvWY/z9/SdPnkwICQsL+/DDD8PCwr744ouwsLCAgADTAY0bNy5WrFjXrl1NCXr27FlK6fbt21mW9ff3DwkJeWOPpA0bNtSqVcvHx6dPnz41a9b87rvvNBrN260tX768IAjPnj2rWLFiXv5XeB8IwoKj8PT2GPld8sH1UXOHuH442q5CbVPXsI0X7XdC2PFYXNaELe6ALZ0AbNCHZRijWEDXfpRvfamOiIiQh8b4+fmtX7/+7Wf5+fldv349ISEhLCxs1apVLVu2NBqNR44c+fzzz//pjWrVqsWyLCGkZMmSpkupWa1atapbt24JCQk6na5ChQqbN2/u16/f2y9i6guad09NBGHBYlin4N4q3yrxv/1PXb2Jc4dBlOWqaum5EO6ba2ItdA0BbFRbL0qI2b7mqlQqg8Fg+ln+4Q329vZ169bdsWNHTExMpUqVatSo8fvvv6ekpFSrVu2fXpbjXiUIpW/u8X7nzp2zZ8/euHFj2rRphBCe51euXPnOIDx+/LhWqy1VqlSO/rK8gc9cM1CVq1Fk7BI+LjJm4Sg+9gX5667hodbcTzfFdgf4F+m4awgAeaZx48Zbtmwx/bx58+Z/OiwgIGD69OnNmjUjhAQGBoaGhrZo0YJhXseEs7Oz6e7gf1q1alXXrl3j//Ls2bNLly7du3fvjcP27t375Zdfjh8/3tSzNBcEoXkwGme3gVMdGgTHLBydfvGIqVhNS8924GroaO1t/P5nyEIAyBvTpk3btm1b27ZtAwMD/2XOQ0BAwMOHDwMDAwkh/v7+8s+yjz76aNSoUc2aNTMNpfkngiCsXbu2V69ecsXJySk4OHj16tWmX4OCgqpUqeLg4DBmzJhp06b953zE/IZLo+ZDqUOjNkrfyvFrvs64fsa1xwhGrVEyZGYdtrUX88ERYUB5OqUmy+CmIQDkjre397Vr1y5fvuzq6pqamtqnT593HtagQYP79+8XL16cEFKtWrX79+97enqaHoqIiChevHivXr06d+4cFRXl4uJCKZWvsrq6ut65c0d+HUmSjh8//sbVzhUrVmRkZLi4uNy/f58QwjCMu7u7g4NDfvy97ws9QjNTFC3lMWoR66SNnj/M8OjV1JzGRejlTlx4lBS4l4/OMG8DAcAW2NvbN2nSpHLlyv9yDMuyvr6+pnErlFJfX1+1Wm16yNvb23T1Uq1We3t7u7i4ODs7u7u7Z32i/Docx/n6+r5xtdPFxcXT09N0pK+vr7e3t4WkIEEQWgKqULp0/tS5wyexy6cl719LJIkQ4mZH9gVzTT1pvR38uWhcJgWAvOHh4dGjRw9zt8KyIAgthbpaoyJjvtffvRqzZKKQFEcIYSmZWov9vhHT8RC/8DrWoAGAPODj4/Neq4wWBghCC8K6uLsPnasqXSX6m8/1966Ziu1LMifbcyvviB8dE9ILakIuAEDhgSC0MAzjFNxb23tc/C+z08L3mmplnOi5DpyKJY128veTcZkUACAvIQgtkapcDffh36ae2J6wcZEk8IQQO5Ys92NHVWUa7+K3PMRlUgDIR2/snWTzEIQWinPz9Bi5UExJiF0yUUx9tQFYn7LM3iBu3HlxxBmBRxoCQP74888/zb41UkFCEFouqlLrBkxRla0e/d1IY+QTU7GWGz3fkbudJAXu4yMxswIAsicmJmbz5s3Hjx9PS0v7888/33j07t27v/322+bNm02T/GSJiYl79uzZtWtXfHy8qXL16tWMjIzDhw9v3bo1ISEh6yts2rTp5MmT8lprf/7559q1a7du3fr06dP8/MvyAILQslHqFNzbqU2fmMXjMq6fNdV0KrIniGtalNbfwV+MxS1DAPgPN2/erF69+p49e5YtWxYcHNy5c+esj27YsKFly5ZXrlw5ceLEyJEj5frt27erVau2bt26TZs2Va1aNSIighASEBAQFBS0bt26zZs316xZ07SF4Zw5c9q3b3/mzJnQ0NCQkBBJkn744YdOnTpFRESEhYVNnDixgP/e94WVZayAfa3mnLZo3KoZfNQTx4DuhBCWkum12TpuYtsD/HI/tn1JfKEBsGhJu1ZIxnevdp3nnNt/TBXKrJXQ0NDBgweHhoYSQkaOHLl3796sj+7YsWPixImDBw9+43UmTZrUq1evr7/+mhAybdq08ePH79u3jxDSqlWrr776ihAybNiw2bNnjxgxYuHChTdv3nRxcSGENG7ceNeuXdu3b58zZ84biWuxEITWQeldwWPUorgV04zPH7h+MMr0f3lIKcbbkYYcFB6nkmGVkIUAlovTeUq8sYDejHnz0+D8+fPyep7BwcFvBGHbtm2HDRt2+vTpNm3ahISE2Nvbm+oXLlwYNWqU6ef27dt///33pp/lHX1bt249ZcqU8PBwtVo9Z84co9GYmpqanJz8xx9/tG/fftCgQfv27WvdunW7du2Uyr8Fs6VBEFoN1lnn/vn8hA3fxfwwXjdgCuvkSgippqXH2rJtDgh3k6QFDbAwKYCFcmjU5r8PyjeSJMm37ih982OiV69eDRs23LVr16JFi6ZNm3b58mX5WfIxlFJRFN+om3Zf4nne3d1dXp67W7du3t7eZcqUadWq1a5du2bOnDlv3rzw8HDmrXi2HJbbMngbVSi1vcaqqzaM/na44emrJW69Henp9lxEvNQtTMjAjHsAeEudOnUOHz5s+vnQoUNvPCqKoq+v74gRI06fPh0ZGfnw4UP5Wfv37zf9vG/fvrp165p+zvpSderUadiw4b1796pWrRr4F19fX1EUK1SoMHbs2PDw8PPnz8tjbSwTeoTWhlLHgO4KT+/YpVNcOn1iX9ufEOKqIgdacx+fEPz38jtbce525m4kAFiSqVOntmzZ8unTpwaD4caNG2882rVrV41G4+vre+3atfLly5crV85UnzFjRqtWrZ49e8ay7O7du/fs2WOq79y58/HjxxkZGYcOHTp58qSPj8+oUaPq16/fpUsXQsiZM2emTp369ddfly5dukSJEufOnQsKCnJzcyvIv/d9vbmtcIHx8vIKDw/38vLKzsE8zxuNRnkddCCEGJ/fj1sx3aFxW9PwGUKIRMi0y8Lae9LeILaccz5eJE1JSXF0dMy/14fsMxgMkiSZtgsAszOdGgsWLHjy5MmCBQvM3Zy/efHixdGjR3U6nYuLS9++fW/fvi0/FBMTc/bs2cjIyFKlSvn7+3Mcd+PGDVdXV09Pz/j4+GPHjomi2Lx5c1OYubm5nTx58t69e6mpqS1btpQT7u7duxcvXhQEoWbNmpUrV3758uXZs2fj4uJKly7drFmz3F8X3b59++rVq7dv357N4zMzM1mWVSgU2TkYPUJrpShe2v2L72KXfiUkxbl0GmK68D+1FltSIzbbzW8K4JoUxQ1DAHilWLFipp1yz50798ZD7u7u7du3z1qpVKmS6QetVvv2yE+O4944nhBStmzZsmXLyr96enp26tQpT1peAHCP0IqxTq7uw+Yan92L/+1/ppXYCCEDyjG/NOe6hPEbH2DtGQB4k7Ozc6NGjXL89GbNmsnDSm0GgtC6MWqN22dzJENm3MoZ8iyllsXp4dbc2PPi3KvIQgD4mwoVKqxatSrHT9+yZYtpC3tbgiC0epRT6Pp9xWpcYn4YJ6Ylm4pVtfRMCPv7A3HEGUHE4jMAAP8MQWgTGMa15xcq3yox3481bepLCClmT0+04+4kSV0OYyNDAIB/hCC0FZQ6hwy0rxsQs2g0H/PcVHNUkB2tOEcFCdrPJxfUohYAANYFQWhTHAO6Owb1ilk8zvj81RLySoasac7WcaP+e/g4vXlbBwBgiRCEtsahXkuXrsNil36lf3DdVKGELGjAtvGiTXfxL9PN2zoAAIuDILRB6qoNtX0mxq2ckRFxWi5Or832Kcv47+Wfp2HwDADAa5hQb5tUZaq5DZ4Z93OomJ7q0CDIVBxfnaGU+O0WDrdhfR0x3R4gv1y+fPmbb74xdytsyvXr1/PvxRGENkvpVdZ9+DexP30ppiU7BnQzFcdVYxwVJGCvcKg1W8YJWQiQ94YNGxYbG/vixQtzN8SmaLXaZs2a5dOLIwhtGefm6T7sf7E/fSllpju17WsqflqR4Shptls40Jqt4oosBMhjCoVi1qxZ5m4FvAfcI7RxrIub++fzM26cS96zRi4OqsDMr88E7ROuxuN+IQAUdghC28c4OLkPnZd562LSzuVy8YPSzI+NmeB9/LloZCEAFGoIwkKBsde4ffq1/u7VrFkYUopZ0ZTrcIgPj0IWAkDhhSAsLF5n4a4VcrGNF/3dn+t8mD/8HFkIAIUUgrAQYew1boNnZt64kDULm3vSTQHch0f5nY+xVQUAFEYIwsKF0Ti7D52beeNC0q6VctGvKN3ZivvklLADWQgAhQ+CsND5KwvPZ83CBh50bxA35JSw+wmukQJA4YIgLIzkLEzev1Yu1nKj+4K5gSf5fU+RhQBQiCAICylTFmZcPZU1C2vo6LaWXL8T/PGXyEIAKCwQhIXX6yw88JtcbOhB17Xguh/hz8cgCwGgUEAQFmqMxtl96JyMP05mzcKAYvTX5lzIQf5yLLIQAGwfgrCwYzQub2dhq+L0x8Zsu4P8jURkIQDYOAQhEEbj4vbZ1+mXj6cc+l0udvJm5tVjA/fyt5OQhQBgyxCEQAghrKOr+9C5aecPpZ7cKRd7l2Fm1WGD9gmPUpCFAGCzEITwCuvk6j50burRLWnnDsrF/uWY0VWZlvuEF+nIQgCwTQhCeI11cXP7dHby3jUZEafl4ueVmU8rMv57hKgMMzYNACC/IAjhbzj34rqPQxM2fq+/d00ujqrKdPOlQfv4eL0ZmwYAkC8QhPAmZclyur4T41bPMjy9Kxdn1GaDStA2B/gUoxmbBgCQ9xCE8A6qstVde34R93MoH/VULs6px9bS0db7+XSBmrFtAAB5C0EI76au0sA5ZGDMT5OE+GhThRKyuBFb2on2OsUZsE0FANgKBCH8I/s6/o7+XWN+nCikJJgqDCUrm7IOHO17XBAxjBQAbAKCEP6Nxi/EvkbT2J8miRmppgpLycqGhthM6fMzgnnbBgCQJ7jcPFmSpPXr158+fbp06dKffPKJRqORH9q3b19ERITpZ5Zlx4wZk6tmgvk4te0rGjLjVkx3GzyTKpSEECVDtgRyzXbz8yLEcdXwXQoArFuuPsWmTZs2a9asmjVrnjhxonXr1lkf2rJly4EDBxL+krtGgpm5dPyE1RaJXzObiK96gU4KsieI/fGmuOoO7hYCgHWjkpTDWz3p6enFixcPCwurVauW0WgsWbLkxo0b/fz8TI8OHDiwfPnyY8eO/aene3l5hYeHe3l5Zee9eJ43Go1qtTpnTYXckwQ+bsV0xl6j7TU2JTXV0dGREHIvWWq2W1jux7b2wjhS8zAYDJIkqVQqczcECCEkJSXFdGqA2WVmZrIsq1AosnNwznuE165do5TWqlWLEKJQKJo2bXrixImsB5w5c2bGjBnr1q3LzMzM8buAhaAsp+v/lZAQk7j1R7lYxolua8n2Pc6ficbIGQCwVjm/R/jy5Ut3d3f5Vw8Pj5cvX8q/+vj4ODo6CoKwcOHCWbNmnT179o0vSnq9fvjw4fb29lmLoaGhJUqUePu9TD3CHHdeIa/Y956Q/PNkcuh3NriXqVLFgfzUgOl0UDrQUijriH+ggmbqEQoCBi5ZhIyMDJZlzd0KICRLj9DOzo5h/qPLl/MgVCqVPM/LvxqNRgcHB/nXSZMmmX6YPHly3bp1ly9fPnLkyKxPZ1m2SZMmWq02a1Gr1b7zIg/LsgzD4PqP+alUik+mxy4aIxQtYV+/lanWwYekilLnY8yJtrQorl4XLEopLo1aDoPBgH8LCyFJkikIKf3vGzc5D8JixYpFRUUZjUbTRdhnz561aNHi7cNYlq1Tp86jR4/efGOO6969ezbvEUqSJIoivmpZAtbFzaHfV6krQhVaD7sKtU3FPuXIkzSx9UHxRDvORWneBhYuLMuaTnhzNwQIIYRlWfxbWAj2L9k5OOf3CKtVq+bh4bFnzx5CSFRU1PHjx0NCQjIzMw8cOKDX69PS0kyHJScnh4WFVa1aNcdvBJaGcS2i+zg04bf5xmf35eJXNZkWnrTTIV6Pq3QAYFVyHoQMw8ybN2/QoEG9evVq3LjxgAEDypYtGxMTExwcHBMTU6pUqTZt2nTt2rV8+fJVqlTp169f3rUZzE9ZspxrzxGxy6fKC7ARQhY0YN3tsOgMAFiZnE+fMHn8+PHFixe9vb1r165NCDEajdevX69SpUpMTMwff/xhMBhKly79zu4gpk9YL3mMeNrpPaknd7gP/5axf7WWgkEkbQ/w5ZzpD41wgaggYPqERcH0CcvxXtMnchuEOYYgtF5Zz/bE7cuMz+66DZlNuVf/wyUbSbPd/AelGSw6UwAQhBYFQWg5CmgeIQAhxKXDIMbBOWHdfPLXNyosOgMA1gVBCLlDqbb3OD4hNnnvL3KtmD3dG8R+eUE49Bx3CwHA0iEIIbeoQuk2MDT96snUU7vlYkUXujmQ++gY/2cCshAALBqCEPIA4+DkNnhmyqH1GdfPysXGReiihmz7g0JUhhmbBgDwHxCEkDc4XVHdwNCE3xcYHt+Si919mX7lmHYH+HT+X54KAGBOCELIM0qvctoPR8etmpl1cuHkmkwlV9oHkwsBwFIhCCEv2VWq59Tqw5ilX4npr3a0p4Qs92MT9NKki1hyBgAsEYIQ8phDozbqSnXjVkyTeKOpomDIpgBu6yPpp5uYUAEAFgdBCHnPOWQg6+Ke8NvryYVaFdkXzM64Ih7GhAoAsDAIQsgHlLr2/EJI+tvkQl9HusGf/fAofx0TKgDAkiAIIV9QhVL3cWj6lePpF4/IxSZF6feN2A6YUAEAlgRBCPmFcXDSDZqatONnw6ObcrGHL9O7DG1/EBMqAMBSIAghHymKlHTtNSZu5XQ+7qVcnFqbreBM+2FCBQBYBgQh5C+7CrUdA7rH/RwqZr7aq5kSsrwpG5sphV7GhAoAMD8EIeQ7TbNOqtLV4tfMIeKr6RNKhmwK5DY8kJbdwoQKADAzBCEUBOfOQyTBmLR7pVzRqci+IHbqZSHsBa6QAoA5IQihIFCW0/WblHHtTFr4XrlY2omua8H1OsrfTkIWAoDZIAihgDD2jm5DZibvX6u/e1UuNvekc+qy7Q8KCXozNg0ACjUEIRQcTuep7Tsx/tc5fOwLudivHNO+JO1+hOdxuxAAzAFBCAVKVbqqU5t+sT+HihmpcvF/9VgVQ8acxyBSADADBCEUNIcGQXbla8Wv+ZqIr5KPoWSdP3f4OQaRAoAZIAjBDFw6DiYsl7h9mVxxUpCtgezkS8Lxlxg4AwAFCkEI5sAwuj4T9Hevpp7eLdfKOVtNq50AACAASURBVNO1zbkeR/gHKchCACg4CEIwD6pS6wZNSzm4Xn/nilxsWZxOqM52PiSkYSVSACgoCEIwG05bRNtnQvyv8/joZ3LxiypMPQ/60TGsRAoABQRBCOakKl3VqV2/2OVT5ZVICSE/NGLj9dL0KxhECgAFAUEIZuZQP8iuQu341bPllUgVDNkcwK29J214gEGkAJDvEIRgfi4dBxNKk/aukStudmRLIDssXLgYiyukAJC/EIRgARhG+9H4jD9Opl8+Jteqa+nSJmzXw9jOHgDyF4IQLAJj76gbODVx60+Gp3fkYmdvpl852vkwr8ftQgDINwhCsBSKoiVduw+PXzVLTE2Ui6G1WC8HOvgUkhAA8guCECyIuloj+zr+sStmSMKriYSUkFVN2T8TpIXXMXAGAPIFghAsi1PrPqyTNmnbUrmi5siWQHZuhLD3KQbOAEDeQxCChaFU++Eo/b2IrFv4ltTQjf7cgBP8vWRkIQDkMQQhWBzT6mvJ+3/V378mF5sUpVNrse0PCslGMzYNAGwQghAsEacr6tprbPwvc4TEWLk4pCLjV5T2PSagVwgAeQhBCBbKrnwtTdOOsctDJYNeLv7QiI3TS3OuYuAMAOQZBCFYLseAbooiJRM2LpQrCoZsDOB+vCFi4AwA5BUEIVg0154j+ehnKce2ypWiarIhgMXAGQDIKwhCsGhUodQNmJx6ZHPmzQtysaEHBs4AQJ5BEIKlY13cdf2/Slj/LR/7Qi5i4AwA5BUEIVgBpU8lp6DecSumS/rXK3D/0IiNxcAZAMg1BCFYB4fGbZXeFeN/+x+RXnUCFQzZhIEzAJBrCEKwGi5dhwopSSlHt8gVDJwBgNxDEILVoCyn6z8p9cT2zFuX5CIGzgBALiEIwZqwTlpd3y8T1s3n4yLlIgbOAEBuIAjByih9KjkG9IhbOeONFWcwcAYAcgZBCNZH06yjsrjvGyvOYOAMAOQMghCskku3z/mop6mndskVDJwBgJxBEIJVogqlrv/klIPrs27VZBo40+mQkIqBMwCQbQhCsFas1kPbe2z8L3OEpDi5OKQiU9+DDjwpmLFhAGBdEIRgxVTlamr8QuJWzZD4133AHxqxD1Kkb69h4AwAZAuCEKybY0B31sU9adtSuaJiyZZAdv414dhL3CwEgP+GIAQrR6n2g1H6+9fSzh6Qa14OdG1zrtdR4VkashAA/gOCEKweVal1H09J3rPa8OSOXPQvRodXZrqGCXrcLgSAf4UgBFvAuRd37TkibvVMMTVRLo6rzpR0oF+cRRICwL9BEIKNsKvcwKFuYNyar4n4KvkoISubsideSituY+AMAPwjBCHYDqfgj6hClbRzhVzRKMi2luyXF4WLsbhZCADvhiAEG0Kptve4jOtn0i8fk2vlnOmyJmzXw0JspvkaBgAWDEEINoWx1+gGTEnc+pPx5SO52KEU07M07XmEx/4UAPA2BCHYGkUxH5dOg+NWzhAz0+Ti7Dosx5DJFzFwBgDehCAEG2Rfu4Vd+Vrxv/6PSK/6gAwl61pwGx5Imx9i4AwA/A2CEGyTc6fBUkZqSthGuaJVkS2B7NBw4UYirpACwGsIQrBNlOW0/SalntqVefOiXKyho/Prs50PCcnYngIA/oIgBJvFOrnq+n6ZsO4bPi5SLn5UhmnqSfsew7gZAHgFQQi2TOlTSePfNW7lDMlokIs/NGJjMqX5EbhZCACEIAjB5jm26MJ5FE/c+qNcUTBkvT+74Lpw5AW6hQCAIIRCQNtzpOHhjbe3p+h9jMf2FACAIATbZ9qeImn3qje2pxhRmcX2FACAIIRCgXMv7tpjRNzqmWJaslwcV53xcqCjziEJAQo1BCEUFuqqDe1rNo9bM5uIr4bJUEJWNWWPvpBW3cHAGYDCC0EIhYhzu/6UYZP3/SpXTNtTTLggXMb2FACFFYIQChNKtR+NT798NCPitFwr70wXNWS7hAlxejO2DADMBkEIhQvj4KTtOzFhw0Jj1BO52MOX6ViK9j3Gi+gWAhQ+uQrCyMjIESNGtG/fftasWXr9375O6/X62bNnt2/ffvjw4S9fvsxdIwHykrJkeee2/eJXz5YMr7co/F89NsVIZv2Bm4UAhU7Og1AUxVatWmVmZg4fPvzQoUNffPFF1kdHjhy5f//+4cOHG43GVq1aiSI+X8CCODRqo/Qqm7BhoVzhGLIxgPv5lrjvKXqFAIULlaQcnvYHDx4cMGDAkydPGIZ5+PBhpUqVnj9/rtVqCSEJCQnFixe/evVq2bJlJUny9vZeunRpcHBw1qd7eXmFh4d7eXll5714njcajWq1OmdNhbyVkpLi6Oho7lbklmQ0xCwabV83UNO0g1w8Ey11OsSHh3C+jtSMbcs+g8EgSZJKpTJ3Q4AQWzk1bENmZibLsgqFIjsH57xHePHixUaNGjEMQwjx8fHRarXXr183PfTnn386OTmVLVuWEEIpbdSo0cWLF//ttQAKHFUodf0npxz+Xf/gulxs6EEnVGc7HxIyeDM2DQAKFJfjZ0ZGRpr6fyZubm7yvcC3H4qMjHzj6RkZGV26dHnjm+ySJUt8fHzefi9Tj1AQMPHZIqSlpVFqHR2m/6C0V3f5PG71bM2QOYzTq/9jB3qT85GKj48JyxpYwV5Nph6h0WgFTS0MbOfUsH5yj9De3t7UYfsXOQ9CBweHlJQU+df09HSNRiM/lHXsTHp6uru7+xtPVyqVY8aM8fDwyFr09fV1cHB4+71wadSiSJIk/1tbveqNmJcPMjd95z5sHmVfnQ4/Nyf1d/DrnykGVbD0YdW4NGpRbOrUsHIcx2X/0mjOg9DLy+vMmTOmnw0Gw8uXL0uWLCk/9PLly8zMTDs7O0LIo0ePatWq9cbTWZZt2LBhNu8RAuQfp6BexucPknetdO74ianiwJFtgWyT3XwNHa3rji/4ADYu5194O3XqdP78+Vu3bhFCfv/995IlS1apUuXevXsbNmyoUqWKj4/P+vXrCSF37tw5e/Zsp06d8qzJAHmLUtcPRmVcP5N+IUyulXWmPzdhu4UJsZn/8kwAsAU5D0JPT8+ZM2c2adLEz89v7NixixcvppSePXt20qRJhJDFixdPmDDBz8+vcePG06ZNK1asWN61GSCPMfYa3YApiTt+NkY+loshpZiepWnPIzw2swewbTmfPmESFRX19OnTihUrmu7tGQyGjIwMZ2dnQkhaWtrNmzdLlChRtGjRt5+I6RPWy1bHiKdfOpK8/zeP0YsYu1c3qgWJtNrHN/SgM+uw5m3bP8E9Qotiq6eGNSqg6RMmRYoUqVOnjjzCRalUmlKQEOLg4FCnTp13piCABbKv7W9XrmbCb/PJX98OWUo2BnDr7ktbH2FFCACbZemD4gAKknPnIUJqcsqRzXJFpyK/+7ODTwk3E3GFFMA2IQgBXqMsp+s/KfXE9sxbl+RiPXc6ozbbPUxIwyx7AFuEIAT4G9ZJq+09Ln7t//j4KLk4pCJTz50OPIklHQBsEIIQ4E2qstUdW3SOX/O1xL9eseWHxuzdJGnhddwsBLA1CEKAd3D078a6uCdtWypX7FiyJZCdc1U4EYmbhQA2BUEI8C6Uaj8cpb8fkXbuoFwrpaG/Nuc+OCI8T0MWAtgOBCHAu1GVWttvUtKuFcbn9+ViYHE6pCLTLUww4BIpgK1AEAL8I0XRUi5dPotbOVNMf72+/Fc1maL2dNx5DJwBsBEIQoB/Y1+zmbpao/g1XxPxVR+QErKyKbvnifTLXfQKAWwBghDgPzi3/5gQknxgrVxxUZKtgezoc8KVONwsBLB6CEKA/8Iw2o/GpZ07lBFxWq5V1dJFDdkeR4REgxlbBgB5AEEI8N8YjYtuwFcJGxYao57IxQ9KM8ElaJ9jgohuIYA1QxACZIuyZHnntv3iVs4QM9Pl4jf12USD9PVV3CwEsGIIQoDscmjURuVTOWHd6+0pFAzZGMD9dFPc9xS9QgBrhSAEeA8uXYcKibEpx7bKlaJqssGf7X+Cf5iCLASwSghCgPdAOYVuwOTUo5szb1+Wi42K0PHV2W5hQibmFgJYIQQhwPthXdy1H41PWPeNkBgrF0dWYaq40kHYngLACiEIAd6bqmwNTdMOcStnZN2e4sfG7I1EackNDJwBsDIIQoCccPTvxuqKJG79Ua6oObIlgJ1xBdtTAFgZBCFAjlCq7TnS8ODPtLP75Zq3I/2lOffhUWxPAWBNEIQAOURVat3HU5J2rzY8uSMXWxann1ZkuoYJetwuBLASCEKAnOPci7v2/CJu9UwxNUkuflmD8XKgI88iCQGsA4IQIFfUVRrY12oRt+ZrIr5KPtP2FCcipRW3MXAGwAogCAFyy7lNX8owSXtWyxWNgmwNZCddFC7E4GYhgKVDEALkGsNo+0zI+ONkxpUTcq2cM13WhO0WJkRnmLFlAPDfEIQAeYBxcNJ9HJqwebHx2T25GFKK+ags7RLGG3GJFMCCIQgB8oaimI9rjxFxq2aKaclycVot1llBJlzAwBkAy4UgBMgz6mqN1dWbxK2ZLQ+cYSj5tTm387G05i56hQAWCkEIkJec239MWUXS7lVyxVVFtgayY84Jl2MxcAbAEiEIAfIUpdqPxmdEhKdfPibXqmrp0iZs1zAhNtN8DQOAf4AgBMhjjL1GN2By4pYlhqd35WJnb6azN/3wKC+gWwhgYRCEAHnPNHAmftXfVpyZW49lKAm9hIEzAJYFQQiQL9TVGqtrNov75fWKMywlvzbn1t6TtjzEwBkAC4IgBMgvzu36U06RtGulXHG3IztasZ+FCxHxuEIKYCkQhAD5hlJt7/EZ186knT8k16pr6ZJGbIdDQgwGzgBYBgQhQD5i7DVug6Yl7VxhePp6q6YuPkxPX9rlMG/AJVIAC4AgBMhfXBEv1x7D41fNyjpwZlYd1llJxp7DwBkA80MQAuQ7ddVG6lrN41ZOlwTeVGEoWdeCO/xc+vkWeoUAZoYgBCgIzm37UTuHpJ0r5IqjgmxryU66KJyMxMAZAHNCEAIUCEq1fSZk3rqUdu6gXCvnTNe24D44KjxLQxYCmA2CEKCAMHb2uv6TknatzDpwplVxOqwS0+GQkM6bsWkAhRqCEKDgKIqW0vYaHbdiupAYKxfHV2fKO9PBpzBwBsA8EIQABcquYl1N43axy6dKBr2pQglZ4cfeTpIWXMfAGQAzQBACFDTHwB6KIiXj180n0qtbg2qObA5g/xch7HuKm4UABQ1BCFDgKHXt+YWQEJMStlGuldTQzQHcgBP8vWRkIUCBQhACmAFVKHUfT0k9vSfj+lm52KgInV2XbX9QSDKYsWkAhQ6CEMA8WCet28dTEn7/zvjykVzsX45p7kn7HhdEdAsBCgqCEMBsFCXKuHYdGrd8atbV1xY1ZBMN0vQrGEQKUEAQhADmpK7hp67ZLOvqawqGbPTnVt+RNmPbQoACgSAEMDPntv0YB6fEjYvkioeabG/JfnZaOB+DK6QA+Q5BCGBulGp7jzM8uZN6apdcq6Gjq5txXQ5j9TWAfIcgBDA/qlLrBk1LOfS7/s4fcrGNF/28MhNyUEjD6msA+QlBCGAROG0Rbd8v49fO5WOey8Vx1Zi67rTPMQwiBchHCEIAS6HyrezUtn/s8qliRqpcXNyITdBLky9hEClAfkEQAlgQh/qt7MrXil/zNRFfDRlVMGRTILfpofTrPQwiBcgXCEIAy+LScTChTNLulXJFpyI7W7JjzgnHX+IKKUDeQxACWBiG0X40PuPambRzB+RaBRf6a3OuxxH+PlYiBchrCEIAi8PYa9wGz0ze+4v+zhW52Ko4/aomViIFyHsIQgBLxLl5avtOjPtlTtaVSIdVYvyL0R5HeB63CwHyDoIQwEKpfKu4dPksbvlUISVBLi5syHKUjD2PQaQAeQZBCGC57Gs2s68bGLcsVN7OnqVknT93+Ln04030CgHyBoIQwKI5BfXiipaM/+31dvZOCrKzFTvzihj2AgNnAPIAghDAslHq2vMLMT0lac9quebjSDf4s72P8neTzdcwAFuBIASwdJTldP0nZUScTj29Wy42KUrn1GPbH6ZxemrGtgHYAAQhgBVg7B3dBs9IObg+8/Zludi3LNOplNTjGNFj6AxALiAIAawDp/PU9v0y/te5WSdUzK5NvBxIn+NYlRsg5xCEAFZD5VvZtcvQrBMqKCFLG0mxmdKEC+gVAuQQghDAmqhrNrWv1ypu2RTJkGmqKBmyNZDb/0xa9CcmVADkBIIQwMo4tfqA8/SOy7JDhbOS7GzJzosQtz5CFgK8NwQhgLWh1LXHCMmoz7pDhbcj3dmSHXJKCI/C3UKA94MgBLA+lOV0AyZn3ryYHr5XLtZyo2uacV0O87eTkIUA7wFBCGCVGDsH3ceh6WEbDFkmVLT2ojPrsG32C1EZZmwagJVBEAJYK87N07nPxOQN3xme3JGLH5dnepWh7Q7wabwZmwZgTRCEAFZMUaqCU/fhcctD+aincnFabbaKlvYIw25NANmSqyBMSkoaMmRI9erVQ0JCbty4kfWhn376qftfPvzww9w1EgD+kbJCHeeQQTE/TRISY0wVSsiyJqxRJEPDMbkQ4L9xuXnyZ599ptfrN27cuH379qCgoPv37yuVStNDFy9edHV17datGyGEYdDvBMhH9nX8hcSY2KWT3T+fz9hrCCEKhmwO5Jru5udFiOOq4QQE+Dc5P0Oio6M3b968YMGC8uXLjx8/3tHRcdeuXVkPKFOmTGBgYGBgoL+/f67bCQD/xjGwh6p8rbiV0ySj4VVFQfYEsUtuiL/ewxVSgH+T8yC8ffu2Vqv18vIy/Vq3bt2IiIisB/z222+tWrUaPHjwrVu3ctVGAMgGlw6DWNci8WtmyxPti9nTvcHsmHMCdi4E+Bc5vzQaHR3t4uIi/+ri4hIdHS3/GhQU1LZtWycnp71799atW/ePP/4oXbp01qenp6c3bNiQZdmsxW3btpUrV+7t9+J53mg0CgJueFiEtLQ0SrH1j0UwGAySJBmNRtOvinYD09fOiV6/wL7DYFOlJEfWNGR6hkm7WhiruCAO8xdODcuRmZnJsqxCobC3t//P23PvF4SXL18eNmwYIaR48eKDBw9OS0uTH0pNTfXw8JB/Nd0dJIQEBATcvHlz7dq1oaGhWV9KrVZv3LjR09Mza7F48eLyXcasTEGoVqvfq7WQTyRJ0mg05m4FEPJXEKpUKrmiGRgas3iceHqHU1AvUyVYQ74TxZ6nmBPt2JIafEznI5waloPjOFMQZuvg93rp8uXLL126lBCiVCo5jouKikpOTnZyciKE3L17t379+u98VtGiRVNSUt4oUkq9vLzkK6sAkCeoSu02eGb0wlGMWqNp2sFU7FWGSTSQgL3CyfZcUXyfBPi797tH6ODgULVq1apVq5YvX7506dK1a9f+8ccfCSEXLly4fPlyly5dkpKSvvzyy+Tk5NOnT0uSZHpo69atGC8DUGAYjbPbkJkpYRszrp6Ui0MrMb3K0KB9fILejE0DsES5Gle9dOnSFStWlCpVKjg4eNmyZTqdLjU19fvvv09NTR04cKCDg4Onp2dwcPDkyZPbtGmTVy0GgP/E6TzdhsxK3LJEf/+aXJxai21ZnLbFojMAf0dN/bYckyQpJibGxcXl7Xt7qampGRkZ7u7u73yil5dXeHh4Ni+N4h6hRUlJSXF0dDR3K4CQd90jzEp/LyJu9Wz3T2cpir8aqiYRMuik8CRV2tWKU7HvfBLkHE4NyyEPlsnOwbmdaUsp9fDweOcIF41G808pCAAFQFWmmmuP4bE/T+Xjo0wVSsjSJqyzkn54VBAwhhSAEIK1RgFsm7pqI8eWPWKXTBRTE00VlpLfWrBpvDTwJKIQgBAEIYDN0zRup67ZNHZZqJiZbqooGbI1kLubJH15AXNzARCEAIWAc5u+ypLlYpd+JelfbVRoz5HdQdy+p9Lcq1iADQo7BCFAIUCpS5fPFJ7esctfL0bqoiR7g9llt8Rlt5CFUKghCAEKB0pdu33OOmvjVkyT+FdLshWzpwdbszOuiL/fRxZC4YUgBCg0KNV+OIax18Sv+VoSXs0lLO1E9wWzI88Ke59i6AwUUghCgMKEYVx7jSWSGP/LHCK+GilTxZXuaMX1P8GfiEQWQmGEIAQoXCjLaft/JRn18esXkL/W06jnTte34LqF8ZdjkYVQ6CAIAQodynK6/pOFxJiE37+Ts9C/GP2xMdv+oHAjEVkIhQuCEKAwogql26DpfOyLxG0/ycXO3syCBkzLvcL1BGQhFCIIQoBCiipVukHTDI9uJm5fJhe7+zILGjBB+5CFUIggCAEKL8bO3u3T2fp7Ecn718rF7r7Mdw2ZoH3CtXhkIRQKCEKAQo1Ra9yHzMr442TK4Q1ysZsP811DJng/shAKBQQhQGHHaJzdhs5Ju3A45chmuWjKwqD9PLIQbB6CEAAI6+jq/unXaaf3pJ7eLRe7+TCLGrLIQrB5CEIAIIQQ1sXN7bPZKYc3pp56nYVdkYVQCCAIAeAVTufp8fn81GNbkw+sk4tyFkYgC8FGIQgB4DVW6+E+/JuMP04k7VopF7v6MN83ZFvu4y9h3RmwRQhCAPgb1snVfeiczFuXEjf/IK8708WHWdKIbXMAWQg2CEEIAG9iNC7uw+Yant1L2PT921l4EVkItgVBCADvwKg17p99zce+iF87T96noosP82Njti2yEGwLghAA3o0q7dw+mSHpM+JWzpD38u3szfzYmG2Pa6RgQxCEAPCPKKfQ9v+KsFzciumS0WAqdvZmljZh2xzgw14gC8EWIAgB4N9QltP1/ZJxdIldOknSZ5iKIaWYHS25j47xGx6I5m0eQO4hCAHgvzCM9oNRCk+fmB/Gi+kpploDD3qoNTf2nLjgOrIQrBuCEACygVKXzp8qfSrFLJkopiWbapVd6Yl27NKb4oQLgnlbB4WFJBlfPMzzV0UQAkD2UOrSaYhdxTox348VkuNNNW9HGh7CHX8pDTgh8OgZQn6SjIa41bOzLvWQVxCEAPAenNv2s6/jH/P9WD4+ylTRqsjhNtyLdKlbmJCJniHkDzE9JfbHiZRldB9PyfMXRxACwPtxDOyhadYpZuEow5PbpooDR3a24lQsab2fTzKYt3Vgg/i4yOiFo5TeFbUfTaCcIs9fH0EIAO9N06Sd6wejYpdNSb9y3FRRMmS9P1tTR/128y/SMa0C8ozxxcOYxeM0fiHOIQMJpfnxFghCAMgJuwq13T+bk7RzRfL+taYKJeTbBmyv0ozfLuFuErIQ8kDm7csxSya4hAzUNGmff++CIASAHFIU8/H4YkHGtTMJGxbKy7CNr86Mrca02CtcxbZNkDtp5w4mrP2fbsAUdc2mpopEyM3EvP//CkEIADnHOus8hs8XkmJjl00RM9NNxSEVme8aMK328cdfIgshRyQpef/alIPr3T//n8q3sqmWzpPuYcLoc3k/IgtBCAC5QlVqt4FTWV3RmEWjhcRYU7GrD7O+Bdf9CL/5ISZVwHsShYSNizJvXPAYuYDzKGGqPUuT/HbzGgXZFsjl+RsiCAEg1xjWtdvnDvWDor/7wvj8vqnmX4weCObGnBOnXhbQMYRskvQZsT9PFVMS3IfNYzQupuL5GKnhTqGnL7OqKati8/5NEYQAkDc0zTq6dBoc8+OXmTcumCo1dPRiR+74SynkIJ9sNG/rwAoIyfHR349hnbW6AZOpUmUqbnootj/IL27EjK2WX4GFIASAPKOu7uc2cGrC7wvSTu8xVdzsyKHWnLeGNtnFP0xBzxD+kTHycfR3X6irNHTtOZIwLCFEImTuVXHMOXF/MNehVD6mFYIQAPKS0rui++f/Szm6JWnPatPu9hxDvm/EDqvENNrJH8HOTfAumTfOxywe79y2v1Nw71cVgfQ+Kux8Ip7vwNXU5cv0QRmCEADyGOde3GPUQsP963GrZ8u7GH5SgdkYwPU+xs+9iuEzkIUkpYRtTNj4vdvHU+xrtzDVXqRLTXfzEiFhbbgi6nxvAoIQAPIeY+/o9ulswtCYxeOEhGhT0a8oPdWeW3tPHHxKMCANgRAxPTX25ymZNy94jF6k9KlkKl6JkxruFNp40XUtWLt8GBrzNgQhAOQLqlDq+ky0r9U8esGIzJuvhs/4OtIzIVxsJvHfw0dlmLeBYGbG5/ejv/mcdXF3+/Rr1tHVVNz0UAzezy9qyEytVSAZSAghJO8nZAAAvEKppllHhVeZ+DVf29cJcG7bjzCMRkE2B7Lzror1dvBbAtk6bvl7+wcsU/qlI4nblrl0+dS+ZjNTRSJk3lVxyU1xf3C+3xR8A3qEAJC/VL5VioxdYnx2L2bJBCE5gRBCCRlfnZlfn2l7gN+EGfeFjSgk7VqZvP8396Fz5BRMMZLuYcLup+LFjgWdggRBCAAFgNE4uw2eqSpTLfrb4YaHf5qK3XyYA8HcuPPilEuCiMGkhYOQFBe9aIwx8onHqIUKT29T8VKsVHs7r1ORw605dzsztApBCAAFgmGcgnu7dv88dsWMlLCNppkVNXT0fAfu+Eup4yEhTm/uFkI+09+/Fv3tcLsKtd0GhjJqDSFEImTBdbHNAX5WHeanJvmyakx2IAgBoODYVapXZMz3GdfC41ZMEzNSCSHuduRwG66iC6mxlT/8HB1Dm5UWvjd+zWzXD0Y5Bfc2bSsYk0naH+A3PhDPhXDdfMwZRghCAChQrIu7+7D/sdoi0d+OMD5/QAhRMGRuPfaX5uyAE8KIM4I+73cXAHOS9Blxa2anhu91H7HArkJtU/HYS6nWNr60Ez3ejvN2NPOAKQQhABQ0yilcOn/qFNw75seJaecPmYotPOm1Llx0Jqmznb+GvQxthfHlo+gFIxiVvceIbzldUUIIL5Kpl4UPj/I/+7ELG7JKC0ghTJ8AAPOwr91CUbx03KqZmNm6gQAAEeRJREFUhgfXXboMpQqls5Ksb8H+clf038t/VYMdXoXB1AorJoopR7ekHN3sHDLQoV5LU+1pmtTrqKBkyKWOCk9787bvNQvIYgAorBRFSxYZtVAy6KO/HW58ds9U7FOWOd+B2/BA7HxIiM00bwMhh/iY59GLRmfeulRk1PdyCm5/LNbZzvsXowdbc5aTggRBCADmRVVqbZ8JTkG9YpdNTtq1UuKNhBAfR3qiHVddR6ptNe57isukVkWS0sL3Rn83Ul21oftnX7NaD0JIpkBGnBFGnhW3BXJTa7GW1tPHpVEAMD91DT9l6SqJm5dEzx/q+sEoZakKHEOm1mKbeTJ9jwmdvOm8emYbWw/ZJ8RHx6//VjJmeoz4Vt5c/lai1POIUMaZXu7IuarM28B3Q48QACwC6+iq6z/JKfijuOXTEjYukgyZhJAWnvRKZ+5ZGmm0i7+ViK6hRUu/cDjq28/tKtTyGP4qBQWJLLwuNt3ND6/CbA5gLTMFCYIQACyKuoZfkYk/E0Ki5n2qvxdBCNGpyJZAdkgFpulufuktLEFjiYTkhLjlU1OObnH/dLZjQHfCMISQS7FS/R38jsfiqfbcgHIWnTUW3TgAKIQYe41r9+EunT+NXzsvYeMiSZ9BCBlUgTnZnlt1R2y6m4/A5ApLkvHHyej/fcYVKekx+ntF8dKEkBQj+eKs0O4AP7wyE9aWK+dsYbcE34IgBABLZFepXpHxPxFT1/DuH4SQ8s70TAg3qDwTvJ8fcUZINpq7iYWemJoYt3JG8v61uk+mO7cfQFmOELLriVh1Cx+fSSK6KPqUtbRhMe+GIAQAC8WoTV3DIfG/zU/c/IOkz6CE9CnL/NlFQQiptJn/5S6ulJqJJKVfPhY17zPOo4THmMVKr7KEkBfpUrcwYfQ5cbkf+0tz1izLZ+cMghAALJpd5QZFxv8kGfVR8z7V37lCCHFVkYUN2R0t2R9uiC328H8mIA0LlOHxrehFo1OPbtV9PMW5XX/KKXiRLLwu1tjKV3Yl1zpzgcWtoh/4GqZPAIClY9Qa1w9GZd68mPD7d4oSZZzb9ec8StR2o2dCuLX3xMC9fHdfZmYd1lFh7obaOiExJmnXKv39COe2/ezrBJjWzr4SJw0+Jdhz5EQ7roKLlUWgCXqEAGAd7CrWKTJphV2F2jGLxyVsXCQkJzCU9CnLXOmsSNCTipv5X+5ij9/8Ihn0KWEbo+Z/zrq4FZ34s33dQEJpGk8mXBCC9vEDyzNH21prChL0CAHAilCWc2jURl2jaUrYxqi5gx0aBDu1+qCoWv1Lc/b4S2louLD+vrioIVvW4ocpWhNJyrh6KmnnckXJch6jFnLaIoQQUSLr74sTLojBJeitbgqtpU4QzCYqSea5vO7l5RUeHu7l5ZWdg3meNxqNarU6v1sF2ZGSkuLo6GjuVgAhhBgMBkmSVCor/xzKESExJvng+sxrZxwDe2j82hOGNYpk4Z/i3KvCpxWZ8dVZhwL/nm97p4bh8a3EbUuJKLp0Gqz0qUQIkQjZ/FCcdll0VpJ59djGRSz0O0dmZibLsgpFti6XIwjhvdne2W69CnMQmhgjHyftXM7HRTq37qOu4UcIeZ4mjTsvhr0QR1Rhh1ZinArwxqEtnRpCQnTS3jX6uxHObfuabgdKhOx4LIZeEu1YMq02G1zCQiPQBEEI+cuWznZrhyA00d+5krhjOVUonEMGqXwrE0IepEhzr4qbH4qDKjDjqxXQ4l62cWpI+oyUo1tST+7U+IU4BnSnCiUh5PBzaeIFQS+SyTWZrj5WMDsQQQj5yzbOdtuAIHxNFNMuHEret1blU8mpXT9O50kIeZgifXddXHdf/LA0M7EGWzSfP0Ks/dQQkuNTT+xIO7tfXbWRU5s+rKMrIeTwc2nSRSGNJ6G1rCMCTRCEkL+s/Wy3JQjCN0gGferxbSnHttpVqqdp2sE00ftxqvTttVdxOL46U8w+vz7MrffUML54mHJsa+b1s/Z1/DXNOpq+RpyKlCZfEqIzyPjqTK8yDGstGUgIQRBCfrPes932IAjfSUxPSTt7IPXUTtbZTdO0g7paY8pyURlkwXVh6U2xuy8zuSZTwiHvP9et79SQpMzbl1OPbjFGPtY07eDQsA1jryGEnI6SQi8JD1LIhOrMx+WtLAJNEISQv6zvbLddCMJ/I0n6u3+kHN9ufHrXvm6gpmkH1lkXmUHmRwir74g9SjPjqzElNXn5GW9Fp4Yk8BmXj6Uc2Uwo1TRpb183kCqUeoFseyQuvSU+SyNTajEflrbKCDR5ryDEPEIAsFGUqsrVVJWryUc/Sz21K2ruYFW5mtrmnefXrzi+OrvgmlBrGx9UgulTlgksTq33E/99iZlp6ecOpRzdzLoWcW7Xz65SfULptXhp+W1h3X2xthv9rBLTqRTDFabVVtAjhPdmRV97bR56hNknpqekndmfenoX6+zu2LSDXbVGSQL32z3x13vi01TSqwztU5b5f3v3HtTUlccB/Nyb8EhIgoQoAQqGJqBFhEUrtbxkK7Rjty5TplNbH4xsZaYd6/iH1ekfzrQz7XT6j3/stDp1O+ugaNdFmcLg0FLTTrW1WtpatHF5P6NJeCQI5H2Te/ePsIyjq4Ia7sX7/fx1E07ILwMn35xzz8nNinuoPBR41wiM3nCe+9J9+ZwsO19RUhGhXTrJkJO97D87WZuHVGVQf8t4xENkHmFqFMJL4L1dVBCEc8ayHtNF5/nGwOgNReFG+TPPS1TqzgnuX71sbTcXLSGV6fT2DDrhgd5shNk1Anab58qPnis/Buw2Rf6GmKK/SpRxv41x/+hgT/ez+QlUZfpjOAREEEJ4CbO3ixOC8IExln7n+UbPHz9J4xNlK5+Nzlor0ep+GuZqe9i6PvbpxdQ2A/1KGi2fy+kjQXWNwOiN6fwbH5Vl58tyCqPTc8YZyal+9uB/WF+QVGXQVRn0ksf0bRVBCOElqN4ucgjCh8WyvoF277WfPX/8RAJM1PLVshV5AcPTX5rpY93s72PcK2l0ZTqdP7svEhNC1wjYrV7Tz+62HwIj5ujMNfI/FUUtf3rQLTFauGYz972VfXkpvWM5/eySx2QK9G4QhBBeQujtEIIgfIRujZCo9BzZimcc+rUnzPJj3awzQEoSqXWJVLGWMqjuGiE8dg3GNuhp+8HTdp71emTZ+bKstc6U7O9s1Lc3OKOFczHc+mT6+WTqZR0tkotVzV8QHjt27Kuvvurp6dm1a1dlZeWtP7p+/Xp1dXVra+sTTzzxySefFBcX3/ZYBOHChSAUDgRhOAQnHV7TJY/pkr/PFKl7KjrrWWty7vde7flhcs7KBTkSSsTiROqpRdStqTifXYMLBgLWAf9Ql/96j6+rjXCsLKeQyiq8FJnxrYUzWrieCa5IS5Um0+uTqCw19ZgPAO8wf9snhoeHi4uLBwYGbDbbbT+qrq7W6/X19fUNDQ0VFRVDQ0NyufxhngsAYH5IVOqY/Bdj8l/k/F5vx29e0yWF8d9/8brKE3URSWmT6rQ2ie5rq+7AH9ETfq5IS4dyMVsd5qxhg4x1wG/u8Zu7GHM3YxuSahIjUzJIcrrlqRe/Dj5pvMH+8gP3p3i2NIn++1r6mSVUxOO1/iV8HsHUaHl5eUFBwb59+2buMZvNer3eYrFoNBpCSG5u7jvvvLNly5ZbH4UR4cKFEaFwYEQ4b1i3k7H0MZZ+xtLPWPoY25BEGRdMSBtU6H6R6M4wab9zCavU7HJ1RKqCSokhKTGUTkm0MuqBv52TY/yBMYt/qIsxd/vN3Yx1gFJrXQnptjhDp8Lwm/TJDldE7yQZ8XLLY6k/J1GlyfQ6LaUQx8znffG/ob6rqysxMTGUgoSQ7Ozszs7OcDwRAMD8oOWKKEN2lCF7+jbLBuxW5kaf0tK3zPLdJutAwDnhVCW7I5U3JTE3idzERX8blI1yMkqmkMllKoU8VimPj5VpFimS4mKig27KPememCSeSdozRbunaM8k5ZnkXFOUZ0LimZJ4JqlgwB2bZFEbOuSGVk3Rd/FpDi5ar6T0CsqgIrkq6hUVpVeSFIWIvg0gTMIShA6HQ6FQzNyMjY0dGxu7rY3L5UpNTb3tzosXL2ZlZd35C0MjwmAw+MhLhQfgcrko0Z1xEKjQiJBhGL4LESVZLDHk0obcKEKiCOF8bslgdwLFEZ+H87k5r5vze4M+h9Nl8Trdvgk36/VyjEfq85CA20FHTUhVrkiVU6qcilBORqimpIlTERnMIqV7idIZEeuMVHglssVRJE3Bpim4LUqyX8ElRPvurMLjmv9XvgDMjAjlcjlN32eOeG5BaDKZXnjhBUJIXFycyWS6WzO1Wj01NTVzc2JiQqfT3dYmJiamra3ttqnRu729YmpUUDiOu/WDDvAIU6MColCQSNmdZw0SeClG9KRSabimRjMzM9vb2wkh9w5Yg8Fgs9nsdnt8fDwhxGQylZaW3tmMoigMLAAAgF9zW1RE07RKpVKpVKEBwejoaF9fn8vlcjgcfX19U1NTFy5c2LNnz9KlS0tKSj788EO/319fX9/f319RURGe+gEAAB7KQ50jPHz4cENDAyHEaDQajcb333+foqiOjg5CyOeff/7GG29oNJqUlJT6+vqYmJhHUy8AAMAjhW+WgTnD9gnhwDlCQUHXEI45bZ9YGPstz549e+jQIb6rgGlVVVVOp5PvKoAQQk6fPn3ixAm+qwBCCAkGg6+99hrfVcC0I0eONDU1zbLxwrgwr9lsvsciVZhnRqMR6/UFore3Fx9KBIJl2W+++YbvKmBae3u7x+OZZeOFMSIEAAAIEwQhAACIGm+LZTZu3Hj58uX7bvgP8fl8gUAAS08F4ubNm7GxsdgDKgRer5fjOKwjE4jx8fG4uDi+qwBCCHG73RKJJCoqqqWlJTMz896NeQvC0O5DXp4aAABEQqvV3nftKG9BCAAAIAQ4RwgAAKKGIAQAAFFbAEHY3NxcVlZWWFj46aef8l2L2J08eXLv3r2vvvoqLjDJuwsXLlRXVxcUFGzYsOHo0aN8lyNqTqdz165dJSUl+fn5b7755tDQEN8VASGE7NmzZ/fu3bNpKfQN9deuXXv99dePHDmSmJi4detWlUpVWVnJd1EixXFcTU3N6tWrm5ubd+/evWzZMr4rErXm5ubMzMyqqiqr1frWW29xHLd9+3a+ixKpQCBgMBg2bdoklUoPHz5cVlbW0dGBZdX8qqmpaWhokEpnlXFCXyzz9ttvBwKBzz77jBBy9OjRgwcPtra28l2U2Gk0msbGxoKCAr4LgWn79+/v6uqqq6vjuxAgw8PDWq12bGwsdBE64IXVai0tLd23b99HH300m+kroU+NXr16NS8vL3Scl5d35coVgSc3wPy7evWqwWDguwqxGxoaam9v//jjj5977jm1Ws13OaK2c+fODz74YNGiRbNsL/Sp0ZGRkZkXo1ar/X7/+Pg4/skAZtTW1v766681NTV8FyJ227Zt6+3t9fv9dXV1mBfl0RdffMEwTEVFRWNj4ywfIvQgVCqVbrc7dOx0OmmaxlVOAGY0NTXt3bu3paUFnw55d+7cOULI2bNnX3rppc7OzuTkZL4rEiOHw/Hee++F/hazJ/Qg1Ol0PT09oeOenp7k5ORZXl8K4LHX0tKyY8eOM2fO5OTk8F0LTCsrK4uLizOZTAhCXgwMDNjt9qysLEIIwzAej0etVnd1dWk0mns8SujnCDdv3nz8+PGJiQmO4w4dOrR582a+KwIQBKPRuHXr1lOnTq1Zs4bvWsRucHBwZGQkdNzU1GS321euXMlvSaK1atUqx/8cP35cr9c7HI57pyARfhCWl5evX79er9enpqZOTk6+++67fFckaitWrKAoym63FxYWUhQ1M1iH+XfgwIGxsbF169ZRFEVRVFFREd8ViVdHR0dmZmZKSopWq925c2dtbW1SUhLfRcEcCH37RMj4+LjP59NqtXwXAgDwfwSDweHh4YiIiMWLF/NdC8zZwghCAACAMBH61CgAAEBYIQgBAEDUEIQAACBqCEIAABA1BCEAAIgaghAAAEQNQQgAAKKGIAQAAFFDEAIAgKghCAEAQNQQhAAAIGoIQgAAEDUEIQAAiNp/AfxO10g03C3sAAAAAElFTkSuQmCC"
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots, ForwardDiff\n",
"gr(fmt = :png);\n",
"\n",
"# operator to get the derivative of this function using AD\n",
"D(f) = x -> ForwardDiff.derivative(f, x)\n",
"\n",
"# compare slopes with AD for sin(x)\n",
"q(x) = sin(x)\n",
"x = 0.0:0.1:4.0\n",
"q_x = q.(x)\n",
"q_slopes_x = slopes(q_x, x)\n",
"\n",
"D_q_x = D(q).(x) # broadcasts AD across vector\n",
"\n",
"plot(x[1:end-1], D_q_x[1:end-1], label = \"q' with AD\")\n",
"plot!(x[1:end-1], q_slopes_x, label = \"q slopes\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider a variation where we pass a function instead of an `AbstractArray`"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(q) <: Function = true\n",
"typeof(x) <: AbstractRange = true\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"q_slopes_x[1] = 0.9983341664682815\n"
]
}
],
"source": [
"slopes(f::Function, x::AbstractRange) = diff(f.(x)) / step(x) # broadcast function\n",
"\n",
"@show typeof(q) <: Function\n",
"@show typeof(x) <: AbstractRange\n",
"q_slopes_x = slopes(q, x) # use slopes(f::Function, x)\n",
"@show q_slopes_x[1];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, if `x` was an `AbstractArray` and not an `AbstractRange` we can no longer use a uniform step.\n",
"\n",
"For this, we add in a version calculating slopes with forward first-differences"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"typeof(x_array) <: AbstractArray = true\n",
"q_slopes_x[1] = "
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.9983341664682815\n"
]
}
],
"source": [
"# broadcasts over the diff\n",
"slopes(f::Function, x::AbstractArray) = diff(f.(x)) ./ diff(x)\n",
"\n",
"x_array = Array(x) # convert range to array\n",
"@show typeof(x_array) <: AbstractArray\n",
"q_slopes_x = slopes(q, x_array)\n",
"@show q_slopes_x[1];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the final example, we see that it is able to use specialized implementations over both the `f` and the `x` arguments.\n",
"\n",
"This is the “multiple” in multiple dispatch."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 1\n",
"\n",
"Explore the package [StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl).\n",
"\n",
"- Describe two abstract types and the hierarchy of three different concrete types. \n",
"- Benchmark the calculation of some simple linear algebra with a static array\n",
" compared to the following for a dense array for `N = 3` and `N = 15`. "
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 82.741 ns (1 allocation: 112 bytes)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 743.714 ns (5 allocations: 1.98 KiB)\n"
]
},
{
"data": {
"text/plain": [
"3×3 Array{Float64,2}:\n",
" -26.7149 -0.887747 22.6895\n",
" 23.9875 1.91814 -20.8351\n",
" 14.9191 0.439589 -11.2581"
]
},
"execution_count": 53,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using BenchmarkTools\n",
"\n",
"N = 3\n",
"A = rand(N, N)\n",
"x = rand(N)\n",
"\n",
"@btime $A * $x # the $ in front of variable names is sometimes important\n",
"@btime inv($A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 2\n",
"\n",
"A key step in the calculation of the Kalman Filter is calculation of the Kalman gain, as can be seen with the following example using dense matrices from [the Kalman lecture](../tools_and_techniques/kalman.html).\n",
"\n",
"Using what you learned from Exercise 1, benchmark this using Static Arrays"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 888.967 ns (10 allocations: 1.94 KiB)\n"
]
},
{
"data": {
"text/plain": [
"2×2 Array{Float64,2}:\n",
" 0.666667 1.11022e-16\n",
" 1.11022e-16 0.666667"
]
},
"execution_count": 54,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Σ = [0.4 0.3;\n",
" 0.3 0.45]\n",
"G = I\n",
"R = 0.5 * Σ\n",
"\n",
"gain(Σ, G, R) = Σ * G' * inv(G * Σ * G' + R)\n",
"@btime gain($Σ, $G, $R)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"How many times faster are static arrays in this example?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 3\n",
"\n",
"The [Polynomial.jl](https://github.com/JuliaMath/Polynomials.jl) provides a package for simple univariate Polynomials."
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"p = Polynomial(2 - 5*x + 2*x^2)\n",
"(p(0.1), p′(0.1)) = "
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"(1.52, -4.6)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"roots(p) = [0.5, 2.0]\n"
]
}
],
"source": [
"using Polynomials\n",
"\n",
"p = Polynomial([2, -5, 2], :x) # :x just gives a symbol for display\n",
"\n",
"@show p\n",
"p′ = derivative(p) # gives the derivative of p, another polynomial\n",
"@show p(0.1), p′(0.1) # call like a function\n",
"@show roots(p); # find roots such that p(x) = 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot both `p(x)` and `p′(x)` for $ x \\in [-2, 2] $."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 4\n",
"\n",
"Use your solution to Exercise 8(a/b) in [Introductory Examples](julia_by_example.html) to\n",
"create a specialized version of Newton’s method for `Polynomials` using the `derivative` function.\n",
"\n",
"The signature of the function should be `newtonsmethod(p::Polynomial, x_0; tolerance = 1E-7, maxiter = 100)`,\n",
"where `p::Polynomial` ensures that this version of the function will be used anytime a polynomial is passed (i.e. dispatch).\n",
"\n",
"Compare the results of this function to the built-in `roots(p)` function.\n",
"\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 5 (Advanced)\n",
"\n",
"The [trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule) approximates an integral with\n",
"\n",
"$$\n",
"\\int_{\\underline{x}}^{\\bar{x}} f(x) \\, dx \\approx \\sum_{n=1}^N \\frac{f(x_{n-1}) + f(x_n)}{2} \\Delta x_n\n",
"$$\n",
"\n",
"where $ x_0 = {\\underline{x}},\\, x_N = \\bar{x} $, and $ \\Delta x_n \\equiv x_{n-1} - x_n $.\n",
"\n",
"Given an `x` and a function `f`, implement a few variations of the trapezoidal rule using multiple dispatch\n",
"\n",
"- `trapezoidal(f, x)` for any `typeof(x) = AbstractArray` and `typeof(f) == AbstractArray` where `length(x) = length(f)` \n",
"- `trapezoidal(f, x)` for any `typeof(x) = AbstractRange` and `typeof(f) == AbstractArray` where `length(x) = length(f)`\n",
" * Exploit the fact that `AbstractRange` has constant step sizes to specialize the algorithm \n",
"- `trapezoidal(f, x̲, x̄, N)` where `typeof(f) = Function`, and the other arguments are `Real`\n",
" * For this, build a uniform grid with `N` points on `[x̲, x̄]` – call the `f` function at those grid points and use the existing `trapezoidal(f, x)` from the implementation \n",
"\n",
"\n",
"With these:\n",
"1. Test each variation of the function with $ f(x) = x^2 $ with $ \\underline{x}=0,\\, \\bar{x} = 1 $.\n",
"2. From the analytical solution of the function, plot the error of `trapezoidal(f, x̲, x̄, N)` relative to the analytical solution for a grid of different `N` values.\n",
"3. Consider trying different functions for $ f(x) $ and compare the solutions for various `N`.\n",
"\n",
"When trying different functions, instead of integrating by hand consider using a high-accuracy\n",
"library for numerical integration such as [QuadGK.jl](https://juliamath.github.io/QuadGK.jl/latest/)"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {
"hide-output": false
},
"outputs": [
{
"data": {
"text/plain": [
"(0.3333333333333333, 5.551115123125783e-17)"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using QuadGK\n",
"\n",
"f(x) = x^2\n",
"value, accuracy = quadgk(f, 0.0, 1.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise 6 (Advanced)\n",
"\n",
"Take a variation of your code in Exercise 5.\n",
"\n",
"Use auto-differentiation to calculate the following derivative for the example functions\n",
"\n",
"$$\n",
"\\frac{d}{d {\\bar{x}}}\\int_{\\underline{x}}^{\\bar{x}} f(x) \\, dx\n",
"$$\n",
"\n",
"Hint: See the following code for the general pattern, and be careful to\n",
"follow the [rules for generic programming](#generic-tips-tricks)."
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"hide-output": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f(0.0, 3.0) = 1.5\n",
"f(0.0, 3.1) = 1.55\n"
]
},
{
"data": {
"text/plain": [
"0.5"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using ForwardDiff\n",
"\n",
"function f(a, b; N = 50)\n",
" r = range(a, b, length=N) # one\n",
"return mean(r)\n",
"end\n",
"\n",
"Df(x) = ForwardDiff.derivative(y -> f(0.0, y), x)\n",
"\n",
"@show f(0.0, 3.0)\n",
"@show f(0.0, 3.1)\n",
"\n",
"Df(3.0)"
]
}
],
"metadata": {
"date": 1591310620.5180163,
"download_nb": 1,
"download_nb_path": "https://julia.quantecon.org/",
"filename": "introduction_to_types.rst",
"filename_with_path": "getting_started_julia/introduction_to_types",
"kernelspec": {
"display_name": "Julia 1.4.2",
"language": "julia",
"name": "julia-1.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.2"
},
"title": "Introduction to Types and Generic Programming"
},
"nbformat": 4,
"nbformat_minor": 2
}
![\"QuantEcon\"](\"https://assets.quantecon.org/img/qe-menubar-logo.svg\"){kind=logo}
\n",
" \n",
"