The Living Thing / Notebooks :

Julia, the programming language

The hippest way to get your IEEE754 on. Hngh.

Usefulness: 🔧 🔧
Novelty: 💡
Uncertainty: 🤪 🤪 🤪
Incompleteness: 🚧 🚧 🚧

Julia: A JIT-compiled language with especial focus on high performance scientific computation.

A Mental Model for Julia: Talking to a Scientist

—Chris Rackauckas, UCI Data Science Initiative Intro to Julia.

Why

tl;dr Not a magic bullet, a handy arrow for the quiver.

Some of Julia’s community makes ambitious claims about Julia being the fastest and bestest thing ever.

Unsurprisingly, Julia is no panacea. It is well designed for numerical computation. In my non-rigorous experiments it seems to do better than other scripting+compilation options, such as cython, on certain tasks. In particular, if you have dynamically defined code in the inner loop it does very well — say, you are doing a Monte Carlo simulation, but don’t know the users’s desired density ahead of time. This is more or less what you’d expect from doing compilation as late as possible rather than shipping a compiled library, but it (at least for my use case) involves less messing around with compilation toolchains, platform libraries, ABIs, makefiles etc.

Julia has its own idiosyncratic frictions. The community process can be problematic (see also giving up on Julia). Library support is patchy, especially at this early stage. It doesn’t run on iphones. In fact it uses lots of memory, so maybe not ideal for embedded controller. Although my colleague Rowan assures me he runs serious Julia code on Raspberry Pi at least all the time, so maybe I need to tighten up my algorithms.

That said, the idea of a science-users-first JIT language is timely, and Julia is that. Python, for example, has clunky legacy issues in the numeric code and a patchy API, and is ill-designed for JIT-compilation, despite various projects that attempt to facilitate it. Matlab is expensive, and nasty for non-numerics, not to mention old and crufty, and code seems to constantly break between MATLAB versions at least as often as it does between Julia versions. Lua has some good science libraries and could likely have filled this niche but for reasons that are perhaps sociological as much as technical doesn’t have the hipness or critical mass of Julia.

The language

Intros

In order of increasing depth

Here’s something that wasn’t obvious to me: What are Symbols?

Exceptions, Errors

What are the Exception types?.

Pitfalls, gotchas

Chris Rackauckas mentions 7 Julia gotchas.

Here are some more.

Implementing methods for custom types is clunky

You want to implement a standard interface on your type so you can, e.g. iterate over it, which commonly looks like this:

for element in iterable
    # body
end

or equivalently

iter_result = iterate(iterable)
while iter_result !== nothing
    (element, state) = iter_result
    # body
    iter_result = iterate(iterable, state)
end

Here is a simple example of that: A range iterator which yields every nth element up to some number of elements could look like

julia> struct EveryNth
           n::Int
           start::Int
           length::Int
       end
julia> function Base.iterate(iter::EveryNth, state=(iter.start, 0))
           element, count = state

           if count >= iter.length
               return nothing
           end

           return (element, (element + iter.n, count + 1))
       end
julia> Base.length(iter::EveryNth) = iter.length
julia> Base.eltype(iter::EveryNth) = Int

(If you are lucky you might be able to inherit from AbstractArray.)

It’s weird for me that this requires injecting your methods into another namespace; in this case, Base. That might feel gross, and it does lead to surprising behaviour that this is how things are done, and some fairly mind-bending namespace resolution rules for methods. Importing one package can magically change the behaviour of another. This monkey patch style (called, in Julia argot, “type piracy”) is everywhere in Julia, and is clearly marked when you write a package, but not when you use the package. Anyway it works fine and I can’t imagine how to handle the multiple dispatch thing better, so deal with it.

The type system is logical, although it’s not obvious if you are used to classical OOP. (Not a criticism.)

It allocates for and copies array slices per default

You are using a chunk of an existing array and don’t want to copy? Consider using views for slices, they say, which means not using slice notation but rather the view function, or the @views macro. Both these are ugly in different ways so I cross my fingers and hope the compiler can optimise away some of this nonsense, but I am not going to check that for now.

Custom containers are scruffy

If you want fast numerics but you are not sure of your float precision, you should not use AbstractFloat for your argument type definitions (although it works fine for simple types). EDIT: You don’t generally need to worry about this in method definitions, as Chris Rackauckas points out.

If you do need well-typed containers, the idiomatic way to do this is using parametric types and parametric methods and so-called orthogonal design.

AFAICT a parametric declaration is the goods if, e.g. you want your method to specialise for arrays of Float32, or of Float64, or of Float16, or of Decimal, but e.g. not an array which mixes types. But also it should work on dense arrays, or view of dense arrays, or sparse arrays thereof, right? Then the correct type would look like one of these:

RealMatrixType = Union{AbstractMatrix{F},SparseMatrixCSC{F}} where {F<:Real}
RealVectorType = Union{AbstractVector{F},SparseVector{F}} where {F<:Real}

