forked from pa-m/sklearn
-
Notifications
You must be signed in to change notification settings - Fork 0
/
matrix.go
477 lines (401 loc) · 12.7 KB
/
matrix.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
package base
import (
"bytes"
"fmt"
"io"
"math"
"math/rand"
"runtime"
"sync"
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
"gonum.org/v1/gonum/mat"
)
// MatConst is a matrix where all cless have the same value
type MatConst struct {
Rows, Columns int
Value float64
}
// Dims for MatConst
func (m MatConst) Dims() (int, int) { return m.Rows, m.Columns }
// At for MatConst
func (m MatConst) At(i, j int) float64 { return m.Value }
// T for MatConst
func (m MatConst) T() mat.Matrix { return MatTranspose{Matrix: m} }
// MatTranspose is a matrix override to transpose a mat.Matrix from its initializer
type MatTranspose struct{ mat.Matrix }
// Dims for MatTranspose
func (m MatTranspose) Dims() (int, int) { r, c := m.Matrix.Dims(); return c, r }
// At for MatTranspose
func (m MatTranspose) At(i, j int) float64 {
return m.Matrix.At(j, i)
}
// Set for MatTranspose
func (m MatTranspose) Set(i, j int, v float64) {
if Mutable, ok := m.Matrix.(mat.Mutable); ok {
Mutable.Set(j, i, v)
} else {
panic("underling Matrix is not Mutable")
}
}
// T for MatTranspose
func (m MatTranspose) T() mat.Matrix { return m.Matrix }
// MatOnesPrepended is a matrix override representing its initializer with an initial column of ones added
type MatOnesPrepended struct{ mat.Matrix }
// Dims for MatOnesPrepended
func (m MatOnesPrepended) Dims() (int, int) { r, c := m.Matrix.Dims(); return r, 1 + c }
// At for MatOnesPrepended
func (m MatOnesPrepended) At(i, j int) float64 {
if j == 0 {
return 1.
}
return m.Matrix.At(i, j-1)
}
// Set for MatOnesPrepended
func (m MatOnesPrepended) Set(i, j int, v float64) {
if Mutable, ok := m.Matrix.(mat.Mutable); ok {
Mutable.Set(i, j-1, v)
} else {
panic("underling Matrix is not Mutable")
}
}
// T for MatOnesPrepended not implemented
func (m MatOnesPrepended) T() mat.Matrix { return MatTranspose{m} }
// MatFirstColumnRemoved is a matrix whose an initial column has been removed respective to its initializer
type MatFirstColumnRemoved struct{ mat.Matrix }
// Dims for MatFirstColumnRemoved
func (m MatFirstColumnRemoved) Dims() (int, int) { r, c := m.Matrix.Dims(); return r, c - 1 }
// At for MatFirstColumnRemoved
func (m MatFirstColumnRemoved) At(i, j int) float64 {
return m.Matrix.At(i, j+1)
}
// Set for MatFirstColumnRemoved
func (m MatFirstColumnRemoved) Set(i, j int, v float64) {
if Mutable, ok := m.Matrix.(mat.Mutable); ok {
Mutable.Set(i, j+1, v)
} else {
panic("underling Matrix is not Mutable")
}
}
// T for MatFirstColumnRemoved
func (m MatFirstColumnRemoved) T() mat.Matrix { return MatTranspose{m} }
// MatFirstRowZeroed is a matrix whose an initial Row has been set to zeros respective to its initializer
type MatFirstRowZeroed struct{ mat.Matrix }
// Dims for MatFirstRowZeroed
func (m MatFirstRowZeroed) Dims() (int, int) { return m.Matrix.Dims() }
// At for MatFirstRowZeroed
func (m MatFirstRowZeroed) At(i, j int) float64 {
if i == 0 {
return 0.
}
return m.Matrix.At(i, j)
}
// Set for MatFirstRowZeroed
func (m MatFirstRowZeroed) Set(i, j int, v float64) {
if Mutable, ok := m.Matrix.(mat.Mutable); ok {
Mutable.Set(i, j, v)
} else {
panic("underling Matrix is not Mutable")
}
}
// T for MatFirstRowZeroed
func (m MatFirstRowZeroed) T() mat.Matrix { return MatTranspose{m} }
// MatRowSlice is a matrix row chunk
type MatRowSlice struct {
mat.Matrix
Start, End int
}
// Dims for MatRowSlice
func (m MatRowSlice) Dims() (int, int) { _, c := m.Matrix.Dims(); return m.End - m.Start, c }
// At for MatRowSlice
func (m MatRowSlice) At(i, j int) float64 {
if i < 0 || i > m.End-m.Start {
panic("indexing error")
}
return m.Matrix.At(i+m.Start, j)
}
// Set for MatRowSlice
func (m MatRowSlice) Set(i, j int, v float64) {
if Mutable, ok := m.Matrix.(mat.Mutable); ok {
Mutable.Set(i-m.Start, j, v)
} else {
panic("underling Matrix is not Mutable")
}
}
// T for MatRowSlice
func (m MatRowSlice) T() mat.Matrix { return MatTranspose{m} }
// MatApply0 is a mat.Matrix override where At returns a func-generated value
type MatApply0 struct {
Rows, Columns int
Func func() float64
}
// Dims for MatApply0
func (m MatApply0) Dims() (int, int) { return m.Rows, m.Columns }
// At for MatApply0
func (m MatApply0) At(i, j int) float64 { return m.Func() }
// T for MatApply0
func (m MatApply0) T() mat.Matrix { return MatTranspose{m} }
// MatApply1 is a mat.Matrix override where At returns a func-trancformed value whose inputs are elements from its initializers
type MatApply1 struct {
mat.Matrix
Func func(float64) float64
}
// Dims for MatApply1
func (m MatApply1) Dims() (int, int) { return m.Matrix.Dims() }
// At for MatApply1
func (m MatApply1) At(i, j int) float64 { return m.Func(m.Matrix.At(i, j)) }
// T for MatApply1 returns a MatTranspose
func (m MatApply1) T() mat.Matrix { return MatTranspose{m} }
// MatApply2 is a mat.Matric overrides returning a function where args are elements from its two Matrix initializers
type MatApply2 struct {
A, B mat.Matrix
Func func(a, b float64) float64
}
// Dims for MatApply2
func (m MatApply2) Dims() (int, int) { return m.A.Dims() }
// At for MatApply2
func (m MatApply2) At(i, j int) float64 { return m.Func(m.A.At(i, j), m.B.At(i, j)) }
// T for MatApply2
func (m MatApply2) T() mat.Matrix { return MatTranspose{m} }
// MatSub is a mat.Matrix override returning difference from its two initializers
type MatSub struct{ A, B mat.Matrix }
// Dims for MatSub
func (m MatSub) Dims() (int, int) { return m.A.Dims() }
// At for MatSub
func (m MatSub) At(i, j int) float64 { return m.A.At(i, j) - m.B.At(i, j) }
// T for MatSub
func (m MatSub) T() mat.Matrix { return MatTranspose{m} }
// MatMulElem is a mat.Matrix override returning elementwize product from its two initializers
type MatMulElem struct{ A, B mat.Matrix }
// Dims for MatMuilElem
func (m MatMulElem) Dims() (int, int) { return m.A.Dims() }
// At for MatMulElem
func (m MatMulElem) At(i, j int) float64 { return m.A.At(i, j) * m.B.At(i, j) }
// T for MatMulElem
func (m MatMulElem) T() mat.Matrix { return MatTranspose{m} }
// MatScaled is a mat.Matrix override returning scaled value from its initializer
type MatScaled struct {
mat.Matrix
Scale float64
}
// Dims for MatScaled
func (m MatScaled) Dims() (int, int) { return m.Matrix.Dims() }
// At for MatScaled
func (m MatScaled) At(i, j int) float64 { return m.Matrix.At(i, j) * m.Scale }
// T for MatScaled
func (m MatScaled) T() mat.Matrix { return MatTranspose{m} }
// MatOneMinus is a mat.Matrix override returning 1.-value from its initializer
type MatOneMinus struct {
mat.Matrix
}
// Dims for MatOnesMinus
func (m MatOneMinus) Dims() (int, int) { return m.Dims() }
// At for MatOneMinus
func (m MatOneMinus) At(i, j int) float64 { return 1. - m.At(i, j) }
// T for MatOneMinus
func (m MatOneMinus) T() mat.Matrix { return MatTranspose{m} }
// MatDimsString returns a string representing Dims of its several Matrix parameters
func MatDimsString(mats ...mat.Matrix) string {
s := ""
for _, m := range mats {
r, c := m.Dims()
s = fmt.Sprintf("%s %d,%d", s, r, c)
}
return s
}
// MatDimsCheck checks compat of operator op and its Matrix parameters Dims.
// R is result of op, X and Y are operands of op.
// "." is dot product. "+","-","*","/" are elementwize opts
func MatDimsCheck(op string, R, X, Y mat.Matrix) {
rx, cx := X.Dims()
ry, cy := Y.Dims()
rr, cr := R.Dims()
switch op {
case "+", "-", "*", "/":
if rx != ry || cx != cy || rr != rx || cr != cx {
panic(fmt.Errorf("%s %s", op, MatDimsString(R, X, Y)))
}
case ".":
if cx != ry || rr != rx || cr != cy {
panic(fmt.Errorf("%s %s", op, MatDimsString(R, X, Y)))
}
}
}
// MatStr return a string from a mat.Matrix
func MatStr(Xs ...mat.Matrix) string {
if len(Xs) == 0 {
return ""
}
nSamples, nFeatures := Xs[0].Dims()
b := bytes.NewBuffer(nil)
for i := 0; i < nSamples; i++ {
for imat, X := range Xs {
_, nFeatures = X.Dims()
for j := 0; j < nFeatures; j++ {
io.WriteString(b, fmt.Sprintf("%g", X.At(i, j)))
if j < nFeatures-1 || imat < len(Xs)-1 {
io.WriteString(b, "\t")
} else {
io.WriteString(b, "\n")
}
}
}
}
return b.String()
}
// MatColStr return the string for a matrix column
func MatColStr(X mat.Matrix, j int) string {
nSamples, _ := X.Dims()
var t = make([]float64, nSamples)
mat.Col(t, j, X)
return fmt.Sprint(t)
}
// MatRowStr returns the string for a matrix row
func MatRowStr(X mat.Matrix, i int) string {
_, nFeatures := X.Dims()
var t = make([]float64, nFeatures)
mat.Row(t, i, X)
return fmt.Sprint(t)
}
// MatShuffle shuffles the rows of X and Y matrices
func MatShuffle(X, Y *mat.Dense) {
nSamples, nFeatures := X.Dims()
_, nOutputs := Y.Dims()
Xrowi := make([]float64, nFeatures, nFeatures)
Yrowi := make([]float64, nOutputs, nOutputs)
for i := nSamples - 1; i > 0; i-- {
j := rand.Intn(i + 1)
copy(Xrowi, X.RawRowView(i))
X.SetRow(i, X.RawRowView(j))
X.SetRow(j, Xrowi)
copy(Yrowi, Y.RawRowView(i))
Y.SetRow(i, Y.RawRowView(j))
Y.SetRow(j, Yrowi)
}
}
// MatSigmoid put emelent-wise sigmoid of X into dst
func MatSigmoid(dst *mat.Dense, X mat.Matrix) *mat.Dense {
if dst == nil {
r, c := X.Dims()
dst = mat.NewDense(r, c, nil)
}
dst.Apply(func(i int, j int, v float64) float64 {
return 1. / (1. + math.Exp(-v))
}, X)
return dst
}
// MatParallelGemm parallelize Gemm on A rows
func MatParallelGemm(tA, tB blas.Transpose, alpha float64, a, b blas64.General, beta float64, c blas64.General) {
var start, end int
// defer func() {
// if r := recover(); r != nil {
// fmt.Printf("a %d,%d start:%d end:%d %s\n", a.Rows, a.Cols, start, end, r)
// }
// }()
nJobs := runtime.NumCPU()
// n is a.Rows if blas.NoTrans, else a.Cols
var n int
if tA == blas.NoTrans {
n = a.Rows
} else {
n = a.Cols
}
sliceRows := (n + nJobs - 1) / nJobs
wg := new(sync.WaitGroup)
fn := func(job, start, end int, wg *sync.WaitGroup) {
//MatDimsCheck(".", dst.Slice(start, end, 0, nOutputs), MatRowSlice{Matrix: A, Start: start, End: end}, B)
if end > n {
end = n
}
if end > start {
//fmt.Printf("start %d end %d aRows %d cRows %d\n", start, end, a.Rows, c.Rows)
var aSlice blas64.General
if tA == blas.NoTrans {
aSlice = MatGeneralSlice(a, start, end, 0, a.Cols)
} else {
aSlice = MatGeneralSlice(a, 0, a.Rows, start, end)
}
blas64.Gemm(tA, tB, alpha,
aSlice,
b,
beta,
MatGeneralSlice(c, start, end, 0, c.Cols))
}
wg.Done()
}
if n < 64 {
wg.Add(1)
fn(0, 0, a.Rows, wg)
} else {
for job := 0; job < nJobs; job++ {
end = start + sliceRows
if end > n {
end = n
}
if end > start {
wg.Add(1)
//fmt.Printf("processing rows %d-%d of a(%d,%d, len:%d)\n", start, end, a.Rows, a.Cols, len(a.Data))
go fn(job, start, end, wg)
}
start = end
}
wg.Wait()
}
}
// MatDenseFirstColumnRemoved returns a *mat.Dense with the same underlaying data as M but 1st column removed
func MatDenseFirstColumnRemoved(src *mat.Dense) *mat.Dense {
nSamples, nOutputs := src.Dims()
return MatDenseSlice(src, 0, nSamples, 1, nOutputs)
}
// MatDenseSlice returns a *mat.Dense with the same underlaying data as src but rows and columns removed
func MatDenseSlice(src mat.RawMatrixer, i, k, j, l int) *mat.Dense {
M := src.RawMatrix()
m := &mat.Dense{}
m.SetRawMatrix(MatGeneralSlice(M, i, k, j, l))
return m
}
// MatGeneralSlice returns a blas64.General with the same underlaying data as M but rows and columns removed
func MatGeneralSlice(M blas64.General, i, k, j, l int) blas64.General {
// defer func() {
// if r := recover(); r != nil {
// fmt.Printf("i*M.Stride+j : %d*%d+%d : %d (k-1)*M.Stride+l: %d*%d+%d : %d len:%d %s\n",
// i, M.Stride, j, i*M.Stride+j,
// k-1, M.Stride, l, (k-1)*M.Stride+l, len(M.Data), r)
// }
// }()
if k <= i {
panic(fmt.Errorf("k<=i %d %d", k, i))
}
return blas64.General{
Rows: k - i,
Cols: l - j,
Stride: M.Stride,
Data: M.Data[i*M.Stride+j : (k-1)*M.Stride+l],
}
}
// MatDenseRowSlice returns a *mat.Dense with the same underlaying data as src but rows and columns removed
func MatDenseRowSlice(src mat.RawMatrixer, i, k int) *mat.Dense {
M := src.RawMatrix()
m := &mat.Dense{}
m.SetRawMatrix(MatGeneralRowSlice(M, i, k))
return m
}
// MatGeneralRowSlice returns a blas64.General with the same underlaying data as M but rows and columns removed
func MatGeneralRowSlice(M blas64.General, i, k int) blas64.General {
// defer func() {
// if r := recover(); r != nil {
// fmt.Printf("i*M.Stride+j : %d*%d+%d : %d (k-1)*M.Stride+l: %d*%d+%d : %d len:%d %s\n",
// i, M.Stride, j, i*M.Stride+j,
// k-1, M.Stride, l, (k-1)*M.Stride+l, len(M.Data), r)
// }
// }()
if k <= i {
panic(fmt.Errorf("k<=i %d %d", k, i))
}
return blas64.General{
Rows: k - i,
Cols: M.Cols,
Stride: M.Stride,
Data: M.Data[i*M.Stride : k*M.Stride],
}
}