<?php
/**
 *    @package JAMA
 *
 *    For an m-by-n matrix A with m >= n, the singular value decomposition is
 *    an m-by-n orthogonal matrix U, an n-by-n diagonal matrix S, and
 *    an n-by-n orthogonal matrix V so that A = U*S*V'.
 *
 *    The singular values, sigma[$k] = S[$k][$k], are ordered so that
 *    sigma[0] >= sigma[1] >= ... >= sigma[n-1].
 *
 *    The singular value decompostion always exists, so the constructor will
 *    never fail.  The matrix condition number and the effective numerical
 *    rank can be computed from this decomposition.
 *
 *    @author  Paul Meagher
 *    @license PHP v3.0
 *    @version 1.1
 */
class SingularValueDecomposition
{
  /**
   * Internal storage of U.
   * 
   * @var array
   */
  private $U = array ();
  
  /**
   * Internal storage of V.
   * 
   * @var array
   */
  private $V = array ();
  
  /**
   * Internal storage of singular values.
   * 
   * @var array
   */
  private $s = array ();
  
  /**
   * Row dimension.
   * 
   * @var int
   */
  private $m;
  
  /**
   * Column dimension.
   * 
   * @var int
   */
  private $n;
  
  /**
   * Construct the singular value decomposition
   *
   * Derived from LINPACK code.
   *
   * @param $A Rectangular
   *          matrix
   * @return Structure to access U, S and V.
   */
  public function __construct($Arg)
  {
    // Initialize.
    $A = $Arg->getArrayCopy ();
    $this->m = $Arg->getRowDimension ();
    $this->n = $Arg->getColumnDimension ();
    $nu = min ( $this->m, $this->n );
    $e = array ();
    $work = array ();
    $wantu = true;
    $wantv = true;
    $nct = min ( $this->m - 1, $this->n );
    $nrt = max ( 0, min ( $this->n - 2, $this->m ) );
    
    // Reduce A to bidiagonal form, storing the diagonal elements
    // in s and the super-diagonal elements in e.
    for($k = 0; $k < max ( $nct, $nrt ); ++ $k)
    {
      if ($k < $nct)
      {
        // Compute the transformation for the k-th column and
        // place the k-th diagonal in s[$k].
        // Compute 2-norm of k-th column without under/overflow.
        $this->s [$k] = 0;
        for($i = $k; $i < $this->m; ++ $i)
        {
          $this->s [$k] = hypo ( $this->s [$k], $A [$i] [$k] );
        }
        if ($this->s [$k] != 0.0)
        {
          if ($A [$k] [$k] < 0.0)
          {
            $this->s [$k] = - $this->s [$k];
          }
          for($i = $k; $i < $this->m; ++ $i)
          {
            $A [$i] [$k] /= $this->s [$k];
          }
          $A [$k] [$k] += 1.0;
        }
        $this->s [$k] = - $this->s [$k];
      }
      
      for($j = $k + 1; $j < $this->n; ++ $j)
      {
        if (($k < $nct) & ($this->s [$k] != 0.0))
        {
          // Apply the transformation.
          $t = 0;
          for($i = $k; $i < $this->m; ++ $i)
          {
            $t += $A [$i] [$k] * $A [$i] [$j];
          }
          $t = - $t / $A [$k] [$k];
          for($i = $k; $i < $this->m; ++ $i)
          {
            $A [$i] [$j] += $t * $A [$i] [$k];
          }
          // Place the k-th row of A into e for the
          // subsequent calculation of the row transformation.
          $e [$j] = $A [$k] [$j];
        }
      }
      
      if ($wantu and ($k < $nct))
      {
        // Place the transformation in U for subsequent back
        // multiplication.
        for($i = $k; $i < $this->m; ++ $i)
        {
          $this->U [$i] [$k] = $A [$i] [$k];
        }
      }
      
      if ($k < $nrt)
      {
        // Compute the k-th row transformation and place the
        // k-th super-diagonal in e[$k].
        // Compute 2-norm without under/overflow.
        $e [$k] = 0;
        for($i = $k + 1; $i < $this->n; ++ $i)
        {
          $e [$k] = hypo ( $e [$k], $e [$i] );
        }
        if ($e [$k] != 0.0)
        {
          if ($e [$k + 1] < 0.0)
          {
            $e [$k] = - $e [$k];
          }
          for($i = $k + 1; $i < $this->n; ++ $i)
          {
            $e [$i] /= $e [$k];
          }
          $e [$k + 1] += 1.0;
        }
        $e [$k] = - $e [$k];
        if (($k + 1 < $this->m) and ($e [$k] != 0.0))
        {
          // Apply the transformation.
          for($i = $k + 1; $i < $this->m; ++ $i)
          {
            $work [$i] = 0.0;
          }
          for($j = $k + 1; $j < $this->n; ++ $j)
          {
            for($i = $k + 1; $i < $this->m; ++ $i)
            {
              $work [$i] += $e [$j] * $A [$i] [$j];
            }
          }
          for($j = $k + 1; $j < $this->n; ++ $j)
          {
            $t = - $e [$j] / $e [$k + 1];
            for($i = $k + 1; $i < $this->m; ++ $i)
            {
              $A [$i] [$j] += $t * $work [$i];
            }
          }
        }
        if ($wantv)
        {
          // Place the transformation in V for subsequent
          // back multiplication.
          for($i = $k + 1; $i < $this->n; ++ $i)
          {
            $this->V [$i] [$k] = $e [$i];
          }
        }
      }
    }
    
    // Set up the final bidiagonal matrix or order p.
    $p = min ( $this->n, $this->m + 1 );
    if ($nct < $this->n)
    {
      $this->s [$nct] = $A [$nct] [$nct];
    }
    if ($this->m < $p)
    {
      $this->s [$p - 1] = 0.0;
    }
    if ($nrt + 1 < $p)
    {
      $e [$nrt] = $A [$nrt] [$p - 1];
    }
    $e [$p - 1] = 0.0;
    // If required, generate U.
    if ($wantu)
    {
      for($j = $nct; $j < $nu; ++ $j)
      {
        for($i = 0; $i < $this->m; ++ $i)
        {
          $this->U [$i] [$j] = 0.0;
        }
        $this->U [$j] [$j] = 1.0;
      }
      for($k = $nct - 1; $k >= 0; -- $k)
      {
        if ($this->s [$k] != 0.0)
        {
          for($j = $k + 1; $j < $nu; ++ $j)
          {
            $t = 0;
            for($i = $k; $i < $this->m; ++ $i)
            {
              $t += $this->U [$i] [$k] * $this->U [$i] [$j];
            }
            $t = - $t / $this->U [$k] [$k];
            for($i = $k; $i < $this->m; ++ $i)
            {
              $this->U [$i] [$j] += $t * $this->U [$i] [$k];
            }
          }
          for($i = $k; $i < $this->m; ++ $i)
          {
            $this->U [$i] [$k] = - $this->U [$i] [$k];
          }
          $this->U [$k] [$k] = 1.0 + $this->U [$k] [$k];
          for($i = 0; $i < $k - 1; ++ $i)
          {
            $this->U [$i] [$k] = 0.0;
          }
        }
        else
        {
          for($i = 0; $i < $this->m; ++ $i)
          {
            $this->U [$i] [$k] = 0.0;
          }
          $this->U [$k] [$k] = 1.0;
        }
      }
    }
    
    // If required, generate V.
    if ($wantv)
    {
      for($k = $this->n - 1; $k >= 0; -- $k)
      {
        if (($k < $nrt) and ($e [$k] != 0.0))
        {
          for($j = $k + 1; $j < $nu; ++ $j)
          {
            $t = 0;
            for($i = $k + 1; $i < $this->n; ++ $i)
            {
              $t += $this->V [$i] [$k] * $this->V [$i] [$j];
            }
            $t = - $t / $this->V [$k + 1] [$k];
            for($i = $k + 1; $i < $this->n; ++ $i)
            {
              $this->V [$i] [$j] += $t * $this->V [$i] [$k];
            }
          }
        }
        for($i = 0; $i < $this->n; ++ $i)
        {
          $this->V [$i] [$k] = 0.0;
        }
        $this->V [$k] [$k] = 1.0;
      }
    }
    
    // Main iteration loop for the singular values.
    $pp = $p - 1;
    $iter = 0;
    $eps = pow ( 2.0, - 52.0 );
    
    while ( $p > 0 )
    {
      // Here is where a test for too many iterations would go.
      // This section of the program inspects for negligible
      // elements in the s and e arrays. On completion the
      // variables kase and k are set as follows:
      // kase = 1 if s(p) and e[k-1] are negligible and k<p
      // kase = 2 if s(k) is negligible and k<p
      // kase = 3 if e[k-1] is negligible, k<p, and
      // s(k), ..., s(p) are not negligible (qr step).
      // kase = 4 if e(p-1) is negligible (convergence).
      for($k = $p - 2; $k >= - 1; -- $k)
      {
        if ($k == - 1)
        {
          break;
        }
        if (abs ( $e [$k] ) <= $eps * (abs ( $this->s [$k] ) + abs ( $this->s [$k + 1] )))
        {
          $e [$k] = 0.0;
          break;
        }
      }
      if ($k == $p - 2)
      {
        $kase = 4;
      }
      else
      {
        for($ks = $p - 1; $ks >= $k; -- $ks)
        {
          if ($ks == $k)
          {
            break;
          }
          $t = ($ks != $p ? abs ( $e [$ks] ) : 0.) + ($ks != $k + 1 ? abs ( $e [$ks - 1] ) : 0.);
          if (abs ( $this->s [$ks] ) <= $eps * $t)
          {
            $this->s [$ks] = 0.0;
            break;
          }
        }
        if ($ks == $k)
        {
          $kase = 3;
        }
        elseif ($ks == $p - 1)
        {
          $kase = 1;
        }
        else
        {
          $kase = 2;
          $k = $ks;
        }
      }
      ++ $k;
      
      // Perform the task indicated by kase.
      switch ($kase)
      {
        // Deflate negligible s(p).
        case 1 :
          $f = $e [$p - 2];
          $e [$p - 2] = 0.0;
          for($j = $p - 2; $j >= $k; -- $j)
          {
            $t = hypo ( $this->s [$j], $f );
            $cs = $this->s [$j] / $t;
            $sn = $f / $t;
            $this->s [$j] = $t;
            if ($j != $k)
            {
              $f = - $sn * $e [$j - 1];
              $e [$j - 1] = $cs * $e [$j - 1];
            }
            if ($wantv)
            {
              for($i = 0; $i < $this->n; ++ $i)
              {
                $t = $cs * $this->V [$i] [$j] + $sn * $this->V [$i] [$p - 1];
                $this->V [$i] [$p - 1] = - $sn * $this->V [$i] [$j] + $cs * $this->V [$i] [$p - 1];
                $this->V [$i] [$j] = $t;
              }
            }
          }
          break;
        // Split at negligible s(k).
        case 2 :
          $f = $e [$k - 1];
          $e [$k - 1] = 0.0;
          for($j = $k; $j < $p; ++ $j)
          {
            $t = hypo ( $this->s [$j], $f );
            $cs = $this->s [$j] / $t;
            $sn = $f / $t;
            $this->s [$j] = $t;
            $f = - $sn * $e [$j];
            $e [$j] = $cs * $e [$j];
            if ($wantu)
            {
              for($i = 0; $i < $this->m; ++ $i)
              {
                $t = $cs * $this->U [$i] [$j] + $sn * $this->U [$i] [$k - 1];
                $this->U [$i] [$k - 1] = - $sn * $this->U [$i] [$j] + $cs * $this->U [$i] [$k - 1];
                $this->U [$i] [$j] = $t;
              }
            }
          }
          break;
        // Perform one qr step.
        case 3 :
          // Calculate the shift.
          $scale = max ( max ( max ( max ( abs ( $this->s [$p - 1] ), abs ( $this->s [$p - 2] ) ), abs ( $e [$p - 2] ) ), abs ( $this->s [$k] ) ), abs ( $e [$k] ) );
          $sp = $this->s [$p - 1] / $scale;
          $spm1 = $this->s [$p - 2] / $scale;
          $epm1 = $e [$p - 2] / $scale;
          $sk = $this->s [$k] / $scale;
          $ek = $e [$k] / $scale;
          $b = (($spm1 + $sp) * ($spm1 - $sp) + $epm1 * $epm1) / 2.0;
          $c = ($sp * $epm1) * ($sp * $epm1);
          $shift = 0.0;
          if (($b != 0.0) || ($c != 0.0))
          {
            $shift = sqrt ( $b * $b + $c );
            if ($b < 0.0)
            {
              $shift = - $shift;
            }
            $shift = $c / ($b + $shift);
          }
          $f = ($sk + $sp) * ($sk - $sp) + $shift;
          $g = $sk * $ek;
          // Chase zeros.
          for($j = $k; $j < $p - 1; ++ $j)
          {
            $t = hypo ( $f, $g );
            $cs = $f / $t;
            $sn = $g / $t;
            if ($j != $k)
            {
              $e [$j - 1] = $t;
            }
            $f = $cs * $this->s [$j] + $sn * $e [$j];
            $e [$j] = $cs * $e [$j] - $sn * $this->s [$j];
            $g = $sn * $this->s [$j + 1];
            $this->s [$j + 1] = $cs * $this->s [$j + 1];
            if ($wantv)
            {
              for($i = 0; $i < $this->n; ++ $i)
              {
                $t = $cs * $this->V [$i] [$j] + $sn * $this->V [$i] [$j + 1];
                $this->V [$i] [$j + 1] = - $sn * $this->V [$i] [$j] + $cs * $this->V [$i] [$j + 1];
                $this->V [$i] [$j] = $t;
              }
            }
            $t = hypo ( $f, $g );
            $cs = $f / $t;
            $sn = $g / $t;
            $this->s [$j] = $t;
            $f = $cs * $e [$j] + $sn * $this->s [$j + 1];
            $this->s [$j + 1] = - $sn * $e [$j] + $cs * $this->s [$j + 1];
            $g = $sn * $e [$j + 1];
            $e [$j + 1] = $cs * $e [$j + 1];
            if ($wantu && ($j < $this->m - 1))
            {
              for($i = 0; $i < $this->m; ++ $i)
              {
                $t = $cs * $this->U [$i] [$j] + $sn * $this->U [$i] [$j + 1];
                $this->U [$i] [$j + 1] = - $sn * $this->U [$i] [$j] + $cs * $this->U [$i] [$j + 1];
                $this->U [$i] [$j] = $t;
              }
            }
          }
          $e [$p - 2] = $f;
          $iter = $iter + 1;
          break;
        // Convergence.
        case 4 :
          // Make the singular values positive.
          if ($this->s [$k] <= 0.0)
          {
            $this->s [$k] = ($this->s [$k] < 0.0 ? - $this->s [$k] : 0.0);
            if ($wantv)
            {
              for($i = 0; $i <= $pp; ++ $i)
              {
                $this->V [$i] [$k] = - $this->V [$i] [$k];
              }
            }
          }
          // Order the singular values.
          while ( $k < $pp )
          {
            if ($this->s [$k] >= $this->s [$k + 1])
            {
              break;
            }
            $t = $this->s [$k];
            $this->s [$k] = $this->s [$k + 1];
            $this->s [$k + 1] = $t;
            if ($wantv and ($k < $this->n - 1))
            {
              for($i = 0; $i < $this->n; ++ $i)
              {
                $t = $this->V [$i] [$k + 1];
                $this->V [$i] [$k + 1] = $this->V [$i] [$k];
                $this->V [$i] [$k] = $t;
              }
            }
            if ($wantu and ($k < $this->m - 1))
            {
              for($i = 0; $i < $this->m; ++ $i)
              {
                $t = $this->U [$i] [$k + 1];
                $this->U [$i] [$k + 1] = $this->U [$i] [$k];
                $this->U [$i] [$k] = $t;
              }
            }
            ++ $k;
          }
          $iter = 0;
          -- $p;
          break;
      } // end switch
    } // end while
  } // end constructor
  