EDIT: Don’t bother with this kind of function declaration; the compiler will work it out so long as the arguments’ types are sensible. In fact, it’s counter productive, since you might be unnecessarily restrictive. The rule is: let the compiler work out the argument types in function definitions, but you should choose the types in variable definitions (i.e. when you are calling said functions)

Argument syntax is only OK

Keyword arguments exist but do not participate in method dispatch. Keyword arguments are second-class citizens and might make things slow or stupid if you need to specialise your code based on them. So… design your functions around that. This is usually OK with careful thought, but small lapses leads to many irritating helper functions to handle default arguments.

Containers of containers need careful thinking about parametric types

It’s hard to work out what went wrong when you hear Type definition: expected UnionAll, got TypeVar.

Traits and dispatch

Dispatching is not always obvious.

I always forget one keyword here is Value Types for these, which is what allows one to choose a method based on the value, rather than type of a thing. Their usage was not obvious (to me) from the manual but explained beautifully by Tim Holy. There are dangers.

Chris Rackaukas explain Holy Traits which are a julia idiomatic duck typing method, explained in depth by Mauro3.

🚧

Going fast isn’t magical

Don’t believe that wild-eyed evangelical gentleman at the Julia meetup. (and so far there has been one at every meetup.) He doesn’t write actual code; he’s too busy talking about it. You still need to optimize and hint to get your code being all that you want.

Here is a convenient overview of timer functions to help you work out what to optimize.

Active Julia contributor Chris Rackauckas, in the comments of this blog, is at pains to point out that type-hinting is rarely required in function definitions, and that often even threaded vectorisation annotation is not needed, which is amazing. (The manual is more ambivalent about that latter point, but anything in this realm is a mild bit of magic.)

Not to minimize how excellent that is, but bold coder, still don’t get lazy! This of course doesn’t eliminate role that performance annotations and careful code structuring typically play, and I also do not regard it as a bug that you occasionally annotate code to inform the computer what you know about what can be done with it. Julia should not be spending CPU cycles proving non-trivial correctness results for you that you already know. That would likely be solving redundant, possibly superpolynomially hard problems to speed up polynomial problems, which is a bad bet. It is a design feature of Julia that it allows you to be clever; it would be a bug if it tried to be clever in a stupid way. You will need to think about array views versus copies. And coding so that, See, say, blog posts on how you structure algorithms for fancy GPU computing.

Profiling

The built-in timing tool is called time. It can even time arbitrary begin/end blocks.

Going deeper, try Profile the statistical profiler:

julia> using Profile
julia> @profile myfunc()

ProfileView is advisable for making sense of the output. It doesn’t work in jupyterlab, but from jupyterlab (or anywhere) )ou can invoke the Gtk interface by setting PROFILEVIEW_USEGTK = true in the Main module before using ProfileView.

It’s kind of the opposite of profiling, but you can give use feedback about what’s happening using a ProgressMeter.

Near-magical performance annotations

@simd
equivocally endorsed in the manual, as it may or may not be needed to squeeze best possible performance out of the inner loop. Caution says try it out and see.
@inbounds
disables bounds checking, but be very careful with your indexing.
@fastmath
violates strict IEEE math setup for performance, in a variety of ways some of which sound benign and some not so much.

Shunting to GPUs or other fancy architecture

Coding on GPUs for Julia is, at least in principle, easy. Thus you can possibly glean some nearly-free acceleration by using, for example, CuArrays to push computation onto GPUs. sdanisch has written the simplest examples of this. AFAICT the hippest interface is the vendor-neutral GPUArrays which, e.g. (I think) for CUDA will transparently use CUArrays and CUDANative, for e.g. an NVIDIA GPU. You can also invoke them in a vendor-specific way, which looks easy. There is even activist support for google cloud TPU compilation. How well this vectorises in practice is not clear to me yet. GPUify loops provides “Support for writing loop-based code that executes both on CPU and GPU”. Simon Danisch writes on GPUs usage in julia.

Installation

Bareback

Distribution packages are generally veeery old in Julia years; if you want to faddish features that are the whole reason you got lured to Julia, download from the main site.

On macOS you could do this to put the binary in your path.

ln -s /Applications/Julia-1.1.app/Contents/Resources/Julia/bin/Julia ~/bin/Julia

JuliaPro

Infrequent but luxurious Julia installs provided through JuliaPro.

I don’t use these since I had early problems with them which turned out to be in the Intel MKL builds. Probably resolved now since Intel MKL is no longer packaged.1 :shrug:

They come with everything packaged up and are a simple way of getting started. I’m not 100% clear on how they do package managmeent and interact with other Julia installations since I don’t use them.

Workflow

General tips: Look out for handy interactive introspection tricks, e.g. @which tells you which method your function invocation would invoke. There are other macros to tell you while file it comes from etc.

Use Revise.jl, and Project environments.

