1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# Install via juliaup (the official version manager, like rustup for Rust)
curl -fsSL https://install.julialang.org | sh

# Install a specific version
juliaup add 1.11
juliaup default 1.11

# Verify
julia --version
# julia version 1.11.2

# Start the REPL
julia

1
2
3
4
5
6
7

# In the REPL:
# Press ] to enter Pkg mode
# pkg> add BenchmarkTools DataFrames Plots

# Press Backspace to return to Julia mode
# Press ; to run shell commands
# shell> ls -la

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# Define multiple methods for the same function name
# Julia picks the right one based on all argument types

function combine(a::Int, b::Int)
    println("Combining two integers: $(a + b)")
    return a + b
end

function combine(a::String, b::String)
    println("Combining two strings: $(a * b)")
    return a * b
end

function combine(a::Int, b::String)
    println("Int and String: $(repeat(b, a))")
    return repeat(b, a)
end

function combine(a::String, b::Int)
    println("String and Int: $(repeat(a, b))")
    return repeat(a, b)
end

combine(3, 4)          # → "Combining two integers: 7", returns 7
combine("hello", "!")  # → "Combining two strings: hello!", returns "hello!"
combine(3, "ha")       # → "Int and String: hahaha", returns "hahaha"

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

# Abstract types define the hierarchy
abstract type Material end
abstract type Metal <: Material end
abstract type Plastic <: Material end

# Concrete types
struct Steel <: Metal
    yield_strength::Float64  # MPa
    density::Float64          # kg/m³
end

struct Aluminum <: Metal
    yield_strength::Float64
    density::Float64
end

struct ABS <: Plastic
    tensile_strength::Float64
    density::Float64
end

# Define interactions based on material combinations
# This is where multiple dispatch shines: behavior from the combination

function friction_coefficient(a::Metal, b::Metal)
    # Metal-metal: depends on similar hardness
    return 0.15 + rand() * 0.05   # roughly 0.15–0.20 for lubricated metal
end

function friction_coefficient(a::Metal, b::Plastic)
    # Metal on plastic: lower friction
    return 0.05 + rand() * 0.03
end

function friction_coefficient(a::Plastic, b::Metal)
    # Symmetric by convention
    return friction_coefficient(b, a)
end

function friction_coefficient(a::Plastic, b::Plastic)
    # Plastic on plastic: varies widely
    return 0.2 + rand() * 0.15
end

# Usage
steel = Steel(250.0, 7850.0)
alum  = Aluminum(270.0, 2700.0)
abs_  = ABS(40.0, 1050.0)

μ1 = friction_coefficient(steel, alum)   # dispatches to Metal-Metal method
μ2 = friction_coefficient(steel, abs_)   # dispatches to Metal-Plastic method
μ3 = friction_coefficient(abs_, steel)   # dispatches to Plastic-Metal (calls Metal-Plastic)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

using InteractiveUtils

# @which tells you exactly which method was selected
@which combine(3, 4)
# combine(a::Int64, b::Int64)
#     @ Main REPL[1]:2

@which friction_coefficient(steel, alum)
# friction_coefficient(a::Metal, b::Metal)
#     @ Main REPL[15]:2

# methods() shows all methods defined for a function
methods(friction_coefficient)
# 4 methods for generic function "friction_coefficient" from Main:
#  [1] friction_coefficient(a::Metal, b::Metal)
#  [2] friction_coefficient(a::Metal, b::Plastic)
#  [3] friction_coefficient(a::Plastic, b::Metal)
#  [4] friction_coefficient(a::Plastic, b::Plastic)

# methodswith() shows all methods that accept a given type
methodswith(Steel)

1
2
3
4
5
6
7

class Matrix:
    def __add__(self, other):
        if isinstance(other, Matrix):
            return self._matrix_add(other)
        elif isinstance(other, (int, float)):
            return self._scalar_add(other)
        return NotImplemented

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

# Core library defines these types and basic operations
struct DenseMatrix
    data::Matrix{Float64}
end

struct SparseMatrix
    rows::Vector{Int}
    cols::Vector{Int}
    vals::Vector{Float64}
    m::Int
    n::Int
end

# Multiplication methods — can be defined anywhere, by anyone
import Base: *

function *(A::DenseMatrix, B::DenseMatrix)
    println("Dense × Dense: using BLAS DGEMM")
    DenseMatrix(A.data * B.data)
end

function *(A::SparseMatrix, B::DenseMatrix)
    println("Sparse × Dense: CSR-Dense multiply")
    # ... sparse matrix multiplication logic
end

function *(A::DenseMatrix, B::SparseMatrix)
    println("Dense × Sparse: transpose trick")
    # ...
end

function *(A::SparseMatrix, B::SparseMatrix)
    println("Sparse × Sparse: CSR-CSR multiply")
    # ...
end

# New type from a new package — extend without modifying core code
struct GPUMatrix
    device_ptr::UInt64
    m::Int
    n::Int
end

function *(A::GPUMatrix, B::DenseMatrix)
    println("GPU × Dense: cudaGemm")
    # ...
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# Abstract types: define the hierarchy, have no data fields
abstract type Number end
abstract type Real <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer <: Real end

# Concrete types: have data fields, are instantiated
struct Point2D{T<:Real}  # parametric concrete type
    x::T
    y::T
end

# Built-in type hierarchy
# Number
#   ├── Complex{T}
#   └── Real
#         ├── AbstractFloat
#         │     ├── Float16
#         │     ├── Float32
#         │     └── Float64     ← default floating point
#         ├── Integer
#         │     ├── Bool
#         │     ├── Signed
#         │     │     ├── Int8, Int16, Int32, Int64, Int128
#         │     └── Unsigned
#         │           └── UInt8, UInt16, UInt32, UInt64, UInt128
#         └── Rational{T}

# Check type relationships
Int64 <: Integer      # true
Float64 <: Real       # true
Float64 <: Integer    # false
Int64 <: Number       # true (transitivity)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

# :: is a type assertion, not a declaration
x::Int = 5           # x is declared as Int
x = 3.14             # ERROR: TypeError: in assignment, expected Int, got Float64