  /**
   * Return the left singular vectors
   *
   * @access public
   * @return U
   */
  public function getU()
  {
    return new Matrix ( $this->U, $this->m, min ( $this->m + 1, $this->n ) );
  }
  
  /**
   * Return the right singular vectors
   *
   * @access public
   * @return V
   */
  public function getV()
  {
    return new Matrix ( $this->V );
  }
  
  /**
   * Return the one-dimensional array of singular values
   *
   * @access public
   * @return diagonal of S.
   */
  public function getSingularValues()
  {
    return $this->s;
  }
  
  /**
   * Return the diagonal matrix of singular values
   *
   * @access public
   * @return S
   */
  public function getS()
  {
    for($i = 0; $i < $this->n; ++ $i)
    {
      for($j = 0; $j < $this->n; ++ $j)
      {
        $S [$i] [$j] = 0.0;
      }
      $S [$i] [$i] = $this->s [$i];
    }
    return new Matrix ( $S );
  }
  
  /**
   * Two norm
   *
   * @access public
   * @return max(S)
   */
  public function norm2()
  {
    return $this->s [0];
  }
  
  /**
   * Two norm condition number
   *
   * @access public
   * @return max(S)/min(S)
   */
  public function cond()
  {
    return $this->s [0] / $this->s [min ( $this->m, $this->n ) - 1];
  }
  
  /**
   * Effective numerical matrix rank
   *
   * @access public
   * @return Number of nonnegligible singular values.
   */
  public function rank()
  {
    $eps = pow ( 2.0, - 52.0 );
    $tol = max ( $this->m, $this->n ) * $this->s [0] * $eps;
    $r = 0;
    for($i = 0; $i < count ( $this->s ); ++ $i)
    {
      if ($this->s [$i] > $tol)
      {
        ++ $r;
      }
    }
    return $r;
  }
}
