Simon Danisch / Aug 28 2018

Julia Workshop

1.
Julia IDEs

1.1.
Atom

https://atom.io/packages/language-julia

1.2.
Visual Studio Code

https://github.com/JuliaEditorSupport/julia-vscode

1.3.
Jupyter

https://github.com/JuliaLang/IJulia.jl

1.4.
Nextjournal

you're looking at it right now!

  • reproducible science in article form
  • articles are interactive, fully versioned and automatically reproducible.
  • supports Python, R, Clojure and Julia
# query the documentation of a function
doc(+)
0.6s
# Install Packages
Pkg.pkg"add Plots"
Base.disable_logging(Base.CoreLogging.Error) # disable warnings
9.2s
# isolate code in Modules
module PlotExample
  using Plots; pyplot()
	function example()
    plot(rand(100) / 3,reg = true,fill = (0, :green))
    scatter!(rand(100), markersize = 6, c = :black)
	end
end

using .PlotExample

PlotExample.example()
77.4s

2.
Julia in one Cell

0.9s
# Function defined for all types T
function func(a::T) where T
  T
end

abstract type SuperType end

# inheritance
struct MyType <: SuperType
	field::Int
end
# mutable type without inheritance but with Type Parameter
mutable struct OtherType{TypeParam}
	field::TypeParam
end

# short function form
func(a::Float64, b::SuperType) = "a string"

func(b::OtherType{T}) where T = T

func(b::OtherType{Int}) = "oh hello there"

# tuple type
func(f, a::Tuple{Any, Any}) = (f(a[1]), f(a[2]))

@show func(1.0, MyType(1))
@show func(OtherType(1im))
@show func(OtherType(1))

# passing function
tuple = (OtherType(0), OtherType(0))
@show func(func, tuple)
# passing closure

closure = x-> x * 2
@show func(closure, (2, 3))
# the same as above, but with do syntax
result = func((2, 3)) do x
  x * 2
end
@show result
# construct array
array = [OtherType(1.0), OtherType(1.0)]
# dot syntax - apply a function to all elements of the argument arrays
@show func.(array)
nothing
0.9s

3.
Julia Core Strengths

1) Multiple Dispatch + Introspection

2) Efficient Duck Typing

3) Efficient Higher order Functions

4) Interprocedural Optimization

5) Dynamism + Meta programing

==> Performance + High Level programing

3.1.
Multiple Dispatch + Introspection

There are three core concepts to multiple dispatch:

  • dispatch on all argument types
  • specialize method body to argument types
  • types are compile time constants
dispatch(x::Int) = 1
dispatch(x::String) = 2
dispatch(x::Complex) = 3

test() = dispatch(1) + dispatch("2") + dispatch(1im)
0.6s
using InteractiveUtils
""" Helper to show the lowered code of a function"""
function show_asts(f, args...)
  argtypes = typeof(args)
  println("$f($args):")
  println(code_typed(f, argtypes, optimize = false)[1][1])
  println("optimized:")
  println(code_typed(f, argtypes)[1][1])
  println("------------------------")
  return
end
1.1s
show_asts(dispatch, 1)
2.7s
show_asts(test)
0.9s
function dispatch(x, ::Nothing)
  isa(x, Int) && return 1
  isa(x, String) && return 2
  isa(x, Complex) && return 3
end
test(x) = dispatch(1, x) + dispatch("2", x) + dispatch(1im, x)
0.7s
show_asts(dispatch, 1, nothing)
0.7s
show_asts(test, nothing)
0.7s
methods(dispatch)
1.7s
# generic base
my_func(a, b) = a + b
# in a package were you have more information:
my_func(a::Float32, b::Float32) = Intrinsics.fmul(a, b)
0.4s

3.1.1.
Enables you to:

  • compiled code with static type information inside function body
  • offer generic fallbacks for all kind of type combinations
  • specialize to type combinations
  • add new methods to existing libraries
  • efficient duck typing

3.2.
Duck Typing

In duck typing, an object's suitability is determined by the presence of certain methods and properties, rather than the type of the object itself.

function do_something(a)
  quack(a)
end

