The Julia Blog

Passing Julia Callback Functions to C

2013-05-10T00:00:00-07:00

One of the great strengths of Julia is that it is so easy to call C code natively, with no special “glue” routines or overhead to marshal arguments and convert return values. For example, if you want to call GNU GSL to compute a special function like a Debye integral, it is as easy as:

debye_1(x) = ccall((:gsl_sf_debye_1,:libgsl), Cdouble, (Cdouble,), x)

at which point you can compute debye_1(2), debye_1(3.7), and so on. (Even easier would be to use Jiahao Chen’s GSL package for Julia, which has already created such wrappers for you.) This makes a vast array of existing C libraries accessible to you in Julia (along with Fortran libraries and other languages with C-accessible calling conventions).

In fact, you can even go the other way around, passing Julia routines to C, so that C code is calling Julia code in the form of callback functions. For example, a C library for numerical integration might expect you to pass the integrand as a function argument, which the library will then call to evaluate the integrand as many times as needed to estimate the integral. Callback functions are also natural for optimization, root-finding, and many other numerical tasks, as well as in many non-numerical problems. The purpose of this blog post is to illustrate the techniques for passing Julia functions as callbacks to C routines, which is straightforward and efficient but requires some lower-level understanding of how functions and other values are passed as arguments.

The code in this post requires Julia 0.2 (or a recent git facsimile thereof); the key features needed for callback functions (especially unsafe_pointer_to_objref) are not available in Julia 0.1.

Sorting with `qsort`

Perhaps the most well-known example of a callback parameter is provided by the qsort function, part of the ANSI C standard library and declared in C as:

void qsort(void *base, size_t nmemb, size_t size,
           int(*compare)(const void *a, const void *b));

The base argument is a pointer to an array of length nmemb, with elements of size bytes each. compare is a callback function which takes pointers to two elements a and b and returns an integer less/greater than zero if a should appear before/after b (or zero if any order is permitted). Now, suppose that we have a 1d array A of values in Julia that we want to sort using the qsort function (rather than Julia’s built-in sort function). Before we worry about calling qsort and passing arguments, we need to write a comparison function that works for some arbitrary type T, e.g.

function mycompare{T}(a_::Ptr{T}, b_::Ptr{T})
    a = unsafe_load(a_)
    b = unsafe_load(b_)
    return a < b ? cint(-1) : a > b ? cint(+1) : cint(0)
end
cint(n) = convert(Cint, n)

Notice that we use the built-in function unsafe_load to fetch the values pointed to by the arguments a_ and b_ (which is “unsafe” because it will crash if these are not valid pointers, but qsort will always pass valid pointers). Also, we have to be a little careful about return values: qsort expects a function returning a C int, so we must be sure to return Cint (the corresponding type in Julia) via a call to convert.

Now, how do we pass this to C? A function pointer in C is essentially just a pointer to the memory location of the machine code implementing that function, whereas a function value mycompare (of type Function) in Julia is quite different. Thanks to Julia’s JIT compilation approach,a Julia function may not even be compiled until the first time it is called, and in general the same Julia function may be compiled into multiple machine-code instantiations, which are specialized for arguments of different types (e.g. different T in this case). So, you can imagine that mycompare must internally point to a rather complicated data structure (a jl_function_t in julia.h, if you are interested), which holds information about the argument types, the compiled versions (if any), and so on. In general, it must store a closure with information about the environment in which the function was defined; we will talk more about this below. In any case, it is a very different object than a simple pointer to machine code for one set of argument types. Fortunately, we can get the latter simply by calling a built-in Julia function called cfunction:

const mycompare_c = cfunction(mycompare, Cint, (Ptr{Cdouble}, Ptr{Cdouble}))

Here, we pass cfunction three arguments: the function mycompare, the return type Cint, and a tuple of the argument types, in this case to sort an array of Cdouble (Float64) elements. Julia compiles a version of mycompare specialized for these argument types (if it has not done so already), and returns a Ptr{Void} holding the address of the machine code, exactly what we need to pass to qsort. We are now ready to call qsort on some sample data:

A = [1.3, -2.7, 4.4, 3.1]
ccall(:qsort, Void, (Ptr{Cdouble}, Csize_t, Csize_t, Ptr{Void}),
      A, length(A), sizeof(eltype(A)), mycompare_c)

After this executes, A is changed to the sorted array [ -2.7, 1.3, 3.1, 4.4]. Note that Julia knows how to convert an array A::Vector{Cdouble} into a Ptr{Cdouble}, how to compute the sizeof a type in bytes (identical to C’s sizeof operator), and so on. For fun, try inserting a println("mycompare($a,$b)") line into mycompare, which will allow you to see the comparisons that qsort is performing (and to verify that it is really calling the Julia function that you passed to it).

The problem with closures

We aren’t done yet, however. If you start passing callback functions to C routines, it won’t be long before you discover that cfunction doesn’t always work. For example, suppose we tried to declare our comparison function inline, via:

mycomp = cfunction((a_,b_) -> unsafe_load(a_) < unsafe_load(b_) ? 
                              cint(-1) : cint(+1),
                   Cint, (Ptr{Cdouble}, Ptr{Cdouble}))

Julia barfs on this, printing ERROR: function is not yet c-callable. In general, cfunction only works for “top-level” functions: named functions defined in the top-level (global or module) scope, but not anonymous (args -> value) functions and not functions defined within other functions (“nested” functions). The reason for this stems from one important concept in computer science: a closure.

To understand the need for closures, and the difficulty they pose for callback functions, suppose that we wanted to provide a nicer interface for qsort, one which permitted the user to simply pass a lessthan function returning true or false while hiding all of the low-level business with pointers, Cint, and so on. We might like to do something of the form:

function qsort!{T}(A::Vector{T}, lessthan::Function)
    function mycompare(a_::Ptr{T}, b_::Ptr{T})
        a = unsafe_load(a_)
        b = unsafe_load(b_)
        return lessthan(a, b) ? cint(-1) : cint(+1)
    end
    mycompare_c = cfunction(mycompare, Cint, (Ptr{T}, Ptr{T}))
    ccall(:qsort, Void, (Ptr{T}, Csize_t, Csize_t, Ptr{Void}),
          A, length(A), sizeof(T), mycompare_c)
    A
end

Then we could simply call qsort!([1.3, -2.7, 4.4, 3.1], <) to sort in ascending order using the built-in < comparison, or any other comparison function we wanted. Unfortunately cfunction will again barf when you try to call qsort!, and it is no longer so difficult to understand why. Notice that the nested mycompare function is no longer self-contained: it uses the variable lessthan from the surrounding scope. This is a common pattern for nested functions and anonymous functions: often, they are parameterized by local variables in the environment where the function is defined. Technically, the ability to have this kind of dependency is provided by lexical scoping in a programming language like Julia, and is typical of any language in which functions are “first-class” objects. In order to support lexical scoping, a Julia Function object needs to internally carry around a pointer to the variables in the enclosing environment, and this encapsulation is called a closure.

In contrast, a C function pointer is not a closure. It doesn’t enclose a pointer to the environment in which the function was defined, or anything else for that matter; it is just the address of a stream of instructions. This makes it hard, in C, to write functions to transform other functions (higher-order functions) or to parameterize functions by local variables. This apparently leaves us with two options, neither of which is especially attractive:

We could store lessthan in a global variable, and reference that from a top-level mycompare function. (This is the traditional solution for C programmers calling qsort with parameterized comparison functions.) The problem with this strategy is that it is not re-entrant: it prevents us from calling qsort! recursively (e.g. if the comparison function itself needs to do a sort, for some complicated datastructure), or from calling qsort! from multiple threads (when a future Julia version supports shared-memory parallelism). Still, this is better than nothing.
Every time qsort! is called, Julia could JIT-compile a new version of mycompare, which hard-codes the reference to the lessthan argument passed on that call. This is technically possible and has been implemented in some languages (e.g. reportedly GNU Guile and Lua do something like this). However, this strategy comes at a price: it requires that callbacks be recompiled every time a parameter in them changes, which is not true of the global-variable strategy. Anyway, it is not implemented yet in Julia.

Fortunately, there is often a third option, because C programmers long ago recognized these limitations of function pointers, and devised a workaround: most modern C callback interfaces allow arbitrary data to be passed through to the callback via a “pass-through” (or “thunk”) pointer parameter. As explained in the next section, we can exploit this technique in Julia to pass a “true” closure as a callback.

Passing closures via pass-through pointers

The qsort interface is nowadays considered rather antiquated. Years ago, it was supplemented on BSD-Unix systems, and eventually in GNU libc, by a function called qsort_r that solves the problem of passing parameters to the callback in a re-entrant way. This is how the BSD (e.g. MacOS) qsort_r function is defined:

void qsort_r(void *base, size_t nmemb, size_t size, void *thunk,
             int (*compare)(void *thunk, const void *a, const void *b));

Compared to qsort, there is an extra thunk parameter, and this is passed through to the compare function as its first argument. In this way, you can pass a pointer to arbitrary data through to your callback, and we can exploit this to pass a closure through for an arbitrary Julia callback.

