416 lines
14 KiB
PHP
416 lines
14 KiB
PHP
|
<?php
|
||
|
|
|
||
|
|
require_once "../Matrix.php";
|
||
|
|
|
||
|
|
class TestMatrix {
|
||
|
|
|
||
|
|
function TestMatrix() {
|
||
|
|
|
||
|
|
// define test variables
|
||
|
|
|
||
|
|
$errorCount = 0;
|
||
|
|
$warningCount = 0;
|
||
|
|
$columnwise = array(1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.);
|
||
|
|
$rowwise = array(1.,4.,7.,10.,2.,5.,8.,11.,3.,6.,9.,12.);
|
||
|
|
$avals = array(array(1.,4.,7.,10.),array(2.,5.,8.,11.),array(3.,6.,9.,12.));
|
||
|
|
$rankdef = $avals;
|
||
|
|
$tvals = array(array(1.,2.,3.),array(4.,5.,6.),array(7.,8.,9.),array(10.,11.,12.));
|
||
|
|
$subavals = array(array(5.,8.,11.),array(6.,9.,12.));
|
||
|
|
$rvals = array(array(1.,4.,7.),array(2.,5.,8.,11.),array(3.,6.,9.,12.));
|
||
|
|
$pvals = array(array(1.,1.,1.),array(1.,2.,3.),array(1.,3.,6.));
|
||
|
|
$ivals = array(array(1.,0.,0.,0.),array(0.,1.,0.,0.),array(0.,0.,1.,0.));
|
||
|
|
$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.));
|
||
|
|
$square = array(array(166.,188.,210.),array(188.,214.,240.),array(210.,240.,270.));
|
||
|
|
$sqSolution = array(array(13.),array(15.));
|
||
|
|
$condmat = array(array(1.,3.),array(7.,9.));
|
||
|
|
$rows = 3;
|
||
|
|
$cols = 4;
|
||
|
|
$invalidID = 5; /* should trigger bad shape for construction with val */
|
||
|
|
$raggedr = 0; /* (raggedr,raggedc) should be out of bounds in ragged array */
|
||
|
|
$raggedc = 4;
|
||
|
|
$validID = 3; /* leading dimension of intended test Matrices */
|
||
|
|
$nonconformld = 4; /* leading dimension which is valid, but nonconforming */
|
||
|
|
$ib = 1; /* index ranges for sub Matrix */
|
||
|
|
$ie = 2;
|
||
|
|
$jb = 1;
|
||
|
|
$je = 3;
|
||
|
|
$rowindexset = array(1,2);
|
||
|
|
$badrowindexset = array(1,3);
|
||
|
|
$columnindexset = array(1,2,3);
|
||
|
|
$badcolumnindexset = array(1,2,4);
|
||
|
|
$columnsummax = 33.;
|
||
|
|
$rowsummax = 30.;
|
||
|
|
$sumofdiagonals = 15;
|
||
|
|
$sumofsquares = 650;
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Test matrix methods
|
||
|
|
*/
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Constructors and constructor-like methods:
|
||
|
|
*
|
||
|
|
* Matrix(double[], int)
|
||
|
|
* Matrix(double[][])
|
||
|
|
* Matrix(int, int)
|
||
|
|
* Matrix(int, int, double)
|
||
|
|
* Matrix(int, int, double[][])
|
||
|
|
* constructWithCopy(double[][])
|
||
|
|
* random(int,int)
|
||
|
|
* identity(int)
|
||
|
|
*/
|
||
|
|
echo "<p>Testing constructors and constructor-like methods...</p>";
|
||
|
|
|
||
|
|
$A = new Matrix($columnwise, 3);
|
||
|
|
if($A instanceof Matrix) {
|
||
|
|
$this->try_success("Column-packed constructor...");
|
||
|
|
} else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "Column-packed constructor...", "Unable to construct Matrix");
|
||
|
|
|
||
|
|
$T = new Matrix($tvals);
|
||
|
|
if($T instanceof Matrix)
|
||
|
|
$this->try_success("2D array constructor...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "2D array constructor...", "Unable to construct Matrix");
|
||
|
|
|
||
|
|
$A = new Matrix($columnwise, $validID);
|
||
|
|
$B = new Matrix($avals);
|
||
|
|
$tmp = $B->get(0,0);
|
||
|
|
$avals[0][0] = 0.0;
|
||
|
|
$C = $B->minus($A);
|
||
|
|
$avals[0][0] = $tmp;
|
||
|
|
$B = Matrix::constructWithCopy($avals);
|
||
|
|
$tmp = $B->get(0,0);
|
||
|
|
$avals[0][0] = 0.0;
|
||
|
|
/** check that constructWithCopy behaves properly **/
|
||
|
|
if ( ( $tmp - $B->get(0,0) ) != 0.0 )
|
||
|
|
$errorCount = $this->try_failure($errorCount,"constructWithCopy... ","copy not effected... data visible outside");
|
||
|
|
else
|
||
|
|
$this->try_success("constructWithCopy... ","");
|
||
|
|
|
||
|
|
$I = new Matrix($ivals);
|
||
|
|
if ( $this->checkMatrices($I,Matrix::identity(3,4)) )
|
||
|
|
$this->try_success("identity... ","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"identity... ","identity Matrix not successfully created");
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Access Methods:
|
||
|
|
*
|
||
|
|
* getColumnDimension()
|
||
|
|
* getRowDimension()
|
||
|
|
* getArray()
|
||
|
|
* getArrayCopy()
|
||
|
|
* getColumnPackedCopy()
|
||
|
|
* getRowPackedCopy()
|
||
|
|
* get(int,int)
|
||
|
|
* getMatrix(int,int,int,int)
|
||
|
|
* getMatrix(int,int,int[])
|
||
|
|
* getMatrix(int[],int,int)
|
||
|
|
* getMatrix(int[],int[])
|
||
|
|
* set(int,int,double)
|
||
|
|
* setMatrix(int,int,int,int,Matrix)
|
||
|
|
* setMatrix(int,int,int[],Matrix)
|
||
|
|
* setMatrix(int[],int,int,Matrix)
|
||
|
|
* setMatrix(int[],int[],Matrix)
|
||
|
|
*/
|
||
|
|
print "<p>Testing access methods...</p>";
|
||
|
|
|
||
|
|
$B = new Matrix($avals);
|
||
|
|
if($B->getRowDimension() == $rows)
|
||
|
|
$this->try_success("getRowDimension...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "getRowDimension...");
|
||
|
|
|
||
|
|
if($B->getColumnDimension() == $cols)
|
||
|
|
$this->try_success("getColumnDimension...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "getColumnDimension...");
|
||
|
|
|
||
|
|
$barray = $B->getArray();
|
||
|
|
if($this->checkArrays($barray, $avals))
|
||
|
|
$this->try_success("getArray...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "getArray...");
|
||
|
|
|
||
|
|
$bpacked = $B->getColumnPackedCopy();
|
||
|
|
if($this->checkArrays($bpacked, $columnwise))
|
||
|
|
$this->try_success("getColumnPackedCopy...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "getColumnPackedCopy...");
|
||
|
|
|
||
|
|
$bpacked = $B->getRowPackedCopy();
|
||
|
|
if($this->checkArrays($bpacked, $rowwise))
|
||
|
|
$this->try_success("getRowPackedCopy...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "getRowPackedCopy...");
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Array-like methods:
|
||
|
|
* minus
|
||
|
|
* minusEquals
|
||
|
|
* plus
|
||
|
|
* plusEquals
|
||
|
|
* arrayLeftDivide
|
||
|
|
* arrayLeftDivideEquals
|
||
|
|
* arrayRightDivide
|
||
|
|
* arrayRightDivideEquals
|
||
|
|
* arrayTimes
|
||
|
|
* arrayTimesEquals
|
||
|
|
* uminus
|
||
|
|
*/
|
||
|
|
print "<p>Testing array-like methods...</p>";
|
||
|
|
|
||
|
|
/**
|
||
|
|
* I/O methods:
|
||
|
|
* read
|
||
|
|
* print
|
||
|
|
* serializable:
|
||
|
|
* writeObject
|
||
|
|
* readObject
|
||
|
|
*/
|
||
|
|
print "<p>Testing I/O methods...</p>";
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Test linear algebra methods
|
||
|
|
*/
|
||
|
|
echo "<p>Testing linear algebra methods...<p>";
|
||
|
|
|
||
|
|
$A = new Matrix($columnwise, 3);
|
||
|
|
if( $this->checkMatrices($A->transpose(), $T) )
|
||
|
|
$this->try_success("Transpose check...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "Transpose check...", "Matrices are not equal");
|
||
|
|
|
||
|
|
if($this->checkScalars($A->norm1(), $columnsummax))
|
||
|
|
$this->try_success("Maximum column sum...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "Maximum column sum...", "Incorrect: " . $A->norm1() . " != " . $columnsummax);
|
||
|
|
|
||
|
|
if($this->checkScalars($A->normInf(), $rowsummax))
|
||
|
|
$this->try_success("Maximum row sum...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "Maximum row sum...", "Incorrect: " . $A->normInf() . " != " . $rowsummax );
|
||
|
|
|
||
|
|
if($this->checkScalars($A->normF(), sqrt($sumofsquares)))
|
||
|
|
$this->try_success("Frobenius norm...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "Frobenius norm...", "Incorrect:" . $A->normF() . " != " . sqrt($sumofsquares));
|
||
|
|
|
||
|
|
if($this->checkScalars($A->trace(), $sumofdiagonals))
|
||
|
|
$this->try_success("Matrix trace...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "Matrix trace...", "Incorrect: " . $A->trace() . " != " . $sumofdiagonals);
|
||
|
|
|
||
|
|
$B = $A->getMatrix(0, $A->getRowDimension(), 0, $A->getRowDimension());
|
||
|
|
if( $B->det() == 0 )
|
||
|
|
$this->try_success("Matrix determinant...");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "Matrix determinant...", "Incorrect: " . $B->det() . " != " . 0);
|
||
|
|
|
||
|
|
$A = new Matrix($columnwise,3);
|
||
|
|
$SQ = new Matrix($square);
|
||
|
|
if ($this->checkMatrices($SQ, $A->times($A->transpose())))
|
||
|
|
$this->try_success("times(Matrix)...");
|
||
|
|
else {
|
||
|
|
$errorCount = $this->try_failure($errorCount, "times(Matrix)...", "Unable to multiply matrices");
|
||
|
|
$SQ->toHTML();
|
||
|
|
$AT->toHTML();
|
||
|
|
}
|
||
|
|
|
||
|
|
$A = new Matrix($columnwise, 4);
|
||
|
|
|
||
|
|
$QR = $A->qr();
|
||
|
|
$R = $QR->getR();
|
||
|
|
$Q = $QR->getQ();
|
||
|
|
if($this->checkMatrices($A, $Q->times($R)))
|
||
|
|
$this->try_success("QRDecomposition...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"QRDecomposition...","incorrect qr decomposition calculation");
|
||
|
|
|
||
|
|
$A = new Matrix($columnwise, 4);
|
||
|
|
$SVD = $A->svd();
|
||
|
|
$U = $SVD->getU();
|
||
|
|
$S = $SVD->getS();
|
||
|
|
$V = $SVD->getV();
|
||
|
|
if ($this->checkMatrices($A, $U->times($S->times($V->transpose()))))
|
||
|
|
$this->try_success("SingularValueDecomposition...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"SingularValueDecomposition...","incorrect singular value decomposition calculation");
|
||
|
|
|
||
|
|
$n = $A->getColumnDimension();
|
||
|
|
$A = $A->getMatrix(0,$n-1,0,$n-1);
|
||
|
|
$A->set(0,0,0.);
|
||
|
|
|
||
|
|
$LU = $A->lu();
|
||
|
|
$L = $LU->getL();
|
||
|
|
if ( $this->checkMatrices($A->getMatrix($LU->getPivot(),0,$n-1), $L->times($LU->getU())) )
|
||
|
|
$this->try_success("LUDecomposition...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"LUDecomposition...","incorrect LU decomposition calculation");
|
||
|
|
|
||
|
|
$X = $A->inverse();
|
||
|
|
if ( $this->checkMatrices($A->times($X),Matrix::identity(3,3)) )
|
||
|
|
$this->try_success("inverse()...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount, "inverse()...","incorrect inverse calculation");
|
||
|
|
|
||
|
|
$DEF = new Matrix($rankdef);
|
||
|
|
if($this->checkScalars($DEF->rank(), min($DEF->getRowDimension(), $DEF->getColumnDimension())-1))
|
||
|
|
$this->try_success("Rank...");
|
||
|
|
else
|
||
|
|
$this->try_failure("Rank...", "incorrect rank calculation");
|
||
|
|
|
||
|
|
$B = new Matrix($condmat);
|
||
|
|
$SVD = $B->svd();
|
||
|
|
$singularvalues = $SVD->getSingularValues();
|
||
|
|
if($this->checkScalars($B->cond(), $singularvalues[0]/$singularvalues[min($B->getRowDimension(), $B->getColumnDimension())-1]))
|
||
|
|
$this->try_success("Condition number...");
|
||
|
|
else
|
||
|
|
$this->try_failure("Condition number...", "incorrect condition number calculation");
|
||
|
|
|
||
|
|
$SUB = new Matrix($subavals);
|
||
|
|
$O = new Matrix($SUB->getRowDimension(),1,1.0);
|
||
|
|
$SOL = new Matrix($sqSolution);
|
||
|
|
$SQ = $SUB->getMatrix(0,$SUB->getRowDimension()-1,0,$SUB->getRowDimension()-1);
|
||
|
|
if ( $this->checkMatrices($SQ->solve($SOL),$O) )
|
||
|
|
$this->try_success("solve()...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"solve()...","incorrect lu solve calculation");
|
||
|
|
|
||
|
|
$A = new Matrix($pvals);
|
||
|
|
$Chol = $A->chol();
|
||
|
|
$L = $Chol->getL();
|
||
|
|
if ( $this->checkMatrices($A, $L->times($L->transpose())) )
|
||
|
|
$this->try_success("CholeskyDecomposition...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"CholeskyDecomposition...","incorrect Cholesky decomposition calculation");
|
||
|
|
|
||
|
|
$X = $Chol->solve(Matrix::identity(3,3));
|
||
|
|
if ( $this->checkMatrices($A->times($X), Matrix::identity(3,3)) )
|
||
|
|
$this->try_success("CholeskyDecomposition solve()...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"CholeskyDecomposition solve()...","incorrect Choleskydecomposition solve calculation");
|
||
|
|
|
||
|
|
$Eig = $A->eig();
|
||
|
|
$D = $Eig->getD();
|
||
|
|
$V = $Eig->getV();
|
||
|
|
if( $this->checkMatrices($A->times($V),$V->times($D)) )
|
||
|
|
$this->try_success("EigenvalueDecomposition (symmetric)...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (symmetric)...","incorrect symmetric Eigenvalue decomposition calculation");
|
||
|
|
|
||
|
|
$A = new Matrix($evals);
|
||
|
|
$Eig = $A->eig();
|
||
|
|
$D = $Eig->getD();
|
||
|
|
$V = $Eig->getV();
|
||
|
|
if ( $this->checkMatrices($A->times($V),$V->times($D)) )
|
||
|
|
$this->try_success("EigenvalueDecomposition (nonsymmetric)...","");
|
||
|
|
else
|
||
|
|
$errorCount = $this->try_failure($errorCount,"EigenvalueDecomposition (nonsymmetric)...","incorrect nonsymmetric Eigenvalue decomposition calculation");
|
||
|
|
|
||
|
|
print("<b>{$errorCount} total errors</b>.");
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Print appropriate messages for successful outcome try
|
||
|
|
* @param string $s
|
||
|
|
* @param string $e
|
||
|
|
*/
|
||
|
|
function try_success($s, $e = "") {
|
||
|
|
print "> ". $s ."success<br />";
|
||
|
|
if ($e != "")
|
||
|
|
print "> Message: ". $e ."<br />";
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Print appropriate messages for unsuccessful outcome try
|
||
|
|
* @param int $count
|
||
|
|
* @param string $s
|
||
|
|
* @param string $e
|
||
|
|
* @return int incremented counter
|
||
|
|
*/
|
||
|
|
function try_failure($count, $s, $e="") {
|
||
|
|
print "> ". $s ."*** failure ***<br />> Message: ". $e ."<br />";
|
||
|
|
return ++$count;
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Print appropriate messages for unsuccessful outcome try
|
||
|
|
* @param int $count
|
||
|
|
* @param string $s
|
||
|
|
* @param string $e
|
||
|
|
* @return int incremented counter
|
||
|
|
*/
|
||
|
|
function try_warning($count, $s, $e="") {
|
||
|
|
print "> ". $s ."*** warning ***<br />> Message: ". $e ."<br />";
|
||
|
|
return ++$count;
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Check magnitude of difference of "scalars".
|
||
|
|
* @param float $x
|
||
|
|
* @param float $y
|
||
|
|
*/
|
||
|
|
function checkScalars($x, $y) {
|
||
|
|
$eps = pow(2.0,-52.0);
|
||
|
|
if ($x == 0 & abs($y) < 10*$eps) return;
|
||
|
|
if ($y == 0 & abs($x) < 10*$eps) return;
|
||
|
|
if (abs($x-$y) > 10 * $eps * max(abs($x),abs($y)))
|
||
|
|
return false;
|
||
|
|
else
|
||
|
|
return true;
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Check norm of difference of "vectors".
|
||
|
|
* @param float $x[]
|
||
|
|
* @param float $y[]
|
||
|
|
*/
|
||
|
|
function checkVectors($x, $y) {
|
||
|
|
$nx = count($x);
|
||
|
|
$ny = count($y);
|
||
|
|
if ($nx == $ny)
|
||
|
|
for($i=0; $i < $nx; ++$i)
|
||
|
|
$this->checkScalars($x[$i],$y[$i]);
|
||
|
|
else
|
||
|
|
die("Attempt to compare vectors of different lengths");
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Check norm of difference of "arrays".
|
||
|
|
* @param float $x[][]
|
||
|
|
* @param float $y[][]
|
||
|
|
*/
|
||
|
|
function checkArrays($x, $y) {
|
||
|
|
$A = new Matrix($x);
|
||
|
|
$B = new Matrix($y);
|
||
|
|
return $this->checkMatrices($A,$B);
|
||
|
|
}
|
||
|
|
|
||
|
|
/**
|
||
|
|
* Check norm of difference of "matrices".
|
||
|
|
* @param matrix $X
|
||
|
|
* @param matrix $Y
|
||
|
|
*/
|
||
|
|
function checkMatrices($X = null, $Y = null) {
|
||
|
|
if( $X == null || $Y == null )
|
||
|
|
return false;
|
||
|
|
|
||
|
|
$eps = pow(2.0,-52.0);
|
||
|
|
if ($X->norm1() == 0. & $Y->norm1() < 10*$eps) return true;
|
||
|
|
if ($Y->norm1() == 0. & $X->norm1() < 10*$eps) return true;
|
||
|
|
|
||
|
|
$A = $X->minus($Y);
|
||
|
|
|
||
|
|
if ($A->norm1() > 1000 * $eps * max($X->norm1(),$Y->norm1()))
|
||
|
|
die("The norm of (X-Y) is too large: ".$A->norm1());
|
||
|
|
else
|
||
|
|
return true;
|
||
|
|
}
|
||
|
|
|
||
|
|
}
|
||
|
|
|
||
|
|
$test = new TestMatrix;
|
||
|
|
?>
|