mat/symmetric_test.go

// Copyright ©2015 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"
	"os"
	"reflect"
	"testing"

	"golang.org/x/exp/rand"

	"gonum.org/v1/gonum/blas"
	"gonum.org/v1/gonum/blas/blas64"
	"gonum.org/v1/gonum/floats/scalar"
)

func TestNewSymmetric(t *testing.T) {
	t.Parallel()
	for i, test := range []struct {
		data []float64
		n    int
		mat  *SymDense
	}{
		{
			data: []float64{
				1, 2, 3,
				4, 5, 6,
				7, 8, 9,
			},
			n: 3,
			mat: &SymDense{
				mat: blas64.Symmetric{
					N:      3,
					Stride: 3,
					Uplo:   blas.Upper,
					Data:   []float64{1, 2, 3, 4, 5, 6, 7, 8, 9},
				},
				cap: 3,
			},
		},
	} {
		sym := NewSymDense(test.n, test.data)
		rows, cols := sym.Dims()

		if rows != test.n {
			t.Errorf("unexpected number of rows for test %d: got: %d want: %d", i, rows, test.n)
		}
		if cols != test.n {
			t.Errorf("unexpected number of cols for test %d: got: %d want: %d", i, cols, test.n)
		}
		if !reflect.DeepEqual(sym, test.mat) {
			t.Errorf("unexpected data slice for test %d: got: %v want: %v", i, sym, test.mat)
		}

		m := NewDense(test.n, test.n, test.data)
		if !reflect.DeepEqual(sym.mat.Data, m.mat.Data) {
			t.Errorf("unexpected data slice mismatch for test %d: got: %v want: %v", i, sym.mat.Data, m.mat.Data)
		}
	}

	panicked, message := panics(func() { NewSymDense(3, []float64{1, 2}) })
	if !panicked || message != ErrShape.Error() {
		t.Error("expected panic for invalid data slice length")
	}
}

func TestSymAtSet(t *testing.T) {
	t.Parallel()
	sym := &SymDense{
		mat: blas64.Symmetric{
			N:      3,
			Stride: 3,
			Uplo:   blas.Upper,
			Data:   []float64{1, 2, 3, 4, 5, 6, 7, 8, 9},
		},
		cap: 3,
	}
	rows, cols := sym.Dims()

	// Check At out of bounds
	for _, row := range []int{-1, rows, rows + 1} {
		panicked, message := panics(func() { sym.At(row, 0) })
		if !panicked || message != ErrRowAccess.Error() {
			t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
		}
	}
	for _, col := range []int{-1, cols, cols + 1} {
		panicked, message := panics(func() { sym.At(0, col) })
		if !panicked || message != ErrColAccess.Error() {
			t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
		}
	}

	// Check Set out of bounds
	for _, row := range []int{-1, rows, rows + 1} {
		panicked, message := panics(func() { sym.SetSym(row, 0, 1.2) })
		if !panicked || message != ErrRowAccess.Error() {
			t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
		}
	}
	for _, col := range []int{-1, cols, cols + 1} {
		panicked, message := panics(func() { sym.SetSym(0, col, 1.2) })
		if !panicked || message != ErrColAccess.Error() {
			t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
		}
	}

	for _, st := range []struct {
		row, col  int
		orig, new float64
	}{
		{row: 1, col: 2, orig: 6, new: 15},
		{row: 2, col: 1, orig: 15, new: 12},
	} {
		if e := sym.At(st.row, st.col); e != st.orig {
			t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.row, st.col, e, st.orig)
		}
		if e := sym.At(st.col, st.row); e != st.orig {
			t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.col, st.row, e, st.orig)
		}
		sym.SetSym(st.row, st.col, st.new)
		if e := sym.At(st.row, st.col); e != st.new {
			t.Errorf("unexpected value for At(%d, %d) after SetSym(%[1]d, %[2]d, %[4]v): got: %[3]v want: %v", st.row, st.col, e, st.new)
		}
		if e := sym.At(st.col, st.row); e != st.new {
			t.Errorf("unexpected value for At(%d, %d) after SetSym(%[2]d, %[1]d, %[4]v): got: %[3]v want: %v", st.col, st.row, e, st.new)
		}
	}
}

