Location: PHPKode > scripts > Matrix new > matrix.php
<?php

/** 
* A class to perform operation on matrices. 
*
* Modified : 19/dec/2009 by Pier-André Bouchard St-Amant pabsta [at] econ.queensu.ca
*
* Changes : 
*	- Introduction of the PLU decomposition of square matrices (through lu() for users, and PLU internally)
*	- Significant improvement in stability and speed for inversion and the computation of the determinant (based on the above).
*	- The "matrix" constructor can now take one or three arguments. The two additionnal arguments are the number of rows and the number of columns.
*	- Bug correction in many functions : dimensions where incorrect after many routines if a matrix had a full column/row of zeros.
*	- Corrections in "get_columns"
*	- See the file "changes.php" for the whole history of changes.
*
*  - Comments :
*		- All mathematical functions are written to ressemble paper algebra : 
*				- AxB ==> $A->times($B);
*				- A-B ==> $A->minus($B);
*				- A^(-1) ==> $A->inv();
*				- A' ==> $A->prime();
*				- det(A) ==> $A->det($A);
*				- etc.
*		- This restricts nested calls, but augment readability.
*		- The algorithm for matrix inversion is based on PLU decomposition. The algorithm performs relatively well numerically and requires an order of 
* 			2*n^2 xsizeof(doubles) of memory space.
*		- Error messages are echoed to the standard output (likely to be your webpage)
*/
class matrix 
{
	//Global vars
	/**
	*	The array containing the numbers of the matrix. 
	*	$numbers[$i][$j] is a two dimensional array where $i represents the depth and $j the width. Hence, $numbers[$i] returns one row of the matrix.
	*	As a convention, zeros are not stored in the array. Thus, the dimension of the array is not necessarely the dimension of the matrix.
	*/
	var $numbers = array();
	
	/**
	* The number of columns in this matrix.
	*/
	var $numColumns = 0;
	
	/**
	* The number of rows in this matrix.
	*/
	var $numRows = 0;
	
	/**
	 * Class contructor.
	 *
	 * The constructor either take one or three arguments. With only one argument, the arrow of numbers, the constructor sets the dimension of the 
	 * matrix as the dimensions of the array. In the three argument version, the constructor sets the dimension to those specified. This may be 		
	 * necessary since some columns/rows full of zeros may not be returned by some internal computations. This may also be useful to creat a matrix of 
	 * zeros by the user. 
	 *
	 * If the number of rows and the number of columns are specified, those must be equal or greater than the actual number of rows ans columns of
	 * the array as otherwise, the constructor will set the array values. 
	 *
	 * @param array $numbers, [int $numRows, int $numColumns] the array of numbers and the dimensions of the matrix.
	 * @return The new matrix
	 */
	function matrix()
	{
		//See which information is passed as arguments.
		$nArgs = func_num_args();
		$numbers = func_get_arg(0);

		if($nArgs > 2) //If the three arguments version, set the values.
		{
			$this->numRows = func_get_arg(1);
			$this->numColumns = func_get_arg(2);
		}
		
		//Routine to initialise the array of numbers.
		$numberColumns = 0;
		$numberRows = 0;

		if(empty($numbers) == false) //If the array is not empty
		{				
			foreach($numbers as $i => $rows) //Check every row...
			{
				foreach($rows as $j => $number) //...and column to remove zeros.
				{
					if($number != 0)
					{
						$this->numbers[$i][$j] = $number; 
					}
					//Find the biggest width of the array.			
					if($j >= $numberColumns) 
					{
						$numberColumns = $j;
					}
				}
				//Find the biggest depth of the array.
				if($i >= $numberRows)
				{
					$numberRows = $i;
				}
			}
			//Arrays starts at zero (so we add one)
			$numberRows++;
			$numberColumns++;
		}
		
		//If the user misspecified the number of rows/columns by a lower value. 
		if($numberRows > $this->numRows) 
		{
			$this->numRows = $numberRows;
		}
		
		if($numberColumns > $this->numColumns)
		{
			$this->numColumns = $numberColumns;
		}
	}
	
	/**
	* Quasi-constructor for an identity matrix.
	*
	* Creates an $n x $m matrix with ones on the main diagonal and zeros elsewhere. 
	* 
	* @param int $n, int $m the depth and width of the matrix.
	* @return The matrix
	*/
	function eye($n, $m)
	{		
		$dim = min($n, $m);
		for($i = 0; $i < $dim; $i++)
		{
			$id_numbers[$i][$i] = 1; 
		}
		
		$the_identity_matrix = new matrix($id_numbers, $n, $m);
		
		return $the_identity_matrix;
	}

	/**
	* Quasi-constructor of a matrix of zeros.
	*
	* Creates an $n x $m matrix with zeros everywhere.
	*
	* @param int $n, int $m the depth and width of the matrix.
	* @return The matrix
	*/
	function zeros($n, $m)
	{
		$the_zero_matrix = new matrix(array(), $n, $m);
		
		return $the_zero_matrix;
	}
	
	/**
	* Quasi-constructor of a matrix of ones
	*
	* Creates an $n x $m matrix with ones everywhere.
	*
	* @param int $n, int $m the depth and width of the matrix.
	* @return The matrix
	*/
	function ones($n, $m)
	{
		$ones = array();
		
		for($i = 0; $i < $n; $i++)
		{
			for($j = 0; $j < $m; $j++)
			{
				$ones[$i][$j] = 1; 
			}
		}
		
		$a_matrix_of_ones = new matrix($ones);
		
		return $a_matrix_of_ones;
	}

	/**
	* Quasi-constructor of a matrix of random numbers
	*
	* Creates an $n x $m matrix with random numbers. The numbers are taken from the inside PHP random number generator and normalised so to have a uniform
	* distribution.
	*
	* @param int $n, int $m the depth and width of the matrix.
	* @return The matrix
	*/
	function random($n, $m)
	{
		$random_numbers = array();
		
		for($i = 0; $i < $n; $i++)
		{
			for($j = 0; $j < $m; $j++)
			{
				$random_numbers[$i][$j] = rand(0,100)/100; 
			}
		}
		
		$the_random_matrix = new matrix($random_numbers, $n, $m);
		return $the_random_matrix;
	}

