php-ai /
php-ml
| Conditions | 71 |
| Paths | > 20000 |
| Total Lines | 376 |
| Code Lines | 270 |
| Lines | 45 |
| Ratio | 11.97 % |
| Changes | 0 | ||
Small methods make your code easier to understand, in particular if combined with a good name. Besides, if your method is small, finding a good name is usually much easier.
For example, if you find yourself adding comments to a method's body, this is usually a good sign to extract the commented part to a new method, and use the comment as a starting point when coming up with a good name for this new method.
Commonly applied refactorings include:
If many parameters/temporary variables are present:
| 1 | <?php |
||
| 409 | private function hqr2() |
||
| 410 | { |
||
| 411 | // Initialize |
||
| 412 | $nn = $this->n; |
||
| 413 | $n = $nn - 1; |
||
| 414 | $low = 0; |
||
| 415 | $high = $nn - 1; |
||
| 416 | $eps = pow(2.0, -52.0); |
||
| 417 | $exshift = 0.0; |
||
| 418 | $p = $q = $r = $s = $z = 0; |
||
| 419 | // Store roots isolated by balanc and compute matrix norm |
||
| 420 | $norm = 0.0; |
||
| 421 | |||
| 422 | for ($i = 0; $i < $nn; ++$i) { |
||
| 423 | if (($i < $low) or ($i > $high)) { |
||
| 424 | $this->d[$i] = $this->H[$i][$i]; |
||
| 425 | $this->e[$i] = 0.0; |
||
| 426 | } |
||
| 427 | for ($j = max($i-1, 0); $j < $nn; ++$j) { |
||
| 428 | $norm = $norm + abs($this->H[$i][$j]); |
||
| 429 | } |
||
| 430 | } |
||
| 431 | |||
| 432 | // Outer loop over eigenvalue index |
||
| 433 | $iter = 0; |
||
| 434 | while ($n >= $low) { |
||
| 435 | // Look for single small sub-diagonal element |
||
| 436 | $l = $n; |
||
| 437 | while ($l > $low) { |
||
| 438 | $s = abs($this->H[$l-1][$l-1]) + abs($this->H[$l][$l]); |
||
| 439 | if ($s == 0.0) { |
||
| 440 | $s = $norm; |
||
| 441 | } |
||
| 442 | if (abs($this->H[$l][$l-1]) < $eps * $s) { |
||
| 443 | break; |
||
| 444 | } |
||
| 445 | --$l; |
||
| 446 | } |
||
| 447 | // Check for convergence |
||
| 448 | // One root found |
||
| 449 | if ($l == $n) { |
||
| 450 | $this->H[$n][$n] = $this->H[$n][$n] + $exshift; |
||
| 451 | $this->d[$n] = $this->H[$n][$n]; |
||
| 452 | $this->e[$n] = 0.0; |
||
| 453 | --$n; |
||
| 454 | $iter = 0; |
||
| 455 | // Two roots found |
||
| 456 | } elseif ($l == $n-1) { |
||
| 457 | $w = $this->H[$n][$n-1] * $this->H[$n-1][$n]; |
||
| 458 | $p = ($this->H[$n-1][$n-1] - $this->H[$n][$n]) / 2.0; |
||
| 459 | $q = $p * $p + $w; |
||
| 460 | $z = sqrt(abs($q)); |
||
| 461 | $this->H[$n][$n] = $this->H[$n][$n] + $exshift; |
||
| 462 | $this->H[$n-1][$n-1] = $this->H[$n-1][$n-1] + $exshift; |
||
| 463 | $x = $this->H[$n][$n]; |
||
| 464 | // Real pair |
||
| 465 | if ($q >= 0) { |
||
| 466 | if ($p >= 0) { |
||
| 467 | $z = $p + $z; |
||
| 468 | } else { |
||
| 469 | $z = $p - $z; |
||
| 470 | } |
||
| 471 | $this->d[$n-1] = $x + $z; |
||
| 472 | $this->d[$n] = $this->d[$n-1]; |
||
| 473 | if ($z != 0.0) { |
||
| 474 | $this->d[$n] = $x - $w / $z; |
||
| 475 | } |
||
| 476 | $this->e[$n-1] = 0.0; |
||
| 477 | $this->e[$n] = 0.0; |
||
| 478 | $x = $this->H[$n][$n-1]; |
||
| 479 | $s = abs($x) + abs($z); |
||
| 480 | $p = $x / $s; |
||
| 481 | $q = $z / $s; |
||
| 482 | $r = sqrt($p * $p + $q * $q); |
||
| 483 | $p = $p / $r; |
||
| 484 | $q = $q / $r; |
||
| 485 | // Row modification |
||
| 486 | View Code Duplication | for ($j = $n-1; $j < $nn; ++$j) { |
|
| 487 | $z = $this->H[$n-1][$j]; |
||
| 488 | $this->H[$n-1][$j] = $q * $z + $p * $this->H[$n][$j]; |
||
| 489 | $this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z; |
||
| 490 | } |
||
| 491 | // Column modification |
||
| 492 | View Code Duplication | for ($i = 0; $i <= $n; ++$i) { |
|
| 493 | $z = $this->H[$i][$n-1]; |
||
| 494 | $this->H[$i][$n-1] = $q * $z + $p * $this->H[$i][$n]; |
||
| 495 | $this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z; |
||
| 496 | } |
||
| 497 | // Accumulate transformations |
||
| 498 | View Code Duplication | for ($i = $low; $i <= $high; ++$i) { |
|
| 499 | $z = $this->V[$i][$n-1]; |
||
| 500 | $this->V[$i][$n-1] = $q * $z + $p * $this->V[$i][$n]; |
||
| 501 | $this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z; |
||
| 502 | } |
||
| 503 | // Complex pair |
||
| 504 | } else { |
||
| 505 | $this->d[$n-1] = $x + $p; |
||
| 506 | $this->d[$n] = $x + $p; |
||
| 507 | $this->e[$n-1] = $z; |
||
| 508 | $this->e[$n] = -$z; |
||
| 509 | } |
||
| 510 | $n = $n - 2; |
||
| 511 | $iter = 0; |
||
| 512 | // No convergence yet |
||
| 513 | } else { |
||
| 514 | // Form shift |
||
| 515 | $x = $this->H[$n][$n]; |
||
| 516 | $y = 0.0; |
||
| 517 | $w = 0.0; |
||
| 518 | if ($l < $n) { |
||
| 519 | $y = $this->H[$n-1][$n-1]; |
||
| 520 | $w = $this->H[$n][$n-1] * $this->H[$n-1][$n]; |
||
| 521 | } |
||
| 522 | // Wilkinson's original ad hoc shift |
||
| 523 | if ($iter == 10) { |
||
| 524 | $exshift += $x; |
||
| 525 | for ($i = $low; $i <= $n; ++$i) { |
||
| 526 | $this->H[$i][$i] -= $x; |
||
| 527 | } |
||
| 528 | $s = abs($this->H[$n][$n-1]) + abs($this->H[$n-1][$n-2]); |
||
| 529 | $x = $y = 0.75 * $s; |
||
| 530 | $w = -0.4375 * $s * $s; |
||
| 531 | } |
||
| 532 | // MATLAB's new ad hoc shift |
||
| 533 | if ($iter == 30) { |
||
| 534 | $s = ($y - $x) / 2.0; |
||
| 535 | $s = $s * $s + $w; |
||
| 536 | if ($s > 0) { |
||
| 537 | $s = sqrt($s); |
||
| 538 | if ($y < $x) { |
||
| 539 | $s = -$s; |
||
| 540 | } |
||
| 541 | $s = $x - $w / (($y - $x) / 2.0 + $s); |
||
| 542 | for ($i = $low; $i <= $n; ++$i) { |
||
| 543 | $this->H[$i][$i] -= $s; |
||
| 544 | } |
||
| 545 | $exshift += $s; |
||
| 546 | $x = $y = $w = 0.964; |
||
| 547 | } |
||
| 548 | } |
||
| 549 | // Could check iteration count here. |
||
| 550 | $iter = $iter + 1; |
||
| 551 | // Look for two consecutive small sub-diagonal elements |
||
| 552 | $m = $n - 2; |
||
| 553 | while ($m >= $l) { |
||
| 554 | $z = $this->H[$m][$m]; |
||
| 555 | $r = $x - $z; |
||
| 556 | $s = $y - $z; |
||
| 557 | $p = ($r * $s - $w) / $this->H[$m+1][$m] + $this->H[$m][$m+1]; |
||
| 558 | $q = $this->H[$m+1][$m+1] - $z - $r - $s; |
||
| 559 | $r = $this->H[$m+2][$m+1]; |
||
| 560 | $s = abs($p) + abs($q) + abs($r); |
||
| 561 | $p = $p / $s; |
||
| 562 | $q = $q / $s; |
||
| 563 | $r = $r / $s; |
||
| 564 | if ($m == $l) { |
||
| 565 | break; |
||
| 566 | } |
||
| 567 | if (abs($this->H[$m][$m-1]) * (abs($q) + abs($r)) < |
||
| 568 | $eps * (abs($p) * (abs($this->H[$m-1][$m-1]) + abs($z) + abs($this->H[$m+1][$m+1])))) { |
||
| 569 | break; |
||
| 570 | } |
||
| 571 | --$m; |
||
| 572 | } |
||
| 573 | for ($i = $m + 2; $i <= $n; ++$i) { |
||
| 574 | $this->H[$i][$i-2] = 0.0; |
||
| 575 | if ($i > $m+2) { |
||
| 576 | $this->H[$i][$i-3] = 0.0; |
||
| 577 | } |
||
| 578 | } |
||
| 579 | // Double QR step involving rows l:n and columns m:n |
||
| 580 | for ($k = $m; $k <= $n-1; ++$k) { |
||
| 581 | $notlast = ($k != $n-1); |
||
| 582 | if ($k != $m) { |
||
| 583 | $p = $this->H[$k][$k-1]; |
||
| 584 | $q = $this->H[$k+1][$k-1]; |
||
| 585 | $r = ($notlast ? $this->H[$k+2][$k-1] : 0.0); |
||
| 586 | $x = abs($p) + abs($q) + abs($r); |
||
| 587 | if ($x != 0.0) { |
||
| 588 | $p = $p / $x; |
||
| 589 | $q = $q / $x; |
||
| 590 | $r = $r / $x; |
||
| 591 | } |
||
| 592 | } |
||
| 593 | if ($x == 0.0) { |
||
| 594 | break; |
||
| 595 | } |
||
| 596 | $s = sqrt($p * $p + $q * $q + $r * $r); |
||
| 597 | if ($p < 0) { |
||
| 598 | $s = -$s; |
||
| 599 | } |
||
| 600 | if ($s != 0) { |
||
| 601 | if ($k != $m) { |
||
| 602 | $this->H[$k][$k-1] = -$s * $x; |
||
| 603 | } elseif ($l != $m) { |
||
| 604 | $this->H[$k][$k-1] = -$this->H[$k][$k-1]; |
||
| 605 | } |
||
| 606 | $p = $p + $s; |
||
| 607 | $x = $p / $s; |
||
| 608 | $y = $q / $s; |
||
| 609 | $z = $r / $s; |
||
| 610 | $q = $q / $p; |
||
| 611 | $r = $r / $p; |
||
| 612 | // Row modification |
||
| 613 | View Code Duplication | for ($j = $k; $j < $nn; ++$j) { |
|
| 614 | $p = $this->H[$k][$j] + $q * $this->H[$k+1][$j]; |
||
| 615 | if ($notlast) { |
||
| 616 | $p = $p + $r * $this->H[$k+2][$j]; |
||
| 617 | $this->H[$k+2][$j] = $this->H[$k+2][$j] - $p * $z; |
||
| 618 | } |
||
| 619 | $this->H[$k][$j] = $this->H[$k][$j] - $p * $x; |
||
| 620 | $this->H[$k+1][$j] = $this->H[$k+1][$j] - $p * $y; |
||
| 621 | } |
||
| 622 | // Column modification |
||
| 623 | View Code Duplication | for ($i = 0; $i <= min($n, $k+3); ++$i) { |
|
| 624 | $p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k+1]; |
||
| 625 | if ($notlast) { |
||
| 626 | $p = $p + $z * $this->H[$i][$k+2]; |
||
| 627 | $this->H[$i][$k+2] = $this->H[$i][$k+2] - $p * $r; |
||
| 628 | } |
||
| 629 | $this->H[$i][$k] = $this->H[$i][$k] - $p; |
||
| 630 | $this->H[$i][$k+1] = $this->H[$i][$k+1] - $p * $q; |
||
| 631 | } |
||
| 632 | // Accumulate transformations |
||
| 633 | View Code Duplication | for ($i = $low; $i <= $high; ++$i) { |
|
| 634 | $p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k+1]; |
||
| 635 | if ($notlast) { |
||
| 636 | $p = $p + $z * $this->V[$i][$k+2]; |
||
| 637 | $this->V[$i][$k+2] = $this->V[$i][$k+2] - $p * $r; |
||
| 638 | } |
||
| 639 | $this->V[$i][$k] = $this->V[$i][$k] - $p; |
||
| 640 | $this->V[$i][$k+1] = $this->V[$i][$k+1] - $p * $q; |
||
| 641 | } |
||
| 642 | } // ($s != 0) |
||
| 643 | } // k loop |
||
| 644 | } // check convergence |
||
| 645 | } // while ($n >= $low) |
||
| 646 | |||
| 647 | // Backsubstitute to find vectors of upper triangular form |
||
| 648 | if ($norm == 0.0) { |
||
| 649 | return; |
||
| 650 | } |
||
| 651 | |||
| 652 | for ($n = $nn-1; $n >= 0; --$n) { |
||
| 653 | $p = $this->d[$n]; |
||
| 654 | $q = $this->e[$n]; |
||
| 655 | // Real vector |
||
| 656 | if ($q == 0) { |
||
| 657 | $l = $n; |
||
| 658 | $this->H[$n][$n] = 1.0; |
||
| 659 | for ($i = $n-1; $i >= 0; --$i) { |
||
| 660 | $w = $this->H[$i][$i] - $p; |
||
| 661 | $r = 0.0; |
||
| 662 | for ($j = $l; $j <= $n; ++$j) { |
||
| 663 | $r = $r + $this->H[$i][$j] * $this->H[$j][$n]; |
||
| 664 | } |
||
| 665 | if ($this->e[$i] < 0.0) { |
||
| 666 | $z = $w; |
||
| 667 | $s = $r; |
||
| 668 | } else { |
||
| 669 | $l = $i; |
||
| 670 | if ($this->e[$i] == 0.0) { |
||
| 671 | if ($w != 0.0) { |
||
| 672 | $this->H[$i][$n] = -$r / $w; |
||
| 673 | } else { |
||
| 674 | $this->H[$i][$n] = -$r / ($eps * $norm); |
||
| 675 | } |
||
| 676 | // Solve real equations |
||
| 677 | } else { |
||
| 678 | $x = $this->H[$i][$i+1]; |
||
| 679 | $y = $this->H[$i+1][$i]; |
||
| 680 | $q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i]; |
||
| 681 | $t = ($x * $s - $z * $r) / $q; |
||
| 682 | $this->H[$i][$n] = $t; |
||
| 683 | if (abs($x) > abs($z)) { |
||
| 684 | $this->H[$i+1][$n] = (-$r - $w * $t) / $x; |
||
| 685 | } else { |
||
| 686 | $this->H[$i+1][$n] = (-$s - $y * $t) / $z; |
||
| 687 | } |
||
| 688 | } |
||
| 689 | // Overflow control |
||
| 690 | $t = abs($this->H[$i][$n]); |
||
| 691 | if (($eps * $t) * $t > 1) { |
||
| 692 | View Code Duplication | for ($j = $i; $j <= $n; ++$j) { |
|
| 693 | $this->H[$j][$n] = $this->H[$j][$n] / $t; |
||
| 694 | } |
||
| 695 | } |
||
| 696 | } |
||
| 697 | } |
||
| 698 | // Complex vector |
||
| 699 | } elseif ($q < 0) { |
||
| 700 | $l = $n-1; |
||
| 701 | // Last vector component imaginary so matrix is triangular |
||
| 702 | if (abs($this->H[$n][$n-1]) > abs($this->H[$n-1][$n])) { |
||
| 703 | $this->H[$n-1][$n-1] = $q / $this->H[$n][$n-1]; |
||
| 704 | $this->H[$n-1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n-1]; |
||
| 705 | } else { |
||
| 706 | $this->cdiv(0.0, -$this->H[$n-1][$n], $this->H[$n-1][$n-1] - $p, $q); |
||
| 707 | $this->H[$n-1][$n-1] = $this->cdivr; |
||
| 708 | $this->H[$n-1][$n] = $this->cdivi; |
||
| 709 | } |
||
| 710 | $this->H[$n][$n-1] = 0.0; |
||
| 711 | $this->H[$n][$n] = 1.0; |
||
| 712 | for ($i = $n-2; $i >= 0; --$i) { |
||
| 713 | // double ra,sa,vr,vi; |
||
| 714 | $ra = 0.0; |
||
| 715 | $sa = 0.0; |
||
| 716 | for ($j = $l; $j <= $n; ++$j) { |
||
| 717 | $ra = $ra + $this->H[$i][$j] * $this->H[$j][$n-1]; |
||
| 718 | $sa = $sa + $this->H[$i][$j] * $this->H[$j][$n]; |
||
| 719 | } |
||
| 720 | $w = $this->H[$i][$i] - $p; |
||
| 721 | if ($this->e[$i] < 0.0) { |
||
| 722 | $z = $w; |
||
| 723 | $r = $ra; |
||
| 724 | $s = $sa; |
||
| 725 | } else { |
||
| 726 | $l = $i; |
||
| 727 | if ($this->e[$i] == 0) { |
||
| 728 | $this->cdiv(-$ra, -$sa, $w, $q); |
||
| 729 | $this->H[$i][$n-1] = $this->cdivr; |
||
| 730 | $this->H[$i][$n] = $this->cdivi; |
||
| 731 | } else { |
||
| 732 | // Solve complex equations |
||
| 733 | $x = $this->H[$i][$i+1]; |
||
| 734 | $y = $this->H[$i+1][$i]; |
||
| 735 | $vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q; |
||
| 736 | $vi = ($this->d[$i] - $p) * 2.0 * $q; |
||
| 737 | if ($vr == 0.0 & $vi == 0.0) { |
||
| 738 | $vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z)); |
||
| 739 | } |
||
| 740 | $this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi); |
||
| 741 | $this->H[$i][$n-1] = $this->cdivr; |
||
| 742 | $this->H[$i][$n] = $this->cdivi; |
||
| 743 | if (abs($x) > (abs($z) + abs($q))) { |
||
| 744 | $this->H[$i+1][$n-1] = (-$ra - $w * $this->H[$i][$n-1] + $q * $this->H[$i][$n]) / $x; |
||
| 745 | $this->H[$i+1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n-1]) / $x; |
||
| 746 | } else { |
||
| 747 | $this->cdiv(-$r - $y * $this->H[$i][$n-1], -$s - $y * $this->H[$i][$n], $z, $q); |
||
| 748 | $this->H[$i+1][$n-1] = $this->cdivr; |
||
| 749 | $this->H[$i+1][$n] = $this->cdivi; |
||
| 750 | } |
||
| 751 | } |
||
| 752 | // Overflow control |
||
| 753 | $t = max(abs($this->H[$i][$n-1]), abs($this->H[$i][$n])); |
||
| 754 | if (($eps * $t) * $t > 1) { |
||
| 755 | for ($j = $i; $j <= $n; ++$j) { |
||
| 756 | $this->H[$j][$n-1] = $this->H[$j][$n-1] / $t; |
||
| 757 | $this->H[$j][$n] = $this->H[$j][$n] / $t; |
||
| 758 | } |
||
| 759 | } |
||
| 760 | } // end else |
||
| 761 | } // end for |
||
| 762 | } // end else for complex case |
||
| 763 | } // end for |
||
| 764 | |||
| 765 | // Vectors of isolated roots |
||
| 766 | for ($i = 0; $i < $nn; ++$i) { |
||
| 767 | if ($i < $low | $i > $high) { |
||
| 768 | for ($j = $i; $j < $nn; ++$j) { |
||
| 769 | $this->V[$i][$j] = $this->H[$i][$j]; |
||
| 770 | } |
||
| 771 | } |
||
| 772 | } |
||
| 773 | |||
| 774 | // Back transformation to get eigenvectors of original matrix |
||
| 775 | for ($j = $nn-1; $j >= $low; --$j) { |
||
| 776 | for ($i = $low; $i <= $high; ++$i) { |
||
| 777 | $z = 0.0; |
||
| 778 | for ($k = $low; $k <= min($j, $high); ++$k) { |
||
| 779 | $z = $z + $this->V[$i][$k] * $this->H[$k][$j]; |
||
| 780 | } |
||
| 781 | $this->V[$i][$j] = $z; |
||
| 782 | } |
||
| 783 | } |
||
| 784 | } // end hqr2 |
||
| 785 | |||
| 891 |
Loading history...
Duplicated code is one of the most pungent code smells. If you need to duplicate the same code in three or more different places, we strongly encourage you to look into extracting the code into a single class or operation.
You can also find more detailed suggestions in the “Code” section of your repository.