func TestSymDenseZero(t *testing.T) {
	t.Parallel()
	// Elements that equal 1 should be set to zero, elements that equal -1
	// should remain unchanged.
	for _, test := range []*SymDense{
		{
			mat: blas64.Symmetric{
				Uplo:   blas.Upper,
				N:      4,
				Stride: 5,
				Data: []float64{
					1, 1, 1, 1, -1,
					-1, 1, 1, 1, -1,
					-1, -1, 1, 1, -1,
					-1, -1, -1, 1, -1,
				},
			},
		},
	} {
		dataCopy := make([]float64, len(test.mat.Data))
		copy(dataCopy, test.mat.Data)
		test.Zero()
		for i, v := range test.mat.Data {
			if dataCopy[i] != -1 && v != 0 {
				t.Errorf("Matrix not zeroed in bounds")
			}
			if dataCopy[i] == -1 && v != -1 {
				t.Errorf("Matrix zeroed out of bounds")
			}
		}
	}
}

func TestSymDiagView(t *testing.T) {
	t.Parallel()
	for cas, test := range []*SymDense{
		NewSymDense(1, []float64{1}),
		NewSymDense(2, []float64{1, 2, 2, 3}),
		NewSymDense(3, []float64{1, 2, 3, 2, 4, 5, 3, 5, 6}),
	} {
		testDiagView(t, cas, test)
	}
}

func TestSymAdd(t *testing.T) {
	t.Parallel()
	rnd := rand.New(rand.NewSource(1))
	for _, test := range []struct {
		n int
	}{
		{n: 1},
		{n: 2},
		{n: 3},
		{n: 4},
		{n: 5},
		{n: 10},
	} {
		n := test.n
		a := NewSymDense(n, nil)
		for i := range a.mat.Data {
			a.mat.Data[i] = rnd.Float64()
		}
		b := NewSymDense(n, nil)
		for i := range a.mat.Data {
			b.mat.Data[i] = rnd.Float64()
		}
		var m Dense
		m.Add(a, b)

		// Check with new receiver
		var s SymDense
		s.AddSym(a, b)
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				want := m.At(i, j)
				if got := s.At(i, j); got != want {
					t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
				}
			}
		}

		// Check with equal receiver
		s.CopySym(a)
		s.AddSym(&s, b)
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				want := m.At(i, j)
				if got := s.At(i, j); got != want {
					t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
				}
			}
		}
	}

	method := func(receiver, a, b Matrix) {
		type addSymer interface {
			AddSym(a, b Symmetric)
		}
		rd := receiver.(addSymer)
		rd.AddSym(a.(Symmetric), b.(Symmetric))
	}
	denseComparison := func(receiver, a, b *Dense) {
		receiver.Add(a, b)
	}
	testTwoInput(t, "AddSym", &SymDense{}, method, denseComparison, legalTypesSym, legalSizeSameSquare, 1e-14)
}

func TestCopy(t *testing.T) {
	t.Parallel()
	rnd := rand.New(rand.NewSource(1))
	for _, test := range []struct {
		n int
	}{
		{n: 1},
		{n: 2},
		{n: 3},
		{n: 4},
		{n: 5},
		{n: 10},
	} {
		n := test.n
		a := NewSymDense(n, nil)
		for i := range a.mat.Data {
			a.mat.Data[i] = rnd.Float64()
		}
		s := NewSymDense(n, nil)
		s.CopySym(a)
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				want := a.At(i, j)
				if got := s.At(i, j); got != want {
					t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
				}
			}
		}
	}
}