	/**
	 * Set all the numbers of this matrix
	 * This changes the actual numbers of this matrix to the specified arguments. The dimensions of the matrix are changed to those of the array passed as an
	 * argument.
	 *
	 * @param array $numbers the numbers of the matrix
	 */
	function set_data($numbers)
	{
		$a = new matrix($numbers);
		$this->numbers = $a->numbers;
		$this->numRows = $a->numRows;
		$this->numColumns = $a->numColumns;
	}
	
	/**
	 * Get the all the numbers of the matrix.
	 *
	 * Get the array of numbers of the matrix, including zeros.
	 *
	 * @return The array of numbers of the matrix.
	 */
	function get_data()
	{
		for($i = 0; $i < $this->numRows; $i++)
		{
			for($j = 0; $j < $this->numColumns; $j++)
			{
				if(empty($this->numbers[$i][$j]) == false)
				{
					$the_numbers[$i][$j] = $this->numbers[$i][$j];
				}
				else
				{
					$the_numbers[$i][$j] = 0;
				}
			}
		}
		return $the_numbers;
	}

	/**
	 * Get a particular element of the matrix
	 *
	 * Pre-requisite : The values of $i and $j must be smaller than the dimensions (rows/columns). $i and $j are indexed from one. Hence get_value(1,1) returns the
	 * value of the first element (in the upper left corner). If the pre-requisite is not met, an error message will be echoed to the standard output.
	 *
	 * @param : int $i, int $j the index of the particular element.
	 * @return The number requested OR zero if the parameters $i, $j are invalid.
	 */
	function get_value($i, $j)
	{
		$the_value = 0;
		if($i-1 < $this->get_num_rows() and $j-1 < $this->get_num_columns())
		{
			if(empty($this->numbers[$i-1][$j-1]) == false)
			{
				$the_value = $this->numbers[$i-1][$j-1];	
			}
		}
		else
		{
			echo "<br><br>\n\n\nWrong parameters in get_value !\n\n\n\<br><br>";
		}
		
		return $the_value;
	}
	
	/**
	 * Set a particular element of the matrix
	 *
	 * Pre-requisite : The values of $i and $j must be smaller than the dimensions (rows/columns). $i and $j are indexed from one. Hence set_value(1,1, 0) sets the
	 * value of the first element (in the upper left corner) to zero. If the pre-requisite is not met, an error message will be echoed to the standard output.
	 *
	 * If there is an error in indexes, no value will be changed. 
	 *	
	 * @param : int $i, int $j, double $value the coordinates of the element and the new value.
	 */
	function set_value($i, $j, $value)
	{
		if($i-1 < $this->get_num_rows() and $j-1 < $this->get_num_columns())
		{
			if($value != 0 )
		  	{
		  		$this->numbers[$i-1][$j-1] = $value;
		  	}
		  	else
		  	{
		  		unset($this->numbers[$i-1][$j-1]);
		  	}
		}
		else
		{
			echo "<br><br>\n\n\nWrong parameters in set_value !\n\n\n\<br><br>";
		}
	}

	/**
	 * Get the number of rows. 
	 *
	 * @return The number of rows of this matrix.
	 */
	function get_num_rows()
	{
		return $this->numRows;
	}

	/**
	 * Get the number of columns
	 *
	 * @return The number of columns of this matrix.
	 */
	function get_num_columns()
	{
		return $this->numColumns;
	}
	
	/**
	 * Get some rows of the matrix
	 *
	 * Returns a matrix containing all the rows inclusively contained in the range $start to $stop.
	 *
	 * Pre-requisite : $start and $stop must be inclusiveley between one an the number of rows of this matrix. If not, an error message will be echoed to the
	 * standard output
	 *
	 * @param int $start, int $stop the range of requested rows
	 * @return The matrix generated by those rows OR a 1x1 matrix with a zero if there is an error in the range of requested rows.
	 */
	function get_rows($start, $stop) 
	{	
		if($start-1 < $this->get_num_rows() and $stop - 1 < $this->get_num_rows())
		{
			
			for($i = $start-1; $i < $stop; $i++)
			{
				if(empty($this->numbers[$i]) == false)
				{
					$the_rows_array[$i - $start + 1] = $this->numbers[$i];
				}
			}
			$the_rows = new matrix($the_rows_array);
		}
		else
		{
			echo "<br><br>Invalid indexes in get_rows<br><br>";
			$the_rows = $this->zeros(1,1);
		}
		return $the_rows;
	}
	
	/**
	 * Get some columns of the matrix
	 *
	 * Returns a matrix containing all the columns inclusively contained in the range $start to $stop.
	 *
	 * Pre-requisite : $start and $stop must be inclusiveley between one an the number of columns of this matrix. If not, an error message will be echoed to the
	 * standard output.
	 *
	 * @param int $start, int $stop the range of requested rows
	 * @return The matrix formed by the requested columns OR a 1x1 matrix with a zero if there is an error in the range of requested columns.
	 */
	function get_columns($start, $stop) 
	{
		if($start-1 < $this->get_num_columns() and $stop - 1 < $this->get_num_columns())
		{
			$toto = $this->prime();
			$the_columns = $toto->get_rows($start, $stop);
			$the_columns = $the_columns->prime();
		}
		else
		{
			echo "<br><br>Invalid indexes in get_columns<br><br>";
			$the_columns = $this->zeros(1,1);
		}
		return $the_columns;
	}
	
	/**
	 * Replace some rows of the matrix
	 *
	 * Changes the rows specified in the inclusive range of $start and $stop by $bloc_matrix
	 *
	 * Pre-requisite : $start and $stop must be inclusiveley between one an the number of columns of this matrix. If not, an error message will be echoed to the
	 * standard output. If there is an error, no values will be changed.
	 * 
	 * Note that no requirements are necessary on $block_matrix. If it is smaller in any dimension, this will be interpreted as filling the missing values with
	 * zeros. If the matrix is bigger, only the necessary values will be taken, starting from the upper-left corner. 
	 *
	 * @param int $start, int $stop, matrix $block_matrix the range and the matrix to use for modifications.
	 */
	function set_rows($start, $stop, $block_matrix)
	{
		if($start-1 < $this->get_num_rows() and $stop < $this->get_num_rows())
		{
			if($bloc_matrix->get_num_columns() == $this->get_num_columns())
			{
				for($i = $start-1; $i< $stop; $i++)
				{
					$this->numbers[$i] = $block_matrix->numbers[$i-$start-1];
				}
			}
		}
		else
		{
			echo "<br><br>Invalid indexes in set_rows<br><br>";
		}
	}
	
	
	/**
	 * Replace some columns of the matrix
	 *
	 * Changes the columns specified in the inclusive range of $start and $stop by $bloc_matrix
	 *
	 * Pre-requisite : $start and $stop must be inclusiveley between one an the number of columns of this matrix. If not, an error message will be echoed to the
	 * standard output. If there is an error, no values will be changed.
	 * 
	 * Note that no requirements are necessary on $block_matrix. If it is smaller in any dimension, this will be interpreted as filling the missing values with
	 * zeros. If the matrix is bigger, only the necessary values will be taken, starting from the upper-left corner. 
	 *
	 * @param int $start, int $stop, matrix $block_matrix the range and the matrix to use for modifications.
	 */
	function set_columns($start, $stop, $block_matrix) 
	{
		if($start-1 < $this->get_num_columns() and $stop < $this->get_num_columns())
		{
			if($bloc_matrix->get_num_rows() == $this->get_num_rows())
			{
				$toto = $this->prime();
				$the_columns = $toto->get_rows($start, $stop, $block_matrix->prime());
				$toto = $the_columns->prime();
				$this->numbers = $thoto->numbers;
			}
		}
		else
		{
			echo "<br><br>Invalid indexes in set_columns<br><br>";
		}
	}
	
	/**
	 * Delete some rows of the matrix
	 *
	 * Deletes the rows specified in the inclusive range of $start and $stop.
	 *
	 * Pre-requisite : $start and $stop must be inclusiveley between one an the number of columns of this matrix. If not, an error message will be echoed to the
	 * standard output. If there is an error, no values will be changed.
	 * 
	 * @param int $start, int $stop the range of rows to be deleted.
	 */
	 function delete_rows($start, $end)
	 {
	 	if($start <= $this->get_num_rows() and $stop <= $this->get_num_rows() and $start > 0 and $start <= $end)
		{
			$p=0;	
			for($i = $start-1; $i < $end; $i++)
			{
				unset($this->numbers[$i]);
			}
			for($i = $start; $i < $this->numRows; $i++)
			{
				if(empty($this->numbers[$i]) == false)
				{
					$this->numbers[$i - ($end - $start + 1)] = $this->numbers[$i];
					unset($this->numbers[$i]);
				}
				
			}
			$this->numRows -= ($end - $start + 1);
		}
		else
		{
			echo "<br><br>Invalid indexes in delete_rows<br><br>";
		}
	 }
	 
	/**
	 * Delete some columns of the matrix
	 *
	 * Deletes the columns specified in the inclusive range of $start and $stop.
	 *
	 * Pre-requisite : $start and $stop must be inclusiveley between one an the number of columns of this matrix. If not, an error message will be echoed to the
	 * standard output. If there is an error, no values will be changed.
	 * 
	 * @param int $start, int $stop the range of the columns to be deleted.
	 */
	 function delete_columns($start, $end)
	 {
	 	if($start <= $this->get_num_columns() and $end <= $this->get_num_columns() and $start > 0 and $start <= $end)
		{
			$toto = $this->prime();
			$toto->delete_rows($start, $end);
			$toto = $toto->prime();

			$this->numbers = $toto->numbers;
			$this->numColumns = $toto->numColumns;
			$this->numRows = $toto->numRows;
		}
		else
		{
			echo "<br><br>Invalid indexes in delete_columns<br><br>";
		}
	 }
	 
	 /**
	 * Obtain the upper triangular part of this matrix.
	 *
	 * This returns a matrix with all the original numbers of the matrix above the principal diagonal and with zeros elsewhere.
	 *
	 * @return The upper diagonal matrix.
	 **/
	 function utrig()
	 {
	 	$utrig = new matrix($this->numbers);
	 	
	 	for($i = 0; $i < $utrig->get_num_rows(); $i++)
	 	{
	 		for($j = 0; $j < $utrig->get_num_columns(); $j++)
	 		{
	 			if($i >= $j)
	 			{
	 				unset($utrig->numbers[$i][$j]);
	 			}
	 		}
	 	}
	 
	 	return $utrig;
	 }
	 
	/**
	 * Obtain the lower triangular part of this matrix.
	 *
	 * This returns a matrix with all the original numbers of the matrix below the principal diagonal and with zeros elsewhere.
	 *
	 * @return The lower diagonal matrix.
	 **/
	 function ltrig()
	 {
	 	$ltrig = new matrix($this->numbers);
	 	
	 	for($i = 0; $i < $ltrig->get_num_rows(); $i++)
	 	{
	 		for($j = 0; $j < $ltrig->get_num_columns(); $j++)
	 		{
	 			if($i <= $j)
	 			{
	 				unset($ltrig->numbers[$i][$j]);
	 			}
	 		}
	 	}
	 
	 	return $ltrig;
	 }
	 
	 
	 
