Why Julia?

• Started in 2009 and first version was released in 2012.
• High-level languages like Python and R let one explore and experiment rapidly, but can run slow.
• Low-level languages like Fortran/C++ tend to take longer to develop, but run fast.
• This is sometimes called the "two language problem" and is something the Julia developers set out to eliminate.
• Julia's promise is to provide a "best of both worlds" experience for programmers who need to develop novel algorithms and bring them into production environments with minimal effort.

Julia: A Fresh Approach to Numerical Computing

Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah

SIAM Rev., 59(1), 65–98. (34 pages) 2012

Julia's Engineering and Design Tradoffs

• Type structures cannot be changed after being created (less dynamism but memory layout can be optimized for)
• All functions are JIT compiled via LLVM (interactive lags but massive runtime improvements)
• All functions specialize on types of arguments (more compilation but give generic programming structures)
• Julia is interactive (use it like Python and R, but makes it harder to get binaries)
• Julia has great methods for handling mutation (more optimization opportunities like C/Fortran, but more cognitive burden)
• Julia's Base library and most packages are written in Julia, (you can understand the source, but learn a new package)
• Julia has expensive tooling for code generation and metaprogramming (concise and more optimizations, but some codes can be for experienced users)

To me, this gives me a language with a lot of depth which works well for computationally-expensive scientific applications.

© ChrisRackaukas

Type-Dispatch Programming

• Centered around implementing the generic template of the algorithm not around building representations of data.
• The data type choose how to efficiently implement the algorithm.
• With this feature share and reuse code is very easy

JuliaCon 2019 | The Unreasonable Effectiveness of Multiple Dispatch | Stefan Karpinski

Time Series : create new type in Julia

"""
    TimeSeries{T,N}( nt )
- `T` : type of values
- `N` : number of variables
- `nt` : number of values
"""
struct TimeSeries{T,N}
   nt  :: Int
   nv  :: Int
   t   :: Vector{T}
   u   :: Vector{Array{T, 1}}
   function TimeSeries{T,N}( nt :: Int) where {T,N}
       t  = zeros(T, nt)
       u  = [zeros(T, N) for i in 1:nt]
       nv = N
       new( nt, nv, t, u)
   end
end

Main.var"ex-index".TimeSeries

@doc TimeSeries

TimeSeries{T,N}( nt )

• T : type of values
• N : number of variables
• nt : number of values

nt, nv = 100, 2
ts = TimeSeries{Float64, nv}(nt)
ts.nv, ts.nt

(2, 100)

Overload Base.length function

Overload Base.length function

import Base:length
length(ts :: TimeSeries) = ts.nt
nt, nv = 100, 2
ts = TimeSeries{Float64, nv}(nt);
@show length(ts) == nt

true

Generate data

ts.t[1] = 0.0
ts.u[1] = [0.0, 1.0]
dt = 0.01
for i in 2:nt
   ts.t[i] = ts.t[i-1] + dt
   ts.u[i][1] = sin(ts.t[i])
   ts.u[i][2] = cos(ts.t[i])
end