(v1.0) pkg> generate PatternMachine
Generating project PatternMachine:
    PatternMachine/Project.toml
    PatternMachine/src/PatternMachine.jl

shell> cd PatternMachine
/media/dan/SchnittPunkt/Source/PatternMachine

(v1.0) pkg> activate .

julia> import PatternMachine
[ Info: Precompiling PatternMachine [93e276e4-e6da-11e8-1ff8-9f9dfed081b5]

(PatternMachine) pkg> add DSP

If you don’t want to manually invoke revise to detect code updates, you can use Revise automatically: To .Julia/config/startup.jl add

atreplinit() do repl
    try
        @eval using Revise
        @async Revise.wait_steal_repl_backend()
    catch
    end
end

and (they claim) that to Julia/config/startup_iJulia.jl you must add

try
    @eval using Revise
catch
end

But for me this leads to the kernel hanging shortly after startup, and the non-IJulia setup is sufficient.

For VS code users .Julia/config/startup.jl should purportedly be

atreplinit() do REPL
    @schedule begin
        sleep(0.1)
        try
            @eval using Revise
        catch err
            warn("Could not load Revise.")
        end
    end
end

Here is Erik Engheim’s workflow walk-through.

IDEs/workbooks

There is a reasonable IDE called Juno, built on Atom. There is jupyter integration through IJulia. I personally just use VS Code to edit code but execute it via the REPL or IJulia.

All these methods have their own joys and frustrations.

Juno

Juno I think is the default. It has magical features - integrated debugger and, I am taold, automtic substitution of LaTeX with actual unicode. Nice.

Juno is single-window only so you can’t use multiple monitors, and thus you end up squinting at tiny windows of code hidden between all the outputs. Atom’s panes aren’t well-designed for this use-case. For me that means there are a million tiny frictions distracting me from writing code in Juno. I can’t stand it.

If you install Juno as an app, but you also already use Jupyter, there is an additional annoyance in that it hijacks your Atom install in a confusing way and mangles your various package preferences. If you love Juno and Atom, I recommend installing Juno from within atom via the uber-juno package.

Possibly you can bypass this using homebrew? I didn’t try. But maybe give this a burl:

brew cask install juno

I personally don’t find juno worth the irritations, and never use it. Instead…

IJulia

IJulia isn’t an IDE er se, it’s a sort-of-light interactive notebook. I don’t want to edit code in Julia; for that I use a text editor. There are loads of those. (I use visual studio code.) However, for presenting experiments and prototyping this is good. These two components have have a less annoying, zippier and more stable workflow for my style by not trying to do everything badly, which is what Juno seems to me to be all about.

A serious plus for the jupyter option is that if you wish to share your results with colleagues you can use juliabox, an simple online julia host.

IJulia is also not flawless. For a start, it’s built on Jupyter, which I don’t love.

For another thing, it does overzealous installs per default, profligately installing another copy of jupyter, which you then have to maintain separately to the one you probably already had. Boring. You can bypass this by commanding it to use the perfectly good other jupyter:

ENV["JUPYTER"] = "/usr/local/bin/jupyter"
Pkg.add("IJulia")

Now Julia appears as a normal kernel in your jupyter setup. (This is only exciting for habitual jupyter users and other anal retentives.)

If you want particular jupyter kernels for particular Julia environments, it is possible: by using the following kernel interpreter template:

Julia -e 'import Pkg; Pkg.activate("myproject")' -L /path/to/kernel.jl

Literate coding

There is a package called Weave.jl which is inspired by R’s knitr but compatible with jupyter, which is also not an IDE, but a good way of invoking reproducible self-documenting code. It could probably be used to fashion a working academic paper if you used the right pandoc hacks.

Pkg.add("Weave")

NB you can also use RMarkdown in Julia mode It’s not clear to me how graphing works in this setup.

See also Literate.jl for some similar functionality.

API, FFIs

See the API list

C

Sort of easy, but there is a tedious need to define the call signature at call time. Survivable.

R

XRJulia:

This package provides an interface from R to Julia, based on the XR structure, as implemented in the XR package, in this repository.

RJulia:

rJulia provides an interface between R and Julia. It allows a user to run a script in Julia from R, and maps objects between the two languages.

Python

tl;dr:

conda create -n conda_jl python
conda activate conda_jl
env CONDA_JL_HOME="$HOME/anaconda3/envs/conda_jl" \
    CONDA_JL_VERSION=3 \
    PYTHON=(which python) \
    JUPYTER=(which jupyter) \
    Julia
using Pkg
Pkg.add("IJulia")
Pkg.add("Conda")
Pkg.resolve()
Pkg.build()

PyCall.jl invokes python. Per default it installs yet another conda python, via Conda.jl, and defaults to the elderly python 2.7. This is irritating for various reasons, such as being ancient, and eating your diskspace with yet more versions of the same shit that you already have installed if you use python, but even more decrepit versions.

Here is how to use an existing version

env PYTHON=(which python3) \
    julia
Pkg.build("PyCall")

Here is how you use Conda, but with python 3:

ENV["CONDA_JL_VERSION"] = "3"
Pkg.build("Conda")

Here is how you use an existing environment

conda create -n conda_jl python
export CONDA_JL_HOME=~/anaconda3/envs/conda_jl
julia -e 'Pkg.build("Conda")'

Either way you should regularly run conda clean to stop your disk filling up with obsolete versions of obscure dependencies for that package you tried out that one time as per standard practice.

conda clean -pt

Pretty printing, formatting

Strings, care and feeding

There are many formatting libraries, because everyone seems to dislike the built-in options, which have various limitations.

Many things seems to be based on Formatting.jl which apes python string interpolation and this is IMO a baseline good idea for formatting numbers. There is a friendly-ish fork, Format.jl which has a slightly different philsophy

This package offers Python-style general formatting and c-style numerical formatting (for speed).

Rich display

Julia has an elaborate display system for types, as illustrated by Simon Dernisch and Henry Shurkus.

tl;dr To display an object, invoke

display(instance)

Say you defined MyType. To ensure MyType displays sanely, define

show(io, instance::MyType)

e.g.

function show(io::IO, inst::MyType)
    print(io, "MyType(length=$(inst))")
end

Matti Pastell’s useful reproducible document rendering system Weave.jl supports magical table displaying for Latex/HTML. TexTables.jl is a specialised table renderer for scientific table output. More specialised yet, RegressionTables does statistical tables with baked-in analysis and a side-order of formatting. (Bonus: its launch announcement is a tour of the statistical ecosystem.)

Mustache.jl also has some traction, being a templating syntax used by some other things.

Latexify marks up certain Julia objects nicely for mathematical TeX display.

latexify("x/y") |> print
$\frac{x}{y}$

Debugging and coding

There is a tracing debug ecosystem based on JuliaInterpreter. This provides such useful commands as @bp to inject breakpoints and @enter to execute a function in a step debugger. Doesn’t work from IJulia, but works fine from the console and also from Juno. (Juno is still, for my tastes, unpleasant enough that this is still not worth it..)

Additionally, there is a somewhat maintained stack debugger, Gallium.jl.

If you want to know how you made some weird mistake in type design, try Cthulhu

Linter, Lint.jl can identify code style issues. (also has an atom linter plugin).

Iteration/Indexing

Julia arrays are 1-indexed, right? No. They have arbitrary indices. Some care is needed to support the (recently introduced) arbitrary indexing. So best practice is not to count from 1:length(A).

If you just need to visit each element, iterate using eachindex.

If you want to iterate by axes… instead of

function mycopy!(dest::AbstractVector, src::AbstractVector)
    length(dest) == length(src) || throw(
        DimensionMismatch("vectors must match"))
    # OK, now we’re safe to use @inbounds, right? (not anymore!)
    for i = 1:length(src)
        @inbounds dest[i] = src[i]
    end
    dest
end

you do

function mycopy!(dest::AbstractVector, src::AbstractVector)
    axes(dest) == axes(src) || throw(
        DimensionMismatch("vectors must match"))
    for i in LinearIndices(src)
        @inbounds dest[i] = src[i]
    end
    dest
end

Generic multidimensional algorithms are not necessarily intuitive but are possible. There are several methods.

Broadcasting is often simplest and fastest but requires some thought.

Dot-fusion/broadcasting

There are some couple issues here; dot syntax in julia provides a means to both compactly express vectorized operations in julia and also to do smart broadcasting and operation fusing (which compiles functions down to vectorised kernels collapsing intermediate operations.) The details are non-obvious, but you can usually ignore them and just use the @. macro to magically invoke it. In most cases it works transparently on your own functions too.

Pro tip: In practice this needs reshape to get the right dimensions lined up.

Note the manual says

In Julia, vectorized functions are not required for performance, and indeed it is often beneficial to write your own loops (see Performance Tips), but they can still be convenient.

This is directly contradicted shortly after by the same manual and also the ephemeral julia discussion on the net. As some recent packages note you probably do want to use this syntax for vector operations because it works on GPUs, whereas using basic iteration syntax does not go fast on GPUs.

This is not generally a huge problem; vectorised syntax is not so opaque.

But what if I want to do custom broadcasting because, e.g. I can do some operation faster using an FFT, and there needs to be an FFT inside the vectorising process?

How do I define custom broadcasting for this kind of case for some type in julia?

I’ll try to work this out from the moving parts in the ApproxFun system.

import Base:  broadcast, broadcast!, similar

import Base.Broadcast: BroadcastStyle, Broadcasted,
        AbstractArrayStyle, broadcastable,
        DefaultArrayStyle, broadcasted

for op in (:+,:-,:*,:/,:^)
    @eval begin
        broadcast(::typeof($op), a::Fun, b::Fun) = $op(a,b)
        broadcast(::typeof($op), a::Fun, b::Number) = $op(a,b)
        broadcast(::typeof($op), a::Number, b::Fun) = $op(a,b)
    end
end

## broadcasting
# for broadcasting, we support broadcasting over `Fun`s, e.g.
#
#    exp.(f) is equivalent to Fun(x->exp(f(x)),domain(f)),
#    exp.(f .+ g) is equivalent to Fun(x->exp(f(x)+g(x)),domain(f) ∪ domain(g)),
#    exp.(f .+ 2) is equivalent to Fun(x->exp(f(x)+2),domain(f)),
#
# When we are broadcasting over arrays and scalar Fun’s together,
# it broadcasts over the Array and treats the scalar Fun’s as constants,
# so will not necessarily call the constructor:
#
#    exp.( x .+ [1,2,3]) is equivalent to [exp(x + 1),exp(x+2),exp(x+3)]

# When array values are mixed with arrays, the Array
# takes precedence:
#
#    exp.([x;x] .+ [x,x]) is equivalent to exp.(Array([x;x]) .+ [x,x])
#
# This precedence is picked by the `promote_containertype` overrides.

struct FunStyle <: BroadcastStyle end

BroadcastStyle(::Type{<:Fun}) = FunStyle()

BroadcastStyle(::FunStyle, ::FunStyle) = FunStyle()
BroadcastStyle(::AbstractArrayStyle{0}, ::FunStyle) = FunStyle()
BroadcastStyle(::FunStyle, ::AbstractArrayStyle{0}) = FunStyle()
BroadcastStyle(A::AbstractArrayStyle, ::FunStyle) = A
BroadcastStyle(::FunStyle, A::AbstractArrayStyle) = A

# Treat Array Fun’s like Arrays when broadcasting with an Array
# note this only gets called when containertype returns Array,
# so will not be used when no argument is an Array
Base.broadcast_axes(::Type{Fun}, A) = axes(A)
Base.broadcastable(x::Fun) = x

broadcastdomain(b) = AnyDomain()
broadcastdomain(b::Fun) = domain(b)
broadcastdomain(b::Broadcasted) = mapreduce(broadcastdomain, ∪, b.args)

broadcasteval(f::Function, x) = f(x)
broadcasteval(c, x) = c
broadcasteval(c::Ref, x) = c.x
broadcasteval(b::Broadcasted, x) = b.f(broadcasteval.(b.args, x)...)

But how do I get custom broadcasting?

🚧

Indexing systems

For some reason there seems to be shame attached to using macros for this, but there is a fairly simple and transparent macro system called [email protected]; See Mike Innes’ explanation.

If you want to avoid macros but you can’t work out how to get your broadcasts to line up, there is the CartesianIndices system, which is compact, general but not especially clear. Here is the “surprisingly simple” (sic) one-pole filter implementation

function expfiltdim(x, dim::Integer, α)
    s = similar(x)
    Rpre = CartesianIndices(size(x)[1:dim-1])
    Rpost = CartesianIndices(size(x)[dim+1:end])
    _expfilt!(s, x, α, Rpre, size(x, dim), Rpost)
end

function _expfilt!(s, x, α, Rpre, n, Rpost)
    for Ipost in Rpost
        # Initialize the first value along the filtered dimension
        for Ipre in Rpre
            s[Ipre, 1, Ipost] = x[Ipre, 1, Ipost]
        end
        # Handle all other entries
        for i = 2:n
            for Ipre in Rpre
                s[Ipre, i, Ipost] = α*x[Ipre, i, Ipost] +
                    (1-α)*s[Ipre, i-1, Ipost]
            end
        end
    end
    s
end

🚧 iscover whether this is as efficient as broadcasting as far as operation fusing etc.

You can avoid any of these for a couple of common cases. EllipsisNotation provides compact intuitive syntactic sugar if you wish to get at the first or last axis.

A[..,1] = [2 1 4 5
           2 2 3 6]

A[..,2] = [3 2 6 5
          3 2 6 6]

A[:,:,1] == [2 1 4 5
             2 2 3 6] #true

For iterating over one index of an arbitrary array, mapslices is convenient, or for extracting a slice along an arbitrary index, selectdim.

There is also AxisAlgorithms, which does batch matrix multiplications and inversions for a couple of handy cases.

Arrays are column-major

The first index should change fastest. That is, walk over the last index in your outer loop, and apply operations to slices comprising the other indices.

Multidimensional arrays in Julia are stored in column-major order. This means that arrays are stacked one column at a time. This can be verified using the vec function or the syntax [:]copy_col_row is much faster than copy_row_col because it follows our rule of thumb that the first element to appear in a slice expression should be coupled with the inner-most loop. This is expected because copy_cols respects the column-based memory layout of the Matrix and fills it one column at a time. Additionally,… it follows our rule of thumb that the first element to appear in a slice expression should be coupled with the inner-most loop.

Einstein summation

Einstein summation is a handy way of thinking about lots of tensor operations which is not as esoteric as it sounds, inclduing, e.g. batch matrix multiply. TensorOperations.jl does some optimised einstein summations. Einsum.jl does reputedly more einstein summations but with possibly less optimistaion.

Structured matrices

Chris says

One of the best scientific computing packages is BandedMatrices.jl. It lets you directly define banded matrices and then overloads things like * and \ to use the right parts of BLAS. Compared to just using a sparse matrix (the standard MATLAB/Python way), this is SCREAMING FAST (the QR factorization difference is pretty big). Banded matrices show up everywhere in scientific computing (pretty much every PDE has a banded matrix operator). If it’s not a banded matrix, then you probably have a block-banded matrix, and thank god that (???) also wrote BlockBandedMatrices.jl which will utilize the same speed tricks but expanded to block banded matrices (so, discretizations of systems of PDEs).

Data loading/saving/exchange

The names are nearly all self explaining

Signal processing

DSP.jl has been split off from core and now needs to be installed separately. Also DirectConvolutions has sensible convolution code.

FFTs are provided by AbstractFFTs, which in-principle wraps many FFT implementations. I don’t know if there is a GPU implementation yet, but there for sure is the classic CPU implementation provided by FFTW.jl which uses FFTW internally.

As for how to use these things, Numerical tours of data scienceshas a Julia edition with lots of signal processing conent.

JuliaAudio processes audio. They recommend PortAudio.jl as a real time soundcard interface, which looks sorta simple. See rkat’s example of how this works. There are useful abstractions like SampledSignals to load audio and keep the data and signal rate bundled together. Although, as SampledSignal maintainer Spencer Russell points out, AxisArrays might be the right data structure for sample signals, and you could use SampledSignals purely for IO, and ignore its data structures thereafter.

Images.jl processes images.

QMC

Low discrepancy and other QMC stuff. Mostly I want low discrepancy sequences. There are two options with near identical interfaces; I’m not sure of the differences.

Sobol.jl claims to have been performance profiled:

] add Sobol
using Sobol
s = SobolSeq(2)
# Then
x = next!(s)