# In function signatures, :: dispatches on type
function scale(v::Vector{Float64}, s::Float64)
    return v .* s   # .* broadcasts element-wise
end

# Without annotation, function accepts any type
function scale_any(v, s)
    return v .* s
end

# isa check
isa(3, Int)           # true
isa(3, Number)        # true
isa(3.0, Int)         # false
isa(3.0, AbstractFloat)  # true

3 isa Number          # alternative syntax

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

# Parametric types are like generics in Java/Rust but more powerful
struct Stack{T}
    data::Vector{T}
end

function push_item!(s::Stack{T}, item::T) where T
    push!(s.data, item)
    return s
end

function pop_item!(s::Stack{T}) where T
    return pop!(s.data)
end

int_stack = Stack{Int}([])
push_item!(int_stack, 1)
push_item!(int_stack, 2)
push_item!(int_stack, 3)
pop_item!(int_stack)   # returns 3

str_stack = Stack{String}([])
push_item!(str_stack, "hello")
push_item!(str_stack, "world")
# push_item!(str_stack, 42)   # ERROR: MethodError, Int64 vs String

# Constrained parametric types
function dot_product(a::Vector{T}, b::Vector{T}) where T<:Number
    return sum(a[i] * b[i] for i in eachindex(a))
end

# Multiple type parameters
struct Pair{A, B}
    first::A
    second::B
end

p = Pair{String, Int}("age", 42)
p.first   # "age"
p.second  # 42

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# Union: a value can be one of several types
function process(x::Union{Int, Float64, String})
    if x isa Int
        return x * 2
    elseif x isa Float64
        return x ^ 2
    else
        return uppercase(x)
    end
end

process(5)      # 10
process(2.5)    # 6.25
process("hi")   # "HI"

# Union{T, Nothing} is the idiomatic Optional/Maybe
function find_first_positive(v::Vector{Float64})::Union{Float64, Nothing}
    for x in v
        if x > 0
            return x
        end
    end
    return nothing
end

result = find_first_positive([-1.0, -2.0, 3.0, 4.0])
if result !== nothing
    println("Found: $result")
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# This feels dynamic — no type annotations
function sum_squares(v)
    total = 0.0
    for x in v
        total += x^2
    end
    return total
end

# First call: Julia compiles a specialized version for Vector{Float64}
v = rand(1_000_000)
@time sum_squares(v)   # includes compile time
# 0.054321 seconds (includes JIT)

# Second call: hits compiled code
@time sum_squares(v)
# 0.002831 seconds (pure execution)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

function add_scaled(x::Float64, y::Float64, scale::Float64)
    return x + scale * y
end

# @code_lowered: Julia's internal IR before type inference
@code_lowered add_scaled(1.0, 2.0, 3.0)
# CodeInfo
#  1 ─ %1 = scale * y
#  │   %2 = x + %1
#  └──      return %2

# @code_typed: after type inference — every value has a known type
@code_typed add_scaled(1.0, 2.0, 3.0)
# CodeInfo for (::typeof(add_scaled))(x::Float64, y::Float64, scale::Float64)
#  1 ─ %1 = Base.mul_float(scale, y)::Float64
#  │   %2 = Base.add_float(x, %1)::Float64
#  └──      return %2
# => Float64

# @code_llvm: LLVM IR — what gets optimized by LLVM
@code_llvm add_scaled(1.0, 2.0, 3.0)
# ; Function Signature: add_scaled(Float64, Float64, Float64)
# define double @julia_add_scaled(...) {
# top:
#   %0 = fmul double %scale, %y
#   %1 = fadd double %x, %0
#   ret double %1
# }

# @code_native: actual machine code (x86-64 assembly)
@code_native add_scaled(1.0, 2.0, 3.0)
# .text
# ; Function Signature: add_scaled(Float64, Float64, Float64)
#   vmulsd  %xmm2, %xmm1, %xmm1
#   vaddsd  %xmm1, %xmm0, %xmm0
#   retq
# Three instructions. Same as a C compiler would emit.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

# TYPE UNSTABLE — bad for performance
function unstable_mean(v)
    total = 0       # Int! This is wrong
    for x in v
        total += x   # Int + Float64 = Float64, type changes mid-loop
    end
    return total / length(v)  # return type depends on runtime values
end

# TYPE STABLE — fast
function stable_mean(v)
    total = zero(eltype(v))   # starts as 0.0 for Float64 arrays
    for x in v
        total += x   # Float64 + Float64 = Float64, type never changes
    end
    return total / length(v)  # always Float64
end

# @code_warntype highlights type instabilities in red
@code_warntype unstable_mean([1.0, 2.0, 3.0])
# Body::Union{Float64, Int64}   ← RED: unstable return type
#  1 ─       total = 0
#  ...

@code_warntype stable_mean([1.0, 2.0, 3.0])
# Body::Float64   ← GREEN: stable return type

# Performance difference is real
using BenchmarkTools
v = rand(1_000_000)

@btime unstable_mean($v)
# 12.4 ms (3 allocations: 48 bytes)  ← extra allocations from type changes

@btime stable_mean($v)
# 1.1 ms (0 allocations)  ← no allocations, pure numeric computation

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

using BenchmarkTools

# @btime: benchmark with warmup, best of N runs, interpolates constants
function naive_dot(a::Vector{Float64}, b::Vector{Float64})
    s = 0.0
    for i in eachindex(a)
        s += a[i] * b[i]
    end
    return s
end

a = rand(10_000)
b = rand(10_000)

# $ interpolates variables to avoid benchmark overhead from global variable lookup
@btime naive_dot($a, $b)
# 6.2 μs (0 allocations: 0 bytes)

@btime dot($a, $b)      # LinearAlgebra.dot — uses BLAS
# 4.8 μs (0 allocations: 0 bytes)

# @benchmark gives detailed statistics
result = @benchmark naive_dot($a, $b)
# BenchmarkTools.Trial: 10000 samples with 1 evaluation.
#  Range (min … max):  5.7 μs … 22.3 μs
#  Time  (median):     6.1 μs               ← use median, not mean
#  Time  (mean ± σ):   6.3 μs ±  0.9 μs
#  Allocs: 0
#  Bytes: 0

