Recently I was invited by Talk Julia to take part in the podcast.
Today it has been released and you can watch it on YouTube.
In the discussion I share my thoughts on DataFrames.jl design principles
and discuss several examples from my upcoming Julia for Data Analysis
book.
One of the things I discuss in the podcast is that in DataFrames.jl development
process we are putting a lot of emphasis on tests so I thought that it is worth
to expand on this topic a bit more in this post.
The codes I use were tested under Julia 1.7.2 and InlineStrings.jl 1.2.2.
Testing and reproducibility of tests
When one develops some software ensuring its proper test coverage is one of the
key practices that help maintaining good quality of code.
My experience with DataFrames.jl is that taking care of proper test coverage
serves three important goals:
Making sure the functionality we provide follows the contract specified in
its documentation. This is particularly important when testing corner case
situations, since they happen rarely in practice (so there is lower
probability that users spot problems and report them) and also allow us to
make sure our design is logically consistent.
Making sure that as we add new functionalities to the package we do not
accidentally break something.
Making sure that upgrading dependencies of DataFrames.jl to newer versions
does not break something (and the biggest such dependency is Base Julia).
A crucial part of writing tests is ensuring their reproducibility. What I mean
by this is that when you run your test suite twice on you should get the same
results. This intends to avoid situations in which you run your tests and get a
bug report. Then you run the tests again and the bug is not present. Such
situation is unwanted, as it is later hard to locate the root cause of the
reported bug.
Randomized tests
When one implements complex algorithms it is hard to cover all possible
hard testing scenarios by writing them down by hand. In such situations,
one of the possible testing methods is to use randomized tests.
An example, old and already resolved, problem that potentially could have been
caught by randomized tests is an issue related to pasting Unicode
characters in Julia REPL. Let me here present a minimal working example of a
code having a similar problem.
Assume I want to write a code that strips last character from a string. Here
is an attempt to implement it:
julia> mychop(s::AbstractString) = isempty(s) ? s : s[1:end-1]
mychop (generic function with 1 method)
Let us write a simple test set for this function:
julia> using Test
julia> @testset "basic test" begin
@test mychop("") == ""
@test mychop("a") == ""
@test mychop("abc") == "ab"
end
Test Summary: | Pass Total
basic test | 3 3
Test.DefaultTestSet("basic test", Any[], 3, false, false)
All looks good so far. However, let us run some more advanced testing
of mychop using randomized tests:
julia> using Random
julia> @testset "advanced tests" begin
Random.seed!(1234)
for _ in 1:40
len = rand(1:10)
input = rand(UInt8, len) |> String
output = join(collect(input)[1:end-1])
@test output == mychop(input)
end
end
advanced tests: Error During Test at REPL[83]:7
Test threw exception
Expression: output == mychop(input)
StringIndexError: invalid index [9], valid nearby indices [8]=>'˄', [10]=>'�'
Test Summary: | Pass Error Total
advanced tests | 39 1 40
ERROR: Some tests did not pass: 39 passed, 0 failed, 1 errored, 0 broken.
What I do in this test is generateing input string using randoom bits and then
use a slow method to get a desired output string by going through individual
characters of the original string. I ensure reproducibility of this test on a
given version of Julia by setting the seed of random number generator using Random.seed!(1234).
From the bug report we can see that something is not right with indexing. The
problem occurs if the second character from the end of the string is not ASCII.
Let us check it:
julia> mychop("∀a")
ERROR: StringIndexError: invalid index [3], valid nearby indices [1]=>'∀', [4]=>'a'
The point of randomized test is that guessing the test scenario second character from the end of the string is not ASCII
is not easy.
Let us fix mychop to resolve this problem:
julia> mychop(s::AbstractString) = isempty(s) ? s : s[1:prevind(s, end)]
mychop (generic function with 1 method)
julia> @testset "basic test" begin
@test mychop("") == ""
@test mychop("a") == ""
@test mychop("abc") == "ab"
end
Test Summary: | Pass Total
basic test | 3 3
Test.DefaultTestSet("basic test", Any[], 3, false, false)
julia> @testset "advanced tests" begin
Random.seed!(1234)
for _ in 1:32768 # pick some round larger number to be sure all works well
len = rand(1:10)
input = rand(UInt8, len) |> String
output = join(collect(input)[1:end-1])
@test output == mychop(input)
end
end
Test Summary: | Pass Total
advanced tests | 32768 32768
Test.DefaultTestSet("advanced tests", Any[], 32768, false, false)
Are we done now? Not yet. We should check our function with some other AbstractString than String. Let me use InlineStrings.jl:
julia> using InlineStrings
julia> @testset "InlineStrings.jl test" begin
@test "ab" == @inferred mychop(InlineString("abc"))
end
InlineStrings.jl test: Error During Test at REPL[110]:2
Test threw exception
Expression: "ab" == #= REPL[110]:2 =# @inferred(mychop(InlineString("abc")))
return type SubString{String3} does not match inferred return type Union{SubString{String3}, String3}
Test Summary: | Error Total
InlineStrings.jl test | 1 1
ERROR: Some tests did not pass: 0 passed, 0 failed, 1 errored, 0 broken.
As you can see there is still some more work to be done, as our mychop
function is not type stable, which is checked by the @inferred macro. The
problem is that indexing into String3 that happens if the passed string is not
empty produces a SubString. Let us fix it by making sure that in every branch
of code we apply the same operation to our source string (in the code, like in
the codes above, we use the fact that AbstractString is guaranteed to use 1-based indexing).
julia> mychop(s::AbstractString) = s[1:prevind(s, max(1, end))]
mychop (generic function with 1 method)
julia> @testset "advanced tests" begin
Random.seed!(1234)
for _ in 1:32768 # pick some round larger number to be sure all works well
len = rand(0:10) # make sure to cover 0-length strings
input = rand(UInt8, len) |> String
output = join(collect(input)[1:end-1])
@test output == @inferred mychop(input)
@test output == @inferred mychop(InlineString(input))
end
end
Test Summary: | Pass Total
advanced tests | 65536 65536
Test.DefaultTestSet("advanced tests", Any[], 65536, false, false)
Now all works as expected and we are done.
Conclusions
I hope you found my thoughts on writing tests useful. As usual, I have some
additional comments regarding the presented codes:
I used Random.seed! explicitly in my tests. However, the @testset macro,
before the execution of its body, makes a call to Random.seed!(seed) where seed is the current seed of the global RNG. Therefore, if you design a
larger test suite you do not have to set the seed in every @testset. It is
enough to set it once per all tests you run.
As I have already commented, the presented codes are reproducible on a given
version of Julia. If you want them to be reproducible across different
versions of Julia I recommend you to use for example the StableRNGs.jl package as a source of randomness in your tests.
I’ve long been interested in exactly how R works – not quite enough for me to learn all the
internals, but I was surprised that I could not find a clear guide towards exactly how vectorization works at the deepest level.
Let’s say we want to add two vectors which we’ve defined as x and y
x <- c(2, 4, 6)
y <- c(1, 3, 2)
One way to do this (the verbose, elementwise way) would be to add each pair of elements
c(x[1] + y[1], x[2] + y[2], x[3] + y[3])
## [1] 3 7 8
but if you are familiar with not repeating yourself, you might write this in a loop. Best practice
involves pre-filling the result to the correct size
res <- c(NA, NA, NA)
for (i in 1:3) {
res[i] <- x[i] + y[i]
}
res
## [1] 3 7 8
There, we wrote a for loop.
Now, if you’ve ever seen R tutorials or discussions, you’ve probably had it drilled into you that
this is “bad” and should be replaced with something else. One of those options is an apply family
function
plus <- function(i) {
x[i] + y[i]
}
sapply(1:3, plus)
## [1] 3 7 8
or something ‘tidier’
purrr::map_dbl(1:3, plus)
## [1] 3 7 8
(yes, yes, these functions are ‘bad’ because they don’t take x or y as arguments) but this stil
performs the operation elementwise. If you’ve seen even more R tutorais/discussions you’ve
probably been seen that vectorization is very handy – The R function + knows what to do with objects that
aren’t just a single value, and does what you might expect
x + y
## [1] 3 7 8
Now, if you’ve really read a lot about R, you’ll know that ‘under the hood’ a for-loop
is involved in every one of these, but it’s “lower down”, “at the C level”. Jenny Bryan makes the
point that “Of course someone has to write loops / It doesn’t have to be you” and for this reason,
vectorization in R is of great benefit.
So, there is a loop, but where exactly does that happen?
At some point, the computer needs to add the elements of x to the elements of y, and the simplest versions of this
happens one element at a time, in a loop. There’s a big sidetrack here about SIMD
which I’ll try to avoid, but I will mention that the Microsoft fork of R (artist, formerly known as Revolution R) running on Intel chips can do SIMD in MKL.
So, let’s start at the operator.
`+`
## function (e1, e2) .Primitive("+")
Digging into primitives is a little tricky, but {pryr} can help
pryr::show_c_source(.Primitive("+"))
+ is implemented by do_arith with op = PLUSOP
We can browse a copy of the source for do_arith (in arithmetic.c) here where we see some
logic paths for scalars and vectors. Let’s assume we’re working with our example which has length(x) == length(y) > 1. With two non-scalar arguments
Depending on the class of the arguments, we need to call different functions, but for the sake of our example let’s say we have non-integer real numbers so we fork to real_binary. This takes a code argument for which type of operation we’re performing, and in our case it’s PLUSOP (noted above). There’s a case branch for this in which case, provided the arguments are of the same length (n1 == n2) we call
R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] + dy[i];);
That’s starting to look a lot like a loop – there’s an iterator i and we’re going to call another function.
This jumps us over to a different file where we see LOOP_WITH_INTERRUPT_CHECK definitely performs some sort of loop. This takes the body above and the argument LOOP_ITERATE_CORE which is finally the actual loop!
#define R_ITERATE_CORE(n, i, loop_body) do { \
for (; i < n; ++i) { loop_body } \
} while (0)
so, that’s where the actual loop in a vectorized R call happens! ALL that sits behind the innocent-looking +.
That was thoroughly satisfying, but I did originally have in mind comparing R to another language – one where loops aren’t frowned upon because of performance, but rather encouraged… How do Julia loops differ?
Julia is not a vectorized language per se, but it has a neat ability to “vectorize” any operation, though in Julia syntax it’s “broadcasting”.
Simple addition can combine scalar values
3+4
## 7
Julia actually has scalar values (in R, even a single value is just a vector of length 1) so a single value could be
typeof(3)
## Int64
whereas several values need to be an Array, even if it only has 1 dimension
Vector{Int64}([1, 2, 3])
## 3-element Array{Int64,1}:
## 1
## 2
## 3
Trying to add two Arrays does work
[1, 2, 3] + [4, 5, 6]
## 3-element Array{Int64,1}:
## 5
## 7
## 9
but only because a specific method has been written for this case, i.e.
methods(+, (Array, Array))
## # 1 method for generic function "+":
## [1] +(A::Array, Bs::Array...) in Base at arraymath.jl:43
One thing I particularly like is that we can see exactly which method was called using the @which macro
@which [1, 2, 3, 4] + [1, 2, 3, 4]
+(A::Array, Bs::Array...) in Base at arraymath.jl:43
something that I really wish was easier to do in R. The @edit macro even jumps us right into the actual code for this dispatched call.
This ‘add vectors’ problem can be solved through broadcasting, which performs an operation elementwise
[1, 2, 3] .+ [4, 5, 6]
## 3-element Array{Int64,1}:
## 5
## 7
## 9
The fun fact about this I recently learned was that broadcasting works on any operation, even if that’s the pipe itself
This blog post serves a couple purposes. The first is to introduce Julia users to GraphBLAS, and the features and performance of the SuiteSparseGraphBLAS.jl package. The second is an update on new features in versions 0.6 and 0.7 as well as a roadmap for the near future. That's a lot for one blogpost so let's dive right in!
What is SuiteSparseGraphBLAS.jl
SuiteSparseGraphBLAS.jl is a sparse linear algebra library with some special features oriented to graph algorithms. It's an interface to the SuiteSparse implementation of the GraphBLAS standard, by Tim Davis. SuiteSparseGraphBLAS.jl serves two purposes: (1) it can replace SparseArrays.jl with faster parallel methods when using the plus-times semiring conventional linear algebra, and (2) it serves as a Julia interface to the full set of GraphBLAS features. GraphBLAS is a powerful tool for writing high-performance graph algorithms at a high level of abstraction, using linear algebraic operations on sparse adjacency matries with different semirings.
The headline feature is taking the matrix multiplication function, and turning it into a higher order function. Instead of the normal arithmetic operations + and * a user could substitute other binary operators, like max or &&.
When we multiply two matrices A and B as normal we use the following index expression:
This expression uses the arithmetic semiring \((+, \times)\). A semiring is an abstract algebra structure consisting of two binary operators that interact in a particular way.
There are plenty of other interesting semirings to choose from, and no reason to limit ourselves to the arithmetic one! We could, for instance, use the max-plus semiring \((\max, +)\) (one of the tropical semirings). Then our index expression above looks like:
This semiring has many applications to dynamic programming and graphs, but to wrap up this short intro, GraphBLAS lets you perform this operation, particularly on sparse matrices:
where \(\oplus\) is some binary operator (a monoid in particular) which takes the place of the \(\sum\) reduction, and \(\otimes\) is some binary operator which takes the place of \(\times\) in typical arithmetic.
Let's see some code!!
What does a simple operation look like? The equivalence between a matrix-vector multiplication and breadth-first search is a key component of linear algebraic graph algorithms.
Adjacency matrix of a directed graph and a single iteration of BFS
julia> using SuiteSparseGraphBLASjulia> A = GBMatrix([1,1,2,2,3,4,4,5,6,7,7,7], [2,4,5,7,6,1,3,6,3,3,4,5], [1:12...])
7x7 GraphBLAS int64_t matrix, bitmap by row
12 entries, memory: 832 bytes (1,2) 1
(1,4) 2
(2,5) 3
(2,7) 4
(3,6) 5
(4,1) 6
(4,3) 7
(5,6) 8
(6,3) 9
(7,3) 10
(7,4) 11
(7,5) 12julia> v = GBVector([4], [10]; nrows = 7)
7x1 GraphBLAS int64_t matrix, bitmap by col
1 entry, memory: 272 bytes
iso value: 10 (4,1) 10julia> A * v
7x1 GraphBLAS int64_t matrix, bitmap by col
2 entries, memory: 328 bytes (1,1) 20
(7,1) 110
This looks exactly like a matrix-vector multiplication with any other array type in Julia. With a different semiring, the result of A*v can return the parent node id of nodes 1 and 7, so that a BFS tree can be constructed.
Where possible, SuiteSparseGraphBLAS.jl will behave exactly as any other array type in Julia. But there are some places where SuiteSparseGraphBLAS.jl has significant extra functionality. I'll illustrate that here by finding the number of triangles in the graph above.
For those that haven't seen triangle counting before, a triangle in a graph is a set of three vertices whose edges form a triangle. We'll pretend like the graph is undirected, and then you can see two triangles: \((2, 5, 7)\) and \((3, 4, 7)\).
julia> using SuiteSparseGraphBLAS: pairjulia> function cohen(A)
U = triu(A)
L = tril(A)
return reduce(+, *(L, U, (+, pair); mask=A)) ÷ 2
end
cohen (generic function with 1 method)julia> function sandia(A)
L = tril(A)
return reduce(+, *(L, L, (+, pair); mask=L))
end
sandia (generic function with 1 method)julia> M = eadd(A, A', +) #Make undirected/symmetric
7x7 GraphBLAS int64_t matrix, bitmap by row
20 entries, memory: 832 bytes (1,2) 1
(1,4) 8
(2,1) 1
(2,5) 3
(2,7) 4
(3,4) 7
(3,6) 14
(3,7) 10
(4,1) 8
(4,3) 7
(4,7) 11
(5,2) 3
(5,6) 8
(5,7) 12
(6,3) 14
(6,5) 8
(7,2) 4
(7,3) 10
(7,4) 11
(7,5) 12julia> cohen(M)
2julia> sandia(M)
2
There are a couple unique features of SuiteSparseGraphBLAS.jl used in the two methods, cohen and sandia above.
The first is the pair function. This function returns 1 whenever both arguments x and y are explicit stored values, otherwise it returns nothing (an implicit zero). To illustrate:
julia> using SuiteSparseGraphBLASjulia> u = GBVector([2,4,5], [4,5,6])
5x1 GraphBLAS int64_t matrix, bitmap by col
3 entries, memory: 328 bytes (2,1) 4
(4,1) 5
(5,1) 6julia> v = GBVector([1,3,5], [1,2,3])
5x1 GraphBLAS int64_t matrix, bitmap by col
3 entries, memory: 328 bytes (1,1) 1
(3,1) 2
(5,1) 3julia> SuiteSparseGraphBLAS.pair.(u, v)
5x1 GraphBLAS int64_t matrix, bitmap by col
1 entry, memory: 272 bytes
iso value: 1 (5,1) 1
This function is primarily a performance enhancement. In an algorithm like triangle counting we don't care about the weight of a particular edge, just that it exists. pair lets us avoid a costly multiplication by just checking that x is a stored value in uandy is a stored value in v.
The second feature is the semirings discussed in What is SuiteSparseGraphBLAS.jl above. The tuple (+, pair) in *(L, U, (+, pair); mask=A) indicates the the * function is using the +–pair semiring:
Finally, *(L, U, (+, pair); mask=A) uses A as a mask. The mask prevents values from being placed in the result where the mask is false (or true if complemented). This is a powerful algorithmic and performance tool.
Summary of Primary Functions
Most of the SuiteSparse:GraphBLAS interface is accessible through the expected Julia functions, alongside some new functions like emul or subassign!.
The complete documentation of supported operations can be found in Operations.
GraphBLAS
Operation
Julia
mxm, mxv, vxm
\(\bf C \langle M \rangle = C \odot AB\)
mul! or *
eWiseMult
\(\bf C \langle M \rangle = C \odot (A \otimes B)\)
emul[!] or . broadcasting
eWiseAdd
\(\bf C \langle M \rangle = C \odot (A \oplus B)\)
eadd[!]
extract
\(\bf C \langle M \rangle = C \odot A(I,J)\)
extract[!], getindex
subassign
\(\bf C (I,J) \langle M \rangle = C(I,J) \odot A\)
subassign[!] or setindex!
assign
\(\bf C \langle M \rangle (I,J) = C(I,J) \odot A\)
assign[!]
apply
\({\bf C \langle M \rangle = C \odot} f{\bf (A)}\)
apply[!], map[!] or . broadcasting
\({\bf C \langle M \rangle = C \odot} f({\bf A},y)\)
\({\bf C \langle M \rangle = C \odot} f(x,{\bf A})\)
select
\({\bf C \langle M \rangle = C \odot} f({\bf A},k)\)
select[!]
reduce
\({\bf w \langle m \rangle = w \odot} [{\oplus}_j {\bf A}(:,j)]\)
reduce[!]
\(s = s \odot [{\oplus}_{ij} {\bf A}(i,j)]\)
transpose
\(\bf C \langle M \rangle = C \odot A^{\sf T}\)
gbtranspose[!], lazy: transpose, '
kronecker
\(\bf C \langle M \rangle = C \odot \text{kron}(A, B)\)
kron[!]
where \(\bf M\) is a GBArray mask, \(\odot\) is a binary operator for accumulating into \(\bf C\), and \(\otimes\) and \(\oplus\) are a binary operation and commutative monoid respectively. \(f\) is either a unary or binary operator.
Show Me the Numbers!
SuiteSparseGraphBLAS.jl has loads of extensions to normal sparse linear algebra, but it's also fast and multithreaded. Let's look at some numbers!
As always, benchmark things yourself. Most operations will be faster in SuiteSparseGraphBLAS.jl, particularly when the matrices are large enough that multithreading kicks in.
However, maintaining good performance can be tricky in any numerical package, and there's plenty of ways to accidentally reduce performance. For instance, below you'll notice that when A is stored in RowMajor format it can be quite a bit faster than operations where A is stored in ColMajor format. This isn't always the case, some operations favor column orientation.
The plots below show the runtime of several operations in SuiteSparseGraphBLAS.jl normalized to the runtime of SparseArrays.SparseMatrixCSC (orange bars) performing the same operations. SuiteSparseGraphBLAS.jl is shown in both column (blue bars) and row (red bars) orientation where necessary, and with 1, 2, and 16 threads.
Matrix multiplication results are discussed briefly after the 3rd plot, all matrices are drawn from the SuiteSparse Matrix Collection unless randomly generated.
When the dense matrix is low-dimensional and only a single thread is used, Julia compares very favorably with, and often beats, SuiteSparseGraphBLAS.jl. However, when 2 threads are available and the sparse matrix is stored in row-major orientation SuiteSparseGraphBLAS.jl begins to pull ahead significantly. Once a full 16 threads are used SuiteSparseGraphBLAS.jl on multiplication of row-major matrices is between 8 and 31 times faster.
SuiteSparse:GraphBLAS uses well over a dozen subalgorithms for matrix multiplication internally to achieve this performance. In particular when A is a sparse row-oriented matrix, and B is a dense column oriented matrix SuiteSparse:GraphBLAS will switch to a highly optimized dot-product algorithm, which is often much faster than the saxpy based algorithm used when A is column oriented. This is noticeable in cases above where the row-major runtimes can be > 2x faster than the column-major ones.
Transpose
SparseArrays.jl maintains performance parity for most problems when SuiteSparseGraphBLAS.jl is restricted to a single thread, but can be as much as 10x slower when 16 threads are used. This test takes advantage of some special SuiteSparseGraphBLAS.jl features like iso-valued matrices, but we restrict some significant optimizations like quick transposition (reinterpreting a by-column matrix as by-row) to remain fair.
Submatrix Assignment
Note the log scale
Benchmark Group
Size (Density) of C
Size (Density) of A
A
\(10,000^{2} \quad (0.001)\)
\(2000^{2} \quad (0.1)\)
B
\(1,000,000^{2} \quad (0.0005)\)
\(5000^{2} \quad (0.005)\)
C
\(25,000,000^{2} \quad (1 \times 10^{-7})\)
\(5000^{2} \quad (0.002)\)
D
\(50,000,000^{2} \quad (1 \times 10^{-7})\)
\(100,000^{2} \quad (0.001)\)
E
\(50,000,000^{2} \quad (1 \times 10^{-7})\)
\(1000^{2} \quad (1.0)\)
The expression C[I,J] = A assigns a matrix A into a submatrix of C, a difficult function to optimize. For a large sparse random matrix C (25 million by 25 million with 62 million entries), and with I and J selected at random, and A of size 5000-by-5000 with 50,000 entries, Julia takes 0.87 seconds using SparseArrays.jl, while using SuiteSparseGraphBLAS.jl with a single thread the time is just 0.04 seconds.
The runtime for a GBMatrix by-row or by-column is equivalent, so only by-column is shown here.
A New Version of SuiteSparseGraphBLAS.jl
Versions 0.6 and 0.7 of SuiteSparseGraphBLAS.jl brought with them several new features and many under the hood changes. The changes under the hood will support faster development in the future and otherwise cleaned up the codebase. The type system was completely overhauled to enable new matrix types with special extensions, the low level wrapper was scrubbed of any human-written contamination, and the operator/user-defined-function system was rewritten.
New Features
User Defined Fill Value
The default GBMatrix returns nothing when indexing into an "implicit-zero" or non-stored value. This better matches the GraphBLAS method of attaching the implicit values to operators like max and *. It also is more natural for graphs, where there is a semantic difference between a stored zero (an edge with weight zero), and an implicit zero (no edge).
Users may now directly set the fill value best suited to their application. This can be done in 3 ways:
On construction: a new keyword argument for constructors, fill.
setfill!: a new mutating function which can be used to change the fill value. This function may only change the fill value to another value within the same domain as the original fill value to maintain type stability.
setfill: a new non-mutating function which returns a new GBMatrix which aliases the original, but with a new fill value. This new value may be of any type.
These changes allow the GBMatrix to be seamlessly used for general scientific applications that expect implicit zeros to be zero, while still adhering to the original design intended for graph algorithms.
mmread Function
Previously reading Matrix Market files was supported only by first reading into a SparseArrays.SparseMatrixCSC before converting, which was slow and memory intensive. We now have a native SuiteSparseGraphBLAS.mmread function. It remains unexported to avoid clashes with other packages. In the future it will likely be superseded (but not replaced) by some function which can automatically detect and read from a number of formats.
Serialization
Using the Serialization Julia standard library users can now easily serialize and deserialize a GraphBLAS data structures. The serialized array is compressed using LZ4 making both read and write operations much faster than operating on a Matrix Market file.
StorageOrders
A new dependency on StorageOrders.jl makes it easier to parametrize new AbstractGBArray types, for instance by restricting the orientation to StorageOrders.RowMajor() or StorageOrders.ColMajor(). The experimental OrientedGBMatrix and GBMatrixC/GBMatrixR subtypes take advantage of this.
Users may now also call gbset(A, :format, RowMajor()), instead of gbset(A, :format, :byrow) avoiding the use of magic symbols. This interface will continue to improve in the future.
apply and Deprecation of Scalar Argument map
The new function apply[!] is now the direct wrapper of the GrB_apply function. map still functions as expected, but no longer accepts a scalar argument. Previously, when map was the direct wrapper of GrB_apply, map(+, A, 3) was legal and equivalent to A .+ 3. This caused many (mostly obscure) ambiguities between map(op::Function, A::GBArray, x), map(op::Function, x, A::GBArray) and several external methods.
apply now performs this function, although the recommended method remains dot-broadcasting.
wait Function Fully Implemented
To fully support calling SuiteSparseGraphBLAS.jl from multiple user threads the wait function has been fully implemented, which allows users to use a matrix from multiple threads simultaneously.
Experimental Functionality
as Function
The new as functions grew out of an internal need to safely and quickly view a GBMatrix as a DenseMatrix/SparseMatrixCSC or vice-versa.
julia> A = GBMatrix([[1, 2] [3,4]])
2x2 GraphBLAS int64_t matrix, full by col
4 entries, memory: 288 bytes (1,1) 1
(2,1) 2
(1,2) 3
(2,2) 4julia> SuiteSparseGraphBLAS.as(Matrix, A) do mat
display(mat)
mat .+= 1
sum(mat)
end
2×2 Matrix{Int64}:
1 3
2 4
14julia> A
2x2 GraphBLAS int64_t matrix, full by col
4 entries, memory: 288 bytes (1,1) 2
(2,1) 3
(1,2) 4
(2,2) 5
Note that this functionality is currently somewhat dangerous. If mat escapes the scope of as in some way, for instance by returning the Transpose of mat, the underlying memory may be freed by SuiteSparseGraphBLAS.jl. If the user attempts to return mat directly the as function will gracefully copy the matrix rather than return an array that may be invalidated in the future.
Automatic Differentiation
Automatic Differentiation support continues to improve, with additional constructor rules, and rules for the second and first families of semirings. The biggest remaining missing pieces are the tropical semirings, user-defined functions and more rigorous testing of user use-cases.
Roadmap
There's lots of upcoming features, extensions, and JuliaLab Compiler Trickery™ on the horizon but here's a few important updates to look out for:
SuiteSparse Support
Supporting the SuiteSparse solvers is 1st on the list, and support will be released sometime in April 2022. This is a relatively simple addition, but involves quite a bit of duplicated code from SuiteSparse.jl, which assumes the use of SparseMatrixCSC.
GBPirate.jl Full Release
GBPirate.jl allows users of SparseArrays.jl to pirate certain methods of that package to use the more performant ones found in SuiteSparseGraphBLAS.jl. Because this involves copying on output it is not always faster, and the remaining development involves finding a heuristic for when this is beneficial.
User Defined Types
User defined types are not currently part of the public interface of SuiteSparseGraphBLAS.jl while the interface is improved. This feature will enable users to use any statically sized isbits type, like dual-numbers, colors, as elements of a GBMatrix. Release is planned for late April 2022.
GBGraphs.jl
A new Graphs.jl backend built on SuiteSparseGraphBLAS.jl will allow users access to the full Graphs.jl ecosystem, while letting them develop fast algorithms in the language of linear algebra!
Conclusion
If you've got any need for sparse matrix operations give SuiteSparseGraphBLAS.jl a try, and feel free to open an issue or contact me directly if you run into problems or have a feature request!
– Will Kimmerer
Acknowledgements
Thanks to Tim Davis for his help on the design of SuiteSparseGraphBLAS.jl and for the underlying C library SuiteSparse:GraphBLAS.
