108 lines
2.1 KiB
Odin
108 lines
2.1 KiB
Odin
package main
|
|
|
|
import "core:hash"
|
|
import "core:math"
|
|
import "core:math/rand"
|
|
|
|
// Gaussian elimination solver
|
|
solve_linear_system :: proc(A: [][]f64, b: []f64) -> []f64 {
|
|
n := len(A)
|
|
|
|
aug := make([][]f64, n)
|
|
for i in 0 ..< n {
|
|
aug[i] = make([]f64, n + 1)
|
|
copy(aug[i][:n], A[i])
|
|
aug[i][n] = b[i]
|
|
}
|
|
defer {
|
|
for row in aug {delete(row)}
|
|
delete(aug)
|
|
}
|
|
|
|
// Forward elimination with partial pivoting
|
|
for col in 0 ..< n {
|
|
max_row := col
|
|
for row in col + 1 ..< n {
|
|
if math.abs(aug[row][col]) > math.abs(aug[max_row][col]) {
|
|
max_row = row
|
|
}
|
|
}
|
|
|
|
if max_row != col {
|
|
aug[col], aug[max_row] = aug[max_row], aug[col]
|
|
}
|
|
|
|
if math.abs(aug[col][col]) < 1e-10 {
|
|
x := make([]f64, n)
|
|
return x
|
|
}
|
|
|
|
for row in col + 1 ..< n {
|
|
factor := aug[row][col] / aug[col][col]
|
|
for j in col ..< n + 1 {
|
|
aug[row][j] -= factor * aug[col][j]
|
|
}
|
|
}
|
|
}
|
|
|
|
// Back substitution
|
|
x := make([]f64, n)
|
|
for i in 0 ..< n {
|
|
row := n - 1 - i
|
|
x[row] = aug[row][n]
|
|
for j in row + 1 ..< n {
|
|
x[row] -= aug[row][j] * x[j]
|
|
}
|
|
x[row] /= aug[row][row]
|
|
}
|
|
|
|
return x
|
|
}
|
|
|
|
train_linear_regression_cached :: proc(selected_features: []int, y: []f64) -> []f64 {
|
|
m := len(selected_features)
|
|
|
|
// Extract submatrix from cache
|
|
XtX := make([][]f64, m)
|
|
defer {
|
|
for row in XtX {delete(row)}
|
|
delete(XtX)
|
|
}
|
|
|
|
for i in 0 ..< m {
|
|
XtX[i] = make([]f64, m)
|
|
for j in 0 ..< m {
|
|
XtX[i][j] = feature_state.XtX_cache[selected_features[i]][selected_features[j]]
|
|
}
|
|
}
|
|
|
|
Xty := make([]f64, m)
|
|
defer delete(Xty)
|
|
for i in 0 ..< m {
|
|
Xty[i] = feature_state.Xty_cache[selected_features[i]]
|
|
}
|
|
|
|
return solve_linear_system(XtX, Xty)
|
|
}
|
|
|
|
predict :: proc(X: [][]f64, beta: []f64) -> []f64 {
|
|
predictions := make([]f64, len(X))
|
|
for i in 0 ..< len(X) {
|
|
sum := 0.0
|
|
for j in 0 ..< len(beta) {
|
|
sum += X[i][j] * beta[j]
|
|
}
|
|
predictions[i] = sum
|
|
}
|
|
return predictions
|
|
}
|
|
|
|
rmse :: proc(predictions: []f64, actual: []f64) -> f64 {
|
|
sum := 0.0
|
|
for i in 0 ..< len(predictions) {
|
|
diff := predictions[i] - actual[i]
|
|
sum += diff * diff
|
|
}
|
|
return math.sqrt(sum / f64(len(predictions)))
|
|
}
|