# Compare functions
suite = BenchmarkGroup()
suite["naive"] = @benchmarkable naive_dot($a, $b)
suite["blas"]  = @benchmarkable dot($a, $b)

tune!(suite)
results = run(suite)
judge(minimum(results["naive"]), minimum(results["blas"]))

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

using Profile

function heavy_computation(n::Int)
    result = zeros(n, n)
    for i in 1:n, j in 1:n
        result[i, j] = sin(i * j / n)
    end
    return result
end

# Profile 100 iterations
@profile for _ in 1:100
    heavy_computation(500)
end

# Text profile: call counts per line
Profile.print()

# Flamegraph with ProfileView.jl
using ProfileView
ProfileView.view()  # Opens interactive flamegraph

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

using LinearAlgebra

# Matrix operations
A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 10.0]
b = [1.0, 2.0, 3.0]

# Basic operations
A'           # Adjoint (transpose for real matrices)
tr(A)        # Trace
det(A)       # Determinant
inv(A)       # Inverse (expensive — usually solve instead)
norm(b)      # 2-norm
norm(A, 2)   # Matrix 2-norm (largest singular value)
rank(A)      # Matrix rank

# Solve Ax = b (better than inv(A) * b)
x = A \ b
@show norm(A * x - b)   # should be ~machine epsilon

# Factorizations
F = lu(A)    # LU decomposition
F.L          # Lower triangular
F.U          # Upper triangular
F.p          # Permutation vector

# Solve multiple right-hand sides efficiently
B = rand(3, 5)
X = F \ B    # Uses precomputed factorization

# QR decomposition
Q, R = qr(A)
norm(Q * R - A)   # ≈ 0