All we need is a way to convert a Julia Function into an opaque Ptr{Void} so that we can pass it through to our callback, and then a way to convert the opaque pointer back into a Function. The former is automatic if we simply declare the ccall argument as type Any (which passes the argument as an opaque Julia object pointer), and the latter is accomplished by the built-in function unsafe_pointer_to_objref. (Technically, we could use type Function or an explicit call to pointer_from_objref instead of Any.) Using these, we can now define a working high-level qsort! function that takes an arbitrary lessthan comparison-function argument:

function qsort!_compare{T}(lessthan_::Ptr{Void}, a_::Ptr{T}, b_::Ptr{T})
    a = unsafe_load(a_)
    b = unsafe_load(b_)
    lessthan = unsafe_pointer_to_objref(lessthan_)::Function
    return lessthan(a, b) ? cint(-1) : cint(+1)
end

function qsort!{T}(A::Vector{T}, lessthan::Function=<)
    compare_c = cfunction(qsort!_compare, Cint, (Ptr{Void}, Ptr{T}, Ptr{T}))
    ccall(:qsort_r, Void, (Ptr{T}, Csize_t, Csize_t, Any, Ptr{Void}),
          A, length(A), sizeof(T), lessthan, compare_c)
    return A
end

qsort!_compare is a top-level function, so cfunction has no problem with it, and it will only be compiled once per type T to be sorted (rather than once per call to qsort! or per lessthan function). We use the explicit ::Function assertion to tell the compiler that we will only pass Function pointers in lessthan_. Note that we gave the lessthan argument a default value of < (default arguments being a recent feature added to Julia).

We can now do qsort!([1.3, -2.7, 4.4, 3.1]) and it will return the array sorted in ascending order, or qsort!([1.3, -2.7, 4.4, 3.1], >) to sort in descending order.

Warning: `qsort_r` is not portable

The example above has one major problem that has nothing to do with Julia: the qsort_r function is not portable. The above example won’t work on Windows, since the Windows C library doesn’t define qsort_r (instead, it has a function called qsort_s, which of course uses an argument order incompatible with both the BSD and GNU qsort_r functions). Worse, it will crash on GNU/Linux systems, which do provide qsort_r but with an incompatible calling convention. And as a result it is difficult to use qsort_r in a way that does not crash either on GNU/Linux or BSD (e.g. MacOS) systems. This is how glibc’s qsort_r is defined:

void qsort_r(void *base, size_t nmemb, size_t size,
             int (*compare)(const void *a, const void *b, void *thunk),
              void *thunk);

Note that the position of the thunk argument is moved, both in qsort_r itself and in the comparison function. So, the corresponding qsort! Julia code on GNU/Linux systems should be:

function qsort!_compare{T}(a_::Ptr{T}, b_::Ptr{T}, lessthan_::Ptr{Void})
    a = unsafe_load(a_)
    b = unsafe_load(b_)
    lessthan = unsafe_pointer_to_objref(lessthan_)::Function
    return lessthan(a, b) ? cint(-1) : cint(+1)
end

function qsort!{T}(A::Vector{T}, lessthan::Function=<)
    compare_c = cfunction(qsort!_compare, Cint, (Ptr{T}, Ptr{T}, Ptr{Void}))
    ccall(:qsort_r, Void, (Ptr{T}, Csize_t, Csize_t, Ptr{Void}, Any),
          A, length(A), sizeof(T), compare_c, lessthan)
    return A
end

If you really needed to call qsort_r from Julia, you could use the above definitions if OS_NAME == :Linux and the BSD definitions otherwise, with a third version using qsort_s on Windows, but fortunately there is not much need as Julia comes with its own perfectly adequate sort and sort! routines.

Passing closures in data structures

As another example that is oriented more towards numerical computations, we’ll examine how we might call the numerical integration routines in the GNU Scientific Library (GSL). There is already a GSL package that handles the wrapper work below for you, but it is instructive to look at how this is implemented because GSL simulates closures in a slightly different way, with data structures.

Like most modern C libraries accepting callbacks, GSL uses a void* pass-through parameter to allow arbitrary data to be passed through to the callback routine, and we can use that to support arbitrary closures in Julia. Unlike qsort_r, however, GSL wraps both the C function pointer and the pass-through pointer in a data structure called gsl_function:

struct {
    double (*function)(double x, void *params);
    void *params;
} gsl_function;

Using the techniques above, we can easily declare a GSL_Function type in Julia that mirrors this C type, and with a constructor GSL_Function(f::Function) that creates a wrapper around an arbitrary Julia function f:

function gsl_function_wrap(x::Cdouble, params::Ptr{Void})
    f = unsafe_pointer_to_objref(params)::Function
    convert(Cdouble, f(x))::Cdouble
end
const gsl_function_wrap_c = cfunction(gsl_function_wrap,
                                      Cdouble, (Cdouble, Ptr{Void}))

type GSL_Function
    func::Ptr{Void}
    params::Any
    GSL_Function(f::Function) = new(gsl_function_wrap_c, f)
end

One subtlety with the above code is that we need to explicitly convert the return value of f to a Cdouble (in case the caller’s code returns some other numeric type for some x, such as an Int). Moreover, we need to explicitly assert (::Cdouble) that the result of the convert was a Cdouble. As with the qsort example, this is because cfunction only works if Julia can guarantee that gsl_function_wrap returns the specified Cdouble type, and Julia cannot infer the return type of convert since it does not know the return type of f(x).

Given the above definitions, it is a simple matter to pass this to the GSL adaptive-integration routines in a wrapper function gsl_integration_qag:

function gsl_integration_qag(f::Function, a::Real, b::Real, epsrel::Real=1e-12,
                             maxintervals::Integer=10^7)
    s = ccall((:gsl_integration_workspace_alloc,:libgsl), Ptr{Void}, (Csize_t,),
              maxintervals)
    result = Array(Cdouble,1)
    abserr = Array(Cdouble,1)
    ccall((:gsl_integration_qag,:libgsl), Cint,
          (Ptr{GSL_Function}, Cdouble,Cdouble, Cdouble, Csize_t, Cint, Ptr{Void}, 
           Ptr{Cdouble}, Ptr{Cdouble}),
          &GSL_Function(f), a, b, epsrel, maxintervals, 1, s, result, abserr)
    ccall((:gsl_integration_workspace_free,:libgsl), Void, (Ptr{Void},), s)
    return (result[1], abserr[1])
end

Note that &GSL_Function(f) passes a pointer to a GSL_Function “struct” containing a pointer to gsl_function_wrap_c and f, corresponding to the gsl_function* argument in C. The return value is a tuple of the estimated integral and an estimated error.

For example, gsl_integration_qag(cos, 0, 1) returns (0.8414709848078965,9.34220461887732e-15), which computes the correct integral sin(1) to machine precision.

Taking out the trash (or not)

In the above examples, we pass an opaque pointer (object reference) to a Julia Function into C. Whenever one passes pointers to Julia data into C code, one has to ensure that the Julia data is not garbage-collected until the C code is done with it, and functions are no exception to this rule. An anonymous function that is no longer referred to by any Julia variable may be garbage collected, at which point any C pointers to it become invalid.

This sounds scary, but in practice you don’t need to worry about it very often, because Julia guarantees that ccall arguments won’t be garbage-collected until the ccall exits. So, in all of the above examples, we are safe: the Function only needs to live as long as the ccall.

The only danger arises when you pass a function pointer to C and the C code saves the pointer in some data structure which it will use in a later ccall. In that case, you are responsible for ensuring that the Function variable lives (is referred to by some Julia variable) as long as the C code might need it.

For example, in the GSL one-dimensional minimization interface, you don’t simply pass your objective function to a minimization routine and wait until it is minimized. Instead, you call a GSL routine to create a “minimizer object”, store your function pointer in this object, call routines to iterate the minimization, and then deallocate the minimizer when you are done. The Julia function must not be garbage-collected until this process is complete. The easiest way to ensure this is to create a Julia wrapper type around the minimizer object that stores an explicit reference to the Julia function, like this:

type GSL_Minimizer
    m::Ptr{Void} # the gsl_min_fminimizer pointer
    f::Any  # explicit reference to objective, to prevent garbage-collection
    function GSL_Minimizer(t)
       m = ccall((:gsl_min_fminimizer_alloc,:libgsl), Ptr{Void}, (Ptr{Void},), t)
       p = new(m, nothing)
       finalizer(p, p -> ccall((:gsl_min_fminimizer_free,:libgsl),
                               Void, (Ptr{Void},), p.m))
       p
    end
end

This wraps around a gsl_min_fminimizer object of type t, with a placeholder f to store a reference to the objective function (once it is set below), including a finalizer to deallocate the GSL object when the GSL_Minimizer is garbage-collected. The parameter t is used to specify the minimization algorithm, which could default to Brent’s algorithm via:

const gsl_brent = unsafe_load(cglobal((:gsl_min_fminimizer_brent,:libgsl), Ptr{Void}))
GSL_Minimizer() = GSL_Minimizer(gsl_brent)

(The call to cglobal yields a pointer to the gsl_min_fminimizer_brent global variable in GSL, which we then dereference to get the actual pointer via unsafe_load.)

Then, when we set the function to minimize (the “objective”), we store an extra reference to it in the GSL_Minimizer to prevent garbage-collection for the lifetime of the GSL_Minimizer, again using the GSL_Function type defined above to wrap the callback:

function gsl_minimizer_set!(m::GSL_Minimizer, f, x0, xmin, xmax)
    ccall((:gsl_min_fminimizer_set,:libgsl), Cint,
          (Ptr{Void}, Ptr{GSL_Function}, Cdouble, Cdouble, Cdouble),
          m.m, &GSL_Function(f), x0, xmin, xmax)
    m.f = f
    m
end

There are then various GSL routines to iterate the minimizer and to check the current x, objective value, or bounds on the minimum, which are convenient to wrap:

gsl_minimizer_iterate!(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_iterate,:libgsl), Cint, (Ptr{Void},), m.m)

gsl_minimizer_x(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_x_minimum,:libgsl), Cdouble, (Ptr{Void},), m.m)

gsl_minimizer_f(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_f_minimum,:libgsl), Cdouble, (Ptr{Void},), m.m)

gsl_minimizer_xmin(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_x_lower,:libgsl), Cdouble, (Ptr{Void},), m.m)
gsl_minimizer_xmax(m::GSL_Minimizer) =
    ccall((:gsl_min_fminimizer_x_upper,:libgsl), Cdouble, (Ptr{Void},), m.m)

Putting all of these together, we can minimize a simple function sin(x) in the interval [-3,1], with a starting guess -1, via:

m = GSL_Minimizer()
gsl_minimizer_set!(m, sin, -1, -3, 1)
while gsl_minimizer_xmax(m) - gsl_minimizer_xmin(m) > 1e-6
    println("iterating at x = $(gsl_minimizer_x(m))")
    gsl_minimizer_iterate!(m)
end
println("found minimum $(gsl_minimizer_f(m)) at x = $(gsl_minimizer_x(m))")

After a few iterations, it prints found minimum -1.0 at x = -1.5707963269964016, which is the correct minimum (−π/2) to about 10 digits.

At this point, I will shamelessly plug my own NLopt package for Julia, which wraps around my free/open-source NLopt library to provide many more optimization algorithms than GSL, with perhaps a nicer interface. However, the techniques used to pass callback functions to NLopt are actually quite similar to those used for GSL.

An even more complicated version of these techniques can be found in the PyCall package to call Python from Julia. In order to pass a Julia function to Python, we again use cfunction on a wrapper function that handles the type conversions and so on, and pass the actual Julia closure through via a pass-through pointer. But in that case, the pass-through pointer consists of a Python object that has been created with a new type that allows it to wrap a Julia object, and garbage-collection is deferred by storing the Julia object in a global dictionary of saved objects (removing it via the Python destructor of the new type). That is all somewhat tricky stuff and beyond the scope of this blog post; I only mention it to illustrate the fact that it is possible to implement quite complex inter-language calling behaviors purely in Julia by building on the above techniques.

Put This In Your Pipe

2013-04-08T00:00:00-07:00

In a previous post, I talked about why “shelling out” to spawn a pipeline of external programs via an intermediate shell is a common cause of bugs, security holes, unnecessary overhead, and silent failures. But it’s so convenient! Why can’t running pipelines of external programs be convenient and safe? Well, there’s no real reason, actually. The shell itself manages to construct and execute pipelines quite well. In principle, there’s nothing stopping high-level languages from doing it at least as well as shells do – the common ones just don’t by default, instead requiring users to make the extra effort to use external programs safely and correctly. There are two major impediments:

Some moderately tricky low-level UNIX plumbing using the pipe, dup2, fork, close, and exec system calls;
The UX problem of designing an easy, flexible programming interface for commands and pipelines.

This post describes the system we designed and implemented for Julia, and how it avoids the major flaws of shelling out in other languages. First, I’ll present the Julia version of the previous post’s example – counting the number of lines in a given directory containing the string “foo”. The fact that Julia provides complete, specific diagnostic error messages when pipelines fail turns out to reveal a surprising and subtle bug, lurking in what appears to be a perfectly innocuous UNIX pipeline. After fixing this bug, we go into details of how Julia’s external command execution and pipeline construction system actually works, and why it provides greater flexibility and safety than the traditional approach of using an intermediate shell to do all the heavy lifting.

Simple Pipeline, Subtle Bug

Here’s how you write the example of counting the number of lines in a directory containing the string “foo” in Julia (you can follow along at home if you have Julia installed from source by changing directories into the Julia source directory and doing cp -a src "source code"; mkdir tmp and then firing up the Julia repl):

julia> dir = "src";

julia> int(readchomp(`find $dir -type f -print0` | `xargs -0 grep foo` | `wc -l`))
5

This Julia command looks suspiciously similar to the naïve Ruby version we started with in the previous post:

`find #{dir} -type f -print0 | xargs -0 grep foo | wc -l`.to_i

However, it isn’t susceptible to the same problems:

julia> dir = "source code";

julia> int(readchomp(`find $dir -type f -print0` | `xargs -0 grep foo` | `wc -l`))
5

julia> dir = "nonexistent";

julia> int(readchomp(`find $dir -type f -print0` | `xargs -0 grep foo` | `wc -l`))
find: `nonexistent': No such file or directory
ERROR: failed processes:
  Process(`find nonexistent -type f -print0`, ProcessExited(1)) [1]
  Process(`xargs -0 grep foo`, ProcessExited(123)) [123]
 in pipeline_error at process.jl:412
 in readall at process.jl:365
 in readchomp at io.jl:172

julia> dir = "foo'; echo MALICIOUS ATTACK; echo '";

julia> int(readchomp(`find $dir -type f -print0` | `xargs -0 grep foo` | `wc -l`))
find: `foo\'; echo MALICIOUS ATTACK; echo \'': No such file or directory
ERROR: failed processes:
  Process(`find "foo'; echo MALICIOUS ATTACK; echo '" -type f -print0`, ProcessExited(1)) [1]
  Process(`xargs -0 grep foo`, ProcessExited(123)) [123]
 in pipeline_error at process.jl:412
 in readall at process.jl:365
 in readchomp at io.jl:172

The default, simplest-to-achieve behavior in Julia is:

not susceptible to any kind of metacharacter breakage,
reliably detects all subprocess failures,
automatically raises an exception if any subprocess fails,
prints error messages including exactly which commands failed.

In the above examples, we can see that even when dir contains spaces or quotes, the expression still behaves exactly as intended – the value of dir is interpolated as a single argument to the find command. When dir is not the name of a directory that exists, find fails – as it should – and this failure is detected and automatically converted into an informative exception, including the fully expanded command-lines that failed.

In the previous post, we observed that using the pipefail option for Bash allows detection of pipeline failures, like this one, occurring before the last process in the pipeline. However, it only allows us to detect that at least one thing in the pipeline failed. We still have to guess at what parts of the pipeline actually failed. In the Julia example, on the other hand, there is no guessing required: when a non-existent directory is given, we can see that both find and xargs fail. While it is unsurprising that find fails in this case, it is unexpected that xargs also fails. Why does xargs fail?

One possibility to check for is that the xargs program fails with no input. We can use Julia’s success predicate to try it out:

julia> success(`cat /dev/null` | `xargs true`)
true

Ok, so xargs seems perfectly happy with no input. Maybe grep doesn’t like not getting any input?

julia> success(`cat /dev/null` | `grep foo`)
false

Aha! grep returns a non-zero status when it doesn’t get any input. Good to know. It turns out that grep indicates whether it matched anything or not with its return status. Most programs use their return status to indicate success or failure, but some, like grep, use it to indicate some other boolean condition – in this case “found something” versus “didn’t find anything”:

julia> success(`echo foo` | `grep foo`)
true

julia> success(`echo bar` | `grep foo`)
false

Now we know why grep is “failing” – and xargs too, since it returns a non-zero status if the program it runs returns non-zero. This means that our Julia pipeline and the “responsible” Ruby version are both susceptible to bogus failures when we search an existing directory that happens not to contain the string “foo” anywhere:

julia> dir = "tmp";

julia> int(readchomp(`find $dir -type f -print0` | `xargs -0 grep foo` | `wc -l`))
ERROR: failed process: Process(`xargs -0 grep foo`, ProcessExited(123)) [123]
 in error at error.jl:22
 in pipeline_error at process.jl:394
 in pipeline_error at process.jl:407
 in readall at process.jl:365
 in readchomp at io.jl:172

Since grep indicates not finding anything using a non-zero return status, the readall function concludes that its pipeline failed and raises an error to that effect. In this case, this default behavior is undesirable: we want the expression to just return 0 without raising an error. The simple fix in Julia is this:

julia> dir = "tmp";

julia> int(readchomp(`find $dir -type f -print0` | ignorestatus(`xargs -0 grep foo`) | `wc -l`))
0

This works correctly in all cases. Next I’ll explain how all of this works, but for now it’s enough to note that the detailed error message provided when our pipeline failed exposed a rather subtle bug that would eventually cause subtle and hard-to-debug problems when used in production. Without such detailed error reporting, this bug would be pretty difficult to track down.

Do-Nothing Backticks

Julia borrows the backtick syntax for external commands form Perl and Ruby, both of which in turn got it from the shell. Unlike in these predecessors, however, in Julia backticks don’t immediately run commands, nor do they necessarily indicate that you want to capture the output of the command. Instead, backticks just construct an object representing a command:

