Author Archives: Andrew Collier

#MonthOfJulia Day 25: Interfacing with Other Languages

Julia-Logo-Other-Languages

Julia has native support for calling C and FORTRAN functions. There are also add on packages which provide interfaces to C++, R and Python. We’ll have a brief look at the support for C and R here. Further details on these and the other supported languages can be found on github.

Why would you want to call other languages from within Julia? Here are a couple of reasons:

  • to access functionality which is not implemented in Julia;
  • to exploit some efficiency associated with another language.

The second reason should apply relatively seldom because, as we saw some time ago, Julia provides performance which rivals native C or FORTRAN code.

C

C functions are called via ccall(), where the name of the C function and the library it lives in are passed as a tuple in the first argument, followed by the return type of the function and the types of the function arguments, and finally the arguments themselves. It’s a bit klunky, but it works!

julia> ccall((:sqrt, "libm"), Float64, (Float64,), 64.0)
8.0

It makes sense to wrap a call like that in a native Julia function.

julia> csqrt(x) = ccall((:sqrt, "libm"), Float64, (Float64,), x);
julia> csqrt(64.0)
8.0

This function will not be vectorised by default (just try call csqrt() on a vector!), but it’s a simple matter to produce a vectorised version using the @vectorize_1arg macro.

julia> @vectorize_1arg Real csqrt;
julia> methods(csqrt)
# 4 methods for generic function "csqrt":
csqrt{T<:Real}(::AbstractArray{T<:Real,1}) at operators.jl:359
csqrt{T<:Real}(::AbstractArray{T<:Real,2}) at operators.jl:360
csqrt{T<:Real}(::AbstractArray{T<:Real,N}) at operators.jl:362
csqrt(x) at none:6

Note that a few extra specialised methods have been introduced and now calling csqrt() on a vector works perfectly.

julia> csqrt([1, 4, 9, 16])
4-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0

R

I’ll freely admit that I don’t dabble in C too often these days. R, on the other hand, is a daily workhorse. So being able to import R functionality into Julia is very appealing. The first thing that we need to do is load up a few packages, the most important of which is RCall. There’s great documentation for the package here.

julia> using RCall
julia> using DataArrays, DataFrames

We immediately have access to R’s builtin data sets and we can display them using rprint().

julia> rprint(:HairEyeColor)
, , Sex = Male

       Eye
Hair    Brown Blue Hazel Green
  Black    32   11    10     3
  Brown    53   50    25    15
  Red      10   10     7     7
  Blond     3   30     5     8

, , Sex = Female

       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

We can also copy those data across from R to Julia.

julia> airquality = DataFrame(:airquality);
julia> head(airquality)
6x6 DataFrame
| Row | Ozone | Solar.R | Wind | Temp | Month | Day |
|-----|-------|---------|------|------|-------|-----|
| 1   | 41    | 190     | 7.4  | 67   | 5     | 1   |
| 2   | 36    | 118     | 8.0  | 72   | 5     | 2   |
| 3   | 12    | 149     | 12.6 | 74   | 5     | 3   |
| 4   | 18    | 313     | 11.5 | 62   | 5     | 4   |
| 5   | NA    | NA      | 14.3 | 56   | 5     | 5   |
| 6   | 28    | NA      | 14.9 | 66   | 5     | 6   |

rcopy() provides a high-level interface to function calls in R.

julia> rcopy("runif(3)")
3-element Array{Float64,1}:
 0.752226
 0.683104
 0.290194

However, for some complex objects there is no simple way to translate between R and Julia, and in these cases rcopy() fails. We can see in the case below that the object of class lm returned by lm() does not diffuse intact across the R-Julia membrane.

julia> "fit <- lm(bwt ~ ., data = MASS::birthwt)" |> rcopy
ERROR: `rcopy` has no method matching rcopy(::LangSxp)
 in rcopy at no file
 in map_to! at abstractarray.jl:1311
 in map_to! at abstractarray.jl:1320
 in map at abstractarray.jl:1331
 in rcopy at /home/colliera/.julia/v0.3/RCall/src/sexp.jl:131
 in rcopy at /home/colliera/.julia/v0.3/RCall/src/iface.jl:35
 in |> at operators.jl:178

