Julia : The Power of Multiple Dispatch

What is a numerical language?

  • Obviously, it's one that specializes in numerical work.

But specialisation means different things

  • Matlab
    • originally, everything was a matrix of complex doubles
    • everything else has been extended from that
  • R
    • data frames as basic type
    • support for "NA" values everywhere
  • Mathematica
    • symbolic evaluation semantics

Every system picks its own angle/niche to specialize.

Julia does NOT have a numerical speciality

  • instead, it is designed so you can define numeric types and their behaviors yourself

Design Overview

  • high-level & dynamic
  • expressive type-system
    • parametric, dependent, invariant
    • concrete types are final
      • but extensive hierarchy of abstract supertypes
      • lots of generic programming with abstract types
    • unobtrusive – don't need to mention types
  • ubiquitous multiple dispatch
    • everything is a generic function
      • even really basic, performance-cricical functions
    • diagonal dispatch
  • metaprogramming
    • homoiconic: code represented as Expr & Symbol objects, etc.
    • can be constructed, manipulated, eval'd
    • macros: @time sleep(1)
  • concurrency & parallelism
    • lightweight coroutine-based I/O
    • distributed global address space
      • easy run code on a cluster of instances
      • first-class remote references

Syntax Basics


In [1]:
function fib(n) 
    if n < 2 
        return n
    else
        return fib(n-1) + fib(n-2)
    end
end


Out[1]:
fib (generic function with 1 method)

In [2]:
fib(n) = n < 2 ? n : fib(n - 1) + fib(n - 2)


Out[2]:
fib (generic function with 1 method)

In [3]:
fib(30)


Out[3]:
832040

In [13]:
@time fib(26)


elapsed time: 0.000987882 seconds (96 bytes allocated)
Out[13]:
121393

In [5]:
A = randn(10,10)


Out[5]:
10x10 Array{Float64,2}:
 -0.574012   0.925585   1.32318     …   0.363387    0.0969455  -1.66204  
 -0.592042  -0.318013  -1.68856         0.236864    0.927813   -0.775747 
  1.18774   -0.171096   1.02094         1.36887    -0.857356   -0.258561 
  0.445789  -0.107505  -0.00775059     -0.173041    1.2877      0.299727 
  0.226309  -0.341116   2.15579         0.925124   -0.215497    0.294889 
 -1.28305    1.40863    0.147011    …   1.24728    -0.127498    0.0244787
  0.876873   1.79837   -0.734348        0.952655    0.363837   -1.44639  
 -1.26743    0.403198  -0.430558        0.283179    0.213105    0.482433 
 -0.474331  -1.71357   -0.209813        1.66785     1.03198     0.6563   
  0.897189   0.763462  -0.566816       -0.0770591  -0.732482    0.395473 

In [6]:
b = randn(10,10)


Out[6]:
10x10 Array{Float64,2}:
 -0.996164    0.773164   -0.509847   …   0.310905   -0.974049   -0.667838
  1.11639     0.465055   -1.30881       -1.14826     0.0343229  -0.830792
  1.9213      1.37146    -1.17559       -0.0292334   1.17679    -0.472685
 -1.76391     1.40594     0.185668      -2.18846    -0.39949    -0.939408
  0.59938     1.23924    -1.58336        0.465168    0.794852   -0.783582
  0.0105152  -0.920143    1.98287    …   1.15895     1.00033    -0.157435
  0.722799    2.07101    -0.641119      -2.00755    -1.25856     0.769873
 -0.807598   -1.60288     1.05169       -0.871025   -0.448713   -1.34616 
  2.14767    -0.541133    0.0462385      1.13291    -0.0402238   0.210212
  0.0435957  -0.0576593   0.0578214      0.319222   -1.8704      0.28868 

In [7]:
x = A \ b


Out[7]:
10x10 Array{Float64,2}:
  0.00879426   1.29178     -0.48704    …  -0.130574  -0.111159  -0.042147 
 -0.472644     0.171434     0.676167      -0.849743   0.232676   0.577437 
 -0.660724     0.553339     0.105551      -0.247019   1.24356   -0.463006 
 -0.783944    -0.396721     0.447256       1.06131   -1.81495   -0.461759 
 -0.720819     0.0153001    0.38196        1.36996   -0.877126  -0.699688 
 -0.330891     0.096513    -0.0387846  …  -0.626891   0.211223  -0.777928 
 -0.0535915   -0.172641    -0.412926       0.923612  -1.78634    0.0643431
  0.35763      0.00456275   0.23226        2.06488   -1.50413   -0.657474 
 -0.397819     0.80955      0.0171044     -1.3113     0.149858   0.570716 
 -0.424861    -0.716788     1.10731        0.364973  -0.243384   0.379121 

In [8]:
A*x - b


