Julia as a CLI Calculator

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.4.2 (2020-05-23)
 _/ |\__'_|_|_|\__'_|  |
|__/                   |

julia>|

I have been using the Julia REPL as my day-to-day calculator on the command line for more than a year now. I am surprised that I haven’t seen more posts talking about how powerful and extensible the language is for simple numerical calculations and scripting In fact, the Julia REPL is much more than a simple calculator, as we shall see later .

I have tried a dozen of alternatives to Julia before, such as MATLAB /Octave, Python, R, Mathematica, and, very briefly, the more old-school bc and dc calculators, but nothing feels more natural to me than Julia.

Of course, Julia can do simple numerical calculations on numbers, just like any other languages. For example, floating point division To open up the Julia REPL , run julia from the command line ,

In [1]

1.75

variables, imaginary numbers,

x = ans + 2im

1.75 + 2.0im

and mathematical functions ( ).

exp(x)

-2.3947596992055047 + 5.232645405696192im

However, who works with plain numbers these days?

Oftentimes, I need to work with lists of objects In Julia’s lingo, it is Array s, but I would prefer to call it lists . For example, vectors, or lists of numbers, and matrices, which can either be viewed as a list of vectors, or a 2 D list of numbers Here, we will use to represent column vectors and to represent row vectors , and sometimes even lists of matrices, lists of higher dimensions, and lists of other objects such as strings or dates.

Julia is a Domain Specific Language ( DSL ) designed for performing numerical calculations on lists (Of course, it’s much more than this). Many people might have been told that Julia is a performant language If you are looking for performace, Julia might disappoints you as a calculator due to its JIT compilation lag , but it is really the syntax of Julia that shines.

Here, I will briefly talk about how I use Julia as calculator on lists. Don’t treat it as a tutorial on how to program , or how to write efficient programs, in Julia. However, hopefully I could show you how to use Julia with syntactically meaningful code. Such skill can be transferred into other languages such as Python or MATLAB , though it might take some work.

Contents

The source code is available as a Jupyter notebook on GitHub , though Jupyter’s capability is too limiting for the code examples to work. I would still suggest you work through the examples in a Julia REPL session run from the command line instead.

If you found any mistakes or have any suggestions, comments, etc., pleaselet me know.

1 A taste of Julia

I originally came to Julia because its MATLAB -like syntax. For example, you could define a matrix just with the exactly the same syntax as MATLAB It’s not exactly the same, but very close ,

A = [1 2 3; 4 5 6]

2×3 Array{Int64,2}:
 1  2  3
 4  5  6

or alternatively,

A = [1 2 3
     4 5 6]

2×3 Array{Int64,2}:
 1  2  3
 4  5  6

and you could do matrix multiplication like this,

A * [3; 2; 1]

2-element Array{Int64,1}:
 10
 28

which calculates the matrix product There would be much more noises in the syntax if you try to do the same thing in numpy .

import numpy as np
A = np.array([[1,2,3],
              [4,5,6]])
A @ np.array([[1,2,3]]).T

array([[10],
       [28]])

If you already know about MATLAB from either academic or work experiences, you shouldn’t have any trouble picking up Julia There are several great cheet­sheets online for transferring from other languages . Many of the MATLAB syntax should just work out of the box. For example, you can write a for loop with the traditional MATLAB syntax.

for i = 1:6
  println(A[i])
end

Out [8]

However, Julia is not a MATLAB clone nor compatible with MATLAB , you can find traces of almost every language mentioned above, plus some more, especially Lisp If you want to desugar the syntax of Julia, knowing Lisp is very important . It might help to think of Julia as an open source, lightweight, extensible MATLAB that looks more like an actual programming language, but keep in mind that it is capable of multiple programming paradigms.

If you are familiar with Python, you can also write a for loop like this Note that Julia is 1 -indexed, unlike Python is 0 -indexed .

for i in 1:3
  println(i)
end

Out [9]

or, alternatively, in Julia’s own syntax.

for i ∈ 1:0.5:2
  println(i)
end

1.0
1.5
2.0

This is one of the highlights of Julia’s syntax. The language is deeply integrated with Unicode characters Because of this, I have to select a different monospace typeface just for this page . In the REPL , you can type the command \in , then press the TAB key to insert the ∈ character.

for i \in|<tab>    ⇒    for i ∈|

If you encounter any character that you don’t know how to spell, for example ℝ , press ? in the REPL . This will lead you to the help mode. Paste the character in, and hit ENTER .

help?> ℝ

"ℝ" can be typed by \bbR<tab>

You can also find get the documentation for a particular function or variable using the same method.

help?> apropos

search: apropos hasproperty

  apropos(string)

  Search through all documentation for a string, ignoring case.

Most of the time, I would prefer to write for loops in the Julia or Python syntax, since the for loop in Julia is more similar to Python— it is based on generators. This is something not possible in MATLAB .

for a ∈ A
  println(a)
end

Out [13]

It is also possible to use the for loop for list comprehensions in Julia Note that it returns a matrix .

[a^2 for a ∈ A]

2×3 Array{Int64,2}:
  1   4   9
 16  25  36

Up until now, we have always been treating a matrix as a 2 D list of numbers. Most of the time, it is more meaningful to treat a matrix as a list of column vectors or a list of row vectors Since Julia has column major arrays, it’s more common to work with column vectors .

2 Rotation: an extended example

Let us look at a practical example that I have encountered several times before. Suppose we want to rotate a triangle around the origin by radians, counterclockwise, in TikZ .

This picture is produced by

\begin{tikzpicture}[>=stealth]
  % axis
  \draw[->] (-2,0) -- (2,0);
  \draw[->] (0,-2) -- (0,2);

  % vertices of the triangle
  \coordinate (v1) at (0.5,0.5);
  \coordinate (v2) at (1.5,0.5);
  \coordinate (v3) at (0.5,1.5);

  % connect edges
  \draw[very thick] (v1) -- (v2) -- (v3) -- cycle;
\end{tikzpicture}

Of course, there are ways to do this directly in TikZ, but rather than searching through the 1300 + pages manual, now that we already know the coordinates of the vertices, it is easy to compute the rotated coordinates using some simple calculations. Let us first input the coordinates into Julia This is easy if you know how to use the visual block mode from Vim .

v1 = [0.5; 0.5]
v2 = [1.5; 0.5]
v3 = [0.5; 1.5]

See? Now we have a list of vectors , or a matrix There is a slight distinction between a matrix and a list of vectors, which will be addressed later .

V = [v1 v2 v3]

2×3 Array{Float64,2}:
 0.5  1.5  0.5
 0.5  0.5  1.5

Our goal is to apply the same operation rotate on all the vectors, i.e.

rotate(v1)
rotate(v2)
rotate(v3)

Before we continue, let us first define the rotate function. To rotate a 2 D vector counterclockwise around the origin by radians, we only need to multiply it with the rotation matrix on the left. If then the rotated vector is Note that the rotation matrix takes an argument as input, i.e. we need to define a function that takes an angle as input and output a matrix.

Julia has a very clean syntax to define such a function.

# ℝ → ℝ^(2×2)
R(θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)]

R (generic function with 1 method)

Since we want to rotate by radians, our rotation matrix is

R(π/3)

2×2 Array{Float64,2}:
 0.5       -0.866025
 0.866025   0.5

Therefore, we can define the rotate function like this Do you notice how easy it is to transfer mathematical ideas to Julia? .

# ℝ^2 → ℝ^2
rotate(v) = R(π/3) * v

