内容简介:I have been using theI have tried a dozen of alternatives to Julia before, such asOf course, Julia can do simple numerical calculations on numbers, just like any other languages. For example, floating point division
_
_ _ _(_)_ | 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]
Out [1]
1.75
variables, imaginary numbers,
In [2]
x = ans + 2im
Out [2]
1.75 + 2.0im
and mathematical functions ( ).
In [3]
exp(x)
Out [3]
-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 ,
In [4]
A = [1 2 3; 4 5 6]
Out [4]
2×3 Array{Int64,2}:
1 2 3
4 5 6
or alternatively,
In [5]
A = [1 2 3
4 5 6]
Out [5]
2×3 Array{Int64,2}:
1 2 3
4 5 6
and you could do matrix multiplication like this,
In [6]
A * [3; 2; 1]
Out [6]
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 .
In [7]
import numpy as np
A = np.array([[1,2,3],
[4,5,6]])
A @ np.array([[1,2,3]]).T
Out [7]
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 cheetsheets 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.
In [8]
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 .
In [9]
for i in 1:3
println(i)
end
Out [9]
or, alternatively, in Julia’s own syntax.
In [10]
for i ∈ 1:0.5:2
println(i)
end
Out [10]
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 .
In [11]
help?> ℝ
Out [11]
"ℝ" can be typed by \bbR<tab>
You can also find get the documentation for a particular function or variable using the same method.
In [12]
help?> apropos
Out [12]
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 .
In [13]
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 .
In [14]
[a^2 for a ∈ A]
Out [14]
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 .
In [15]
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 .
In [16]
V = [v1 v2 v3]
Out [16]
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.
In [17]
# ℝ → ℝ^(2×2)
R(θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)]
Out [17]
R (generic function with 1 method)
Since we want to rotate by radians, our rotation matrix is
In [18]
R(π/3)
Out [18]
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? .
In [19]
# ℝ^2 → ℝ^2
rotate(v) = R(π/3) * v
Out [19]
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
In [20]
rotate(v1)
Out [20]
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,
In [21]
# ℝ → ℝ^2 (?)
dup = x -> [x;x]
Out [21]
#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,
In [22]
map(dup, V)
Out [22]
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.
In [23]
map(dup, map(dup, V))
Out [23]
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.
In [24]
dup([1 2; 3 4])
Out [24]
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.
In [25]
@which [1;1]
Out [25]
vcat(X::T...) where T<:Number in Base at abstractarray.jl:1277
In [26]
@which [v1;v1]
Out [26]
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 .
In [27]
[V v1; v2 V]
Out [27]
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
.
In [28]
cols = eachcol(V)
Out [28]
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
In [29]
collect(cols)
Out [29]
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 '
.
In [30]
V′ = map(rotate, eachcol(V))
Out [30]
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 .
In [31]
V′ = rotate.(eachcol(V))
Out [31]
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.
In [32]
exp.([1,2,3,4])
Out [32]
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.
In [33]
[1,2,3,4] .* [4,3,2,1]
Out [33]
4-element Array{Int64,1}:
4
6
6
4
It is equivalent to
In [34]
.*([1,2,3,4], [4,3,2,1])
Out [34]
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.
In [35]
# Vector → String
cleanup(v) = join(round.(v, digits=5), ", ")
Out [35]
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 .
In [36]
1 |> exp |> sin
Out [36]
0.41078129050290885
This means the same as
In [37]
sin(exp(1))
Out [37]
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.
In [38]
[1,2,3] .|> exp .|> sin
Out [38]
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.
In [39]
[1,2,3] .|> [exp sin]
Out [39]
3×2 Array{Float64,2}:
2.71828 0.841471
7.38906 0.909297
20.0855 0.14112
This is equivalent to
In [40]
[exp(1) sin(1); exp(2) sin(2); exp(3) sin(3)]
Out [40]
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 .
In [41]
# (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 ->
,
In [42]
(1, 2) .|> x -> x^2 |> sum
Out [42]
(1, 4)
In [43]
(1, 2) .|> (x -> x^2) |> sum
Out [43]
we need to group anonymous functions with parentheses Thanks u/multipleattackers for pointing this out , though for our encrypt
function, this doesn’t matter.
We can encrypt a string by
In [44]
encrypt("hellojulia", 20)
Out [44]
"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.
In [45]
V′ .|>
(x -> round.(x, digits=5)) |>
(x -> join(x, ", ")) |>
println
Out [45]
-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.
In [46]
V′ .|>
(x -> round.(x, digits=5)) |>
(x -> join(x, ", ")) |>
println;
Out [46]
-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
In [47]
map (3*) [1,2,3]
Out [47]
[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,
In [48]
curry = Base.Fix2
map(curry(*,3), [1,2,3])
Out [48]
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.
In [49]
eachcol(V) .|> rotate .|> cleanup .|> println;
Out [49]
-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}
,
In [50]
A = [1 2 3; 4 5 6]
Out [50]
2×3 Array{Int64,2}:
1 2 3
4 5 6
In [51]
Array{Integer,2} == Matrix{Integer}
Out [51]
true
and a list of vectors, which has type Array{Array{T,1},1} = Vector{Vector{T}}
,
In [52]
B = [[1, 4], [2, 5], [3, 6]]
Out [52]
3-element Array{Array{Int64,1},1}:
[1, 4]
[2, 5]
[3, 6]
In [53]
Array{Array{Integer,1},1} == Vector{Vector{Integer}}
Out [53]
true
They are, of course, not equivalent.
In [54]
A == B
Out [54]
false
In [55]
Matrix{Integer} == Vector{Vector{Integer}}
Out [55]
false
To understand this, it would be easier to look at the distinction between the following three syntax.
In [56]
# List (Array)
[1,2,3]
Out [56]
3-element Array{Int64,1}:
1
2
3
In [57]
# Column vector
[1;2;3]
Out [57]
3-element Array{Int64,1}:
1
2
3
In [58]
# Row vector
[1 2 3]
Out [58]
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]
.
In [59]
[1,2,3] == [1;2;3]
Out [59]
true
The range syntax also produces a list, and thus a column vector.
In [60]
1:3 == [1;2;3]
Out [60]
true
This is why Julia is not compatible with the matrix syntax of MATLAB . This is valid in MATLAB ,
In [61]
[1,2,3; 4,5,6]
Out [61]
ans =
1 2 3
4 5 6
but not valid in Julia.
In [62]
[1,2,3; 4,5,6]
Out [62]
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.
In [63]
[1 2; 3 4] * [1,1]
Out [63]
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 .
In [64]
[[1,2], [3,4], [5,6]]
Out [64]
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,
In [65]
[[1,2]; [3,4]; [5,6]]
Out [65]
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
).
In [66]
[[1 2] [3 4] [5 6]]
Out [66]
1×6 Array{Int64,2}:
1 2 3 4 5 6
Compare with the result using the list syntax,
In [67]
[[1 2], [3 4], [5 6]]
Out [67]
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.
In [68]
[[1;2] [3;4] [5;6]]
Out [68]
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.
In [69]
# A list of a single number
a = [1]
Out [69]
1-element Array{Int64,1}:
1
In [70]
# More lists of a single number
b = [2]
c = [3]
d = [4]
In [71]
# Vertical concatenation of horizontal concatenations of lists
# of lists of a single number
[[a] [b]; [c] [d]]
Out [71]
2×2 Array{Array{Int64,1},2}:
[1] [2]
[3] [4]
In [72]
# Vertical concatenation of lists of a single number
[a; b]
Out [72]
2-element Array{Int64,1}:
1
2
In [73]
# 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]
Out [73]
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.
In [74]
Meta.show_sexpr(:([1,2,3]))
Out [74]
(:vect, 1, 2, 3)
In [75]
Meta.show_sexpr(:([1;2;3]))
Out [75]
(:vcat, 1, 2, 3)
In [76]
Meta.show_sexpr(:([1 2 3]))
Out [76]
(: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.
In [77]
'vcat
Out [77]
'vcat
In [78]
'1
Out [78]
and in Julia
In [79]
:vcat
Out [79]
:vcat
In [80]
Out [80]
The notation :(1 + 1)
denotes an expression , where the :()
syntax is similar to quasiquotes in Lisps,
In [81]
`(+ 1 2)
Out [81]
'(+ 1 2)
In [82]
`(+ ,(+ 1 3) 2)
Out [82]
'(+ 4 2)
and in Julia The splicing syntax works on strings as well
In [83]
:(1 + 2)
Out [83]
:(1 + 2)
In [84]
:($(1 + 3) + 2)
Out [84]
:(4 + 2)
We can evaluate an expression using the eval
function.
In [85]
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.
In [86]
@. [1,2,3] |> exp |> sin
Out [86]
3-element Array{Float64,1}:
0.41078129050290885
0.8938549549128102
0.9444710089262849
If we take a look at the syntax for defining matrices,
In [87]
Meta.show_sexpr(:([1 2 3; 4 5 6]))
Out [87]
(: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
In [88]
@edit [1 2 3; 4 5 6]
Out [88]
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,
In [89]
Meta.show_sexpr(:([[1,4], [2,5], [3,6]]))
Out [89]
(:vect, (:vect, 1, 4), (:vect, 2, 5), (:vect, 3, 6))
and the function being executed is vect
.
In [90]
@edit [[1,4], [2,5], [3,6]]
Out [90]
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 .
In [91]
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
.
In [92]
eachslice(V, dims=2) |> collect
Out [92]
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 .
In [93]
V
Out [93]
2×3 Array{Float64,2}:
0.5 1.5 0.5
0.5 0.5 1.5
You can query this directly.
In [94]
row, col = size(V)
Out [94]
(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,
In [95]
V[1,2]
Out [95]
1.5
and we could select multiple elements by passing a list as index You can also use 1:2
or 1:end
here .
In [96]
V[[1,2], 2]
Out [96]
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.
In [97]
V[[1 2], 2]
Out [97]
1×2 Array{Float64,2}:
1.5 0.5
or even a matrix.
In [98]
V[[1 2; 2 1], 2]
Out [98]
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.
In [99]
map(i -> V[i,2], [1 2; 2 1])
Out [99]
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.
In [100]
V[:,2]
Out [100]
2-element Array{Float64,1}:
1.5
0.5
This is exactly our . Hence, we could write the matrix V
as
In [101]
[V[:,1] V[:,2] V[:,3]] == V
Out [101]
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
In [102]
map(i -> V[:,i], 1:3)
Out [102]
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.
In [103]
mapslices(rotate, V, dims=1)
Out [103]
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 .
In [104]
mapslices(sum, V, dims=2)
Out [104]
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.
In [105]
T = rand(3,3,3)
map(i -> T[:,i,:], 1:3) == eachslice(T, dims=2) |> collect
Out [105]
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.
In [106]
for v ∈ eachcol(V)
v |> rotate |> cleanup |> println
end
Out [106]
-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.
In [107]
[v |> rotate |> cleanup for v ∈ eachcol(V)]
Out [107]
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,
In [108]
for i ∈ 1:size(V,2)
V[:,i] |> rotate |> cleanup |> println
end
Out [108]
-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.
In [109]
# 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
.
In [110]
# 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.
In [111]
# 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
In [112]
for i ∈ 1:size(V,2)
view(V,:,i) |> rotate |> cleanup |> println
end
Out [112]
-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.
In [113]
@views for i ∈ 1:size(V,2)
V[:,i] |> rotate |> cleanup |> println
end
Out [113]
-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.
In [114]
R(π/3) * V
Out [114]
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 .
In [115]
ans'
Out [115]
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
In [116]
using LinearAlgebra
det(R(π/3))
Out [116]
1.0
or you could show off the pipe operator by writing
In [117]
π/3 |> R |> det
Out [117]
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 .
In [118]
V .- V[:,1]
Out [118]
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.
In [119]
V .+ V[1,:]'
Out [119]
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
In [120]
V[1,:]'
Out [120]
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
In [121]
V .+ V[1,1]
Out [121]
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.
In [122]
R(π/3) * (V .- V[:,1]) .+ V[:,1]
Out [122]
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.
In [123]
c = V * [1;1;1]/3
Out [123]
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
In [124]
reduce(+, V / 3, dims=2)
Out [124]
2×1 Array{Float64,2}:
0.8333333333333333
0.8333333333333333
The argument dims=2
means the reduce
operation smashes the second dimension
In [125]
V[1,:]
Out [125]
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
In [126]
reduce(*, V, dims=2).^(1/size(V,2))
Out [126]
2×1 Array{Float64,2}:
0.7211247851537042
0.7211247851537042
To finish our example, let us calculate the rotated vertices,
In [127]
R(π/3) * (V .- c) .+ c
Out [127]
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,
In [128]
translate(t) = v -> v + t
rotate(θ) = v -> R(θ)*v
eachcol(V) .|> translate(-v1) .|> rotate(π/3) .|> translate(v1)
Out [128]
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
.
In [129]
(@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 .
In [130]
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.
In [131]
@basis S"∞∅++" G31 e
Out [131]
(⟨∞∅++⟩, 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∞
.
In [132]
e1 + 2*e2
Out [132]
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
.
In [133]
[e1 e2] * V
Out [133]
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
In [134]
W = ans |> transpose
Out [134]
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
In [135]
↑(0.5*e1 + 0.5*e2)
Out [135]
-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.
In [136]
= e12
Out [136]
v₁₂
In [137]
^2
Out [137]
-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
In [138]
↓(↑(0.5*e1 + 0.5*e2) ⊘ exp( *(π/3)/2))
Out [138]
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
In [139]
# rotation angle
θ = π/3
# displacement vector
t = W[1]
.↑(W) .⊘ exp(-t*e∞/2) .⊘ exp( *θ/2) .⊘ exp(t*e∞/2) .|> ↓
Out [139]
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
In [140]
@. ↑(W) ⊘ exp(-t*e∞/2) ⊘ exp( *θ/2) ⊘ exp(t*e∞/2) |> ↓
Out [140]
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
In [141]
T = exp(-t*e∞/2) * exp( *θ/2) * exp(t*e∞/2)
Out [141]
0.8660254037844387 - 0.24999999999999997v∞₁ +
0.24999999999999994v∞₂ + 0.49999999999999994v₁₂
and apply it to each element in W
using versor product directly.
In [142]
@. ↑(W) ⊘ T |> ↓
Out [142]
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.
In [143]
(@v1.4) pkg> add UnicodePlots
In [144]
using UnicodePlots
Now, we can draw the same diagram as before directly in the terminal using the lineplot
function
In [145]
# 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
)
Out [145]
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢸⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢸⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⢸⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠸⠤⠤⠤⠤⠤⠤⠵⠄⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⡧⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
We can apply the rotation in any of the approaches described before, but for simplicity, we will just use the linear algebra approach.
In [146]
U′ = R(π/3) * U
Out [146]
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
In [147]
lineplot( U′[1,:], U′[2,:] # x, y
, xlim=(-2,2), ylim=(-2,2)
, border=:none, labels=false
, width=30
)
Out [147]
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣇⣀⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⢀⣀⡠⠤⠒⠒⠉⠉⣇⠜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠉⠒⢄⡀⠀⠀⠀⢀⡏⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠑⠤⣠⠃⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⡧⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
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 Windows
~/.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
,
In [148]
@vars θ x y real=true
Out [148]
(θ, x, y)
and perform operations on them directly,
In [149]
w = R(θ) * [x;y]
Out [149]
2-element Array{Sym,1}:
x*cos(θ) - y*sin(θ)
x*sin(θ) + y*cos(θ)
In [150]
using LinearAlgebra
norm(w) |> simplify |> print
Out [150]
sqrt(x^2 + y^2)
3.1.2 Latexify
The package Latexify.jl is for printing the commands for a Julia expression directly.
In [151]
V |> latexify |> print
Out [151]
\begin{equation}
\left[
\begin{array}{ccc}
0.5 & 1.5 & 0.5 \\
0.5 & 0.5 & 1.5 \\
\end{array}
\right]
\end{equation}
In [152]
:(sqrt(θ/3 + exp(ϕ^2))) |> latexraw |> print
Out [152]
\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.
In [153]
bin(0xce - 5)
Out [153]
In [154]
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,
In [155]
$ time julia hello.jl
Out [155]
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
In [156]
$ time julia --startup-file=no hello.jl
Out [156]
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
In [157]
@time println("hello, julia")
Out [157]
hello, julia
0.002414 seconds (24 allocations: 1.734 KiB)
In [158]
@time println("hello, julia")
Out [158]
hello, julia
0.000045 seconds (4 allocations: 80 bytes)
In [159]
@time println("bye, julia")
Out [159]
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.
In [160]
@time println([1,2,3,4])
Out [160]
[1, 2, 3, 4]
0.110054 seconds (243.47 k allocations: 11.987 MiB)
In [161]
@time println([1,2,3,4])
Out [161]
[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
In [162]
$ time julia hello.jl
Out [162]
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 .
In [163]
hex(0b11110) |> print
Out [163]
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
In [164]
@vars θ x y real=true
Out [164]
(θ, x, y)
similarly for Latexify,
In [165]
:(sqrt(θ/3 + exp(ϕ^2))) |> latexraw |> print
Out [165]
\sqrt{\frac{\theta}{3} + e^{\phi^{2}}}
and LinearAlgebra.
In [166]
eigen([1 2; 3 4])
Out [166]
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 PackageCompiler documentation for details .
In [167]
exit()
This session should generate a precompile.jl
file in the directory you run the julia
command.
In [168]
$ head -n 5 precompile.jl
Out [168]
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
In [169]
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
In [170]
$ 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 PackageCompiler 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.
In [171]
$ jl
Out [171]
_
_ _ _(_)_ | 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
In [172]
:ℝ |> 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.
In [173]
"∃" |> clipboard
We could write a simple macro called @clip
to save some typing. You could put it into your startup.jl
file.
In [174]
macro clip(symb)
:(clipboard($symb))
end
Out [174]
@clip (macro with 1 method)
Then we can type
In [175]
@clip :ℂ
to copy a ℂ
into the clipboard.
By the way, Julia also supports emoji. They usually have prefix \:
. For example, \:heart:
is ❤
.
In [176]
@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,
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.
We only covered a little bit of the dot syntax, for more on this, check out the blog post on Julia blog
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.
Of course, the official documentation for Julia is also quite friendly and complete
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.
以上就是本文的全部内容,希望本文的内容对大家的学习或者工作能带来一定的帮助,也希望大家多多支持 码农网
猜你喜欢:本站部分资源来源于网络,本站转载出于传递更多信息之目的,版权归原作者或者来源机构所有,如转载稿涉及版权问题,请联系我们。