Category Archives: Julia

How much do collections of allocated objects cost?

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2021/06/11/vecvec.html

Introduction

In my recent comment on StackOverflow I have said that using a vector of
vectors is: a) slow, b) uses more memory than needed, and c) puts much more
stress on garbage collection. Having read it one of my students asked me to
expand on this issue. In this post I want to give you some examples that were
designed to clarify this issue.

The post was tested under Julia 1.6.1 on Linux with a machine having 16 GB of
RAM (the last point affects the frequency of triggering GC).

In the post we will compare the same operations for vector of vectors and
vector of tuples defined using the following functions:

test1(n) =[rand(2) for _ in 1:n]
test2(n) =[(rand(), rand()) for _ in 1:n]

The difference is that test1 creates a vector of references to dynamically
allocated objects, while in test2 the tuples are stored directly in the vector.

Let us test all three claims. In the comparisons I use a fresh session in
each code block (as the examples given use a lot of memory). This means that
the timings will include compilation time, but the size of the computations is
large enough so that this is relatively negligible.

Examples

Start with a vector of vectors:

julia> test1(n) =[rand(2) for _ in 1:n]
test1 (generic function with 1 method)

julia> @time t1 = test1(10^8);
 21.907014 seconds (100.00 M allocations: 9.686 GiB, 67.93% gc time)

julia> @time sum(x -> x[1], t1);
  0.768779 seconds (179.18 k allocations: 11.105 MiB, 7.44% compilation time)

julia> @time Base.summarysize(test1(10^7))
262.832256 seconds (50.06 M allocations: 2.864 GiB, 4.58% gc time, 0.02% compilation time)
640000040

julia> @time GC.gc()
  4.631474 seconds (100.00% gc time)

julia> @time GC.gc()
  2.647404 seconds (100.00% gc time)

julia> @time GC.gc(false)
  0.004821 seconds (99.89% gc time)

julia> @time GC.gc(false)
  0.005526 seconds (99.89% gc time)

And now vector of tuples (fresh Julia session):

julia> test2(n) =[(rand(), rand()) for _ in 1:n]
test2 (generic function with 1 method)

julia> @time t2 = test2(10^8);
  1.458052 seconds (25 allocations: 1.490 GiB, 0.36% gc time)

julia> @time sum(x -> x[1], t2);
  0.164072 seconds (178.23 k allocations: 11.036 MiB, 35.07% compilation time)

julia> @time Base.summarysize(test2(10^7))
  0.190496 seconds (24.27 k allocations: 153.937 MiB, 3.25% gc time, 17.60% compilation time)
160000040

julia> @time GC.gc()
  0.070696 seconds (99.99% gc time)

julia> @time GC.gc()
  0.102222 seconds (99.99% gc time)

julia> @time GC.gc(false)
  0.000517 seconds (99.13% gc time)

julia> @time GC.gc(false)
  0.000523 seconds (98.97% gc time)

As you can see:

  • creation of vector of vectors is much slower; in particular a lot of small
    allocations happens (which is expensive) and in total also more memory is
    allocated.
  • a simple aggregation with sum is also slower because for vector of vectors
    we have to go through references to objects (which takes time) which also
    means that this is less CPU cache friendly.
  • With Base.summarysize we can check that using vectors also uses up much more
    memory; also as a side issue we learn that functions like Base.summarysize
    which traverse the tree of object references are much, much slower for a
    vector of vectors.
  • Finally both full sweep GC.gc() and incremental sweep GC.gc(false) are
    slower with vector of vectors; this is especially visible for full sweep case
    (fortunately it is triggered less often in normal usage). The important thing
    to note here is that for garbage collection time to be affected it is enough
    that the vector of vectors is somewhere in the memory; it does not have to be
    used in some operation you do.

Conclusions

The conclusion is something that seasoned Julia developers know very well:
avoid having many small allocated objects in your Julia programs. Having read
this post I hope you now have a better understanding what aspects of the
performance of your code can be affected when you have to use such data.

Before I finish let me add that collections of any mutable objects will be
affected by the same issue, and even some special immutable objects (like
String, or to some extent e.g. Symbol) can cause issues like presented
in this post.

Calling R From Julia

By: Josh Day

Re-posted from: https://www.juliafordatascience.com/calling-r-from-julia/