QMC.jl:

] add https://github.com/PieterjanRobbe/QMC.jl
using QMC
lat = LatSeq(2)
#then
next(lat)

Matrix Factorisation

NMF.jl contains reference implementations of non-negative matrix factorisation.

Laplacians.jl by Dan Spielman et al is a matrix factorisation toolkit especially for Laplacian (graph adjacency) matrices.

Plotting

Due to the explosion of options, segmented off into a separate Julia plotting page.

Destructuring assignments

a.k.a. splats. Not built in, but very fancy ones are provided by the macro Destructure.jl.

Differentiating, optimisation

Optimizing

JuMP support many types of optimisation, including over non-continuous domains, and is part of the JuliaOpt family of confusingly diverse optimizers, which invoke various sub-families of optimizers. The famous NLOpt solvers comprise one such class, and they can additionally be invoked separately.

Unlike NLOpt and the JuMP family, Optim.jl (part of JuliaNLSolvers, a different family entirely) solves optimisation problems purely inside Julia. It has nieces and nephews such as LsqFit for Levenberg-Marquardt non-linear least squares fits. Optim will automatically invoke ForwardDiff. Assumes mostly unconstrained problems.

Krylov.jl is a collection of Krylov-type iterative method for large iterative linear and least-squares objectives.

Autodiff

