Category Archives: Julia

Drawing Vector Graphics with Julia can be Awesome

Re-posted from: https://medium.com/coffee-in-a-klein-bottle/drawing-vector-graphics-with-julia-can-be-awesome-298169609032?source=rss-8bd6ec95ab58------2

A quick tutorial on Luxor.jl, the go-to package for static vector graphics

We are all very used to plotting packages, yet, from time to time, we want to create visualizations that require more freedom. When that happens, it might be better to go straight to vector graphics. This is an area where people tend to be less experienced, and it’s not always clear what package to use.

Vector Graphics Packages

When it comes to vector graphics, you’ve probably heard the name Cairo, since it is perhaps on of the most famous 2D graphics libraries, and it’s commonly used by plotting packages as a backend. Thus, a first option for drawing your vector graphics is using Cairo.jl, which is Julia’s api for the Cairo library. Yet, it might not be a very pleasant experience, as Cairo can be fairly complicated for newcomers.

Here is where Luxor.jl comes in.

Luxor is a Julia package for drawing simple static vector graphics. It provides basic drawing functions and utilities for working with shapes, polygons, clipping masks, PNG and SVG images, turtle graphics, and simple animations.

The focus of Luxor is on simplicity and ease of use: it should be easier to use than plain Cairo.jl, with shorter names, fewer underscores, default contexts, and simplified functions. — Luxor’s documentation.

I could actually stop here and tell you to go read the docs, cause this package is insanely well documented. But since you are already here, I might as well take some time to quickly spill out some of the basics.

Quick Introduction to Luxor.jl

Perhaps the first thing to note about Luxor is that it’s procedural and static, which means that you code a series of drawing commands until you say it’s finished, which will in turn produce either a PNG, SVG, PDF or EPS based on what you defined at the start.

If you are used to dynamic programming (which you probably are, since you are using Julia), then this static way of drawing with Luxor will be a bit strange at first. Yet, there is a flip side. Luxor is fast! What do I mean? Just try it and you’ll see that not only the precompilation time is minimal, but every redrawing is also quite fast (compared specially to plotting packages such as VegaLite.jl and Plotly.jl).

Let’s do our first drawing in Luxor.

using Luxor
Drawing(500, 500, "my-drawing.svg")
origin()
setcolor("red")
circle(Point(0, 0), 100, :fill)
finish()

The code above encapsulates very well the workflow in Luxor. The very first line is just importing the package. In the second line, we create a drawing with size of 500 by 500 pixels, and once we finish the drawing, Luxor will create a file called “my-drawing.svg” in the folder we are working at. All good, let’s proceed to the third line of the code.

You might think that the orgin() line is “starting” the drawing, but what it’s actually doing in altering the origin of the coordinate system to the center of the image… What do I mean? Well, in Luxor the default origin is in the top left of the figure, and the y-axis actually points down. Thus, the origin() function is responsible to altering the location of origin coordinate. And if we call this function without arguments, it places the origin at the very center of the drawing (you might be a bit confused, unless you’ve already used other vector drawing packages, in which case, this convention might actually make sense to you).

Moving on, the next line is setcolor("red") and it sets the current color to red. The important thing to note here is that this is similar to choosing the current color of the drawing pen. Hence, every drawing that we do after calling setcolor("red") will be red, until we change the color of the pen.

Let’s finally draw something in our canvas. And the command circle(Point(0,0),100,:fill) does exactly that. This command is drawing a circle with center at point (0,0) , which is the center of our figure, since we used the command origin() . The next argument sets the radius, which in our case is 100 pixels. Finally, the :fill argument indicates that we must paint the inside of the circle. Again, this is all very common if you have worked with vector drawings before.

We conclude our drawing with the finish() command, which then saves the drawing file. Now, if you want to see the drawing right after finishing it, you can just use the function preview() , which will show the most current finished drawing.

Another Simple Example

Let’s change a bit of our first drawing.

