Resources for evangelism

Community Teaching & Outreach
1

I’m looking for resources for introducing Julia to people as a possible solution to challenges experienced in a python-based scientific-computing project. The challenges arise from inefficiency and lack of composability (this is actually part of my question, how to characterize and articulate the challenges of the current system as well as the opportunities offered by Julia). The format is an hour long presentation. The audience are relatively experienced scientists/engineers. My plan is to scroll through a Jupyter notebook. But, if people know of other notebooks or slides I can use, please let me know.

I don’t plan to dwell too much on the execution speed of Julia (no plots of microbenchmarks). I guess this is the easiest thing for people to understand and appreciate. I’ll integrate execution times in the story of other features.

The harder part is understanding why Julia is more productive (in many situations, at least). Regarding composability and multiple dispatch: Is the best thing to pull examples out of “The unreasonable effectiveness of multiple dispatch.” ? IIRC Chris R and Lyndon may have some related blog posts or presentations.

I want to use the situation in Python frequently as a counterpoint. I know we support the notion that all languages are beautiful in their own special way, which is laudable. My goal here is not to be antagonistic, but rather to explain specific technical benefits to using Julia in this kind of project, as well as the social benefit of happier, more enthusiastic developers (not least myself); some are quite happy using Python, others not so much. This is important, because you need to give a very compelling reason to justify adding another language to your stack.

Some topics are:

Specialization on input types, inference, optimization.

Heavy use of functions as first-class objects, devirtualizing/inlining functions passed as arguments because each has it’s own type. For example, adding 1 and then calling a passed function that subtracts one will be compiled into a no-op.

