My quest to responsive visualizations with Julia

By: Simon Danisch

Re-posted from: http://randomfantasies.com/2016/09/my-quest-to-responsive-visualizations-with-julia/

Responsive and fun interactivity is notoriously hard! But it’s also the key to less frustration and more patience when working on a project.

There are endless important projects that humanity needs to solve. To become less of a burden to the environment and have a better society we constantly need to advance the state of the art.

Still, people give up on doing this when the initial amount of patience, motivation and curiosity is depleted without getting new incentives.

This actually happened to me with my 3D modeling hobby a few years ago. Not only was the 3D software difficult to learn, it also turned unresponsive fairly quickly, even on a powerful machine.

Just as an example, rendering a 30 second clip in high quality took ~2 weeks with my PC running day and night. This slowness drained my patience and at some point I decided that this is too much hassle for a hobby.

Half a year later, I ran into a similar problem when working on a computer vision pipeline. We simply couldn’t find frameworks which were fast enough for real time processing while maintaining high programming productivity. This time I didn’t give up, but the job could have been much more fun and productive!

At this point I decided that I want to work on software that turns really tough problems like object recognition of 3D models into weekend projects.

I started working on GLVisualize, a library for fast and interactive graphics written in Julia.

Why graphics?

Visualizing something is the first step to understanding and it allows us to explore huge problem spaces that were invisible before.

But current tools seem to be divided into 2D and 3D libraries, high performance libraries, which are kind of a hassle to use, and easily usable libraries which turn into an unresponsive mess quickly.

So with my background, this seemed like a worthy start!

Why Julia?

While there are several reasons for choosing Julia, today I want to show you its impressive speed despite being an easy-to-use dynamic language.

Let me show you the benefit with two toy examples (all graphics including the plots were produced with GLVisualize!).

Consider the following function, the famous Lorenz Attractor:

It takes a point t0 and parameters a to d and returns a new point. Could you imagine what kind of line this will draw, when recursively applying it to the new result and changing parameters?

I guess no one has such a powerful brain, which is why we visualize these kind of things, together with sliders to explore the function:

lorenz

code

These are 100,000 points, and as you can see GLVisualize can deal with that amount easily.

Lets see how long we could have stayed at pleasant interactive speeds in another language.

Lets take for example Python, a language comparable in usability:

 

pyjulo

 

minimal speedup 70.4680453961563
maximal speedup 148.02061279172827
max seconds Julia 0.193422334
max seconds Python 13.630093812942505

As you can see, computation times jump beyond one second at around 10⁵ points for Python, while Julia stays interactive with up to 10⁷ points.

So if you’re interested in a higher resolution in Python you better crank up your patience or call out to C, eliminating all convenience of Python abruptly!

Next example is a 3D mandelbulb visualization:

mandelbulb on the gpu

code

One step through the function takes around 24 seconds with Julia on the CPU, which makes it fairly painful to explore the parameters.

So why is the shown animation still smooth? Well, when choosing Julia, I’ve been betting on the fact that it should be fairly simple

to run Julia code on the GPU. This bet is now starting to become reality.

I’m using CUDAnative and GPUArrays to speed up the calculation of the mandelbulb.

CUDAnative allows you to run Julia code on the GPU while GPUArrays offers a simpler array interface for this task.

Here are the benchmark results:

cujlmandel

minimal speedup 37.05708590791309
maximal speedup 401.59165062897495
max seconds cpu 24.275534426
max seconds cuda 0.128474155

This means that by running Julia on the GPU we can still enjoy a smooth interaction with this function.

But the problem is actually so hard, that I can’t run this interactively at the resolution I would like to.

As you can see, the iso surface visualization still looks very coarse. So even when using state of the art software, you still run into problems that won’t compute under one second.

But with Julia you can at least squeeze out the last bit of performance of your hardware, be it CPU or GPU, while enjoying the comfort of a dynamic high level language!

I’m sure that with Julia and advances in hardware, algorithms and optimizations, we can soon crack even the hardest problem with ease!

 


Benchmarking system:

RAM: 32GB

CPU: Intel® Core™ i7-6700 CPU @ 3.40GHz × 8

