index.qmd

---
title: "Why Julia?"
subtitle: "A gentle pitch"
author: "Jose Storopoli, PhD"
format:
  revealjs: 
    slide-number: true
    transition: slide
    chalkboard: 
      buttons: false
    preview-links: auto
    footer: <https://storopoli.github.io/Why-Julia>
    logo: images/julia-dots.svg
    callout-appearance: minimal
execute:
  echo: true
  cache: true
---

```{julia}
#| echo: false
#| output: false
ENV["PYCALL_JL_RUNTIME_PYTHON"] = Sys.which("python")
using PyCall
```

## Agenda

</br> </br>

::: incremental

1. speed
2. ease-of-use
3. composability
:::

## What I Assume?

</br> </br>

::: incremental

- Python background
- scientific computing background
:::

## So let's dive in?

![](images/julia_python_meme.jpg){fig-align="center"}

## Julia is past beyond "experimental" {.smaller .scrollable}

:::incremental

- NASA uses Julia in a supercomputer to analyze the
   ["Largest Batch of Earth-Sized Planets Ever Found"](https://exoplanets.nasa.gov/news/1669/seven-rocky-trappist-1-planets-may-be-made-of-similar-stuff/) and achieve a whopping **1,000x speedup** to catalog 188 million astronomical objects in 15 minutes.
- [The Climate Modeling Alliance (CliMa)](https://clima.caltech.edu/)
   is using mostly Julia to **model climate in the GPU and CPU**.
Launched in 2018 in collaboration with researchers at Caltech, the NASA Jet Propulsion Laboratory, and the Naval Postgraduate School, CliMA is utilizing recent progress in computational science to develop an Earth system model that can predict droughts, heat waves, and rainfall with unprecedented precision and speed.
- [US Federal Aviation Administration (FAA) is developing an **Airborne Collision Avoidance System (ACAS-X)** using Julia](https://youtu.be/19zm1Fn0S9M).
This is a nice example of the "Two-Language Problem".
Previous solutions used Matlab to develop the algorithms and C++ for a fast implementation.
Now, FAA is using one language to do all this: Julia.
- [**175x speedup** for Pfizer's pharmacology models using GPUs in Julia](https://juliahub.com/case-studies/pfizer/).
It was presented as a [poster](https://chrisrackauckas.com/assets/Posters/ACoP11_Poster_Abstracts_2020.pdf) in the 11th American Conference of Pharmacometrics (ACoP11) and [won a quality award](https://web.archive.org/web/20210121164011/https://www.go-acop.org/abstract-awards).
- [The Attitude and Orbit Control Subsystem (AOCS) of the Brazilian satellite Amazonia-1 is **written 100% in Julia**](https://discourse.julialang.org/t/julia-and-the-satellite-amazonia-1/57541) by Ronan Arraes Jardim Chagas (<https://ronanarraes.com/>).
- [Brazil's national development bank (BNDES) ditched a paid solution and opted for open-source Julia modeling and gained a **10x speedup**.](https://youtu.be/NY0HcGqHj3g)
:::

::: footer
If this is not enough, there are more case studies in [JuliaHub website](https://juliahub.com/case-studies/).
:::

## Speed

</br>

**Julia is fast!**

</br>

::: {.fragment .fade-in}
Two examples:
:::

::: incremental

- Data Wrangling: `pandas` versus `DataFrames.jl`
- ODE solving: `scipy` versus `DifferentialEquations.jl`
:::

## Benchmarking --- Data Wrangling

Common data wrangling scenario doing "split-apply-combine" operations.

::: incremental

- 10,000 observations
- 1 categorical variable `x` $\in \{\mathrm{A}, \mathrm{B}, \mathrm{C}, \mathrm{D}\}$
- 2 continuous variables:
  - `y` $\in [0, 1]$
  - `z` $\text{Normal}(0, 1)$
:::

## Benchmarking --- Data Wrangling (Python)

```{julia}
using BenchmarkTools
py"""
import pandas as pd
import numpy as np

n = 10000

df = pd.DataFrame({'x': np.random.choice(['A', 'B', 'C', 'D'], n, replace=True),
                   'y': np.random.randn(n),
                   'z': np.random.rand(n)})
"""
@btime py"df.groupby('x').agg({'y': 'median', 'z': 'mean'})";
```

## Benchmarking --- Data Wrangling (Julia)

```{julia}
using Random
using DataFrames
using BenchmarkTools
using Chain
Random.seed!(123)

n = 10_000

df = DataFrame(
    x=rand('A':'D', n),
    y=rand(n),
    z=randn(n),
)

@btime @chain $df begin
    groupby(:x)
    combine(:y => median, :z => mean)
end;
```

## Benchmarking --- ODE Solver

Second order non-linear ODE example with a **simple pendulum**

:::: {.columns}

::: {.column width="50%"}

$$
\begin{align*}
&\dot{\theta} = d{\theta} \\
&\dot{d\theta} = - \frac{g}{L}{\sin(\theta)}
\end{align*}
$$

:::

::: {.column width="50%"}

![](images/pendulum.png){width=60%}
:::

::::

## Benchmarking --- ODE Solver (Julia)

```{julia}
using DifferentialEquations

# Constants
const g = 9.81
L = 1.0

# Initial Conditions
u₀ = [0, π/2]
tspan = (0.0, 6.3)

# Define the problem
function simplependulum(du, u, p, t)
    θ, dθ = u
    du[1] = dθ
    du[2] = -(g/L)*sin(θ)
end

# Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
# RK 4/5th order solver (Tsitouras)
@btime solve(prob, Tsit5(); saveat=range(tspan...; length=1_000));
```

## Benchmarking --- ODE Solver (Python)

```{julia}
py"""
import numpy as np
from scipy.integrate import odeint

# Constants
g = 9.81
L = 1.0

# Initial Conditions
u0 = [0, np.pi/2]
tspan = np.linspace(0.0, 6.3, 1000)

def simplependulum(u, t, g, L):
    theta, dtheta = u
    dydt = [dtheta, -(g/L)*np.sin(theta)]
    return dydt
"""

# RK 4/5th order solver (Dormand-Prince)
@btime py"odeint(simplependulum, u0, tspan, args=(g, L))";
```

## Why Julia is so Fast?

![LLVM](images/LLVM_logo.png)

::: incremental

- *just-in-time* compilation for the LLVM compiler
- exposes everything in *intermediate representation* code
- then LLVM does what does best: **OPTIMIZE**
- including `for`-loops
:::

## Why Julia is so Fast? --- LLVM code {.scrollable}

```{julia}
#| output-location: slide
using Statistics: mean
@code_llvm mean(1:10)
```

::: footer
output in next slide
:::

## Ease of Use

The syntax is quite similar to Python.

. . .

But, no indentation and every keyword needs an `end`.

. . .

:::: {.columns}

::: {.column width="50%"}

Julia:

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

:::

::: {.column width="50%"}

Python:

```python
for i in range(10):
    print(i)
```

:::

::::

## It's Julia all the way down

If you need to find something just use the
`@which` macro on a type or a function signature.

</br>

. . .

```{julia}
@which DataFrame # type
```

. . .

</br>

```{julia}
@which mean(1:10) # function signature.
```

## Composability

::: incremental

- It is very easy to create new packages that have types and functions.

- You can extend other package's functions,
  including Julia's `Base` standard library to you new types.

- And you can also create new functions for other Package's types.
:::

## Composability --- Example with `Point`

```{julia}
#| output: false
struct Point
  x::Float64
  y::Float64
end

function Base.:+(x::Point, y::Point)
  return Point(x.x + y.x, x.y + y.y)
end

function distance(x::Point, y::Point)
  return sqrt( (x.x - y.x)^2 + (x.y - y.y)^2 )
end

p1 = Point(1, 1); p2 = Point(2, 2)
```

. . .

```{julia}
p1 + p2
```

. . .

```{julia}
distance(p1, p2)
```

## Composability --- Example with Autodiff {.smaller}

Suppose you are creating a new sort of graph structure that allows for differentiation and integration, i.e you can take gradients, Jacobians, Hessians and so on.

. . .

</br>

Imagine having to code the whole API in `libtorch` (`PyTorch` C++ backend). Including:

. . .

::: incremental

- types
- constructors
- linear algebra functions
- autodiff rules
:::

. . .

And in the end you can *only* use PyTorch.
You would have to do the whole thing again for JAX or any other autodiff backend.

## Composability --- Example with Autodiff (Julia) {.smaller}

Now let's see how we do this in Julia?

::: incremental

- We can create a package `DifferentialGraph.jl`.
- Add [`ChainRulesCore.jl`](https://github.com/JuliaDiff/ChainRulesCore.jl) as a dependency.
- Create forward- and reverse-mode derivative rules: `rrules` or `frules`
:::

. . .

Now we can use you differential graphs with all of these backends:

::: incremental

- [`ForwardDiff.jl`](https://github.com/JuliaDiff/ForwardDiff.jl): forward-mode AD
- [`ReverseDiff.jl`](https://github.com/JuliaDiff/ReverseDiff.jl): tape-based reverse-mode AD
- [`Zygote.jl`](https://github.com/FluxML/Zygote.jl): source-to-source reverse-mode AD
- [`Enzyme.jl`](https://github.com/EnzymeAD/Enzyme.jl): Julia bindings for Enzyme which ADs LLVM (low-level)
- [`Diffractor.jl`](https://github.com/JuliaDiff/Diffractor.jl): experimental mixed-mode AD meant to replace Zygote.jl
:::

::: footer
Since your graph has derivatives you can use gradient-based solvers to perform optimization.
:::

## Composability --- Examples from the Julia Ecosystem

::: incremental

- Bayesian Neural Nets: [`Flux.jl`](https://fluxml.ai/) neural network inside a [`Turing.jl`](https://turing.ml) model.
- Bayesian COVID modeling: [`DifferentialEquations.jl`](https://github.com/SciML/DifferentialEquations.jl) ODE inside a [`Turing.jl`](https://turing.ml) model.
- Quaternion ODE solver in the GPU: [`Quaternions.jl`](https://github.com/JuliaGeometry/Quaternions.jl) types in a [`DifferentialEquations.jl`](https://github.com/SciML/DifferentialEquations.jl)
   ODE running in `CuArrays` from [`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl).
:::

## My Pitch {.scrollable}

::: incremental

- It is fast.
- It is easy to use.
- Learning the basics of Julia will make your life so much easier
   in all other packages.
   You don't need to learn specific package syntax to be effective in
   using a certain package.
- A bliss to install in Windows, Mac, and Linux (even in clusters).
- Very good community, check the [discourse](https://discourse.julialang.org).
- Very "nerdy", "mathy", and "geeky" userbase.
- If you are creating new stuff, like research or algorithms,
   you don't want to have to stumble upon FORTRAN or C code
   (`scipy`, `numpy`, `pytorch` etc.).
   In Julia everything is in Julia.
- You can easily mix-and-match types and functions from different packages,
   as you saw in the previous slide.
- Good language interop:
  - C: [standard library](https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/)
  - FORTRAN: [standard library](https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/)
  - Python: [`PyCall.jl`](https://github.com/JuliaPy/PyCall.jl)
  - R: [`RCall.jl`](https://github.com/JuliaInterop/RCall.jl)
:::

## It's not all rainbows

::: incremental

- Hard to onboard people. Sometimes they don't want to learn new stuff
   (I mean we still have FORTRAN around ...).
- Not widely used in the marketplace (but tons of academic usage).
- Some package ecosystems are not mature enough, e.g. survival analysis.
   But, differential equations is way more mature than other
   scientific computing languages.
- In my point of view, Julia's strength is in **scientific computing**.
   For all other things, you might not have additional benefits.
:::

## Some nice packages {.scrollable}

::: incremental

- The whole [standard library](https://docs.julialang.org/en/v1/), especially [`LinearAlgebra`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/) module
- [`DifferentialEquations.jl`](https://docs.sciml.ai/DiffEqDocs/latest/), [`NeuralPDE.jl`](https://docs.sciml.ai/NeuralPDE/dev/)
   and the whole [SciML package ecosystem](https://sciml.ai)
- [`Flux.jl`](https://fluxml.ai/) and [`MLJ.jl`](https://alan-turing-institute.github.io/MLJ.jl/dev/)
- [`DataFrames.jl`](https://dataframes.juliadata.org/stable/) and [`DataFramesMeta.jl`](https://juliadata.github.io/DataFramesMeta.jl/stable/)
- [`Makie.jl`](https://makie.org) and [`AlgebraOfGraphics.jl`](https://aog.makie.org)
- [`Turing.jl`](https://turing.ml)
- [`Pluto.jl`](https://plutojl.org)
- [`JuMP`](https://jump.dev)
- [`Distributions.jl`](https://juliastats.org/Distributions.jl/latest/)
:::

::: footer
these are all clickable links
:::

## Conclusions

:::incremental

- Julia is pretty darn awesome.
- Easy to get going, and you can always make it faster by
   just optimizing your Julia code.
- No need to drop down to C++.
- Buuuut it can't beat Python at deep learning.
- Otherwise, it's worth a try.
- Godspeed to you.
:::

## Packages Used

:::: {.columns}

::: {.column width="50%"}

```{julia}
#| echo: false
using Pkg
deps = [pair.second for pair in Pkg.dependencies()]
deps = filter(p -> p.is_direct_dep, deps)
deps = filter(p -> !isnothing(p.version), deps)
list = ["$(p.name) $(p.version)" for p in deps]
sort!(list)
println("Julia: $(VERSION)")
println(join(list, '\n'))
```

:::

::: {.column width="50%"}

```{julia}
#| echo: false
py"""
import sys
import numpy as np
import scipy
import pandas as pd

print(f"Python: {sys.version}")
print(f"numpy: {np.__version__}")
print(f"scipy: {scipy.__version__}")
print(f"pandas: {pd.__version__}")
"""
```

:::
::::

## System Information

```{julia}
#| echo: false
versioninfo()
```