text: /**
   * @file tests/pca_test.cpp
   * @author Ajinkya Kale
   * @author Marcus Edel
   *
   * Test file for PCA class.
   *
   * mlpack is free software; you may redistribute it and/or modify it under the
   * terms of the 3-clause BSD license.  You should have received a copy of the
   * 3-clause BSD license along with mlpack.  If not, see
   * http://www.opensource.org/licenses/BSD-3-Clause for more information.
   */
  #include <mlpack/core.hpp>
  #include <mlpack/methods/pca/pca.hpp>
  #include <mlpack/methods/pca/decomposition_policies/exact_svd_method.hpp>
  #include <mlpack/methods/pca/decomposition_policies/quic_svd_method.hpp>
  #include <mlpack/methods/pca/decomposition_policies/randomized_svd_method.hpp>
  #include <mlpack/methods/pca/decomposition_policies/randomized_block_krylov_method.hpp>
  
  #include <boost/test/unit_test.hpp>
  #include "test_tools.hpp"
  
  BOOST_AUTO_TEST_SUITE(PCATest);
  
  using namespace arma;
  using namespace mlpack;
  using namespace mlpack::pca;
  using namespace mlpack::distribution;
  
  /*
   * Compare the output of the our PCA implementation with Armadillo's using the
   * specified decomposition policy.
   */
  template<typename DecompositionPolicy>
  void ArmaComparisonPCA(
      const bool scaleData = false,
      const DecompositionPolicy& decomposition = DecompositionPolicy())
  {
    arma::mat coeff, coeff1, score, score1;
    arma::vec eigVal, eigVal1;
  
    arma::mat data = arma::randu<arma::mat>(3, 1000);
  
    PCA<DecompositionPolicy> pcaType(scaleData, decomposition);
    pcaType.Apply(data, score1, eigVal1, coeff1);
  
    princomp(coeff, score, eigVal, trans(data));
  
    // Verify the PCA results based on the eigenvalues.
    for (size_t i = 0; i < eigVal.n_elem; ++i)
    {
      if (eigVal[i] == 0.0)
        BOOST_REQUIRE_SMALL(eigVal1[i], 1e-15);
      else
        BOOST_REQUIRE_CLOSE(eigVal[i], eigVal1[i], 0.0001);
    }
  }
  
  /*
   * Test that dimensionality reduction with PCA works the same way MATLAB does
   * (which should be correct!) using the specified decomposition policy.
   */
  template<typename DecompositionPolicy>
  void PCADimensionalityReduction(
      const bool scaleData = false,
      const DecompositionPolicy& decomposition = DecompositionPolicy())
  {
    // Fake, simple dataset.  The results we will compare against are from MATLAB.
    mat data("1 0 2 3 9;"
             "5 2 8 4 8;"
             "6 7 3 1 8");
  
    // Now run PCA to reduce the dimensionality.
    size_t trial = 0;
    bool success = false;
    double varRetained = 0.0;
    while (trial < 3 && !success)
    {
      // In some cases the LU decomposition may fail.
      try
      {
        PCA<DecompositionPolicy> p(scaleData, decomposition);
        varRetained = p.Apply(data, 2); // Reduce to 2 dimensions.
        success = true;
      }
      catch (std::logic_error&) { }
  
      ++trial;
    }
  
    BOOST_REQUIRE_EQUAL(success, true);
  
    // Compare with correct results.
    mat correct("-1.53781086 -3.51358020 -0.16139887 -1.87706634  7.08985628;"
                " 1.29937798  3.45762685 -2.69910005 -3.15620704  1.09830225");
  
    BOOST_REQUIRE_EQUAL(data.n_rows, correct.n_rows);
    BOOST_REQUIRE_EQUAL(data.n_cols, correct.n_cols);
  
    // If the eigenvectors are pointed opposite directions, they will cancel
    // each other out in this summation.
    for (size_t i = 0; i < data.n_rows; ++i)
    {
      if (accu(abs(correct.row(i) + data.row(i))) < 0.001 /* arbitrary */)
      {
        // Flip Armadillo coefficients for this column.
        data.row(i) *= -1;
      }
    }
  
    for (size_t row = 0; row < 2; row++)
      for (size_t col = 0; col < 5; col++)
        BOOST_REQUIRE_CLOSE(data(row, col), correct(row, col), 1e-3);
  
    // Check that the amount of variance retained is right.
    BOOST_REQUIRE_CLOSE(varRetained, 0.904876047045906, 1e-5);
  }
  
  /**
   * Test that setting the variance retained parameter to perform dimensionality
   * reduction works using the specified decomposition policy.
   */
  template<typename DecompositionPolicy>
  void PCAVarianceRetained()
  {
      // Fake, simple dataset.
    mat data("1 0 2 3 9;"
             "5 2 8 4 8;"
             "6 7 3 1 8");
  
    // The normalized eigenvalues:
    //   0.616237391936100
    //   0.288638655109805
    //   0.095123952954094
    // So if we keep one dimension, the actual variance retained is
    //   0.616237391936100
    // and if we keep two, the actual variance retained is
    //   0.904876047045906
    // and if we keep three, the actual variance retained is 1.
    PCA<DecompositionPolicy> p;
    arma::mat origData = data;
    double varRetained = p.Apply(data, 0.1);
  
    BOOST_REQUIRE_EQUAL(data.n_rows, 1);
    BOOST_REQUIRE_EQUAL(data.n_cols, 5);
    BOOST_REQUIRE_CLOSE(varRetained, 0.616237391936100, 1e-5);
  
    data = origData;
    varRetained = p.Apply(data, 0.5);
  
    BOOST_REQUIRE_EQUAL(data.n_rows, 1);
    BOOST_REQUIRE_EQUAL(data.n_cols, 5);
    BOOST_REQUIRE_CLOSE(varRetained, 0.616237391936100, 1e-5);
  
    data = origData;
    varRetained = p.Apply(data, 0.7);
  
    BOOST_REQUIRE_EQUAL(data.n_rows, 2);
    BOOST_REQUIRE_EQUAL(data.n_cols, 5);
    BOOST_REQUIRE_CLOSE(varRetained, 0.904876047045906, 1e-5);
  
    data = origData;
    varRetained = p.Apply(data, 0.904);
  
    BOOST_REQUIRE_EQUAL(data.n_rows, 2);
    BOOST_REQUIRE_EQUAL(data.n_cols, 5);
    BOOST_REQUIRE_CLOSE(varRetained, 0.904876047045906, 1e-5);
  
    data = origData;
    varRetained = p.Apply(data, 0.905);
  
    BOOST_REQUIRE_EQUAL(data.n_rows, 3);
    BOOST_REQUIRE_EQUAL(data.n_cols, 5);
    BOOST_REQUIRE_CLOSE(varRetained, 1.0, 1e-5);
  
    data = origData;
    varRetained = p.Apply(data, 1.0);
  
    BOOST_REQUIRE_EQUAL(data.n_rows, 3);
    BOOST_REQUIRE_EQUAL(data.n_cols, 5);
    BOOST_REQUIRE_CLOSE(varRetained, 1.0, 1e-5);
  }
  
  /**
   * Compare the output of our exact PCA implementation with Armadillo's.
   */
  BOOST_AUTO_TEST_CASE(ArmaComparisonExactPCATest)
  {
    ArmaComparisonPCA<ExactSVDPolicy>();
  }
  
  /**
   * Compare the output of our randomized block krylov PCA implementation with
   * Armadillo's.
   */
  BOOST_AUTO_TEST_CASE(ArmaComparisonRandomizedBlockKrylovPCATest)
  {
    RandomizedBlockKrylovSVDPolicy decomposition(5);
    ArmaComparisonPCA<RandomizedBlockKrylovSVDPolicy>(false, decomposition);
  }
  
  /**
   * Compare the output of our randomized-SVD PCA implementation with Armadillo's.
   */
  BOOST_AUTO_TEST_CASE(ArmaComparisonRandomizedPCATest)
  {
    ArmaComparisonPCA<RandomizedSVDPolicy>();
  }
  
  /**
   * Test that dimensionality reduction with exact-svd PCA works the same way
   * MATLAB does (which should be correct!).
   */
  BOOST_AUTO_TEST_CASE(ExactPCADimensionalityReductionTest)
  {
    PCADimensionalityReduction<ExactSVDPolicy>();
  }
  
  /**
   * Test that dimensionality reduction with randomized block krylov PCA works the
   * same way MATLAB does (which should be correct!).
   */
  BOOST_AUTO_TEST_CASE(RandomizedBlockKrylovPCADimensionalityReductionTest)
  {
    RandomizedBlockKrylovSVDPolicy decomposition(5);
    PCADimensionalityReduction<RandomizedBlockKrylovSVDPolicy>(false,
        decomposition);
  }
  
  /**
   * Test that dimensionality reduction with randomized-svd PCA works the same way
   * MATLAB does (which should be correct!).
   */
  BOOST_AUTO_TEST_CASE(RandomizedPCADimensionalityReductionTest)
  {
    PCADimensionalityReduction<RandomizedSVDPolicy>();
  }
  
  /**
   * Test that dimensionality reduction with QUIC-SVD PCA works the same way
   * as the Exact-SVD PCA method.
   */
  BOOST_AUTO_TEST_CASE(QUICPCADimensionalityReductionTest)
  {
    arma::mat data, data1;
    data::Load("test_data_3_1000.csv", data);
    data1 = data;
  
    arma::mat backupData(data);
  
    // It isn't guaranteed that the QUIC-SVD will match with the exact SVD method,
    // starting with random samples. If this works 1 of 5 times, I'm fine with
    // that. All I want to know is that the QUIC-SVD method is able to solve the
    // task and is at least as good as the exact method (plus a little bit for
    // noise).
    size_t successes = 0;
    for (size_t trial = 0; trial < 5; ++trial)
    {
      if (trial > 0)
      {
        data = backupData;
        data1 = backupData;
      }
  
      PCA<ExactSVDPolicy> exactPCA;
      const double varRetainedExact = exactPCA.Apply(data, 1);
  
      PCA<QUICSVDPolicy> quicPCA;
      const double varRetainedQUIC = quicPCA.Apply(data1, 1);
  
      if (std::abs(varRetainedExact - varRetainedQUIC) < 0.2)
      {
        ++successes;
        break;
      }
    }
  
    BOOST_REQUIRE_GE(successes, 1);
    BOOST_REQUIRE_EQUAL(data.n_rows, data1.n_rows);
    BOOST_REQUIRE_EQUAL(data.n_cols, data1.n_cols);
  }
  
  /**
   * Test that setting the variance retained parameter to perform dimensionality
   * reduction works using the exact svd PCA method.
   */
  BOOST_AUTO_TEST_CASE(ExactPCAVarianceRetainedTest)
  {
    PCAVarianceRetained<ExactSVDPolicy>();
  }
  
  /**
   * Test that scaling PCA works.
   */
  BOOST_AUTO_TEST_CASE(PCAScalingTest)
  {
    // Generate an artificial dataset in 3 dimensions.
    arma::mat data(3, 5000);
  
    arma::vec mean("1.0 3.0 -12.0");
    arma::mat cov("1.0 0.9 0.0;"
                  "0.9 1.0 0.0;"
                  "0.0 0.0 12.0");
    GaussianDistribution g(mean, cov);
  
    for (size_t i = 0; i < 5000; ++i)
      data.col(i) = g.Random();
  
    // Now get the principal components when we are scaling.
    PCA<> p(true);
    arma::mat transData;
    arma::vec eigval;
    arma::mat eigvec;
  
    p.Apply(data, transData, eigval, eigvec);
  
    // The first two components of the eigenvector with largest eigenvalue should
    // be somewhere near sqrt(2) / 2.  The third component should be close to
    // zero.  There is noise, of course...
    BOOST_REQUIRE_CLOSE(std::abs(eigvec(0, 0)), sqrt(2) / 2, 0.35);
    BOOST_REQUIRE_CLOSE(std::abs(eigvec(1, 0)), sqrt(2) / 2, 0.35);
    BOOST_REQUIRE_SMALL(eigvec(2, 0), 0.1); // Large tolerance for noise.
  
    // The second component should be focused almost entirely in the third
    // dimension.
    BOOST_REQUIRE_SMALL(eigvec(0, 1), 0.1);
    BOOST_REQUIRE_SMALL(eigvec(1, 1), 0.1);
    BOOST_REQUIRE_CLOSE(std::abs(eigvec(2, 1)), 1.0, 0.35);
  
    // The third component should have the same absolute value characteristics as
    // the first (plus 20% tolerance).
    BOOST_REQUIRE_CLOSE(std::abs(eigvec(0, 0)), sqrt(2) / 2, 0.35);
    BOOST_REQUIRE_CLOSE(std::abs(eigvec(1, 0)), sqrt(2) / 2, 0.35);
    BOOST_REQUIRE_SMALL(eigvec(2, 0), 0.1); // Large tolerance for noise.
  
    // The eigenvalues should sum to three.
    BOOST_REQUIRE_CLOSE(accu(eigval), 3.0, 0.1); // 10% tolerance.
  }
  
  
  BOOST_AUTO_TEST_SUITE_END();
  