struct Duck end
quack(::Duck) = println("quoaak")
# by overloading a function, we get quacking & hence can `do_something`
quack(x::Int) = println("quack")

do_something(Duck())
do_something(1)
0.6s

3.2.1.
Practical Examples

3.2.1.1.
Algorithm Visualization
# untyped generic code: perfect to inject our own types :)
function magic(y::AbstractArray, func)
    l = length(y)
    k = ceil(Int, log2(l))
    for j=1:k, i=2^j:2^j:min(l, 2^k)
        y[i] = func(y[i-2^(j-1)], y[i])
    end
    for j=(k-1):-1:1, i=3*2^(j-1):2^j:min(l, 2^k)
        y[i] = func(y[i-2^(j-1)], y[i])
    end
    y
end
X = rand(10)
Y = magic(copy(X), +)
1.3s
import Base: getindex, setindex!, length, size
mutable struct AccessArray <: AbstractArray{Nothing, 1}
    length::Int
    read::Vector
    history::Vector
end
function AccessArray(length, read = [], history = [])
  AccessArray(length, read, history)
end
length(A::AccessArray) = A.length
size(A::AccessArray) = (A.length,)
function getindex(A::AccessArray, i)
    push!(A.read, i)
    nothing
end
function setindex!(A::AccessArray, x, i)
    push!(A.history, (A.read, [i]))
    A.read = []
end   
import Base.+
+(a::Nothing, b::Nothing) = a
0.4s
# we can now simply feed this new arrya type into magic - and record what it does
A = magic(AccessArray(8), +)
1.0s
140.7s
using GeometryTypes, Makie

function render(A::AccessArray)
    olast = depth = 0
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        olast = maximum(y[2])
    end
    maxdepth = depth
    olast = depth = 0
    C = []
    for y in A.history
        (any(y[1] .≤ olast)) && (depth += 1)
        push!(C, ((y...,), A.length, maxdepth, depth))
        olast = maximum(y[2])
    end
    msize = 0.1
    outsize = 0.15
    x1 = Point2f0.(first.(first.(first.(C))), last.(C) .+ outsize .+ 0.05)
    x2 = Point2f0.(last.(first.(first.(C))), last.(C) .+ outsize .+ 0.05)
    x3 = Point2f0.(first.(last.(first.(C))), last.(C) .+ 1)
    connections = Point2f0[]

    yoff = Point2f0(0, msize / 2)
    ooff = Point2f0(0, outsize / 2 + 0.05)
    for i = 1:length(x3)
        push!(connections, x3[i] .- ooff, x1[i] .+ yoff, x3[i] .- ooff, x2[i] .+ yoff)
    end
    node_theme = NT(
        markersize = msize, strokewidth = 3,
        strokecolor = :black, color = (:white, 0.0),
        axis = NT(
            ticks = NT(ranges = (1:8, 1:5)),
            names = NT(axisnames = ("Array Index", "Depth")),
            frame = NT(axis_position = :none)
        )
    )
  	s = Scene(resolution = (500, 500))
    s = scatter!(s, node_theme, x1)
    s = scatter!(s, node_theme, x2)
    scatter!(x3, color = :white, markersize = 0.2, strokewidth = 4, strokecolor = :black)
    scatter!(x3, color = :red, marker = '+', markersize = outsize)
    linesegments!(connections, color = :red)
    s
end
render(A)
140.7s
# Check our theory:
Y - [0; Y[1:end-1];] ≈ X
1.0s

[example adapted from Jiahao Che n]

3.2.1.2.
Automatic Differentiation
using ForwardDiff

f(x) = x^2 + 2.0x + 3.0
f_grad(x) = 2x + 2 # manual gradient function

val = 5
dual = ForwardDiff.Dual(val, 1)
dx = f(dual)
@show ForwardDiff.value(dx) == f(val)
@show ForwardDiff.partials(dx)[] == f_grad(val)
2.1s
using BenchmarkTools
vals = fill(val, 10^6)
duals = ForwardDiff.Dual.(vals, 1)
@btime f($val)
@btime f($dual)

@btime f.($vals)
@btime f.($duals)
nothing
30.2s