Julia is a hotbed of autodiff for technical and community reasons. Such a hotbed that it’s worth discussing in the autodiff notebook.

Closely related, projects like ModelingToolkit.jl blur the lines between equations and coding, and allow easy definition of differentiable or probabilistic programming.

Statistics, probability and data analysis

Hayden Klok and Yoni Nazarathy are writing a free Julia Statistics textbook which seems a thorough introduction to statistics as well as Julia, albeit statistics in a classical frame that won’t be fashionable with either your learning theory or Bayesian types.

A good starting point for doing stuff is JuliaStats which organisation produces many statistics megapackages, for kernel density estimates, generalised linear models etc. Install them all using Statskit:

using StatsKit

Data frames

The workhorse data structure of statistics.

Data frames are provided by DataFrames.jl.

There are some older ones you might encoutner such as DataTables.jl which are subtly incompatible in tedious ways which we can hopefully ignore soon. For now, use IterableTables.jl to translate where needed between these and many other data sources.

One can access DataFrames (And DataTables and SQL databases and streaming data sources) using Query.jl. DataFramesMeta has also been recommended. You can load a lot of the R standard datasets using RDatasets.

using RDatasets
iris = dataset("datasets", "iris")
neuro = dataset("boot", "neuro")

Frequentist statistics

Lasso and other sparse regressions are available in Lasso.jl which reimplements the lasso algorithm in pure Julia, GLMNET.jl which wrap the classic Friedman FORTAN code for same. There is also (functionality unattested) an orthogonal matching pursuit one called OMP.jl but that algorithm is simple enough to bang out oneself in an afternoon, so no stress if it doesn’t work. Incremental/online versions of (presumably exponential family) statistics are in OnlineStats. MixedModels.jl