julia> `echo Hello`
`echo Hello`

julia> typeof(ans)
Cmd

(In the Julia repl, ans is automatically bound to the value of the last evaluated input.) In order to actually run a command, you have to do something with a command object. To run a command and capture its output into a string – what other languages do with backticks automatically – you can apply the readall function:

julia> readall(`echo Hello`)
"Hello\n"

Since it’s very common to want to discard the trailing line break at the end of a command’s output, Julia provides the readchomp(x) command which is equivalent to writing chomp(readall(x)):

julia> readchomp(`echo Hello`)
"Hello"

To run a command without capturing its output, letting it just print to the same stdout stream as the main process – i.e. what the system function does when given a command as a string in other languages – use the run function:

julia> run(`echo Hello`)
Hello

The "Hello\n" after the readall command is a returned value, whereas the Hello after the run command is printed output. (If your terminal supports color, these are colored differently so that you can easily distinguish them visually.) Nothing is returned by the run command, but if something goes wrong, an exception is raised:

julia> run(`false`)
ERROR: failed process: Process(`false`, ProcessExited(1)) [1]
 in error at error.jl:22
 in pipeline_error at process.jl:394
 in run at process.jl:384

julia> run(`notaprogram`)
execvp(): No such file or directory
ERROR: failed process: Process(`notaprogram`, ProcessExited(-1)) [-1]
 in error at error.jl:22
 in pipeline_error at process.jl:394
 in run at process.jl:384

As with xargs and grep above, this may not always be desirable. In such cases, you can use ignorestatus to indicate that the command returning a non-zero value should not be considered an error:

julia> run(ignorestatus(`false`))

Although Julia’s backtick syntax intentionally mimics the shell as closely as possible, there is an important distinction: the command string is never passed to a shell to be interpreted and executed; instead it is parsed in Julia code, using the same rules the shell uses to determine what the command and arguments are. Command objects allow you to see what the program and arguments were determined to be by accessing the .exec field:

julia> cmd = `perl -e 'print "Hello\n"'`
`perl -e 'print "Hello\n"'`

julia> cmd.exec
3-element Union(UTF8String,ASCIIString) Array:
 "perl"
 "-e"
 "print \"Hello\\n\""

This field is a plain old array of strings that can be manipulated like any other Julia array.

Constructing Commands

The purpose of the backtick notation in Julia is to provide a familiar, shell-like syntax for making objects representing commands with arguments. To that end, quotes and spaces work just as they do in the shell. The real power of backtick syntax doesn’t emerge, however, until we begin constructing commands programmatically. Just as in the shell (and in Julia strings), you can interpolate values into commands using the dollar sign ($):

julia> dir = "src";

julia> `find $dir -type f`.exec
4-element Union(UTF8String,ASCIIString) Array:
 "find"
 "src"
 "-type"
 "f"

Unlike in the shell, however, Julia values interpolated into commands are interpolated as a single verbatim argument – no characters inside the value are interpreted as special after the value has been interpolated:

julia> dir = "two words";

julia> `find $dir -type f`.exec
4-element Union(UTF8String,ASCIIString) Array:
 "find"
 "two words"
 "-type"
 "f"

julia> dir = "foo'bar";

julia> `find $dir -type f`.exec
4-element Union(UTF8String,ASCIIString) Array:
 "find"
 "foo'bar"
 "-type"
 "f"

This works no matter what the contents of the interpolated value is, allowing simple interpolation of characters that are quite difficult to pass as parts of command-line arguments even in the shell (for the following examples, tmp/a.tsv and tmp/b.tsv can be created in the shell with echo -e "foo\tbar\nbaz\tqux" > tmp/a.tsv; echo -e "foo\t1\nbaz\t2" > tmp/b.tsv):

julia> tab = "\t";

julia> cmd = `join -t$tab tmp/a.tsv tmp/b.tsv`;

julia> cmd.exec
4-element Union(UTF8String,ASCIIString) Array:
 "join"
 "-t\t"
 "tmp/a.tsv"
 "tmp/b.tsv"

julia> run(cmd)
foo     bar     1
baz     qux     2

Moreover, what comes after the $ can actually be any valid Julia expression, not just a variable name:

julia> `join -t$"\t" tmp/a.tsv tmp/b.tsv`.exec
4-element Union(UTF8String,ASCIIString) Array:
 "join"
 "-t\t"
 "a.tsv"
 "b.tsv"

A tab character is somewhat harder to pass in the shell, requiring command interpolation and some tricky quoting:

bash-3.2$ join -t"$(printf '\t')" tmp/a.tsv tmp/b.tsv
foo     bar     1
baz     qux     2

While interpolating values with spaces and other strange characters is great for non-brittle construction of commands, there was a reason why the shell split values on spaces in the first place: to allow interpolation of multiple arguments. Most modern shells have first-class array types, but older shells used space-separation to simulate arrays. Thus, if you interpolate a value like “foo bar” into a command in the shell, it’s treated as two separate words by default. In languages with first-class array types, however, there’s a much better option: consistently interpolate single values as single arguments and interpolate arrays as multiple values. This is precisely what Julia’s backtick interpolation does:

julia> dirs = ["foo", "bar", "baz"];

julia> `find $dirs -type f`.exec
6-element Union(UTF8String,ASCIIString) Array:
 "find"
 "foo"
 "bar"
 "baz"
 "-type"
 "f"

And of course, no matter how strange the strings contained in an interpolated array are, they become verbatim arguments, without any shell interpretation. Julia’s backticks have one more fancy trick up their sleeve. We saw earlier (without really remarking on it) that you could interpolate single values into a larger argument:

julia> x = "bar";

julia> `echo foo$x`
`echo foobar`

What happens if x is an array? Only one way to find out:

julia> x = ["bar", "baz"];

julia> `echo foo$x`
`echo foobar foobaz`

Julia does what the shell would do if you wrote echo foo{bar,baz}. This even works correctly for multiple values interpolated into the same shell word:

julia> dir = "/data"; names = ["foo","bar"]; exts=["csv","tsv"];

julia> `cat $dir/$names.$exts`
`cat /data/foo.csv /data/foo.tsv /data/bar.csv /data/bar.tsv`

This is the same Cartesian product expansion that the shell does if multiple {...} expressions are used in the same word.

Distributed Numerical Optimization

2013-04-05T00:00:00-07:00

This post walks through the parallel computing functionality of Julia to implement an asynchronous parallel version of the classical cutting-plane algorithm for convex (nonsmooth) optimization, demonstrating the complete workflow including running on both Amazon EC2 and a large multicore server. I will quickly review the cutting-plane algorithm and will be focusing primarily on parallel computation patterns, so don’t worry if you’re not familiar with the optimization side of things.

Cutting-plane algorithm

The cutting-plane algorithm is a method for solving the optimization problem

$$\text{minimize} _ {x \in \mathbb{R}^d} \sum_{i=1}^n f_i(x) $$

where the functions $ f_i $ are convex but not necessarily differentiable. The absolute value function $ |x| $ and the 1-norm $ ||x|| _ 1 $ are typical examples. Important applications also arise from Lagrangian relaxation. The idea of the algorithm is to approximate the functions $ f_i $ with piecewise linear models $ m_i $ which are built up from information obtained by evaluating $ f_i $ at different points. We iteratively minimize over the models to generate candidate solution points.

We can state the algorithm as

Choose starting point $ x $.
For $i = 1,\ldots,n$, evaluate $ f_i(x) $ and update corresponding model $ m_i $.
Let the next candidate $ x $ be the minimizer of $ \sum_{i=1}^n m_i(x) $.
If not converged, goto step 2.

If it is costly to evaluate $ f_i(x) $, then the algorithm is naturally parallelizable at step 2. The minimization in step 3 can be computed by solving a linear optimization problem, which is usually very fast. (Let me point out here that Julia has interfaces to the linear programming solvers Clp, Gurobi, and GLPK.)

Abstracting the math, we can write the algorithm using the following Julia code.

# functions initialize, isconverged, solvesubproblem, and process implemented elsewhere
state, subproblems = initialize()
while !isconverged(state)
    results = map(solvesubproblem,subproblems)
    state, subproblems = process(state, results)
end

The function solvesubproblem corresponds to evaluating $ f_i(x) $ for a given $ i $ and $ x $ (the elements of subproblems could be tuples (i,x)). The function process corresponds to minimizing the model in step 3, and it produces a new state and a new set of subproblems to solve.

Note that the algorithm looks much like a map-reduce that would be easy to parallelize using many existing frameworks. Indeed, in Julia we can simply replace map with pmap (parallel map). Let’s consider a twist that makes the parallelism not so straightforward.

Asynchronous variant

Variability in the time taken by the solvesubproblem function can lead to load imbalance and limit parallel efficiency as workers sit idle waiting for new tasks. Such variability arises naturally if solvesubproblem itself requires solving a optimization problem, or if the workers and network are shared, as is often the case with cloud computing.

We can consider a new variant of the cutting-plane algorithm to address this issue. The key point is

When proportion $0 < \alpha \le 1 $ of subproblems for a given candidate have been solved, generate a new candidate and corresponding set of subproblems by using whatever information is presently available.

