2 Julia programming

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.

2.1 Interacting with Julia

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.

2.1.1 The Julia prompt

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 $π$ (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 $π$ .
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.

2.1.2 Jupyter notebooks

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
IJulia.notebook()

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.

2.1.3 Pluto notebooks

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.

2.2 Simple expressions

We start our deeper dive into Julia with a treatment of simple expressions, the very simplest of which are literals (Table 2.1).

Table 2.1: Examples of literals in Julia

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

2.2.1 Assignment, scopes, and variable types

Like many other languages, the Julia assignment operation = represents binding a name to a value. For example,

x = 10

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.

L, U, p = lu(A)  # Assign multiple outputs

There is also a form of destructuring assignment for when the right hand side has named quantities, e.g.

p = (x = 1.0, y = 2.0, z = 3.0)
(; y, z) = p  # Sets y = p.y, z = p.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.

2.2.1.1 Scope

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
    p = 0
    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:

A series of assignment statements on the same line as the let, which result in bindings of variable names that are purely local to the let;
A code block that by default makes use of the bindings in the assignment block and any bindings in the outer scope.

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
    p = 0
    let
        p = 1      # Reassign outer version of 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
    p = 0
    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()
    x = 0
    for k = 1:5
        x = k    # x is the same variable as before
    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
        x = k    # x is local to this iteration
    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()
    xyz = 0
    let
        global xyz
        xyz = 1
    end
end

modify_global_xyz()
println(xyz) # Prints 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.

2.2.1.2 Variable types

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:

x :: Float64 = 0  # Implicitly converts the integer 0 to the floating point 0.0

Once the type of x has been specified, we cannot assign an incompatible value to x:

x :: Float64 = 0  # OK so far, x is now floating point 0.0
x = "Oops"        # Error!

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:

x :: 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!

We will have much more to say about types in ?sec-tbd.

2.2.2 Arithmetic operations

Table 2.2: Standard arithmetic operations

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 $2 i$ .

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,

x += 1  # Equivalent to x = x + 1
x *= 2  # Equivalent to x = x * 2

2.2.3 Logical operations

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
    x = true || begin println("Hello 1"); true end
    y = false || begin println("Hello 2"); true end
    x, y
end

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:

z1 = p ? x : y
z2 = if p then x else y end

2.2.4 Comparisons

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

2.2.5 Simple calls

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

C = 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

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

a = (1.0, 1.00000001)
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

b = (a, 2)          # Yields ((1.0, 1.00000001), 2)
bunpack = (a..., 2) # Yields (1.0, 1.00000001, 2)

2.2.6 Broadcasting operations

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} + 2 x + 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.

2.2.7 Indexing and slicing

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,

lst = [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]

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

lst = [5; 2; 2; 3]
sublst = lst[1:2]    # New storage for sublst
sublist[:] = [3; 4]  # Changes contents of sublst, not lst

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:

lst = [5; 2; 2; 3]
sublst = view(lst, 1:2)  # Provides a view into lst
sublist[:] = [3; 4]      # lst is now [3; 4; 2; 3]

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:

lst = [5; 2; 2; 3]
idx_big = (lst .> 2)  # [true, false, false, true]
lst[idx_big]          # Results in [5; 3]

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:

xyz = (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)

Our examples so far have been for single-index arrays (vectors), but everything generalizes natural to multi-dimensional arrays.

2.2.8 Structure access

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:

xyz = (x=4, y=5, z=6)  # Create a named tuple
xyz.z                  # Returns 6 (structure access syntax)

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.

2.2.9 Strings

One particularly convenient feature of Julia is string interpolation. For example, consider the following Julia code:

let
  sqrt2approx = 99.0/70.0
  "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 ".

2.3 Control flow

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.

2.3.1 Conditional statements

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.

2.3.2 Loops

A while loop executes while the loop condition is true

function gcd(a,b)
    while b != 0
        a, b = b, a%b
    end
    a
end

A for loop iterates over a collection. Perhaps most frequently, this is an integer range:

