VCCRI / CIDR

Clustering through Imputation and Dimensionality Reduction
GNU General Public License v2.0
62 stars 22 forks source link

why cpp_dist and cpp_dist_weighted function can work? #8

Open eleozzr opened 7 years ago

eleozzr commented 7 years ago

When using the pakcage CIDR, I don't know why the functions cpp_dist and cpp_dist_weighted defined in RcppExports.cpp can work. There is only the statement of cpp_dist and cpp_dist_weighted, how the parameters in cpp_dist and cpp_dist_weighted, like dist, truth, counts, ncol et al dalculated.

The definitions of dist_cpp and dist_cpp_weighted function as following,

// cpp_dist

  | NumericMatrix cpp_dist(NumericMatrix dist, IntegerMatrix truth, NumericMatrix counts, int ncol, double threshold);   | RcppExport SEXP _cidr_cpp_dist(SEXP distSEXP, SEXP truthSEXP, SEXP countsSEXP, SEXP ncolSEXP, SEXP thresholdSEXP) {   | BEGIN_RCPP   | Rcpp::RObject rcpp_result_gen;   | Rcpp::RNGScope rcpp_rngScope_gen;   | Rcpp::traits::input_parameter< NumericMatrix >::type dist(distSEXP);   | Rcpp::traits::input_parameter< IntegerMatrix >::type truth(truthSEXP);   | Rcpp::traits::input_parameter< NumericMatrix >::type counts(countsSEXP);   | Rcpp::traits::input_parameter< int >::type ncol(ncolSEXP);   | Rcpp::traits::input_parameter< double >::type threshold(thresholdSEXP);   | rcpp_result_gen = Rcpp::wrap(cpp_dist(dist, truth, counts, ncol, threshold));   | return rcpp_result_gen;   | END_RCPP   | }   | // cpp_dist_weighted   | NumericMatrix cpp_dist_weighted(NumericMatrix dist, IntegerMatrix truth, NumericMatrix counts, int ncol, double a, double b);   | RcppExport SEXP _cidr_cpp_dist_weighted(SEXP distSEXP, SEXP truthSEXP, SEXP countsSEXP, SEXP ncolSEXP, SEXP aSEXP, SEXP bSEXP) {   | BEGIN_RCPP   | Rcpp::RObject rcpp_result_gen;   | Rcpp::RNGScope rcpp_rngScope_gen;   | Rcpp::traits::input_parameter< NumericMatrix >::type dist(distSEXP);   | Rcpp::traits::input_parameter< IntegerMatrix >::type truth(truthSEXP);   | Rcpp::traits::input_parameter< NumericMatrix >::type counts(countsSEXP);   | Rcpp::traits::input_parameter< int >::type ncol(ncolSEXP);   | Rcpp::traits::input_parameter< double >::type a(aSEXP);   | Rcpp::traits::input_parameter< double >::type b(bSEXP);   | rcpp_result_gen = Rcpp::wrap(cpp_dist_weighted(dist, truth, counts, ncol, a, b));   | return rcpp_result_gen;   | END_RCPP   | }

hojwk commented 7 years ago

CIDR utilises the package Rcpp, which automatically generates C & C++ code based on the Rcpp code (slightly modified C++ code) that we write.

This is the case when Rcpp generates _cpp_dist from our cpp_dist code.

For more information, you could look into the Rcpp documentation.

eleozzr commented 7 years ago

Thank you very much. But what's the your cpp_dist and cpp_dist_weighted code? I know the function eigen_centre and getEigenSpace have been defined explictly in the scr directory. But I still could not understand why the funtion cpp_dist and cpp_dist_weighted can be called directly in your scDissim function. Can you give me some notations?

hojwk commented 7 years ago

We can call the C++ functions (that use Rcpp libraries) cpp_dist & cpp_dist_weighted directly from the R function scDissim because we are using the R package Rcpp. This is the beauty of using Rcpp - it abstracts the need for us to know how this process works.

You may notice in the source code located in src/scPCA.cpp, which contains the functions cpp_dist & cpp_dist_weighted, there is a line in the code directly before each function: // [[Rcpp::export]] This is part of the required Rcpp procedure to allow the C++ function to become available in the R code.

eleozzr commented 7 years ago

Thank you very much. I guess I see now. It is very interesting and useful, which is a great work. By the way, I have another question. At the end of your code, thre is some lines to compute the optimal number of clusters, But I don't fully understand why you calculate like that.
image

hojwk commented 6 years ago

The typical shape for the CH index may have one peak that is clearly defined - see example(cidr) CH plot for one such example. In this case it may be reasonable to select the peak, which may correspond to the @Best.nc. However, we have found with a number of data sets, where we may know the real number of clusters, that this single peak shape does not hold. In some of these cases, the CH plot is monotinically decreasing, and @Best.nc is a lower number of clusters than we know to be true. Hence we empirically determined the displayed code via experimentation with a number of data sets, which is a form of the scree method (see CIDR paper for reference). The idea of the code is to find an "elbow" in the CH plot where the gradient may change significantly from the previous part of the plot. The user is encouraged to examine the CH plot, and the default number of clusters suggested by CIDR. If not happy with what CIDR produces, the user is able to override the value.