Test Network, Test Changes, Algorithm Updates

Author: allywarner

src/Core/Algorithms/Math/ComputePCA.cc

@@ -21,67 +21,86 @@
  DEALINGS IN THE SOFTWARE.
  */
 
+
+//ComputePCA Algorithm: Computes Principal Component Analysis.
+
 #include <Core/Algorithms/Base/AlgorithmPreconditions.h>
 #include <Core/Algorithms/Math/ComputePCA.h>
 #include <Core/Datatypes/DenseMatrix.h>
 #include <Core/Datatypes/DenseColumnMatrix.h>
 #include <Core/Datatypes/MatrixTypeConversions.h>
 #include <Eigen/SVD>
-
 #include <Core/Algorithms/Base/AlgorithmVariableNames.h>
 
 using namespace SCIRun;
 using namespace SCIRun::Core::Algorithms;
 using namespace SCIRun::Core::Datatypes;
 using namespace SCIRun::Core::Algorithms::Math;
 
-void ComputePCAAlgo::run(MatrixHandle input, DenseMatrixHandle& LeftPrinMat, DenseMatrixHandle& PrinVals, DenseMatrixHandle& RightPrinMat) const
-{
+//Let's do some math.
+//Algorithm:
+void ComputePCAAlgo::run(MatrixHandle input, DenseMatrixHandle& LeftPrinMat, DenseMatrixHandle& PrinVals, DenseMatrixHandle& RightPrinMat) const{
+    
+    //Throws an error if one or both of the input matrix dimensions is zero.
+    if (input->nrows() == 0 || input->ncols() == 0){
+    
+        THROW_ALGORITHM_INPUT_ERROR("Input has a zero dimension.");
+    }
+    
+    //Input matrix: nxm
     if (matrix_is::dense(input))
     {
+        //First, we have to center the data.
+        auto denseInputCentered = centerData(input);
 
-        auto denseInput = matrix_cast::as_dense(input);
-        
-        auto rows = denseInput->rows();
-        
-        auto centerMatrix = Eigen::MatrixXd::Identity(rows,rows) - (1/rows)*Eigen::MatrixXd::Constant(rows,rows,1);
-        
-        auto denseInputCentered = centerMatrix * *denseInput;
-        
+        //After the data is centered, then we compute SVD on the centered matrix.
+        //Centered Matrix = U*S*Vt, Vt = V transpose
         Eigen::JacobiSVD<DenseMatrix::EigenBase> svd_mat(denseInputCentered, Eigen::ComputeFullU | Eigen::ComputeFullV);
 
+        //U: Left principal matrix, nxn, orthogonal
         LeftPrinMat = boost::make_shared<DenseMatrix>(svd_mat.matrixU());
 
+        //S: Principal values nxm, diagonal
         PrinVals = boost::make_shared<DenseMatrix>(svd_mat.singularValues());
 
+        //V: Right singular mxm, orthognol
         RightPrinMat = boost::make_shared<DenseMatrix>(svd_mat.matrixV());
     }
     else
     {
+        //Throw an error if the matrix is not dense.
+        //Sparse matrices not supported at this time.
         THROW_ALGORITHM_INPUT_ERROR("ComputePCA works for dense matrix input only.");
     }
 }
 
-DenseMatrix ComputePCAAlgo::centerData(DenseMatrixHandle input_matrix)
+//Centers input matrix.
+DenseMatrix ComputePCAAlgo::centerData(MatrixHandle input_matrix)
 {
+    //Casts the matrix as dense.
     auto denseInput = matrix_cast::as_dense(input_matrix);
 
+    //Counts the number of rows in the input matrix.
     auto rows = denseInput->rows();
 
+    //Calulates the centering matrix (C).
+    // C = Identity(nxn) - 1/n * matrix of ones(nxn)
     auto centerMatrix = Eigen::MatrixXd::Identity(rows,rows) - (1.0/rows)*Eigen::MatrixXd::Constant(rows,rows,1);
 
+    //Multiplying the input matrix by the centering matrix.
     auto denseInputCentered = centerMatrix * *denseInput;
-        
-    return denseInputCentered.eval();
+    
+    return denseInputCentered;
 }
 
+//Run the algorithm.
 AlgorithmOutput ComputePCAAlgo::run_generic(const AlgorithmInput& input) const
 {
     auto input_matrix = input.get<Matrix>(Variables::InputMatrix);
 
     DenseMatrixHandle LeftPrinMat;
-    DenseMatrixHandle RightPrinMat;
     DenseMatrixHandle PrinVals;
+    DenseMatrixHandle RightPrinMat;
 
     run(input_matrix, LeftPrinMat, PrinVals, RightPrinMat);
 
@@ -94,6 +113,7 @@ AlgorithmOutput ComputePCAAlgo::run_generic(const AlgorithmInput& input) const
     return output;
 }
 
+//Outputs:
 AlgorithmOutputName ComputePCAAlgo::LeftPrincipalMatrix("LeftPrincipalMatrix");
 AlgorithmOutputName ComputePCAAlgo::PrincipalValues("PrincipalValues");
 AlgorithmOutputName ComputePCAAlgo::RightPrincipalMatrix("RightPrincipalMatrix");

src/Core/Algorithms/Math/ComputePCA.h

@@ -21,6 +21,7 @@
  DEALINGS IN THE SOFTWARE.
  */
 
+//ComputePCA Algortithm Header
 #ifndef ComputePCA_ComputePCA_h_1
 #define ComputePCA_ComputePCA_h_1
 
@@ -43,7 +44,7 @@ namespace SCIRun {
                     static AlgorithmOutputName RightPrincipalMatrix;
                     void run(Datatypes::MatrixHandle input_matrix, Datatypes::DenseMatrixHandle& LeftPrinMat, Datatypes::DenseMatrixHandle& PrinVals, Datatypes::DenseMatrixHandle& RightPrinMat) const;
                     virtual AlgorithmOutput run_generic(const AlgorithmInput& input) const;
-                    static Datatypes::DenseMatrix centerData(Datatypes::DenseMatrixHandle input_matrix);
+                    static Datatypes::DenseMatrix centerData(Datatypes::MatrixHandle input_matrix);
                 };
 
             }}}}