Out[8]:
10x10 Array{Float64,2}:
 -5.55112e-16   5.55112e-16   1.11022e-15  …   8.88178e-16  -5.55112e-16
  1.11022e-15  -7.21645e-16   0.0             -2.46331e-15   3.33067e-16
 -4.44089e-16   0.0          -2.22045e-16      8.88178e-16  -2.22045e-16
 -2.22045e-16   0.0          -1.66533e-16      2.22045e-16   5.55112e-16
  4.44089e-16   2.22045e-16   2.22045e-16     -1.22125e-15   1.11022e-16
  7.80626e-17   1.11022e-16   2.22045e-16  …   0.0           1.11022e-16
  1.11022e-16   0.0          -3.33067e-16      0.0          -5.55112e-16
  0.0          -2.22045e-16   2.22045e-16     -7.21645e-16  -6.66134e-16
  0.0          -1.11022e-16   3.95517e-16     -6.8695e-16   -2.77556e-16
  3.33067e-16   6.59195e-16  -4.57967e-16     -4.44089e-16  -3.33067e-16

In [9]:
eps()


Out[9]:
2.220446049250313e-16

In [10]:
round(A*x-b, 12)


Out[10]:
10x10 Array{Float64,2}:
 -0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0  -0.0
  0.0  -0.0   0.0  -0.0   0.0  -0.0  -0.0   0.0  -0.0   0.0
 -0.0   0.0  -0.0  -0.0   0.0   0.0   0.0  -0.0   0.0  -0.0
 -0.0   0.0  -0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0
  0.0   0.0   0.0  -0.0   0.0   0.0  -0.0  -0.0  -0.0   0.0
  0.0   0.0   0.0  -0.0  -0.0   0.0  -0.0  -0.0   0.0   0.0
  0.0   0.0  -0.0  -0.0  -0.0   0.0   0.0   0.0   0.0  -0.0
  0.0  -0.0   0.0  -0.0  -0.0   0.0  -0.0  -0.0  -0.0  -0.0
  0.0  -0.0   0.0  -0.0   0.0  -0.0   0.0  -0.0  -0.0  -0.0
  0.0   0.0  -0.0  -0.0  -0.0  -0.0   0.0   0.0  -0.0  -0.0

In [11]:
fft(A)


Out[11]:
10x10 Array{Complex{Float64},2}:
  11.7415+0.0im       -4.66568+9.90026im  …   -4.66568-9.90026im
 -5.06181-3.36135im   -2.61363-1.3522im      -0.443209+1.47246im
 -0.21693+4.5584im    -1.97886-3.90164im      0.358735+9.43976im
  9.23658+2.08115im    2.50709+3.97186im      -14.5272+15.2404im
 -1.38978+2.5823im     4.61229-3.92095im      -3.21567+6.01804im
 -5.12975+0.0im        13.2683-2.34219im  …    13.2683+2.34219im
 -1.38978-2.5823im    -3.21567-6.01804im       4.61229+3.92095im
  9.23658-2.08115im   -14.5272-15.2404im       2.50709-3.97186im
 -0.21693-4.5584im    0.358735-9.43976im      -1.97886+3.90164im
 -5.06181+3.36135im  -0.443209-1.47246im      -2.61363+1.3522im 

Multiple Dispatch

Basic Dispatch


In [12]:
f(a, b) = "fallback"
f(a::Number, b::Number) = "a and b are both numbers"
f(a::Number, b) = "a is a number"
f(a, b::Number) = "b is a number"
f(a::Integer, b::Integer) = "a and b are both integers"


Out[12]:
f (generic function with 5 methods)

In [13]:
f(1.5,2)


Out[13]:
"a and b are both numbers"

In [14]:
f(1, "bar")


Out[14]:
"a is a number"

In [15]:
f(1,2)


Out[15]:
"a and b are both integers"

In [16]:
f("foo", [1,2])


Out[16]:
"fallback"

In [17]:
methods(f)


Out[17]:
5 methods for generic function f:
  • f(a::Integer,b::Integer) at In[12]:5
  • f(a::Number,b::Number) at In[12]:2
  • f(a::Number,b) at In[12]:3
  • f(a,b::Number) at In[12]:4
  • f(a,b) at In[12]:1

Diagonal dispatch


In [18]:
f{T<:Number}(a::T, b::T) = "a and b are both $(T)s"


Out[18]:
f (generic function with 6 methods)

In [19]:
f(big(1.5),big(2.5))


Out[19]:
"a and b are both BigFloats"

In [20]:
f(big(1),big(2)) #<== integer rule is more specific


Out[20]:
"a and b are both integers"

In [21]:
f("foo","bar") #<== still doesn't apply to non-numbers


Out[21]:
"fallback"

Varargs Methods


In [22]:
f(args::Number...) = "$(length(args))-ary heterogeneous call"
f{T<:Number}(args::T...) = "$(length(args))-ary homogeneous call"


Out[22]:
f (generic function with 8 methods)

In [23]:
f(1)


Out[23]:
"1-ary homogeneous call"

In [24]:
f(1,2,3)


Out[24]:
"3-ary homogeneous call"

In [25]:
f(1,1.5,2)


Out[25]:
"3-ary heterogeneous call"

In [26]:
f() #==> heterogeneous because we can't bind T


Out[26]:
"0-ary heterogeneous call"

In [27]:
f(1,2) #<== previous 2-arg method is more specific


Out[27]:
"a and b are both integers"

