“Looks like Python, feels like lisp, runs like C.”

Introduction

As discussed in the previous article, I plan to explore the capability of Julia in scientific computing, esp. in applications like computational mechanics and optimization. I have gotten even more attracted to Julia when I heard that it has joined the Petaflop Club: achieved peak performance exceeding one petaflop per second on supercomputers. Julia appears to be the only “high-level” language in the Petaflop club - the others are C, C++ and Fortran.

This article is a digest of the basics of Julia that I have learned recently. The syntax of Julia will not be discussed in detail, as it could be mastered by actaully writing the code - So there will not be many examples. Instead, in this article, I would focus more on the features and highlights of Julia. In addition, due to my background, I tend to highlight the similarities and differences between Julia and Python/C++.

Overview

Two-Language Problem

One of the goals of Julia is to solve the two-language problem: one can either have a programming language that’s easy to use, or fast.

There are indeed some existing solutions, e.g. those from the Python ecosystem. In the standard CPython implementation, the Python code is interpreted as bytecodes that run on the virtual machine, which is of course slower than compiled machine code that runs on the actual machine. Three representative approaches to accelerate Python include PyPy, Cython, and Numba. PyPy is basically an implementation different from CPython using a just-in-time (JIT) compiler. Cython generates C code from Python code - or wraps C libraries for Python ¹. Numba compiles Python snippets in a JIT manner. I myself have mainly used Python with Cython-wrapped C++ library. So I am really using two languages to “solve” the two-language problem. The three approaches have been discussed and compared in this post. These approaches are not perfect. To accelerate a given Python code, some extra work are involved, e.g. annotating the code using certain “peculiar” syntaxes.

The two-language problem had been even more challenging when the Julia project started in 2009, since the Python ecosystem has not matured enough to solve the problem. As one of the developers discussed in this post: The initial version of Cython was released in 2007; The first working verion of PyPy was released in 2007; Numpy was yet not mature at the time; Numba would not be released until 2012. ²

Solution by Julia

Julia is mainly written in Julia with some dependencies on C, Fortran, and assembly. However, in the core, there is a Scheme parser FemtoLisp and a LLVM compiler in C++. They are responsible for gluing Julia and the underlying operating system together.

The Julia solution to the two-language problem is best described by the quote at the beginning of this article:

1. Writes like Python: The user writes Julia like a high-level dynamic language.
2. Feels like Lisp: The Julia code is parsed into abstract syntax trees (ASTs) using FemtoLisp (lowered, typed).
3. Runs like C: The LLVM generates the intermediate representation (IR) of the ASTs, which would run on the virtual machine (llvm), and compiles the IR into the native code optimized for the local machine (native) ³.

The user essentially only needs to care about writing their applications. The compiler takes care of the rest. Of course, the code would be faster if it is written with care taken for the parsing and compiling stages.

Key Features

As I understand, Julia as a dynamic programming language is made distinct from the others by its type system combined with multiple dispatch. Well, I am not a computer scientist, I could be wrong in this assertion and some of the concepts discussed in this section.

Dynamic languages are usually considered “typeless”. However, more precisely speaking, types do exist, but what lacks is type declaration. The user cannot explicitly instruct the compiler to specify certain types for certain values. On the other hand, in static languages, type declaration is mandatory. And that is one reason why static languages are usually faster than dynamic languages, as some optimizations of the machine code can be done based on the types. Cython and Numba are like patches of partial type declaration for accelerating Python. Julia takes a step further, as it has the complete type system and the type inference mechanism. Like static language, all the values in Julia can be assigned a type by the user. Like dynamic language, not all values require a type - the most suitable type can be inferred by the compiler.

Based on the type system, the multiple dispatch mechanism is enabled with parametric polymorphism. In polymorphism, a single interface, or method, is provided to values of different types ⁴. Typically, the polymorphic methods are stored in a virtual method table, a.k.a. vtable. In the run-time, the specific method is chosen from the vtable. In single dispatch, the vtable is stored with the classes. The calling of a method can be written as object.action(args), because the choice of polymorphic method is solely determined by the type of the calling object. In multiple dispatch, the vtable is stored with the method itself. The calling of a method should be better written as action(objects, args), as is done in Julia, the choice of polymorphic method is determined by the types of all the arguments.

