# Who am I ? * My name is *Pierre Navaro* * **Fortran 77 + PVM** : during my PhD 1998-2002 (Université du Havre) * **Fortran 90-2003 + OpenMP-MPI** : Engineer in Strasbourg (2003-2015) at IRMA * **Numpy + Cython, R + Rcpp** : Engineer in Rennes (2015-now) at IRMAR * **Julia v1.0** since July 2018 ## Instructions to open the notebooks https://gitlab.inria.fr/navarop/JuliaInriaTech --- # Why Julia? * Started in 2009 and first version was released in 2012. * High-level languages like Python and R let one explore and experiment rapidly, but can run slow. * Low-level languages like Fortran/C++ tend to take longer to develop, but run fast. * This is sometimes called the "two language problem" and is something the Julia developers set out to eliminate. * Julia's promise is to provide a "best of both worlds" experience for programmers who need to develop novel algorithms and bring them into production environments with minimal effort. **Julia: A Fresh Approach to Numerical Computing** *Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah* SIAM Rev., 59(1), 65–98. (34 pages) 2012 --- # Julia's Engineering and Design Tradoffs * Type structures cannot be changed after being created (less dynamism but memory layout can be optimized for) * All functions are JIT compiled via LLVM (interactive lags but massive runtime improvements) * All functions specialize on types of arguments (more compilation but give generic programming structures) * Julia is interactive (use it like Python and R, but makes it harder to get binaries) * Julia has great methods for handling mutation (more optimization opportunities like C/Fortran, but more cognitive burden) * Julia's Base library and most packages are written in Julia, (you can understand the source, but learn a new package) * Julia has expensive tooling for code generation and metaprogramming (concise and more optimizations, but some codes can be for experienced users) To me, this gives me a language with a lot of depth which works well for computationally-expensive scientific applications. [© ChrisRackaukas](https://www.youtube.com/watch?v=zJ3R6vOhibA&feature=em-uploademail) --- # Type-Dispatch Programming * Centered around implementing the generic template of the algorithm not around building representations of data. * The data type choose how to efficiently implement the algorithm. * With this feature share and reuse code is very easy [JuliaCon 2019 | The Unreasonable Effectiveness of Multiple Dispatch | Stefan Karpinski](https://youtu.be/kc9HwsxE1OY) --- # Simple gravity pendulum ```julia using DifferentialEquations, Plots g = 9.79 # Gravitational constants L = 1.00 # Length of the pendulum #Initial Conditions u₀ = [0, π / 60] # Initial speed and initial angle tspan = (0.0, 6.3) # time domain #Define the problem function simplependulum(du, u, p, t) θ = u[1] dθ = u[2] du[1] = dθ du[2] = -(g/L)*θ end #Pass to solvers prob = ODEProblem(simplependulum, u₀, tspan) sol = solve(prob, Tsit5(), reltol = 1e-6) ``` --- Analytic and computed solution ```julia u = u₀[2] .* cos.(sqrt(g / L) .* sol.t) scatter(sol.t, getindex.(sol.u, 2), label = "Numerical") plot!(sol.t, u, label = "Analytic") ``` ![](pendulum1.svg) --- [Numbers with Uncertainties](http://tutorials.juliadiffeq.org/html/type_handling/02-uncertainties.html) ```julia using Measurements g = 9.79 ± 0.02; # Gravitational constants L = 1.00 ± 0.01; # Length of the pendulum #Initial Conditions u₀ = [0 ± 0, π / 60 ± 0.01] # Initial speed and initial angle #Define the problem function simplependulum(du, u, p, t) θ = u[1] dθ = u[2] du[1] = dθ du[2] = -(g/L)*θ end #Pass to solvers prob = ODEProblem(simplependulum, u₀, tspan) sol = solve(prob, Tsit5(), reltol = 1e-6); ``` --- Analytic solution ```julia u = u₀[2] .* cos.(sqrt(g / L) .* sol.t) plot(sol.t, getindex.(sol.u, 2), label = "Numerical") plot!(sol.t, u, label = "Analytic") ``` ![](pendulum2.svg) --- ## GPU Computing https://github.com/JuliaGPU ```julia using Plots, BenchmarkTools, FFTW, LinearAlgebra ``` ### Advection equation for a rotation in two dimensional domain $$ \\frac{d f}{dt} + (y \\frac{d f}{dx} - x \\frac{d f}{dy}) = 0 $$ $x \in [-π, π], y \in [-π, π]$ and $t \in [0, 200π]$ --- # Composite type for mesh definition ```julia struct Mesh nx :: Int64 ny :: Int64 x :: Vector{Float64} y :: Vector{Float64} kx :: Vector{Float64} ky :: Vector{Float64} function Mesh( xmin, xmax, nx, ymin, ymax, ny) # periodic boundary condition, we remove the end point. x = LinRange(xmin, xmax, nx+1)[1:end-1] y = LinRange(ymin, ymax, ny+1)[1:end-1] kx = 2π ./ (xmax-xmin) .* [0:nx÷2-1;nx÷2-nx:-1] ky = 2π ./ (ymax-ymin) .* [0:ny÷2-1;ny÷2-ny:-1] new( nx, ny, x, y, kx, ky) end end ``` --- # Exact computation of solution ```julia function exact!(f, time, mesh :: Mesh; shift=1.0) for (i, x) in enumerate(mesh.x), (j, y) in enumerate(mesh.y) xn = cos(time)*x - sin(time)*y - shift yn = sin(time)*x + cos(time)*y - shift f[i,j] = exp(-(xn^2+yn^2)/0.1) end end ``` ``` exact! (generic function with 1 method) ``` --- ```julia """ exact( time, mesh; shift=1.0) Computes the solution of the rotation problem """ function exact( time, mesh :: Mesh; shift=1.0) f = zeros(Float64, (mesh.nx, mesh.ny)) exact!(f, time, mesh, shift = shift) return f end ``` ``` Main.ex-index.exact ``` -- ```julia @doc exact ``` ``` exact( time, mesh; shift=1.0) ``` Computes the solution of the rotation problem --- # Create animation to show what we compute ```julia using Plots function animation( tf, nt) mesh = Mesh( -π, π, 64, -π, π, 64) dt = tf / nt t = 0 f = zeros(Float64, (mesh.nx, mesh.ny)) anim = @animate for n=1:nt exact!(f, t, mesh) t += dt p = contour(mesh.x, mesh.y, f, axis=[], framestyle=:none) plot!(p[1]; clims=(0.,1.), aspect_ratio=:equal, colorbar=false, show=false) plot!(√2 .* cos.(-π:0.1:π+0.1), √2 .* sin.(-π:0.1:π+0.1), label="", show=false) xlims!(-π,π) ylims!(-π,π) end anim end ``` ``` animation (generic function with 1 method) ``` --- ```julia anim = animation( 2π, 100) gif(anim, "rotation2d.gif", fps = 20); ``` ``` ┌ Info: Saved animation to └ fn = "/home/gitlab-runner/builds/AvRy3Hze/0/navarop/JuliaInriaTech/docs/build/rotation2d.gif" ``` ![](rotation2d.gif) --- ```julia function rotation_on_cpu( mesh :: Mesh, nt :: Int64, tf :: Float64) dt = tf / nt f = zeros(ComplexF64,(mesh.nx,mesh.ny)) exact!( f, 0.0, mesh ) exky = exp.( 1im*tan(dt/2) .* mesh.x .* mesh.ky') ekxy = exp.(-1im*sin(dt) .* mesh.y' .* mesh.kx ) for n = 1:nt fft!(f, 2) f .= exky .* f ifft!(f,2) fft!(f, 1) f .= ekxy .* f ifft!(f, 1) fft!(f, 2) f .= exky .* f ifft!(f, 2) end real(f) end ``` ``` rotation_on_cpu (generic function with 1 method) ``` --- Run the simulation and test error. ```julia mesh = Mesh( -π, π, 1024, -π, π, 1024) nt, tf = 100, 20. rotation_on_cpu(mesh, 1, 0.1) # trigger building etime = @time norm( rotation_on_cpu(mesh, nt, tf) .- exact( tf, mesh)) println(etime) ``` ``` 13.141572 seconds (496.71 k allocations: 96.259 MiB, 0.33% gc time) 4.346515881369408e-12 ``` --- # Test if GPU packages are installed ```julia using CUDAdrv GPU_ENABLED = CUDAdrv.functional() if GPU_ENABLED using CuArrays, CuArrays.CUFFT println(CUDAdrv.name(CuDevice(0))) end ``` ``` ┌ Warning: CuArrays.jl only supports CUDNN v7.6 or higher └ @ CuArrays ~/.julia/packages/CuArrays/HE8G6/src/CuArrays.jl:121 ┌ Warning: You are using CUDNN 7.0.5 for CUDA 9.0.0 with CUDA toolkit 10.2.89; these might be incompatible. └ @ CuArrays ~/.julia/packages/CuArrays/HE8G6/src/CuArrays.jl:127 Tesla K20m ``` **JuliaGPU** GPU Computing in Julia https://juliagpu.org/ --- ```julia if GPU_ENABLED function rotation_on_gpu( mesh :: Mesh, nt :: Int64, tf :: Float64) dt = tf / nt f = zeros(ComplexF64,(mesh.nx, mesh.ny)) exact!( f, 0.0, mesh) d_f = CuArray(f) # allocate f and create fft plans on GPU p_x, pinv_x = plan_fft!(d_f, [1]), plan_ifft!(d_f, [1]) p_y, pinv_y = plan_fft!(d_f, [2]), plan_ifft!(d_f, [2]) d_exky = CuArray(exp.( 1im*tan(dt/2) .* mesh.x .* mesh.ky')) d_ekxy = CuArray(exp.(-1im*sin(dt) .* mesh.y' .* mesh.kx )) for n = 1:nt p_y * d_f d_f .*= d_exky pinv_y * d_f p_x * d_f d_f .*= d_ekxy pinv_x * d_f p_y * d_f d_f .*= d_exky pinv_y * d_f end real(collect(d_f)) # Transfer f from GPU to CPU end end ``` ``` rotation_on_gpu (generic function with 1 method) ``` --- ```julia if GPU_ENABLED nt, tf = 100, 20. rotation_on_gpu(mesh, 1, 0.1) etime = @time norm( rotation_on_gpu(mesh, nt, tf) .- exact( tf, mesh)) println(etime) end ``` ``` [ Info: Building the CUDAnative run-time library for your sm_35 device, this might take a while... 0.971772 seconds (15.89 k allocations: 136.546 MiB, 3.62% gc time) 4.631189950924857e-12 ``` --- # Time Series ```julia struct TimeSeries{T,N} nt :: Int nv :: Int t :: Vector{T} u :: Vector{Array{T, 1}} function TimeSeries{T,N}( nt :: Int) where {T,N} t = zeros(T, nt) u = [zeros(T, N) for i in 1:nt] nv = N new( nt, nv, t, u) end end ``` --- ## Overload `Base.length` function -- ```julia import Base:length length(ts :: TimeSeries) = ts.nt nt, nv = 100, 2 ts = TimeSeries{Float64, nv}(nt); @show length(ts) == nt ``` ``` true ``` Generate data ```julia ts.t[1] = 0.0 ts.u[1] = [0.0, 1.0] dt = 0.01 for i in 2:nt ts.t[i] = ts.t[i-1] + dt ts.u[i][1] = sin(ts.t[i]) ts.u[i][2] = cos(ts.t[i]) end ``` --- ```julia using Plots plot(ts.t, vcat(ts.u'...)) ``` ![](plot1.svg) --- ```julia plot(ts.t, [getindex.(ts.u, i) for i in 1:nv]) ``` ![](plot2.svg) --- ## Overload the `[]` operator we want `ts[i]` equal to `ts.u[:][i]` values -- ```julia import Base: getindex ``` -- ```julia getindex( ts :: TimeSeries, i ) = getindex.(ts.u, i) ``` ``` getindex (generic function with 395 methods) ``` --- ```julia plot(ts[1], ts[2]) ``` ![](plot3.svg) --- ## Overload the `+` operator to add noise ```julia import Base:+ ``` -- ```julia function +(ts :: TimeSeries, ϵ ) for n in 1:ts.nt, d in 1:ts.nv ts.u[n][d] += ϵ[n,d] end return ts end ``` ``` + (generic function with 267 methods) ``` --- ```julia ts = ts + 0.1*randn((nt,nv)); ``` -- ```julia scatter(ts.t, [ts[1],ts[2]]) ``` ![](plot4.svg) --- # Linear regression with obvious operation ```julia using LinearAlgebra X = hcat(ones(nt), ts.t, ts[1]) y = ts[2] @show β = inv(X'X) * X'y ``` ``` 3-element Array{Float64,1}: 1.084666295350754 -0.47803032839141224 0.002909660504045064 ``` --- # Version with QR factorisation ```julia @show β = X \ y ``` ``` 3-element Array{Float64,1}: 1.0846662953507533 -0.4780303283914258 0.0029096605040541063 ``` The `\` operator is the short-hand for ```julia Q, R = qr(X) @show β = inv(factorize(R))Q'y ``` ``` 3-element Array{Float64,1}: 1.0846662953507542 -0.4780303283914261 0.0029096605040542997 ``` --- # Version with singular values decomposition ```julia U, S, V = svd(X) @show β = V * diagm(1 ./ S) * U' * y ``` ``` 3-element Array{Float64,1}: 1.0846662953507535 -0.47803032839142634 0.0029096605040541193 ``` -- ```julia @show β = pinv(X, atol=1e-6) * y ``` ``` 3-element Array{Float64,1}: 1.0846662953507538 -0.47803032839142656 0.002909660504054175 ``` --- ## With GLM.jl ```julia using GLM fitted = lm(X, y) ``` ``` GLM.LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}}: Coefficients: ─────────────────────────────────────────────────────────────────────── Estimate Std. Error t value Pr(>|t|) Lower 95% Upper 95% ─────────────────────────────────────────────────────────────────────── x1 1.08467 0.0211186 51.3607 <1e-71 1.04275 1.12658 x2 -0.47803 0.106364 -4.49431 <1e-4 -0.689132 -0.266928 x3 0.00290966 0.113689 0.0255932 0.9796 -0.222731 0.228551 ─────────────────────────────────────────────────────────────────────── ``` --- # Metaprogramming ## The ParticleGroup example ```julia """ ParticleGroup{D,V}(n_particles, charge, mass) - `D` : number of dimension in physical space - `V` : number of dimension in phase space - `n` : number of particles """ struct ParticleGroup{D,V} n_particles :: Int64 data :: Array{Float64, 2} function ParticleGroup{D,V}(n) where {D, V} data = zeros( Float64, (D+V, n)) new( n, data) end end ``` ``` Main.ex-index.ParticleGroup ``` --- Set position of ith particle of p to x ```julia @generated function set_x!( p :: ParticleGroup{D,V}, i, x :: Float64 ) where {D, V} :(p.data[1, i] = x) end ``` ``` set_x! (generic function with 1 method) ``` -- Set position of ith particle of p to x when x is a vector ```julia @generated function set_x!( p :: ParticleGroup{D,V}, i, x :: Vector{Float64} ) where {D, V} :(for j in 1:$D p.data[j, i] = x[j] end) end ``` ``` set_x! (generic function with 2 methods) ``` --- Set velocity of ith particle of p to v ```julia @generated function set_v!( p :: ParticleGroup{D,V}, i, v :: Float64 ) where {D, V} :(p.data[$D+1, i] = v) end ``` ``` set_v! (generic function with 1 method) ``` -- Set velocity of ith particle of p to v ```julia @generated function set_v!( p :: ParticleGroup{D,V}, i, v :: Vector{Float64} ) where {D, V} :(for j in 1:$V p.data[$D+j, i] = v[j] end) end ``` ``` set_v! (generic function with 2 methods) ``` --- Get position of ith particle of p ```julia @generated function get_x( p :: ParticleGroup{D,V}, i ) where {D, V} :(p.data[1:$D, i]) end ``` ``` get_x (generic function with 1 method) ``` Get velocity of ith particle of p ```julia @generated function get_v( p :: ParticleGroup{D,V}, i ) where {D, V} :(p.data[$D+1:$D+$V, i]) end ``` ``` get_v (generic function with 1 method) ``` --- ```julia import Sobol ``` f(x,v) = 1/2π (1 + α cos (kx * x) * exp(-(vx^2+vy^2)) ```julia function landau_sampling!( pg :: ParticleGroup{1,2}, alpha, kx ) function newton(r) x0, x1 = 0.0, 1.0 r *= 2π / kx while (abs(x1-x0) > 1e-12) p = x0 + alpha * sin( kx * x0) / kx f = 1 + alpha * cos( kx * x0) x0, x1 = x1, x0 - (p - r) / f end x1 end s = Sobol.SobolSeq(2) nbpart = pg.n_particles for i=1:nbpart v = sqrt(-2 * log( (i-0.5)/nbpart)) r1, r2 = Sobol.next!(s) θ = r1 * 2π set_x!(pg, i, newton(r2)) set_v!(pg, i, [v * cos(θ), v * sin(θ)]) end end ``` ``` landau_sampling! (generic function with 1 method) ``` --- ```julia n_particles = 10000 pg = ParticleGroup{1,2}( n_particles) alpha, kx = 0.1, 0.5 landau_sampling!(pg, alpha, kx) ``` -- ```julia xp = vcat([get_x(pg, i) for i in 1:pg.n_particles]...) vp = vcat([get_v(pg, i) for i in 1:pg.n_particles]'...) ``` ``` 10000×2 Array{Float64,2}: -4.4505 5.45029e-16 -7.70866e-16 -4.1964 2.4939e-16 4.07285 -2.82092 2.82092 2.77602 -2.77602 -2.73963 -2.73963 2.70897 2.70897 1.45172 3.50477 -1.43904 -3.47415 3.44671 -1.42768 ⋮ 0.0146796 0.0385387 -0.0137887 -0.0361999 0.0336995 -0.0128363 -0.0309982 0.0118074 -0.0122742 0.0273778 0.0108246 -0.0241444 -0.0204052 -0.0091482 0.0158054 0.00708599 0.00815159 0.00579259 ``` --- ```julia using Plots pp = plot(layout=(3,1)) histogram!(pp[1,1], xp, normalize=true, bins = 100, lab=:x) plot!(pp[1,1], x -> (1+alpha*cos(kx*x))/(2π/kx), 0., 2π/kx, lab="") histogram!(pp[2,1], vp[:,1], normalize=true, bins = 100, lab=:vx) plot!(pp[2,1], v -> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="") histogram!(pp[3,1], vp[:,2], normalize=true, bins = 100, lab=:vy) plot!(pp[3,1], v -> exp( - v^2 / 2) * 4 / π^2 , -6, 6, lab="") ``` ![](particles.svg) --- ```julia histogram2d(vp[:,1], vp[:,2], normalize=true, bins=100) ``` ![](hist2d.svg) --- ### Optimizing Julia code is often done at the expense of transparency ```julia using Random, LinearAlgebra, BenchmarkTools A = rand(1024, 256); B = rand(256, 1024); C = rand(1024, 1024) function test1(A, B, C) C = C - A * B return C end @btime test1(A, B, C); #C, A and B are matrices. ``` ``` 5.128 ms (4 allocations: 16.00 MiB) ``` -- ```julia function test2(A, B, C) C .-= A * B return C end @btime test2(A, B, C); #C, A and B are matrices. ``` ``` 4.608 ms (2 allocations: 8.00 MiB) ``` --- ```julia function test_opt(A, B, C) BLAS.gemm!('N','N', -1., A, B, 1., C) return C end @btime test_opt(A, B, C) # avoids taking two unnecessary copies of the matrix C. ``` ``` 2.250 ms (0 allocations: 0 bytes) ``` -- ```julia C = rand(1024, 1024) all(test1(A, B, C) .== test2(A, B, C)) ``` ``` true ``` -- ```julia C = rand(1024, 1024) all(test1(A, B, C) .== test_opt(A, B, C)) ``` ``` true ``` --- ### Derivative computation with FFT ```julia using FFTW xmin, xmax, nx = 0, 4π, 1024 ymin, ymax, ny = 0, 4π, 1024 x = LinRange(xmin, xmax, nx+1)[1:end-1] y = LinRange(ymin, ymax, ny+1)[1:end-1] ky = 2π ./ (ymax-ymin) .* [0:ny÷2-1;ny÷2-ny:-1] exky = exp.( 1im .* ky' .* x) function df_dy( f, exky ) ifft(exky .* fft(f, 2), 2) end f = sin.(x) .* cos.(y') # f is a 2d array created by broadcasting @btime df_dy(f, exky) ``` ``` 84.282 ms (118 allocations: 64.01 MiB) ``` --- ### Memory alignement, and inplace computation. ```julia f = zeros(ComplexF64, (nx,ny)) fᵗ = zeros(ComplexF64, reverse(size(f))) f̂ᵗ = zeros(ComplexF64, reverse(size(f))) f .= sin.(x) .* cos.(y') fft_plan = plan_fft(fᵗ, 1, flags=FFTW.PATIENT) function df_dy!( f, fᵗ, f̂ᵗ, exky ) transpose!(fᵗ,f) mul!(f̂ᵗ, fft_plan, fᵗ) f̂ᵗ .= f̂ᵗ .* exky ldiv!(fᵗ, fft_plan, f̂ᵗ) transpose!(f, fᵗ) end @btime df_dy!(f, fᵗ, f̂ᵗ, exky ) ``` ``` 33.638 ms (3 allocations: 176 bytes) ``` --- # Why use Julia language! * **You develop in the same language in which you optimize.** * Packaging system is very efficient (3173 registered packages) * PyPi (198,360 projects) R (14993 packages) * It is very easy to create a package (easier than R and Python) * It is very easy to use GPU device. * Nice interface for Linear Algebra and Differential Equations * Easy access to BLAS and LAPACK * Julia talks to all major Languages - mostly without overhead! --- # What's bad * It is still hard to build shared library or executable from Julia code. * Compilation times can be unbearable. * Plotting takes time (20 seconds for the first plot) * OpenMP is better than the Julia multithreading library but it is progressing. * For parallelization, The Julia community seems to prefer the distributed processing approach. * Does not work well with vectorized code, you need to do a lot of inplace computation to avoid memory allocations and use explicit views to avoid copy. * Julia website proclaims that it is faster than Fortran but this is not true. But it is very close and it is progressing. [What's Bad About Julia by Jeff Bezanson](https://www.youtube.com/watch?v=TPuJsgyu87U) --- ## So when should i use Julia? * Now! If you need performance and plan to write your own libraries. * In ~1-2 Years if you want a smooth deploy. * In ~3-5 Years if you want a 100% smooth experience. ## Think Julia: How to Think Like a Computer Scientist https://github.com/BenLauwens/ThinkJulia.jl --- ## Python-Julia benchmarks by Thierry Dumont https://github.com/Thierry-Dumont/BenchmarksPythonJuliaAndCo/wiki ## Mailing Lists * Rennes : https://listes.univ-rennes1.fr/wws/info/math-julia * France : https://listes.services.cnrs.fr/wws/info/julia * World : https://discourse.julialang.org --- # Julia is a language made for Science. [Some State of the Art Packages](http://www.stochasticlifestyle.com/some-state-of-the-art-packages-in-julia-v1-0) * JuliaDiffEq – Differential equation solving and analysis. * JuliaDiff – Differentiation tools. * JuliaGeometry – Computational Geometry. * JuliaGraphs – Graph Theory and Implementation. * JuliaIntervals - Rigorous numerics with interval arithmetic & applications. * JuliaMath – Mathematics made easy in Julia. * JuliaOpt – Optimization. * JuliaPolyhedra – Polyhedral computation. * JuliaSparse – Sparse matrix solvers. * JuliaStats – Statistics and Machine Learning. * JuliaPlots - powerful convenience for visualization. * JuliaGPU - GPU Computing for Julia. * FluxML - The Elegant Machine Learning Stack. --- *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*