In other words, we generate new tasks to feed to workers without needing to wait for all current tasks to complete, making the algorithm asynchronous. The algorithm remains convergent, although the total number of iterations may increase. For more details, see this paper by Jeff Linderoth and Stephen Wright.

By introducing asynchronicity we can no longer use a nice black-box pmap function and have to dig deeper into the parallel implementation. Fortunately, this is easy to do in Julia.

Parallel implementation in Julia

Julia implements distributed-memory parallelism based on one-sided message passing, where process push work onto others (via remotecall) and the results are retrieved (via fetch) by the process which requires them. Macros such as @spawn and @parallel provide pretty syntax around this low-level functionality. This model of parallelism is very different from the typical SIMD style of MPI. Both approaches are useful in different contexts, and I expect an MPI wrapper for Julia will appear in the future (see also here).

Reading the manual on parallel computing is highly recommended, and I won’t try to reproduce it in this post. Instead, we’ll dig into and extend one of the examples it presents.

The implementation of pmap in Julia is

function pmap(f, lst)
    np = nprocs()  # determine the number of processors available
    n = length(lst)
    results = cell(n)
    i = 1
    # function to produce the next work item from the queue.
    # in this case it's just an index.
    next_idx() = (idx=i; i+=1; idx)
    @sync begin
        for p=1:np
            if p != myid() || np == 1
                @spawnlocal begin
                    while true
                        idx = next_idx()
                        if idx > n
                            break
                        end
                        results[idx] = remotecall_fetch(p, f, lst[idx])
                    end
                end
            end
        end
    end
    results
end

On first sight, this code is not particularly intuitive. The @spawnlocal macro creates a task on the master process (e.g. process 1). Each task feeds work to a corresponding worker; the call remotecall_fetch(p, f, lst[idx]) function calls f on process p and returns the result when finished. Tasks are uninterruptable and only surrender control at specific points such as remotecall_fetch. Tasks cannot directly modify variables from the enclosing scope, but the same effect can be achieved by using the next_idx function to access and mutate i. The task idiom functions in place of using a loop to poll for results from each worker process.

Implementing our asynchronous algorithm is not much more than a modification of the above code:

# given constants n and 0 < alpha <= 1
# functions initialize and solvesubproblem defined elsewhere
np = nprocs() 
state, subproblems = initialize()
converged = false
isconverged() = converged
function updatemodel(mysubproblem, result)
    # store result
    ...
    # decide whether to generate new subproblems
    state.numback[mysubproblem.parent] += 1
    if state.numback[mysubproblem.parent] >= alpha*n && !state.didtrigger[mysubproblem.parent]
        state.didtrigger[mysubproblem.parent] = true
        # generate newsubproblems by solving linear optimization problem
        ...
        if ... # convergence test
            converged = true
        else
            append!(subproblems, newsubproblems)
            push!(state.didtrigger, false)
            push!(state.numback, 0)
            # ensure that for s in newsubproblems, s.parent == length(state.numback)
        end
    end
end

@sync begin
    for p=1:np
        if p != myid() || np == 1
            @spawnlocal begin
                while !isconverged()
                    if length(subproblems) == 0
                        # no more subproblems but haven't converged yet
                        yield()
                        continue
                    end
                    mysubproblem = shift!(subproblems) # pop subproblem from queue
                    result = remotecall_fetch(p, solvesubproblem, mysubproblem)
                    updatemodel(mysubproblem, result)
                end
            end
        end
    end
end

where state is an instance of a type defined as

type State
    didtrigger::Vector{Bool}
    numback::Vector{Int}
    ...
end

There is little difference in the structure of the code inside the @sync blocks, and the asynchronous logic is encapsulated in the local updatemodel function which conditionally generates new subproblems. A strength of Julia is that functions like pmap are implemented in Julia itself, so that it is particularly straightforward to make modifications like this.

Running it

Now for the fun part. The complete cutting-plane algorithm (along with additional variants) is implemented in JuliaBenders. The code is specialized for stochastic programming where the cutting-plane algorithm is known as the L-shaped method or Benders decomposition and is used to decompose the solution of large linear optimization problems. Here, solvesubproblem entails solving a relatively small linear optimization problem. Test instances are taken from the previously mentioned paper.

We’ll first run on a large multicore server. The runals.jl (asynchronous L-shaped) file contains the algorithm we’ll use. Its usage is

julia runals.jl [data source] [num subproblems] [async param] [block size]

where [num subproblems] is the $n$ as above and [async param] is the proportion $\alpha$. By setting $\alpha = 1$ we obtain the synchronous algorithm. For the asynchronous version we will take $\alpha = 0.6$. The [block size] parameter controls how many subproblems are sent to a worker at once (in the previous code, this value was always 1). We will use 4000 subproblems in our experiments.

To run multiple Julia processes on a shared-memory machine, we pass the -p N option to the julia executable, which will start up N system processes. To execute the asynchronous version with 10 workers, we run

julia -p 12 runals.jl Data/storm 4000 0.6 30

Note that we start 12 processes. These are the 10 workers, the master (which distributes tasks), and another process to perform the master’s computations (an additional refinement which was not described above). Results from various runs are presented in the table below.

	Synchronous		Asynchronous
No. Workers	Speed	Efficiency	Speed	Efficiency
10	154	Baseline	166	Baseline
20	309	100.3%	348	105%
40	517	84%	654	98%
60	674	73%	918	92%

Table: Results on a shared-memory 8x Xeon E7-8850 server. Workers correspond to individual cores. Speed is the rate of subproblems solved per second. Efficiency is calculated as the percent of ideal parallel speedup obtained. The superlinear scaling observed with 20 workers is likely a system artifact.

There are a few more hoops to jump through in order to run on EC2. First we must build a system image (AMI) with Julia installed. Julia connects to workers over ssh, so I found it useful to put my EC2 ssh key on the AMI and also set StrictHostKeyChecking no in /etc/ssh/ssh_config to disable the authenticity prompt when connecting to new workers. Someone will likely correct me on if this is the right approach.

Assuming we have an AMI in place, we can fire up the instances. I used an m3.xlarge instance for the master and m1.medium instances for the workers. (Note: you can save a lot of money by using the spot market.)

To add remote workers on startup, Julia accepts a file with a list of host names through the --machinefile option. We can generate this easily enough by using the EC2 API Tools (Ubuntu package ec2-api-tools) with the command

ec2-describe-instances | grep running | awk '{ print $5; }' > mfile

On the master instance we can then run

julia --machinefile mfile runatr.jl Data/storm 4000 0.6 30

Results from various runs are presented in the table below.

	Synchronous		Asynchronous
No. Workers	Speed	Efficiency	Speed	Efficiency
10	149	Baseline	151	Baseline
20	289	97%	301	99.7%
40	532	89%	602	99.5%

Table: Results on Amazon EC2. Workers correspond to individual m1.medium instances. The master process is run on an m3.xlarge instance.

On both architectures the asynchronous version solves subproblems at a higher rate and has significantly better parallel efficiency. Scaling is better on EC2 than on the shared-memory server likely because the subproblem calculation is memory bound, and so performance is better on the distributed-memory architecture. Anyway, with Julia we can easily experiment on both.

Videos from the Julia tutorial at MIT

2013-03-30T00:00:00-07:00

We held a two day Julia tutorial at MIT in January 2013, which included 10 sessions. MIT Open Courseware and MIT-X graciously provided support for recording of these lectures, so that the wider Julia community can benefit from these sessions.

Julia Lightning Round (slides)

This session is a rapid introduction to julia, using a number of lightning rounds. It uses a number of short examples to demonstrate syntax and features, and gives a quick feel for the language.

Rationale behind Julia and the Vision (slides)

The rationale and vision behind julia, and its design principles are discussed in this session.

Data Analysis with DataFrames (slides)

DataFrames is one of the most widely used Julia packages. This session is an introduction to data analysis with Julia using DataFrames.

Statistical Models in Julia (slides)

This session demonstrates Julia’s statistics capabilities, which are provided by these packages: Distributions, GLM, and LM.

Fast Fourier Transforms

Julia provides a built-in interface to the FFTW library. This session demonstrates the Julia’s signal processing capabilities, such as FFTs and DCTs. Also see the Hadamard package.

Optimization (slides)

This session focuses largely on using Julia for solving linear programming problems. The algebraic modeling language discussed was later released as MathProg. Benchmarks are shown evaluating the performance of Julia for implementing low-level optimization code. Linear programming solvers are available through the CLP, GLPK, and Gurobi packages. The Optim and NLopt packages provide more general optimization algorithms.

Metaprogramming and Macros

Julia is homoiconic: it represents its own code as a data structure of the language itself. Since code is represented by objects that can be created and manipulated from within the language, it is possible for a program to transform and generate its own code. Metaprogramming is described in detail in the Julia manual.

Parallel and Distributed Computing (Lab, Solution)

Parallel and distributed computing have been an integral part of Julia’s capabilities from an early stage. This session describes existing basic capabilities, which can be used as building blocks for higher level parallel libraries.

Networking

Julia provides asynchronous networking I/O using the libuv library. Libuv is a portable networking library created as part of the Node.js project.

Grid of Resistors (Lab, Solution)