Many ways to compute over an array that are roughly equivalent for efficiency. Explicit loops (with nice iterators, like enumerate). Array comprehensions. Fused broadcast with dot notation: No functions are vectorized, and all functions are vectorized. (How much engineering goes into vectorization in the Python world ? I don’t have a detailed understanding: elementwise_grad, numpy.vectorize, etc.

Interop. @ccall, PyCall, @cfunction, PackageCompiler.jl, pyjulia.

Ways to illustrate the two-language problem. One story is the following about AD (which I know nothing about) I wrote LambertW.jl, which I chose only because I know it. Maybe it’s a good example in a way because the LambertW function is a bit obscure and is less likely to get the library support that cos does. With very little, almost no, effort, I ensured that lambertw is compatible with FowardDiff and Zygote. Not by supporting them specifically, but by making sure it’s written in a generic way. In the Python world, there is mpmath.lambertw which is orders of magnitude slower, and returns numbers that are not python floats or complexes, but something else. The complex one mpmath.ctx_mp_python.mpc fails isinstance(., complex) There is also scipy.special.lambertw which is written in cython. You can find the source online after some digging, it’s not installed by pip. It can accept a real number, but always returns a complex type. The type is, again, not a Python complex, but something else numpy.complex128, which nonetheless does satisfy isinstance(., complex). The latter function is as efficient as the Julia lambertw, but only when operating on a numpy array (it is vectorized). mpmath.lambertw will not accept an array. As far as I can tell, neither of the Python lambertw implementations is compatible with any AD package, because that would require explicit support. The wrapper (or reimplementation, or whatever) autograd.scipy.special does not implement lambertw. I don’t know about JAX. Furthermore, even if autograd did support lambertw, it would be as a separate function, in a hierarchy separated from other implementations of lambertw. There is more to the story, but you get the idea. In contrast, a single source LambertW.jl works efficiently with AD (and IntervalArithmetic.jl, etc) with no engineering effort and, in contrast to Python, incurring no complexity in multiple libraries. I don’t understand all the details of this story; maybe someone would like to clarify. But, that’s part of the point. I ensured LambertW works (and efficiently) with AD even though I know very little about it.

I should add that if someone has another story illustrating the same point, but with a less-obscure function, I’d be happy to hear it. Of course the lesson is not about lambertw, and only partly about AD. It’s that similar engineering hurdles will appear repeatedly in one language and either not appear at all or be much easier to clear in the other.

10 Likes
2

Here’s a notebook that I used for a general intro to Julia a few months ago: Jupyter Notebook Viewer

I spent a bit more time on performance than it sounds like you want to, but I also covered some of the magic of multiple dispatch and generics in Julia. You’re welcome to use any of that material if you find it useful.

By the way, I’d also suggest mentioning package management in Julia if your audience are scientists. Reproducibility is critical in science, and a Project.toml and Manifest.toml are a very nice and portable way to create a reproducible Julia environment. It’s possible to create a similarly reproducible environment in Python, but it is rarely actually done (in my experience).

11 Likes
3

1000% this. At some point I really need to write up my julia / conda #apprecigripes…

But more generally, I think the LambertW thing is perfect, because you wrote it. Your challenge will be too make it generic, for people that don’t know about our care about LambertW. The point is that scientists want to work on the cutting edge. Sometimes, you gotta write your own thing. Do you want that thing to be easy and slow, difficult and fast, or easy and fast. Maybe show the table of DiffEq algorithms (if Chris had an updated one) - why does Julia have everything, and why is it all fast? It’s not (only) because Chris is a wizard.

3 Likes
4

It would be interesting to compare the lines of code of various Julia libraries vs. their Python (or R) counterparts (including their C / Cython parts) of roughly equal functionality, e.g.

  • Julia base Arrays vs. Numpy
  • (DataFrames.jl + CSV.jl + Tables.jl) vs. Pandas
  • Genie.jl vs. Django
  • Flux vs. PyTorch
  • Pluto.jl vs. Jupyter Notebook

My bet would be that the Julia libraries are significantly smaller.

1 Like
5

This might not be the best tactic since the obvious response is the (sometimes true) “that’s because Julia libraries aren’t nature enough for our needs”

2 Likes
6

I don’t know how C code compares, but I’m 100% sure the Julia code would be substantially clearer for the average user. How many people using numpy understand C? How many people using Julia arrays understand Julia?

2 Likes
7

Thanks. I’ll put something in about that. My post got kind of long and I lost steam.

2 Likes
8

Oh. I’ll look for that table. I’ll also show (well, I’ll have to cut something out somewhere) the story of ODEs, uncertainty, and error bars.

9

At my institution and among the research groups that I interact with the most (applied mathematicians and engineers), the languages most used are Fortran and Matlab.

Most of the Fortran users do not get impressed by multiple dispatch, iterators, etc. They are convinced to try Julia when we show them that translating a Fortran code to Julia is something almost trivial, because you can write loops etc, that the syntax is simple and that, on top of that, there is a huge library of numerical packages to use (which are easier to fetch and use than libraries in fortran) and (yes, very important), Plots!*

The Matlab users here are convinced to try the first time they have to pay for the license again.

*Fortran users are sometimes reluctant because Fortran is there since ever. If Julia was called Fortran, they would have already moved enthusiastically :slight_smile:

9 Likes
10

Once the Fortran code being translated to Julia, you may show how it can be automatically differentiated : I guess that it can be pretty convincing.

I think that unique features are strong motivation for language adoption.

2 Likes
11

Sorry if this is already provided in some of the given resources, but I saw a pretty concise and convincing example of the composability and “everything is first class” in an interview with Jeff (sorry, I don’t remember where I saw the link, will add if I find it heres the link, thanks for posting it) which went something like this:

julia> using LinearAlgebra

julia> dmat = Diagonal(1:4) # Only stores the start and stop indices
4×4 Diagonal{Int64,UnitRange{Int64}}:
 1  ⋅  ⋅  ⋅
 ⋅  2  ⋅  ⋅
 ⋅  ⋅  3  ⋅
 ⋅  ⋅  ⋅  4

julia> dmat * randn(4,4) # Works with normal matrices ofc
4×4 Array{Float64,2}:
 -1.36422  -0.507138  -0.238564  -0.635628
  2.37113  -5.29489    1.42308    3.46297
 -4.95212   1.61776   -3.96482   -3.83821
  3.38688   2.9376     2.73302    4.6389

julia> struct OneHot{T} <: AbstractVector{T} # Lets make our own vector
       length::Int
       hot::Int
       end

julia> Base.size(v::OneHot) = (v.length,)

julia> Base.getindex(v::OneHot{T}, i) where T = i == v.hot ? one(T) : zero(T)

julia> Diagonal(OneHot{Float64}(4, 2)) # Also composes nicely
4×4 Diagonal{Float64,OneHot{Float64}}:
 0.0   ⋅    ⋅    ⋅
  ⋅   1.0   ⋅    ⋅
  ⋅    ⋅   0.0   ⋅
  ⋅    ⋅    ⋅   0.0

I was thinking one could expand on this and show how to specialize e.g. matrix multiplication with OneHot to be faster than the default.

1 Like
12

Note that hopefully this example will work even better with 1.6 as then Diagonal(OneHot) should automatically preserve the sparsity of OneHot

13

I think it was

Taking advantage of the structure of the one-hot vector is also discussed in

2 Likes
14

I had already worked in the same example from “The unreasonable …”. I added the example you mentioned, which makes an additional point. I think I have seen Jeff bring that up in yet another talk as well.

Also, I found the notebook that Stefan used to get the plot of the ODE with errors. I copied this (with attribution to the original authors) into my notebook.

1 Like
15

I encountered this example here, at around 13:00:

The video also shows how to define multiplication between a matrix and OneHotVector (it is just picking a row/column).

3 Likes
16

Contrasting it with Python seems sensible, that’s how I think many people would set up a modern simulation code nowadays – if they’re not aware of Julia. I can share a story in that regards: a year ago, we had the choice of starting a new project in either Python (making heavy use of Numba (JIT) and Dask (parallel/out-of-core scheduling)) or in Julia.

I started another Python-Numba-Dask another year before, and I ran into many hurdles. My issues were all due to having very poor availability for data structures, basically I needed an array of unequally sized structs holding some physical states. Every struct contains a number of arrays, which might grow or shrink. The math was rather simple.

Essentially, numba supports:

  • nD numpy arrays
  • numpy structured arrays (all equally sized per “row”)
  • classes holding arrays (still labelled experimental)

I ended up re-implementing my code two times along these data structures, each time running into limitations (on the upside, I’ve gotten a fairly good feel for the limitations of numba). The class-based approach still has a great deal of boilerplate to get numba’s type-inferencing to work properly. Performance also seems to be only around 25-30 times better than interpreted Python, which seems on the low side for numba – I’m pretty sure these abstractions are not without cost in this case.

Alternatively, I could’ve chosen C++, Fortran, Cython, or plain C. These are a much greater pain to distribute (especially with compilers and Python on Windows!), and I’d still have to write/generate a good amount of boilerplate. I’d also have to instruct colleagues who mostly only know Python.

In Julia:

  • I could’ve taken the obvious approach (array of structs)
  • I would’ve probably cut the number of lines half
  • I would’ve gotten better performance to boot

So even without quality-of-life features such as multiple dispatch, you can market Julia as a “better” Python which gives good performance without having to jump through many hoops: you can start a project with full confidence that you won’t run into (almost) crippling limitations which require large rewrites (possibly in a different language). Maybe you could find a colleague who started a non-trivial project in Python, experienced similar difficulties, and contrast it with how pleasant life is using Julia instead.

3 Likes