// TODO(kortschak) Roll this into testOneInput when it exists.
// https://github.com/gonum/matrix/issues/171
func TestSymCopyPanic(t *testing.T) {
	t.Parallel()
	var (
		a SymDense
		n int
	)
	m := NewSymDense(1, nil)
	panicked, message := panics(func() { n = m.CopySym(&a) })
	if panicked {
		t.Errorf("unexpected panic: %v", message)
	}
	if n != 0 {
		t.Errorf("unexpected n: got: %d want: 0", n)
	}
}

func TestSymRankOne(t *testing.T) {
	t.Parallel()
	rnd := rand.New(rand.NewSource(1))
	const tol = 1e-15

	for _, test := range []struct {
		n int
	}{
		{n: 1},
		{n: 2},
		{n: 3},
		{n: 4},
		{n: 5},
		{n: 10},
	} {
		n := test.n
		alpha := 2.0
		a := NewSymDense(n, nil)
		for i := range a.mat.Data {
			a.mat.Data[i] = rnd.Float64()
		}
		x := make([]float64, n)
		for i := range x {
			x[i] = rnd.Float64()
		}

		xMat := NewDense(n, 1, x)
		var m Dense
		m.Mul(xMat, xMat.T())
		m.Scale(alpha, &m)
		m.Add(&m, a)

		// Check with new receiver
		s := NewSymDense(n, nil)
		s.SymRankOne(a, alpha, NewVecDense(len(x), x))
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				want := m.At(i, j)
				if got := s.At(i, j); !scalar.EqualWithinAbsOrRel(got, want, tol, tol) {
					t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
				}
			}
		}

		// Check with reused receiver
		copy(s.mat.Data, a.mat.Data)
		s.SymRankOne(s, alpha, NewVecDense(len(x), x))
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				want := m.At(i, j)
				if got := s.At(i, j); !scalar.EqualWithinAbsOrRel(got, want, tol, tol) {
					t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want)
				}
			}
		}
	}

	alpha := 3.0
	method := func(receiver, a, b Matrix) {
		type SymRankOner interface {
			SymRankOne(a Symmetric, alpha float64, x Vector)
		}
		rd := receiver.(SymRankOner)
		rd.SymRankOne(a.(Symmetric), alpha, b.(Vector))
	}
	denseComparison := func(receiver, a, b *Dense) {
		var tmp Dense
		tmp.Mul(b, b.T())
		tmp.Scale(alpha, &tmp)
		receiver.Add(a, &tmp)
	}
	legalTypes := func(a, b Matrix) bool {
		_, ok := a.(Symmetric)
		if !ok {
			return false
		}
		_, ok = b.(Vector)
		return ok
	}
	legalSize := func(ar, ac, br, bc int) bool {
		if ar != ac {
			return false
		}
		return br == ar
	}
	testTwoInput(t, "SymRankOne", &SymDense{}, method, denseComparison, legalTypes, legalSize, 1e-14)
}

func TestIssue250SymRankOne(t *testing.T) {
	t.Parallel()
	x := NewVecDense(5, []float64{1, 2, 3, 4, 5})
	var s1, s2 SymDense
	s1.SymRankOne(NewSymDense(5, nil), 1, x)
	s2.SymRankOne(NewSymDense(5, nil), 1, x)
	s2.SymRankOne(NewSymDense(5, nil), 1, x)
	if !Equal(&s1, &s2) {
		t.Error("unexpected result from repeat")
	}
}