In [28]:
f("foo") #<== doesn't apply to non-numbers


no method f(ASCIIString)
while loading In[28], in expression starting on line 1

The Expression Problem

Mads Torgersen's paper The Expression Problem Revisited:

To which degree can your application be structured in such a way that both the data model and the set of virtual operations over it can be extended without the need to modify existing code, without the need for code repetition and without runtime type errors.

Basically you want to be able to add:

  1. new types to which you can apply existing operations → easy in o.o., hard in functional
  2. new operations which you can apply to existing types → easy in functional, hard in o.o.

In [29]:
module ModInts
export ModInt
import Base: convert, promote_rule, show, showcompact

immutable ModInt{n} <: Integer
    k::Int
    ModInt(k) = new(k % n)
end

-{n2}(a::ModInt{n2}) = ModInt{n2}(-a.k)
+{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k+b.k)
-{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k-b.k)
*{n}(a::ModInt{n}, b::ModInt{n}) = ModInt{n}(a.k*b.k)

convert{n}(::Type{ModInt{n}}, i::Int) = ModInt{n}(i)
promote_rule{n}(::Type{ModInt{n}}, ::Type{Int}) = ModInt{n}

show{n}(io::IO, k::ModInt{n}) = print(io, "$(k.k) mod $n")
showcompact(io::IO, k::ModInt) = print(io, k.k)

end # module

In [30]:
using ModInts

In [31]:
ModInt


Out[31]:
ModInt{n} (constructor with 0 methods)

In [32]:
ModInt{11}


Out[32]:
ModInt{11} (constructor with 1 method)

In [33]:
a=ModInt{11}(12345532425)


Out[33]:
6 mod 11

In [34]:
b=ModInt{11}(23445134156)


Out[34]:
4 mod 11

In [35]:
a+b


Out[35]:
10 mod 11

In [36]:
a+b+b


Out[36]:
3 mod 11

In [37]:
a+1


Out[37]:
7 mod 11

In [38]:
2a


Out[38]:
1 mod 11

In [39]:
2a+1


Out[39]:
2 mod 11

In [40]:
A = map(ModInt{11}, rand(1:1000, 5,5))


Out[40]:
5x5 Array{ModInt{11},2}:
 10  7   3  8  10
  7  0   7  0   1
  9  5  10  8   8
  2  4  10  7   5
  8  1   1  9   0

In [41]:
A^2


Out[41]:
5x5 Array{ModInt{11},2}:
 8   6   1  8  6
 9   8   4  0  5
 9  10   8  5  6
 5   9   0  3  7
 4   9  10  3  2

In [52]:
A^10000000


Out[52]:
5x5 Array{ModInt{11},2}:
 4  1  9  8  3
 8  9  5  8  8
 9  2  7  6  1
 3  1  2  3  5
 1  1  3  2  9

In [43]:
B = map(ModInt{11}, rand(1:1000, 5,5))


Out[43]:
5x5 Array{ModInt{11},2}:
 2   8  8  5  7
 1  10  5  8  4
 0   5  3  6  6
 1   3  5  3  8
 1   3  4  7  8

In [44]:
A+B


Out[44]:
5x5 Array{ModInt{11},2}:
 1   4  0   2  6
 8  10  1   8  5
 9  10  2   3  3
 3   7  4  10  2
 9   4  5   5  8

In [45]:
foo{n}(a::ModInt{n}, b::ModInt{n}) = a^2 + 2b - 1
foo{n}(a::ModInt{n}, b::Int) = foo(a, ModInt{n}(b))
foo{n}(a::Int, b::ModInt{n}) = foo(ModInt{n}(a), b)
foo(a::Int, b::Int) = foo(ModInt{11}(a), ModInt{11}(b)).k


Out[45]:
foo (generic function with 4 methods)

In [46]:
foo(3,4)


Out[46]:
5

In [47]:
code_llvm(foo,(Int,Int))


define i64 @julia_foo15929(i64, i64) {
top:
  %2 = srem i64 %0, 11, !dbg !2947
  %3 = insertvalue %ModInt undef, i64 %2, 0, !dbg !2947, !julia_type !2948
  %4 = call %ModInt @julia_power_by_squaring15922(%ModInt %3, i64 2), !dbg !2947, !julia_type !2948
  %5 = extractvalue %ModInt %4, 0, !dbg !2947
  %6 = srem i64 %1, 11, !dbg !2947
  %7 = shl i64 %6, 1, !dbg !2947
  %8 = srem i64 %7, 11, !dbg !2947
  %9 = add i64 %8, %5, !dbg !2947
  %10 = srem i64 %9, 11, !dbg !2947
  %11 = add i64 %10, -1, !dbg !2947
  %12 = srem i64 %11, 11, !dbg !2947
  ret i64 %12, !dbg !2947
}

In [48]:
using Gadfly
@plot(cos(x)/x, 5, 25)


Out[48]:

Acknowledgements

Stefan Karpinski for inspiration for most of this presentation

Stefan, with Jeff Bezanson, Viral B Shah and Alan Edelman for creating Julia

Daniel C Jones for creating Gadfly


In [ ]: