/PHPExcel/Shared/JAMA/tests/TestMatrix.php
PHP | 415 lines | 261 code | 53 blank | 101 comment | 56 complexity | 19723aca94cb1db21a8ee2346d65c4bd MD5 | raw file
1<?php 2 3require_once "../Matrix.php"; 4 5class TestMatrix { 6 7 function TestMatrix() { 8 9 // define test variables 10 11 $errorCount = 0; 12 $warningCount = 0; 13 $columnwise = array(1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.); 14 $rowwise = array(1.,4.,7.,10.,2.,5.,8.,11.,3.,6.,9.,12.); 15 $avals = array(array(1.,4.,7.,10.),array(2.,5.,8.,11.),array(3.,6.,9.,12.)); 16 $rankdef = $avals; 17 $tvals = array(array(1.,2.,3.),array(4.,5.,6.),array(7.,8.,9.),array(10.,11.,12.)); 18 $subavals = array(array(5.,8.,11.),array(6.,9.,12.)); 19 $rvals = array(array(1.,4.,7.),array(2.,5.,8.,11.),array(3.,6.,9.,12.)); 20 $pvals = array(array(1.,1.,1.),array(1.,2.,3.),array(1.,3.,6.)); 21 $ivals = array(array(1.,0.,0.,0.),array(0.,1.,0.,0.),array(0.,0.,1.,0.)); 22 $evals = array(array(0.,1.,0.,0.),array(1.,0.,2.e-7,0.),array(0.,-2.e-7,0.,1.),array(0.,0.,1.,0.)); 23 $square = array(array(166.,188.,210.),array(188.,214.,240.),array(210.,240.,270.)); 24 $sqSolution = array(array(13.),array(15.)); 25 $condmat = array(array(1.,3.),array(7.,9.)); 26 $rows = 3; 27 $cols = 4; 28 $invalidID = 5; /* should trigger bad shape for construction with val */ 29 $raggedr = 0; /* (raggedr,raggedc) should be out of bounds in ragged array */ 30 $raggedc = 4; 31 $validID = 3; /* leading dimension of intended test Matrices */ 32 $nonconformld = 4; /* leading dimension which is valid, but nonconforming */ 33 $ib = 1; /* index ranges for sub Matrix */ 34 $ie = 2; 35 $jb = 1; 36 $je = 3; 37 $rowindexset = array(1,2); 38 $badrowindexset = array(1,3); 39 $columnindexset = array(1,2,3); 40 $badcolumnindexset = array(1,2,4); 41 $columnsummax = 33.; 42 $rowsummax = 30.; 43 $sumofdiagonals = 15; 44 $sumofsquares = 650; 45 46 /** 47 * Test matrix methods 48 */ 49 50 /** 51 * Constructors and constructor-like methods: 52 * 53 * Matrix(double[], int) 54 * Matrix(double[][]) 55 * Matrix(int, int) 56 * Matrix(int, int, double) 57 * Matrix(int, int, double[][]) 58 * constructWithCopy(double[][]) 59 * random(int,int) 60 * identity(int) 61 */ 62 echo "<p>Testing constructors and constructor-like methods...</p>"; 63 64 $A = new Matrix($columnwise, 3); 65 if($A instanceof Matrix) { 66 $this->try_success("Column-packed constructor..."); 67 } else 68 $errorCount = $this->try_failure($errorCount, "Column-packed constructor...", "Unable to construct Matrix"); 69 70 $T = new Matrix($tvals); 71 if($T instanceof Matrix) 72 $this->try_success("2D array constructor..."); 73 else 74 $errorCount = $this->try_failure($errorCount, "2D array constructor...", "Unable to construct Matrix"); 75 76 $A = new Matrix($columnwise, $validID); 77 $B = new Matrix($avals); 78 $tmp = $B->get(0,0); 79 $avals[0][0] = 0.0; 80 $C = $B->minus($A); 81 $avals[0][0] = $tmp; 82 $B = Matrix::constructWithCopy($avals); 83 $tmp = $B->get(0,0); 84 $avals[0][0] = 0.0; 85 /** check that constructWithCopy behaves properly **/ 86 if ( ( $tmp - $B->get(0,0) ) != 0.0 ) 87 $errorCount = $this->try_failure($errorCount,"constructWithCopy... ","copy not effected... data visible outside"); 88 else 89 $this->try_success("constructWithCopy... ",""); 90 91 $I = new Matrix($ivals); 92 if ( $this->checkMatrices($I,Matrix::identity(3,4)) ) 93 $this->try_success("identity... ",""); 94 else 95 $errorCount = $this->try_failure($errorCount,"identity... ","identity Matrix not successfully created"); 96 97 /** 98 * Access Methods: 99 * 100 * getColumnDimension() 101 * getRowDimension() 102 * getArray() 103 * getArrayCopy() 104 * getColumnPackedCopy() 105 * getRowPackedCopy() 106 * get(int,int) 107 * getMatrix(int,int,int,int) 108 * getMatrix(int,int,int[]) 109 * getMatrix(int[],int,int) 110 * getMatrix(int[],int[]) 111 * set(int,int,double) 112 * setMatrix(int,int,int,int,Matrix) 113 * setMatrix(int,int,int[],Matrix) 114 * setMatrix(int[],int,int,Matrix) 115 * setMatrix(int[],int[],Matrix) 116 */ 117 print "<p>Testing access methods...</p>"; 118 119 $B = new Matrix($avals); 120 if($B->getRowDimension() == $rows) 121 $this->try_success("getRowDimension..."); 122 else 123 $errorCount = $this->try_failure($errorCount, "getRowDimension..."); 124 125 if($B->getColumnDimension() == $cols) 126 $this->try_success("getColumnDimension..."); 127 else 128 $errorCount = $this->try_failure($errorCount, "getColumnDimension..."); 129 130 $barray = $B->getArray(); 131 if($this->checkArrays($barray, $avals)) 132 $this->try_success("getArray..."); 133 else 134 $errorCount = $this->try_failure($errorCount, "getArray..."); 135 136 $bpacked = $B->getColumnPackedCopy(); 137 if($this->checkArrays($bpacked, $columnwise)) 138 $this->try_success("getColumnPackedCopy..."); 139 else 140 $errorCount = $this->try_failure($errorCount, "getColumnPackedCopy..."); 141 142 $bpacked = $B->getRowPackedCopy(); 143 if($this->checkArrays($bpacked, $rowwise)) 144 $this->try_success("getRowPackedCopy..."); 145 else 146 $errorCount = $this->try_failure($errorCount, "getRowPackedCopy..."); 147 148 /** 149 * Array-like methods: 150 * minus 151 * minusEquals 152 * plus 153 * plusEquals 154 * arrayLeftDivide 155 * arrayLeftDivideEquals 156 * arrayRightDivide 157 * arrayRightDivideEquals 158 * arrayTimes 159 * arrayTimesEquals 160 * uminus 161 */ 162 print "<p>Testing array-like methods...</p>"; 163 164 /** 165 * I/O methods: 166 * read 167 * print 168 * serializable: 169 * writeObject 170 * readObject 171 */ 172 print "<p>Testing I/O methods...</p>"; 173 174 /** 175 * Test linear algebra methods 176 */ 177 echo "<p>Testing linear algebra methods...<p>"; 178 179 $A = new Matrix($columnwise, 3); 180 if( $this->checkMatrices($A->transpose(), $T) ) 181 $this->try_success("Transpose check..."); 182 else 183 $errorCount = $this->try_failure($errorCount, "Transpose check...", "Matrices are not equal"); 184 185 if($this->checkScalars($A->norm1(), $columnsummax)) 186 $this->try_success("Maximum column sum..."); 187 else 188 $errorCount = $this->try_failure($errorCount, "Maximum column sum...", "Incorrect: " . $A->norm1() . " != " . $columnsummax); 189 190 if($this->checkScalars($A->normInf(), $rowsummax)) 191 $this->try_success("Maximum row sum..."); 192 else 193 $errorCount = $this->try_failure($errorCount, "Maximum row sum...", "Incorrect: " . $A->normInf() . " != " . $rowsummax ); 194 195 if($this->checkScalars($A->normF(), sqrt($sumofsquares))) 196 $this->try_success("Frobenius norm..."); 197 else 198 $errorCount = $this->try_failure($errorCount, "Frobenius norm...", "Incorrect:" . $A->normF() . " != " . sqrt($sumofsquares)); 199 200 if($this->checkScalars($A->trace(), $sumofdiagonals)) 201 $this->try_success("Matrix trace..."); 202 else 203 $errorCount = $this->try_failure($errorCount, "Matrix trace...", "Incorrect: " . $A->trace() . " != " . $sumofdiagonals); 204 205 $B = $A->getMatrix(0, $A->getRowDimension(), 0, $A->getRowDimension()); 206 if( $B->det() == 0 ) 207 $this->try_success("Matrix determinant..."); 208 else 209 $errorCount = $this->try_failure($errorCount, "Matrix determinant...", "Incorrect: " . $B->det() . " != " . 0); 210 211 $A = new Matrix($columnwise,3); 212 $SQ = new Matrix($square); 213 if ($this->checkMatrices($SQ, $A->times($A->transpose()))) 214 $this->try_success("times(Matrix)..."); 215 else { 216 $errorCount = $this->try_failure($errorCount, "times(Matrix)...", "Unable to multiply matrices"); 217 $SQ->toHTML(); 218 $AT->toHTML(); 219 } 220 221 $A = new Matrix($columnwise, 4); 222 223 $QR = $A->qr(); 224 $R = $QR->getR(); 225 $Q = $QR->getQ(); 226 if($this->checkMatrices($A, $Q->times($R))) 227 $this->try_success("QRDecomposition...",""); 228 else 229 $errorCount = $this->try_failure($errorCount,"QRDecomposition...","incorrect qr decomposition calculation"); 230 231 $A = new Matrix($columnwise, 4); 232 $SVD = $A->svd(); 233 $U = $SVD->getU(); 234 $S = $SVD->getS(); 235 $V = $SVD->getV(); 236 if ($this->checkMatrices($A, $U->times($S->times($V->transpose())))) 237 $this->try_success("SingularValueDecomposition...",""); 238 else 239 $errorCount = $this->try_failure($errorCount,"SingularValueDecomposition...","incorrect singular value decomposition calculation"); 240 241 $n = $A->getColumnDimension(); 242 $A = $A->getMatrix(0,$n-1,0,$n-1); 243 $A->set(0,0,0.); 244 245 $LU = $A->lu(); 246 $L = $LU->getL(); 247 if ( $this->checkMatrices($A->getMatrix($LU->getPivot(),0,$n-1), $L->times($LU->getU())) ) 248 $this->try_success("LUDecomposition...",""); 249 else 250 $errorCount = $this->try_failure($errorCount,"LUDecomposition...","incorrect LU decomposition calculation"); 251 252 $X = $A->inverse(); 253 if ( $this->checkMatrices($A->times($X),Matrix::identity(3,3)) ) 254 $this->try_success("inverse()...",""); 255 else 256 $errorCount = $this->try_failure($errorCount, "inverse()...","incorrect inverse calculation"); 257 258 $DEF = new Matrix($rankdef); 259 if($this->checkScalars($DEF->rank(), min($DEF->getRowDimension(), $DEF->getColumnDimension())-1)) 260 $this->try_success("Rank..."); 261 else 262 $this->try_failure("Rank...", "incorrect rank calculation"); 263 264 $B = new Matrix($condmat); 265 $SVD = $B->svd(); 266 $singularvalues = $SVD->getSingularValues(); 267 if($this->checkScalars($B->cond(), $singularvalues[0]/$singularvalues[min($B->getRowDimension(), $B->getColumnDimension())-1])) 268 $this->try_success("Condition number..."); 269 else 270 $this->try_failure("Condition number...", "incorrect condition number calculation"); 271 272 $SUB = new Matrix($subavals); 273 $O = new Matrix($SUB->getRowDimension(),1,1.0); 274 $SOL = new Matrix($sqSolution); 275 $SQ = $SUB->getMatrix(0,$SUB->getRowDimension()-1,0,$SUB->getRowDimension()-1); 276 if ( $this->checkMatrices($SQ->solve($SOL),$O) ) 277 $this->try_success("solve()...",""); 278 else 279 $errorCount = $this->try_failure($errorCount,"solve()...","incorrect lu solve calculation"); 280 281 $A = new Matrix($pvals); 282 $Chol = $A->chol(); 283 $L = $Chol->getL(); 284 if ( $this->checkMatrices($A, $L->times($L->transpose())) ) 285 $this->try_success("CholeskyDecomposition...",""); 286 else 287 $errorCount = $this->try_failure($errorCount,"CholeskyDecomposition...","incorrect Cholesky decomposition calculation"); 288 289 $X = $Chol->solve(Matrix::identity(3,3)); 290 if ( $this->checkMatrices($A->times($X), Matrix::identity(3,3)) ) 291 $this->try_success("CholeskyDecomposition solve()...",""); 292 else 293 $errorCount = $this->try_failure($errorCount,"CholeskyDecomposition solve()...","incorrect Choleskydecomposition solve calculation"); 294 295 $Eig = $A->eig(); 296 $D = $Eig->getD(); 297 $V = $Eig->getV(); 298 if( $this->checkMatrices($A->times($V),$V->times($D)) ) 299 $this->try_success("EigenvalueDecomposition (symmetric)...",""); 300 else 301 $errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (symmetric)...","incorrect symmetric Eigenvalue decomposition calculation"); 302 303 $A = new Matrix($evals); 304 $Eig = $A->eig(); 305 $D = $Eig->getD(); 306 $V = $Eig->getV(); 307 if ( $this->checkMatrices($A->times($V),$V->times($D)) ) 308 $this->try_success("EigenvalueDecomposition (nonsymmetric)...",""); 309 else 310 $errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (nonsymmetric)...","incorrect nonsymmetric Eigenvalue decomposition calculation"); 311 312 print("<b>{$errorCount} total errors</b>."); 313 } 314 315 /** 316 * Print appropriate messages for successful outcome try 317 * @param string $s 318 * @param string $e 319 */ 320 function try_success($s, $e = "") { 321 print "> ". $s ."success<br />"; 322 if ($e != "") 323 print "> Message: ". $e ."<br />"; 324 } 325 326 /** 327 * Print appropriate messages for unsuccessful outcome try 328 * @param int $count 329 * @param string $s 330 * @param string $e 331 * @return int incremented counter 332 */ 333 function try_failure($count, $s, $e="") { 334 print "> ". $s ."*** failure ***<br />> Message: ". $e ."<br />"; 335 return ++$count; 336 } 337 338 /** 339 * Print appropriate messages for unsuccessful outcome try 340 * @param int $count 341 * @param string $s 342 * @param string $e 343 * @return int incremented counter 344 */ 345 function try_warning($count, $s, $e="") { 346 print "> ". $s ."*** warning ***<br />> Message: ". $e ."<br />"; 347 return ++$count; 348 } 349 350 /** 351 * Check magnitude of difference of "scalars". 352 * @param float $x 353 * @param float $y 354 */ 355 function checkScalars($x, $y) { 356 $eps = pow(2.0,-52.0); 357 if ($x == 0 & abs($y) < 10*$eps) return; 358 if ($y == 0 & abs($x) < 10*$eps) return; 359 if (abs($x-$y) > 10 * $eps * max(abs($x),abs($y))) 360 return false; 361 else 362 return true; 363 } 364 365 /** 366 * Check norm of difference of "vectors". 367 * @param float $x[] 368 * @param float $y[] 369 */ 370 function checkVectors($x, $y) { 371 $nx = count($x); 372 $ny = count($y); 373 if ($nx == $ny) 374 for($i=0; $i < $nx; ++$i) 375 $this->checkScalars($x[$i],$y[$i]); 376 else 377 die("Attempt to compare vectors of different lengths"); 378 } 379 380 /** 381 * Check norm of difference of "arrays". 382 * @param float $x[][] 383 * @param float $y[][] 384 */ 385 function checkArrays($x, $y) { 386 $A = new Matrix($x); 387 $B = new Matrix($y); 388 return $this->checkMatrices($A,$B); 389 } 390 391 /** 392 * Check norm of difference of "matrices". 393 * @param matrix $X 394 * @param matrix $Y 395 */ 396 function checkMatrices($X = null, $Y = null) { 397 if( $X == null || $Y == null ) 398 return false; 399 400 $eps = pow(2.0,-52.0); 401 if ($X->norm1() == 0. & $Y->norm1() < 10*$eps) return true; 402 if ($Y->norm1() == 0. & $X->norm1() < 10*$eps) return true; 403 404 $A = $X->minus($Y); 405 406 if ($A->norm1() > 1000 * $eps * max($X->norm1(),$Y->norm1())) 407 die("The norm of (X-Y) is too large: ".$A->norm1()); 408 else 409 return true; 410 } 411 412} 413 414$test = new TestMatrix; 415?>