Drawing(200, 200, "my-drawing.png")
background(0.9, 0.2, 0.2, .5)
setcolor("blue")
setline(10)
setdash("dash")
circle(Point(0, 0), 100, :stroke)
finish()

This is very similar to our initial drawing, but we’ve added some new functions and removed the origin() . Note that now our circle is actually centered in the top left corner of our figure, which can now be seen because of the background() function. These numbers inside the background function is just the color specification in hsla.

Another difference is that we used :stroke instead of :fill in the circle function. Thus, instead of filling the circle, it just draws the border. The setline() function alters the thickness of the stroke, and the setdash() changes the line style. There are many more variations in line styles, which you can check in the function docstring.

Luxor for the Mathematically Inclined

All I showed up until now is very basic, and perhaps not very useful, so perhaps this next example will grab your interest, that is, if you have a more mathematical inclination.

Let’s use Luxor to draw some diagrams, which occur in a discipline called Category Theory ?.

using MathTeXEngine

Drawing(200, 200, "my-drawing.png")

# Drawing the borders of our drawing. Note that the `O` is a
# shortcut to Point(0,0).
rect(O,200,200,:stroke)

# Let's change the origin to the center of the drawing
origin()

# We now draw the our diagram
setline(10)
arrow(Point(-40,0),Point(40,0))

# Let's define the font size and write some text
fontsize(15)
text(L"f",Point(0,20), halign = :center)

circle(Point(-50,0),5,:fill)
text(L"ℕ",Point(-50,20),halign = :center)

circle(Point(50,0),5,:fill)
text(L"ℚ",Point(50,20),halign = :center)

#### From here, we are drawing the looping arrows ####

loopx = 30
loopy = 40
adjx = 6
adjy = -2

arrow(
    Point(-50,0) + Point(-adjx,adjy),
    Point(-50,0) + Point(-loopx,-loopy),
    Point(-50,0)+ Point(loopx,-loopy),
    Point(-50,0)+Point(adjx,adjy)
)

arrow(
    Point(50,0) + Point(-adjx,adjy),
    Point(50,0) + Point(-loopx,-loopy),
    Point(50,0) + Point(loopx,-loopy),
    Point(50,0) + Point(adjx,adjy)
)

finish()

Pretty neat, no? You might be a bit overwhelmed by the code above, but just read it through and you will see it’s quite straightforward. Here are some comments on things that might be obscure.

First, you might have noticed that we can write LaTeX in Luxor (I know, it’s awesome)! Just writeL”\sum^n_{i=10}”, i.e. place “L” before the string. The MathTeX is a package that is required in order to render the LaTeX string without actually needing to use a LaTeX version installed in your computer.

The final part is the code is more complicated, but the only thing going on there is that one can pass control points in order to draw curves. Hence, instead of passing only the start and finish of the arrow, I’m actually passing four points in order to make the arrow loop, and I’m also altering the start and finish position so that they do not touch the circle. The loopx and the other defined variables are just some auxiliary variables that I used in order to figure out the control points to give me a nice looking curve.

Final Words and a Hot Tip

There are many other things I could say about Luxor, since the package is quite extensive. Fortunately, as I’ve already said, the docs are very good. Not only it has many examples, but it also has many explanations about the many functionalities.

Yet, the fact that the docs are so comprehensive also makes a bit harder to find some specific things one might be looking for. Hence, here are some hot tips that might come in handy.

If you are working in a notebook environment such as Jupyter, you might wish to store the drawings in variables, instead of saving it to files right away. Also, you might wish to do small drawings, and then place them inside another drawing… Is this possible? YES! Here is how to do both.

d1 = Drawing(200,200,:svg)
origin()
circle(O+Point(30,0),10,:fill)
finish()

d2 = Drawing(200,200,:svg)
origin()
setcolor("blue")
circle(O,10,:fill)
finish()

d = Drawing(200,200,:svg)
placeimage(d1)
placeimage(d2)
finish()

