Julia Language

ENPC Course — Scientific Computing in Practice

Course Organization (1/2)

Session Structure

Each session will include

  • a theoretical part with slides,

  • a practical part on computer.

Assessment

  • 3 assignments (3 × 10%) as Jupyter notebooks:

    • Assignment 1 due today,

    • Assignment 2 due in the 3rd session (April 29),

    • Assignment 3 due in the 5th session (May 13).

  • Final exam on computer, without AI (70%)

Course Organization (2/2)

Detailed Schedule

Date Topic
Tuesday, April 8 Introduction to Julia (assignment 1 \(\rightarrow\) Monday April 14)
Tuesday, April 15 Interpolation and integration
Tuesday, April 29 Linear systems (assignment 2 \(\rightarrow\) Monday May 12)
Tuesday, May 6 Eigenvalue problems
Tuesday, May 13 Nonlinear systems (assignment 3 \(\rightarrow\) Monday May 26)
Tuesday, May 20 Ordinary differential equations
Tuesday, May 27 Final exam

Important Notes

  • The content of the first notebook is essential for the rest of the course.

  • The course moves quickly and requires regular work.

  • Attendance is mandatory (−2/20 per unjustified absence).

  • Programming language created in 2009 at MIT

  • By Jeff Bezanson, Stefan Karpinski, Viral B. Shah, Alan Edelman

  • Current version 1.12.5

  • General documentation for Julia: https://docs.julialang.org/en/v1/

Installing Julia → two options

  1. Version manager juliaup
  1. Individual versions from https://julialang.org/

Development Editor (IDE)

The REPL and some important commands to know

  • REPL = read–eval–print loop → Julia console

  • The REPL is accessible

    • as a standalone executable or by typing julia in a terminal (if PATH is set)

    • integrated in VSCode by pressing Shift+Ctrl+P and searching for Julia: Start REPL

      In VSCode, you can type your code lines in a *.jl file and run them one by one in the integrated console with Shift+Enter.

  • From the REPL, type

    • ] gives access to the package manager

    • ; gives access to shell mode

    • ? gives access to online help

    • ⟵ returns to normal mode

    • Tab for automatic completion

    • the up and down arrows to navigate history (optionally filtered by the beginning of the typed line)

The Package Manager

To call an already installed package

  • with import
import LinearAlgebra
LinearAlgebra.inv([4 5 ; 3 4])
2×2 Matrix{Float64}:
  4.0  -5.0
 -3.0   4.0
  • or with using (no need to prefix the package name afterwards)
using LinearAlgebra
inv([4 5 ; 3 4])
2×2 Matrix{Float64}:
  4.0  -5.0
 -3.0   4.0

To install a new package: 2 options

  1. in the package manager accessible from the REPL: type ] then
pkg> add Plots
  1. or directly in the code
import Pkg
Pkg.add("Plots")

Important

Installing a package in an environment is only done once.

Any attempt to reinstall is not only useless but also time-consuming.

Why Julia?

Python Matlab Julia C
Free/open source? ✔️ ❌ ✔️ ✔️
Fast? ❔ ❔ ✔️ ✔️
Easy to learn? ✔️ ✔️ ✔️ ❌
Suited for scientific computing? ❔ ✔️ ✔️ ❌
Rich ecosystem? ✔️ ❔ ❔ ✔️

  • Julia is easy to learn, read and maintain (comparable to Python)

  • Julia is fast (comparable to C/C++, Fortran)

  • Julia is efficient (multiple dispatch, macros, …)

  • Julia is well suited for scientific computing

Julia is easy to learn, read and maintain

  • Dynamic typing
x = 1 ; y = π ; z = 1//2 ; t = x + y + z
4.641592653589793
for v ∈ (x,y,z,t) println("$v is a $(typeof(v))") end
1 is a Int64
π is a Irrational{:π}
1//2 is a Rational{Int64}
4.641592653589793 is a Float64
  • But static typing also possible
