# 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 (revisited)

[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))
```

* Keep function argument types generic if possible
* Avoid explicit typing if possible


```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
```


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

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


```julia
using BenchmarkTools
```


```julia
x = rand(Float64, 100);
@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(Int, 100);
@btime vander_naive($x);
@btime vander_generic($x);
```


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

On the other hand, sometimes it is worth converting to a different type to dispatch to a faster method or to utilize magic like compiler optimizations.


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

# Generic Programming + Multiple Dispatch + JIT

The possibility to write generic algorithms that compile to fast machine code 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 generic algorithms**
 * Examples: StatsBase.jl, SortingAlgorithms.jl, GenericLinearAlgebra.jl,
 * **Methods are selected based on all argument types**
2. **Sharing types**
 * Examples: DataStructures.jl, OrderedCollections.jl, StaticArrays.jl, Colors.jl, Measurements.jl
 * **You can define methods on types after the type is defined!**

 

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)

# Let's get a bit more fancy!


```julia
using Interact
```


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


```julia
function insert_block(A::AbstractMatrix, i, j, what=7)
    B = copy(A)
    B[i:i+2, j:j+2] .= what
    B
end
```


```julia
A = fill(0, 9, 9)
insert_block(A, 3, 5)  # this returns the new matrix
```


```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)
end
```

### Let's add some color!

Our function `insert_block` is generic. Since the first argument `A isa AbstractArray`, we can index into it and set new values. Pretty much every value type is fine!


```julia
using Colors
```


```julia
@manipulate for n in 1:80
    distinguishable_colors(n)
end
```


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


```julia
colors[1]
```


```julia
A = fill(colors[1], 10, 10)
n = size(A, 1)

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

## Emergent features: Example DifferentiaEquations.jl

$$\frac{du}{dt} = f(u,p,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.


```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])
```

# 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 + JIT = **lots of code sharing and emergent features**

c4557548771c9786c7d7ab8507a6d95ad63ac8fa
ipynb
Jupyter Notebook
Day1/5_generic_programming.ipynb
michaellass/JuliaNRW21
337eacb40950ea37c5e0de485fe7ecb9192efb68
MIT
2021-02-24T05:08:46.000Z
2021-05-28T11:13:52.000Z
Day1/5_generic_programming.ipynb
michaellass/JuliaNRW21
337eacb40950ea37c5e0de485fe7ecb9192efb68
MIT
2021-03-02T19:10:00.000Z
2021-03-03T20:28:17.000Z
Day1/5_generic_programming.ipynb
michaellass/JuliaNRW21
337eacb40950ea37c5e0de485fe7ecb9192efb68
MIT
2021-03-01T07:14:26.000Z
2021-06-21T08:50:25.000Z
Qwen/Qwen-72B
1. YES
2. YES