Merge pull request #66 from FMKerckhof/master

Potential fixes for #48 (assisted by copilot)

Author: mlampros

.Rbuildignore

@@ -5,3 +5,5 @@
 ^_pkgdown\.yml$
 ^docs$
 ^pkgdown$
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -1,3 +1,4 @@
 docs/
 .Rhistory
 docs
+.Rproj.user

DESCRIPTION

@@ -3,7 +3,7 @@ Type: Package
 Title: Gaussian Mixture Models, K-Means, Mini-Batch-Kmeans, K-Medoids and Affinity Propagation Clustering
 Version: 1.3.4
 Date: 2025-09-14
-Authors@R: c( person(given = "Lampros", family = "Mouselimis", email = "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-8024-1546")), person(given = "Conrad", family = "Sanderson", role = "cph", comment = "Author of the C++ Armadillo library"), person(given = "Ryan", family = "Curtin", role = "cph", comment = "Author of the C++ Armadillo library"), person(given = "Siddharth", family = "Agrawal", role = "cph", comment = "Author of the C code of the Mini-Batch-Kmeans algorithm (https://github.com/siddharth-agrawal/Mini-Batch-K-Means)"), person(given = "Brendan", family = "Frey", email = "[email protected]", role = "cph", comment = "Author of the matlab code of the Affinity propagation algorithm (for commercial use please contact the author of the matlab code)"), person(given = "Delbert", family = "Dueck", role = "cph", comment = "Author of the matlab code of the Affinity propagation algorithm"), person(given = "Vitalie", family = "Spinu", email = "[email protected]", role = "ctb", comment = "Github Contributor") )  
+Authors@R: c( person(given = "Lampros", family = "Mouselimis", email = "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-8024-1546")), person(given = "Conrad", family = "Sanderson", role = "cph", comment = "Author of the C++ Armadillo library"), person(given = "Ryan", family = "Curtin", role = "cph", comment = "Author of the C++ Armadillo library"), person(given = "Siddharth", family = "Agrawal", role = "cph", comment = "Author of the C code of the Mini-Batch-Kmeans algorithm (https://github.com/siddharth-agrawal/Mini-Batch-K-Means)"), person(given = "Brendan", family = "Frey", email = "[email protected]", role = "cph", comment = "Author of the matlab code of the Affinity propagation algorithm (for commercial use please contact the author of the matlab code)"), person(given = "Delbert", family = "Dueck", role = "cph", comment = "Author of the matlab code of the Affinity propagation algorithm"), person(given = "Vitalie", family = "Spinu", email = "[email protected]", role = "ctb", comment = "Github Contributor"),person(given = "Frederiek - Maarten", family = "Kerckhof", email = "[email protected]", role = "ctb", comment = "Github Contributor") )  
 BugReports: https://github.com/mlampros/ClusterR/issues
 URL: https://github.com/mlampros/ClusterR, https://mlampros.github.io/ClusterR/
 Description: Gaussian mixture models, k-means, mini-batch-kmeans, k-medoids and affinity propagation clustering with the option to plot, validate, predict (new data) and estimate the optimal number of clusters. The package takes advantage of 'RcppArmadillo' to speed up the computationally intensive parts of the functions. For more information, see (i) "Clustering in an Object-Oriented Environment" by Anja Struyf, Mia Hubert, Peter Rousseeuw (1997), Journal of Statistical Software, <doi:10.18637/jss.v001.i04>; (ii) "Web-scale k-means clustering" by D. Sculley (2010), ACM Digital Library, <doi:10.1145/1772690.1772862>; (iii) "Armadillo: a template-based C++ library for linear algebra" by Sanderson et al (2016), The Journal of Open Source Software, <doi:10.21105/joss.00026>; (iv) "Clustering by Passing Messages Between Data Points" by Brendan J. Frey and Delbert Dueck, Science 16 Feb 2007: Vol. 315, Issue 5814, pp. 972-976, <doi:10.1126/science.1136800>.
@@ -30,4 +30,4 @@ Suggests:
     knitr,
     rmarkdown
 VignetteBuilder: knitr
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3

R/RcppExports.R

@@ -49,6 +49,10 @@ predict_MGausDPDF <- function(data, CENTROIDS, COVARIANCE, WEIGHTS, eps = 1.0e-8
     .Call(`_ClusterR_predict_MGausDPDF`, data, CENTROIDS, COVARIANCE, WEIGHTS, eps)
 }
 
