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).

For some time, we will use the console and work interactively in what is known as the REPL (read–eval–print loop) which is a sort of super pocket-calculator where the computer attends your every command on the spot, e.g.,

julia> 2+2
4

So let's start talking to it! To do so we need a common language. Computers natively don't speak English so it is for us to learn their language. The lowest-level language (assembly language) is not very convenient for us, although that's what the computer understands best and execute most efficiently. Instead, we use intermediate languages which need to be "compiled" into the computer's mother tongue (assembly). There are many languages, some very famous, like C or Java, others less so like golang or Brainfuck (you can have fun reading at the list of the most esoteric programming languages, such as Whitespace which consists only of whitespace, i.e., spaces, tabs and linebreaks, making it particularly difficult to read!). Some are very high-level and straightforward like Mathematica (which we have used last time) or Matlab, others are more cryptic and technically involved, like C++. We'll use Julia as a great compromise between simplicity and power. This is very close to other popular languages like python so you'll be able to go from one to the other easily (would python still exist in the future as Julia might take over).

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...

Although sometimes you have to look for them. This is the golden ratio $\phi$:

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

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 with 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). One can also copy the "structure" only, i.e., the same format:

julia> c=similar(a)
10-element Array{Int64,1}:
 140612253716432
 140612233906160
 140612233906160
 140612327663648
 140612253721936
 140612363299584
 140612327663648
 140612327664096
 140612233904704
 140612328019200

The content is just what happened to be there in memory.

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

One can also make 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.

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.

We have seen how matrices multiply themselves in a particular way. If $A=BC$ then

$$A_{ij}=\sum_{k=1}^lB_{ik}C_{jk}\,,$$

demanding in particular that the sizes of the corresponding matrices obey:

$$\label{eq:dimfits}(m\times n)=(m\times l)(l\times m)$$

Julia supports matrix multiplication as products of arrays. Let us check:

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

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

julia> A=B*C
3×4 Array{Int64,2}:
  9  19  29  39
 12  26  40  54
 15  33  51  69

Multiplying a vector by a matrix becomes a simple particular case.

If we try the other way around, we'll make our first acquaintance with a recurrent theme of computer programming: the computer error:

julia> C*B
ERROR: DimensionMismatch("matrix A has dimensions (2,4), matrix B has dimensions (3,2)")
Stacktrace:
 [1] _generic_matmatmul!(::Array{Int64,2}, ::Char, ::Char, ::Array{Int64,2}, ::Array{Int64,2}) at /build/julia-WtZlcm/julia-1.0.4+dfsg/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:591
 [2] generic_matmatmul!(::Array{Int64,2}, ::Char, ::Char, ::Array{Int64,2}, ::Array{Int64,2}) at /build/julia-WtZlcm/julia-1.0.4+dfsg/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:581
 [3] mul! at /build/julia-WtZlcm/julia-1.0.4+dfsg/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:175 [inlined]
 [4] *(::Array{Int64,2}, ::Array{Int64,2}) at /build/julia-WtZlcm/julia-1.0.4+dfsg/usr/share/julia/stdlib/v1.0/LinearAlgebra/src/matmul.jl:142
 [5] top-level scope at none:0

It is instructive to get hints from the nasty soup of details that is echoed in the process. What happened of course is that we did not satisfy the criterion Eq. \ref{eq:dimfits}.

If we want to make elementwise multiplication (Hadamard product), then we need to broadcast the multiplication so that it distributes over the arrays:

julia> A=reshape([i^2 for i in 1:9],3,3)
3×3 Array{Int64,2}:
 1  16  49
 4  25  64
 9  36  81

julia> B=reshape([1/i for i in 1:9],3,3)
3×3 Array{Float64,2}:
 1.0       0.25      0.142857
 0.5       0.2       0.125   
 0.333333  0.166667  0.111111

julia> A*B
3×3 Array{Float64,2}:
 25.3333  11.6167   7.5873
 37.8333  16.6667  10.8075
 54.0     22.95    14.7857

julia> A.*B
3×3 Array{Float64,2}:
 1.0  4.0  7.0
 2.0  5.0  8.0
 3.0  6.0  9.0

Here it is easy to check one of the most important results of matrix algebra:

julia> A*B==B*A
false

julia> A.*B==B.*A
true

