PageRenderTime 32ms CodeModel.GetById 11ms app.highlight 17ms RepoModel.GetById 1ms app.codeStats 0ms

/PHPExcel/Shared/JAMA/tests/TestMatrix.php

https://bitbucket.org/nfredricks/wp-employee-time
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?>