Category Archives: Julia

Everyone’s Favorite Blogpost: CSV Benchmarks

By: Jacob Quinn

Re-posted from: https://quinnj.home.blog/2019/08/24/everyones-favorite-blogpost-csv-benchmarks/

We all know ’em, we all hate ’em, let’s get a good benchmarking blogpost up in here. Most people who groggily glance at their phone at 5:30 AM just roll over and go back to sleep. For some of us, you also open up the email just to see if there’s anything interesting really quick. And for the very small minority out there, we have a “performance issue” opened on one of our precious darling open-source libraries and THAT’S IT! no more sleep until this internet stranger can be proven wrong! (for the record, the issue in question is here and xiaodiagh isn’t an internet stranger, but a great buddy I got to meet at JuliaCon 2019 in Baltimore this year who is doing some really cool work on grouping performance in Julia; but you know, “internet stranger” is a lot funnier).

Some of you may heard/seen that multithreaded csv parsing support recently landed in CSV.jl, my aforementioned precious darling. So naturally, it’s a good time to round up some benchmarks and show how competitive we are in the csv parsing landscape. I apologize for the lack of fancy graphics and pretty charts, but I’m more interested in the numbers; pretty chart PRs welcome to CSV.jl!

Hey! Look at that! CSV.jl is basically on par with R’s fread!

Now, I’ll add my personal opinion on these kind of benchmark comparisons: always take them with a grain of salt. Benchmarking tends to rely on contrived data that can sometimes be biased one way or another; it can be system dependent in a bunch of ways, they are only accurate for a given amount of time while packages continue to develop/decay, caveat, caveat, caveat. BUT, they also tend to be directionally accurate, and that’s what I’m most pleased with here, particularly with regards to fread‘s one remaining advantage over CSV.jl in terms of multithreading support. (For full feature comparison between CSV.jl, fread, and pandas, see the CSV.jl 0.5 release announcement).

The files + benchmark script can be found here; the single-typed files and mixed.csv are derived from this benchmark site. The benchmark numbers shown above were run on my system: 2019 MacBook Pro, 2.4 GHz Intel Core i9, 32 GB 2400 MHz DDR4 RAM. The benchmarks were run using CSV.jl#master branch, as the last few things get ironed out before a release which will include multithreading support for Julia versions 1.3+.

As always, hit me up on Twitter to chat or follow my #JuliaLang tweets.

GSoC’19: Duckietown.jl Summary

By: Tejan Karmali

Re-posted from: https://tejank10.github.io/jekyll/update/2019/08/24/GSoC-2019-Duckietown.jl-Summary.html

Hello there,

Over the past year, I continued my streak with Julia by contributing to some interesting experiments with differentiable programming. That got me super-excited about the paradigm of differentiable learning. The main idea being that if we know the system, it could be used to simplify and accelerate the training process. Together with Mike and Avik, we planned on a mission: a self-driving car simulator using differentiable programming.

We chose duckietown environment to test our approach. Duckietown is a project started by Prof. Liam Paull. It is a miniature model of a town having buildings, vehicles, traffic signals, and pedestrians. Maxime Chevalier-Boisvert et al, from MILA, have built an awesome simulator of the duckietown, called gym-duckietown. It is used for testing algorithms before deploying it in real duckietown environment. But since it is written using python, it is not differentiable. Hence to make it differentiable we had to build it in julia.

Duckietown.jl creation is spread over three parts:

  • The Simulator
  • Rendering with RayTracer
  • Training

Let’s get started!

The Environment

Duckietown environment contains maps for different tasks – straight road, loop, zigzag turns, UdeM, etc. There are also some variants of these maps with dynamic objects, like traffic signal and pedestrians. The maps are encoded in a .yaml files and can be parsed with YAML.jl. Each environment contains a map in the form of a grid. Each element of the grid is allocated a tile: for eg., a road, an asphalt, an office floor or it can be a grassy surface. A logical arrangement of these tiles is all that is required for a bare-bones version of the map to set up. That’s all the straight road and loopy road maps have. To make these maps more challenging, we need to add objects to it. Objects can be static or dynamic. Static objects include house, tree, traffic sign, bus, truck, traffic cone, etc. Dynamic objects are traffic signals and duckies which are the pedestrians in the duckietown. The positions of these objects are defined in the map. An object is represented as meshes, with texture wrapped around it.

UdeM

The Simulator

using Duckietown, Flux, Zygote