But the call to lm() was successful and we can still look at the results.

julia> rprint(:fit)

Call:
lm(formula = bwt ~ ., data = MASS::birthwt)

Coefficients:
(Intercept)          low          age          lwt         race  
    3612.51     -1131.22        -6.25         1.05      -100.90  
      smoke          ptl           ht           ui          ftv  
    -174.12        81.34      -181.95      -336.78        -7.58 

You can use R to generate plots with either the base functionality or that provided by libraries like ggplot2 or lattice.

julia> reval("plot(1:10)");                # Will pop up a graphics window...
julia> reval("library(ggplot2)");
julia> rprint("ggplot(MASS::birthwt, aes(x = age, y = bwt)) + geom_point() + theme_classic()")
julia> reval("dev.off()")                  # ... and close the window.

Watch the videos below for some other perspectives on multi-language programming with Julia. Also check out the complete code for today (including examples with C++, FORTRAN and Python) on github.



The post #MonthOfJulia Day 25: Interfacing with Other Languages appeared first on Exegetic Analytics.

#MonthOfJulia Day 24: Graphs

Julia-Logo-Graphs

If you’re not too familiar Graph Theory, then it might be an idea to take a moment to get the basics. Graphs are an extremely versatile data structure for storing data consisting of linked entities. I’m going to look at two packages for managing graphs in Julia: LightGraphs and Graphs.

LightGraphs

As usual, the first step is to load the package.

julia> using LightGraphs

LightGraphs has methods which generate a selection of standard graphs like StarGraph(), WheelGraph() and FruchtGraph(). There are also functions for random graphs, for example, erdos_renyi() and watts_strogatz(). We’ll start off by creating two small graphs. One will have 10 nodes connected by 20 random edges. The other will be a directed star graph consisting of four nodes, the central node being connected to every other node.

julia> g1 = Graph(10, 20)
{10, 20} undirected graph
julia> g2 = StarDiGraph(4)
{4, 3} directed graph
julia> edges(g2)
Set{Pair{Int64,Int64}}({edge 1 - 2,edge 1 - 4,edge 1 - 3})

It’s simple to find the degree and neighbours of a given node.

julia> degree(g1, 4)                       # How many neighbours for vertex 4?
6
julia> neighbors(g1, 4)                    # Find neighbours of vertex 4
6-element Array{Int64,1}:
 1
 3
 6
 2
 9
 7

There’s a straightforward means to add and remove edges from the graph.

julia> add_edge!(g1, 4, 8)                 # Add edge between vertices 4 and 8
edge 4 - 8
julia> rem_edge!(g1, 4, 6)                 # Remove edge between vertices 4 and 6
edge 6 - 4

The package has functionality for performing high level tests on the graph (checking, for instance, whether it is cyclic or connected). There’s also support for path based algorithms, but we’ll dig into those when we look at the Graphs package.

Graphs

Before we get started with the Graphs package you might want to restart your Julia session to purge all of that LightGraphs goodness. Take a moment to browse the Graphs.jl documentation, which is very comprehensive.

julia> using Graphs

As with LightGraphs, there are numerous options for generating standard graphs.

julia> g1a = simple_frucht_graph()
Undirected Graph (20 vertices, 18 edges)
julia> g1b = simple_star_graph(8)
Directed Graph (8 vertices, 7 edges)
julia> g1c = simple_wheel_graph(8)
Directed Graph (8 vertices, 14 edges)

Graphs uses the GraphViz library to generate plots.

julia> plot(g1a)

sample-graph

Of course, a graph can also be constructed manually.

julia> g2 = simple_graph(4)
Directed Graph (4 vertices, 0 edges)
julia> add_edge!(g2, 1, 2)
edge [1]: 1 -- 2
julia> add_edge!(g2, 1, 3)
edge [2]: 1 -- 3
julia> add_edge!(g2, 2, 3)
edge [3]: 2 -- 3

