From f4862a65d4a0ebd2933c06341001624ea3f2ab71 Mon Sep 17 00:00:00 2001 From: Ulises Jeremias Cornejo Fandos Date: Sat, 30 Apr 2022 00:47:38 -0300 Subject: [PATCH] vsl.la: Reimplement vsl.la and vsl.ml using generics (#82) * Initial iteration for la with generics * Updated la * Updated CI * Updated functions * Updated functions --- .github/workflows/ci.yml | 2 + blas/README.md | 2 +- la/blas.v | 173 +++++++++++++++++++++++++++------------ la/matrix.v | 18 ++-- la/matrix_ops.v | 6 +- ml/data.v | 58 +++++++------ ml/kmeans.v | 12 +-- ml/knn.v | 4 +- ml/linreg.v | 46 +++++------ ml/linreg_test.v | 2 +- ml/paramsreg.v | 78 ++++++++++-------- ml/paramsreg_test.v | 2 +- ml/workspace.v | 64 +++++++-------- ml/workspace_test.v | 2 +- 14 files changed, 275 insertions(+), 194 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 74c69ced6..914a0b990 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,6 +12,8 @@ on: - cron: "31 1,12 * * *" push: + branches: + - master pull_request: branches: diff --git a/blas/README.md b/blas/README.md index a38c83f9d..0df4f086c 100644 --- a/blas/README.md +++ b/blas/README.md @@ -22,7 +22,7 @@ import vsl.la // // c := alpha⋅a⋅b ⇒ cij := alpha * aik * bkj // -pub fn matrix_matrix_mul(mut c la.Matrix, alpha f64, a la.Matrix, b la.Matrix) { +pub fn matrix_matrix_mul(mut c la.Matrix, alpha f64, a &la.Matrix, b &la.Matrix) { if c.m < 6 && c.n < 6 && a.n < 30 { for i := 0; i < c.m; i++ { for j := 0; j < c.n; j++ { diff --git a/la/blas.v b/la/blas.v index 34f09228d..0a4c8f16e 100644 --- a/la/blas.v +++ b/la/blas.v @@ -3,6 +3,15 @@ module la import vsl.blas import math +// @todo: @ulises-jeremias to remove this once https://github.com/vlang/v/issues/14047 is finished +fn arr_to_f64arr(arr []T) []f64 { + mut ret := []f64{cap: arr.len} + for v in arr { + ret << f64(v) + } + return ret +} + /* * vector_rms_error returns the scaled root-mean-square of the difference between two vectors * with components normalised by a scaling factor @@ -16,57 +25,81 @@ import math * * scale[i] = a + m*|s[i]| */ -pub fn vector_rms_error(u []f64, v []f64, a f64, m f64, s []f64) f64 { - mut rms := 0.0 - for i := 0; i < u.len; i++ { +pub fn vector_rms_error(u []T, v []T, a T, m T, s []T) T { + mut rms := T(0) + for i in 0 .. u.len { scale := a + m * math.abs(s[i]) err := math.abs(u[i] - v[i]) rms += err * err / (scale * scale) } - return math.sqrt(rms / f64(u.len)) + return T(math.sqrt(f64(rms) / f64(u.len))) } // vector_dot returns the dot product between two vectors: // s := u・v -pub fn vector_dot(u []f64, v []f64) f64 { - mut res := 0.0 - cutoff := 150 - if u.len <= cutoff { - for i := 0; i < u.len; i++ { +pub fn vector_dot(u []T, v []T) T { + $if T is f64 { + mut res := T(0) + cutoff := 150 + if u.len <= cutoff { + for i in 0 .. u.len { + res += u[i] * v[i] + } + return res + } + return blas.ddot(u.len, arr_to_f64arr(u), 1, arr_to_f64arr(v), 1) + } $else { + mut res := T(0) + for i in 0 .. u.len { res += u[i] * v[i] } return res } - return blas.ddot(u.len, u, 1, v, 1) } // vector_add adds the scaled components of two vectors // res := alpha⋅u + beta⋅v ⇒ result[i] := alpha⋅u[i] + beta⋅v[i] -pub fn vector_add(alpha f64, u []f64, beta f64, v []f64) []f64 { - mut res := []f64{len: v.len} - n := u.len - cutoff := 150 - if beta == 1 && n > cutoff { - res = v.clone() - blas.daxpy(n, alpha, u, 1, mut res, 1) +pub fn vector_add(alpha T, u []T, beta T, v []T) []T { + $if T is f64 { + mut res := []f64{len: v.len} + n := u.len + cutoff := 150 + if beta == 1 && n > cutoff { + res = v.clone() + blas.daxpy(n, alpha, arr_to_f64arr(u), 1, mut res, 1) + return res + } + m := n % 4 + for i in 0 .. m { + res[i] = alpha * u[i] + beta * v[i] + } + for i := m; i < n; i += 4 { + res[i + 0] = alpha * u[i + 0] + beta * v[i + 0] + res[i + 1] = alpha * u[i + 1] + beta * v[i + 1] + res[i + 2] = alpha * u[i + 2] + beta * v[i + 2] + res[i + 3] = alpha * u[i + 3] + beta * v[i + 3] + } + return res + } $else { + mut res := []T{len: v.len} + n := u.len + m := n % 4 + for i in 0 .. m { + res[i] = alpha * u[i] + beta * v[i] + } + for i := m; i < n; i += 4 { + res[i + 0] = alpha * u[i + 0] + beta * v[i + 0] + res[i + 1] = alpha * u[i + 1] + beta * v[i + 1] + res[i + 2] = alpha * u[i + 2] + beta * v[i + 2] + res[i + 3] = alpha * u[i + 3] + beta * v[i + 3] + } return res } - m := n % 4 - for i in 0 .. m { - res[i] = alpha * u[i] + beta * v[i] - } - for i := m; i < n; i += 4 { - res[i + 0] = alpha * u[i + 0] + beta * v[i + 0] - res[i + 1] = alpha * u[i + 1] + beta * v[i + 1] - res[i + 2] = alpha * u[i + 2] + beta * v[i + 2] - res[i + 3] = alpha * u[i + 3] + beta * v[i + 3] - } - return res } // vector_max_diff returns the maximum absolute difference between two vectors // maxdiff = max(|u - v|) -pub fn vector_max_diff(u []f64, v []f64) f64 { +pub fn vector_max_diff(u []T, v []T) T { mut maxdiff := math.abs(u[0] - v[0]) for i := 1; i < u.len; i++ { diff := math.abs(u[i] - v[i]) @@ -79,9 +112,9 @@ pub fn vector_max_diff(u []f64, v []f64) f64 { // vector_scale_abs creates a "scale" vector using the absolute value of another vector // scale := a + m ⋅ |x| ⇒ scale[i] := a + m ⋅ |x[i]| -pub fn vector_scale_abs(a f64, m f64, x []f64) []f64 { - mut scale := []f64{len: x.len} - for i := 0; i < x.len; i++ { +pub fn vector_scale_abs(a T, m T, x []T) []T { + mut scale := []T{len: x.len} + for i in 0 .. x.len { scale[i] = a + m * math.abs(x[i]) } return scale @@ -91,10 +124,24 @@ pub fn vector_scale_abs(a f64, m f64, x []f64) []f64 { // // v = alpha⋅a⋅u ⇒ vi = alpha * aij * uj // -pub fn matrix_vector_mul(alpha f64, a &Matrix, u []f64) []f64 { - mut v := []f64{len: a.m} - if a.m < 9 && a.n < 9 { - for i := 0; i < a.m; i++ { +pub fn matrix_vector_mul(alpha T, a &Matrix, u []T) []T { + $if T is f64 { + mut v := []f64{len: a.m} + if a.m < 9 && a.n < 9 { + for i in 0 .. a.m { + v[i] = 0.0 + for j := 0; j < a.n; j++ { + v[i] += alpha * a.get(i, j) * u[j] + } + } + return v + } + blas.dgemv(false, a.m, a.n, alpha, arr_to_f64arr(a.data), a.m, arr_to_f64arr(u), + 1, 0.0, mut v, v.len) + return v + } $else { + mut v := []T{len: a.m} + for i in 0 .. a.m { v[i] = 0.0 for j := 0; j < a.n; j++ { v[i] += alpha * a.get(i, j) * u[j] @@ -102,18 +149,30 @@ pub fn matrix_vector_mul(alpha f64, a &Matrix, u []f64) []f64 { } return v } - blas.dgemv(false, a.m, a.n, alpha, a.data, a.m, u, 1, 0.0, mut v, v.len) - return v } // matrix_tr_vector_mul returns the transpose(matrix)-vector multiplication // // v = alpha⋅aᵀ⋅u ⇒ vi = alpha * aji * uj = alpha * uj * aji // -pub fn matrix_tr_vector_mul(alpha f64, a &Matrix, u []f64) []f64 { - mut v := []f64{len: a.n} - if a.m < 9 && a.n < 9 { - for i := 0; i < a.n; i++ { +pub fn matrix_tr_vector_mul(alpha T, a &Matrix, u []T) []T { + $if T is f64 { + mut v := []f64{len: a.n} + if a.m < 9 && a.n < 9 { + for i in 0 .. a.n { + v[i] = 0.0 + for j := 0; j < a.m; j++ { + v[i] += alpha * a.get(j, i) * u[j] + } + } + return v + } + blas.dgemv(true, a.m, a.n, alpha, arr_to_f64arr(a.data), a.n, arr_to_f64arr(u), + 1, 0.0, mut v, v.len) + return v + } $else { + mut v := []T{len: a.n} + for i in 0 .. a.n { v[i] = 0.0 for j := 0; j < a.m; j++ { v[i] += alpha * a.get(j, i) * u[j] @@ -121,8 +180,6 @@ pub fn matrix_tr_vector_mul(alpha f64, a &Matrix, u []f64) []f64 { } return v } - blas.dgemv(true, a.m, a.n, alpha, a.data, a.n, u, 1, 0.0, mut v, v.len) - return v } // vector_vector_tr_mul returns the matrix = vector-transpose(vector) multiplication @@ -130,9 +187,24 @@ pub fn matrix_tr_vector_mul(alpha f64, a &Matrix, u []f64) []f64 { // // a = alpha⋅u⋅vᵀ ⇒ aij = alpha * ui * vj // -pub fn vector_vector_tr_mul(alpha f64, u []f64, v []f64) &Matrix { - mut m := new_matrix(u.len, v.len) - if m.m < 9 && m.n < 9 { +pub fn vector_vector_tr_mul(alpha T, u []T, v []T) &Matrix { + $if T is f64 { + mut m := new_matrix(u.len, v.len) + if m.m < 9 && m.n < 9 { + for i in 0 .. m.m { + for j in 0 .. m.n { + m.set(i, j, alpha * u[i] * v[j]) + } + } + return m + } + mut a := []f64{len: u.len * v.len} + blas.dger(m.m, m.n, alpha, arr_to_f64arr(u), 1, arr_to_f64arr(v), 1, mut + a, int(math.max(m.m, m.n))) + return matrix_raw(u.len, v.len, a) + } $else { + mut m := new_matrix(u.len, v.len) + for i in 0 .. m.m { for j in 0 .. m.n { m.set(i, j, alpha * u[i] * v[j]) @@ -140,9 +212,6 @@ pub fn vector_vector_tr_mul(alpha f64, u []f64, v []f64) &Matrix { } return m } - mut a := []f64{len: u.len * v.len} - blas.dger(m.m, m.n, alpha, u, 1, v, 1, mut a, int(math.max(m.m, m.n))) - return matrix_raw(u.len, v.len, a) } // matrix_vector_mul_add returns the matrix-vector multiplication with addition @@ -161,7 +230,7 @@ pub fn matrix_vector_mul_add(alpha f64, a &Matrix, u []f64) []f64 { // pub fn matrix_matrix_mul(mut c Matrix, alpha f64, a &Matrix, b &Matrix) { if c.m < 6 && c.n < 6 && a.n < 30 { - for i := 0; i < c.m; i++ { + for i in 0 .. c.m { for j := 0; j < c.n; j++ { c.set(i, j, 0.0) for k := 0; k < a.n; k++ { @@ -181,7 +250,7 @@ pub fn matrix_matrix_mul(mut c Matrix, alpha f64, a &Matrix, b &Matrix // pub fn matrix_tr_matrix_mul(mut c Matrix, alpha f64, a &Matrix, b &Matrix) { if c.m < 6 && c.n < 6 && a.m < 30 { - for i := 0; i < c.m; i++ { + for i in 0 .. c.m { for j := 0; j < c.n; j++ { c.set(i, j, 0.0) for k := 0; k < a.m; k++ { diff --git a/la/matrix.v b/la/matrix.v index 807143486..727043b2d 100644 --- a/la/matrix.v +++ b/la/matrix.v @@ -1,8 +1,8 @@ module la +import math import strconv import vsl.errors -import math [heap] pub struct Matrix { @@ -48,7 +48,7 @@ pub fn matrix_raw(m int, n int, rawdata []T) &Matrix { pub fn (mut o Matrix) set_from_deep2(a [][]T) { mut k := 0 for i in 0 .. o.m { - for j := 0; j < o.n; j++ { + for j in 0 .. o.n { o.data[k] = a[i][j] k++ } @@ -58,7 +58,7 @@ pub fn (mut o Matrix) set_from_deep2(a [][]T) { // set_diag sets diagonal matrix with diagonal components equal to val pub fn (mut o Matrix) set_diag(val T) { for i in 0 .. o.m { - for j := 0; j < o.n; j++ { + for j in 0 .. o.n { if i == j { o.data[i * o.n + j] = val } else { @@ -82,7 +82,7 @@ pub fn (o &Matrix) get(i int, j int) T { pub fn (o &Matrix) get_deep2() [][]T { mut m := [][]T{len: o.m, init: []T{len: o.n}} for i in 0 .. o.m { - for j := 0; j < o.n; j++ { + for j in 0 .. o.n { m[i][j] = o.data[i * o.n + j] } } @@ -137,7 +137,7 @@ pub fn (mut o Matrix) fill(val T) { */ pub fn (mut o Matrix) clear_rc(rows []int, cols []int, diag T) { for r in rows { - for j := 0; j < o.n; j++ { + for j in 0 .. o.n { if r == j { o.set(r, j, diag) } else { @@ -280,7 +280,7 @@ pub fn (o &Matrix) norm_inf() T { mut sumrow := 0.0 for i := 1; i < o.m; i++ { sumrow = 0.0 - for j := 0; j < o.n; j++ { // sum the other rows + for j in 0 .. o.n { // sum the other rows sumrow += math.abs(o.data[i * o.n + j]) if sumrow > nrm { nrm = sumrow @@ -314,7 +314,7 @@ pub fn (o &Matrix) print(nfmt_ string) string { if i > 0 { l += '\n' } - for j := 0; j < o.n; j++ { + for j in 0 .. o.n { l += safe_print(nfmt, o.get(i, j)) } } @@ -330,7 +330,7 @@ pub fn (o &Matrix) print_v(nfmt_ string) string { mut l := '[][]$T.name{\n' for i in 0 .. o.m { l += ' {' - for j := 0; j < o.n; j++ { + for j in 0 .. o.n { if j > 0 { l += ',' } @@ -351,7 +351,7 @@ pub fn (o &Matrix) print_py(nfmt_ string) string { mut l := 'np.matrix([\n' for i in 0 .. o.m { l += ' [' - for j := 0; j < o.n; j++ { + for j in 0 .. o.n { if j > 0 { l += ',' } diff --git a/la/matrix_ops.v b/la/matrix_ops.v index a0cb5ccda..7a6cc8bdc 100644 --- a/la/matrix_ops.v +++ b/la/matrix_ops.v @@ -124,8 +124,8 @@ pub fn matrix_inv(mut ai Matrix, mut a Matrix, calc_det bool) T { } // singular value decomposition mut s := []T{len: int(math.min(a.m, a.n))} - mut u := new_matrix(a.m, a.m) - mut vt := new_matrix(a.n, a.n) + mut u := new_matrix(a.m, a.m) + mut vt := new_matrix(a.n, a.n) matrix_svd(mut s, mut u, mut vt, mut a, true) // pseudo inverse tol_s := 1e-8 // TODO: improve this tolerance with a better estimate @@ -149,7 +149,7 @@ pub fn matrix_inv(mut ai Matrix, mut a Matrix, calc_det bool) T { // "F" or "" (default) => Frobenius pub fn matrix_cond_num(mut a Matrix, normtype string) T { mut res := 0.0 - mut ai := new_matrix(a.m, a.n) + mut ai := new_matrix(a.m, a.n) matrix_inv(mut ai, mut a, false) if normtype == 'I' { res = a.norm_inf() * ai.norm_inf() diff --git a/ml/data.v b/ml/data.v index b4ca82f6a..d79f5c9d0 100644 --- a/ml/data.v +++ b/ml/data.v @@ -5,7 +5,7 @@ import vsl.la import vsl.errors /* -* Data holds data in matrix format; e.g. for regression computations +* Data holds data in matrix format; e.g. for regression computations * * Example: * _ _ _ _ @@ -19,13 +19,13 @@ import vsl.errors * */ [heap] -pub struct Data { - util.Observable +pub struct Data { pub mut: + observable util.Observable nb_samples int // number of data points (samples). number of rows in x and y nb_features int // number of features. number of columns in x - x &la.Matrix // [nb_samples][nb_features] x values - y []f64 // [nb_samples] y values [optional] + x &la.Matrix // [nb_samples][nb_features] x values + y []T // [nb_samples] y values [optional] } // new_data returns a new object to hold ML data @@ -37,17 +37,13 @@ pub mut: // x and y must be set using set() method // Output: // new object -pub fn new_data(nb_samples int, nb_features int, use_y bool, allocate bool) ?&Data { - x := if allocate { - la.new_matrix(nb_samples, nb_features) - } else { - la.new_matrix(0, 0) - } - mut y := []f64{} +pub fn new_data(nb_samples int, nb_features int, use_y bool, allocate bool) ?&Data { + x := if allocate { la.new_matrix(nb_samples, nb_features) } else { la.new_matrix(0, 0) } + mut y := []T{} if allocate && use_y { - y = []f64{len: nb_samples} + y = []T{len: nb_samples} } - return &Data{ + return &Data{ x: x y: y nb_samples: nb_samples @@ -59,10 +55,10 @@ pub fn new_data(nb_samples int, nb_features int, use_y bool, allocate bool) ?&Da // Input: // x -- x values // y -- y values [optional] -pub fn (mut o Data) set(x &la.Matrix, y []f64) { +pub fn (mut o Data) set(x &la.Matrix, y []T) { o.x = x o.y = y - o.notify_update() + o.observable.notify_update() } // data_from_raw_x returns a new object with data set from raw x values @@ -70,7 +66,7 @@ pub fn (mut o Data) set(x &la.Matrix, y []f64) { // xraw -- [nb_samples][nb_features] table with x values (NO y values) // Output: // new object -pub fn data_from_raw_x(xraw [][]f64) ?&Data { +pub fn data_from_raw_x(xraw [][]T) ?&Data { // check nb_samples := xraw.len if nb_samples < 1 { @@ -78,7 +74,7 @@ pub fn data_from_raw_x(xraw [][]f64) ?&Data { } // allocate new object nb_features := xraw[0].len - mut o := new_data(nb_samples, nb_features, true, true) ? + mut o := new_data(nb_samples, nb_features, true, true) ? // copy data from raw table to x matrix for i := 0; i < nb_samples; i++ { for j := 0; j < nb_features; j++ { @@ -88,11 +84,11 @@ pub fn data_from_raw_x(xraw [][]f64) ?&Data { return o } -// data_from_raw_xy_sep accepts two parameters: xraw [][]f64 and -// yraw []f64. It acts similarly to data_from_raw_xy, but instead +// data_from_raw_xy_sep accepts two parameters: xraw [][]T and +// yraw []T. It acts similarly to data_from_raw_xy, but instead // of using the last column of xraw as the y data, it uses yraw // instead. -pub fn data_from_raw_xy_sep(xraw [][]f64, yraw []f64) ?&Data { +pub fn data_from_raw_xy_sep(xraw [][]T, yraw []T) ?&Data { // check nb_samples := xraw.len if nb_samples < 1 { @@ -100,7 +96,7 @@ pub fn data_from_raw_xy_sep(xraw [][]f64, yraw []f64) ?&Data { } // allocate new object nb_features := xraw[0].len - mut o := new_data(nb_samples, nb_features, false, true) ? + mut o := new_data(nb_samples, nb_features, false, true) ? // copy data from raw table to x matrix for i := 0; i < nb_samples; i++ { for j := 0; j < nb_features; j++ { @@ -119,7 +115,7 @@ pub fn data_from_raw_xy_sep(xraw [][]f64, yraw []f64) ?&Data { // where the last column contains y-values // Output: // new object -pub fn data_from_raw_xy(xyraw [][]f64) ?&Data { +pub fn data_from_raw_xy(xyraw [][]T) ?&Data { // check nb_samples := xyraw.len if nb_samples < 1 { @@ -127,7 +123,7 @@ pub fn data_from_raw_xy(xyraw [][]f64) ?&Data { } // allocate new object nb_features := xyraw[0].len - 1 // -1 because of y column - mut o := new_data(nb_samples, nb_features, true, true) ? + mut o := new_data(nb_samples, nb_features, true, true) ? // copy data from raw table to x and y arrays for i := 0; i < nb_samples; i++ { for j := 0; j < nb_features; j++ { @@ -139,12 +135,22 @@ pub fn data_from_raw_xy(xyraw [][]f64) ?&Data { } // clone returns a deep copy of this object -pub fn (o &Data) clone() ?&Data { +pub fn (o &Data) clone() ?&Data { use_y := o.y.len > 0 - mut p := new_data(o.nb_samples, o.nb_features, use_y, true) ? + mut p := new_data(o.nb_samples, o.nb_features, use_y, true) ? o.x.copy_into(mut p.x, 1) if use_y { p.y = o.y.clone() } return p } + +// add_observer adds an object to the list of interested observers +pub fn (mut o Data) add_observer(obs util.Observer) { + o.observable.add_observer(obs) +} + +// notify_update notifies observers of updates +pub fn (mut o Data) notify_update() { + o.observable.notify_update() +} diff --git a/ml/kmeans.v b/ml/kmeans.v index 83bcbad67..2862c408a 100644 --- a/ml/kmeans.v +++ b/ml/kmeans.v @@ -8,11 +8,11 @@ import vsl.la [heap] pub struct Kmeans { mut: - name string // name of this "observer" - data &Data // x data - stat &Stat // statistics about x (data) - nb_classes int // expected number of classes - bins &gm.Bins // "bins" to speed up searching for data points given their coordinates (2D or 3D only at the moment) + name string // name of this "observer" + data &Data // x data + stat &Stat // statistics about x (data) + nb_classes int // expected number of classes + bins &gm.Bins // "bins" to speed up searching for data points given their coordinates (2D or 3D only at the moment) pub mut: classes []int // [nb_samples] indices of classes of each sample centroids [][]f64 // [nb_classes][nb_features] coordinates of centroids @@ -20,7 +20,7 @@ pub mut: } // new_kmeans returns a new K-means model -pub fn new_kmeans(mut data Data, nb_classes int, name string) &Kmeans { +pub fn new_kmeans(mut data Data, nb_classes int, name string) &Kmeans { // classes classes := []int{len: data.nb_samples} centroids := [][]f64{len: nb_classes} diff --git a/ml/knn.v b/ml/knn.v index 9d82e0899..20a0e3792 100644 --- a/ml/knn.v +++ b/ml/knn.v @@ -8,7 +8,7 @@ import vsl.errors pub struct KNN { mut: name string // name of this "observer" - data &Data + data &Data weights map[f64]f64 // weights[class] = weight pub mut: neighbors []Neighbor @@ -29,7 +29,7 @@ mut: // ```mut knn := new_knn(mut data_from_raw_xy_sep([[0.0, 0.0], [10.0, 10.0]], [0.0, 1.0]))``` // If you predict with `knn.predict(1, [9.0, 9.0])`, it should return 1.0 as it is the closest // to [10.0, 10.0] (which is class 1.0). -pub fn new_knn(mut data Data, name string) ?&KNN { +pub fn new_knn(mut data Data, name string) ?&KNN { if data.x.data.len == 0 { return errors.error('with name $name expects `data.x` to have at least one element.', .einval) diff --git a/ml/linreg.v b/ml/linreg.v index c9af0b103..8f64987aa 100644 --- a/ml/linreg.v +++ b/ml/linreg.v @@ -5,12 +5,12 @@ import vsl.la // LinReg implements a linear regression model [heap] pub struct LinReg { - ParamsReg mut: + params &ParamsReg // main - name string // name of this "observer" - data &Data // x-y data - stat &Stat // statistics + name string // name of this "observer" + data &Data // x-y data + stat &Stat // statistics // workspace e []f64 // vector e = b⋅o + x⋅theta - y [nb_samples] } @@ -19,16 +19,17 @@ mut: // Input: // data -- x,y data // name -- unique name of this (observer) object -pub fn new_lin_reg(mut data Data, name string) &LinReg { +pub fn new_lin_reg(mut data Data, name string) &LinReg { mut stat := stat_from_data(mut data, 'stat_' + name) stat.update() + params := new_params_reg(data.nb_features) mut reg := &LinReg{ name: name data: data stat: stat e: []f64{len: data.nb_samples} + params: params } - reg.init(data.nb_features) return reg } @@ -43,8 +44,8 @@ pub fn (o &LinReg) name() string { // Output: // y -- model prediction y(x) pub fn (o &LinReg) predict(x []f64) f64 { - theta := o.access_thetas() - b := o.get_bias() + theta := o.params.access_thetas() + b := o.params.get_bias() return b + la.vector_dot(x, theta) // b + xᵀtheta } @@ -58,8 +59,8 @@ pub fn (o &LinReg) predict(x []f64) f64 { pub fn (mut o LinReg) cost() f64 { // auxiliary m_1 := 1.0 / f64(o.data.nb_samples) - lambda := o.get_lambda() - theta := o.access_thetas() + lambda := o.params.get_lambda() + theta := o.params.access_thetas() // cost o.calce() // e := b⋅o + x⋅theta - y mut c := (0.5 * m_1) * la.vector_dot(o.e, o.e) // C := (0.5/m) eᵀe @@ -76,8 +77,8 @@ pub fn (mut o LinReg) cost() f64 { pub fn (mut o LinReg) gradients() ([]f64, f64) { // auxiliary m_1 := 1.0 / f64(o.data.nb_samples) - lambda := o.get_lambda() - theta := o.access_thetas() + lambda := o.params.get_lambda() + theta := o.params.access_thetas() x := o.data.x // dcdtheta o.calce() // e := b⋅o + x⋅theta - y @@ -96,7 +97,7 @@ pub fn (mut o LinReg) gradients() ([]f64, f64) { // params -- theta and b pub fn (mut o LinReg) train() { // auxiliary - lambda := o.get_lambda() + lambda := o.params.get_lambda() x, y := o.data.x, o.data.y s, t := o.stat.sum_vars() // r vector @@ -117,17 +118,17 @@ pub fn (mut o LinReg) train() { } } // solve system - mut theta := o.access_thetas() + mut theta := o.params.access_thetas() la.den_solve(mut theta, k, r, false) b_ := (t - la.vector_dot(s, theta)) * m_1 - o.set_bias(b_) + o.params.set_bias(b_) } // calce calculates e vector (save into o.e) // Output: e = b⋅o + x⋅theta - y pub fn (mut o LinReg) calce() { - theta := o.access_thetas() - b := o.get_bias() + theta := o.params.access_thetas() + b := o.params.get_bias() x, y := o.data.x, o.data.y o.e = [b] { @@ -142,15 +143,8 @@ pub fn (mut o LinReg) calce() { pub fn (o &LinReg) str() string { mut res := []string{} res << 'vsl.ml.LinReg{' - res << ' name: $o.name' - res << ' theta: $o.theta' - res << ' bias: $o.bias' - res << ' lambda: $o.lambda' - res << ' degree: $o.degree' - res << ' bkp_theta: $o.bkp_theta' - res << ' bkp_bias: $o.bkp_bias' - res << ' bkp_lambda: $o.bkp_lambda' - res << ' bkp_degree: $o.bkp_degree' + res << ' name: $o.name' + res << ' params: $o.params' res << ' stat: $o.stat' res << ' e: $o.e' res << '}' diff --git a/ml/linreg_test.v b/ml/linreg_test.v index 21e8fa1d9..f13c19562 100644 --- a/ml/linreg_test.v +++ b/ml/linreg_test.v @@ -66,7 +66,7 @@ fn test_lin_reg_pred() { mut reg := new_lin_reg(mut data, 'linear regression') // set regularization parameter - reg.set_lambda(1e12) // very high bias => constant line + reg.params.set_lambda(1e12) // very high bias => constant line reg.train() diff --git a/ml/paramsreg.v b/ml/paramsreg.v index d7ecbf571..5e9826b70 100644 --- a/ml/paramsreg.v +++ b/ml/paramsreg.v @@ -4,39 +4,39 @@ import vsl.la import vsl.util [heap] -pub struct ParamsReg { - util.Observable +pub struct ParamsReg { pub mut: + observable util.Observable // main - theta []f64 // theta parameter [nb_features] - bias f64 // bias parameter - lambda f64 // regularization parameter - degree int // degree of polynomial + theta []T // theta parameter [nb_features] + bias T // bias parameter + lambda T // regularization parameter + degree int // degree of polynomial // backup - bkp_theta []f64 // copy of theta - bkp_bias f64 // copy of b - bkp_lambda f64 // copy of lambda - bkp_degree int // copy of degree + bkp_theta []T // copy of theta + bkp_bias T // copy of b + bkp_lambda T // copy of lambda + bkp_degree int // copy of degree } // new_params_reg returns a new object to hold regression parameters -pub fn new_params_reg(nb_features int) &ParamsReg { - theta := []f64{len: nb_features} - bkp_theta := []f64{len: nb_features} - return &ParamsReg{ +pub fn new_params_reg(nb_features int) &ParamsReg { + theta := []T{len: nb_features} + bkp_theta := []T{len: nb_features} + return &ParamsReg{ theta: theta bkp_theta: bkp_theta } } // init initializes ParamsReg with nb_features (number of features) -pub fn (mut o ParamsReg) init(nb_features int) { - o.theta = []f64{len: nb_features} - o.bkp_theta = []f64{len: nb_features} +pub fn (mut o ParamsReg) init(nb_features int) { + o.theta = []T{len: nb_features} + o.bkp_theta = []T{len: nb_features} } // backup creates an internal copy of parameters -pub fn (mut o ParamsReg) backup() { +pub fn (mut o ParamsReg) backup() { o.bkp_theta = o.theta.clone() o.bkp_bias = o.bias o.bkp_lambda = o.lambda @@ -44,7 +44,7 @@ pub fn (mut o ParamsReg) backup() { } // restore restores an internal copy of parameters and notifies observers -pub fn (mut o ParamsReg) restore(skip_notification bool) { +pub fn (mut o ParamsReg) restore(skip_notification bool) { o.theta = o.bkp_theta.clone() o.bias = o.bkp_bias o.lambda = o.bkp_lambda @@ -55,7 +55,7 @@ pub fn (mut o ParamsReg) restore(skip_notification bool) { } // set_params sets theta and b and notifies observers -pub fn (mut o ParamsReg) set_params(theta []f64, b f64) { +pub fn (mut o ParamsReg) set_params(theta []T, b T) { o.theta = theta.clone() o.bias = b o.notify_update() @@ -63,7 +63,7 @@ pub fn (mut o ParamsReg) set_params(theta []f64, b f64) { // set_param sets either theta or b (use negative indices for b). Notifies observers // i -- index of theta or -1 for bias -pub fn (mut o ParamsReg) set_param(i int, value f64) { +pub fn (mut o ParamsReg) set_param(i int, value T) { defer { o.notify_update() } @@ -76,7 +76,7 @@ pub fn (mut o ParamsReg) set_param(i int, value f64) { // get_param returns either theta or b (use negative indices for b) // i -- index of theta or -1 for bias -pub fn (o &ParamsReg) get_param(i int) f64 { +pub fn (o &ParamsReg) get_param(i int) T { if i < 0 { return o.bias } @@ -84,66 +84,76 @@ pub fn (o &ParamsReg) get_param(i int) f64 { } // set_thetas sets the whole vector theta and notifies observers -pub fn (mut o ParamsReg) set_thetas(theta []f64) { +pub fn (mut o ParamsReg) set_thetas(theta []T) { la.vector_apply(mut o.theta, 1.0, theta) o.notify_update() } // get_thetas gets a copy of theta -pub fn (o &ParamsReg) get_thetas() []f64 { +pub fn (o &ParamsReg) get_thetas() []T { return o.theta.clone() } // access_thetas returns access (slice) to theta -pub fn (o &ParamsReg) access_thetas() []f64 { +pub fn (o &ParamsReg) access_thetas() []T { return o.theta } // access_bias returns access (pointer) to b -pub fn (o &ParamsReg) access_bias() &f64 { +pub fn (o &ParamsReg) access_bias() &T { return &o.bias } // set_theta sets one component of theta and notifies observers -pub fn (mut o ParamsReg) set_theta(i int, thetai f64) { +pub fn (mut o ParamsReg) set_theta(i int, thetai T) { o.theta[i] = thetai o.notify_update() } // get_theta returns the value of theta[i] -pub fn (o &ParamsReg) get_theta(i int) f64 { +pub fn (o &ParamsReg) get_theta(i int) T { return o.theta[i] } // set_bias sets b and notifies observers -pub fn (mut o ParamsReg) set_bias(b f64) { +pub fn (mut o ParamsReg) set_bias(b T) { o.bias = b o.notify_update() } // get_bias gets a copy of b -pub fn (o &ParamsReg) get_bias() f64 { +pub fn (o &ParamsReg) get_bias() T { return o.bias } // set_lambda sets lambda and notifies observers -pub fn (mut o ParamsReg) set_lambda(lambda f64) { +pub fn (mut o ParamsReg) set_lambda(lambda T) { o.lambda = lambda o.notify_update() } // get_lambda gets a copy of lambda -pub fn (o &ParamsReg) get_lambda() f64 { +pub fn (o &ParamsReg) get_lambda() T { return o.lambda } // set_degree sets p and notifies observers -pub fn (mut o ParamsReg) set_degree(p int) { +pub fn (mut o ParamsReg) set_degree(p int) { o.degree = p o.notify_update() } // get_degree gets a copy of p -pub fn (o &ParamsReg) get_degree() int { +pub fn (o &ParamsReg) get_degree() int { return o.degree } + +// add_observer adds an object to the list of interested observers +pub fn (mut o ParamsReg) add_observer(obs util.Observer) { + o.observable.add_observer(obs) +} + +// notify_update notifies observers of updates +pub fn (mut o ParamsReg) notify_update() { + o.observable.notify_update() +} diff --git a/ml/paramsreg_test.v b/ml/paramsreg_test.v index 092635813..845bafdde 100644 --- a/ml/paramsreg_test.v +++ b/ml/paramsreg_test.v @@ -4,7 +4,7 @@ import vsl.float.float64 fn test_params_reg() { nb_features := 3 - mut params := new_params_reg(nb_features) + mut params := new_params_reg(nb_features) params.theta[0] = 1 params.theta[1] = 2 params.theta[2] = 3 diff --git a/ml/workspace.v b/ml/workspace.v index 6fb7976a5..8927f1bce 100644 --- a/ml/workspace.v +++ b/ml/workspace.v @@ -9,52 +9,52 @@ import vsl.la // NOTE: Stat is an Observer of Data; thus, data.notify_update() will recompute stat // [heap] -pub struct Stat { +pub struct Stat { pub mut: - data &Data // data - name string // name of this object - min_x []f64 // [n_features] min x values - max_x []f64 // [n_features] max x values - sum_x []f64 // [n_features] sum of x values - mean_x []f64 // [n_features] mean of x values - sig_x []f64 // [n_features] standard deviations of x - del_x []f64 // [n_features] difference: max(x) - min(x) - min_y f64 // min of y values - max_y f64 // max of y values - sum_y f64 // sum of y values - mean_y f64 // mean of y values - sig_y f64 // standard deviation of y - del_y f64 // difference: max(y) - min(y) + data &Data // data + name string // name of this object + min_x []T // [n_features] min x values + max_x []T // [n_features] max x values + sum_x []T // [n_features] sum of x values + mean_x []T // [n_features] mean of x values + sig_x []T // [n_features] standard deviations of x + del_x []T // [n_features] difference: max(x) - min(x) + min_y T // min of y values + max_y T // max of y values + sum_y T // sum of y values + mean_y T // mean of y values + sig_y T // standard deviation of y + del_y T // difference: max(y) - min(y) } // stat returns a new Stat object -pub fn stat_from_data(mut data Data, name string) &Stat { - mut o := &Stat{ +pub fn stat_from_data(mut data Data, name string) &Stat { + mut o := &Stat{ name: name data: data - min_x: []f64{len: data.nb_features} - max_x: []f64{len: data.nb_features} - sum_x: []f64{len: data.nb_features} - mean_x: []f64{len: data.nb_features} - sig_x: []f64{len: data.nb_features} - del_x: []f64{len: data.nb_features} + min_x: []T{len: data.nb_features} + max_x: []T{len: data.nb_features} + sum_x: []T{len: data.nb_features} + mean_x: []T{len: data.nb_features} + sig_x: []T{len: data.nb_features} + del_x: []T{len: data.nb_features} } data.add_observer(o) return o } // name returns the name of this stat object (thus defining the Observer interface) -pub fn (o &Stat) name() string { +pub fn (o &Stat) name() string { return o.name } // update compute statistics for given data (an Observer of Data) -pub fn (mut o Stat) update() { +pub fn (mut o Stat) update() { // constants m := o.data.x.m // number of samples n := o.data.x.n // number of features // x values - mf := f64(m) + mf := T(m) for j in 0 .. n { o.min_x[j] = o.data.x.get(0, j) o.max_x[j] = o.min_x[j] @@ -89,8 +89,8 @@ pub fn (mut o Stat) update() { // Output: // t -- scalar t = oᵀy sum of columns of the y vector: t = Σ_i^m o_i y_i // s -- vector s = Xᵀo sum of columns of the X matrix: s_j = Σ_i^m o_i X_ij [n_features] -pub fn (mut o Stat) sum_vars() ([]f64, f64) { - one := []f64{len: o.data.x.m, init: 1.0} +pub fn (mut o Stat) sum_vars() ([]T, T) { + one := []T{len: o.data.x.m, init: T(1)} s := la.matrix_tr_vector_mul(1, o.data.x, one) mut t := 0.0 if o.data.y.len > 0 { @@ -100,7 +100,7 @@ pub fn (mut o Stat) sum_vars() ([]f64, f64) { } // copy_into copies stat into p -pub fn (o &Stat) copy_into(mut p Stat) { +pub fn (o &Stat) copy_into(mut p Stat) { p.min_x = o.min_x.clone() p.max_x = o.max_x.clone() p.sum_x = o.sum_x.clone() @@ -116,10 +116,10 @@ pub fn (o &Stat) copy_into(mut p Stat) { } // str is a custom str function for observers to avoid printing data -pub fn (o &Stat) str() string { +pub fn (o &Stat) str() string { mut res := []string{} - res << 'vsl.ml.Stat{' - res << ' name: $o.name' + res << 'vsl.ml.Stat<$T.name>{' + res << ' name: $o.name' res << ' min_x: $o.min_x' res << ' max_x: $o.max_x' res << ' sum_x: $o.sum_x' diff --git a/ml/workspace_test.v b/ml/workspace_test.v index eacc58e15..952d0208f 100644 --- a/ml/workspace_test.v +++ b/ml/workspace_test.v @@ -2,7 +2,7 @@ module ml import vsl.float.float64 -fn sample_01_check_sat(o &Stat) { +fn sample_01_check_sat(o &Stat) { assert float64.tolerance(o.min_x[0], 0.87, 1e-15) assert float64.tolerance(o.max_x[0], 1.55, 1e-15) assert float64.tolerance(o.mean_x[0], 1.1960, 1e-15)