Multiple dispatch has two advantages. First, it is the direct reason for the high performance of Julia. For each polymorphic method, specialized low-level code is generated with all possible optimizations under the host machine. In the run-time, the most suitable/optimized method is called from the vtable. Second, multiple dispatch provides greater flexibility and reusability, as a method is not bound to a type. For example, suppose there is an existing method for LU decomposition for matrices of Float64. To do LU decomposition of matrices of rationals, one only need to define the elementary arithmetic operations for rationals, and call the same LU method. However, in single dispatch languages, one will need to not only define the elementary arithmetic operations, but also rewriting the LU method using the new rational type.

Besides the type system and multiple dispatch, another important feature of Julia is its metaprogramming capability: a program can be designed to read, generate, analyse or transform other programs, and even modify itself while running. The roles of code and data are more or less interchangeable. In Julia, everything is an expression that returns a value and an expression can take values and/or expressions as input. The metaprogramming capability directly stems from the lisp core of Julia, that enables the manipulation of the ASTs generated from the Julia code, e.g. using macros ⁵ - functions that acts on expressions. The metaprogramming capability allows the users to minimize the efforts and computational time for their applications.

Learning materials

The best learning material for Julia is probably its documentation, which is always (expected to be) update-to-date, and contains both introductory and advanced topics. The official website also provides a large list of learning materials, including videos (e.g. this and this), online documents, and books. Among these are two books that I recommend: (1) Getting Started with Julia by Ivo Balbaert, indeed good for getting started; (2) “Julia High Performance” by Avik Sengupta, a little advanced but very useful for beginners aiming for performant programs. Another book, “Mastering Julia” by Malcolm Sherrington, might be good for people who need immediate hands-on experiences. However, a notable issue is that, with the advent of the more stable 1.0 version, some of the syntaxes and claims in the books are outdated. So make sure to practice and verify everything read from the books.

Another good source of learning materials are the Julia code itself. As one of the books has suggested, there are several interesting folders in Julia to look at:

1. base: contains a great portion of the standard library and the coding style exemplary.
2. test: has some code that illustrates how to write test scripts and use the Base.Test system.
3. examples: gives Julia’s take on some well-known old computing chestnuts.

This article is mainly based on materials from Getting Started with Julia by Ivo Balbaert and the official documentation.

More on Features

Type System

Type hierarchy

Every value has a type and there are a lot of built-in types. These types are organized in tree, with the root, type Any. The types themselves are of type DataType. A type is either an abstract type or a concrete type. An abstract type can have a lot of children or subtypes but cannot be instantiated, while a concrete type is on the contrary. Essentially, an abstract type is only a name that groups multiple subtypes together, which works in a manner similar to Union that gives a type/alias to a user-specified group of types. Naturally, Any is the supertype of all types (including itself). The opposite of Any is a special type called Union{}, i.e. union of nothing. It is subtype of all other types (even the concrete ones). There is also a special type Nothing.

Here, it is stressed again: In the context of multiple dispatch, a subtype do not inherit any methods from its supertype. This is on the contrary to the object-oriented thinking.

Conversion, promotion, and annotation

Variables of different types can be converted to each other, if the corresponding convert method exists. There is also a built-in system called automatic type promotion to promote arguments of functions so as to match those of the existing functions. This is a key device that facilitates the multi-dispatch mechanism. Additionally, the combination of type promotion and Union could bring some convenience to the user. The user can define one function using data types from a user-defined union and save the effort of defining multiple similar functions for each data type from the union.

Variables can be associated to types using two operators for type annotation: :: and <:. Essentially, var::T asserts if the variable var is of type T. This is also used in the sense of a type declaration, but only in local scopes such as in functions. var<:T asserts if the variable var is of the subtype of T. This can be used in, e.g. the restriction of types in composite types.

Composite and parametric types

The composite types are the types composed of a set of named fields with an optional type annotation (by default Any). A composite type is similar to a struct in C. In fact, the keyword is struct too in Julia.