function mydot(x, y)
    @assert length(x) == length(y)
    result = 0.0
    for i = 1:length(x)
        result += x[i]*conj(y[i])
    end
    result
end

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)
    result = 0.0
    for (xi, yi) in zip(x, y)
        result += xi*conj(yi)
    end
    result
end

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.

2.4 Functions and methods

2.4.1 Defining functions

Julia provides three primary syntactic variants for defining functions:

Using the function keyword
Using an equational definition
Using the -> operator

The 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
        a, b = b, a%b
    end
    a
end

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:

f1 = function(x) x^2 + 2x + 1 end
f2 = x -> x^2 + 2x + 1

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.

2.4.2 Defaults and keywords

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)
    l = sort(l, by=by)
    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.

2.4.3 Closures

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)
    fx = f(x)
    for _ = 1:maxiter
        if abs(fx) < ftol
            return x
        end
        x -= fx/df(x)
        fx = f(x)
    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)
        fx = f(x)
        push!(fs, fx)
        fx
    end
  
    x = naive_newton(frecord, df, x, 
                     maxiter=maxiter, ftol=ftol)
    plot(filter(fx->fx > 0, abs.(fs)),
         xlabel="k",
         yscale=:log10, label="\$f(x_k)\$")
end

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)

2.4.4 Methods

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]
        [c*x-s*y; s*x+c*y]
    end

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)
    result
end

let
    t = ((1,2), (3,(4,5)), 6)
    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.

2.4.5 Operator overloading

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
    coeffs :: Vector{Float64}
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)
      n = max(length(p.coeffs), length(q.coeffs))
      c = zeros(n)
      c[1:length(p.coeffs)] = p.coeffs
      c[1:length(q.coeffs)] += q.coeffs
      Polynomial(c)
end

function Base.:*(p :: Polynomial, q :: Polynomial)
      n = length(p.coeffs) + length(q.coeffs) - 1
      c = zeros(n)
      for (i, pi) in enumerate(p.coeffs)
            c[i:length(q.coeffs)+i-1] += pi*q.coeffs
      end
      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)
    px = 0*x
    for c in reverse(p.coeffs)
        px = px*x+c
    end
    px
end

With these definitions in place, we can write code like

let
    x = Polynomial([0; 1])  # x
    p = x*(2*x+1)
    p(2)
end

10.0

The ambitious reader might also try writing subtraction, polynomial division and remainders.

2.4.6 Higher-order functions

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.

2.4.7 do syntax

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)
           ispositive = x > 0
           action = ispositive ? "Keep it" : "Toss it"
           println("Testing x = $x: $action")
           ispositive
        end, [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
    ispositive = x > 0
    action = ispositive ? "Keep it" : "Toss it"
    println("Testing x = $x: $action")
    ispositive
end

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
    body
end

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.

2.4.8 Composition and pipes

In addition to nested calls (e.g. f(g(x))), Julia supports two special syntax constructs for function composition:

The composition operator (∘) acts like composition in mathematics.
The pipeline operator (|>) 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)))
norm2b = sqrt ∘ sum ∘ sqr
norm2c(v) = v |> sqr |> sum |> sqrt

2.4.9 Generators and comprehensions

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

g = ((i,j) for i=1:2, j=1:2)

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
    g = ((i,j) for i=1:2, j=1:2)
    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)

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:

xgrid = range(0.0, 1.0, length=20)
kxx = [k(x,y) for x in xgrid, y in xgrid]

2.5 Types

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 time¹. 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
        a, b = b, a%b
    end
    a
end

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.

2.5.1 Type hierarchy

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.

2.5.2 Numeric types

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.²
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 + 1 i$ of type Complex{Float64}.

Julia also supports arbitrary-precision BigInt and BigFloat types, though we will rarely make use of these.

2.5.3 Strings and symbols

Strings in Julia inherit from the AbstractString type, and are logically sequences of Unicode characters³. 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.

2.5.4 Tuples

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

2.5.5 Structs and parameterization

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
    q
end

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
    p :: Integer
    q :: Integer
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
    p :: T
    q :: T
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
    p :: T
    q :: T
    MyRational(p :: Integer, q :: Integer) = 
        if q == zero(q)
            error("Divide by zero")
        else
            p, q = promote(p, q)
            r = gcd(p, q)
            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.

