Tag Archives: Julia

Julia iFEM1: Porting Mesh Generation

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/julia-ifem-1-mesh-generation-and-julias-type-system/

My first project on the quest for a Julia finite element method is a simple homework problem. Just some background, this is for UC Irvine’s graduate Computational PDEs 226B course where in the first quarter we did all sorts of finite difference methods and now is our first foray into finite element methods. The purpose of the project is to grasp the data structure enough to use simple tools (i.e. mesh creation and plotting) to create a finite element solver for Poisson’s equation in 2D and check the performance differences. For those who haven’t programmed much this is a great learning exercise, but being a pretty standard exercise the MATLAB code took no time to writeup and so I started thinking: how does Julia compare to MATLAB for solving simple PDEs with the finite element method?

Testing this would take a few steps. The code base is all in MATLAB, so there is a step of (partial) porting MATLAB to Julia. Then there is a fiasco as how to plot the data structures given that our display commands are different (i.e. is there a way to plot in Julia that’s almost exactly like MATLAB?). Then we have to optimize the Julia code and see how it fares against the MATLAB code. This post will focus on how I did the first two steps (and the problems I encountered), the last will be its own post later.

Semi-Porting to Julia: MATLAB.jl and number casting

The moment I looked at Julia code I was struck by how similar it looked to MATLAB code. The first part of the problem is to make a uniform square grid, which is easy enough to be just a few MATLAB commands, so my first response to porting it into Julia was to simply copy-paste it, find out what caused errors, and beat through the errors until it worked. So I downloaded the Julia+Juno IDE bundle and started hacking away at it. The first thing I noticed was the Julia did not have a command for meshgrid in its base library. However, not surprisingly enough people want this command that someone already implemented it and ndgrid. So I took that file and it played as a template for the syntax as well.

With meshgrid in hand, things went really smooth. Array dereferencing in Julia uses square brackets, so I had to go through and change X(:,1) to X[:,1] and the like, but that was easy enough. Another issue is that in MATLAB you can drop parts of the array using

t2nidxMap(topNode) = [];

whereas in Julia I used the deleteat function:

t2nidxMap = deleteat!(collect(t2nidxMap),collect(topNode));

You may be asking, what is that peculiar “collect” function? Well in Julia arrays defined by ranges (i.e. A=1:100) return a range object and not an array. In most cases it seems that they can be used just like an array (with better performance and less memory usage), but this seems like a case where they could not, and so the collect function builds the array from the range object. This gave me the code to make the square mesh as:

function squaremesh(square,h)
  #square = [x0 x1 y0 y1], h is stepsize
  x0 = square[1]; x1= square[2]
  y0 = square[3]; y1= square[4]
  x,y = meshgrid(x0:h:x1,y0:h:y1)
  node = [x[:] y[:]];
 
  ni = size(x,1); # number of rows
  N = size(node,1);
  t2nidxMap = 1:N-ni;
  topNode = ni:ni:N-ni;
  t2nidxMap = deleteat!(collect(t2nidxMap),collect(topNode))
  k = t2nidxMap;
  elem = [k+ni k+ni+1 k ; k+1 k k+ni+1]
  return node,elem
end

whereas the MATLAB code was

%% Generate nodes
x0 = square(1); x1 = square(2); 
y0 = square(3); y1 = square(4);
[x,y] = meshgrid(x0:h:x1,y0:h:y1);
node = [x(:),y(:)];
 
%% Generate elements
ni = size(x,1); % number of rows
N = size(node,1);
t2nidxMap = 1:N-ni;
topNode = ni:ni:N-ni;
t2nidxMap(topNode) = [];
k = (t2nidxMap)';
elem = [k+ni k+ni+1 k; k+1 k k+ni+1];

As you can see, porting that function over took just a few minutes and very few things changed. It would take less than a minute if I knew about the range and deleteat issues. In fact, for next many things I did in Julia, those (and square brackets) were really the only differences, other than the plots.

Plotting via matplotlib

The next part of the project is to plot the triangulation. The MATLAB code for the “showmesh” throws this directly into MATLAB’s trisurf function. Thus I Google’d Julia trisurf and the first hit was Julia’s PyPlot package. It didn’t tell me how to do it, but the package simply calls matplotlib, so I went over to the matplotlib page and and found their function conveniently named plot_trisurf. All I had to do was pick a better color scheme than the default and I was good to go. So while the MATLAB code called

trisurf(elem(:,1:3),node(:,1),node(:,2),zeros(size(node,1),1));

to get the code on the project webpage, the Julia code called

plot_trisurf(node[:,1],node[:,2],zeros(size(node,1)),cmap=get_cmap("ocean"))

to give me this nifty 3D plot which shows the triangalization:

juliaTriangulation

Simple enough.

Getting lazy: MATLAB time?

Okay, things are pretty easy. Now we have to make a circle mesh. So I go into the code only to find that it calls a much larger MATLAB package for mesh generation. This is the part where most people give up on the new language: this part is done in MATLAB and looks like it would take more than the next 15 minutes to port… could I just call MATLAB? If you’ve ever dealt with MEX or similar things that “glue together languages”, you they are a mess. But I found the MATLAB.jl package so I decided to give it a try. The directions are dead simple: just stick the folder from the GitHub repository into the spot they tell you and you’re done. So I fire it up and… it didn’t work. But recall that in MATLAB numbers always promote to the highest state, i.e.