I recommend always working with svg if your drawing is not too heavy, since you’ll always have crystal clear quality. With the code above, you can just run d in a cell in order to visualize the drawing stored in the variable.

Finally, you might wish to see the actual svg specification. Here is a possible solution.

dsvg = String(copy(d.bufferdata))

This will give you the svg string, which you can just save into a mydrawing.svg file.

That’s all. Now go read the docs!

Drawing Vector Graphics with Julia can be Awesome was originally published in Coffee in a Klein Bottle on Medium, where people are continuing the conversation by highlighting and responding to this story.

Triangle frenzy

By: Can Candan's Blog

Re-posted from: https://cancandan.github.io/julia/graphics/cuda/2022/05/07/triangles.html

Suppose we want to draw a batch of images, where each image is made up of randomly positioned and colored triangles, blending, so it will look like this:

triangles

and then find the euclidean distance of each such image to a given target image. Why on earth am I doing this? Well this turns into an interesting optimization problem of finding the closest triangle image and also an excuse to practise Julia. The inspiration is from this repo.

Towards framing this as an optimization problem, we will represent a triangle as a vector of size 10, made up of floating point numbers between 0 and 1. Four numbers of this vector are for the color of the triangle; r,g,b and alpha, and for the three vertices of the triangle we need 6 numbers, (x1,y1), (x2,y2), (x3,y3). Hence, if we want to draw M images, each image having N triangles, we need a matrix of size (10,N,M), which will be our parameters matrix. I want to randomly create such a matrix and min-max scale it along the triangle dimension, by which I mean, for each image, I first find the minimum and maximum of a triangle parameter among the N triangles, and then subtract from the parameter this minimum and then divide the result by the difference between the maximum and the minimum. I want to end up with an array of size (3,w,h,M) for the images, where w is width and h is height, and an array of size M for the distances. Let’s see how fast we can do this.

First order of business is setting this up, note that I am scaling the numbers to the given width and height:

using Images
using Statistics
using Cairo
using BenchmarkTools
using Random: seed!

function prepare()
    seed!(123)
    w,h=128,128
    num_params=10
    num_triangles=50
    num_images=256
    target = rand(Float32, 3, h, w)
    prms = rand(Float32, num_params, num_triangles, num_images);
    prms .= (prms .- minimum(prms, dims=2)) ./ (maximum(prms, dims=2) .- minimum(prms, dims=2))   
    prms[collect(1:2:6),:,:] .*= w
    prms[collect(2:2:6),:,:] .*= h         
    return prms, target, num_images, num_triangles, w, h
end

The first thing that comes to mind is to use a 2d graphics library for drawing, and since the Cairo lib is available, let’s try that, the function below is drawing the triangles on a blank white Cairo canvas, and copying it to the img array at the end:

@views function renderCairo(img, prms, num_triangles, w, h)                  
    # blank white canvas
    buffer = ones(UInt32, w, h) * typemax(UInt32)    
    c = Cairo.CairoImageSurface(buffer, Cairo.FORMAT_ARGB32, flipxy=false)
    cr = CairoContext(c);        
    
    for i in 1:num_triangles
        q = prms[:,i]
        set_source_rgba(cr,q[7], q[8], q[9], q[10])        
        move_to(cr, q[1],q[2]);
        line_to(cr, q[3],q[4]);
        line_to(cr, q[5],q[6]);
        close_path(cr);
        fill(cr);            
    end        
            
    resultimg = reinterpret(ARGB32, permutedims(c.data, (2, 1)))
    resultchn = Float32.(channelview(RGB.(resultimg)))                
    img .= resultchn
    Cairo.finish(c)
    Cairo.destroy(c)        
end

Now let’s draw each image in this fashion:

