---
title: Introduction to Julia
subtitle: Biostat/Biomath M257
author: Dr. Hua Zhou @ UCLA
date: today
format:
  html:
    theme: cosmo
    embed-resources: true
    number-sections: true
    toc: true
    toc-depth: 4
    toc-location: left
    code-fold: false
jupyter:
  jupytext:
    formats: 'ipynb,qmd'
    text_representation:
      extension: .qmd
      format_name: quarto
      format_version: '1.0'
      jupytext_version: 1.14.5
  kernelspec:
    display_name: Julia 1.8.5
    language: julia
    name: julia-1.8
---

<img src="./julia_logo.png" align="center" width="400"/>

System information (for reproducibility)

```{julia}
versioninfo()
```

Load packages:

```{julia}
using Pkg
Pkg.activate(pwd())
Pkg.instantiate()
Pkg.status()
```

```{julia}
using BenchmarkTools, Distributions, RCall
using LinearAlgebra, Profile, SparseArrays
```

## What's Julia?

> Julia is a high-level, high-performance dynamic programming language for technical computing, with syntax that is familiar to users of other technical computing environments

* History:
  - Project started in 2009. First public release in 2012 
  - Creators: Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral Shah
  - First major release v1.0 was released on Aug 8, 2018
  - Current stable release v1.8.5 (as of Apr 6, 2023)

* Aim to solve the notorious **two language problem**: Prototype code goes into high-level languages like R/Python, production code goes into low-level language like C/C++. 

    Julia aims to:
> Walks like Python. Runs like C.

<img src="./julia_vs_otherlang.png" align="center" width="600"/>

See <https://julialang.org/benchmarks/> for the details of benchmark.

* Write high-level, abstract code that closely resembles mathematical formulas
    - yet produces fast, low-level machine code that has traditionally only been generated by static languages.

* Julia is more than just "Fast R" or "Fast Matlab"
    - Performance comes from features that work well together.  
    - You can't just take the magic dust that makes Julia fast and sprinkle it on [language of choice]

## Learning resources

