Julia and Julia sets

You can read this article on my own blog, if you prefer that to medium.

If I were forced to make a bet, on which programming language will revolutionize scientific software in the next 5 years, I would put my money on Julia.

The system programming revolution that post-metaclasses C++ and Rust will bring to systems engineering, Julia will bring to scientific computing.

What is Julia and why should you care about it ?

Despite looking very similar to Python, from the end user’s perspective, Julia has a few tricks up its sleeves.

Once Julia code starts executing, it’s a match for any well-written C, C++ or Fortran implementation of that same program.

Julia makes parallelism easy and allows you to write pretty and concise code that runs fast. You can read more about what makes Julia so fast and easy to integrate & deploy on its official website.

Julia is still a dynamically typed language, which means values are typed, but variables aren’t. But, it provides a very easy syntax for enforcing the type of a value at its creation. Unlike Python, its type system is very advanced, meaning you can do something like this:

arr = Array{Float16,1}()push!(arr, 12.9) # Workspush!(arr, "my string") # Give a type errorval = UInt8(11)push!(arr, val) # `val` is cast to Float16 then appended

It also provides native support for things like multidimensional arrays and array operations using SIMD optimizations. Further more, its arrays are actual arrays, not lists, but large chunks of memory, allocated continuously, for fast iteration and lookup.

All this, means that Julia can achieve amazing speeds for vector computations, without needing 3rd party containers or linear algebra libraries.

Generating a Julia set

Since it’s such a low hanging fruit, I might as well show you Julia by generating and plotting some Julia sets.

For the rest of this blog, I will try and present to you some Julia code, to showcase some of its features. Hopefully I’ll convince you of its efficiency.

To follow through you will need to install Julia and open up a .jl file in a text editor of your choice. You can run the code using julia my_file.jl.

One thing I like about Julia, is that it takes most of the clutter out of writing fast procedural code. Take, for example, these functions that generate our Julia set matrix.

function julia(x, y, width, height, c)    # Scale size the fractal within our frame    z = ((y/width)*2.7 - 1.3) + ((x/height)*4.5 - 2.5)im    for i = 1:255        z = z^2 + c        if abs(z) >= 4            return i        end    end    return -255endfunction julia_set(height, width, c)    [julia(x, y, width, height, c) for x = 1:height, y = 1:width]end

We can invoke julia_set with a given image size and c parameter, then plot a heatmap of the resulting matrix to view the fractal.

Plotting a heatmap

Another nice thing about Julia, are its community maintained libraries. The ease, and practically non existent overhead, with which C and Python libraries can be mixed within Julia code, gives us the ability to use many mature libraries where a native Julia alternative doesn’t exist or isn’t needed.

Plots is a very powerful, all be it less than ideally documented, plotting library. It provides a Juila interface under which different back ends, such as pyplot or plotyjs, will generate our images. I will be using the GR backend.

# Install the packages we needPkg.add("Plots")Pkg.add("GR")# Import the Plots library in the current runtimeusing Plots# Initialize the gr backendgr()

Alright, now the only thing left is to chose our c term, generate the julia set and plot it.

# Set the size of our imagesheight = 800width = 800

c = 1.2e^(1.1π*im)

# Generate the setdata = julia_set(height, width, c)

# Generate the heatmap and save it as a PNGfile_name = "/tmp/julia_$(height)_$(width)_$(c.re)_$(c.im)"png(heatmap(data, size=(width,height), color=:ice), file_name)

This should be enough for us to generate a pretty fractal, however, I think we should go further.

Animating our fractal

Another nice thing we can do with the Plots package is create animations, by simply putting multiple pngs back to back via ffmpeg.

This will be done using the @animate macro to create an animation from our for loop. Then, the gif function to save it to a file and give it our desired frame rate.

Now, depending on the value of c, framerate, color scheme and resolution you select, you could get wildly different animations here. I’ve create one that I’m quite happy with.

Making it parallel

Generating the images and computing the sets could take quite a long while. However, Julia provides tools to easily parallelize for loops without side effects. This can be done for the loop that generates our sets using the @parallel macro.

One thing that should be note here, is that you need to run Julia with the -p argument to benefit from this. For example julia -p 10, uses 9 additional processes.

Once you do this, functions need to be declared with the @everywhere macro to be accessible from all processes.

@everywhere julia = julia@everywhere julia_set = julia_set# Re-write our code to generate the sets in parallelnr_frames = 40sets = @sync @parallel append! for n=0:nr_frames    c = e ^ (2*π*n/nr_frames*im)    set = julia_set(height, width, c)    [set]endanim = @animate for set in sets    heatmap(set, size=(width,height), color=:speed, leg=false)end

file_name = "/tmp/julia_$(height)_$(width)_$(nr_frames).gif"gif(anim, file_name, fps=4)

Parallelizing the animation may seem trivial, but, due to the way @animate works, it’s actually a bit more complicated. I will leave that as an exercise, to the reader.

I’d encourage anyone interested to try using this language for some of their scientific computing needs.

I’ve come to prefer Julia over Python in a very short while and I’m convinced some of my readers may feel the same way

Note: If you liked this article I have a followup article on learning Julia (this one deals with *hype noises* machine learning): Adversarially generated Julia sets