rotate (generic function with 1 method)

The rotate function takes an vector and returns another vector, i.e. its type is Julia does have a type system , but it is dynamic Because of this, the rotate function above would happily accept a scalar value as well and, because the (optional) type signature is declared inline, it can get quite messy if you want to use Julia as a simple calculator. Hence, I’ll use implicit type signature as comments.

Here is how we can use this function

rotate(v1)

2-element Array{Float64,1}:
 -0.18301270189221924
  0.6830127018922194

To apply the function rotate on V , there are quite a few approaches, and each of them will have a different flavor.

2.1 The functional approach

If you are familiar with functional programming, you might immediately recognize that rotate is a pure function and might want to do something like this.

map(rotate, V)

Yes, Julia does have a map function, but it will treat the matrix as a 2 D list of numbers. For example, if we have a function dup , defined using lambda notation,

# ℝ → ℝ^2 (?)
dup = x -> [x;x]

#3 (generic function with 1 method)

that duplicates a number into a 2 D vector. Applying dup on V using map will give us a 2 D list of vectors,

map(dup, V)

2×3 Array{Array{Float64,1},2}:
 [0.5, 0.5]  [1.5, 1.5]  [0.5, 0.5]
 [0.5, 0.5]  [0.5, 0.5]  [1.5, 1.5]

and we can try something even crazier.

map(dup, map(dup, V))

2×3 Array{Array{Float64,1},2}:
 [0.5, 0.5, 0.5, 0.5]  [1.5, 1.5, 1.5, 1.5]  [0.5, 0.5, 0.5, 0.5]
 [0.5, 0.5, 0.5, 0.5]  [0.5, 0.5, 0.5, 0.5]  [1.5, 1.5, 1.5, 1.5]

Well, this doesn’t look like our intended type signature at all.

This is because the syntax [v1; v2] will concatenate vectors vertically if v1 and v2 are not scalars. dup can be applied to any vector and it will return a vector, it can even be applied to matrices.

dup([1 2; 3 4])

4×2 Array{Int64,2}:
 1  2
 3  4
 1  2
 3  4

The compiler will dispatch the underlying function ( vcat ) of the syntax [v1;v2] to its corresponding implementations based on the type of the argument A more important example of this is that special matrices, e.g. symmetric ones, are tagged with types to dispatch them into the most efficient implementation . The @which macro can be used to see which implementation is used.

@which [1;1]

vcat(X::T...) where T<:Number in Base at abstractarray.jl:1277

@which [v1;v1]

vcat(arrays::Array{T,1}...) where T in Base at array.jl:1594

The multiple dispatch of the underlying function hvcat of [a b; c d] allows you to write block matrices easily We will take a look at this in detail .

[V v1; v2 V]

4×4 Array{Float64,2}:
 0.5  1.5  0.5  0.5
 0.5  0.5  1.5  0.5
 1.5  0.5  1.5  0.5
 0.5  0.5  0.5  1.5

Let us go back to our main topic. To apply a function to each column vector of a matrix, there are also multiple ways. Let us first look at a simpler one.

2.1.1 Approach 1 : eachcol

To map over each column of a matrix, Julia provides a function conveniently named eachcol .

cols = eachcol(V)

Base.Generator{Base.OneTo{Int64},Base.var"#179#180"{Array{Float64,2}}}
(Base.var"#179#180"{Array{Float64,2}}([0.5 1.5 0.5; 0.5 0.5 1.5]),
Base.OneTo(3))

The exact meaning of the output is not important, we just need to know that it is a generator . It is mostly the same as the generators in Python. To convert a generator into a list, we could use the collect function

collect(cols)