Composite types declared with struct are immutable; they cannot be modified after construction. This design has some advantages, including efficiency in storage and memory allocation, thread safety, and readability. Nevertheless, the contents of a mutable field in an immutable type can be changed, because only the reference is stored. It could be a good practice to define immutable types without mutable fields. When necessary, mutable composite types can be declared with the keywords mutable struct.

Composite types can be generalized via parametrization. A parametric type is conceptually similar to templates in C++ and creates a whole family of new possible concrete types. However, the new types are only compiled as needed at runtime. Finally, for better performance, always annotate the type of the fields in a composite type when possible.

As a side note, to see whether the two objects x and y are identical, they must be compared with the is function, The is(x, y) function can also be written with the triple-equal operator as x === y. The addresses in memory are compared to check whether they point to the same memory location

Constructors

Constructors creates and initializes an instance of a type. There is always a default (inner) constructor that comes with a type, if not otherwise specified. One can, of course, define a specific inner constructor in the type declaration, which overrides the default one. The inner constructors are the only constructors that has access to a special locally existent function called new that creates an object of the declared type. Once a type is declared, there is no way to add more inner constructor methods. Instead, one can only add outer constructors, that are defined just like regular methods (due to multiple dispatch). The outer constructors, ultimately, rely on calling one of the inner constructors to complete the instantiation.

Special inner constructors are needed in certain cases, e.g. when argument checking is necessary. A good practice is to combine outer constructors with a limited set of inner constructors.

Metaprogramming

This is the part of Julia that feel like Lisp: “code is data and data is code”. Every piece of code is internally represented as an expression, Expr, which can be manipulated by other expressions. This characteristic of Julia, homoiconicity, enables the transformation and generation of the code. This is not unfamiliar, if one works with Lisp (e.g. in Emacs) and Mathematica.

The homoiconicity facilitates the reflection capability of Julia. The structure of a program and its types can be explored programmatically just like any other data. This means that a running program can dynamically discover its own properties. Reflection is indispensable for tools that need to inspect the internals of the code objects programmatically, such as debugging and profiling.

Expressions

The Julia code is parsed by FemtoLisp as an AST, a tree representation of the abstract syntactic structure of the source code. The nodes of AST are all of type Expr. A formal definition of the Expr type is

1
2
3
4

type Expr
    head::Symbol
    args::Array{Any,1}
end

The head indicates the type of the expression, e.g. :call as a function call and = as assignment. The args is an array of Symbols and (maybe) Exprs, which recursively represents the AST. One can get a formatted output of the AST of an expression using the dump function.

Expressions can be built directly using the constructor for Expr. However, to make life easier, one can create expressions using the quote operator : that treats its argument as data instead of code. It means “prevented evaluation” in a sense. The quote preceding a symbol makes a Symbol, e.g :sym. The quote preceding a bracketed expression makes an Expr, e.g :(1+2).

A compound expression consists of multiple sub-expressions separated by ;. The value of a compound expression is the value of the last sub-expression. For a larger compound expression, it is better to use the quote-end block instead of quote operator. Note that the evaluation of an expression is done in the global scope and the internal variables will be visible to the outside code.

Writing expressions is made even simpler with the interpolation operator $. The interpolation is evaluated once an expression is constructed at parse time. This is unlike the quotation itself, which is evaluated by eval at runtime.

Macros

Macros are like functions, but instead of values, they take expressions, symbols or literals as input arguments. The evaluation of a macro means the expansion of the input expression. This expansion occurs at parse time when the AST is being built. This is different from a function, which takes the input values and returns the computed values at run-time.

A macro is hygenic, when it differentiates between the macro context and the calling context. The hygienity is important because the macros inject the code directly in the namespace in which they are called, which may clash with the existing methods in the same module. There are several rules for writing hygienic macros:

1. Declare the variables used in the macro as local, so as not to conflict with the outer variables.
2. Use the escape function esc to make sure that an interpolated expression is not expanded, but instead is used literally.
3. Don’t call eval inside a macro.

Note that escaping is similar to Hold in Mathematica. An expression wrapped in this manner is left alone by the macro expander and simply pasted into the output verbatim. Therefore it will be resolved when the macro is called, instead of when it is defined.

Notes on the Usage

This section highlights some basic functionalities of Julia.

General Behaviors and Conventions

