<commit_before>/** \file CompressedMatrix.h
 *
 *  This file defines the templated class CompressedMatrix.
 *
 *  \copyright Copyright (c) 2014,
 *             Universita della Svizzera italiana (USI) &
 *             Centro Svizzero di Calcolo Scientifico (CSCS).
 *             All rights reserved.
 *             This software may be modified and distributed under the terms
 *             of the BSD license. See the LICENSE file for details.
 *
 *  \author <NAME> (CSCS)
 *  \author <NAME> (CSCS)
 *  \author <NAME> (CSCS)
 */

#pragma once
#include <cstdlib>  // std::srand (Eigen only)
#include <iostream> // std::cout, std::endl
#include <vector>   // std::vector
#include <limits>   // std::numeric_limits
#include <algorithm>// std::max
#include <mpi.h>    // MPI_Comm_size, MPI_Comm_rank, MPI_Allreduce, MPI_Allgather
#include "mpi_type_helper.h"
#include "matrices.h"

#if defined(USE_EIGEN)
#include "lanczos_correlation_eigen.h"
#elif defined(USE_MINLIN)
#include "lanczos_correlation_minlin.h"
#endif


/**
 * The class CompressedMatrix stores all the information needed to represent
 * the data in compressed form. The compression is done in the constructor
 * and there is a member function for reconstructing the full matrix. In
 * addition, there are two public variables 'original_size' and
 * 'compressed_size' with the number of elements stored for the full or
 * compressed matrix respectively.
 */
template<class Scalar>
class CompressedMatrix
{

public:

  int original_size;    ///< \brief The number of elements (scalar numbers) in
                        ///< the original (uncompressed) matrix.
  int compressed_size;  ///< \brief The number of elements (scalar numbers)
                        ///< used for the compressed representation of the
                        ///< matrix.

  /**
   * \brief Construct a compressed representation of a matrix.
   *
   * The constructor of CompressedMatrix does all the compression work. We
   * set up the data structures and call member functions to find a good
   * clustering and compress the data.
   *
   * \see initialize_data(), do_iterative_clustering(), do_final_pca(),
   *      calculate_reduced_form()
   *
   * \param[in] X           The matrix that is to be compressed.
   * \param[in] K           The number of clusters used for compression.
   * \param[in] M           The number of eigenvectors used per cluster.
   * \param[in] column_ids  A vector of unique IDs for the column of the
   *                        matrix, independent of the number of processes,
   *                        used for assigning the initial clusters.
   */
  CompressedMatrix(const DeviceMatrix<Scalar> &X, const int K, const int M, std::vector<int> column_ids) {

    Nc_ = X.rows(); // number of entries along compressed direction
    Nd_ = X.cols(); // number of entries along distributed direction
    K_  = K;        // number of clusters
    M_  = M;        // number of eigenvectors per cluster

    initialize_data(column_ids);
    do_iterative_clustering(X);
    do_final_pca(X);
    calculate_reduced_form(X);

  }

  /**
   * \brief Reconstruct the full matrix.
   *
   * This function reconstructs the full matrix using the compressed data.
   * This matrix should be similar to the original matrix, except for the
   * compression losses.
   *
   * \return      The reconstructed matrix.
   */
  DeviceMatrix<Scalar> reconstruct() {

    DeviceMatrix<Scalar> X_reconstructed(Nc_, Nd_);

    // This loop can be multithreaded, if all threads have a separate
    // copy of X_translated
    for(int k = 0; k < K_; k++) {
      std::vector<int> current_cluster_indices = indices_for_cluster(k);
      for (int m = 0; m < current_cluster_indices.size() ; m++ ) {
        GET_COLUMN(X_reconstructed, current_cluster_indices[m]) = eigenvectors_[k]
            * GET_COLUMN(X_reduced_[k], current_cluster_indices[m]);
        GET_COLUMN(X_reconstructed, current_cluster_indices[m]) += GET_COLUMN(cluster_means_, k);
      }
    }

    return X_reconstructed;
  }


private:

  /// \brief The reduced representation of the data. For each column of the
  /// original matrix, this contains the coefficients for the eigenvectors of
  /// its cluster. This corresponds to the projection of the original column
  /// vector into the subspace spanned by the eigenvectors.
  std::vector< DeviceMatrix<Scalar> > X_reduced_;
  /// The mean vector for each cluster.
  DeviceMatrix<Scalar> cluster_means_;
  /// The M largest eigenvectors for each cluster.
  std::vector< DeviceMatrix<Scalar> > eigenvectors_;
  /// The index of the cluster each column of the original matrix belongs to.
  std::vector<int> cluster_indices_;

  // data sizes
  int K_;         ///< The number of clusters.
  int M_;         ///< The number of eigenvectors calculated fo each cluster.
  int Nc_;        ///< The number of rows of the original matrix.
  int Nd_;        ///< \brief The number of columns of the original matrix
                  ///< that are assigned to the current process.
  int Nd_total_;  ///< The total number of columns of the original matrix.

