diff --git a/lib/scikit/cluster.flow b/lib/scikit/cluster.flow index 3aaae37..3731f80 100644 --- a/lib/scikit/cluster.flow +++ b/lib/scikit/cluster.flow @@ -5,6 +5,13 @@ import "lib/scikit/matrix.flow" import "lib/scikit/prng.flow" import "lib/scikit/blas.flow" +extern { + # clang lowers both to a single instruction, so a loop written with them + # holds no branch and the vectorizer can take it. + function fmaxf(a: f32, b: f32) -> f32 + function fminf(a: f32, b: f32) -> f32 +} + export struct KMeans { centroids: ptr >, n_clusters: i32, @@ -2007,17 +2014,34 @@ export function affinity_propagation_fit(X: Matrix, damping: f32, max_iter: i32, } } - for k in 0 to n { - let mut max_other: f32 = best - if k == best_k { max_other = second } - let v: f32 = take * (si3[k] - max_other) + keep * ri2[k] + # Every element of the row takes the row maximum, and the two + # elements that do not are corrected afterwards. Deciding it per + # element put a branch on k against best_k and another on k + # against i inside a loop that is otherwise a fused multiply add + # and a max, which is the difference between a vector loop and a + # scalar one. + let old_best: f32 = ri2[best_k] + # Written as a while loop with a constant step. A `for k in 0 to n` + # emits an increment chosen by a runtime test of the bounds, which + # is not a stride the LLVM vectorizer will accept, and this loop + # was scalar with it. + let mut k: i32 = 0 + while k < n { + let v: f32 = take * (si3[k] - best) + keep * ri2[k] ri2[k] = v - if k == i { - rkk[i] = v - } elif v > 0.0 { - col_pos[k] = col_pos[k] + v - } + col_pos[k] = col_pos[k] + fmaxf(v, 0.0) + k = k + 1 } + + # The largest entry of the row competes against the runner up. + let fixed: f32 = take * (si3[best_k] - second) + keep * old_best + col_pos[best_k] = col_pos[best_k] + fmaxf(fixed, 0.0) - fmaxf(ri2[best_k], 0.0) + ri2[best_k] = fixed + + # The diagonal is the self responsibility, and it is not part of + # the column sum. + rkk[i] = ri2[i] + col_pos[i] = col_pos[i] - fmaxf(ri2[i], 0.0) } # The availability update. The term that depends only on the column is @@ -2029,21 +2053,16 @@ export function affinity_propagation_fit(X: Matrix, damping: f32, max_iter: i32, for i in 0 to n { let ai2: ptr = A + i * n let ri4: ptr = R + i * n - for k in 0 to i { - let mut val: f32 = base[k] - let r_ik: f32 = ri4[k] - if r_ik > 0.0 { val = val - r_ik } - if val > 0.0 { val = 0.0 } - ai2[k] = take * val + keep * ai2[k] - } - ai2[i] = take * base[i] + keep * ai2[i] - for k in i + 1 to n { - let mut val2: f32 = base[k] - let r_ik2: f32 = ri4[k] - if r_ik2 > 0.0 { val2 = val2 - r_ik2 } - if val2 > 0.0 { val2 = 0.0 } - ai2[k] = take * val2 + keep * ai2[k] - } + let old_diag: f32 = ai2[i] + let mut k2: i32 = 0 + while k2 < n { + let val: f32 = fminf(base[k2] - fmaxf(ri4[k2], 0.0), 0.0) + ai2[k2] = take * val + keep * ai2[k2] + k2 = k2 + 1 + } + # The diagonal keeps the whole column sum rather than the clipped + # one, and is written after the loop for the reason above. + ai2[i] = take * base[i] + keep * old_diag } # Convergence is the set of exemplars holding still, which is what diff --git a/lib/scikit/linear.flow b/lib/scikit/linear.flow index 976f312..c8b5193 100644 --- a/lib/scikit/linear.flow +++ b/lib/scikit/linear.flow @@ -7252,8 +7252,12 @@ function _glm_gradient_fit(X: Matrix, y: ptr, weights: ptr, bias_init: for iter in 0 to max_iter { blas_matvec(xdata, n, p, weights, eta, 1.0, 0.0) + # While loops with a constant step. Flow emits a range loop with a + # stride chosen at runtime, which LLVM will not vectorize. Flow issue + # #1004. let mut grad_b: f32 = 0.0 - for i in 0 to n { + let mut i: i32 = 0 + while i < n { let e: f32 = eta[i] + bias let mut mu: f32 = 0.0 if kind == 0 { mu = _poisson_inverse_link(e) } @@ -7262,11 +7266,14 @@ function _glm_gradient_fit(X: Matrix, y: ptr, weights: ptr, bias_init: let r: f32 = y[i] - mu resid[i] = r grad_b = grad_b + r + i = i + 1 } blas_matvec_trans(xdata, n, p, resid, grad_w, 1.0, 0.0) - for j in 0 to p { + let mut j: i32 = 0 + while j < p { let g: f32 = grad_w[j] * inv_n - alpha * weights[j] weights[j] = weights[j] + lr * g + j = j + 1 } bias = bias + lr * grad_b * inv_n } diff --git a/lib/scikit/manifold.flow b/lib/scikit/manifold.flow index d323c36..1a5bd5c 100644 --- a/lib/scikit/manifold.flow +++ b/lib/scikit/manifold.flow @@ -8,6 +8,9 @@ import "lib/scikit/blas.flow" extern { function rand() -> i32 function srand(seed: u32) -> void + # clang lowers this to a single instruction, so a loop written with it + # holds no branch and the vectorizer can take it. + function fminf(a: f32, b: f32) -> f32 } export struct TSNE { @@ -306,125 +309,132 @@ export struct Isomap { fitted: bool } +# Modified Gram-Schmidt on the columns of a row-major n by k block. +# A column that collapses is left at zero, and the remaining columns stay +# orthonormal. +function _orthonormalize_block(V: ptr, n: i32, k: i32) -> void { + for c in 0 to k { + for prev in 0 to c { + let mut dot: f32 = 0.0 + for i in 0 to n { dot = dot + V[i * k + c] * V[i * k + prev] } + for i in 0 to n { V[i * k + c] = V[i * k + c] - dot * V[i * k + prev] } + } + let mut norm: f32 = 0.0 + for i in 0 to n { + let v: f32 = V[i * k + c] + norm = norm + v * v + } + norm = sqrt((norm) as f64) as f32 + if norm > 0.0000000001 { + let inv: f32 = 1.0 / norm + for i in 0 to n { V[i * k + c] = V[i * k + c] * inv } + } + } +} + export function isomap_fit(X: Matrix, n_components: i32, n_neighbors: i32) -> Isomap { let n: i32 = X.rows let p: i32 = X.cols + let kc: i32 = n_components - let D: ptr > = malloc((n as i64) * 8) as ptr > + # Distances, then geodesics, then the double centering, all in contiguous + # blocks. Every inner loop below is a while loop with a constant step, + # because Flow emits a range loop with a stride chosen at runtime and LLVM + # will not vectorize that. Flow issue #1004. + let xdata: ptr = X.data + let geo: ptr = array_new_f32(n * n) for i in 0 to n { - D[i] = array_new_f32(n) - for j in 0 to n { + let ri: ptr = xdata + i * p + let gi: ptr = geo + i * n + for j in i + 1 to n { + let rj: ptr = xdata + j * p let mut s: f32 = 0.0 - for k in 0 to p { - let d: f32 = matrix_at(X, i, k) - matrix_at(X, j, k) + let mut m: i32 = 0 + while m < p { + let d: f32 = ri[m] - rj[m] s = s + d * d + m = m + 1 } - D[i][j] = sqrt((s as f64)) as f32 + let dist: f32 = sqrt((s) as f64) as f32 + gi[j] = dist + geo[j * n + i] = dist } } - let geo: ptr > = malloc((n as i64) * 8) as ptr > - for i in 0 to n { - geo[i] = array_copy_f32(D[i], n) - } - + # Floyd-Warshall. The term that depends only on i and k is lifted out of + # the inner loop, and the comparison is a minimum, so the loop holds no + # branch and vectorizes. for k in 0 to n { + let rowk: ptr = geo + k * n for i in 0 to n { - for j in 0 to n { - if geo[i][k] + geo[k][j] < geo[i][j] { - geo[i][j] = geo[i][k] + geo[k][j] - } + let rowi: ptr = geo + i * n + let dik: f32 = rowi[k] + let mut j: i32 = 0 + while j < n { + rowi[j] = fminf(rowi[j], dik + rowk[j]) + j = j + 1 } } } let mut mean_d: f32 = 0.0 for i in 0 to n { - for j in 0 to n { - geo[i][j] = geo[i][j] * geo[i][j] - mean_d = mean_d + geo[i][j] + let gi2: ptr = geo + i * n + let mut j2: i32 = 0 + while j2 < n { + let v: f32 = gi2[j2] * gi2[j2] + gi2[j2] = v + mean_d = mean_d + v + j2 = j2 + 1 } } mean_d = mean_d / ((n * n) as f32) - let B: ptr > = malloc((n as i64) * 8) as ptr > + let B: ptr = array_new_f32(n * n) for i in 0 to n { - B[i] = array_new_f32(n) - for j in 0 to n { - B[i][j] = -0.5 * (geo[i][j] - mean_d) + let gi3: ptr = geo + i * n + let bi: ptr = B + i * n + let mut j3: i32 = 0 + while j3 < n { + bi[j3] = -0.5 * (gi3[j3] - mean_d) + j3 = j3 + 1 } } + array_free_f32(geo) - let embedding: ptr = array_new_f32(n * n_components) - # Track found eigenvectors for deflation - let found_vecs: ptr > = malloc((n_components as i64) * 8) as ptr > - for c in 0 to n_components { - let v: ptr = array_new_f32(n) - srand(42 + c) - for i in 0 to n { v[i] = ((rand() % 1000) as f32) / 500.0 - 1.0 } - - # Orthogonalize against previously found eigenvectors - for prev in 0 to c { - let mut dot: f32 = 0.0 - for i in 0 to n { dot = dot + v[i] * found_vecs[prev][i] } - for i in 0 to n { v[i] = v[i] - dot * found_vecs[prev][i] } - } - - let mut iter: i32 = 0 - while iter < 200 { - let Bv: ptr = array_new_f32(n) - for i in 0 to n { - let mut s: f32 = 0.0 - for j in 0 to n { s = s + B[i][j] * v[j] } - Bv[i] = s - } - - # Deflate against previous eigenvectors - for prev in 0 to c { - let mut dot: f32 = 0.0 - for i in 0 to n { dot = dot + Bv[i] * found_vecs[prev][i] } - for i in 0 to n { Bv[i] = Bv[i] - dot * found_vecs[prev][i] } - } - - let mut norm: f32 = 0.0 - for i in 0 to n { norm = norm + Bv[i] * Bv[i] } - norm = sqrt((norm as f64)) as f32 - if norm < 0.0000000001 { norm = 0.0000000001 } - - # Convergence check - let mut max_change: f32 = 0.0 - for i in 0 to n { - let new_v: f32 = Bv[i] / norm - let diff: f32 = new_v - v[i] - if diff < 0.0 { diff = -diff } - if diff > max_change { max_change = diff } - v[i] = new_v - } - array_free_f32(Bv) + # Block power iteration. One sgemm per iteration advances every component, + # where the old loop ran a scalar matvec per component per iteration and + # allocated its result buffer inside the iteration. + let V: ptr = array_new_f32(n * kc) + let Z: ptr = array_new_f32(n * kc) + srand(42) + for i in 0 to n * kc { V[i] = ((rand() % 1000) as f32) / 500.0 - 1.0 } + _orthonormalize_block(V, n, kc) - if max_change < 0.000001 { break } - iter = iter + 1 - } + let mut iter: i32 = 0 + while iter < 200 { + blas_matmat(B, V, Z, n, kc, n, 1.0, 0.0) + _orthonormalize_block(Z, n, kc) - for i in 0 to n { - embedding[i * n_components + c] = v[i] + let mut max_change: f32 = 0.0 + for i in 0 to n * kc { + let mut diff: f32 = Z[i] - V[i] + if diff < 0.0 { diff = 0.0 - diff } + if diff > max_change { max_change = diff } + V[i] = Z[i] } - found_vecs[c] = v + if max_change < 0.000001 { break } + iter = iter + 1 } - for c in 0 to n_components { array_free_f32(found_vecs[c]) } - free(found_vecs as ptr) - - for i in 0 to n { array_free_f32(D[i]); array_free_f32(geo[i]); array_free_f32(B[i]) } - free(D as ptr) - free(geo as ptr) - free(B as ptr) + array_free_f32(B) + array_free_f32(Z) return Isomap { - embedding: embedding, + embedding: V, n_samples: n, n_features: p, - n_components: n_components, + n_components: kc, n_neighbors: n_neighbors, fitted: true } @@ -456,29 +466,6 @@ export struct LLE { fitted: bool } -# Modified Gram-Schmidt on the columns of a row-major n by k block. -# A column that collapses is left at zero, and the remaining columns stay -# orthonormal. -function _orthonormalize_block(V: ptr, n: i32, k: i32) -> void { - for c in 0 to k { - for prev in 0 to c { - let mut dot: f32 = 0.0 - for i in 0 to n { dot = dot + V[i * k + c] * V[i * k + prev] } - for i in 0 to n { V[i * k + c] = V[i * k + c] - dot * V[i * k + prev] } - } - let mut norm: f32 = 0.0 - for i in 0 to n { - let v: f32 = V[i * k + c] - norm = norm + v * v - } - norm = sqrt((norm) as f64) as f32 - if norm > 0.0000000001 { - let inv: f32 = 1.0 / norm - for i in 0 to n { V[i * k + c] = V[i * k + c] * inv } - } - } -} - export function lle_fit(X: Matrix, n_components: i32, n_neighbors: i32) -> LLE { let n: i32 = X.rows let p: i32 = X.cols diff --git a/lib/scikit/neighbors.flow b/lib/scikit/neighbors.flow index ceeb059..b0444da 100644 --- a/lib/scikit/neighbors.flow +++ b/lib/scikit/neighbors.flow @@ -212,39 +212,61 @@ export function knn_regressor_predict(model: KNNRegressor, X: Matrix) -> ptr = array_new_f32(X.rows) let neighbors: ptr = malloc((model.k as i64) * 12) as ptr + # Row pointers into both designs, and the worst of the k held in hand. + # The distance loop read every element through matrix_at, which copies the + # Matrix struct per element, and the loop that finds the slot to replace + # ran over all k for every training point when it only has to run when a + # point actually displaces one. The feature loop is a while loop with a + # constant step, because Flow emits a range loop with a stride chosen at + # runtime and LLVM will not vectorize that. Flow issue #1004. + let cols: i32 = X.cols + let xdata: ptr = X.data + let tdata: ptr = model.X_train.data + let n_train: i32 = model.X_train.rows + let k: i32 = model.k + for i in 0 to X.rows { - for j in 0 to model.k { + let xrow: ptr = xdata + i * cols + for j in 0 to k { neighbors[j].distance = 9999999999.0 neighbors[j].label = 0.0 } + let mut worst_idx: i32 = 0 + let mut worst_dist: f32 = 9999999999.0 - for j in 0 to model.X_train.rows { - # Use f64 for distance to match sklearn precision + for j in 0 to n_train { + let trow: ptr = tdata + j * cols + # f64 for the distance, to match scikit-learn's precision. let mut dist: f64 = 0.0 - for f in 0 to X.cols { - let d: f64 = (matrix_at(X, i, f) - matrix_at(model.X_train, j, f)) as f64 + let mut f: i32 = 0 + while f < cols { + let d: f64 = (xrow[f] - trow[f]) as f64 dist = dist + d * d + f = f + 1 } let dist_f32: f32 = dist as f32 - let mut worst_idx: i32 = 0 - for n in 1 to model.k { - if neighbors[n].distance > neighbors[worst_idx].distance { - worst_idx = n - } - } - - if dist_f32 < neighbors[worst_idx].distance { + if dist_f32 < worst_dist { neighbors[worst_idx].distance = dist_f32 neighbors[worst_idx].label = model.y_train[j] + let mut w: i32 = 0 + let mut wd: f32 = neighbors[0].distance + for m in 1 to k { + if neighbors[m].distance > wd { + wd = neighbors[m].distance + w = m + } + } + worst_idx = w + worst_dist = wd } } let mut sum: f32 = 0.0 - for n in 0 to model.k { + for n in 0 to k { sum = sum + neighbors[n].label } - preds[i] = sum / (model.k as f32) + preds[i] = sum / (k as f32) } free(neighbors as ptr) @@ -597,10 +619,15 @@ export function local_outlier_factor_fit(X: Matrix, n_neighbors: i32) -> LocalOu di[i] = 0.0 for j in i + 1 to n { let rj: ptr = xdata + j * p + # A while loop with a constant step, because Flow emits a range + # loop with a stride chosen at runtime and LLVM will not vectorize + # that. Flow issue #1004. let mut acc: f32 = 0.0 - for m in 0 to p { + let mut m: i32 = 0 + while m < p { let d: f32 = ri[m] - rj[m] acc = acc + d * d + m = m + 1 } let dist: f32 = sqrt((acc as f64)) as f32 di[j] = dist diff --git a/lib/scikit/semi_supervised.flow b/lib/scikit/semi_supervised.flow index 376d6f7..05f618c 100644 --- a/lib/scikit/semi_supervised.flow +++ b/lib/scikit/semi_supervised.flow @@ -2,6 +2,7 @@ # Implements LabelPropagation, LabelSpreading, SelfTrainingClassifier. import "lib/scikit/matrix.flow" +import "lib/scikit/blas.flow" extern { function malloc(size: i64) -> ptr @@ -50,10 +51,15 @@ function _build_affinity(X: Matrix, gamma: f32) -> ptr > { row_i[i] = self_affinity for j in i + 1 to n { let base_j: i32 = j * p + # A while loop with a constant step, because Flow emits a range + # loop with a stride chosen at runtime and LLVM will not vectorize + # that. Flow issue #1004. let mut s: f32 = 0.0 - for k in 0 to p { + let mut k: i32 = 0 + while k < p { let d: f32 = X.data[base_i + k] - X.data[base_j + k] s = s + d * d + k = k + 1 } let w: f32 = exp((-gamma * s) as f64) as f32 row_i[j] = w @@ -259,58 +265,69 @@ export function label_spreading_fit(X: Matrix, y: ptr, n_classes: i32, gamm } } - let S: ptr > = malloc((n as i64) * 8) as ptr > + # The propagation is new_Y = S * Y, a matrix product, so it is one sgemm + # against contiguous storage. It used to be a scalar triple loop that + # walked Y down a column of an array of row allocations, which is a cache + # miss per element, and it allocated a fresh n by n_classes set of rows on + # every iteration. + let Sc: ptr = array_new_f32(n * n) for i in 0 to n { - S[i] = array_new_f32(n) - for j in 0 to n { - S[i][j] = alpha * W[i][j] / deg[i] + let srow: ptr = Sc + i * n + let wrow: ptr = W[i] + let scale: f32 = alpha / deg[i] + let mut j: i32 = 0 + while j < n { + srow[j] = scale * wrow[j] + j = j + 1 } } - let Y: ptr > = malloc((n as i64) * 8) as ptr > - let clamp: ptr > = malloc((n as i64) * 8) as ptr > + let nc: i32 = n_classes + let Yc: ptr = array_new_f32(n * nc) + let Ync: ptr = array_new_f32(n * nc) + let clampc: ptr = array_new_f32(n * nc) for i in 0 to n { - Y[i] = array_new_f32(n_classes) - clamp[i] = array_new_f32(n_classes) - if y[i] >= 0 && y[i] < n_classes { - clamp[i][y[i]] = 1.0 - Y[i][y[i]] = 1.0 + if y[i] >= 0 && y[i] < nc { + clampc[i * nc + y[i]] = 1.0 + Yc[i * nc + y[i]] = 1.0 } } + let hold: f32 = 1.0 - alpha + let total: i32 = n * nc let mut n_iter: i32 = 0 for iter in 0 to max_iter { n_iter = iter + 1 let mut max_change: f32 = 0.0 - let new_Y: ptr > = malloc((n as i64) * 8) as ptr > - for i in 0 to n { - new_Y[i] = array_new_f32(n_classes) - } - - for i in 0 to n { - for c in 0 to n_classes { - let mut s: f32 = 0.0 - for j in 0 to n { - s = s + S[i][j] * Y[j][c] - } - new_Y[i][c] = s + (1.0 - alpha) * clamp[i][c] - } - } + blas_matmat(Sc, Yc, Ync, n, nc, n, 1.0, 0.0) - for i in 0 to n { - for c in 0 to n_classes { - let change: f32 = fabs((new_Y[i][c] - Y[i][c]) as f64) as f32 - if change > max_change { max_change = change } - Y[i][c] = new_Y[i][c] - } - array_free_f32(new_Y[i]) + let mut idx: i32 = 0 + while idx < total { + let v: f32 = Ync[idx] + hold * clampc[idx] + let change: f32 = fabs((v - Yc[idx]) as f64) as f32 + if change > max_change { max_change = change } + Yc[idx] = v + idx = idx + 1 } - free(new_Y as ptr) if max_change < tol { break } } + # The model hands back its distributions a row at a time, so the rows are + # built once here. + let Y: ptr > = malloc((n as i64) * 8) as ptr > + for i in 0 to n { + let row: ptr = array_new_f32(nc) + let src: ptr = Yc + i * nc + for c in 0 to nc { row[c] = src[c] } + Y[i] = row + } + array_free_f32(Sc) + array_free_f32(Yc) + array_free_f32(Ync) + array_free_f32(clampc) + let labels: ptr = malloc((n as i64) * 4) as ptr for i in 0 to n { let mut best_c: i32 = 0 @@ -326,12 +343,8 @@ export function label_spreading_fit(X: Matrix, y: ptr, n_classes: i32, gamm for i in 0 to n { array_free_f32(W[i]) - array_free_f32(S[i]) - array_free_f32(clamp[i]) } free(W as ptr) - free(S as ptr) - free(clamp as ptr) array_free_f32(deg) return LabelSpreading { diff --git a/lib/scikit/svm.flow b/lib/scikit/svm.flow index aa9dec8..f98b54a 100644 --- a/lib/scikit/svm.flow +++ b/lib/scikit/svm.flow @@ -2202,61 +2202,68 @@ function _svm_tanh_f32(x: f32) -> f32 { return val } -function _svm_kernel_rows(X1: Matrix, i1: i32, X2: Matrix, i2: i32, n_features: i32, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> f32 { - # The two rows are addressed once. matrix_at copies the Matrix struct by - # value for every element, and this is called once per pair of rows in - # every kernel prediction, so it was two struct copies per feature per - # pair. - let r1: ptr = X1.data + i1 * X1.cols - let r2: ptr = X2.data + i2 * X2.cols - if kernel == KERNEL_LINEAR { - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + r1[k] * r2[k] - } - return s - } - if kernel == KERNEL_POLY { - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + r1[k] * r2[k] - } - let base: f32 = gamma * s + coef0 - let val: f32 = _svm_pow_f32(base, degree) - return val - } +# The kernel between two rows, taken as pointers. +# +# Every loop here is a while loop with a constant step. Flow's `for i in a to b` +# emits an increment chosen by a runtime test of the bounds, which LLVM will not +# accept as a stride, so a `for` loop over an array is scalar code. Flow issue +# #1004. +function _svm_kernel_vec(r1: ptr, r2: ptr, n_features: i32, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> f32 { if kernel == KERNEL_RBF { let mut s: f32 = 0.0 - for k in 0 to n_features { + let mut k: i32 = 0 + while k < n_features { let d: f32 = r1[k] - r2[k] s = s + d * d + k = k + 1 } - let val: f32 = exp((0.0 - gamma * s) as f64) as f32 - return val - } - if kernel == KERNEL_SIGMOID { - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + matrix_at(X1, i1, k) * matrix_at(X2, i2, k) - } - let val: f32 = _svm_tanh_f32(gamma * s + coef0) - return val + return exp((0.0 - gamma * s) as f64) as f32 } - let mut s: f32 = 0.0 - for k in 0 to n_features { - s = s + matrix_at(X1, i1, k) * matrix_at(X2, i2, k) + + let mut dot: f32 = 0.0 + let mut k2: i32 = 0 + while k2 < n_features { + dot = dot + r1[k2] * r2[k2] + k2 = k2 + 1 } - return s + if kernel == KERNEL_LINEAR { return dot } + if kernel == KERNEL_POLY { return _svm_pow_f32(gamma * dot + coef0, degree) } + if kernel == KERNEL_SIGMOID { return _svm_tanh_f32(gamma * dot + coef0) } + return dot } +function _svm_kernel_rows(X1: Matrix, i1: i32, X2: Matrix, i2: i32, n_features: i32, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> f32 { + # The two rows are addressed once. matrix_at copies the Matrix struct by + # value for every element, and this is called once per pair of rows in + # every kernel prediction, so it was two struct copies per feature per + # pair. + let r1: ptr = X1.data + i1 * X1.cols + let r2: ptr = X2.data + i2 * X2.cols + return _svm_kernel_vec(r1, r2, n_features, kernel, gamma, degree, coef0) +} + +# Every kernel here is symmetric, so one triangle is computed and mirrored. +# The entries are written through row pointers, where matrix_set copied the +# Matrix struct by value for each of the n^2 of them. function _svm_kernel_matrix(X: Matrix, kernel: i32, gamma: f32, degree: i32, coef0: f32) -> Matrix { let n: i32 = X.rows + let p: i32 = X.cols let K: Matrix = matrix_new(n, n) - for i in 0 to n { - for j in 0 to n { - let val: f32 = _svm_kernel_rows(X, i, X, j, X.cols, kernel, gamma, degree, coef0) - matrix_set(K, i, j, val) + let kdata: ptr = K.data + let xdata: ptr = X.data + let mut i: i32 = 0 + while i < n { + let ri: ptr = xdata + i * p + let ki: ptr = kdata + i * n + let mut j: i32 = i + while j < n { + let rj: ptr = xdata + j * p + let val: f32 = _svm_kernel_vec(ri, rj, p, kernel, gamma, degree, coef0) + ki[j] = val + kdata[j * n + i] = val + j = j + 1 } + i = i + 1 } return K } @@ -2682,9 +2689,8 @@ export function svr_fit(X: Matrix, y: ptr, C: f32, epsilon: f32, kernel: i3 let adelta: f32 = fabs((delta) as f64) as f32 if adelta > 0.0 { alpha_diff[i] = a_new - for k in 0 to n { - f[k] = f[k] + row[k] * delta - } + # f = f + delta * row is a saxpy, and BLAS has one. + blas_axpy(n, delta, row, f) if adelta > max_change { max_change = adelta } } }