f(x) = 2x
f(x::Float64) = 3x
@show f(1), f(1.) ;
(f(1), f(1.0)) = (2, 3.0)
  • Simplified definition of short functions (equivalent to Python’s lambda) or anonymous functions
print( (x -> 2x^2)(5) );
50

  • The standard library contains arrays and common operations (+, -, *, \), then LinearAlgebra for det, eigen…

  • Array comprehension (for inside the array)

A = [63/(i+2j) for i ∈ 1:2, j ∈ 1:2]
2×2 Matrix{Float64}:
 21.0   12.6
 15.75  10.5
x = [1.2, 3.4]
b = A*x
A\b
2-element Vector{Float64}:
 1.2000000000000008
 3.399999999999998
@show A\b == x
@show A\b ≈ x ;
A \ b == x = false
A \ b ≈ x = true
  • Indexing starts at 1 and filling is column-major
@show x[1]
@show x[2] ;
x[1] = 1.2
x[2] = 3.4
  • Unicode characters
𝒢(x,μ=0.,σ=1.) = 1/σ√(2π) * exp(-(x-μ)^2/2σ^2)
𝒢 (generic function with 3 methods)
@show Σaᵢⱼ² = sum(A.^2)
@show √Σaᵢⱼ²
ω = 1  # ⚠️ even emojis!
f(t) = sin(ω*t)
∇f(t) = ω * cos(ω*t)
f²(t) = f(t) * f(t) ;
Σaᵢⱼ² = sum(A .^ 2) = 958.0725
√Σaᵢⱼ² = 30.952746243265718
  • New binary operators
⊕(α,β) = α + β
5 ⊕ 7
12
⊕(α::Array{T,2}, β::Array{T,2}) where {T} = α + β'
⊕ (generic function with 2 methods)
display(A + A)
display(A ⊕ A)
2×2 Matrix{Float64}:
 42.0  25.2
 31.5  21.0
2×2 Matrix{Float64}:
 42.0   28.35
 28.35  21.0

Julia is fast

⚠ These results do not account for compilation time.

The vertical axis shows each benchmark time normalized against the C implementation. The benchmark data shown above were computed with Julia v1.0.0, SciLua v1.0.0-b12, Rust 1.27.0, Go 1.9, Java 1.8.0_17, Javascript V8 6.2.414.54, Matlab R2018a, Anaconda Python 3.6.3, R 3.5.0, and Octave 4.2.2. C and Fortran are compiled with gcc 7.3.1, taking the best timing from all optimization levels (-O0 through -O3). C, Fortran, Go, Julia, Lua, Python, and Octave use OpenBLAS v0.2.20 for matrix operations; Mathematica uses Intel MKL. The Python implementations of matrix_statistics and matrix_multiply use NumPy v1.14.0 and OpenBLAS v0.2.20 functions; the rest are pure Python implementations. Raw benchmark numbers in CSV format are available under https://github.com/JuliaLang/Microbenchmarks.

These micro-benchmark results were obtained on a single core (serial execution) on an Intel Core i7-3960X 3.30GHz CPU with 64GB of 1600MHz DDR3 RAM, running openSUSE LEAP 15.0 Linux.

Julia is efficient

A major advantage is multiple dispatch

Example taken from the talk The Unreasonable Effectiveness of Multiple Dispatch by Stefan Karpinski

abstract type Animal end
struct Chien <: Animal; nom::String end
struct Chat <: Animal; nom::String end

function rencontre(a::Animal,b::Animal)
    verb = agit(a,b)
    println("$(a.nom) rencontre $(b.nom) et $verb")
end

agit(a::Chien,b::Chien) = "le renifle" ; agit(a::Chien,b::Chat) = "le chasse"
agit(a::Chat,b::Chien) = "s'enfuit" ; agit(a::Chat,b::Chat) = "miaule"