Calling R From Julia

Since Julia has a younger ecosystem, you won't always find the functionality you need for every task.  In those cases, you're left with the choice to either:

  • Implement it yourself.
    • This is great for the Julia community in the long run (especially if you release your work as a package! Please do this!), but you don't always have the time/energy to do so.
  • Use a package in another language.
    • Language wars are boring. You should use all the tools at your disposable. You don't need to stick with Julia for everything (especially since interop is easy!).

Particularly for niche models in the field of statistics, you'll often find R packages that do not have Julia counterparts.  Thankfully, Julia and R work together rather seamlessly through the help of one package.

Introducing RCall.jl

To get started, let's install RCall:

] add RCall

Running R Code from Julia

  • After installing and typing using RCall, you'll have access to the R REPL Mode by typing $.  Your prompt will change from julia> to R>  and now all of your commands will run in R instead of Julia.  You can use it just like a normal R session.
julia> using RCall

# type `$`

R> install.packages("ggplot2")

R> library(ggplot2)

R> data(diamonds)

R> ggplot(diamonds, aes(x=carat, y=price)) + geom_point()
Calling R From Julia
  • If you want to call R from Julia in a non-interactive manner (not from the REPL), you can use @R_str macro:
julia> R"y = 2"
RObject{RealSxp}
[1] 2

Sending Julia Variables to R

  • The @rput macro sends a variable to R and uses the same name.
julia> x = 1
1

julia> @rput x
1

R> x
[1] 1
  • You can interpolate Julia values in @R_str commands as well as the R REPL Mode:
julia> x = 1
1

julia> R"y = $x"
RObject{IntSxp}
[1] 1

R> 1 + $x
[1] 2

Retrieving R Variables in Julia

  • The @rget macro sends a variable from R to Julia and uses the same name.
julia> R"z = 5"
RObject{RealSxp}
[1] 5


julia> @rget z
5.0

julia> z
5.0
  • You can also convert an RObject into the appropriate Julia counterpart with rcopy.
julia> robj = R"z"
RObject{RealSxp}
[1] 5


julia> rcopy(robj)
5.0

🚀 That's It!

You now know the basics of working in R from Julia.  There are deeper depths to dive into on the topic (such as type conversions between R and Julia), but this should give you a start on using your favorite R package together with Julia.

Enjoying Julia For Data Science?  Please share us with a friend and follow us on Twitter at @JuliaForDataSci.

Additional Resources

CUDA.jl 3.3

By: Tim Besard

Re-posted from: https://juliagpu.org/post/2021-06-10-cuda_3.3/index.html

There have been several releases of CUDA.jl in the past couple of months, with many bugfixes and many exciting new features to improve GPU programming in Julia: CuArray now supports isbits Unions, CUDA.jl can emit debug info for use with NVIDIA tools, and changes to the compiler make it even easier to use the latest version of the CUDA toolkit.

CuArray support for isbits Unions

Unions are a way to represent values of one type or another, e.g., a value that can be an integer or a floating point. If all possible element types of a Union are so-called bitstypes, which can be stored contiguously in memory, the Union of these types can be stored contiguously too. This kind of optimization is implemented by the Array type, which can store such "isbits Unions" inline, as opposed to storing a pointer to a heap-allocated box. For more details, refer to the Julia documentation.

With CUDA.jl 3.3, the CuArray GPU array type now supports this optimization too. That means you can safely allocate CuArrays with isbits union element types and perform GPU-accelerated operations on then:

julia> a = CuArray([1, nothing, 3])
3-element CuArray{Union{Nothing, Int64}, 1}:
 1
  nothing
 3julia> findfirst(isnothing, a)
2

It is also safe to pass these CuArrays to a kernel and use unions there:

julia> function kernel(a)
         i = threadIdx().x
         if a[i] !== nothing
           a[i] += 1
         end
         return
       endjulia> @cuda threads=3 kernel(a)julia> a
3-element CuArray{Union{Nothing, Int64}, 1}:
 2
  nothing
 4

This feature is especially valuable to represent missing values, and is an important step towards GPU support for DataFrames.jl.

Debug and location information

