let
= 0
p let p = 1 # New version of p
println(p) # Prints 1 (p from inner scope)
end
println(p) # Prints 0 (p from outer scope)
end
1
0
Entering this chapter, we assume you have written code in some other languages, but not programmed much (or at all) in Julia, the language that we will use for the code in this book. Julia is a relatively young language initially released in 2012. In comparison the first releases of MATLAB and Python were 1984 and 1991, respectively. Despite its relative youth, Julia has become popular in scientific computing for its speed, simple MATLAB-like matrix syntax, and support for a variety of programming paradigms. Among other reasons, we will use Julia for our codes because
Julia has strong support for interactive computation, similar to Python, R, or MATLAB. I can try things out on a regular command line, or inside a code editor like VSCode, or in a notebook environment (using Jupyter or Pluto).
Carefully written Julia code can run as fast as codes written in C, C++, or Fortran.
Like MATLAB or modern Fortran, Julia provides a concise and convenient syntax for linear algebra computations. Julia also provides native support for Unicode, which is nice for writing mathematics in code the same as one would write it on paper.
Julia has an easy-to-use package manager and a growing ecosystem of useful package. Because so many users of Julia come from scientific computing and data science, there tend to be many good packages for numerical computing.
Julia treats functions (and closures) as first-class objects and supports a variety of functional programming features.
Julia has an expressive type system that is useful for code correctness, for writing expressive code via multiple dispatch (different specialized methods with the same name but different type signatures), and for fast execution.
Julia provides a rich metaprogramming facility, explicitly modeled on Lisp macros. This makes it easy to write Julia programs that write programs, and can be a powerful way to extend the language.
The rest of this chapter is intended to get you started with Julia. It is not comprehensive, and we recommend that the curious reader explore further with the excellent Julia documentation. But we also includes some treatment of more advanced corners of the language.
We can use Julia non-interactively: typing julia filename.jl
at the command line will execute any Julia code in the text file filename.jl
. But both for developing new code and for learning the language, it is very useful to work with Julia interactively. The two most common approaches are
Typing julia
at the command line will start an interactive Julia session, sometimes known as a read-evaluate-print loop (REPL). One can also start an interactive Julia session from within in some editors, such as Visual Studio.
One can interact with Julia using a notebook interface, such as Jupyter (likely the most popular) or Pluto. We view notebooks using a web browser.
Running julia
with no arguments takes us to the Julia prompt. From here, we can type and run general Julia code, though for anything too long we would usually write the code in a file. But even typing simple commands can have some surprises!
The Julia prompt supports many of the amentities of modern command prompts, including:
History: One can scroll through past commands using the up and down arrows.
Emacs-style editing: Code typed into the command prompt can be edited using navigation key bindings similar to the Emacs editor. For example, Ctrl-A and Ctrl-E jump to the beginning and end of the current line, Ctrl-K “kills” text (putting it into a buffer called the “kill ring”), and Ctrl-Y retrieves the last line from the kill ring.
Tab completion: Typing a partial command into Julia and then pressing the tab key will lead to Julia either finishing the command (if it is unambiguous). If the completion is ambiguous, pressing tab a second time will produce a list of suggested possible completions (if the completion is ambiguous).
Julia allows us to use Unicode characters in code, including code written at the command line. The simplest way to enter these characters is to type a LaTeX special character command and then tab; for example \pi
followed by tab produces the character \(\pi\) (and in Julia, as in mathematics, this represents \(3.14159...\) unless superceded by another local definition). While it is possible to use Julia without taking advantage of Unicode, in many cases it is the most concise way of writing expressions. For example, writing isapprox(x, y)
or x ≈ y
(rendering the approximately equal symbol by tab-completion on \approx
are equivalent in functionality. But the latter is certainly shorter and arguably more idiomatic. Several popular editors also have Julia modes that support Unicode via tab-completion of LaTeX commands, including Visual Studio, Emacs, and Vim.
From the Julia prompt, one can also access
Help: We access the help system by typing the ?
character and then the name of a function (or module, global variable, etc). Tab completion works in the help system as well, both for completing command names and for entering Unicode. For example, typing ?\pi
followed by a tab will get system help on the constant \(\pi\).
Shell: If we would like to list directories or quickly run other programs, we can access the shell by typing the ;
character. To exit shell mode, we type the backspace character at an otherwise-empty shell prompt.
Packages: To enter the package manager, we type the ]
character; when done, we can exit the package manager by typing a backspace, as with shell mode. We will have more to say about the package manager later in this chapter.
Pressing Ctrl-C lets us break out of a running Julia command. We can exit from Julia by typing exit()
at the prompt or pressing Ctrl-D.
The Jupyter notebook system was first released in 2015 as a multi-language extension to the earlier IPython notebooks. The Jupyter system supports “computational notebooks” built from cells containing text documentation written in Markdown; code written in Julia or another supported language; and code outputs, which may include both text and graphics. Early versions of Jupyter supported Julia, Python, and R (hence the name),though support for many other languages has been added since. In 2018, a more full-featured environment called JupyterLab was released, and today lives alongside the “classic” Jupyter notebook interface. Both classic Jupyter and JupyterLab use the same file formats .ipynb
files), which may contain all cell types, including computed outputs.
Users typically interact with Jupyter through a web browser interface, though there is also support for working with Jupyter notebooks in editors like Visual Studio Code and Emacs. An example of a running notebook is shown in ?fig-tbd. The actual computations for a Jupyter noteook run in a separate “kernel” that communicates with the web browser via a network protocol. The kernel and the brower interface are separate; they may run on the same computer, but are not required to. The kernel responds as though individual cells were typed into the command prompt in the order that the user chooses to execute them, with output then redirected to be shown in the web browser. Because the order of cell execution is chosen by the user, it is not always possible to tell how a notebook arrived at some result just by looking at the output: the user might have executed cells in some unintuitive order, or perhaps a cell was run and then deleted. For this reason, when preparing to save a notebook (whether to share with collaborators or because it is time to quit for the day), we recommend restarting the kernel and re-running all cells from the start so that “what you see is what you get.”
The Julia kernel for Julia is provided in the IJulia
package. This can be installed from the Julia command prompt by entering the package mode (by typing ]
at the prompt) and typing
add IJulia
Once the IJulia
package is installed and we have exited back to the normal Julia prompt, we can run Jupyter with a Julia kernel from the Julia prompt by typing
using IJulia
notebook() IJulia.
Running Jupyter in this way results in a new installation of the Jupyter system, independent of any installations of Jupyter that may already exist on your computer. If you already have Jupyter installed (e.g. as part of a Python installation), you can tell it about the IJulia
kernel by making sure the JUPYTER
environment variable points to the desired jupyter
executable (e.g. by assigning the path to ENV["JUPYTER"]
from the Julia prompt), and either running build IJulia
at the package prompt, or at the Julia prompt running
using IJulia
installkernel()
Once the kernel is installed, one can start Jupyter as normal (e.g. by typing jupyter lab
at the shell), and then select the Julia kernel when starting a new notebook or loading an existing one.
The Pluto notebook is a Julia-specific take on the computational notebook concept. As with Jupyter, the user interacts with Julia computations through a web interface, with computations executing in a separate Julia process. However, there are several key differences between the two systems:
Pluto does not distinguish between “code cells” and “documentation cells.” The Pluto equivalent of a documentation cell is a code cell that renders Markdown text, which is then displayed as the cell output.
A single Pluto cell can only contain one statement, though this includes things like function definitions or begin
/end
blocks.
Pluto maintains a dependency graph between cells that is used to determine execution order. If cell A defines a symbol that is used within cell B, then executing any update to cell A will also trigger an update of cell B, along with any other cells that depend on cell A, whether directly or indirectly. This mostly addresses the issue of uncertain execution order that sometimes causes headaches for Jupyter users. However, the dependency graph construction is not perfect, and can sometimes be confused by too-clever macros. When this happens, or when there is a change in the packages used by a Pluto notebook, it is usually a good idea to recompute the whole notebook from scratch.
Pluto notebooks are saved as regular Julia files, with special comments to delimit the cells and indicate the order in which they are displayed. Within the file, the cells are ordered according to the depencency graph, so that a Pluto notebook can also be run as an ordinary Julia program.
We start our deeper dive into Julia with a treatment of simple expressions, the very simplest of which are literals (Table 2.1).
Type | Examples |
---|---|
Int |
123 , 1_234 |
UInt |
0x7B , 0o173 0b1111011 |
Float64 |
10.1 , 10.1e0 , 1.01e1 |
Float32 |
10.1f0 , 1.01f1 |
Logical |
true , false |
Char |
'c' |
String |
"string" |
Symbol |
:symbol |
Ordinary decimal integer expressions (e.g. 123
) are interpreted as type Int
, the system signed integer type. In principle, this will be 32-bit or 64-bit depending on the preferred default integer size for the system. In practice, most current systems are 64-bit. Unsigned integer literals can be written in hexadecimal, octal, or binary; these listerals default to the smallest type that will contain them. Floating point literals default to 64-bit, unless an exponent is included. For both integer and floating point types, underscores can be included on input as spacers; that is, 1_234
and 1234
are interpreted the same way.
In addition to standard numerical literals, Julia provides a few constants, including \(\pi\) and \(e\); we type Unicode symbols for these with tab completion on \pi
and \euler
. There are also special constants for exceptional floating point values. For example, Inf
or Inf64
denote double-precision floating point infinity and NaN
or NaN64
denote double-precision not-a-number encodings; the single-precision (32-bit) analogues are Inf32
and NaN32
.
Julia supports string literals surrounded by either single quotes or triple quotes. Triple-quoted strings can include regular quote characters and newlines with no special characters, e.g.
= """
box_quote "All models are wrong, but some models are useful."
-- George Box"""
Julia characters correspond to Unicode codepoints, and strings likewise use Unicode encodings. In particular, Julia string literals are represented in program code and in memory using UTF-8, which defines a variable-length encoding for characters. This leads to some complexities in Julia string processing, which we address in Section 2.5.3.
A symbol in Julia is represented in memory as a special type of string that is entered in a global dictionary for fast comparisons. Symbols can be used like strings that admit particularly fast comparison, but they are more important for the role that they play in metaprogramming, where evaluating a symbol corresponds to evaluating a variable named in the program.
Like many other languages, the Julia assignment operation =
represents binding a name to a value. For example,
= 10 x
assigns a name x
to the value 10
. Julia also allows destructuring assignments in which an ordered collection (a tuple or vector) is unpacked into multiple outputs, e.g.
= lu(A) # Assign multiple outputs L, U, p
There is also a form of destructuring assignment for when the right hand side has named quantities, e.g.
= (x = 1.0, y = 2.0, z = 3.0)
p = p # Sets y = p.y, z = p.z (; y, z)
This format works for named tuples and for structures.
What is more complicated is the scope, i.e. the part of the code in which this assignment holds.
Global variables in Julia belong to a “top-level” namespace called the global scope. Technically, there is a separate global scope for each module, as we will describe in more detail in ?sec-tbd. But when we start a Julia session, we are working with a new (anonymous) module with its own global scope, and when we refer to “the” global scope, this is the namespace that we mean.
Certain types of code blocks create new local scopes, or variable namespaces with bindings that are only meaningful inside the code block. Julia is lexically scoped, so the visibility of variable bindings depends on the code structure, not the execution path through the code. Local scopes can (and often are) nested. The let
statement in Julia serves the sole purpose of creating a new scope with local bindings; for example:
let
= 0
p let p = 1 # New version of p
println(p) # Prints 1 (p from inner scope)
end
println(p) # Prints 0 (p from outer scope)
end
1
0
The let
statement actually has two parts:
let
, which result in bindings of variable names that are purely local to the let
;The difference can be somewhat subtle. For example, a single newline can change the semantics of the previous example by moving the assignment p=0
from the assignment block of the let
(where it creates a new version of p
) to the code block of the let
(where it re-binds the symbol p
from the enclosing scope):
let
= 0
p let
= 1 # Reassign outer version of p
p println(p) # Prints 1
end
println(p) # Prints 1 again
end
1
1
At the same time, we can use the local
keyword to make a version of the variable that is explicitly local to the let
code block:
let
= 0
p let
local p = 1 # Create local version of p
println(p) # Prints 1
end
println(p) # Prints 0
end
1
0
Julia distinguishes between constructs that create “hard” local scopes (such as function and macro definitions or the assignment part of let
blocks) and those that create “soft” local scopes (such as loops, conditionals, the code part of let
statements, and try
blocks), with the two differing in their treatment of assignments to existing names. For example, a for
statement establishes a soft local scope, in which assignments to existing variables in enclosing scopes do not create a shadow variable:
function do_more_stuff()
= 0
x for k = 1:5
= k # x is the same variable as before
x end
println(x) # Prints 5
end
However, if x
were not declared before the for
loop, the version inside the loop body would be purely local to the loop:
function do_more_stuff()
for k = 1:5
= k # x is local to this iteration
x end
println(x) # Error -- x doesn't exist here!
end
In any local scope, we can explicitly declare that we want to use a global variable. For example:
function modify_global_xyz()
= 0
xyz let
global xyz
= 1
xyz end
end
modify_global_xyz()
println(xyz) # Prints 1
1
However, if we try to assign to a global variable name in a local scope without using the global
keyword, different things will happen depending not only on whether the local scope is “hard” or “soft,” but also on whether we are in an interactive session or running code non-interactively:
Hard local scope: Always create a shadow of the global variable.
Soft local scope: If interactive, assign to the global variable. Otherwise, create a shadow variable and issue a warning to the user.
All values in Julia have concrete types. When a variable in Julia is given a type, it restricts the types of values that can be assigned to the variable to only those that are compatible. For example, if we would like x
to represent a floating point number, we could write:
:: Float64 = 0 # Implicitly converts the integer 0 to the floating point 0.0 x
Once the type of x
has been specified, we cannot assign an incompatible value to x
:
:: Float64 = 0 # OK so far, x is now floating point 0.0
x = "Oops" # Error! x
Though values must have concrete types, we can give variables abstract types. For example, Julia has many numeric types, all of which are subtypes of type Number
. We can assign any numeric value to a variable declared to have type Number
, but non-numeric values are not allowed:
:: Number = 0 # OK, x is now Int(0), and Int is a subtype of Number
x = 0.0 # Still OK, Float64 is a subtype of Number
x = "Oops" # Error! x
We will have much more to say about types in ?sec-tbd.
Syntax | Description |
---|---|
+x |
Unary plus (identity) |
-x |
Unary minus |
x + y |
Addition |
x - y |
Subtraction |
x * y |
Multiplication |
x / y |
Division |
x ÷ y |
Integer division |
x \ y |
Inverse divide (\(\equiv\) y/x ) |
p // q |
Rational construction |
x ^ y |
Power |
x % y |
Remainder |
We show some of the standard Julia arithmetic operations in Table 2.2. These mostly behave similarly to other languages, but with some subtleties around division. In particular:
The ordinary division operator produces a floating point result even when the inputs are integers and the result is exactly representable as an integer. For example, 6/3
yields 2.0
.
Integer division always rounds toward zero; for example, 14 ÷ 5
yields 2
, while -14 ÷ 5
yields -2
. The remainder operation is consistent with integer division.
Integer division and remainder can also be used with non-integer types; for example, 14.1 ÷ 5
yields 2.0
, and 14.1 % 5
yields 4.1
.
The backslash (inverse divide) operation can be used with any numeric types, but is most frequently used in the context of linear algebra. When A
and b
are a matrix and a vector, we use A\b
to compute \(A^{-1} b\) without computing an explicit inverse.
The construct p // q
is used to construct a rational number. The numerator and denominator must be integers.
Julia also allows for implicit multiplication when a literal number (integer or floating point) precedes a symbol with no space. For example, 2.0im
produces a complex floating point value equal to \(2i\).
There are updating operations of all the arithmetic operations (except for rational construction), formed by combining the operation name with an equal sign. These are essentially equivalent to a binary operation followed by an assignment; for example,
+= 1 # Equivalent to x = x + 1
x *= 2 # Equivalent to x = x * 2 x
Julia supports the usual logical and p && q
, or p || q
, and !p
constructions. The arguments to these operations must be logical; unlike some languages, we cannot interpret integers or lists as logical values, for example. The and and or operators are short-circuited, meaning that the second argument is only executed if the first argument is insufficient to determine the truth value. For example, consider the code fragment
let
= true || begin println("Hello 1"); true end
x = false || begin println("Hello 2"); true end
y
x, yend
Hello 2
(true, true)
Both x
and y
are assigned to be true, but only the string "Hello 2"
is printed.
Julia also supports a ternary conditional operator: p ? x : y
evaluates p
and returns either x
or y
depending whether p
is true or false, respectively. Like the logical and and or, the ternary conditional operator is short-circuiting; only one of the x
and y
clauses will be evaluated. The conditional operator is equivalent to a Julia if
statment; that is, these two statements do the same thing:
= p ? x : y
z1 = if p then x else y end z2
Julia provides several different comparison operations. Most types support tests for structural equality (==
) or inequality (!=
or ≠
). This is a forgiving sort of test, allowing for things like promotion of the arguments to a common type. For example, all the following expressions evaluate to true
:
1 == 1 &&
1 != 2 &&
2 == 2.0 &&
1; 2] == [1.0; 2.0] [
true
The substitutional equality operation (===
) and inequality operation (!==
) test whether two entities are indistinguishable. This is usually a more restrictive notion of equality. For example, all the following expressions evaluate to true
:
1 === 1 &&
1 !== 2 &&
2 !== 2.0
true
Values that correspond to mutable objects never satisfy ===
, unless they are exactly the same object. For example, Julia allows us to change the contents of an array (it is mutable), but not the contents of a string (it is immutable). Therefore, both the following expressions are true
:
"12" === "12" &&
1; 2] !== [1; 2] [
true
For the most part, substitutional equality implies structural equality, but not vice-versa. We say “for the most part” because of the example of the default not-a-number (NaN
). According to the standard floating point semantics (Chapter 9), not-a-number is not equal to anything in the floating point system, including itself. But the NaN
symbol does satisfy substitutional equality. Therefore, both the following are true
:
NaN != NaN &&
NaN === NaN
true
We will have more to say about the vagaries of floating point equality when we get to Chapter 9. For now we simply mention that Julia also supports an approximate equality operation that can be used for comparing whether floating point quantities are within some tolerance of each other. For example, the following statements are both true
:
1.0 ≈ 1.00000001 &&
1.0 ≠ 1.00000001
true
In addition to testing for equality or inequality, we can test ordered types with <
, >
, <=
(or ≤
), and >=
(or ≥
). Ordered types include integers, floating point numbers, strings, symbols, and logicals. For strings and symbols, we use lexicographic ordering (aka alphabetic ordering); for logicals, we have false < true
.
:a < :b &&
"a" < "b" &&
false < true
true
Simple calls in Julia have the form function_name(args)
. The arguments are comma-separated, starting with any positional arguments followed by any keyword arguments. Julia uses the types of the provided positional arguments to distinguish between different methods with the same name and number of arguments. Most operators in Julia are just a thin layer of “syntactic sugar” that hide ordinary method calls. For example, writing 1+2
and +(1,2)
are equivalent; but 1+2
and 1.0+2.0
invoke different methods, since the types are different (both Int
for the former expression, both Float64
in the latter).
A call may provide only the leading positional arguments if the remaining arguments are assigned default values. Keyword arguments always have a default value. For example, the linear algebra routine cholesky
takes positional arguments: a matrix A
and a flag variable that indicates whether the computation should use pivoting or not (defaulting to NoPivot()
). In addition, Cholesky takes a logical keyword argument check
that indicates whether the routine should produce an exception on failure (check == true
, the default) or should produce no exception (check == false
). Therefore, all the following are valid calls to cholesky
= cholesky(A) # without pivoting, check result
C = cholesky(A, check=false) # without pivot, no check
C = cholesky(A, RowMaximum()) # pivoted, check result
C = cholesky(A, RowMaximum(), check=false), # pivot, no check C
Julia passes arguments using a “call-by-sharing” convention, in which names within the function are bound to the values passed in. This means that if a mutable value (e.g. an array) is passed as an argument to a function, the function is allowed to modify that argument. By convention, functions that modify their arguments end with an exclamation mark. For example, calling cholesky(A)
function produces a new output matrix and leaves A
alone, while cholesky!(A)
modifies the matrix A
.
Sometimes, it is convenient to package the positional arguments to a function into a tuple. A tuple can be unpacked into an argument list by following it with ...
, sometimes called the splat operation. This is sometimes useful in function calls; for example
= (1.0, 1.00000001)
a isapprox(a...) # Same as isapprox(1.0, 1.00000001)
While it is mostly used for calling, Julia does also allow this sort of unpacking in other contexts than function calls; for example
= (a, 2) # Yields ((1.0, 1.00000001), 2)
b = (a..., 2) # Yields (1.0, 1.00000001, 2) bunpack
Function calls and operations can be applied elementwise to arrays using broadcasting syntax. This involves putting a dot before the operator name or after the function name. For example: {.julia} θs = [0; π/4; π/2; 3π/4; π] cos.(θs) # [1; sqrt(2)/2; 0; -sqrt(2)/2; -1] θs .+ 1 # [1; π/4+1; π/2+1; 3π/4+1; π+1]
The @.
macro can be used to “dot-ify” everything in an expression, e.g. {.julia} xs = range(0, 1, length=100) y = @. xs^2 + 2xs + 1.0
evaluates the map \(x \mapsto x^2 + 2x + 1\) to every element of the vector xs
.
Vectorized operations are critical to performance in some languages. In Julia, they are not so performance critical. However, they are convenient for concise code.
Square brackets are used to access elements of an array or similar object,or to refer to sub-arrays of an array (slices). Indices are one-based; and in the context of indexing, the symbol end
refers to the final index. For example,
= [5; 2; 2; 3]
lst 2] # Evaluates to 2
lst[end] # Evaluates to 3
lst[1:2] # Evaluates to [5; 2]
lst[2; 1]] # Evaluates to [2; 5]
lst[[2:end] # Evaluates to [2; 2; 3]
lst[2:3] .= 4 # Now lst is [5; 4; 4; 3] lst[
We can index into an array to either get elements or subarrays or to assign to them. When getting elements or subarrays, the default behavior is to get a copy of the extracted data; so, for example
= [5; 2; 2; 3]
lst = lst[1:2] # New storage for sublst
sublst :] = [3; 4] # Changes contents of sublst, not lst sublist[
If we want to modify a subarray in the main array from which it was extracted, we can use a view. A view has its own indexing, but uses the storage of another object. For example:
= [5; 2; 2; 3]
lst = view(lst, 1:2) # Provides a view into lst
sublst :] = [3; 4] # lst is now [3; 4; 2; 3] sublist[
We frequently provide indices that are integers (to get individual elements) or arrays of integers (whether ranges like 1:2
or concrete index lists like [1 4 2]
). But we can also use logical indexing, in which an iterable boolean list tells us whether the corresponding entry should be kept in the result or excluded. For example:
= [5; 2; 2; 3]
lst = (lst .> 2) # [true, false, false, true]
idx_big # Results in [5; 3] lst[idx_big]
In some circumstances, we also have index spaces that are fundamentally non-integer (or can be). An example is with named tuples, which can be accessed either by indexing positionally or by name. For example, consider the following tuple:
= (x=4, y=5, z=6) # Create a named tuple
xyz 1] # Returns 4 (the value of x)
xyz[:y] # Returns 5 (the value of y) xyz[
Our examples so far have been for single-index arrays (vectors), but everything generalizes natural to multi-dimensional arrays.
Several compound data types contain named fields, including record types (defined with the struct
keyword) and named tuples (mentioned in the previous section). These fields can be accessed with the syntax val.field
. For example, consider again the xyz
tuple from the previous example:
= (x=4, y=5, z=6) # Create a named tuple
xyz # Returns 6 (structure access syntax) xyz.z
Many compound data types are immutable, so that fields cannot be assigned. For those that are mutable, though, we can also use val.field
on the left hand side of an assignment statement.
One particularly convenient feature of Julia is string interpolation. For example, consider the following Julia code:
let
= 99.0/70.0
sqrt2approx "The error in sqrt(2)-$sqrt2approx is $(sqrt(2)-sqrt2approx)"
end
"The error in sqrt(2)-1.4142857142857144 is -7.215191261922271e-5"
In evaluating the string, Julia replaces $sqrt2approx
with the value of the variable sqrt2approx
, and it replaces $(sqrt(2)-sqrt2approx)
with the string representation of the number sqrt(2)-sqrt2approx
. For more controlled formatting (e.g. printing only a few digits of the error), the standard recommendation is to use the @sprintf
macro from the Printf
package.
In addition to building strings by interpolation, Julia lets us build strings by concatentation. The operator *
concatenates two strings; for example "Hello " * "world"
produces "Hello world"
. The ^
operation on strings can also be used for repetition; for example "Hello "^2
produces "Hello Hello "
.
We have already discussed one form of control flow, with short-circuit evaluation of logical and conditional operators. Beyond this, we have conditionals and loops. We treat exception handling in Section 2.6.
Julia has an if
/elseif
/else
/end
statement structure. The statement is actually an expression, which evaluates to whatever branch is taking. For example
is_leap_year(year) =
if year % 400 == 0
true # Multiples of 400 years are leap years
elseif year % 100 == 0
false # Multiples of 100 years otherwise are not
elseif year % 4 == 0
true # Which is an exception to "every fourth year"
else
false # Otherwise, not a leap year
end
If none of the conditions are satisfied and there is no else
clause, the result evaluates to nothing
.
Unlike loops, if
statements do not introduce a local scope.
A while
loop executes while the loop condition is true
function gcd(a,b)
while b != 0
= b, a%b
a, b end
aend
A for
loop iterates over a collection. Perhaps most frequently, this is an integer range:
function mydot(x, y)
@assert length(x) == length(y)
= 0.0
result for i = 1:length(x)
+= x[i]*conj(y[i])
result end
resultend
However, we can also iterate over other collections, whether they are concrete or abstract. For example, the zip
function creates an abstract collection of tuples constructed from elements of parallel lists:
function mydot2(x, y)
@assert length(x) == length(y)
= 0.0
result for (xi, yi) in zip(x, y)
+= xi*conj(yi)
result end
resultend
In the for
loop syntax, we can use =
, in
, or ∈
between the loop variable name and the collection to be iterated over.
The break
statement jumps out of a loop in progress.
Both while
loops and for
loops produce “soft” local scope: references to already-used names from surrounding scopes will not be shadowed, but any new names will be bound in the local scope.
Julia provides three primary syntactic variants for defining functions:
function
keyword->
operatorThe function
keyword is the standard choice for defining complicated functions that might involve several statements. We have seen a few examples already, but for now let’s consider again the Euclidean algorithm for the greatest common denominator, seen in Section 2.3.2:
function gcd(a, b)
while b != 0
= b, a%b
a, b end
aend
The defined function gcd
takes two positional arguments (a
and b
). On a function call, the concrete arguments to the call are bound to the names a
and b
in a local scope, and then the body of the function is then evaluated to obtain a return value. We note that as with other compound statements in Julia (e.g. begin/end
blocks and let/end
blocks), the value of the function body is the last expression in the body – in this case, the final value of a
. We could also write return a
, but this is not strictly necessary in this case. In general, return
statements are only really needed when we want to exit before reaching the end of a function body.
For short expressions, it is more concise to define functions using Julia’s equational syntax:
N0(x) = (x-1.0)*(x-2.0)*(x-3.0)/6.0
The function body on the right-hand side of the function definition can be any valid Julia expression, including a compound statement inside a begin
/end
constrution. Conventionally, though, we use simple expressions for the function bodies when using this syntax.
Our examples thus far have all involved named functions, but Julia also supports anonymous functions (sometimes called lambdas). We can either write anonymous functions with the function
keyword, or using the ->
syntax:
= function(x) x^2 + 2x + 1 end
f1 = x -> x^2 + 2x + 1 f2
The body in the function
version is a compound expression ending with end
. The body in the ->
version is a single expression, similar to what we saw with the equational syntax.
Most Julia functions take arguments. We can define arguments later in the list to have default values; for example:
increment(x, amount=1) = x + amount
Default argument values may depend on earlier arguments.
Julia functions can take keyword arguments as well as positional arguments. Keyword arguments must have a default value. In the function definition, a semicolon separates the positional argument list from the keyword argument list. For example, consider the following routine to print a list in sorted order:
function print_sorted(l; by=identity)
= sort(l, by=by)
l for v in l
println(v)
end
end
The keyword argument by
indicates a function that returns the sort key, and is passed through to the call to the sort
function in the first line. We need the semicolon separator in the function definition to distinguish keyword arguments from positional arguments with defaults. However, there is no such requirement when calling the function, so the call to sort
only uses commas to separate arguments.
When a Julia function is created inside an enclosing local scope, it may refer to variables that were created inside that scope. In this case, Julia creates a closure of the function definition together with variable bindings from the surrounding environment.
For example, suppose we have a naive Newton iteration for solving a 1D equation \(f(x) = 0\):
""" naive_newton(f, df, x; maxiter=10, ftol=1e-8)
Run Newton iteration to find a zero of `f`
(with derivative `df`) from the starting
point `x`. Returns when the magnitude of `f`
is less than `ftol`, otherwise errors out
after running for `maxiter` steps.
"""
function naive_newton(f, df, x; maxiter=10, ftol=1e-8)
= f(x)
fx for _ = 1:maxiter
if abs(fx) < ftol
return x
end
-= fx/df(x)
x = f(x)
fx end
if abs(fx) < ftol
return x
end
error("Did not converge in $maxiter steps")
end
We would like to evaluate the list of function values in order to plot the convergence. We could do this by modifying the naive_newton
routine, or we could record each call to the function.
function plot_naive_newton(f, df, x; maxiter=10, ftol=1e-8)
# Create a list of values and a closure that
# records function evaluations to the list
= []
fs function frecord(x)
= f(x)
fx push!(fs, fx)
fxend
= naive_newton(frecord, df, x,
x =maxiter, ftol=ftol)
maxiterplot(filter(fx->fx > 0, abs.(fs)),
="k",
xlabel=:log10, label="\$f(x_k)\$")
yscaleend
Now a call to plot_naive_newton(sin, cos, 0.5)
(for example) will show a plot of the function value magnitudes at each step of the iteration.
plot_naive_newton(sin, cos, 0.5)
Julia allows a single function name to provide a common interface to different implementations (methods) based on the number and type of the arguments. The process of looking up the appropriate method from the invocation can take into account the types of more than one of the arguments (multiple dispatch). For example, consider the definitions
myrotate(z :: Complex, θ) = exp(1im*θ)*z
myrotate(xy :: Vector, θ) =
let c = cos(θ), s = sin(θ), x = xy[1], y = xy[2]
*x-s*y; s*x+c*y]
[cend
Both versions of myrotate
correspond to a rotation of a plane, but we have different implementations depending on whether the first argument is a complex number or a vector.
println(myrotate(1.0 + 0.0im, π/4))
println(myrotate([1.0; 0.0], π/4))
0.7071067811865476 + 0.7071067811865475im
[0.7071067811865476, 0.7071067811865475]
Methods are also very useful for defining structurally recursive functions to process data structures. For example, consider “flattening” a structure of nested tuples into a vector:
function flatten_tuple(t)
= []
result process(x :: Tuple) = process.(x)
process(x) = push!(result, x)
process(t)
resultend
let
= ((1,2), (3,(4,5)), 6)
t flatten_tuple(t)
end
6-element Vector{Any}:
1
2
3
4
5
6
We will have more to say about Julia’s type system in Section 2.5.
As mentioned previously, operators in Julia are just a special type of function. We can add methods to these functions in order to implement operators acting on new types. For example, consider a new type Polynomial
representing polynomials expressed in the monomial basis (we will return to this in Section 2.5):
struct Polynomial
:: Vector{Float64}
coeffs end
We would like to be able to add, subtract, and multiply polynomials. For example, to add implementation of the +
and *
operators from Base
that works with a Polynomial
, we would write
function Base.:+(p :: Polynomial, q :: Polynomial)
= max(length(p.coeffs), length(q.coeffs))
n = zeros(n)
c 1:length(p.coeffs)] = p.coeffs
c[1:length(q.coeffs)] += q.coeffs
c[Polynomial(c)
end
function Base.:*(p :: Polynomial, q :: Polynomial)
= length(p.coeffs) + length(q.coeffs) - 1
n = zeros(n)
c for (i, pi) in enumerate(p.coeffs)
:length(q.coeffs)+i-1] += pi*q.coeffs
c[iend
Polynomial(c)
end
If we would like to treat constants as a type of polynomial, we need to add a conversion routine and some additional methods for the +
and *
functions.
Base.convert(::Type{Polynomial}, c :: Number) = Polynomial([c])
Base.:+(p :: Polynomial, q) = p + convert(Polynomial, q)
Base.:+(p, q :: Polynomial) = convert(Polynomial, p) + q
Base.:*(p :: Polynomial, q) = p * convert(Polynomial, q)
Base.:*(p, q :: Polynomial) = convert(Polynomial, p) * q
We might also want to be able to evaluate a polynomial, since it represents a function:
function (p :: Polynomial)(x)
= 0*x
px for c in reverse(p.coeffs)
= px*x+c
px end
pxend
With these definitions in place, we can write code like
let
= Polynomial([0; 1]) # x
x = x*(2*x+1)
p p(2)
end
10.0
The ambitious reader might also try writing subtraction, polynomial division and remainders.
A higher-order function is a function that takes another function as an argument. Higher-order functions like map
, filter
, and reduce
are the bread-and-butter of functional programming languages, and Julia supports these as well:
map(x->x+1, [1; 2; 3]) # Yields [2; 3; 4]
filter(x->x > 0, [1; -1; 2; 3]) # Yields [1; 2; 3]
reduce(+, [1; 2; 3]) # Yields 6
Higher-order functions are quite common in numerical methods, as many of the operations that we want to perform are naturally expressed in terms of operations on functions (computing integrals and derivatives, optimizing, finding zeros). Indeed, we already saw a simple example of a higher-order function with the naive_newton
example in Section 2.4.3.
A common pattern in higher-order functions is to have a function that takes one function argument and potentially a few other operations. When the argument function is small, it is convenient to use the ->
operation to define a short anonymous function, as we did in Section 2.4.6. But sometimes when the argument function is not small, it gets clunky to write it in place as an anonymous function. For example, suppose we wanted a verbose version of the filter
example above, one that prints out a description of what is going on at each step. Using the notation we have written so far, we could write something like
filter(function(x)
= x > 0
ispositive = ispositive ? "Keep it" : "Toss it"
action println("Testing x = $x: $action")
ispositiveend, [1; -1; 2; 3])
Testing x = 1: Keep it
Testing x = -1: Toss it
Testing x = 2: Keep it
Testing x = 3: Keep it
3-element Vector{Int64}:
1
2
3
Because this pattern occurs so often in Julia, however, there is a special syntax that makes it easier to write:
filter([1; -1; 2; 3]) do x
= x > 0
ispositive = ispositive ? "Keep it" : "Toss it"
action println("Testing x = $x: $action")
ispositiveend
Testing x = 1: Keep it
Testing x = -1: Toss it
Testing x = 2: Keep it
Testing x = 3: Keep it
3-element Vector{Int64}:
1
2
3
Behind the scenes, Julia rewrites constructions of the form
f(fargs) do args
bodyend
to
f(function(args) body end, fargs)
Beyond being used with functions like map
and reduce
, the do
syntax is often used in Julia for tasks like I/O that involve a required setup and teardown phase. For example, writing to a file in Julia is often done using syntax like the following:
open("foo.txt", "w") do f
println(f, "Hello, world")
end
The version of the open
method that takes a function argument in the terse syntax will call the function with the file handle, then close the file afterward.
In addition to nested calls (e.g. f(g(x))
), Julia supports two special syntax constructs for function composition:
∘
) acts like composition in mathematics.|>
) chains together inputs and outputs like the pipe symbol in a Unix shell command or pipes in R.For example:
sqr(x) = x.^2
norm2a(v) = sqrt(sum(sqr(v)))
= sqrt ∘ sum ∘ sqr
norm2b norm2c(v) = v |> sqr |> sum |> sqrt
Generators and comprehensions are a final type of specialized use of anonymous functions. A generator expression has the form expr for var in collection
, with the option of additional comma-separated var in collection
expressions. The expression evaluates to an iterable collection that is lazily evaluated. For example, consider the generator expression
= ((i,j) for i=1:2, j=1:2) g
We can loop over the generator the same as we would over any other collection; with each iteration, the expression (i,j)
is evaluated with a different binding of i
and j
from the iteration space, visited in column-major order (early indices change fastest). For example, consider the output of the following code fragment:
let
= ((i,j) for i=1:2, j=1:2)
g for v in g
println(v)
end
end
(1, 1)
(2, 1)
(1, 2)
(2, 2)
Iterable expressions can be used in many places where Julia would accept a list or vector. For example:
sum(i^2 for i = 1:10)
385
An array comprehension is similar to a generator, but with surrounding square brackets. Rather than producing a lazily-evaluated collection, it produces a concrete array of results. For example, if \(k(x,y)\) is a bivariate function of interest, we can use a generator to concisely produce a matrix of function samples over a regular grid:
= range(0.0, 1.0, length=20)
xgrid = [k(x,y) for x in xgrid, y in xgrid] kxx
At a certain level, computers only manipulate sequences of bits. But programmers usually work at a more abstract level, where those bit sequences are interpreted as different types of values: these 64 bits represent an integer, those 64 bits represent a floating point number, that other set of 64 bits represents the string "program"
(with a null terminator). One of the classic distinctions drawn among programming languages is how types are treated. In statically-typed languages like C++ and Fortran, every variable and expression is associated with a particular type that can be determined at compile time1. In dynamically-typed languages like MATLAB, Python, and R, values have types but variables are not restricted in the types of values to which they may refer. But, of course, things are not so simple — C++ is mostly statically typed, but run-time polymorphism allows some elements of dynamic typing; while modern Python is mostly dynamically typed, but allows optional type annotations.
Julia is a dynamically-typed language, in which variables are just names that can be bound to values of different types. At the same time, Julia has a rich type system with optional type declarations using the ::
syntax as described in Section 2.2.1.2. These type declarations serve three purposes:
Correctness: Sometimes an operation only makes sense with certain types of inputs, and declaring that explicitly makes it possible for the system to provide an early error rather than an error at runtime.
Expressivity: Declaring the types of arguments to functions allows the same function name to dispatch to multiple methods, as discussed in Section 2.4.4.
Performance: When the Julia system can infer types in a function call, the just-in-time compiler can create faster specialized versions of the code, much like what we would find in C or Fortran.
While type declarations is often helpful, Julia also benefits from not requiring them. As an example, consider the gcd
function from Section 2.3.2. We might argue that it is appropriate to restrict the arguments to be Integer
types, i.e.
function gcd(a :: Integer, b :: Integer)
while b != 0
= b, a%b
a, b end
aend
However, the Euclidean algorithm also makes sense in more general Euclidean domains, such as polynomials. So we might extend the operations on the Polynomial
type from Section 2.4.5 so that comparisons with zero and the remainder operator are both allowed. With this extension, the original gcd
function (without the :: Integer
annotations) would produce the correct result.
Julia types form a tree-shaped hierarchy, where each abstract type may have multiple subtypes. For example, Integer
is a subtype of Number
, and concrete integer types like Int64
are subtypes of Integer
. The syntax <:
is used to denote a subtype relation; for example, we would write Integer <: Number
to test that Integer
is indeed a subtype of Number
. At the root of the type hierarchy is the type Any
.
Only abstract types may serve as supertypes in Julia. Concrete types, whether primitive types (like Logical
, UInt16
, or Float64
), or various compound types, are always leaves of the hierarchy.
The primitive concrete numeric types in Julia are
Floating point numbers Float16
, Float32
, and Float64
. These are all subtypes of AbstractFloat
. The 32-bit and 64-bit formats are supported in hardware, but operations with the 16-bit format are software only. We will have much more to say about floating point in Chapter 9.
Signed integer types Int8
, Int16
, Int32
, Int64
, and Int128
. These are all subtypes of Signed
, which is a subtype of Integer
. The type Int
is an alias for Int64
or Int32
, depending on whether Julia is being run on a 64-bit or a 32-bit system.2
Unsigned integer types UInt8
, UInt16
, UInt32
, UInt64
, and UInt128
. These are all subtypes of Unsigned
, which is a subtype of Integer
.
Technically, type Logical
is also a subtype of type Integer
, and it can be used in arithmetic expressions like an integer – for example, true + true
evaluates to 2
, while 1.0 * false
evaluates to 0.0
. However, we do not recommend using logical variables in this way.
In addition to the primitive types, Julia supports parameterized types for rational and complex numbers. For example, 3 // 4
produces a representation of \(3/4\) of type Rational{Int64}
, while 1.0 + 1.0im
produces a representation of \(1+1i\) of type Complex{Float64}
.
Julia also supports arbitrary-precision BigInt
and BigFloat
types, though we will rarely make use of these.
Strings in Julia inherit from the AbstractString
type, and are logically sequences of Unicode characters3. Unicode characters (or “code points”) in turn inherit from the AbstractChar
type, and are associated with 32-bit unsigned integers. The character A
, for example, corresponds to code point U+0041; this is character 65 in decimal, but Unicode codepoints are typically written in hexadecimal. More generally, code points U+0000 through U+007F correspond to the classic ASCII character set, which includes Latin characters and punctuation used in standard English.
The default character type Char
in Julia takes 32 bits in memory. The default String
type is not an array of 32-bit Char
entities; instead, Julia uses UTF-8 strings, which use a variable-length encoding. Unicode codepoint numbers corresponding to ASCII are encoded with one byte per character, and other characters may take more than one byte to represent. Hence, the length
of a string (in characters) and the sizeof
a string (in bytes) mean different things in Julia. These are also different complexity: computing the length (in codeunits) in a UTF-8 string requires that we scan the string from start to end. To avoid the time cost of a scan, indexing into a standard string is done by bytes. For example, in the string s = "sπot"
, we can safely refer to s[1]
(s), s[2]
(π), s[4]
(o), and s[5]
(t); but accessing s[3]
will cause a runtime exception.
Symbols inherit from type Symbol
. They are encoded as “interned” UTF-8 strings, which means that the string is stored in a hash table, and the data passed around the Julia environment is the location in that hash table.
The type of a tuple is a Cartesian product of the tuple elements. In Julia, this is represented as Tuple
, a parameterized type whose parameters are the component types. For example, the type of (1.0, 2)
is Tuple{Float64,Int64}
. Named tuples are similar; the type of (x=1.0, y=2)
is @NamedTuple{x :: Float64, y :: Int64}
.
We declare a structure type (sometimes called a record type) with the struct
keyword; for example, we might define a type for rational numbers as
# Version 0
struct MyRational <: Number
p
qend
Then MyRational(1,2)
will give us an instance of MyRational
with p
set to 1 and q
set to 2. However, this definition also allows us to legitimately write MyRational("bad", "idea")
, and we probably do not want to allow a ratio of two strings. A second attempt at a definition might be
# Version 1
struct MyRational <: Number
:: Integer
p :: Integer
q end
With this definition, we are restricted to integer values for p
and q
. However, we may still be unhappy with this if we read about Julia performance, and see the advice to make fields be concrete types to produce efficient code (Integer
is an abstract type). This brings us to our third attempt, which involves making MyRational
a parametric type:
struct MyRational{T<:Integer} <: Number
:: T
p :: T
q end
Writing {T<:Integer}
after the structure name says that we want to have a type parameter T
which is a subtype of Integer
. With this definition, the constructor MyRational(1,2)
will give us an instance of type MyRational{Int64}
(assuming a 64-bit machine).
We may decide that we would like to make p
and q
have reduced form, and disallow a zero denominator. To do this, we would create a new inner constructor function:
struct MyRational{T<:Integer} <: Number
:: T
p :: T
q MyRational(p :: Integer, q :: Integer) =
if q == zero(q)
error("Divide by zero")
else
= promote(p, q)
p, q = gcd(p, q)
r new{typeof(p)}(p÷r, q÷r)
end
end
There are several things to unpack here:
The MyRational
takes two abstract Integers
. These could be different types of integers! In order to gracefully deal with this, we promote the arguments to a common type with the promote
command (see Section 2.5.11)
Inside the MyRational
constructor, we cannot just call another constructor called MyRational
. Instead, we call new
to create the object. This is only used inside of inner constructors.
The new
call needs a type parameter, which we get from the type of p
.
In addition to inner constructors, we can also define an outer constructor, so named because the definition is outside the struct
statement. For example, we might want to define a constructor that takes just one argument (with the second being an implicit 1):
MyRational(p :: Integer) = MyRational(p, one(p))
As described earlier, we might also decide that we wanted to overload the operators. Note that this does not require any special magic with the type parameters to MyRational
:
Base.:+(a :: MyRational, b :: MyRational) =
MyRational(a.p*b.q + a.q*b.p, a.q*b.q)
Base.:-(a :: MyRational, b :: MyRational) =
MyRational(a.p*b.q - a.q*b.p, a.q*b.q)
Base.:*(a :: MyRational, b :: MyRational) =
MyRational(a.p*b.p, a.q*b.q)
Base.:/(a :: MyRational, b :: MyRational) =
MyRational(a.p*b.q, a.q*b.p)
What if we want to get the type associated with p
and q
? We can, of course, use typeof(p)
to do this, but we can also write a method to extract the type. Unlike the overloaded operators, this method does explicitly refer to the type parameter, and so we need some additional syntax:
inttype(:: MyRational{T}) where {T} = T
Unpacking this construction, we have
An anonymous argument of type MyRational{T}
. We do not need to give the argument a name, because we do not care about the value of the variable; we only care about its type.
The where {T}
clause to specify that this method definition is parameterized by the type T
.
By default, the entries of a struct are immutable. Using the example from the previous section, once we have created a MyRational
with a given p
and q
, we cannot change those values. Note that if any field of a structure is a writeable container (like an array), we can write into that container; the container is immutable, but its contents are not. Still, sometimes we would genuinely like to be able to update a structure, and for this we have mutable structures, which are declared with the mutable
keyword. For example, consider the following counter struct:
mutable struct FooCounter
:: Int
p FooCounter() = new(0)
end
Now we can write an increment operator that updates this structure (using an exclamation mark in the name, as is the convention for functions that change their inputs):
increment!(f :: FooCounter) = (f.p += 1)
= FooCounter() # foo.p initialized to 0
foo println(increment!(foo)) # Prints 1
println(increment!(foo)) # Prints 2
An Array{T,d}
is a \(d\)-index array of type T
. The type aliases Vector{T}
and Matrix{T}
correspond to arrays with d
set to 1 and 2, respectively. It is possible to make T
an abstract type, but it is generally much more efficient to make T
a concrete type.
In addition to arrays, there are several constructs that are subtypes of the AbstractArray
interface. These includes:
Range
: An abstraction for a regularly-spaced range, which can be constructed by the colon operator (e.g. 1:3
) or with the range
command.
Adjoint
: (Conjugate) transpose of a matrix.
BitArray
: A compact storage format for an array of logical values.
SparseMatrixCSC
: A sparse matrix format.
SubArray
: A submatrix of an existing matrix, which can be constructied using the view
macros.
StaticArray
: An fast type for small arrays whose dimensions are part of the type.
These constructs can all be indexed, sliced, iterated over, and otherwise treated much like ordinary arrays (though not all are writeable).
In addition to tuples and array-like collections, Julia supports dictionary and set types. A dictionary contains key-value pairs, written as key => value
in Julia; the =>
operator constructs a Pair
object. The structure supports addition and deletion of pairs, testing for keys, iteration, and lookup (either with get
, which supports default values, or with dict[key]
, which does not). There are a few variants, the most frequently used of which is the ordinary dictionary (Dict
).
Julia also supports a Set
type, which supports addition and removal, testing for inclusion or subset relations, and operations like union, intersection, and set differences.
Julia allows Union
types whose values may be one of a selection of types. Two specific examples have to do with data that may be “null” or missing4. In Julia, the type Nothing
(with instance nothing
) is used the way that programmers use “null” in some other languages, while the type Missing
(with instance missing
) is used to indicate missing data. These types are most often used in the context of Union{T, Nothing}
or Union{T, Missing}
types. For example, values of type Union{T, Missing}
can either be something of type T or the value missing
.
Types are themselves things that can be represented and manipulated in Julia. Each type has type DataType
, which is a subtype of the parametric abstract type Type{T}
. In general, a type T
(and no other type) is an instance of Type{T}
. This can be useful for certain special types of method dispatch situations. Similarly, value types (Val{T}
) are parameterized types that take a value as a parameter, and can be used for specialized method dispatch or for helping the compiler with performant code.
The convert
function in Julia is used to convert values of one type to another type (when possible). The system calls convert
automatically when assigning a value of one type to a variable or field of another type. The first argument is the target type, and the second is the value to be converted. Defining new methods for convert
is a good use of the Type
type; for example, we might define a conversion from an Integer
to the MyRational
type defined earlier as:
Base.convert(::Type{MyRational}, p :: Integer) = MyRational(p)
A more elaborate example might involve converting MyRational
values with different integer parameter types:
Base.convert(::Type{MyRational{T}}, r :: MyRational) where {T} =
MyRational(T(r.p), T(r.q))
The promote
function in Julia takes a list of values and converts them both to a compatible greatest type. We typically would not write new methods for promote
; instead, we write new methods for promote_rule
. Again using MyRational
as an example, we might have
Base.promote_rule(::Type{MyRational{T}}, ::Type{S}) where {T<:Integer,S<:Integer} =
promote_type(T,S)} MyRational{
Julia has throw/catch system for handling exceptional conditions, similar to many other modern languages (Python, MATLAB, R, Java, C++). Information about exceptions is encoded in via a struct that is a subtype of Exception
, which we throw
when the error is detected:
""" imlog(x)
Compute the imaginary part of the principal branch
of the log of negative number x.
"""
imlog(x) =
>= zero(x)) ?
(x throw(DomainError("Argument must be negative")) :
log(-x)
The error
function throws an instance of an ErrorException
.
Exceptions that are thrown at some point in a call chain can be caught by callers using a try
/catch
statement.
""" generous_log(x :: Number)
Computes the log of a number; if a domain error is
thrown for the log of a negative real number, tries
with a complex version of the number.
"""
generous_log(x :: Number) =
trylog(x)
e
catch if e isa DomainError
log(Complex(x))
else
throw(e)
end
end
These statements may also have an optional else
clause that is executed when the try
statement succeeded without exceptions; and an optional finally
clause that is executed when exiting the try
/catch
, independent of how the exit occurred.
Functions, types, and modules can all have documentation strings just before them, which are referenced by the Julia help system. Documentation strings are formatted using Markdown. A fenced code block beginning with jldoctest
can be used to produce a usage example that is suitable for automatic testing.
We have already discussed some of Julia’s array support. However, one of the attractive features of Julia is the support not only for arrays, but for numerical linear algebra computations. Many of these facilities are in the LinearAlgebra
package, which we will almost always use throughout this book.
By default, we think of a one-dimensional array as a column vector, and a two-dimensional array as a matrix. We can do standard linear algebra operations like scaling (2*A
), summing like types of objects (v1+v2
), and matrix multiplication (A*v
). The expression
= v' w
represents the adjoint of the vector v
with respect to the standard inner product (i.e. the conjugate transpose). The tick operator also gives the (conjugate) transpose of a matrix. We note that the tick operator in Julia does not actually copy any storage; it just gives us a re-interpretation of the argument. This shows up, for example, if we write
let
= [1, 2] # v is a 2-element Vector{Int64} containing [1, 2]
v = v' # w is a 1-2 adjoint(::Vector{Int64}) with eltype Int64
w 2] = 3 # Now v contains [1, 3] and w is the adjoint [1, 3]'
v[end
3
Julia gives us several standard matrix and vector construction functions.
= zeros(n) # Length n vector of zeros
Z = zeros(n,n) # n-by-n matrix of zeros
Z = rand(n) # Length n random vector of U[0,1] entries (from Random)
b e = ones(n) # Length n vector of ones
= diagm(e) # Construct a diagonal matrix
D = diag(D) # Extract a matrix diagonal e2
The identity matrix in Julia is simply I
. This is an abstract matrix with a size that can usually be inferred from context. In the rare cases when you need a concrete instantiation of an identity matrix, you can use Matrix(I, n, n)
.
In addition to functions for constructing specific types of matrices and vectors, Julia lets us put together matrices and vectors by horizontal and vertical concatenation. This works with matrices just as well as with vectors! Spaces are used for horizontal concatenation and semicolons for vertical concatenation.
= [1; 2] # Length-2 vector
y = [1 2] # 1-by-2 matrix
y = [1 2; 3 4] # 2-by-2 matrix
M = [I A] # Horizontal matrix concatenation
M = [I; A] # Vertical matrix concatenation M
Julia uses commas to separate elements of a list-like data type or an array. So [1, 2]
and [1; 2]
give us the same thing (a length 2 vector), but [I, A]
gives us a list consisting of a uniform scaling object and a matrix — not quite the same as horizontal matrix concatenation.
Julia lets us rearrange the data inside a matrix or vector in a variety of ways. In addition to the usual transposition operation, we can also do “reshape” operations that let us interpret the same data layout in computer memory in different ways.
# Reshape A to a vector, then back to a matrix
# Note: Julia is column-major
= reshape(A, prod(size(A)));
avec = reshape(avec, n, n)
A
= randperm(n) # Random permutation of indices (need to use Random)
idx = A[:,idx] # Permute columns of A
Ac = A[idx,:] # Permute rows of A
Ar = A[idx,idx] # Permute rows and columns Ap
Julia lets us extract specific parts of a matrix, like the diagonal entries or the upper or lower triangle. Some operations make separate copies of the data referenced:
= randn(6,6) # 6-by-6 random matrix
A 1:3,1:3] # Leading 3-by-3 submatrix
A[1:2:end,:] # Rows 1, 3, 5
A[:,3:end] # Columns 3-6
A[
= diag(A) # Diagonal of A (as vector)
Ad = diag(A,1) # First superdiagonal
A1 = triu(A) # Upper triangle
Au = tril(A) # Lower triangle Al
Other operations give a view of the matrix without making a copy of the contents, which can be much faster:
= randn(6,6) # 6-by-6 random matrix
A view(A,1:3,1:3) # View of leading 3-by-3 submatrix
view(A,:,3:end) # View of columns 3-6
= UpperTriangular(A) # View of upper triangle
Au = LowerTriangular(A) # View of lower triangle Al
Julia provides a variety of elementwise operations as well as linear algebraic operations. To distinguish elementwise multiplication or division from matrix multiplication and linear solves or least squares, we put a dot in front of the elementwise operations.
= d.*x # Elementwise multiplication of vectors/matrices
y = x./d # Elementwise division
y = x + y # Add vectors/matrices
z = x .+ 1 # Add scalar to every element of a vector/matrix
z
= A*x # Matrix times vector
y = x'*A # Vector times matrix
y = A*B # Matrix times matrix
C
# Don't use inv!
= A\b # Solve Ax = b *or* least squares
x = b/A # Solve yA = b or least squares y
There are few good reasons to compute explicit matrix inverses or determinants in numerical computations. Julia does provide these operations. But if you find yourself typing inv
or det
in Julia, think long and hard. Is there an alternate formulation? Could you use the forward slash or backslash operations for solving a linear system?
Beyond the LinearAlgebra
package mentioned in Section 2.8, there are several other packages that we will use or recommend in this book:
The Statistics
package provides a variety of standard statistics computations (mean, median, etc). There is also a StatsBase
package that provides some additional statistics. There has been discussion of refactoring these.
The StatsFuns
package provides a variety of special functions that occur in statistical computations (e.g. the log gamma function or various common pdf, cdf, and inverse cdf functions). The Distributions
package is recommended as providing a more convenient interface, but sometimes we will want the low-level interface.
The SparseArrays
package is generally useful for setting up sparse matrix and vectors. In addition, Julia provides packages for a variety of sparse direct and iterative solvers, including the SuiteSparse
package (direct solvers) and the IterativeSolvers
package (iterative solvers).
The DataFrames
package provides functionality for manipulating tabular data similar to the Python Pandas
package.
Though we will only make limited use of these within the book, there are several Julia packages that support automatic differentiation, including ForwardDiff
, ReverseDiff
, Zygote
, and Enzyme
.
The Plots
package provides a standard interface for producing plots with a variety of plotting backends. There are other plotting packages in Julia (the second-most popular of which is probably Makie
), but the Plots
package makes a good starting point.
For text output of tabular data, we recommend the PrettyTable
package. This can also write to a variety of output formats, including HTML, Markdown, and LaTeX.
The FileIO
package provides a common interface for reading (and writing) a wide variety of data file types, including text, tabular, image, sound, and video formats.
The Test
package in Julia provides a variety of useful tools for writing unit tests. When tests fail and it is necessary to step through code, the Debugger
package is useful. For performance evaluation, we recommend the BenchmarkTools
and Profile
packages.
Julia provides a LISP-like metaprogramming facility that allows us to manipulate Julia code as data. This is a very powerful tool, allowing us to (for example) write macros, i.e. code that generates new code on our behalf. Metaprogramming requires care to get right, and the usual advice is to only use it when other tools do not suffice. However, there are some times when it really is the right tool, and so we will touch on it here at least lightly.
The heart of metaprogramming is the understanding that code is a type of data, and can be manipulated as such. We distinguish code that we want to treat as data by quoting it, which involves bracketing the code by :(
and )
or by quote
and end
. For example, the expression
= :(1 + 1) quoted_sum
:(1 + 1)
assigns quoted_sum
to the unevaluated code for 1 + 1
. To evaluate the code, we can call the eval
function
eval(quoted_sum)
2
Quoted expressions are printed by the Julia interpreter in natural Julia syntax. Internally, though, these expressions are not represented by strings, but by expression trees. We can see the structure of the expresion tree with the dump
command, e.g.
dump(quoted_sum)
Expr
head: Symbol call
args: Array{Any}((3,))
1: Symbol +
2: Int64 1
3: Int64 1
The expression tree data structure involves internal nodes of type expr
that have a head
field indicating the type of operation represented, and an args
array indicating the arguments of the operation. The leaves in the expression tree are literals, including symbols.
Just as we can use string interpolation to build strings that include computations, we can also use interpolation with quoted code; for example
let
= 1
x :($x + 1)
end
:(1 + 1)
We can produce the same expression tree without quoting by explicitly using the constructor for Expr
objects, e.g.
let
= 1
x Expr(:call, :+, x, 1)
end
:(1 + 1)
For the reader who really wants to learn to write macros, we recommend the documentation together with the book On Lisp (Graham 1993). We will write small macros a few times in order to simplify our life; and in any case, we also use macros written by others often enough that we are obliged to say a few words from a user’s perspective.
A macro maps a tuple of input expressions to an output expression, which is then compiled. The actual work of the macro code occurs at compile time; at run time, the system evaluates whatever output expression was generated by the macro. Macros are frequently used for various types of programming utilities. Debugging, timing, assertions, and formatted output in Julia are all conventionally done with macros.
The argument expressions can be parenthesized or not; for example, the following two lines are functionally identical invocations of the @assert
macro:
@assert(f(), "Call to f failed")
@assert f() "Call to f failed"
As a toy example, of a macro definition, consider a macro that prints the input expression
macro printarg(e)
= "$e = "
s quote
let v = $e
println($s, v)
vend
end
end
@printarg (macro with 1 method)
When calling a macro, we preface the name with an @
sign; for example
@printarg 1+1
1 + 1 = 2
2
This is essentially the functionality of the built-in @show
macro, except that our printarg
macro fails when the argument involves an expression with a variable. For example, @show
works fine in this context, while @printarg
would fail:
let
= 1
x @show x + x
end
x + x = 2
2
The problem has to do with the fact that Julia macros and rename variables internally when there is a possibility of accidental conflict with names in an outer environment in a local scope. To deliberately refer to a variable in the outer environment, that variable must be escaped. We can see the difference by expanding both the @printarg
and the @show
macros on the same input using the @macroexpand
macro:
@macroexpand @printarg x + x
quote #= In[40]:4 =# let var"#97#v" = Main.x + Main.x #= In[40]:5 =# Main.println("x + x = ", var"#97#v") #= In[40]:6 =# var"#97#v" end end
We can see that all the variable names in this expression are explicitly scoped to the Main
module. In contrast, for @show
, we have
@macroexpand @show x + x
quote
Base.println("x + x = ", Base.repr(begin
#= show.jl:1232 =#
local var"#99#value" = x + x
end))
var"#99#value"
end
Here we see that a local variable (with a name that will not conflict with any other name in the system) is assigned to the expression x+x
.
As an example of a more significant small utility macro, consider a closure defined inside of an outer function. Such closures sometimes suffer5 from a performance issue because the compiler analysis does not determine that variables from the outer context (captured variables) are “type stable” (see Chapter 3). Putting the function into a let
block can help address the issue, e.g. replacing {.julia} f = x -> C*x
with {.julia} f = let C=C x -> C*x end
This does make a semantic change to the program; for example, the code {.julia} C = 10 f = x -> C*X C = 20 f(5)
will produce 100, while introducing an enclosing let
block would produce 50. Nonetheless, in many cases captured variables remain unchanging for the useful lifetime of a closure, and this is a valid transformation.
Let us assume we are interested in automatically producing a surrounding let
block. In this case, we first need a list of variables referenced in an expression. We do this by recursing through the call tree and adding any symbols we encounter (except for function names) to a set. Note that we can make this code very concise by defining different methods associated with the type of the first argument (dispatching based on an Expr
, a Symbol
, or another type)
symbols!(e, s = Set{Symbol}()) = nothing
symbols!(e :: Symbol, s = Set{Symbol}()) = push!(s, e)
function symbols!(e :: Expr, s = Set{Symbol}())
= e.head == :call ? e.args[2:end] : e.args
args for arg in args
symbols!(arg, s)
end
send
Now we are in a position to define the macro:
macro letclosure(e)
# Check that this is a function definition with ->
@assert e.head == :(->)
# Get the symbols that are not in the argument list
= Set{Symbol}()
s = Set{Symbol}()
sarg symbols!(e.args[2], s)
symbols!(e.args[1], sarg)
setdiff!(s, sarg)
# Bind local variables to the outer version
= [quote $v = $(esc(v)) end for v in s]
bindings
# Put the expression into the let block
quote
let
$(bindings...)
$e
end
end
end
There are a few points to note about our definition:
For local variables, the macro
system produces “local” versions of names that do not coincide with any other names in the module. When we want such a coincidence, we need to explicitly escape a symbol with the esc
function.
As in other contexts, the “splatting” operator ($(bindings...)
) expands the contents of a list in place.
When we evaluate the macro, the code is transformed and the result is inserted into our program; for example:
let
= 10
C = @letclosure x -> C*x
f = 20
C f(5)
end
50
If we want to see how a macro expannds out to regular code, we can use the @macroexpand
macro. Continuing with our example above, we have:
@macroexpand @letclosure x -> C*x
quote #= In[46]:18 =# let #= In[46]:19 =# begin #= In[46]:14 =# var"#103#C" = C end #= In[46]:20 =# (var"#105#x",)->begin #= In[48]:1 =# var"#103#C" * var"#105#x" end end end
Note in particular that in the variable bindings, the regular C
appears on the right hand side, but a local version of C
is assigned to. All the other variables in the generated code likewise have purely local names that are guaranteed not to coincide with other names in the module.
Beyond writing language utilities, Julia macros are useful for writing embedded domain-specific languages (DSLs) for accomplishing particular tasks. In this setting, we are really writing a language interpreter embedded in Julia, using the Julia parse tree as an input and producing Julia code as output.
As an example of both Julia macro design and Julia programming more generally, we will design a small language extension for writing functions that transform code based on structural pattern matching; for example, a rule written as:
0-x => x -> -x
corresponds to a function something like
e ->
# Check that e is a binary minus and the first argument is zero.
if (e isa Expr && e.head == :call &&
e.args[1] == :- && length(e.args) == 3 &&
e.args[2] == 0)
# Bind the name "x" to the second argument
let x = e.args[3]
= :(-$x) # Create a "-x" expression
result true, result # Yes it matched, return result
end
elsefalse, nothing # No, it didn't match
end
That is, we want a function that takes an input
-
operation at the top?x
to the second term in the sum, produce a new expression :(-$x)
, and assign to result
. Return the pair (true, result)
.(false, nothing)
.The one-line description is already more concise for this simple example; if we wanted to match against a more complicated pattern, the Julia code corresponding to a rule in the syntax
=> (argument_symbols) -> result_code pattern
becomes that much more complex. Our key task is therefore take rules written with the first syntax and to convert them automatically to standard functions.
This is a more ambitious example, and can be safely skipped over. However, it is both a useful example of a variety of features in Julia and introduces concrete tools we will use in later chapters (e.g. when discussing automatic differentiation).
The parse tree for Julia code includes not only expressions, but also line numbers, e.g.
:(x->x+1)
:(x->begin
#= In[49]:1 =#
x + 1
end)
In the code shown, the comment line immediately after the begin
statement corresponds to a LineNumberNode
. Such LineNumberNode
nodes are useful for debugging, but are a nuisance for our attempts at pattern matching. Therefore, before doing anything else, we will write a utility to get rid of such nodes. Because LineNumberNode
is a distinct type, we can use Julia’s multiple dispatch as an easy way to accomplish this.
filter_line_numbers(e :: Expr) =
let args = filter(a -> !(a isa LineNumberNode), e.args)
Expr(e.head, filter_line_numbers.(args)...)
end
filter_line_numbers(e) = e
At the heart of our matching algorithm will be a function match_gen
that generates code to check whether an expression matches a pattern. The inputs to the match_gen
function are
In the simplest case, for non-symbol leaf nodes, we declare a match when the pattern agrees with the expression.
match_gen!(bindings, e, pattern) = :($e == $pattern)
For example, the following code matches the expression named expr
to the pattern 0
:
match_gen!(Dict(), :expr, 0)
:(expr == 0)
Things are more complicated when we match a symbol. If the symbol is not in the bindings dictionary, we just check that the expression equals the (quoted) symbol. Otherwise, there are two cases:
expr
matches the previous binding.match_gen!(bindings, e, s :: Symbol) =
if !(s in keys(bindings))
= QuoteNode(s)
qs :($e == $qs)
elseif bindings[s] == nothing
= gensym()
binding = binding
bindings[s] quote $binding = $e; true end
else= bindings[s]
binding :($e == $binding)
end
We illustrate all three cases below.
let
= Dict{Symbol,Any}(:x => nothing)
bindings println( match_gen!(bindings, :expr, :x) )
println( match_gen!(bindings, :expr, :x) )
println( match_gen!(bindings, :expr, :y) )
end
begin
#= In[53]:8 =#
var"##231" = expr
#= In[53]:8 =#
true
end
expr == var"##231"
expr == :y
The most complex case is matching expressions. The basics are not complicated: an item e
matches an Expr
pattern if e
is an expression with the same head as the pattern and if the arguments match.
match_gen!(bindings, e, pattern :: Expr) =
let head = QuoteNode(pattern.head),
= match_gen_args!(bindings, e, pattern.args)
argmatch :( $e isa Expr && $e.head == $head && $argmatch )
end
The building block for checking the arguments will be to check that a list of expressions matches a list of patterns.
match_gen_lists!(bindings, exprs, patterns) =
foldr((x,y) -> Expr(:(&&), x, y),
match_gen!(bindings, e, p)
[e,p) in zip(exprs, patterns)]) for (
The more complicated case is ensuring that the arguments match. This is in part because we want to accomodate the possibility that the last argument in the list is “splatted”; that is, a pattern like f(args...)
should match f(1, 2, 3)
with args
bound to the tuple [1, 2, 3]
. In order to do this, we would first like to make sure that we can sensibly identify a “splatted argument.”
is_splat_arg(bindings, e) =
e isa Expr &&
e.head == :(...) &&
e.args[1] isa Symbol &&
e.args[1] in keys(bindings)
Note that we only consider something a “splatted argument” if the argument to the splat operator is a symbol in the bindings table.
To implement the check, we create a let
statement to bind a local name to each argument. If the last pattern is splatted, we make sure the last term in the tuple on the left-hand side of the statement is splatted (and then remove the pattern splat). Finally, we generate checks to make sure each of the locally-named arguments matches with the associated term in the pattern list.
match_gen_args!(bindings, e, patterns) =
if isempty(patterns)
:(length($e.args) == 0)
else
= length(patterns)
nargs = :(length($e.args) == $nargs)
lencheck = Vector{Any}([gensym() for j = 1:length(patterns)])
args = Expr(:tuple, args...)
argstuple
# De-splat pattern / splat arg assignment and
# adjust the length check
if is_splat_arg(bindings, patterns[end])
= copy(patterns)
patterns end] = patterns[end].args[1]
patterns[end] = Expr(:(...), argstuple.args[end])
argstuple.args[= :(length($e.args) >= $(nargs-1))
lencheck end
= match_gen_lists!(bindings, args, patterns)
argchecks :($lencheck && let $argstuple = $e.args; $argchecks end)
end
We strongly recommend the reader trace through this code for some examples until enlightenment strikes.
Given a list of symbols and a pattern in which they appear, we can produce code to generate a mapping from expressions to either (true,bindings)
where bindings
is a list of all the subexpressions bound to the name list; or (false, nothing)
if there is no match.
function compile_matcher(symbols, pattern)
= Dict{Symbol,Any}(s => nothing for s in symbols)
bindings
# Input expression symbol and matching code
# (symbol bindings are indicated by bindings table)
= gensym()
expr = match_gen!(bindings, expr, pattern)
test
# Get result vals (symbol/nothing) and associated variable names
= [bindings[s] for s in symbols]
result_vals = filter(x -> x != nothing, result_vals)
declarations
# Produce the matching code
= Expr(:tuple, result_vals...)
results :($expr ->
let $(declarations...)
if $test
true, $results)
(else
false, nothing)
(end
end)
end
It is convenient to compile this into a macro. The macro version also filters line numbers out of the input pattern and out of any exprssion we are trying to match.
macro match(symbols, pattern)
@assert(symbols.head == :tuple &&
all(isa.(symbols.args, Symbol)),
"Invalid input symbol list")
= filter_line_numbers(pattern)
pattern = compile_matcher(symbols.args, pattern)
matcher :($matcher ∘ filter_line_numbers)
end
A rule has the form
=> args -> code pattern
where args
is the list of names that are bound in the pattern, and these names can be used in the subsequent code. We want to allow the code to potentially use not only the matched subexpression, but also to refer to the symbol as a whole; we use the named argument feature in tuples for that. So, for example, in the rule
- y => (x,y;e) -> process(e) x
we mean to match any subtraction, bind the operands to x
and y
and the expression as a whole to e
, and call process(e)
to process the expression.
Now that we have a macro for computing matches, we can use it to help with parsing the rule declarations into argument names, expression name (if present), the pattern, and the code to be called on a match.
function parse_rule(rule)
= @match (pattern, args, code) pattern => args -> code
match_rule = match_rule(rule)
isok, rule_parts if !isok
error("Syntax error in rule $rule")
end
= rule_parts
pattern, args, code
println("Pattern: $pattern")
println("Args: $args")
println("Code: $code")
= @match (args, expr) (args..., ; expr)
match_arg1 = @match (args, expr) (args..., )
match_arg2 println("Did match 1: $match_arg1")
println("Did match 2: $match_arg2")
begin println("m1"); ismatch, bindings = match_arg1(args); ismatch end ||
begin println("m2"); ismatch, bindings = match_arg2(args); ismatch end ||
begin bindings = (Vector{Any}([args]), nothing) end
println("Bindings: $bindings")
= bindings
symbols, expr_name
if !all(isa.(symbols, Symbol))
error("Arguments should all be symbols in $symbols")
end
if !(expr_name == nothing || expr_name isa Symbol)
error("Expression parameter should be a symbol")
end
symbols, expr_name, pattern, codeend
Matching a pattern on an input produces a boolean variable (whether there was a match or not) and a table of bindings from names in the pattern to symbols generated during the match. In order to safely access the right-hand-side symbols that were generated during the match, we need to declare them (in a let
statement). If there is a match, we bind them to the ordinary input names (things like :x
) in a second let statement, then call the code in that second let statement and assign the result to the result
symbol. The resulting code must be used in a context where the expr
and result
symbols are already set up.
function compile_rule(rule, expr, result)
= parse_rule(rule)
symbols, expr_name, pattern, code = Dict{Symbol,Any}(s => nothing for s in symbols)
bindings = match_gen!(bindings, expr, pattern)
test
# Get list of match symbols and associated declarations
= [bindings[s] for s in symbols]
result_vals = filter(x -> x != nothing, result_vals)
declarations
# Set up local bindings of argument names in code
= [:($s = $(r == nothing ? (:nothing) : r))
binding_code in zip(symbols, result_vals)]
for (s,r) if expr_name != nothing
push!(binding_code, :($expr_name = $expr))
end
# Produce the rule
= gensym()
ismatch quote
let $(declarations...)
$ismatch = $test
if $ismatch
$result = let $(binding_code...); $code end
end
$ismatch
end
end
end
Now that we are able to compile a rule, we set up a macro for compiling a rule into a standalone function. The input expression symbol (expr
) is the named argument to this standalone function, and at the end the function returns both the match condition (which is what the compiled code evaluates to) and the result (which the compiled code produces as a side effect).
macro rule(r)
= gensym(), gensym()
expr, result = compile_rule(filter_line_numbers(r), expr, result)
code esc(quote
$expr ->
let $result = nothing
$code, $result
end
end)
end
Finally, we often want to combine rules, executing the first one that matches a given input. The @rules
macro takes a set of rules packaged into a block and tries them in sequence, returning the result of whichever rule was executed (or nothing
if no rule was executed).
macro rules(rblock :: Expr)
= filter_line_numbers(rblock)
rblock if rblock.head != :block
error("Rules must be in a begin/end block")
end
# Set up input name, output name, and rule list
= gensym(), gensym()
expr, result = rblock.args
rules
# Or together all the tests
=
rule_calls foldr((x,y) -> Expr(:(||), x, y),
compile_rule(r, expr, result) for r in rules])
[
# Call all the rules, return the computed result
esc(quote
$expr ->
let $result = nothing
$rule_calls
$result
end
end)
end
We conclude our matching example by giving a few examples of the process in practice.
FIXME: Examples mostly commented out while looking for a scoping fix from Julia 1.11 breaking change.
Julia parses chains of additions and multiplications into one long list. For example, 1+2+3
parses to (1+2)+3
rather than to +(+(1,2), 3)
. For some operations, it is simpler to only allow binary addition and multiplication nodes. The following function converts these nodes.
#| output: false
function reassoc_addmul(e)
# Apply folding to an op(args...) node (op = +, *)
fold_args(args, op) =
foldl((x,y) -> Expr(:call, op, x, y), args)
# Rules for processing one expression node
= @rules begin
r +(args...) => args -> fold_args(args, :+)
*(args...) => args -> fold_args(args, :*)
e => e -> e
end
# Recursively process the graph
process(e :: Expr) = r(Expr(e.head, process.(e.args)...))
process(e) = e
# Process the input expression
process(e)
end
For example, the following expression has both long sums and long products that we re-associate into binary operations.
let
e = :(sin(1+2+3*4*f))
println("Original: $e")
println("Transformed: $(reassoc_addmul(e))")
end
If we wanted to, we could also “de-associate” by collapsing sums and multiplications into single nodes.
#| output: false
function deassociate(e)
= @rules begin
r * (*(args...)) => (x, args) -> Expr(:call, :*, x, args...)
x *(args...)) * x => (args, x) -> Expr(:call, :*, args..., x)
(+ (+(args...)) => (x, args) -> Expr(:call, :+, x, args...)
x +(args...)) + x => (args, x) -> Expr(:call, :+, args..., x)
(e => e -> e
end
process(e :: Expr) = r(Expr(e.head, process.(e.args)...))
process(e) = e
process(e)
end
We take the output of our “reassociate” code as an example input.
deassociate(:(sin((1+2) + (3*4)*f)))
Automatically-generated code of the sort that comes out of naive automatic differentiation often involves operations that can be easily simplified by removing additions with zeros, multiplications with one, and so forth. For example, we write simplify_sum
and simplify_mul
routines to handle simplification of sums and products involving zeros and ones.
#| output: false
simplify_sum(args) =
let args = filter(x -> x != 0, args)
if length(args) == 1
1]
args[else
Expr(:call, :+, args...)
end
end
simplify_mul(args) =
let args = filter(x -> x != 1, args)
if any(args .== 0)
0
elseif length(args) == 1
1]
args[else
Expr(:call, :*, args...)
end
end
Building on simplifying sums and products, we can recurse up and down an expression tree to apply these and other similar rules.
#| output: false
function simplify(e)
= @rules begin
r +(args...) => args -> simplify_sum(args)
*(args...) => args -> simplify_mul(args)
- 0 => x -> x
x 0 - x => x -> :(-$x)
-(0) => (;x) -> 0
/ 1 => x -> x
x 0 / x => x -> 0
^ 1 => x -> x
x ^ 0 => x -> 1
x e => e -> e
end
process(e :: Expr) = r(Expr(e.head, process.(e.args)...))
process(e) = e
process(e)
end
We illustrate our simplifier with some operations
let
compare(e) = println("$e simplifies to $(simplify(e))")
compare(:(0*x + C*dx))
compare(:(2*x^1*dx))
compare(:(1*x^0*dx))
end
We write programs to be executed by a computer, an unforgivingly literal-minded audience. But the humans who read our code are a more challenging audience by far. We humans get bored. We get sidetracked. We misunderstand the simple things, and declare that we understand complicated things based on only-partly-correct models of how the world works. We read, forget, read again, and forget again. In the face of such an audience, what’s a programmer to do?
We seek to write code that accomplishes something worthwhile, and to write in a simple, direct style that can be easily understood by anyone. This is easier said than done. There are some classic books with practical guidance that transcends specific programming languages, and we highly recommend the books of Kernighan and Plauger (1978), Knuth (1984), Kernighan and Pike (1999), and Hunt and Thomas (1999). The style guide for the Julia programming language also provides useful guidance. But since we are here, we will provide our own two cents on what we think are some of the key issues to writing “good” Julia code.
Julia code tends to have some common formatting conventions that make it easier to read code written by varying authors:
CamelCase
lowercase
and omit separators when feasible (underscores can be used otherwise)We also recommend limiting globally visible names (modules, structures, or functions) to Latin characters, or providing a Latin alternative when a larger set of Unicode characters is used (for example, pi
as an alternative to π
).
There are a number of style guides with more detailed opinions on the proper formatting of Julia code. If you want to adhere tightly to such a guide, we recommend a tool like JuliaFormatter
.
Correctness is a tricky concept in numerical computing. Approximation is at the heart of most numerical methods, whether for data science or otherwise. Hence, we rarely get to decide whether we have “the right answer”; instead, we have to reason about whether the result of a computation is accurate enough for a downstream task. To make matters even more complicated, we usually face tradeoffs between speed and accuracy, and so need to design simultaneously for “right enough” and “fast enough.”
Neither correctness nor performance are easy to reason about. However, it is easier to tell when numerical code is slow than when it is wrong. Our mental models of why code is slow may often be wrong, but our observation that code is slow are not. Therefore, we usually want to start with simple code together with analysis and tests to ensure that it is accurate, and revise the code to tune performance later, informed by profiling (Chapter 3).
The earlier we catch an error, the easier it is to correct. We would rather:
Tests, type annotations, and careful assertions all help us with finding errors early.
For a human audience, functions are a natural unit for code. A well-designed function gives name to a concept, and ideally the name, comments, and code are all in obvious agreement – it “does what it says on the cover.” Small functions are also relatively easy to test. And assigning appropriate type annotations to the arguments and return values of a function also helps us catch errors in how the function is used.
Apart from readibility, there are good performance reasons for putting code in functions. The Julia compiler works with functions. Code written at a command line or in a script style will not be processed in the same way, and will typically run more slowly. It is also easier to use Julia’s analysis tools to understand the performance implications of design choices made in short functions.
Functions and control flow constructs in Julia create new scopes, so Julia variable names often have short lives. This is convenient for human readers, with our limited memories, particularly when we want to save typing and use short names. Variable names with limited scope (and limited opportunity for reassignment) are also easy for the compiler to analyze, and result in more efficient code.
Unlike some languages (e.g. Python, Java, and C++), Julia does not have a strong data hiding mechanism, though it is possible to code in a way that hides data structures (e.g. by tucking them inside functions). Nonetheless, if the option exists, it is usually preferable to use getter and setter functions rather than directly accessing the contents of a Julia structure. This makes it easier to monitor access to the data structures for debugging, and it makes it easier to change data structures as the code evolves.
A general design principle for any data system is that code and data should ideally have a single, concise, authoritative representation.
Compared to some languages, Julia code is fairly concise, and does not require much boilerplate6. Several aspects of the language contribute to this conciseness:
Peer to DRY in our opinions on design style are two other acronyms:
That is, we want to start with simple, working code that solves a known problem, and only get fancier if there is a measurable reason to do so. Simple code and simple data structures are easier for people to read, and preferable absent other real concerns. For example, any number of Julia constructs that we have discussed (promotion rules, macros, and expansion of tuples into parameter lists) are arguably harder to read than alternative constructions. It is nice to have these features in the language in situations to avoid writing redundant code or to address performance bottlenecks during profiling. But if we have not observed redundancy or performance issues in practice, it is usually best to write the simplest thing possible. In the event we decide that we do need a less simple code construct, whether to avoid repetition or to deal with a performance issue, we have the experience of writing the short and simple thing first. And, as with writing English prose, it is usually easier to revise what exists than to start from a blank page.
In the beginning, the creators of the C language allowed the same piece of memory to be treated as belonging to different types. This has made a lot of people very angry and been widely regarded as a bad move.↩︎
At the time of this writing, the default Int
type in Julia is almost always 64 bits. It is not so easy to find 32-bit systems any more.↩︎
In 2003, Joel Spolsky wrote a blog post on “The Absolute Minimum Every Software Developer Absolutely, Positively Must Know About Unicode and Character Sets (No Excuses!)”, which includes the memorably line “if you are a programmer working in 2003 and you don’t know the basics of characters, character sets, encodings, and Unicode, and I catch you, I’m going to punish you by making you peel onions for 6 months in a submarine. I swear I will.” That was two decades ago at the time of this writing. Unicode as a concept and UTF-8 as an encoding are both ubiquitous. There is no excuse for not knowing a little about it.↩︎
The notion of “null” is surprisingly subtle. This can refer to something that is “deliberately left blank” (e.g. a default value for an optional argument), something that is “well-defined but unknown” (e.g. a record of someone’s age) or something that is inappropriate to the task at hand (e.g. spouses’ name when not married).↩︎
As the Julia documentation points out, improvements to the compiler may mean this will no longer be an issue soon. This does not change the value of this as an example.↩︎
Per Wikipedia, the metal printing plates used for typesetting widely-distributed ads or syndicated columns were called “boilerplate” by analogy with the rolled steel plates used to make water biolers.↩︎