The Grid of Resistors is a classic numerical problem to compute the voltages and the effective resistance of a 2n+1 by 2n+2 grid of 1 ohm resistors if a battery is connected to the two center points. As part of this lab, the problem is solved in Julia in a number of different ways such as a vectorized implementation, a devectorized implementation, and using comprehensions, in order to study the performance characteristics of various methods.

Efficient Aggregates in Julia

2013-03-05T00:00:00-08:00

We recently introduced an exciting feature that has been in planning for some time: immutable aggregate types. In fact, we have been planning to do this for so long that this feature is the subject of our issue #13 on GitHub, out of more than 2400 total issues so far.

Essentially, this feature drastically reduces the overhead of user-defined types that represent small number-like values, or that wrap a small number of other objects. Consider an RGB pixel type:

immutable Pixel
    r::Uint8
    g::Uint8
    b::Uint8
end

Instances of this type can now be packed efficiently into arrays, using exactly 3 bytes per object. In all other respects, these objects continue to act like normal first-class objects. To see how we might use this, here is a function that converts an RGB image in standard 24-bit framebuffer format to grayscale:

function rgb2gray!(img::Array{Pixel})
    for i=1:length(img)
        p = img[i]
        v = uint8(0.30*p.r + 0.59*p.g + 0.11*p.b)
        img[i] = Pixel(v,v,v)
    end
end

This code will run blazing fast, performing no memory allocation. We have not done thorough benchmarking, but this is in fact likely to be the fastest way to write this function in Julia from now on.

The key to this behavior is the new immutable keyword, which means instances of the type cannot be modified. At first this sounds like a mere restriction — how come I’m not allowed to modify one? — but what it really means is that the object is identified with its contents, rather than its memory address. A mutable object has “behavior”; it changes over time, and there may be many references to the object, all of which can observe those changes. An immutable object, on the other hand, has only a value, and no time-varying behavior. Its location does not matter. It is “just some bits”.

Julia has always had some immutable values, in the form of bits types, which are used to represent fixed-bit-width numbers. It is highly intuitive that numbers are immutable. If x equals 2, you might later change the value of x, but it is understood that the value of 2 itself does not change. The immutable keyword generalizes this idea to structured data types with named fields. Julia variables and containers, including arrays, are all still mutable. While a Pixel object itself can’t change, a new Pixel can be written over an old one within an array, since the array is mutable.

Let’s take a look at the benefits of this feature.

The compiler and GC have a lot of freedom to move and copy these objects around. This flexibility can be used to store data more efficiently, for example keeping the real and imaginary parts of a complex number in separate registers, or keeping only one part in a register.
Immutable objects are easy to reason about. Some languages, such as C++ and C#, provide “value types”, which have many of the benefits of immutable objects. However, their behavior can be confusing. Consider code like the following:
```
item = lookup(collection, index)
modify!(item)
```
The question here is whether we have modified the same item that is in the collection, or if we have modified a local copy. In Julia there are only two possibilities: either item is mutable, in which case we modified the one and only copy of it, or it is immutable, in which case modifying it is not allowed.
No-overhead data abstractions become possible. It is often useful to define a new type that simply wraps a single value, and modifies its behavior in some way. Our favorite modular integer example type fits this description:
```
immutable ModInt{n} <: Integer
    k::Int
    ModInt(k) = new(mod(k,n))
end
```
Since a given ModInt doesn’t need to exist at a particular address, it can be passed to functions, stored in arrays, and so on, as efficiently as a single Int, with no wrapping overhead. But, in Julia, the overhead will not always be zero. The ModInt type information will “follow the data around” at compile time to the extent possible, but heap-allocated wrappers will be added as needed at run time. Typically these wrappers will be short-lived; if the final destination of a ModInt is in a ModInt array, for example, the wrapper can be discarded when the value is assigned. But if the value is only used locally inside a function, there will most likely be no wrappers at all.
Abstractions are fully enforced. If a custom constructor is written for an immutable type, then all instances will be created by it. Since the constructed objects are never modified, the invariants provided by the constructor cannot be violated. At this time, uninitialized arrays are an exception to this rule. New arrays of “plain data” immutable types have unspecified contents, so it is possible to obtain an invalid value from one. This is usually harmless in practice, since arrays must be initialized anyway, and are often created through functions like zeros that do so.
We can automatically type-specialize fields. Since field values at construction time are final, their types are too, so we learn everything about the type of an immutable object when it is constructed.

There are many potential optimizations here, and we have not implemented all of them yet. But having this feature in place provides another lever to help us improve performance over time.

For now though, we at least have a much simpler implementation of complex numbers, and will be able to take advantage of efficient rational matrices and other similar niceties.

Addendum: Under the hood

For purposes of calling C and writing reflective code, it helps to know a bit about how immutable types are implemented. Before this change, we had types AbstractKind, BitsKind, and CompositeKind, for separating which types are abstract, which are represented by immutable bit strings, and which are mutable aggregates. It was sometimes convenient that the type system reflected these differences, but also a bit unwarranted since all these types participate in the same hierarchy and follow the same subtyping rules.

Now, the type landscape is both simpler and more complex. The three Kinds have been merged into a single kind called DataType. The type of every value in Julia is now either a DataType, or else a tuple type (union types still exist, but of course are always abstract). To find out the details of a DataType’s physical representation, you must query its properties. DataTypes have three boolean properties abstract, mutable, and pointerfree, and an integer property size. The CompositeKind properties names and types are still there to describe fields.

The abstract property indicates that the type was declared with the abstract keyword and has no direct instances. mutable indicates, for concrete types, whether instances are mutable. pointerfree means that instances contain “just data” and no references to other Julia values. size gives the size of an instance in bytes.

What used to be BitsKinds are now DataTypes that are immutable, concrete, have no fields, and have non-zero size. The former CompositeKinds are mutable and concrete, and either have fields or are zero size if they have zero fields. Clearly, new combinations are now possible. We have already mentioned immutable types with fields. We could have the equivalent of mutable BitsKinds, but this combination is not exposed in the language, since it is easily emulated using mutable fields. Another new combination is abstract types with fields, which would allow you to declare that all subtypes of some abstract type should have certain fields. That one is definitely useful, and we plan to provide syntax for it.

Typically, the only time you need to worry about these things is when calling native code, when you want to know whether some array or struct has C-compatible data layout. This is handled by the type predicate isbits(T).

Design and implementation of Julia

2012-08-16T00:00:00-07:00

We describe the design and implementation of Julia in our first paper – Julia: A Fast Dynamic Language for Technical Computing. This is work in progress and comments are appreciated.

New York Open Stats Meetup

2012-04-18T00:00:00-07:00

I’ll be giving a talk on Julia at the New York Open Statistical Programming Meetup on May 1st. After my presentation, John Myles White and Shane Conway are going to give followup demos of statistical applications using Julia. Then we’re going to hang out and grab drinks nearby. Thanks to Harlan Harris and Drew Conway for setting the whole thing up!

Announcement:

After a brief hiatus, we are very excited to announce our May meetup will feature one of the hottest new languages in statistical computing: Julia. We are delighted to welcome Stefan Karpinski, one of the creators of Julia, to give an introduction to the language and his perspective on statistical computing.

Julia is a general-purpose, high-level, dynamic language in the tradition of Lisp, Perl, Python and Ruby. It is designed to take advantage of modern techniques for executing dynamic languages with statically-compiled performance. As part of this design, the language has an expressive type system, which programmers may leverage for dispatch and error checking — incidentally providing the compiler with useful type information. Using types is entirely optional, however: “typeless Julia” is a valid and useful subset of the language, similar to traditional dynamic languages, which nevertheless runs at statically compiled speeds.\

Julia is especially good at running Matlab and R-style programs. Given its level of performance, we envision a new era of technical computing where libraries can be developed in a high-level language instead of C or Fortran. We have also experimented with cloud API integration, and begun to develop a web-based interactive computing environment. The ultimate goal is to make cloud-based supercomputing as easy and accessible as Google Docs.

We will also hear from a mix of people who have already started developing in Julia and see some examples of what they have developed.

The meetup will follow our typical schedule: pizza will begin at 6:15pm, Stefan will begin promptly at 7pm, and we will head to The Central Bar around 8:30pm.

Update: You can see the slides for the talk here. There was no video of the talk, but hopefully the slides are informative — there are, among other things, a lot of code examples that should just work if pasted into the Julia repl.

Lang.NEXT Announcement

2012-03-24T00:00:00-07:00

Jeff and I will be giving a presentation on Julia at the upcoming Lang.NEXT conference, a gathering of “programming language design experts and enthusiasts” featuring “talks, panels and discussion on leading programming language work from industry and research.” We are honored and excited to have been invited to speak at an event alongside so many programming language luminaries.

Abstract:

Julia is a dynamic language in the tradition of Lisp, Perl, Python and Ruby. It aims to advance expressiveness and convenience for scientific and technical computing beyond that of environments like Matlab and NumPy, while simultaneously closing the performance gap with compiled languages like C, C++, Fortran and Java.

Most high-performance dynamic language implementations have taken an existing interpreted language and worked to accelerate its execution. In creating Julia, we have reconsidered the basic language design, taking into account the capabilities of modern JIT compilers and the specific needs of technical computing. Our design includes:

Multiple dispatch as the core language paradigm.
Exposing a sophisticated type system including parametric dependent types.
Dynamic type inference to generate fast code from programs with no declarations.
Aggressive specialization of generated code for types encountered at run-time.

