Subversion Repositories eFlore/Applications.cel

<?php
/**
 *      @package JAMA
 *
 *      For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
 *      orthogonal matrix Q and an n-by-n upper triangular matrix R so that
 *      A = Q*R.
 *
 *      The QR decompostion always exists, even if the matrix does not have
 *      full rank, so the constructor will never fail.  The primary use of the
 *      QR decomposition is in the least squares solution of nonsquare systems
 *      of simultaneous linear equations.  This will fail if isFullRank()
 *      returns false.
 *
 *      @author  Paul Meagher
 *      @license PHP v3.0
 *      @version 1.1
 */
class PHPExcel_Shared_JAMA_QRDecomposition {

        const MatrixRankException       = "Can only perform operation on full-rank matrix.";

        /**
         *      Array for internal storage of decomposition.
         *      @var array
         */
        private $QR = array();

        /**
         *      Row dimension.
         *      @var integer
         */
        private $m;

        /**
        *       Column dimension.
        *       @var integer
        */
        private $n;

        /**
         *      Array for internal storage of diagonal of R.
         *      @var  array
         */
        private $Rdiag = array();


        /**
         *      QR Decomposition computed by Householder reflections.
         *
         *      @param matrix $A Rectangular matrix
         *      @return Structure to access R and the Householder vectors and compute Q.
         */
        public function __construct($A) {
                if($A instanceof PHPExcel_Shared_JAMA_Matrix) {
                        // Initialize.
                        $this->QR = $A->getArrayCopy();
                        $this->m  = $A->getRowDimension();
                        $this->n  = $A->getColumnDimension();
                        // Main loop.
                        for ($k = 0; $k < $this->n; ++$k) {
                                // Compute 2-norm of k-th column without under/overflow.
                                $nrm = 0.0;
                                for ($i = $k; $i < $this->m; ++$i) {
                                        $nrm = hypo($nrm, $this->QR[$i][$k]);
                                }
                                if ($nrm != 0.0) {
                                        // Form k-th Householder vector.
                                        if ($this->QR[$k][$k] < 0) {
                                                $nrm = -$nrm;
                                        }
                                        for ($i = $k; $i < $this->m; ++$i) {
                                                $this->QR[$i][$k] /= $nrm;
                                        }
                                        $this->QR[$k][$k] += 1.0;
                                        // Apply transformation to remaining columns.
                                        for ($j = $k+1; $j < $this->n; ++$j) {
                                                $s = 0.0;
                                                for ($i = $k; $i < $this->m; ++$i) {
                                                        $s += $this->QR[$i][$k] * $this->QR[$i][$j];
                                                }
                                                $s = -$s/$this->QR[$k][$k];
                                                for ($i = $k; $i < $this->m; ++$i) {
                                                        $this->QR[$i][$j] += $s * $this->QR[$i][$k];
                                                }
                                        }
                                }
                                $this->Rdiag[$k] = -$nrm;
                        }
                } else {
                        throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::ArgumentTypeException);
                }
        }       //      function __construct()


        /**
         *      Is the matrix full rank?
         *
         *      @return boolean true if R, and hence A, has full rank, else false.
         */
        public function isFullRank() {
                for ($j = 0; $j < $this->n; ++$j) {
                        if ($this->Rdiag[$j] == 0) {
                                return false;
                        }
                }
                return true;
        }       //      function isFullRank()


        /**
         *      Return the Householder vectors
         *
         *      @return Matrix Lower trapezoidal matrix whose columns define the reflections
         */
        public function getH() {
                for ($i = 0; $i < $this->m; ++$i) {
                        for ($j = 0; $j < $this->n; ++$j) {
                                if ($i >= $j) {
                                        $H[$i][$j] = $this->QR[$i][$j];
                                } else {
                                        $H[$i][$j] = 0.0;
                                }
                        }
                }
                return new PHPExcel_Shared_JAMA_Matrix($H);
        }       //      function getH()


        /**
         *      Return the upper triangular factor
         *
         *      @return Matrix upper triangular factor
         */
        public function getR() {
                for ($i = 0; $i < $this->n; ++$i) {
                        for ($j = 0; $j < $this->n; ++$j) {
                                if ($i < $j) {
                                        $R[$i][$j] = $this->QR[$i][$j];
                                } elseif ($i == $j) {
                                        $R[$i][$j] = $this->Rdiag[$i];
                                } else {
                                        $R[$i][$j] = 0.0;
                                }
                        }
                }
                return new PHPExcel_Shared_JAMA_Matrix($R);
        }       //      function getR()


        /**
         *      Generate and return the (economy-sized) orthogonal factor
         *
         *      @return Matrix orthogonal factor
         */
        public function getQ() {
                for ($k = $this->n-1; $k >= 0; --$k) {
                        for ($i = 0; $i < $this->m; ++$i) {
                                $Q[$i][$k] = 0.0;
                        }
                        $Q[$k][$k] = 1.0;
                        for ($j = $k; $j < $this->n; ++$j) {
                                if ($this->QR[$k][$k] != 0) {
                                        $s = 0.0;
                                        for ($i = $k; $i < $this->m; ++$i) {
                                                $s += $this->QR[$i][$k] * $Q[$i][$j];
                                        }
                                        $s = -$s/$this->QR[$k][$k];
                                        for ($i = $k; $i < $this->m; ++$i) {
                                                $Q[$i][$j] += $s * $this->QR[$i][$k];
                                        }
                                }
                        }
                }
                /*
                for($i = 0; $i < count($Q); ++$i) {
                        for($j = 0; $j < count($Q); ++$j) {
                                if(! isset($Q[$i][$j]) ) {
                                        $Q[$i][$j] = 0;
                                }
                        }
                }
                */
                return new PHPExcel_Shared_JAMA_Matrix($Q);
        }       //      function getQ()


        /**
         *      Least squares solution of A*X = B
         *
         *      @param Matrix $B A Matrix with as many rows as A and any number of columns.
         *      @return Matrix Matrix that minimizes the two norm of Q*R*X-B.
         */
        public function solve($B) {
                if ($B->getRowDimension() == $this->m) {
                        if ($this->isFullRank()) {
                                // Copy right hand side
                                $nx = $B->getColumnDimension();
                                $X  = $B->getArrayCopy();
                                // Compute Y = transpose(Q)*B
                                for ($k = 0; $k < $this->n; ++$k) {
                                        for ($j = 0; $j < $nx; ++$j) {
                                                $s = 0.0;
                                                for ($i = $k; $i < $this->m; ++$i) {
                                                        $s += $this->QR[$i][$k] * $X[$i][$j];
                                                }
                                                $s = -$s/$this->QR[$k][$k];
                                                for ($i = $k; $i < $this->m; ++$i) {
                                                        $X[$i][$j] += $s * $this->QR[$i][$k];
                                                }
                                        }
                                }
                                // Solve R*X = Y;
                                for ($k = $this->n-1; $k >= 0; --$k) {
                                        for ($j = 0; $j < $nx; ++$j) {
                                                $X[$k][$j] /= $this->Rdiag[$k];
                                        }
                                        for ($i = 0; $i < $k; ++$i) {
                                                for ($j = 0; $j < $nx; ++$j) {
                                                        $X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
                                                }
                                        }
                                }
                                $X = new PHPExcel_Shared_JAMA_Matrix($X);
                                return ($X->getMatrix(0, $this->n-1, 0, $nx));
                        } else {
                                throw new PHPExcel_Calculation_Exception(self::MatrixRankException);
                        }
                } else {
                        throw new PHPExcel_Calculation_Exception(PHPExcel_Shared_JAMA_Matrix::MatrixDimensionException);
                }
        }       //      function solve()

}       //      PHPExcel_Shared_JAMA_class QRDecomposition