Natarajan_likelihood_summaryStats.hpp

Go to the documentation of this file.

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

class Natarajan_likelihood_summaryStats : public Likelihood {
private:
  const Natarajan_params &params; 
  const utils_params &utils;      
  // Precomputed values for efficiency
  const double lgamma_delta1; 
  const double
      log_beta_alpha;         
  const double lgamma_delta2; 
  const double log_gamma_zeta; 
  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_summaryStats(const Data &data,
                                    const Natarajan_params &params,
                                    const utils_params &utils)
      : Likelihood(data), params(params), utils(utils),
        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(utils.D.cols()) {
    // Precompute log distances manually
    const int n = utils.D.rows() * utils.D.cols();
    log_D_data.resize(n);
 
    // Log distance precomputation
    const double *D_ptr = utils.D.data();
    for (int i = 0; i < n; ++i) {
      log_D_data[i] = std::log(std::max(D_ptr[i], 1e-12));
    }
 
    // Cache lgamma values for efficiency
    lgamma_alpha_mh_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 = lgamma(params.zeta + params.delta2);
  }

  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 clusters_loglikelihood(int c1, int c2) const override final;

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

  double pairwise_repulsion(int cluster_1, int cluster_2) const;
};