2.5.6 Mutable structs

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
    p :: Int
    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)

foo = FooCounter()        # foo.p initialized to 0
println(increment!(foo))  # Prints 1
println(increment!(foo))  # Prints 2

2.5.7 Array types

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

2.5.8 Other collections

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.

2.5.9 Unions, Nothing, and Missing

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

2.5.10 Types of types and value types

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.

2.5.11 Type conversions and promotions

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} =
    MyRational{promote_type(T,S)}

2.6 Exception handling

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) = 
    (x >= zero(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) =
    try
        log(x)
    catch e
        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.

2.7 Documentation

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.

2.8 Vectors, matrices, and linear algebra

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.

2.8.1 Building matrices and vectors

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

w = v'

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
    v = [1, 2]  # v is a 2-element Vector{Int64} containing [1, 2]
    w = v'      # w is a 1-2 adjoint(::Vector{Int64}) with eltype Int64
    v[2] = 3    # Now v contains [1, 3] and w is the adjoint [1, 3]'
end

Julia gives us several standard matrix and vector construction functions.

Z = zeros(n)   # Length n vector of zeros
Z = zeros(n,n) # n-by-n matrix of zeros
b = rand(n)    # Length n random vector of U[0,1] entries (from Random)
e = ones(n)    # Length n vector of ones
D = diagm(e)   # Construct a diagonal matrix
e2 = diag(D)   # Extract a matrix diagonal

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

2.8.2 Concatenating matrices and vectors

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.

y = [1; 2]        # Length-2 vector
y = [1 2]         # 1-by-2 matrix
M = [1 2; 3 4]    # 2-by-2 matrix
M = [I A]         # Horizontal matrix concatenation
M = [I; A]        # Vertical matrix concatenation

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.

2.8.3 Transpose and rearrangemenent

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
avec = reshape(A, prod(size(A)));
A = reshape(avec, n, n)
  
idx = randperm(n)   # Random permutation of indices (need to use Random)
Ac = A[:,idx]       # Permute columns of A
Ar = A[idx,:]       # Permute rows of A
Ap = A[idx,idx]     # Permute rows and columns

2.8.4 Submatrices, diagonals, and triangles

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:

A = 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

Ad = diag(A)      # Diagonal of A (as vector)
A1 = diag(A,1)    # First superdiagonal
Au = triu(A)      # Upper triangle
Al = tril(A)      # Lower triangle

Other operations give a view of the matrix without making a copy of the contents, which can be much faster:

A = randn(6,6)          # 6-by-6 random matrix
view(A,1:3,1:3)         # View of leading 3-by-3 submatrix
view(A,:,3:end)         # View of columns 3-6
Au = UpperTriangular(A) # View of upper triangle
Al = LowerTriangular(A) # View of lower triangle

2.8.5 Matrix and vector operations

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.

y = d.*x   # Elementwise multiplication of vectors/matrices
y = x./d   # Elementwise division
z = x + y  # Add vectors/matrices
z = x .+ 1 # Add scalar to every element of a vector/matrix
  
y = A*x    # Matrix times vector
y = x'*A   # Vector times matrix
C = A*B    # Matrix times matrix

# Don't use inv!
x = A\b    # Solve Ax = b *or* least squares
y = b/A    # Solve yA = b or least squares

2.8.6 Things best avoided

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?

2.9 Useful packages

Beyond the LinearAlgebra package mentioned in Section 2.8, there are several other packages that we will use or recommend in this book:

2.9.1 Statistics

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.

2.9.2 Sparse matrices

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

2.9.3 DataFrames

The DataFrames package provides functionality for manipulating tabular data similar to the Python Pandas package.

2.9.4 Automatic differentiation

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.

2.9.5 Plots and tables

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.

2.9.6 I/O

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.

2.9.7 Development tools

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.

2.10 Macros and metaprogramming

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.

2.10.1 Code as data

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

quoted_sum = :(1 + 1)

:(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)

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
    x = 1
    :($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
    x = 1
    Expr(:call, :+, x, 1)
end

:(1 + 1)

2.10.2 Macros

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"

2.10.2.1 A toy definition

As a toy example, of a macro definition, consider a macro that prints the input expression

macro printarg(e)
    s = "$e = "
    quote
        let v = $e
            println($s, v)
            v
        end
    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

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
    x = 1
    @show x + x
end

x + x = 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.

2.10.2.2 A modest example

As an example of a more significant small utility macro, consider a closure defined inside of an outer function. Such closures sometimes suffer⁵ 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}())
    args = e.head == :call ? e.args[2:end] : e.args
    for arg in args
        symbols!(arg, s)
    end
    s