3-element Array{SubArray{Float64,1,Array{Float64,2},
Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1}:
 [0.5, 0.5]
 [1.5, 0.5]
 [0.5, 1.5]

Again, the type signature is not important here, we only need to know that it converts a matrix to a list of column vectors.

In fact, we don’t need to convert the generator into a list, map can map over a generator directly, so we can write apply our rotate like this As a side note on the notation, \prime ( ′ ) is not the same as ' .

V′ = map(rotate, eachcol(V))

3-element Array{Array{Float64,1},1}:
 [-0.18301270189221924, 0.6830127018922194]
 [0.3169872981077809, 1.549038105676658]
 [-1.049038105676658, 1.1830127018922196]

It does exactly what it claims— map the function rotate to each column of V. This is exactly what we needed. In fact, there is a simpler way to write it using the dot syntax .

V′ = rotate.(eachcol(V))

3-element Array{Array{Float64,1},1}:
 [-0.18301270189221924, 0.6830127018922194]
 [0.3169872981077809, 1.549038105676658]
 [-1.049038105676658, 1.1830127018922196]

For generators and 1 D lists of objects with functions that take only 1 argument, the dot syntax means exactly the same as map But for higher dimensions, it is broadcast under the hood . MATLAB has a similar syntax, but it only applies to a few operators. In Julia, this syntax applies to any functions.

exp.([1,2,3,4])

4-element Array{Float64,1}:
  2.718281828459045
  7.38905609893065
 20.085536923187668
 54.598150033144236

The dot syntax also applies to binary operators and functions with more arguments. For infix operators, the dot is in front of the operator.

[1,2,3,4] .* [4,3,2,1]

4-element Array{Int64,1}:
 4
 6
 6
 4

It is equivalent to

.*([1,2,3,4], [4,3,2,1])

4-element Array{Int64,1}:
 4
 6
 6
 4

In order to use the rotated vectors V′ in our TikZ plot, we have to do some clean up.

For each vector v , we will first round each element to 5 digits, then concatenate with comma.

# Vector → String
cleanup(v) = join(round.(v, digits=5), ", ")

cleanup (generic function with 1 method)

Here, the dot syntax means exactly the same as

map(x -> round(x, digits=5), v)

For function compositions like this, Julia also borrows the forward pipe operator |> popularized by F# R also has a similar operator via the magrittr package .

1 |> exp |> sin

0.41078129050290885

This means the same as

sin(exp(1))

0.41078129050290885

In particular, x |> f is equivalent to f(x) . Because of the dot syntax, we also have the .|> operator that allows us to apply functions on lists.

[1,2,3] .|> exp .|> sin

3-element Array{Float64,1}:
 0.41078129050290885
 0.8938549549128102
 0.9444710089262849

Or, with a little bit of trick It’s broadcasting under the hood , apply multiple functions on each element of the list.

[1,2,3] .|> [exp sin]

3×2 Array{Float64,2}:
  2.71828  0.841471
  7.38906  0.909297
 20.0855   0.14112

This is equivalent to

[exp(1) sin(1); exp(2) sin(2); exp(3) sin(3)]

3×2 Array{Float64,2}:
  2.71828  0.841471
  7.38906  0.909297
 20.0855   0.14112

Don’t worry if you don’t understand the example above right now, there is a distinction, but also similarity, between a matrix and list of vectors that we will address in the next section.

Just for fun, let us implement a Caesar cipher , which simply shifts all the characters in a string by a fixed amount, using the pipe operator You can use ESC + ENTER key to insert a line break in the REPL .

# (String, ℤ) → String
encrypt(str, key) = (
                      str                     |>
                      collect                .|>
                     (x -> x - 'a')          .|>
                     (x -> mod(x + key, 26)) .|>
                     (x -> x + 'a')
                    )                         |>
                      join

Note that the .|> operator has a higher precedence than |> , so we have to use a parenthesis here. Also, since both |> and .|> have higher precedence than -> ,

(1, 2) .|> x -> x^2 |> sum

(1, 4)

(1, 2) .|> (x -> x^2) |> sum

Out [43]

we need to group anonymous functions with parentheses Thanks u/multiple­attackers for pointing this out , though for our encrypt function, this doesn’t matter.

We can encrypt a string by

encrypt("hellojulia", 20)

"byffidofcu"

Try come up with the decryption function using pipe operators.

Using the forward pipe operator, we don’t even need an extra cleanup function to map over V′ . We can directly apply the clean up using pipes.

V′                         .|>
(x -> round.(x, digits=5))  |>
(x -> join(x, ", "))        |>
println

-0.18301, 0.68301
0.31699, 1.54904
-1.04904, 1.18301
3-element Array{Nothing,1}:
 nothing
 nothing
 nothing

Don’t worry about the nothing in the output, it is the return value of println , similar to IO () in Haskell. We could suppress it by adding a semicolon ; at the end.

V′                         .|>
(x -> round.(x, digits=5))  |>
(x -> join(x, ", "))        |>
println;

-0.18301, 0.68301
0.31699, 1.54904
-1.04904, 1.18301

One of the complaints I have about Julia, though, is that it doesn’t support partial application and currying by default. For example, in Haskell, you could write something like

map (3*) [1,2,3]

[3,6,9]

but it is not possible in Julia. There is a function in the standard library called Base.Fix2 , which you could rename it curry to mimic the syntax,

curry = Base.Fix2
map(curry(*,3), [1,2,3])

3-element Array{Int64,1}:
 3
 6
 9

but it is still not as convenient. There is an ongoing proposal to solve this issue, but it’s not merged yet. Once it’s been merged, we could write something much cleaner.

V′                  .|>
round.(_, digits=5)  |>
join(_, ", ")        |>
println

However, now that we have already defined the cleanup function, which takes only one argument, there is a much slicker way to do what we have done so far. The trick is that the pipe operator |> happily takes generators as well.

eachcol(V) .|> rotate .|> cleanup .|> println;

-0.18301, 0.68301
0.31699, 1.54904
-1.04904, 1.18301

Now, we can finally use the coordinates to draw the rotated triangle,

where the picture is produced using

\begin{tikzpicture}[>=stealth]
  % axis
  \draw[->] (-2,0) -- (2,0);
  \draw[->] (0,-2) -- (0,2);

  % vertices of the triangle
  \coordinate (v1) at (-0.18301, 0.68301);
  \coordinate (v2) at ( 0.31699, 1.54904);
  \coordinate (v3) at (-1.04904, 1.18301);

  % connect edges
  \draw[very thick] (v1) -- (v2) -- (v3) -- cycle;
\end{tikzpicture}

2.1.2 Desugaring the array syntax

As you might have noticed, there is a slight distinction between a matrix, which has type Array{T,2} = Matrix{T} ,

A = [1 2 3; 4 5 6]

2×3 Array{Int64,2}:
 1  2  3
 4  5  6

Array{Integer,2} == Matrix{Integer}

true

and a list of vectors, which has type Array{Array{T,1},1} = Vector{Vector{T}} ,

B = [[1, 4], [2, 5], [3, 6]]

3-element Array{Array{Int64,1},1}:
 [1, 4]
 [2, 5]
 [3, 6]

Array{Array{Integer,1},1} == Vector{Vector{Integer}}

true

They are, of course, not equivalent.

A == B

false

Matrix{Integer} == Vector{Vector{Integer}}

false

To understand this, it would be easier to look at the distinction between the following three syntax.

# List (Array)
[1,2,3]

3-element Array{Int64,1}:
 1
 2
 3

# Column vector
[1;2;3]

3-element Array{Int64,1}:
 1
 2
 3

# Row vector
[1 2 3]

1×3 Array{Int64,2}:
 1  2  3

Looking at the output, we should expect that [1,2,3] to be equivalent to [1;2;3] .

[1,2,3] == [1;2;3]

true

The range syntax also produces a list, and thus a column vector.

1:3 == [1;2;3]

true

This is why Julia is not compatible with the matrix syntax of MATLAB . This is valid in MATLAB ,

[1,2,3; 4,5,6]

ans =
     1     2     3
     4     5     6

but not valid in Julia.

[1,2,3; 4,5,6]

ERROR: syntax: unexpected semicolon in array expression

It is indeed true, though, that in Julia a column vector (of numbers) is equivalent to a list of numbers , so we could apply matrix multiplication to a list of numbers as well.

[1 2; 3 4] * [1,1]

2-element Array{Int64,1}:
 3
 7

However, there is a significant difference between the list syntax and the column vector syntax for a list of vectors .

[[1,2], [3,4], [5,6]]

3-element Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]
 [5, 6]

In mathematical notation, this object is However, if we define the matrix using the ; syntax,

[[1,2]; [3,4]; [5,6]]

6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6

the result is a flat matrix As we can see, it’s much appropriate to call the ; symbol vertical concatenation ( vcat ), instead of creating a column. Similarly, the space character is called horizontal concatenation ( hcat ).

[[1 2] [3 4] [5 6]]

1×6 Array{Int64,2}:
 1  2  3  4  5  6

Compare with the result using the list syntax,

[[1 2], [3 4], [5 6]]

3-element Array{Array{Int64,2},1}:
 [1 2]
 [3 4]
 [5 6]

This object should be visualized as If we hcat column vectors, we would get an matrix.

[[1;2] [3;4] [5;6]]

2×3 Array{Int64,2}:
 1  3  5
 2  4  6

Let us revise our comments.

# List (Array)
[a,b,c]
# Vertical concatenation
[a;b;c]
# Horizontal concatenation
[a b c]
# A list of lists
[[a],[b],[c]]
# Vertical concatenation of lists
[[a];[b];[c]]
# Horizontal concatenation of lists
[[a] [b] [c]]
# Vertical concatenation of horizontal concatenations of lists
[[a] [b]; [c] [d]]

Remember this. Now let us have some fun.

# A list of a single number
a = [1]

1-element Array{Int64,1}:
 1

# More lists of a single number
b = [2]
c = [3]
d = [4]

# Vertical concatenation of horizontal concatenations of lists
# of lists of a single number
[[a] [b]; [c] [d]]

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

# Vertical concatenation of lists of a single number
[a; b]

2-element Array{Int64,1}:
 1
 2

# Matrix multiplication of
# (vertical concatenation of horizontal concatenations of lists
# of lists of a single number)
# and
# (vertical concatenation of lists of a single number)
[[a] [b]; [c] [d]] * [a; b]

2-element Array{Array{Int64,1},1}:
 [5]
 [11]

Visually, Every step can be justified mathematically, though it involves two vector spaces, one is , the vector space of all matrices, and the other one is .

Try it in MATLAB and you would be disappointed.

Intermezzo: The Julia AST

In fact, there is a much easier way to understand this. We can take a look at the abstract syntax tree ( AST ) of each notation.

Meta.show_sexpr(:([1,2,3]))

(:vect, 1, 2, 3)

Meta.show_sexpr(:([1;2;3]))

(:vcat, 1, 2, 3)

Meta.show_sexpr(:([1 2 3]))

(:hcat, 1, 2, 3)

Here, the function Meta.show_sexpr can show any expression as a Lisp-style S-expression . It is very useful if you want to desugar the syntax of Julia. As we can see, [a;b] is parsed as :vcat and [a b] is parsed as :hcat . They really just mean concatenations.