	 /**
	  * Obtain the diagonal of the matrix OR form a diagonal matrix with a supplied vector.
	  *
	  * Two possible usages : 
	  *  - Without any arguments : Obtain the elements on the main diagonal of the matrix.
	  *  - With a vertical matrix argument : Form a matrix with the vertical matrix on the diagonal.
	  *
	  * If any other number of arguments is supplied, the function will echo an error message.
	  * 
	  *	Pre-requisite : 
	  *  - Without any arguments : the matrix must be square
	  *  - With a vertical matrix argument : the parameter must be an 1 x n matrix (a vertical matrix).
	  *
	  * If the pre-requisites are not met, an error message will be echoed to the standard output.
	  *
	  * @return 
	  *  - Without any arguments : an n x 1 matrix with the diagonal elements of this matrix. 
	  *	 - With a vertical matrix : a square matrix with the 1 x n matrix on the principal diagonal.
	  *	 - If any error, a 1 x 1 matrix with a zero.
	  */
	 function diag()
	 {
	 	if(func_num_args() == 1)
	 	{
	 		$vector = func_get_arg(0);
	 		if($vector->get_num_columns() == 1 or $vector->get_num_rows() == 1)
	 		{
	 			if($vector->get_num_columns() == 1)
	 			{
	 				$vector->prime();
	 			}
	 			
	 			for($i = 0; $i < $vector->get_num_columns(); $i++)
	 			{
	 				if($i != 0)
	 				{	
	 					$vector->numbers[$i][$i] = $vector->numbers[0][$i];
	 					unset($vector->numbers[0][$i]);
	 				}
	 			}
	 			$vector->numRows = $vector->numColumns;
	 			
	 			$diag_matrix = $vector->prime();
	 		}
	 		else
	 		{
	 			echo "Wrong vector argument in diag function !<br><br>";
	 			$diag_matrix = $this->zeros(1,1);
	 		}
	 	}
	 	elseif(func_num_args() == 0)
	 	{
	 		if($this->is_square())
	 		{
	 			for($i = 0; $i < $this->get_num_columns(); $i++)
	 			{
	 				if(empty($this->numbers[$i][$i]) == false)
	 				{
	 					$diag_data[$i][0] = $this->numbers[$i][$i];
	 				}
	 				
	 			}
	 			$diag_matrix = new matrix($diag_data, $this->get_num_rows(), 1);
	 		}
	 		else
	 		{
		 		echo "This matrix is not square for diag function !<br><br>";
		 		$diag_matrix = $this->zeros(1,1);
	 		}
	 	}
	 	else
	 	{
	 		echo "Wrong arguments in the diag function !<br><br>";
	 		$diag_matrix = $this->zeros(1,1);
	 	}
	 	
	 	return $diag_matrix;
	 }
	 

	/**
	 * Check if two matrices are of the same size.
	 *
	 * E.g. : checks if they both have the same number of rows and columns.
	 *
	 * @param matrix $someMatrix the matrix to compare with this matrix.
	 * @return bool true or false depending on if they are of the same size or not.
	 */
	function size_eq($someMatrix)
	{
		$return = false;
			
		if ($someMatrix->get_num_rows() == $this->get_num_rows() and $someMatrix->get_num_columns() == $this->get_num_columns())
		{
			$return = true;
		}
			
		return $return;
	}
		

	/**
	 * Check if a matrix is square.
	 *
	 * E.g. : if the number of rows is equal to the number of columns.
	 *
	 * @return bool true or false, depending on the result of the test.
	 */
	function is_square()
	{
		$return = false;
		if ($this->get_num_rows() == $this->get_num_columns())
		{
			 $return = true;
		}
		
			 
		return $return;
	}


	/**
	 * Substract $someMatrix from this matrix.
	 *
	 * A - B is $A->minus($B)
	 *
	 * Pre-requisite: Both matrices must have the same size. If not, an error message will be echoed to standard output.
	 *
	 * @param matrix $someMatrix the matrix to substract.
	 * @return The matrix resulting of the substraction OR a matrix of zeros if the pre-requisite is not met.
	 */
	function minus($someMatrix)
	{
		$substract = $this->zeros($this->get_num_rows(), $this->get_num_columns());
		if($this->size_eq($someMatrix))
		{
			for($i = 0; $i < $this->get_num_rows(); $i++)
			{   
				for($j = 0; $j < $this->get_num_columns(); $j++)
				{
					if(empty($this->numbers[$i][$j]) == false)
					{
						if(empty($someMatrix->numbers[$i][$j]) == false)
						{
							$substract->numbers[$i][$j] = $this->numbers[$i][$j] - $someMatrix->numbers[$i][$j];
						}
						else
						{
							$substract->numbers[$i][$j] = $this->numbers[$i][$j];
						}
					}
					else
					{
						if(empty($someMatrix->numbers[$i][$j]) == false)
						{
							$substract->numbers[$i][$j] = -1*$someMatrix->numbers[$i][$j];
							
						}
					}
					
				}
			}
		}
		else
		{
			echo "\n\n\n Wrong dimensions in matrix substraction \n\n\n";
		}
		return $substract;
	}
	
	/**
	 * Add $someMatrix to this matrix
	 *
	 * A + B is $A->plus($B)
	 *
	 * Pre-requisite: Both must have the same size. Otherwise, an error message is echoed to the standard output.
	 *
	 * @param matrix $someMatrix, the matrix to add.
	 * @return The matrix resulting of the addition OR a matrix of zeros if the pre-requisite is not met.
	 */
	function plus($someMatrix)
	{	
		$add = new matrix($this->numbers, $this->get_num_rows(), $this->get_num_columns());
		 
		if($this->size_eq($someMatrix))
		{

			$someMatrix = $someMatrix->s_times(-1);
			$add = $this->minus($someMatrix);
			
		}
		else
		{
			echo "\n\n\n Wrong dimensions in matrix addition \n\n\n";
		}
		
		return $add;
	}	
	


	/**
	 * Multiply this matrix by $someMatrix
	 *
	 * A x B is $A->times($B)
	 *
	 * Pre-requiste : the number of columns of this matrix $this must equal the number of rows of $someMatrix. Otherwise, an error message is echoed to the standard
	 * output. 
	 *
	 * @param numbers $someMatrix the matrix to multiply.
	 * @return The result of the matrix product OR a matrix of zeros if the pre-requisite is not met.
	 */
	function times($someMatrix) 
	{
		$easier = $someMatrix->prime();
		
		if($this->get_num_columns() == $someMatrix->get_num_rows())
		{	
			foreach($this->numbers as $i => $row) 
			{	
				foreach($easier->numbers as $j => $column)
				{
					$total = 0;
					foreach($row as $k => $number)
					{
						if(empty($column[$k]) == false)
						{
							$total += $number * $column[$k];
						}				
					}
					$the_product_data[$i][$j] = $total;
				}
			}
			$the_product = new matrix($the_product_data);
		}
		else
		{
			echo "\n\n\n Wrong dimensions in matrix multiplication \n\n\n";
			$the_product = $this->zeros(1,1);
		}
		return $the_product;
	}