src/Core/Algorithms/Math/Tests/ComputePCAtest.cc

@@ -27,25 +27,26 @@
  */
 
 #include <gtest/gtest.h>
-
 #include <Core/Datatypes/DenseMatrix.h>
 #include <Core/Datatypes/MatrixComparison.h>
 #include <Testing/Utils/MatrixTestUtilities.h>
 #include <Core/Algorithms/Math/ComputePCA.h>
+#include <Eigen/SVD>
 
 using namespace SCIRun::Core::Datatypes;
 using namespace SCIRun::Core::Algorithms::Math;
 using namespace SCIRun::TestUtils;
 
-
 namespace
 {
-    //Matrix for input
+    //Matrix for input.
     DenseMatrixHandle inputMatrix()
     {
+        //Hard coded for testing.
         int column1 [12] = {5,4,7,8,10,10,4,6,8,7,9,6};
         int column2 [12] = {7,6,9,8,5,7,9,3,4,5,5,9};
 
+        //Puts the columns into a matrix.
         DenseMatrixHandle inputM(boost::make_shared<DenseMatrix>(12,2));
         for (int i = 0; i < inputM->rows(); i++){
             (*inputM)(i,0) = column1[i];
@@ -54,34 +55,24 @@ namespace
         return inputM;
     }
 
-    //A after it is centered
+    //The input matrix after centering.
     DenseMatrixHandle centeredInputMatrix()
     {
-    
+        //Hard coded for testing.
         double centeredColumn1 [12] = {-2.0,-3.0,-0.0,1.0,3.0,3.0,-3.0,-1.0,1.0,-0.0,2.0,-1.0};
         double centeredColumn2 [12] = {0.583333,-0.416667,2.583333,1.583333,-1.416667,0.583333,2.583333,-3.416667,-2.416667,-1.416667,-1.416667,2.583333};
 
+        //Puts the columns into a matrix.
         DenseMatrixHandle centeredM(boost::make_shared<DenseMatrix>(12,2));
         for (int i = 0; i < centeredM->rows(); i++){
             (*centeredM)(i,0) = centeredColumn1[i];
             (*centeredM)(i,1) = centeredColumn2[i];
         }
         return centeredM;
-                                
     }
-                            
-    //Outputs: U,S,V
-    DenseMatrixHandle outputU()
-    {}
-    DenseMatrixHandle outputS()
-    {}
-    DenseMatrixHandle outputV()
-    {}
-    
 }
 
-
-
+//Checks if the algorithm centers the data correctly.
 TEST(ComputePCAtest, CenterData)
 {
     ComputePCAAlgo algo;
@@ -92,33 +83,77 @@ TEST(ComputePCAtest, CenterData)
 
     auto expected = *centeredInputMatrix();
 
+    //Checking every element with a tolerance of 1e-5.
     for (int i = 0; i < centered.rows(); ++i) {
         for (int j = 0; j < centered.cols(); ++j)
-            EXPECT_NEAR(expected(i,j), centered(i,j), 1e-5);
+            ASSERT_NEAR(expected(i,j), centered(i,j), 1e-5);
+    }
+}
+
+//Checks if the outputs are correct.
+//U,S,V: U: Left principal matrix, S: principal values, V: Right principal matrix.
+TEST(ComputePCAtest, checkOutputs)
+{
+    ComputePCAAlgo algo;
+    
+    DenseMatrixHandle m1(inputMatrix());
+    DenseMatrixHandle LeftPrinMat_U;
+    DenseMatrixHandle PrinVals_S;
+    DenseMatrixHandle RightPrinMat_V;
+    
+    //Runs algorithm.
+    algo.run(m1,LeftPrinMat_U,PrinVals_S,RightPrinMat_V);
+    
+    //Testing if the results are null.
+    ASSERT_NE(nullptr,LeftPrinMat_U);
+    ASSERT_NE(nullptr,PrinVals_S);
+    ASSERT_NE(nullptr,RightPrinMat_V);
+    
+    //Check the dimensions of the matrices that were created for output.
+    
+    //Rows
+    ASSERT_EQ(12,LeftPrinMat_U->rows());
+    ASSERT_EQ(2,PrinVals_S->rows());
+    ASSERT_EQ(2,RightPrinMat_V->rows());
+    
+    //Columns
+    ASSERT_EQ(12,LeftPrinMat_U->cols());
+    ASSERT_EQ(1,PrinVals_S->cols());
+    ASSERT_EQ(2,RightPrinMat_V->cols());
+    
+    //Eigen does not create a diagonal matrix when it computes SVD, it just has a column with the principal values, so we must put it into a diagonal matrix to be able to do some matrix multiplication later.
+    DenseMatrix sDiag = Eigen::MatrixXd::Constant(12,2,0);
+    sDiag.diagonal() = PrinVals_S->col(0);
+    
+    //Multiplying back together and comparing to the centered matrix. They should be equal to each other with some tolerance.
+    auto product = (*LeftPrinMat_U) * sDiag * (*RightPrinMat_V).transpose();
+    
+    auto expected = *centeredInputMatrix();
+    
+    //Comparing each element in the matrices.
+    for (int i = 0; i < product.rows(); ++i) {
+        for (int j = 0; j < product.cols(); ++j)
+            ASSERT_NEAR(expected(i,j), product(i,j), 1e-5);
     }
+
 }
 
-//TEST(ComputePCAtest, output1)
-//{
-//    ComputePCAAlgorithm algo;
-//    
-//    DenseMatrixHandle m1(testMatrix());
-//    ComputePCAAlgorithm::Outputs result = algo.run(ComputePCAAlgorithm::Inputs(m1),ComputePCAAlgorithm);
-//    
-//}
-//
-//TEST(ComputePCAtest, output2)
-//{
-//    ComputePCAAlgorithm algo;
-//    
-//    DenseMatrixHandle m1(testMatrix());
-//    
-//}
-//
-//TEST(ComputePCAtest, output3)
-//{
-//    ComputePCAAlgorithm algo;
-//    
-//    DenseMatrixHandle m1(testMatrix());
-//    
-//}
+//Tests for input with a dimension of zero.
+TEST(ComputePCAtest, ThrowsForZeroDimensionInput)
+{
+    ComputePCAAlgo algo;
+    
+    DenseMatrixHandle m1(new DenseMatrix(5, 0));
+    DenseMatrixHandle m2(new DenseMatrix(0, 5));
+    DenseMatrixHandle m3(new DenseMatrix(0, 0));
+    
+    DenseMatrixHandle LeftPrinMat_U;
+    DenseMatrixHandle PrinVals_S;
+    DenseMatrixHandle RightPrinMat_V;
+    
+    //Runs algorithm and expects an error.
+    EXPECT_ANY_THROW(algo.run(m1,LeftPrinMat_U,PrinVals_S,RightPrinMat_V));
+    EXPECT_ANY_THROW(algo.run(m2,LeftPrinMat_U,PrinVals_S,RightPrinMat_V));
+    EXPECT_ANY_THROW(algo.run(m3,LeftPrinMat_U,PrinVals_S,RightPrinMat_V));
+
+}