The notation :vcat denotes a symbol in Julia, which is similar to atoms/symbols in Lisp and the syntax : is similar to quotes in Lisps, e.g.

'vcat

'vcat

'1

Out [78]

and in Julia

:vcat

:vcat

In [80]

Out [80]

The notation :(1 + 1) denotes an expression , where the :() syntax is similar to quasiquotes in Lisps,

`(+ 1 2)

'(+ 1 2)

`(+ ,(+ 1 3) 2)

'(+ 4 2)

and in Julia The splicing syntax works on strings as well

:(1 + 2)

:(1 + 2)

:($(1 + 3) + 2)

:(4 + 2)

We can evaluate an expression using the eval function.

eval(:($(1 + 3) + 2))

Out [85]

One of the reasons Julia has such first class support for these is because it supports metaprogramming , so we could transform the AST directly The other one is that part of the Julia compiler is written in a custom Lisp dialect called femtolisp . For example, there is a @. macro for converting every function call into the dot syntax.

@. [1,2,3] |> exp |> sin

3-element Array{Float64,1}:
 0.41078129050290885
 0.8938549549128102
 0.9444710089262849

If we take a look at the syntax for defining matrices,

Meta.show_sexpr(:([1 2 3; 4 5 6]))

(:vcat, (:row, 1, 2, 3), (:row, 4, 5, 6))

we can notice that it’s interpreted as a vertical concatenation of two row s. We already know we can use the @which macro to see which function is responsible for creating the array, but we can go even further and use the @edit macro to jump to the source code of the implementation

@edit [1 2 3; 4 5 6]

function hvcat(rows::Tuple{Vararg{Int}}, xs::T...) where T<:Number
  nr = length(rows)
  nc = rows[1]

  a = Matrix{T}(undef, nr, nc)
  if length(a) != length(xs)
    throw(ArgumentError("argument count does not match specified shape (expected $(length(a)), got $(length(xs)))"))
  end
  k = 1
  @inbounds for i=1:nr
    if nc != rows[i]
      throw(ArgumentError("row $(i) has mismatched number of columns (expected $nc, got $(rows[i]))"))
    end
    for j=1:nc
      a[i,j] = xs[k]
      k += 1
    end
  end
  a
end

As we can see, even though the parser interprets matrices as vertical concatenation of rows, this particular instance of the function hvcat is responsible for actually creating the matrix.

The counterpart, a list of lists, is parsed differently,

Meta.show_sexpr(:([[1,4], [2,5], [3,6]]))

(:vect, (:vect, 1, 4), (:vect, 2, 5), (:vect, 3, 6))

and the function being executed is vect .

@edit [[1,4], [2,5], [3,6]]

vect(X::T...) where {T} = T[ X[i] for i = 1:length(X) ]

If you encounter any weird Julia syntax in the future, checking the AST and the source code might be a good bet. For those of you who are getting hungry, how about some snacks?

$ julia --lisp

2.1.3 Approach 2 : eachslice and mapslices

In the previous section, we have been talking about mapping a function over a matrix, or 2 D list, which has dimension 2 .

ndims(V)

Out [91]

For mapping a function to lists of higher dimensions, there is also a function called eachslice . The following is equivalent to eachcol .

eachslice(V, dims=2) |> collect

3-element Array{SubArray{Float64,1,Array{Float64,2},
Tuple{Base.Slice{Base.OneTo{Int64}},Int64},true},1}:
 [0.5, 0.5]
 [1.5, 0.5]
 [0.5, 1.5]

However, to understand it, we have to talk about slices first. It is mostly the same as the slices in MATLAB and Python.

Back when we defined V , the REPL tells us that the dimension of V is 2 × 3 .

V

2×3 Array{Float64,2}:
 0.5  1.5  0.5
 0.5  0.5  1.5

You can query this directly.

row, col = size(V)

(2, 3)

Each element in the matrix is indexed by Note that even though Julia is column major, the first index is still the row number We could select a single element from the matrix using square bracket,

V[1,2]

1.5

and we could select multiple elements by passing a list as index You can also use 1:2 or 1:end here .

V[[1,2], 2]

2-element Array{Float64,1}:
 1.5
 0.5

This selects the element corresponding to and , returned as a list, or column vector. The returned vector will have the same shape as the index, so you could also index with a row vector.

V[[1 2], 2]

1×2 Array{Float64,2}:
 1.5  0.5

or even a matrix.

V[[1 2; 2 1], 2]

2×2 Array{Float64,2}:
 1.5  0.5
 0.5  1.5

You could understand the result of this particular case using map , but do notice that this interpretation is not true if the index in the second slot is not a scalar.

map(i -> V[i,2], [1 2; 2 1])

2×2 Array{Float64,2}:
 1.5  0.5
 0.5  1.5

The symbol : can be used to to select all the indices in the given dimension, and return the result with the same shape as a list, or equivalently, a column vector , as mentioned before.

V[:,2]

2-element Array{Float64,1}:
 1.5
 0.5

This is exactly our . Hence, we could write the matrix V as

[V[:,1] V[:,2] V[:,3]] == V

true

The eachslice function essentially replaces every slot of the index with : , except the slot denoted by dims . The case above is equivalent to Not exactly, the underlying implementation uses view under the hood to avoid copying and allow modification

map(i -> V[:,i], 1:3)

3-element Array{Array{Float64,1},1}:
 [0.5, 0.5]
 [1.5, 0.5]
 [0.5, 1.5]

but do notice that eachslice returns a generator, not a list.

There is also a more flexible function called mapslices , which can be used to map a function over a matrix as it is, without converting it to a list of vector.

mapslices(rotate, V, dims=1)

2×3 Array{Float64,2}:
 -0.183013  0.316987  -1.04904
  0.683013  1.54904    1.18301

where the dims=1 here means the arguments to rotate are V[:,j] for each j∈1:3 , where the colon : is in the first slot.

Similarly, if we set dims=2 , the arguments would be V[i,:] for i∈1:2 , where : is in the second slot A caveat, though, is that remember V[i,:] is also a column vector because it has the same shape as a list .

mapslices(sum, V, dims=2)

2×1 Array{Float64,2}:
 2.5
 2.5

This is different from eachslice . For eachslice , the dims=2 means the colon : is in all the slots, except the second slot.

T = rand(3,3,3)
map(i -> T[:,i,:], 1:3) == eachslice(T, dims=2) |> collect

true

If this sounds confusing, there are also other approaches to try.

2.2 The imperative approach

Unlike Python and MATLAB , because of JIT compilation, it is perfectly fine, and sometimes preferred , to use loops to iterate over large matrices.

Recall that eachcol returns an generator, so we could use it in for loops as well.

for v ∈ eachcol(V)
  v |> rotate |> cleanup |> println
end

-0.18301, 0.68301
0.31699, 1.54904
-1.04904, 1.18301

Notice how clean and descriptive is the syntax, though not as clean as our functional approach :) . You can also perform the same operation using list comprehensions.

[v |> rotate |> cleanup for v ∈ eachcol(V)]

3-element Array{String,1}:
 "-0.18301, 0.68301"
 "0.31699, 1.54904"
 "-1.04904, 1.18301"

I would usually prefer to use iterators and generators over the raw way,

for i ∈ 1:size(V,2)
  V[:,i] |> rotate |> cleanup |> println
end

-0.18301, 0.68301
0.31699, 1.54904
-1.04904, 1.18301

but the raw way is sometimes more flexible on how you could index the list.