# SVD
U, S, V = svd(A)
norm(U * Diagonal(S) * V' - A)   # ≈ 0

# Eigendecomposition
vals, vecs = eigen(A)
# Verify: A * vecs[:, 1] ≈ vals[1] * vecs[:, 1]

# Cholesky (for symmetric positive definite matrices)
Σ = [4.0 2.0; 2.0 3.0]   # covariance matrix
C = cholesky(Σ)
C.L * C.L'   # ≈ Σ

# BLAS operations (direct BLAS calls for maximum performance)
using LinearAlgebra.BLAS

# DGEMM: C = α * A * B + β * C
A_big = rand(1000, 1000)
B_big = rand(1000, 1000)
C_big = zeros(1000, 1000)
BLAS.gemm!('N', 'N', 1.0, A_big, B_big, 0.0, C_big)

# Set number of BLAS threads
BLAS.set_num_threads(8)
BLAS.get_num_threads()

# Specialized matrix types — Julia uses dispatch to pick optimal algorithms
D = Diagonal([1.0, 2.0, 3.0])        # Diagonal matrix
T = Tridiagonal([1,1,1], [4,4,4,4], [1,1,1])  # Tridiagonal
S = Symmetric(A + A')                 # Symmetric (uses symmetric algorithms)
H = Hermitian(A + A')                 # Hermitian (complex matrices)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

using Plots

# Backend selection
# gr()       — GR backend, fast, good for exploration (default)
# pyplot()   — Matplotlib, familiar to Python users
# plotlyjs() — Interactive in browser
# pgfplotsx() — LaTeX-quality for publications

# Basic plotting
x = range(0, 2π, length=200)
y = sin.(x)   # .() broadcasts over the array

plot(x, y, label="sin(x)", linewidth=2, color=:blue)
plot!(x, cos.(x), label="cos(x)", linewidth=2, color=:red, linestyle=:dash)
xlabel!("x")
ylabel!("f(x)")
title!("Trigonometric Functions")
savefig("trig.png")

# Subplots
p1 = plot(x, sin.(x), title="sin(x)", legend=false)
p2 = plot(x, cos.(x), title="cos(x)", legend=false)
p3 = plot(x, sin.(x) .* exp.(-x/4), title="Damped sin", legend=false)
p4 = scatter(sin.(x), cos.(x), title="Unit Circle", aspect_ratio=:equal, legend=false)

plot(p1, p2, p3, p4, layout=(2,2), size=(900, 600))
savefig("subplot_demo.png")

# Contour and surface plots
f(x, y) = sin(x) * cos(y)
xs = range(-π, π, length=100)
ys = range(-π, π, length=100)
z  = [f(x, y) for y in ys, x in xs]

contourf(xs, ys, z, levels=20, color=:viridis)
surface(xs, ys, z, xlabel="x", ylabel="y", zlabel="f(x,y)")

# CairoMakie for publication-quality figures
using CairoMakie

fig = Figure(resolution=(800, 600), fontsize=14)
ax = Axis(fig[1, 1],
    xlabel="Time (s)",
    ylabel="Amplitude",
    title="Signal Processing",
    xgridvisible=false)

t = LinRange(0, 10, 1000)
signal = @. sin(2π * 2 * t) + 0.3 * sin(2π * 5 * t)
noise  = 0.1 * randn(length(t))

lines!(ax, t, signal, color=:steelblue, linewidth=1.5, label="Signal")
lines!(ax, t, signal + noise, color=(:red, 0.4), linewidth=0.8, label="Noisy")
axislegend(ax)

save("publication_figure.pdf", fig)  # Vector PDF, LaTeX-ready

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

using DataFrames, CSV, Statistics

# Create a DataFrame
df = DataFrame(
    id       = 1:1000,
    name     = ["Person_$i" for i in 1:1000],
    age      = rand(18:80, 1000),
    salary   = round.(rand(30000:200000, 1000) ./ 1000) .* 1000,
    dept     = rand(["Engineering", "Sales", "HR", "Finance"], 1000),
    active   = rand(Bool, 1000),
)

# Load from CSV
df = CSV.read("employees.csv", DataFrame)

# Basic exploration
size(df)           # (1000, 6)
names(df)          # column names
describe(df)       # summary statistics
first(df, 5)       # head
last(df, 5)        # tail

# Selecting columns
df.salary          # column as Vector
df[!, :salary]     # same (! means copy not needed)
df[:, :salary]     # makes a copy
df[!, [:name, :salary]]   # multiple columns

# Filtering rows
subset(df, :age => ByRow(>(30)))          # age > 30
filter(:dept => ==("Engineering"), df)    # Engineering only
df[df.active .& (df.salary .> 100_000), :]  # Boolean indexing

# Transformations with transform
transform(df,
    :salary => (s -> s ./ 1000) => :salary_k,   # salary in thousands
    :age => ByRow(x -> x >= 65 ? "senior" : "regular") => :age_group
)

# GroupBy and aggregation
gdf = groupby(df, :dept)

# combine: apply aggregation functions to each group
summary = combine(gdf,
    :salary => mean => :mean_salary,
    :salary => std  => :std_salary,
    :salary => minimum => :min_salary,
    :salary => maximum => :max_salary,
    nrow => :headcount,
    :active => sum => :active_count,
)

# Sort
sort!(summary, :mean_salary, rev=true)

# Joins (similar to SQL)
managers = DataFrame(
    dept = ["Engineering", "Sales", "HR"],
    manager = ["Alice", "Bob", "Carol"]
)

# Inner join
innerjoin(df, managers, on=:dept)

# Left join — keeps all rows from left, fills missing for unmatched
leftjoin(df, managers, on=:dept)

# Anti-join — rows in left NOT in right
antijoin(df, managers, on=:dept)

# Reshape: wide to long and back
using DataFrames: stack, unstack

# Stack (wide → long / melt)
wide = DataFrame(id=1:5, q1=rand(5), q2=rand(5), q3=rand(5))
long = stack(wide, [:q1, :q2, :q3], variable_name=:quarter, value_name=:value)

# Unstack (long → wide / pivot)
wide_again = unstack(long, :id, :quarter, :value)

# Missing data handling
df_with_missing = DataFrame(x=[1, missing, 3, missing, 5], y=[1.0, 2.0, missing, 4.0, 5.0])
dropmissing(df_with_missing)          # drop rows with any missing
dropmissing(df_with_missing, :x)      # drop only where x is missing
coalesce.(df_with_missing.x, 0)       # replace missing with 0

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

using DifferentialEquations, Plots

# --- Ordinary Differential Equations ---

# The Lorenz system: classic chaotic ODE
# dx/dt = σ(y - x)
# dy/dt = x(ρ - z) - y
# dz/dt = xy - βz

function lorenz!(du, u, p, t)
    σ, ρ, β = p
    du[1] = σ * (u[2] - u[1])
    du[2] = u[1] * (ρ - u[3]) - u[2]
    du[3] = u[1] * u[2] - β * u[3]
end

u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 100.0)
p = (10.0, 28.0, 8/3)

prob = ODEProblem(lorenz!, u0, tspan, p)

# Tsit5: 4th/5th order Runge-Kutta, good default for non-stiff ODEs
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

# Solution interpolates — access any time point
sol(50.0)       # state at t=50
sol.t           # time points chosen by adaptive stepper
sol.u           # states at those time points

# Plot directly
plot(sol, idxs=(1, 2, 3), title="Lorenz Attractor")

# --- Stiff ODEs ---

# Robertson chemical kinetics — a classic stiff problem
function robertson!(du, u, p, t)
    k1, k2, k3 = p
    du[1] = -k1 * u[1] + k3 * u[2] * u[3]
    du[2] =  k1 * u[1] - k2 * u[2]^2 - k3 * u[2] * u[3]
    du[3] =  k2 * u[2]^2
end

u0 = [1.0, 0.0, 0.0]
tspan = (0.0, 1e11)
p = (0.04, 3e7, 1e4)

prob_stiff = ODEProblem(robertson!, u0, tspan, p)

# Rodas5: stiff solver, much faster than Tsit5 on stiff problems
sol_stiff = solve(prob_stiff, Rodas5(), reltol=1e-8, abstol=1e-8)

# --- Stochastic Differential Equations (SDEs) ---

# Geometric Brownian Motion: dS = μS dt + σS dW
function gbm_drift!(du, u, p, t)
    μ, σ = p
    du[1] = μ * u[1]
end

function gbm_diffusion!(du, u, p, t)
    μ, σ = p
    du[1] = σ * u[1]
end

u0 = [100.0]   # Initial price
tspan = (0.0, 1.0)   # 1 year
p = (0.05, 0.20)     # 5% drift, 20% volatility

prob_sde = SDEProblem(gbm_drift!, gbm_diffusion!, u0, tspan, p)

# Simulate 1000 paths
ensembleprob = EnsembleProblem(prob_sde)
sol_ensemble = solve(ensembleprob, EnsembleThreads(), trajectories=1000)

# Summary statistics across paths
using DifferentialEquations: EnsembleSummary
summary = EnsembleSummary(sol_ensemble, 0:0.01:1.0)
plot(summary, title="GBM: 1000 Paths (5th-95th percentile)")

# --- Parameter Estimation ---

using DifferentialEquations, Optimization, OptimizationOptimJL

# Generate synthetic data
true_params = [1.5, 1.0, 3.0, 1.0]
function lotka_volterra!(du, u, p, t)
    α, β, γ, δ = p
    du[1] =  α * u[1] - β * u[1] * u[2]
    du[2] = -γ * u[2] + δ * u[1] * u[2]
end

prob_lv = ODEProblem(lotka_volterra!, [1.0, 1.0], (0.0, 10.0), true_params)
sol_true = solve(prob_lv, Tsit5())
t_obs = 0.0:0.5:10.0
data = Array(sol_true(t_obs)) + 0.1 * randn(2, length(t_obs))

# Loss function: mean squared error between model and data
function loss(params, _)
    prob = ODEProblem(lotka_volterra!, [1.0, 1.0], (0.0, 10.0), params)
    sol = solve(prob, Tsit5(), saveat=0.5)
    if sol.retcode != ReturnCode.Success
        return Inf, sol
    end
    pred = Array(sol)
    l = sum(abs2, data .- pred)
    return l, sol
end

# Optimize with BFGS
optfunc = OptimizationFunction(loss, Optimization.AutoZygote())
optprob = OptimizationProblem(optfunc, [1.0, 1.0, 2.0, 1.5])
result = solve(optprob, BFGS())
println("Recovered parameters: ", result.u)
println("True parameters:      ", true_params)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84

using Flux
using Flux: DataLoader, train!
using Statistics

# Flux.jl: ML framework entirely in Julia — no Python backend
# This matters: you can differentiate through your ODE solver,
# your custom loss, and your model all in the same AD pass.

# --- Simple feedforward network ---
model = Chain(
    Dense(784 => 256, relu),
    Dropout(0.3),
    Dense(256 => 128, relu),
    Dropout(0.2),
    Dense(128 => 10),
)

# Flux uses Zygote.jl for automatic differentiation
params = Flux.params(model)

# Loss
loss_fn(x, y) = Flux.crossentropy(softmax(model(x)), y)

# Optimizer
opt = Adam(0.001)

# Training loop
function train_epoch!(model, loader, opt)
    total_loss = 0.0
    n_batches = 0
    for (x, y) in loader
        gs = gradient(params) do
            loss_fn(x, y)
        end
        Flux.update!(opt, params, gs)
        total_loss += loss_fn(x, y)
        n_batches += 1
    end
    return total_loss / n_batches
end

# --- Neural ODE: combine DiffEq and ML ---
using DifferentialEquations, Flux, SciMLSensitivity

# A neural network as the right-hand side of an ODE
neural_ode_rhs = Chain(Dense(2, 32, tanh), Dense(32, 2))

function neural_ode_forward(u0, p, tspan)
    function dudt!(du, u, p, t)
        du .= neural_ode_rhs(u)
    end
    prob = ODEProblem(dudt!, u0, tspan)
    return solve(prob, Tsit5(), saveat=0.1, sensealg=BacksolveAdjoint())
end

# Gradient flows through the ODE solver — this is what makes Julia special
# In Python, you'd need a specialized library (torchdiffeq); in Julia, it's built in.

# --- Lux.jl: explicit state, better for research ---
using Lux, Random, Zygote

rng = Random.default_rng()

# Lux separates model structure, parameters, and state explicitly
# (like PyTorch's functional API, or Flax in JAX)
model_lux = Lux.Chain(
    Lux.Dense(10 => 32, relu),
    Lux.Dense(32 => 1),
)

ps, st = Lux.setup(rng, model_lux)

x_input = randn(Float32, 10, 32)  # 10 features, batch of 32

# Forward pass: explicit params and state
y_pred, st_new = model_lux(x_input, ps, st)

# Gradient w.r.t. parameters
loss_lux(ps, st, x, y) = begin
    y_pred, st_new = model_lux(x, ps, st)
    return mean(abs2, y_pred .- y), st_new
end

gs = gradient(p -> first(loss_lux(p, st, x_input, zeros(Float32, 1, 32))), ps)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

using Turing, StatsPlots, Distributions

# Turing.jl: PPL built on Julia's multiple dispatch and AD
# Competes with Stan but in pure Julia (no external DSL)

# --- Bayesian Linear Regression ---
@model function bayesian_regression(x, y)
    # Priors
    α  ~ Normal(0, 10)           # intercept
    β  ~ Normal(0, 10)           # slope
    σ  ~ Exponential(1)          # noise standard deviation

    # Likelihood
    for i in eachindex(y)
        y[i] ~ Normal(α + β * x[i], σ)
    end
end

# Generate synthetic data
true_α, true_β, true_σ = 2.0, 3.0, 0.5
x_data = randn(100)
y_data = true_α .+ true_β .* x_data .+ true_σ .* randn(100)

# Sample with NUTS (No-U-Turn Sampler — same as Stan's default)
model = bayesian_regression(x_data, y_data)
chain = sample(model, NUTS(0.65), 2000; progress=false)

# Summary statistics
summarize(chain)
# parameters    mean       std    naive_se  mcse      ess   r̂
#   α           2.02      0.05     0.001    0.001   1752   1.00
#   β           2.98      0.05     0.001    0.001   1823   1.00
#   σ           0.51      0.04     0.001    0.001   1891   1.00

plot(chain)   # Trace plots and marginal distributions

# --- Hierarchical Model ---
@model function hierarchical_schools(J, sigma, y)
    μ ~ Normal(0, 10)
    τ ~ truncated(Cauchy(0, 2.5), lower=0)
    θ ~ filldist(Normal(μ, τ), J)

    for i in 1:J
        y[i] ~ Normal(θ[i], sigma[i])
    end
end

# Eight Schools data (classic hierarchical example)
J = 8
sigma = [15, 10, 16, 11, 9, 11, 10, 18]
y_schools = [28, 8, -3, 7, -1, 1, 18, 12]

chain_schools = sample(
    hierarchical_schools(J, sigma, y_schools),
    NUTS(0.65),
    MCMCThreads(),  # parallel chains using Julia threads
    1000, 4         # 1000 samples, 4 chains
)

# --- Variational Inference (faster for large models) ---
q = vi(model, ADVI(10, 1000))  # 10 samples, 1000 iterations
posterior_samples = rand(q, 1000)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

using Symbolics

# Declare symbolic variables
@variables x y z t

# Symbolic expressions
expr = sin(x)^2 + cos(x)^2
simplify(expr)   # → 1

poly = (x + y)^3
expand(poly)     # → x³ + 3x²y + 3xy² + y³

# Differentiation
D = Differential(x)

f = x^3 + 2*x^2 + 5*x + 1
D(f)           # symbolic: 3x² + 4x + 5
expand_derivatives(D(f))  # → 3(x^2) + 4x + 5

# Higher-order derivatives
D2 = Differential(x)^2
expand_derivatives(D2(f))  # → 6x + 4

# Partial derivatives
g = x^2 * y + sin(x * y)
expand_derivatives(Differential(x)(g))  # → 2xy + y*cos(x*y)
expand_derivatives(Differential(y)(g))  # → x^2 + x*cos(x*y)

# Jacobian and Hessian
@variables q[1:3]   # vector of symbolic variables
f_vec = [q[1]^2 + q[2], q[2]*q[3], sin(q[1])]

jac = Symbolics.jacobian(f_vec, collect(q))
# 3×3 symbolic Jacobian matrix

hess = Symbolics.hessian(q[1]^2 + q[2]^2 + q[3]^2, collect(q))
# 3×3 symbolic Hessian (diagonal: [2, 2, 2])

# Code generation: symbolic → optimized Julia code
fn, fn! = build_function(f_vec, collect(q))
# fn is a callable function — fast runtime code generated from symbolic expression

# Solve symbolic equations
@variables a b c
roots = Symbolics.solve_for([a*x^2 + b*x + c ~ 0], [x])

# Integration with ModelingToolkit.jl for physical systems
using ModelingToolkit

@variables t
D_t = Differential(t)
@variables x(t) y(t)   # time-varying symbolic variables

# Define an ODE system symbolically
eqs = [
    D_t(x) ~ -x + y,
    D_t(y) ~ x - y^2,
]

@named sys = ODESystem(eqs, t)
sys_simplified = structural_simplify(sys)

# Convert to DifferentialEquations.jl problem automatically
prob_sym = ODEProblem(sys_simplified, [x => 1.0, y => 0.5], (0.0, 10.0))
sol_sym = solve(prob_sym, Tsit5())

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

using Base.Threads
println("Threads: ", nthreads())   # Set with JULIA_NUM_THREADS env var

# Start Julia with: julia --threads 8
# or: JULIA_NUM_THREADS=8 julia

# Simple parallel loop
function parallel_sum_squares(n::Int)
    result = Atomic{Float64}(0.0)
    @threads for i in 1:n
        atomic_add!(result, Float64(i)^2)
    end
    return result[]
end

# Better pattern: thread-local accumulators (avoid atomic contention)
function parallel_sum_squares_fast(n::Int)
    partial = zeros(Float64, nthreads())
    @threads for i in 1:n
        partial[threadid()] += Float64(i)^2
    end
    return sum(partial)
end

# @threads on realistic work: parallel matrix row operations
function parallel_normalize_rows!(A::Matrix{Float64})
    @threads for i in axes(A, 1)
        row_norm = norm(view(A, i, :))
        if row_norm > 0
            A[i, :] ./= row_norm
        end
    end
    return A
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

# @spawn: create asynchronous tasks, more flexible than @threads
function parallel_map(f, items)
    tasks = [Threads.@spawn f(item) for item in items]
    return [fetch(t) for t in tasks]
end

# Useful for IO-heavy or heterogeneous workloads
function fetch_and_process(urls::Vector{String})
    tasks = map(urls) do url
        Threads.@spawn begin
            # Each URL fetched concurrently on a thread pool
            data = download_data(url)   # IO-bound work
            process(data)               # CPU-bound work
        end
    end
    return fetch.(tasks)
end

# Recursive parallel algorithms
function parallel_fib(n::Int)
    n <= 20 && return fib_sequential(n)   # Don't spawn tiny tasks
    t = Threads.@spawn parallel_fib(n - 1)
    r = parallel_fib(n - 2)
    return fetch(t) + r
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

using Distributed

# Add worker processes
addprocs(4)   # 4 additional processes (workers)
# addprocs([("worker1.example.com", 2), ("worker2.example.com", 4)])  # Remote workers
nprocs()      # total (1 main + 4 workers)
nworkers()    # 4

# Load code on all workers
@everywhere using LinearAlgebra

# @everywhere: run on all processes
@everywhere function expensive_work(x::Float64)
    # Simulate heavy computation
    return sum(sin(x * k) for k in 1:10_000)
end

# pmap: parallel map (auto-distributes work)
inputs = rand(1000)
results = pmap(expensive_work, inputs)

# @distributed: parallel for loop with reduction
total = @distributed (+) for i in 1:1_000_000
    rand()^2
end
# Approximates π/4 via Monte Carlo

# SharedArray: shared memory between processes on same machine
using SharedArrays
A = SharedArray{Float64}(1000, 1000)
@distributed for i in 1:1000
    A[i, :] = rand(1000)
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

using CUDA

# Check GPU
CUDA.versioninfo()
device()

# Move arrays to GPU — operations are JIT-compiled for GPU
x_cpu = rand(Float32, 10_000_000)
x_gpu = CuArray(x_cpu)   # transfer to GPU memory

# Standard Julia operations work on GPU arrays
y_gpu = x_gpu .^ 2 + 2 .* x_gpu   # element-wise, runs on GPU
z_gpu = sum(y_gpu)                 # GPU reduction

# Move back
y_cpu = Array(y_gpu)

# GPU matrix multiplication (CUBLAS)
A_gpu = CUDA.randn(Float32, 1000, 1000)
B_gpu = CUDA.randn(Float32, 1000, 1000)
C_gpu = A_gpu * B_gpu   # Uses CUBLAS SGEMM automatically

# Write custom GPU kernels in Julia
function vector_add_kernel!(c, a, b)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(c)
        c[i] = a[i] + b[i]
    end
    return nothing
end

n = 10_000
a_gpu = CUDA.ones(Float32, n)
b_gpu = CUDA.ones(Float32, n) .* 2
c_gpu = CUDA.zeros(Float32, n)

threads = 256
blocks = cld(n, threads)   # ceil(n / threads)
@cuda threads=threads blocks=blocks vector_add_kernel!(c_gpu, a_gpu, b_gpu)
CUDA.synchronize()
all(Array(c_gpu) .== 3.0f0)   # true

# @turbo: CPU SIMD vectorization with LoopVectorization.jl
using LoopVectorization

function dot_turbo(a::Vector{Float64}, b::Vector{Float64})
    s = 0.0
    @turbo for i in eachindex(a)
        s += a[i] * b[i]
    end
    return s
end

# @turbo often approaches or exceeds BLAS for simple operations
a = rand(10_000); b = rand(10_000)
@btime dot_turbo($a, $b)   # often < 5μs on modern hardware

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# @simd: hint to compiler that loop iterations are independent, enable SIMD
function sum_simd(v::Vector{Float64})
    s = 0.0
    @simd for x in v
        s += x
    end
    return s
end

# @inbounds + @simd: disable bounds checking and enable SIMD together
function inner_product_fast(a::Vector{Float64}, b::Vector{Float64})
    s = 0.0
    @inbounds @simd for i in eachindex(a)
        s += a[i] * b[i]
    end
    return s
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37

using PythonCall

# PythonCall.jl is the modern Python interop (replaces PyCall.jl)
# Uses a separate Python environment managed by CondaPkg.jl

np = pyimport("numpy")
pd = pyimport("pandas")

# NumPy arrays
arr = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
julia_arr = pyconvert(Vector{Float64}, arr)

# Go the other way: Julia → Python
jl_data = [1.0, 2.0, 3.0]
py_data = Py(jl_data)   # wraps without copy when types are compatible

# Call pandas
df_py = pd.DataFrame(pydict(Dict("a" => [1,2,3], "b" => [4,5,6])))
py_result = df_py.groupby("a").sum()

# Use scikit-learn from Julia
sklearn = pyimport("sklearn.linear_model")
LinearRegression = sklearn.LinearRegression

X = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])
y = np.array([1.0, 2.0, 3.0])