Julia behaves similarly to Python in many senses, e.g. all names are just aliases bound to values. However, there are differences. For example, in Julia, arrays are 1-indexed. Negative index is not supported and the last index is end. Refer to this for more differences between Julia and other languages.

Entities in Julia must start with a letter or an underscore _. Conventionally variable names consist of lowercase letters with long names separated by underscores rather than using camel case. Numbers can be separated by underscores for readability. It is a good practice to write constant variables in uppercase letters. The package names are kept upper camel case. ! indicates a mutating function that modifies its the first argument.

Julia does not have line-continuation operators. It will read until the expression is finished.

Julia supports unicode, so that functions and variables can be named using special characters, e.g. ⊗ instead of * for outer product, and α instead of alpha for a Greek symbol.

Data structures

Concepts like strings, arrays, dictionaries, tuples, and sets all behave in a manner that is more or less similar to those in Python. Note that except string, all the rest are parametric types.

Strings

Julia differentiates between single characters and strings. Strings can contain any number of characters and are specified using double quotes " or triple double quotes """. The character literals are specified using single quotes '. """ is used (esp.) when there are double quotes in the string itself.

The strings are immutable, which means that they cannot be altered once they have been defined

Julia has an elegant string interpolation mechanism for constructing strings using the $ operator (again). $var inside a string is replaced by the value of var, and $(expr) is replaced by the computed value of expression expr. The @printf macro takes a format string and one or more variables to substitute into this string while being formatted. It works in a manner similar to printf in C. In string interpolation, another @sprintf from Printf package is used.

Strings are joined by *, unlike + in Python. As a result, ^ applied on a string indicates an implicit multiplication.

Unicode characters in a string contains multiple bytes. Therefore, there are two ways to define the length of a string: (1) Number of bytes in string, endof(str); (2) Number of characters in string, length().

There are several special types of string.

1. v"" for version number, which can be compared.
2. r"" for regular expressions, which can be used for pattern matching.
3. b"" for byte array literal, which represents a string using Uint8 values.

Arrays

An array is of parametric type: Array{Type,Dimension}. For example, Array{Int64,2} means a 2-D array containing values of type Int64. Additionally, there are short-hands for vectors and matrices: Array{Type,1} == Vector{Type} and Array{Type,2} == Matrix{Type}.

The arrays are by default “real” arrays, meaning that the data is stored continuously in memory (column-major). This is in contrast to the list of lists in Python and vectors of vectors in C++. It is possible to create array of arrays in Julia, of course. However, the array would only contain the references to the underlying arrays.

There are multiple ways to create an array,

1. Direct specification, e.g. a=[1,2,3].
2. Using some methods, e.g. collect(range(...))
3. List comprehension, e.g. [i for i in 1:3]. Multiple for’s result in multi-D arrays.

Note that due to the column-major storage, 1-D array is logically column vectors. Logically row vectors are essentially 1-by-N matrices. To create row vectors, separate the values by space instead of comma, e.g. a=[1 2 3]. Furthermore, a semicolon ; indicates a new row, leading to a matrix, e.g. a=[1 2 3;4 5 6].

In the creation of an array, constraining the type is often helpful for the performance. Using typed comprehensions everywhere for explicitness and safety in production code is certainly a best practice. For example, use a=Int64[1,2,3] instead of a=[1,2,3]. When dealing with large arrays, it is better to indicate the final number of elements for pre-allocation using the sizehint method to improve the performance. The memory of a large array can be freed by setting the array to nothing.

Dictionaries and enumerations

A dictionary is of parametric type: Dict{Type_key,Type_value}. It is created by by Dict(key=>val, ...), where => is a pair operator. All the keys must have the same type, and the same is true for the values. The types can be specified at creation. A good practice is to use immutable types for the key, e.g. a Symbol ⁶.

Dictionaries are mutable. New key-value pairs can be appended to a dictionary without calling method like update in Python. Also, one can directly loop over keys, values, or key-value pairs of a dict. In Python, one need to call iteritems to loops over key-value pairs.

Julia does not provide a standard enumeration type. However, there is a @enum macro that creates a type that behaves like a Enum,

1

@enum EnumName[::BaseType] key[=val] ...