As a side note, even though this is irrelevant if you want to use Julia as a calculator, do notice that, unlike numpy, but similar to MATLAB , the slicing operation will create a copy of the slice in the memory. We can use the @allocated to see the number of bytes allocated for an expression.

# a large matrix
A = rand(1000, 1000);
# trigger JIT compilation, as will be explained later
A[:,1];

@allocated A[:,1]

Out [109]

To avoid such copying, we need to use the view function, which creates a mutable reference to the corresponding elements in A .

# trigger JIT compilation, as will be explained later
view(A,:,1);

@allocated view(A,:,1)

Out [110]

or equivalently there is a macro called @view for converting a slicing operation to view.

# trigger JIT compilation, as will be explained later
@view A[:,1];

@allocated @view A[:,1]

Out [111]

Therefore, the more memory-efficient way to write the raw iteration above should be

for i ∈ 1:size(V,2)
  view(V,:,i) |> rotate |> cleanup |> println
end

-0.18301, 0.68301
0.31699, 1.54904
-1.04904, 1.18301

There is also a @views macro for converting every slicing operation in a block with views, so we could also equivalently write.

@views for i ∈ 1:size(V,2)
         V[:,i] |> rotate |> cleanup |> println
       end

-0.18301, 0.68301
0.31699, 1.54904
-1.04904, 1.18301

2.3 The linear algebra approach

If you are a mathematician, thank you for your patience.

I deliberately left this approach at the end, because there would be no place for other approaches if I talk about it first. Applying a (linear) operation on a list of vectors is so fundamental to linear algebra that it is how matrix multiplication is defined and similarly, for linear operations on row vectors (dual vectors), The problem of applying the rotation matrix to each column of is essentially a simple matrix product. Therefore, to calculate the rotated vector, we simply use the * symbol to perform a matrix product.

R(π/3) * V

2×3 Array{Float64,2}:
 -0.183013  0.316987  -1.04904
  0.683013  1.54904    1.18301

and we can possibly take a transpose (adjoint) for easier copy-pasting For real transpose, pipe ans to the transpose function, but it doesn’t matter here .

ans'

3×2 Adjoint{Float64,Array{Float64,2}}:
 -0.183013  0.683013
  0.316987  1.54904
 -1.04904   1.18301

For more linear algebra related functions, you can import the LinearAlgebra package. For example, if we want to check that our rotation matrix is indeed a rotation, i.e. has determinant 1 , we could use the det function

using LinearAlgebra
det(R(π/3))

1.0

or you could show off the pipe operator by writing

π/3 |> R |> det

1.0

However, the matrix multiplication trick only works for linear transformations, though luckily, even though the world is often filled non-linear problems, usually they can be well approximated via linear transformations, or with linear combination of simple nonlinear transformations, at least for those we could comprehend.

To deal with simple non-linear transformations, Julia has broadcasting syntax to apply a non-linear (or linear) transformation to a list of vectors.

Let say we are no longer satisfied with rotating the vertices around the origin. We want to use , which is the bottom-left corner of the triangle as the pivot for rotation.

This transformation can be done by first translating all the vertices by , shifting to origin.

perform the rotation

and then shift everything back.

There is a problem if we try to come up with a matrix for the transformation. Because is the only point in the entire plane kept fixed, the origin is not mapped to the origin under this transformation, so there won’t be a matrix.

However, as we have already decomposed the transformation into three smaller pieces. We can try to do it step by step.

First, we need to subtract from all the vectors in . This operation is not linear, but we can use the dot syntax mentioned before to broadcast subtraction by to all the column vectors of .

V .- V[:,1]

2×3 Array{Float64,2}:
 0.0  1.0  0.0
 0.0  0.0  1.0

It will basically extend to the same dimension as the matrix, then perform subtraction. However, since Julia’s broadcasting also applies to ordinary functions, it is sometimes better to think with the following mental model We could also perform operation to each row by broadcasting with a row matrix.

V .+ V[1,:]'

2×3 Array{Float64,2}:
 1.0  3.0  1.0
 1.0  2.0  2.0

Recall that V[1,:] has the same shape as a list, or a column vector, so we have to apply a transpose (adjoint) to get the corresponding row vector Again, for complex matrices, pipe it into the transpose function

V[1,:]'

1×3 Adjoint{Float64,Array{Float64,1}}:
 0.5  1.5  0.5

The mental model in this case should be Or, we could apply operation to each individual entry of the matrix by supplying a number

V .+ V[1,1]

2×3 Array{Float64,2}:
 1.0  2.0  1.0
 1.0  1.0  2.0

and the mental picture should be Now, to continue our main topic, we now perform the remaining two steps to V.

R(π/3) * (V .- V[:,1]) .+ V[:,1]

2×3 Array{Float64,2}:
 0.5  1.0      -0.366025
 0.5  1.36603   1.0

Now, let us consider another scenario, what if we want to rotate around the center, i.e. the average of the vertices, of the triangle?

This require us to calculate the mean of all the column vectors in V . Of course, for people familiar with linear algebra, such folding pattern can be inherently represented as linear functionals. By treating the column vector as a dual vector (something that eats an regular vector and spit out a number) to V ’s row vectors, we could directly calculate the mean via matrix product.

c = V * [1;1;1]/3

2-element Array{Float64,1}:
 0.8333333333333334
 0.8333333333333334

2.4 Back to functional constructs

However, this approach will require us to construct an additional vector for summing across columns, and more importantly, linear combination, i.e. weighted sum of each element in a list, is only a special case of the more generalized folding pattern. We already know how to do the weighted step via broadcasting, but instead of + , we could apply other binary operations, such as * , min , max to each pair of elements and combine the result.

Julia has a function called reduce and two more called foldl and foldr for explicit associativity , for such general patterns. To compute the mean of each row, we could write

reduce(+, V / 3, dims=2)

2×1 Array{Float64,2}:
 0.8333333333333333
 0.8333333333333333

The argument dims=2 means the reduce operation smashes the second dimension

V[1,:]

3-element Array{Float64,1}:
 0.5
 1.5
 0.5

into 1 number.

It is sometimes useful to think about dims=1 as smashing row vectors (smash downward), and dims=2 as smashing column vectors (smash to the right).

For no particular reason, if we want to calculate the geometric mean of each row we could similarly do

reduce(*, V, dims=2).^(1/size(V,2))

2×1 Array{Float64,2}:
 0.7211247851537042
 0.7211247851537042

To finish our example, let us calculate the rotated vertices,

R(π/3) * (V .- c) .+ c

2×3 Array{Float64,2}:
 0.955342  1.45534  0.0893164
 0.377992  1.24402  0.877992

and hence we have

Intermezzo: Packages

Although Julia provides a nice syntax for handling nonlinear transformations, and we could even define functions to make them composable with linear ones,

translate(t) = v -> v + t
rotate(θ) = v -> R(θ)*v

eachcol(V) .|> translate(-v1) .|> rotate(π/3) .|> translate(v1)

3-element Array{Array{Float64,1},1}:
 [0.5, 0.5]
 [1.0, 1.3660254037844386]
 [-0.3660254037844386, 1.0]

we still don’t consider rotation and translation composable mathematically, since matrix multiplication and vector addition are very different operations.

Mathematicians have come up with other models of describing geometric transformations that make rotation and translation composable, notably the homogeneous model, where, by bringing in one extra dimension, translations can be represented as linear transformations, and the conformal model, where, by bringing in two extra dimensions, translations can be represented as orthogonal transformations (or versor/sandwich product), just like rotations.