	/**
	 * Multiply this matrix by some scalar.
	 *
	 * A x scalar is $A->s_times($scalar)
	 *
	 * @param double $value
	 */
	function s_times($value) 
	{	
		$the_matrix = new matrix($this->numbers, $this->get_num_rows(), $this->get_num_columns());
		foreach($this->numbers as $i => $column)
		{
			foreach($column as $j => $number)
			{
				$the_matrix->numbers[$i][$j] *= $value;
			}
		}
		
		return $the_matrix;
	}
	
	/**
	 * Compute the determinant of some matrix
	 *
	 * Pre-requisite: the matrix must be square
	 *
	 * This function computes the determinant of the matrix through a PLU decomposition. Hence, it computes the product of the main diagonal of the U matrix.
	 *
	 * @param matrix $someMatrix
	 * @return int $the_determinant 
	 */
	function det($someMatrix) 
	{
		$the_determinant = 0;		
		if($someMatrix->is_square())
		{
			$the_determinant = 1;
		  	$PLU = $someMatrix->PLU();
		  	for($j = 0; $j < $this->get_num_rows(); $j++)
		  	{
		  		$the_determinant =	$the_determinant * $PLU['LU'][$j][$j];
		  	}
		  	$the_determinant = $the_determinant * $PLU['P_set'];
		}		
		return $the_determinant;
	}


	/**
	 * Finds a minor matrix. 
	 *
	 * A "minor matrix" is the matrix which corresponds to deleting a specified row and column of the original matrix.
	 *
	 * Pre-requisite : the $column and $row values must be within the existing range of this matrix. Otherwise, an error message will be echoed to the standard
	 * output.
	 *
	 * @param matrix $someMatrix, int $column, int $row that is, the matrix from which we seek the minor matrix, the column and the row to delete.
	 * @return The minor matrix OR a matrix of zeros if the pre-requisite is not met.
	 */
	function minor_matrix($someMatrix, $column, $row) 
	{
		$the_result = new matrix($someMatrix->numbers);
		$the_result->delete_rows($row, $row);
		$the_result->delete_columns($column, $column);
		return $the_result;
	}


	/**
	 * Compute the transpose of this matrix
	 *
	 * E.g. : A_ij becomes A_ji 
	 *
	 * A' is $A->prime() 
	 *
	 * @return The transposed matrix.
	 */
	function prime() 
	{	
		foreach($this->numbers as $i => $row)
		{
			foreach($row as $j => $number)
			{
				$the_transpose_data[$j][$i] = $number;
			}
		}
		$the_transpose = new matrix($the_transpose_data, $this->get_num_columns(), $this->get_num_rows());
		
		return $the_transpose;
	}

	/**
	 * Smoothe the elements of this matrix to zero
	 *
	 * All elements within (-$tolerance, $tolerance) are set to zero.
	 *
	 * @param : $tolerance
	 *	
	 */
	function smoothe($tolerance)
	{
		foreach($this->numbers as $i => $row)
		{
			foreach($row as $j => $number)
			{
				if(abs($number)<$tolerance)
				{
					unset($this->numbers[$i][$j]);
				}
			}
		}
	}

	/**
	 * Compute the inverse of this matrix
	 *
	 * This is done in a few steps : 
	 *	- Find the PLU decomposition of the matrix. That is find P, L and U where P is a permutation matrix, L is a lower diagonal matrix with ones on the diagonal
	 *  and U is an upper triangular matrix such that P*L*U is the original matrix. (See the private function PLU)
	 *	- Check if the determinant is non-null. An error is echoed to the standard output otherwise.
	 *	- Perform the inversion of U, L, and P (which is much faster).
	 *  - Multiply the inverted results to find the inverse of the original matrix.
	 *
	 * Pre-requisite : this matrix is invertible.
	 *	
	 * @return The inverse matrix OR a matrix of zeros if the pre-requisite is not met.
	 */
	function inv() 
	{	
		if($this->is_square())
		{
			$n = $this->get_num_rows();
			$the_determinant = 1;
		  	$PLU = $this->PLU();
		  	$inv = $PLU['LU'];
		  	for($j = 0; $j < $this->get_num_rows(); $j++)
		  	{
		  		$the_determinant *=	$PLU['LU'][$j][$j];
		  	}
		  	$the_determinant = $the_determinant * $PLU['P_set'];
		  	
			if($the_determinant != 0)
			{
			
				//Perform the inversion algorithm. 
				for($i = 0; $i < $n; $i++)
				{
					//Inverse of the main diagonal of U
					$inv[$i][$i]  = 1/$PLU['LU'][$i][$i];				
				}
				//For each diagonal except the principal...
				for($i = 1; $i < $n; $i++)
				{	
					//Now, for each element of this diagonal... (there are $n - $i of them)
					for($j = 0; $j < $n - $i; $j++)
					{
						//For all elements of this row from this element to the principal diagonal
						//Except the first one
						//(there are $i-1 of them)
						$temp1 = 0; 
						$temp2 = 0;
						for($k = 0; $k <= $n; $k++)
						{
							//We start with the elements of L
							if($k > $j and $k <= $i + $j)
							{
								if($k == $i + $j)
								{
									$temp1 -= $PLU['LU'][$k][$j];
								}
								else
								{	
									$temp1 -= $inv[$i+$j][$k]*$PLU['LU'][$k][$j];
								}
							}
							
							//And then U
							if($k >= $j and $k < $i + $j)
							{
								$temp2 -= $inv[$j][$k] * $PLU['LU'][$k][$i + $j];
							}
							
							
						}
						$inv[$j + $i][$j] = $temp1;
						$inv[$j][$i+$j] = $temp2 * $inv[$i + $j][$i + $j];  
					}				
				}
				$PLU['LU'] = $inv;
				unset($inv);
				//Now, we perform the multiplication of U_inv * L_inv
				for($i = 0; $i < $n; $i++)
				{
					for($j = 0; $j < $n; $j++)
					{
						for($k = max($j, $i); $k <= $n; $k++)
						{
							if($k == max($j, $i))
							{
								if($k == $i)
								{
									$inv[$j][$i] = $PLU['LU'][$j][$k];						
								}
								else
								{
									$inv[$j][$i] = $PLU['LU'][$k][$i] * $PLU['LU'][$j][$k];						
								}
								
							}
							else
							{
								if($k == $i)
								{
									$inv[$j][$i] += $PLU['LU'][$j][$k];						
								}
								else
								{
									$inv[$j][$i] += $PLU['LU'][$k][$i]*$PLU['LU'][$j][$k];						
								}
							}					
						}
					}
				}
				//Finally, we multiply UL_inv * P_inv to get the inverse
				$PLU['LU'] = $inv;
				unset($inv);
				for($i = 0; $i < $n; $i++)
				{
					for($j = 0; $j < $n; $j++)
					{
						for($k = 0; $k < $n; $k++)
						{
							if($k == 0)
							{
								$inv[$i][$j] = $PLU['LU'][$i][$k] * $PLU['P'][$k][$j];						
							}
							else
							{
								$inv[$i][$j] += $PLU['LU'][$i][$k] * $PLU['P'][$k][$j];						
							}
							
						}
					}
				}
				$the_inverse = new matrix($inv);
			}
			else
			{
				echo "\n\n\n Matrix has a zero determinant !! \n\n\n";
				$the_inverse = $this->zeros(1,1);
			}
		}
		else
		{
			echo "\n\n\n Matrix is not square !! \n\n\n";
			$the_inverse = $this->zeros(1,1);
		}
		return $the_inverse;
	}
	