@views function withCairo() 
    prms, target, num_images, num_triangles, w, h = prepare() 

    imgs = Array{Float32}(undef, 3, h, w, num_images)    
    for i in 1:num_images
        img = imgs[:,:,:,i]        
        renderCairo(img, prms[:,:,i], num_triangles, w, h)        
    end
    dists = reshape(mean((imgs .- target) .^2, dims=(1,2,3)), num_images)
    return imgs, 1 .- dists
end

Benchmarking this with @btime withCairo();

I see 428.101 ms (3157 allocations: 205.09 MiB).

The cross product method

Now something really cool. Move your mouse inside and outside of the triangle below, you will see a bar chart depicting the magnitude and direction of the 3rd components of the cross products of vectors AB with AP (reds), BC with BP (greens) and CA with CP (blues). Observe that all those bars point same direction only inside the triangle!

Whats great about this is that cross products are just multiplications and subtractions, perfect job to parallelize with a GPU.
What needs to be done is clear, for each of the M images, for each of the N triangles, our operation is to update a pixel color to blend with the current triangles color if that pixel is inside the triangle. We will parallelize this operation with a CUDA kernel, shown below:

function puttri(prms, imgs, tri, ins)    
    idx = (blockIdx().x - 1) * blockDim().x + threadIdx().x  
    idy = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    
    abx = prms[3,tri,ins] - prms[1,tri,ins]
    aby = prms[4,tri,ins] - prms[2,tri,ins]
    apx = idx - prms[1,tri,ins]
    apy = idy - prms[2,tri,ins]
    cr1 = apx * aby - apy * abx

    bcx = prms[5,tri,ins] - prms[3,tri,ins]
    bcy = prms[6,tri,ins] - prms[4,tri,ins]
    bpx = idx - prms[3,tri,ins]
    bpy = idy - prms[4,tri,ins]
    cr2 = bpx * bcy - bpy * bcx

    cax = prms[1,tri,ins] - prms[5,tri,ins]
    cay = prms[2,tri,ins] - prms[6,tri,ins]
    cpx = idx - prms[5,tri,ins]
    cpy = idy - prms[6,tri,ins]
    cr3 = cpx * cay - cpy * cax

    if ((cr1>=0) & (cr2>=0) & (cr3>=0)) | ((cr1<=0) & (cr2<=0) & (cr3<=0))
        oneMinusAlpha = (1.0f0-prms[10,tri,ins])        
        imgs[1,idx,idy,ins] = imgs[1,idx,idy,ins] * oneMinusAlpha + prms[7,tri,ins] * prms[10,tri,ins]
        imgs[2,idx,idy,ins] = imgs[2,idx,idy,ins] * oneMinusAlpha + prms[8,tri,ins] * prms[10,tri,ins]
        imgs[3,idx,idy,ins] = imgs[3,idx,idy,ins] * oneMinusAlpha + prms[9,tri,ins] * prms[10,tri,ins]
    end
    return
end

We will need to pass our parameters and target array to the GPU, and then call the kernel with @cuda. We can create white canvases with CUDA.ones here, so no need to pass it.

function withGpu()
    prms, target, num_images, num_triangles, w, h = prepare()     

    prms = CuArray(prms)    
    imgs = CUDA.ones(3, h, w, num_images)
    target = CuArray(target)
    for tri in 1:num_triangles
        for i in 1:num_images
            @cuda threads=(32,32) blocks=(8,8) puttri(prms, imgs, tri, i)        
        end                                
    end
    gpufitnesses = 1.0f0 .- reshape(mean((imgs .- target) .^ 2, dims=(1,2,3)),num_images)
    return Array(imgs), Array(gpufitnesses)
end

Benchmarking this I see 120.315 ms (38689 allocations: 52.53 MiB)

Awesome, that’s about 4x speedup. Note that I benchmarked with --check-bounds=no, which is a startup option that you pass when launching julia that disables the bounds checking. In the next post, I will talk about the very cool PGPE algorithm used in evojax to steer these images towards a target image. You can see one example of this here.

I am pretty new to Julia, please let me know if you have any comments, suggestions.

