When I was reading about visualizing 2.7 billion points in a minute on a , I got curious how fast this would be in . this blog post macbook Julia Not only for fun, but also to integrate it into my plotting library ( ). Since I’ve been very happy at how quickly I was able to create a very fast solution, I decided to share my experience! MakiE To get started, one can download the data from http://planet.osm.org/gps/ Parsing the CSV I'm assuming that one wants to do lots of operations on the data, so the most sensible thing to do is to convert the 55gb CSV to a 22gb binary blob, which can get and doesn't need parsing! memory mapped There wasn't any tool around that could do this out of the box, but with the help of .jl I was able to whip up a fairly fast custom solution - which also demonstrates how nicely one can extend existing Julia libraries. TextParser The first problem is, that TextParser expects a string to parse the CSV. Since I can't just load the whole 55gb dataset into RAM, our first trick is to create a memory mapped string type in Julia: struct MMapString{T}# using Type parameter to not write out the result type of mmap# but I also don't want to loose performancex::Tend # Constructor from a pathMMapString(path::String) = MMapString(Mmap.mmap(open(path), Vector{UInt8}))# we only need to overload the iteration protocol for TextParse to work!Base.getindex(x::MMapString, i) = Char(x.x[i])Base.next(x::MMapString, i) = Char(x.x[i]), i + 1Base.done(x::MMapString, i) = i > length(x.x)Base.start(x::MMapString, i) = 1Base.length(x::MMapString) = length(x.x) This is all one needs to satisfy the basic interface of most libraries expecting a string like type. Now we can write a function to read the 55gb dataset, parse it line by line and write it into a binary blob. One of the problems that occurred was, that the function to write to an IO file stream was allocating quite a bit. creates a mutable container (Ref) for its argument to safely pass it to C. Even that directly takes a pointer still has a check that allocates 128 bytes. Be ready for , short . This saves us 40 seconds and allocates 3.5 gb less memory (which ends up as ~1.2x faster)!This is actually a great example of how Julia currently fails, and how one can still claw back the performance.We just call the C library directly and use the almost deprecated syntax (`&`) to take a pointer from the stack allocated argument.in Julia 0.7 this is likely no issue and we should be able to just use ` `, since the compiler got a lot better at stack allocating and eliminating mutable containers. Base.write Base.unsafe_write unsafe_unsafe_write uuwrite write By the way I was figuring this out with which shows the allocations, and then using (or alternatively to jump directly into the editor) to figure out where the function is defined in Julia and where the allocations come from. @time write(io, 1f0) @which write(io, 1f0) @edit using TextParseusing TextParse: Record, Field, Numeric, tryparsenext function uuwrite(io, ptr::T) where TInt(ccall(:ios_write,Csize_t, (Ptr{Void}, Ptr{T}, Csize_t),io.ios, &ptr, sizeof(T)))end function save2bin(path, n = typemax(Int))str = MMapString(path)# Descriptor of 'num,num\n' which is the format in the csvrec = Record((Field(Numeric(Float32), delim = ','),Field(Numeric(Float32), eoldelim = true)))# skip the header! Nice is, that Julia's findfirst works with# any iteratorpos = findfirst(str, '\n') + 1io = open(homedir()*"/gpspoints.bin", "w")i = 0while !done(str, pos) && i <= n# tryparsenext takes roughly 35ns so it's fairly optimalp_or_null, pos = tryparsenext(rec, str, pos, length(str))isnull(p_or_null) && continue # continue when parsing failsp = get(p_or_null)uuwrite(io, p)i += 1endclose(io)iendtic() save2bin(homedir() * "/Downloads/gps-points.csv");toc() @time This rewrites the file and saves it as a binary in 190 seconds, which isn't too bad for such a simple implementation and an operation that only needs to be done one time. Drawing the Image Now we can memory map the file as a which holds the raw gps coordinates. We can then directly loop over the points and accumulate each position in a 2D array. Vector{NTuple{2, Float32}} """Transforms from longitude/latitude to pixel on screen, with `dims` refering tothe dimensions of the screen in pixel""" function gps2pixel(point, dims)lon, lat = point[1], point[2]x = ((dims[1] / 180.0) * (90.0 + (lon / 10^7)))y = ((dims[2] / 360.0) * (180.0 - (lat / 10^7)))(x, y)end @inline function to_image_inner!(img, points, start, stop)dims = size(img)for i = start:stop# don't check bounds when indexing into an array beginp0 = gps2pixel(points[i], dims)idx = Int.(round.(p0))xidx, yidx = dims[1] - idx[1], dims[2] - idx[2]if checkbounds(Bool, img, xidx, yidx)# we should give each point a radius and then add# the coverage to each pixel for a smoother image# But this does well enough for this short example:img[xidx, yidx] += 0.001f0 # accumulateendendendendfunction to_image!(img, points, range)N = length(range)# Number of threads Julia has available# Use the environment variable JULIA_NUM_THREADS=8 to start# Julia with 8 threadsNT = Threads.nthreads()slices = floor(Int, N / NT)offset = minimum(range)Threads. for i = 1:NT# creates a closure, which sometimes# confuses type inference leading to slow code.# this is why it's a good practice to move the loop body# behind a function barrier# ( )to_image_inner!(img, points,offset + ((i-1) * slices + 1), offset + (i * slices))endreturn imgend @inbounds @threads @threads https://docs.julialang.org/en/latest/manual/performance-tips/#kernal-functions-1 We can use the above function in the following way to get a simple gray scale image: img = zeros(Float32, 600, 960)io = open(homedir() * "/gpspoints.bin")points = Mmap.mmap(io, Vector{NTuple{2, Float32}})tic()to_image!(img, points, 1:length(points))toc()using FileIO, Colors # save it with FileIO as a Gray imageFileIO.save("gps.png", Gray.(clamp.(1f0 .- img, 0, 1))) Et voila! This took 10 seconds to draw 22gb of points on a normal desktop computer. One thing to point out is the great dot call syntax you can see in the example: Gray.(clamp.(1f0 .- img, 0, 1))) This is actually a wonderful mechanism to . So the above will turn into a single for loopover , clamping, inverting and converting it to a color in one go! apply a function per element on an array while fusing all operations img We can also use my new to create a nice animation of the progress: interactive plotting library using MakiE resolution = (600, 960)scene = Scene(resolution = reverse(resolution))img = zeros(Float32, resolution)imviz = heatmap(img, colornorm = (0, 1))center!(scene)# slice `points` into multiple batches which we will updateslice = 10^7range = slice:slice:length(points)stop = 1vio = VideoStream(scene, homedir()*"/Desktop/", "gps")while truestart = stopstop = min(stop + slice, length(points))to_image!(img, points, start:stop)imviz[:heatmap] = img # update heatmap plotrecordframe!(vio) # add frame to streamstop == length(points) && breakendfinish(vio, "gif") # finish streaming and export as gif! There are of course various improvements we could apply, but I’m very happy with this as a first result. It is also not clear how comparable this solution is to the initially mentioned blog post, since I haven’t run that on my computer yet. If you want to try this on your machine, check out the full ! code
