Mathieu Besançon
Zuse Institute Berlin
Run the notebooks (local or online): https://github.com/matbesancon/Julia_JuMP_workshop
"but we already have a programming language at home"
Technical consequence: workflow interruption
Data exploration $\Rightarrow$ pre-processing $\Rightarrow$ computation $\Rightarrow$ post-processing $\Rightarrow$ visualization
Human consequence: culture gap between science and development
# scalars
x = 3
y = 2x
6
# Vector of integers
v = [3, 4, 5]
3-element Array{Int64,1}: 3 4 5
d = 3.0 # Float64
f = 3.14f0 # Float32
3.14f0
π
π = 3.1415926535897...
# complex numbers
z = 2.0 + 3im
2.0 + 3.0im
typeof(z)
Complex{Float64}
# arbitrary precision
b = big"3.4020496842"
3.402049684199999999999999999999999999999999999999999999999999999999999999999983
typeof(b)
BigFloat
# matrices and multiplication
M = [
0 0 1
0 1 0
1 0 0
]
3×3 Array{Int64,2}: 0 0 1 0 1 0 1 0 0
M * v
3-element Array{Int64,1}: 5 4 3
1.5 * M * v
3-element Array{Float64,1}: 7.5 6.0 4.5
M * M * v
3-element Array{Int64,1}: 3 4 5
# N-dimensional arrays
x_scalar = rand()
x_vec = rand(3)
x_matrix = rand(3, 2)
x_tensor = rand(3,2,5)
size(x_tensor)
(3, 2, 5)
vec_int = [3, 4, 5] # Array{Int64,1}
vec_float = [3.0, 4, 5] # Array{Float64,1}
vec_comp = [3.0, 4 + 2im, 5] # Array{Complex{Float64},1}
vec_comp = [3, 4, "hello"]; # Array{Any,1}
Central building block in Julia.
Multiple dispatch $\Rightarrow$ secret ingredient for programs looking like maths.
function f1(x, y)
if x < y
return x * y
end
return y / x
end
f2(x, y) = 3x + 5y
f2 (generic function with 1 method)
# can be typed with \rightdotarrow<tab>
⤑(a, b) = b[a]
⤑ (generic function with 1 method)
3 ⤑ [1, 2, 42, 10]
42
"""
Given an exponent `α`, returns a function `x -> x^α`
"""
function exponentiate(α)
function(x)
x^α
end
end
exponentiate
square = exponentiate(2)
square(3)
9
?exponentiate
search: exponentiate ExponentialBackOff
Given an exponent α
, returns a function x -> x^α
?exponenti
search: exponentiate ExponentialBackOff exponent Couldn't find exponenti Perhaps you meant exponent, exponentiate or export
No documentation found.
Binding exponenti
does not exist.
using LinearAlgebra: dot
function quadratic(x::AbstractVector)
return dot(x, x)
end
quadratic (generic function with 1 method)
quadratic([1, 1])
2
function quadratic(x::AbstractVector, y::AbstractVector)
return dot(x, y)
end
quadratic (generic function with 2 methods)
quadratic([1, 2], [1, 1])
3
function quadratic(A::AbstractMatrix, x::AbstractVector)
return x' * A * x
end
quadratic (generic function with 3 methods)
M = rand(3, 3)
quadratic(M, [1, 0, 1])
1.7776375865298464
Look at Stefan Karpinski's talk at Juliacon 2019 for comparison with other systems.
"""
A wrapper type for an Int
"""
struct IntWrapper
x::Int
end
IntWrapper
v = IntWrapper(3)
v.x
3
import Base: +
# we define + for our new type and an integer
+(iw::IntWrapper, v::Integer) = IntWrapper(iw.x + v)
+ (generic function with 185 methods)
v + 33
IntWrapper(36)
+(iw::IntWrapper, s::String) = IntWrapper(iw.x + length(s))
+ (generic function with 186 methods)
v + "hello"
IntWrapper(8)
Brings best software engineering practice to scientists. Source: https://unsplash.com/photos/IClZBVw5W5A
(from julialang.org)
"Basic unit" of Julia code.
├── LICENSE
├── Manifest.toml
├── Project.toml
├── README.md
├── src
│ ├── BilevelOptimization.jl
│ ├── file1.jl
│ └── file2.jl
└── test
└── runtests.jl
Project.toml defines the project name, dependencies...
───────┬───────────────────────────────────────────────
│ File: Project.toml
───────┼───────────────────────────────────────────────
1 │ name = "BilevelOptimization"
2 │ uuid = "98803d92-2a2a-5e5c-9642-fb423c87776e"
3 │ authors = ["Mathieu Besançon <mathieu.besancon@gmail.com>"]
4 │ version = "0.2.2"
5 │
6 │ [deps]
7 │ JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
8 │ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
9 │ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
10 │
11 │ [compat]
12 │ JuMP = "0.21"
13 │ julia = "1"
14 │
15 │ [extras]
16 │ Cbc = "9961bab8-2fa3-5c5a-9d89-47fab24efd76"
17 │ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
18 │
19 │ [targets]
20 │ test = ["Test", "Cbc"]
Manifest.toml shows how the project was run (versions for all dependencies).
Useful to reproduce exactly an environment.
───────┬───────────────────────────────────
│ File: Manifest.toml
───────┼───────────────────────────────────
1 │ # This file is machine-generated - editing it directly is not advised
2 │
3 │ [[Base64]]
4 │ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
5 │
6 │ [[BenchmarkTools]]
7 │ deps = ["JSON", "Printf", "Statistics"]
8 │ git-tree-sha1 = "90b73db83791c5f83155016dd1cc1f684d4e1361"
9 │ uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
10 │ version = "0.4.3"
Source: Cormullion
Some cool tricks to convince you.
Multiple dispatch $\rightarrow$ functions can be extended by a 3rd party.
# package 1 develops a new number-like thing
struct Zero
end
import Base: *, /
*(::Zero, x) = Zero()
*(x, ::Zero) = Zero()
+(::Zero, x) = x
+(x, ::Zero) = x
/(::Zero, x) = Zero()
/ (generic function with 119 methods)
# package 2 develops an algorithm using a "number-like" object
function complex_algorithm(x, y)
return x + 2y + x/y
end
# user plugs the new type when using the algorithm
@show complex_algorithm(Zero(), 3.0);
@show complex_algorithm(0, 3.0);
complex_algorithm(Zero(), 3.0) = 6.0 complex_algorithm(0, 3.0) = 6.0
# symbolic maths engine
using SymEngine
@vars a, b
@show complex_algorithm(a, b);
@show complex_algorithm(a, π);
complex_algorithm(a, b) = a + 2*b + a/b complex_algorithm(a, π) = 6.28318530717959 + a + a/pi
Differentiation of almost arbitrary programs with respect to their input.
using ForwardDiff
function sqrt_babylonian(s)
x = s / 2
while abs(x^2 - s) > 0.001
x = (x + s/x) / 2
end
x
end
sqrt_babylonian (generic function with 1 method)
sqrt_babylonian(2) - sqrt(2)
2.123901414519125e-6
@show ForwardDiff.derivative(sqrt_babylonian, 2);
@show ForwardDiff.derivative(sqrt, 2);
ForwardDiff.derivative(sqrt_babylonian, 2) = 0.353541906958862 ForwardDiff.derivative(sqrt, 2) = 0.35355339059327373
Physicists' dreams finally made true.
using Unitful
using Unitful: J, kg, m, s
3J + 1kg * (1m / 1s)^2
4.0 kg m^2 s^-2
complex_algorithm(3kg, 2s)
DimensionError: 3 kg and 4 s are not dimensionally compatible. Stacktrace: [1] +(::Quantity{Int64,𝐌,Unitful.FreeUnits{(kg,),𝐌,nothing}}, ::Quantity{Int64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}}) at /home/matbesancon/.julia/packages/Unitful/KE9TK/src/quantities.jl:105 [2] +(::Quantity{Int64,𝐌,Unitful.FreeUnits{(kg,),𝐌,nothing}}, ::Quantity{Int64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}}, ::Quantity{Float64,𝐌*𝐓^-1,Unitful.FreeUnits{(kg, s^-1),𝐌*𝐓^-1,nothing}}) at ./operators.jl:538 [3] complex_algorithm(::Quantity{Int64,𝐌,Unitful.FreeUnits{(kg,),𝐌,nothing}}, ::Quantity{Int64,𝐓,Unitful.FreeUnits{(s,),𝐓,nothing}}) at ./In[35]:4 [4] top-level scope at In[42]:1
Getting help: https://docs.julialang.org/
Reaching out: https://julialang.org/community
Relevant links:
Is your research software correct: https://mikecroucher.github.io/MLPM_talk/