  // information about MPI
  int mpi_processes_;
  int my_rank_;   ///< \brief The rank (ID) of the current MPI process. This
                  ///< is mainly used to limit the console output to one
                  ///< process.

  /**
   * This function is called by the constructor and sets up the member
   * variables with the correct sizes. It initializes the cluster indicies
   * based on the column ids that are passed to it and it calculates the
   * original and compressed matrix sizes based on the matrix dimensions
   * and parameters.
   *
   * \param[in] column_ids  A vector with an ID for each column of the data
   *                        matrix that is passed to the current process.
   */
  void initialize_data(std::vector<int> column_ids) {
    
    // collect number of columns each process has
    MPI_Comm_size(MPI_COMM_WORLD, &mpi_processes_);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank_);
    int *Nd_global = new int[mpi_processes_];
    MPI_Allgather( &Nd_, 1, MPI_INT, Nd_global, 1, MPI_INT,
        MPI_COMM_WORLD);
    Nd_total_ = 0;
    for (int i=0; i<mpi_processes_; i++) {Nd_total_ += Nd_global[i];}

    // set up data
    cluster_indices_ = std::vector<int>(Nd_, 0);
    cluster_means_ = DeviceMatrix<Scalar>(Nc_, K_);
    eigenvectors_ = std::vector< DeviceMatrix<Scalar> >(K_, DeviceMatrix<Scalar>(Nc_, M_));
    X_reduced_ = std::vector< DeviceMatrix<Scalar> >(K_, DeviceMatrix<Scalar>(M_, Nd_));
    
    // initialize clusters independent from number of processes
    for (int i=0; i<Nd_; i++) cluster_indices_[i] = column_ids[i]%K_;

    // set up sizes
    original_size = Nc_ * Nd_total_;
    compressed_size = Nd_total_ + Nc_ * K_ + K_ * Nc_ * M_ + K_ * M_ * Nd_total_;

    delete[] Nd_global;

  }

  /**
   * This function iterates until a near optimal clustering is found. For each
   * iteration, it does a PCA for each cluster finding the largest eigenvector
   * only. It then reassigns all columns to the cluster where they are best
   * represented by the eigenvector. This is repeated until the norm measuring
   * the difference between the columns and the eigenvectors doesn't change<commit_msg>don’t save column numbers for all processes in CompressedMatrix initialization, as it is no longer needed for initial clustering<commit_after>/** \file CompressedMatrix.h
 *
 *  This file defines the templated class CompressedMatrix.
 *
 *  \copyright Copyright (c) 2014,
 *             Universita della Svizzera italiana (USI) &
 *             Centro Svizzero di Calcolo Scientifico (CSCS).
 *             All rights reserved.
 *             This software may be modified and distributed under the terms
 *             of the BSD license. See the LICENSE file for details.
 *
 *  \author <NAME> (CSCS)
 *  \author <NAME> (CSCS)
 *  \author <NAME> (CSCS)
 */

#pragma once
#include <cstdlib>  // std::srand (Eigen only)
#include <iostream> // std::cout, std::endl
#include <vector>   // std::vector
#include <limits>   // std::numeric_limits
#include <algorithm>// std::max
#include <mpi.h>    // MPI_Comm_size, MPI_Comm_rank, MPI_Allreduce, MPI_Allgather
#include "mpi_type_helper.h"
#include "matrices.h"

#if defined(USE_EIGEN)
#include "lanczos_correlation_eigen.h"
#elif defined(USE_MINLIN)
#include "lanczos_correlation_minlin.h"
#endif


/**
 * The class CompressedMatrix stores all the information needed to represent
 * the data in compressed form. The compression is done in the constructor
 * and there is a member function for reconstructing the full matrix. In
 * addition, there are two public variables 'original_size' and
 * 'compressed_size' with the number of elements stored for the full or
 * compressed matrix respectively.
 */
template<class Scalar>
class CompressedMatrix
{

public:

  int original_size;    ///< \brief The number of elements (scalar numbers) in
                        ///< the original (uncompressed) matrix.
  int compressed_size;  ///< \brief The number of elements (scalar numbers)
                        ///< used for the compressed representation of the
                        ///< matrix.

  /**
   * \brief Construct a compressed representation of a matrix.
   *
   * The constructor of CompressedMatrix does all the compression work. We
   * set up the data structures and call member functions to find a good
   * clustering and compress the data.
   *
   * \see initialize_data(), do_iterative_clustering(), do_final_pca(),
   *      calculate_reduced_form()
   *
   * \param[in] X           The matrix that is to be compressed.
   * \param[in] K           The number of clusters used for compression.
   * \param[in] M           The number of eigenvectors used per cluster.
   * \param[in] column_ids  A vector of unique IDs for the column of the
   *                        matrix, independent of the number of processes,
   *                        used for assigning the initial clusters.
   */
  CompressedMatrix(const DeviceMatrix<Scalar> &X, const int K, const int M, std::vector<int> column_ids) {

    Nc_ = X.rows(); // number of entries along compressed direction
    Nd_ = X.cols(); // number of entries along distributed direction
    K_  = K;        // number of clusters
    M_  = M;        // number of eigenvectors per cluster

    initialize_data(column_ids);
    do_iterative_clustering(X);
    do_final_pca(X);
    calculate_reduced_form(X);

  }