end

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
    s = Set{Symbol}()
    sarg = Set{Symbol}()
    symbols!(e.args[2], s)
    symbols!(e.args[1], sarg)
    setdiff!(s, sarg)

    # Bind local variables to the outer version
    bindings = [quote $v = $(esc(v)) end for v in s]

    # 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
    C = 10
    f = @letclosure x -> C*x
    C = 20
    f(5)
end

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.

2.11 A matching example

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]
            result = :(-$x)   # Create a "-x" expression
            true, result      # Yes it matched, return result
        end
    else
        false, nothing        # No, it didn't match
    end

That is, we want a function that takes an input

Is this an expression with a binary - operation at the top?
Is the second term in the sum a zero?
If both are true, bind x to the second term in the sum, produce a new expression :(-$x), and assign to result. Return the pair (true, result).
Otherwise, return (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

pattern => (argument_symbols) -> result_code

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

2.11.1 Preprocessing

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

2.11.2 Matching

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

A dictionary of bindings, mapping symbols to values.
A pattern
A symbol naming the current expression being matched.

2.11.2.1 Literals

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)

2.11.2.2 Symbols

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:

First use: we create a new binding for it.
Repeat use: Check if expr matches the previous binding.

match_gen!(bindings, e, s :: Symbol) =
    if !(s in keys(bindings))
        qs = QuoteNode(s)
        :($e == $qs)
    elseif bindings[s] == nothing
        binding = gensym()
        bindings[s] = binding
        quote $binding = $e; true end
    else
        binding = bindings[s]
        :($e == $binding)
    end

We illustrate all three cases below.

let
    bindings = Dict{Symbol,Any}(:x => nothing)
    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

2.11.2.3 Expressions

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),
        argmatch = match_gen_args!(bindings, e, pattern.args)
        :( $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)
           for (e,p) in zip(exprs, patterns)])

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
        nargs = length(patterns)
        lencheck = :(length($e.args) == $nargs)
        args = Vector{Any}([gensym() for j = 1:length(patterns)])
        argstuple = Expr(:tuple, args...)

        # De-splat pattern / splat arg assignment and
        #   adjust the length check
        if is_splat_arg(bindings, patterns[end])
            patterns = copy(patterns)
            patterns[end] = patterns[end].args[1]
            argstuple.args[end] = Expr(:(...), argstuple.args[end])
            lencheck = :(length($e.args) >= $(nargs-1))
        end

        argchecks = match_gen_lists!(bindings, args, patterns)
        :($lencheck && let $argstuple = $e.args; $argchecks end)
    end

We strongly recommend the reader trace through this code for some examples until enlightenment strikes.