is a Julia package providing capabilities for fitting and examining linear and generalized linear mixed-effect models. It is similar in scope to the lme4 package for R.

Sampling-based inference

Another notable product of JuliaStats organisation is Distributions.jl, a probability distribution toolkit providing densities, sampling etc. It is complemented by the nifty Bridge.jl which simulates random processes with univariate indices.

Turing.jl does simulation-based posterior inference, in a compositional model setting.

Turing.jl is a Julia library for (universal) probabilistic programming. Current features include:

More recently and with a bigger splash Gen

Gen simplifies the use of probabilistic modeling and inference, by providing modeling languages in which users express models, and high-level programming constructs that automate aspects of inference.

Like some probabilistic programming research languages, Gen includes universal modeling languages that can represent any model, including models with stochastic structure, discrete and continuous random variables, and simulators. However, Gen is distinguished by the flexibility that it affords to users for customizing their inference algorithm.

Gen’s flexible modeling and inference programming capabilities unify symbolic, neural, probabilistic, and simulation-based approaches to modeling and inference, including causal modeling, symbolic programming, deep learning, hierarchical Bayesiam modeling, graphics and physics engines, and planning and reinforcement learning.

Hmm.

DynamicHMC does Hamiltonian/NUTS sampling in a raw likelihood setting.

