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

@@ -12,6 +12,8 @@ on:
     - cron: "31 1,12 * * *"
 
   push:
+    branches:
+      - master
 
   pull_request:
     branches:

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<f64>, alpha f64, a la.Matrix<f64>, b la.Matrix<f64>) {
+pub fn matrix_matrix_mul(mut c la.Matrix<f64>, alpha f64, a &la.Matrix<f64>, b &la.Matrix<f64>) {
 	if c.m < 6 && c.n < 6 && a.n < 30 {
 		for i := 0; i < c.m; i++ {
 			for j := 0; j < c.n; j++ {

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<T>(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<T>(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<T>(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<T>(u), 1, arr_to_f64arr<T>(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<T>(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<T>(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<T>(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,58 +124,94 @@ 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<f64>, 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<T>(alpha T, a &Matrix<T>, 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<T>(a.data), a.m, arr_to_f64arr<T>(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]
 			}
 		}
 		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<f64>, 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<T>(alpha T, a &Matrix<T>, 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<T>(a.data), a.n, arr_to_f64arr<T>(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]
 			}
 		}
 		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
 // (e.g. dyadic product)
 //
 // a = alpha⋅u⋅vᵀ    ⇒    aij = alpha * ui * vj
 //
-pub fn vector_vector_tr_mul(alpha f64, u []f64, v []f64) &Matrix<f64> {
-	mut m := new_matrix<f64>(u.len, v.len)
-	if m.m < 9 && m.n < 9 {
+pub fn vector_vector_tr_mul<T>(alpha T, u []T, v []T) &Matrix<T> {
+	$if T is f64 {
+		mut m := new_matrix<f64>(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<T>(u), 1, arr_to_f64arr<T>(v), 1, mut
+			a, int(math.max(m.m, m.n)))
+		return matrix_raw(u.len, v.len, a)
+	} $else {
+		mut m := new_matrix<T>(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])
 			}
 		}
 		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<f64>, u []f64) []f64 {
 //
 pub fn matrix_matrix_mul(mut c Matrix<f64>, alpha f64, a &Matrix<f64>, b &Matrix<f64>) {
 	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<f64>, alpha f64, a &Matrix<f64>, b &Matrix
 //
 pub fn matrix_tr_matrix_mul(mut c Matrix<f64>, alpha f64, a &Matrix<f64>, b &Matrix<f64>) {
 	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++ {

la/matrix.v

@@ -1,8 +1,8 @@
 module la
 
+import math
 import strconv
 import vsl.errors
-import math
 
 [heap]
 pub struct Matrix<T> {
@@ -48,7 +48,7 @@ pub fn matrix_raw<T>(m int, n int, rawdata []T) &Matrix<T> {
 pub fn (mut o Matrix<T>) 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<T>) set_from_deep2(a [][]T) {
 // set_diag sets diagonal matrix with diagonal components equal to val
 pub fn (mut o Matrix<T>) 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<T>) get(i int, j int) T {
 pub fn (o &Matrix<T>) 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<T>) fill(val T) {
 */
 pub fn (mut o Matrix<T>) 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<T>) 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<T>) 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<T>) 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<T>) 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 += ','
 			}

la/matrix_ops.v

@@ -124,8 +124,8 @@ pub fn matrix_inv<T>(mut ai Matrix<T>, mut a Matrix<T>, 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<T>(a.m, a.m)
+	mut vt := new_matrix<T>(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<T>(mut ai Matrix<T>, mut a Matrix<T>, calc_det bool) T {
 // "F" or "" (default) => Frobenius
 pub fn matrix_cond_num<T>(mut a Matrix<T>, normtype string) T {
 	mut res := 0.0
-	mut ai := new_matrix(a.m, a.n)
+	mut ai := new_matrix<T>(a.m, a.n)
 	matrix_inv(mut ai, mut a, false)
 	if normtype == 'I' {
 		res = a.norm_inf() * ai.norm_inf()