If we want to do more advanced matrix operations, we need to import them as they are not always needed or are not by everybody and it would be unpractical to have the language cluttered with advanced functionalities. The idea of empowering the language with further tools brings us to "libraries". In the cases of matrices, we would "import" the LinearAlgebra library:

julia> using LinearAlgebra

julia> det(A)
-216.00000000000034

julia> inv(A)
3×3 Array{Float64,2}:
  1.29167  -2.16667   0.930556
 -1.16667   1.66667  -0.611111
  0.375    -0.5       0.180556

That's computing the inverse and determinant. Be careful with computers, they can be missing the point:

julia> det(A*B)==det(B*A)
false

which is unexpected because we have proven in class that the determinant commutes as a consequence of $\det(AB)=\det(A)\det(B)$. Well, this is not even true from the computer's perspective, but here's why:

julia> detA=det(A)
-216.00000000000034

julia> detB=det(B)
0.00029761904761904656

julia> det(A*B)
-0.06428571428571417

julia> detA*detB
-0.06428571428571415

Do you see where it fails? Remember: the computer is literal, no subtlety, no understanding of the approximations and of course, it is cluttered with numerical errors. So be wary of them when dealing with machine precision. Of course in principle you can always instruct them better and have them work in so-called infinite precision (symbolically):

