# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Julia 1.2.0
#     language: julia
#     name: julia-1.2
# ---

# # 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}

# ### Keep function argument types generic if possible

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

x = rand(3)
vander_naive(x)

vander_naive(1:3)

# The annotation `::Vector` in the the function signature is unnecessarily specific.

function vander_naive(x::AbstractVector)
    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

vander_naive(1:3)

# ### Avoid explicit typing if possible

vander_naive([1,2,3])

# Why is the result a matrix of floating point numbers....?
#
# Even worse:

vander_naive(rand(ComplexF64, 3))

# We can easily cover those cases as well by only slightly modifying our code.

function vander_almost_generic(x::AbstractVector)
    T = eltype(x)
    m = length(x)
    V = Matrix{T}(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

vander_almost_generic([1,2,3])

vander_almost_generic(rand(ComplexF64, 3))

vander_almost_generic(["Stadt", "Land", "Fluss"])

function vander_generic(x::AbstractVector{T}) where T # this is the same as just x::AbstractVector
    m = length(x)
    V = Matrix{T}(undef, m, m)
    for j = 1:m
        V[j,1] = one(T)
    end
    for i= 2:m
        for j = 1:m
            V[j,i] = x[j] * V[j,i-1]
            end
        end
    return V
end

vander_generic(["Stadt", "Land", "Fluss"])

# One more level of generality, just because we can. :)

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

function vander_supergeneric(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

vander_supergeneric([3, "Stadt", 4 + 5im])

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

using BenchmarkTools

x = rand(Float64, 100);
@btime vander_naive($x);
@btime vander_supergeneric($x);

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

x = rand(Int, 100);
@btime vander_naive($x);
@btime vander_supergeneric($x);

x = rand(Bool, 100);
@btime vander_naive($x);
@btime vander_supergeneric($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.

x = rand(Float32, 100);
@btime vander_naive($x);
@btime vander_supergeneric($x);

# # Arbitrary precision computations

# Let's say you have implemented the following crazily complex physics code as part of your thesis project which takes in a number and spits out the answer to life, the universe, and everything.

function answer_to_life_universe_and_everything(x::Integer)
    m = sin((2*x)^100)
    c = 42
    E = m*c^2
    a = sqrt(abs(E))
    b = atan(m)
    c = a^2 + b^2
    answer = sqrt(1764)/(1+exp(-c))
end

answer_to_life_universe_and_everything(2)

# Alright, apparently the answer is 21!
#
# The author: "I checked the code multiple times, it is correct. So let's publish."

# Without changing a line of code we can check the *correctness* (in the numerical sense) of our result using [arbitrary precision arithmetics](https://docs.julialang.org/en/latest/manual/integers-and-floating-point-numbers/#Arbitrary-Precision-Arithmetic-1).

big(2)

typeof(big(2))

answer_to_life_universe_and_everything(big(2))

# # Let's get a bit more fancy!

using Interact

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

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

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

# +
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!

using Colors

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

colors = distinguishable_colors(10)

colors[1]

# +
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
# -

# # 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://pretalx.com/juliacon2019/talk/BCYWZJ/)) amount of code reuse. This sharing of code comes in two forms:
# 1. **Sharing types** among a wide variety of packages implementing different algorithms;
# 2. **Sharing generic algorithms** that work for different package-defined types implementing common abstractions.

#  <img src="revdeps.png" alt="drawing" width="800"/>

# 1. **Sharing types:** DataStructures.jl, OrderedCollections.jl, StaticArrays.jl, Colors.jl, Measurements.jl ...

# As of the time of this writing, **804 packages** depend on the data types provided in [DataStructures.jl](https://juliacollections.github.io/DataStructures.jl/latest/).
#
# **851 packages** reuse type implementations in [OrderedCollections.jl](https://github.com/JuliaCollections/OrderedCollections.jl).
#
# **That's about every third package!**

# 2. **Sharing generic algorithms:** StatsBase.jl, SortingAlgorithms.jl, GenericLinearAlgebra.jl, ...

# # Emergent features

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

# +
using OrdinaryDiffEq, Measurements, PyPlot

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

#Setup
u0 = 1.0 ± 0.1
tspan = (0.0 ± 0.0, 1.0 ± 0.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)

# +
# analytic solution
u = u0 .* exp.(-c .* sol.t);

# plot solution
ts = getfield.(sol.t, :val)
solvals = getfield.(sol, :val)
solerrs = getfield.(sol, :err);

errorbar(ts, solvals, yerr=solerrs)
plot(ts, getfield.(u, :val), color="red", lw=2)
ylabel("u(t)")
xlabel("t");
# -

# 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**.

# # 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**

StarcoderdataJupyterScriptsDedupFiltered