2.11.2.4 Compiling the match

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)
    bindings = Dict{Symbol,Any}(s => nothing for s in symbols)

    # Input expression symbol and matching code
    # (symbol bindings are indicated by bindings table)
    expr = gensym()
    test = match_gen!(bindings, expr, pattern)

    # Get result vals (symbol/nothing) and associated variable names
    result_vals = [bindings[s] for s in symbols]
    declarations = filter(x -> x != nothing, result_vals)

    # Produce the matching code
    results = Expr(:tuple, result_vals...)
    :($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")

    pattern = filter_line_numbers(pattern)
    matcher = compile_matcher(symbols.args, pattern)
    :($matcher ∘ filter_line_numbers)
end

2.11.3 Rules

A rule has the form

pattern => args -> code

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

x - y => (x,y;e) -> process(e)

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_rule = @match (pattern, args, code) pattern => args -> code
    isok, rule_parts = match_rule(rule)
    if !isok
        error("Syntax error in rule $rule")
    end
    pattern, args, code = rule_parts

    println("Pattern: $pattern")
    println("Args: $args")
    println("Code: $code")
    match_arg1 = @match (args, expr) (args..., ; expr)
    match_arg2 = @match (args, expr) (args..., )
    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")
    symbols, expr_name = bindings
    
    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, code
end

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)
    symbols, expr_name, pattern, code = parse_rule(rule)
    bindings = Dict{Symbol,Any}(s => nothing for s in symbols)
    test = match_gen!(bindings, expr, pattern)

    # Get list of match symbols and associated declarations
    result_vals = [bindings[s] for s in symbols]
    declarations = filter(x -> x != nothing, result_vals)

    # Set up local bindings of argument names in code
    binding_code = [:($s = $(r == nothing ? (:nothing) : r))
                    for (s,r) in zip(symbols, result_vals)]
    if expr_name != nothing
        push!(binding_code, :($expr_name = $expr))
    end

    # Produce the rule
    ismatch = gensym()
    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)
    expr, result = gensym(), gensym()
    code = compile_rule(filter_line_numbers(r), expr, result)
    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)
    rblock = filter_line_numbers(rblock)
    if rblock.head != :block
        error("Rules must be in a begin/end block")
    end

    # Set up input name, output name, and rule list
    expr, result = gensym(), gensym()
    rules = rblock.args

    # 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

2.11.4 Examples

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.

2.11.4.1 Re-associating operations

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
    r = @rules begin
        +(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)
    r = @rules begin
        x * (*(args...)) => (x, args) -> Expr(:call, :*, x, args...)
        (*(args...)) * x => (args, x) -> Expr(:call, :*, args..., x)
        x + (+(args...)) => (x, args) -> Expr(:call, :+, x, args...)
        (+(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)))

2.11.4.2 Simplifying expressions

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
            args[1]
        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
            args[1]
        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)
    r = @rules begin
        +(args...) => args -> simplify_sum(args)
        *(args...) => args -> simplify_mul(args)
        x - 0 => x -> x
        0 - x => x -> :(-$x)
        -(0)  => (;x) -> 0
        x / 1 => x -> x
        0 / x => x -> 0
        x ^ 1 => x -> x
        x ^ 0 => x -> 1
        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

2.12 Elements of Julia style

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.

2.12.1 Formatting conventions

Julia code tends to have some common formatting conventions that make it easier to read code written by varying authors:

Indentation by four spaces
Functions that modify (mutate) an input end with an exclamation mark
Module and type names use CamelCase
Function names are 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.

2.12.2 Correctness first

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

2.12.3 Catch errors early

The earlier we catch an error, the easier it is to correct. We would rather:

Catch conceptual errors before we start writing code.
Catch implementation errors at compile time than at run time.
Analyze errors in an individual function than have to perform a root cause analysis.

Tests, type annotations, and careful assertions all help us with finding errors early.

2.12.4 Put it in a (small) function

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.

2.12.5 Keep the scope small

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.

2.12.6 Respect interfaces

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.

2.12.7 DRY (Don’t Repeat Yourself)

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 boilerplate⁶. Several aspects of the language contribute to this conciseness:

The ability of the system to dispatch to different methods based on argument types simplifies the process of writing generic code.
Functions are first-class entities that can be passed around with context, so we can pass around specialized functions for a specific subtask rather than coding it everywhere. For example, we can pass logging routines into our functions as callbacks rather than adding logging infrastructure to all our routines.
Conversion methods and promotion rules limit the number of variants of methods that we might have to write for different combinations of argument types.
When other methods fail, we can remove redundancy by writing macros to generate repetitive code for us.

2.12.8 KISS and YAGNI

Peer to DRY in our opinions on design style are two other acronyms:

KISS: Keep It Simple, Stupid
YAGNI: You Ain’t Gonna Need It

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