--- /dev/null
+// Copyright ©2013 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mat
+
+import (
+ "fmt"
+ "math"
+ "reflect"
+ "testing"
+
+ "gonum.org/v1/gonum/blas"
+ "gonum.org/v1/gonum/blas/blas64"
+ "gonum.org/v1/gonum/floats"
+)
+
+func panics(fn func()) (panicked bool, message string) {
+ defer func() {
+ r := recover()
+ panicked = r != nil
+ message = fmt.Sprint(r)
+ }()
+ fn()
+ return
+}
+
+func flatten(f [][]float64) (r, c int, d []float64) {
+ r = len(f)
+ if r == 0 {
+ panic("bad test: no row")
+ }
+ c = len(f[0])
+ d = make([]float64, 0, r*c)
+ for _, row := range f {
+ if len(row) != c {
+ panic("bad test: ragged input")
+ }
+ d = append(d, row...)
+ }
+ return r, c, d
+}
+
+func unflatten(r, c int, d []float64) [][]float64 {
+ m := make([][]float64, r)
+ for i := 0; i < r; i++ {
+ m[i] = d[i*c : (i+1)*c]
+ }
+ return m
+}
+
+// eye returns a new identity matrix of size n×n.
+func eye(n int) *Dense {
+ d := make([]float64, n*n)
+ for i := 0; i < n*n; i += n + 1 {
+ d[i] = 1
+ }
+ return NewDense(n, n, d)
+}
+
+func TestCol(t *testing.T) {
+ for id, af := range [][][]float64{
+ {
+ {1, 2, 3},
+ {4, 5, 6},
+ {7, 8, 9},
+ },
+ {
+ {1, 2, 3},
+ {4, 5, 6},
+ {7, 8, 9},
+ {10, 11, 12},
+ },
+ {
+ {1, 2, 3, 4},
+ {5, 6, 7, 8},
+ {9, 10, 11, 12},
+ },
+ } {
+ a := NewDense(flatten(af))
+ col := make([]float64, a.mat.Rows)
+ for j := range af[0] {
+ for i := range col {
+ col[i] = float64(i*a.mat.Cols + j + 1)
+ }
+
+ if got := Col(nil, j, a); !reflect.DeepEqual(got, col) {
+ t.Errorf("test %d: unexpected values returned for dense col %d: got: %v want: %v",
+ id, j, got, col)
+ }
+
+ got := make([]float64, a.mat.Rows)
+ if Col(got, j, a); !reflect.DeepEqual(got, col) {
+ t.Errorf("test %d: unexpected values filled for dense col %d: got: %v want: %v",
+ id, j, got, col)
+ }
+ }
+ }
+
+ denseComparison := func(a *Dense) interface{} {
+ r, c := a.Dims()
+ ans := make([][]float64, c)
+ for j := range ans {
+ ans[j] = make([]float64, r)
+ for i := range ans[j] {
+ ans[j][i] = a.At(i, j)
+ }
+ }
+ return ans
+ }
+
+ f := func(a Matrix) interface{} {
+ _, c := a.Dims()
+ ans := make([][]float64, c)
+ for j := range ans {
+ ans[j] = Col(nil, j, a)
+ }
+ return ans
+ }
+ testOneInputFunc(t, "Col", f, denseComparison, sameAnswerF64SliceOfSlice, isAnyType, isAnySize)
+
+ f = func(a Matrix) interface{} {
+ r, c := a.Dims()
+ ans := make([][]float64, c)
+ for j := range ans {
+ ans[j] = make([]float64, r)
+ Col(ans[j], j, a)
+ }
+ return ans
+ }
+ testOneInputFunc(t, "Col", f, denseComparison, sameAnswerF64SliceOfSlice, isAnyType, isAnySize)
+}
+
+func TestRow(t *testing.T) {
+ for id, af := range [][][]float64{
+ {
+ {1, 2, 3},
+ {4, 5, 6},
+ {7, 8, 9},
+ },
+ {
+ {1, 2, 3},
+ {4, 5, 6},
+ {7, 8, 9},
+ {10, 11, 12},
+ },
+ {
+ {1, 2, 3, 4},
+ {5, 6, 7, 8},
+ {9, 10, 11, 12},
+ },
+ } {
+ a := NewDense(flatten(af))
+ for i, row := range af {
+ if got := Row(nil, i, a); !reflect.DeepEqual(got, row) {
+ t.Errorf("test %d: unexpected values returned for dense row %d: got: %v want: %v",
+ id, i, got, row)
+ }
+
+ got := make([]float64, len(row))
+ if Row(got, i, a); !reflect.DeepEqual(got, row) {
+ t.Errorf("test %d: unexpected values filled for dense row %d: got: %v want: %v",
+ id, i, got, row)
+ }
+ }
+ }
+
+ denseComparison := func(a *Dense) interface{} {
+ r, c := a.Dims()
+ ans := make([][]float64, r)
+ for i := range ans {
+ ans[i] = make([]float64, c)
+ for j := range ans[i] {
+ ans[i][j] = a.At(i, j)
+ }
+ }
+ return ans
+ }
+
+ f := func(a Matrix) interface{} {
+ r, _ := a.Dims()
+ ans := make([][]float64, r)
+ for i := range ans {
+ ans[i] = Row(nil, i, a)
+ }
+ return ans
+ }
+ testOneInputFunc(t, "Row", f, denseComparison, sameAnswerF64SliceOfSlice, isAnyType, isAnySize)
+
+ f = func(a Matrix) interface{} {
+ r, c := a.Dims()
+ ans := make([][]float64, r)
+ for i := range ans {
+ ans[i] = make([]float64, c)
+ Row(ans[i], i, a)
+ }
+ return ans
+ }
+ testOneInputFunc(t, "Row", f, denseComparison, sameAnswerF64SliceOfSlice, isAnyType, isAnySize)
+}
+
+func TestCond(t *testing.T) {
+ for i, test := range []struct {
+ a *Dense
+ condOne float64
+ condTwo float64
+ condInf float64
+ }{
+ {
+ a: NewDense(3, 3, []float64{
+ 8, 1, 6,
+ 3, 5, 7,
+ 4, 9, 2,
+ }),
+ condOne: 16.0 / 3.0,
+ condTwo: 4.330127018922192,
+ condInf: 16.0 / 3.0,
+ },
+ {
+ a: NewDense(4, 4, []float64{
+ 2, 9, 3, 2,
+ 10, 9, 9, 3,
+ 1, 1, 5, 2,
+ 8, 4, 10, 2,
+ }),
+ condOne: 1 / 0.024740155174938,
+ condTwo: 34.521576567075087,
+ condInf: 1 / 0.012034465570035,
+ },
+ {
+ a: NewDense(3, 3, []float64{
+ 5, 6, 7,
+ 8, -2, 1,
+ 7, 7, 7}),
+ condOne: 30.769230769230749,
+ condTwo: 21.662689498448440,
+ condInf: 31.153846153846136,
+ },
+ } {
+ orig := DenseCopyOf(test.a)
+ condOne := Cond(test.a, 1)
+ if !floats.EqualWithinAbsOrRel(test.condOne, condOne, 1e-13, 1e-13) {
+ t.Errorf("Case %d: one norm mismatch. Want %v, got %v", i, test.condOne, condOne)
+ }
+ if !Equal(test.a, orig) {
+ t.Errorf("Case %d: unexpected mutation of input matrix for one norm. Want %v, got %v", i, orig, test.a)
+ }
+ condTwo := Cond(test.a, 2)
+ if !floats.EqualWithinAbsOrRel(test.condTwo, condTwo, 1e-13, 1e-13) {
+ t.Errorf("Case %d: two norm mismatch. Want %v, got %v", i, test.condTwo, condTwo)
+ }
+ if !Equal(test.a, orig) {
+ t.Errorf("Case %d: unexpected mutation of input matrix for two norm. Want %v, got %v", i, orig, test.a)
+ }
+ condInf := Cond(test.a, math.Inf(1))
+ if !floats.EqualWithinAbsOrRel(test.condInf, condInf, 1e-13, 1e-13) {
+ t.Errorf("Case %d: inf norm mismatch. Want %v, got %v", i, test.condInf, condInf)
+ }
+ if !Equal(test.a, orig) {
+ t.Errorf("Case %d: unexpected mutation of input matrix for inf norm. Want %v, got %v", i, orig, test.a)
+ }
+ }
+
+ for _, test := range []struct {
+ name string
+ norm float64
+ }{
+ {
+ name: "CondOne",
+ norm: 1,
+ },
+ {
+ name: "CondTwo",
+ norm: 2,
+ },
+ {
+ name: "CondInf",
+ norm: math.Inf(1),
+ },
+ } {
+ f := func(a Matrix) interface{} {
+ return Cond(a, test.norm)
+ }
+ denseComparison := func(a *Dense) interface{} {
+ return Cond(a, test.norm)
+ }
+ testOneInputFunc(t, test.name, f, denseComparison, sameAnswerFloatApproxTol(1e-12), isAnyType, isAnySize)
+ }
+}
+
+func TestDet(t *testing.T) {
+ for c, test := range []struct {
+ a *Dense
+ ans float64
+ }{
+ {
+ a: NewDense(2, 2, []float64{1, 0, 0, 1}),
+ ans: 1,
+ },
+ {
+ a: NewDense(2, 2, []float64{1, 0, 0, -1}),
+ ans: -1,
+ },
+ {
+ a: NewDense(3, 3, []float64{
+ 1, 2, 0,
+ 0, 1, 2,
+ 0, 2, 1,
+ }),
+ ans: -3,
+ },
+ {
+ a: NewDense(3, 3, []float64{
+ 1, 2, 3,
+ 5, 7, 9,
+ 6, 9, 12,
+ }),
+ ans: 0,
+ },
+ } {
+ a := DenseCopyOf(test.a)
+ det := Det(a)
+ if !Equal(a, test.a) {
+ t.Errorf("Input matrix changed during Det. Case %d.", c)
+ }
+ if !floats.EqualWithinAbsOrRel(det, test.ans, 1e-14, 1e-14) {
+ t.Errorf("Det mismatch case %d. Got %v, want %v", c, det, test.ans)
+ }
+ }
+ // Perform the normal list test to ensure it works for all types.
+ f := func(a Matrix) interface{} {
+ return Det(a)
+ }
+ denseComparison := func(a *Dense) interface{} {
+ return Det(a)
+ }
+ testOneInputFunc(t, "Det", f, denseComparison, sameAnswerFloatApproxTol(1e-12), isAnyType, isSquare)
+
+ // Check that it gives approximately the same answer as Cholesky
+ // Ensure the input matrices are wider than tall so they are full rank
+ isWide := func(ar, ac int) bool {
+ return ar <= ac
+ }
+ f = func(a Matrix) interface{} {
+ ar, ac := a.Dims()
+ if !isWide(ar, ac) {
+ panic(ErrShape)
+ }
+ var tmp Dense
+ tmp.Mul(a, a.T())
+ return Det(&tmp)
+ }
+ denseComparison = func(a *Dense) interface{} {
+ ar, ac := a.Dims()
+ if !isWide(ar, ac) {
+ panic(ErrShape)
+ }
+ var tmp SymDense
+ tmp.SymOuterK(1, a)
+ var chol Cholesky
+ ok := chol.Factorize(&tmp)
+ if !ok {
+ panic("bad chol test")
+ }
+ return chol.Det()
+ }
+ testOneInputFunc(t, "DetVsChol", f, denseComparison, sameAnswerFloatApproxTol(1e-10), isAnyType, isWide)
+}
+
+type basicVector struct {
+ m []float64
+}
+
+func (v *basicVector) AtVec(i int) float64 {
+ if i < 0 || i >= v.Len() {
+ panic(ErrRowAccess)
+ }
+ return v.m[i]
+}
+
+func (v *basicVector) At(r, c int) float64 {
+ if c != 0 {
+ panic(ErrColAccess)
+ }
+ return v.AtVec(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.(Vector), b.(Vector))
+ }
+ denseComparison := func(a, b *Dense) interface{} {
+ ra, ca := a.Dims()
+ rb, cb := b.Dims()
+ if ra != rb || ca != cb {
+ panic(ErrShape)
+ }
+ var sum float64
+ for i := 0; i < ra; i++ {
+ for j := 0; j < ca; j++ {
+ sum += a.At(i, j) * b.At(i, j)
+ }
+ }
+ return sum
+ }
+ testTwoInputFunc(t, "Dot", f, denseComparison, sameAnswerFloatApproxTol(1e-12), legalTypesVectorVector, legalSizeSameVec)
+}
+
+func TestEqual(t *testing.T) {
+ f := func(a, b Matrix) interface{} {
+ return Equal(a, b)
+ }
+ denseComparison := func(a, b *Dense) interface{} {
+ return Equal(a, b)
+ }
+ testTwoInputFunc(t, "Equal", f, denseComparison, sameAnswerBool, legalTypesAll, isAnySize2)
+}
+
+func TestMax(t *testing.T) {
+ // A direct test of Max with *Dense arguments is in TestNewDense.
+ f := func(a Matrix) interface{} {
+ return Max(a)
+ }
+ denseComparison := func(a *Dense) interface{} {
+ return Max(a)
+ }
+ testOneInputFunc(t, "Max", f, denseComparison, sameAnswerFloat, isAnyType, isAnySize)
+}
+
+func TestMin(t *testing.T) {
+ // A direct test of Min with *Dense arguments is in TestNewDense.
+ f := func(a Matrix) interface{} {
+ return Min(a)
+ }
+ denseComparison := func(a *Dense) interface{} {
+ return Min(a)
+ }
+ testOneInputFunc(t, "Min", f, denseComparison, sameAnswerFloat, isAnyType, isAnySize)
+}
+
+func TestNorm(t *testing.T) {
+ for i, test := range []struct {
+ a [][]float64
+ ord float64
+ norm float64
+ }{
+ {
+ a: [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}},
+ ord: 1,
+ norm: 30,
+ },
+ {
+ a: [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}},
+ ord: 2,
+ norm: 25.495097567963924,
+ },
+ {
+ a: [][]float64{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}, {10, 11, 12}},
+ ord: math.Inf(1),
+ norm: 33,
+ },
+ {
+ a: [][]float64{{1, -2, -2}, {-4, 5, 6}},
+ ord: 1,
+ norm: 8,
+ },
+ {
+ a: [][]float64{{1, -2, -2}, {-4, 5, 6}},
+ ord: math.Inf(1),
+ norm: 15,
+ },
+ } {
+ a := NewDense(flatten(test.a))
+ if math.Abs(Norm(a, test.ord)-test.norm) > 1e-14 {
+ t.Errorf("Mismatch test %d: %v norm = %f", i, test.a, test.norm)
+ }
+ }
+
+ for _, test := range []struct {
+ name string
+ norm float64
+ }{
+ {"NormOne", 1},
+ {"NormTwo", 2},
+ {"NormInf", math.Inf(1)},
+ } {
+ f := func(a Matrix) interface{} {
+ return Norm(a, test.norm)
+ }
+ denseComparison := func(a *Dense) interface{} {
+ return Norm(a, test.norm)
+ }
+ testOneInputFunc(t, test.name, f, denseComparison, sameAnswerFloatApproxTol(1e-12), isAnyType, isAnySize)
+ }
+}
+
+func TestNormZero(t *testing.T) {
+ for _, a := range []Matrix{
+ &Dense{},
+ &SymDense{},
+ &SymDense{mat: blas64.Symmetric{Uplo: blas.Upper}},
+ &TriDense{},
+ &TriDense{mat: blas64.Triangular{Uplo: blas.Upper, Diag: blas.NonUnit}},
+ &VecDense{},
+ } {
+ for _, norm := range []float64{1, 2, math.Inf(1)} {
+ panicked, message := panics(func() { Norm(a, norm) })
+ if !panicked {
+ t.Errorf("expected panic for Norm(&%T{}, %v)", a, norm)
+ }
+ if message != ErrShape.Error() {
+ t.Errorf("unexpected panic string for Norm(&%T{}, %v): got:%s want:%s",
+ a, norm, message, ErrShape.Error())
+ }
+ }
+ }
+}
+
+func TestSum(t *testing.T) {
+ f := func(a Matrix) interface{} {
+ return Sum(a)
+ }
+ denseComparison := func(a *Dense) interface{} {
+ return Sum(a)
+ }
+ testOneInputFunc(t, "Sum", f, denseComparison, sameAnswerFloatApproxTol(1e-12), isAnyType, isAnySize)
+}
+
+func TestTrace(t *testing.T) {
+ for _, test := range []struct {
+ a *Dense
+ trace float64
+ }{
+ {
+ a: NewDense(3, 3, []float64{1, 2, 3, 4, 5, 6, 7, 8, 9}),
+ trace: 15,
+ },
+ } {
+ trace := Trace(test.a)
+ if trace != test.trace {
+ t.Errorf("Trace mismatch. Want %v, got %v", test.trace, trace)
+ }
+ }
+ f := func(a Matrix) interface{} {
+ return Trace(a)
+ }
+ denseComparison := func(a *Dense) interface{} {
+ return Trace(a)
+ }
+ testOneInputFunc(t, "Trace", f, denseComparison, sameAnswerFloat, isAnyType, isSquare)
+}
+
+func TestDoer(t *testing.T) {
+ type MatrixDoer interface {
+ Matrix
+ NonZeroDoer
+ RowNonZeroDoer
+ ColNonZeroDoer
+ }
+ ones := func(n int) []float64 {
+ data := make([]float64, n)
+ for i := range data {
+ data[i] = 1
+ }
+ return data
+ }
+ for i, m := range []MatrixDoer{
+ NewTriDense(3, Lower, ones(3*3)),
+ NewTriDense(3, Upper, ones(3*3)),
+ NewBandDense(6, 6, 1, 1, ones(3*6)),
+ NewBandDense(6, 10, 1, 1, ones(3*6)),
+ NewBandDense(10, 6, 1, 1, ones(7*3)),
+ NewSymBandDense(3, 0, ones(3)),
+ NewSymBandDense(3, 1, ones(3*(1+1))),
+ NewSymBandDense(6, 1, ones(6*(1+1))),
+ NewSymBandDense(6, 2, ones(6*(2+1))),
+ } {
+ r, c := m.Dims()
+
+ want := Sum(m)
+
+ // got and fn sum the accessed elements in
+ // the Doer that is being operated on.
+ // fn also tests that the accessed elements
+ // are within the writable areas of the
+ // matrix to check that only valid elements
+ // are operated on.
+ var got float64
+ fn := func(i, j int, v float64) {
+ got += v
+ switch m := m.(type) {
+ case MutableTriangular:
+ m.SetTri(i, j, v)
+ case MutableBanded:
+ m.SetBand(i, j, v)
+ case MutableSymBanded:
+ m.SetSymBand(i, j, v)
+ default:
+ panic("bad test: need mutable type")
+ }
+ }
+
+ panicked, message := panics(func() { m.DoNonZero(fn) })
+ if panicked {
+ t.Errorf("unexpected panic for Doer test %d: %q", i, message)
+ continue
+ }
+ if got != want {
+ t.Errorf("unexpected Doer sum: got:%f want:%f", got, want)
+ }
+
+ // Reset got for testing with DoRowNonZero.
+ got = 0
+ panicked, message = panics(func() {
+ for i := 0; i < r; i++ {
+ m.DoRowNonZero(i, fn)
+ }
+ })
+ if panicked {
+ t.Errorf("unexpected panic for RowDoer test %d: %q", i, message)
+ continue
+ }
+ if got != want {
+ t.Errorf("unexpected RowDoer sum: got:%f want:%f", got, want)
+ }
+
+ // Reset got for testing with DoColNonZero.
+ got = 0
+ panicked, message = panics(func() {
+ for j := 0; j < c; j++ {
+ m.DoColNonZero(j, fn)
+ }
+ })
+ if panicked {
+ t.Errorf("unexpected panic for ColDoer test %d: %q", i, message)
+ continue
+ }
+ if got != want {
+ t.Errorf("unexpected ColDoer sum: got:%f want:%f", got, want)
+ }
+ }
+}