Choose which topic you want to present.
If you want to discuss suitable topics, use the GitHub forum! (Or send me an email.)
Comments?
Here are two examples when programming languages try to be too smart.
diag(c(8,9,10))
gives the result $$ \begin{pmatrix} 8 & 0 & 0 \\ 0 & 9 & 0 \\ 0 & 0 & 10 \\ \end{pmatrix}. $$ and
diag(c(8,9))
gives $$ \begin{pmatrix} 8 & 0 \\ 0 & 9 \\ \end{pmatrix}. $$ But
diag(c(8))
gives $$ \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ \end{pmatrix}. $$
diag([1 2 3; 4 5 6])gives $$ \begin{pmatrix} 1 \\ 5 \\ \end{pmatrix}. $$ but
diag([1 2 3])gives $$ \begin{pmatrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 3 \\ \end{pmatrix}. $$
Both of these behaviors are quite ok when using R/Matlab interactively. But when the code is hidden deep inside layers of functions, it causes much more trouble than it's worth.
Julia solves this particular problem by having two different functions. Each defining a clear idea or protocol.
In [1]:
?diag
Out[1]:
In [2]:
?diagm
Out[2]:
1-dimensional arrays are column vectors.
In [3]:
v = [9, 8, 7]
Out[3]:
Matrices can be constructed in many ways.
In [4]:
A = [1 2 3; 4 5 6; 7 8 9]
Out[4]:
In [5]:
zeros(3,3)
Out[5]:
In [6]:
ones(2,4)
Out[6]:
In [7]:
eye(3)
Out[7]:
In [8]:
rand(1,2)
Out[8]:
In [9]:
rand(1:10,3,4)
Out[9]:
In [10]:
randn(4,2)
Out[10]:
In [11]:
A1 = [1 2; 3 4]
Out[11]:
In [12]:
A2 = [1.0 2.0; 3.0 4.0]
Out[12]:
In [13]:
zeros(A1)
Out[13]:
In [14]:
zeros(A2)
Out[14]:
Syntax:
[ F(x) for x=... ]
or
[ F(x,y) for y=..., x=... ]
and likewise for higher dimensions.
In [15]:
C = [ sin(i)+cos(j) for j=-5:5, i=-10:10 ]
Out[15]:
In [16]:
B1 = ones(2,3)
Out[16]:
In [17]:
B2 = 200*ones(2,3)
Out[17]:
Build matrix by blocks:
In [18]:
[B1 -B2; -B1 B2]
Out[18]:
where each $I_k$ may be:
In [19]:
A = [1 2 3; 4 5 6; 7 8 9]
Out[19]:
In [20]:
A[2,3]
Out[20]:
In [21]:
A[1:2,2]
Out[21]:
In [22]:
A[2,:]
Out[22]:
In [23]:
A[[3, 2], :]
Out[23]:
In [24]:
A[:, [false true true]]
Out[24]:
Elementwise:
min(A, B)
max(A, B)
In [25]:
B = rand(1:10,size(A))
Out[25]:
In [26]:
max(A,B)
Out[26]:
Find the minimum/maximum in an array:
minimum(A, dims)
maximum(A, dims)
In [27]:
maximum(B,1)
Out[27]:
In [28]:
maximum(B,2)
Out[28]:
In [29]:
A = rand(1:10,3,3)
Out[29]:
In [30]:
A * eye(A)
Out[30]:
Putting a . in front of the operator will apply it elementwise.
In [31]:
A .* eye(A)
Out[31]:
Let
In [32]:
A = rand(0:1,3,3)
Out[32]:
In [33]:
B = rand(0:1,3,3)
Out[33]:
And compare
In [34]:
A==B
Out[34]:
with
In [35]:
A.==B
Out[35]:
Broadcasting extends the idea of working elementwise by expanding singleton dimensions.
In [36]:
a = [1 2 3]
Out[36]:
In [37]:
b = [1,3,5]
Out[37]:
In [38]:
a+b
In [39]:
a.+b
Out[39]:
In [40]:
a.>=b
Out[40]:
In summary, the clear distinction between normal operators
and broadcasting operators
makes the code easier to read and less prone to unexpected side effects.
Why are types important in Julia?
Read :: as "is an instance of".
In [41]:
1+2::Integer
Out[41]:
In [42]:
1+2::AbstractFloat
Most commonly used for function arguments.
In [43]:
f(a, x::Number) = length(a)*x
Out[43]:
In [44]:
f("abc",0.5)
Out[44]:
In [45]:
f([1,2,3],2)
Out[45]:
In [46]:
f("abc","def")
Julia will compile different versions of the function f for every combination of arguments that is encountered!
This minimizes the (slow) type inference that the program needs to do at runtime.
In [47]:
typeof(f("abc",0.5))
Out[47]:
In [48]:
typeof(f("abc",2))
Out[48]:
Compare
for i=1:5
print(i)
end
and
for i=[1 2 3 4 5]
print(i)
end
The types are different
In [49]:
typeof(1:5)
Out[49]:
In [50]:
typeof([1 2 3 4 5])
Out[50]:
In [51]:
rng = 1:5
Out[51]:
In [52]:
length(rng)
Out[52]:
In [53]:
rng[2]
Out[53]:
The range object behaves (mostly) like an array, but it represents a the common abstraction of arrays that are lists of consecutive integers. This improves performance while giving readable code.
Multiple dispatch is used by Julia to achieve this behavior.
Code is said to be type-stable if variables do not change type.
A function is said to be type-stable if the return type can be deduced from the input types.
A motivating example:
In [54]:
function array_sum(a)
s = 0
for x in a
s += x
end
s
end
Out[54]:
In [55]:
myArray = randn(100000000);
In [56]:
array_sum(myArray)
Out[56]:
In [57]:
@time array_sum(myArray)
Out[57]:
In [58]:
function array_sum2(a)
s = zero(eltype(a))
for x in a
s += x
end
s
end
Out[58]:
In [59]:
array_sum2(myArray)
Out[59]:
In [60]:
@time array_sum2(myArray)
Out[60]:
Which of the following functions are type-stable?
In [61]:
e(x) = x<0 ? 0 : 1
f(x) = x<0 ? -x : x
g(x) = x==0 ? 1 : sin(x)/x
h(x,y) = x<y ? x : y
Out[61]:
In [62]:
y = e(-1)
y, typeof(y)
Out[62]:
In [63]:
y = e(1.0)
y, typeof(y)
Out[63]:
In [64]:
typeof(f(-3)), typeof(f(3))
Out[64]:
In [65]:
y = g(1)
y, typeof(y)
Out[65]:
In [66]:
y = g(1.0)
y, typeof(y)
Out[66]:
In [67]:
y = g(0)
y, typeof(y)
Out[67]:
In [68]:
y = h(1,2.0)
y, typeof(y)
Out[68]:
In [69]:
y = h(2,1.0)
y, typeof(y)
Out[69]:
Why should we write type-stable functions?
When a type-stable function is called, the compiler will know the return type and can avoid type checks. It also helps the outer function becoming type-stable.
This works.
In [70]:
2^3
Out[70]:
But this doesn't.
In [71]:
2^-3
However, we can make it work:
In [72]:
2.0^-3
Out[72]:
In [73]:
float(2)^-3
Out[73]:
In [74]:
2^-3.0
Out[74]:
To make 2^-3 return 0.125, either
or
Similar behavior for $\sqrt{}$.
In [75]:
sqrt(4)
Out[75]:
In [76]:
sqrt(-1)
In [77]:
sqrt(-1 + 0im)
Out[77]:
In [78]:
sqrt(complex(-1))
Out[78]:
This can also be helpful to track down bugs.
It's better to get an error message than getting unexpected complex numbers later in the program!
More on type-stability: https://www.youtube.com/watch?list=PLP8iPy9hna6Sdx4soiGrSefrmOPdUWixM&v=L0rx_Id8EKQ
Julia has a type hierarchy. Only leaf types are concrete. (As opposed to object-oritented languages like C++ or Java.)
Read <: as "is a subtype of"
abstract Number
abstract Real <: Number
abstract AbstractFloat <: Real
abstract Integer <: Real
abstract Signed <: Integer
abstract Unsigned <: Integer
Useful when defining functions.
In [79]:
square(x::Number) = x*x
Out[79]:
In [80]:
square(2)
Out[80]:
In [81]:
square(1.5)
Out[81]:
In [82]:
square("abc")
We can define our own types. They will be just as fast as the built-in types in Julia (if implemented correctly).
In [83]:
type Person
name::AbstractString
age::Int
end
Specifying the types of the fields is not necessary, but can improve performance for the same reasons as we have seen above.
In [84]:
p = Person("Anna Svensson", 42)
Out[84]:
In [85]:
p.name
Out[85]:
In [86]:
p.age
Out[86]:
Many simple types in Julia are immutable. It means that they are not allowed to be modified after construction. It greatly helps the compiler create fast code in many cases.
Julia uses pass-by-sharing for function arguments.
In [87]:
a = [1,2,3,4]
Out[87]:
In [88]:
function dostuff!(v)
v[1] = 100
end
Out[88]:
In [89]:
dostuff!(a)
Out[89]:
In [90]:
a
Out[90]:
However
In [91]:
function dostuff2(v)
v = -v
end
Out[91]:
In [92]:
dostuff2(a)
Out[92]:
In [93]:
a
Out[93]:
Here -v creates a new array and the variable v is rebound to the new array. Hence we do not change the original array a.
In [94]:
import Base: +, -, *
Definition of type.
In [95]:
immutable ModInt{n} <: Integer
k::Int
ModInt(k) = new(mod(k,n))
end
Arithmetic operations.
In [96]:
-{n}(a::ModInt{n}) = ModInt{n}(-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)
Out[96]:
Conversion and promotions (more on these later).
In [97]:
Base.convert{n}(::Type{ModInt{n}}, i::Int) = ModInt{n}(i)
Base.promote_rule{n}(::Type{ModInt{n}}, ::Type{Int}) = ModInt{n}
Out[97]:
How to print it.
In [98]:
Base.show{n}(io::IO, k::ModInt{n}) = print(io, "$(k.k) mod $n")
Base.showcompact(io::IO, k::ModInt) = print(io, k.k)
Out[98]:
Inversion.
In [99]:
Base.inv{n}(a::ModInt{n}) = ModInt{n}(invmod(a.k, n))
Out[99]:
Usage:
In [100]:
a = ModInt{11}(120)
Out[100]:
In [101]:
b = ModInt{11}(987)
Out[101]:
In [102]:
a+b
Out[102]:
In [103]:
a+2
Out[103]:
Since a lot of functionality is built on top of the definitions above, we can do many more things!
Create matrices:
In [104]:
A = map(ModInt{13}, rand(1:1000,5,5))
Out[104]:
In [105]:
B = map(ModInt{13}, rand(1:1000,5,5))
Out[105]:
Operate on matrices:
In [106]:
A+B
Out[106]:
In [107]:
A*B
Out[107]:
In [108]:
A.*B
Out[108]:
In [109]:
2A^3 - 3B + 2I
Out[109]:
We can also use common functions:
In [110]:
sum(A)
Out[110]:
In [111]:
cumsum(A[1,:])
Out[111]:
Remember that we think of functions as ideas or protocols.
In 14 lines of code, we got a new type that is already very useful!
There are rules for converting and promoting types in many languages. The difference is that in Julia this is exposed to the programmer.
In [112]:
1 + 3.2
Out[112]:
Let's break down what Julia is doing here.
+ is defined between for Int+Int and for Float64+Float64. And then there is this general definition for Numbers (that can be of different types):
+(x::Number, y::Number) = +(promote(x,y)...)
which is invoked above.
The function promote converts the arguments to a common type:
In [113]:
promote(1,3.2)
Out[113]:
Which in turn relies on
In [114]:
promote_type(Int,Float64)
Out[114]:
to find the common type and
In [115]:
convert(Float64,1)
Out[115]:
to do the conversion.
For our own types, we need to define
convert()
and
promote_rule()
as in the ModInt example above.
Let's make a positive definite symmetric matrix.
In [116]:
N = 5
Out[116]:
In [117]:
A = randn(N,N)
A = A'A
Out[117]:
And compute the Cholesky factorization. (Every positive definite symmetric matrix $A$ has a decomposition $A=U^TU$.)
In [118]:
C = cholfact(A)
Out[118]:
Note that variable $C$ represents the decomposition $U^TU$.
In [119]:
Y = randn(N)
Out[119]:
Now let's solve for $X$ in $AX=Y$.
In [120]:
A\Y
Out[120]:
Now let's solve for $X$ in $U^TUX=Y$.
In [121]:
C\Y
Out[121]:
This way, we can work with a more efficient representation while keeping the code easy to read.
Types are necessary to make the code run fast.
However, the Julia JIT compiler does much of the work for us. Many times, we don't need to specify types in our code and Julia will make it fast anyway.
High-level code doesn't need bother too much with types.
Packages and Julia library functions tend to be well-written and fast. They will be fast even if there are some type problems in your high-level code.
If you write low-level, performance critical code.
Type stability is key.
Types are necessary for Multiple Dispatch.
Which is absoletely central for code abstraction in Julia.
We should be aware of how types work in Julia.
To know when we need to care.
To understand why certain things are designed the way they are.