3.2.2.
Enables you to

  • compose functionality across packages
  • low effort to extend packages
  • use your own types that do extra work with generic APIs

3.3.
Efficient higher order function

using BenchmarkTools
# we annotate mutating functions usually with '!'
function map1!(f::F, dest, A) where F
    for i in 1:length(dest)
        @inbounds dest[i] = f(A[i])
    end
		return dest
end
# This is usually how other languages implement higher order functions
function map2!(@nospecialize(f), dest, A)
    for i in 1:length(dest)
        @inbounds dest[i] = f(A[i])
    end
		return dest
end
# some heavy functions
@inline test1(x) = x / 10
@noinline test2(x) = test1(x)# some heavy functions
A, dest = rand(10^6), zeros(10^6);
@btime map1!($test1, $dest, $A)
@btime map1!($test2, $dest, $A)
@btime map2!($test1, $dest, $A)
nothing
36.6s

3.3.1.
Enables you to:

  • implement mapreduce workflows with optimal performance
  • compose functionality freely

3.4.
Interprocedural optimization (IPO)

3.4.1.
What would Python do?

import Solver # Probably uses Numpy(C) + Numba
import Arithmetic # Maybe CPython + Python?

class MyType:
    ...
    
# CPython vs Numba, 
# Numpy & CPython vs MyClass
# some_function written in pure python? -> no specialization
Solver.solve(Arithmetic.some_function, convert(MyType, data))
Julia

Making a single Python algorithm fast in isolation is not that hard

making those algorithms generic and interface with other packages is hard!

  • use cython - loose generics and most higher level constructs & needs compilation
  • use numba - loose substantial subset of language, no IPO
  • use C & Python - loose everything 😛
using Solver, Arithmetic, GPUArrays, PrecisionFloats

data = convert(PrecisionFloats.Float, some_float_data)

# fast higher order functions, duck typing, 
# constant propagation + inlining across packages
Solver.solve(Arithmetic.some_function, GPUArray(data))
# Solver doesn't even need to know Arithmetic nor GPUArrays,
Julia

3.5.
Dynamics

str = "1 + 1"
expr = Meta.parse(str)
eval(expr)
0.3s
expr = quote 
  1 + 1
end
# or the short form
expr = :(1 + 1)
dump(expr) # dump = detailed output of fields
1.0s
eval(Expr(:call, +, 1, 1))
0.1s
macro replace_ab(expr, val)
  map!(expr.args, expr.args) do arg
    arg in (:a, :b) && return val
    arg
  end
  return expr
end
@replace_ab(a * b, "hi ")
0.4s

3.5.1.
Performance of dynamic constructs

3.5.1.1.
Macros - no cost

Macros are just a code transformation on the syntax level

3.5.1.2.
Eval - scoped performance penalty
using BenchmarkTools
function worst(a, expr::Expr)
  f = eval(expr)
  map!(a, a) do val
    Base.invokelatest(f, val)
  end
end

function best(f, a)
  map!(f, a, a)
end

x = rand(10^6)
expr = :((a)-> cos(a))
@btime worst($x, expr)

efunc = eval(:((a)-> cos(a)))
@btime best($efunc, $x)
nothing
25.3s
test(x) = eval(x)
show_asts(test, :(1 + 1))
2.5s

3.5.2.
Generated Functions

@generated function test(a, b)
  # Println is not really allowed here - Core.println works though
  Core.println(a, " ", b) 
  return :() # Need to return an Expr
end
0.3s
test(1, 2)
0.8s
@generated function test(a, b)
  Core.println(a, " ", b)
  return :(a + b) # Need to return an Expr
end
test(1, 2)
test(1, 2)
0.8s
expressions = []
@generated function map_unrolled(f, arg::NTuple{N, Any}) where N
  args = map(1:N) do i
    :(f(arg[$i]))
  end
  expr = Expr(:tuple, args...)
  push!(expressions, expr) # Core.println is really ugly for Expr
  expr
end
0.3s
map_unrolled(sin, (1.0, 2.0, 3.0))
0.7s
expressions[1]
0.9s

4.
Workshop