RcppbdCorr_matrix_cross

C++ Function Reference

1 Signature

corr_result BigDataStatMeth::RcppbdCorr_matrix_cross(const Eigen::MatrixXd &X, const Eigen::MatrixXd &Y, const std::string &method="pearson", bool use_complete_obs=true, bool compute_pvalues=true, bool trans_x=false, bool trans_y=false, Rcpp::Nullable< int > threads=R_NilValue)

2 Parameters

  • X (const Eigen::MatrixXd &)
  • Y (const Eigen::MatrixXd &)
  • method (const std::string &)
  • use_complete_obs (bool)
  • compute_pvalues (bool)
  • trans_x (bool)
  • trans_y (bool)
  • threads (Rcpp::Nullable< int >)

3 Returns

Type: corr_result

4 Call Graph

Function dependencies

5 Source Code

File: inst/include/hdf5Algebra/matrixCorrelation.hppLines 1044-1153

inline corr_result RcppbdCorr_matrix_cross(const Eigen::MatrixXd& X,
                                                         const Eigen::MatrixXd& Y,
                                                         const std::string& method = "pearson",
                                                         bool use_complete_obs = true,
                                                         bool compute_pvalues = true,
                                                         bool trans_x = false,
                                                         bool trans_y = false,
                                                         Rcpp::Nullable<int> threads = R_NilValue) {
        
        corr_result result;
        result.trans_x = trans_x;
        result.trans_y = trans_y;
        
        try {
            
            // Determine effective dimensions after transposition
            int n_obs_x, n_vars_x, n_obs_y, n_vars_y;
            
            if (trans_x) {
                n_obs_x = X.cols();  // Original variables become observations
                n_vars_x = X.rows(); // Original observations become variables
            } else {
                n_obs_x = X.rows();  // Standard layout
                n_vars_x = X.cols();
            }
            
            if (trans_y) {
                n_obs_y = Y.cols();  // Original variables become observations
                n_vars_y = Y.rows(); // Original observations become variables
            } else {
                n_obs_y = Y.rows();  // Standard layout
                n_vars_y = Y.cols();
            }
            
            // Validate dimensions after transposition
            if (n_obs_x != n_obs_y) {
                Rcpp::Rcerr << "Matrices must have the same number of observations after transposition" << std::endl;
                result.bcomputed = false;
                return result;
            }
            
            result.correlation_matrix = Eigen::MatrixXd::Zero(n_vars_x, n_vars_y);
            result.method = method;
            result.n_obs = n_obs_x;
            result.n_vars_x = n_vars_x;
            result.n_vars_y = n_vars_y;
            result.has_pvalues = compute_pvalues; 
            
            // Optimized threading for cross-correlation
            
    #ifdef _OPENMP
            int num_threads = 1;
            if (!threads.isNull()) {
                num_threads = Rcpp::as<int>(threads);
            } else {
                // Scale with matrix size but be conservative
                int total_pairs = n_vars_x * n_vars_y;
                num_threads = (total_pairs > 1000) ? std::min(6, omp_get_max_threads()) : 1;
            }
            num_threads = std::max(1, std::min(num_threads, omp_get_max_threads()));
    // #endif
    //         
    //         // Single level parallelization with good cache locality
    // #ifdef _OPENMP
    #pragma omp parallel for num_threads(num_threads) schedule(dynamic, 1)
    #endif
            for (int i = 0; i < n_vars_x; ++i) {
                for (int j = 0; j < n_vars_y; ++j) {
                    
                    Eigen::VectorXd vec_x, vec_y;
                    double corr_val;//, pval;
                    
                    // Get vectors with logical transposition
                    if (trans_x) {
                        vec_x = X.row(i).transpose(); // Row becomes variable
                    } else {
                        vec_x = X.col(i);            // Column is variable
                    }
                    
                    if (trans_y) {
                        vec_y = Y.row(j).transpose(); // Row becomes variable
                    } else {
                        vec_y = Y.col(j);            // Column is variable
                    }
                    
                    // double corr_val;
                    if (method == "spearman") {
                        corr_val = spearman_correlation(vec_x, vec_y, use_complete_obs);
                    } else {
                        corr_val = pearson_correlation(vec_x, vec_y, use_complete_obs);
                    }
                    
                    result.correlation_matrix(i, j) = corr_val;
                }
            }
            
            // Compute p-values using optimized vectorized function
            if (compute_pvalues) {
                result.pvalues = compute_pvalues_optimized(result.correlation_matrix, n_obs_x, false);
            }
            
            result.bcomputed = true;
            
        } catch (std::exception &ex) {
            Rcpp::Rcerr << "C++ exception RcppbdCorr_matrix_cross: " << ex.what() << std::endl;
            result.bcomputed = false;
        }
        
        return result;
    }

6 Usage Example

#include "BigDataStatMeth.hpp"

// Example usage
auto result = RcppbdCorr_matrix_cross(...);