Merge 00b6a4c into 20644db

AUTHORS

@@ -19,6 +19,7 @@ Dan Kortschak <dan.kortschak@adelaide.edu.au> <dan@kortschak.io>
 Daniel Fireman <danielfireman@gmail.com>
 David Samborski <bloggingarrow@gmail.com>
 Davor Kapsa <davor.kapsa@gmail.com>
+DeepMind Technologies
 Egon Elbre <egonelbre@gmail.com>
 Ekaterina Efimova <katerina.efimova@gmail.com>
 Ethan Burns <burns.ethan@gmail.com>

mat/diagonal.go

@@ -0,0 +1,119 @@
+// Copyright ©2018 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 (
+	"gonum.org/v1/gonum/blas"
+	"gonum.org/v1/gonum/blas/blas64"
+)
+
+var (
+	diagDense *DiagDense
+	_         Matrix          = diagDense
+	_         Diagonal        = diagDense
+	_         MutableDiagonal = diagDense
+	_         Triangular      = diagDense
+	_         Symmetric       = diagDense
+	_         Banded          = diagDense
+	_         RawBander       = diagDense
+	_         RawSymBander    = diagDense
+)
+
+// Diagonal represents a diagonal matrix, that is a square matrix that only
+// has non-zero terms on the diagonal.
+type Diagonal interface {
+	Matrix
+	// Diag returns the number of rows/columns in the matrix
+	Diag() int
+}
+
+// MutableDiagonal is a Diagonal matrix whose elements can be set.
+type MutableDiagonal interface {
+	Diagonal
+	SetDiag(i int, v float64)
+}
+
+// DiagDense represents a diagonal matrix in dense storage format.
+type DiagDense struct {
+	data []float64
+}
+
+// NewDiagonal creates a new Diagonal matrix with n rows and n columns.
+// The length of data must be n or data must be nil, otherwise NewDiagonal
+// will panic.
+func NewDiagonal(n int, data []float64) *DiagDense {
+	if n < 0 {
+		panic("mat: negative dimension")
+	}
+	if data == nil {
+		data = make([]float64, n)
+	}
+	if len(data) != n {
+		panic(ErrShape)
+	}
+	return &DiagDense{
+		data: data,
+	}
+}
+
+// Diag returns the dimension of the receiver.
+func (d *DiagDense) Diag() int {
+	return len(d.data)
+}
+
+// Dims returns the dimensions of the matrix.
+func (d *DiagDense) Dims() (r, c int) {
+	return len(d.data), len(d.data)
+}
+
+// T returns the transpose of the matrix.
+func (d *DiagDense) T() Matrix {
+	return d
+}
+
+// TTri returns the transpose of the matrix. Note that Diagonal matrices are
+// Upper by default
+func (d *DiagDense) TTri() Triangular {
+	return TransposeTri{d}
+}
+
+func (d *DiagDense) TBand() Banded {
+	return TransposeBand{d}
+}
+
+func (d *DiagDense) Bandwidth() (kl, ku int) {
+	return 0, 0
+}
+
+// Symmetric implements the Symmetric interface.
+func (d *DiagDense) Symmetric() int {
+	return len(d.data)
+}
+
+// Triangle implements the Triangular interface.
+func (d *DiagDense) Triangle() (int, TriKind) {
+	return len(d.data), Upper
+}
+
+func (d *DiagDense) RawBand() blas64.Band {
+	return blas64.Band{
+		Rows:   len(d.data),
+		Cols:   len(d.data),
+		KL:     0,
+		KU:     0,
+		Stride: 1,
+		Data:   d.data,
+	}
+}
+
+func (d *DiagDense) RawSymBand() blas64.SymmetricBand {
+	return blas64.SymmetricBand{
+		N:      len(d.data),
+		K:      0,
+		Stride: 1,
+		Uplo:   blas.Upper,
+		Data:   d.data,
+	}
+}

mat/diagonal_test.go