An exercise in DataFrames.jl transformation minilanguage

By: Blog by Bogumił Kamiński

Re-posted from: https://bkamins.github.io/julialang/2022/05/06/minilanguage.html

Introduction

Recently I answered an interesting question about transformation of a data
frame. I thought that the problem and solution are instructive enough to
warrant writing a blog post about them.

This post was written under Julia 1.7.0 and DataFrames.jl 1.3.4.

The problem

Assume you are given the following data frame:

julia> using DataFrames

julia> df1 = DataFrame(reshape(1:30, 5, 6), vec(string.(["x", "y"], [1 2 3])))
5×6 DataFrame
 Row │ x1     y1     x2     y2     x3     y3
     │ Int64  Int64  Int64  Int64  Int64  Int64
─────┼──────────────────────────────────────────
   1 │     1      6     11     16     21     26
   2 │     2      7     12     17     22     27
   3 │     3      8     13     18     23     28
   4 │     4      9     14     19     24     29
   5 │     5     10     15     20     25     30

Before we move forward let me comment a bit about this code.
The reshape(1:30, 5, 6) part creates a 5×6 matrix filled with integers
ranging from 1 to 30:

julia> reshape(1:30, 5, 6)
5×6 reshape(::UnitRange{Int64}, 5, 6) with eltype Int64:
 6  11  16  21  26
 7  12  17  22  27
 8  13  18  23  28
 9  14  19  24  29
10  15  20  25  30

Next the string.(["x", "y"], [1 2 3]) part creates a matrix of column names:

julia> string.(["x", "y"], [1 2 3])
2×3 Matrix{String}:
 "x1"  "x2"  "x3"
 "y1"  "y2"  "y3"

I use vec on it since the DataFrame constructor requires column names to
be passed as a vector.

We want to create a new data frame having the following four columns:

x_minimum: storing for each row minimum value stored in
the columns containing "x" in their name;
x_maximum: storing for each row maximum value stored in
the columns containing "x" in their name;
y_minimum: storing for each row minimum value stored in
the columns containing "y" in their name;
y_maximum: storing for each row maximum value stored in
the columns containing "y" in their name.

The question is how to do it in DataFrames.jl. Below I will discuss several
options you can consider.

Using a loop

Here is a simple approach for performing this operation which relies
on knowledge of Base Julia:

julia> df2 = DataFrame()
0×0 DataFrame

julia> for n in ["x", "y"]
           mat = Matrix(df1[:, Regex(n)])
           for fun in [minimum, maximum]
               df2[:, string(n, "_", fun)] = fun.(eachrow(mat))
           end
       end

julia> df2
5×4 DataFrame
 Row │ x_minimum  x_maximum  y_minimum  y_maximum
     │ Int64      Int64      Int64      Int64
─────┼────────────────────────────────────────────
   1 │         1         21          6         26
   2 │         2         22          7         27
   3 │         3         23          8         28
   4 │         4         24          9         29
   5 │         5         25         10         30

What we do in this code can be explained as follows. First we create an empty
target data frame df2. Next we iteratively add columns to it. To be able to
use Base Julia functionality we select from the data frame the columns
respectively having "x" or "y" in their name using a regular expression and
convert the result to a Matrix. Finally we apply the minimum or maximum
function to rows of this matrix with the fun.(eachrow(mat)) expression and
assign the result to a new column in the df2 data frame.

Using broadcasting in transformation minilanguage

Now let us turn to using the DataFrames.jl transformation minilanguage
(if you do not have much experience with it I recommend you to first read
this post before proceeding):

julia> df2 = select(df1, AsTable.([r"x" r"y"]) .=>
                         ByRow.([minimum, maximum]) .=>
                         string.(["x_" "y_"], [minimum, maximum]))
5×4 DataFrame
 Row │ x_minimum  x_maximum  y_minimum  y_maximum
     │ Int64      Int64      Int64      Int64