The data is stored as Vector{Tuple{Symbol,Integer}}.

Tuples

The type of a tuple is just a tuple of the types of the values it contains. Tuples are immutable, just like other languages. A tuple can be unpacked or deconstructed, e.g. a, b = tuple1. Notice that an error will not be raised if the left-hand side cannot take all the values of the tuple.

Finally, note that there is a type NamedTuples, which can be considered an ordered dictionary. Its underlying structure consists of a tuple of all the keys as Symbol and a tuple of all the values.

Control flow

Basic Blocks

One thing to be kept in mind is that, everything is an expression that returns a value. That means, structs like if-else-end and try-except-end all return a value, which can be assigned to a variable.

The familiar if-elseif-else-end syntax is used for conditionals. There are two short-hand versions: (1) Ternary operators: a ? b : c, and (2) Short-circuit evaluation ⁷,

1. cond && expr, if cond=true, do expr
2. cond || expr, if cond=false, do expr

Like Python, Julia has no switch/case statement, and the language provides no built-in pattern matching

The familiar while-end and for-end syntaxes are used for loops, which can be combined with break and continue. For the for loop, the iteration can be done on a range, a string, an array, or any other iterable collection, just like Python. Furthermore, enumerate also exists in Julia. One thing special, though, is that multiple for loops can be combined into one block - this is how multi-for list comprehension works. Note that for loops are typically faster than vectorization in languages like Python. That is because (1) in Python, the implementation of vectorization (e.g. by Numpy) is faster than naive loops generated by the interpreter; (2) in Julia, vectorization requires extra memory allocation and is thus slower than naive loops. Nevertheless, list comprehension can be far slower than the for loops and vectorizations.

The exceptions can be handled using the familiar try-catch-finally-end construct. However, the try-catch should not be used in performance bottlenecks. Whenever possible, all possible exceptions should be handled via conditionals. There are 24 predefined exceptions that Julia can generate, as of writing. These exceptions can be throw’d to interrupt the execution. User-defined exceptions can be created by deriving from the base type, Exception. Finally, Exceptions with error/warning/info messages are handled by the error, warn, and info methods.

Scopes

The for, while, and try blocks (but not the if blocks) all introduce a new scope. There are several things to note.

To create a new local binding for a variable, one should use the let block. This is similar to the cell-var-from-loop in Python discussed in a previous article, where a local variable is used in the definition of a function that would be used outside the current block.

The for loops and comprehensions differ in the way they scope an iteration variable. When i is initialized before a for loop, after executing a for loop iterating on i, the variable i will be updated. However, this is not true for comprehension.

Tasks/Coroutines

Julia also supports the yield mechanism, like Python. This can done by the combination of put! and take! methods associated with the Channel type. There is also a more fundamental Task type to define the task to be executed by put! and take!. There is also a convenient @task macro for creating a task from a function. ⁸

Finally, note that currently Julia tasks are not scheduled to run on separate CPU cores (yet?).

Functions

There are several ways for creating a function:

1. Code block of function-end
2. One-liner: function(variables) = expression
3. Anonymous function: function = variables -> expression

When the performance is important, try to use named functions instead of anonymous ones, because calling anonymous functions involves a huge overhead. Anonymous functions are mostly used when passing a function as an argument to another function, or as output of a function. That latter case is related to closures.

The input and output of a function are essentially specified using tuples. There are three types of arguments: (1) Normal arguments; (2) Optional positional arguments; (3) Optional keyword arguments. A typical arglist would be

1

func(arg, pos_arg=val1; kw_arg=val2)

Only the first two types of arguments, the positional ones, are considered in the multiple dispatch mechanism. The optional keyword arguments are mainly introduced for code clarity. The last positional argument can be “spliced”: arg..., meaning that it would take all the extra positional arguments into one tuple. Note that the splice operator ..., or splat, can be also used to pass an array to a function as individual arguments rather than an array as a whole. This is similar to the * operator in Python.

It is important to realize that in Julia, all the arguments to functions (with the exception of plain data such as numbers and chars) are passed by reference. Their values are not copied when they are passed, which means they can be changed from inside the function, and the changes will be visible to the calling code.