Possibly it is a competitor of Klara.jl, the Juliastats MCMC.

If finance is your game, the Julia team has manufactured a fintech toolbox Juliafin. This possibly pays the bills for them.

Machine learning

Let’s put the automatic differntiation, the optimizers and the samplers together to do differentiable learning!

The deep learning toolkits have shorter feature lists than the lengthy ones of those fancy python/C++ libraries (e.g. mobile app building, cuDNN-backed optimisations are all less present in julia libraries) But maybe elegance/performance of Julia makes some of those features irrelevant? I for one don’t care about most of those because I’m a researcher not a deployer.

Having said that, Tensorflow.jl gets all the features, becaues it invokes C++ tensorflow. Surely one misses the benefit of Juliathgis way, since there are two different array-processing infrastructures to pass between, and a very different approach to JIT versus pre-compiled execution. Or no?

Flux.jl sounds like a reimplementation of Tensorflow-style differentiable programming inside Julia, which strikes me as the right way to do this to benefit from the end-to-end-optimised design philosophy of Julia.

Flux is a library for machine learning. It comes “batteries-included” with many useful tools built in, but also lets you use the full power of the Julia language where you need it. The whole stack is implemented in clean Julia code (right down to the GPU kernels) and any part can be tweaked to your liking.

It’s missing some features of Tensorflow, but has simple dynamism like Pytorch, and includes compensatory suprising/unique feature combinationss. Its GPU support is real, and elegant, but AFAICS more proof-of-concept than production-ready, relying upon the supposition that CuArrays can represent all the operations I need and will perform them optimally, and that I don’t need any fancy DNN-specific GPU optimizations. This is dubious – for example, it does not have all the FFT operations I want, such as the Discrete Cosine Transform. Its end-to-end Julia philosophy is justifed, IMO, by twists like DiffEqFlux – see below – which makes Neural ODEs sort-of simple to create.

Alternatively, Mocha.jl is a belt-and-braces deep learning thing, with a library of pre-defined layers. Swap out some of the buzzwords of Flux with newer ones, and skip some older ones, and there you are. It doesn’t appear, on brief perusal, to be as flexible as Flux, missing state filters and recurrent nets and many other models that are made less pure by their distance from the platonic ideal of deep learning, the convnet. As such it is not much use to me, despite promising clever things about using CUDA etc. YMMV.

Knet.jl is another deep learning library that claims to show the ease of implementing deep learning frameworks in Julia. Their RNN support looks a little limp, and they are using the less-popular autodiff frameworks, so I won’t be using ’em but point is well taken.

If one were aiming to do that, why not do something left-field like use the dynamical systems approach to deep learning? This neat trick was popularised by Haber and Ruthotto et al, who have released some of their models as Meganet.jl. I’m curious to see how they work.

A straight, if not fancy package for Gaussian Process is Stheno. It aspires to play well with Turing.jl for non-Gaussian likelihood and Flux.jl for deep Gaussian processes. Currently has automatic differentiation problems.

MLJ is a scikit-learn-like pipeline for data analysis in Julia which standardises model composition automates some of the training etc. It has various adaptors for other ML systems via MLJModels.

ODEs, PDEs, SDEs

Chris Rauckackas is a veritable wizard with this stuff; just read his blog.