func TestRankTwo(t *testing.T) {
	t.Parallel()
	rnd := rand.New(rand.NewSource(1))
	for _, test := range []struct {
		n int
	}{
		{n: 1},
		{n: 2},
		{n: 3},
		{n: 4},
		{n: 5},
		{n: 10},
	} {
		n := test.n
		alpha := 2.0
		a := NewSymDense(n, nil)
		for i := range a.mat.Data {
			a.mat.Data[i] = rnd.Float64()
		}
		x := make([]float64, n)
		y := make([]float64, n)
		for i := range x {
			x[i] = rnd.Float64()
			y[i] = rnd.Float64()
		}

		xMat := NewDense(n, 1, x)
		yMat := NewDense(n, 1, y)
		var m Dense
		m.Mul(xMat, yMat.T())
		var tmp Dense
		tmp.Mul(yMat, xMat.T())
		m.Add(&m, &tmp)
		m.Scale(alpha, &m)
		m.Add(&m, a)

		// Check with new receiver
		s := NewSymDense(n, nil)
		s.RankTwo(a, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y))
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				if !scalar.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
					t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j))
				}
			}
		}

		// Check with reused receiver
		copy(s.mat.Data, a.mat.Data)
		s.RankTwo(s, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y))
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				if !scalar.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) {
					t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j))
				}
			}
		}
	}
}

func TestSymRankK(t *testing.T) {
	t.Parallel()
	alpha := 3.0
	method := func(receiver, a, b Matrix) {
		type SymRankKer interface {
			SymRankK(a Symmetric, alpha float64, x Matrix)
		}
		rd := receiver.(SymRankKer)
		rd.SymRankK(a.(Symmetric), alpha, b)
	}
	denseComparison := func(receiver, a, b *Dense) {
		var tmp Dense
		tmp.Mul(b, b.T())
		tmp.Scale(alpha, &tmp)
		receiver.Add(a, &tmp)
	}
	legalTypes := func(a, b Matrix) bool {
		_, ok := a.(Symmetric)
		return ok
	}
	legalSize := func(ar, ac, br, bc int) bool {
		if ar != ac {
			return false
		}
		return br == ar
	}
	testTwoInput(t, "SymRankK", &SymDense{}, method, denseComparison, legalTypes, legalSize, 1e-14)
}

func TestSymOuterK(t *testing.T) {
	t.Parallel()
	for _, f := range []float64{0.5, 1, 3} {
		method := func(receiver, x Matrix) {
			type SymOuterKer interface {
				SymOuterK(alpha float64, x Matrix)
			}
			rd := receiver.(SymOuterKer)
			rd.SymOuterK(f, x)
		}
		denseComparison := func(receiver, x *Dense) {
			receiver.Mul(x, x.T())
			receiver.Scale(f, receiver)
		}
		testOneInput(t, "SymOuterK", &SymDense{}, method, denseComparison, isAnyType, isAnySize, 1e-14)
	}
}

func TestIssue250SymOuterK(t *testing.T) {
	t.Parallel()
	x := NewVecDense(5, []float64{1, 2, 3, 4, 5})
	var s1, s2 SymDense
	s1.SymOuterK(1, x)
	s2.SymOuterK(1, x)
	s2.SymOuterK(1, x)
	if !Equal(&s1, &s2) {
		t.Error("unexpected result from repeat")
	}
}

func TestScaleSym(t *testing.T) {
	t.Parallel()
	for _, f := range []float64{0.5, 1, 3} {
		method := func(receiver, a Matrix) {
			type ScaleSymer interface {
				ScaleSym(f float64, a Symmetric)
			}
			rd := receiver.(ScaleSymer)
			rd.ScaleSym(f, a.(Symmetric))
		}
		denseComparison := func(receiver, a *Dense) {
			receiver.Scale(f, a)
		}
		testOneInput(t, "ScaleSym", &SymDense{}, method, denseComparison, legalTypeSym, isSquare, 1e-14)
	}
}