Higher-order function in Julia just works like that in Python. There are methods like map, filter, and broadcast. Note that the short-hand for broadcast is the dot .. This is the underlying mechanism for element-wise operations in arrays, such as .*. Also, wise applications of the dot operator can avoid extra memory allocations in vectorized loops, or making a devectorized loop look like a vectorized one - See example #6 here.

Here are some general tips for writing performant functions/program:

0. Refrain from using global variables. If unavoidable, make them constant, or at least annotate the types.
1. Split the complete program in functions, small functions, that work on local variables.
2. Types should be stable: Avoid changing the types of variables during execution.
3. Types could be inferred: Always type-annotate the arguments using concrete types or type unions.
4. The return type should only depend on the types of the arguments.
5. Avoid using the splat operator for dynamic lists of keyword arguments.

Advanced Topics

Two advanced topics in Julia are briefly touched here and will be investigated more in the future.

Calling other languages from Julia

Shell

Julia offers an efficient shell integration through the run function, which creates an object of type Cmd that is defined by enclosing a command string in backticks `. Julia forks commands as child processes from the Julia process. The Cmd object is not executed immediately. It can be run, connected to other commands via pipes, and read or write to it. The pipe operator |> can be used to redirect the output of a Cmd as the input to the following Cmd, just like that in shell environment ⁹.

C and Fortran

Julia’s LLVM compiler generates native code, the C functions can be called directly by Julia without any glue code. Calling a C function from Julia has exactly the same overhead as calling the same function from C code itself. However, care should be taken when dealing with arguments. First, one needs to work with pointer types, i.e. a native pointer Ptr{T} representing the memory address for a variable of type T. Second, the data consists of bits, e.g. Int8, Uint8, Int32, Float64, Bool, and Char, are considered bitstype. These will be viewed as contiguous byte arrays from C. The calling basically takes the following form,

1

result = ccall( (:function, "library"), return_type, (tuple_of_argument_types,), arguments)

Arguments to C functions are, in general, automatically converted, and the returned values in C types are also converted to Julia types. Arrays of Booleans are handled differently in C and Julia and cannot be passed directly, so they must be manually converted. The ccall function will also automatically ensure that all of its arguments will be preserved from garbage collection until the call returns.

Calling Fortran functions is similar. Just note that all inputs must be passed by reference. The support for C++ is more complex, see packages like Cxx and CxxWrap.

Python

Calling Python from Julia is done easily using the PyCall package. With the @pyimport macro, one can easily import any Python library, whose functions are called with the familiar dot notation.

Parallel Computing

Julia’s model for building a large parallel application works by means of a global distributed address space. This means that you can hold a reference to an object that lives on another machine participating in a computation. These references are easily manipulated and passed around between machines, making it simple to keep track of what’s being computed where. Also, machines can be added in mid computation when needed.

Besides manually receiving and sending messages between processes, or workers, Julia also provides some useful wrappers.

Parallel macro

The @parallel macro acts on a for loop, splitting the range, and distributing it to each process. It optionally takes a “reducer” as its first argument. If a reducer is specified, the results from each remote procedure will be aggregated using the reducer.

If the computational task is to apply a function to all elements in some collection, the parallel map operation can be applied using the pmap function pmap(f, coll). It applies a function on each element of a collection in parallel, while preserving the order of the collection in the result.

Distributed array

When computations have to be done on a very large array (or arrays), the array can be distributed, so that each process works in parallel on a different part of the array. In this way, we can make use of the memory resources of multiple machines, and allow the manipulation of arrays that would be too large to fit on one machine.

The specific data type used here is called a distributed array or DArray; most operations behave exactly as on the normal Array type, so the parallelism is invisible. With DArray, each process has local access to just a part of the data, and no two processes share the same data.

More Topics for Later

There are yet many more topics to talk about in Julia, e.g.

1. Packages: Julia has a built-in package system with git capability.
2. Graphics: Plotting capability is typically enabled by interfacing to external packages, e.g. Matplotlib from Python.
3. Benchmark: There are some tools/packages for benchmarking, which are critical for developing real performant Julia code.
4. Performance notes: The whole book “Julia High Performance” by Avik Sengupta is devoted to this topic.

These topics (and more) will probably be discussed in future posts.