model_py = LinearRegression()
model_py.fit(X, y)
println("Coefficients: ", pyconvert(Vector{Float64}, model_py.coef_))

# Use matplotlib for plotting in Julia notebooks
plt = pyimport("matplotlib.pyplot")
x_py = np.linspace(0, 2*3.14159, 100)
plt.plot(x_py, np.sin(x_py))
plt.savefig("from_julia.png")
plt.close()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

using RCall

# R expressions
R"summary(c(1, 2, 3, 4, 5))"

# Transfer data
x = randn(100)
@rput x   # send Julia variable to R
R"hist(x, main='Julia data in R')"
R"shapiro_result <- shapiro.test(x)"
@rget shapiro_result   # get R result back

# Run R scripts
R"""
library(ggplot2)
data(mtcars)
p <- ggplot(mtcars, aes(x=wt, y=mpg)) +
    geom_point() +
    geom_smooth(method='lm') +
    labs(title='Weight vs MPG')
ggsave('ggplot_from_julia.png', p)
"""

# Use R packages not available in Julia
R"""
library(forecast)
data <- ts($x, frequency=12)
model <- auto.arima(data)
forecast_result <- forecast(model, h=24)
"""
@rget forecast_result

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

# ccall has zero overhead — no marshaling, no FFI wrapper layer
# The JIT sees through ccall and generates the same code as a direct call