medor = Chien("Médor") ; 🐶 = Chien("🐶")
felix = Chat("Félix") ; 🐱 = Chat("🐱")

rencontre(medor,🐶)
rencontre(🐶,felix)
rencontre(felix,🐶)
rencontre(🐱,felix)
Médor rencontre 🐶 et le renifle
🐶 rencontre Félix et le chasse
Félix rencontre 🐶 et s'enfuit
🐱 rencontre Félix et miaule

Example inspired by Mosè Giordano’s blog

abstract type HandShape end
struct Rock     <: HandShape end
struct Paper    <: HandShape end
struct Scissors <: HandShape end
play(::Type{Paper}, ::Type{Rock}) = println("Paper VS Rock ⇒ Paper wins")
play(::Type{Scissors}, ::Type{Paper}) = println("Scissors VS Paper ⇒ Scissors wins")
play(::Type{Rock}, ::Type{Scissors}) = println("Rock VS Scissors ⇒ Rock wins")
play(::Type{T}, ::Type{T}) where {T<: HandShape} = println("Tie between $(T), try again")
play(a::Type{<:HandShape}, b::Type{<:HandShape}) = play(b, a) # Commutativity
play(Paper, Rock)
play(Scissors, Scissors)
play(Paper, Scissors)
Paper VS Rock ⇒ Paper wins
Tie between Main.Notebook.Scissors, try again
Scissors VS Paper ⇒ Scissors wins

It is easy to add new types after the fact

struct Well <: HandShape end
play(::Type{Well}, ::Type{Rock})     = "Well VS Rock ⇒ Well wins"
play(::Type{Well}, ::Type{Scissors}) = "Well VS Scissors ⇒ Well wins"
play(::Type{Paper}, ::Type{Well})    = "Paper VS Well ⇒ Paper wins"
play(Scissors, Well)
"Well VS Scissors ⇒ Well wins"

A more mathematical example of multiple dispatch

Example inspired by the talk The Unreasonable Effectiveness of Multiple Dispatch by Stefan Karpinski

\[\sum_{i=1}^p\mathbf{\boldsymbol{v}}_i^\mathrm{T} \mathsf{A} \mathbf{\boldsymbol{v}}_i, \quad \mathsf{A}∈{\mathbb{{R}}}^{n×n},\;\mathbf{\boldsymbol{v}}_i∈{\mathbb{{R}}}^n\]

  • Standard computation
using LinearAlgebra

inner(v, A, w) = v'A*w

# function inner_sum(A, vs)
#     t = zero(eltype(A))
#     for v ∈ vs t += inner(v, A, v) end
#     return t
# end

inner_sum(A, vs) = sum(inner(v, A, v) for v ∈ vs)

n = 10_000 ; p = 100
A = rand(n, n)
vs = [rand(n) for _ ∈ 1:p]

@time inner_sum(A, vs)
  6.530119 seconds (227.91 k allocations: 18.700 MiB, 1.29% compilation time)
1.2492995156471162e9

Special case of vectors \(\mathbf{\boldsymbol{v}}_i=(0,\cdots,0,1,0,\cdots,0)^\mathrm{T}\)

\[\mathbf{\boldsymbol{v}}^\mathrm{T} \mathsf{A} \mathbf{\boldsymbol{w}} = A_{\mathrm{ind}(\mathbf{\boldsymbol{v}}), \mathrm{ind}(\mathbf{\boldsymbol{w}})}\]

  • Construction of a new vector type (size and getindex must be defined for standard matrix computations)
import Base: size, getindex
struct OneVector <: AbstractVector{Bool}
  len::Int
  ind::Int
end
size(v::OneVector) = (v.len,)
getindex(v::OneVector, i::Integer) = i == v.ind
print(OneVector(5, rand(1:5)))
Bool[0, 1, 0, 0, 0]
  • Standard computation
vs = [OneVector(n, rand(1:n)) for _ ∈ 1:p]
@time inner_sum(A, vs)
 20.919272 seconds (409.25 k allocations: 28.009 MiB, 1.45% compilation time)