	/**
	 * Compute each element of this matrix with each element of some matrix
	 *
	 * E.g. this_ij * someMatrix_ij
	 *
	 * Pre-requisite : both matrix must have the same size. Otherwise, an error message is echoed to the standard output.
	 *
	 * @param matrix $someMatrix the matrix used to perform the element by element multiplication.
	 * @return The element by element product matrix OR a matrix of zeros if the pre-requisite is not met.
	 */
	function p_times($someMatrix)
	{
		
		if($this->size_eq($someMatrix))
		{
			foreach($this->numbers as $i => $row)
			{
				foreach($row as $j => $number)
				{
					if(empty($someMatrix->numbers[$i][$j]) == false)
					{
						$the_point_to_point[$i][$j] = $number * $someMatrix->numbers[$i][$j];
					}
				}
			}
			$the_point_to_point_product = new matrix($the_point_to_point, $this->get_num_rows(), $this->get_num_columns());
		}
		else
		{
			echo "\n\n\n wrong dimensions in point product !! \n\n\n";
			$the_point_to_point_product = $this->zeros(1,1);
		}
		
		return $the_point_to_point_product;
	}

	/**
	 * Compute element by element division
	 *	 
	 * E.g. this_ij / someMatrix2_ij
	 *
	 * Pre-requisites : both matrix must have the same size and no elements of someMatrix is zero. Otherwise an error message is echoed to the standard output.
	 *
	 * @param matrix $someMatrix the matrix of divisors.
	 * @return The element by element division matrix OR a matrix of zeros if the pre-requisites are not met.
	 */
	function p_div($someMatrix)
	{
		$noZeros = false;
		foreach($someMatrix->numbers as $key => $rows)
		{
			if(count($rows) < $someMatrix->get_num_rows())
			{
				$noZeros = true;
			}
		}	
		if($this->size_eq($someMatrix) and $noZeros == false)
		{
			foreach($this->numbers as $i => $rows)
			{
				foreach($rows as $j => $numbers)
				{
					$the_point_to_point[$i][$j] = $numbers / $someMatrix->numbers[$i][$j];
				}
			}
		}
		else
		{
			echo "\n\n\n Wrong dimensions in point division and/or the matrix two has some zeros !! \n\n\n";
		}
		$the_point_to_point_division = new matrix($the_point_to_point, $someMatrix->get_num_rows(), $someMatrix->get_num_columns());
		
		return $the_point_to_point_division;
	}
	
	/**
	 * Finds the max of each column of this matrix
	 *
	 * @return An 1 x m matrix with the maximum of each columns.
	 */
	function mat_max()
	{
		$i = 0;
		$position = -1; 
		while($i < $this->get_num_rows())
		{
			if(empty($this->numbers[$i]) == false)
			{
				$position = $i;
				$i = $this->get_num_rows();
			}
			else
			{
				$i++;
			}
		}
		if($position > -1)
		{
			$max[0] = $this->numbers[$position];
			foreach($this->numbers as $i => $rows)
			{
				foreach($rows as $j => $number)
				{
					if($number > $max[0][$j])
					{
						$max[0][$j] = $number;
					}
				}
			}	
		}
		$the_max = new matrix($max, 1, $this->get_num_columns());

		return $the_max;
	}
	
	
	/**
	 * Compute the min of each column of this matrix
	 *
	 * @return An 1 x m matrix with the minimum of each columns.
	 */
	function mat_min()
	{
		$i = 0;
		$position = -1; 
		while($i < $this->get_num_rows())
		{
			if(empty($this->numbers[$i]) == false)
			{
				$position = $i;
				$i = $this->get_num_rows();
			}
			else
			{
				$i++;
			}
		}
		if($position > -1)
		{
			$min[0] = $this->numbers[$position];
			foreach($this->numbers as $i => $rows)
			{
				foreach($rows as $j => $number)
				{
					if($number < $min[0][$j])
					{
						$min[0][$j] = $number;
					}
				}
			}	
		}
		$the_min = new matrix($min, 1, $this->get_num_columns());

		return $the_min;
	}