# Call C standard library
# ccall((function_name, library), return_type, (arg_types...,), args...)

# strlen from libc
s = "Hello, World!"
len = ccall(:strlen, Csize_t, (Cstring,), s)
println("Length: $len")   # 13

# Math functions (often already inlined by Julia, but demonstrates the syntax)
val = ccall((:sin, "libm"), Cdouble, (Cdouble,), 1.5)

# Call custom C function
# Suppose you have: double dot_product(double *a, double *b, int n);
a = [1.0, 2.0, 3.0]
b = [4.0, 5.0, 6.0]
n = length(a)

result = ccall(
    (:dot_product, "libmymath"),  # function name and library
    Float64,                       # return type
    (Ptr{Float64}, Ptr{Float64}, Cint),  # argument types
    a, b, n                        # arguments (Julia arrays pass as pointers)
)

# Call Fortran (Fortran mangles names with underscore)
# DGEMM from BLAS (Fortran subroutine — no return value → Cvoid)
using LinearAlgebra.BLAS: @blasfunc

# Julia's LinearAlgebra calls BLAS via ccall internally.
# Example of what that looks like:
function my_dgemm!(C, A, B, α=1.0, β=0.0)
    m, k = size(A)
    k2, n = size(B)
    @assert k == k2
    @assert size(C) == (m, n)

    ccall((@blasfunc(dgemm_), BLAS.libblastrampoline),
        Cvoid,
        (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
         Ref{Float64}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
         Ref{Float64}, Ptr{Float64}, Ref{BlasInt}),
        'N', 'N', m, n, k, α, A, m, B, k, β, C, m
    )
    return C
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