int8(5)+4.4

returns 9. We will return to this later when talking about efficiency, but the problem I was having boiled down this problem. All I had to do was pass in the same array explicitly telling it the numbers were floats by putting a . after it (i.e. 1 vs 1.) and all of the MATLAB codes worked. I was able to mix and matching generating meshes in Julia/MATLAB and plotting them in Julia/MATLAB. An example looks like this. Here I use the iFEM package functions and some of their Julia translations I described before, and mix and match calling them in Julia and MATLAB.

# Testing Square Mesh Generation
node,elem = squaremesh([0 1 0 1],.1)
showmesh(node,elem);
 
# Testing MATLAB call
node2,elem2 = mxcall(:squaremesh,2,[0. 1. 0. 1.],.1) #Array has periods to turn to Float
showmesh(node2,elem2)
 
# Testing MATLAB showmesh
mxcall(:showmesh,0,node,elem) #surpress return of function handle, cannot be saved STDIO for Julia return
mxcall(:showmesh,0,node2,elem2)
 
# Testing Quadralateral showmesh
C = .5*zeros(length(node[:,1]),length(node[:,2]))
D = zeros(length(node[:,1]))
p = pcolormesh(node[:,1],node[:,2],C,edgecolor=D,cmap="ocean")
 
# Doing circle and refine via MATLAB, pyplot plotting
node,elem = mxcall(:circlemesh,2,0.,0.,1.,0.2)
showmesh(node,elem)
node,elem = mxcall(:uniformrefine,2,node,elem)
showmesh(node,elem)
 
## 3D Mesh Generation and Plotting
node3,elem3 = mxcall(:cubemesh,2,[0. 1. 0. 1. 0. 1.],.25)
showmesh(node3,elem3)

You notice that to do MATLAB calls you simply use mxcall, tell it the function (with a colon in front, in Julia this means you converted it to the symbol type), the number of outputs, and give it the inputs. Be careful about integer vs floating point problems and everything works.

Conclusion: Moving to Julia is dead simple, and you don’t have to move too fast

This really appealed to me in two ways. First of all, Julia’s syntax was very natural to pick up and within minutes I was fine. In fact, the only errors that were hard to diagnose were the ones due to MATLAB! That leads me to the second point, Julia’s language interfacing is so good that I didn’t even have to stray too far. If I wanted to use old MATLAB code, it took 2 minutes to get it setup and I was calling it from Julia. The same seems to be true for calling Python and calling R as well. Not only that, but using C and Fortran code is built right into the core of Julia and has the same syntax structure as mxcall (instead you use ccall), meaning it’s easier to use than MATLAB’s MEX interface. Julia melds together different languages so nicely that you can act like you have the packages of all of them! All of this ease of use comes down Julia’s type system, which I will start looking into in my next post.

So was it worth it? As of right now I can’t say that, but I can say it didn’t cost me more than 30 minutes to get up and running in Julia. But will it be better? We will look at Julia’s type system and performance next time.

The post Julia iFEM1: Porting Mesh Generation appeared first on Stochastic Lifestyle.

WordPress Julia Syntax Highlighting

By: Christopher Rackauckas

Re-posted from: http://www.stochasticlifestyle.com/wordpress-julia-syntax-highlighting/

using Julia
print("Test for Julia")
A = [0 1 ; 0 1];

Syntax highlighting in Julia exists? Yes, you can make it work even though no plugins currently support it. The way I was able to get it to work was to use the WP-GeSHi-Highlight WordPress plugin which uses GeSHi to perform the syntax highlighting. Sure enough, an unofficial Julia syntax file for GeSHi is hosted on GitHub. I manually plopped this file into the GeSHi folder in my WordPress install, used the command lang=”julia” and you can see at the top of this page that it worked!

The post WordPress Julia Syntax Highlighting appeared first on Stochastic Lifestyle.

Why Julia’s DataFrames are Still Slow

By: John Myles White

Re-posted from: http://www.johnmyleswhite.com/notebook/2015/11/28/why-julias-dataframes-are-still-slow/

Introduction

Although I’ve recently decided to take a break from working on OSS for a little while, I’m still as excited as ever about Julia as a language.

That said, I’m still unhappy with the performance of Julia’s core data analysis infrastructure. The performance of code that deals with missing values has been substantially improved thanks to the beta release of the NullableArrays package, which David Gold developed during this past Julia Summer of Code. But the DataFrames package is still a source of performance problems.

The goal of this post is to explain why Julia’s DataFrames are still unacceptably slow in many important use cases — and will remain slow even after the current dependency on the DataArrays package is replaced with a dependency on NullableArrays.

Problematic Interactions with Julia’s Compiler

The core problem with the DataFrames library is that a DataFrame is, at its core, a black-box container that could, in theory, contain objects of arbitrary types. In practice, a DataFrame contains highly constrained objects, but those constraints are (a) hard to express to the compiler and (b) still too weak to allow the compiler to produce the most efficient machine code.