sim = Simulator(map_name="straight_road", camera_width = 128, camera_height = 128)

Simulator manages the subtasks involved in running the duckietown and maintains related statistics. The subtasks include updating the states and positions of different objects involved, running an action on duckiebot, maintaining data such as velocity, position, the action performed on the duckiebot, rendering what the bit sees, etc. The parameters of the simulator are defined in a FixedParams object. These are the parameters that define the behavior of the simulator.

Rendering the view

rende_obs(...) is used to render was duckiebot sees. We use differentiable RayTracer.jl for this purpose. For rendering, we first need to define a camera model. The camera needs to know where the bot is looking from & at, dimensions of the image, field of view, focal length, and up vector.

x, y, z = sim.cur_pos
# get the direction in which bot is looking
dx, dy, dz = get_dir_vec(sim.cur_angle)

## Define camera model
# Looking from
eye = Vec3([x], [y], [z])
# Looking at
target = Vec3([x + dx], [y + dy], [z + dz])
# vup is vector pointing in upward direction
vup = Vec3([0f0], [1f0], [0f0])
cam = Camera(eye, target, vup, cam_fov_y, focal_length, cam_width, cam_height)

A scene is generated containing all the objects in the environments. These objects are decomposed into triangles.

## Scene generation
scene = Vector{Triangle}()
# Decompose the objects into triangles
obj_Δs = map(obj->render(obj, fp.draw_bbox), objs)

for  in obj_Δs
    scene = vcat(scene, )
end

Light source and its position is defined, which is then used to raytrace the scene.

# Define light source
light_pos = Vec3([-40f0], [200f0], [100f0])
# PointLight takes color, intensity and position of light source as args
light = PointLight(Vec3([1f0]), 5f15, light_pos)
origin, direction = get_primary_rays(cam)

# Rendering what duckiebot sees
im = raytrace(origin, direction, observation, light, origin, 2)

Taking action

step!(...) is used to take action on the duckiebot. action is a vector of length 2. It specifies the speed of the left and the right wheel. each element belongs to [-1, 1], where positive speed implies moving in the forward direction. Velocities of both the wheels give us the information about the steering direction of the robot. For eg: to move on a straight road both the velocities should be equal, whereas to take a left turn velocity of the left wheel should be less than that of the right wheel. Using this, the robot’s new position and direction is calculated.

A path to be followed is determined by bezier curves. Each tile has its bezier curve defined. For example, a straight road tile would have a straight line as its curve whereas that for a left turn would be approximately circular. Based on this curve, two kinds of rewards are defined. The distance from the curve and angular distance from the tangent of the curve. There is also a penalty to prevent a collision. Each object has a safety circle surrounding itself. Collision penalty is the degree of overlap between the safety circle of the bot and that of an object.

With these details, we are now equipped to train a model!

Training a model

We define a simple Flux model, which extracts features from the image using Conv and passes them onto the FC layers.

# model: Input- Rendering of what duckiebot sees
#        Output- Action to be taken
model = Chain(
           Conv((3, 3), 3=>8, relu, pad = 1),
           MaxPool((2, 2)),
           Conv((3, 3), 8=>16, relu, pad = 1),
           MaxPool((2, 2)),
           Conv((3, 3), 16=>32, relu, pad = 1),
           x -> reshape(x, :, 1),
           Dense((32 * 32 * 32), 64, relu),
           Dense(64, 16, relu),
           Dense(16, 2),
           x -> reshape(x, 2))

opt = ADAM(0.001f0)

We begin with dividing an episode into sequences. Let’s call a sequence as μEpisode. In each μEpisode, actions are performed for a short number of timesteps. We take the loss as negative of reward and add an action penalty. Recall that actions should lie in [-1, 1]. Since for initial few timesteps actions could arbitrarily lie anywhere in the real domain, this penalty is required. Also, the reward is proportional to speed. If Very high action is chosen, then it should also set very high speed and in turn very high reward, which is not expected ideally. Action penalty is somewhat similar to the regularisation loss.

function μEpisode(model, sim, initial_render, μEp_len)
    obs, action, reward, done, info = step!(sim, model(initial_render))
    loss = -reward + action_penalty(action)

    done && return loss

    for iter in 2:μEp_len
        obs, action, reward, done, info = step!(sim, model(obs))
        loss += -reward + action_penalty(action)
        done && return loss
    end

    return loss
end