# C++ interop via CxxWrap.jl (requires a C++ wrapper module)
# In C++:
# #include <jlcxx/jlcxx.hpp>
# JLCXX_MODULE define_julia_module(jlcxx::Module& mod) {
#     mod.add_type<MyClass>("MyClass")
#        .constructor<int, double>()
#        .method("compute", &MyClass::compute);
# }

using CxxWrap
@wrapmodule(() -> "libmymodule")
@initcxx

obj = MyClass(42, 3.14)
result = compute(obj)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

# In the Julia REPL, press ] to enter Pkg mode

# Add packages
pkg> add DataFrames CSV Plots BenchmarkTools

# Add a specific version
pkg> add DataFrames@1.6.0

# Add from GitHub (unregistered package)
pkg> add https://github.com/username/MyPackage.jl

# Add a local path
pkg> add /path/to/MyPackage

# Remove a package
pkg> remove Plots

# Update all packages
pkg> update

# Update a specific package
pkg> update DataFrames

# Show package status
pkg> status
# Project lunarOps_demo v0.1.0
#   [a93c6f00] DataFrames v1.6.1
#   [91a5bcdd] Plots v1.40.0
#   [6f49c342] CSV v0.10.12

# Show outdated packages
pkg> outdated

# Pin a package version (prevent updates)
pkg> pin DataFrames

# Resolve and precompile all dependencies
pkg> precompile

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

# Activate an environment (creates if not exists)
using Pkg
Pkg.activate("./MyProject")
# or from the REPL: pkg> activate ./MyProject

# Project.toml: human-edited, version constraints
# [deps]
# DataFrames = "a93c6f00-e57d-5684-b890-1f30009b9b0e"
# CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
#
# [compat]
# julia = "1.9"
# DataFrames = "1"
# CSV = "0.10"

# Manifest.toml: auto-generated, exact versions of all deps including transitive
# Never edit Manifest.toml by hand.
# Commit both for reproducible research; commit only Project.toml for libraries.