using Plots
plot(ts.t, vcat(ts.u'...))

plot(ts.t, [getindex.(ts.u, i) for i in 1:nv])

Overload the [] operator

we want ts[i] equal to ts.u[:][i] values

Overload the [] operator

we want ts[i] equal to ts.u[:][i] values

import Base: getindex

Overload the [] operator

we want ts[i] equal to ts.u[:][i] values

import Base: getindex

getindex( ts :: TimeSeries, i ) = getindex.(ts.u, i)

getindex (generic function with 284 methods)

plot(ts[1])
plot!(ts[2])

Plot recipe

using RecipesBase
@recipe f(ts::TimeSeries) = vcat(ts.u'...)
plot(ts)

Overload the + operator to add noise

import Base:+

Overload the + operator to add noise

import Base:+

function +(ts :: TimeSeries, ϵ )
    for n in 1:ts.nt, d in 1:ts.nv
       ts.u[n][d] += ϵ[n,d]
    end
    return ts
end

+ (generic function with 266 methods)

ts = ts + 0.1*randn((nt,nv));

ts = ts + 0.1*randn((nt,nv));

scatter(ts.t, [ts[1],ts[2]])

Linear regression with obvious operation

using LinearAlgebra
X = hcat(ones(nt), ts.t, ts[1])
y = ts[2]
@show β = inv(X'X) * X'y

3-element Vector{Float64}:
  1.0997092611654475
 -0.4491752626990362
 -0.043450126831331204

Version with QR factorisation

@show β = X \ y

3-element Vector{Float64}:
  1.0997092611654484
 -0.4491752626990573
 -0.04345012683131561

The \ operator is the short-hand for

Q, R = qr(X)
@show β = (inv(factorize(R))Q') * y

3-element Vector{Float64}:
  1.0997092611654484
 -0.4491752626990569
 -0.043450126831315806

Version with singular values decomposition

U, S, V = svd(X)
@show β = V * Diagonal(1 ./ S) * U' * y

3-element Vector{Float64}:
  1.0997092611654486
 -0.44917526269905755
 -0.04345012683131474

Version with singular values decomposition

U, S, V = svd(X)
@show β = V * Diagonal(1 ./ S) * U' * y

3-element Vector{Float64}:
  1.0997092611654486
 -0.44917526269905755
 -0.04345012683131474

@show β = pinv(X, atol=1e-8) * y

3-element Vector{Float64}:
  1.0997092611654489
 -0.44917526269905783
 -0.043450126831314974

With GLM.jl

using GLM
fitted = lm(X, y)

GLM.LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}, Vector{Int64}}}}:
Coefficients:
─────────────────────────────────────────────────────────────────
         Coef.  Std. Error      t  Pr(>|t|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────
x1   1.09971     0.0220369  49.90    <1e-70   1.05597    1.14345
x2  -0.449175    0.100115   -4.49    <1e-04  -0.647876  -0.250475
x3  -0.0434501   0.106612   -0.41    0.6845  -0.255045   0.168145
─────────────────────────────────────────────────────────────────

Clustering

using CluGen, Plots
n, k = 500, 4
o = clugen(2, k, n, [1, 0], π / 2, [10, 10], 10, 1, 1)
plot(o.points[:, 1], o.points[:, 2], seriestype = :scatter, group=o.clusters)

using BenchmarkTools
@btime omin = minimum(view(o.points,:,1))
@btime omin = minimum(o.points[:,1])
@btime @views omin = minimum(o.points[:,1])

  559.805 ns (4 allocations: 160 bytes)
  1.226 μs (4 allocations: 4.17 KiB)
  665.025 ns (5 allocations: 176 bytes)

omin_x = minimum(o.points[:,1])
o.points[:, 1] .-= omin_x
omax_x = maximum(o.points[:,1])
o.points[:, 1] ./= omax_x
omin_y = minimum(o.points[:,2])
o.points[:, 2] .-= omin_y
omax_y = maximum(o.points[:,2])
o.points[:, 2] ./= omax_y

500-element view(::Matrix{Float64}, :, 2) with eltype Float64:
 0.05150913591429948
 0.05076195000993125
 0.08211837757379832
 0.07432054289550463
 0.09184561936232206
 0.1026934742700536
 0.13366228807596342
 0.07408076732799851
 0.07130896101332615
 0.06047619216682257
 ⋮
 0.5832383569369922
 0.5837470305697766
 0.5822081564156026
 0.5875960498028272
 0.5826117734565991
 0.5465090690792758
 0.6068451125993966
 0.5921324601450353
 0.6096462939860531

plot(o.points[:, 1], o.points[:, 2], seriestype = :scatter, group=o.clusters)

Bregman divergence in R

using RCall
R"""
bregman <- function(x,y){
    distance <- function(x,y){
        if (x==0) {
            return(y)
        }
        else {
            return(x * log(x) - x + y - x * log(y))
        }
    }
    return(sum(divergences = mapply(distance, x, y)))
}
"""

RCall.RObject{RCall.ClosSxp}
function (x, y) 
{
    distance <- function(x, y) {
        if (x == 0) {
            return(y)
        }
        else {
            return(x * log(x) - x + y - x * log(y))
        }
    }
    return(sum(divergences = mapply(distance, x, y)))
}

Bregman divergence in Julia

"""
    bregman(x, y)
Compute Bregman divergence between vectors `x` and `y`
"""
function bregman(x, y)
    function distance(x, y)
        if x == 0
            return y
        else
            return x * log(x) - (x - y) - x * log(y)
        end
    end
    return sum(distance.(x, y))
end

Main.var"ex-index".bregman

R"""
compute_divergence_min <- function( x, centers ) {
    n = nrow(x)
    k = ncol(centers)
    divergence_min = rep(Inf,n)
    cluster = rep(0,n)
    for(i in 1:k){
      divergence = apply(x,1,bregman, y = centers[,i])
      improvement = (divergence < divergence_min)
      divergence_min[improvement] = divergence[improvement]
      cluster[improvement] = i
    }
    return(list(divergence_min, cluster))
}
"""
function compute_divergence_min( x, centers )
    n = size(x, 1)
    k = size(centers, 2)
    divergence = zeros(n)
    divergence_min = zeros(n)
    cluster = zeros(Int, k)
    for i = 1:k
       for j = 1:n
           divergence[j] = bregman(x[j,:], centers[:,i])
           if divergence[j] < divergence_min[j]
               divergence_min[j] = divergence[j]
               cluster[j] = i
           end
       end
    end
    return divergence_min, cluster
end

compute_divergence_min (generic function with 1 method)

using Random
centers = transpose(o.points[first(randperm(n),k), :])
x = o.points
@rput x
@rput centers
R"""
system.time({compute_divergence_min(x, centers)})
"""

RCall.RObject{RCall.RealSxp}
   user  system elapsed 
  0.068   0.005   0.073

@time divergence_min, cluster = compute_divergence_min( x, centers )

  0.107422 seconds (261.76 k allocations: 13.196 MiB, 99.32% compilation time)

@time divergence_min, cluster = compute_divergence_min( x, centers )

  0.000450 seconds (6.01 k allocations: 477.062 KiB)

function compute_divergence_min!(divergence_min, cluster, x, centers )
    n = size(x, 1)
    k = size(centers, 2)
    divergence = zeros(n)
    for i = 1:k
       for j = 1:n
           divergence[j] = bregman(x[j,:], centers[:,i])
           if divergence[j] < divergence_min[j]
               divergence_min[j] = divergence[j]
               cluster[j] = i
           end
       end
    end
end

compute_divergence_min! (generic function with 1 method)

divergence_min = zeros(n)
cluster = zeros(Int, k)
@time compute_divergence_min!( divergence_min, cluster, x, centers )

  0.038574 seconds (37.91 k allocations: 2.074 MiB, 98.96% compilation time)

@time compute_divergence_min!( divergence_min, cluster, x, centers )

  0.000303 seconds (6.00 k allocations: 472.812 KiB)

Application to authors texts clustering

"textes.txt"

using DelimitedFiles
table = readdlm("textes.txt")

51×210 Matrix{Any}:
 "Mark Twain"     "Mark Twain"  …     "Obama"     "Obama"     ""
 "be"          435                 274         276         283
 "have"        112                 196         191         277
 "say"         153                  10          16          15
 "do"          136                  84          73         106
 "go"          108              …   14          13          26
 "come"         26                  16          20          17
 "see"          42                   7          15          16
 "know"         35                  21          29          23
 "make"         48                  58          43          68
 ⋮                              ⋱                          
 "young"         0                   7          12           9
 "night"        11                   0           3           0
 "cry"           9                   1           0           0
 "word"          7              …    1           5           4
 "mind"          7                   0           4           3
 "house"         8                   8           4           0
 "sit"           1                   2           1           3
 "sir"           3                   0           0           0
 "last"         10              …   18          15          12

Data from texts are stored in some variable df. The commands used for displaying data are the following.

import DataFrames: DataFrame, rename!, ByRow, transform!
import MultivariateStats
import MultivariateStats: PCA, MulticlassLDA
import CategoricalArrays: recode

table = readdlm(joinpath(@__DIR__, "textes.txt"))
df = DataFrame(
    hcat(table[2:end, 1], table[2:end, 2:end]),
    vec(vcat("authors", table[1, 1:end-1])),
    makeunique = true,
)