Natarajan_likelihood.hpp

Go to the documentation of this file.

 
#pragma once
 
#include "../utils/Likelihood.hpp"
#include <vector>

class Natarajan_likelihood : public Likelihood {
private:
  // Precomputed values for efficiency
  const double lgamma_delta1; 
  const double log_beta_alpha; 
  const double lgamma_delta2; 
  const double log_gamma_zeta; 
  std::vector<double> lgamma_zeta_mt_cache; 
  std::vector<double> lgamma_alpha_mh_cache; 
 
  std::vector<double> log_D_data;   
  const int D_cols; 

  double compute_cohesion(int point_index, int cluster_index,
                          const Eigen::Ref<const Eigen::VectorXi> &cls_ass_k,
                          int n_k) const;

  double compute_repulsion(int point_index, int cluster_index,
                           const Eigen::Ref<const Eigen::VectorXi> &cls_ass_k,
                           int n_k) const;
 
public:
  Natarajan_likelihood(const Data &data, const Params &param)
      : Likelihood(data, param), lgamma_delta1(lgamma(params.delta1)),
        log_beta_alpha(log(params.beta) * params.alpha - lgamma(params.alpha)),
        lgamma_delta2(lgamma(params.delta2)),
        log_gamma_zeta(log(params.gamma) * params.zeta - lgamma(params.zeta)),
        D_cols(params.D.cols()) {
    // Precompute log distances manually
    const int n = params.D.rows() * params.D.cols();
    log_D_data.resize(n);
 
    // Log distance precomputation
    const double *D_ptr = params.D.data();
    for (int i = 0; i < n; ++i) {
      log_D_data[i] = std::log(D_ptr[i]);
    }
 
    // Cache lgamma values for efficiency
    lgamma_alpha_mh_cache.resize(data.get_n() + 1, 0.0);
    lgamma_zeta_mt_cache.resize(data.get_n() + 1, 0.0);
 
    for (int val = 0; val <= data.get_n(); ++val) {
      lgamma_alpha_mh_cache[val] = lgamma(params.alpha + params.delta1 * val);
      lgamma_zeta_mt_cache[val] = lgamma(params.zeta + params.delta2 * val);
    }
 
  }

  double cluster_loglikelihood(int cluster_index) const override final;

  double cluster_loglikelihood(
      int cluster_index,
      const Eigen::Ref<const Eigen::VectorXi> &cls_ass_k) const override final __attribute__((hot));

  double point_loglikelihood_cond(int point_index, int cluster_index) const override final __attribute__((hot));
};