julia> A=reshape([i^2//1 for i in 1:9],3,3)
3×3 Array{Rational{Int64},2}:
 1//1  16//1  49//1
 4//1  25//1  64//1
 9//1  36//1  81//1

julia> B=reshape([1//i for i in 1:9],3,3)
3×3 Array{Rational{Int64},2}:
 1//1  1//4  1//7
 1//2  1//5  1//8
 1//3  1//6  1//9

julia> A*B
3×3 Array{Rational{Int64},2}:
  76//3  697//60   478//63 
 227//6   50//3   5447//504
  54//1  459//20   207//14 

julia> B*A
3×3 Array{Rational{Int64},2}:
 23//7   767//28   536//7 
 97//40   35//2   1897//40
  2//1    27//2     36//1 

julia> det(A*B)
-9//140

julia> [det(A), det(B)]
2-element Array{Rational{Int64},1}:
 -216//1   
    1//3360

julia> det(A*B)==det(A)det(B)
true

We can access the elements of an array by specifying its indices:

julia> A=reshape([i^2 for i=1:9], 3, 3)
3×3 Array{Int64,2}:
 1  16  49
 4  25  64
 9  36  81

julia> A[1,2]
16

julia> A[2,1]
4

We use a colon to retrieve every row or column. For example, here's "every row, second column" and vice-versa:

julia> A[:, 2]
3-element Array{Int64,1}:
 16
 25
 36

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

You can also make more sophisticated selections:

julia> A[[1, 3], [3, 2, 1]]
2×3 Array{Int64,2}:
 49  16  1
 81  36  9

julia> A[[1, 3], [1, 2, 3]]
2×3 Array{Int64,2}:
 1  16  49
 9  36  81

We'll leave further explorations and usage to such things for the future, introducing them as they'll be needed, and now turn to the most famous feature of computer programming: control-flow. This is deciding on what to do depending on the state of the system, mainly, of the variables. This is a possible implementation of the Fibonacci sequence:

julia> ngen=25;

julia> fibo=zeros(ngen);

julia> fibo[1]=fibo[2]=1;

julia> for i=3:ngen
         fibo[i]=fibo[i-1]+fibo[i-2]
       end

julia> fibo
25-element Array{Float64,1}:
     1.0
     1.0
     2.0
     3.0
     5.0
     8.0
    13.0
    21.0
    34.0
    55.0
    89.0
   144.0
   233.0
   377.0
   610.0
   987.0
  1597.0
  2584.0
  4181.0
  6765.0
 10946.0
 17711.0
 28657.0
 46368.0
 75025.0

We are merely iterating over the variable which serves as a counter in this case. Euler noticed something, the ratio of subsequent Fibonacci numbers tends to a famous constant:

julia> fibo./circshift(fibo,1)
25-element Array{Float64,1}:
 1.3328890369876708e-5
 1.0                  
 2.0                  
 1.5                  
 1.6666666666666667   
 1.6                  
 1.625                
 1.6153846153846154   
 1.619047619047619    
 1.6176470588235294   
 1.6181818181818182   
 1.6179775280898876   
 1.6180555555555556   
 1.6180257510729614   
 1.6180371352785146   
 1.618032786885246    
 1.618034447821682    
 1.6180338134001253   
 1.618034055727554    
 1.6180339631667064   
 1.6180339985218033   
 1.618033985017358    
 1.6180339901755971   
 1.618033988205325    
 1.618033988957902 

where we used circshift to rotate the array. Let us now look at Collatz sequence which relies even more importantly on control flow, and for which we can use an if-else statement:

julia> n=5;

julia> if n%2==0
           n/=2
       else
           n=3n+1
       end
16

and further iterations in the REPL give us: 8, 4, 2, 1. This can be iterated by the computer itself until the conjecture is fulfilled. We need to define variables within the loops as "global" so that they are identified by the machine as part of the loop:

Julia> n=15;counter=0

julia> while n!=1
           global counter+=1
           print("$(n)\u2192")
           if n%2==0
               global n÷=2
           else
               n=3n+1
           end
       end
       println("1.")

15→46→23→70→35→106→53→160→80→40→20→10→5→16→8→4→2→1.

Note the use of print where piece of code enter in the escaping $() sequence.

To run over all instances, it is useful to package our piece of code in a so-called function, while defines a piece of code that can take some input and actuates on it, returning a value that we can further process. This is a possible implementation of such a function:

julia> function collatz(n)
           global counter=0;
           while n!=1
               counter+=1;
               print("$(n)\u2192")
               if n%2==0
                   n÷=2
               else
                   n=3n+1
               end
           end
           println("1. ($(counter) steps)")
       end
collatz (generic function with 1 method)

julia> collatz(24)
24→12→6→3→10→5→16→8→4→2→1. (10 steps)

which is now easy to iterate:

julia> for i=1:30
           collatz(i)
       end

1. (0 steps)
2→1. (1 steps)
3→10→5→16→8→4→2→1. (7 steps)
4→2→1. (2 steps)
5→16→8→4→2→1. (5 steps)
6→3→10→5→16→8→4→2→1. (8 steps)
7→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (16 steps)
8→4→2→1. (3 steps)
9→28→14→7→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (19 steps)
10→5→16→8→4→2→1. (6 steps)
11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (14 steps)
12→6→3→10→5→16→8→4→2→1. (9 steps)
13→40→20→10→5→16→8→4→2→1. (9 steps)
14→7→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (17 steps)
15→46→23→70→35→106→53→160→80→40→20→10→5→16→8→4→2→1. (17 steps)
16→8→4→2→1. (4 steps)
17→52→26→13→40→20→10→5→16→8→4→2→1. (12 steps)
18→9→28→14→7→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (20 steps)
19→58→29→88→44→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (20 steps)
20→10→5→16→8→4→2→1. (7 steps)
21→64→32→16→8→4→2→1. (7 steps)
22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (15 steps)
23→70→35→106→53→160→80→40→20→10→5→16→8→4→2→1. (15 steps)
24→12→6→3→10→5→16→8→4→2→1. (10 steps)
25→76→38→19→58→29→88→44→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (23 steps)
26→13→40→20→10→5→16→8→4→2→1. (10 steps)
27→82→41→124→62→31→94→47→142→71→214→107→322→161→484→242→121→364→182→91→274→137→412→206→103→310→155→466→233→700→350→175→526→263→790→395→1186→593→1780→890→445→1336→668→334→167→502→251→754→377→1132→566→283→850→425→1276→638→319→958→479→1438→719→2158→1079→3238→1619→4858→2429→7288→3644→1822→911→2734→1367→4102→2051→6154→3077→9232→4616→2308→1154→577→1732→866→433→1300→650→325→976→488→244→122→61→184→92→46→23→70→35→106→53→160→80→40→20→10→5→16→8→4→2→1. (111 steps)
28→14→7→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (18 steps)
29→88→44→22→11→34→17→52→26→13→40→20→10→5→16→8→4→2→1. (18 steps)
30→15→46→23→70→35→106→53→160→80→40→20→10→5→16→8→4→2→1. (18 steps)

Notice that something happens with 27!

We conclude with Conway's game of Life and this time we'll define functions to keep things handy. cleanspace initializes (or reset, if already created) the space of our cellular automaton.

julia> function cleanspace(size)
           global space=[0 for i=1:size, j=1:size]
       end
cleanspace (generic function with 1 method)

julia> cleanspace(12)
12×12 Array{Int64,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

julia> function randomizespace()
           for i=1:size(space)[1]
               for j=1:size(space)[2]
                   space[i,j]=rand([0,1])
               end
           end
       end
randomizespace (generic function with 1 method)

We now have a space to work on:

julia> space
12×12 Array{Int64,2}:
 0  0  1  0  1  1  1  0  1  0  1  1
 1  1  0  1  0  1  0  0  1  0  0  1
 0  1  1  1  1  0  0  0  1  1  0  1
 0  0  0  1  1  0  1  1  0  1  1  0
 1  1  1  1  0  0  1  0  1  1  0  0
 1  0  1  1  0  1  1  1  1  1  1  1
 0  0  0  1  1  0  0  0  0  0  1  0
 1  1  0  0  1  1  1  0  1  0  1  0
 1  1  1  0  0  0  0  0  1  0  0  1
 0  0  1  0  0  0  0  1  1  1  1  1
 0  0  1  0  0  0  1  0  0  1  0  0
 0  1  1  1  1  0  1  1  1  0  0  1

Or we could prefer to work with predefinite patterns:

julia> function putglider(i, j)
           space[i-2:i, j-2:j]=[0 1 0; 0 0 1; 1 1 1]
       end

Defining some particular i and j, which by the way can be done in one shot:

julia> i, j = 2, 5
(2, 5)

we can count how many live cells surround this point with:

sum(vcat(space[i-1, [j-1,j,j+1]], space[i, [j-1,j+1]], space[i+1, [j-1,j,j+1]]))

where the only point needing explanation here is vcat which concatenate the three arrays. The programming of the updating rule is simple enough and follow previous' lecture implementation:

julia> function nextgen()
               nextspace=copy(space);
                     for i=2:size(space)[1]-1
                         for j=2:size(space)[2]-1
                             nsurround=sum(vcat(space[i-1, [j-1,j,j+1]], space[i, [j-1,j+1]], space[i+1, [j-1,j,j+1]]))
                             if space[i,j]==1
                                 if nsurround<2 || nsurround>3
                                     nextspace[i,j]=0
                                 end
                             else
                                 if nsurround==3
                                     nextspace[i,j]=1
                                 end
                             end
                         end
                     end
                     global space=copy(nextspace);
               end

we can now have fun:

julia> cleanspace(12)
12×12 Array{Int64,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

julia> putglider(2,2)
3×3 Array{Int64,2}:
 0  1  0
 0  0  1
 1  1  1

julia> nextgen()
12×12 Array{Int64,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  1  0  1  0  0  0  0  0  0  0  0
 0  0  1  1  0  0  0  0  0  0  0  0
 0  0  1  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

julia> nextgen()
12×12 Array{Int64,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  1  0  0  0  0  0  0  0  0
 0  1  0  1  0  0  0  0  0  0  0  0
 0  0  1  1  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

julia> nextgen()
12×12 Array{Int64,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  1  0  0  0  0  0  0  0  0  0
 0  0  0  1  1  0  0  0  0  0  0  0
 0  0  1  1  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

julia> nextgen()
12×12 Array{Int64,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  1  0  0  0  0  0  0  0  0
 0  0  0  0  1  0  0  0  0  0  0  0
 0  0  1  1  1  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

Did you see what happened there? We can't blame you if not, it's not very easy to visualize streams of 0 & 1. Next Lecture, we'll start to look at another thing that computers are very good at doing, namely, plotting graphics! Then the bit of magic you may have already perceived will be revealed in its entirety!

To conclude today's lecture, one last useful comment on the REPL: the history of commands are stored in a file (typically ~/.julia/logs/repl_history.jl) so if you are looking for something you did a long time ago, you can grep through this. Do you remember where we started today?!

shell> head ~/.julia/logs/repl_history.jl
# time: 2021-02-12 11:21:03 CET
# mode: julia
	2+2
# time: 2021-02-12 11:21:09 CET
# mode: julia
	x=3
# time: 2021-02-12 11:22:27 CET
# mode: julia
	x=4
# time: 2021-02-12 11:22:29 CET

Links

• Learn Julia in x minutes.