48.54209890814333
  • Multiple dispatch: specialization of inner
inner(v::OneVector, A, w::OneVector) = A[v.ind, w.ind]
@time inner_sum(A, vs)
@time inner_sum(A, vs)
  0.048355 seconds (113.32 k allocations: 5.383 MiB, 99.91% compilation time: 100% of which was recompilation)
  0.000013 seconds (1 allocation: 16 bytes)
48.54209890814333

Julia is well suited for scientific computing

  • Automatic differentiation
using Zygote
f(x) = log(x)
f'(2), f''(2)
(0.5, -0.25)
  • Dimensional computations
using Unitful
m = u"m" ; cm = u"cm" ;
@show 1.0m + 1.0cm
println()
@show 1.0u"MPa" + 2.0u"bar" + 3.0u"daN/cm^2"
println()
@show uconvert(u"MPa",1u"bar")
;
1.0m + 1.0cm = 1.01 m

1.0 * u"MPa" + 2.0 * u"bar" + 3.0 * u"daN/cm^2" = 1.5e6 kg m^-1 s^-2

uconvert(u"MPa", 1 * u"bar") = 1//10 MPa
  • Differential equations
using ModelingToolkit, DifferentialEquations, Plots
@independent_variables t
∂_∂(x) = y -> expand_derivatives(Differential(x)(y))
d_dt = ∂_∂(t)
q = @variables θ(t) ϕ(t) ψ(t)
θ̇, ϕ̇, ψ̇ = q̇ = d_dt.(q)
q̈ = d_dt.(q̇)
@parameters g R M m₁ m₂ ℓ₁ ℓ₂ L ;

… Euler-Lagrange eq. \(\frac{{\mathrm{d}}}{{\mathrm{d}}t}\left(\frac{∂ℒ}{∂\dot{q_i}}\right) -\frac{∂ℒ}{∂q_i}=0\)

K = Jᴵ * θ̇^2 / 2 + m₁ * (ẋ₁^2 + ẏ₁^2) / 2 + m₂ * (ẋ₂^2 + ẏ₂^2) / 2
V = -M * g * yᴳ * cos(θ) + m₁ * g * y₁ + m₂ * g * y₂
ℒ = K - V
EL = [d_dt(∂_∂(v̇)(ℒ)) - ∂_∂(v)(ℒ) for (v̇, v) ∈ zip(q̇, q)] ;

…

  • Tensor calculus
using SymPy, LinearAlgebra, TensND
Spherical = coorsys_spherical()
θ, ϕ, r = getcoords(Spherical)
𝐞ᶿ, 𝐞ᵠ, 𝐞ʳ = unitvec(Spherical)
@set_coorsys Spherical
𝕀, 𝕁, 𝕂 = ISO() ; 𝟏 = tensId2()
k, μ = symbols("k μ", positive = true)
λ = k -2μ/3
u = SymFunction("u", real = true)
𝛆 = SYMGRAD(u(r) * 𝐞ʳ)
𝛆 |> intrinsic
(u(r)/r)𝐞ᶿ⊗𝐞ᶿ + (u(r)/r)𝐞ᵠ⊗𝐞ᵠ + (Derivative(u(r), r))𝐞ʳ⊗𝐞ʳ
𝛔 = λ * tr(𝛆) * 𝟏 + 2μ * 𝛆
eq = DIV(𝛔) ⋅ 𝐞ʳ
sol = dsolve(eq, u(r))

\(u{\left(r \right)} = \frac{C_{1}}{r^{2}} + C_{2} r\)

T̂ = tfactor(tsimplify(tsubs(𝐞ʳ ⋅ 𝛔 ⋅ 𝐞ʳ, u(r) => sol.rhs())))

\(\frac{- 4 C_{1} μ + 3 C_{2} k r^{3}}{r^{3}}\)

References for self-study