+predict_MGausDPDF_full <- function(data, CENTROIDS, COVARIANCE, WEIGHTS, eps = 1.0e-8) {
+    .Call(`_ClusterR_predict_MGausDPDF_full`, data, CENTROIDS, COVARIANCE, WEIGHTS, eps)
+}
+
 GMM_arma_AIC_BIC <- function(data, max_clusters, dist_mode, seed_mode, km_iter, em_iter, verbose, var_floor = 1e-10, criterion = "AIC", seed = 1L) {
     .Call(`_ClusterR_GMM_arma_AIC_BIC`, data, max_clusters, dist_mode, seed_mode, km_iter, em_iter, verbose, var_floor, criterion, seed)
 }

R/clustering_functions.R

@@ -131,12 +131,13 @@ GMM = function(data,
 #'
 #' @param data matrix or data frame
 #' @param CENTROIDS matrix or data frame containing the centroids (means), stored as row vectors
-#' @param COVARIANCE matrix or data frame containing the diagonal covariance matrices, stored as row vectors
+#' @param COVARIANCE matrix or data frame (for diagonal covariance) or 3D array (for full covariance matrices)
 #' @param WEIGHTS vector containing the weights
 #' @return a list consisting of the log-likelihoods, cluster probabilities and cluster labels.
 #' @author Lampros Mouselimis
 #' @details
 #' This function takes the centroids, covariance matrix and weights from a trained model and returns the log-likelihoods, cluster probabilities and cluster labels for new data.
+#' The function handles both diagonal covariance matrices (2D matrix) and full covariance matrices (3D array/cube).
 #' @export
 #' @examples
 #'
@@ -156,18 +157,34 @@ predict_GMM = function(data, CENTROIDS, COVARIANCE, WEIGHTS) {
   if (!inherits(data, 'matrix')) stop('data should be either a matrix or a data frame')
   if ('data.frame' %in% class(CENTROIDS)) CENTROIDS = as.matrix(CENTROIDS)
   if (!inherits(CENTROIDS, 'matrix')) stop('CENTROIDS should be either a matrix or a data frame')
-  if ('data.frame' %in% class(COVARIANCE)) COVARIANCE = as.matrix(COVARIANCE)
-  if (!inherits(COVARIANCE, 'matrix')) stop('COVARIANCE should be either a matrix or a data frame')
-  if (ncol(data) != ncol(CENTROIDS) || ncol(data) != ncol(COVARIANCE) || length(WEIGHTS) != nrow(CENTROIDS) || length(WEIGHTS) != nrow(COVARIANCE))
-    stop('the number of columns of the data, CENTROIDS and COVARIANCE should match and the number of rows of the CENTROIDS AND COVARIANCE should be equal to the length of the WEIGHTS vector')
   if (!inherits(WEIGHTS, 'numeric') || !is.vector(WEIGHTS))
     stop('WEIGHTS should be a numeric vector')
 
   flag_non_finite = check_NaN_Inf(data)
 
   if (!flag_non_finite) stop("the data includes NaN's or +/- Inf values")
 
-  res = predict_MGausDPDF(data, CENTROIDS, COVARIANCE, WEIGHTS, eps = 1.0e-8)
+  # Check if COVARIANCE is a 3D array (full covariance) or 2D matrix (diagonal covariance)
+  is_full_covariance = length(dim(COVARIANCE)) == 3
+
+  if (is_full_covariance) {
+    # Full covariance matrices - COVARIANCE is a 3D array (cube)
+    if (!is.array(COVARIANCE)) stop('COVARIANCE should be a 3D array for full covariance matrices')
+    if (dim(COVARIANCE)[1] != ncol(data) || dim(COVARIANCE)[2] != ncol(data) || dim(COVARIANCE)[3] != length(WEIGHTS))
+      stop('for full covariance: dim(COVARIANCE)[1] and dim(COVARIANCE)[2] should equal ncol(data), and dim(COVARIANCE)[3] should equal length(WEIGHTS)')
+    if (ncol(data) != ncol(CENTROIDS) || length(WEIGHTS) != nrow(CENTROIDS))
+      stop('the number of columns of the data and CENTROIDS should match and the number of rows of CENTROIDS should equal the length of the WEIGHTS vector')
+
+    res = predict_MGausDPDF_full(data, CENTROIDS, COVARIANCE, WEIGHTS, eps = 1.0e-8)
+  } else {
+    # Diagonal covariance matrices - COVARIANCE is a 2D matrix
+    if ('data.frame' %in% class(COVARIANCE)) COVARIANCE = as.matrix(COVARIANCE)
+    if (!inherits(COVARIANCE, 'matrix')) stop('COVARIANCE should be either a matrix or a data frame for diagonal covariance')
+    if (ncol(data) != ncol(CENTROIDS) || ncol(data) != ncol(COVARIANCE) || length(WEIGHTS) != nrow(CENTROIDS) || length(WEIGHTS) != nrow(COVARIANCE))
+      stop('the number of columns of the data, CENTROIDS and COVARIANCE should match and the number of rows of the CENTROIDS AND COVARIANCE should be equal to the length of the WEIGHTS vector')
+
+    res = predict_MGausDPDF(data, CENTROIDS, COVARIANCE, WEIGHTS, eps = 1.0e-8)
+  }
 
   # I've added 1 to the output cluster labels to account for the difference in indexing between R and C++
   list(log_likelihood = res$Log_likelihood_raw,

inst/include/ClusterRHeader.h

@@ -1749,6 +1749,80 @@ namespace clustR {
 
 
 
+      // predict function for full covariance matrices (3D cube)
+      // This handles the case when GMM is fitted with full_covariance_matrices=TRUE
+      //
+
+      Rcpp::List predict_MGausDPDF_full(arma::mat data, arma::mat CENTROIDS, arma::cube COVARIANCE, arma::vec WEIGHTS, double eps = 1.0e-8) {
+
+        arma::mat gaus_mat(data.n_rows, WEIGHTS.n_elem, arma::fill::zeros);
+
+        arma::mat gaus_mat_log_lik(data.n_rows, WEIGHTS.n_elem, arma::fill::zeros);
+
+        for (unsigned int i = 0; i < WEIGHTS.n_elem; i++) {
+
+          arma::vec gaus_vec(data.n_rows, arma::fill::zeros);
+
+          arma::vec gaus_vec_log(data.n_rows, arma::fill::zeros);
+
+          // Extract the full covariance matrix for this Gaussian component
+          arma::mat tmp_cov_mt = COVARIANCE.slice(i);
+
+          double tmp_determinant = arma::det(tmp_cov_mt);
+
+          if (tmp_determinant == 0.0) {
+
+            Rcpp::stop("the determinant is zero or approximately zero. The data might include highly correlated variables or variables with low variance");
+          }
+
+          // Compute inverse of covariance matrix
+          arma::mat inv_cov_mt = arma::inv(tmp_cov_mt);
+
+          for (unsigned int j = 0; j < data.n_rows; j++) {
+
+            double n = data.n_cols;
+
+            arma::vec tmp_vec = (arma::conv_to< arma::vec >::from(data.row(j)) - arma::conv_to< arma::vec >::from(CENTROIDS.row(i)));
+
+            double tmp_val = 1.0 / std::sqrt(std::pow(2.0 * arma::datum::pi, n) * tmp_determinant);
+
+            double inner_likelih = 0.5 * (arma::as_scalar(tmp_vec.t() * inv_cov_mt * tmp_vec));
+
+            gaus_vec_log(j) = -(n / 2.0) * std::log(2.0 * arma::datum::pi) - (1.0 / 2.0) * (std::log(tmp_determinant)) - inner_likelih;
+
+            gaus_vec(j) = tmp_val * std::exp(-inner_likelih);
+          }
+
+          gaus_mat.col(i) = arma::as_scalar(WEIGHTS(i)) * gaus_vec;
+
+          gaus_mat_log_lik.col(i) = gaus_vec_log;
+        }
+
+        arma::mat loglik1(data.n_rows, WEIGHTS.n_elem, arma::fill::zeros);
+
+        arma::rowvec loglik2(data.n_rows, arma::fill::zeros);
+
+        for (unsigned int j = 0; j < loglik1.n_rows; j++) {
+
+          arma::rowvec tmp_vec = arma::conv_to< arma::rowvec >::from(gaus_mat.row(j)) + eps;
+
+          tmp_vec /= arma::sum(tmp_vec);                                                               // normalize row-data to get probabilities
+
+          loglik1.row(j) = tmp_vec;                                                                    // assign probabilities
+
+          arma::uvec log_lik_label = arma::find(tmp_vec == arma::max(tmp_vec));
+
+          loglik2(j) = arma::as_scalar(log_lik_label(0));                                              // assign labels
+        }
+
+        return Rcpp::List::create( Rcpp::Named("Log_likelihood_raw") = gaus_mat_log_lik,
+                                   Rcpp::Named("cluster_proba") = loglik1,
+                                   Rcpp::Named("cluster_labels") = loglik2 );
+      }
+
+
+
+
       // function to calculate bic-aic
       //
 

man/predict_GMM.Rd

@@ -199,6 +199,21 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// predict_MGausDPDF_full
+Rcpp::List predict_MGausDPDF_full(arma::mat data, arma::mat CENTROIDS, arma::cube COVARIANCE, arma::vec WEIGHTS, double eps);
+RcppExport SEXP _ClusterR_predict_MGausDPDF_full(SEXP dataSEXP, SEXP CENTROIDSSEXP, SEXP COVARIANCESEXP, SEXP WEIGHTSSEXP, SEXP epsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type data(dataSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type CENTROIDS(CENTROIDSSEXP);
+    Rcpp::traits::input_parameter< arma::cube >::type COVARIANCE(COVARIANCESEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type WEIGHTS(WEIGHTSSEXP);
+    Rcpp::traits::input_parameter< double >::type eps(epsSEXP);
+    rcpp_result_gen = Rcpp::wrap(predict_MGausDPDF_full(data, CENTROIDS, COVARIANCE, WEIGHTS, eps));
+    return rcpp_result_gen;
+END_RCPP
+}
 // GMM_arma_AIC_BIC
 arma::rowvec GMM_arma_AIC_BIC(arma::mat& data, arma::rowvec max_clusters, std::string dist_mode, std::string seed_mode, int km_iter, int em_iter, bool verbose, double var_floor, std::string criterion, int seed);
 RcppExport SEXP _ClusterR_GMM_arma_AIC_BIC(SEXP dataSEXP, SEXP max_clustersSEXP, SEXP dist_modeSEXP, SEXP seed_modeSEXP, SEXP km_iterSEXP, SEXP em_iterSEXP, SEXP verboseSEXP, SEXP var_floorSEXP, SEXP criterionSEXP, SEXP seedSEXP) {

src/RcppExports.cpp

@@ -294,6 +294,21 @@ Rcpp::List predict_MGausDPDF(arma::mat data, arma::mat CENTROIDS, arma::mat COVA
 
 
 
+// predict function for full covariance matrices (3D cube)
+// This handles the case when GMM is fitted with full_covariance_matrices=TRUE
+//
+
+// [[Rcpp::export]]
+Rcpp::List predict_MGausDPDF_full(arma::mat data, arma::mat CENTROIDS, arma::cube COVARIANCE, arma::vec WEIGHTS, double eps = 1.0e-8) {
+
+  ClustHeader CRH;
+
+  return CRH.predict_MGausDPDF_full(data, CENTROIDS, COVARIANCE, WEIGHTS, eps);
+}
+
+
+
+
 // function to calculate bic-aic
 //
 

src/export_inst_folder_headers.cpp

@@ -25,6 +25,7 @@ extern SEXP _ClusterR_opt_clust_fK(SEXP, SEXP, SEXP);
 extern SEXP _ClusterR_OptClust(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP _ClusterR_predict_medoids(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP _ClusterR_predict_MGausDPDF(SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP _ClusterR_predict_MGausDPDF_full(SEXP, SEXP, SEXP, SEXP, SEXP);
 extern SEXP _ClusterR_Predict_mini_batch_kmeans(SEXP, SEXP, SEXP, SEXP);
 extern SEXP _ClusterR_preferenceRange(SEXP, SEXP, SEXP);
 extern SEXP _ClusterR_SCALE(SEXP, SEXP, SEXP);
@@ -50,6 +51,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_ClusterR_OptClust",                    (DL_FUNC) &_ClusterR_OptClust,                    12},
     {"_ClusterR_predict_medoids",             (DL_FUNC) &_ClusterR_predict_medoids,              7},
     {"_ClusterR_predict_MGausDPDF",           (DL_FUNC) &_ClusterR_predict_MGausDPDF,            5},
+    {"_ClusterR_predict_MGausDPDF_full",      (DL_FUNC) &_ClusterR_predict_MGausDPDF_full,       5},
     {"_ClusterR_Predict_mini_batch_kmeans",   (DL_FUNC) &_ClusterR_Predict_mini_batch_kmeans,    4},
     {"_ClusterR_preferenceRange",             (DL_FUNC) &_ClusterR_preferenceRange,              3},
     {"_ClusterR_SCALE",                       (DL_FUNC) &_ClusterR_SCALE,                        3},