# Abstraction

We have seen **duck typing** as a convient abstraction tool: 
* As a **user**, we don't have to care about how a specific type, say an `Array`, is implemented. By *being an array* (i.e. a subtype of `AbstractArray`) it behaves like we expect and we can just use it.
* As a **developer**, as long as we make our objects behave like, say, and `AbstractArray`, we can implement it in whatever way we deem appropriate and it will work with all kinds of algorithms.

Building upon this principle, we also want to formulate our **algorithms** in an abstract way such that it works will all kinds of data types irrespective of their precise implementation (or even meaning). This is generally known as **generic programming**.

From [Wikipedia](https://en.wikipedia.org/wiki/Generic_programming):
> **Generic programming** is a style of computer programming in which algorithms are written in terms of types *to-be-specified-later* that are then *instantiated* when needed for specific types provided as parameters.

## Example: Vandermonde matrix

[Vandermonde matrix:](https://en.wikipedia.org/wiki/Vandermonde_matrix)
\begin{align}V=\begin{bmatrix}1&\alpha _{1}&\alpha _{1}^{2}&\dots &\alpha _{1}^{n-1}\\1&\alpha _{2}&\alpha _{2}^{2}&\dots &\alpha _{2}^{n-1}\\1&\alpha _{3}&\alpha _{3}^{2}&\dots &\alpha _{3}^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&\alpha _{m}&\alpha _{m}^{2}&\dots &\alpha _{m}^{n-1}\end{bmatrix}\end{align}


```julia
function vander_naive(x::Vector)
    m = length(x)
    V = Matrix{Float64}(undef, m, m)
    for j = 1:m
        V[j,1] = 1.0
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end
```


```julia
vander_naive(rand(3))
```


```julia
vander_naive(1:3)
```


```julia
vander_naive([1,2,3])
```

Why is the result a matrix of floating point numbers....?

Even worse:


```julia
vander_naive(rand(ComplexF64, 3))
```

Our algorithm `vander_naive` is not generic! It assumes too much about the input type (it only takes `Vector`s) and also special cases `Float64`. Instead we should try to
* keep function argument types generic if possible
* avoid explicit typing (which is rarely necessary)


```julia
function vander_generic(x::AbstractVector{T}) where T
    m = length(x)
    V = Matrix{T}(undef, m, m)
    for j = 1:m
        V[j,1] = one(x[j])
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end
```

Not only does our code now work in the obvious cases above:


```julia
vander_generic(1:3)
```


```julia
vander_generic(rand(ComplexF64, 3))
```

It also extends to data types that we probably didn't even have in mind when writing down the algorithm:


```julia
vander_generic([3, "Stadt", 4 + 5im])
```

(Note that `*` is string concatenation and `one(String) == ""`, representing the identity element under multiplication.)

### And all of this genericity comes at no performance penality


```julia
using BenchmarkTools
```


```julia
x = rand(Float64, 10);
@btime vander_naive($x);
@btime vander_generic($x);
```

Actually, for this specific example **our generic code is faster** in a few cases inasmuch as type conversions are unnecessary.


```julia
x = rand(Bool, 10);
@btime vander_naive($x);
@btime vander_generic($x);
```


```julia
x = rand(Float32, 10);
@btime vander_naive($x);
@btime vander_generic($x);
```

# Generic Programming + Multiple Dispatch

The possibility to write generic algorithms that compile to fast machine code (a point we will investigate later) in combination with multiple dispatch leads to an ([unreasonable](https://www.youtube.com/watch?v=kc9HwsxE1OY)) amount of code reuse. This sharing of code comes in two forms:
1. **Sharing types**
 * Examples: DataStructures.jl, OrderedCollections.jl, StaticArrays.jl, Colors.jl, Measurements.jl
 * **You can define methods on types after the type is defined!**
 
 
2. **Sharing generic algorithms**
 * Examples: StatsBase.jl, SortingAlgorithms.jl, GenericLinearAlgebra.jl,
 * **Methods are selected based on all argument types**

 

As of the time of this writing, **~ 800 packages** depend on the data types provided in [DataStructures.jl](https://juliacollections.github.io/DataStructures.jl/latest/).

**~ 900 packages** reuse type implementations in [OrderedCollections.jl](https://github.com/JuliaCollections/OrderedCollections.jl).

**~ 550 packages** reuse [SortingAlgorithms.jl](https://github.com/JuliaCollections/SortingAlgorithms.jl).

**~ 300 packages** reuse [DiffRules.jl](https://github.com/JuliaDiff/DiffRules.jl).

**~ 50 packages** reuse [GenericSVD.jl](GenericSVD.jl).

(~ 5000 packages overall)

# Fancy example of generic programming


```julia
using Interact
```


```julia
@manipulate for n in 1:20
    [i*j for i in 1:n, j in 1:n]
end
```


```julia
# generic algorithm
function insert_block(A, i, j, blockvalue)
    B = copy(A)
    B[i:i+2, j:j+2] .= blockvalue
    return B
end
```


```julia
# generic algorithm
function insert_block(A::AbstractMatrix{T}, i::Integer, j::Integer, blockvalue::T) where T
    B = copy(A)
    B[i:i+2, j:j+2] .= blockvalue
    return B
end
```


```julia
A = fill(0, 9, 9)
```


```julia
insert_block(A, 3, 5, 1)
```


```julia
A = fill(0, 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j, 1)
end
```

### Let's add some color!

Our function `insert_block` is generic. We don't care about the element type of `A`. Pretty much every type is fine!


```julia
using Colors
```


```julia
colors = distinguishable_colors(10)
```


```julia
basecolor = colors[1]
```


```julia
accentcolor = colors[4]
```


```julia
A = fill(basecolor, 10, 10)
```


```julia
A = fill(basecolor, 10, 10)
n = size(A, 1)

@manipulate for i in 1:n-2, j in 1:n-2
    insert_block(A, i, j, accentcolor)
end
```

## Emergent features: Example DifferentiaEquations.jl

$$\frac{du(t)}{dt} = -cu(t)$$


```julia
using OrdinaryDiffEq

#Half-life of Carbon-14 is 5730 years.
c = 5.730

#Setup
u0 = 1.0
tspan = (0.0, 1.0)

#Define the problem
radioactivedecay(u,p,t) = -c*u

#Pass to solver
prob = ODEProblem(radioactivedecay,u0,tspan)
```


```julia
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8);
```


```julia
using Plots

plot(sol.t, sol.u, ylabel="u(t)", xlabel="t", lw=2, legend=false)
```

### Arbitrary precision


```julia
using OrdinaryDiffEq, Measurements, Plots

#Half-life of Carbon-14 is 5730 years.
c = big(5.730) # now a BigFloat

#Setup
u0 = big(1.0) # now a BigFloat
tspan = (0.0, 1.0)

#Define the problem
radioactivedecay(u,p,t) = -c*u

#Pass to solver
prob = ODEProblem(radioactivedecay,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8);

plot(sol.t, sol.u, ylabel="u(t)", xlabel="t", lw=2, legend=false)
```


```julia
sol.u
```

### With uncertainties: Measurements.jl

`Measurement` type from [Measurements.jl]() and differential equation solver from [OrdinaryDiffEq.jl](https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl) (i.e. [DifferentialEquations.jl](https://github.com/JuliaDiffEq/DifferentialEquations.jl))


```julia
using Measurements
```


```julia
1 ± 0.1
```


```julia
(1 ± 0.1) * (2 ± 0.05)
```


```julia
using OrdinaryDiffEq, Measurements, Plots

#Half-life of Carbon-14 is 5730 years.
c = 5.730 ± 2

#Setup
u0 = 1.0 ± 0.1
tspan = (0.0, 1.0)

#Define the problem
radioactivedecay(u,p,t) = -c*u

#Pass to solver
prob = ODEProblem(radioactivedecay,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8);

plot(sol.t, sol.u, ylabel="u(t)", xlabel="t", lw=2, legend=false)
```

Note that, in some sense, **Julia implemented that feature by itself**.

The authors of Measurements.jl and DifferentialEquations.jl never had any collabration on this.

It **just works**.

## Symbolic computations

Let's compute the Vandermonde matrix symbolically.


\begin{align}V=\begin{bmatrix}1&\alpha _{1}&\alpha _{1}^{2}&\dots &\alpha _{1}^{n-1}\\1&\alpha _{2}&\alpha _{2}^{2}&\dots &\alpha _{2}^{n-1}\\1&\alpha _{3}&\alpha _{3}^{2}&\dots &\alpha _{3}^{n-1}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&\alpha _{m}&\alpha _{m}^{2}&\dots &\alpha _{m}^{n-1}\end{bmatrix}\end{align}


```julia
using SymPy
```


```julia
@vars a b c d e
```


```julia
sin(a)sin(b)+cos(a)*cos(b)
```


```julia
simplify(ans)
```


```julia
vander_generic([a,b,c,d,e])
```

(See [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) for a pure-Julia SymPy replacement.)

# Core messages of this Notebook

* It is simple to write **type-generic code** in Julia and you should do it.
* Generally, **generic code is just as fast as specific code**.
* Generic Programming + Multiple Dispatch = **lots of code sharing and emergent features**

410f5dcbfaa464e25e647f6574f5b39291c51e60
ipynb
Jupyter Notebook
Day1/4_abstraction.ipynb
tbdrstl/JuliaCologne21
0e2778cd8d0601291dcc8573c20469b3fb84ebe7
MIT
2021-03-09T19:05:00.000Z
2021-08-11T05:00:07.000Z
Day1/4_abstraction.ipynb
tbdrstl/JuliaCologne21
0e2778cd8d0601291dcc8573c20469b3fb84ebe7
MIT
2021-03-18T09:27:40.000Z
2021-04-08T09:37:47.000Z
Day1/4_abstraction.ipynb
tbdrstl/JuliaCologne21
0e2778cd8d0601291dcc8573c20469b3fb84ebe7
MIT
2021-03-18T09:30:02.000Z
2021-06-18T06:50:58.000Z
Qwen/Qwen-72B
1. YES
2. YES