splitmerge_LSS_SDDS.hpp

Go to the documentation of this file.

 
#pragma once
 
#include "../params/utils-params.hpp"
#include "../utils/Sampler.hpp"

class SplitMerge_LSS_SDDS : public Sampler {
private:
  const utils_params &params;
 
  std::vector<double> row_max_distances;
  double temperature;
 
  // ========== Random Number Generation ==========

  mutable std::mt19937 gen;
 
  // ========== Move Selection Variables ==========

  int idx_i;

  int idx_j;

  int ci;

  int cj;
 
  // ========== Algorithm Configuration ==========

  bool shuffle_bool = false;
 
  // ========== State Management ==========

  Eigen::VectorXi launch_state;

  Eigen::VectorXi S;

  int launch_state_size;

  const Eigen::VectorXi &original_allocations;
 
  // ========== Proposal Probabilities ==========

  std::uniform_real_distribution<> dis_real;

  std::uniform_int_distribution<> dis_int;

  Eigen::Vector2d log_probs;

  Eigen::Vector2d probs;

  double log_split_gibbs_prob = 0;

  double log_merge_gibbs_prob = 0;

  const double rand_split_prob = log(0.5);
 
  // ========== Debug variables ==========
  int accepted_split = 0;
  int accepted_merge = 0;
  int accepted_shuffle = 0;
  int split_moves = 0;
  int merge_moves = 0;
  int shuffle_moves = 0;
 
  // ========== Move Selection Methods ==========

  void choose_indeces(bool similarity = false);

  void choose_clusters_shuffle();
 
  // ========== SAMS-Specific Proposal Generation ==========

  void sequential_allocation(int iterations, bool only_probabilities = false,
                             bool sequential = true);
 
  // ========== Split Move Implementation ==========

  void smart_split_move();

  void dumb_split_move();

  double compute_acceptance_ratio_split(double likelihood_old_cluster);
 
  // ========== Merge Move Implementation ==========

  void smart_merge_move();

  void dumb_merge_move();

  double compute_acceptance_ratio_merge(double likelihood_old_clusters);
 
  // ========== Shuffle Move Implementation ==========

  void shuffle();

  double compute_acceptance_ratio_shuffle(double likelihood_old_clusters,
                                          int old_ci_size, int old_cj_size);
 
public:
  // ========== Constructor ==========

  SplitMerge_LSS_SDDS(Data &d, utils_params &p, Likelihood &l, Process &pr,
                      bool shuffle)
      : Sampler(d, l, pr), params(p), shuffle_bool(shuffle),
        original_allocations(pr.old_allocations_view()), gen(rd()),
        dis_real(0.0, 1.0), dis_int(0, 1) {
 
    // Precompute row max distances for efficiency
    row_max_distances.resize(data.get_n());
    for (int i = 0; i < data.get_n(); ++i) {
      row_max_distances[i] = params.D.row(i).maxCoeff();
    }
    const int n = data.get_n();
 
    // Take all pairwise distances (i < j)
    std::vector<double> pairwise_dists;
    pairwise_dists.reserve(n * (n - 1) / 2);
    for (int i = 0; i < n; ++i) {
      for (int j = i + 1; j < n; ++j) {
        pairwise_dists.push_back(params.D(i, j));
      }
    }
 
    // Compute the 50th percentile of pairwise distances
    if (!pairwise_dists.empty()) {
      auto q_it = pairwise_dists.begin() +
                  static_cast<size_t>(pairwise_dists.size() * 0.50);
      std::nth_element(pairwise_dists.begin(), q_it, pairwise_dists.end());
      temperature = *q_it;
      if (temperature < 1e-12)
        temperature = 1.0;
    } else {
      temperature = 1.0;
    }
  };
 
  // ========== MCMC Interface ==========

  void step() override final;
 
  // ========== Accessor Methods ==========
  double get_accepted_split() const {
    return split_moves > 0 ? static_cast<double>(accepted_split) / split_moves : 0.0;
  };

  double get_accepted_merge() const {
    return merge_moves > 0 ? static_cast<double>(accepted_merge) / merge_moves : 0.0;
  };

  double get_accepted_shuffle() const {
    return shuffle_moves > 0 ? static_cast<double>(accepted_shuffle) / shuffle_moves : 0.0;
  };
};