Julia feels light and natural for data exploration and algorithm prototyping, but has performance that lets you deploy your prototypes.

Update: You can see the slides for our talk here. Video of the presentation is available here.

Shelling Out Sucks

2012-03-11T00:00:00-08:00

Spawning a pipeline of connected programs via an intermediate shell — a.k.a. “shelling out” — is a really convenient and effective way to get things done. It’s so handy that some “glue languages,” like Perl and Ruby, even have special syntax for it (backticks). However, shelling out is also a common source of bugs, security holes, unnecessary overhead, and silent failures. Here are the three reasons why shelling out is problematic:

Metacharacter brittleness. When commands are constructed programmatically, the resulting code is almost always brittle: if a variable used to construct the command contains any shell metacharacters, including spaces, the command will likely break and do something very different than what was intended — potentially something quite dangerous.
Indirection and inefficiency. When shelling out, the main program forks and execs a shell process just so that the shell can in turn fork and exec a series of commands with their inputs and outputs appropriately connected. Not only is starting a shell an unnecessary step, but since the main program is not the parent of the pipeline commands, it cannot be notified when they terminate — it can only wait for the pipeline to finish and hope the shell indicates what happened.
Silent failures by default. Errors in shelled out commands don’t automatically become exceptions in most languages. This default leniency leads to code that fails silently when shelled out commands don’t work. Worse still, because of the indirection problem, there are many cases where the failure of a process in a spawned pipeline cannot be detected by the parent process, even if errors are fastidiously checked for.

In the rest of this post, I’ll go over examples demonstrating each of these problems. At the end, I’ll talk about better alternatives to shelling out, and in a followup post, I’ll demonstrate how Julia makes these better alternatives dead simple to use. Examples below are given in Ruby which shells out to Bash, but the same problems exist no matter what language one shells out from: it’s the technique of using an intermediate shell process to spawn external commands that’s at fault, not the language.

Metacharacter Brittleness

Let’s start with a simple example of shelling out from Ruby. Suppose you want to count the number of lines containing the string “foo” in all the files under a directory given as an argument. One option is to write Ruby code that reads the contents of the given directory, finds all the files, opens them and iterates through them looking for the string “foo”. However, that’s a lot of work and it’s going to be much slower than using a pipeline of standard UNIX commands, which are written in C and heavily optimized. The most natural and convenient thing to do in Ruby is to shell out, using backticks to capture output:

`find #{dir} -type f -print0 | xargs -0 grep foo | wc -l`.to_i

This expression interpolates the dir variable into a command, spawns a Bash shell to execute the resulting command, captures the output into a string, and then converts that string to an integer. The command uses the -print0 and -0 options to correctly handle strange characters in file names piped from find to xargs (these options cause file names to be delimited by NULs instead of whitespace). Even with extra-careful options, this code for shelling out is simple and clear. Here it is in action:

irb(main):001:0> dir="src"
=> "src"
irb(main):002:0> `find #{dir} -type f -print0 | xargs -0 grep foo | wc -l`.to_i
=> 5

Great. However, this only works as expected if the directory name dir doesn’t contain any characters that the shell considers special. For example, the shell decides what constitutes a single argument to a command using whitespace. Thus, if the value of dir is a directory name containing a space, this will fail:

irb(main):003:0> dir="source code"
=> "source code"
irb(main):004:0> `find #{dir} -type f -print0 | xargs -0 grep foo | wc -l`.to_i
find: `source': No such file or directory
find: `code': No such file or directory
=> 0

The simple solution to the problem of spaces is to surround the interpolated directory name in quotes, telling the shell to treat spaces inside as normal characters:

irb(main):005:0> `find '#{dir}' -type f -print0 | xargs -0 grep foo | wc -l`.to_i
=> 5

Excellent. So what’s the problem? While this solution addresses the issue of file names with spaces in them, it is still brittle with respect to other shell metacharacters. What if a file name has a quote character in it? Let’s try it. First, let’s create a very weirdly named directory:

bash-3.2$ mkdir "foo'bar"
bash-3.2$ echo foo > "foo'bar"/test.txt
bash-3.2$ ls -ld foo*bar
drwxr-xr-x 3 stefan staff 102 Feb  3 16:17 foo'bar/

That’s an admittedly strange directory name, but it’s perfectly legal in UNIXes of all flavors. Now back to Ruby:

irb(main):006:0> dir="foo'bar"
=> "foo'bar"
irb(main):007:0> `find '#{dir}' -type f -print0  | xargs -0 grep foo | wc -l`.to_i
sh: -c: line 0: unexpected EOF while looking for matching `''
sh: -c: line 1: syntax error: unexpected end of file
=> 0

Doh. Although this may seem like an unlikely corner case that one needn’t realistically worry about, there are serious security ramifications. Suppose the name of the directory came from an untrusted source — like a web submission, or an argument to a setuid program from an untrusted user. Suppose an attacker could arrange for any value of dir they wanted:

irb(main):008:0> dir="foo'; echo MALICIOUS ATTACK 1>&2; echo '"
=> "foo'; echo MALICIOUS ATTACK 1>&2; echo '"
irb(main):009:0> `find '#{dir}' -type f -print0  | xargs -0 grep foo | wc -l`.to_i
find: `foo': No such file or directory
MALICIOUS ATTACK
grep:  -type f -print0
: No such file or directory
=> 0

Your box is now owned. Of course, you could sanitize the value of the dir variable, but there’s a fundamental tug-of-war between security (as limited as possible) and flexibility (as unlimited as possible). The ideal behavior is to allow any directory name, no matter how bizarre, as long as it actually exists, but “defang” all shell metacharacters.

The only two way to fully protect against these sorts of metacharacter attacks — whether malicious or accidental — while still using an external shell to construct the pipeline, is to do full shell metacharacter escaping:

irb(main):010:0> require 'shellwords'
=> true
irb(main):011:0> `find #{Shellwords.shellescape(dir)} -type f -print0  | xargs -0 grep foo | wc -l`.to_i
find: `foo\'; echo MALICIOUS ATTACK 1>&2; echo \'': No such file or directory
=> 0

With shell escaping, this safely attempts to search a very oddly named directory instead of executing the malicious attack. Although shell escaping does work (assuming that there aren’t any mistakes in the shell escaping implementation), realistically, no one actually bothers — it’s too much trouble. Instead, code that shells out with programmatically constructed commands is typically riddled with potential bugs in the best case and massive security holes in the worst case.

Indirection and Inefficiency

If we were using the above code to count the number of lines with the string “foo” in a directory, we would want to check to see if everything worked and respond appropriately if something went wrong. In Ruby, you can check if a shelled out command was successful using the bizarrely named $?.success? indicator:

irb(main):012:0> dir="src"                                                              
=> "src"
irb(main):013:0> `find #{Shellwords.shellescape(dir)} -type f -print0  | xargs -0 grep foo | wc -l`.to_i
=> 5
irb(main):014:0> $?.success?                                                                
=> true

Ok, that correctly indicates success. Let’s make sure that it can detect failure:

irb(main):015:0> dir="nonexistent"                                                              
=> "nonexistent"
irb(main):016:0> `find #{Shellwords.shellescape(dir)} -type f -print0  | xargs -0 grep foo | wc -l`.to_i
find: `nonexistent': No such file or directory
=> 0
irb(main):017:0> $?.success?
=> true

Wait. What?! That wasn’t successful. What’s going on?

The heart of the problem is that when you shell out, the commands in the pipeline are not immediate children of the main program, but rather its grandchildren: the program spawns a shell, which makes a bunch of UNIX pipes, forks child processes, connects inputs and outputs to pipes using the dup2 system call, and then execs the appropriate commands. As a result, your main program is not the parent of the commands in the pipeline, but rather, their grandparent. Therefore, it doesn’t know their process IDs, nor can it wait on them or get their exit statuses when they terminate. The shell process, which is their parent, has to do all of that. Your program can only wait for the shell to finish and see if that was successful. If the shell is only executing a single command, this is fine:

irb(main):018:0> `cat /dev/null`
=> ""
irb(main):019:0> $?.success?
=> true
irb(main):020:0> `cat /dev/nada`
cat: /dev/nada: No such file or directory
=> ""
irb(main):021:0> $?.success?
=> false

Unfortunately, by default the shell is quite lenient about what it considers to be a successful pipeline:

irb(main):022:0> `cat /dev/nada | sort`
cat: /dev/nada: No such file or directory
=> ""
irb(main):023:0> $?.success?
=> true

As long as the last command in a pipeline succeeds — in this case sort — the entire pipeline is considered a success. Thus, even when one or more of the earlier programs in a pipeline fails spectacularly, the last command may not, leading the shell to consider the entire pipeline to be successful. This is probably not what you meant by success.

Bash’s notion of pipeline success can fortunately be made stricter with the pipefail option. This option causes the shell to consider a pipeline successful only if all of its commands are successful:

irb(main):024:0> `set -o pipefail; cat /dev/nada | sort`
cat: /dev/nada: No such file or directory
=> ""
irb(main):025:0> $?.success?
=> false

Since shelling out spawns a new shell every time, this option has to be set for every multi-command pipeline in order to be able to determine its true success status. Of course, just like shell-escaping every interpolated variable, setting pipefail at the start of every command is simply something that no one actually does. Moreover, even with the pipefail option, your program has no way of determining which commands in a pipeline were unsuccessful — it just knows that something somewhere went wrong. While that’s better than silently failing and continuing as if there were no problem, its not very helpful for postmortem debugging: many programs are not as well-behaved as cat and don’t actually identify themselves or the specific problem when printing error messages before going belly up.

Given the other problems caused by the indirection of shelling out, it seems like a barely relevant afterthought to mention that execing a shell process just to spawn a bunch of other processes is inefficient. However, it is a real source of unnecessary overhead: the main process could just do the work the shell does itself. Asking the kernel to fork a process and exec a new program is a non-trivial amount of work. The only reason to have the shell do this work for you is that it’s complicated and hard to get right. The shell makes it easy. So programming languages have traditionally relied on the shell to setup pipelines for them, regardless of the additional overhead and problems caused by indirection.

Silent Failures by Default

Let’s return to our example of shelling out to count “foo” lines. Here’s the total expression we need to use in order to shell out without being susceptible to metacharacter breakage and so we can actually tell whether the entire pipeline succeeded:

`set -o pipefail; find #{Shellwords.shellescape(dir)} -type f -print0  | xargs -0 grep foo | wc -l`.to_i