func TestSubsetSym(t *testing.T) {
	t.Parallel()
	for _, test := range []struct {
		a    *SymDense
		dims []int
		ans  *SymDense
	}{
		{
			a: NewSymDense(3, []float64{
				1, 2, 3,
				0, 4, 5,
				0, 0, 6,
			}),
			dims: []int{0, 2},
			ans: NewSymDense(2, []float64{
				1, 3,
				0, 6,
			}),
		},
		{
			a: NewSymDense(3, []float64{
				1, 2, 3,
				0, 4, 5,
				0, 0, 6,
			}),
			dims: []int{2, 0},
			ans: NewSymDense(2, []float64{
				6, 3,
				0, 1,
			}),
		},
		{
			a: NewSymDense(3, []float64{
				1, 2, 3,
				0, 4, 5,
				0, 0, 6,
			}),
			dims: []int{1, 1, 1},
			ans: NewSymDense(3, []float64{
				4, 4, 4,
				0, 4, 4,
				0, 0, 4,
			}),
		},
	} {
		var s SymDense
		s.SubsetSym(test.a, test.dims)
		if !Equal(&s, test.ans) {
			t.Errorf("SubsetSym mismatch dims %v\nGot:\n% v\nWant:\n% v\n", test.dims, s, test.ans)
		}
	}

	dims := []int{0, 2}
	maxDim := dims[0]
	for _, v := range dims {
		if maxDim < v {
			maxDim = v
		}
	}
	method := func(receiver, a Matrix) {
		type SubsetSymer interface {
			SubsetSym(a Symmetric, set []int)
		}
		rd := receiver.(SubsetSymer)
		rd.SubsetSym(a.(Symmetric), dims)
	}
	denseComparison := func(receiver, a *Dense) {
		*receiver = *NewDense(len(dims), len(dims), nil)
		sz := len(dims)
		for i := 0; i < sz; i++ {
			for j := 0; j < sz; j++ {
				receiver.Set(i, j, a.At(dims[i], dims[j]))
			}
		}
	}
	legalSize := func(ar, ac int) bool {
		return ar == ac && ar > maxDim
	}

	testOneInput(t, "SubsetSym", &SymDense{}, method, denseComparison, legalTypeSym, legalSize, 0)
}

func TestViewGrowSquare(t *testing.T) {
	t.Parallel()
	// n is the size of the original SymDense.
	// The first view uses start1, span1. The second view uses start2, span2 on
	// the first view.
	for _, test := range []struct {
		n, start1, span1, start2, span2 int
	}{
		{10, 0, 10, 0, 10},
		{10, 0, 8, 0, 8},
		{10, 2, 8, 0, 6},
		{10, 2, 7, 4, 2},
		{10, 2, 6, 0, 5},
	} {
		n := test.n
		s := NewSymDense(n, nil)
		for i := 0; i < n; i++ {
			for j := i; j < n; j++ {
				s.SetSym(i, j, float64((i+1)*n+j+1))
			}
		}

		// Take a subset and check the view matches.
		start1 := test.start1
		span1 := test.span1
		v := s.sliceSym(start1, start1+span1)
		for i := 0; i < span1; i++ {
			for j := i; j < span1; j++ {
				if v.At(i, j) != s.At(start1+i, start1+j) {
					t.Errorf("View mismatch")
				}
			}
		}

		start2 := test.start2
		span2 := test.span2
		v2 := v.SliceSym(start2, start2+span2).(*SymDense)

		for i := 0; i < span2; i++ {
			for j := i; j < span2; j++ {
				if v2.At(i, j) != s.At(start1+start2+i, start1+start2+j) {
					t.Errorf("Second view mismatch")
				}
			}
		}

		// Check that a write to the view is reflected in the original.
		v2.SetSym(0, 0, 1.2)
		if s.At(start1+start2, start1+start2) != 1.2 {
			t.Errorf("Write to view not reflected in original")
		}

		// Grow the matrix back to the original view
		gn := n - start1 - start2
		g := v2.GrowSym(gn - v2.Symmetric()).(*SymDense)
		g.SetSym(1, 1, 2.2)

		for i := 0; i < gn; i++ {
			for j := 0; j < gn; j++ {
				if g.At(i, j) != s.At(start1+start2+i, start1+start2+j) {
					t.Errorf("Grow mismatch")

					fmt.Printf("g=\n% v\n", Formatted(g))
					fmt.Printf("s=\n% v\n", Formatted(s))
					os.Exit(1)
				}
			}
		}

		// View g, then grow it and make sure all the elements were copied.
		gv := g.SliceSym(0, gn-1).(*SymDense)

		gg := gv.GrowSym(2)
		for i := 0; i < gn; i++ {
			for j := 0; j < gn; j++ {
				if g.At(i, j) != gg.At(i, j) {
					t.Errorf("Expand mismatch")
				}
			}
		}

		s.Reset()
		rg := s.GrowSym(n).(*SymDense)
		if rg.mat.Stride < n {
			t.Errorf("unexpected stride after GrowSym on empty matrix: got:%d want >= %d", rg.mat.Stride, n)
		}
	}
}