GPU: GeForce GTX 950

 

An Introduction to Machine Learning in Julia

Introduction

Machine learning is now pervasive in every field of inquiry and has lead to
breakthroughs in various fields from medical diagnoses to online advertising.
Practical machine learning is quite computationally intensive, whether it involves
millions of repetitions of simple mathematical methods such as
Euclidian Distance or more
intricate optimizers
or backpropagation algorithms.
Such computationally intensive techniques need a fast and expressive language –
one that enables scientists to write simple, readable code that
performs well. Enter Julia.

In this post, we introduce a simple machine learning algorithm
called K Nearest Neighbors,
and demonstrate certain Julia features that allow for its easy and efficient
implementation. We will demonstrate that the code we write is inherently
generic, and show the use of the same code to run on GPUs via the
ArrayFire package.

Problem statement

We need to come up with a character recognition software that can identify images
of textual characters. The data we’ll be using for this problem are a
set of images from Google Street View.
The algorithm we’ll be using is called k-Nearest Neighbors (kNN), which is a
useful starting point for working on problems of this nature.

Solution

Let us begin by taking a good look at our data. First, we read the data from a set of
CSV files from disk, and plot a histogram that represents the distribution of the training data:

# Read training set 
trainLabels = readtable("$(path)/trainLabels.csv")
counts=by(trainLabels, :Class, nrow);
plot(x = counts[:Class], y=counts[:x1], color = counts[:Class], Theme(key_position = :none), Guide.xlabel("Characters")

Training data distribution

When developing this code within an integrared environment such as IJulia Notebooks, we can explore this data interactively. For example, we can scroll through the images using the
Interact package:

@manipulate for i = 1:size(trainLabels,1)
    load("$(path)/trainResized/$i.Bmp")
end

kNN

The solution of labelling the image with its textual character involves finding the “distances” of each image in the training
set to every other image. Then, the image under consideration is classified
as having the most common label among its k nearest images.
We use Euclidean distance
to measure the proximity of the images:

Let X be a set of m training samples, Y be the labels for each of those samples, and x be
the unknown sample to be classified. We perform the following steps:

  • First, calculate all the distances:
for i = 1 to m
    distance = d(X(i), x)
end

We can also represent this generally in matrix notation:

function get_all_distances(imageI::AbstractArray, x::AbstractArray)
    diff = imageI .- x
    distances = vec(sum(diff .* diff,1))
end
  • Then sort the distances in ascending order and select first k indices
    from the sorted list.
function get_k_nearest_neighbors(x::AbstractArray, imageI::AbstractArray, k::Int = 3, train = true)
    nRows, nCols = size(x)
    distances = get_all_distances(imageI, x)
    sortedNeighbors = sortperm(distances)
    if train
        return kNearestNeighbors = Array(sortedNeighbors[2:k+1])
    else
        return kNearestNeighbors = Array(sortedNeighbors[1:k])  
    end
end
  • Finally, classify x as the most common element among the k indices.
function assign_label(x::AbstractArray, y::AbstractArray{Int64}, imageI::AbstractArray, k, train::Bool)
    kNearestNeighbors = get_k_nearest_neighbors(x, imageI, k, train)
    counts = Dict{Int, Int}()
    highestCount = 0
    mostPopularLabel = 0
    for n in kNearestNeighbors
        labelOfN = y[n]
        if !haskey(counts, labelOfN)
            counts[labelOfN] = 0
        end
        counts[labelOfN] += 1
        if counts[labelOfN] > highestCount
            highestCount = counts[labelOfN]
            mostPopularLabel = labelOfN
        end
     end
    mostPopularLabel
end

The objective is to find the best value of k such that we reach optimal
accuracy.

Training

We use the labeled data to find the optimal value of k. The following code
will output the accuracy for values of k from 3 to 5, for example.

for k  = 3:5
    checks = trues(size(train_data, 2))
    @time for i = 1:size(train_data, 2)
        imageI = train_data[:,i]
        checks[i]  = (assign_label(train_data, trainLabelsInt, imageI, k, true) == trainLabelsInt[i])
    end
    accuracy = length(find(checks)) / length(checks)
    @show accuracy
end

Testing

Now, from our training, we should have identified a value of k that gives us
the best accuracy. Of course, this varies depending on the data. Suppose our runs
tell us the best value of k is 3. It’s time to now test our model. For this,
we can actually use the same code that we used for training, but we needn’t check
for accuracy or loop over values of k.

x = test_data
xT = train_data
yT = trainLabelsInt
k = 3, #say
prediction = zeros(Int,size(x,2))
@time for i = 1:size(x,2)
    imageI = x[:,i]
    prediction[i]  = assign_label(xT, yT, imageI, k, false)
end

Improving Performance

There are a couple of ways we can speed up our code. Notice that each iteration
of the for-loop in the main compute loop is independent. This looks like a great
candidate for parallelisation. So, we can spawn a couple
of Julia processes and execute the for loop in parallel.

Parallel for-loops

The macro @parallel will split the iteration space amongst the available Julia
processes and then perform the computation in each process simultaneously.

Let us first add n processes by using the addprocs function. We
then execute the following for loop in parallel.

addprocs(n)
@time sumValues = @parallel (+) for i in 1:size(xTrain, 2)
    assign_label(xTrain, yTrain, k, i) == yTrain[i, 1]
end

The following scaling plot shows performance variation as you spawn more processes:

GPU Acceleration

The ArrayFire package enables us to run
computations on various general purpose GPUs from within Julia programs.
The AFArray constructor from this package
converts Julia vectors to ArrayFire Arrays, which can be transferred to the GPU.
Converting Arrays to AFArrays allows us to trivially use the same code as above to run our
computation on the GPU.

checks = trues(AFArray{Bool}, size(train_data, 2))
train_data = AFArray(train_data)
trainLabelsInt = AFArray(trainLabelsInt)
for k  = 3:5
    @time for i = 1:size(train_data, 2)
        imageI = train_data[:,i]
        checks[i]  = (assign_label(train_data, trainLabelsInt, imageI, k, true) == trainLabelsInt[i])
    end
    accuracy = length(find(checks)) / length(checks)
    @show accuracy
end

Benchmarking against the serial version shows a significant speedup, which has been achieved without much extra code!

Serial Time ArrayFire Time
34 sec 20 sec

Conclusion

  • Julia allows for easy prototyping and deployment of machine learning models.
  • Parallelization also very easy to implement.
  • Built-in parallelism gives the programmer a variety of options to speed up their code.
  • Generic code can be run on GPUs using the package ArrayFire

Future Work

Congratulations! You now know how to write a simple image recognition model using
kNN. While this represents a good starting point for problems such as these, it is by
no means the state of the art. We shall explore more sophisticated methods, such as deep
neural networks
, in a future blog post.

Julia calling C: A minimal example

This blog is a “Hello World” example of Julia calling C.

We start of by at bit of C code we want to call from Julia. We write the following in calc_mean.c

double mean(double a, double b) {
  return (a+b) / 2;
}

To build the library, we need to create a Makefile

CC=gcc 
 
CFLAGS=-c -Wall -fPIC
 
SOURCES=calc_mean.c 
OBJECTS=$(SOURCES:.c=.o)
 
.c.o:
    $(CC) $(CFLAGS) $< -o $@ 
 
lib: $(OBJECTS)
    $(CC) -shared -fPIC -o libmean.so $(OBJECTS)
 
clean:
    rm *.o *.so

The option fPIC and -shared are essential for Julia to be able to resolve the function in our library. Now we are almost ready to build our library. From the bash terminal we invoke:

make lib

This will generate a libmean.so file.

In Julia we call the function in our c library by

x=ccall((:mean,"libmean"),Float64,(Float64,Float64),2.0,5.0)
println(x)
3.5

For this to work,

  • Julia must be running either on the same path where libmean.so resides,
  • the path to libmean.so is in LD_LIBRARY_PATH, or
  • the path to the library is pushed to Libdl.DL_LOAD_PATH via

push!(Libdl.DL_LOAD_PATH,"path_to_libmean.so")

P.S. Thanks to Christopher Rackauckas for tips on Julia highlighting.