Individual vertices (a vertex is the same as a node) can be interrogated. Since we are considering a directed graph we look separately at the edges exiting and entering a node.

julia> num_vertices(g2)
4
julia> vertices(g2)
1:4
julia> out_degree(1, g2)
2
julia> out_edges(1, g2)
2-element Array{Edge{Int64},1}:
 edge [1]: 1 -- 2
 edge [2]: 1 -- 3
julia> in_degree(2, g2)
1
julia> in_edges(2, g2)
1-element Array{Edge{Int64},1}:
 edge [1]: 1 -- 2

Vertices can be created with labels and attributes.

julia> V1 = ExVertex(1, "V1");
julia> V1.attributes["size"] = 5.0
5.0
julia> V2 = ExVertex(2, "V2");
julia> V2.attributes["size"] = 3.0
3.0
julia> V3 = ExVertex(3, "V3")
vertex [3] "V3"

Those vertices can then be used to define edges, which in turn can have labels and attributes.

julia> E1 = ExEdge(1, V1, V2)
edge [1]: vertex [1] "V1" -- vertex [2] "V2"
julia> E1.attributes["distance"] = 50
50
julia> E1.attributes["color"] = "green"
"green"

Finally the collection of vertices and edges can be gathered into a graph.

julia> g3 = edgelist([V1, V2], [E1], is_directed = true)
Directed Graph (2 vertices, 1 edges)

It’s possible to systematically visit all connected vertices in a graph, applying an operation at every vertex. traverse_graph() performs the graph traversal using either a depth first or breadth first algorithm. In the sample code below the operation applied at each vertex is LogGraphVisitor(), which is a simple logger.

julia> traverse_graph(g1c, DepthFirst(), 1, LogGraphVisitor(STDOUT))
discover vertex: 1
examine neighbor: 1 -> 2 (vertexcolor = 0, edgecolor= 0)
discover vertex: 2
open vertex: 2
examine neighbor: 2 -> 3 (vertexcolor = 0, edgecolor= 0)
discover vertex: 3
open vertex: 3
examine neighbor: 3 -> 4 (vertexcolor = 0, edgecolor= 0)
discover vertex: 4
open vertex: 4
examine neighbor: 4 -> 5 (vertexcolor = 0, edgecolor= 0)
discover vertex: 5
open vertex: 5
examine neighbor: 5 -> 6 (vertexcolor = 0, edgecolor= 0)
discover vertex: 6
open vertex: 6
examine neighbor: 6 -> 7 (vertexcolor = 0, edgecolor= 0)
discover vertex: 7
open vertex: 7
examine neighbor: 7 -> 8 (vertexcolor = 0, edgecolor= 0)
discover vertex: 8
open vertex: 8
examine neighbor: 8 -> 2 (vertexcolor = 1, edgecolor= 0)
close vertex: 8
close vertex: 7
close vertex: 6
close vertex: 5
close vertex: 4
close vertex: 3
close vertex: 2
examine neighbor: 1 -> 3 (vertexcolor = 2, edgecolor= 0)
examine neighbor: 1 -> 4 (vertexcolor = 2, edgecolor= 0)
examine neighbor: 1 -> 5 (vertexcolor = 2, edgecolor= 0)
examine neighbor: 1 -> 6 (vertexcolor = 2, edgecolor= 0)
examine neighbor: 1 -> 7 (vertexcolor = 2, edgecolor= 0)
examine neighbor: 1 -> 8 (vertexcolor = 2, edgecolor= 0)
close vertex: 1

We can use Dijkstra’s Algorithm to calculate the distance from a given vertex to all other vertices in the graph. We see, for instance, that the distance from vertex 1 to vertex 4 is three steps. Since vertex 1 and vertex 20 are not connected, the distance between them is infinite. There are a couple of other algorithms available for calculating shortest paths.

