vendor/gonum.org/v1/gonum/mat/triangular_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         "math"
   9         "reflect"
  10         "testing"
  11
  12         "golang.org/x/exp/rand"
  13
  14         "gonum.org/v1/gonum/blas"
  15         "gonum.org/v1/gonum/blas/blas64"
  16 )
  17
  18 func TestNewTriangular(t *testing.T) {
  19         for i, test := range []struct {
  20                 data []float64
  21                 n    int
  22                 kind TriKind
  23                 mat  *TriDense
  24         }{
  25                 {
  26                         data: []float64{
  27                                 1, 2, 3,
  28                                 4, 5, 6,
  29                                 7, 8, 9,
  30                         },
  31                         n:    3,
  32                         kind: Upper,
  33                         mat: &TriDense{
  34                                 mat: blas64.Triangular{
  35                                         N:      3,
  36                                         Stride: 3,
  37                                         Uplo:   blas.Upper,
  38                                         Data:   []float64{1, 2, 3, 4, 5, 6, 7, 8, 9},
  39                                         Diag:   blas.NonUnit,
  40                                 },
  41                                 cap: 3,
  42                         },
  43                 },
  44         } {
  45                 tri := NewTriDense(test.n, test.kind, test.data)
  46                 rows, cols := tri.Dims()
  47
  48                 if rows != test.n {
  49                         t.Errorf("unexpected number of rows for test %d: got: %d want: %d", i, rows, test.n)
  50                 }
  51                 if cols != test.n {
  52                         t.Errorf("unexpected number of cols for test %d: got: %d want: %d", i, cols, test.n)
  53                 }
  54                 if !reflect.DeepEqual(tri, test.mat) {
  55                         t.Errorf("unexpected data slice for test %d: got: %v want: %v", i, tri, test.mat)
  56                 }
  57         }
  58
  59         for _, kind := range []TriKind{Lower, Upper} {
  60                 panicked, message := panics(func() { NewTriDense(3, kind, []float64{1, 2}) })
  61                 if !panicked || message != ErrShape.Error() {
  62                         t.Errorf("expected panic for invalid data slice length for upper=%t", kind)
  63                 }
  64         }
  65 }
  66
  67 func TestTriAtSet(t *testing.T) {
  68         tri := &TriDense{
  69                 mat: blas64.Triangular{
  70                         N:      3,
  71                         Stride: 3,
  72                         Uplo:   blas.Upper,
  73                         Data:   []float64{1, 2, 3, 4, 5, 6, 7, 8, 9},
  74                         Diag:   blas.NonUnit,
  75                 },
  76                 cap: 3,
  77         }
  78
  79         rows, cols := tri.Dims()
  80
  81         // Check At out of bounds
  82         for _, row := range []int{-1, rows, rows + 1} {
  83                 panicked, message := panics(func() { tri.At(row, 0) })
  84                 if !panicked || message != ErrRowAccess.Error() {
  85                         t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
  86                 }
  87         }
  88         for _, col := range []int{-1, cols, cols + 1} {
  89                 panicked, message := panics(func() { tri.At(0, col) })
  90                 if !panicked || message != ErrColAccess.Error() {
  91                         t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
  92                 }
  93         }
  94
  95         // Check Set out of bounds
  96         for _, row := range []int{-1, rows, rows + 1} {
  97                 panicked, message := panics(func() { tri.SetTri(row, 0, 1.2) })
  98                 if !panicked || message != ErrRowAccess.Error() {
  99                         t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row)
 100                 }
 101         }
 102         for _, col := range []int{-1, cols, cols + 1} {
 103                 panicked, message := panics(func() { tri.SetTri(0, col, 1.2) })
 104                 if !panicked || message != ErrColAccess.Error() {
 105                         t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col)
 106                 }
 107         }
 108
 109         for _, st := range []struct {
 110                 row, col int
 111                 uplo     blas.Uplo
 112         }{
 113                 {row: 2, col: 1, uplo: blas.Upper},
 114                 {row: 1, col: 2, uplo: blas.Lower},
 115         } {
 116                 tri.mat.Uplo = st.uplo
 117                 panicked, message := panics(func() { tri.SetTri(st.row, st.col, 1.2) })
 118                 if !panicked || message != ErrTriangleSet.Error() {
 119                         t.Errorf("expected panic for %+v", st)
 120                 }
 121         }
 122
 123         for _, st := range []struct {
 124                 row, col  int
 125                 uplo      blas.Uplo
 126                 orig, new float64
 127         }{
 128                 {row: 2, col: 1, uplo: blas.Lower, orig: 8, new: 15},
 129                 {row: 1, col: 2, uplo: blas.Upper, orig: 6, new: 15},
 130         } {
 131                 tri.mat.Uplo = st.uplo
 132                 if e := tri.At(st.row, st.col); e != st.orig {
 133                         t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.row, st.col, e, st.orig)
 134                 }
 135                 tri.SetTri(st.row, st.col, st.new)
 136                 if e := tri.At(st.row, st.col); e != st.new {
 137                         t.Errorf("unexpected value for At(%d, %d) after SetTri(%[1]d, %d, %v): got: %v want: %[3]v", st.row, st.col, st.new, e)
 138                 }
 139         }
 140 }
 141
 142 func TestTriDenseCopy(t *testing.T) {
 143         for i := 0; i < 100; i++ {
 144                 size := rand.Intn(100)
 145                 r, err := randDense(size, 0.9, rand.NormFloat64)
 146                 if size == 0 {
 147                         if err != ErrZeroLength {
 148                                 t.Fatalf("expected error %v: got: %v", ErrZeroLength, err)
 149                         }
 150                         continue
 151                 }
 152                 if err != nil {
 153                         t.Fatalf("unexpected error: %v", err)
 154                 }
 155
 156                 u := NewTriDense(size, true, nil)
 157                 l := NewTriDense(size, false, nil)
 158
 159                 for _, typ := range []Matrix{r, (*basicMatrix)(r)} {
 160                         for j := range u.mat.Data {
 161                                 u.mat.Data[j] = math.NaN()
 162                                 l.mat.Data[j] = math.NaN()
 163                         }
 164                         u.Copy(typ)
 165                         l.Copy(typ)
 166                         for m := 0; m < size; m++ {
 167                                 for n := 0; n < size; n++ {
 168                                         want := typ.At(m, n)
 169                                         switch {
 170                                         case m < n: // Upper triangular matrix.
 171                                                 if got := u.At(m, n); got != want {
 172                                                         t.Errorf("unexpected upper value for At(%d, %d) for test %d: got: %v want: %v", m, n, i, got, want)
 173                                                 }
 174                                         case m == n: // Diagonal matrix.
 175                                                 if got := u.At(m, n); got != want {
 176                                                         t.Errorf("unexpected upper value for At(%d, %d) for test %d: got: %v want: %v", m, n, i, got, want)
 177                                                 }
 178                                                 if got := l.At(m, n); got != want {
 179                                                         t.Errorf("unexpected diagonal value for At(%d, %d) for test %d: got: %v want: %v", m, n, i, got, want)
 180                                                 }
 181                                         case m < n: // Lower triangular matrix.
 182                                                 if got := l.At(m, n); got != want {
 183                                                         t.Errorf("unexpected lower value for At(%d, %d) for test %d: got: %v want: %v", m, n, i, got, want)
 184                                                 }
 185                                         }
 186                                 }
 187                         }
 188                 }
 189         }
 190 }
 191
 192 func TestTriTriDenseCopy(t *testing.T) {
 193         for i := 0; i < 100; i++ {
 194                 size := rand.Intn(100)
 195                 r, err := randDense(size, 1, rand.NormFloat64)
 196                 if size == 0 {
 197                         if err != ErrZeroLength {
 198                                 t.Fatalf("expected error %v: got: %v", ErrZeroLength, err)
 199                         }
 200                         continue
 201                 }
 202                 if err != nil {
 203                         t.Fatalf("unexpected error: %v", err)
 204                 }
 205
 206                 ur := NewTriDense(size, true, nil)
 207                 lr := NewTriDense(size, false, nil)
 208
 209                 ur.Copy(r)
 210                 lr.Copy(r)
 211
 212                 u := NewTriDense(size, true, nil)
 213                 u.Copy(ur)
 214                 if !equal(u, ur) {
 215                         t.Fatal("unexpected result for U triangle copy of U triangle: not equal")
 216                 }
 217
 218                 l := NewTriDense(size, false, nil)
 219                 l.Copy(lr)
 220                 if !equal(l, lr) {
 221                         t.Fatal("unexpected result for L triangle copy of L triangle: not equal")
 222                 }
 223
 224                 zero(u.mat.Data)
 225                 u.Copy(lr)
 226                 if !isDiagonal(u) {
 227                         t.Fatal("unexpected result for U triangle copy of L triangle: off diagonal non-zero element")
 228                 }
 229                 if !equalDiagonal(u, lr) {
 230                         t.Fatal("unexpected result for U triangle copy of L triangle: diagonal not equal")
 231                 }
 232
 233                 zero(l.mat.Data)
 234                 l.Copy(ur)
 235                 if !isDiagonal(l) {
 236                         t.Fatal("unexpected result for L triangle copy of U triangle: off diagonal non-zero element")
 237                 }
 238                 if !equalDiagonal(l, ur) {
 239                         t.Fatal("unexpected result for L triangle copy of U triangle: diagonal not equal")
 240                 }
 241         }
 242 }
 243
 244 func TestTriInverse(t *testing.T) {
 245         for _, kind := range []TriKind{Upper, Lower} {
 246                 for _, n := range []int{1, 3, 5, 9} {
 247                         data := make([]float64, n*n)
 248                         for i := range data {
 249                                 data[i] = rand.NormFloat64()
 250                         }
 251                         a := NewTriDense(n, kind, data)
 252                         var tr TriDense
 253                         err := tr.InverseTri(a)
 254                         if err != nil {
 255                                 t.Errorf("Bad test: %s", err)
 256                         }
 257                         var d Dense
 258                         d.Mul(a, &tr)
 259                         if !equalApprox(eye(n), &d, 1e-8, false) {
 260                                 var diff Dense
 261                                 diff.Sub(eye(n), &d)
 262                                 t.Errorf("Tri times inverse is not identity. Norm of difference: %v", Norm(&diff, 2))
 263                         }
 264                 }
 265         }
 266 }
 267
 268 func TestTriMul(t *testing.T) {
 269         method := func(receiver, a, b Matrix) {
 270                 type MulTrier interface {
 271                         MulTri(a, b Triangular)
 272                 }
 273                 receiver.(MulTrier).MulTri(a.(Triangular), b.(Triangular))
 274         }
 275         denseComparison := func(receiver, a, b *Dense) {
 276                 receiver.Mul(a, b)
 277         }
 278         legalSizeTriMul := func(ar, ac, br, bc int) bool {
 279                 // Need both to be square and the sizes to be the same
 280                 return ar == ac && br == bc && ar == br
 281         }
 282
 283         // The legal types are triangles with the same TriKind.
 284         // legalTypesTri returns whether both input arguments are Triangular.
 285         legalTypes := func(a, b Matrix) bool {
 286                 at, ok := a.(Triangular)
 287                 if !ok {
 288                         return false
 289                 }
 290                 bt, ok := b.(Triangular)
 291                 if !ok {
 292                         return false
 293                 }
 294                 _, ak := at.Triangle()
 295                 _, bk := bt.Triangle()
 296                 return ak == bk
 297         }
 298         legalTypesLower := func(a, b Matrix) bool {
 299                 legal := legalTypes(a, b)
 300                 if !legal {
 301                         return false
 302                 }
 303                 _, kind := a.(Triangular).Triangle()
 304                 r := kind == Lower
 305                 return r
 306         }
 307         receiver := NewTriDense(3, Lower, nil)
 308         testTwoInput(t, "TriMul", receiver, method, denseComparison, legalTypesLower, legalSizeTriMul, 1e-14)
 309
 310         legalTypesUpper := func(a, b Matrix) bool {
 311                 legal := legalTypes(a, b)
 312                 if !legal {
 313                         return false
 314                 }
 315                 _, kind := a.(Triangular).Triangle()
 316                 r := kind == Upper
 317                 return r
 318         }
 319         receiver = NewTriDense(3, Upper, nil)
 320         testTwoInput(t, "TriMul", receiver, method, denseComparison, legalTypesUpper, legalSizeTriMul, 1e-14)
 321 }
 322
 323 func TestScaleTri(t *testing.T) {
 324         for _, f := range []float64{0.5, 1, 3} {
 325                 method := func(receiver, a Matrix) {
 326                         type ScaleTrier interface {
 327                                 ScaleTri(f float64, a Triangular)
 328                         }
 329                         rd := receiver.(ScaleTrier)
 330                         rd.ScaleTri(f, a.(Triangular))
 331                 }
 332                 denseComparison := func(receiver, a *Dense) {
 333                         receiver.Scale(f, a)
 334                 }
 335                 testOneInput(t, "ScaleTriUpper", NewTriDense(3, Upper, nil), method, denseComparison, legalTypeTriUpper, isSquare, 1e-14)
 336                 testOneInput(t, "ScaleTriLower", NewTriDense(3, Lower, nil), method, denseComparison, legalTypeTriLower, isSquare, 1e-14)
 337         }
 338 }
 339
 340 func TestCopySymIntoTriangle(t *testing.T) {
 341         nan := math.NaN()
 342         for tc, test := range []struct {
 343                 n     int
 344                 sUplo blas.Uplo
 345                 s     []float64
 346
 347                 tUplo TriKind
 348                 want  []float64
 349         }{
 350                 {
 351                         n:     3,
 352                         sUplo: blas.Upper,
 353                         s: []float64{
 354                                 1, 2, 3,
 355                                 nan, 4, 5,
 356                                 nan, nan, 6,
 357                         },
 358                         tUplo: Upper,
 359                         want: []float64{
 360                                 1, 2, 3,
 361                                 0, 4, 5,
 362                                 0, 0, 6,
 363                         },
 364                 },
 365                 {
 366                         n:     3,
 367                         sUplo: blas.Lower,
 368                         s: []float64{
 369                                 1, nan, nan,
 370                                 2, 3, nan,
 371                                 4, 5, 6,
 372                         },
 373                         tUplo: Upper,
 374                         want: []float64{
 375                                 1, 2, 4,
 376                                 0, 3, 5,
 377                                 0, 0, 6,
 378                         },
 379                 },
 380                 {
 381                         n:     3,
 382                         sUplo: blas.Upper,
 383                         s: []float64{
 384                                 1, 2, 3,
 385                                 nan, 4, 5,
 386                                 nan, nan, 6,
 387                         },
 388                         tUplo: Lower,
 389                         want: []float64{
 390                                 1, 0, 0,
 391                                 2, 4, 0,
 392                                 3, 5, 6,
 393                         },
 394                 },
 395                 {
 396                         n:     3,
 397                         sUplo: blas.Lower,
 398                         s: []float64{
 399                                 1, nan, nan,
 400                                 2, 3, nan,
 401                                 4, 5, 6,
 402                         },
 403                         tUplo: Lower,
 404                         want: []float64{
 405                                 1, 0, 0,
 406                                 2, 3, 0,
 407                                 4, 5, 6,
 408                         },
 409                 },
 410         } {
 411                 n := test.n
 412                 s := NewSymDense(n, test.s)
 413                 // For the purpose of the test, break the assumption that
 414                 // symmetric is stored in the upper triangle (only when S is
 415                 // RawSymmetricer).
 416                 s.mat.Uplo = test.sUplo
 417
 418                 t1 := NewTriDense(n, test.tUplo, nil)
 419                 copySymIntoTriangle(t1, s)
 420
 421                 equal := true
 422         loop1:
 423                 for i := 0; i < n; i++ {
 424                         for j := 0; j < n; j++ {
 425                                 if t1.At(i, j) != test.want[i*n+j] {
 426                                         equal = false
 427                                         break loop1
 428                                 }
 429                         }
 430                 }
 431                 if !equal {
 432                         t.Errorf("Case %v: unexpected T when S is RawSymmetricer", tc)
 433                 }
 434
 435                 if test.sUplo == blas.Lower {
 436                         continue
 437                 }
 438
 439                 sb := (basicSymmetric)(*s)
 440                 t2 := NewTriDense(n, test.tUplo, nil)
 441                 copySymIntoTriangle(t2, &sb)
 442                 equal = true
 443         loop2:
 444                 for i := 0; i < n; i++ {
 445                         for j := 0; j < n; j++ {
 446                                 if t1.At(i, j) != test.want[i*n+j] {
 447                                         equal = false
 448                                         break loop2
 449                                 }
 450                         }
 451                 }
 452                 if !equal {
 453                         t.Errorf("Case %v: unexpected T when S is not RawSymmetricer", tc)
 454                 }
 455         }
 456 }