Merge pull request #23 from jonreiter/mulelemvec

add MulElemVec for element-wise multiplication of sparse vectors

Author: james-bowman

cholesky.go

@@ -41,7 +41,7 @@ func (ch *Cholesky) At(i, j int) float64 {
 	ri := ch.chol.RowView(i).(*Vector)
 	rj := ch.chol.RowView(j).(*Vector)
 	// FIXME: check types
-	val = dotSparseSparseNoSortBefore(ri, rj, min(i, j)+1)
+	val = dotSparseSparseNoSortBefore(ri, rj, nil, min(i, j)+1)
 	return val
 }
 
@@ -212,10 +212,6 @@ func cholCSR(matrix *CSR, lower *CSR) {
 		if matrix.RowNNZ(i) == 0 {
 			continue
 		}
-		//		rowDotSum := 0.0
-		//		aPos := 0
-		//		bPos := 0
-		//		thisSum := 0.0
 		for j := 0; j <= i; j++ {
 			iRow := lower.RowView(i)
 			iRowS, iRowIsSparse := iRow.(*Vector)
@@ -231,9 +227,7 @@ func cholCSR(matrix *CSR, lower *CSR) {
 				}
 				lower.Set(j, j, math.Sqrt(matrix.At(i, i)-sum))
 			} else {
-				//				thisSum, _, _ = dotSparseSparseNoSortBeforeWithStart(iRowS, jRowS, j, aPos, bPos)
-				//				rowDotSum = thisSum
-				rowDotSum := dotSparseSparseNoSort(iRowS, jRowS)
+				rowDotSum := dotSparseSparseNoSort(iRowS, jRowS, nil)
 				if rowDotSum == 0.0 && matrix.At(i, j) == 0.0 {
 					continue
 				}

vector.go

@@ -349,7 +349,7 @@ func Dot(a, b mat.Vector) float64 {
 
 	if aIsSparse {
 		if bIsSparse {
-			return dotSparseSparse(as, bs)
+			return dotSparseSparse(as, bs, nil)
 		}
 		if bdense, bIsDense := b.(mat.RawVectorer); bIsDense {
 			raw := bdense.RawVector()
@@ -557,23 +557,23 @@ func (v *Vector) IsSorted() bool {
 
 // dotSparseSparse computes the dot product of two sparse vectors.
 // This will be called by Dot if both entered vectors are Sparse.
-func dotSparseSparse(a, b *Vector) float64 {
+func dotSparseSparse(a, b, c *Vector) float64 {
 	a.Sort()
 	b.Sort()
-	return dotSparseSparseNoSort(a, b)
+	return dotSparseSparseNoSort(a, b, c)
 }
 
-func dotSparseSparseNoSort(a, b *Vector) float64 {
+func dotSparseSparseNoSort(a, b, c *Vector) float64 {
 	n := a.Len()
-	return dotSparseSparseNoSortBefore(a, b, n)
+	return dotSparseSparseNoSortBefore(a, b, c, n)
 }
 
-func dotSparseSparseNoSortBefore(a, b *Vector, n int) float64 {
-	v, _, _ := dotSparseSparseNoSortBeforeWithStart(a, b, n, 0, 0)
+func dotSparseSparseNoSortBefore(a, b, c *Vector, n int) float64 {
+	v, _, _ := dotSparseSparseNoSortBeforeWithStart(a, b, c, n, 0, 0)
 	return v
 }
 
-func dotSparseSparseNoSortBeforeWithStart(a, b *Vector, n, aStart, bStart int) (float64, int, int) {
+func dotSparseSparseNoSortBeforeWithStart(a, b, c *Vector, n, aStart, bStart int) (float64, int, int) {
 	tot := 0.0
 	aPos := aStart
 	bPos := bStart
@@ -583,7 +583,11 @@ func dotSparseSparseNoSortBeforeWithStart(a, b *Vector, n, aStart, bStart int) (
 		aIndex = a.ind[aPos]
 		bIndex = b.ind[bPos]
 		if aIndex == bIndex {
-			tot += a.data[aPos] * b.data[bPos]
+			val := a.data[aPos] * b.data[bPos]
+			tot += val
+			if c != nil {
+				c.SetVec(aIndex, val)
+			}
 			aPos++
 			bPos++
 		} else {
@@ -596,3 +600,21 @@ func dotSparseSparseNoSortBeforeWithStart(a, b *Vector, n, aStart, bStart int) (
 	}
 	return tot, aPos, bPos
 }
+
+// MulElemVec does element-by-element multiplication of a and b
+// and puts the result in the receiver.
+func (v *Vector) MulElemVec(a, b *Vector) {
+	ar := a.Len()
+	br := b.Len()
+	if ar != br {
+		panic(mat.ErrShape)
+	}
+	aNNZ := a.NNZ()
+	bNNZ := b.NNZ()
+	minNNZ := aNNZ
+	if bNNZ < minNNZ {
+		minNNZ = bNNZ
+	}
+	v.reuseAs(ar, minNNZ, true)
+	dotSparseSparse(a, b, v)
+}

vector_test.go

@@ -156,6 +156,11 @@ func TestDot(t *testing.T) {
 			b: NewVector(6, []int{0, 1, 3}, []float64{1, 1, 2}),
 			r: 5,
 		},
+		{
+			a: NewVector(6, []int{4, 3, 1}, []float64{1, 2, 1}),
+			b: NewVector(6, []int{1, 0, 3}, []float64{1, 1, 2}),
+			r: 5,
+		},
 		{
 			a: mat.NewVecDense(6, []float64{0, 1, 0, 2, 1, 0}),
 			b: NewVector(6, []int{0, 1, 3}, []float64{1, 1, 2}),
@@ -171,13 +176,40 @@ func TestDot(t *testing.T) {
 	for ti, test := range tests {
 
 		result := Dot(test.a, test.b)
-
 		if result != test.r {
 			t.Errorf("Test %d: Incorrect result for Dot - expected:\n%v\nbut received:\n%v\n", ti+1, test.r, result)
 		}
 	}
 }
 
+func TestMulElemVec(t *testing.T) {
+	v := NewVector(6, nil, nil)
+	tests := []struct {
+		a *Vector
+		b *Vector
+		r *Vector
+	}{
+		{
+			a: NewVector(6, []int{1, 3, 4}, []float64{1, 2, 1}),
+			b: NewVector(6, []int{0, 1, 3}, []float64{1, 1, 2}),
+			r: NewVector(6, []int{1, 3}, []float64{1, 4}),
+		},
+		{
+			a: NewVector(6, []int{4, 3, 1}, []float64{1, 2, 1}),
+			b: NewVector(6, []int{1, 0, 3}, []float64{1, 1, 2}),
+			r: NewVector(6, []int{1, 3}, []float64{1, 4}),
+		},
+	}
+
+	for ti, test := range tests {
+
+		v.MulElemVec(test.a, test.b)
+		if !mat.Equal(v, test.r) {
+			t.Errorf("Test %d: Incorrect result for MulElemVec - expected:\n%v\nbut received:\n%v\n", ti+1, test.r, v)
+		}
+	}
+}
+
 func TestVectorNorm(t *testing.T) {
 	tests := []struct {
 		a      mat.Vector
@@ -538,7 +570,7 @@ func TestVectorSet(t *testing.T) {
 		},
 	}
 
-	for ti, test := range tests {      
+	for ti, test := range tests {
 		act := new(Vector)
 		act.CloneVec(test.source)
 		act.Set(test.idx, 0, test.val)
@@ -635,9 +667,9 @@ func TestVecSetPanic(t *testing.T) {
 			act.Set(2, col, 1.1)
 
 		}(ti, test.col, test.source)
-  }
+	}
 }
-    
+
 func TestMulMatSparseVec(t *testing.T) {
 	permsB := []struct {
 		name   string
@@ -687,7 +719,7 @@ func TestMulMatSparseVec(t *testing.T) {
 		},
 	}
 
-	for ti, test := range tests {    
+	for ti, test := range tests {
 		for _, b := range permsB {
 			for _, rawa := range matrixPermutationsForA {
 				var matPair = []struct {