/**
 * @file exact_svd_method.hpp
 * @author Ajinkya Kale
 * @author Ryan Curtin
 * @author Marcus Edel
 *
 * Implementation of the exact svd method for use in the Principal Components
 * Analysis method.
 */

#ifndef MLPACK_METHODS_PCA_DECOMPOSITION_POLICIES_ITERATIVE_PCA_METHOD_HPP
#define MLPACK_METHODS_PCA_DECOMPOSITION_POLICIES_ITERATIVE_PCA_METHOD_HPP

#include <mlpack/core.hpp>

#include <chrono>

namespace mlpack {
namespace pca {

/**
 * Implementation of the iterative PCA policy.
 */
class IterativePCAPolicy
{
public:
  /**
   * @param dimension desired dimension of the transformed data
   * @param iterateTime upper bound iteration times of each eigen vector
   */

  explicit IterativePCAPolicy(size_t dimension, size_t iterateTime = 50) :
      dimension(dimension),
      iterateTime(iterateTime)
  {}

  /**
   * Apply Iterative PCA to the provided data set using the
   * exact SVD method.
   *
   * @param data Data matrix.
   * @param centeredData Centered data matrix.
   * @param transformedData Matrix to put results of PCA into.
   * @param eigVal Vector to put eigenvalues into.
   * @param eigvec Matrix to put eigenvectors (loadings) into.
   * @param rank Rank of the decomposition.
   */
  void Apply(const arma::mat& /* data */,
             const arma::mat& centeredData,
             arma::mat& transformedData,
             arma::vec& /* eigval */,
             arma::mat& eigvec,
             const size_t /* rank */)
  {
    std::chrono::time_point<std::chrono::system_clock> start, end;
    start = std::chrono::system_clock::now();
    const arma::mat covMat = arma::cov(arma::trans(centeredData));
    end = std::chrono::system_clock::now();
    std::cout<<"covMat elapsed : "<<
               std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()<<std::endl;

    start = std::chrono::system_clock::now();
    eigvec = arma::randu<arma::mat>(covMat.n_rows, dimension);
    end = std::chrono::system_clock::now();
    std::cout<<"randu elapsed : "<<
               std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()<<std::endl;

    start = std::chrono::system_clock::now();
    GramSchmidtIter(covMat, eigvec);
    end = std::chrono::system_clock::now();
    std::cout<<"GramSchmidtIter elapsed : "<<
               std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()<<std::endl;

    start = std::chrono::system_clock::now();
    transformedData = arma::trans(eigvec) * centeredData;
    end = std::chrono::system_clock::now();
    std::cout<<"transformedData elapsed : "<<
               std::chrono::duration_cast<std::chrono::microseconds>(end - start).count()<<std::endl;
  }

  void Dimension(size_t dimension)
  {
    dimension = dimension;
  }

  size_t Dimension()
  {
    return dimension;
  }

  size_t Dimension() const
  {
    return dimension;
  }

  void IterateTime(size_t iterateTime)
  {
    iterateTime = iterateTime;
  }

  size_t IterateTime()
  {
    return iterateTime;
  }

  size_t IterateTime() const
  {
    return iterateTime;
  }

private:
  void GramSchmidtIter(const arma::mat &covMat, arma::mat &eigvec) const
  {
    arma::mat temp(eigvec.n_rows, 1);
    arma::mat innerProd;
    for(size_t i = 0; i != dimension; ++i)
    {
      size_t iterate = 0;
      do{
        temp.zeros();
        eigvec.col(i) = covMat * eigvec.col(i);
        const arma::mat eivecTrans = arma::trans(eigvec.col(i));
        for(size_t j = 0; j <= i; ++j){
          innerProd = eivecTrans * eigvec.col(j);
          temp += innerProd(0) * eigvec.col(j);
        }
        //const arma::mat temp = arma::repmat(arma::trans(eigvec.col(i)), i + 1, 1) *
        //    eigvec.submat(0,0,eigvec.n_rows-1, i);
        eigvec.col(i) -= temp;
        eigvec.col(i) /= std::sqrt(arma::sum(eigvec.col(i) % eigvec.col(i)));
        ++iterate;
      }while(iterate < iterateTime);
    }
  }

  size_t dimension;
  size_t iterateTime;
};

} // namespace pca
} // namespace mlpack

#endif