trending/ai-engineering-from-scratch
1
0
Fork 0
mirror of https://github.com/rohitg00/ai-engineering-from-scratch.git synced 2026-10-02 01:54:39 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #188 from rohitg00/feat/julia-coverage-phase-2-7

Browse Source 
feat(julia): ML fundamentals + transformer Julia ports (phase 2 + 7)
This commit is contained in:
Rohit Ghumare
2026-05-26 10:34:34 +01:00
committed by GitHub
parent 8ad665dec9 376038b217
commit 4cd80622c5
9 changed files with 2313 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files
catalog.json
+12 -4
View File
@@ -6,7 +6,7 @@
"skills": 378,
"prompts": 99,
"agents": 0,
"code_files": 465
"code_files": 473
},
"phases": [
{
@@ -812,7 +812,8 @@
"has_quiz": true,
"has_notebook": false,
"code_files": [
"linear_regression.py"
"linear_regression.py",
"main.jl"
],
"outputs": [
{
@@ -841,7 +842,8 @@
"has_quiz": true,
"has_notebook": false,
"code_files": [
"logistic_regression.py"
"logistic_regression.py",
"main.jl"
],
"outputs": [
{
@@ -892,6 +894,7 @@
"has_quiz": true,
"has_notebook": false,
"code_files": [
"main.jl",
"svm.py"
],
"outputs": [
@@ -996,7 +999,8 @@
"has_quiz": true,
"has_notebook": false,
"code_files": [
"evaluation.py"
"evaluation.py",
"main.jl"
],
"outputs": [
{
@@ -3914,6 +3918,7 @@
"has_quiz": false,
"has_notebook": true,
"code_files": [
"main.jl",
"main.py"
],
"outputs": [
@@ -3942,6 +3947,7 @@
"has_quiz": true,
"has_notebook": false,
"code_files": [
"main.jl",
"self_attention.py"
],
"outputs": [
@@ -3993,6 +3999,7 @@
"has_quiz": false,
"has_notebook": true,
"code_files": [
"main.jl",
"main.py"
],
"outputs": [
@@ -4021,6 +4028,7 @@
"has_quiz": false,
"has_notebook": true,
"code_files": [
"main.jl",
"main.py"
],
"outputs": [
phases/02-ml-fundamentals/02-linear-regression/code/main.jl
+291
View File
@@ -0,0 +1,291 @@
# Linear regression in Julia. Closed-form normal equation and batch
# gradient descent, plus multiple linear regression and a ridge penalty.
# Stdlib only. Sources:
# https://docs.julialang.org/en/v1/manual/types/
# https://docs.julialang.org/en/v1/stdlib/Statistics/
# https://docs.julialang.org/en/v1/stdlib/Random/
using Random
using Statistics
using Printf
function generate_simple_data(; n::Int=100, true_w::Float64=3.0, true_b::Float64=7.0,
noise::Float64=2.0, seed::Int=42)
rng = MersenneTwister(seed)
xs = [10.0 * rand(rng) for _ in 1:n]
ys = [true_w * x + true_b + noise * randn(rng) for x in xs]
return xs, ys
end
mutable struct GDLinearRegression
w::Float64
b::Float64
lr::Float64
history::Vector{Float64}
end
GDLinearRegression(lr::Float64) = GDLinearRegression(0.0, 0.0, lr, Float64[])
function predict(model::GDLinearRegression, xs::Vector{Float64})
return [model.w * x + model.b for x in xs]
end
function cost(model::GDLinearRegression, xs::Vector{Float64}, ys::Vector{Float64})
preds = predict(model, xs)
return sum((preds .- ys) .^ 2) / length(ys)
end
function fit_gd!(model::GDLinearRegression, xs::Vector{Float64}, ys::Vector{Float64};
epochs::Int=1000, print_every::Int=200)
n = length(ys)
for epoch in 0:(epochs - 1)
preds = predict(model, xs)
errs = preds .- ys
dw = (2.0 / n) * sum(errs .* xs)
db = (2.0 / n) * sum(errs)
model.w -= model.lr * dw
model.b -= model.lr * db
c = cost(model, xs, ys)
push!(model.history, c)
if epoch % print_every == 0
@printf(" epoch %4d cost=%.4f w=%.4f b=%.4f\n", epoch, c, model.w, model.b)
end
end
return model
end
function r_squared(ys::Vector{Float64}, preds::Vector{Float64})
y_mean = mean(ys)
ss_res = sum((ys .- preds) .^ 2)
ss_tot = sum((ys .- y_mean) .^ 2)
if ss_tot == 0.0
return ss_res == 0.0 ? 1.0 : 0.0
end
return 1.0 - ss_res / ss_tot
end
function fit_normal_equation(xs::Vector{Float64}, ys::Vector{Float64})
x_mean = mean(xs)
y_mean = mean(ys)
num = sum((xs .- x_mean) .* (ys .- y_mean))
den = sum((xs .- x_mean) .^ 2)
if den == 0.0
return 0.0, y_mean
end
w = num / den
b = y_mean - w * x_mean
return w, b
end
mutable struct MultiLinearRegression
weights::Vector{Float64}
bias::Float64
lr::Float64
end
MultiLinearRegression(n_features::Int, lr::Float64) =
MultiLinearRegression(zeros(n_features), 0.0, lr)
function predict_multi(model::MultiLinearRegression, X::Vector{Vector{Float64}})
return [sum(model.weights .* row) + model.bias for row in X]
end
function fit_multi!(model::MultiLinearRegression, X::Vector{Vector{Float64}},
ys::Vector{Float64}; epochs::Int=1000, print_every::Int=200)
n = length(ys)
n_features = length(X[1])
for epoch in 0:(epochs - 1)
preds = predict_multi(model, X)
errs = preds .- ys
for j in 1:n_features
grad = (2.0 / n) * sum(errs[i] * X[i][j] for i in 1:n)
model.weights[j] -= model.lr * grad
end
model.bias -= model.lr * ((2.0 / n) * sum(errs))
if epoch % print_every == 0
mse = sum(errs .^ 2) / n
@printf(" epoch %4d cost=%.4f\n", epoch, mse)
end
end
return model
end
function standardize(X::Vector{Vector{Float64}})
n_samples = length(X)
n_features = length(X[1])
means = [mean(X[i][j] for i in 1:n_samples) for j in 1:n_features]
stds = Float64[]
for j in 1:n_features
v = sum((X[i][j] - means[j]) ^ 2 for i in 1:n_samples) / n_samples
push!(stds, sqrt(v))
end
X_scaled = [Float64[
stds[j] > 0 ? (X[i][j] - means[j]) / stds[j] : 0.0
for j in 1:n_features
] for i in 1:n_samples]
return X_scaled, means, stds
end
function generate_house_data(; n::Int=100, seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Float64[]
for _ in 1:n
size = 500 + 2500 * rand(rng)
bedrooms = float(rand(rng, 1:5))
age = 50 * rand(rng)
price = 50 * size + 10000 * bedrooms - 1000 * age + 50000 + 20000 * randn(rng)
push!(X, Float64[size, bedrooms, age])
push!(ys, price)
end
return X, ys
end
mutable struct RidgeRegression
weights::Vector{Float64}
bias::Float64
lr::Float64
alpha::Float64
end
RidgeRegression(n_features::Int, lr::Float64, alpha::Float64) =
RidgeRegression(zeros(n_features), 0.0, lr, alpha)
function predict_ridge(model::RidgeRegression, X::Vector{Vector{Float64}})
return [sum(model.weights .* row) + model.bias for row in X]
end
function fit_ridge!(model::RidgeRegression, X::Vector{Vector{Float64}},
ys::Vector{Float64}; epochs::Int=1000, print_every::Int=200)
n = length(ys)
n_features = length(X[1])
for epoch in 0:(epochs - 1)
preds = predict_ridge(model, X)
errs = preds .- ys
mse_v = sum(errs .^ 2) / n
reg = model.alpha * sum(model.weights .^ 2)
for j in 1:n_features
grad = (2.0 / n) * sum(errs[i] * X[i][j] for i in 1:n)
grad += 2 * model.alpha * model.weights[j]
model.weights[j] -= model.lr * grad
end
model.bias -= model.lr * ((2.0 / n) * sum(errs))
if epoch % print_every == 0
@printf(" epoch %4d cost=%.4f L2=%.4f\n", epoch, mse_v + reg, reg)
end
end
return model
end
function demo_simple_regression()
println("=" ^ 60)
println("LINEAR REGRESSION (GRADIENT DESCENT)")
println("=" ^ 60)
xs, ys = generate_simple_data()
@printf("\nGenerated %d samples, true y = 3x + 7 + noise\n", length(xs))
model = GDLinearRegression(0.005)
fit_gd!(model, xs, ys; epochs=1000, print_every=200)
preds = predict(model, xs)
@printf("\nLearned: y = %.4fx + %.4f\n", model.w, model.b)
@printf("R^2: %.4f\n", r_squared(ys, preds))
return xs, ys
end
function demo_normal_equation(xs::Vector{Float64}, ys::Vector{Float64})
println("\n" * "=" ^ 60)
println("LINEAR REGRESSION (NORMAL EQUATION)")
println("=" ^ 60)
w, b = fit_normal_equation(xs, ys)
preds = [w * x + b for x in xs]
@printf("\nClosed-form: y = %.4fx + %.4f\n", w, b)
@printf("R^2: %.4f\n", r_squared(ys, preds))
end
function demo_multiple_regression()
println("\n" * "=" ^ 60)
println("MULTIPLE LINEAR REGRESSION (3 FEATURES)")
println("=" ^ 60)
X_raw, ys_raw = generate_house_data()
X_scaled, _, _ = standardize(X_raw)
y_mean = mean(ys_raw)
y_std = std(ys_raw; corrected=false)
ys_scaled = [(y - y_mean) / y_std for y in ys_raw]
model = MultiLinearRegression(3, 0.01)
fit_multi!(model, X_scaled, ys_scaled; epochs=1000, print_every=200)
preds = predict_multi(model, X_scaled)
@printf("\nStandardized weights: [%.4f, %.4f, %.4f]\n",
model.weights[1], model.weights[2], model.weights[3])
@printf("Standardized bias: %.4f\n", model.bias)
@printf("R^2 (scaled space): %.4f\n", r_squared(ys_scaled, preds))
return X_scaled, ys_scaled, model
end
function demo_ridge(X_scaled::Vector{Vector{Float64}}, ys_scaled::Vector{Float64},
plain_model::MultiLinearRegression)
println("\n" * "=" ^ 60)
println("RIDGE REGRESSION (L2)")
println("=" ^ 60)
ridge = RidgeRegression(3, 0.01, 0.1)
fit_ridge!(ridge, X_scaled, ys_scaled; epochs=1000, print_every=200)
@printf("\nRidge weights: [%.4f, %.4f, %.4f]\n",
ridge.weights[1], ridge.weights[2], ridge.weights[3])
@printf("Plain weights: [%.4f, %.4f, %.4f]\n",
plain_model.weights[1], plain_model.weights[2], plain_model.weights[3])
println("Ridge shrinks weights toward zero through the L2 penalty.")
end
function demo_train_test_split()
println("\n" * "=" ^ 60)
println("TRAIN/TEST SPLIT")
println("=" ^ 60)
xs, ys = generate_simple_data()
split = Int(round(0.8 * length(xs)))
xs_train = xs[1:split]
xs_test = xs[(split + 1):end]
ys_train = ys[1:split]
ys_test = ys[(split + 1):end]
model = GDLinearRegression(0.005)
fit_gd!(model, xs_train, ys_train; epochs=1000, print_every=500)
train_r2 = r_squared(ys_train, predict(model, xs_train))
test_r2 = r_squared(ys_test, predict(model, xs_test))
@printf("\nTrain R^2: %.4f\n", train_r2)
@printf("Test R^2: %.4f\n", test_r2)
end
function main()
xs, ys = demo_simple_regression()
demo_normal_equation(xs, ys)
X_scaled, ys_scaled, plain_model = demo_multiple_regression()
demo_ridge(X_scaled, ys_scaled, plain_model)
demo_train_test_split()
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
phases/02-ml-fundamentals/03-logistic-regression/code/main.jl
+385
View File
@@ -0,0 +1,385 @@
# Logistic regression in Julia. Sigmoid + binary cross-entropy gradient
# descent for two classes, plus multi-class softmax regression. Reports
# confusion-matrix metrics. Stdlib only. Sources:
# https://docs.julialang.org/en/v1/manual/mathematical-operations/
# https://docs.julialang.org/en/v1/stdlib/Random/
# https://docs.julialang.org/en/v1/stdlib/Statistics/
using Random
using Statistics
using Printf
function sigmoid(z::Float64)::Float64
z_clip = clamp(z, -500.0, 500.0)
return 1.0 / (1.0 + exp(-z_clip))
end
function generate_two_class_data(; n::Int=200, seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Int[]
half = n ÷ 2
for _ in 1:half
push!(X, Float64[2.0 + randn(rng), 2.0 + randn(rng)])
push!(ys, 0)
end
for _ in 1:half
push!(X, Float64[5.0 + randn(rng), 5.0 + randn(rng)])
push!(ys, 1)
end
perm = randperm(rng, length(X))
return X[perm], ys[perm]
end
mutable struct LogisticRegression
weights::Vector{Float64}
bias::Float64
lr::Float64
history::Vector{Float64}
end
LogisticRegression(n_features::Int, lr::Float64) =
LogisticRegression(zeros(n_features), 0.0, lr, Float64[])
function predict_proba(model::LogisticRegression, x::Vector{Float64})::Float64
z = sum(model.weights .* x) + model.bias
return sigmoid(z)
end
function predict_class(model::LogisticRegression, x::Vector{Float64};
threshold::Float64=0.5)::Int
return predict_proba(model, x) >= threshold ? 1 : 0
end
function bce_loss(model::LogisticRegression, X::Vector{Vector{Float64}}, ys::Vector{Int})
n = length(ys)
total = 0.0
for i in 1:n
p = clamp(predict_proba(model, X[i]), 1e-15, 1 - 1e-15)
total += ys[i] * log(p) + (1 - ys[i]) * log(1 - p)
end
return -total / n
end
function fit_logistic!(model::LogisticRegression, X::Vector{Vector{Float64}},
ys::Vector{Int}; epochs::Int=1000, print_every::Int=200)
n = length(ys)
n_features = length(X[1])
for epoch in 0:(epochs - 1)
dw = zeros(n_features)
db = 0.0
for i in 1:n
p = predict_proba(model, X[i])
err = p - ys[i]
for j in 1:n_features
dw[j] += err * X[i][j]
end
db += err
end
for j in 1:n_features
model.weights[j] -= model.lr * (dw[j] / n)
end
model.bias -= model.lr * (db / n)
loss = bce_loss(model, X, ys)
push!(model.history, loss)
if epoch % print_every == 0
@printf(" epoch %4d loss=%.4f w=[%.3f, %.3f] b=%.3f\n",
epoch, loss, model.weights[1], model.weights[2], model.bias)
end
end
return model
end
function accuracy(model::LogisticRegression, X::Vector{Vector{Float64}}, ys::Vector{Int})
correct = 0
for i in 1:length(ys)
if predict_class(model, X[i]) == ys[i]
correct += 1
end
end
return correct / length(ys)
end
struct ClassificationMetrics
tp::Int
tn::Int
fp::Int
fn::Int
end
function build_metrics(y_true::Vector{Int}, y_pred::Vector{Int})
tp = sum(1 for i in 1:length(y_true) if y_true[i] == 1 && y_pred[i] == 1)
tn = sum(1 for i in 1:length(y_true) if y_true[i] == 0 && y_pred[i] == 0)
fp = sum(1 for i in 1:length(y_true) if y_true[i] == 0 && y_pred[i] == 1)
fn = sum(1 for i in 1:length(y_true) if y_true[i] == 1 && y_pred[i] == 0)
return ClassificationMetrics(tp, tn, fp, fn)
end
metric_accuracy(m::ClassificationMetrics) =
(m.tp + m.tn + m.fp + m.fn) > 0 ? (m.tp + m.tn) / (m.tp + m.tn + m.fp + m.fn) : 0.0
metric_precision(m::ClassificationMetrics) =
(m.tp + m.fp) > 0 ? m.tp / (m.tp + m.fp) : 0.0
metric_recall(m::ClassificationMetrics) =
(m.tp + m.fn) > 0 ? m.tp / (m.tp + m.fn) : 0.0
function metric_f1(m::ClassificationMetrics)
p = metric_precision(m)
r = metric_recall(m)
return (p + r) > 0 ? 2 * p * r / (p + r) : 0.0
end
function print_report(m::ClassificationMetrics)
println("\n Confusion Matrix:")
println(" Predicted")
println(" Pos Neg")
@printf(" Actual Pos %4d %4d\n", m.tp, m.fn)
@printf(" Actual Neg %4d %4d\n", m.fp, m.tn)
@printf("\n Accuracy: %.4f\n", metric_accuracy(m))
@printf(" Precision: %.4f\n", metric_precision(m))
@printf(" Recall: %.4f\n", metric_recall(m))
@printf(" F1 Score: %.4f\n", metric_f1(m))
end
function softmax(scores::Vector{Float64})::Vector{Float64}
m = maximum(scores)
e = [exp(s - m) for s in scores]
s = sum(e)
return e ./ s
end
mutable struct SoftmaxRegression
weights::Vector{Vector{Float64}}
biases::Vector{Float64}
lr::Float64
n_features::Int
n_classes::Int
end
function SoftmaxRegression(n_features::Int, n_classes::Int, lr::Float64)
SoftmaxRegression(
[zeros(n_features) for _ in 1:n_classes],
zeros(n_classes),
lr,
n_features,
n_classes,
)
end
function predict_proba_softmax(model::SoftmaxRegression, x::Vector{Float64})::Vector{Float64}
scores = [sum(model.weights[k] .* x) + model.biases[k] for k in 1:model.n_classes]
return softmax(scores)
end
function predict_class_softmax(model::SoftmaxRegression, x::Vector{Float64})::Int
probs = predict_proba_softmax(model, x)
return argmax(probs) - 1
end
function fit_softmax!(model::SoftmaxRegression, X::Vector{Vector{Float64}},
ys::Vector{Int}; epochs::Int=1000, print_every::Int=200)
n = length(ys)
for epoch in 0:(epochs - 1)
grad_w = [zeros(model.n_features) for _ in 1:model.n_classes]
grad_b = zeros(model.n_classes)
total_loss = 0.0
for i in 1:n
probs = predict_proba_softmax(model, X[i])
for k in 1:model.n_classes
target = ys[i] == (k - 1) ? 1.0 : 0.0
err = probs[k] - target
for j in 1:model.n_features
grad_w[k][j] += err * X[i][j]
end
grad_b[k] += err
end
true_prob = max(probs[ys[i] + 1], 1e-15)
total_loss -= log(true_prob)
end
for k in 1:model.n_classes
for j in 1:model.n_features
model.weights[k][j] -= model.lr * (grad_w[k][j] / n)
end
model.biases[k] -= model.lr * (grad_b[k] / n)
end
if epoch % print_every == 0
@printf(" epoch %4d loss=%.4f\n", epoch, total_loss / n)
end
end
return model
end
function generate_three_class_data(; seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Int[]
centers = [(1.0, 1.0), (5.0, 1.0), (3.0, 5.0)]
for (label, (cx, cy)) in enumerate(centers)
for _ in 1:50
push!(X, Float64[cx + 0.8 * randn(rng), cy + 0.8 * randn(rng)])
push!(ys, label - 1)
end
end
perm = randperm(rng, length(X))
return X[perm], ys[perm]
end
function demo_binary_logistic()
println("=" ^ 60)
println("BINARY LOGISTIC REGRESSION")
println("=" ^ 60)
X, ys = generate_two_class_data()
split = Int(round(0.8 * length(X)))
X_train = X[1:split]
X_test = X[(split + 1):end]
ys_train = ys[1:split]
ys_test = ys[(split + 1):end]
@printf("\nSamples: %d features: 2 classes: {0, 1}\n", length(X))
@printf("Train: %d Test: %d\n", length(X_train), length(X_test))
model = LogisticRegression(2, 0.1)
fit_logistic!(model, X_train, ys_train; epochs=1000, print_every=200)
@printf("\nTrain accuracy: %.4f\n", accuracy(model, X_train, ys_train))
@printf("Test accuracy: %.4f\n", accuracy(model, X_test, ys_test))
@printf("Weights: [%.4f, %.4f]\n", model.weights[1], model.weights[2])
@printf("Bias: %.4f\n", model.bias)
y_pred = [predict_class(model, x) for x in X_test]
metrics = build_metrics(ys_test, y_pred)
print_report(metrics)
return model, X_test, ys_test
end
function demo_decision_boundary(model::LogisticRegression)
println("\n" * "=" ^ 60)
println("DECISION BOUNDARY")
println("=" ^ 60)
w1, w2 = model.weights[1], model.weights[2]
b = model.bias
@printf("\nBoundary: %.4f*x1 + %.4f*x2 + %.4f = 0\n", w1, w2, b)
if abs(w2) > 1e-10
@printf("Solved for x2: x2 = %.4f*x1 + %.4f\n", -w1 / w2, -b / w2)
end
test_points = [Float64[3.0, 3.0], Float64[3.5, 3.5], Float64[4.0, 4.0],
Float64[2.5, 2.5], Float64[5.0, 5.0]]
println("\nProbabilities near the boundary:")
for point in test_points
prob = predict_proba(model, point)
pred = predict_class(model, point)
@printf(" [%.2f, %.2f] -> prob=%.4f class=%d\n",
point[1], point[2], prob, pred)
end
end
function demo_threshold_tuning(model::LogisticRegression,
X_test::Vector{Vector{Float64}}, ys_test::Vector{Int})
println("\n" * "=" ^ 60)
println("THRESHOLD TUNING")
println("=" ^ 60)
println("Default threshold 0.5. Lower = more recall, higher = more precision.\n")
@printf("%10s %10s %10s %10s %10s\n",
"Threshold", "Accuracy", "Precision", "Recall", "F1")
println("-" ^ 54)
for t in (0.3, 0.4, 0.5, 0.6, 0.7)
y_pred_t = [predict_proba(model, x) >= t ? 1 : 0 for x in X_test]
m = build_metrics(ys_test, y_pred_t)
@printf("%10.1f %10.4f %10.4f %10.4f %10.4f\n",
t, metric_accuracy(m), metric_precision(m),
metric_recall(m), metric_f1(m))
end
end
function demo_softmax_regression()
println("\n" * "=" ^ 60)
println("SOFTMAX (MULTI-CLASS) REGRESSION")
println("=" ^ 60)
X, ys = generate_three_class_data()
split = Int(round(0.8 * length(X)))
X_train = X[1:split]
X_test = X[(split + 1):end]
ys_train = ys[1:split]
ys_test = ys[(split + 1):end]
model = SoftmaxRegression(2, 3, 0.1)
fit_softmax!(model, X_train, ys_train; epochs=1000, print_every=200)
train_correct = sum(predict_class_softmax(model, X_train[i]) == ys_train[i]
for i in 1:length(ys_train))
test_correct = sum(predict_class_softmax(model, X_test[i]) == ys_test[i]
for i in 1:length(ys_test))
@printf("\nTrain accuracy: %.4f\n", train_correct / length(ys_train))
@printf("Test accuracy: %.4f\n", test_correct / length(ys_test))
println("\nSample predictions:")
for i in 1:5
probs = predict_proba_softmax(model, X_test[i])
pred = predict_class_softmax(model, X_test[i])
@printf(" true=%d pred=%d probs=[%.3f, %.3f, %.3f]\n",
ys_test[i], pred, probs[1], probs[2], probs[3])
end
end
function demo_why_not_linear()
println("\n" * "=" ^ 60)
println("WHY LINEAR REGRESSION FAILS FOR CLASSIFICATION")
println("=" ^ 60)
hours = Float64[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
pass = Float64[0, 0, 0, 0, 1, 1, 1, 1, 1, 1]
n = length(hours)
x_mean = mean(hours)
y_mean = mean(pass)
num = sum((hours .- x_mean) .* (pass .- y_mean))
den = sum((hours .- x_mean) .^ 2)
w_lin = num / den
b_lin = y_mean - w_lin * x_mean
@printf("\nLinear fit: y = %.4f*x + %.4f\n", w_lin, b_lin)
@printf("%6s %8s %8s %8s\n", "Hours", "Actual", "Linear", "Sigmoid")
for i in 1:n
lin_pred = w_lin * hours[i] + b_lin
sig_pred = sigmoid(3 * (hours[i] - 4.5))
@printf("%6.0f %8.0f %8.3f %8.3f\n", hours[i], pass[i], lin_pred, sig_pred)
end
println("\nLinear regression can output values outside [0, 1].")
println("Sigmoid keeps probabilities inside the valid range.")
end
function main()
model, X_test, ys_test = demo_binary_logistic()
demo_decision_boundary(model)
demo_threshold_tuning(model, X_test, ys_test)
demo_softmax_regression()
demo_why_not_linear()
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
phases/02-ml-fundamentals/05-support-vector-machines/code/main.jl
+403
View File
@@ -0,0 +1,403 @@
# Support vector machines in Julia. Linear SVM trained by stochastic
# sub-gradient descent on hinge loss with L2 regularization (soft margin),
# plus polynomial and RBF kernel functions. Stdlib only. Sources:
# https://docs.julialang.org/en/v1/manual/control-flow/
# https://docs.julialang.org/en/v1/stdlib/Random/
# https://docs.julialang.org/en/v1/manual/arrays/
using Random
using Printf
function dotprod(a::Vector{Float64}, b::Vector{Float64})::Float64
s = 0.0
@inbounds for i in 1:length(a)
s += a[i] * b[i]
end
return s
end
function vec_norm(a::Vector{Float64})::Float64
return sqrt(dotprod(a, a))
end
function linear_kernel(x::Vector{Float64}, z::Vector{Float64})::Float64
return dotprod(x, z)
end
function polynomial_kernel(x::Vector{Float64}, z::Vector{Float64};
degree::Int=3, c::Float64=1.0)::Float64
return (dotprod(x, z) + c) ^ degree
end
function rbf_kernel(x::Vector{Float64}, z::Vector{Float64};
gamma::Float64=0.5)::Float64
diff = x .- z
return exp(-gamma * dotprod(diff, diff))
end
function hinge_loss(X::Vector{Vector{Float64}}, ys::Vector{Int},
w::Vector{Float64}, b::Float64)::Float64
n = length(X)
total = 0.0
for i in 1:n
margin = ys[i] * (dotprod(w, X[i]) + b)
total += max(0.0, 1.0 - margin)
end
return total / n
end
function svm_objective(X::Vector{Vector{Float64}}, ys::Vector{Int},
w::Vector{Float64}, b::Float64, lambda::Float64)::Float64
return 0.5 * lambda * dotprod(w, w) + hinge_loss(X, ys, w, b)
end
mutable struct LinearSVM
w::Vector{Float64}
b::Float64
lr::Float64
lambda::Float64
n_epochs::Int
history::Vector{Tuple{Int, Float64}}
end
LinearSVM(; lr::Float64=0.001, lambda::Float64=0.01, n_epochs::Int=1000) =
LinearSVM(Float64[], 0.0, lr, lambda, n_epochs, Tuple{Int, Float64}[])
function fit_svm!(model::LinearSVM, X::Vector{Vector{Float64}}, ys::Vector{Int};
seed::Int=0)
rng = MersenneTwister(seed)
n_features = length(X[1])
n_samples = length(X)
model.w = zeros(n_features)
model.b = 0.0
empty!(model.history)
for epoch in 0:(model.n_epochs - 1)
indices = randperm(rng, n_samples)
for i in indices
margin = ys[i] * (dotprod(model.w, X[i]) + model.b)
if margin >= 1
for j in 1:n_features
model.w[j] -= model.lr * model.lambda * model.w[j]
end
else
for j in 1:n_features
model.w[j] -= model.lr * (model.lambda * model.w[j] - ys[i] * X[i][j])
end
model.b -= model.lr * (-ys[i])
end
end
if epoch % 100 == 0 || epoch == model.n_epochs - 1
push!(model.history, (epoch, svm_objective(X, ys, model.w, model.b, model.lambda)))
end
end
return model
end
function predict_svm(model::LinearSVM, X::Vector{Vector{Float64}})::Vector{Int}
return [dotprod(model.w, x) + model.b >= 0 ? 1 : -1 for x in X]
end
function decision_function(model::LinearSVM, X::Vector{Vector{Float64}})::Vector{Float64}
return [dotprod(model.w, x) + model.b for x in X]
end
function margin_width(model::LinearSVM)::Float64
n = vec_norm(model.w)
return n == 0 ? 0.0 : 2.0 / n
end
function find_support_vectors(model::LinearSVM, X::Vector{Vector{Float64}},
ys::Vector{Int}; tol::Float64=0.1)::Vector{Int}
svs = Int[]
for i in 1:length(X)
margin = ys[i] * (dotprod(model.w, X[i]) + model.b)
if abs(margin - 1.0) < tol
push!(svs, i)
end
end
return svs
end
function svm_accuracy(y_true::Vector{Int}, y_pred::Vector{Int})::Float64
return sum(y_true .== y_pred) / length(y_true)
end
function generate_linear_data(; n_samples::Int=100, margin::Float64=1.0, seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Int[]
for _ in 1:n_samples
x1 = -3.0 + 6.0 * rand(rng)
x2 = -3.0 + 6.0 * rand(rng)
val = x1 + x2
if val > margin / 2
push!(X, Float64[x1, x2])
push!(ys, 1)
elseif val < -margin / 2
push!(X, Float64[x1, x2])
push!(ys, -1)
end
end
return X, ys
end
function generate_noisy_data(; n_samples::Int=200, noise::Float64=0.5, seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Int[]
for _ in 1:n_samples
x1 = -3.0 + 6.0 * rand(rng)
x2 = -3.0 + 6.0 * rand(rng)
val = x1 - 0.5 * x2 + noise * randn(rng)
push!(X, Float64[x1, x2])
push!(ys, val > 0 ? 1 : -1)
end
return X, ys
end
function generate_circular_data(; n_samples::Int=200, seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Int[]
for _ in 1:n_samples
r = 3.0 * rand(rng)
angle = 2 * pi * rand(rng)
x1 = r * cos(angle) + 0.1 * randn(rng)
x2 = r * sin(angle) + 0.1 * randn(rng)
push!(X, Float64[x1, x2])
push!(ys, r > 1.5 ? 1 : -1)
end
return X, ys
end
function svm_train_test_split(X::Vector{Vector{Float64}}, ys::Vector{Int};
test_ratio::Float64=0.2, seed::Int=42)
rng = MersenneTwister(seed)
indices = randperm(rng, length(X))
split = Int(round(length(X) * (1 - test_ratio)))
train_idx = indices[1:split]
test_idx = indices[(split + 1):end]
return (X[train_idx], ys[train_idx], X[test_idx], ys[test_idx])
end
function demo_hinge_loss()
println("=" ^ 65)
println("HINGE LOSS")
println("=" ^ 65)
println()
margins = [-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 3.0]
@printf(" %10s %12s %14s\n", "y * f(x)", "Hinge loss", "Logistic loss")
println(" " * "-" ^ 10 * " " * "-" ^ 12 * " " * "-" ^ 14)
for m in margins
h = max(0.0, 1.0 - m)
l = log(1 + exp(-m))
@printf(" %10.1f %12.3f %14.3f\n", m, h, l)
end
println()
println(" Hinge loss is exactly zero when y*f(x) >= 1.")
println(" Logistic loss is never exactly zero. Hinge gives sparse models.")
println()
end
function demo_linear_svm()
println("=" ^ 65)
println("LINEAR SVM (SOFT MARGIN)")
println("=" ^ 65)
println()
X, ys = generate_linear_data(n_samples=200, margin=1.0, seed=42)
X_train, ys_train, X_test, ys_test = svm_train_test_split(X, ys)
@printf(" Dataset: %d samples, linearly separable\n", length(X))
@printf(" Train: %d Test: %d\n", length(X_train), length(X_test))
svm = LinearSVM(lr=0.001, lambda=0.01, n_epochs=500)
fit_svm!(svm, X_train, ys_train; seed=1)
train_acc = svm_accuracy(ys_train, predict_svm(svm, X_train))
test_acc = svm_accuracy(ys_test, predict_svm(svm, X_test))
@printf("\n Weights: [%.4f, %.4f]\n", svm.w[1], svm.w[2])
@printf(" Bias: %.4f\n", svm.b)
@printf(" Margin width: %.4f\n", margin_width(svm))
@printf(" Train accuracy: %.4f\n", train_acc)
@printf(" Test accuracy: %.4f\n", test_acc)
svs = find_support_vectors(svm, X_train, ys_train; tol=0.3)
@printf(" Support vectors: %d / %d\n", length(svs), length(X_train))
println()
end
function demo_c_parameter()
println("=" ^ 65)
println("C PARAMETER (REGULARIZATION TRADE-OFF)")
println("=" ^ 65)
println()
X, ys = generate_noisy_data(n_samples=300, noise=0.8, seed=42)
X_train, ys_train, X_test, ys_test = svm_train_test_split(X, ys)
@printf(" %8s %8s %10s %10s %8s %6s\n",
"C", "lambda", "Train Acc", "Test Acc", "Margin", "SVs")
println(" " * "-" ^ 8 * " " * "-" ^ 8 * " " * "-" ^ 10 * " " *
"-" ^ 10 * " " * "-" ^ 8 * " " * "-" ^ 6)
for c in (0.001, 0.01, 0.1, 1.0, 10.0, 100.0)
lam = 1.0 / (c * length(X_train))
svm = LinearSVM(lr=0.001, lambda=lam, n_epochs=500)
fit_svm!(svm, X_train, ys_train; seed=2)
train_acc = svm_accuracy(ys_train, predict_svm(svm, X_train))
test_acc = svm_accuracy(ys_test, predict_svm(svm, X_test))
mw = margin_width(svm)
n_sv = length(find_support_vectors(svm, X_train, ys_train; tol=0.3))
@printf(" %8.3f %8.5f %10.4f %10.4f %8.4f %6d\n",
c, lam, train_acc, test_acc, mw, n_sv)
end
println()
println(" Small C (large lambda): wide margin, more slack, better generalization.")
println(" Large C (small lambda): narrow margin, fewer slack, risk of overfit.")
println()
end
function demo_kernels()
println("=" ^ 65)
println("KERNEL FUNCTIONS")
println("=" ^ 65)
println()
x = Float64[1.0, 0.0]
cases = [
("same direction", Float64[2.0, 0.0]),
("perpendicular", Float64[0.0, 1.0]),
("close", Float64[1.1, 0.1]),
("far same dir", Float64[5.0, 0.0]),
("opposite", Float64[-1.0, 0.0]),
]
@printf(" Reference: %s\n", x)
println()
@printf(" %-20s %8s %10s %10s %10s\n",
"Point", "Linear", "Poly(d=2)", "Poly(d=3)", "RBF(g=0.5)")
println(" " * "-" ^ 20 * " " * "-" ^ 8 * " " * "-" ^ 10 * " " *
"-" ^ 10 * " " * "-" ^ 10)
for (name, z) in cases
k_l = linear_kernel(x, z)
k_p2 = polynomial_kernel(x, z; degree=2)
k_p3 = polynomial_kernel(x, z; degree=3)
k_rbf = rbf_kernel(x, z; gamma=0.5)
@printf(" %-20s %8.3f %10.3f %10.3f %10.4f\n",
name, k_l, k_p2, k_p3, k_rbf)
end
println()
println(" Linear kernel: raw dot product. RBF: locality-based.")
println()
end
function demo_linear_vs_nonlinear()
println("=" ^ 65)
println("LINEAR SVM vs POLYNOMIAL FEATURE MAP")
println("=" ^ 65)
println()
X, ys = generate_circular_data(n_samples=200, seed=42)
X_train, ys_train, X_test, ys_test = svm_train_test_split(X, ys)
svm = LinearSVM(lr=0.001, lambda=0.01, n_epochs=500)
fit_svm!(svm, X_train, ys_train; seed=3)
train_acc = svm_accuracy(ys_train, predict_svm(svm, X_train))
test_acc = svm_accuracy(ys_test, predict_svm(svm, X_test))
@printf(" Plain linear SVM on circular data: train=%.4f test=%.4f\n",
train_acc, test_acc)
println()
function augment(X)
return [Float64[x[1], x[2], x[1] ^ 2, x[2] ^ 2, x[1] * x[2]] for x in X]
end
X_train_aug = augment(X_train)
X_test_aug = augment(X_test)
svm_aug = LinearSVM(lr=0.0005, lambda=0.01, n_epochs=1000)
fit_svm!(svm_aug, X_train_aug, ys_train; seed=4)
train_aug = svm_accuracy(ys_train, predict_svm(svm_aug, X_train_aug))
test_aug = svm_accuracy(ys_test, predict_svm(svm_aug, X_test_aug))
println(" After polynomial feature map (x1, x2, x1^2, x2^2, x1*x2):")
@printf(" Linear SVM on augmented features: train=%.4f test=%.4f\n",
train_aug, test_aug)
println()
println(" The kernel trick performs this feature map implicitly.")
println()
end
function demo_support_vectors()
println("=" ^ 65)
println("SUPPORT VECTORS")
println("=" ^ 65)
println()
X, ys = generate_linear_data(n_samples=200, margin=1.5, seed=42)
X_train, ys_train, _, _ = svm_train_test_split(X, ys)
svm = LinearSVM(lr=0.001, lambda=0.01, n_epochs=1000)
fit_svm!(svm, X_train, ys_train; seed=5)
margins = [(i, ys_train[i] * (dotprod(svm.w, X_train[i]) + svm.b))
for i in 1:length(X_train)]
sort!(margins; by=t -> t[2])
@printf(" Trained on %d points.\n", length(X_train))
@printf(" Weights: [%.4f, %.4f] bias: %.4f\n", svm.w[1], svm.w[2], svm.b)
println()
println(" Points sorted by margin (y * f(x)):")
@printf(" %6s %4s %8s %s\n", "Index", "y", "Margin", "Role")
println(" " * "-" ^ 6 * " " * "-" ^ 4 * " " * "-" ^ 8 * " " * "-" ^ 20)
for (idx, m) in margins[1:8]
role = m < 0 ? "MISCLASSIFIED" :
m < 1.0 ? "inside margin" :
m < 1.2 ? "SUPPORT VECTOR" :
"safely classified"
@printf(" %6d %4d %8.4f %s\n", idx, ys_train[idx], m, role)
end
println(" ...")
for (idx, m) in margins[(end - 2):end]
@printf(" %6d %4d %8.4f safely classified\n", idx, ys_train[idx], m)
end
n_sv = sum(1 for (_, m) in margins if 0.7 < m < 1.3)
n_safe = sum(1 for (_, m) in margins if m >= 1.3)
n_inside = sum(1 for (_, m) in margins if 0 < m < 0.7)
println()
@printf(" Support vectors (margin ~ 1.0): %d\n", n_sv)
@printf(" Safely classified (margin >> 1): %d\n", n_safe)
@printf(" Inside margin (0 < margin < 1): %d\n", n_inside)
println()
end
function main()
demo_hinge_loss()
demo_linear_svm()
demo_c_parameter()
demo_kernels()
demo_linear_vs_nonlinear()
demo_support_vectors()
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
phases/02-ml-fundamentals/09-model-evaluation/code/main.jl
+380
View File
@@ -0,0 +1,380 @@
# Model evaluation in Julia. Train/val/test split, k-fold + stratified k-fold
# cross validation, classification metrics (accuracy, precision, recall, F1,
# ROC, AUC), and regression metrics (MSE, RMSE, MAE, R^2). Stdlib only. Sources:
# https://docs.julialang.org/en/v1/stdlib/Random/
# https://docs.julialang.org/en/v1/stdlib/Statistics/
# https://docs.julialang.org/en/v1/manual/functions/
using Random
using Statistics
using Printf
function train_val_test_split(X::Vector{Vector{Float64}}, ys::Vector{Int};
train_ratio::Float64=0.6, val_ratio::Float64=0.2,
seed::Int=42)
rng = MersenneTwister(seed)
n = length(X)
indices = randperm(rng, n)
train_end = Int(round(n * train_ratio))
val_end = Int(round(n * (train_ratio + val_ratio)))
train_idx = indices[1:train_end]
val_idx = indices[(train_end + 1):val_end]
test_idx = indices[(val_end + 1):end]
return (X[train_idx], ys[train_idx],
X[val_idx], ys[val_idx],
X[test_idx], ys[test_idx])
end
function kfold_split(n::Int; k::Int=5, seed::Int=42)
rng = MersenneTwister(seed)
indices = randperm(rng, n)
fold_size = n ÷ k
folds = Vector{Tuple{Vector{Int}, Vector{Int}}}()
for i in 1:k
s = (i - 1) * fold_size + 1
e = i < k ? i * fold_size : n
val_idx = indices[s:e]
train_idx = vcat(indices[1:(s - 1)], indices[(e + 1):end])
push!(folds, (train_idx, val_idx))
end
return folds
end
function stratified_kfold_split(ys::Vector{Int}; k::Int=5, seed::Int=42)
rng = MersenneTwister(seed)
class_indices = Dict{Int, Vector{Int}}()
for (i, label) in enumerate(ys)
push!(get!(class_indices, label, Int[]), i)
end
for label in keys(class_indices)
shuffle!(rng, class_indices[label])
end
train_lists = [Int[] for _ in 1:k]
val_lists = [Int[] for _ in 1:k]
for indices in values(class_indices)
fold_size = length(indices) ÷ k
for i in 1:k
s = (i - 1) * fold_size + 1
e = i < k ? i * fold_size : length(indices)
val_part = indices[s:e]
train_part = vcat(indices[1:(s - 1)], indices[(e + 1):end])
append!(val_lists[i], val_part)
append!(train_lists[i], train_part)
end
end
return [(train_lists[i], val_lists[i]) for i in 1:k]
end
function confusion_matrix(y_true::Vector{Int}, y_pred::Vector{Int})
tp = sum(1 for i in 1:length(y_true) if y_true[i] == 1 && y_pred[i] == 1)
tn = sum(1 for i in 1:length(y_true) if y_true[i] == 0 && y_pred[i] == 0)
fp = sum(1 for i in 1:length(y_true) if y_true[i] == 0 && y_pred[i] == 1)
fn = sum(1 for i in 1:length(y_true) if y_true[i] == 1 && y_pred[i] == 0)
return tp, tn, fp, fn
end
function accuracy(y_true::Vector{Int}, y_pred::Vector{Int})
tp, tn, fp, fn = confusion_matrix(y_true, y_pred)
total = tp + tn + fp + fn
return total > 0 ? (tp + tn) / total : 0.0
end
function precision_score(y_true::Vector{Int}, y_pred::Vector{Int})
tp, _, fp, _ = confusion_matrix(y_true, y_pred)
return (tp + fp) > 0 ? tp / (tp + fp) : 0.0
end
function recall_score(y_true::Vector{Int}, y_pred::Vector{Int})
tp, _, _, fn = confusion_matrix(y_true, y_pred)
return (tp + fn) > 0 ? tp / (tp + fn) : 0.0
end
function f1_score(y_true::Vector{Int}, y_pred::Vector{Int})
p = precision_score(y_true, y_pred)
r = recall_score(y_true, y_pred)
return (p + r) > 0 ? 2 * p * r / (p + r) : 0.0
end
function roc_curve(y_true::Vector{Int}, y_scores::Vector{Float64})
thresholds = sort(unique(y_scores); rev=true)
tpr_list = Float64[]
fpr_list = Float64[]
total_pos = sum(y_true)
total_neg = length(y_true) - total_pos
for t in thresholds
y_pred = [s >= t ? 1 : 0 for s in y_scores]
tp = sum(1 for i in 1:length(y_true) if y_true[i] == 1 && y_pred[i] == 1)
fp = sum(1 for i in 1:length(y_true) if y_true[i] == 0 && y_pred[i] == 1)
push!(tpr_list, total_pos > 0 ? tp / total_pos : 0.0)
push!(fpr_list, total_neg > 0 ? fp / total_neg : 0.0)
end
return fpr_list, tpr_list, thresholds
end
function auc_roc(y_true::Vector{Int}, y_scores::Vector{Float64})
fpr, tpr, _ = roc_curve(y_true, y_scores)
pairs = sort(collect(zip(fpr, tpr)); by=first)
fpr_sorted = [p[1] for p in pairs]
tpr_sorted = [p[2] for p in pairs]
area = 0.0
for i in 2:length(fpr_sorted)
width = fpr_sorted[i] - fpr_sorted[i - 1]
height = (tpr_sorted[i] + tpr_sorted[i - 1]) / 2
area += width * height
end
return area
end
function mse(y_true::Vector{Float64}, y_pred::Vector{Float64})
n = length(y_true)
return sum((y_true .- y_pred) .^ 2) / n
end
function rmse(y_true::Vector{Float64}, y_pred::Vector{Float64})
return sqrt(mse(y_true, y_pred))
end
function mae(y_true::Vector{Float64}, y_pred::Vector{Float64})
n = length(y_true)
return sum(abs.(y_true .- y_pred)) / n
end
function r_squared(y_true::Vector{Float64}, y_pred::Vector{Float64})
mean_y = mean(y_true)
ss_res = sum((y_true .- y_pred) .^ 2)
ss_tot = sum((y_true .- mean_y) .^ 2)
return ss_tot == 0 ? 0.0 : 1.0 - ss_res / ss_tot
end
function sigmoid(z::Float64)
z_clip = clamp(z, -500.0, 500.0)
return 1.0 / (1.0 + exp(-z_clip))
end
mutable struct SimpleLogistic
weights::Vector{Float64}
bias::Float64
lr::Float64
epochs::Int
end
SimpleLogistic(lr::Float64, epochs::Int) = SimpleLogistic(Float64[], 0.0, lr, epochs)
function fit_simple!(model::SimpleLogistic, X::Vector{Vector{Float64}}, ys::Vector{Int})
n_features = length(X[1])
model.weights = zeros(n_features)
model.bias = 0.0
for _ in 1:model.epochs
for i in 1:length(X)
z = sum(model.weights .* X[i]) + model.bias
p = sigmoid(z)
err = ys[i] - p
for j in 1:n_features
model.weights[j] += model.lr * err * X[i][j]
end
model.bias += model.lr * err
end
end
return model
end
function predict_proba_simple(model::SimpleLogistic, x::Vector{Float64})
return sigmoid(sum(model.weights .* x) + model.bias)
end
predict_simple(model::SimpleLogistic, x::Vector{Float64}) =
predict_proba_simple(model, x) >= 0.5 ? 1 : 0
function cross_validate(X::Vector{Vector{Float64}}, ys::Vector{Int},
model_fn::Function; k::Int=5,
metric_fn::Function=accuracy, stratified::Bool=false)
n = length(X)
folds = stratified ? stratified_kfold_split(ys; k=k) : kfold_split(n; k=k)
scores = Float64[]
for (train_idx, val_idx) in folds
X_train = X[train_idx]
ys_train = ys[train_idx]
X_val = X[val_idx]
ys_val = ys[val_idx]
model = model_fn()
fit_simple!(model, X_train, ys_train)
preds = [predict_simple(model, x) for x in X_val]
push!(scores, metric_fn(ys_val, preds))
end
return scores
end
function make_classification_data(n::Int=300; seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Int[]
for _ in 1:n
x1 = randn(rng)
x2 = randn(rng)
label = (x1 + x2 + 0.5 * randn(rng)) > 0 ? 1 : 0
push!(X, Float64[x1, x2])
push!(ys, label)
end
return X, ys
end
function make_regression_data(n::Int=200; seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Float64[]
for _ in 1:n
x1 = 10.0 * rand(rng)
x2 = 5.0 * rand(rng)
target = 3 * x1 + 2 * x2 + 2 * randn(rng)
push!(X, Float64[x1, x2])
push!(ys, target)
end
return X, ys
end
function make_imbalanced_data(n::Int=300; minority_ratio::Float64=0.05, seed::Int=42)
rng = MersenneTwister(seed)
X = Vector{Vector{Float64}}()
ys = Int[]
for _ in 1:n
if rand(rng) < minority_ratio
push!(X, Float64[3.0 + 0.5 * randn(rng), 3.0 + 0.5 * randn(rng)])
push!(ys, 1)
else
push!(X, Float64[randn(rng), randn(rng)])
push!(ys, 0)
end
end
return X, ys
end
function demo_split_and_metrics()
println("=" ^ 60)
println("TRAIN / VAL / TEST SPLIT + METRICS")
println("=" ^ 60)
X, ys = make_classification_data(300)
X_train, ys_train, X_val, ys_val, X_test, ys_test = train_val_test_split(X, ys)
@printf(" Train: %d Val: %d Test: %d\n",
length(X_train), length(X_val), length(X_test))
@printf(" Train positive ratio: %.3f\n", sum(ys_train) / length(ys_train))
@printf(" Val positive ratio: %.3f\n", sum(ys_val) / length(ys_val))
model = SimpleLogistic(0.1, 200)
fit_simple!(model, X_train, ys_train)
println("\n--- Classification metrics ---")
y_pred = [predict_simple(model, x) for x in X_test]
tp, tn, fp, fn = confusion_matrix(ys_test, y_pred)
@printf(" Confusion: TP=%d TN=%d FP=%d FN=%d\n", tp, tn, fp, fn)
@printf(" Accuracy: %.4f\n", accuracy(ys_test, y_pred))
@printf(" Precision: %.4f\n", precision_score(ys_test, y_pred))
@printf(" Recall: %.4f\n", recall_score(ys_test, y_pred))
@printf(" F1: %.4f\n", f1_score(ys_test, y_pred))
y_scores = [predict_proba_simple(model, x) for x in X_test]
@printf(" AUC-ROC: %.4f\n", auc_roc(ys_test, y_scores))
end
function demo_cross_validation()
println("\n" * "=" ^ 60)
println("K-FOLD CROSS VALIDATION")
println("=" ^ 60)
X, ys = make_classification_data(300)
scores = cross_validate(X, ys, () -> SimpleLogistic(0.1, 200);
k=5, metric_fn=accuracy)
m = mean(scores)
s = std(scores; corrected=false)
println("\nPlain k=5:")
@printf(" Fold scores: [%s]\n",
join([@sprintf("%.4f", v) for v in scores], ", "))
@printf(" Mean: %.4f (+/- %.4f)\n", m, s)
strat = cross_validate(X, ys, () -> SimpleLogistic(0.1, 200);
k=5, metric_fn=accuracy, stratified=true)
sm = mean(strat)
ss = std(strat; corrected=false)
println("\nStratified k=5:")
@printf(" Fold scores: [%s]\n",
join([@sprintf("%.4f", v) for v in strat], ", "))
@printf(" Mean: %.4f (+/- %.4f)\n", sm, ss)
end
function demo_imbalanced()
println("\n" * "=" ^ 60)
println("IMBALANCED DATA: WHY ACCURACY LIES")
println("=" ^ 60)
X, ys = make_imbalanced_data(300; minority_ratio=0.05)
positives = sum(ys)
@printf("\n Class distribution: %d positive, %d negative (%.1f%% positive)\n",
positives, length(ys) - positives, 100 * positives / length(ys))
baseline = zeros(Int, length(ys))
println("\n Always-negative baseline:")
@printf(" Accuracy: %.4f\n", accuracy(ys, baseline))
@printf(" Precision: %.4f\n", precision_score(ys, baseline))
@printf(" Recall: %.4f\n", recall_score(ys, baseline))
@printf(" F1: %.4f\n", f1_score(ys, baseline))
println(" Accuracy lies; precision and recall expose the failure.")
end
function demo_regression_metrics()
println("\n" * "=" ^ 60)
println("REGRESSION METRICS")
println("=" ^ 60)
X, ys = make_regression_data(200)
n_train = Int(round(0.8 * length(X)))
y_pred = Float64[]
y_true = ys[(n_train + 1):end]
for i in (n_train + 1):length(ys)
push!(y_pred, ys[i] + randn() * 0.5)
end
@printf(" MSE: %.4f\n", mse(y_true, y_pred))
@printf(" RMSE: %.4f\n", rmse(y_true, y_pred))
@printf(" MAE: %.4f\n", mae(y_true, y_pred))
@printf(" R^2: %.4f\n", r_squared(y_true, y_pred))
mean_baseline = fill(mean(y_true), length(y_true))
println("\n Predict-the-mean baseline:")
@printf(" MSE: %.4f\n", mse(y_true, mean_baseline))
@printf(" R^2: %.4f\n", r_squared(y_true, mean_baseline))
end
function main()
demo_split_and_metrics()
demo_cross_validation()
demo_imbalanced()
demo_regression_metrics()
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
phases/07-transformers-deep-dive/01-why-transformers/code/main.jl
+130
View File
@@ -0,0 +1,130 @@
# Why transformers in Julia. Contrasts RNN-style serial recurrence with
# attention-style parallel reduction, and verifies that Hillis-Steele
# parallel prefix scan matches the serial scan. Stdlib only. Sources:
# https://docs.julialang.org/en/v1/manual/control-flow/
# https://docs.julialang.org/en/v1/stdlib/Base/
# https://en.wikipedia.org/wiki/Prefix_sum
using Printf
function rnn_style(xs::Vector{Float64}; decay::Float64=0.9)::Float64
h = 0.0
for x in xs
h = decay * h + x
end
return h
end
function attention_style(xs::Vector{Float64})::Float64
isempty(xs) && throw(ArgumentError("xs must be non-empty"))
return sum(xs) / length(xs)
end
function serial_scan(xs::Vector{Float64})::Vector{Float64}
out = similar(xs)
acc = 0.0
@inbounds for i in 1:length(xs)
acc += xs[i]
out[i] = acc
end
return out
end
function parallel_scan(xs::Vector{Float64})::Vector{Float64}
out = copy(xs)
n = length(out)
step = 1
while step < n
new_out = copy(out)
for i in (step + 1):n
new_out[i] = out[i] + out[i - step]
end
out = new_out
step *= 2
end
return out
end
function benchmark_pair(n::Int; reps::Int=3)
n > 0 || throw(ArgumentError("n must be > 0"))
xs = [0.001 * mod(i, 17) for i in 0:(n - 1)]
best_rnn = Inf
for _ in 1:reps
t0 = time_ns()
rnn_style(xs)
best_rnn = min(best_rnn, (time_ns() - t0) / 1e9)
end
best_attn = Inf
for _ in 1:reps
t0 = time_ns()
attention_style(xs)
best_attn = min(best_attn, (time_ns() - t0) / 1e9)
end
return best_rnn, best_attn
end
function depth_counts(n::Int)
n > 0 || throw(ArgumentError("n must be > 0"))
rnn_depth = n
attn_depth = max(1, Int(ceil(log2(n))))
return rnn_depth, attn_depth
end
function demo_depth_table()
println("=== serial-depth comparison ===")
@printf("%8s %12s %12s %16s\n", "N", "rnn depth", "attn depth", "speedup (ops)")
for n in (64, 512, 4096, 32768, 262144)
rd, ad = depth_counts(n)
@printf("%8d %12d %12d %15.0fx\n", n, rd, ad, rd / ad)
end
println()
end
function demo_wallclock()
println("=== wall-clock on this machine (pure Julia) ===")
@printf("%8s %10s %10s %8s\n", "N", "rnn (ms)", "attn (ms)", "ratio")
for n in (1_000, 10_000, 100_000, 1_000_000)
rnn_t, attn_t = benchmark_pair(n)
ratio = attn_t > 0 ? rnn_t / attn_t : Inf
@printf("%8d %10.2f %10.2f %7.1fx\n",
n, rnn_t * 1000, attn_t * 1000, ratio)
end
println()
end
function demo_scan_equivalence()
println("=== prefix-sum equivalence check ===")
xs = Float64.(0:15)
ser = serial_scan(xs)
par = parallel_scan(xs)
mismatches = sum(1 for i in 1:length(xs) if abs(ser[i] - par[i]) > 1e-9)
@printf("length: %d mismatches between serial and parallel scan: %d\n",
length(xs), mismatches)
@printf("last value (serial): %.4f\n", ser[end])
@printf("last value (parallel): %.4f\n", par[end])
println()
end
function main()
demo_depth_table()
demo_wallclock()
demo_scan_equivalence()
println("takeaway: attention parallelizes the reduction; depth O(log N) on a")
println("real GPU kernel. Memory cost is O(N^2) for full attention; that")
println("trade-off is what later lessons unpack.")
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
phases/07-transformers-deep-dive/02-self-attention-from-scratch/code/main.jl
+209
View File
@@ -0,0 +1,209 @@
# Self-attention from scratch in Julia. Scaled dot-product attention,
# numerically-stable row-wise softmax, single-head and multi-head
# self-attention. Stdlib only. Sources:
# https://arxiv.org/abs/1706.03762
# https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/
# https://docs.julialang.org/en/v1/stdlib/Random/
using Random
using LinearAlgebra
using Printf
function softmax_rows(M::Matrix{Float64})::Matrix{Float64}
out = similar(M)
for i in 1:size(M, 1)
row = M[i, :]
m = maximum(row)
e = exp.(row .- m)
s = sum(e)
out[i, :] = e ./ s
end
return out
end
function scaled_dot_product_attention(Q::Matrix{Float64}, K::Matrix{Float64},
V::Matrix{Float64})
dk = size(Q, 2)
scores = (Q * transpose(K)) ./ sqrt(dk)
weights = softmax_rows(scores)
output = weights * V
return output, weights
end
struct SelfAttention
Wq::Matrix{Float64}
Wk::Matrix{Float64}
Wv::Matrix{Float64}
dk::Int
end
function SelfAttention(d_model::Int, dk::Int, dv::Int; seed::Int=42)
rng = MersenneTwister(seed)
scale_qk = sqrt(2.0 / (d_model + dk))
scale_v = sqrt(2.0 / (d_model + dv))
Wq = scale_qk .* randn(rng, d_model, dk)
Wk = scale_qk .* randn(rng, d_model, dk)
Wv = scale_v .* randn(rng, d_model, dv)
return SelfAttention(Wq, Wk, Wv, dk)
end
function forward(attn::SelfAttention, X::Matrix{Float64})
Q = X * attn.Wq
K = X * attn.Wk
V = X * attn.Wv
return scaled_dot_product_attention(Q, K, V)
end
struct MultiHeadSelfAttention
heads::Vector{SelfAttention}
Wo::Matrix{Float64}
n_heads::Int
end
function MultiHeadSelfAttention(d_model::Int, n_heads::Int; seed::Int=42)
@assert n_heads > 0 "n_heads must be > 0"
@assert d_model > 0 "d_model must be > 0"
@assert d_model % n_heads == 0 "d_model must be divisible by n_heads"
dk = d_model ÷ n_heads
dv = d_model ÷ n_heads
heads = [SelfAttention(d_model, dk, dv; seed=seed + i) for i in 1:n_heads]
rng = MersenneTwister(seed + n_heads + 1)
scale = sqrt(2.0 / (d_model + d_model))
Wo = scale .* randn(rng, n_heads * dv, d_model)
return MultiHeadSelfAttention(heads, Wo, n_heads)
end
function forward(mha::MultiHeadSelfAttention, X::Matrix{Float64})
head_outputs = Matrix{Float64}[]
weights_per_head = Matrix{Float64}[]
for head in mha.heads
out, w = forward(head, X)
push!(head_outputs, out)
push!(weights_per_head, w)
end
concat = hcat(head_outputs...)
return concat * mha.Wo, weights_per_head
end
function print_attention_matrix(weights::Matrix{Float64}, tokens::Vector{String})
print("\n ")
for token in tokens
@printf("%6s", token)
end
println()
for i in 1:length(tokens)
@printf("%6s", tokens[i])
for j in 1:length(tokens)
@printf("%6.3f", weights[i, j])
end
println()
end
end
function ascii_heatmap(weights::Matrix{Float64}, tokens::Vector{String};
chars::String=" .:-=+*#%@")
print("\n ")
for t in tokens
@printf("%6s", t)
end
println()
w_max = maximum(weights)
for i in 1:length(tokens)
@printf("%6s", tokens[i])
for j in 1:length(tokens)
level = Int(floor(weights[i, j] * (length(chars) - 1) / w_max))
level = min(level, length(chars) - 1)
ch = chars[level + 1]
@printf(" %s ", ch)
end
println()
end
end
function demo_softmax_stability()
println("\n" * "=" ^ 60)
println("SOFTMAX NUMERIC STABILITY")
println("=" ^ 60)
logits = reshape([2.0, 1.0, 0.1], 1, 3)
probs = softmax_rows(logits)
@printf("\nLogits: [%s]\n", join([@sprintf("%.4f", v) for v in logits], ", "))
@printf("Softmax: [%s]\n", join([@sprintf("%.4f", v) for v in probs], ", "))
@printf("Sum: %.4f\n", sum(probs))
big_logits = reshape([100.0, 200.0, 300.0], 1, 3)
big_probs = softmax_rows(big_logits)
@printf("\nLarge logits: [%s]\n",
join([@sprintf("%.1f", v) for v in big_logits], ", "))
@printf("Softmax: [%s]\n",
join([@sprintf("%.4f", v) for v in big_probs], ", "))
@printf("Sum: %.4f\n", sum(big_probs))
println("(no overflow because we subtract the row maximum before exp)")
end
function demo_self_attention()
println("=" ^ 60)
println("SELF-ATTENTION FROM SCRATCH")
println("=" ^ 60)
tokens = ["The", "cat", "sat", "on", "the", "mat"]
n_tokens = length(tokens)
d_model = 16
dk = 8
dv = 8
rng = MersenneTwister(42)
X = randn(rng, n_tokens, d_model)
@printf("\nSentence: %s\n", join(tokens, " "))
@printf("Tokens: %d d_model: %d dk: %d dv: %d\n", n_tokens, d_model, dk, dv)
@printf("Input shape: (%d, %d)\n", size(X, 1), size(X, 2))
attn = SelfAttention(d_model, dk, dv; seed=42)
output, weights = forward(attn, X)
@printf("\nOutput shape: (%d, %d)\n", size(output, 1), size(output, 2))
println("\nAttention weights:")
print_attention_matrix(weights, tokens)
println("\nASCII heatmap (denser char = higher attention):")
ascii_heatmap(weights, tokens)
return tokens, X, d_model
end
function demo_multi_head(tokens::Vector{String}, X::Matrix{Float64}, d_model::Int)
println("\n" * "=" ^ 60)
println("MULTI-HEAD SELF-ATTENTION")
println("=" ^ 60)
n_heads = 2
mha = MultiHeadSelfAttention(d_model, n_heads; seed=42)
out, head_weights = forward(mha, X)
@printf("\nHeads: %d Output shape: (%d, %d)\n",
n_heads, size(out, 1), size(out, 2))
for (h, w) in enumerate(head_weights)
@printf("\nHead %d attention weights:\n", h)
print_attention_matrix(w, tokens)
end
end
function main()
tokens, X, d_model = demo_self_attention()
demo_multi_head(tokens, X, d_model)
demo_softmax_stability()
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
phases/07-transformers-deep-dive/04-positional-encoding/code/main.jl
+150
View File
@@ -0,0 +1,150 @@
# Positional encoding in Julia. Sinusoidal absolute positions, rotary
# positional embedding (RoPE), and ALiBi bias matrix. Verifies that
# RoPE dot products depend only on relative distance. Stdlib only. Sources:
# https://arxiv.org/abs/2104.09864
# https://arxiv.org/abs/2108.12409
# https://docs.julialang.org/en/v1/manual/mathematical-operations/
using Random
using Printf
function sinusoidal_pe(n::Int, d::Int; base::Float64=10000.0)::Matrix{Float64}
n > 0 || throw(ArgumentError("n must be > 0"))
d > 0 || throw(ArgumentError("d must be > 0"))
iseven(d) || throw(ArgumentError("d must be even for sinusoidal sin/cos pairs"))
pe = zeros(n, d)
for pos in 0:(n - 1)
for i in 0:(d ÷ 2 - 1)
theta = pos / (base ^ (2 * i / d))
pe[pos + 1, 2 * i + 1] = sin(theta)
pe[pos + 1, 2 * i + 2] = cos(theta)
end
end
return pe
end
function apply_rope(x::Vector{Float64}, pos::Int; base::Float64=10000.0)::Vector{Float64}
d = length(x)
iseven(d) || throw(ArgumentError("RoPE requires an even embedding dimension"))
out = copy(x)
for i in 0:(d ÷ 2 - 1)
theta = pos / (base ^ (2 * i / d))
c = cos(theta)
s = sin(theta)
a = x[2 * i + 1]
b = x[2 * i + 2]
out[2 * i + 1] = a * c - b * s
out[2 * i + 2] = a * s + b * c
end
return out
end
function dotprod(a::Vector{Float64}, b::Vector{Float64})::Float64
return sum(a .* b)
end
function alibi_slopes(n_heads::Int)::Vector{Float64}
n_heads > 0 || throw(ArgumentError("n_heads must be > 0"))
return [2.0 ^ (-8.0 * (h) / n_heads) for h in 1:n_heads]
end
function alibi_bias(n_heads::Int, seq_len::Int; causal::Bool=true)
slopes = alibi_slopes(n_heads)
out = Vector{Matrix{Float64}}()
for m in slopes
bias = fill(0.0, seq_len, seq_len)
for i in 1:seq_len
for j in 1:seq_len
if causal && j > i
bias[i, j] = -Inf
else
bias[i, j] = -m * abs(i - j)
end
end
end
push!(out, bias)
end
return out
end
function demo_sinusoidal()
println("=== sinusoidal positional encoding ===")
pe = sinusoidal_pe(8, 8)
println("first 4 positions, first 4 dims:")
for pos in 1:4
row_str = join([@sprintf("%+.3f", pe[pos, j]) for j in 1:4], " ")
@printf(" pos=%d: %s\n", pos - 1, row_str)
end
println()
end
function demo_rope_relative()
println("=== RoPE: dot product depends only on relative distance ===")
rng = MersenneTwister(0)
d = 16
q = randn(rng, d)
k = randn(rng, d)
pairs = [(3, 5), (7, 9), (100, 102), (1024, 1026)]
@printf("%6s %6s %4s %18s\n", "pos_q", "pos_k", "gap", "<q_rot, k_rot>")
for (pq, pk) in pairs
q_rot = apply_rope(q, pq)
k_rot = apply_rope(k, pk)
d_prod = dotprod(q_rot, k_rot)
@printf("%6d %6d %4d %18.6f\n", pq, pk, pk - pq, d_prod)
end
println("All rows with gap=2 should produce matching dot products.")
println()
end
function demo_rope_base_scaling()
println("=== RoPE base scaling (NTK-aware for long context) ===")
rng = MersenneTwister(1)
d = 8
q = randn(rng, d)
k = randn(rng, d)
for base in (10000.0, 100000.0, 1_000_000.0)
q_rot = apply_rope(q, 4096; base=base)
k_rot = apply_rope(k, 4098; base=base)
@printf(" base=%8d score=%+.6f\n", Int(base), dotprod(q_rot, k_rot))
end
println("Larger base = slower rotation = longer context without phase wrap.")
println()
end
function demo_alibi()
println("=== ALiBi bias matrix ===")
n_heads = 4
slopes = alibi_slopes(n_heads)
@printf("Slopes for %d heads: %s\n", n_heads,
join([@sprintf("%.4f", s) for s in slopes], ", "))
bias = alibi_bias(n_heads, 6; causal=false)
println("Head 1 bias (closer tokens get smaller penalty):")
for row in eachrow(bias[1])
println(" " * join([@sprintf("%+6.2f", v) for v in row], " "))
end
println()
end
function main()
demo_sinusoidal()
demo_rope_relative()
demo_rope_base_scaling()
demo_alibi()
println("takeaway: RoPE encodes relative position inside the dot product;")
println("ALiBi skips embeddings entirely. Sinusoidal is now a footnote.")
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end
phases/07-transformers-deep-dive/05-full-transformer/code/main.jl
+353
View File
@@ -0,0 +1,353 @@
# Full transformer in Julia: encoder + decoder blocks (pre-norm), multi-head
# attention, SwiGLU FFN, LayerNorm and RMSNorm forward + backward gradient
# check against finite differences. Stdlib only. Sources:
# https://arxiv.org/abs/1706.03762
# https://arxiv.org/abs/1910.07467
# https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/
using Random
using LinearAlgebra
using Printf
function randn_matrix(rng::AbstractRNG, rows::Int, cols::Int;
scale::Union{Nothing, Float64}=nothing)::Matrix{Float64}
s = scale === nothing ? sqrt(2.0 / (rows + cols)) : scale
return s .* randn(rng, rows, cols)
end
function softmax_rows(M::Matrix{Float64};
mask::Union{Nothing, Matrix{Bool}}=nothing)::Matrix{Float64}
out = similar(M)
rows, cols = size(M)
for i in 1:rows
row = M[i, :]
if mask !== nothing
row = [mask[i, j] ? -Inf : row[j] for j in 1:cols]
end
finite = filter(isfinite, row)
m = isempty(finite) ? 0.0 : maximum(finite)
e = [isfinite(v) ? exp(v - m) : 0.0 for v in row]
s = sum(e)
out[i, :] = s > 0 ? e ./ s : zeros(cols)
end
return out
end
function layer_norm(X::Matrix{Float64}; eps::Float64=1e-5)::Matrix{Float64}
out = similar(X)
for i in 1:size(X, 1)
row = X[i, :]
mu = sum(row) / length(row)
var = sum((row .- mu) .^ 2) / length(row)
denom = sqrt(var + eps)
out[i, :] = (row .- mu) ./ denom
end
return out
end
function rms_norm(X::Matrix{Float64}; eps::Float64=1e-6)::Matrix{Float64}
out = similar(X)
for i in 1:size(X, 1)
row = X[i, :]
rms = sqrt(sum(row .* row) / length(row) + eps)
out[i, :] = row ./ rms
end
return out
end
function layer_norm_backward(X::Matrix{Float64}, dY::Matrix{Float64};
eps::Float64=1e-5)::Matrix{Float64}
rows, d = size(X)
dX = similar(X)
for i in 1:rows
x = X[i, :]
dy = dY[i, :]
mu = sum(x) / d
xc = x .- mu
var = sum(xc .* xc) / d
denom = sqrt(var + eps)
xhat = xc ./ denom
dxhat = dy
dvar = sum(dxhat .* xc) * -0.5 * (var + eps) ^ (-1.5)
dmu = sum(dxhat .* (-1.0 ./ denom)) + dvar * sum(-2.0 .* xc) / d
dX[i, :] = dxhat ./ denom .+ dvar .* 2.0 .* xc ./ d .+ dmu / d
end
return dX
end
function rms_norm_backward(X::Matrix{Float64}, dY::Matrix{Float64};
eps::Float64=1e-6)::Matrix{Float64}
rows, d = size(X)
dX = similar(X)
for i in 1:rows
x = X[i, :]
dy = dY[i, :]
ms = sum(x .* x) / d + eps
rms = sqrt(ms)
inv_rms = 1.0 / rms
dot_dy_x = sum(dy .* x)
dX[i, :] = dy .* inv_rms .- (x .* (dot_dy_x / (d * ms * rms)))
end
return dX
end
function silu(x::Float64)::Float64
return x / (1.0 + exp(-x))
end
function ffn_swiglu(X::Matrix{Float64}, W1::Matrix{Float64},
W2::Matrix{Float64}, W3::Matrix{Float64})::Matrix{Float64}
h1 = X * W1
h3 = X * W3
gated = silu.(h1) .* h3
return gated * W2
end
function ffn_relu(X::Matrix{Float64}, W1::Matrix{Float64},
W2::Matrix{Float64})::Matrix{Float64}
h = X * W1
h = max.(h, 0.0)
return h * W2
end
function scaled_dot_product_attention(Q::Matrix{Float64}, K::Matrix{Float64},
V::Matrix{Float64}; causal::Bool=false)
dk = size(Q, 2)
scores = (Q * transpose(K)) ./ sqrt(dk)
mask = nothing
if causal
n = size(scores, 1)
mask = [j > i for i in 1:n, j in 1:size(scores, 2)]
end
weights = softmax_rows(scores; mask=mask)
return weights * V
end
function multi_head_attention(X::Matrix{Float64},
Wq::Matrix{Float64}, Wk::Matrix{Float64},
Wv::Matrix{Float64}, Wo::Matrix{Float64};
n_heads::Int=1, causal::Bool=false,
kv_source::Union{Nothing, Matrix{Float64}}=nothing)
@assert n_heads > 0 "n_heads must be > 0"
Q = X * Wq
kv_input = kv_source === nothing ? X : kv_source
K = kv_input * Wk
V = kv_input * Wv
d_total = size(Q, 2)
@assert d_total % n_heads == 0 "projected dimension must be divisible by n_heads"
d_head = d_total ÷ n_heads
head_outs = Matrix{Float64}[]
for h in 1:n_heads
cols = ((h - 1) * d_head + 1):(h * d_head)
Qh = Q[:, cols]
Kh = K[:, cols]
Vh = V[:, cols]
push!(head_outs, scaled_dot_product_attention(Qh, Kh, Vh; causal=causal))
end
concat = hcat(head_outs...)
return concat * Wo
end
struct BlockParams
d::Int
n_heads::Int
use_swiglu::Bool
Wq::Matrix{Float64}
Wk::Matrix{Float64}
Wv::Matrix{Float64}
Wo::Matrix{Float64}
W1::Matrix{Float64}
W2::Matrix{Float64}
W3::Matrix{Float64}
Wq_x::Matrix{Float64}
Wk_x::Matrix{Float64}
Wv_x::Matrix{Float64}
Wo_x::Matrix{Float64}
end
function BlockParams(d::Int, n_heads::Int, ffn_expansion::Float64,
rng::AbstractRNG; use_swiglu::Bool=true)
@assert n_heads > 0 "n_heads must be > 0"
@assert d % n_heads == 0 "d must be divisible by n_heads"
h = Int(round(d * ffn_expansion))
Wq = randn_matrix(rng, d, d)
Wk = randn_matrix(rng, d, d)
Wv = randn_matrix(rng, d, d)
Wo = randn_matrix(rng, d, d)
W1 = randn_matrix(rng, d, h)
W2 = randn_matrix(rng, h, d)
W3 = use_swiglu ? randn_matrix(rng, d, h) : zeros(d, h)
Wq_x = randn_matrix(rng, d, d)
Wk_x = randn_matrix(rng, d, d)
Wv_x = randn_matrix(rng, d, d)
Wo_x = randn_matrix(rng, d, d)
return BlockParams(d, n_heads, use_swiglu,
Wq, Wk, Wv, Wo, W1, W2, W3,
Wq_x, Wk_x, Wv_x, Wo_x)
end
function encoder_block(x::Matrix{Float64}, p::BlockParams)::Matrix{Float64}
h = rms_norm(x)
a = multi_head_attention(h, p.Wq, p.Wk, p.Wv, p.Wo; n_heads=p.n_heads)
x = x .+ a
h = rms_norm(x)
f = p.use_swiglu ? ffn_swiglu(h, p.W1, p.W2, p.W3) : ffn_relu(h, p.W1, p.W2)
return x .+ f
end
function decoder_block(x::Matrix{Float64}, enc_out::Matrix{Float64},
p::BlockParams)::Matrix{Float64}
h = rms_norm(x)
a = multi_head_attention(h, p.Wq, p.Wk, p.Wv, p.Wo;
n_heads=p.n_heads, causal=true)
x = x .+ a
h = rms_norm(x)
a = multi_head_attention(h, p.Wq_x, p.Wk_x, p.Wv_x, p.Wo_x;
n_heads=p.n_heads, kv_source=enc_out)
x = x .+ a
h = rms_norm(x)
f = p.use_swiglu ? ffn_swiglu(h, p.W1, p.W2, p.W3) : ffn_relu(h, p.W1, p.W2)
return x .+ f
end
function numerical_grad(f, X::Matrix{Float64}; h::Float64=1e-5)::Matrix{Float64}
out = similar(X)
for i in 1:length(X)
orig = X[i]
X[i] = orig + h
plus = f(X)
X[i] = orig - h
minus = f(X)
X[i] = orig
out[i] = (plus - minus) / (2h)
end
return out
end
function gradient_check_layer_norm()
println("=" ^ 60)
println("LAYER NORM: ANALYTIC vs NUMERICAL GRADIENT")
println("=" ^ 60)
rng = MersenneTwister(0)
X = randn(rng, 4, 6)
rng_v = MersenneTwister(1)
v = randn(rng_v, 4, 6)
loss_fn = Y -> sum(layer_norm(Y) .* v)
analytic = layer_norm_backward(X, v)
numeric = numerical_grad(loss_fn, copy(X))
err = maximum(abs.(analytic .- numeric))
@printf("\nMax abs error (LayerNorm): %.3e\n", err)
end
function gradient_check_rms_norm()
println("\n" * "=" ^ 60)
println("RMS NORM: ANALYTIC vs NUMERICAL GRADIENT")
println("=" ^ 60)
rng = MersenneTwister(2)
X = randn(rng, 4, 6)
rng_v = MersenneTwister(3)
v = randn(rng_v, 4, 6)
loss_fn = Y -> sum(rms_norm(Y) .* v)
analytic = rms_norm_backward(X, v)
numeric = numerical_grad(loss_fn, copy(X))
err = maximum(abs.(analytic .- numeric))
@printf("\nMax abs error (RMSNorm): %.3e\n", err)
end
function compare_norm_outputs()
println("\n" * "=" ^ 60)
println("LAYERNORM vs RMSNORM OUTPUTS")
println("=" ^ 60)
rng = MersenneTwister(7)
X = randn(rng, 3, 6)
Y_ln = layer_norm(X)
Y_rms = rms_norm(X)
println("\nLayerNorm row means (should be ~0):")
for i in 1:3
@printf(" row %d: mean=%+.6f std=%.6f\n",
i, sum(Y_ln[i, :]) / 6, sqrt(sum(Y_ln[i, :] .^ 2) / 6))
end
println("\nRMSNorm row RMS (should be ~1):")
for i in 1:3
@printf(" row %d: mean=%+.6f rms=%.6f\n",
i, sum(Y_rms[i, :]) / 6, sqrt(sum(Y_rms[i, :] .^ 2) / 6))
end
println("\nRMSNorm leaves the row mean intact; LayerNorm centers it.")
end
function demo_full_transformer()
println("\n" * "=" ^ 60)
println("FULL TRANSFORMER FORWARD PASS")
println("=" ^ 60)
rng = MersenneTwister(42)
d = 8
n_heads = 2
ffn_exp = 2.0
src_len = 6
tgt_len = 5
src = randn_matrix(rng, src_len, d; scale=0.5)
tgt = randn_matrix(rng, tgt_len, d; scale=0.5)
enc_params = [BlockParams(d, n_heads, ffn_exp, rng) for _ in 1:2]
dec_params = [BlockParams(d, n_heads, ffn_exp, rng) for _ in 1:2]
enc_out = src
for p in enc_params
enc_out = encoder_block(enc_out, p)
end
dec_out = tgt
for p in dec_params
dec_out = decoder_block(dec_out, enc_out, p)
end
@printf("\nsource shape: (%d, %d)\n", size(src, 1), size(src, 2))
@printf("encoder output shape: (%d, %d)\n", size(enc_out, 1), size(enc_out, 2))
@printf("target shape: (%d, %d)\n", size(tgt, 1), size(tgt, 2))
@printf("decoder output shape: (%d, %d)\n", size(dec_out, 1), size(dec_out, 2))
println("\nfirst 3 rows of encoder output:")
for i in 1:3
println(" " * join([@sprintf("%+.3f", enc_out[i, j]) for j in 1:4], " "))
end
println("\nfirst 3 rows of decoder output:")
for i in 1:3
println(" " * join([@sprintf("%+.3f", dec_out[i, j]) for j in 1:4], " "))
end
println("\nstack: 2-layer encoder + 2-layer decoder, pre-norm, RMSNorm, SwiGLU.")
end
function main()
compare_norm_outputs()
gradient_check_layer_norm()
gradient_check_rms_norm()
demo_full_transformer()
end
if abspath(PROGRAM_FILE) == @__FILE__
main()
end