1
#include <Rcpp.h>
2
// [[Rcpp::plugins(cpp11)]]
3
// [[Rcpp::depends(RcppEigen)]]
4

5
using namespace Rcpp;
6

7
#include <vector>
8
#include <RcppEigen.h>
9
#include <Rcpp.h>
10

11
// [[Rcpp::export]]
12 4
double rcpp_proportion_held(Rcpp::NumericMatrix x, Rcpp::NumericMatrix y, Rcpp::NumericVector y_weights) {
13
  //// Initialization
14 4
  Eigen::MatrixXd x_MTX(Rcpp::as<Eigen::MatrixXd>(x));
15 4
  Eigen::MatrixXd y_MTX(Rcpp::as<Eigen::MatrixXd>(y));
16
    
17
  //// Preliminary processing
18
  // declare variables
19 4
  std::size_t N_x_INT = x_MTX.rows();
20 4
  std::size_t N_y_INT = y_MTX.rows(); 
21 4
  std::size_t N_k_INT = y_MTX.cols(); 
22 4
  Eigen::ArrayXd tmp_AXD;
23

24
  //// Main processing  
25
  /// calculate tss
26
  // compute centroid
27 4
  Eigen::RowVectorXd y_centroid_MTX(N_k_INT);
28 4
  y_centroid_MTX = y_MTX.colwise().sum();
29 4
  y_centroid_MTX /= static_cast<double>(N_y_INT);
30
  // tss
31 4
  double tss_DBL=0.0;
32 4
  for (std::size_t i=0; i<N_y_INT; ++i) {
33 4
    tmp_AXD = y_centroid_MTX - y_MTX.row(i);
34 4
    tss_DBL += (tmp_AXD.square().sum() * y_weights[i]);
35
  }
36
  
37
  /// calculate wss
38 4
  double wss_DBL = 0.0;
39
  double curr_lowest_DBL;
40 4
  for (std::size_t j=0; j < N_y_INT; ++j) {
41 4
    curr_lowest_DBL = std::numeric_limits<double>::max();
42 4
    for (std::size_t i=0; i < N_x_INT; ++i) {
43 4
      tmp_AXD = x_MTX.row(i) - y_MTX.row(j);
44 4
      curr_lowest_DBL = std::min(tmp_AXD.square().sum() * y_weights[i], curr_lowest_DBL);
45
    }
46 4
    wss_DBL += curr_lowest_DBL;
47
  }
48

49
  //// Exports
50
  // return proportion held
51 4
  return (1.0 - (wss_DBL / tss_DBL));
52
}

Read our documentation on viewing source code .

Loading