introduction-to-julia

Let's start with a quick overview of the basic syntax, emphasising differences with Python.

Why Jula (slides show)

  • David P. Sanders
  • Department of Physics, Faculty of Sciences
  • National University of Mexico (UNAM)
  • Twitter: @DavidPSanders
  • GitHub: dpsanders

Variables and arithmetic

Numeric values in Julia work similarly to Python:

Variables in Julia are created as in Python, with a simple assignment operation; variable names can be arbitrary unicode characters. Many may be entered in LaTeX notation, using tab substitution: type \alpha<TAB>. There is also tab completion on partial names: \alp<TAB>

In [1]:
x = 3
Out[1]:
3
In [2]:
y = 5
Out[2]:
5
In [3]:
α = 3;  = 10
Out[3]:
10

Functions use parentheses (round brackets, ()) around the arguments being passed. println prints its arguments, followed by a new line. [print omits the new line.]

In [4]:
println("α = ", α)
α = 3

Simple functions may be defined with a nice mathematical syntax; * is not needed in simple expressions:

In [5]:
f(x) = 2x^2 + 3x + 1
g(x) = f(x) - (2x+1)*(x+1)
Out[5]:
g (generic function with 1 method)
In [6]:
f(3)
Out[6]:
28
In [7]:
g(3.5)
Out[7]:
0.0

Variable substitution

The values of variables may be substituted into strings in a simple way using the $ operator:

In [8]:
# Variable substitution with $:
name = "David"
greeting = "Hello, $name"
println(greeting)
Hello, David

More complicated expressions are wrapped in parentheses:

In [9]:
μ = 3
println("The sine of $μ is $(sin(μ))")
The sine of 3 is 0.1411200080598672

Numerical types

There are numerical types with different precisions: typing Float<TAB> or Int<TAB> will provide a list. Currently, in arithmetic calculations types are promoted to the machine type. (This looks likely to change soon.)

Machine integers!

In [10]:
a = int(1e16)
a * 10
Out[10]:
100000000000000000
In [11]:
a = int8(1)
b = int8(2)
a + b
Out[11]:
3
In [12]:
typeof(ans) # ans is the last result
Out[12]:
Int8

These promotion rules are defined in int.jl.

Arbitrary-precision arithmetic

Arbitrary-precision integers and floating points are available through the types BigInt and BigFloat. The function big converts a number into the corresponding Big type:

In [13]:
big(10)
Out[13]:
10
In [14]:
typeof(ans)
Out[14]:
BigInt (constructor with 10 methods)

Note that, unlike in Python, integers are not automatically promoted to arbitrary-precision integers.


Exercise: Calculate powers of 10 using standard integers and BigInts


In [15]:
10^18
Out[15]:
1000000000000000000
In [16]:
10**5
syntax: use "^" instead of "**"
while loading In[16], in expression starting on line 1
In [17]:
10^19
Out[17]:
-8446744073709551616
In [18]:
ten = big(10)
Out[18]:
10
In [19]:
ten^109
Out[19]:
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
In [20]:
i = int8(10)
Out[20]:
10
In [21]:
i * big(11)
Out[21]:
110
In [22]:
typeof(ans)
Out[22]:
BigInt (constructor with 10 methods)
In [23]:
typemax(Int64)
Out[23]:
9223372036854775807
In [24]:
2^63-1
Out[24]:
9223372036854775807
In [25]:
typemin(Int64)
Out[25]:
-9223372036854775808
In [26]:
-2^63
Out[26]:
-9223372036854775808

Complex numbers

Complex numbers are written using im for the imaginary part:

In [27]:
a = 7
c = (1+3.5im) * a
Out[27]:
7.0 + 24.5im
In [28]:
c.im
Out[28]:
24.5
In [29]:
c.re, c.im
Out[29]:
(7.0,24.5)
In [30]:
c * conj(c)  # conj is a function that returns the conjugate of a complex number
Out[30]:
649.25 + 0.0im

(Tuples behave similarly to Python.)

Rational numbers

Rational numbers are also built into Julia; they are created using the // operator:

