diff --git a/AUTHORS b/AUTHORS index 41ac8c22..6770c0ce 100644 --- a/AUTHORS +++ b/AUTHORS @@ -23,6 +23,7 @@ Fazlul Shahriar Iakov Davydov Jalem Raj Rohit James Bell +James Bowman Jeff Juozapaitis Jonathan J Lawlor Jonathan Schroeder diff --git a/CONTRIBUTORS b/CONTRIBUTORS index f7eb0c85..13bf3f92 100644 --- a/CONTRIBUTORS +++ b/CONTRIBUTORS @@ -30,6 +30,7 @@ Fazlul Shahriar Iakov Davydov Jalem Raj Rohit James Bell +James Bowman Jeff Juozapaitis Jonathan J Lawlor Jonathan Schroeder diff --git a/diff/fd/jacobian.go b/diff/fd/jacobian.go index 658d1183..d27b99e3 100644 --- a/diff/fd/jacobian.go +++ b/diff/fd/jacobian.go @@ -136,13 +136,14 @@ func jacobianConcurrent(dst *mat.Dense, f func([]float64, []float64), x, origin xcopy := make([]float64, n) y := make([]float64, m) yVec := mat.NewVecDense(m, y) + var col mat.VecDense for job := range jobs { copy(xcopy, x) xcopy[job.j] += job.pt.Loc * step f(y, xcopy) - col := dst.ColView(job.j) + col.ColViewOf(dst, job.j) mu[job.j].Lock() - col.AddScaledVec(col, job.pt.Coeff, yVec) + col.AddScaledVec(&col, job.pt.Coeff, yVec) mu[job.j].Unlock() } } @@ -184,9 +185,10 @@ func jacobianConcurrent(dst *mat.Dense, f func([]float64, []float64), x, origin if pt.Loc != 0 { continue } + var col mat.VecDense for j := 0; j < n; j++ { - col := dst.ColView(j) - col.AddScaledVec(col, pt.Coeff, originVec) + col.ColViewOf(dst, j) + col.AddScaledVec(&col, pt.Coeff, originVec) } } } diff --git a/mat/dense.go b/mat/dense.go index 1fab99c2..9e24008d 100644 --- a/mat/dense.go +++ b/mat/dense.go @@ -201,20 +201,13 @@ func (m *Dense) T() Matrix { return Transpose{m} } -// ColView returns a VecDense reflecting the column j, backed by the matrix data. +// ColView returns a Vector reflecting the column j, backed by the matrix data. // // See ColViewer for more information. -func (m *Dense) ColView(j int) *VecDense { - if j >= m.mat.Cols || j < 0 { - panic(ErrColAccess) - } - return &VecDense{ - mat: blas64.Vector{ - Inc: m.mat.Stride, - Data: m.mat.Data[j : (m.mat.Rows-1)*m.mat.Stride+j+1], - }, - n: m.mat.Rows, - } +func (m *Dense) ColView(j int) Vector { + var v VecDense + v.ColViewOf(m, j) + return &v } // SetCol sets the values in the specified column of the matrix to the values @@ -250,17 +243,10 @@ func (m *Dense) SetRow(i int, src []float64) { // backed by the matrix data. // // See RowViewer for more information. -func (m *Dense) RowView(i int) *VecDense { - if i >= m.mat.Rows || i < 0 { - panic(ErrRowAccess) - } - return &VecDense{ - mat: blas64.Vector{ - Inc: 1, - Data: m.rawRowView(i), - }, - n: m.mat.Cols, - } +func (m *Dense) RowView(i int) Vector { + var v VecDense + v.RowViewOf(m, i) + return &v } // RawRowView returns a slice backed by the same array as backing the diff --git a/mat/dense_test.go b/mat/dense_test.go index 34b37cc3..609764ba 100644 --- a/mat/dense_test.go +++ b/mat/dense_test.go @@ -1212,7 +1212,7 @@ func TestCopyVecDenseAlias(t *testing.T) { 4, 5, 6, 7, 8, 9, }) - var src *VecDense + var src Vector var want *Dense if horiz { src = a.RowView(si) diff --git a/mat/hogsvd.go b/mat/hogsvd.go index 2d3c0403..9da1ddc5 100644 --- a/mat/hogsvd.go +++ b/mat/hogsvd.go @@ -102,9 +102,10 @@ func (gsvd *HOGSVD) Factorize(m ...Matrix) (ok bool) { return false } v := eig.Vectors() + var cv VecDense for j := 0; j < c; j++ { - cv := v.ColView(j) - cv.ScaleVec(1/blas64.Nrm2(c, cv.mat), cv) + cv.ColViewOf(v, j) + cv.ScaleVec(1/blas64.Nrm2(c, cv.mat), &cv) } b := make([]Dense, len(m)) @@ -159,9 +160,10 @@ func (gsvd *HOGSVD) UTo(dst *Dense, n int) *Dense { dst.reuseAs(gsvd.b[n].Dims()) } dst.Copy(&gsvd.b[n]) + var v VecDense for j, f := range gsvd.Values(nil, n) { - v := dst.ColView(j) - v.ScaleVec(1/f, v) + v.ColViewOf(dst, j) + v.ScaleVec(1/f, &v) } return dst } @@ -187,8 +189,10 @@ func (gsvd *HOGSVD) Values(s []float64, n int) []float64 { } else if len(s) != c { panic(ErrSliceLengthMismatch) } + var v VecDense for j := 0; j < c; j++ { - s[j] = blas64.Nrm2(r, gsvd.b[n].ColView(j).mat) + v.ColViewOf(&gsvd.b[n], j) + s[j] = blas64.Nrm2(r, v.mat) } return s } diff --git a/mat/list_test.go b/mat/list_test.go index 5f0cdafa..ae2d535e 100644 --- a/mat/list_test.go +++ b/mat/list_test.go @@ -148,8 +148,19 @@ func legalTypeVec(v Matrix) bool { return ok } -// legalTypesVecVec returns whether both inputs are *VecDense. +// legalTypesVecVec returns whether both inputs are Vector func legalTypesVecVec(a, b Matrix) bool { + if _, ok := a.(Vector); !ok { + return false + } + if _, ok := b.(Vector); !ok { + return false + } + return true +} + +// legalTypesVecDenseVecDense returns whether both inputs are *VecDense. +func legalTypesVecDenseVecDense(a, b Matrix) bool { if _, ok := a.(*VecDense); !ok { return false } @@ -183,7 +194,7 @@ func legalDims(a Matrix, m, n int) bool { return false } return true - case *VecDense: + case *VecDense, *basicVector: if m < 0 || n < 0 { return false } @@ -281,6 +292,20 @@ func makeRandOf(a Matrix, m, n int) Matrix { mat.SetVec(i, rand.NormFloat64()) } return mat + case *basicVector: + if m == 0 && n == 0 { + return &basicVector{} + } + if n != 1 { + panic(fmt.Sprintf("bad vector size: m = %v, n = %v", m, n)) + } + mat := &basicVector{ + m: make([]float64, m), + } + for i := 0; i < m; i++ { + mat.m[i] = rand.NormFloat64() + } + return mat case *SymDense, *basicSymmetric: if m != n { panic("bad size") @@ -373,6 +398,12 @@ func makeCopyOf(a Matrix) Matrix { } copy(m.mat.Data, t.mat.Data) return m + case *basicVector: + m := &basicVector{ + m: make([]float64, t.Len()), + } + copy(m.m, t.m) + return m } } @@ -506,6 +537,7 @@ var testMatrices = []Matrix{ NewTriDense(3, true, nil), NewTriDense(3, false, nil), NewVecDense(0, nil), + &basicVector{}, &VecDense{mat: blas64.Vector{Inc: 10}}, &basicMatrix{}, &basicSymmetric{}, diff --git a/mat/matrix.go b/mat/matrix.go index 6ba51eaf..cc3648a8 100644 --- a/mat/matrix.go +++ b/mat/matrix.go @@ -98,10 +98,10 @@ type Mutable interface { Matrix } -// A RowViewer can return a VecDense reflecting a row that is backed by the matrix -// data. The VecDense returned will have length equal to the number of columns. +// A RowViewer can return a Vector reflecting a row that is backed by the matrix +// data. The Vector returned will have length equal to the number of columns. type RowViewer interface { - RowView(i int) *VecDense + RowView(i int) Vector } // A RawRowViewer can return a slice of float64 reflecting a row that is backed by the matrix @@ -110,10 +110,10 @@ type RawRowViewer interface { RawRowView(i int) []float64 } -// A ColViewer can return a VecDense reflecting a column that is backed by the matrix -// data. The VecDense returned will have length equal to the number of rows. +// A ColViewer can return a Vector reflecting a column that is backed by the matrix +// data. The Vector returned will have length equal to the number of rows. type ColViewer interface { - ColView(j int) *VecDense + ColView(j int) Vector } // A RawColViewer can return a slice of float64 reflecting a column that is backed by the matrix @@ -330,13 +330,22 @@ func Det(a Matrix) float64 { // Dot returns the sum of the element-wise product of a and b. // Dot panics if the matrix sizes are unequal. -func Dot(a, b *VecDense) float64 { +func Dot(a, b Vector) float64 { la := a.Len() lb := b.Len() if la != lb { panic(ErrShape) } - return blas64.Dot(la, a.mat, b.mat) + if arv, ok := a.(RawVectorer); ok { + if brv, ok := b.(RawVectorer); ok { + return blas64.Dot(la, arv.RawVector(), brv.RawVector()) + } + } + var sum float64 + for i := 0; i < la; i++ { + sum += a.At(i, 0) * b.At(i, 0) + } + return sum } // Equal returns whether the matrices a and b have the same size diff --git a/mat/matrix_test.go b/mat/matrix_test.go index 7a576529..14b5f831 100644 --- a/mat/matrix_test.go +++ b/mat/matrix_test.go @@ -367,9 +367,35 @@ func TestDet(t *testing.T) { testOneInputFunc(t, "DetVsChol", f, denseComparison, sameAnswerFloatApproxTol(1e-10), isAnyType, isWide) } +type basicVector struct { + m []float64 +} + +func (v *basicVector) At(r, c int) float64 { + if c != 0 { + panic(ErrColAccess) + } + if r < 0 || r >= v.Len() { + panic(ErrRowAccess) + } + return v.m[r] +} + +func (v *basicVector) Dims() (r, c int) { + return v.Len(), 1 +} + +func (v *basicVector) T() Matrix { + return Transpose{v} +} + +func (v *basicVector) Len() int { + return len(v.m) +} + func TestDot(t *testing.T) { f := func(a, b Matrix) interface{} { - return Dot(a.(*VecDense), b.(*VecDense)) + return Dot(a.(Vector), b.(Vector)) } denseComparison := func(a, b *Dense) interface{} { ra, ca := a.Dims() diff --git a/mat/symmetric.go b/mat/symmetric.go index 98802d94..547a3c89 100644 --- a/mat/symmetric.go +++ b/mat/symmetric.go @@ -495,8 +495,11 @@ func (s *SymDense) PowPSD(a Symmetric, pow float64) error { u.EigenvectorsSym(&eigen) s.SymOuterK(values[0], u.ColView(0)) + + var v VecDense for i := 1; i < dim; i++ { - s.SymRankOne(s, values[i], u.ColView(i)) + v.ColViewOf(&u, i) + s.SymRankOne(s, values[i], &v) } return nil } diff --git a/mat/vector.go b/mat/vector.go index 110581ad..203a97ca 100644 --- a/mat/vector.go +++ b/mat/vector.go @@ -482,3 +482,39 @@ func (v *VecDense) asGeneral() blas64.General { Data: v.mat.Data, } } + +// ColViewOf reflects the column j of the RawMatrixer m, into the receiver +// backed by the same underlying data. The length of the receiver must either be +// zero or match the number of rows in m. +func (v *VecDense) ColViewOf(m RawMatrixer, j int) { + rm := m.RawMatrix() + + if j >= rm.Cols || j < 0 { + panic(ErrColAccess) + } + if !v.IsZero() && v.n != rm.Rows { + panic(ErrShape) + } + + v.mat.Inc = rm.Stride + v.mat.Data = rm.Data[j : (rm.Rows-1)*rm.Stride+j+1] + v.n = rm.Rows +} + +// RowViewOf reflects the row i of the RawMatrixer m, into the receiver +// backed by the same underlying data. The length of the receiver must either be +// zero or match the number of columns in m. +func (v *VecDense) RowViewOf(m RawMatrixer, i int) { + rm := m.RawMatrix() + + if i >= rm.Rows || i < 0 { + panic(ErrRowAccess) + } + if !v.IsZero() && v.n != rm.Cols { + panic(ErrShape) + } + + v.mat.Inc = 1 + v.mat.Data = rm.Data[i*rm.Stride : i*rm.Stride+rm.Cols] + v.n = rm.Cols +} diff --git a/mat/vector_test.go b/mat/vector_test.go index 4f39f71f..23d9a720 100644 --- a/mat/vector_test.go +++ b/mat/vector_test.go @@ -205,7 +205,7 @@ func TestVecDenseMul(t *testing.T) { func TestVecDenseScale(t *testing.T) { for i, test := range []struct { - a *VecDense + a Vector alpha float64 want *VecDense }{ @@ -250,12 +250,12 @@ func TestVecDenseScale(t *testing.T) { }, } { var v VecDense - v.ScaleVec(test.alpha, test.a) + v.ScaleVec(test.alpha, test.a.(*VecDense)) if !reflect.DeepEqual(v.RawVector(), test.want.RawVector()) { t.Errorf("test %d: unexpected result for v = alpha * a: got: %v want: %v", i, v.RawVector(), test.want.RawVector()) } - v.CopyVec(test.a) + v.CopyVec(test.a.(*VecDense)) v.ScaleVec(test.alpha, &v) if !reflect.DeepEqual(v.RawVector(), test.want.RawVector()) { t.Errorf("test %d: unexpected result for v = alpha * v: got: %v want: %v", i, v.RawVector(), test.want.RawVector()) @@ -291,13 +291,13 @@ func TestVecDenseAddScaled(t *testing.T) { sb.Scale(alpha, b) receiver.Add(a, &sb) } - testTwoInput(t, "AddScaledVec", &VecDense{}, method, denseComparison, legalTypesVecVec, legalSizeSameVec, 1e-14) + testTwoInput(t, "AddScaledVec", &VecDense{}, method, denseComparison, legalTypesVecDenseVecDense, legalSizeSameVec, 1e-14) } } func TestVecDenseAdd(t *testing.T) { for i, test := range []struct { - a, b *VecDense + a, b Vector want *VecDense }{ { @@ -317,7 +317,7 @@ func TestVecDenseAdd(t *testing.T) { }, } { var v VecDense - v.AddVec(test.a, test.b) + v.AddVec(test.a.(*VecDense), test.b.(*VecDense)) if !reflect.DeepEqual(v.RawVector(), test.want.RawVector()) { t.Errorf("unexpected result for test %d: got: %v want: %v", i, v.RawVector(), test.want.RawVector()) } @@ -326,7 +326,7 @@ func TestVecDenseAdd(t *testing.T) { func TestVecDenseSub(t *testing.T) { for i, test := range []struct { - a, b *VecDense + a, b Vector want *VecDense }{ { @@ -346,7 +346,7 @@ func TestVecDenseSub(t *testing.T) { }, } { var v VecDense - v.SubVec(test.a, test.b) + v.SubVec(test.a.(*VecDense), test.b.(*VecDense)) if !reflect.DeepEqual(v.RawVector(), test.want.RawVector()) { t.Errorf("unexpected result for test %d: got: %v want: %v", i, v.RawVector(), test.want.RawVector()) } @@ -355,7 +355,7 @@ func TestVecDenseSub(t *testing.T) { func TestVecDenseMulElem(t *testing.T) { for i, test := range []struct { - a, b *VecDense + a, b Vector want *VecDense }{ { @@ -375,7 +375,7 @@ func TestVecDenseMulElem(t *testing.T) { }, } { var v VecDense - v.MulElemVec(test.a, test.b) + v.MulElemVec(test.a.(*VecDense), test.b.(*VecDense)) if !reflect.DeepEqual(v.RawVector(), test.want.RawVector()) { t.Errorf("unexpected result for test %d: got: %v want: %v", i, v.RawVector(), test.want.RawVector()) } @@ -384,7 +384,7 @@ func TestVecDenseMulElem(t *testing.T) { func TestVecDenseDivElem(t *testing.T) { for i, test := range []struct { - a, b *VecDense + a, b Vector want *VecDense }{ { @@ -404,7 +404,7 @@ func TestVecDenseDivElem(t *testing.T) { }, } { var v VecDense - v.DivElemVec(test.a, test.b) + v.DivElemVec(test.a.(*VecDense), test.b.(*VecDense)) if !reflect.DeepEqual(v.RawVector(), test.want.RawVector()) { t.Errorf("unexpected result for test %d: got: %v want: %v", i, v.RawVector(), test.want.RawVector()) }