# Instantiate from Project.toml (recreate exact environment from Manifest.toml)
Pkg.instantiate()

# Example: set up a project from scratch
Pkg.activate("my_analysis")
Pkg.add(["DataFrames", "CSV", "Plots", "Statistics", "BenchmarkTools"])
Pkg.status()

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# Generate package structure
julia -e 'using Pkg; Pkg.generate("MyPackage")'

# Creates:
# MyPackage/
# ├── Project.toml
# └── src/
#     └── MyPackage.jl

# Better: use PkgTemplates.jl for a full setup

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

using PkgTemplates

t = Template(
    user="yourusername",
    license="MIT",
    plugins=[
        Git(ssh=true),
        GitHubActions(; x86=false),
        Codecov(),
        Documenter{GitHubActions}(),
        Tests(; project=true),
    ]
)

t("MyPackage")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# src/MyPackage.jl
module MyPackage

export greet, PhysicsModel, simulate

include("physics.jl")
include("io.jl")

function greet(name::String)
    println("Hello from MyPackage, $name!")
end

end  # module

# test/runtests.jl
using MyPackage
using Test

@testset "MyPackage" begin
    @testset "greet" begin
        @test greet("world") === nothing   # returns nothing, doesn't throw
    end

    @testset "PhysicsModel" begin
        m = PhysicsModel(mass=1.0, stiffness=10.0)
        sol = simulate(m, tspan=(0.0, 10.0))
        @test length(sol.t) > 1
        @test all(isfinite, sol.u)
    end
end

# Run tests
# pkg> test MyPackage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# The General registry is the default — contains thousands of packages
# Check what's available at: https://juliahub.com/ui/Packages

# Add a private registry (for internal packages)
pkg> registry add https://github.com/myorg/MyRegistry

# Register your own package
# 1. Create the package with proper Project.toml
# 2. Tag a version: git tag v0.1.0
# 3. Use Registrator.jl or JuliaHub to register
# 4. Or for local use: keep in the local path, addprefs

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

# Check precompilation cache status
using Pkg
Pkg.precompile()   # Force precompilation of all installed packages

# Packages can declare precompile workloads with PrecompileTools.jl
# (good packages already do this; check their source)

using PrecompileTools

@setup_workload begin
    # Code here runs during package precompilation, not at load time
    x = collect(1.0:10.0)
    y = sin.(x)

    @compile_workload begin
        # This call is compiled and cached at precompile time
        plot(x, y)       # ← "time to first plot" is now paid during ]precompile
    end
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

using PackageCompiler

# Create a custom system image that bakes in precompiled code
# This eliminates startup latency entirely for the included packages

create_sysimage(
    ["DataFrames", "CSV", "Plots", "DifferentialEquations"],
    sysimage_path = "my_sysimage.so",
    precompile_execution_file = "warmup.jl"  # Run this to determine what to bake in
)

# warmup.jl: a script that exercises the code you want fast
# It runs during sysimage creation, and the compiled code is baked in.
# using DataFrames; df = DataFrame(a=1:100, b=randn(100)); combine(df, :a => sum)
# using Plots; x = 1:10; plot(x, sin.(x))

# Use the sysimage
# julia --sysimage my_sysimage.so my_script.jl

# Create a deployable app (no Julia installation needed on target machine)
create_app(
    "MyProject",           # directory with Project.toml and src/
    "build/MyApp",         # output directory
    precompile_execution_file = "warmup.jl",
    include_lazy_artifacts = true,
)
# Creates: build/MyApp/bin/MyApp — a self-contained binary

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

# Julia: type-stable tight loop
function mandelbrot_count(c::ComplexF64, max_iter::Int)::Int
    z = zero(ComplexF64)
    for i in 1:max_iter
        z = z^2 + c
        abs2(z) > 4.0 && return i
    end
    return max_iter
end

function mandelbrot_julia(n::Int, max_iter::Int=100)
    xs = LinRange(-2.5, 1.0, n)
    ys = LinRange(-1.25, 1.25, n)
    result = Matrix{Int}(undef, n, n)
    @threads for j in 1:n, i in 1:n
        result[i, j] = mandelbrot_count(complex(xs[i], ys[j]), max_iter)
    end
    return result
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

# Python: equivalent NumPy/Numba version
import numpy as np

def mandelbrot_numpy(n, max_iter=100):
    xs = np.linspace(-2.5, 1.0, n)
    ys = np.linspace(-1.25, 1.25, n)
    X, Y = np.meshgrid(xs, ys)
    C = X + 1j * Y
    Z = np.zeros_like(C)
    result = np.zeros(C.shape, dtype=int)
    for i in range(max_iter):
        mask = np.abs(Z) <= 2
        Z[mask] = Z[mask]**2 + C[mask]
        result[mask] += 1
    return result
# NumPy version is vectorized but still ~5–10x slower for n=1000
# Julia threaded version: ~40ms. NumPy: ~200ms. Pure Python: ~8000ms.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

# Python: ~50ms for basic script
time python -c "import numpy; print('hello')"
# real    0m0.089s

# Julia: 2–4 seconds for basic script (JIT warmup)
time julia -e "println(\"hello\")"
# real    0m2.341s

# Julia with sysimage: ~0.2s
time julia --sysimage my_sysimage.so -e "println(\"hello\")"
# real    0m0.198s

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

# Julia: parallelism is first-class
using Base.Threads

# This just works — no GIL, no multiprocessing overhead
@threads for i in 1:N
    results[i] = expensive_computation(data[i])
end

# Distributed across machines is one addprocs() call
using Distributed
addprocs(4)
@distributed for i in 1:N
    # runs on workers
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# Python: GIL prevents true thread parallelism for CPU work
# Must use multiprocessing, which has significant overhead
from multiprocessing import Pool
with Pool(4) as p:
    results = p.map(expensive_computation, data)
# Good performance but: fork overhead, pickle serialization, memory copying

# Numba with parallel=True works for some cases
@nb.njit(parallel=True)
def parallel_loop(data):
    for i in nb.prange(len(data)):
        # only works if the kernel is Numba-compatible