Tag Archives: Programming

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.

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.

What’s Wrong with Statistics in Julia?

By: John Myles White

Re-posted from: http://www.johnmyleswhite.com/notebook/2014/11/29/whats-wrong-with-statistics-in-julia/

Introduction

Several months ago, I promised to write an updated version of my old post, “The State of Statistics in Julia”, that would describe how Julia’s support for statistical computing has evolved since December 2012.

I’ve kept putting off writing that post for several reasons, but the most important reason is that all of my attention for the last few months has been focused on what’s wrong with how Julia handles statistical computing. As such, the post I’ve decided to write isn’t a review of what’s already been done in Julia, but a summary of what’s being done right now to improve Julia’s support for statistical computing.

In particular, this post focuses on several big changes to the core data structures that are used in Julia to represent statistical data. These changes should all ship when Julia 0.4 is released.

What’s Wrong with Statistics in Julia Today?

The primary problem with statistical computing in Julia is that the current tools were all designed to emulate R. Unfortunately, R’s approach to statistical computing isn’t amenable to the kinds of static analysis techniques that Julia uses to produce efficient machine code.

In particular, the following differences between R and Julia have repeatedly created problems for developers:

  • In Julia, computations involving scalars are at least as important as computations involving vectors. In particular, iterative computations are first-class citizens in Julia. This implies that statistical libraries must allow developers to write efficient code that iterates over the elements of a vector in pure Julia. Because Julia’s compiler can only produce efficient machine code for computations that are type-stable, the representations of missing values, categorical values and ordinal values in Julia programs must all be type-stable. Whether a value is missing or not, its type must remain the same.
  • In Julia, almost all end-users will end up creating their own types. As such, any tools for statistical computing must be generic enough that they can be extended to arbitrary types with little to no effort. In contrast to R, which can heavily optimize its algorithms for a very small number of primitive types, Julia developers must ensure that their libraries are both highly performant and highly abstract.
  • Julia, like most mainstream languages, eagerly evaluates the arguments passed to functions. This implies that idioms from R which depend upon non-standard evaluation are not appropriate for Julia, although it is possible to emulate some forms of non-standard evaluation using macros. In addition, Julia doesn’t allow programmers to reify scope. This implies that idioms from R that require access to the caller’s scope are not appropriate for Julia.

The most important way in which these issues came up in the first generation of statistical libraries was in the representation of a single scalar missing value. In Julia 0.3, this concept is represented by the value NA, but that representation will be replaced when 0.4 is released. Most of this post will focus on the problems created by NA.

In addition to problems involving NA, there were also problems with how expressions were being passed to some functions. These problems have been resolved by removing the function signatures for statistical functions that involved passing expressions as arguments to those functions. A prototype package called DataFramesMeta, which uses macros to emulate some kinds of non-standard evaluation, is being developed by Tom Short.

Representing Missing Values

In Julia 0.3, missing values are represented by a singleton object, NA, of type NAtype. Thus, a variable x, which might be either a Float64 value or a missing value encoded as NA, will end up with type Union(Float64, NAtype). This Union type is a source of performance problems because it defeats Julia’s compiler’s attempts to assign a unique concrete type to every variable.

We could remove this type-instability by ensuring that every type has a specific value, such as NaN, that signals missingness. This is the approach that both R and pandas take. It offers acceptable performance, but does so at the expense of generic handling of non-primitive types. Given Julia’s rampant usage of custom types, the sentinel values approach is not viable.

As such, we’re going to represent missing values in Julia 0.4 by borrowing some ideas from functional languages. In particular, we’ll be replacing the singleton object NA with a new parametric type Nullable{T}. Unlike NA, a Nullable object isn’t a direct scalar value. Rather, a Nullable object is a specialized container type that either contains one value or zero values. An empty Nullable container is taken to represent a missing value.

The Nullable approach to representing a missing scalar value offers two distinct improvements:

  • Nullable{T} provides radically better performance than Union(T, NA). In some benchmarks, I find that iterative constructs can be as much as 100x faster when using Nullable{Float64} instead of Union(Float64, NA). Alternatively, I’ve found that Nullable{Float64} is about 60% slower than using NaN to represent missing values, but involves a generic approach that trivially extends to arbitrary new types, including integers, dates, complex numbers, quaternions, etc…
  • Nullable{T} provides more type safety by requiring that all attempts to interact with potentially missing values explicitly indicate how missing values should be treated.

In a future blog post, I’ll describe how Nullable works in greater detail.

Categorical Values

In addition to revising the representation of missing values, I’ve also been working on revising our representation of categorical values. Working with categorical data in Julia has always been a little strange, because the main tool for representing categorical data, the PooledDataArray, has always occupied an awkward intermediate position between two incompatible objectives:

  • A container that keeps track of the unique values present in the container and uses this information to efficiently represent values as pointers to a pool of unique values.
  • A container that contains values of a categorical variable drawn from a well-defined universe of possible values. The universe can include values that are not currently present in the container.

These two goals come into severe tension when considering subsets of a PooledDataArray. The uniqueness constraint suggests that the pool should shrink, whereas the categorical variable definition suggests that the pool should be maintained without change. In Julia 0.4, we’re going to commit completely to the latter behavior and leave the problem of efficiently representing highly compressible data for another data structure.

We’ll also begin representing scalar values of categorical variables using custom types. The new CategoricalVariable and OrdinalVariable types that will ship with Julia 0.4 will further the efforts to put scalar computations on an equal footing with vector computations. This will be particularly notable for dealing with ordinal variables, which are not supported at all in Julia 0.3.

Metaprogramming

Many R functions employ non-standard evaluation as a mechanism for augmenting the current scope with the column names of a data.frame. In Julia, it’s often possible to emulate this behavior using macros. The in-progress DataFramesMeta package explores this alternative to non-standard evaluation. We will also be exploring other alternatives to non-standard evaluation in the future.

What’s Next

In the long-term future, I’m hoping to improve several other parts of Julia’s core statistical infrastructure. In particular, I’d like to replace DataFrames with a new type that no longer occupies a strange intermediate position between matrices and relational tables. I’ll write another post about those issues later.