Revision as of 00:08, 12 February 2021

Crash Course in Scientific Computing

IV. Julia

Julia is a powerful/efficient/high-level computer programming language. You can get into interacting mode right-away with:

julia

Exit with exit(). One may also also use Wolfram's notebook concept:

jupyter-notebook

We will often need to install packages, which, back to julia, is done as follows:

using Pkg
Pkg.add("IJulia")

Then, from a terminal, jupyter-notebook opens a session in a browser, in which one can start a Julia session with a New invocation (and choosing julia there).

There are some conventions in computer programming which become natural as you get used to them. For instance, $x=3$ takes a different meaning in most languages than in Mathematics. In the latter case, it means that the variable $x$ has the same value as the number 3. For a computer x=3 assigns the value 3 to the variable x. In Julia that becomes:

julia> x=3
3

The output 3 is the result of our line of code, which reads: assigns the value 3 to x. The result of this is 3. In maths, if $x=3$, then $x\neq 4$, while in computer code, following x=3, we can have x=4 meaning we now assign the value 4 to x. The computer code equivalent of the mathematical $=$ is actually ==. Indeed:

julia> x=4
4

julia> x==3
false

A typical line of code is to increment the value of a variable, which in many cases we would refer to as a counter:

x=x+1

Note that in Maths, this equality is trivially wrong. In computer code, on the other hand, it is so useful that we even have a special notation for it:

x+=1

Incrementing the value of x and storing the result:

julia> x=4
4

julia> x+=1
5

A very early computer language which defined much of the conventions, notations and overall syntax of computer programming is C. An improved (and modernized) version implementing object-programming is known as C++, implying the joke it's an added version of C.

Most conventions are natural or familiar already

julia> 3>2
true

julia> 17<=5
false

julia> 1+5!=2+3
true

but some others, we're not even sure we really want them:

julia> 7/2==2\7
true

(this is called inverse division). More useful is the remainder $r$ from Euclidean division $a=bq+r$ which is obtained as a%b, for instance, $17=3\times 5+2$ so

julia> 17%5
2

julia> 17%3
2

so for instance, when we found that it was taking 205s for the computer of a dumb programmer to compute a Fibonacci number, we could give a better idea of the poor performance by giving the time in minutes:

julia> 205÷60
3

julia> 205%60
25

