Julia: A numerical programming language

Overview

  • Very recent: v1.0. released August 2018
  • Compiled scripting language
  • JIT with LLVM backend
  • High-level syntax, dynamical
  • Garbage collected
  • Strongly typed
  • A bit like python with integrated numpy
  • Multiple dispatch
  • Rich interoperability: Fortran, C, C++, python, R, ...
  • Tries to facilitate modern developments:
    • Machine learning
    • Large-scala data analysis
    • Automatic differentiation)

Functions and multiple dispatch

In [ ]:
function mysquare(x)
    x * x
end

# or alternatively just:
mysquare(x) = x * x
In [ ]:
mysquare(3)    # Integer arithmetic
In [ ]:
mysquare(2im)  # Complex arithmetic

Lists == Arrays == Vectorisation

In [ ]:
typeof(["abc"])
In [ ]:
typeof([1, 2, 3])
In [ ]:
typeof([1.1, 2, 3])
In [ ]:
typeof(["abc", 1, 2.3])
In [ ]:
arr = [1.1, 2, 3]

mysquare.(arr)   # Vectorised call
In [ ]:
sin.(exp.(mysquare.(arr)))

# or:
@. sin(exp(mysquare(arr)))
In [ ]:
using Printf

function printarray(array::Array{T, 1}) where T
    println("General array passed, elements:")
    for element in array
        println("   $element")
    end
end

function printarray(array::Array{T, 1}) where T <: Number
    println("Numeric array passed, elements:")
    for element in array
        println("   $element")
    end
end

function printarray(array::Array{Float64, 1})
    println("Float array passed, elements:")
    for element in array
        @printf "  %10.5f\n" element
    end
end

printarray(["a", "b", 1])
println()
printarray([1, 2, 3])
println()
printarray([1., 2., 3.])

Changing precision

In [ ]:
Type = Float16

a = rand(Type, 3, 4)
sum(mysquare.(a))

Changing computational backend

In [ ]:
using StaticArrays

# Array operations on static-sized arrays:
a = @SArray[1,2,3,4]
mysquare.(a)
In [ ]:
using SparseArrays

# random sparse array
a = sprandn(6, 6, 0.3)
mysquare.(a)
In [ ]:
using Distributed

# Add 4 processes and setup environment at each of them
#addprocs(4)

@everywhere begin
    mysquare(x) = x * x
    using DistributedArrays
end

a = drandn(100,100,10)
println(sum(mysquare.(a)))
close(a)

Structs and inheritance

  • A struct is a collection of data
In [ ]:
struct MyPoint
    x::Float64
    y::Float64
end

Note: A struct contains no methods and is immutable by default.

In [ ]:
p = MyPoint(1, 2)
p.x = 6

Functionality is provided by external functions only:

In [ ]:
pointlength(p::MyPoint) = sqrt(p.x^2 + p.y^2)
pointlength(p)

Abstract types are only names to "group" concrete structs:

In [ ]:
using StaticArrays

@show AbstractArray{Float64} <: AbstractArray
@show SArray                 <: AbstractArray

@show Array{Float64}         <: AbstractArray{Float64}
@show Array{Int64}           <: AbstractArray{Float64};
In [ ]:
# Function, which works on all Abstract Arrays of element type Float64
function strangesum(a::AbstractArray{Float64})
    a[1] + sum(a[2:end])
end
In [ ]:
# ok:
strangesum([1., 2., 3.])
strangesum(@SArray[1., 2., 3.])
In [ ]:
# not ok:
strangesum([1, 2, 3])

Automatic differentiation (Adjoint mode)

  • FD and Forward well-established
  • Adjoint mode: Many packages, Zygote very promising, but experimental
In [ ]:
using Zygote
using SpecialFunctions

f(x) = 2x*sin(x) - x^2/2
In [ ]:
@time f'(1)
In [ ]:
@time f'(2)

Integration with existing scientific codes

  • e.g. python, C++, Fortan, R
In [ ]:
# Call to lapack to compute dot product
v = [1.0, 2.0, 3.0, 4.0]
w = [2.0, 3.0, 4.0, 5.0]

VdotW = ccall((:ddot_, "liblapack"), Float64,
              (Ref{Int32}, Ptr{Float64}, Ref{Int32}, Ptr{Float64}, Ref{Int32}),
              length(v), v, 1, w, 1)
println("v ⋅ w = $VdotW")
In [ ]:
# Plot function in matplotlib
using PyPlot

# f defined as above
x = collect(-4:0.2:7)
PyPlot.plot(x, f.(x), label="f(x) = 2x*sin(x) - x^2/2 ")
PyPlot.plot(x, f'.(x), label="f'(x)")
PyPlot.axhline(0, color="grey", linewidth=0.5)
PyPlot.legend()
PyPlot.show()

A word about performance

  • Best out of five
  • Heat equation example (based on code samples by Antoine Levitt)
  • Used software: gcc 8.3, clang 7.0.1, python 3.7, numpy 1.16.2, julia 1.0.3
         |                           |   duration (s)
-------- | ------------------------- | ---------------
python   | array operation (numpy)   |           11.8
         |                           |                 
C        | gcc                       |            8.1
C        | gcc -O3                   |            1.1
C        | gcc -O3 -march=native     |            0.5
         |                           |                
C        | clang                     |            8.0
C        | clang -O3                 |            1.1
C        | clang -O3 -march=native   |            0.8
         |                           |                
julia    | array operations          |            3.3
julia    | loops                     |            1.9
julia    | loops, no bounds check    |            0.5

Wrapping up

  • Stronger emphasis on functional-style programming
  • Multiple dispatch based on JIT-compiling and strong typing
  • Two-way interoperability with Fortran, C / C++, python, R
  • Threading, vectorisation and distributed memory in high-level syntax
  • Uses cases:
    • Great language for algorithmic developments
    • High-level alternative for Fortran and C in high-performance computing
    • Helpful for one-shot scripts (easy parallelisation)
    • Not so much for classical scripting (glue code)