However, an error isn’t raised by default when a shelled out command fails. To avoid silent errors, we need to explicitly check $?.success? after every time we shell out and raise an exception if it indicates failure. Of course, doing this manually is tedious, and as a result, it largely isn’t done. The default behavior — and therefore the easiest and most common behavior — is to assume that shelled out commands worked and completely ignore failures. To make our “foo” counting example well-behaved, we would have to wrap it in a function like so:

def foo_count(dir)
  n = `set -o pipefail;
       find #{Shellwords.shellescape(dir)} -type f -print0  | xargs -0 grep foo | wc -l`.to_i
  raise("pipeline failed") unless $?.success?
  return n
end

This function behaves the way we would like it to:

irb(main):026:0> foo_count("src")
=> 5
irb(main):027:0> foo_count("source code")
=> 5
irb(main):028:0> foo_count("nonexistent")
find: `nonexistent': No such file or directory
RuntimeError: pipeline failed
    from (irb):5:in `foo_count'
    from (irb):13
    from :0
irb(main):029:0> foo_count("foo'; echo MALICIOUS ATTACK; echo '")
find: `foo\'; echo MALICIOUS ATTACK; echo \'': No such file or directory
RuntimeError: pipeline failed
    from (irb):5:in `foo_count'
    from (irb):14
    from :0

However, this 6-line, 200-character function is a far cry from the clarity and brevity we started with:

`find #{dir} -type f -print0 | xargs -0 grep foo | wc -l`.to_i

If most programmers saw the longer, safer version of this in a program, they’d probably wonder why someone was writing such verbose, cryptic code to get something so simple and straightforward done.

Summary and Remedy

To sum it up, shelling out is great, but making code that shells out bug-free, secure, and not prone to silent failures requires three things that typically aren’t done:

Shell-escaping all values used to construct commands
Prefixing each multi-command pipeline with “set -o pipefail;”
Explicitly checking for failure after each shelled out command.

The trouble is that after doing all of these things, shelling out is no longer terribly convenient, and the code becomes annoyingly verbose. In short, shelling out responsibly kind of sucks.

As is so often the case, the root of all of these problems is relying on a middleman rather than doing things yourself. If a program constructs and executes pipelines itself, it remains in control of all the subprocesses, can determine their individual exit conditions, automatically handle errors appropriately, and give accurate, comprehensive diagnostic messages when things go wrong. Moreover, without a shell to interpret commands, there is also no shell to treat metacharacters specially, and therefore no danger of metacharacter brittleness. Python gets this right: using os.popen to shell out is officially deprecated, and the recommended way to call external programs is to use the subprocess module, which spawns external programs without using a shell. Constructing pipelines using subprocess can be a little verbose, but it is safe and avoids all the problems that shelling out is prone to. In my followup post, I will describe how Julia makes constructing and executing pipelines of external commands as safe as Python’s subprocess and as convenient as shelling out.

Stanford Talk Video

2012-03-01T00:00:00-08:00

Jeff gave his previously announced, invited talk at Stanford yesterday and the video is available here. Congrats, Jeff!

Stanford Talk Announcement

2012-02-27T00:00:00-08:00

I will be speaking about Julia at the Stanford EE Computer Systems Colloquium on Wednesday, February 29 at 4:15PM PST. The title of the talk is Julia: A Fast Dynamic Language For Technical Computing.

Abstract:

Julia is a general-purpose, high-level, dynamic language, designed from the start to take advantage of techniques for executing dynamic languages at statically-compiled language speeds. As a result the language has a more powerful type system, and generally provides better type information to the compiler.

Julia is especially good at running MATLAB and R-style programs. Given its level of performance, we envision a new era of technical computing where libraries can be developed in a high-level language instead of C or FORTRAN. We have also experimented with cloud API integration, and begun to develop a web-based, language-neutral platform for visualization and collaboration. The ultimate goal is to make cloud-based supercomputing as easy and accessible as Google Docs.

Speaker Bio:

Jeff Bezanson has been developing the Julia language for two and a half years with a small distributed team of collaborators. Previously, he worked as a software engineer at Interactive Supercomputing, which developed the Star-P parallel extension to MATLAB. At the company, Jeff was a principal developer of “M#”, an implementation of the MATLAB language running on .NET. He is now a second-year graduate student at MIT. Jeff received an A.B. in Computer Science from Harvard University in 2004, and has experience with applications of technical computing in medical imaging.

The talk will be webcast live.

Edit: the video of the talk can be found here.

Why We Created Julia

2012-02-14T00:00:00-08:00

In short, because we are greedy.

We are power Matlab users. Some of us are Lisp hackers. Some are Pythonistas, others Rubyists, still others Perl hackers. There are those of us who used Mathematica before we could grow facial hair. There are those who still can’t grow facial hair. We’ve generated more R plots than any sane person should. C is our desert island programming language.

We love all of these languages; they are wonderful and powerful. For the work we do — scientific computing, machine learning, data mining, large-scale linear algebra, distributed and parallel computing — each one is perfect for some aspects of the work and terrible for others. Each one is a trade-off.

We are greedy: we want more.

We want a language that’s open source, with a liberal license. We want the speed of C with the dynamism of Ruby. We want a language that’s homoiconic, with true macros like Lisp, but with obvious, familiar mathematical notation like Matlab. We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell. Something that is dirt simple to learn, yet keeps the most serious hackers happy. We want it interactive and we want it compiled.

(Did we mention it should be as fast as C?)

While we’re being demanding, we want something that provides the distributed power of Hadoop — without the kilobytes of boilerplate Java and XML; without being forced to sift through gigabytes of log files on hundreds of machines to find our bugs. We want the power without the layers of impenetrable complexity. We want to write simple scalar loops that compile down to tight machine code using just the registers on a single CPU. We want to write A*B and launch a thousand computations on a thousand machines, calculating a vast matrix product together.

We never want to mention types when we don’t feel like it. But when we need polymorphic functions, we want to use generic programming to write an algorithm just once and apply it to an infinite lattice of types; we want to use multiple dispatch to efficiently pick the best method for all of a function’s arguments, from dozens of method definitions, providing common functionality across drastically different types. Despite all this power, we want the language to be simple and clean.

All this doesn’t seem like too much to ask for, does it?

Even though we recognize that we are inexcusably greedy, we still want to have it all. About two and a half years ago, we set out to create the language of our greed. It’s not complete, but it’s time for a 1.0 release — the language we’ve created is called Julia. It already delivers on 90% of our ungracious demands, and now it needs the ungracious demands of others to shape it further. So, if you are also a greedy, unreasonable, demanding programmer, we want you to give it a try.

The Julia Blog

Passing Julia Callback Functions to C

Sorting with qsort

The problem with closures

Passing closures via pass-through pointers

Warning: qsort_r is not portable

Passing closures in data structures

Taking out the trash (or not)

Put This In Your Pipe

Simple Pipeline, Subtle Bug

Do-Nothing Backticks

Constructing Commands

Further Reading

Distributed Numerical Optimization

Cutting-plane algorithm

Asynchronous variant

Parallel implementation in Julia

Running it

Further reading

Videos from the Julia tutorial at MIT

Julia Lightning Round (slides)

Rationale behind Julia and the Vision (slides)

Data Analysis with DataFrames (slides)

Statistical Models in Julia (slides)

Fast Fourier Transforms

Optimization (slides)

Metaprogramming and Macros

Parallel and Distributed Computing (Lab, Solution)

Networking

Grid of Resistors (Lab, Solution)

Efficient Aggregates in Julia

Addendum: Under the hood

Design and implementation of Julia

New York Open Stats Meetup

Lang.NEXT Announcement

Shelling Out Sucks

Metacharacter Brittleness

Indirection and Inefficiency

Silent Failures by Default

Summary and Remedy

Stanford Talk Video

Stanford Talk Announcement

Why We Created Julia

Sorting with `qsort`

Warning: `qsort_r` is not portable