	/**
	 * Element by element greater than comparison.
	 *
	 * E.g. 1 if this_ij > someMatrix_ij and 0 otherwise. 
	 *
	 * Pre-requiste : both matrices must have the same dimensions. Otherwise an error message is echoed to the standart output.
	 *
	 * @param matrix $someMatrix
	 * @return A matrix of zeros and ones OR a matrix of zeros if the pre-requisite is not met. 
	 */
	function gt($someMatrix)
	{
		if($this->size_eq($someMatrix) == true)
		{
			$what_is_greater = $this->zeros($this->get_num_rows(), $this->get_num_columns());
			
			for($i = 0; $i < $this->get_num_rows(); $i++)
			{
				for($j = 0; $j < $this->get_num_columns(); $j++)
				{
					if(empty($this->numbers[$i][$j])== false)
					{
						if(empty($someMatrix->numbers[$i][$j]) == false)
						{
							if($this->numbers[$i][$j] > $someMatrix->numbers[$i][$j])
							{
								$what_is_greater->numbers[$i][$j] = 1;
							}
						}
						else
						{
							if($this->numbers[$i][$j] > 0)
							{
								$what_is_greater->numbers[$i][$j] = 1;
							}
						}
					}
					else
					{
						if(0 > $someMatrix->numbers[$i][$j])
						{
							$what_is_greater->numbers[$i][$j] = 1;
						}
					}
					
				}
			}
		}
		else
		{
			echo "Wrong dimensions in gt !!!<br><br>";
			$what_is_greater = $this->zeros(1,1);
		};
		return $what_is_greater;
	}
	
	/**
	 * Element by element greater than or equal comparison.
	 *
	 * E.g. 1 if this_ij >= someMatrix_ij and 0 otherwise. 
	 *
	 * Pre-requiste : both matrices must have the same dimensions. Otherwise an error message is echoed to the standart output.
	 *
	 * @param matrix $someMatrix
	 * @return A matrix of zeros and ones OR a matrix of zeros if the pre-requisite is not met. 
	 */
	function gteq($someMatrix)
	{
		
		if($this->size_eq($someMatrix))
		{
			$what_is_greater = $this->zeros($this->get_num_rows(), $this->get_num_columns());
			for($i = 0; $i < $this->get_num_rows(); $i++)
			{
				for($j = 0; $j < $this->get_num_columns(); $j++)
				{
					if(empty($this->numbers[$i][$j])== false)
					{
						if(empty($someMatrix->numbers[$i][$j]) == false)
						{
							if($this->numbers[$i][$j] >= $someMatrix->numbers[$i][$j])
							{
								$what_is_greater->numbers[$i][$j] = 1;
							}
						}
						else
						{
							if($this->numbers[$i][$j] >= 0)
							{
								$what_is_greater->numbers[$i][$j] = 1;
							}
						}
					}
					else
					{
						if(0 >= $someMatrix->numbers[$i][$j])
						{
							$what_is_greater->numbers[$i][$j] = 1;
						}
					}
				}
			}
		}
		else
		{
			echo "Wrong dimensions in gteq !!!<br><br>";
			$what_is_greater = $this->zeros(1,1);
		}
		
		return $what_is_greater;
	}
	
	/**
	 * Element by element less than comparison.
	 *
	 * E.g. 1 if this_ij < someMatrix_ij and 0 otherwise. 
	 *
	 * Pre-requiste : both matrices must have the same dimensions. Otherwise an error message is echoed to the standart output.
	 *
	 * @param matrix $someMatrix
	 * @return A matrix of zeros and ones OR a matrix of zeros if the pre-requisite is not met. 
	 */
	function lt($someMatrix)
	{
		
		if($this->size_eq($someMatrix))
		{
			$what_is_smaller = $this->ones($this->get_num_rows(), $this->get_num_columns());
			$what_is_smaller = $what_is_smaller->minus($this->gt($someMatrix));
		}
		else
		{
			echo "Wrong dimensions in lt !!!<br><br>";
			$what_is_smaller = $this->zeros(1,1);
		}
		
		return $what_is_smaller;
	}
	
	/**
	 * Element by element greater than comparison.
	 *
	 * E.g. 1 if this_ij <= someMatrix_ij and 0 otherwise. 
	 *
	 * Pre-requiste : both matrices must have the same dimensions. Otherwise an error message is echoed to the standart output.
	 *
	 * @param matrix $someMatrix
	 * @return A matrix of zeros and ones OR a matrix of zeros if the pre-requisite is not met. 
	 */
	function lteq($someMatrix)
	{
		
		if($this->size_eq($someMatrix))
		{
			$what_is_smaller = $this->ones($this->get_num_rows(), $this->get_num_columns());
			$what_is_smaller = $what_is_smaller->minus($this->gteq($someMatrix));
		}
		else
		{
			echo "Wrong dimensions in lteq !!!<br><br>";
			$what_is_smaller = $this->zeros(1,1);
		}
		
		return $what_is_smaller;
		return $what_is_smaller;
	}

	
	/**
	 * Get the sum of each column.
	 * 	
	 * @return An 1 x m matrix with the sum of each column.
	 */
	function sum()
	{
		foreach($this->numbers as $i => $rows)
		{
			foreach($rows as $j => $number)
			{
				$sum[0][$j] += $number;
			}
		}
		$the_sum = new matrix($sum, 1, $this->get_num_columns());
		
		return $the_sum;
	}
	