─────┼────────────────────────────────────────────
   1 │         1         21          6         26
   2 │         2         22          7         27
   3 │         3         23          8         28
   4 │         4         24          9         29
   5 │         5         25         10         30

To understand what is going on in this expression we first need to inspect
what is passed to select as a second argument:

julia> AsTable.([r"x" r"y"]) .=>
       ByRow.([minimum, maximum]) .=>
       string.(["x_" "y_"], [minimum, maximum])
2×2 Matrix{Pair{AsTable}}:
 AsTable(r"x")=>(ByRow{typeof(minimum)}(minimum)=>"x_minimum")  AsTable(r"y")=>(ByRow{typeof(minimum)}(minimum)=>"y_minimum")
 AsTable(r"x")=>(ByRow{typeof(maximum)}(maximum)=>"x_maximum")  AsTable(r"y")=>(ByRow{typeof(maximum)}(maximum)=>"y_maximum")

As you can see Julia broadcasting mechanism magically created four operation
specification expressions. The trick is what since we wanted function names to
change faster I passed them in a vector [minimum, maximum] twice, while I
wanted column names to change slower, so I passed them to broadcasting as one
row matrices with [r"x" r"y"] and ["x_" "y_"] expressions respectively.

The second part is understanding what each of the operation specification
expression means. Let us concentrate on the first one:

AsTable(r"x")=>(ByRow{typeof(minimum)}(minimum)=>"x_minimum")

The decomposition is:

AsTable(r"x") means: select all columns that contain "x" and pass them
to the transformation function as a single positional argument (this is what
AsTable serves for here);
the ByRow(minimum) part means that we want to apply the minimum function
to each row of the data passed to it;
finally "x_minimum" part means that we want to store the result in the column
having this name.

An alternative way to write this transformation would be to replace minimum
and maximum with min and max. The difference is that min and max
take multiple positional arguments. It means that we would need to drop the
AsTable part in the transformation specification like this:

julia> df2 = select(df1, [r"x" r"y"] .=>
                         ByRow.([min, max]) .=>
                         string.(["x_" "y_"], [minimum, maximum]))
5×4 DataFrame
 Row │ x_minimum  x_maximum  y_minimum  y_maximum
     │ Int64      Int64      Int64      Int64
─────┼────────────────────────────────────────────
   1 │         1         21          6         26
   2 │         2         22          7         27
   3 │         3         23          8         28
   4 │         4         24          9         29
   5 │         5         25         10         30

As you can see we get the same result. You might ask why I have not used this
style initially? The reason is that r"x" potentially could have selected
thousands of columns from the source data frame. In Julia, in general, it is
not a good idea to pass very many positional arguments to functions as in some
cases it might put too much strain on the Julia compiler. In such cases
AsTable wrapper is preferred as it guarantees to pass a single argument to
the function(a collection of passed columns).

Using a comprehension in transformation minilanguage

Above I have shown you how to use broadcasting to achieve the desired result.
Let me show below that it is equally easy to use a comprehension to achieve
the same:

julia> df2 = select(df1, [AsTable(Regex(n)) => ByRow(fun) => string(n, "_", fun)
                    for n in ["x", "y"] for fun in [minimum, maximum]])
5×4 DataFrame
 Row │ x_minimum  x_maximum  y_minimum  y_maximum
     │ Int64      Int64      Int64      Int64
─────┼────────────────────────────────────────────
   1 │         1         21          6         26
   2 │         2         22          7         27
   3 │         3         23          8         28
   4 │         4         24          9         29
   5 │         5         25         10         30

The choice between using broadcasting and a comprehension is mostly a personal
preference.

Conclusions

I hope you will find the presented examples useful to better understand how to
write complex transformations in DataFrames.jl.

The codes might look scary to you at a first glance. However, in my experience,
after having some practice with broadcasting or writing comprehensions in
Julia they become natural.

juliabloggers.com

A Julia Language Blog Aggregator