From 40e217a4945292a90fcb3d1c0d1686abe52001f7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E7=8E=8B=E5=B8=83=E8=A1=A3?= Date: Sat, 25 Feb 2023 01:57:55 +0000 Subject: [PATCH 1/2] =?UTF-8?q?!74=20#I6HES3=20=E5=AE=9E=E7=8E=B0numpy.dot?= =?UTF-8?q?=E5=87=BD=E6=95=B0=20*=20#I6HES3=20=E6=96=B0=E5=A2=9Edot?= =?UTF-8?q?=E5=87=BD=E6=95=B0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- stat/dot.go | 169 ++++++++++++++++++++++++++++++++++++++++++++ stat/dot_test.go | 33 +++++++++ stat/matrix.go | 99 ++++++++++++++++++++++++++ stat/matrix_test.go | 27 +++++++ stat/ones.go | 6 ++ stat/pow.go | 18 +++++ stat/range.go | 18 +++++ stat/shape.go | 23 ++++++ stat/shape_test.go | 12 ++++ stat/util.go | 16 +++++ stat/util_test.go | 22 ++++++ 11 files changed, 443 insertions(+) create mode 100644 stat/dot.go create mode 100644 stat/dot_test.go create mode 100644 stat/matrix.go create mode 100644 stat/matrix_test.go create mode 100644 stat/ones.go create mode 100644 stat/pow.go create mode 100644 stat/range.go create mode 100644 stat/shape.go create mode 100644 stat/shape_test.go create mode 100644 stat/util.go create mode 100644 stat/util_test.go diff --git a/stat/dot.go b/stat/dot.go new file mode 100644 index 0000000..966e26a --- /dev/null +++ b/stat/dot.go @@ -0,0 +1,169 @@ +package stat + +import ( + "fmt" +) + +func Dot1D[T Number](a, b []T) T { + return __dot1d_go(a, b) +} + +func __dot1d_go[T Number](x, y []T) T { + res := T(0) + for i := 0; i < len(x); i++ { + res += x[i] * y[i] + } + return res +} + +// 3 x 3 +func __dot2d_go[T Number](a, b [][]T) [][]T { + A := a + B := b + rLen := len(B) + cLen := len(B[0]) + c := make([][]T, rLen) + for i := 0; i < rLen; i++ { + col := make([]T, cLen) + for j := 0; j < cLen; j++ { + for k := 0; k < rLen; k++ { + col[j] += A[i][k] * B[k][j] + } + } + c[i] = col + } + return c +} + +func __align_go[T Number](x [][]T, n int) [][]T { + //if n < 3 { + // panic("lt 3") + //} + d := make([][]T, n) + rLen := len(x) + cLen := len(x[0]) + for i := 0; i < n; i++ { + col := make([]T, n) + for j := 0; j < n; j++ { + col[j] = x[rLen-n+i][cLen-n+j] + } + d[i] = col + } + fmt.Println(d) + return d +} + +func Dot2D_V1[T Number](a, b [][]T) [][]T { + W := 3 + A := __align_go(a, W) + B := __align_go(b, W) + + return __dot2d_go(A, B) +} + +// Dot2D 二维矩阵点积 +// +// 点积(dot)运算及简单应用 https://www.jianshu.com/p/482abac8798c +func Dot2D[T Number](a, b [][]T) [][]T { + A := a + B := b + rLen := len(A[0]) + cLen := len(B[0]) + xLen := __min_n_go(rLen, cLen) + x := make([][]T, xLen) + // 行 + for i := 0; i < xLen; i++ { + col := make([]T, cLen) + // 列 + for j := 0; j < cLen; j++ { + for k := 0; k < rLen; k++ { + col[j] += A[i][k] * B[k][j] + } + } + x[i] = col + } + return x +} + +// Dot 二维点积 +func Dot[T Number](a, b [][]T) [][]T { + m, n := Shape[T](a) + k, l := Shape[T](b) + if n != k { + panic("dot 2d a.rows<>b.cols") + } + //fmt.Println("m, n:", m, n) + //fmt.Println("k, l:", k, l) + x := make([][]T, m) + // 行 + for i := 0; i < m; i++ { + col := make([]T, l) + // 列 + for c := 0; c < l; c++ { + for r := 0; r < k; r++ { + col[c] += a[i][r] * b[r][c] + } + } + x[i] = col + } + return x +} + +func Dot_v1[T Number](a, b [][]T) [][]T { + m, n := Shape[T](a) + k, l := Shape[T](b) + if n != k { + panic("dot 2d vs 1d a.rows<>b.cols") + } + //fmt.Println("m, n:", m, n) + //fmt.Println("k, l:", k, l) + x := make([][]T, m) + // 行 + for i := 0; i < m; i++ { + col := make([]T, l) + // 列 + for c := 0; c < l; c++ { + for r := 0; r < k; r++ { + col[c] += a[i][r] * b[r][c] + } + } + x[i] = col + } + return x +} + +// Dot2D1 二维矩阵和一维矩阵计算点积 +func Dot2D1[T Number](a [][]T, b []T) []T { + B := [][]T{b} + b1 := Transpose2D(B) + //fmt.Println("Dot2D1: b1 =", b1) + x1 := Dot[T](a, b1) + //fmt.Println("Dot2D1: x1 =", x1) + x2 := Transpose2D(x1) + //fmt.Println("Dot2D1: x2 =", x2) + return x2[0] +} + +func Dot2D1_v2[T Number](a [][]T, b []T) []T { + m, n := Shape[T](a) + k, l := Shape[T](b) + if l < 1 { + l = 1 + } + + fmt.Println("m, n:", m, n) + fmt.Println("k, l:", k, l) + x := make( /*[]*/ []T, m) + // 行 + for i := 0; i < m; i++ { + col := T(0) + // 列 + for c := 0; c < l; c++ { + for r := 0; r < k; r++ { + col += a[i][r] * b[r] + } + } + x[i] = col + } + return x +} diff --git a/stat/dot_test.go b/stat/dot_test.go new file mode 100644 index 0000000..9335fda --- /dev/null +++ b/stat/dot_test.go @@ -0,0 +1,33 @@ +package stat + +import ( + "fmt" + "testing" +) + +func TestDot2D(t *testing.T) { + // https://blog.csdn.net/llittleSun/article/details/115045660 + // + //a := [][]int{{1, 2}, {3, 7}} + //b := [][]int{{4, 3}, {5, 0}} + // + //c := Dot2D(a, b) + //fmt.Println(c) + + A := [][]int{{1, 4, 9}, {1, 2, 3}, {1, 1, 1}} + fmt.Println("A =", A) + B := [][]int{{1, 1, 1}, {4, 2, 1}, {9, 3, 1}} + fmt.Println("B =", B) + C := Dot2D(A, B) + fmt.Println("C =", C) +} + +func TestDot(t *testing.T) { + A := [][]int{{1, 4, 9}, {1, 2, 3}, {1, 1, 1}} + A = [][]int{{1, 4, 9, 16, 25}, {1, 2, 3, 4, 5}, {1, 1, 1, 1, 1}} + B := Transpose2D(A) + fmt.Println("A =", A) + fmt.Println("B =", B) + C := Dot(A, B) + fmt.Println("C =", C) +} diff --git a/stat/matrix.go b/stat/matrix.go new file mode 100644 index 0000000..522a800 --- /dev/null +++ b/stat/matrix.go @@ -0,0 +1,99 @@ +package stat + +// Concat1D +// +// Translates slice objects to concatenation along the second axis. +// 沿第二个轴将切片对象转换为串联 +func Concat1D[T Number](a, b, c []T) [][]T { + length := len(a) + cLen := 3 // a,b,c TODO:这个位置需要扩展, 如果传入的是2或者3个以上的数组怎么处理? + rows := make([][]T, length) + for i := 0; i < length; i++ { + col := make([]T, cLen) + col[0] = a[i] + col[1] = b[i] + col[2] = c[i] + rows[i] = col + } + return rows +} + +// Transpose2D 矩阵转置 +func Transpose2D[T Number](x [][]T) [][]T { + length := len(x[0]) + cLen := len(x) + rows := make([][]T, length) + for i := 0; i < length; i++ { + col := make([]T, cLen) + for j := 0; j < cLen; j++ { + col[j] = x[j][i] + } + rows[i] = col + } + return rows +} + +// Inverse 计算矩阵的（乘法）逆 +// +// Compute the (multiplicative) inverse of a matrix. +// +// Given a square matrix `a`, return the matrix `ainv` satisfying +// ``dot(a, ainv) = dot(ainv, a) = eye(a.shape[0])``. +func Inverse(a [][]float64) [][]float64 { + var n = len(a) + + // Create augmented matrix + var augmented = make([][]float64, n) + for i := range augmented { + augmented[i] = make([]float64, 2*n) + for j := 0; j < n; j++ { + augmented[i][j] = a[i][j] + } + } + for i := 0; i < n; i++ { + augmented[i][i+n] = 1 + } + + // Perform Gaussian elimination + for i := 0; i < n; i++ { + var pivot = augmented[i][i] + for j := i + 1; j < n; j++ { + var factor = augmented[j][i] / pivot + for k := i; k < 2*n; k++ { + augmented[j][k] -= factor * augmented[i][k] + } + } + } + + // Perform back-substitution + for i := n - 1; i >= 0; i-- { + var pivot = augmented[i][i] + for j := i - 1; j >= 0; j-- { + var factor = augmented[j][i] / pivot + for k := i; k < 2*n; k++ { + augmented[j][k] -= factor * augmented[i][k] + } + } + } + + // Normalize rows + for i := 0; i < n; i++ { + var pivot = augmented[i][i] + for j := 0; j < 2*n; j++ { + augmented[i][j] /= pivot + } + } + + // Extract inverse from augmented matrix + var inverse = make([][]float64, n) + for i := range inverse { + inverse[i] = make([]float64, n) + } + for i := 0; i < n; i++ { + for j := 0; j < n; j++ { + inverse[i][j] = augmented[i][j+n] + } + } + + return inverse +} diff --git a/stat/matrix_test.go b/stat/matrix_test.go new file mode 100644 index 0000000..827a344 --- /dev/null +++ b/stat/matrix_test.go @@ -0,0 +1,27 @@ +package stat + +import ( + "fmt" + "testing" +) + +func TestConcat1D(t *testing.T) { + a := []int{1, 2, 3, 4, 5} + b := Mul(a, a) + c := []int{1, 1, 1, 1, 1} + fmt.Println(Concat1D[int](b, a, c)) +} + +func TestTranspose2D(t *testing.T) { + N := 5 + x := Arange[float64](1, float64(N)+1, 1) + + t1 := Pow[float64](x, 2) + t2 := x + t3 := Ones[float64](x) + + A := Concat1D(t1, t2, t3) + T := Transpose2D(A) + fmt.Println("A =", A) + fmt.Println("T =", T) +} diff --git a/stat/ones.go b/stat/ones.go new file mode 100644 index 0000000..3001398 --- /dev/null +++ b/stat/ones.go @@ -0,0 +1,6 @@ +package stat + +// Ones v -> shape +func Ones[T Number](v []T) []T { + return Repeat[T](T(1), len(v)) +} diff --git a/stat/pow.go b/stat/pow.go new file mode 100644 index 0000000..6f42b15 --- /dev/null +++ b/stat/pow.go @@ -0,0 +1,18 @@ +package stat + +import ( + "math" +) + +func Pow[T Number](v []T, n int) []T { + x := make([]T, len(v)) + for i := 0; i < len(v); i++ { + x[i] = __pow_go(v[i], n) + } + return x +} + +func __pow_go[T Number](x T, n int) T { + y := math.Pow(float64(x), float64(n)) + return T(y) +} diff --git a/stat/range.go b/stat/range.go new file mode 100644 index 0000000..6abaeb6 --- /dev/null +++ b/stat/range.go @@ -0,0 +1,18 @@ +package stat + +// Arange Return evenly spaced values within a given interval. +// +// 返回给定间隔内的等间距值 +func Arange[T Number](start T, end T, argv ...T) []T { + step := T(1) + if len(argv) > 0 { + step = argv[0] + } + x := make([]T, 0) + for i := start; i < end; i += step { + x = append(x, start) + start += T(step) + } + + return x +} diff --git a/stat/shape.go b/stat/shape.go new file mode 100644 index 0000000..32de271 --- /dev/null +++ b/stat/shape.go @@ -0,0 +1,23 @@ +package stat + +func Shape[T Number](x any) (r, c int) { + return __slice_shape_go[T](x) +} + +func __slice_shape_go[T Number](x any) (r, c int) { + switch vs := x.(type) { + case T: + return 0, 0 + case []T: + return len(vs), 0 + case [][]T: + r = len(vs) + if r > 0 { + c = len(vs[0]) + } + return + default: + return -1, -1 + } + +} diff --git a/stat/shape_test.go b/stat/shape_test.go new file mode 100644 index 0000000..facfa5d --- /dev/null +++ b/stat/shape_test.go @@ -0,0 +1,12 @@ +package stat + +import ( + "fmt" + "testing" +) + +func Test___slice_shape_go(t *testing.T) { + fmt.Println(__slice_shape_go[int](0)) + fmt.Println(__slice_shape_go[int]([]int{1, 2, 3})) + fmt.Println(__slice_shape_go[int]([][]int{{1, 2, 3}, {4, 5, 6}})) +} diff --git a/stat/util.go b/stat/util.go new file mode 100644 index 0000000..f3135db --- /dev/null +++ b/stat/util.go @@ -0,0 +1,16 @@ +package stat + +import ( + "math" +) + +func __min_n_go[T Number](a, b T) T { + if a < b { + return a + } + return b +} + +func __float_exp_go(f float64) (frac float64, exp int) { + return math.Frexp(f) +} diff --git a/stat/util_test.go b/stat/util_test.go new file mode 100644 index 0000000..d688bdf --- /dev/null +++ b/stat/util_test.go @@ -0,0 +1,22 @@ +package stat + +import ( + "fmt" + "strconv" + "testing" +) + +func Test___index_go(t *testing.T) { + f1 := 1.7763568394002505e-15 + fmt.Println(f1) + fmt.Printf("1.7763568394002505e-15: %f\n", f1) + fmt.Println(__float_exp_go(f1)) + f2 := -6.21724894e-15 + fmt.Println(__float_exp_go(f2)) + fmt.Printf("-6.21724894e-15: %f\n", f2) + + f2 = -0.07142857142857112 + fmt.Println(__float_exp_go(f2)) + s2 := strconv.FormatFloat(f2, 'E', -1, 64) + fmt.Println(s2) +} -- Gitee From e287b84f9442180879603f37c21fa8510aac0c52 Mon Sep 17 00:00:00 2001 From: wangfeng Date: Sat, 25 Feb 2023 10:00:48 +0800 Subject: [PATCH 2/2] =?UTF-8?q?#I6HES2=20=E4=BF=AE=E8=AE=A2shape=E5=87=BD?= =?UTF-8?q?=E6=95=B0?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- stat/dot.go | 5 +---- stat/shape.go | 1 + 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/stat/dot.go b/stat/dot.go index 966e26a..1f0bcdc 100644 --- a/stat/dot.go +++ b/stat/dot.go @@ -136,11 +136,8 @@ func Dot_v1[T Number](a, b [][]T) [][]T { func Dot2D1[T Number](a [][]T, b []T) []T { B := [][]T{b} b1 := Transpose2D(B) - //fmt.Println("Dot2D1: b1 =", b1) x1 := Dot[T](a, b1) - //fmt.Println("Dot2D1: x1 =", x1) x2 := Transpose2D(x1) - //fmt.Println("Dot2D1: x2 =", x2) return x2[0] } @@ -153,7 +150,7 @@ func Dot2D1_v2[T Number](a [][]T, b []T) []T { fmt.Println("m, n:", m, n) fmt.Println("k, l:", k, l) - x := make( /*[]*/ []T, m) + x := make([]T, m) // 行 for i := 0; i < m; i++ { col := T(0) diff --git a/stat/shape.go b/stat/shape.go index 32de271..f0819f6 100644 --- a/stat/shape.go +++ b/stat/shape.go @@ -1,5 +1,6 @@ package stat +// Shape 返回一维或2维数组的行数和列数 func Shape[T Number](x any) (r, c int) { return __slice_shape_go[T](x) } -- Gitee