vendor/gonum.org/v1/gonum/mat/list_test.go

   1 // Copyright ©2015 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package mat
   6
   7 import (
   8         "fmt"
   9         "math"
  10         "reflect"
  11         "testing"
  12
  13         "golang.org/x/exp/rand"
  14
  15         "gonum.org/v1/gonum/blas"
  16         "gonum.org/v1/gonum/blas/blas64"
  17         "gonum.org/v1/gonum/floats"
  18 )
  19
  20 // legalSizeSameRectangular returns whether the two matrices have the same rectangular shape.
  21 func legalSizeSameRectangular(ar, ac, br, bc int) bool {
  22         if ar != br {
  23                 return false
  24         }
  25         if ac != bc {
  26                 return false
  27         }
  28         return true
  29 }
  30
  31 // legalSizeSameSquare returns whether the two matrices have the same square shape.
  32 func legalSizeSameSquare(ar, ac, br, bc int) bool {
  33         if ar != br {
  34                 return false
  35         }
  36         if ac != bc {
  37                 return false
  38         }
  39         if ar != ac {
  40                 return false
  41         }
  42         return true
  43 }
  44
  45 // legalSizeSameHeight returns whether the two matrices have the same number of rows.
  46 func legalSizeSameHeight(ar, _, br, _ int) bool {
  47         return ar == br
  48 }
  49
  50 // legalSizeSameWidth returns whether the two matrices have the same number of columns.
  51 func legalSizeSameWidth(_, ac, _, bc int) bool {
  52         return ac == bc
  53 }
  54
  55 // legalSizeSolve returns whether the two matrices can be used in a linear solve.
  56 func legalSizeSolve(ar, ac, br, bc int) bool {
  57         return ar == br
  58 }
  59
  60 // legalSizeSameVec returns whether the two matrices are column vectors.
  61 func legalSizeVector(_, ac, _, bc int) bool {
  62         return ac == 1 && bc == 1
  63 }
  64
  65 // legalSizeSameVec returns whether the two matrices are column vectors of the
  66 // same dimension.
  67 func legalSizeSameVec(ar, ac, br, bc int) bool {
  68         return ac == 1 && bc == 1 && ar == br
  69 }
  70
  71 // isAnySize returns true for all matrix sizes.
  72 func isAnySize(ar, ac int) bool {
  73         return true
  74 }
  75
  76 // isAnySize2 returns true for all matrix sizes.
  77 func isAnySize2(ar, ac, br, bc int) bool {
  78         return true
  79 }
  80
  81 // isAnyColumnVector returns true for any column vector sizes.
  82 func isAnyColumnVector(ar, ac int) bool {
  83         return ac == 1
  84 }
  85
  86 // isSquare returns whether the input matrix is square.
  87 func isSquare(r, c int) bool {
  88         return r == c
  89 }
  90
  91 // sameAnswerFloat returns whether the two inputs are both NaN or are equal.
  92 func sameAnswerFloat(a, b interface{}) bool {
  93         if math.IsNaN(a.(float64)) {
  94                 return math.IsNaN(b.(float64))
  95         }
  96         return a.(float64) == b.(float64)
  97 }
  98
  99 // sameAnswerFloatApproxTol returns a function that determines whether its two
 100 // inputs are both NaN or within tol of each other.
 101 func sameAnswerFloatApproxTol(tol float64) func(a, b interface{}) bool {
 102         return func(a, b interface{}) bool {
 103                 if math.IsNaN(a.(float64)) {
 104                         return math.IsNaN(b.(float64))
 105                 }
 106                 return floats.EqualWithinAbsOrRel(a.(float64), b.(float64), tol, tol)
 107         }
 108 }
 109
 110 func sameAnswerF64SliceOfSlice(a, b interface{}) bool {
 111         for i, v := range a.([][]float64) {
 112                 if same := floats.Same(v, b.([][]float64)[i]); !same {
 113                         return false
 114                 }
 115         }
 116         return true
 117 }
 118
 119 // sameAnswerBool returns whether the two inputs have the same value.
 120 func sameAnswerBool(a, b interface{}) bool {
 121         return a.(bool) == b.(bool)
 122 }
 123
 124 // isAnyType returns true for all Matrix types.
 125 func isAnyType(Matrix) bool {
 126         return true
 127 }
 128
 129 // legalTypesAll returns true for all Matrix types.
 130 func legalTypesAll(a, b Matrix) bool {
 131         return true
 132 }
 133
 134 // legalTypeSym returns whether a is a Symmetric.
 135 func legalTypeSym(a Matrix) bool {
 136         _, ok := a.(Symmetric)
 137         return ok
 138 }
 139
 140 // legalTypeTri returns whether a is a Triangular.
 141 func legalTypeTri(a Matrix) bool {
 142         _, ok := a.(Triangular)
 143         return ok
 144 }
 145
 146 // legalTypeTriLower returns whether a is a Triangular with kind == Lower.
 147 func legalTypeTriLower(a Matrix) bool {
 148         t, ok := a.(Triangular)
 149         if !ok {
 150                 return false
 151         }
 152         _, kind := t.Triangle()
 153         return kind == Lower
 154 }
 155
 156 // legalTypeTriUpper returns whether a is a Triangular with kind == Upper.
 157 func legalTypeTriUpper(a Matrix) bool {
 158         t, ok := a.(Triangular)
 159         if !ok {
 160                 return false
 161         }
 162         _, kind := t.Triangle()
 163         return kind == Upper
 164 }
 165
 166 // legalTypesSym returns whether both input arguments are Symmetric.
 167 func legalTypesSym(a, b Matrix) bool {
 168         if _, ok := a.(Symmetric); !ok {
 169                 return false
 170         }
 171         if _, ok := b.(Symmetric); !ok {
 172                 return false
 173         }
 174         return true
 175 }
 176
 177 // legalTypeVector returns whether v is a Vector.
 178 func legalTypeVector(v Matrix) bool {
 179         _, ok := v.(Vector)
 180         return ok
 181 }
 182
 183 // legalTypeVec returns whether v is a *VecDense.
 184 func legalTypeVecDense(v Matrix) bool {
 185         _, ok := v.(*VecDense)
 186         return ok
 187 }
 188
 189 // legalTypesVectorVector returns whether both inputs are Vector
 190 func legalTypesVectorVector(a, b Matrix) bool {
 191         if _, ok := a.(Vector); !ok {
 192                 return false
 193         }
 194         if _, ok := b.(Vector); !ok {
 195                 return false
 196         }
 197         return true
 198 }
 199
 200 // legalTypesVecDenseVecDense returns whether both inputs are *VecDense.
 201 func legalTypesVecDenseVecDense(a, b Matrix) bool {
 202         if _, ok := a.(*VecDense); !ok {
 203                 return false
 204         }
 205         if _, ok := b.(*VecDense); !ok {
 206                 return false
 207         }
 208         return true
 209 }
 210
 211 // legalTypesMatrixVector returns whether the first input is an arbitrary Matrix
 212 // and the second input is a Vector.
 213 func legalTypesMatrixVector(a, b Matrix) bool {
 214         _, ok := b.(Vector)
 215         return ok
 216 }
 217
 218 // legalTypesMatrixVecDense returns whether the first input is an arbitrary Matrix
 219 // and the second input is a *VecDense.
 220 func legalTypesMatrixVecDense(a, b Matrix) bool {
 221         _, ok := b.(*VecDense)
 222         return ok
 223 }
 224
 225 // legalDims returns whether {m,n} is a valid dimension of the given matrix type.
 226 func legalDims(a Matrix, m, n int) bool {
 227         switch t := a.(type) {
 228         default:
 229                 panic("legal dims type not coded")
 230         case Untransposer:
 231                 return legalDims(t.Untranspose(), n, m)
 232         case *Dense, *basicMatrix:
 233                 if m < 0 || n < 0 {
 234                         return false
 235                 }
 236                 return true
 237         case *SymDense, *TriDense, *basicSymmetric, *basicTriangular:
 238                 if m < 0 || n < 0 || m != n {
 239                         return false
 240                 }
 241                 return true
 242         case *VecDense, *basicVector:
 243                 if m < 0 || n < 0 {
 244                         return false
 245                 }
 246                 return n == 1
 247         }
 248 }
 249
 250 // returnAs returns the matrix a with the type of t. Used for making a concrete
 251 // type and changing to the basic form.
 252 func returnAs(a, t Matrix) Matrix {
 253         switch mat := a.(type) {
 254         default:
 255                 panic("unknown type for a")
 256         case *Dense:
 257                 switch t.(type) {
 258                 default:
 259                         panic("bad type")
 260                 case *Dense:
 261                         return mat
 262                 case *basicMatrix:
 263                         return asBasicMatrix(mat)
 264                 }
 265         case *SymDense:
 266                 switch t.(type) {
 267                 default:
 268                         panic("bad type")
 269                 case *SymDense:
 270                         return mat
 271                 case *basicSymmetric:
 272                         return asBasicSymmetric(mat)
 273                 }
 274         case *TriDense:
 275                 switch t.(type) {
 276                 default:
 277                         panic("bad type")
 278                 case *TriDense:
 279                         return mat
 280                 case *basicTriangular:
 281                         return asBasicTriangular(mat)
 282                 }
 283         }
 284 }
 285
 286 // retranspose returns the matrix m inside an Untransposer of the type
 287 // of a.
 288 func retranspose(a, m Matrix) Matrix {
 289         switch a.(type) {
 290         case TransposeTri:
 291                 return TransposeTri{m.(Triangular)}
 292         case Transpose:
 293                 return Transpose{m}
 294         case Untransposer:
 295                 panic("unknown transposer type")
 296         default:
 297                 panic("a is not an untransposer")
 298         }
 299 }
 300
 301 // makeRandOf returns a new randomly filled m×n matrix of the underlying matrix type.
 302 func makeRandOf(a Matrix, m, n int) Matrix {
 303         var rMatrix Matrix
 304         switch t := a.(type) {
 305         default:
 306                 panic("unknown type for make rand of")
 307         case Untransposer:
 308                 rMatrix = retranspose(a, makeRandOf(t.Untranspose(), n, m))
 309         case *Dense, *basicMatrix:
 310                 mat := NewDense(m, n, nil)
 311                 for i := 0; i < m; i++ {
 312                         for j := 0; j < n; j++ {
 313                                 mat.Set(i, j, rand.NormFloat64())
 314                         }
 315                 }
 316                 rMatrix = returnAs(mat, t)
 317         case *VecDense:
 318                 if m == 0 && n == 0 {
 319                         return &VecDense{}
 320                 }
 321                 if n != 1 {
 322                         panic(fmt.Sprintf("bad vector size: m = %v, n = %v", m, n))
 323                 }
 324                 length := m
 325                 inc := 1
 326                 if t.mat.Inc != 0 {
 327                         inc = t.mat.Inc
 328                 }
 329                 mat := &VecDense{
 330                         mat: blas64.Vector{
 331                                 Inc:  inc,
 332                                 Data: make([]float64, inc*(length-1)+1),
 333                         },
 334                         n: length,
 335                 }
 336                 for i := 0; i < length; i++ {
 337                         mat.SetVec(i, rand.NormFloat64())
 338                 }
 339                 return mat
 340         case *basicVector:
 341                 if m == 0 && n == 0 {
 342                         return &basicVector{}
 343                 }
 344                 if n != 1 {
 345                         panic(fmt.Sprintf("bad vector size: m = %v, n = %v", m, n))
 346                 }
 347                 mat := &basicVector{
 348                         m: make([]float64, m),
 349                 }
 350                 for i := 0; i < m; i++ {
 351                         mat.m[i] = rand.NormFloat64()
 352                 }
 353                 return mat
 354         case *SymDense, *basicSymmetric:
 355                 if m != n {
 356                         panic("bad size")
 357                 }
 358                 mat := NewSymDense(n, nil)
 359                 for i := 0; i < m; i++ {
 360                         for j := i; j < n; j++ {
 361                                 mat.SetSym(i, j, rand.NormFloat64())
 362                         }
 363                 }
 364                 rMatrix = returnAs(mat, t)
 365         case *TriDense, *basicTriangular:
 366                 if m != n {
 367                         panic("bad size")
 368                 }
 369
 370                 // This is necessary because we are making
 371                 // a triangle from the zero value, which
 372                 // always returns upper as true.
 373                 var triKind TriKind
 374                 switch t := t.(type) {
 375                 case *TriDense:
 376                         triKind = t.triKind()
 377                 case *basicTriangular:
 378                         triKind = (*TriDense)(t).triKind()
 379                 }
 380
 381                 mat := NewTriDense(n, triKind, nil)
 382                 if triKind == Upper {
 383                         for i := 0; i < m; i++ {
 384                                 for j := i; j < n; j++ {
 385                                         mat.SetTri(i, j, rand.NormFloat64())
 386                                 }
 387                         }
 388                 } else {
 389                         for i := 0; i < m; i++ {
 390                                 for j := 0; j <= i; j++ {
 391                                         mat.SetTri(i, j, rand.NormFloat64())
 392                                 }
 393                         }
 394                 }
 395                 rMatrix = returnAs(mat, t)
 396         }
 397         if mr, mc := rMatrix.Dims(); mr != m || mc != n {
 398                 panic(fmt.Sprintf("makeRandOf for %T returns wrong size: %d×%d != %d×%d", a, m, n, mr, mc))
 399         }
 400         return rMatrix
 401 }
 402
 403 // makeCopyOf returns a copy of the matrix.
 404 func makeCopyOf(a Matrix) Matrix {
 405         switch t := a.(type) {
 406         default:
 407                 panic("unknown type in makeCopyOf")
 408         case Untransposer:
 409                 return retranspose(a, makeCopyOf(t.Untranspose()))
 410         case *Dense, *basicMatrix:
 411                 var m Dense
 412                 m.Clone(a)
 413                 return returnAs(&m, t)
 414         case *SymDense, *basicSymmetric:
 415                 n := t.(Symmetric).Symmetric()
 416                 m := NewSymDense(n, nil)
 417                 m.CopySym(t.(Symmetric))
 418                 return returnAs(m, t)
 419         case *TriDense, *basicTriangular:
 420                 n, upper := t.(Triangular).Triangle()
 421                 m := NewTriDense(n, upper, nil)
 422                 if upper {
 423                         for i := 0; i < n; i++ {
 424                                 for j := i; j < n; j++ {
 425                                         m.SetTri(i, j, t.At(i, j))
 426                                 }
 427                         }
 428                 } else {
 429                         for i := 0; i < n; i++ {
 430                                 for j := 0; j <= i; j++ {
 431                                         m.SetTri(i, j, t.At(i, j))
 432                                 }
 433                         }
 434                 }
 435                 return returnAs(m, t)
 436         case *VecDense:
 437                 m := &VecDense{
 438                         mat: blas64.Vector{
 439                                 Inc:  t.mat.Inc,
 440                                 Data: make([]float64, t.mat.Inc*(t.n-1)+1),
 441                         },
 442                         n: t.n,
 443                 }
 444                 copy(m.mat.Data, t.mat.Data)
 445                 return m
 446         case *basicVector:
 447                 m := &basicVector{
 448                         m: make([]float64, t.Len()),
 449                 }
 450                 copy(m.m, t.m)
 451                 return m
 452         }
 453 }
 454
 455 // sameType returns true if a and b have the same underlying type.
 456 func sameType(a, b Matrix) bool {
 457         return reflect.ValueOf(a).Type() == reflect.ValueOf(b).Type()
 458 }
 459
 460 // maybeSame returns true if the two matrices could be represented by the same
 461 // pointer.
 462 func maybeSame(receiver, a Matrix) bool {
 463         rr, rc := receiver.Dims()
 464         u, trans := a.(Untransposer)
 465         if trans {
 466                 a = u.Untranspose()
 467         }
 468         if !sameType(receiver, a) {
 469                 return false
 470         }
 471         ar, ac := a.Dims()
 472         if rr != ar || rc != ac {
 473                 return false
 474         }
 475         if _, ok := a.(Triangular); ok {
 476                 // They are both triangular types. The TriType needs to match
 477                 _, aKind := a.(Triangular).Triangle()
 478                 _, rKind := receiver.(Triangular).Triangle()
 479                 if aKind != rKind {
 480                         return false
 481                 }
 482         }
 483         return true
 484 }
 485
 486 // equalApprox returns whether the elements of a and b are the same to within
 487 // the tolerance. If ignoreNaN is true the test is relaxed such that NaN == NaN.
 488 func equalApprox(a, b Matrix, tol float64, ignoreNaN bool) bool {
 489         ar, ac := a.Dims()
 490         br, bc := b.Dims()
 491         if ar != br {
 492                 return false
 493         }
 494         if ac != bc {
 495                 return false
 496         }
 497         for i := 0; i < ar; i++ {
 498                 for j := 0; j < ac; j++ {
 499                         if !floats.EqualWithinAbsOrRel(a.At(i, j), b.At(i, j), tol, tol) {
 500                                 if ignoreNaN && math.IsNaN(a.At(i, j)) && math.IsNaN(b.At(i, j)) {
 501                                         continue
 502                                 }
 503                                 return false
 504                         }
 505                 }
 506         }
 507         return true
 508 }
 509
 510 // equal returns true if the matrices have equal entries.
 511 func equal(a, b Matrix) bool {
 512         ar, ac := a.Dims()
 513         br, bc := b.Dims()
 514         if ar != br {
 515                 return false
 516         }
 517         if ac != bc {
 518                 return false
 519         }
 520         for i := 0; i < ar; i++ {
 521                 for j := 0; j < ac; j++ {
 522                         if a.At(i, j) != b.At(i, j) {
 523                                 return false
 524                         }
 525                 }
 526         }
 527         return true
 528 }
 529
 530 // isDiagonal returns whether a is a diagonal matrix.
 531 func isDiagonal(a Matrix) bool {
 532         r, c := a.Dims()
 533         for i := 0; i < r; i++ {
 534                 for j := 0; j < c; j++ {
 535                         if a.At(i, j) != 0 && i != j {
 536                                 return false
 537                         }
 538                 }
 539         }
 540         return true
 541 }
 542
 543 // equalDiagonal returns whether a and b are equal on the diagonal.
 544 func equalDiagonal(a, b Matrix) bool {
 545         ar, ac := a.Dims()
 546         br, bc := a.Dims()
 547         if min(ar, ac) != min(br, bc) {
 548                 return false
 549         }
 550         for i := 0; i < min(ar, ac); i++ {
 551                 if a.At(i, i) != b.At(i, i) {
 552                         return false
 553                 }
 554         }
 555         return true
 556 }
 557
 558 // underlyingData extracts the underlying data of the matrix a.
 559 func underlyingData(a Matrix) []float64 {
 560         switch t := a.(type) {
 561         default:
 562                 panic("matrix type not implemented for extracting underlying data")
 563         case Untransposer:
 564                 return underlyingData(t.Untranspose())
 565         case *Dense:
 566                 return t.mat.Data
 567         case *SymDense:
 568                 return t.mat.Data
 569         case *TriDense:
 570                 return t.mat.Data
 571         case *VecDense:
 572                 return t.mat.Data
 573         }
 574 }
 575
 576 // testMatrices is a list of matrix types to test.
 577 // The TriDense types have actual sizes because the return from Triangular is
 578 // only valid when n == 0.
 579 var testMatrices = []Matrix{
 580         &Dense{},
 581         &SymDense{},
 582         NewTriDense(3, true, nil),
 583         NewTriDense(3, false, nil),
 584         NewVecDense(0, nil),
 585         &basicVector{},
 586         &VecDense{mat: blas64.Vector{Inc: 10}},
 587         &basicMatrix{},
 588         &basicSymmetric{},
 589         &basicTriangular{cap: 3, mat: blas64.Triangular{N: 3, Stride: 3, Uplo: blas.Upper}},
 590         &basicTriangular{cap: 3, mat: blas64.Triangular{N: 3, Stride: 3, Uplo: blas.Lower}},
 591
 592         Transpose{&Dense{}},
 593         Transpose{NewTriDense(3, true, nil)},
 594         TransposeTri{NewTriDense(3, true, nil)},
 595         Transpose{NewTriDense(3, false, nil)},
 596         TransposeTri{NewTriDense(3, false, nil)},
 597         Transpose{NewVecDense(0, nil)},
 598         Transpose{&VecDense{mat: blas64.Vector{Inc: 10}}},
 599         Transpose{&basicMatrix{}},
 600         Transpose{&basicSymmetric{}},
 601         Transpose{&basicTriangular{cap: 3, mat: blas64.Triangular{N: 3, Stride: 3, Uplo: blas.Upper}}},
 602         Transpose{&basicTriangular{cap: 3, mat: blas64.Triangular{N: 3, Stride: 3, Uplo: blas.Lower}}},
 603 }
 604
 605 var sizes = []struct {
 606         ar, ac int
 607 }{
 608         {1, 1},
 609         {1, 3},
 610         {3, 1},
 611
 612         {6, 6},
 613         {6, 11},
 614         {11, 6},
 615 }
 616
 617 func testOneInputFunc(t *testing.T,
 618         // name is the name of the function being tested.
 619         name string,
 620
 621         // f is the function being tested.
 622         f func(a Matrix) interface{},
 623
 624         // denseComparison performs the same operation, but using Dense matrices for
 625         // comparison.
 626         denseComparison func(a *Dense) interface{},
 627
 628         // sameAnswer compares the result from two different evaluations of the function
 629         // and returns true if they are the same. The specific function being tested
 630         // determines the definition of "same". It may mean identical or it may mean
 631         // approximately equal.
 632         sameAnswer func(a, b interface{}) bool,
 633
 634         // legalType returns true if the type of the input is a legal type for the
 635         // input of the function.
 636         legalType func(a Matrix) bool,
 637
 638         // legalSize returns true if the size is valid for the function.
 639         legalSize func(r, c int) bool,
 640 ) {
 641         for _, aMat := range testMatrices {
 642                 for _, test := range sizes {
 643                         // Skip the test if the argument would not be assignable to the
 644                         // method's corresponding input parameter or it is not possible
 645                         // to construct an argument of the requested size.
 646                         if !legalType(aMat) {
 647                                 continue
 648                         }
 649                         if !legalDims(aMat, test.ar, test.ac) {
 650                                 continue
 651                         }
 652                         a := makeRandOf(aMat, test.ar, test.ac)
 653
 654                         // Compute the true answer if the sizes are legal.
 655                         dimsOK := legalSize(test.ar, test.ac)
 656                         var want interface{}
 657                         if dimsOK {
 658                                 var aDense Dense
 659                                 aDense.Clone(a)
 660                                 want = denseComparison(&aDense)
 661                         }
 662                         aCopy := makeCopyOf(a)
 663                         // Test the method for a zero-value of the receiver.
 664                         aType, aTrans := untranspose(a)
 665                         errStr := fmt.Sprintf("%v(%T), size: %#v, atrans %t", name, aType, test, aTrans)
 666                         var got interface{}
 667                         panicked, err := panics(func() { got = f(a) })
 668                         if !dimsOK && !panicked {
 669                                 t.Errorf("Did not panic with illegal size: %s", errStr)
 670                                 continue
 671                         }
 672                         if dimsOK && panicked {
 673                                 t.Errorf("Panicked with legal size: %s: %v", errStr, err)
 674                                 continue
 675                         }
 676                         if !equal(a, aCopy) {
 677                                 t.Errorf("First input argument changed in call: %s", errStr)
 678                         }
 679                         if !dimsOK {
 680                                 continue
 681                         }
 682                         if !sameAnswer(want, got) {
 683                                 t.Errorf("Answer mismatch: %s", errStr)
 684                         }
 685                 }
 686         }
 687 }
 688
 689 var sizePairs = []struct {
 690         ar, ac, br, bc int
 691 }{
 692         {1, 1, 1, 1},
 693         {6, 6, 6, 6},
 694         {7, 7, 7, 7},
 695
 696         {1, 1, 1, 5},
 697         {1, 1, 5, 1},
 698         {1, 5, 1, 1},
 699         {5, 1, 1, 1},
 700
 701         {5, 5, 5, 1},
 702         {5, 5, 1, 5},
 703         {5, 1, 5, 5},
 704         {1, 5, 5, 5},
 705
 706         {6, 6, 6, 11},
 707         {6, 6, 11, 6},
 708         {6, 11, 6, 6},
 709         {11, 6, 6, 6},
 710         {11, 11, 11, 6},
 711         {11, 11, 6, 11},
 712         {11, 6, 11, 11},
 713         {6, 11, 11, 11},
 714
 715         {1, 1, 5, 5},
 716         {1, 5, 1, 5},
 717         {1, 5, 5, 1},
 718         {5, 1, 1, 5},
 719         {5, 1, 5, 1},
 720         {5, 5, 1, 1},
 721         {6, 6, 11, 11},
 722         {6, 11, 6, 11},
 723         {6, 11, 11, 6},
 724         {11, 6, 6, 11},
 725         {11, 6, 11, 6},
 726         {11, 11, 6, 6},
 727
 728         {1, 1, 17, 11},
 729         {1, 1, 11, 17},
 730         {1, 11, 1, 17},
 731         {1, 17, 1, 11},
 732         {1, 11, 17, 1},
 733         {1, 17, 11, 1},
 734         {11, 1, 1, 17},
 735         {17, 1, 1, 11},
 736         {11, 1, 17, 1},
 737         {17, 1, 11, 1},
 738         {11, 17, 1, 1},
 739         {17, 11, 1, 1},
 740
 741         {6, 6, 1, 11},
 742         {6, 6, 11, 1},
 743         {6, 11, 6, 1},
 744         {6, 1, 6, 11},
 745         {6, 11, 1, 6},
 746         {6, 1, 11, 6},
 747         {11, 6, 6, 1},
 748         {1, 6, 6, 11},
 749         {11, 6, 1, 6},
 750         {1, 6, 11, 6},
 751         {11, 1, 6, 6},
 752         {1, 11, 6, 6},
 753
 754         {6, 6, 17, 1},
 755         {6, 6, 1, 17},
 756         {6, 1, 6, 17},
 757         {6, 17, 6, 1},
 758         {6, 1, 17, 6},
 759         {6, 17, 1, 6},
 760         {1, 6, 6, 17},
 761         {17, 6, 6, 1},
 762         {1, 6, 17, 6},
 763         {17, 6, 1, 6},
 764         {1, 17, 6, 6},
 765         {17, 1, 6, 6},
 766
 767         {6, 6, 17, 11},
 768         {6, 6, 11, 17},
 769         {6, 11, 6, 17},
 770         {6, 17, 6, 11},
 771         {6, 11, 17, 6},
 772         {6, 17, 11, 6},
 773         {11, 6, 6, 17},
 774         {17, 6, 6, 11},
 775         {11, 6, 17, 6},
 776         {17, 6, 11, 6},
 777         {11, 17, 6, 6},
 778         {17, 11, 6, 6},
 779 }
 780
 781 func testTwoInputFunc(t *testing.T,
 782         // name is the name of the function being tested.
 783         name string,
 784
 785         // f is the function being tested.
 786         f func(a, b Matrix) interface{},
 787
 788         // denseComparison performs the same operation, but using Dense matrices for
 789         // comparison.
 790         denseComparison func(a, b *Dense) interface{},
 791
 792         // sameAnswer compares the result from two different evaluations of the function
 793         // and returns true if they are the same. The specific function being tested
 794         // determines the definition of "same". It may mean identical or it may mean
 795         // approximately equal.
 796         sameAnswer func(a, b interface{}) bool,
 797
 798         // legalType returns true if the types of the inputs are legal for the
 799         // input of the function.
 800         legalType func(a, b Matrix) bool,
 801
 802         // legalSize returns true if the sizes are valid for the function.
 803         legalSize func(ar, ac, br, bc int) bool,
 804 ) {
 805         for _, aMat := range testMatrices {
 806                 for _, bMat := range testMatrices {
 807                         // Loop over all of the size combinations (bigger, smaller, etc.).
 808                         for _, test := range sizePairs {
 809                                 // Skip the test if the argument would not be assignable to the
 810                                 // method's corresponding input parameter or it is not possible
 811                                 // to construct an argument of the requested size.
 812                                 if !legalType(aMat, bMat) {
 813                                         continue
 814                                 }
 815                                 if !legalDims(aMat, test.ar, test.ac) {
 816                                         continue
 817                                 }
 818                                 if !legalDims(bMat, test.br, test.bc) {
 819                                         continue
 820                                 }
 821                                 a := makeRandOf(aMat, test.ar, test.ac)
 822                                 b := makeRandOf(bMat, test.br, test.bc)
 823
 824                                 // Compute the true answer if the sizes are legal.
 825                                 dimsOK := legalSize(test.ar, test.ac, test.br, test.bc)
 826                                 var want interface{}
 827                                 if dimsOK {
 828                                         var aDense, bDense Dense
 829                                         aDense.Clone(a)
 830                                         bDense.Clone(b)
 831                                         want = denseComparison(&aDense, &bDense)
 832                                 }
 833                                 aCopy := makeCopyOf(a)
 834                                 bCopy := makeCopyOf(b)
 835                                 // Test the method for a zero-value of the receiver.
 836                                 aType, aTrans := untranspose(a)
 837                                 bType, bTrans := untranspose(b)
 838                                 errStr := fmt.Sprintf("%v(%T, %T), size: %#v, atrans %t, btrans %t", name, aType, bType, test, aTrans, bTrans)
 839                                 var got interface{}
 840                                 panicked, err := panics(func() { got = f(a, b) })
 841                                 if !dimsOK && !panicked {
 842                                         t.Errorf("Did not panic with illegal size: %s", errStr)
 843                                         continue
 844                                 }
 845                                 if dimsOK && panicked {
 846                                         t.Errorf("Panicked with legal size: %s: %v", errStr, err)
 847                                         continue
 848                                 }
 849                                 if !equal(a, aCopy) {
 850                                         t.Errorf("First input argument changed in call: %s", errStr)
 851                                 }
 852                                 if !equal(b, bCopy) {
 853                                         t.Errorf("First input argument changed in call: %s", errStr)
 854                                 }
 855                                 if !dimsOK {
 856                                         continue
 857                                 }
 858                                 if !sameAnswer(want, got) {
 859                                         t.Errorf("Answer mismatch: %s", errStr)
 860                                 }
 861                         }
 862                 }
 863         }
 864 }
 865
 866 // testOneInput tests a method that has one matrix input argument
 867 func testOneInput(t *testing.T,
 868         // name is the name of the method being tested.
 869         name string,
 870
 871         // receiver is a value of the receiver type.
 872         receiver Matrix,
 873
 874         // method is the generalized receiver.Method(a).
 875         method func(receiver, a Matrix),
 876
 877         // denseComparison performs the same operation as method, but with dense
 878         // matrices for comparison with the result.
 879         denseComparison func(receiver, a *Dense),
 880
 881         // legalTypes returns whether the concrete types in Matrix are valid for
 882         // the method.
 883         legalType func(a Matrix) bool,
 884
 885         // legalSize returns whether the matrix sizes are valid for the method.
 886         legalSize func(ar, ac int) bool,
 887
 888         // tol is the tolerance for equality when comparing method results.
 889         tol float64,
 890 ) {
 891         for _, aMat := range testMatrices {
 892                 for _, test := range sizes {
 893                         // Skip the test if the argument would not be assignable to the
 894                         // method's corresponding input parameter or it is not possible
 895                         // to construct an argument of the requested size.
 896                         if !legalType(aMat) {
 897                                 continue
 898                         }
 899                         if !legalDims(aMat, test.ar, test.ac) {
 900                                 continue
 901                         }
 902                         a := makeRandOf(aMat, test.ar, test.ac)
 903
 904                         // Compute the true answer if the sizes are legal.
 905                         dimsOK := legalSize(test.ar, test.ac)
 906                         var want Dense
 907                         if dimsOK {
 908                                 var aDense Dense
 909                                 aDense.Clone(a)
 910                                 denseComparison(&want, &aDense)
 911                         }
 912                         aCopy := makeCopyOf(a)
 913
 914                         // Test the method for a zero-value of the receiver.
 915                         aType, aTrans := untranspose(a)
 916                         errStr := fmt.Sprintf("%T.%s(%T), size: %#v, atrans %v", receiver, name, aType, test, aTrans)
 917                         zero := makeRandOf(receiver, 0, 0)
 918                         panicked, err := panics(func() { method(zero, a) })
 919                         if !dimsOK && !panicked {
 920                                 t.Errorf("Did not panic with illegal size: %s", errStr)
 921                                 continue
 922                         }
 923                         if dimsOK && panicked {
 924                                 t.Errorf("Panicked with legal size: %s: %v", errStr, err)
 925                                 continue
 926                         }
 927                         if !equal(a, aCopy) {
 928                                 t.Errorf("First input argument changed in call: %s", errStr)
 929                         }
 930                         if !dimsOK {
 931                                 continue
 932                         }
 933                         if !equalApprox(zero, &want, tol, false) {
 934                                 t.Errorf("Answer mismatch with zero receiver: %s.\nGot:\n% v\nWant:\n% v\n", errStr, Formatted(zero), Formatted(&want))
 935                                 continue
 936                         }
 937
 938                         // Test the method with a non-zero-value of the receiver.
 939                         // The receiver has been overwritten in place so use its size
 940                         // to construct a new random matrix.
 941                         rr, rc := zero.Dims()
 942                         neverZero := makeRandOf(receiver, rr, rc)
 943                         panicked, _ = panics(func() { method(neverZero, a) })
 944                         if panicked {
 945                                 t.Errorf("Panicked with non-zero receiver: %s", errStr)
 946                         }
 947                         if !equalApprox(neverZero, &want, tol, false) {
 948                                 t.Errorf("Answer mismatch non-zero receiver: %s", errStr)
 949                         }
 950
 951                         // Test with an incorrectly sized matrix.
 952                         switch receiver.(type) {
 953                         default:
 954                                 panic("matrix type not coded for incorrect receiver size")
 955                         case *Dense:
 956                                 wrongSize := makeRandOf(receiver, rr+1, rc)
 957                                 panicked, _ = panics(func() { method(wrongSize, a) })
 958                                 if !panicked {
 959                                         t.Errorf("Did not panic with wrong number of rows: %s", errStr)
 960                                 }
 961                                 wrongSize = makeRandOf(receiver, rr, rc+1)
 962                                 panicked, _ = panics(func() { method(wrongSize, a) })
 963                                 if !panicked {
 964                                         t.Errorf("Did not panic with wrong number of columns: %s", errStr)
 965                                 }
 966                         case *TriDense, *SymDense:
 967                                 // Add to the square size.
 968                                 wrongSize := makeRandOf(receiver, rr+1, rc+1)
 969                                 panicked, _ = panics(func() { method(wrongSize, a) })
 970                                 if !panicked {
 971                                         t.Errorf("Did not panic with wrong size: %s", errStr)
 972                                 }
 973                         case *VecDense:
 974                                 // Add to the column length.
 975                                 wrongSize := makeRandOf(receiver, rr+1, rc)
 976                                 panicked, _ = panics(func() { method(wrongSize, a) })
 977                                 if !panicked {
 978                                         t.Errorf("Did not panic with wrong number of rows: %s", errStr)
 979                                 }
 980                         }
 981
 982                         // The receiver and the input may share a matrix pointer
 983                         // if the type and size of the receiver and one of the
 984                         // arguments match. Test the method works properly
 985                         // when this is the case.
 986                         aMaybeSame := maybeSame(neverZero, a)
 987                         if aMaybeSame {
 988                                 aSame := makeCopyOf(a)
 989                                 receiver = aSame
 990                                 u, ok := aSame.(Untransposer)
 991                                 if ok {
 992                                         receiver = u.Untranspose()
 993                                 }
 994                                 preData := underlyingData(receiver)
 995                                 panicked, err = panics(func() { method(receiver, aSame) })
 996                                 if panicked {
 997                                         t.Errorf("Panics when a maybeSame: %s: %v", errStr, err)
 998                                 } else {
 999                                         if !equalApprox(receiver, &want, tol, false) {
1000                                                 t.Errorf("Wrong answer when a maybeSame: %s", errStr)
1001                                         }
1002                                         postData := underlyingData(receiver)
1003                                         if !floats.Equal(preData, postData) {
1004                                                 t.Errorf("Original data slice not modified when a maybeSame: %s", errStr)
1005                                         }
1006                                 }
1007                         }
1008                 }
1009         }
1010 }
1011
1012 // testTwoInput tests a method that has two input arguments.
1013 func testTwoInput(t *testing.T,
1014         // name is the name of the method being tested.
1015         name string,
1016
1017         // receiver is a value of the receiver type.
1018         receiver Matrix,
1019
1020         // method is the generalized receiver.Method(a, b).
1021         method func(receiver, a, b Matrix),
1022
1023         // denseComparison performs the same operation as method, but with dense
1024         // matrices for comparison with the result.
1025         denseComparison func(receiver, a, b *Dense),
1026
1027         // legalTypes returns whether the concrete types in Matrix are valid for
1028         // the method.
1029         legalTypes func(a, b Matrix) bool,
1030
1031         // legalSize returns whether the matrix sizes are valid for the method.
1032         legalSize func(ar, ac, br, bc int) bool,
1033
1034         // tol is the tolerance for equality when comparing method results.
1035         tol float64,
1036 ) {
1037         for _, aMat := range testMatrices {
1038                 for _, bMat := range testMatrices {
1039                         // Loop over all of the size combinations (bigger, smaller, etc.).
1040                         for _, test := range sizePairs {
1041                                 // Skip the test if any argument would not be assignable to the
1042                                 // method's corresponding input parameter or it is not possible
1043                                 // to construct an argument of the requested size.
1044                                 if !legalTypes(aMat, bMat) {
1045                                         continue
1046                                 }
1047                                 if !legalDims(aMat, test.ar, test.ac) {
1048                                         continue
1049                                 }
1050                                 if !legalDims(bMat, test.br, test.bc) {
1051                                         continue
1052                                 }
1053                                 a := makeRandOf(aMat, test.ar, test.ac)
1054                                 b := makeRandOf(bMat, test.br, test.bc)
1055
1056                                 // Compute the true answer if the sizes are legal.
1057                                 dimsOK := legalSize(test.ar, test.ac, test.br, test.bc)
1058                                 var want Dense
1059                                 if dimsOK {
1060                                         var aDense, bDense Dense
1061                                         aDense.Clone(a)
1062                                         bDense.Clone(b)
1063                                         denseComparison(&want, &aDense, &bDense)
1064                                 }
1065                                 aCopy := makeCopyOf(a)
1066                                 bCopy := makeCopyOf(b)
1067
1068                                 // Test the method for a zero-value of the receiver.
1069                                 aType, aTrans := untranspose(a)
1070                                 bType, bTrans := untranspose(b)
1071                                 errStr := fmt.Sprintf("%T.%s(%T, %T), sizes: %#v, atrans %v, btrans %v", receiver, name, aType, bType, test, aTrans, bTrans)
1072                                 zero := makeRandOf(receiver, 0, 0)
1073                                 panicked, err := panics(func() { method(zero, a, b) })
1074                                 if !dimsOK && !panicked {
1075                                         t.Errorf("Did not panic with illegal size: %s", errStr)
1076                                         continue
1077                                 }
1078                                 if dimsOK && panicked {
1079                                         t.Errorf("Panicked with legal size: %s: %v", errStr, err)
1080                                         continue
1081                                 }
1082                                 if !equal(a, aCopy) {
1083                                         t.Errorf("First input argument changed in call: %s", errStr)
1084                                 }
1085                                 if !equal(b, bCopy) {
1086                                         t.Errorf("Second input argument changed in call: %s", errStr)
1087                                 }
1088                                 if !dimsOK {
1089                                         continue
1090                                 }
1091                                 wasZero, zero := zero, nil // Nil-out zero so we detect illegal use.
1092                                 // NaN equality is allowed because of 0/0 in DivElem test.
1093                                 if !equalApprox(wasZero, &want, tol, true) {
1094                                         t.Errorf("Answer mismatch with zero receiver: %s", errStr)
1095                                         continue
1096                                 }
1097
1098                                 // Test the method with a non-zero-value of the receiver.
1099                                 // The receiver has been overwritten in place so use its size
1100                                 // to construct a new random matrix.
1101                                 rr, rc := wasZero.Dims()
1102                                 neverZero := makeRandOf(receiver, rr, rc)
1103                                 panicked, message := panics(func() { method(neverZero, a, b) })
1104                                 if panicked {
1105                                         t.Errorf("Panicked with non-zero receiver: %s: %s", errStr, message)
1106                                 }
1107                                 // NaN equality is allowed because of 0/0 in DivElem test.
1108                                 if !equalApprox(neverZero, &want, tol, true) {
1109                                         t.Errorf("Answer mismatch non-zero receiver: %s", errStr)
1110                                 }
1111
1112                                 // Test with an incorrectly sized matrix.
1113                                 switch receiver.(type) {
1114                                 default:
1115                                         panic("matrix type not coded for incorrect receiver size")
1116                                 case *Dense:
1117                                         wrongSize := makeRandOf(receiver, rr+1, rc)
1118                                         panicked, _ = panics(func() { method(wrongSize, a, b) })
1119                                         if !panicked {
1120                                                 t.Errorf("Did not panic with wrong number of rows: %s", errStr)
1121                                         }
1122                                         wrongSize = makeRandOf(receiver, rr, rc+1)
1123                                         panicked, _ = panics(func() { method(wrongSize, a, b) })
1124                                         if !panicked {
1125                                                 t.Errorf("Did not panic with wrong number of columns: %s", errStr)
1126                                         }
1127                                 case *TriDense, *SymDense:
1128                                         // Add to the square size.
1129                                         wrongSize := makeRandOf(receiver, rr+1, rc+1)
1130                                         panicked, _ = panics(func() { method(wrongSize, a, b) })
1131                                         if !panicked {
1132                                                 t.Errorf("Did not panic with wrong size: %s", errStr)
1133                                         }
1134                                 case *VecDense:
1135                                         // Add to the column length.
1136                                         wrongSize := makeRandOf(receiver, rr+1, rc)
1137                                         panicked, _ = panics(func() { method(wrongSize, a, b) })
1138                                         if !panicked {
1139                                                 t.Errorf("Did not panic with wrong number of rows: %s", errStr)
1140                                         }
1141                                 }
1142
1143                                 // The receiver and an input may share a matrix pointer
1144                                 // if the type and size of the receiver and one of the
1145                                 // arguments match. Test the method works properly
1146                                 // when this is the case.
1147                                 aMaybeSame := maybeSame(neverZero, a)
1148                                 bMaybeSame := maybeSame(neverZero, b)
1149                                 if aMaybeSame {
1150                                         aSame := makeCopyOf(a)
1151                                         receiver = aSame
1152                                         u, ok := aSame.(Untransposer)
1153                                         if ok {
1154                                                 receiver = u.Untranspose()
1155                                         }
1156                                         preData := underlyingData(receiver)
1157                                         panicked, err = panics(func() { method(receiver, aSame, b) })
1158                                         if panicked {
1159                                                 t.Errorf("Panics when a maybeSame: %s: %v", errStr, err)
1160                                         } else {
1161                                                 if !equalApprox(receiver, &want, tol, false) {
1162                                                         t.Errorf("Wrong answer when a maybeSame: %s", errStr)
1163                                                 }
1164                                                 postData := underlyingData(receiver)
1165                                                 if !floats.Equal(preData, postData) {
1166                                                         t.Errorf("Original data slice not modified when a maybeSame: %s", errStr)
1167                                                 }
1168                                         }
1169                                 }
1170                                 if bMaybeSame {
1171                                         bSame := makeCopyOf(b)
1172                                         receiver = bSame
1173                                         u, ok := bSame.(Untransposer)
1174                                         if ok {
1175                                                 receiver = u.Untranspose()
1176                                         }
1177                                         preData := underlyingData(receiver)
1178                                         panicked, err = panics(func() { method(receiver, a, bSame) })
1179                                         if panicked {
1180                                                 t.Errorf("Panics when b maybeSame: %s: %v", errStr, err)
1181                                         } else {
1182                                                 if !equalApprox(receiver, &want, tol, false) {
1183                                                         t.Errorf("Wrong answer when b maybeSame: %s", errStr)
1184                                                 }
1185                                                 postData := underlyingData(receiver)
1186                                                 if !floats.Equal(preData, postData) {
1187                                                         t.Errorf("Original data slice not modified when b maybeSame: %s", errStr)
1188                                                 }
1189                                         }
1190                                 }
1191                                 if aMaybeSame && bMaybeSame {
1192                                         aSame := makeCopyOf(a)
1193                                         receiver = aSame
1194                                         u, ok := aSame.(Untransposer)
1195                                         if ok {
1196                                                 receiver = u.Untranspose()
1197                                         }
1198                                         // Ensure that b is the correct transpose type if applicable.
1199                                         // The receiver is always a concrete type so use it.
1200                                         bSame := receiver
1201                                         u, ok = b.(Untransposer)
1202                                         if ok {
1203                                                 bSame = retranspose(b, receiver)
1204                                         }
1205                                         // Compute the real answer for this case. It is different
1206                                         // from the initial answer since now a and b have the
1207                                         // same data.
1208                                         zero = makeRandOf(wasZero, 0, 0)
1209                                         method(zero, aSame, bSame)
1210                                         wasZero, zero = zero, nil // Nil-out zero so we detect illegal use.
1211                                         preData := underlyingData(receiver)
1212                                         panicked, err = panics(func() { method(receiver, aSame, bSame) })
1213                                         if panicked {
1214                                                 t.Errorf("Panics when both maybeSame: %s: %v", errStr, err)
1215                                         } else {
1216                                                 if !equalApprox(receiver, wasZero, tol, false) {
1217                                                         t.Errorf("Wrong answer when both maybeSame: %s", errStr)
1218                                                 }
1219                                                 postData := underlyingData(receiver)
1220                                                 if !floats.Equal(preData, postData) {
1221                                                         t.Errorf("Original data slice not modified when both maybeSame: %s", errStr)
1222                                                 }
1223                                         }
1224                                 }
1225                         }
1226                 }
1227         }
1228 }