@@ -0,0 +1,91 @@
+// Copyright ©2018 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 (
+	"reflect"
+	"testing"
+)
+
+func TestNewDiagonal(t *testing.T) {
+	for i, test := range []struct {
+		data  []float64
+		n     int
+		mat   *DiagDense
+		dense *Dense
+	}{
+		{
+			data: []float64{1, 2, 3, 4, 5, 6},
+			n:    6,
+			mat: &DiagDense{
+				data: []float64{1, 2, 3, 4, 5, 6},
+			},
+			dense: NewDense(6, 6, []float64{
+				1, 0, 0, 0, 0, 0,
+				0, 2, 0, 0, 0, 0,
+				0, 0, 3, 0, 0, 0,
+				0, 0, 0, 4, 0, 0,
+				0, 0, 0, 0, 5, 0,
+				0, 0, 0, 0, 0, 6,
+			}),
+		},
+	} {
+		band := NewDiagonal(test.n, test.data)
+		rows, cols := band.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(band, test.mat) {
+			t.Errorf("unexpected value via reflect for test %d: got: %v want: %v", i, band, test.mat)
+		}
+		if !Equal(band, test.mat) {
+			t.Errorf("unexpected value via mat.Equal for test %d: got: %v want: %v", i, band, test.mat)
+		}
+		if !Equal(band, test.dense) {
+			t.Errorf("unexpected value via mat.Equal(band, dense) for test %d:\ngot:\n% v\nwant:\n% v", i, Formatted(band), Formatted(test.dense))
+		}
+	}
+}
+
+func TestDiagonalAtSet(t *testing.T) {
+	for _, n := range []int{1, 3, 8} {
+		for _, nilstart := range []bool{true, false} {
+			var diag *DiagDense
+			if nilstart {
+				diag = NewDiagonal(n, nil)
+			} else {
+				data := make([]float64, n)
+				diag = NewDiagonal(n, data)
+				// Test the data is used.
+				for i := range data {
+					data[i] = -float64(i) - 1
+					v := diag.At(i, i)
+					if v != data[i] {
+						t.Errorf("Diag shadow mismatch. Got %v, want %v", v, data[i])
+					}
+				}
+			}
+			for i := 0; i < n; i++ {
+				for j := 0; j < n; j++ {
+					if i != j {
+						if diag.At(i, j) != 0 {
+							t.Errorf("Diag returned non-zero off diagonal element at %d, %d", i, j)
+						}
+					}
+					v := float64(i) + 1
+					diag.SetDiag(i, v)
+					v2 := diag.At(i, i)
+					if v2 != v {
+						t.Errorf("Diag at/set mismatch. Got %v, want %v", v, v2)
+					}
+				}
+			}
+		}
+	}
+}

mat/errors.go

@@ -130,6 +130,7 @@ var (
 	ErrTriangle            = Error{"matrix: triangular storage mismatch"}
 	ErrTriangleSet         = Error{"matrix: triangular set out of bounds"}
 	ErrBandSet             = Error{"matrix: band set out of bounds"}
+	ErrDiagSet             = Error{"matrix: diagonal set out of bounds"}
 	ErrSliceLengthMismatch = Error{"matrix: input slice length mismatch"}
 	ErrNotPSD              = Error{"matrix: input not positive symmetric definite"}
 	ErrFailedEigen         = Error{"matrix: eigendecomposition not successful"}

mat/index_bound_checks.go

@@ -231,3 +231,34 @@ func (s *SymBandDense) set(i, j int, v float64) {
 	}
 	s.mat.Data[i*s.mat.Stride+pj] = v
 }
+
+// At returns the element at row i, column j.
+func (d *DiagDense) At(i, j int) float64 {
+	return d.at(i, j)
+}
+
+func (d *DiagDense) at(i, j int) float64 {
+	if uint(i) >= uint(len(d.data)) {
+		panic(ErrRowAccess)
+	}
+	if uint(j) >= uint(len(d.data)) {
+		panic(ErrColAccess)
+	}
+	if i != j {
+		return 0
+	}
+	return d.data[i]
+}
+
+// SetDiag sets the element at row i, column i to the value v.
+// It panics if the location is outside the appropriate region of the matrix.
+func (d *DiagDense) SetDiag(i int, v float64) {
+	d.setDiag(i, v)
+}
+
+func (d *DiagDense) setDiag(i int, v float64) {
+	if uint(i) >= uint(len(d.data)) {
+		panic(ErrRowAccess)
+	}
+	d.data[i] = v
+}

mat/index_no_bound_checks.go

@@ -235,3 +235,34 @@ func (s *SymBandDense) set(i, j int, v float64) {
 	}
 	s.mat.Data[i*s.mat.Stride+pj] = v
 }
+
+// At returns the element at row i, column j.
+func (d *DiagDense) At(i, j int) float64 {
+	if uint(i) >= uint(len(d.data)) {
+		panic(ErrRowAccess)
+	}
+	if uint(j) >= uint(len(d.data)) {
+		panic(ErrColAccess)
+	}
+	return d.at(i, j)
+}
+
+func (d *DiagDense) at(i, j int) float64 {
+	if i != j {
+		return 0
+	}
+	return d.data[i]
+}
+
+// SetDiag sets the element at row i, column i to the value v.
+// It panics if the location is outside the appropriate region of the matrix.
+func (d *DiagDense) SetDiag(i int, v float64) {
+	if uint(i) >= uint(len(d.data)) {
+		panic(ErrRowAccess)
+	}
+	d.setDiag(i, v)
+}
+
+func (d *DiagDense) setDiag(i int, v float64) {
+	d.data[i] = v
+}

mat/list_test.go

@@ -234,7 +234,7 @@ func legalDims(a Matrix, m, n int) bool {
 			return false
 		}
 		return true
-	case *SymDense, *TriDense, *basicSymmetric, *basicTriangular:
+	case *SymDense, *TriDense, *basicSymmetric, *basicTriangular, *DiagDense:
 		if m < 0 || n < 0 || m != n {
 			return false
 		}
@@ -280,6 +280,14 @@ func returnAs(a, t Matrix) Matrix {
 		case *basicTriangular:
 			return asBasicTriangular(mat)
 		}
+	case *DiagDense:
+		switch t.(type) {
+		default:
+			panic("bad type")
+		case *DiagDense:
+			return mat
+			// TODO(btracey): basic diagonaler
+		}
 	}
 }
 
@@ -393,6 +401,15 @@ func makeRandOf(a Matrix, m, n int) Matrix {
 			}
 		}
 		rMatrix = returnAs(mat, t)
+	case *DiagDense:
+		if m != n {
+			panic("bad size")
+		}
+		diag := NewDiagonal(m, nil)
+		for i := 0; i < n; i++ {
+			diag.SetDiag(i, rand.NormFloat64())
+		}
+		rMatrix = returnAs(diag, t)
 	}
 	if mr, mc := rMatrix.Dims(); mr != m || mc != n {
 		panic(fmt.Sprintf("makeRandOf for %T returns wrong size: %d×%d != %d×%d", a, m, n, mr, mc))
@@ -449,6 +466,12 @@ func makeCopyOf(a Matrix) Matrix {
 		}
 		copy(m.m, t.m)
 		return m
+	case *DiagDense:
+		d := &DiagDense{
+			data: make([]float64, len(t.data)),
+		}
+		copy(d.data, t.data)
+		return d
 	}
 }
 
@@ -582,6 +605,7 @@ var testMatrices = []Matrix{
 	NewTriDense(3, true, nil),
 	NewTriDense(3, false, nil),
 	NewVecDense(0, nil),
+	&DiagDense{},
 	&basicVector{},
 	&VecDense{mat: blas64.Vector{Inc: 10}},
 	&basicMatrix{},
@@ -596,6 +620,8 @@ var testMatrices = []Matrix{
 	TransposeTri{NewTriDense(3, false, nil)},
 	Transpose{NewVecDense(0, nil)},
 	Transpose{&VecDense{mat: blas64.Vector{Inc: 10}}},
+	Transpose{&DiagDense{}},
+	TransposeTri{&DiagDense{}},
 	Transpose{&basicMatrix{}},
 	Transpose{&basicSymmetric{}},
 	Transpose{&basicTriangular{cap: 3, mat: blas64.Triangular{N: 3, Stride: 3, Uplo: blas.Upper}}},

mat/symband.go

@@ -91,13 +91,6 @@ func NewSymBandDense(n, k int, data []float64) *SymBandDense {
 	}
 }
 
-// NewDiagonal is a convenience function that returns a diagonal matrix represented by a
-// SymBandDense. The length of data must be n or data must be nil, otherwise NewDiagonal
-// will panic.
-func NewDiagonal(n int, data []float64) *SymBandDense {
-	return NewSymBandDense(n, 0, data)
-}
-
 // Dims returns the number of rows and columns in the matrix.
 func (s *SymBandDense) Dims() (r, c int) {
 	return s.mat.N, s.mat.N