func TestPowPSD(t *testing.T) {
	t.Parallel()
	for cas, test := range []struct {
		a   *SymDense
		pow float64
		ans *SymDense
	}{
		// Comparison with Matlab.
		{
			a:   NewSymDense(2, []float64{10, 5, 5, 12}),
			pow: 0.5,
			ans: NewSymDense(2, []float64{3.065533767740645, 0.776210486171016, 0.776210486171016, 3.376017962209052}),
		},
		{
			a:   NewSymDense(2, []float64{11, -1, -1, 8}),
			pow: 0.5,
			ans: NewSymDense(2, []float64{3.312618742210524, -0.162963396980939, -0.162963396980939, 2.823728551267709}),
		},
		{
			a:   NewSymDense(2, []float64{10, 5, 5, 12}),
			pow: -0.5,
			ans: NewSymDense(2, []float64{0.346372134547712, -0.079637515547296, -0.079637515547296, 0.314517128328794}),
		},
		{
			a:   NewSymDense(3, []float64{15, -1, -3, -1, 8, 6, -3, 6, 14}),
			pow: 0.6,
			ans: NewSymDense(3, []float64{
				5.051214323034288, -0.163162161893975, -0.612153996497505,
				-0.163162161893976, 3.283474884617009, 1.432842761381493,
				-0.612153996497505, 1.432842761381494, 4.695873060862573,
			}),
		},
	} {
		var s SymDense
		err := s.PowPSD(test.a, test.pow)
		if err != nil {
			panic("bad test")
		}
		if !EqualApprox(&s, test.ans, 1e-10) {
			t.Errorf("Case %d, pow mismatch", cas)
			fmt.Println(Formatted(&s))
			fmt.Println(Formatted(test.ans))
		}
	}

	// Compare with Dense.Pow
	rnd := rand.New(rand.NewSource(1))
	for dim := 2; dim < 10; dim++ {
		for pow := 2; pow < 6; pow++ {
			a := NewDense(dim, dim, nil)
			for i := 0; i < dim; i++ {
				for j := 0; j < dim; j++ {
					a.Set(i, j, rnd.Float64())
				}
			}
			var mat SymDense
			mat.SymOuterK(1, a)

			var sym SymDense
			err := sym.PowPSD(&mat, float64(pow))
			if err != nil {
				t.Errorf("unexpected error: %v", err)
			}

			var dense Dense
			dense.Pow(&mat, pow)

			if !EqualApprox(&sym, &dense, 1e-10) {
				t.Errorf("Dim %d: pow mismatch", dim)
			}
		}
	}
}

func BenchmarkSymSum1000(b *testing.B) { symSumBench(b, 1000) }

var symSumForBench float64

func symSumBench(b *testing.B, size int) {
	src := rand.NewSource(1)
	a := randSymDense(size, src)
	b.ResetTimer()
	for i := 0; i < b.N; i++ {
		symSumForBench = Sum(a)
	}
}

func randSymDense(size int, src rand.Source) *SymDense {
	rnd := rand.New(src)
	backData := make([]float64, size*size)
	for i := 0; i < size; i++ {
		backData[i*size+i] = rnd.Float64()
		for j := i + 1; j < size; j++ {
			v := rnd.Float64()
			backData[i*size+j] = v
			backData[j*size+i] = v
		}
	}
	s := NewSymDense(size, backData)
	return s
}