In [31]:
3//4
Out[31]:
3//4
In [32]:
typeof(ans)
Out[32]:
Rational{Int64} (constructor with 1 method)
In [33]:
(big(3)//4)^50
Out[33]:
717897987691852588770249//1267650600228229401496703205376
In [34]:
typeof(ans)
Out[34]:
Rational{BigInt} (constructor with 1 method)
In [35]:
3//4 + 5//6
Out[35]:
19//12

Operators are a convenient way of writing functions:

In [36]:
+(3, 4)
Out[36]:
7
In [37]:
//(3, 4)
Out[37]:
3//4
In [38]:
//
Out[38]:
// (generic function with 8 methods)

We see that // is a function, implemented as a series of methods. We can see what these methods are:

In [39]:
methods(//)
Out[39]:
8 methods for generic function //:
In [40]:
3 + -34
Out[40]:
-31

The expression n::Integer is a type annotation that specifies that the method applies when its first argument is of type Integer.

Clicking on the file name takes us directly to the Julia standard library source code on GitHub where these functions are defined!

Vectors: an "equivalent" of Python lists and numpy arrays

To store several values "in one variable", we can try to imitate using a "list" as we would in Python:

In [41]:
l = [3, 4, 5]
Out[41]:
3-element Array{Int64,1}:
 3
 4
 5
In [42]:
typeof(l)
Out[42]:
Array{Int64,1}

In Julia these objects are called Arrays. The curly braces indicate type parameters of the Array type. The first is the type of element contained in the Array (all must be of the same type) and the second the number of dimensions.


Exercise: Try to create an array in this way with elements of different types. What happens?

Exercise: What does the following syntax do?: l = {3, 4, 7.5}


In [43]:
l = [3, 4, 7.5]
Out[43]:
3-element Array{Float64,1}:
 3.0
 4.0
 7.5
In [44]:
l = [3., 4, 5]
Out[44]:
3-element Array{Float64,1}:
 3.0
 4.0
 5.0
In [45]:
l = [3., "a"]
Out[45]:
2-element Array{Any,1}:
 3.0 
  "a"
In [46]:
l = [3., 'a']
Out[46]:
2-element Array{Float64,1}:
  3.0
 97.0
In [47]:
l = {3., 4, "hello", [3, 4]}
Out[47]:
4-element Array{Any,1}:
 3.0     
 4       
  "hello"
  [3,4]  

Indexing

The indices of Julia arrays are numbered starting at 1, unlike Python (where they are numbered starting at 0).

In [48]:
l[1]  
Out[48]:
3.0

The syntax for ranges is similar to that for Python:

In [49]:
l[1:2]
Out[49]:
2-element Array{Any,1}:
 3.0
 4  

However, the limits must be explicitly specified:

In [50]:
l[2:end]   # Use `end` explicitly
Out[50]:
3-element Array{Any,1}:
 4       
  "hello"
  [3,4]  
In [51]:
l[1:end-1]
Out[51]:
3-element Array{Any,1}:
 3.0     
 4       
  "hello"
In [52]:
l[-1]
BoundsError()
while loading In[52], in expression starting on line 1

 in getindex_3B_2413 at /home/tj2/Documents/julia/usr/bin/../lib/julia/sys.so

Julia Arrays, like Python lists, but unlike numpy arrays, are dynamic. However, the syntax is rather different from Python -- to add an element at the end of the list, we write

In [53]:
l = [3,4,5]
l + l
Out[53]:
3-element Array{Int64,1}:
  6
  8
 10
In [54]:
names(l)
Out[54]:
0-element Array{Any,1}
In [55]:
l = [3, 4, 5]

push!(l, 7)
Out[55]:
4-element Array{Int64,1}:
 3
 4
 5
 7
In [56]:
methods(sizehint)
Out[56]:
11 methods for generic function sizehint:
In [57]:
push!
Out[57]:
push! (generic function with 22 methods)
In [58]:
methods(push!)
Out[58]:
22 methods for generic function push!:
In [60]:
l = [3, 4, 5]
push!(l, 12.001)
InexactError()
while loading In[60], in expression starting on line 2

 in push! at array.jl:451
In [61]:
l = [3.0, 4, 5]
push!(l, 12.001)
Out[61]:
4-element Array{Float64,1}:
  3.0  
  4.0  
  5.0  
 12.001
In [62]:
12.0 == 12
Out[62]:
true
In [63]:
push!(l, 12.1)
Out[63]:
5-element Array{Float64,1}:
  3.0  
  4.0  
  5.0  
 12.001
 12.1  
In [64]:
l = [3, 4, 5, 12]
append!(l, [10., 11., 12])
Out[64]:
7-element Array{Int64,1}:
  3
  4
  5
 12
 10
 11
 12

push! replaces append in Python. There are no methods of objects as in Python; rather, we use functions and send the object as an argument of the function.

The exclamation mark, or bang, (!) indicates that the function modifies its argument; this is a standard convention in Julia.

Arrays which have been defined with a certain type cannot acquire elements of a different type:

In [65]:
l = [3, 4, 5]
push!(l, "hello")
`convert` has no method matching convert(::Type{Int64}, ::ASCIIString)
while loading In[65], in expression starting on line 2

 in push! at array.jl:451

Arrays work as mathematical vectors, with the sum of two vectors and scalar multiplication being defined:

In [66]:
a = [1.1, 2.2, 3.3]
b = [4.4, 5.5, 6.6]
Out[66]:
3-element Array{Float64,1}:
 4.4
 5.5
 6.6
In [67]:
a + b
Out[67]:
3-element Array{Float64,1}:
 5.5
 7.7
 9.9
In [68]:
3.5 * a
Out[68]:
3-element Array{Float64,1}:
  3.85
  7.7 
 11.55

However, operators are, in general, not treated in an elementwise fashion (as they would be e.g. in numpy):

In [69]:
a * b
`*` has no method matching *(::Array{Float64,1}, ::Array{Float64,1})
while loading In[69], in expression starting on line 1

Rather, elementwise operations use a Matlab-like syntax, with an extra . before the symbol for the operator:

In [70]:
a .* b
Out[70]:
3-element Array{Float64,1}:
  4.84
 12.1 
 21.78

There are many useful operations on vectors predefined, without needing to explicitly import them.

In [71]:
dot(a,b)
Out[71]:
38.72
In [72]:
?dot
dot (generic function with 7 methods)
In [73]:
cross(a, b)
Out[73]:
3-element Array{Float64,1}:
 -3.63
  7.26
 -3.63
In [74]:
?cross
cross (generic function with 1 method)
In [75]:
norm(a)
Out[75]:
4.115823125451335
In [76]:
?norm
Base.norm(A[, p])

   Compute the "p"-norm of a vector or the operator norm of a matrix
   "A", defaulting to the "p=2"-norm.

   For vectors, "p" can assume any numeric value (even though not
   all values produce a mathematically valid vector norm). In
   particular, "norm(A, Inf)" returns the largest value in
   "abs(A)", whereas "norm(A, -Inf)" returns the smallest.

   For matrices, valid values of "p" are "1", "2", or "Inf".
   (Note that for sparse matrices, "p=2" is currently not
   implemented.) Use "vecnorm()" to compute the Frobenius norm.

Use help or ? (before the command) to obtain help:

In [77]:
help(dot)
dot (generic function with 7 methods)
In [78]:
?dot
dot (generic function with 7 methods)
In [79]:
transpose(a)
Out[79]:
1x3 Array{Float64,2}:
 1.1  2.2  3.3
In [80]:
a'
Out[80]:
1x3 Array{Float64,2}:
 1.1  2.2  3.3
In [81]:
transpose(a) == a'
Out[81]:
true
In [82]:
M = [[2,1], [1,1]] # 1-d
Out[82]:
4-element Array{Int64,1}:
 2
 1
 1
 1
In [83]:
push!(l, 12.0)
Out[83]:
4-element Array{Int64,1}:
  3
  4
  5
 12
In [84]:
M = [2. 1; 1 1] # 2-d
Out[84]:
2x2 Array{Float64,2}:
 2.0  1.0
 1.0  1.0
In [85]:
M = reshape([1:8], (2,2,2))
Out[85]:
2x2x2 Array{Int64,3}:
[:, :, 1] =
 1  3
 2  4

[:, :, 2] =
 5  7
 6  8
In [86]:
a  b
Out[86]:
38.72
In [87]:
dot(a,b) == a  b
Out[87]:
true
In [88]:
a × b
Out[88]:
3-element Array{Float64,1}:
 -3.63
  7.26
 -3.63
In [89]:
cross(a,b) == a × b
Out[89]:
true

I used \cdot

In [90]:

Out[90]:
dot (generic function with 7 methods)

[Note that in the Julia command-line REPL, typing ? puts it immediately into a special help mode. Similarly, ; puts it into shell mode, in which commands are sent straight to the shell.]

Control flow: drop the colon (:) and add end

White space in Julia is not significant. Commands on one line can be separated by ;. Blocks must finish with end

In [91]:
i = 0
while i < 5 
    print("$i\t")
    i += 1
end
0	1	2	3	4	
In [92]:
total = 0
for i = 1:26
    total += 2^i
end
println("Sum is $total")
Sum is 134217726

Here, 1:26 is a range object which may be iterated over.

In [93]:
typeof(1:10)
Out[93]:
UnitRange{Int64} (constructor with 1 method)

We can construct an array from this by enclosing it in square brackets:

In [94]:
collect(1:10)
Out[94]:
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
In [95]:
?collect
Base.collect(collection)

   Return an array of all items in a collection. For associative
   collections, returns (key, value) tuples.

Base.collect(element_type, collection)

   Return an array of type "Array{element_type,1}" of all items in a
   collection.
In [96]:
[1:2:10, 17]
Out[96]:
6-element Array{Int64,1}:
  1
  3
  5
  7
  9
 17
In [97]:
help("dot")
Base.LinAlg.BLAS.dot(n, X, incx, Y, incy)

   Dot product of two vectors consisting of "n" elements of array
   "X" with stride "incx" and "n" elements of array "Y" with
   stride "incy".

Exercise: Implement the Babylonian method for calculating the square root of a positive number $y$, via the iteration $$x_{n+1} = \textstyle \frac{1}{2} (x_n + \frac{y}{x_n})$$

In [98]:
x = 7
# Babylonian method
r1 = big(x)
for i=1:10
    r2 = (r1+x//r1)//2
    r1 = r2
    println("$i: $r2 is $(float(r2))")
end
sqrt(x)
1: 4//1 is 4e+00
2: 23//8 is 2.875e+00
3: 977//368 is 2.654891304347826086956521739130434782608695652173913043478260869565217391304333e+00
4: 1902497//719072 is 2.645767044190289706733122691469004494682034622402207289395220506430510435672636e+00
5: 7238946623297//2736064645568 is 2.645751311111369320474788060415261270876928274529823742547961659814056687694087e+00
6: 104804696428033056657448577//39612451854313553433195392 is 2.645751311064590590502029293938515493859933284910652657274139597282698886360655e+00
7: 21968048786744329890159858491687242656637852949560577//8303141982722814210287293056980408998691454066714368 is 2.645751311064590590501615753639260425710259215401356797324429072042618334785781e+00
8: 965190334993558048952747224707737662553546023969581717336893753948305178090306054973804711697153676048897//364807656319439656893574197813762741889129729119092813240115363343939412817707718620985339231655544540672 is 2.645751311064590590501615753639260425710259183082450180368334459201068823230298e+00
9: 1863184765529953614442676969410516768705348877426018592874619682164981370927579559773987864606096009580887559742891544147556039040815814753743851100010309178655530447444083262328821461724034060383304637370397697//704217648022349512771442189018194725339511495157190929516136319260161990261007461575776992568187313708188185637682782310132281240170939460318854442124039911242297405920199037531309472261446561026954943354477568 is 2.645751311064590590501615753639260425710259183082450180368334459201068823230264e+00
10: 6942914941005816452906820447809898528401948252945007180207641686978710969653763380353810662539739407094605662843293836754945636431635508103618963852112274901648239421682834044604657380796983687411621541616911160282084512012755071820846262756649514791264502834906610162581301591431279046989968494539999498942229056113659738389773073293970820133886115413882850641756087247746883347629340780021440048884451568277591449266177//2624175186825153359249378049906502688043810272014845844292434248994264560667710669982033339453245573538963306071655131184848403463172113026723021373675261241203122087530398269881154662495920113030441863771241633914463530742637259717641425928569498576675875622441689707127692085122998906191651056070668909263006767739390431278990942834266700554990212884068483635241604906439758589382529929519521217527464936293523650721792 is 2.645751311064590590501615753639260425710259183082450180368334459201068823230264e+00
Out[98]:
2.6457513110645907
In [99]:
 = norm
Out[99]:
norm (generic function with 15 methods)
In [100]:
(a)
Out[100]:
4.115823125451335
In [101]:
snowman(x) = x^3
Out[101]:
snowman (generic function with 1 method)
In [102]:
 = snowman
(3)
Out[102]:
27

Short-circuit evaluation

In [103]:
a = 3
a < 5 && println("Small")   # evaluate the second statement only if the first is true;  semantics of if-then

a > 10 || println("Small")  # semantics of if not-then

# equal to

if a < 5
    println("small, if-then")
end

if a > 10
    println("small, if not-then")
end
Small
Small
small, if-then
In [104]:
a == 3 ? println("Hello") : println("Not true")

# equal to

if a < 5
    println("small, if-then")
else
    println("small, if not-then")
end
Hello
small, if-then

Array comprehensions

There is an equivalent of list comprehensions in Python, as follows. Note that the array construction syntax is quite flexible.

In [105]:
squares = [i^2 for i in [1:2:10, 7]]
Out[105]:
6-element Array{Any,1}:
  1
  9
 25
 49
 81
 49
In [106]:
sums = [i+j for i=1:5, j=1:5] # boardcasting
Out[106]:
5x5 Array{Int64,2}:
 2  3  4  5   6
 3  4  5  6   7
 4  5  6  7   8
 5  6  7  8   9
 6  7  8  9  10
In [107]:
sums = [i+j+k for i=1:5, j=1:5, k=1:5]
Out[107]:
5x5x5 Array{Int64,3}:
[:, :, 1] =
 3  4  5   6   7
 4  5  6   7   8
 5  6  7   8   9
 6  7  8   9  10
 7  8  9  10  11

[:, :, 2] =
 4  5   6   7   8
 5  6   7   8   9
 6  7   8   9  10
 7  8   9  10  11
 8  9  10  11  12

[:, :, 3] =
 5   6   7   8   9
 6   7   8   9  10
 7   8   9  10  11
 8   9  10  11  12
 9  10  11  12  13

[:, :, 4] =
  6   7   8   9  10
  7   8   9  10  11
  8   9  10  11  12
  9  10  11  12  13
 10  11  12  13  14

[:, :, 5] =
  7   8   9  10  11
  8   9  10  11  12
  9  10  11  12  13
 10  11  12  13  14
 11  12  13  14  15

Matrices

Square brackets with commas gives a one-dimensional vector. This is printed in a way that treats it as if it were a column vector (although there is in fact no difference between a one-dimensional row vector and column vector).

In [108]:
v = [3, 4, 5]
Out[108]:
3-element Array{Int64,1}:
 3
 4
 5

To create explicit matrices, Matlab-style notation is used. If we omit the commas, something different happens: we now obtain a two-dimensional Array, i.e. a matrix, of size $1 \times n$. [Recall that in the standard notation for matrices, an $m \times n$ matrix has $m$ rows and $n$ columns.]

In [109]:
row_vec = [3 4 5]
Out[109]:
1x3 Array{Int64,2}:
 3  4  5
In [110]:
v == row_vec
Out[110]:
false

We can also use the transpose operator, '. [This is actually the conjugate-transpose operator, which also takes the complex conjugate of complex numbers. Transpose without conjugate is denoted .']

In [111]:
row_vec = [2im, 2]'
Out[111]:
1x2 Array{Complex{Int64},2}:
 0-2im  2+0im
In [112]:
row_vec = [2im, 2].'
Out[112]:
1x2 Array{Complex{Int64},2}:
 0+2im  2+0im

A complete matrix may be constructed using a semicolon (;) to separate rows:

In [113]:
M = [1 2im; 3 4]
Out[113]:
2x2 Array{Complex{Int64},2}:
 1+0im  0+2im
 3+0im  4+0im

As in numpy, it may also be created using a reshape:

In [114]:
M = reshape([1 2 3 4], (2,2))

# equal to

N = reshape([1, 2, 3, 4], (2,2))
Out[114]:
2x2 Array{Int64,2}:
 1  3
 2  4
In [115]:
M == N
Out[115]:
true

Here, as in Python, (2,2) denotes an (immutable) tuple:

In [116]:
t = (2, 2)
typeof(t)
Out[116]:
(Int64,Int64)

There is an important difference in the way that Python and Julia treat slices of matrices. While in Python a one-dimensional slice in either direction returns a 1-dimensional vector, in Julia there is a difference. A vertical one-dimensional slice gives a 1-dimensional vector (a "column vector"):

In [117]:
M[:,1]
Out[117]:
2-element Array{Int64,1}:
 1
 2

However, a horizontal one-dimensional slice produces a $1 \times n$ matrix:

In [118]:
M[1,:]
Out[118]:
1x2 Array{Int64,2}:
 1  3

This is the same result that is produced using the following Matlab-like syntax:

In [119]:
[1 2]
Out[119]:
1x2 Array{Int64,2}:
 1  2

Random numbers

The Mersenne Twister (pseudo-)random number generator is built-in to Julia:

In [120]:
rand()
Out[120]:
0.40186768641976944
In [121]:
rand(5)
Out[121]:
5-element Array{Float64,1}:
 0.344519
 0.617192
 0.829609
 0.981215
 0.477052
In [122]:
x = rand(5, 5)
Out[122]:
5x5 Array{Float64,2}:
 0.634627  0.121441   0.793247  0.78431   0.965882
 0.989307  0.616908   0.353485  0.710715  0.635409
 0.352762  0.118025   0.185011  0.280691  0.28924 
 0.748569  0.195795   0.279713  0.831955  0.888495
 0.765     0.0837347  0.537341  0.429161  0.147381
In [123]:
y = rand(Int32,3,3)
Out[123]:
3x3 Array{Int32,2}:
  -113324951   196613133  -1363228423
 -2088630349  1597727404   -358295295
 -1080816523  1898363541    571418313
In [124]:
?rand
Base.rand() -> Float64

   Generate a "Float64" random number uniformly in [0,1)

Base.rand(rng::AbstractRNG[, dims...])

   Generate a random "Float64" number or array of the size specified
   by dims, using the specified RNG object. Currently,
   "MersenneTwister" is the only available Random Number Generator
   (RNG), which may be seeded using srand.

Base.rand(dims or [dims...])

   Generate a random "Float64" array of the size specified by dims

Base.rand(t::Type[, dims...])

   Generate a random number or array of random numbes of the given
   type.

Base.rand(r[, dims...])

   Pick a random element or array of random elements from range "r"
   (for example, "1:n" or "0:2:10").
In [125]:
methods(rand)
Out[125]:
33 methods for generic function rand:

Matrix multiplication

In [126]:
v = [1, 2]
Out[126]:
2-element Array{Int64,1}:
 1
 2
In [127]:
v*v
`*` has no method matching *(::Array{Int64,1}, ::Array{Int64,1})
while loading In[127], in expression starting on line 1
In [128]:
v^2
`*` has no method matching *(::Array{Int64,1}, ::Array{Int64,1})
while loading In[128], in expression starting on line 1

 in power_by_squaring at intfuncs.jl:56
 in ^ at intfuncs.jl:86
In [129]:
dot(v, v)
Out[129]:
5
In [130]:
M = [2 1; 1 1]
Out[130]:
2x2 Array{Int64,2}:
 2  1
 1  1
In [131]:
dot(M, v)
`dot` has no method matching dot(::Array{Int64,2}, ::Array{Int64,1})
while loading In[131], in expression starting on line 1

Matrix multiplication uses the * operator:

In [132]:
M * v
Out[132]:
2-element Array{Int64,1}:
 4
 3
In [133]:
@which M*v
Out[133]:
*{T,S}(A::AbstractArray{T,2},x::AbstractArray{S,1}) at linalg/matmul.jl:71

Exercise: Use the power method to calculate the largest eigenvalue $\lambda_1$ of the matrix $M = \begin{pmatrix} 2 & 1 \\ 1 & 1 \end{pmatrix}$. In this method, we start from an arbitrary non-zero vector $\mathbf{w}$, and repeatedly apply $M$ to it, thus calculating powers of the matrix $M$ applied to $\mathbf{w}$. The resulting vector converges to the eigenvector $\mathbf{v}_1$ corresponding to $\lambda_1$.

In [134]:
w = [1., 1]
M = reshape([2., 1, 1, 1], (2,2))

M, w
Out[134]:
(
2x2 Array{Float64,2}:
 2.0  1.0
 1.0  1.0,

[1.0,1.0])
In [135]:
eig(M)
Out[135]:
([0.381966,2.61803],
2x2 Array{Float64,2}:
  0.525731  -0.850651
 -0.850651  -0.525731)
In [136]:
?eig
Base.eig(A,[irange,][vl,][vu,][permute=true,][scale=true]) -> D, V

   Computes eigenvalues and eigenvectors of "A". See "eigfact()"
   for details on the "balance" keyword argument.

      julia> eig([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
      ([1.0,3.0,18.0],
      3x3 Array{Float64,2}:
       1.0  0.0  0.0
       0.0  1.0  0.0
       0.0  0.0  1.0)

   "eig" is a wrapper around "eigfact()", extracting all parts of
   the factorization to a tuple; where possible, using "eigfact()"
   is recommended.

Base.eig(A, B) -> D, V

   Computes generalized eigenvalues and vectors of "A" with respect
   to "B".

   "eig" is a wrapper around "eigfact()", extracting all parts of
   the factorization to a tuple; where possible, using "eigfact()"
   is recommended.
In [137]:
w0 = [1., 1]
w = copy(w0)

for i in 1:10
    w_new = M*w
    println(w_new)
end
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
[3.0,2.0]
In [138]:
for i =1:10
    t = 3
end
t
Out[138]:
3

Linear algebra

Julia has built-in linear algebra, not only using LAPACK, but now also generic routines that work for arbitrary element types, implemented completely in Julia.

For example, given a matrix $A$, the LU-decomposition of $A$ is equivalent to Gaussian elimination; it expresses $A$ as the product $A = LU$, with $L$ a lower-triangular and $U$ an upper-triangular matrix.

This is implemented in pure Julia for arbitrary element types. When the elements are standard floating-point numbers, it uses the corresponding fast LAPACK implementation.

In [139]:
M = rand(100, 100)
eig(M)
Out[139]:
(Complex{Float64}[50.0183+0.0im,1.84743+2.30961im,1.84743-2.30961im,-2.72534+0.880116im,-2.72534-0.880116im,-2.86948+0.0im,2.64375+0.964703im,2.64375-0.964703im,2.79553+0.0im,0.196765+2.72857im  …  -0.867929-0.457136im,0.478782+0.882849im,0.478782-0.882849im,-0.120951+0.762254im,-0.120951-0.762254im,-0.54004+0.0im,0.471338+0.0im,0.165157+0.118936im,0.165157-0.118936im,-0.288467+0.0im],
100x100 Array{Complex{Float64},2}:
  -0.101511+0.0im   -0.0115356-0.039612im    …   -0.0524037+0.0im
 -0.0952353+0.0im   0.00489796-0.00599155im      -0.0682674+0.0im
  -0.108918+0.0im    0.0161753-0.092386im         -0.114214+0.0im
  -0.111416+0.0im     0.102416+0.0351175im        -0.103224+0.0im
 -0.0977978+0.0im   -0.0198729+0.0722435im      -0.00218148+0.0im
 -0.0836507+0.0im     0.082544+0.0627261im   …    -0.036751+0.0im
 -0.0963353+0.0im   -0.0388372+0.0032277im        0.0954342+0.0im
 -0.0929018+0.0im     0.029581+0.0575031im        0.0794094+0.0im
  -0.108383+0.0im   -0.0832873-0.160539im         0.0634216+0.0im
 -0.0902722+0.0im   -0.0375187+0.0593948im       -0.0071044+0.0im
   -0.10374+0.0im    -0.102945+0.00231717im  …   -0.0538666+0.0im
  -0.106464+0.0im    -0.104489-0.00207164im      -0.0406458+0.0im
 -0.0960566+0.0im    0.0159221-0.049508im          -0.16083+0.0im
           ⋮                                 ⋱                   
 -0.0946263+0.0im   -0.0375218+0.0672119im          -0.1078+0.0im
 -0.0908183+0.0im  -0.00540686+0.142578im           0.12745+0.0im
  -0.104126+0.0im   -0.0670472+0.0553909im   …    -0.158563+0.0im
 -0.0999715+0.0im    0.0382579-0.0947829im       -0.0696019+0.0im
 -0.0986926+0.0im    -0.084612-0.12682im           0.262889+0.0im
 -0.0949404+0.0im   -0.0400257+0.023941im          -0.11199+0.0im
  -0.107486+0.0im   -0.0344824-0.0759473im         0.087478+0.0im
 -0.0944187+0.0im    -0.021001+0.0461968im   …    0.0110733+0.0im
 -0.0999898+0.0im    0.0392409-0.0264316im       0.00222877+0.0im
 -0.0991767+0.0im  -0.00947725+0.0518241im        0.0205667+0.0im
  -0.108745+0.0im    0.0184766-0.0401366im         0.232166+0.0im
 -0.0993031+0.0im    0.0796098-0.0182621im        0.0712234+0.0im)
In [140]:
M = rand(100, 100)
M2 = map(big, M)
Out[140]:
100x100 Array{BigFloat,2}:
 2.225506282223690401878002376179210841655731201171875e-01  …  6.56851621473254265737296009319834411144256591796875e-01 
 7.319635804294619685350653526256792247295379638671875e-01     9.672824011227445506477806702605448663234710693359375e-01
 5.340947103998476830355457423138432204723358154296875e-01     9.836664954220417644847884730552323162555694580078125e-01
 3.162661116152143403468244287068955600261688232421875e-01     6.32452881486862406035243111546151340007781982421875e-01 
 2.223050262746237848432429018430411815643310546875e-01        9.4728250126541180264894137508235871791839599609375e-02  
 9.375878604020722040246482720249332487583160400390625e-01  …  5.216604894602110409351780617726035416126251220703125e-01
 2.82637383973149969307314677280373871326446533203125e-02      7.1043767348334707634194273850880563259124755859375e-01  
 9.31414288307105042719058474176563322544097900390625e-01      6.843501947718078071147829177789390087127685546875e-01   
 9.08755575606785992448521938058547675609588623046875e-01      7.7541670796339889903947550919838249683380126953125e-02  
 9.90974279386318102069708402268588542938232421875e-01         3.446265247750661675496530733653344213962554931640625e-01
 3.62011301348933667298979344195686280727386474609375e-01   …  5.674886600100188527306954711093567311763763427734375e-01
 9.70427871977250333657138980925083160400390625e-01            1.845040577422647487537687993608415126800537109375e-01   
 9.0041572068816577711913851089775562286376953125e-02          2.86408155710278666816748227574862539768218994140625e-01 
 ⋮                                                          ⋱                                                           
 5.51151260187439984150614691316150128841400146484375e-01      3.306386894959143507577437048894353210926055908203125e-01
 3.497623439347330442927841431810520589351654052734375e-01     8.4233377002282594503412838093936443328857421875e-01     
 2.454565323093833573153688121237792074680328369140625e-01  …  3.436356112336635959536579321138560771942138671875e-01   
 9.3207496351350638263966175145469605922698974609375e-01       7.007260372885506516382747577154077589511871337890625e-01
 8.144309263423663569625432501197792589664459228515625e-01     3.333819991414610495183978855493478477001190185546875e-01
 7.57929320416632634760389919392764568328857421875e-02         6.67666265761114541277265743701718747615814208984375e-01 
 6.366199276186090827422958682291209697723388671875e-01        3.03926419592160268479119622497819364070892333984375e-01 
 5.72819224497360668379997150623239576816558837890625e-01   …  6.6274981934506360659042911720462143421173095703125e-01  
 9.56949160642551976962977278162725269794464111328125e-01      8.29108303980233163343882551998831331729888916015625e-01 
 7.970001203346741558419807915925048291683197021484375e-01     1.888462575949938493380386717035435140132904052734375e-01
 9.0568890895109799288320573396049439907073974609375e-01       4.966056054863334789928330792463384568691253662109375e-01
 2.048115155536647247203063670895062386989593505859375e-01     3.74601565881409026559367703157477080821990966796875e-01 
In [141]:
?map
Base.map(f, c...) -> collection

   Transform collection "c" by applying "f" to each element. For
   multiple collection arguments, apply "f" elementwise.

      julia> map((x) -> x * 2, [1, 2, 3])
      3-element Array{Int64,1}:
       2
       4
       6

      julia> map(+, [1, 2, 3], [10, 20, 30])
      3-element Array{Int64,1}:
       11
       22
       33
In [142]:
lu(M2)
Out[142]:
(
100x100 Array{BigFloat,2}:
 1e+00                                                                                 …  0e+00
 1.838523316026038069423294824104585502107622762920382948468461321049017107287582e-01     0e+00
 2.226652672640778428153578974542223057899272696712333628373710024305295178288089e-01     0e+00
 8.800271373546953526824682655928247523713656239133309029931077636194740548240004e-01     0e+00
 9.969711511643861050008990963774833601635694233722855443788571591033609498720476e-01     0e+00
 8.979138730555442894457974595685742086077359813056139849438793026237408391767279e-01  …  0e+00
 6.807588600102166833365124716208762430663471006732089118695641612340596096693672e-01     0e+00
 2.455829705129951553304067727592573195211437393497849255275277512239547428646326e-01     0e+00
 7.733884474442739333431587703752807939448278524827717594937087861232859025738641e-01     0e+00
 6.369478597346308029171973720850203624729145464284684008762430813593740043172259e-01     0e+00
 2.827829745702584616783062241061194670528384475163341040075224413128537647704774e-02  …  0e+00
 6.8916911836155407719497226896317385232855759890401728866826844091322063237215e-02       0e+00
 9.914847445449848621627477787268697917524046066273956341540063379389126840704449e-01     0e+00
 ⋮                                                                                     ⋱       
 5.917560795509736662262312292818648981295636865290101585946113867783160165303854e-01     0e+00
 2.730343465699633349451394482373398269247854760360330831523299677598397581623309e-01     0e+00
 9.039945797670267226102309624603231261343146722095674272138480977753157246173324e-01  …  0e+00
 9.462017626922501325170145941450688057814574666304378855494052732254706598746676e-01     0e+00
 7.922089379792833548468864680973299515938762464750620693002149619341420960071084e-01     0e+00
 8.91312804253419886178066534318762609663940375583837777475766448455788684217046e-01      0e+00
 2.792642705237073137742461219913623698283706391737491380867611670762794730952418e-05     0e+00
 7.573650294135219377704801284400430317640792346704970287945898791774355476300744e-01  …  0e+00
 4.340348981096360281012666931631459479552193700232478278658625428081540612698508e-01     0e+00
 8.425932003310413942243844696491498599262760249784090231673454808048725464374257e-01     0e+00
 9.574420989713105503105950862187424587135981099834153833267817582826759295038591e-01     0e+00
 9.937128828323402433568347457741518522388968395097111212977949553756464564214356e-01     1e+00,

100x100 Array{BigFloat,2}:
 9.994851507686071112601666754926554858684539794921875e-01  …   5.00779152589674314555168166407383978366851806640625e-02                            
 0e+00                                                          5.957614820766671895663417312682526372739723842260882343182715867242722586063138e-01
 0e+00                                                          4.469478342275836995433682687239315802830255734475567159873597681467845826134468e-01
 0e+00                                                          1.594225612955579740227438886545435994704893579284040834538915577684978835022275e-01
 0e+00                                                          3.51418102474667774543916344035850933636933755274101969213535739730486430086481e-01 
 0e+00                                                      …   3.853995991912961375481683208700128985955666356716160345936599990779107394373079e-01
 0e+00                                                         -1.388928850457794922700969846901469586577590670720174042695508511464725140648597e-01
 0e+00                                                         -1.963736253744149588710187786364289757425844420781275911168560619527357939135026e-01
 0e+00                                                          3.178968402329430072315582253427327119635802392272492696354540300940096734477185e-01
 0e+00                                                          2.176214988150537681248929128626937855865936488765697272966015579849791295536862e-01
 0e+00                                                      …   5.119997557073573321754864325074991553008298165056610101108332375190160965916915e-01
 0e+00                                                         -8.199420519824984131304944439473283677509641603892215267311260421539949314147843e-02
 0e+00                                                         -1.100757926861985762949819942253220494124734219500225416371630465493477225668268e+00
 ⋮                                                          ⋱                                                                                       
 0e+00                                                          3.119374583443352091823412074951378032702116370133229255033911885197311498449721e-01
 0e+00                                                         -3.963473902224899280333406215892276084125185726495321949218801950999434601486697e-01
 0e+00                                                      …   1.460022051064841242988743472828596892824793640959853693084700574539034836575632e+00
 0e+00                                                          9.460562261692564028602511319807472584567333376119182633775041922970959644819617e-01
 0e+00                                                          1.916100769270754349291340364996423402303124329778281922547647681942010509253265e+00
 0e+00                                                          2.453777141130431399782611832707715011759461721288957685494547893432872020437905e+00
 0e+00                                                          3.092241770787052786375610695487741355942081214332538992688100491955019679622272e-01
 0e+00                                                      …  -7.690085518435571170493252506836527889003361725238261885731324700815770035046439e-01
 0e+00                                                          7.057077397497909132059254043595925169793816268455912227899415362222622859241975e-01
 0e+00                                                          2.11408055146617577269295905987997866661481241577267775907335892526116912022093e-01 
 0e+00                                                          4.72135706225358566629488251105981838448104777282753072505691201455089121862885e-01 
 0e+00                                                         -1.411605504403213576549449087599331540683320578509496538094218931547453928045177e+00,

[41,78,1,52,36,63,62,91,21,95  …  67,34,76,69,31,77,28,82,97,58])
In [143]:
lu(M)
Out[143]:
(
100x100 Array{Float64,2}:
 1.0         0.0         0.0         0.0        …   0.0       0.0       0.0
 0.183852    1.0         0.0         0.0            0.0       0.0       0.0
 0.222665    0.333612    1.0         0.0            0.0       0.0       0.0
 0.880027    0.908157    0.485849    1.0            0.0       0.0       0.0
 0.996971    0.0121974   0.880719    0.482328       0.0       0.0       0.0
 0.897914    0.688716   -0.173968   -0.535249   …   0.0       0.0       0.0
 0.680759    0.256184    0.805169    0.0615926      0.0       0.0       0.0
 0.245583    0.122523    0.852098    0.230656       0.0       0.0       0.0
 0.773388    0.783385    0.294589    0.813216       0.0       0.0       0.0
 0.636948    0.621283   -0.464593    0.0427875      0.0       0.0       0.0
 0.0282783   0.895956   -0.28519    -0.723122   …   0.0       0.0       0.0
 0.0689169   0.529095    0.623271    0.567317       0.0       0.0       0.0
 0.991485    0.723194   -0.0744293   0.0926657      0.0       0.0       0.0
 ⋮                                              ⋱                          
 0.591756    0.553754    0.585172    0.792699       0.0       0.0       0.0
 0.273034    0.946579    0.312432    0.252827       0.0       0.0       0.0
 0.903995    0.190557    0.113748   -0.253964   …   0.0       0.0       0.0
 0.946202    0.644426    0.227142    0.306999       0.0       0.0       0.0
 0.792209    0.616671    0.267607    0.875396       0.0       0.0       0.0
 0.891313    0.895161    0.31949     0.0199976      0.0       0.0       0.0
 2.79264e-5  0.187737    0.926009    0.460432       0.0       0.0       0.0
 0.757365    0.538407    0.300644    0.140578   …   0.0       0.0       0.0
 0.434035    0.293604    0.46559     0.452693       0.0       0.0       0.0
 0.842593    0.325906    0.090703    0.620878       1.0       0.0       0.0
 0.957442    0.628242    0.0057032  -0.198115      -0.163927  1.0       0.0
 0.993713    0.465608   -0.431142   -0.125049      -0.617413  0.640633  1.0,

100x100 Array{Float64,2}:
 0.999485  0.0617704  0.296413  …   0.536738     0.579976    0.0500779
 0.0       0.964559   0.333478     -0.0592369   -0.0351511   0.595761 
 0.0       0.0        0.756701      0.751446     0.616732    0.446948 
 0.0       0.0        0.0          -0.572233    -0.401953    0.159423 
 0.0       0.0        0.0          -0.097152    -0.229394    0.351418 
 0.0       0.0        0.0       …   0.171652     0.103321    0.3854   
 0.0       0.0        0.0          -0.823456    -0.31806    -0.138893 
 0.0       0.0        0.0           0.346634     0.117767   -0.196374 
 0.0       0.0        0.0           0.114205     0.694677    0.317897 
 0.0       0.0        0.0           0.18621      0.454778    0.217621 
 0.0       0.0        0.0       …   0.875338     1.06396     0.512    
 0.0       0.0        0.0          -0.666298    -0.761167   -0.0819942
 0.0       0.0        0.0          -0.463261    -0.622826   -1.10076  
 ⋮                              ⋱                                     
 0.0       0.0        0.0           1.10803     -2.71365     0.311937 
 0.0       0.0        0.0           1.55302     -4.63257    -0.396347 
 0.0       0.0        0.0       …   0.00811131  -1.63731     1.46002  
 0.0       0.0        0.0           0.959823    -1.79163     0.946056 
 0.0       0.0        0.0           2.08064      0.0531956   1.9161   
 0.0       0.0        0.0           0.803113    -1.73451     2.45378  
 0.0       0.0        0.0          -1.11381      0.340109    0.309224 
 0.0       0.0        0.0       …   1.40266      0.277447   -0.769009 
 0.0       0.0        0.0          -1.36105      0.370272    0.705708 
 0.0       0.0        0.0          -1.64671      0.297193    0.211408 
 0.0       0.0        0.0           0.0          1.00489     0.472136 
 0.0       0.0        0.0           0.0          0.0        -1.41161  ,

[41,78,1,52,36,63,62,91,21,95  …  67,34,76,69,31,77,28,82,97,58])
In [144]:
?lu
Base.lu(A) -> L, U, p

   Compute the LU factorization of "A", such that "A[p,:] = L*U".
In [145]:
methods(lu)
Out[145]:
2 methods for generic function lu:
In [146]:
# @edit lu(M)
In [147]:
?edit
Base.edit(file::String[, line])

   Edit a file optionally providing a line number to edit at. Returns
   to the julia prompt when you quit the editor.

Base.edit(function[, types])

   Edit the definition of a function, optionally specifying a tuple of
   types to indicate which method to edit.

Interacting with the system

Command-line arguments

A Julia script, similar to a Python script, is a sequence of Julia commands placed in a file, with the termination .jl.

From the command line, a script script.jl can be run as

julia script.jl arg1 arg2 

where arg1 and arg2 are command-line arguments.

These command-line arguments to Julia scripts are placed in the variable ARGS as an array of strings.

Files

Simple file input and output is easy:

In [148]:
outfile = open("test.txt", "w")
Out[148]:
IOStream(<file test.txt>)
In [149]:
for i in 1:10
    println(outfile, "The value of i is $i")
end

close(outfile)
In [150]:
;cat test.txt
The value of i is 1
The value of i is 2
The value of i is 3
The value of i is 4
The value of i is 5
The value of i is 6
The value of i is 7
The value of i is 8
The value of i is 9
The value of i is 10
In [151]:
infile = open("test.txt", "r")
Out[151]:
IOStream(<file test.txt>)
In [152]:
lines = readlines(infile)
Out[152]:
10-element Array{Union(UTF8String,ASCIIString),1}:
 "The value of i is 1\n" 
 "The value of i is 2\n" 
 "The value of i is 3\n" 
 "The value of i is 4\n" 
 "The value of i is 5\n" 
 "The value of i is 6\n" 
 "The value of i is 7\n" 
 "The value of i is 8\n" 
 "The value of i is 9\n" 
 "The value of i is 10\n"
In [153]:
map(split, lines)
Out[153]:
10-element Array{Array{SubString{ASCIIString},1},1}:
 SubString{ASCIIString}["The","value","of","i","is","1"] 
 SubString{ASCIIString}["The","value","of","i","is","2"] 
 SubString{ASCIIString}["The","value","of","i","is","3"] 
 SubString{ASCIIString}["The","value","of","i","is","4"] 
 SubString{ASCIIString}["The","value","of","i","is","5"] 
 SubString{ASCIIString}["The","value","of","i","is","6"] 
 SubString{ASCIIString}["The","value","of","i","is","7"] 
 SubString{ASCIIString}["The","value","of","i","is","8"] 
 SubString{ASCIIString}["The","value","of","i","is","9"] 
 SubString{ASCIIString}["The","value","of","i","is","10"]
In [154]:
[float(line[6]) for line in map(split, lines)]
Out[154]:
10-element Array{Any,1}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0
In [155]:
x = rand(5,5)
Out[155]:
5x5 Array{Float64,2}:
 0.494374  0.548784    0.265116  0.14175    0.225112 
 0.610879  0.863226    0.90218   0.905241   0.0066329
 0.352839  0.163394    0.514312  0.46365    0.750558 
 0.72891   0.539688    0.37556   0.0172088  0.0674222
 0.069424  0.00747654  0.951597  0.984837   0.350702 
In [156]:
writedlm("random.txt", x)
In [157]:
;cat random.txt
.49437356567549995	.5487835147411582	.2651160558475778	.1417496876038895	.22511205695160474
.6108787919798604	.863226095224578	.9021796958394908	.9052411849611	.006632897854429665
.3528390706690816	.16339401668213926	.5143118341066373	.4636496583658838	.7505578231005474
.7289100209460473	.5396877767037889	.37555992341554356	.017208754116336156	.06742216219992869
.069423968121745	.007476542969398592	.9515968215572792	.9848367987309561	.35070195619498623
In [158]:
y = readdlm("random.txt")  # note that tab completion works for files
Out[158]:
5x5 Array{Float64,2}:
 0.494374  0.548784    0.265116  0.14175    0.225112 
 0.610879  0.863226    0.90218   0.905241   0.0066329
 0.352839  0.163394    0.514312  0.46365    0.750558 
 0.72891   0.539688    0.37556   0.0172088  0.0674222
 0.069424  0.00747654  0.951597  0.984837   0.350702 
In [159]:
?writedlm
Base.writedlm(f, A, delim='t')

   Write "A" (either an array type or an iterable collection of
   iterable rows) as text to "f" (either a filename string or an
   "IO" stream) using the given delimeter "delim" (which defaults
   to tab, but can be any printable Julia object, typically a "Char"
   or "String").

   For example, two vectors "x" and "y" of the same length can be
   written as two columns of tab-delimited text to "f" by either
   "writedlm(f, [x y])" or by "writedlm(f, zip(x, y))".
In [160]:
?readdlm
Base.readdlm(source, delim::Char, T::Type, eol::Char; header=false, skipstart=0, skipblanks=true, use_mmap, ignore_invalid_chars=false, quotes=true, dims, comments=true, comment_char='#')

   Read a matrix from the source where each line (separated by
   "eol") gives one row, with elements separated by the given
   delimeter. The source can be a text file, stream or byte array.
   Memory mapped files can be used by passing the byte array
   representation of the mapped segment as source.

   If "T" is a numeric type, the result is an array of that type,
   with any non-numeric elements as "NaN" for floating-point types,
   or zero. Other useful values of "T" include "ASCIIString",
   "String", and "Any".

   If "header" is "true", the first row of data will be read as
   header and the tuple "(data_cells, header_cells)" is returned
   instead of only "data_cells".

   Specifying "skipstart" will ignore the corresponding number of
   initial lines from the input.

   If "skipblanks" is "true", blank lines in the input will be
   ignored.

   If "use_mmap" is "true", the file specified by "source" is
   memory mapped for potential speedups. Default is "true" except on
   Windows. On Windows, you may want to specify "true" if the file
   is large, and is only read once and not written to.

   If "ignore_invalid_chars" is "true", bytes in "source" with
   invalid character encoding will be ignored. Otherwise an error is
   thrown indicating the offending character position.

   If "quotes" is "true", column enclosed within double-quote (``)
   characters are allowed to contain new lines and column delimiters.
   Double-quote characters within a quoted field must be escaped with
   another double-quote.

   Specifying "dims" as a tuple of the expected rows and columns
   (including header, if any) may speed up reading of large files.

   If "comments" is "true", lines beginning with "comment_char"
   and text following "comment_char" in any line are ignored.

Base.readdlm(source, delim::Char, eol::Char; options...)

   If all data is numeric, the result will be a numeric array. If some
   elements cannot be parsed as numbers, a cell array of numbers and
   strings is returned.

Base.readdlm(source, delim::Char, T::Type; options...)

   The end of line delimiter is taken as "\n".

Base.readdlm(source, delim::Char; options...)

   The end of line delimiter is taken as "\n". If all data is
   numeric, the result will be a numeric array. If some elements
   cannot be parsed as numbers, a cell array of numbers and strings is
   returned.

Base.readdlm(source, T::Type; options...)

   The columns are assumed to be separated by one or more whitespaces.
   The end of line delimiter is taken as "\n".

Base.readdlm(source; options...)

   The columns are assumed to be separated by one or more whitespaces.
   The end of line delimiter is taken as "\n". If all data is
   numeric, the result will be a numeric array. If some elements
   cannot be parsed as numbers, a cell array of numbers and strings is
   returned.

Shelling out

In [161]:
;ls
Example of defining a new type.ipynb
fractal_example
helllo_julia.png
Index.ipynb
Install IPython Notebook.md
Install Julia.md
Install Julia Packages.md
Interoperability with Python and C.ipynb
Introduction to Julia.ipynb
intro.txt
julia_tutorial.md
Metaprogramming.ipynb
Packages and graphics in Julia.ipynb
Performance in Julia.ipynb
Plan of the tutorial.ipynb
random.txt
README.md
stuff.pdf
test.txt
Vector2d.jl
What's next.ipynb
Why Julia.ipynb
Why Julia.slides.html
In [162]:
run(`echo Hello`)
Hello

Scientific computing

Linear algebra, FFT, random numbers, special functions. Packages for optimization, ODEs etc.

Functions

Functions may be defined using the short syntax f(x) = 3x + 1 or using a longer form:

In [163]:
dup(x) = 2x
Out[163]:
dup (generic function with 1 method)
In [164]:
function duplicate(x)
    2x   # no explicit "return" needed
end
Out[164]:
duplicate (generic function with 1 method)

The last value computed in the function is automatically returned; no explicit return statement is required.

In [165]:
duplicate(x) = x^2
Out[165]:
duplicate (generic function with 1 method)

Every operator in Julia is a function. Functions are implemented by specifying their action on different types. Until now, we have written only functions that are generic, in the sense that they do not specify which type they accept, and as in Python they will work as long as the operations performed in them make sense for the input value:

In [166]:
duplicate(3), duplicate(3.5), duplicate(1+3im)
Out[166]:
(9,12.25,-8 + 6im)
In [167]:
duplicate("Hola")
Out[167]:
"HolaHola"
In [168]:
2 * "Hola"
`*` has no method matching *(::Int64, ::ASCIIString)
while loading In[168], in expression starting on line 1

Note that string concatenation uses the * operator in Julia, instead of the + operator as in Python. Repeating a string is thus done by raising to an integer power:

In [169]:
"Hello"^2
Out[169]:
"HelloHello"

As a simple example, suppose that we wish to concatenate two strings. In Python we would write:

In [170]:
s1 = "Hello, "
s2 = "David"

s1 + s2
`+` has no method matching +(::ASCIIString, ::ASCIIString)
while loading In[170], in expression starting on line 4

However, we see that in Julia, summation is not defined for strings. What is it defined for?

In [171]:
+
Out[171]:
+ (generic function with 148 methods)

We see that + is treated as a function, and that it has a multitude of methods, which, in Julia, are specialised versions of the function that act on different types:

If we were unaware of the * operator for string concatenation, we could just define our own + for the concatenation of two strings:

In [173]:
+(s1::String, s2::String) = string(s1, s2)
Out[173]:
+ (generic function with 149 methods)
In [174]:
"First" + " second"
Out[174]:
"First second"

However, we cannot add a number to a string, since we have not (yet) defined it:

In [175]:
"The value of x is " + 3
`+` has no method matching +(::ASCIIString, ::Int64)
while loading In[175], in expression starting on line 1

This we can also define, using the previous new definition:

In [176]:
+(s::String, x::Number) = s + "$(2x)"
Out[176]:
+ (generic function with 150 methods)
In [177]:
"The value of x is " + 3
Out[177]:
"The value of x is 6"
In [178]:
x = 3.5
"The value of x is " + x
Out[178]:
"The value of x is 7.0"
In [179]:
3 + "hello"
`+` has no method matching +(::Int64, ::ASCIIString)
while loading In[179], in expression starting on line 1

In fact, we can define the summation of a string with any other object:

In [180]:
+(s::String, x) = s + string(x)
Out[180]:
+ (generic function with 151 methods)
In [181]:
"Complex " + [3,4,5]
Out[181]:
"Complex [3,4,5]"
In [182]:
"a" + 3
Out[182]:
"a6"
In [183]:
Number
Out[183]:
Number
In [184]:
typeof(Number)
Out[184]:
DataType
In [185]:
super(Int64)
Out[185]:
Signed
In [186]:
super(Signed)
Out[186]:
Integer
In [187]:
super(Integer)
Out[187]:
Real
In [188]:
super(Real)
Out[188]:
Number
In [189]:
super(Number)
Out[189]:
Any

In this way, the concept of "function" is replaced by a "patchwork" of different definitions for objects of different types, easily modifiable by the user. This is also exactly the way to define "operator overloading" for user-defined types.

In the above, we also begin to see the power of multiple dispatch: we defined two methods (versions) of the function +, both with the same number but different types of arguments.

User-defined types

A user-defined "composite type" is a collection of data. Unlike in Python, types do not "own" methods (functions internal to the type).

Rather, methods are defined separately, and are characterised by the types of all of their arguments; this is known as multiple dispatch. (Dispatch is the process of choosing which "version" of a given function to execute.)

A simple, but useful, example, is that of defining a 2D vector type. (See also the ImmutableArrays.jl package; fixed-size arrays will later be incorporated into base Julia.)

In [190]:
@which 3//4
Out[190]:
//(n::Integer,d::Integer) at rational.jl:15
In [191]:
Rational(3)
Out[191]:
3//1
In [192]:
6//4
Out[192]:
3//2
In [193]:
im*im
Out[193]:
-1 + 0im
In [194]:
immutable Vector2D   # type
    x::Float64
    y::Float64
end 
In [195]:
v = Vector2D(3, 4)
w = Vector2D(5, 6)
Out[195]:
Vector2D(5.0,6.0)
In [196]:
v + w
`+` has no method matching +(::Vector2D, ::Vector2D)
while loading In[196], in expression starting on line 1
In [197]:
+(v::Vector2D, w::Vector2D) = Vector2D(v.x+w.x, v.y+w.y)
Out[197]:
+ (generic function with 152 methods)
In [198]:
v + w
Out[198]:
Vector2D(8.0,10.0)
In [199]:
*(v::Vector2D, α::Number) = Vector2D(v.x*α, v.y*α)
*(α::Number, v::Vector2D) = Vector2D(v.x*α, v.y*α)
Out[199]:
* (generic function with 126 methods)
In [200]:
v * 3.5
Out[200]:
Vector2D(10.5,14.0)
In [201]:
3.5 * v
Out[201]:
Vector2D(10.5,14.0)

Exercise: Define mathematical operations on Vector2D. Define a particle with position and velocity in 2D. Define function move that acts on a particle to move it over a time $\delta t$.

In [202]:
immutable VectorT   # type
    p::Vector2D
    v::Vector2D
end
In [203]:
a = VectorT(Vector2D(0,0),Vector2D(2,3))
Out[203]:
VectorT(Vector2D(0.0,0.0),Vector2D(2.0,3.0))
In [204]:
move(a::VectorT, t::Number) = a.p+a.v*t
Out[204]:
move (generic function with 1 method)
In [205]:
println("new position of particle $a with origin position $(a.p) and velocity $(a.v) at time 3 is $(move(a, 3))")
new position of particle VectorT(Vector2D(0.0,0.0),Vector2D(2.0,3.0)) with origin position Vector2D(0.0,0.0) and velocity Vector2D(2.0,3.0) at time 3 is Vector2D(6.0,9.0)

Here, we have used immutable instead of type for efficiency: the object is stored in an efficient packed form.

The equivalent of the Python __repr__ method for an object is to extend the show method:

In [206]:
?show
Base.show(x)

   Write an informative text representation of a value to the current
   output stream. New types should overload "show(io, x)" where the
   first argument is a stream. The representation used by "show"
   generally includes Julia-specific formatting and type information.
In [207]:
import Base.show

show(io::IO, v::Vector2D) = print(io, "[$(v.x), $(v.y)]")
Out[207]:
show (generic function with 94 methods)
In [208]:
v
Out[208]:
[3.0, 4.0]
In [209]:
+(v1::Vector2D, v2::Vector2D) = Vector2D(v1.x+v2.x, v1.y+v2.y)

*(v::Vector2D, lamb::Number)  = Vector2D(lamb*v.x, lamb*v.y)
Out[209]:
* (generic function with 126 methods)

We can confirm that the new method for the function + has indeed been defined:

In [211]:
x = Vector2D(3, 4)
y = Vector2D(5, 6)

x + y
Out[211]:
[8.0, 10.0]

Parametrised types

Types may have a parameter, for example:

In [212]:
immutable Vector2D{T <: Real} #must same type T
    x::T
    y::T
end 
invalid redefinition of constant Vector2D
while loading In[212], in expression starting on line 1

T is a type parameter. The expression T <: Real means that T must be a subtype of the abstract type Real. We can investigate the type hierarchy with the super function:

In [213]:
Integer
Out[213]:
Integer
In [214]:
super(Integer)
Out[214]:
Real
In [215]:
super(Real)
Out[215]:
Number
In [216]:
super(Number)
Out[216]:
Any
In [217]:
v = Vector2D(3., 4.)
Out[217]:
[3.0, 4.0]
In [218]:
v = Vector2D(3, 4.)
Out[218]:
[3.0, 4.0]

Here, the types of the two arguments were different, so there is no match for the type signature.

In [219]:
v = Vector2D(3//4, 5//6)
Out[219]:
[0.75, 0.8333333333333334]
In [220]:
import Base.show
show{T}(io::IO, v::Vector2D{T}) = print(io, "[$(v.x), $(v.y)]")
too many parameters for type Vector2D
while loading In[220], in expression starting on line 2
In [221]:
v
Out[221]:
[0.75, 0.8333333333333334]

We can define outer constructors, defined outside the type definition itself, which allow other ways of constructing the object:

In [222]:
Vector2D{T}(x::T) = Vector2D(x, x)
Out[222]:
Vector2D (constructor with 3 methods)
In [223]:
Vector2D(3)
Out[223]:
[3.0, 3.0]

Example: a simple type for a collection of particles

Let's define a particle:

In [224]:
type Particle
    position::Vector2D{Float64}
    velocity::Vector2D{Float64}
end
move(p::Particle, dt::Real) = p.position += p.velocity * dt
too many parameters for type Vector2D
while loading In[224], in expression starting on line 1
In [225]:
show(io::IO, p::Particle) = print(io, "pos: $(p.position); vel: $(p.velocity)")
Particle not defined
while loading In[225], in expression starting on line 1
In [226]:
p = Particle(Vector2D(0.,0.), Vector2D(1.,1.))
+{T}(v1::Vector2D{T}, v2::Vector2D{T}) = Vector2D{T}(v1.x+v2.x, v1.y+v2.y)

*{T}(v::Vector2D{T}, lamb::Number)  = Vector2D{T}(lamb*v.x, lamb*v.y)
Particle not defined
while loading In[226], in expression starting on line 1
In [227]:
move(p, 0.1)
p not defined
while loading In[227], in expression starting on line 1
In [228]:
p
p not defined
while loading In[228], in expression starting on line 1

Now we can define a gas as a collection of particle:

In [229]:
Int
Out[229]:
Int64
In [230]:
type Gas
    particles::Vector{Particle}  # Array{Particle, 1}
    
    function Gas(N::Integer)
        parts = [Particle(Vector2D(rand(2)...), Vector2D(rand(2)...)) for i in 1:N]
        new(parts)
    end
end
Particle not defined
while loading In[230], in expression starting on line 1
In [231]:
show(io::IO, g::Gas) = for i in 1:length(g.particles); \
    println(io, "Particle $i: $(g.particles[i])"); end
Gas not defined
while loading In[231], in expression starting on line 1
In [232]:
g = Gas(10)
Gas not defined
while loading In[232], in expression starting on line 1
In [233]:
move(g, 1)
g
`move` has no method matching move(::Function, ::Int64)
while loading In[233], in expression starting on line 1
In [234]:
function move(g::Gas, dt::Number)
    for particle in g.particles
        move(particle, dt)
    end
end
Gas not defined
while loading In[234], in expression starting on line 1
In [235]:
move(g, 1)
`move` has no method matching move(::Function, ::Int64)
while loading In[235], in expression starting on line 1
In [236]:
g
Out[236]:
g (generic function with 1 method)

Matrix multplication

In [237]:
M = [1 1; 0 1]
Out[237]:
2x2 Array{Int64,2}:
 1  1
 0  1
In [238]:
M * M
Out[238]:
2x2 Array{Int64,2}:
 1  2
 0  1
In [239]:
M^10
Out[239]:
2x2 Array{Int64,2}:
 1  10
 0   1
In [240]:
M^(-1)
Out[240]:
2x2 Array{Float64,2}:
 1.0  -1.0
 0.0   1.0
In [241]:
inv(M)
Out[241]:
2x2 Array{Float64,2}:
 1.0  -1.0
 0.0   1.0
In [242]:
?inv
Base.inv(M)

   Matrix inverse
In [243]:
det(M)
Out[243]:
1
In [244]:
norm(M)
Out[244]:
1.618033988749895
In [245]:
?norm
Base.norm(A[, p])

   Compute the "p"-norm of a vector or the operator norm of a matrix
   "A", defaulting to the "p=2"-norm.

   For vectors, "p" can assume any numeric value (even though not
   all values produce a mathematically valid vector norm). In
   particular, "norm(A, Inf)" returns the largest value in
   "abs(A)", whereas "norm(A, -Inf)" returns the smallest.

   For matrices, valid values of "p" are "1", "2", or "Inf".
   (Note that for sparse matrices, "p=2" is currently not
   implemented.) Use "vecnorm()" to compute the Frobenius norm.
In [246]:
lamb, vv = eig(M)
Out[246]:
([1.0,1.0],
2x2 Array{Float64,2}:
 1.0  -1.0        
 0.0   2.22045e-16)
In [247]:
Symmetric(M)
Out[247]:
2x2 Symmetric{Int64,Array{Int64,2}}:
 1  1
 1  1
In [248]:
?Symmetric
DataType   : Symmetric{T,S<:AbstractArray{T,2}} (constructor with 2 methods)
  supertype: AbstractArray{T,2}
  fields   : (:data,:uplo)
In [249]:
methods(Symmetric)
Out[249]:
2 methods for generic function Symmetric:

Eigenvalues of a random matrix


Exercise: Generate a random matrix of $1000 \times 1000$ gaussian random variates. Calculate its eigenvalues and plot them. Calculate the differences between consecutive eigenvalues and plot a histogram of them.


In [250]:
N = 1000
M = randn(N, N)
M = Symmetric(M);
In [251]:
typeof(M)
Out[251]:
Symmetric{Float64,Array{Float64,2}} (constructor with 1 method)
In [252]:
@which dot([3], [4])
Out[252]:
dot(x::AbstractArray{T,1},y::AbstractArray{T,1}) at linalg/matmul.jl:49
In [253]:
help("dot")
Base.LinAlg.BLAS.dot(n, X, incx, Y, incy)

   Dot product of two vectors consisting of "n" elements of array
   "X" with stride "incx" and "n" elements of array "Y" with
   stride "incy".
In [254]:
@time lamb, vv = eig(M);
elapsed time: 0.730055745 seconds (24964608 bytes allocated, 6.74% gc time)
In [255]:
using PyPlot
In [256]:
plot(lamb, "o-")
Out[256]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7ff8f8bd3390>
In [257]:
PyPlot.figure(figsize=(8,4))
plot(lamb, "o-")
Out[257]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7ff8f8a19310>
In [258]:
differences = diff(lamb);
In [259]:
plot(differences, "o-")
Out[259]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7ff8f89d3650>
In [260]:
h = hist(differences, 100)
Out[260]:
(0.0:0.01:0.93,[8,16,23,48,47,54,59,56,54,66  …  0,0,0,0,0,0,0,0,0,1])
In [261]:
plot(collect(h[1][1:end-1]), h[2])
Out[261]:
1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7ff8f88ed850>

requires vs include vs using vs import

In [262]:
g.particles
type Function has no field particles
while loading In[262], in expression starting on line 1
In [263]:
for particle in g.particles
    move(particle, 1)
end

g.particles
type Function has no field particles
while loading In[263], in expression starting on line 1

 in anonymous at no file