1. The (free) online course [Introduction to Julia](https://juliaacademy.com/p/intro-to-julia), by Jane Herriman. 

2. Cheat sheet: [The Fast Track to Julia](https://juliadocs.github.io/Julia-Cheat-Sheet/).  

3. Browse the Julia [documentation](https://docs.julialang.org/en).  
bm
4. For Matlab users, read [Noteworthy Differences From Matlab](https://docs.julialang.org/en/v1/manual/noteworthy-differences/#Noteworthy-differences-from-MATLAB-1).  

    For R users, read [Noteworthy Differences From R](https://docs.julialang.org/en/v1/manual/noteworthy-differences/#Noteworthy-differences-from-R-1).  

    For Python users, read [Noteworthy Differences From Python](https://docs.julialang.org/en/v1/manual/noteworthy-differences/?highlight=matlab#Noteworthy-differences-from-Python-1).  

5. The [Learning page](http://julialang.org/learning/) on Julia's website has pointers to many other learning resources.  

## Julia REPL (Read-Evaluation-Print-Loop)

The `Julia` REPL, or `Julia` shell, has at least five modes.

1. **Default mode** is the Julian prompt `julia>`. Type backspace in other modes to enter default mode.    

2. **Help mode** `help?>`. Type `?` to enter help mode. `?search_term` does a fuzzy search for `search_term`.  

3. **Shell mode** `shell>`. Type `;` to enter shell mode.  

4. **Package mode** `(@v1.8) pkg>`. Type `]` to enter package mode for managing Julia packages (install, uninstall, update, ...).

5. **Search mode** `(reverse-i-search)`. Press `ctrl+R` to enter search model. 

6. With `RCall.jl` package installed, we can enter the **R mode** by typing `$` (shift+4) at Julia REPL.

Some survival commands in Julia REPL:  
1. `exit()` or `Ctrl+D`: exit Julia.

2. `Ctrl+C`: interrupt execution.

3. `Ctrl+L`: clear screen.

0. Append `;` (semi-colon) to suppress displaying output from a command. Same as Matlab.

0. `include("filename.jl")` to source a Julia code file.

## Seek help

* Online help from REPL: `?function_name`.

* Google.

* Julia documentation: <https://docs.julialang.org/en/>.

* Look up source code: `@edit sin(π)`.

* Discourse: <https://discourse.julialang.org>.

* Friends.

## Which IDE?

* Julia homepage lists many choices: Juno, VS Code, Vim, ...

* Unfortunately at the moment there are no mature RStudio- or Matlab-like IDE for Julia yet.

* For dynamic document, e.g., homework, I recommend [Jupyter Notebook](https://jupyter.org/install.html) or [JupyterLab](http://jupyterlab.readthedocs.io/en/stable/). 

* For extensive Julia coding, myself has been happily using the editor [VS Code](https://code.visualstudio.com) with extensions `Julia` and `VS Code Jupyter Notebook Previewer` installed. 

## Julia package system

* Each Julia package is a Git repository. Each Julia package name ends with `.jl`. E.g., `Distributions.jl` package lives at <https://github.com/JuliaStats/Distributions.jl>.   
Google search with `PackageName.jl` usually leads to the package on github.com. 

* The package ecosystem is rapidly maturing; a complete list of **registered** packages (which are required to have a certain level of testing and documentation) is at [http://pkg.julialang.org/](http://pkg.julialang.org/).

* For example, the package called `Distributions.jl` is added with
```julia
# in Pkg mode
(@v1.8) pkg> add Distributions
```
and "removed" (although not completely deleted) with
```julia
# in Pkg mode
(@v1.8) pkg> rm Distributions
```
* The package manager provides a dependency solver that determines which packages are actually required to be installed.

* **Non-registered** packages are added by cloning the relevant Git repository. E.g.,
```julia
# in Pkg mode
(@v1.8) pkg> add https://github.com/OpenMendel/SnpArrays.jl
```

* A package needs only be added once, at which point it is downloaded into your local `.julia/packages` directory in your home directory. 

```{julia}
readdir(Sys.islinux() ? ENV["JULIA_PATH"] * "/pkg/packages" : ENV["HOME"] * "/.julia/packages")
```

* Directory of a specific package can be queried by `pathof()`:

```{julia}
pathof(Distributions)
```

* If you start having problems with packages that seem to be unsolvable, you may try just deleting your .julia directory and reinstalling all your packages. 

* Periodically, one should run `update` in Pkg mode, which checks for, downloads and installs updated versions of all the packages you currently have installed.

* `status` lists the status of all installed packages.

* Using functions in package.
```julia
using Distributions
```
This pulls all of the *exported* functions in the module into your local namespace, as you can check using the `whos()` command. An alternative is
```julia
import Distributions
```
Now, the functions from the Distributions package are available only using 
```julia
Distributions.<FUNNAME>
```
All functions, not only exported functions, are always available like this.

## Calling R from Julia

* The [`RCall.jl`](https://github.com/JuliaInterop/RCall.jl) package allows us to embed R code inside of Julia.

* There are also `PyCall.jl`, `MATLAB.jl`, `JavaCall.jl`, `CxxWrap.jl` packages for interfacing with other languages.

```{julia}
# x is in Julia workspace
x = randn(1000)

# $ is the interpolation operator
R"""
hist($x, main = "I'm plotting a Julia vector")
""";
```

```{julia}
R"""
library(ggplot2)

qplot($x)
"""
```

```{julia}
x = R"""
rnorm(10)
"""
```

```{julia}
# collect R variable into Julia workspace
y = collect(x)
```

* Access Julia variables in R REPL mode:
```julia
julia> x = rand(5) # Julia variable
R> y <- $x
```

* Pass Julia expression in R REPL mode:
```julia
R> y <- $(rand(5))
```

* Put Julia variable into R environment:
```julia
julia> @rput x
R> x
```

* Get R variable into Julia environment:
```julia
R> r <- 2
Julia> @rget r
```

* If you want to call Julia within R, check out the [`JuliaCall`](https://non-contradiction.github.io/JuliaCall//index.html) package.

## Some basic Julia code

```{julia}
# an integer, same as int in R
y = 1
```

```{julia}
# query type of a Julia object
typeof(y)
```

```{julia}
# a Float64 number, same as double in R
y = 1.0
```

```{julia}
typeof(y) 
```

```{julia}
# Greek letters:  `\pi<tab>`
π
```

```{julia}
typeof(π)
```

```{julia}
# Greek letters:  `\theta<tab>`
θ = y + π
```

```{julia}
# emoji! `\:kissing_cat:<tab>`
😽 = 5.0
😽 + 1
```

```{julia}
# `\alpha<tab>\hat<tab>`
α̂ = π
```

For a list of unicode symbols that can be tab-completed, see <https://docs.julialang.org/en/v1/manual/unicode-input/>. Or in the help mode, type `?` followed by the unicode symbol you want to input.

```{julia}
# vector of Float64 0s
x = zeros(5)
```

```{julia}
# vector Int64 0s
x = zeros(Int, 5)
```

```{julia}
# matrix of Float64 0s
x = zeros(5, 3)
```

```{julia}
# matrix of Float64 1s
x = ones(5, 3)
```

```{julia}
# define array without initialization
x = Matrix{Float64}(undef, 5, 3)
```

```{julia}
# fill a matrix by 0s
fill!(x, 0)
```

```{julia}
# initialize an array to be constant 2.5
fill(2.5, (5, 3))
```

```{julia}
# rational number
a = 3//5
```

```{julia}
typeof(a)
```

```{julia}
b = 3//7
```

```{julia}
a + b
```

```{julia}
# uniform [0, 1) random numbers
x = rand(5, 3)
```

```{julia}
# uniform random numbers (in single precision)
x = rand(Float16, 5, 3)
```

```{julia}
# random numbers from {1,...,5}
x = rand(1:5, 5, 3)
```

```{julia}
# standard normal random numbers
x = randn(5, 3)
```

```{julia}
# range
1:10
```

```{julia}
typeof(1:10)
```

```{julia}
1:2:10
```

```{julia}
typeof(1:2:10)
```

```{julia}
# integers 1-10
x = collect(1:10)
```

```{julia}
# or equivalently
[1:10...]
```

```{julia}
# Float64 numbers 1-10
x = collect(1.0:10)
```

```{julia}
# convert to a specific type
convert(Vector{Float64}, 1:10)
```

## Matrices and vectors

### Dimensions

```{julia}
x = randn(5, 3)
```

```{julia}
size(x)
```

```{julia}
size(x, 1) # nrow() in R
```

```{julia}
size(x, 2) # ncol() in R
```

```{julia}
# total number of elements
length(x)
```

### Indexing

```{julia}
# 5 × 5 matrix of random Normal(0, 1)
x = randn(5, 5)
```

```{julia}
# first column
x[:, 1]
```

```{julia}
# first row
x[1, :]
```

```{julia}
# sub-array
x[1:2, 2:3]
```

```{julia}
# getting a subset of a matrix creates a copy, but you can also create "views"
z = view(x, 1:2, 2:3)
```

```{julia}
# same as
@views z = x[1:2, 2:3]
```

```{julia}
# change in z (view) changes x as well
z[2, 2] = 0.0
x
```

```{julia}
# y points to same data as x
y = x
```

```{julia}
# x and y point to same data
pointer(x), pointer(y)
```

```{julia}
# changing y also changes x
y[:, 1] .= 0
x
```

```{julia}
# create a new copy of data
z = copy(x)
```

```{julia}
pointer(x), pointer(z)
```

### Concatenate matrices

```{julia}
# 3-by-1 vector
[1, 2, 3]
```

```{julia}
# 1-by-3 array
[1 2 3]
```

```{julia}
# multiple assignment by tuple
x, y, z = randn(5, 3), randn(5, 2), randn(3, 5)
```

```{julia}
[x y] # 5-by-5 matrix
```

```{julia}
[x y; z] # 8-by-5 matrix
```

### Dot operation (broadcasting)

Dot operation in Julia is elementwise operation, similar to Matlab.

```{julia}
x = randn(5, 3)
```

```{julia}
y = ones(5, 3)
```

```{julia}
x .* y # same x * y in R
```

```{julia}
x .^ (-2)
```

```{julia}
sin.(x)
```

### Basic linear algebra

```{julia}
x = randn(5)
```

```{julia}
# vector L2 norm
norm(x)
```

```{julia}
# same as
sqrt(sum(abs2, x))
```

```{julia}
y = randn(5) # another vector
# dot product
dot(x, y) # x' * y
```

```{julia}
# same as
x'y
```

```{julia}
x, y = randn(5, 3), randn(3, 2)
# matrix multiplication, same as %*% in R
x * y
```

```{julia}
x = randn(3, 3)
```

```{julia}
# conjugate transpose
x'
```

```{julia}
b = rand(3)
x'b # same as x' * b
```

```{julia}
# trace
tr(x)
```

```{julia}
det(x)
```

```{julia}
rank(x)
```

### Sparse matrices

```{julia}
# 10-by-10 sparse matrix with sparsity 0.1
X = sprandn(10, 10, .1)
```

```{julia}
# dump() in Julia is like str() in R
dump(X)
```

```{julia}
# convert to dense matrix; be cautious when dealing with big data
Xfull = convert(Matrix{Float64}, X)
```

```{julia}
# convert a dense matrix to sparse matrix
sparse(Xfull)
```

```{julia}
# syntax for sparse linear algebra is same as dense linear algebra
β = ones(10)
X * β
```

```{julia}
# many functions apply to sparse matrices as well
sum(X)
```

## Control flow and loops

* if-elseif-else-end

```julia
if condition1
    # do something
elseif condition2
    # do something
else
    # do something
end
```

* `for` loop

```julia
for i in 1:10
    println(i)
end
```

* Nested `for` loop:

```julia
for i in 1:10
    for j in 1:5
        println(i * j)
    end
end
```
Same as

```julia
for i in 1:10, j in 1:5
    println(i * j)
end
```

* Exit loop:

```julia
for i in 1:10
    # do something
    if condition1
        break # skip remaining loop
    end
end
```

* Exit iteration:  

```julia
for i in 1:10
    # do something
    if condition1
        continue # skip to next iteration
    end
    # do something
end
```

## Functions 

* In Julia, all arguments to functions are **passed by reference**, in contrast to R and Matlab.

* Function names ending with `!` indicates that function mutates at least one argument, typically the first.
```julia
sort!(x) # vs sort(x)
```

* Function definition

```julia
function func(req1, req2; key1=dflt1, key2=dflt2)
    # do stuff
    return out1, out2, out3
end
```
**Required arguments** are separated with a comma and use the positional notation.  
**Optional arguments** need a default value in the signature.  
**Semicolon** is not required in function call.  
**return** statement is optional.  
Multiple outputs can be returned as a **tuple**, e.g., `return out1, out2, out3`.  

* Anonymous functions, e.g., `x -> x^2`, is commonly used in collection function or list comprehensions.
```julia
map(x -> x^2, y) # square each element in x
```

* Functions can be nested:

```julia
function outerfunction()
    # do some outer stuff
    function innerfunction()
        # do inner stuff
        # can access prior outer definitions
    end
    # do more outer stuff
end
```

* Functions can be vectorized using the Dot syntax:

```{julia}
# defined for scalar
function myfunc(x)
    return sin(x^2)
end

x = randn(5, 3)
myfunc.(x)
```

* **Collection function** (think this as the series of `apply` functions in R).

    Apply a function to each element of a collection:

```julia
map(f, coll) # or
map(coll) do elem
    # do stuff with elem
    # must contain return
end
```

```{julia}
map(x -> sin(x^2), x)
```

```{julia}
map(x) do elem
    elem = elem^2
    return sin(elem)
end
```

```{julia}
# Mapreduce
mapreduce(x -> sin(x^2), +, x)
```

```{julia}
# same as
sum(x -> sin(x^2), x)
```

* List **comprehension**

```{julia}
[sin(2i + j) for i in 1:5, j in 1:3] # similar to Python
```

## Type system

* Every variable in Julia has a type.

* When thinking about types, think about sets.

* Everything is a subtype of the abstract type `Any`.

* An abstract type defines a set of types
    - Consider types in Julia that are a `Number`:

<img src="tree.png" width="800" align="center"/>

* We can explore type hierarchy with `typeof()`, `supertype()`, and `subtypes()`.

```{julia}
# 1.0: double precision, 1: 64-bit integer
typeof(1.0), typeof(1)
```

```{julia}
supertype(Float64)
```

```{julia}
subtypes(AbstractFloat)
```

```{julia}
# Is Float64 a subtype of AbstractFloat?
Float64 <: AbstractFloat
```

```{julia}
# On 64bit machine, Int == Int64
Int == Int64
```

```{julia}
# convert to Float64
convert(Float64, 1)
```

```{julia}
# same as casting
Float64(1)
```

```{julia}
# Float32 vector
x = randn(Float32, 5)
```

```{julia}
# convert to Float64
convert(Vector{Float64}, x)
```

```{julia}
# same as broadcasting (dot operatation)
Float64.(x)
```

```{julia}
# convert Float64 to Int64
convert(Int, 1.0)
```

```{julia}
convert(Int, 1.5) # should use round(1.5)
```

```{julia}
round(Int, 1.5)
```

## Multiple dispatch

* Multiple dispatch lies in the core of Julia design. It allows built-in and user-defined functions to be overloaded for different combinations of argument types.

* Let's consider a simple "doubling" function:

```{julia}
g(x) = x + x
```

```{julia}
g(1.5)
```

This definition is too broad, since some things, e.g., strings, can't be added 

```{julia}
g("hello world")
```

* This definition is correct but too restrictive, since any `Number` can be added.

```{julia}
g(x::Float64) = x + x
```

* This definition will automatically work on the entire type tree above!

```{julia}
g(x::Number) = x + x
```

This is a lot nicer than 
```julia
function g(x)
    if isa(x, Number)
        return x + x
    else
        throw(ArgumentError("x should be a number"))
    end
end
```

* `methods(func)` function display all methods defined for `func`.

```{julia}
methods(g)
```

* When calling a function with multiple definitions, Julia will search from the narrowest signature to the broadest signature.

* `@which func(x)` marco tells which method is being used for argument signature `x`.

```{julia}
# an Int64 input
@which g(1)
```

```{julia}
# a Vector{Float64} input
@which g(randn(5))
```

## Just-in-time compilation (JIT)

Following figures are taken from Arch D. Robinson's slides [Introduction to Writing High Performance Julia](https://www.youtube.com/watch?v=szE4txAD8mk).


| <img src="./julia_toolchain.png" alt="Julia toolchain" style="width: 400px;"/> | <img src="./julia_introspect.png" alt="Julia toolchain" style="width: 500px;"/> |
|----------------------------------|------------------------------------|
|||

* `Julia`'s efficiency results from its capability to infer the types of **all** variables within a function and then call LLVM compiler to generate optimized machine code at run-time. 

Consider the `g` (doubling) function defined earlier. This function will work on **any** type which has a method for `+`.

```{julia}
g(2), g(2.0)
```

**Step 1**: Parse Julia code into [abstract syntax tree (AST)](https://en.wikipedia.org/wiki/Abstract_syntax_tree).

```{julia}
@code_lowered g(2)
```

**Step 2**: Type inference according to input type.

```{julia}
# type inference for integer input
@code_warntype g(2)
```

```{julia}
# type inference for Float64 input
@code_warntype g(2.0)
```

**Step 3**: Compile into **LLVM bitcode** (equivalent of R bytecode generated by the `compiler` package).

```{julia}
# LLVM bitcode for integer input
@code_llvm g(2)
```

```{julia}
# LLVM bitcode for Float64 input
@code_llvm g(2.0)
```

We didn't provide a type annotation. But different LLVM bitcodes were generated depending on the argument type!

In R or Python, `g(2)` and `g(2.0)` would use the same code for both.

In Julia, `g(2)` and `g(2.0)` dispatches to optimized code for `Int64` and `Float64`, respectively.

For integer input `x`, LLVM compiler is smart enough to know `x + x` is simply shifting `x` by 1 bit, which is faster than addition.
 
* **Step 4**: Lowest level is the **assembly code**, which is machine dependent.

```{julia}
# Assembly code for integer input
@code_native g(2)
```

```{julia}
# Assembly code for Float64 input
@code_native g(2.0)
```

## Profiling Julia code

Julia has several built-in tools for profiling. The `@time` marco outputs run time and heap allocation. Note the first call of a function incurs (substantial) compilation time.

```{julia}
# a function defined earlier
function tally(x::Array)
    s = zero(eltype(x))
    for v in x
        s += v
    end
    s
end

using Random
Random.seed!(257)
a = rand(20_000_000)
@time tally(a) # first run: include compile time
```

```{julia}
@time tally(a)
```

For more robust benchmarking, the [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) package is highly recommended.

```{julia}
@benchmark tally($a)
```

The `Profile` module gives line by line profile results.

```{julia}
Profile.clear()
@profile tally(a)
Profile.print(format=:flat)
```

One can use [`ProfileView`](https://github.com/timholy/ProfileView.jl) package for better visualization of profile data:

```julia
using ProfileView

ProfileView.view()
```

```{julia}
# check type stability
@code_warntype tally(a)
```

```{julia}
# check LLVM bitcode
@code_llvm tally(a)
```

```{julia}
@code_native tally(a)
```

**Exercise:** Annotate the loop in `tally` function by `@simd` and look for the difference in LLVM bitcode and machine code.

## Memory profiling

Detailed memory profiling requires a detour. First let's write a script `bar.jl`, which contains the workload function `tally` and a wrapper for profiling.

```{julia}
;cat bar.jl
```

Next, in terminal, we run the script with `--track-allocation=user` option.

```{julia}
#;julia --track-allocation=user bar.jl
```

The profiler outputs a file `bar.jl.51116.mem`.

```{julia}
;cat bar.jl.51116.mem
```
