cluster.h

Go to the documentation of this file.

#ifndef CLUSTER_H
#define CLUSTER_H
 
#include <cmp_defines.h>
#include <distribution.h>
#include <grid.h>
 
namespace cmp::cluster {
 
class GeometricCluster {
  private:
    Eigen::VectorXs labels_;    
    Eigen::MatrixXd centroids_; 
 
    size_t nClusters_;          
    size_t nPoints_;            
    size_t dim_;                
 
  public:
    GeometricCluster() = default;
 
    bool fit(const Eigen::Ref<const Eigen::MatrixXd> &points, size_t nClusters, std::default_random_engine &rng, size_t max_iter = 1000);
 
    const size_t &operator[](size_t i) const;
 
    Eigen::VectorXd centroid(size_t i) const;
 
    size_t nClusters() const;
 
    size_t nPoints() const;
 
    size_t dim() const;
 
    const Eigen::VectorXs &getLabels() const {
        return labels_;
    }
};
 
class DirichletProcessMixtureModel {
 
  private:
    struct Cluster {
        size_t id_;                 
        size_t nPoints_;            
 
        // Sufficient statistics
        Eigen::VectorXd sumOfX_;       
        Eigen::MatrixXd sumOfXXT_;     
 
        // Default constructor
        Cluster() : id_(0), nPoints_(0) {}
 
        // Constructor with id and dimension
        Cluster(size_t id, size_t dim) : id_(id), nPoints_(0),
            sumOfX_(Eigen::VectorXd::Zero(dim)),
            sumOfXXT_(Eigen::MatrixXd::Zero(dim, dim)) {}
    };
 
  private:
 
    // Workspaces for the Gibbs sampler to avoid heap allocations
    std::vector<double> logp_workspace_;        
    std::vector<double> probs_workspace_;       
    std::vector<size_t> clusterIDs_workspace_;  
 
    Eigen::MatrixXd xObs_;                      
    Eigen::VectorXs labels_;                    
    size_t nPoints_{0};                         
    int dim_{0};                                
    std::vector<size_t>  clusterIDs_;           
 
    // This prior is conjugate to the Gaussian likelihood
    cmp::distribution::NormalInverseWishartDistribution hyper_; 
    cmp::distribution::GammaDistribution alphaPrior_;           
    double alpha_{1.0};                                         
 
    std::map<size_t, Cluster> clusters_;        
    size_t nextClusterId_;                      
 
    // Random number generator for sampling
    std::default_random_engine rng_;            
    std::uniform_real_distribution<double> distU_{0.0, 1.0}; 
 
  public:
 
// public API
    DirichletProcessMixtureModel(
        double alpha,
        const cmp::distribution::GammaDistribution &alphaPrior,
        const cmp::distribution::NormalInverseWishartDistribution &hyper,
        unsigned int seed = 12345
    );
 
    void condition(const Eigen::Ref<const Eigen::MatrixXd> &data, const Eigen::VectorXs& init_labels);
 
    void step();
 
    void remapLabels();
 
    Eigen::VectorXs getLabels() const {
        return labels_;
    }
 
    size_t nClusters() const {
        return clusters_.size();
    }
 
    double getAlpha() const {
        return alpha_;
    }
 
  private:
 
    void init(const Eigen::VectorXs& init_labels);
 
    void removePointFromCluster(const size_t &pointIndex, const size_t &clusterID);
 
    void addPointToCluster(const size_t &pointIndex, const size_t &clusterID);
 
    double logNIWPosterior(const Eigen::VectorXd& x, const Cluster& c) const;
 
    void updateAlpha() {
        double K = static_cast<double>(clusters_.size());
        double N = static_cast<double>(nPoints_);
 
        // Extract prior parameters from your distribution object
        double a = alphaPrior_.getAlpha();
        double b = alphaPrior_.getBeta();
 
        // 1. Instantiate the Beta distribution and sample eta
        cmp::distribution::BetaDistribution beta_dist(alpha_ + 1.0, N);
        double eta = beta_dist.sample(rng_);
 
        // 2. Calculate the weights for the Gamma mixture
        double weight1 = (a + K - 1.0) / (N * (b - std::log(eta)));
        double pi_eta = weight1 / (weight1 + 1.0);
 
        // 3. Sample the new alpha using your custom cmp::distribution classes
        double u = distU_(rng_);
        double updated_beta = b - std::log(eta); // The rate parameter
 
        if(u < pi_eta) {
            // Gamma(a + K, b - ln(eta))
            cmp::distribution::GammaDistribution new_gamma(a + K, updated_beta);
            alpha_ = new_gamma.sample(rng_);
        } else {
            // Gamma(a + K - 1, b - ln(eta))
            cmp::distribution::GammaDistribution new_gamma(a + K - 1.0, updated_beta);
            alpha_ = new_gamma.sample(rng_);
        }
 
        // Safety bounds just in case of severe numerical underflow
        if(alpha_ < 1e-5) alpha_ = 1e-5;
    }
 
};
 
 
class DummyCluster {
  private:
    size_t nClusters_;               
    std::default_random_engine rng_;  
    Eigen::VectorXs labels_;         
  public:
    DummyCluster(size_t nClusters, unsigned int seed = 42)
        : nClusters_(nClusters), rng_(seed) {}
 
    void fit(const Eigen::Ref<const Eigen::MatrixXd> &points) {
        size_t nPoints = points.rows();
        std::uniform_int_distribution<size_t> dist(0, nClusters_ - 1);
        labels_ = Eigen::VectorXs::Zero(nPoints);
        for(size_t i = 0; i < nPoints; i++) {
            labels_(i) = dist(rng_);
        }
    }
 
    const Eigen::VectorXs &getLabels() const {
        return labels_;
    }
 
 
};
 
} // namespace cmp::cluster
 
 
#endif