There is a nice package for Julia called Grassmann.jl that let us do conformal geometric algebra ( CGA ) in Julia, so we can compose translations and rotations mathematically. I know this is like killing a mosquito with a cannon, but why not?

To install it, press ] in the REPL to enter the Package mode, then run the command add Grassmann .

(@v1.4) pkg> add Grassmann

Press backspace to go back to Julian mode (i.e. the default mode). We can import a package with the using keyword There is also an import keyword if you want to avoid conflicts with existing names .

using Grassmann

It might take a little bit of time to load, since Julia is based on JIT compilation and it needs to compile every expression before it can run it. We will talk about a way to mitigate this issue later.

We won’t go into the details of how geometric algebra or CGA works. Alan Macdonald’s “ A Survey of Geometric Algebra and Geometric Calculus ” and Ben Lynn’s blog posts on geometric algebra in Haskell provide a good overview with some nice animations.

The basic idea is to introduce two extra basis vectors and to the usual orthonormal basis of , giving us a new 4 D space called , sometimes called a Minkowski space in physics, and then generate a geometric algebra with dot product It is conventional to define two more vector, the origin Some sources might have different signs or normalizations in the defition, but they are equivalent and the point at infinity We can import the basis of using the @basis macro.

@basis S"∞∅++" G31 e

(⟨∞∅++⟩, v, v∞, v∅, v₁, v₂, v∞∅, v∞₁, v∞₂, v∅₁, v∅₂,
v₁₂, v∞∅₁, v∞∅₂, v∞₁₂, v∅₁₂, v∞∅₁₂)

Note that we are using a custom prefix e for the basis vectors, so we can reference them using e1 , e2 , e∞ , etc., but they will still be printed as v1 , v2 , v∞ .

e1 + 2*e2

0v∞ + 0v∅ + 1v₁ + 2v₂

To convert our vector to geometric algebra notation , we can use some matrix algebra trick. Recall that matrix-vector product is defined as a linear combination Because is nothing limiting us to use other objects in a matrix, we can construct a row vector [e1 e2] and apply it to each column of V via matrix multiplication. This will convert each column of V to a linear combination of e1 and e2 .

[e1 e2] * V

1×3 Array{Any,2}:
 0.0 + 0.5v₁ + 0.5v₂  0.0 + 1.5v₁ + 0.5v₂  0.0 + 0.5v₁ + 1.5v₂

Let’s take a transpose of the matrix to use it as a column vector The usual adjoint is slightly different here

W = ans |> transpose

3×1 Transpose{Any,Array{Any,2}}:
 0.0 + 0.5v₁ + 0.5v₂
 0.0 + 1.5v₁ + 0.5v₂
 0.0 + 0.5v₁ + 1.5v₂

To apply geometric transformation on vectors, we first need to convert a vector to its conformal representation . The formula for the conversion is less of our focus here. We only need to know that this can be done using the ↑ function. For example, the corresponding conformal vector of is

↑(0.5*e1 + 0.5*e2)

-0.0 + 0.25v∞ + 1.0v∅ + 0.5v₁ + 0.5v₂

Rotations in geometric algebra are defined on planes (bivectors), for , there is only one plane in the entire space, which is , we will denote it by since it squares to under the geometric product.

   = e12

v₁₂

  ^2

-1v

The v here represents a scalar, or simply 1 .

The rotation of a vector by radians counterclockwise can be done using a versor (sandwich) product on the conformal vector If this reminds you of quaternion rotations, quaternions are embedded in where the multiplications are geometric products . In Grassmann.jl’s notation, this can also be written as For example, to rotate the vector by radians, we could write

↓(↑(0.5*e1 + 0.5*e2) ⊘ exp(  *(π/3)/2))

0.0 - 0.1830127018922192v₁ + 0.6830127018922194v₂

which is exactly the same as the result we get previously. The ↓ function converts a conformal vector back to an Euclidean vector.

The translation of a vector by a displacement vector can be performed using versor product as well, Using the construction of CGA , we can now do the rotation around by

# rotation angle
θ = π/3
# displacement vector
t = W[1]

.↑(W) .⊘ exp(-t*e∞/2) .⊘ exp(  *θ/2) .⊘ exp(t*e∞/2) .|> ↓

3×1 Array{MultiVector{⟨∞∅++⟩,Float64,16},2}:
                               0.0 + 0.5v₁ + 0.5v₂
 0.0 + 1.0000000000000002v₁ + 1.3660254037844386v₂
 0.0 - 0.3660254037844386v₁ + 1.0000000000000002v₂

or with the @. macro

@. ↑(W) ⊘ exp(-t*e∞/2) ⊘ exp(  *θ/2) ⊘ exp(t*e∞/2) |> ↓

3×1 Array{MultiVector{⟨∞∅++⟩,Float64,16},2}:
                               0.0 + 0.5v₁ + 0.5v₂
 0.0 + 1.0000000000000002v₁ + 1.3660254037844386v₂
 0.0 - 0.3660254037844386v₁ + 1.0000000000000002v₂

The nice thing about CGA is that it is perfectly composable, we could define a new versor

T = exp(-t*e∞/2) * exp(  *θ/2) * exp(t*e∞/2)

0.8660254037844387 - 0.24999999999999997v∞₁ +
0.24999999999999994v∞₂ + 0.49999999999999994v₁₂

and apply it to each element in W using versor product directly.

@. ↑(W) ⊘ T |> ↓

3×1 Array{MultiVector{⟨∞∅++⟩,Float64,16},2}:
 0.0 + 0.49999999999999983v₁ + 0.5000000000000001v₂
                 0.0 + 1.0v₁ + 1.3660254037844388v₂
                0.0 - 0.36602540378443854v₁ + 1.0v₂

There does exist some issue with precisions, but the answer we get is same. Anyway, at least I won’t use this package for such a simple problem, but I had some fun.

Intermezzo: Plotting

In fact, we don’t actually need to use TikZ for drawing our simple triangle. There are plenty of libraries for making nice looking plots, and even interactive plots, for Julia. Some examples include Plots.jl , Makie.jl and Gadfly.jl .

However, since we are using Julia as a command line calculator, there is a more lightweight library called UnicodePlots.jl to let us do the plotting directly in the terminal. It can be installed the same way in the Package mode.

(@v1.4) pkg> add UnicodePlots

using UnicodePlots

Now, we can draw the same diagram as before directly in the terminal using the lineplot function

# of course, there is no `cycle` command
# so we have to add an extra vertex
U = [V v1]
lineplot( U[1,:], U[2,:]             # x, y
        , xlim=(-2,2), ylim=(-2,2)
        , border=:none, labels=false
        , width=30
        )

⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢸⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢸⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢸⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠸⠤⠤⠤⠤⠤⠤⠵⠄⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⡧⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀

We can apply the rotation in any of the approaches described before, but for simplicity, we will just use the linear algebra approach.

U′ = R(π/3) * U

2×4 Array{Float64,2}:
 -0.183013  0.316987  -1.04904  -0.183013
  0.683013  1.54904    1.18301   0.683013

Then, we can draw the rotated triangle by

lineplot( U′[1,:], U′[2,:]           # x, y
        , xlim=(-2,2), ylim=(-2,2)
        , border=:none, labels=false
        , width=30
        )

⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣇⣀⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⢀⣀⡠⠤⠒⠒⠉⠉⣇⠜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠉⠒⢄⡀⠀⠀⠀⢀⡏⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠤⣠⠃⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⡧⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀

3 Customizations

An advantage of using Julia as a command calculator is that it can be easily customized. Whenever you start a REPL session, a custom script located at Check the documentation if you on Win­dows

~/.julia/config/startup.jl

will be loaded before you can interact with the terminal. You can add helper functions or load common packages, such as LinearAlgebra and Statistics, in the startup.jl file so that you can use them every time you open up Julia.

3.1 My startup.jl

For example, my startup.jl for using Julia as calculator roughly looks like this

ENV["JULIA_EDITOR"] = "nvim"

# packages
using LinearAlgebra
using SymPy
using Latexify

# convert number to hex string
function hex(n :: Number) :: String
  string(n, base=16)
end

# convert number to binary string
function bin(n :: Number) :: String
  string(n, base=2)
end

# convert number to decimal string
function dec(n :: Number) :: String
  string(n, base=10)
end

# --- and some common functions such as rotation matrix, etc.

It is very minimal and does not load many packages, since my Julia calculator session is always running in the background. I would load another file manually if I want to develop packages.

I set my editor to nvim , so that the @edit macro will open up the source code in neovim . I’m also importing the LinearAlgebra package by default since I use linear algebra functions quite often.

3.1.1 SymPy

The SymPy package is for doing symbolic operations. For example, we could define symbolic variables using @vars ,

@vars θ x y real=true

(θ, x, y)

and perform operations on them directly,

w = R(θ) * [x;y]

2-element Array{Sym,1}:
 x*cos(θ) - y*sin(θ)
 x*sin(θ) + y*cos(θ)

using LinearAlgebra
norm(w) |> simplify |> print

sqrt(x^2 + y^2)

3.1.2 Latexify

The package Latexify.jl is for printing the commands for a Julia expression directly.

V |> latexify |> print

\begin{equation}
\left[
\begin{array}{ccc}
0.5 & 1.5 & 0.5 \\
0.5 & 0.5 & 1.5 \\
\end{array}
\right]
\end{equation}

:(sqrt(θ/3 + exp(ϕ^2))) |> latexraw |> print

\sqrt{\frac{\theta}{3} + e^{\phi^{2}}}

3.1.3 Utilities

I have quite a few custom functions defined for Julia, but the most frequently used ones are hex , bin and dec , which can be used to convert between hexadecimal, binary, and decimal numbers.

bin(0xce - 5)

Out [153]

dec(0xf1) |> print

Out [154]

3.2 The JIT overhead

One thing to keep in mind is that, everything you have in the startup.jl file will be executed when you run the julia command from the command line, even when executing a simple hello world script. For example, if we have a file named hello.jl with content

println("hello, " * "world")

and we use the startup.jl file from before, executing this script takes around 8 seconds,

$ time julia hello.jl

hello, world
julia hello.jl  7.85s user 0.17s system 99% cpu 8.052 total

because the SymPy package takes quite a bit of time to load. We will take a look at how you can accelerate the loading time for SymPy later, but to prevent Julia from loading the startup.jl file for simple scripts, we could add a command line argument --startup-file=no

$ time julia --startup-file=no hello.jl

hello, world
julia --startup-file=no hello.jl  0.08s user 0.06s system 88% cpu 0.156 total

It is still fairly slow compared to other languages. This is one of Julia’s flaws, since it only compiles when absolutely necessary, and the overhead of JIT compilation only pays off in the long term. It does run fast, but the compilation is slow. We can use the @time macro in the REPL to display the time it takes to run a command. Let us open up a fresh Julia REPL session and run

@time println("hello, julia")

hello, julia
  0.002414 seconds (24 allocations: 1.734 KiB)

@time println("hello, julia")

hello, julia
  0.000045 seconds (4 allocations: 80 bytes)

@time println("bye, julia")

bye, julia
  0.000042 seconds (4 allocations: 80 bytes)

See? It takes much longer and more memory for Julia to run the first call to the println function, because Julia has to compile code before it can run it. In contrast, the function calls following it are very fast.

However, the pre-compilation only accounts for this particular instance of println . If we pass in another type, the compiler has to compile again.

@time println([1,2,3,4])

[1, 2, 3, 4]
  0.110054 seconds (243.47 k allocations: 11.987 MiB)

@time println([1,2,3,4])

[1, 2, 3, 4]
  0.000106 seconds (34 allocations: 976 bytes)

This is sometimes called the time-to-first-plot issue , because it is very prominent in plotting libraries such as Plots.jl , but you might have noticed that Grassmann.jl also has this issue.

The time-to-first-plot issue is still unsolved, but it is critical if we want to use Julia as a command line calculator, since most of the operations will be run only once. Luckily, there are ways to mitigate this phenomenon, and we will dedicate the rest of this section to address it.

3.3 atreplinit

We can use the atreplinit function to register a function to run before only at a REPL session. It will be ignored when the julia command is not run in the REPL mode, so our hello world example won’t take as long even if we don’t pass in the --startup-file=no option.

3.3.1 My startup.jl , take 2

Using the atreplinit function, my new startup.jl would look like this

function initrepl(repl)
  @eval ENV["JULIA_EDITOR"] = "nvim"

  # packages
  @eval using LinearAlgebra
  @eval using SymPy
  @eval using Latexify

  # convert number to hex string
  @eval function hex(n :: Number) :: String
          string(n, base=16)
        end

  # convert number to binary string
  @eval function bin(n :: Number) :: String
          string(n, base=2)
        end

  # convert number to decimal string
  @eval function dec(n :: Number) :: String
          string(n, base=10)
        end

  # --- and some common functions such as rotation matrix, etc.
end

atreplinit(initrepl)

Notice that I have used the @eval macro for evaluating expressions, if you don’t use the @eval macro, the functions defined in the initrepl function would only be available within the scope of the initrepl function, instead of defined in the global scope. Furthermore, the using statement is not allowed in a function, so we have to use @eval to execute it in the global scope as well.

3.3.2 My startup.jl , take 3

There is an equivalent do -block syntax for defining temporary functions meant to be used only once, so that our temporary function initrepl doesn’t take up a name in the global scope.

atreplinit() do repl
  @eval ENV["JULIA_EDITOR"] = "nvim"

  # packages
  @eval using LinearAlgebra
  @eval using SymPy
  @eval using Latexify

  # convert number to hex string
  @eval function hex(n :: Number) :: String
          string(n, base=16)
        end

  # convert number to binary string
  @eval function bin(n :: Number) :: String
          string(n, base=2)
        end

  # convert number to decimal string
  @eval function dec(n :: Number) :: String
          string(n, base=10)
        end

  # --- and some common functions such as rotation matrix, etc.
end

Now the hello world script will run fast without the flag

$ time julia hello.jl

hello, world
julia hello.jl  0.09s user 0.06s system 99% cpu 0.152 total

and we can still use our defined functions in the REPL .

hex(0b11110) |> print

1e

3.4 PackageCompiler.jl

Using the atreplinit function does reduce the time it takes to run the julia command in non- REPL sessions. However, we are still facing a long start-up time due to SymPy. To solve this issue, we need to create a sysimage for our REPL session using the PackageCompiler.jl package. This will essentially eliminate the JIT overhead during the initial start-up.

To install this package, press ] to enter the package mode and run add PackageCompiler

(@v1.4) pkg> add PackageCompiler

Now, exit the current Julia REPL session, customize your startup.jl file to your liking, and start a new session with flags

$ julia --trace-compile=precompile.jl