  /**
   * \brief Reconstruct the full matrix.
   *
   * This function reconstructs the full matrix using the compressed data.
   * This matrix should be similar to the original matrix, except for the
   * compression losses.
   *
   * \return      The reconstructed matrix.
   */
  DeviceMatrix<Scalar> reconstruct() {

    DeviceMatrix<Scalar> X_reconstructed(Nc_, Nd_);

    // This loop can be multithreaded, if all threads have a separate
    // copy of X_translated
    for(int k = 0; k < K_; k++) {
      std::vector<int> current_cluster_indices = indices_for_cluster(k);
      for (int m = 0; m < current_cluster_indices.size() ; m++ ) {
        GET_COLUMN(X_reconstructed, current_cluster_indices[m]) = eigenvectors_[k]
            * GET_COLUMN(X_reduced_[k], current_cluster_indices[m]);
        GET_COLUMN(X_reconstructed, current_cluster_indices[m]) += GET_COLUMN(cluster_means_, k);
      }
    }

    return X_reconstructed;
  }


private:

  /// \brief The reduced representation of the data. For each column of the
  /// original matrix, this contains the coefficients for the eigenvectors of
  /// its cluster. This corresponds to the projection of the original column
  /// vector into the subspace spanned by the eigenvectors.
  std::vector< DeviceMatrix<Scalar> > X_reduced_;
  /// The mean vector for each cluster.
  DeviceMatrix<Scalar> cluster_means_;
  /// The M largest eigenvectors for each cluster.
  std::vector< DeviceMatrix<Scalar> > eigenvectors_;
  /// The index of the cluster each column of the original matrix belongs to.
  std::vector<int> cluster_indices_;

  // data sizes
  int K_;         ///< The number of clusters.
  int M_;         ///< The number of eigenvectors calculated fo each cluster.
  int Nc_;        ///< The number of rows of the original matrix.
  int Nd_;        ///< \brief The number of columns of the original matrix
                  ///< that are assigned to the current process.
  int Nd_total_;  ///< The total number of columns of the original matrix.

  // information about MPI
  int mpi_processes_;
  int my_rank_;   ///< \brief The rank (ID) of the current MPI process. This
                  ///< is mainly used to limit the console output to one
                  ///< process.

  /**
   * This function is called by the constructor and sets up the member
   * variables with the correct sizes. It initializes the cluster indicies
   * based on the column ids that are passed to it and it calculates the
   * original and compressed matrix sizes based on the matrix dimensions
   * and parameters.
   *
   * \param[in] column_ids  A vector with an ID for each column of the data
   *                        matrix that is passed to the current process.
   */
  void initialize_data(std::vector<int> column_ids) {
    
    // collect number of columns each process has
    MPI_Comm_size(MPI_COMM_WORLD, &mpi_processes_);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank_);
    MPI_Allreduce(&Nd_, &Nd_total_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

    // set up data
    cluster_indices_ = std::vector<int>(Nd_, 0);
    cluster_means_ = DeviceMatrix<Scalar>(Nc_, K_);
    eigenvectors_ = std::vector< DeviceMatrix<Scalar> >(K_, DeviceMatrix<Scalar>(Nc_, M_));
    X_reduced_ = std::vector< DeviceMatrix<Scalar> >(K_, DeviceMatrix<Scalar>(M_, Nd_));
    
    // initialize clusters independent from number of processes
    for (int i=0; i<Nd_; i++) cluster_indices_[i] = column_ids[i]%K_;

    // set up sizes
    original_size = Nc_ * Nd_total_;
    compressed_size = Nd_total_ + Nc_ * K_ + K_ * Nc_ * M_ + K_ * M_ * Nd_total_;
  }

  /**
   * This function iterates until a near optimal clustering is found. For each
   * iteration, it does a PCA for each cluster finding the largest eigenvector
   * only. It then reassigns all columns to the cluster where they are best
   * represented by the eigenvector. This is repeated until the norm measuring
   * the difference between the columns and the eigenvectors doesn't change
{'fraction_non_alphanumeric': 0.08913899054063384, 'fraction_numerical': 0.0014499758337361043, 'mean_word_length': 2.8675567423230977, 'pattern_counts': {'":': 0, '<': 98, '<?xml version=': 0, '>': 63, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7211782', 'n_tokens_mistral': 4280, 'n_tokens_neox': 4069, 'n_words': 1562}