	/**
	 * 	Compute the mean of each column. 
	 *
	 *	@return An 1 x m matrix with the mean of each column.
	 */
	function mean()
	{
		$the_mean = $this->sum();
		
		$columns = $the_mean->get_num_columns();
		foreach($the_mean->numbers[0] as $j => $number)
		{
			$the_mean->numbers[0][$j] = $number / $columns;
		}

		return $the_mean;
	}
	
	
	/**
	 * 	Performs a PLU decomposition of this matrix.
	 * 	
	 * 	Pre-requisite : this matrix must be square. Otherwise an error message is echoed to the standard output.
	 *
	 *	This function sets two arrays containing the data for matrices : "P" "L", "U".
	 *
	 *	@return : The array $PLU which consists of $PLU['P_det'], $PLU['P'], $PLU['LU']. 
	 *
	 *  Example : if the matrix is : 2    1    1
	 *								 1    2    1
	 *								 1    1    2
	 *
	 *  The decomposition is then : 
	 *
	 *	L			 				 1    0    0
	 *								 1/2  1    0
	 *								 1/2  4/3  1
	 *
	 *	U			 				 2    1    1
	 *								 0    3/2  1/2
	 *								 0    0    4/3
	 *
	 *  P (determinant = 1)			 1    0    0 
	 *								 0    1    0
	 *								 0    0    1
	 *
	 *
	 *	So the return arrays will be : 
	 * 	$PLU['P'] 					 1    0    0 
	 *								 0    1    0
	 *								 0    0    1
	 *
	 *	$PLU['P_det']	: 1
	 *
	 *	And the we omit the ones in the L matrix to get a "stacked" array with LU : 
	 *
	 *	$PLU['LU'] together :		2    1    1
	 *								1/2  3/2  1/2
	 *								1/2  4/3  1							
	 *
	 */
	private function PLU()
	{
	 	if($this->is_square())
		{
			$PLU['LU'] = $this->numbers;
			$PLU['P'] = array();
			// Set the permutation matrix. 
			$PLU['P_set'] = 1;
			for($k = 0; $k < $this->get_num_rows(); $k++)
			{
				$PLU['P'][$k][$k] = 1; 
			}
		
			//For all rows of the matrix....
			for($k = 0; $k < $this->get_num_rows(); $k++)
			{
	    		//Find the highest value of the actual column, starting at $k and switch it.
	    		$row = $k;
	    		for($j = $k; $j < $this->get_num_rows(); $j++)
	    		{
	    			if($PLU['LU'][$j][$k] > $PLU['LU'][$row][$k])
	    			{
	    				$row = $j;
	    			}	
	    		}
	    		
	    		if($row != $k) //If they are not equal, we have to switch rows
	    		{
	    			//Row switching on L
	    			$temp = $PLU['LU'][$k];
	    			$PLU['LU'][$k] = $PLU['LU'][$row];
	    			$PLU['LU'][$row] = $temp;
	    				    			
	    			//Row switching on P
	    			$temp = $PLU['P'][$k];
	    			$PLU['P'][$k] = $PLU['P'][$row];
	    			$PLU['P'][$row] = $temp;
	    			$PLU['P_set'] *= -1;
	    		}
	    		
	    		//Then, we perform the algorithm on L 
	    		//(note : you gotta admire this. Such a small piece of code for a so powerful result)
	    		
	    		//For all the values in the column $k
	    		for($i = $k+1; $i < $this->get_num_rows(); $i++)
	    		{	
	    			//If we are below the diagonal, then we work on L...
			        if(empty($PLU['LU'][$i][$k]) == false and $i > $k)
			        {
			        	$PLU['LU'][$i][$k] = $PLU['LU'][$i][$k]/$PLU['LU'][$k][$k];
			        }
			        
			        //Now for all elements of the U matrix...
			    	for($j = $k+1; $j < $this->get_num_rows(); $j++)
			    	{
	       				if(empty($PLU['LU'][$i][$j]) == false)
	           			{
	           				if(empty($PLU['LU'][$k][$j]) == false)
	           				{
	           					$PLU['LU'][$i][$j] = $PLU['LU'][$i][$j] - $PLU['LU'][$i][$k] * $PLU['LU'][$k][$j];
	           				}
	           			}
	           			else
	           			{
	           				if(empty($PLU['LU'][$k][$j]) == false)
	           				{
	           					$PLU['LU'][$i][$j] = - $PLU['LU'][$i][$k] * $PLU['LU'][$k][$j];
	           				}
	           			}
	        		}
	    		}
			}	
			//P needs to be transposed
	    	foreach($PLU['P'] as $i => $row)
	    	{
	    		foreach($row as $j => $value)
	    		{
	    			$new_P[$j][$i] = $value; 
	    		}
	    	}
		}
		else
		{
			echo "Wrong dimensions in LU !!!<br><br>";
			$PLU['LU'] = $this->numbers;
		}
		return $PLU;
	}
	
	/**
	* Performs a PLU decomposition of this matrix.
	* 	
	* Pre-requisite :this matrix must be square. Otherwise an error message is echoed to the standard output.
	*
	* @return : A 2n x n matrix where :
	*			- the first n x n matrix is the permutation matrix.
	*			- the second n x n matrix is the LU decomposition "stacked" where "ones" of the lower diagonal matrix
	*			  are omitted (by convention).
	*
	*  Example : if the matrix is the 3 x 3 :   
	* 
	* 2    1    1
	*
	* 1    2    1
	*
	* 1    1    2
	*
	*  The decomposition is then : 
	*
	*	L:
	*
	* 1    0    0
	*
	* 1/2  1    0
	*
	* 1/2  4/3  1
	*
	*	U:
	*
	* 2    1    1
	*
	* 0    3/2  1/2
	*
	* 0    0    4/3
	*
	*   P:
	*
	* 1    0    0 
	*
	* 0    1    0
	*
	* 0    0    1
	*
	*
	*	So the return matrix will be : 
	*
	* 1    0    0   2    1    1
	*
	* 0    1    0   1/2  3/2  1/2
	*
	* 0    0    1   1/2  4/3  1			
	*
	*/
	function lu()
	{
		$PLU = $this->PLU();
		$n = $this->get_num_rows();
		for($i = 0; $i < $n; $i++)
		{
			for($j = $n; $j < $n*2; $j++)
			{
				$PLU['P'][$i][$j] = $PLU['LU'][$i][$j - $n];
			}
		}
		$the_lu = new matrix($PLU['P']);
		return $the_lu;
	}

		
	/**
	* Prints the matrix
	*
	* Prints the matrix in an ugly table with borders in the standard output. This is intended for debugging purpose.
	*/
	function print_matrix()
	{
		$numb = $this->numbers;
		echo '<table border="1">'."\n";
		for($i = 0; $i < $this->get_num_rows();$i++)
		{
			echo "<tr>";
			for($j = 0; $j < $this->get_num_columns(); $j++)
			{
				if(empty($numb[$i][$j]) == false)
				{
					echo "<td>".$numb[$i][$j]."</td>";
				}
				else
				{
					echo "<td>0</td>";
				}
				
			}
			echo "<tr>\n";
		}
		echo "</table>\n";
	}

}
?>
Return current item: Matrix new