<commit_before>#include "clustering.hpp"

#include "constant_series.hpp"
#include "monomials.hpp"
#include "nchoosek.hpp"
#include "p_max_total.hpp"


namespace fgt {


namespace {


double ddist(const arma::rowvec& x, const arma::rowvec& y) {
    return arma::accu(arma::pow(x - y, 2));
}
}


Clustering::Clustering(const arma::mat& source, arma::uword K, double bandwidth,
                       double epsilon)
    : m_source(source),
      m_K(K),
      m_indices(arma::zeros<arma::uvec>(source.n_rows)),
      m_centers(arma::zeros<arma::mat>(K, source.n_cols)),
      m_num_points(arma::zeros<arma::uvec>(K)),
      m_radii(K),
      m_rx(),
      m_bandwidth(bandwidth),
      m_epsilon(epsilon),
      m_p_max(0),
      m_constant_series() {}


Clustering::~Clustering() {}


arma::uword Clustering::choose_truncation_number(int dimensions,
                                                 double bandwidth,
                                                 double epsilon, double rx) {
    double r = std::min(std::sqrt(dimensions),
                        bandwidth * std::sqrt(std::log(1 / epsilon)));
    double rx2 = rx * rx;
    double h2 = bandwidth * bandwidth;
    double error = 1;
    double temp = 1;
    int p = 0;

    while ((error > epsilon) and (p <= TruncationNumberUpperLimit)) {
        ++p;
        double b = std::min((rx + std::sqrt(rx2 + 2 * p * h2)) / 2, rx + r);
        double c = rx - b;
        temp *= 2 * rx * b / h2 / p;
        error = temp * std::exp(-c * c / h2);
    }

    return p;
}


void Clustering::cluster() {
    cluster_impl();
    m_p_max =
        choose_truncation_number(m_source.n_cols, m_bandwidth, m_epsilon, m_rx);
    m_constant_series = compute_constant_series(m_source.n_cols, m_p_max);
}


arma::mat Clustering::compute_C(const arma::vec& q) const {
    arma::mat C = arma::zeros<arma::mat>(
        m_centers.n_rows, get_p_max_total(m_source.n_cols, m_p_max));
    std::vector<double> monomials(C.n_cols);
    double h2 = m_bandwidth * m_bandwidth;

    for (arma::uword i = 0; i < m_source.n_rows; ++i) {
        arma::uword k = m_indices(i);
        arma::rowvec dx = m_source.row(i) - m_centers.row(k);
        double distance2 = arma::accu(arma::pow(dx, 2));
        compute_monomials(dx / m_bandwidth, m_p_max, monomials);
        double f = q(i) * std::exp(-distance2 / h2);
        // TODO remove rowvec
        arma::rowvec center_monomials(monomials);
        C.row(k) += f * center_monomials;
    }

    arma::rowvec constant_series =
        compute_constant_series(m_source.n_cols, m_p_max);
    for (arma::uword i = 0; i < C.n_rows; ++i) {
        C.row(i) = C.row(i) % constant_series;
    }

    return C;
}


const optional_arma_uword_t GonzalezClustering::DefaultStartingIndex = {false,
                                                                        0};


GonzalezClustering::GonzalezClustering(const arma::mat& source, int K,
                                       double bandwidth, double epsilon,
                                       optional_arma_uword_t starting_index)
    : Clustering(source, K, bandwidth, epsilon),
      m_starting_index(starting_index) {}


void GonzalezClustering::cluster_impl() {
    arma::uword N = get_n_rows();
    arma::uword K = get_K();
    arma::uvec centers(K);
    arma::uvec cprev(N);
    arma::uvec cnext(N);<commit_msg>Don't allocate a new arma rowvec every monomial<commit_after>#include "clustering.hpp"

#include "constant_series.hpp"
#include "monomials.hpp"
#include "nchoosek.hpp"
#include "p_max_total.hpp"


namespace fgt {


namespace {


double ddist(const arma::rowvec& x, const arma::rowvec& y) {
    return arma::accu(arma::pow(x - y, 2));
}
}


Clustering::Clustering(const arma::mat& source, arma::uword K, double bandwidth,
                       double epsilon)
    : m_source(source),
      m_K(K),
      m_indices(arma::zeros<arma::uvec>(source.n_rows)),
      m_centers(arma::zeros<arma::mat>(K, source.n_cols)),
      m_num_points(arma::zeros<arma::uvec>(K)),
      m_radii(K),
      m_rx(),
      m_bandwidth(bandwidth),
      m_epsilon(epsilon),
      m_p_max(0),
      m_constant_series() {}


Clustering::~Clustering() {}


arma::uword Clustering::choose_truncation_number(int dimensions,
                                                 double bandwidth,
                                                 double epsilon, double rx) {
    double r = std::min(std::sqrt(dimensions),
                        bandwidth * std::sqrt(std::log(1 / epsilon)));
    double rx2 = rx * rx;
    double h2 = bandwidth * bandwidth;
    double error = 1;
    double temp = 1;
    int p = 0;

    while ((error > epsilon) and (p <= TruncationNumberUpperLimit)) {
        ++p;
        double b = std::min((rx + std::sqrt(rx2 + 2 * p * h2)) / 2, rx + r);
        double c = rx - b;
        temp *= 2 * rx * b / h2 / p;
        error = temp * std::exp(-c * c / h2);
    }

    return p;
}


void Clustering::cluster() {
    cluster_impl();
    m_p_max =
        choose_truncation_number(m_source.n_cols, m_bandwidth, m_epsilon, m_rx);
    m_constant_series = compute_constant_series(m_source.n_cols, m_p_max);
}


arma::mat Clustering::compute_C(const arma::vec& q) const {
    arma::mat C = arma::zeros<arma::mat>(
        m_centers.n_rows, get_p_max_total(m_source.n_cols, m_p_max));
    std::vector<double> monomials(C.n_cols);
    double h2 = m_bandwidth * m_bandwidth;

    for (arma::uword i = 0; i < m_source.n_rows; ++i) {
        arma::uword k = m_indices(i);
        arma::rowvec dx = m_source.row(i) - m_centers.row(k);
        double distance2 = arma::accu(arma::pow(dx, 2));
        compute_monomials(dx / m_bandwidth, m_p_max, monomials);
        double f = q(i) * std::exp(-distance2 / h2);
        for (arma::uword j = 0; j < monomials.size(); ++j) {
            C(k, j) += f * monomials[j];
        }
    }

    arma::rowvec constant_series =
        compute_constant_series(m_source.n_cols, m_p_max);
    for (arma::uword i = 0; i < C.n_rows; ++i) {
        C.row(i) = C.row(i) % constant_series;
    }

    return C;
}


const optional_arma_uword_t GonzalezClustering::DefaultStartingIndex = {false,
                                                                        0};


GonzalezClustering::GonzalezClustering(const arma::mat& source, int K,
                                       double bandwidth, double epsilon,
                                       optional_arma_uword_t starting_index)
    : Clustering(source, K, bandwidth, epsilon),
      m_starting_index(starting_index) {}


void GonzalezClustering::cluster_impl() {
    arma::uword N = get_n_rows();
    arma::uword K = get_K();
    arma::uvec centers(K);
    arma::uvec cprev(N);
    arma::uvec cnext(N);
{'fraction_non_alphanumeric': 0.1365121806298277, 'fraction_numerical': 0.006981580510992276, 'mean_word_length': 2.407388663967611, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1025062', 'n_tokens_mistral': 2465, 'n_tokens_neox': 2279, 'n_words': 461}