The use of any black-box container creates the potential for performance problems in Julia because of the way that Julia’s compiler works. In particular, Julia’s compiler is able to execute code quickly because it can generate custom machine code for every function call — and this custom machine code is specialized for the specific run-time types of the function’s arguments.

This run-time generation of custom machine code is called specialization. When working with black-box containers, Julia’s approach to specialization is not used to full effect because machine code specialization based on run-time types only occurs at function call sites. If you access objects from a black-box container and then perform extended computations on the results, those computations will not be fully specialized because there is no function call between (a) the moment at which type uncertainty about the contents of the black-box container is removed and (b) the moment at which code that could benefit from type information is executed.

A Minimal Example

To see this concern in practice, consider the following minimal example of a hot loop being executed on values that are extracted from a black-box container:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
function g1(black_box_container)
    x, y = black_box_container[1], black_box_container[2]
    n = length(x)
    s = 0.0
    for i in 1:n
        s += x[i] * y[i]
    end
    s
end
 
function hot_loop(x, y)
    n = length(x)
    s = 0.0
    for i in 1:n
        s += x[i] * y[i]
    end
    s
end
 
function g2(black_box_container)
    x, y = black_box_container[1], black_box_container[2]
    hot_loop(x, y)
end
 
container = Any[randn(10_000_000), randn(10_000_000)];
 
@time g1(container)
# 2.258571 seconds (70.00 M allocations: 1.192 GB, 5.03% gc time)
 
@time g2(container)
# 0.015286 seconds (5 allocations: 176 bytes)

g1 is approximately 150x slower than g2 on my machine. But g2 is, at a certain level of abstraction, exactly equivalent to g1 — the only difference is that the hot loop in g1 has been put inside of a function call. To convince yourself that the function call boundary is the only important difference between these two functions, consider the following variation of g2 and hot_loop:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
@inline function hot_loop_alternative(x, y)
    n = length(x)
    s = 0.0
    for i in 1:n
        s += x[i] * y[i]
    end
    s
end
 
function g3(black_box_container)
    x, y = black_box_container[1], black_box_container[2]
    hot_loop_alternative(x, y)
end
 
@time g1(container)
# 2.290116 seconds (70.00 M allocations: 1.192 GB, 4.90% gc time)
 
@time g2(container)
# 0.017835 seconds (5 allocations: 176 bytes)
 
@time g3(container)
# 2.250301 seconds (70.00 M allocations: 1.192 GB, 5.08% gc time)

On my system, forcing the hot loop code to be inlined removes all of the performance difference between g1 and g2. Somewhat ironically, by inlining the hot loop, we’ve prevented the compiler from generating machine code that’s specialized on the types of the x and y values we pull out of our black_box_container. Inlining removes a function call site — and function call sites are the only times when machine code can be fully specialized based on run-time type information.

This problem is the core issue that needs to be resolved to make Julia’s DataFrames as efficient as they should be. Below I outline three potential solutions to this problem. I do not claim that these three are the only solutions; I offer them only to illustrate important issues that need to be addressed.

Potential Solutions to the Under-Specialization Problem

One possible solution to the problem of under-specialization is to change Julia’s compiler. I think that work on that front could be very effective, but the introduction of specialization strategies beyond Julia’s current "specialize at function call sites" would make Julia’s compiler much more complex — and could, in theory, make some code slower if the compiler were to spend more time performing compilation and less time performing the actual computations that a user wants to perform.

A second possible solution is to generate custom DataFrame types for every distinct DataFrame object. This could convert DataFrames from black-box containers that contain objects of arbitrary type into fully typed containers that can only contain objects of types that are fully known to the compiler.

The danger with this strategy is that you could generate an excessively large number of different specializations — which would again run the risk of spending more time inside the compiler than inside of the code you actually want to execute. It could also create excessive memory pressure as an increasing number of specialized code paths are stored in memory. Despite these concerns, a more aggressively typed DataFrame might be a powerful tool for doing data analysis.

The last possible solution I know of is the introduction of a high-level API that ensures that operations on DataFrames always reduce down to operations on objects whose types are known when hot loops execute. This is essentially the computational model used in traditional databases: take in a SQL specification of a computation, make use of knowledge about the data actually stored in existing tables to formulate an optimized plan for performing that computation, and then perform that optimized computation.

I think this third option is the best because it will also solve another problem Julia’s data infrastructure will hit eventually: the creation of code that is insufficiently generic and not portable to other backends. If people learn to write code that only works efficiently for a specific implementation of DataFrames, then their code will likely not work when they try to apply it to data stored in alternative backends (e.g. traditional databases). This would trap users into data structures that may not suit their needs. The introduction of a layer of appropriate abstractions (as in dplyr) would resolve both issues at once.

Take-Aways

  • Making Julia’s DataFrames better is still a work-in-progress.
  • The core issue is still the usage of data structures that are not amenable to Julia’s type inference machinery. One of the two main issues is now resolved; another must be addressed before things function smoothly.
  • Several solutions to this remaining are possible; we will probably see one or more of these solutions gain traction in the near-term future.