In the episode!(...) function, the gradient of the loss wrt to the parameters of the model is taken using Zygote.jl a source-to-source AD package. Gradients are clamped to prevent the overflow due to gradient explosion.

function episode!(sim)    
    for μEp in 1:NUM_μEPISODES
        # Get the gradients of μEpisode
        initial_render = render_obs(sim)
        gs = Zygote.gradient(params(model)) do
            μEpisode(model, sim, initial_render, μEPISODE_LENGTH)
        end

        # Update the weights
        for p in params(model)
            clamp!(gs[p], -0.01f0, 0.01f0)
            Flux.Optimise.update!(opt, p, gs[p])
        end

        sim.done && break
    end

    reset!(sim)
end

And after a while, you should be able to see the bot guiding itself on the lane!
Straight road

What’s next?

What a productive summer it was! With Duckietown.jl you can now research autonomous driving in Julia, and also leverage the differentiability aspect of it. I believe this is just a start for differentiable programming. By knowing the system, we can speed up the training of a model on it by leaps and bounds. In the future, I plan to:

  • Transfer learning: Evaluating the performance model trained on one map by testing it on other maps.
  • Defining tasks over different maps
  • There has been some advances in terms of the packages for physical environements for deep learning. I plan to do some experiments on that using diffferentiable programming.

Acknowledgments

I am extremely grateful to my mentor Mike Innes for posing faith in me for this ambitious project. A huge thanks to my fellow GSoC’er Avik Pal for his amazing RayTracer, and helping me out from time to time. I would also like to thank Dhairya Gandhi for his valuable inputs, Julia Computing Bengaluru for hosting me, and Julia Computing for providing machines for training. Finally, I thank Google for providing me this amazing opportunity in being part of the mission to drive open-source culture.s

B-splines

About three month ago I published my post about Bézier curves and mentioned a follow up post about B-splines. Here it is 😉
I finished the course on Geometric Modelling and Animation in university so I thought it’s time now. Maybe there will be also a post on surfaces.

The blog post contains:

  • An explanation of B-splines
  • How to draw B-splines with different techniques
    • Normal way
    • De Boor Algorithm
  • Animations using Julia

What are B-splines?

The idea of splines in general is to combine several curves together to obtain one single curve. The curve itself should have some nice properties like continuity and it shouldn’t behave in a strange way like making weird curves as a polynomial function of degree 15 might do. Another important point is local control which means if we move one control point a little bit we don’t want to change the whole curve but only the small area around the point. For Bézier curves for example the whole curve can be changed even though the curve is changed mostly around the point that moved.

Additionally having a convex hull property is quite useful for collision detection as checking directly against the curve is quite complex which means it takes time whereas checking whether an object is inside a convex hull is easy. Later the tests can be refined if it at least has the possibility of being part of a collision with the actual curve.

In comparison to Bézier curves we gain local control, a better convex hull property and we can control the degree of the curve by combining curves together. Compared to Bézier splines we don’t have to worry about A-frames anymore which was the last part of the previous post.

Let’s have a look at a specific example:

Bezier spline vs B-spline

I have drawn this using my Bézier curve code of the last post and simple \(C^1\) continuity but for \(C^2\) this would be harder. In this example you can see the normal bezier points \(b_i\) and the so called de Boor points \(d_i\) which given the knot vector which I roughly explained last time (and the continuity) perfectly define the curve and the inner points \(b_2,b_4, b_6\) are defined implicitly by this.

I’ll come back to the Knot vector later which defines the continuity and whether it is interpolating or approximating the point. For the comparison to Bézier splines it should be end point interpolating and approximating everywhere else which means that the curve touches \(d_0\) and \(d_{end}\) but not the ones in between.

I think we need a bit more definitions for the next steps:

We have \(m+1\) de Boor points, \(d_0,\dots,d_m\) which means in the example above \(m=5\)
and degree \(n\) and the Knot Vector \(\mathbf{T} = (u_0, \dots, u_{n+m+1})\)

The spline is defined by:

\[
s(t) = \sum_{i=0}^{m} \mathbf{d}_iN_i^n(t)\]

where \(N_i^n(t)\) are the B-spline basis functions which are defined by:

\[
N_i^0(t) = \begin{cases}
1 & \text{if } u_i \leq t

\[
N_i^r(t) = \frac{t-u_i}{u_{i+r}-u_i}N_i^{r-1}(t) +
\frac{u_{i+1+r}-t}{u_{i+1+r}-u_{i+1}} N_{i+1}^{r-1}(t) \quad 1…