This flag will let Julia trace all the JIT compilations that would happen in the new REPL session, including the using statements in the startup.jl file.

If you only want to reduce the time it takes to start up a REPL session, run the exit() command to exit the current session. However, if you also want to reduce the time it takes to run common functions, try simulate your usual workflow in the current REPL session and try to trigger as many pre-compilations as possible.

For example, I would always use the @vars macro to define new symbolic variables in SymPy, so I could run

@vars θ x y real=true

(θ, x, y)

similarly for Latexify,

:(sqrt(θ/3 + exp(ϕ^2))) |> latexraw |> print

\sqrt{\frac{\theta}{3} + e^{\phi^{2}}}

and LinearAlgebra.

eigen([1 2; 3 4])

Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
values:
2-element Array{Float64,1}:
 -0.3722813232690143
  5.372281323269014
vectors:
2×2 Array{Float64,2}:
 -0.824565  -0.415974
  0.565767  -0.909377

Continue like this if you have want to trigger more pre-compilations, but we will exit for now You can also load the functions to run in a file. Refer to the Package­Compiler documentation for details .

exit()

This session should generate a precompile.jl file in the directory you run the julia command.

$ head -n 5 precompile.jl

precompile(Tuple{typeof(Base.atreplinit), Function})
precompile(Tuple{typeof(REPL.Terminals.hascolor), REPL.Terminals.TTYTerminal})
precompile(Tuple{getfield(Main, Symbol("#3#4")), REPL.LineEditREPL})
precompile(Tuple{typeof(MacroTools.__init__)})
precompile(Tuple{typeof(Base.nextpow), Int64, Int64})

This is all the pre-compilation instructions generated from our last session. Now, we could go back to Julia REPL , and run

using PackageCompiler
create_sysimage( [:SymPy, :LinearAlgebra, :Latexify]
               , sysimage_path="sys_calc.so"
               , precompile_statements_file="precompile.jl"
               )

This will use the precompile.jl file we created earlier to create a system image named sys_calc.so in the current directory.

Creating a system image can take some time, so be patient. After the compilation finishes, exit the current REPL session, and start up again using the flag

$ julia -Jsys_calc.so

If everything is done correctly, the REPL should start-up almost instantly, at least much faster than before. If you don’t want to add this flag every time you run Julia, put this file in a fixed directory, for example ~/.julia/lib/sys_calc.so , and add an alias in your .bashrc or .zshrc .

alias jl="julia -J'$HOME/.julia/lib/sys_calc.so'"

or you could instead create a shell script called jl with content

#!/bin/sh

julia -J"$HOME/.julia/lib/sys_calc.so"

and put it into your $PATH You could also use the replace_default option in Package­Compiler to replace the default sysimage, but I would prefer to have separate sysimages for calculator and development .

Now you can use the alias jl to run Julia without the start-up overhead of JIT compilation.

$ jl

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.4.2 (2020-05-23)
 _/ |\__'_|_|_|\__'_|  |
|__/                   |

julia>

4 Tips

The Julia REPL is a beast. There are a ton of small features quite useful for daily usage. Some of my findings are documented here, and I will update the list as I discover more. If you have any other tips, pleasesend me an email and I might consider include some of them.

4.1 Unicode IME

The Unicode input and TAB completion in Julia REPL is quite powerful. I frequently need to type some Unicode math symbols in daily communications when is missing.

There is a nice command called clipboard that basically transformed the Julia REPL as an IME for Unicode characters. We could pipe a symbol or a string into clipboard

:ℝ |> clipboard

and the character ℝ magically appears in your clipboard. Not every character can be used as a symbol, though, so sometimes we need to use strings.

"∃" |> clipboard

We could write a simple macro called @clip to save some typing. You could put it into your startup.jl file.

macro clip(symb)
  :(clipboard($symb))
end

@clip (macro with 1 method)

Then we can type

@clip :ℂ

to copy a ℂ into the clipboard.

By the way, Julia also supports emoji. They usually have prefix \: . For example, \:heart: is ❤ .

@clip :❤

4.2 Scratchpad window

If you are using a window manager such as i 3 , xmonad , or dwm etc., you could set up a scratchpad window for Julia so you can spawn a Julia window with a single key press Mine is bind to Mod+A , and it will keep running in the background so you can hide it and bring it back in with a keyboard shortcut.

This has some added benefits that the JIT pre-compilation of functions are kept in the background, which might save you some time in the long run.

Most window managers has either built-in or third-party support for scratchpad windows, for example,

• i 3
• xmonad
• dwm

Check the documentation for your window manager if you are not sure.

5 Further reading

Since our focus is to use Julia as a simple calculator, we only covered a tiny bit of Julia’s full capabilities. There are so much more, such as data frames , type system , FFI , creating packages , etc. If you want a quick overview of what Julia can do and the syntax we have missed, these cheatsheets could be a good start.

• The Fast Track to Julia
• MATLAB –Python–Julia cheatsheet

We only covered a little bit of the dot syntax, for more on this, check out the blog post on Julia blog

• More Dots: Syntactic Loop Fusion in Julia

I didn’t spend too much time on imperative programming in Julia, because I personally prefer thinking about abstraction and patterns instead of caring about the details of indexing when using Julia as a scratchpad. Many people don’t like Julia’s 1 -based indexing, but we barely used indices in this entire discussion. We are saved by the abstractions of linear algebra, iterators, and functional paradigm to worry about them.

However, if you want to have granular control over performance, manual iterations can be necessary. The following posts might be helpful.

• Multidimensional algorithms and iteration
• Performance Tips
• The Relationship between Vectorized and Devectorized Code by John Myles White

There is also a quite complete course on Julia by Jesse Perla, Thomas J. Sargent, and John Stachurski.

• Quantitative Economics with Julia

Of course, the official documentation for Julia is also quite friendly and complete

• The Julia Language

And for avid Lisp programmers,

• LispSyntax.jl
• Switching from Common Lisp to Julia by Tamás K. Papp

6 Colophon

Besides my usualcolophon, the monospace typeface is Iosevka by Belleve Invis. The Unicode support of this type is quite complete and thus it is a good choice if you want to program in Julia. In particular, I’m using the Extended variant, Fira Mono style ( ss05 ), and modified * ( cv19 ), with ligature off. I don’t want programming ligatures to confuse people, and they are not cool.

The typeface in images is newpx by Michael Sharpe. It is based on pxfonts and Gyre Pagella, both of which are based on URW Palladio L. that is derived from Palatino by Hermann Zapf.

This is my first experiment with the Web layout. It is not written in Jupyter, since the typesetting capability of Jupyter is too limiting for me, but I miss the reproducibility provided by Jupyter. Pleaselet me know if you know an alternative capable of handling my narrative style.

The page shares pretty much the same code base withkarasu, except a different set of filters from fl.tr is applied here. I think the Web layout works quite nicely, and since I can cache all the rendered images and syntax highlighting, the compilation is so much faster than working with real . However, I am still not satisfied with the cross-referencing and widow/orphan control capability on the Web. Maybe I’ll write more filters in the future to fix these problems.

7 License, etc.

The words are licensed under CC BY-NC-SA 4 . 0 , so any non-commercial derivations or citations must be attributed to either Krasjet or my real name, which is readily available if you are willing tocontact me.

Let me know if you want to produce any derived works, e.g. translations, since I think the format and typography are quite tricky to handle, so the source document might be helpful.

The code is licensed separately under the MIT License, so commercial and private use are all allowed.