Here is a tour of fun tricks with stochastic PDEs. There is a lot of tooling for this; DiffEqOperators … does something. DiffEqFlux (EZ neural ODEs works with Flux and claims to make Neural ODEs simple. The implementation of these things in python, for the award-winning NeurIPS paper that made them famous was a nightmare. +1 for Julia here. The neural SDE section is mostly julia; Go check that out.

Graphs/networks

JuliaGraphs supports a number of graph algorithms. It seems to be mostly simple, matrix-ish graph representations.

Laplacians.jl by Dan Spielman et al is an advanced matrix factorisation toolkit for Laplacian (graph adjacency) matrices, i.e. networks viewed as matrices.

A little less vibrant but maybe still good is JuNet which supports a similar set of algorithms, but seems, AFAICT, to concentrate on node-edge representations over matrix representations.

Approximating and interpolating

ApproxFun.jl does Chebychev and Fourier approximations of given functions This is not, at least primarily, a tool for data analysis, but for solving eigenfunction problems and such like using computational Hilbert space methods for functions which are otherwise difficult. In particular, we must be able to evaluate the target functions at arbitrary points to construct the interpolant, rather than at, say, provided sample points. Useful companion package FastTransforms converts between different basis function representaitons. Interpolations.jl does arbitrary order spline interpolations of mathematical functions, but also data. This enables some very clever tricks, e.g. approximate random sampling of tricky distributions.

Caching

Doing heavy calculations? Want to save time? Caching notionally helps, but of course is a quagmire of bikeshedding.

The Anamemsis creator(s) have opinion about this:

Julia is a scientific programming language which was designed from the ground up with computational efficiency as one of its goals. As such, even ad hoc and “naive” functions written in Julia tend to be extremely efficient. Furthermore, any memoization implementation is liable to have some performance pain points: in general they contain type unstable code (even if they have type stable output) and they must include a value cache which is accessible from the same scope as the function itself so that global objects are potentially involved. Because of this, memoization comes with a significant overhead, even with Julia’s highly efficient AbstractDict implementations.

My own opinion of caching is that I’m for it. Despite the popularity of the “memoize” pattern I’m against that as a rule. I usually end up having to think carefully about caching, and magically making functions cache never seems to help with the hard part of that. Caching logic should, for projects I’ve worked on at least, be explicit, opt-in, and that should be an acceptable design overhead because I should only be doing it for a few expensive functions, rather than creating magic cache functions everywhere.

LRUCache is the most commonly recommended cache. It provides simple, explicit, opt-in, in-memory cache structure. I am unclear how this would perform with multiple processes or threads.

Memoize provides an automatic caching macro but primitive in-memory caching using plain old dictionaries. AFAICT it is possibly compatible with LRUCache, but weirdly I could not find discussion of that. However, I don’t like memoisation, so that’s the end of that story.

Anamemsis has advanced cache decorators and lots of levers to pull to alter caching behaviour.

Caching is especially full-featured with macros, disk persistence and configurable evictions etc. It wants type annotations to be used in the function definitions.

Cleverman Simon Byrne notes that you can roll your own fancy caching using wild Julia introspection stunts via Cassette.jl, but I don’t know how battle-tested this is. Possibly this is the correct amount of magic?

Units

People sometimes mention that Julia can handle physical units in a syntactically pleasant fashion, but do not actually explain it. It’s not clear how, since all the documentation assumes I already know. From what I can tell, I first need to actually install Unitful, which includes some useful units, and in particular Quantity types.

Then I can use units in the following fashion:

$ using Unitful
$ const s = u"s"
$ const minute = u"minute"
$ 5.3s + 2.0minute
125.3 s

To know know how this works, and also how I can invent your own units, I read Erik Engheim who explains it in depth.

See SampledSignals for a concrete example of how to do things with units, such as the following method definitions to handle time-to-index conversion.

toindex(buf::SampleBuf, t::Number) = t
toindex(buf::SampleBuf, t::Unitful.Time) = inframes(Int, t, samplerate(buf)) + 1
toindex(buf::SampleBuf, t::Quantity) = throw(Unitful.DimensionError(t, s))

UIs and servers

Escher/Interact ecosystem

Interact.jl is a UI/widget framework that is reasonably front-end agnostic. A web-app-style interface that can deploy to IDE, web server, or desktop electron via Bink.jl. How much tedious javascript management is required I do not know.

(If you really wanted to have lots of control and a proper web framework, you could use the web framework Genie.jl.)

Immediate mode

You can plug into classic C++ immediate mode Dear Imgui via CImGui.jl.

Reactive programming

Reactive.jl is widely used to manage signals in UI toolkits. It has a competitor Signals.jl, whose creator describes the differences:

Signals.jl, while offering the same functionality as Reactive, is different on some key factors

The last point is, IMO, a minus, but the second one is a plus.

Traditional toolkits

QML.jl magically plugs into Qt5. Gtk.jl plugs into GTK.


  1. In general if you are a research you should avoid Intel MKL because you are unlikely to truly want to make your code many gigabytes larger, complicate the licensing, and crash more often to maybe shave a few points off the execution time for some hypothetical end user who you will never contact.↩︎