i.e., 205s=3min25s. We entered the integer division ÷ with \div+TAB. If this is too esoteric, we could have used floor(204/60) where floor extracts the integer part (to lowest value; we'd use ceil for the highest value and round for the closest value).

Julia also has basic symbolic abilities, beyond numerical ones:

julia> 24/9
2.6666666666666665

julia> 24//9
8//3

julia> numerator(1389//234)
463

Julia knows about important constants:

julia> pi
π = 3.1415926535897...

It even supports unicode. We could use such characters using LaTeX encoding and tab:

\alpha+[tab]→α

that we can use as a variable. In some cases it could be useful, compare for instance:

julia> 5≥2
true

julia> 5>=2
true

julia> 5=>2
5 => 2

where the top version was input as \ge+[tab]. If you know the unicode code, you can enter it as "\uXXXX". For instance:

julia> c="\u2020"
"†"

By the way, we have just introduced a new type of "variable", that we don't find in maths. So-called strings:

julia> hello="Hi there!"
"Hi there!"

julia> hello
"Hi there!"

You can somehow "operate" on such things as well. Here's concatenation, or joining strings:

julia> hello*" Hey, you!"
"Hi there! Hey, you!"

Strings are clearly important in computer programming, as we deal a lot with texts. We shall see concepts such as regular expressions that allow powerful text processing, e.g.,

julia> replace("Julia", "a" => "us")
"Julius"

But for now, we'll be more concerned on the numerical aspect. In this context, another type of variable that is very useful is the concept of arrays, which is a list or vector. Here's an array of 10 consecutive integers:

julia> a = [1:10;]
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

You can construct arrays explictely:

julia> a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

or have it filled-in with a given value:

julia> z=zeros(10)
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0

or ones(n) for n-sized array full of 1.0 (you can multiply this by something to be constant to some other value or uses fill instead).

Things can quickly go from completely obvious to quite tricky. Consider the following:

julia> b=a
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

julia> b[5]+=1
6

where b[5] looks at the 5th value of b, so we have simply incremented the 5th element of b by one, so that:

julia> b
10-element Array{Int64,1}:
  1
  2
  3
  4
  6
  6
  7
  8
  9
 10

Now for the tricky part. what do you think happens with:

a==b

The answer is true. This is counterintuitive because with "variables" rather than arrays, we get a different results:

julia> α=5
5

julia> β=α
5

julia> β+=1
6

julia> α==β
false

This is because the command b=a does not copy the array a into b (it does for the variables), but it copied the "address" of the array in memory. If we want to duplicate an independent array, we must do:

ulia> b=copy(a)
10-element Array{Int64,1}:
  1
  2
  3
  4
  7
  6
  7
  8
  9
 10

Now changing b[5] does not change a[5]. Be careful of these kinds of subtle things, they are the source of many mistakes (or "bugs", as we call them).

Another subtlety is that [1:10] (without the ;) returns still something else:

julia> [1:10]
1-element Array{UnitRange{Int64},1}:
 1:10

which is an array with one element, namely, 1:10, the latter being a so-called range, which appears in other constructs, such as:

julia> collect(1:10)
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

that builds a proper array! collect and ranges are good to make more interesting grids, e.g.:

julia> collect(0:10:100)
11-element Array{Int64,1}:
   0
  10
  20
  30
  40
  50
  60
  70
  80
  90
 100

julia> collect(0:.25:1)
5-element Array{Float64,1}:
 0.0 
 0.25
 0.5 
 0.75
 1.0 

Let us start bringing in the famous control-flow commands. Here's the list of the first 10 squares:

julia> [i^2 for i in 1:10]
10-element Array{Int64,1}:
   1
   4
   9
  16
  25
  36
  49
  64
  81
 100

The result of the last evaluation is ans. If we want to know which is the 7th square, we could do:

julia> ls=ans;

julia> ls[7]
49

where we called ls the "list of squares". Note that ; suppress printing the output, in case you're not interested in seeing it.

One can add a if statement to impose a test in the array construction. This returns the squares less than 100 divisible by 3:

julia> [i^2 for i=1:10 if i % 3 == 0]
3-element Array{Int64,1}:
  9
 36
 81

Arrays are basically what we call "vectors" in Mathematics. The corresponding matrices are two-dimensional arrays. Compare:

julia> [1, 2, 3, 4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> [1 2 3 4]
1×4 Array{Int64,2}:
 1  2  3  4

julia> [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

2D arrays can be built with two iterators. Here is how to input the matrix $M$ such that $M_{ij}=ij$:

julia> [i*j for i in 1:5, j in 1:5]
5×5 Array{Int64,2}:
 1   2   3   4   5
 2   4   6   8  10
 3   6   9  12  15
 4   8  12  16  20
 5  10  15  20  25

A zero $2\times 3$ matrix:

julia> zeros(2,3)
2×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0

A powerful command allows to reshape arrays from one type into another:

julia> reshape([1:9;],3,3)
3×3 Array{Int64,2}:
 1  4  7
 2  5  8
 3  6  9

julia> reshape([1:18;],3,6)
3×6 Array{Int64,2}:
 1  4  7  10  13  16
 2  5  8  11  14  17
 3  6  9  12  15  18

julia> reshape([1:18;],6,3)
6×3 Array{Int64,2}:
 1   7  13
 2   8  14
 3   9  15
 4  10  16
 5  11  17
 6  12  18

We'll keep the following generalization of the idea (tensors) for later:

julia> reshape([1:27;],3,3,3)
3×3×3 Array{Int64,3}:
[:, :, 1] =
 1  4  7
 2  5  8
 3  6  9

[:, :, 2] =
 10  13  16
 11  14  17
 12  15  18

[:, :, 3] =
 19  22  25
 20  23  26
 21  24  27

even though you surely have a pretty good idea of what is going on.

One can also work with arrays of arrays:

julia> A=[[1, 2], [3, 4]]
2-element Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]

When dealing with arrays, we are often interested in accessing parts of the array. We have seen already a[5] for the 5th value. We can also access from the end:

julia> a=[1:10;]
10-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10

julia> a[end]
10

julia> a[end-3]
7

To do beyond one value at a time, we asks for ranges:

julia> a[3:2:end-2]
3-element Array{Int64,1}:
 3
 5
 7

since 3:2:end-2 refers to the range from 3 till the 2nd-to-last items by steps of 2. Consecutive numbers are boring, so let's work with their square. Instead of creating the array directly (like ls above), let's keep using a but apply a command to it to process each input from the array. We would think of something like a^2 but we need something else, known as "broadcasting" and that involves a . to force the operation to be distributed to all elements of the array, so we'd go:

julia> a.^2
10-element Array{Int64,1}:
   1
   4
   9
  16
  25
  36
  49
  64
  81
 100

Now we can query parts of this, for instance, if we want the 3rd, 5th and 9th elements, we'd use ...:

julia> (a.^2)[[3,5,9]]
3-element Array{Int64,1}:
  9
 25
 81

We can use a length-matching sequence of true & false:

julia> (a.^2)[[true, true, true, false, false, false, false, true, true, true]]
6-element Array{Int64,1}:
   1
   4
   9
  64
  81
 100

with which we can achieve more interest results such as:

julia> (a.^2)[[i%4==0 for i in 1:10]]
2-element Array{Int64,1}:
 16
 64

that finds all the squares less than 100 that are divisible by 4.

julia> Base.MathConstants.golden
φ = 1.6180339887498...

Incrementing a counter is a good way to keep track of what is going, along with control flow. The most famous are:

Let us now play with things computers are good at, to start with, plotting: