github.com/gopherd/gonum@v0.0.4/mat/symmetric_test.go (about) 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 "os" 10 "reflect" 11 "testing" 12 13 "math/rand" 14 15 "github.com/gopherd/gonum/blas" 16 "github.com/gopherd/gonum/blas/blas64" 17 "github.com/gopherd/gonum/floats/scalar" 18 ) 19 20 func TestNewSymmetric(t *testing.T) { 21 t.Parallel() 22 for i, test := range []struct { 23 data []float64 24 n int 25 mat *SymDense 26 }{ 27 { 28 data: []float64{ 29 1, 2, 3, 30 4, 5, 6, 31 7, 8, 9, 32 }, 33 n: 3, 34 mat: &SymDense{ 35 mat: blas64.Symmetric{ 36 N: 3, 37 Stride: 3, 38 Uplo: blas.Upper, 39 Data: []float64{1, 2, 3, 4, 5, 6, 7, 8, 9}, 40 }, 41 cap: 3, 42 }, 43 }, 44 } { 45 sym := NewSymDense(test.n, test.data) 46 rows, cols := sym.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(sym, test.mat) { 55 t.Errorf("unexpected data slice for test %d: got: %v want: %v", i, sym, test.mat) 56 } 57 58 m := NewDense(test.n, test.n, test.data) 59 if !reflect.DeepEqual(sym.mat.Data, m.mat.Data) { 60 t.Errorf("unexpected data slice mismatch for test %d: got: %v want: %v", i, sym.mat.Data, m.mat.Data) 61 } 62 } 63 64 panicked, message := panics(func() { NewSymDense(3, []float64{1, 2}) }) 65 if !panicked || message != ErrShape.Error() { 66 t.Error("expected panic for invalid data slice length") 67 } 68 } 69 70 func TestSymAtSet(t *testing.T) { 71 t.Parallel() 72 sym := &SymDense{ 73 mat: blas64.Symmetric{ 74 N: 3, 75 Stride: 3, 76 Uplo: blas.Upper, 77 Data: []float64{1, 2, 3, 4, 5, 6, 7, 8, 9}, 78 }, 79 cap: 3, 80 } 81 rows, cols := sym.Dims() 82 83 // Check At out of bounds 84 for _, row := range []int{-1, rows, rows + 1} { 85 panicked, message := panics(func() { sym.At(row, 0) }) 86 if !panicked || message != ErrRowAccess.Error() { 87 t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row) 88 } 89 } 90 for _, col := range []int{-1, cols, cols + 1} { 91 panicked, message := panics(func() { sym.At(0, col) }) 92 if !panicked || message != ErrColAccess.Error() { 93 t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col) 94 } 95 } 96 97 // Check Set out of bounds 98 for _, row := range []int{-1, rows, rows + 1} { 99 panicked, message := panics(func() { sym.SetSym(row, 0, 1.2) }) 100 if !panicked || message != ErrRowAccess.Error() { 101 t.Errorf("expected panic for invalid row access N=%d r=%d", rows, row) 102 } 103 } 104 for _, col := range []int{-1, cols, cols + 1} { 105 panicked, message := panics(func() { sym.SetSym(0, col, 1.2) }) 106 if !panicked || message != ErrColAccess.Error() { 107 t.Errorf("expected panic for invalid column access N=%d c=%d", cols, col) 108 } 109 } 110 111 for _, st := range []struct { 112 row, col int 113 orig, new float64 114 }{ 115 {row: 1, col: 2, orig: 6, new: 15}, 116 {row: 2, col: 1, orig: 15, new: 12}, 117 } { 118 if e := sym.At(st.row, st.col); e != st.orig { 119 t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.row, st.col, e, st.orig) 120 } 121 if e := sym.At(st.col, st.row); e != st.orig { 122 t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", st.col, st.row, e, st.orig) 123 } 124 sym.SetSym(st.row, st.col, st.new) 125 if e := sym.At(st.row, st.col); e != st.new { 126 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) 127 } 128 if e := sym.At(st.col, st.row); e != st.new { 129 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) 130 } 131 } 132 } 133 134 func TestSymDenseZero(t *testing.T) { 135 t.Parallel() 136 // Elements that equal 1 should be set to zero, elements that equal -1 137 // should remain unchanged. 138 for _, test := range []*SymDense{ 139 { 140 mat: blas64.Symmetric{ 141 Uplo: blas.Upper, 142 N: 4, 143 Stride: 5, 144 Data: []float64{ 145 1, 1, 1, 1, -1, 146 -1, 1, 1, 1, -1, 147 -1, -1, 1, 1, -1, 148 -1, -1, -1, 1, -1, 149 }, 150 }, 151 }, 152 } { 153 dataCopy := make([]float64, len(test.mat.Data)) 154 copy(dataCopy, test.mat.Data) 155 test.Zero() 156 for i, v := range test.mat.Data { 157 if dataCopy[i] != -1 && v != 0 { 158 t.Errorf("Matrix not zeroed in bounds") 159 } 160 if dataCopy[i] == -1 && v != -1 { 161 t.Errorf("Matrix zeroed out of bounds") 162 } 163 } 164 } 165 } 166 167 func TestSymDiagView(t *testing.T) { 168 t.Parallel() 169 for cas, test := range []*SymDense{ 170 NewSymDense(1, []float64{1}), 171 NewSymDense(2, []float64{1, 2, 2, 3}), 172 NewSymDense(3, []float64{1, 2, 3, 2, 4, 5, 3, 5, 6}), 173 } { 174 testDiagView(t, cas, test) 175 } 176 } 177 178 func TestSymAdd(t *testing.T) { 179 t.Parallel() 180 rnd := rand.New(rand.NewSource(1)) 181 for _, test := range []struct { 182 n int 183 }{ 184 {n: 1}, 185 {n: 2}, 186 {n: 3}, 187 {n: 4}, 188 {n: 5}, 189 {n: 10}, 190 } { 191 n := test.n 192 a := NewSymDense(n, nil) 193 for i := range a.mat.Data { 194 a.mat.Data[i] = rnd.Float64() 195 } 196 b := NewSymDense(n, nil) 197 for i := range a.mat.Data { 198 b.mat.Data[i] = rnd.Float64() 199 } 200 var m Dense 201 m.Add(a, b) 202 203 // Check with new receiver 204 var s SymDense 205 s.AddSym(a, b) 206 for i := 0; i < n; i++ { 207 for j := i; j < n; j++ { 208 want := m.At(i, j) 209 if got := s.At(i, j); got != want { 210 t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want) 211 } 212 } 213 } 214 215 // Check with equal receiver 216 s.CopySym(a) 217 s.AddSym(&s, b) 218 for i := 0; i < n; i++ { 219 for j := i; j < n; j++ { 220 want := m.At(i, j) 221 if got := s.At(i, j); got != want { 222 t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want) 223 } 224 } 225 } 226 } 227 228 method := func(receiver, a, b Matrix) { 229 type addSymer interface { 230 AddSym(a, b Symmetric) 231 } 232 rd := receiver.(addSymer) 233 rd.AddSym(a.(Symmetric), b.(Symmetric)) 234 } 235 denseComparison := func(receiver, a, b *Dense) { 236 receiver.Add(a, b) 237 } 238 testTwoInput(t, "AddSym", &SymDense{}, method, denseComparison, legalTypesSym, legalSizeSameSquare, 1e-14) 239 } 240 241 func TestCopy(t *testing.T) { 242 t.Parallel() 243 rnd := rand.New(rand.NewSource(1)) 244 for _, test := range []struct { 245 n int 246 }{ 247 {n: 1}, 248 {n: 2}, 249 {n: 3}, 250 {n: 4}, 251 {n: 5}, 252 {n: 10}, 253 } { 254 n := test.n 255 a := NewSymDense(n, nil) 256 for i := range a.mat.Data { 257 a.mat.Data[i] = rnd.Float64() 258 } 259 s := NewSymDense(n, nil) 260 s.CopySym(a) 261 for i := 0; i < n; i++ { 262 for j := i; j < n; j++ { 263 want := a.At(i, j) 264 if got := s.At(i, j); got != want { 265 t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want) 266 } 267 } 268 } 269 } 270 } 271 272 // TODO(kortschak) Roll this into testOneInput when it exists. 273 // https://github.com/gonum/matrix/issues/171 274 func TestSymCopyPanic(t *testing.T) { 275 t.Parallel() 276 var ( 277 a SymDense 278 n int 279 ) 280 m := NewSymDense(1, nil) 281 panicked, message := panics(func() { n = m.CopySym(&a) }) 282 if panicked { 283 t.Errorf("unexpected panic: %v", message) 284 } 285 if n != 0 { 286 t.Errorf("unexpected n: got: %d want: 0", n) 287 } 288 } 289 290 func TestSymRankOne(t *testing.T) { 291 t.Parallel() 292 rnd := rand.New(rand.NewSource(1)) 293 const tol = 1e-15 294 295 for _, test := range []struct { 296 n int 297 }{ 298 {n: 1}, 299 {n: 2}, 300 {n: 3}, 301 {n: 4}, 302 {n: 5}, 303 {n: 10}, 304 } { 305 n := test.n 306 alpha := 2.0 307 a := NewSymDense(n, nil) 308 for i := range a.mat.Data { 309 a.mat.Data[i] = rnd.Float64() 310 } 311 x := make([]float64, n) 312 for i := range x { 313 x[i] = rnd.Float64() 314 } 315 316 xMat := NewDense(n, 1, x) 317 var m Dense 318 m.Mul(xMat, xMat.T()) 319 m.Scale(alpha, &m) 320 m.Add(&m, a) 321 322 // Check with new receiver 323 s := NewSymDense(n, nil) 324 s.SymRankOne(a, alpha, NewVecDense(len(x), x)) 325 for i := 0; i < n; i++ { 326 for j := i; j < n; j++ { 327 want := m.At(i, j) 328 if got := s.At(i, j); !scalar.EqualWithinAbsOrRel(got, want, tol, tol) { 329 t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want) 330 } 331 } 332 } 333 334 // Check with reused receiver 335 copy(s.mat.Data, a.mat.Data) 336 s.SymRankOne(s, alpha, NewVecDense(len(x), x)) 337 for i := 0; i < n; i++ { 338 for j := i; j < n; j++ { 339 want := m.At(i, j) 340 if got := s.At(i, j); !scalar.EqualWithinAbsOrRel(got, want, tol, tol) { 341 t.Errorf("unexpected value for At(%d, %d): got: %v want: %v", i, j, got, want) 342 } 343 } 344 } 345 } 346 347 alpha := 3.0 348 method := func(receiver, a, b Matrix) { 349 type SymRankOner interface { 350 SymRankOne(a Symmetric, alpha float64, x Vector) 351 } 352 rd := receiver.(SymRankOner) 353 rd.SymRankOne(a.(Symmetric), alpha, b.(Vector)) 354 } 355 denseComparison := func(receiver, a, b *Dense) { 356 var tmp Dense 357 tmp.Mul(b, b.T()) 358 tmp.Scale(alpha, &tmp) 359 receiver.Add(a, &tmp) 360 } 361 legalTypes := func(a, b Matrix) bool { 362 _, ok := a.(Symmetric) 363 if !ok { 364 return false 365 } 366 _, ok = b.(Vector) 367 return ok 368 } 369 legalSize := func(ar, ac, br, bc int) bool { 370 if ar != ac { 371 return false 372 } 373 return br == ar 374 } 375 testTwoInput(t, "SymRankOne", &SymDense{}, method, denseComparison, legalTypes, legalSize, 1e-14) 376 } 377 378 func TestIssue250SymRankOne(t *testing.T) { 379 t.Parallel() 380 x := NewVecDense(5, []float64{1, 2, 3, 4, 5}) 381 var s1, s2 SymDense 382 s1.SymRankOne(NewSymDense(5, nil), 1, x) 383 s2.SymRankOne(NewSymDense(5, nil), 1, x) 384 s2.SymRankOne(NewSymDense(5, nil), 1, x) 385 if !Equal(&s1, &s2) { 386 t.Error("unexpected result from repeat") 387 } 388 } 389 390 func TestRankTwo(t *testing.T) { 391 t.Parallel() 392 rnd := rand.New(rand.NewSource(1)) 393 for _, test := range []struct { 394 n int 395 }{ 396 {n: 1}, 397 {n: 2}, 398 {n: 3}, 399 {n: 4}, 400 {n: 5}, 401 {n: 10}, 402 } { 403 n := test.n 404 alpha := 2.0 405 a := NewSymDense(n, nil) 406 for i := range a.mat.Data { 407 a.mat.Data[i] = rnd.Float64() 408 } 409 x := make([]float64, n) 410 y := make([]float64, n) 411 for i := range x { 412 x[i] = rnd.Float64() 413 y[i] = rnd.Float64() 414 } 415 416 xMat := NewDense(n, 1, x) 417 yMat := NewDense(n, 1, y) 418 var m Dense 419 m.Mul(xMat, yMat.T()) 420 var tmp Dense 421 tmp.Mul(yMat, xMat.T()) 422 m.Add(&m, &tmp) 423 m.Scale(alpha, &m) 424 m.Add(&m, a) 425 426 // Check with new receiver 427 s := NewSymDense(n, nil) 428 s.RankTwo(a, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y)) 429 for i := 0; i < n; i++ { 430 for j := i; j < n; j++ { 431 if !scalar.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) { 432 t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j)) 433 } 434 } 435 } 436 437 // Check with reused receiver 438 copy(s.mat.Data, a.mat.Data) 439 s.RankTwo(s, alpha, NewVecDense(len(x), x), NewVecDense(len(y), y)) 440 for i := 0; i < n; i++ { 441 for j := i; j < n; j++ { 442 if !scalar.EqualWithinAbsOrRel(s.At(i, j), m.At(i, j), 1e-14, 1e-14) { 443 t.Errorf("unexpected element value at (%d,%d): got: %f want: %f", i, j, m.At(i, j), s.At(i, j)) 444 } 445 } 446 } 447 } 448 } 449 450 func TestSymRankK(t *testing.T) { 451 t.Parallel() 452 alpha := 3.0 453 method := func(receiver, a, b Matrix) { 454 type SymRankKer interface { 455 SymRankK(a Symmetric, alpha float64, x Matrix) 456 } 457 rd := receiver.(SymRankKer) 458 rd.SymRankK(a.(Symmetric), alpha, b) 459 } 460 denseComparison := func(receiver, a, b *Dense) { 461 var tmp Dense 462 tmp.Mul(b, b.T()) 463 tmp.Scale(alpha, &tmp) 464 receiver.Add(a, &tmp) 465 } 466 legalTypes := func(a, b Matrix) bool { 467 _, ok := a.(Symmetric) 468 return ok 469 } 470 legalSize := func(ar, ac, br, bc int) bool { 471 if ar != ac { 472 return false 473 } 474 return br == ar 475 } 476 testTwoInput(t, "SymRankK", &SymDense{}, method, denseComparison, legalTypes, legalSize, 1e-14) 477 } 478 479 func TestSymOuterK(t *testing.T) { 480 t.Parallel() 481 for _, f := range []float64{0.5, 1, 3} { 482 method := func(receiver, x Matrix) { 483 type SymOuterKer interface { 484 SymOuterK(alpha float64, x Matrix) 485 } 486 rd := receiver.(SymOuterKer) 487 rd.SymOuterK(f, x) 488 } 489 denseComparison := func(receiver, x *Dense) { 490 receiver.Mul(x, x.T()) 491 receiver.Scale(f, receiver) 492 } 493 testOneInput(t, "SymOuterK", &SymDense{}, method, denseComparison, isAnyType, isAnySize, 1e-14) 494 } 495 } 496 497 func TestIssue250SymOuterK(t *testing.T) { 498 t.Parallel() 499 x := NewVecDense(5, []float64{1, 2, 3, 4, 5}) 500 var s1, s2 SymDense 501 s1.SymOuterK(1, x) 502 s2.SymOuterK(1, x) 503 s2.SymOuterK(1, x) 504 if !Equal(&s1, &s2) { 505 t.Error("unexpected result from repeat") 506 } 507 } 508 509 func TestScaleSym(t *testing.T) { 510 t.Parallel() 511 for _, f := range []float64{0.5, 1, 3} { 512 method := func(receiver, a Matrix) { 513 type ScaleSymer interface { 514 ScaleSym(f float64, a Symmetric) 515 } 516 rd := receiver.(ScaleSymer) 517 rd.ScaleSym(f, a.(Symmetric)) 518 } 519 denseComparison := func(receiver, a *Dense) { 520 receiver.Scale(f, a) 521 } 522 testOneInput(t, "ScaleSym", &SymDense{}, method, denseComparison, legalTypeSym, isSquare, 1e-14) 523 } 524 } 525 526 func TestSubsetSym(t *testing.T) { 527 t.Parallel() 528 for _, test := range []struct { 529 a *SymDense 530 dims []int 531 ans *SymDense 532 }{ 533 { 534 a: NewSymDense(3, []float64{ 535 1, 2, 3, 536 0, 4, 5, 537 0, 0, 6, 538 }), 539 dims: []int{0, 2}, 540 ans: NewSymDense(2, []float64{ 541 1, 3, 542 0, 6, 543 }), 544 }, 545 { 546 a: NewSymDense(3, []float64{ 547 1, 2, 3, 548 0, 4, 5, 549 0, 0, 6, 550 }), 551 dims: []int{2, 0}, 552 ans: NewSymDense(2, []float64{ 553 6, 3, 554 0, 1, 555 }), 556 }, 557 { 558 a: NewSymDense(3, []float64{ 559 1, 2, 3, 560 0, 4, 5, 561 0, 0, 6, 562 }), 563 dims: []int{1, 1, 1}, 564 ans: NewSymDense(3, []float64{ 565 4, 4, 4, 566 0, 4, 4, 567 0, 0, 4, 568 }), 569 }, 570 } { 571 var s SymDense 572 s.SubsetSym(test.a, test.dims) 573 if !Equal(&s, test.ans) { 574 t.Errorf("SubsetSym mismatch dims %v\nGot:\n% v\nWant:\n% v\n", test.dims, s, test.ans) 575 } 576 } 577 578 dims := []int{0, 2} 579 maxDim := dims[0] 580 for _, v := range dims { 581 if maxDim < v { 582 maxDim = v 583 } 584 } 585 method := func(receiver, a Matrix) { 586 type SubsetSymer interface { 587 SubsetSym(a Symmetric, set []int) 588 } 589 rd := receiver.(SubsetSymer) 590 rd.SubsetSym(a.(Symmetric), dims) 591 } 592 denseComparison := func(receiver, a *Dense) { 593 *receiver = *NewDense(len(dims), len(dims), nil) 594 sz := len(dims) 595 for i := 0; i < sz; i++ { 596 for j := 0; j < sz; j++ { 597 receiver.Set(i, j, a.At(dims[i], dims[j])) 598 } 599 } 600 } 601 legalSize := func(ar, ac int) bool { 602 return ar == ac && ar > maxDim 603 } 604 605 testOneInput(t, "SubsetSym", &SymDense{}, method, denseComparison, legalTypeSym, legalSize, 0) 606 } 607 608 func TestViewGrowSquare(t *testing.T) { 609 t.Parallel() 610 // n is the size of the original SymDense. 611 // The first view uses start1, span1. The second view uses start2, span2 on 612 // the first view. 613 for _, test := range []struct { 614 n, start1, span1, start2, span2 int 615 }{ 616 {10, 0, 10, 0, 10}, 617 {10, 0, 8, 0, 8}, 618 {10, 2, 8, 0, 6}, 619 {10, 2, 7, 4, 2}, 620 {10, 2, 6, 0, 5}, 621 } { 622 n := test.n 623 s := NewSymDense(n, nil) 624 for i := 0; i < n; i++ { 625 for j := i; j < n; j++ { 626 s.SetSym(i, j, float64((i+1)*n+j+1)) 627 } 628 } 629 630 // Take a subset and check the view matches. 631 start1 := test.start1 632 span1 := test.span1 633 v := s.sliceSym(start1, start1+span1) 634 for i := 0; i < span1; i++ { 635 for j := i; j < span1; j++ { 636 if v.At(i, j) != s.At(start1+i, start1+j) { 637 t.Errorf("View mismatch") 638 } 639 } 640 } 641 642 start2 := test.start2 643 span2 := test.span2 644 v2 := v.SliceSym(start2, start2+span2).(*SymDense) 645 646 for i := 0; i < span2; i++ { 647 for j := i; j < span2; j++ { 648 if v2.At(i, j) != s.At(start1+start2+i, start1+start2+j) { 649 t.Errorf("Second view mismatch") 650 } 651 } 652 } 653 654 // Check that a write to the view is reflected in the original. 655 v2.SetSym(0, 0, 1.2) 656 if s.At(start1+start2, start1+start2) != 1.2 { 657 t.Errorf("Write to view not reflected in original") 658 } 659 660 // Grow the matrix back to the original view 661 gn := n - start1 - start2 662 g := v2.GrowSym(gn - v2.SymmetricDim()).(*SymDense) 663 g.SetSym(1, 1, 2.2) 664 665 for i := 0; i < gn; i++ { 666 for j := 0; j < gn; j++ { 667 if g.At(i, j) != s.At(start1+start2+i, start1+start2+j) { 668 t.Errorf("Grow mismatch") 669 670 fmt.Printf("g=\n% v\n", Formatted(g)) 671 fmt.Printf("s=\n% v\n", Formatted(s)) 672 os.Exit(1) 673 } 674 } 675 } 676 677 // View g, then grow it and make sure all the elements were copied. 678 gv := g.SliceSym(0, gn-1).(*SymDense) 679 680 gg := gv.GrowSym(2) 681 for i := 0; i < gn; i++ { 682 for j := 0; j < gn; j++ { 683 if g.At(i, j) != gg.At(i, j) { 684 t.Errorf("Expand mismatch") 685 } 686 } 687 } 688 689 s.Reset() 690 rg := s.GrowSym(n).(*SymDense) 691 if rg.mat.Stride < n { 692 t.Errorf("unexpected stride after GrowSym on empty matrix: got:%d want >= %d", rg.mat.Stride, n) 693 } 694 } 695 } 696 697 func TestPowPSD(t *testing.T) { 698 t.Parallel() 699 for cas, test := range []struct { 700 a *SymDense 701 pow float64 702 ans *SymDense 703 }{ 704 // Comparison with Matlab. 705 { 706 a: NewSymDense(2, []float64{10, 5, 5, 12}), 707 pow: 0.5, 708 ans: NewSymDense(2, []float64{3.065533767740645, 0.776210486171016, 0.776210486171016, 3.376017962209052}), 709 }, 710 { 711 a: NewSymDense(2, []float64{11, -1, -1, 8}), 712 pow: 0.5, 713 ans: NewSymDense(2, []float64{3.312618742210524, -0.162963396980939, -0.162963396980939, 2.823728551267709}), 714 }, 715 { 716 a: NewSymDense(2, []float64{10, 5, 5, 12}), 717 pow: -0.5, 718 ans: NewSymDense(2, []float64{0.346372134547712, -0.079637515547296, -0.079637515547296, 0.314517128328794}), 719 }, 720 { 721 a: NewSymDense(3, []float64{15, -1, -3, -1, 8, 6, -3, 6, 14}), 722 pow: 0.6, 723 ans: NewSymDense(3, []float64{ 724 5.051214323034288, -0.163162161893975, -0.612153996497505, 725 -0.163162161893976, 3.283474884617009, 1.432842761381493, 726 -0.612153996497505, 1.432842761381494, 4.695873060862573, 727 }), 728 }, 729 } { 730 var s SymDense 731 err := s.PowPSD(test.a, test.pow) 732 if err != nil { 733 panic("bad test") 734 } 735 if !EqualApprox(&s, test.ans, 1e-10) { 736 t.Errorf("Case %d, pow mismatch", cas) 737 fmt.Println(Formatted(&s)) 738 fmt.Println(Formatted(test.ans)) 739 } 740 } 741 742 // Compare with Dense.Pow 743 rnd := rand.New(rand.NewSource(1)) 744 for dim := 2; dim < 10; dim++ { 745 for pow := 2; pow < 6; pow++ { 746 a := NewDense(dim, dim, nil) 747 for i := 0; i < dim; i++ { 748 for j := 0; j < dim; j++ { 749 a.Set(i, j, rnd.Float64()) 750 } 751 } 752 var mat SymDense 753 mat.SymOuterK(1, a) 754 755 var sym SymDense 756 err := sym.PowPSD(&mat, float64(pow)) 757 if err != nil { 758 t.Errorf("unexpected error: %v", err) 759 } 760 761 var dense Dense 762 dense.Pow(&mat, pow) 763 764 if !EqualApprox(&sym, &dense, 1e-10) { 765 t.Errorf("Dim %d: pow mismatch", dim) 766 } 767 } 768 } 769 } 770 771 func BenchmarkSymSum1000(b *testing.B) { symSumBench(b, 1000) } 772 773 var symSumForBench float64 774 775 func symSumBench(b *testing.B, size int) { 776 src := rand.NewSource(1) 777 a := randSymDense(size, src) 778 b.ResetTimer() 779 for i := 0; i < b.N; i++ { 780 symSumForBench = Sum(a) 781 } 782 } 783 784 func randSymDense(size int, src rand.Source) *SymDense { 785 rnd := rand.New(src) 786 backData := make([]float64, size*size) 787 for i := 0; i < size; i++ { 788 backData[i*size+i] = rnd.Float64() 789 for j := i + 1; j < size; j++ { 790 v := rnd.Float64() 791 backData[i*size+j] = v 792 backData[j*size+i] = v 793 } 794 } 795 s := NewSymDense(size, backData) 796 return s 797 }