julia> distances = ones(num_edges(g1a));   # Assign distance of 1 to each edge.
julia> d = dijkstra_shortest_paths(g1a, distances, 1);
julia> d.dists                             # Vector of distances to all other vertices.
20-element Array{Float64,1}:
   0.0
   1.0
   2.0
   3.0
   3.0
   2.0
   1.0
   1.0
   3.0
   4.0
   2.0
   2.0
 Inf  
 Inf  
 Inf  
 Inf  
 Inf  
 Inf  
 Inf  
 Inf  

As with the most of the packages that I have looked at already, the functionality summarised above is just a small subset of what’s available. Have a look at the home pages for these packages and check out the full code for today (which looks at a number of other features) on github. Some time in the future I plan on looking at the EvolvingGraphs which caters for graphs where the structure changes with time.

493-drawing-stars-0

The post #MonthOfJulia Day 24: Graphs appeared first on Exegetic Analytics.

#MonthOfJulia Day 23: Data Structures

Julia-Logo-DataStructure

Although Julia has integrated support for various data structures (arrays, tuples, dictionaries, sets), it doesn’t exhaust the full gamut of ptions. More exotic structures (like queues and deques, stacks, counters, heaps, tries and variations on sets and dictionaries) are implemented in the DataStructures package.

As always we start by loading the required package.

julia> using DataStructures

I won’t attempt to illustrate all structures offered by the package (that would make for an absurdly dull post), but focus instead on queues and counters. The remaining types are self-explanatory and well illustrated in the package documentation.

Queue

Let’s start off with a queue. The data type being queued must be specified at instantiation. We’ll make a queue which can hold items of Any type. Can’t get more general than that.

julia> queue = Queue(Any);

The rules of a queue are such that new items are always added to the back. Adding items is done with enqueue!().

julia> enqueue!(queue, "First in.");
julia> for n in [2:4]; enqueue!(queue, n); end
julia> enqueue!(queue, "Last in.")
Queue{Deque{Any}}(Deque [{"First in.",2,3,4,"Last in."}])
julia> length(queue)
5

The queue now holds five items. We can take a look at the items at the front and back of the queue using front() and back(). Note that indexing does not work on a queue (that would violate the principles of queuing!).

julia> front(queue)
"First in."
julia> back(queue)
"Last in."

Finally we can remove items from the front of the queue using dequeue!(). The queue implements FIFO (which is completely different from the other form of FIFO, which I only discovered today).

julia> dequeue!(queue)
"First in."

Counter

The counter() function returns an Accumulator object, which is used to assemble item counts.

julia> cnt = counter(ASCIIString)
Accumulator{ASCIIString,Int64}(Dict{ASCIIString,Int64}())

Using a Noah’s Ark example we’ll count the instances of different types of domestic animals.

julia> push!(cnt, "dog")                   # Add 1 dog
1
julia> push!(cnt, "cat", 3)                # Add 3 cats
3
julia> push!(cnt, "cat")                   # Add another cat (returns current count)
4
julia> push!(cnt, "mouse", 5)              # Add 5 mice
5

Let’s see what the counter looks like now.

julia> cnt
Accumulator{ASCIIString,Int64}(["mouse"=>5,"cat"=>4,"dog"=>1])

We can return (and remove) the count for a particular item using pop!().

julia> pop!(cnt, "cat")
4
julia> cnt["cat"]                          # How many cats do we have now? All gone.
0

And simply accessing the count for an item is done using [] indexing notation.

julia> cnt["mouse"]                        # But we still have a handful of mice.
5

I’ve just finished reading through the second early access version of Learn Julia by Chris von Csefalvay. In the chapter on Strings the author present a nice example in which he counts the times each character speaks in Shakespeare’s Hamlet. I couldn’t help but think that this would’ve been even more elegant using an Accumulator.

Tomorrow we’ll take a look at an extremely useful data structure: a graph. Until then, feel free to check out the full code for today on github.

The post #MonthOfJulia Day 23: Data Structures appeared first on Exegetic Analytics.