mat: make RowView and ColView return Vector type and add RowViewOf and ColViewOf

Also change Dot signature to accept Vector parameters.

AUTHORS

@@ -23,6 +23,7 @@ Fazlul Shahriar <fshahriar@gmail.com>
 Iakov Davydov <iakov.davydov@unil.ch>
 Jalem Raj Rohit <jrajrohit33@gmail.com>
 James Bell <james@stellentus.com>
+James Bowman <james.edward.bowman@gmail.com>
 Jeff Juozapaitis <jjjuozap@email.arizona.edu>
 Jonathan J Lawlor <jonathan.lawlor@gmail.com>
 Jonathan Schroeder <jd.schroeder@gmail.com>

CONTRIBUTORS

@@ -30,6 +30,7 @@ Fazlul Shahriar <fshahriar@gmail.com>
 Iakov Davydov <iakov.davydov@unil.ch>
 Jalem Raj Rohit <jrajrohit33@gmail.com>
 James Bell <james@stellentus.com>
+James Bowman <james.edward.bowman@gmail.com>
 Jeff Juozapaitis <jjjuozap@email.arizona.edu>
 Jonathan J Lawlor <jonathan.lawlor@gmail.com>
 Jonathan Schroeder <jd.schroeder@gmail.com>

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.ColViewOf(dst, j)
-				col.AddScaledVec(col, pt.Coeff, originVec)
+				col.AddScaledVec(&col, pt.Coeff, originVec)
 			}
 		}
 	}

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 {
+func (m *Dense) ColView(j int) Vector {
-	if j >= m.mat.Cols || j < 0 {
+	var v VecDense
-		panic(ErrColAccess)
+	v.ColViewOf(m, j)
-	}
+	return &v
-	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,
-	}
 }
 
 // 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 {
+func (m *Dense) RowView(i int) Vector {
-	if i >= m.mat.Rows || i < 0 {
+	var v VecDense
-		panic(ErrRowAccess)
+	v.RowViewOf(m, i)
-	}
+	return &v
-	return &VecDense{
-		mat: blas64.Vector{
-			Inc:  1,
-			Data: m.rawRowView(i),
-		},
-		n: m.mat.Cols,
-	}
 }
 
 // RawRowView returns a slice backed by the same array as backing the

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)

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.ColViewOf(v, j)
-		cv.ScaleVec(1/blas64.Nrm2(c, cv.mat), cv)
+		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.ColViewOf(dst, j)
-		v.ScaleVec(1/f, v)
+		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
 }

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{},

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
+// A RowViewer can return a Vector reflecting a row that is backed by the matrix
-// data. The VecDense returned will have length equal to the number of columns.
+// 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
+// A ColViewer can return a Vector reflecting a column that is backed by the matrix
-// data. The VecDense returned will have length equal to the number of rows.
+// 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

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()

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
 }

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
+}

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())
 		}