Another noteworthy addition is the support for emitting debug and location information. The debug level, set by passing -g <level> to the julia executable, determines how much info is emitted. The default of level 1 only enables location information instructions which should not impact performance. Passing -g0 disables this, while passing -g2 also enables the output of DWARF debug information and compiles in debug mode.

Location information is useful for a variety of reasons. Many tools, like the NVIDIA profilers, use it corelate instructions to source code:

NVIDIA Visual Profiler with source-code location information

Debug information can be used to debug compiled code using cuda-gdb:

$ cuda-gdb --args julia -g2 examples/vadd.jl
(cuda-gdb) set cuda break_on_launch all
(cuda-gdb) run
[Switching focus to CUDA kernel 0, grid 1, block (0,0,0), thread (0,0,0), device 0, sm 0, warp 0, lane 0]
macro expansion () at .julia/packages/LLVM/hHQuD/src/interop/base.jl:74
74                  Base.llvmcall(($ir,$fn), $rettyp, $argtyp, $(args.args...))(cuda-gdb) bt
#0  macro expansion () at .julia/packages/LLVM/hHQuD/src/interop/base.jl:74
#1  macro expansion () at .julia/dev/CUDA/src/device/intrinsics/indexing.jl:6
#2  _index () at .julia/dev/CUDA/src/device/intrinsics/indexing.jl:6
#3  blockIdx_x () at .julia/dev/CUDA/src/device/intrinsics/indexing.jl:56
#4  blockIdx () at .julia/dev/CUDA/src/device/intrinsics/indexing.jl:76
#5  julia_vadd<<<(1,1,1),(12,1,1)>>> (a=..., b=..., c=...) at .julia/dev/CUDA/examples/vadd.jl:6(cuda-gdb) f 5
#5  julia_vadd<<<(1,1,1),(12,1,1)>>> (a=..., b=..., c=...) at .julia/dev/CUDA/examples/vadd.jl:6
6           i = (blockIdx().x-1) * blockDim().x + threadIdx().x(cuda-gdb) l
1       using Test
2
3       using CUDA
4
5       function vadd(a, b, c)
6           i = (blockIdx().x-1) * blockDim().x + threadIdx().x
7           c[i] = a[i] + b[i]
8           return
9       end
10

Improved CUDA compatibility support

As always, new CUDA.jl releases come with updated support for the CUDA toolkit. CUDA.jl is now compatible with CUDA 11.3, as well as CUDA 11.3 Update 1. Users don't have to do anything to update to these versions, as CUDA.jl will automatically select and download the latest supported version.

Of course, for CUDA.jl to use the latest versions of the CUDA toolkit, a sufficiently recent version of the NVIDIA driver is required. Before CUDA 11.0, the driver's CUDA compatibility was a strict lower bound, and every minor CUDA release required a driver update. CUDA 11.0 comes with an enhanced compatibility option that follows semantic versioning, e.g., CUDA 11.3 can be used on an NVIDIA driver that only supports up to CUDA 11.0. CUDA.jl now follows semantic versioning when selecting a compatible toolkit, making it easier to use the latest version of the CUDA toolkit in Julia.

For those interested: Implementing semantic versioning required the CUDA.jl compiler to use ptxas instead of the driver's embedded JIT to generate GPU machine code. At the same time, many parts of CUDA.jl still use the CUDA driver APIs, so it's always recommended to keep your NVIDIA driver up-to-date.

High-level graph APIs

To overcome the cost of launching kernels, CUDA makes it possible to build computational graphs, and execute those graphs with less overhead than the underlying operations. In CUDA.jl we provide easy access to the APIs to record and execute these graphs:

A = CUDA.zeros(Int, 1)# ensure the operation is compiled
A .+= 1# capture
graph = capture() do
    A .+= 1
end
@test Array(A) == [1]   # didn't change anything# instantiate and launch
exec = instantiate(graph)
CUDA.launch(exec)
@test Array(A) == [2]# update and instantiate/launch again
graph′ = capture() do
    A .+= 2
end
update(exec, graph′)
CUDA.launch(exec)
@test Array(A) == [4]

This sequence of operations is common enough that we provide a high-level @captured macro wraps that automatically records, updates, instantiates and launches the graph:

A = CUDA.zeros(Int, 1)for i in 1:2
    @captured A .+= 1
end
@test Array(A) == [2]

Minor changes and features